Scipy Optimize Minimize For Combinatorics - Integer only solutions

Question

I have the below code that tries to assign the best product (a or b) to each factory (1,2,3,4 or 5).
How do I set constraints/bounds so that it only uses the integer values 0 or 1? e.g I have set x0 = [1,0,1,1,0] and the solution for this sample data is [0,0,0,0,0] Also, Is this even this correct approach to solve this problem? I have also tried pulp, but it does not seem to support the dynamic emission charges, as these are calculated in aggregate not at a factory level.

import numpy as np
from scipy.optimize import minimize

data = {
    '1': (1.2,),
    '2': (1.0,),
    '3': (1.7,),
    '4': (1.8,),
    '5': (1.6,)
}

unit_cost_b = 8
limit_b = 100
echarge_b = 0.004

unit_cost_a = 5
limit_a = 60
echarge_a = 0.007


def objective(x):
    x_a = [data[str(i+1)][0] for i in range(len(data)) if x[i] == 0]
    x_b = [data[str(i+1)][0] for i in range(len(data)) if x[i] == 1]

    total_unit_cost = (len(x_a) * unit_cost_a) + (len(x_b) * unit_cost_b)
    total_usage_a = np.sum(x_a)  
    total_usage_b = np.sum(x_b)
    
    emission_cost_a = max(total_usage_a - (len(x_a) * limit_a), 0) * echarge_a * 31
    emission_cost_b = max(total_usage_b - (len(x_b) * limit_b), 0) * echarge_b * 31

    return total_unit_cost + emission_cost_a + emission_cost_b

bounds = [(0, 1)] * len(data) 
eq_cons = {'type': 'eq', 'fun': lambda x: np.sum(x) - 5}


x0 = [1,0,1,1,0]
print(x0)
result = minimize(objective, x0, method='SLSQP', bounds=bounds, constraints=[eq_cons])

Answer 1

This is a pulp implementation below. You could do this in scipy using the scipy.optimize.MILP pkg, but I'm pretty sure you would need to develop the problem in matrix format, which is doable, but packages such as pulp or pyomo do it for you and I think they are much more approachable.

Realize, you cannot use if or max() in problem formulation as those things are nonlinear and the value of the variables is not known when the problem is handed over to a solver. And the bulk of what pulp is doing is making the problem to hand over....

On your model, the emissions data is odd and seems to always be zero because of the large values of limit, but perhaps it makes sense in larger context or you can adjust from the below.

CODE:

# factory picker

import pulp

data = {
    'f1': (1.2,),
    'f2': (1.0,),
    'f3': (1.7,),
    'f4': (1.8,),
    'f5': (1.6,)
}

# unit_cost_b = 8
# limit_b = 100
# echarge_b = 0.004

# unit_cost_a = 5
# limit_a = 60
# echarge_a = 0.007

unit_cost = {'a': 5, 'b': 8}
limit = {'a': 60, 'b': 100}
echarge = {'a': 0.007, 'b': 0.004}

# conveniences....
factories = list(data.keys())
products = ['a', 'b']
fp = [(f, p) for f in factories for p in products]

### PROB SETUP
prob = pulp.LpProblem('factory_assignments', pulp.LpMinimize)

### VARS

# make[factory, product] == 1 if making that product, else 0...
make = pulp.LpVariable.dicts('make', fp, cat=pulp.LpBinary)

emission = pulp.LpVariable.dicts('emission', products, lowBound=0)

### OBJECTIVE

# a convenience expression...
tot_unit_cost = sum(unit_cost[p] * make[f, p] for f, p in fp)

prob += tot_unit_cost + sum(emission[p] for p in products)


### CONSTRAINTS

# constrain the emission to GTE the expression.  It already has a lowBound=0, so that
# handles the max() expression in a linear expression
for p in products:
    prob += emission[p] >= sum(make[f, p] * (data[f][0] - limit[p]) for f in factories) * echarge[p] * 31

# each factory must produce *something*  (and only 1 thing)  (this is implied in the orig problem)
for f in factories:
    prob += sum(make[f, p] for p in products) == 1

### QA it...
print(prob)

### SOLVE
soln = prob.solve()

for f, p in fp:
    if make[f, p].varValue:   # True if the variable is 1
        print(f'make {p} in factory {f}')

print(f'tot unit cost: {pulp.value(tot_unit_cost)}')

for p in products:
    print(f'emission cost for {p} is: {emission[p].varValue}')

OUTPUT (truncated):

Result - Optimal solution found

Objective value:                25.00000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             0.00
Time (Wallclock seconds):       0.00

Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.00   (Wallclock seconds):       0.00

make a in factory f1
make a in factory f2
make a in factory f3
make a in factory f4
make a in factory f5
tot unit cost: 25.0
emission cost for a is: 0.0
emission cost for b is: 0.0

[Finished in 100ms]

Answer 2

The simplest approach would be to use brute force. If you can't express your objective as a linear function, that might be your best bet.

Example:

dims = 5
ranges = ((slice(0, 2, 1),) * dims)
x = scipy.optimize.brute(objective, ranges, finish=None)