Use this Skill for operations research and optimization: LP/QP/MILP with cvxpy and OR-Tools, convex relaxation, sensitivity analysis, and Gurobi/GLPK solver interface.
Use this Skill for operations research and optimization: LP/QP/MILP with cvxpy and OR-Tools, convex relaxation, sensitivity analysis, and Gurobi/GLPK solver interface.
Switch between open-source solvers (GLPK, ECOS, SCS) and commercial ones (Gurobi).
Model constraint satisfaction problems with OR-Tools CP-SAT.
Problem class
Recommended tool
LP (continuous)
cvxpy + GLPK / ECOS
QP (quadratic objective)
cvxpy + OSQP / ECOS
MILP (integer variables)
cvxpy + GLPK_MI, OR-Tools CP-SAT
Large-scale LP/MILP
OR-Tools linear solver (pywraplp)
Commercial (fastest)
Gurobi via cvxpy
Background & Key Concepts
Standard Form LP
A linear program in standard form:
minimize c^T x
subject to A_eq x = b_eq
A_ub x <= b_ub
x >= 0
The dual problem attaches a multiplier (shadow price) to each constraint.
Shadow price = marginal value of relaxing that constraint by one unit.
Mixed-Integer LP (MILP)
Replace some variables with cp.Variable(integer=True) or cp.Variable(boolean=True).
Solved by branch-and-bound: the LP relaxation at each node is solved and branched on
fractional integer variables.
Quadratic Programming (QP)
Markowitz portfolio optimization:
minimize (1/2) w^T Sigma w (portfolio variance)
subject to mu^T w >= r_min (return floor)
1^T w = 1 (fully invested)
w >= 0 (long-only)
Convex Relaxation
MILP is NP-hard in general. Convex (LP) relaxation of integer constraints provides
a lower bound and guides branch-and-bound. Tight relaxations lead to fast solves.
import cvxpy as cp
# List all solvers cvxpy can find on this machineprint("Available solvers:", cp.installed_solvers())
# Expected output includes: ['CLARABEL', 'ECOS', 'ECOS_BB', 'GLPK', 'GLPK_MI',# 'OSQP', 'SCS', 'SCIPY']# Gurobi will appear only if gurobipy is installed and licensed.
Core Workflow
Step 1 — Resource Allocation LP with cvxpy
Classic production planning: allocate limited resources across products to maximise profit.
import cvxpy as cp
import numpy as np
# ── Problem data ────────────────────────────────────────────────────────────────# 4 products, 3 resources
n_products = 4
n_resources = 3# Profit per unit of each product
profit = np.array([25.0, 30.0, 15.0, 20.0])
# Resource consumption matrix A[i, j] = units of resource i per unit of product j
A = np.array([
[1.0, 2.0, 1.0, 3.0], # Labour (hours)
[3.0, 1.0, 2.0, 1.0], # Material (kg)
[2.0, 2.0, 1.0, 2.0], # Machine time (hrs)
])
# Available resource capacities
b = np.array([240.0, 300.0, 200.0])
# ── Decision variables ──────────────────────────────────────────────────────────
x = cp.Variable(n_products, name="production")
# ── Objective ───────────────────────────────────────────────────────────────────
objective = cp.Maximize(profit @ x)
# ── Constraints ─────────────────────────────────────────────────────────────────
constraints = [
A @ x <= b, # Resource capacity
x >= 0, # Non-negativity
]
# ── Solve ────────────────────────────────────────────────────────────────────────
problem = cp.Problem(objective, constraints)
problem.solve(solver=cp.GLPK)
print(f"Status : {problem.status}")
print(f"Optimal profit: {problem.value:.2f}")
print("Production plan:")
for j, xj inenumerate(x.value):
print(f" Product {j+1}: {xj:.3f} units")
# ── Dual variables (shadow prices) ──────────────────────────────────────────────# Each shadow price = marginal profit gain per extra unit of that resource
shadow_prices = constraints[0].dual_value
resource_names = ["Labour", "Material", "Machine time"]
print("\nShadow prices (marginal value of relaxing each constraint):")
for name, lam inzip(resource_names, shadow_prices):
print(f" {name}: {lam:.4f}")
Step 2 — Facility Location MILP with OR-Tools CP-SAT
Decide which warehouses to open and which customers to assign to minimise total cost.
from ortools.sat.python import cp_model
import numpy as np
np.random.seed(0)
# ── Problem data ────────────────────────────────────────────────────────────────
n_facilities = 5# candidate warehouse sites
n_customers = 10# demand points# Fixed opening cost per facility
fixed_cost = np.array([100, 120, 80, 150, 90], dtype=int)
# Shipping cost[i][j] = cost to serve customer j from facility i
ship_cost = np.random.randint(5, 30, size=(n_facilities, n_customers))
# Demand (units) per customer
demand = np.random.randint(10, 50, size=n_customers)
# Capacity per facility
capacity = np.array([200, 180, 220, 160, 240], dtype=int)
# ── Scale to integers (CP-SAT requires integer coefficients) ────────────────────
SCALE = 1# already integer; increase if using fractional costs# ── Model ────────────────────────────────────────────────────────────────────────
model = cp_model.CpModel()
# y[i] = 1 if facility i is opened
y = [model.NewBoolVar(f"y_{i}") for i inrange(n_facilities)]
# x[i][j] = fraction of customer j's demand served from facility i (scaled to int)# Here we use full-integer assignment: x[i][j] in {0,1} (each customer assigned to 1 fac)
x = [[model.NewBoolVar(f"x_{i}_{j}") for j inrange(n_customers)]
for i inrange(n_facilities)]
# ── Constraints ──────────────────────────────────────────────────────────────────# 1. Each customer assigned to exactly one facilityfor j inrange(n_customers):
model.Add(sum(x[i][j] for i inrange(n_facilities)) == 1)
# 2. Customers can only be assigned to open facilitiesfor i inrange(n_facilities):
for j inrange(n_customers):
model.Add(x[i][j] <= y[i])
# 3. Capacity constraintsfor i inrange(n_facilities):
model.Add(
sum(x[i][j] * int(demand[j]) for j inrange(n_customers)) <= int(capacity[i])
)
# ── Objective: minimise fixed + shipping cost ─────────────────────────────────
total_cost = (
sum(int(fixed_cost[i]) * y[i] for i inrange(n_facilities))
+ sum(int(ship_cost[i][j]) * x[i][j]
for i inrange(n_facilities)
for j inrange(n_customers))
)
model.Minimize(total_cost)
# ── Solve ─────────────────────────────────────────────────────────────────────
solver = cp_model.CpSolver()
solver.parameters.max_time_in_seconds = 30.0
status = solver.Solve(model)
status_name = solver.StatusName(status)
print(f"Solver status : {status_name}")
if status in (cp_model.OPTIMAL, cp_model.FEASIBLE):
print(f"Total cost : {solver.ObjectiveValue():.0f}")
open_facs = [i for i inrange(n_facilities) if solver.Value(y[i]) == 1]
print(f"Open facilities: {open_facs}")
for j inrange(n_customers):
assigned = next(i for i inrange(n_facilities) if solver.Value(x[i][j]) == 1)
print(f" Customer {j:2d} -> Facility {assigned}")
Step 3 — Portfolio QP (Markowitz) with cvxpy
import cvxpy as cp
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
# ── Synthetic market data ────────────────────────────────────────────────────────
n_assets = 8
n_days = 252# one year of daily returns# Random expected returns and covariance matrix
mu_true = np.random.uniform(0.0005, 0.002, n_assets) # daily expected return# Build a random positive-definite covariance matrix
F = np.random.randn(n_assets, n_assets) * 0.01
Sigma = F @ F.T + np.diag(np.random.uniform(0.0001, 0.001, n_assets))
# ── Efficient frontier ──────────────────────────────────────────────────────────
w = cp.Variable(n_assets, name="weights")
r_targets = np.linspace(mu_true.min(), mu_true.max(), 40)
frontier_risk = []
frontier_return = []
frontier_weights = []
for r_target in r_targets:
objective = cp.Minimize(cp.quad_form(w, Sigma))
constraints = [
mu_true @ w >= r_target, # return floor
cp.sum(w) == 1, # fully invested
w >= 0, # long-only
]
prob = cp.Problem(objective, constraints)
prob.solve(solver=cp.ECOS, verbose=False)
if prob.status in ("optimal", "optimal_inaccurate") and w.value isnotNone:
frontier_risk.append(float(cp.sqrt(cp.quad_form(w, Sigma)).value))
frontier_return.append(float(mu_true @ w.value))
frontier_weights.append(w.value.copy())
# ── Maximum Sharpe ratio portfolio ──────────────────────────────────────────────
rf = 0.0001# daily risk-free rate
sharpe_ratios = [
(ret - rf) / risk if risk > 0else -np.inf
for ret, risk inzip(frontier_return, frontier_risk)
]
best_idx = int(np.argmax(sharpe_ratios))
print(f"Max-Sharpe portfolio:")
print(f" Daily return : {frontier_return[best_idx]*252:.4f} (annualised)")
print(f" Daily vol : {frontier_risk[best_idx]*np.sqrt(252):.4f} (annualised)")
print(f" Sharpe ratio : {sharpe_ratios[best_idx]*np.sqrt(252):.4f} (annualised)")
print(" Weights:", np.round(frontier_weights[best_idx], 4))
# ── Plot efficient frontier ──────────────────────────────────────────────────────
fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(
[r * np.sqrt(252) for r in frontier_risk],
[r * 252for r in frontier_return],
"b-o", markersize=3, label="Efficient frontier",
)
ax.scatter(
frontier_risk[best_idx] * np.sqrt(252),
frontier_return[best_idx] * 252,
color="red", zorder=5, s=80, label="Max Sharpe",
)
ax.set_xlabel("Annualised volatility")
ax.set_ylabel("Annualised expected return")
ax.set_title("Markowitz Efficient Frontier")
ax.legend()
fig.tight_layout()
fig.savefig("efficient_frontier.png", dpi=150)
print("Saved efficient_frontier.png")
Advanced Usage
Sensitivity Analysis with scipy.optimize.linprog
from scipy.optimize import linprog
import numpy as np
# Minimise c^T x# s.t. A_ub x <= b_ub# lb <= x <= ub
c = np.array([-25.0, -30.0, -15.0, -20.0]) # negate for maximisation
A_ub = np.array([
[1, 2, 1, 3],
[3, 1, 2, 1],
[2, 2, 1, 2],
])
b_ub = np.array([240.0, 300.0, 200.0])
bounds = [(0, None)] * 4
result = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method="highs")
print(f"scipy linprog optimal value: {-result.fun:.2f}")
print(f"Solution: {result.x}")
# Sensitivity information (available with HiGHS solver)ifhasattr(result, "ineqlin"):
print("\nShadow prices (HiGHS marginals):", result.ineqlin.marginals)
print("Slack values:", result.ineqlin.residual)
Solver Selection Guide and Gurobi Interface
import cvxpy as cp
import numpy as np
# ── Build a simple LP ────────────────────────────────────────────────────────────
n = 100
np.random.seed(7)
c = np.random.randn(n)
A = np.random.randn(50, n)
b = np.abs(A).sum(axis=1)
x = cp.Variable(n)
prob = cp.Problem(cp.Minimize(c @ x), [A @ x <= b, x >= 0])
# ── Try different solvers ────────────────────────────────────────────────────────
solver_map = {
"GLPK" : cp.GLPK,
"ECOS" : cp.ECOS,
"SCS" : cp.SCS,
"CLARABEL" : cp.CLARABEL,
# "GUROBI" : cp.GUROBI, # uncomment if Gurobi is installed and licensed
}
for name, solver in solver_map.items():
try:
prob.solve(solver=solver, warm_start=True)
status = prob.status
val = prob.value if prob.value isnotNoneelsefloat("nan")
print(f" {name:12s}: status={status:8s} value={val:.6f}")
except cp.SolverError as e:
print(f" {name:12s}: UNAVAILABLE ({e})")
