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 xsubject to:
b_l \leq A x \leq b_u, \quad l \leq x \leq u, \quad x_i \in \mathbb{Z} \; \forall i \in Iwhere 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.
Examples
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
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.
"""
def to2d(x):
"""Normalize scalar inputs (Excel passes [[x]] as x) to 2D lists."""
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 "Invalid input: c must be a 2D list of numeric values."
if c_vec.size == 0:
return "Invalid input: c must contain at least one coefficient."
n_vars = c_vec.size
# Helper to extract 1D arrays from potential 2D lists
def _to_array(values, name):
if values is None:
return None
try:
arr = np.array(values, dtype=float)
except Exception:
return f"Invalid input: {name} must contain numeric values."
return arr
# Integrality handling
integrality_arr = None
if integrality is not None:
try:
integrality_arr = np.array(integrality, dtype=int).flatten()
except Exception:
return "Invalid input: integrality must be integers."
if integrality_arr.size != n_vars:
return "Invalid input: integrality length must match number of variables."
if not np.all(np.isin(integrality_arr, [0, 1, 2, 3])):
return "Invalid input: integrality values must be 0 (continuous), 1 (integer), 2 (semi-continuous), or 3 (semi-integer)."
# Bounds handling
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 "Invalid input: bounds_lower must provide one value per variable."
if upper_arr is not None and upper_arr.size != n_vars:
return "Invalid input: 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 "Invalid input: 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:
# SciPy requires both vectors; if one missing, fill with defaults
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)
# Constraint processing helper
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"Invalid input: {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"Invalid input: {sense} matrix must have {n_vars} columns."
if vec.size != mat.shape[0]:
return f"Invalid input: {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
try:
result = scipy_milp(
c_vec,
integrality=integrality_arr,
bounds=bounds_obj,
constraints=constraints,
)
except Exception as exc:
return f"Error during MILP solving: {exc}"
if not result.success:
message = result.message if hasattr(result, "message") else "MILP solver did not succeed."
return f"MILP solving failed: {message}"
solution = result.x
if solution is None:
return "MILP solving failed: no solution returned."
try:
solution_list = list(map(float, solution))
except Exception:
return "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]]