32

I'm using scipy.optimize.minimize to optimize a real-world problem for which the answers can only be integers. My current code looks like this: 

from scipy.optimize import minimize

def f(x):
    return (481.79/(5+x[0]))+(412.04/(4+x[1]))+(365.54/(3+x[2]))+(375.88/(3+x[3]))+(379.75/(3+x[4]))+(632.92/(5+x[5]))+(127.89/(1+x[6]))+(835.71/(6+x[7]))+(200.21/(1+x[8]))

def con(x):
    return sum(x)-7

cons = {'type':'eq', 'fun': con}

print scipy.optimize.minimize(f, [1,1,1,1,1,1,1,0,0], constraints=cons, bounds=([0,7],[0,7],[0,7],[0,7],[0,7],[0,7],[0,7],[0,7],[0,7]))

This yields: 

x: array([  2.91950510e-16,   2.44504019e-01,   9.97850733e-01,
     1.05398840e+00,   1.07481251e+00,   2.60570253e-01,
     1.36470363e+00,   4.48527831e-02,   1.95871767e+00]

But I want it optimized with integer values (rounding all x to the nearest whole number doesn't always give the minimum).

Is there a way to use scipy.optimize.minimize with only integer values?

(I guess I could create an array with all possible permutations of x and evaluate f(x) for each combination, but that doesn't seem like a very elegant or quick solution.)

nicmet
  • 510
  • 1
  • 5
  • 9
  • 4
    That's not possible. There is no solver for **(Mixed-)Integer-Programming** available within numpy/scipy. You may want to use [pulp](https://github.com/coin-or/pulp) or some alternatives (pyomo, cvxpy, ...). Or if you are crazy: write your own branch-and-bound procedure. – sascha Aug 30 '16 at 20:56
  • Couldn't you round `x` to integer values within `f`, as well as round the optimised result to integers? – oulenz Feb 24 '20 at 10:55
  • [scipy.optimize import milp](https://scipy.github.io/devdocs/tutorial/optimize.html) – JeeyCi Aug 02 '23 at 04:30

4 Answers4

12

pulp solution

After some research, I don't think your objective function is linear. I recreated the problem in the Python pulp library but pulp doesn't like that we're dividing by a float and 'LpAffineExpression'. This answer suggests that linear programming "doesn't understand divisions" but that comment is in context of adding constraints, not the objective function. That comment pointed me to "Mixed Integer Linear Fractional Programming (MILFP)" and on Wikipedia.

Here's how you could do it in pulp if it actually worked (maybe someone can figure out why):

import pulp

data = [(481.79, 5), (412.04, 4), (365.54, 3)] #, (375.88, 3), (379.75, 3), (632.92, 5), (127.89, 1), (835.71, 6), (200.21, 1)]
x = pulp.LpVariable.dicts('x', range(len(data)), lowBound=0, upBound=7, cat=pulp.LpInteger)

numerator = dict((i,tup[0]) for i,tup in enumerate(data))
denom_int = dict((i,tup[1]) for i,tup in enumerate(data))

problem = pulp.LpProblem('Mixed Integer Linear Programming', sense=pulp.LpMinimize)

# objective function (doesn't work)
# TypeError: unsupported operand type(s) for /: 'float' and 'LpAffineExpression'
problem += sum([numerator[i] / (denom_int[i] + x[i]) for i in range(len(data))])

problem.solve()

for v in problem.variables():
  print(v.name, "=", v.varValue)

brute solution with scipy.optimize

You can use brute and ranges of slices for each x in your function. If you have 3 xs in your function, you'll also have 3 slices in your ranges tuple. The key to all of this is to add the step size of 1 to the slice(start, stop,step) so slice(#, #, 1).

from scipy.optimize import brute
import itertools

def f(x):
  return (481.79/(5+x[0]))+(412.04/(4+x[1]))+(365.54/(3+x[2]))

ranges = (slice(0, 9, 1),) * 3
result = brute(f, ranges, disp=True, finish=None)
print(result)

itertools solution

Or you can use itertools to generate all combinations:

combinations = list(itertools.product(*[[0,1,2,3,4,5,6,7,8]]*3))

values = []
for combination in combinations:
  values.append((combination, f(combination)))

best = [c for c,v in values if v == min([v for c,v in values])]
print(best)

Note: this is a scaled-down version of your original function for example purposes.

Jarad
  • 17,409
  • 19
  • 95
  • 154
  • 2
    Note that this is the brute-force try-every-possibility option already mentioned in the question. – user2357112 Oct 19 '17 at 22:40
  • 3
    The crux of the question is how to use something in `scipy.optimize` to return integer answers under a minimization strategy. Just because the concept was mentioned in the question doesn't imply that someone would know to use `brute` or think to use itertools - especially a beginner. It's also not obvious that that step size should be 1 initially. If you have a better answer, definitely post it! – Jarad Oct 19 '17 at 22:52
  • Thanks so much - will definitely look into pulp for future optimization problems, somehow didn't know it existed before this! – nicmet Nov 16 '17 at 20:26
8

One thing that might help your problem you could have a constraint as:

max([x-int(x)])=0

This is not going to completely solve your problem, the algorithm will still try and cheat and you will get values with some level of error ~±5e-10 that it will still try and optimize towards just by the error in scipy's algorithm but it's better than nothing.

cons = ({'type':'eq', 'fun': con},
        {'type':'eq','fun': lambda x : max([x[i]-int(x[i]) for i in range(len(x))])})

having tested this process on some optimizations I know the solution to, this process is more sensitive to the initial values than the unconstrained search, it gets fairly accurate answers however the solution may actually not find the true value, you are basically requiring the large jump of the optimization process (what it uses to make sure it's not optimizing to a local minimum) to search the sample space as the smaller increments are usually not strong enough to move to the next number over.

Simon Jones
  • 89
  • 1
  • 2
6

Here is a way to solve the Mixed Integer Nonlinear Programming problem with Python Gekko (a package that I maintain):

from gekko import GEKKO

m = GEKKO(remote=False)
x = m.Array(m.Var,9,lb=0,ub=7,integer=True)

def f(x):
    return (481.79/(5+x[0]))+(412.04/(4+x[1]))\
           +(365.54/(3+x[2]))+(375.88/(3+x[3]))\
           +(379.75/(3+x[4]))+(632.92/(5+x[5]))\
           +(127.89/(1+x[6]))+(835.71/(6+x[7]))\
           +(200.21/(1+x[8]))

m.Minimize(f(x))
m.Equation(sum(x)==7)
m.options.SOLVER=1
m.solve()
print(x)

This gives the solution:

 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :  0.0529 sec
 Objective      :  859.5269999999999
 Successful solution
 ---------------------------------------------------


[[0.0] [0.0] [1.0] [1.0] [1.0] [0.0] [1.0] [0.0] [3.0]]
John Hedengren
  • 12,068
  • 1
  • 21
  • 25
0

Generic function for bruteforce solution. Does somewhat a better job than brute in scipy, since scipy actually runs function with float numbers, not integers only, though range explicitly says so, as Jarad stated

def brute(func, arg_ranges, finish=min):
if isinstance(arg_ranges, dict):
    args = {k:np.unique(np.hstack([a for r in rs for a in r]) if isinstance(rs, list) else [a for a in rs]) for k,rs in arg_ranges.items()}
    print(args)
    return finish([(dict(zip(args.keys(), vs)), func(**dict(zip(args.keys(), vs)))) for vs in itertools.product(*args.values())], key=lambda x: x[1])
elif isinstance(arg_ranges, list):
    return finish([(i, func(i)) for r in arg_ranges for i in r], key=lambda x: x[1])
else:
    return finish([(i, func(i)) for i in arg_ranges], key=lambda x: x[1])

print(brute(lambda x,y: x / (y + 2), {'x':[range(1,5,2), range(0,6,1)], 'y':range(2,5,1)}))
Konstantin Samkowskiy
  • 21
  • 4