Binary and Integer Program in Python

I have a binary/integer program that I am trying to solve using milp in Python. Below the image and below is the code I have so far. When I run it, I get

`bounds.lb` and `bounds.ub` must contain reals and be broadcastable to `c.shape`.

I believed I accounted for the upper and lower bounds in the "bounds" variable. But apparently not. What is missing?

from scipy.optimize import linprog
import numpy as np
from scipy.optimize import milp, LinearConstraint, Bounds 
c = [3, 5, 7, 9, 2, 8, 6, 4, 1]  # x1 to x9 (office setup cost multipliers)
c = c + [0]*36     


A = [
    [1700, 3600, 2100, 2500, 3100, 2700, 4100, 3400, 4100] + [24.3]*9 + [32.4]*9 + [28.9]*9+[26]*9,  # Setup and transfer costs
    [0, 0, 0, -1, 0, -1, 1, 0, 0] + [0]*36,  # x7 <= x4 + x6
    [1, 0, 0, 0, 0, 1, 0, 0, 0] + [0]*36,  # x1 + x6 <= 1
    [0]*9 + [1]*9 + [0]*27,  # Chicago
    [0]*18 + [1]*9 + [0]*18,  # Charlotte
    [0]*27 + [1]*9 + [0]*9,  # Pittsburgh
    [0]*36 + [1]*9   # Houston
]

b = [
    14000,  
    0,  
    1,  
    29, 
    22, 
    19,  
    27,  
]

# Define bounds for each variable (0 <= xi <= 1 for binary variables, 0 <= yij)
x_bounds = [(0, 1)] * 9  # Bounds for the x variables
y_bounds = [(0, None)] * 36  # Bounds for the y variables (unbounded upper limit)

# Combine bounds
bounds = x_bounds + y_bounds

integrality = ([1]) * 45

milp(c, integrality = integrality,
    bounds = Bounds(bounds),
     constraints = LinearConstraint(A, b))

There are two issues with your code. First, to disable bounds or define "unbounded bounds", one should use np.inf or -np.inf instead of None.

Second, you are currently giving a list of tuples to scipy.optimize.Bounds which does not work. Instead you should provide a list of the lower bounds and a list for the upper bounds.

To resolve these issues, you can change the last part of your code to:

# Define bounds for each variable (0 <= xi <= 1 for binary variables, 0 <= yij)
x_bounds = [(0, 1)] * 9  # Bounds for the x variables
y_bounds = [(0, np.inf)] * 36  # Bounds for the y variables (unbounded upper limit)

# Combine bounds
bounds = x_bounds + y_bounds

integrality = ([1]) * 45

l_bounds = [x for (x,y) in bounds]
u_bounds = [y for (x,y) in bounds]

milp(c, integrality = integrality,
    bounds = Bounds(l_bounds, u_bounds),
     constraints = LinearConstraint(A, b))

Instead of the list comprehension you could also use the zip function for a more concise solution:

l_bounds, u_bounds = zip(*bounds)

or directly

bounds = Bounds(*zip(*bounds))

Personally, I prefer the solution using list comprehension as it is clearer to the reader what is happening and with lists of this length, performance is not really an issue.