MILP
Overview
The MILP function solves mixed-integer linear programming problems, a class of mathematical optimization where some or all decision variables are constrained to take integer values. MILP extends standard linear programming by adding integrality constraints, making it suitable for problems involving discrete decisions such as scheduling, resource allocation, network design, and capital budgeting.
This implementation uses scipy.optimize.milp, which is a wrapper around the HiGHS optimization solver—a high-performance, open-source solver developed at the University of Edinburgh. HiGHS employs a branch-and-cut algorithm for MIP problems, combining the dual revised simplex method with cutting planes and presolve techniques to efficiently find optimal or near-optimal solutions.
The function solves problems of the form:
\min_{x} \quad c^T x subject to: b_l \leq A x \leq b_u, \quad l \leq x \leq u, \quad x_i \in \mathbb{Z} \; \forall i \in I
where c is the vector of objective coefficients, A is the constraint matrix, b_l and b_u are constraint bounds, l and u are variable bounds, and I is the set of indices for integer-constrained variables.
The function supports four types of variable integrality: continuous (0), integer (1), semi-continuous (2), and semi-integer (3). Semi-continuous and semi-integer variables must either equal zero or fall within their specified bounds—useful for modeling minimum lot sizes or on/off decisions.
For more details on the underlying algorithm, see the SciPy MILP documentation and the HiGHS GitHub repository. For background on integer programming concepts, see Integer programming on Wikipedia.
This example function is provided as-is without any representation of accuracy.
Excel Usage
=MILP(c, integrality, bounds_lower, bounds_upper, a_ub, b_ub, a_eq, b_eq)c(list[list], required): Objective coefficients as a single row or column vector.integrality(list[list], optional, default: null): Variable integrality codes (0=continuous, 1=integer, 2=semi-continuous, 3=semi-integer).bounds_lower(list[list], optional, default: null): Lower bounds for each decision variable.bounds_upper(list[list], optional, default: null): Upper bounds for each decision variable.a_ub(list[list], optional, default: null): Inequality constraint coefficient matrix (A_ub x <= b_ub).b_ub(list[list], optional, default: null): Right-hand side values for inequality constraints.a_eq(list[list], optional, default: null): Equality constraint coefficient matrix (A_eq x = b_eq).b_eq(list[list], optional, default: null): Right-hand side values for equality constraints.
Returns (list[list]): 2D list [[x1, x2, …, objective]], or error message string.
Example 1: Integer variables with multiple inequality constraints
Inputs:
| c | integrality | bounds_lower | bounds_upper | a_ub | b_ub | |||||
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | -1 | 1 | 1 | 0 | 0 | 5 | 5 | -1 | 1 | 1 |
| 3 | 2 | 12 | ||||||||
| 2 | 3 | 12 |
Excel formula:
=MILP({0,-1}, {1,1}, {0,0}, {5,5}, {-1,1;3,2;2,3}, {1;12;12})Expected output:
| Result | ||
|---|---|---|
| 2 | 2 | -2 |
Example 2: Bounded integer optimization with production constraint
Inputs:
| c | integrality | bounds_lower | bounds_upper | a_ub | b_ub | |||||
|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 2 | 1 | 1 | 0 | 0 | 4 | 3 | -1 | -1 | -3 |
Excel formula:
=MILP({1,2}, {1,1}, {0,0}, {4,3}, {-1,-1}, -3)Expected output:
| Result | ||
|---|---|---|
| 3 | 0 | 3 |
Example 3: Mixed continuous and integer variables with equality constraint
Inputs:
| c | integrality | a_eq | b_eq | bounds_lower | bounds_upper | |||||
|---|---|---|---|---|---|---|---|---|---|---|
| 3 | 2 | 0 | 1 | 1 | 1 | 5 | 0 | 0 | 5 | 5 |
Excel formula:
=MILP({3,2}, {0,1}, {1,1}, 5, {0,0}, {5,5})Expected output:
| Result | ||
|---|---|---|
| 0 | 5 | 10 |
Example 4: Continuous LP with default bounds and inequality constraints
Inputs:
| c | a_ub | b_ub | ||||
|---|---|---|---|---|---|---|
| 1 | 1 | 1 | 1 | 2 | 0 | 4 |
| 0 | 1 | 1 | 3 |
Excel formula:
=MILP({1,1,1}, {1,2,0;0,1,1}, {4;3})Expected output:
| Result | |||
|---|---|---|---|
| 0 | 0 | 0 | 0 |
Python Code
Show Code
import math
import numpy as np
from scipy.optimize import Bounds, LinearConstraint
from scipy.optimize import milp as scipy_milp
def milp(c, integrality=None, bounds_lower=None, bounds_upper=None, a_ub=None, b_ub=None, a_eq=None, b_eq=None):
"""
Solve a mixed-integer linear program using scipy.optimize.milp.
See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.milp.html
This example function is provided as-is without any representation of accuracy.
Args:
c (list[list]): Objective coefficients as a single row or column vector.
integrality (list[list], optional): Variable integrality codes (0=continuous, 1=integer, 2=semi-continuous, 3=semi-integer). Default is None.
bounds_lower (list[list], optional): Lower bounds for each decision variable. Default is None.
bounds_upper (list[list], optional): Upper bounds for each decision variable. Default is None.
a_ub (list[list], optional): Inequality constraint coefficient matrix (A_ub x <= b_ub). Default is None.
b_ub (list[list], optional): Right-hand side values for inequality constraints. Default is None.
a_eq (list[list], optional): Equality constraint coefficient matrix (A_eq x = b_eq). Default is None.
b_eq (list[list], optional): Right-hand side values for equality constraints. Default is None.
Returns:
list[list]: 2D list [[x1, x2, ..., objective]], or error message string.
"""
try:
def to2d(x):
return [[x]] if not isinstance(x, list) else x
c = to2d(c)
if integrality is not None:
integrality = to2d(integrality)
if bounds_lower is not None:
bounds_lower = to2d(bounds_lower)
if bounds_upper is not None:
bounds_upper = to2d(bounds_upper)
if a_ub is not None:
a_ub = to2d(a_ub)
if b_ub is not None:
b_ub = to2d(b_ub)
if a_eq is not None:
a_eq = to2d(a_eq)
if b_eq is not None:
b_eq = to2d(b_eq)
try:
c_vec = np.array(c, dtype=float).flatten()
except Exception:
return "Error: c must be a 2D list of numeric values."
if c_vec.size == 0:
return "Error: c must contain at least one coefficient."
n_vars = c_vec.size
def _to_array(values, name):
if values is None:
return None
try:
arr = np.array(values, dtype=float)
except Exception:
return f"Error: {name} must contain numeric values."
return arr
integrality_arr = None
if integrality is not None:
try:
integrality_arr = np.array(integrality, dtype=int).flatten()
except Exception:
return "Error: integrality must be integers."
if integrality_arr.size != n_vars:
return "Error: integrality length must match number of variables."
if not np.all(np.isin(integrality_arr, [0, 1, 2, 3])):
return "Error: integrality values must be 0 (continuous), 1 (integer), 2 (semi-continuous), or 3 (semi-integer)."
lower_arr = None
upper_arr = None
if bounds_lower is not None:
lower_arr = _to_array(bounds_lower, "bounds_lower")
if isinstance(lower_arr, str):
return lower_arr
lower_arr = lower_arr.flatten()
if bounds_upper is not None:
upper_arr = _to_array(bounds_upper, "bounds_upper")
if isinstance(upper_arr, str):
return upper_arr
upper_arr = upper_arr.flatten()
if lower_arr is not None and lower_arr.size != n_vars:
return "Error: bounds_lower must provide one value per variable."
if upper_arr is not None and upper_arr.size != n_vars:
return "Error: bounds_upper must provide one value per variable."
if lower_arr is not None and upper_arr is not None:
if np.any(lower_arr > upper_arr):
return "Error: each lower bound must be less than or equal to the corresponding upper bound."
bounds_obj = None
if lower_arr is not None or upper_arr is not None:
if lower_arr is None:
lower_arr = np.zeros(n_vars)
if upper_arr is None:
upper_arr = np.full(n_vars, math.inf)
bounds_obj = Bounds(lower_arr, upper_arr)
def _build_linear_constraint(matrix, vector, sense):
if matrix is None and vector is None:
return None
if matrix is None or vector is None:
return f"Error: {sense} constraints require both matrix and vector."
mat = _to_array(matrix, f"{sense} matrix")
if isinstance(mat, str):
return mat
vec = _to_array(vector, f"{sense} vector")
if isinstance(vec, str):
return vec
vec = vec.flatten()
if mat.ndim == 1:
mat = mat.reshape(1, -1)
if mat.shape[1] != n_vars:
return f"Error: {sense} matrix must have {n_vars} columns."
if vec.size != mat.shape[0]:
return f"Error: {sense} vector length must equal number of rows in the matrix."
if sense == "inequality":
return LinearConstraint(mat, -np.inf, vec)
return LinearConstraint(mat, vec, vec)
constraints = []
for sense, matrix, vector in (
("inequality", a_ub, b_ub),
("equality", a_eq, b_eq),
):
constraint = _build_linear_constraint(matrix, vector, sense)
if isinstance(constraint, str):
return constraint
if constraint is not None:
constraints.append(constraint)
if len(constraints) == 0:
constraints = None
result = scipy_milp(
c_vec,
integrality=integrality_arr,
bounds=bounds_obj,
constraints=constraints,
)
if not result.success:
message = result.message if hasattr(result, "message") else "MILP solver did not succeed."
return f"Error: {message}"
solution = result.x
if solution is None:
return "Error: No solution returned."
try:
solution_list = list(map(float, solution))
except Exception:
return "Error: Error converting MILP solution to floats."
objective = float(result.fun) if result.fun is not None else float(np.dot(c_vec, solution))
return [solution_list + [objective]]
except Exception as e:
return f"Error: {str(e)}"