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.
Install with Codex or Claude Copy this prompt, paste it into Codex, Claude, or another assistant, and let it review the skill page and install it for you.
A direct command skips the review prompt. Inspect the source before running it.
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})")
