13. Computational Tools and Practice

Every model in this document ends as input to software. A modeling layer turns the algebra and the data into the numbers a solver understands, and the solver's answer has to be read, checked and trusted before anyone acts on it. This chapter covers both halves of that work: the tools, which are modeling languages, Python libraries, commercial and open-source solvers, and spreadsheets (13.1), and the practical craft of making them work on real problems, which means numerics, solver logs, infeasible models, performance tuning and deployed applications (13.2).

<a id="131-modeling-languages-and-solvers"></a>

13.1 Modeling Languages and Solvers

A modern optimization system separates three things: the model (the algebra), the data (sets and numbers), and the solver (the algorithm). The modeling layer combines model and data into a concrete instance, passes it to any solver, and maps the answer back to meaningful names:

   model (algebra)          data (sets, parameters)
          \                        /
           v                      v
     +--------------------------------------+
     |  modeling layer: instance generation |   AMPL, GAMS, Pyomo, PuLP, OR-Tools, gurobipy
     +--------------------------------------+
                       |   c, A, b, bounds, integrality (in memory, or as an MPS / LP file)
                       v
     +--------------------------------------+
     |  solver                              |   HiGHS, CBC, SCIP, Gurobi, CPLEX, Xpress
     +--------------------------------------+
                       |   status, x*, duals, reduced costs, bound and gap
                       v
       results mapped back to names, e.g. Ship['seattle', 'chicago'] = 300

The running example of 13.1 is the classic transportation problem trnsport from the GAMS model library, which goes back to Dantzig's 1963 textbook. Two canning plants (Seattle with 350 cases, San Diego with 600) supply three markets (New York 325, Chicago 300, Topeka 275). Freight costs 90 dollars per case per thousand miles. The known optimal cost of 153.675 thousand dollars makes it easy to check every tool against the same answer. Later topics use a harder capacitated facility location model.

Algebraic Modeling Languages (AMPL, GAMS)

Why it matters. Before modeling languages existed, practitioners wrote matrix generators: custom Fortran programs that emitted the constraint matrix column by column. The model existed only implicitly in that code, so a wrong coefficient in row 4,217 was practically invisible. Algebraic modeling languages (AMLs) let you write the model the way you would write it on paper, with indexed sets and summations, and generate the matrix automatically. The same 15 lines describe an instance with 6 variables or 6 million; only the data file changes. AMLs remain the standard in energy-system planning, economics, chemical engineering and many industrial OR teams, and they shaped the design of every Python library in the next topic.

Theory. The transportation model in mathematical form is

$$ \min \sum_{i \in I} \sum_{j \in J} c_{ij} x_{ij} \quad \text{s.t.} \quad \sum_{j \in J} x_{ij} \le a_i \ \ \forall i \in I, \qquad \sum_{i \in I} x_{ij} \ge b_j \ \ \forall j \in J, \qquad x_{ij} \ge 0, $$

with $c_{ij} = f , d_{ij} / 1000$. An AML has one construct for each element of this notation:

Concept	Mathematics	AMPL	GAMS
Set	$i \in I$	set PLANTS;	Set i / seattle, san-diego /;
Parameter	$a_i$	param supply {PLANTS};	Parameter a(i);
Variable	$x_{ij} \ge 0$	var Ship {PLANTS, MARKETS} >= 0;	Positive Variable x(i,j);
Constraint family	$\sum_j x_{ij} \le a_i \ \forall i$	subject to Supply {i in PLANTS}: ...;	supply(i).. sum(j, x(i,j)) =l= a(i);
Objective	$\min \sum c_{ij} x_{ij}$	minimize Total_Cost: ...;	cost.. z =e= ...; and minimizing z
Solve		solve;	solve transport using lp minimizing z;

The key ideas are:

• Declarative modeling. You state what the model is, not how to build the matrix. A constraint family such as Supply {i in PLANTS} is a single statement that generates one row per plant.
• Model/data separation. The model file contains no numbers. The same model runs on test data, last month's data and a what-if scenario.
• Solver independence. Changing option solver highs to option solver gurobi is a one-line change. The AML writes the instance in a solver-neutral form and translates the answer back.
• Automatic differentiation. For nonlinear models, the AML computes exact gradients and Hessians and passes them to NLP solvers such as Ipopt or Knitro. This is a major advantage over hand-coded derivatives.
• Presolve and checking. AMLs check the model against its declarations (an index outside a set, a negative value for a parameter declared >= 0) before any solver runs.

The two classic AMLs differ in style more than in power:

	AMPL	GAMS
Origin	Bell Labs (Fourer, Gay, Kernighan), late 1980s	World Bank, 1970s and 1980s
Syntax	close to mathematical notation	compact, equation-oriented (=l=, =g=, =e=)
Objective	named minimize or maximize statement	a free variable z defined by an equation
Model type	inferred from the model	declared in the solve statement (lp, mip, nlp, minlp, ...)
Data	separate .dat files, databases, spreadsheets	inline tables, or GDX binary files
Python access	amplpy	GAMS Python API, GAMSPy
Typical users	industrial OR, teaching	energy and economic models, process engineering

Other AMLs follow the same pattern: AIMMS, LINGO, Xpress Mosel, IBM's OPL, and the open-source ZIMPL and GNU MathProg (a subset of AMPL used by GLPK). Both AMPL and GAMS are commercial, with free size-limited editions for learning.

Example: the same model in AMPL and GAMS. In AMPL, the model, the data and the commands live in three files:

# transport.mod: the model, with no numbers in it
set PLANTS;
set MARKETS;

param supply {PLANTS} >= 0;                 # cases available at each plant
param demand {MARKETS} >= 0;                # cases required at each market
param dist {PLANTS, MARKETS} >= 0;          # thousands of miles
param freight >= 0;                         # dollars per case per thousand miles
param cost {i in PLANTS, j in MARKETS} := freight * dist[i,j] / 1000;   # thousands of $ per case

var Ship {PLANTS, MARKETS} >= 0;            # cases shipped from i to j

minimize Total_Cost: sum {i in PLANTS, j in MARKETS} cost[i,j] * Ship[i,j];

subject to Supply {i in PLANTS}:  sum {j in MARKETS} Ship[i,j] <= supply[i];
subject to Demand {j in MARKETS}: sum {i in PLANTS}  Ship[i,j] >= demand[j];

# transport.dat: one data instance
set PLANTS  := seattle san_diego;
set MARKETS := new_york chicago topeka;

param supply := seattle 350  san_diego 600;
param demand := new_york 325  chicago 300  topeka 275;

param dist:     new_york  chicago  topeka :=
  seattle          2.5      1.7      1.8
  san_diego        2.5      1.8      1.4 ;

param freight := 90;

# transport.run: the session
model transport.mod;
data transport.dat;
option solver highs;          # or gurobi, cplex, cbc, ...: the model does not change
solve;
display Total_Cost, Ship;
display Demand.dual;          # shadow price of each market's demand

The GAMS version is the original trnsport model, slightly abridged:

Set
   i 'canning plants' / seattle, san-diego /
   j 'markets'        / new-york, chicago, topeka /;

Parameter
   a(i) 'capacity of plant i in cases'
        / seattle 350, san-diego 600 /
   b(j) 'demand at market j in cases'
        / new-york 325, chicago 300, topeka 275 /;

Table d(i,j) 'distance in thousands of miles'
              new-york  chicago  topeka
   seattle       2.5      1.7     1.8
   san-diego     2.5      1.8     1.4;

Scalar f 'freight in dollars per case per thousand miles' / 90 /;

Parameter c(i,j) 'transport cost in thousands of dollars per case';
c(i,j) = f*d(i,j)/1000;

Variable          z      'total transportation cost in thousands of dollars';
Positive Variable x(i,j) 'shipment quantities in cases';

Equation
   cost      'define objective function'
   supply(i) 'observe supply limit at plant i'
   demand(j) 'satisfy demand at market j';

cost..      z =e= sum((i,j), c(i,j)*x(i,j));
supply(i).. sum(j, x(i,j)) =l= a(i);
demand(j).. sum(i, x(i,j)) =g= b(j);

Model transport / all /;
solve transport using lp minimizing z;
* levels (.l) and marginals (.m) of the shipments
display x.l, x.m;

What does the AML do between solve and the solver? The following Python script performs the same instance generation by hand. It expands the indexed constraints into the matrix that the solver actually receives:

import numpy as np
from scipy.optimize import linprog

# The "data section": sets and parameters, exactly as an AML data file would hold them
plants = ["seattle", "san-diego"]
markets = ["new-york", "chicago", "topeka"]
supply = {"seattle": 350, "san-diego": 600}
demand = {"new-york": 325, "chicago": 300, "topeka": 275}
dist = {("seattle", "new-york"): 2.5, ("seattle", "chicago"): 1.7, ("seattle", "topeka"): 1.8,
        ("san-diego", "new-york"): 2.5, ("san-diego", "chicago"): 1.8, ("san-diego", "topeka"): 1.4}
freight = 90                                            # $ per case per thousand miles

# "Model generation": expand the indexed model into c, A_ub, b_ub, the step an AML automates
arcs = [(i, j) for i in plants for j in markets]        # one column per (i, j) pair
col = {a: k for k, a in enumerate(arcs)}
c = np.array([freight * dist[a] / 1000 for a in arcs])  # cost in thousands of $ per case
A, b, names = [], [], []
for i in plants:                                        # supply[i]:  sum_j x[i,j] <= a[i]
    row = np.zeros(len(arcs)); row[[col[i, j] for j in markets]] = 1
    A.append(row); b.append(supply[i]); names.append(f"supply[{i}]")
for j in markets:                                       # demand[j]:  sum_i x[i,j] >= b[j]
    row = np.zeros(len(arcs)); row[[col[i, j] for i in plants]] = -1   # written as <= with sign flip
    A.append(row); b.append(-demand[j]); names.append(f"demand[{j}]")
A = np.array(A)

print("generated matrix (rows = constraints, columns = x[i,j]):")
for n, row, rhs in zip(names, A, b):
    print(f"  {n:20s} {' '.join(f'{v:+.0f}' for v in row)}  <= {rhs}")
res = linprog(c, A_ub=A, b_ub=b, method="highs")       # the solver only ever sees numbers
print(f"z = {res.fun:.3f} thousand $")
for a, v in zip(arcs, res.x):
    if v > 1e-9:
        print(f"  x[{a[0]},{a[1]}] = {v:.0f}")
# generated matrix (rows = constraints, columns = x[i,j]):
#   supply[seattle]      +1 +1 +1 +0 +0 +0  <= 350
#   supply[san-diego]    +0 +0 +0 +1 +1 +1  <= 600
#   demand[new-york]     -1 +0 +0 -1 +0 +0  <= -325
#   demand[chicago]      +0 -1 +0 +0 -1 +0  <= -300
#   demand[topeka]       +0 +0 -1 +0 +0 -1  <= -275
# z = 153.675 thousand $
#   x[seattle,chicago] = 300
#   x[san-diego,new-york] = 325
#   x[san-diego,topeka] = 275

The solver receives a $5 \times 6$ matrix and three vectors. Names, sets and units exist only in the modeling layer, which is why it must map the solution back. The optimal cost matches the published 153.675. The shipping plan, however, differs from the one printed in the GAMS tutorial, which sends 50 cases from Seattle to New York. Both plans cost the same: the problem has alternative optima, and different solvers stop at different optimal vertices (more in Section 13.2). Doing this expansion by hand for a model with twenty index sets is exactly the error-prone work that AMLs eliminate.

Key takeaways

• An AML states the model algebraically, with sets, parameters, variables, constraint families and an objective, and generates the solver's matrix automatically.
• Model/data separation and solver independence let one model serve many instances and many solvers.
• AMPL reads like mathematics, while GAMS is equation-oriented and declares the model type. Both provide automatic derivatives for nonlinear models.
• Solvers see only numbers; the modeling layer is responsible for names, units and mapping results back.

Practice

1. (Beginner) Each plant now has a production cost per case (Seattle 0.10, San Diego 0.12 thousand dollars). Write the new parameter and the modified objective in AMPL. Which of the three files change, and which stay the same?
2. (Intermediate) Add a warehouse in Denver through which goods may be transshipped, with distances of your choice. Rewrite the AMPL model as a minimum-cost flow over a set of NODES and a set of ARCS within {NODES, NODES} with flow-balance constraints (Section 3.5), so that the network can change without touching the model.
3. (Interview) "Our team is fluent in Python. Why would anyone keep a production model in AMPL or GAMS?" Hint: think about readability for non-programmers, automatic differentiation, decades of validated models, generation speed for very large models, and the existence of Python APIs (amplpy, GAMSPy) that combine both worlds.

Python Modeling Libraries (Pyomo, PuLP, OR-Tools)

Why it matters. Most optimization today runs inside larger software systems: a model reads data from a database, solves, and writes a schedule back through an API, often in a container that runs every night. Python modeling libraries embed the model in a general-purpose language, so data handling, solving and integration happen in one program, with no licence for the modeling layer. Knowing their differences saves a painful rewrite when a project outgrows its first choice.

Theory. All Python libraries follow the same five steps: create a model, create variables indexed by Python collections, add constraints, set the objective, call a solver and read the values. They rely on operator overloading. In x + y <= 3, the <= does not compute a boolean; it builds a constraint object that the library stores. This is the same trick SQL query builders use, where Python expressions describe a query that an engine later executes.

Library	Problem classes	Solvers	Strengths	Limitations
PuLP (COIN-OR)	LP, MILP	CBC, HiGHS, SCIP, Gurobi, CPLEX, Xpress, ...	very small API, quick to learn	linear models only, few advanced features
Pyomo (Sandia National Laboratories)	LP, MILP, NLP, MINLP, disjunctive, dynamic (DAE), stochastic	dozens through plugins (HiGHS, Ipopt, Gurobi, CPLEX, ...)	the most general, many extensions	more verbose, slower model building
OR-Tools (Google)	LP, MILP, constraint programming (CP-SAT), vehicle routing, graph algorithms	GLOP, PDLP, CP-SAT, plus CBC, SCIP, HiGHS, Gurobi through wrappers	best-in-class CP-SAT and routing library	no general nonlinear modeling

Other options exist for specific needs. Vendor APIs (gurobipy, docplex) give full access to one solver. CVXPY handles convex problems with automatic convexity checks. linopy builds very large models with vectorized arrays, and JuMP is the equivalent ecosystem in Julia.

Three practical points apply to all of them:

• Tolerances. A value that should be 0 may come back as 6e-14. Compare against a threshold such as > 1e-6, never == 0.
• Names. Libraries sanitize names for the solver files: PuLP turns demand_new-york into demand_new_york. Look constraints up by their sanitized names.
• Versions. APIs change. PuLP 4.0 replaced pulp.LpVariable.dicts(...) with methods on the problem (prob.add_variable_dicts, prob.add_variable_dict) and returns a statistics object from solve. Pin library versions in production.

Example: one model, three libraries. All three scripts read the same data module:

# transport_data.py
plants = ["seattle", "san-diego"]
markets = ["new-york", "chicago", "topeka"]
supply = {"seattle": 350, "san-diego": 600}
demand = {"new-york": 325, "chicago": 300, "topeka": 275}
dist = {("seattle", "new-york"): 2.5, ("seattle", "chicago"): 1.7, ("seattle", "topeka"): 1.8,
        ("san-diego", "new-york"): 2.5, ("san-diego", "chicago"): 1.8, ("san-diego", "topeka"): 1.4}
cost = {a: 90 * d / 1000 for a, d in dist.items()}      # thousands of $ per case

PuLP (version 4.0) with its default solver, CBC:

import pulp
from transport_data import plants, markets, supply, demand, cost

m = pulp.LpProblem("transport", pulp.LpMinimize)
x = m.add_variable_dicts("x", (plants, markets), lowBound=0)        # PuLP < 4: pulp.LpVariable.dicts(...)
m += pulp.lpSum(cost[i, j] * x[i][j] for i in plants for j in markets)   # first expression = objective
for i in plants:
    m += pulp.lpSum(x[i][j] for j in markets) <= supply[i], f"supply_{i}"
for j in markets:
    m += pulp.lpSum(x[i][j] for i in plants) >= demand[j], f"demand_{j}"

stats = m.solve(pulp.COIN_CMD(msg=False))                                # CBC; or pulp.HiGHS(), pulp.GUROBI()
print(stats.status_str, round(pulp.value(m.objective), 3))
print({f"{i}->{j}": x[i][j].value() for i in plants for j in markets if x[i][j].value() > 1e-6})
print("dual of demand_new-york:", m.get_constraint_by_name("demand_new_york").pi)
# Optimal 153.675
# {'seattle->new-york': 50.0, 'seattle->chicago': 300.0, 'san-diego->new-york': 275.0, 'san-diego->topeka': 275.0}
# dual of demand_new-york: 0.225

Pyomo, solved with HiGHS:

import pyomo.environ as pyo
from transport_data import plants, markets, supply, demand, cost

m = pyo.ConcreteModel()
m.I, m.J = pyo.Set(initialize=plants), pyo.Set(initialize=markets)
m.x = pyo.Var(m.I, m.J, domain=pyo.NonNegativeReals)
m.cost = pyo.Objective(expr=sum(cost[i, j] * m.x[i, j] for i in m.I for j in m.J))
m.supply = pyo.Constraint(m.I, rule=lambda m, i: sum(m.x[i, j] for j in m.J) <= supply[i])
m.demand = pyo.Constraint(m.J, rule=lambda m, j: sum(m.x[i, j] for i in m.I) >= demand[j])
m.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT)       # ask the solver to return duals

res = pyo.SolverFactory("appsi_highs").solve(m)         # or "gurobi_direct", "cplex", "cbc", ...
print(res.solver.termination_condition, round(pyo.value(m.cost), 3))
print({f"{i}->{j}": m.x[i, j].value for i in m.I for j in m.J if m.x[i, j].value > 0})
print("dual of demand[new-york]:", m.dual[m.demand["new-york"]])
# optimal 153.675
# {'seattle->chicago': 300.0, 'san-diego->new-york': 325.0, 'san-diego->topeka': 275.0}
# dual of demand[new-york]: 0.225

OR-Tools, solved with GLOP, Google's simplex code:

from ortools.linear_solver import pywraplp
from transport_data import plants, markets, supply, demand, cost

s = pywraplp.Solver.CreateSolver("GLOP")                # Google's LP simplex; "CBC", "SCIP", "CP_SAT" for MIPs
x = {(i, j): s.NumVar(0, s.infinity(), f"x[{i},{j}]") for i in plants for j in markets}
s.Minimize(sum(cost[a] * x[a] for a in x))
for i in plants:
    s.Add(sum(x[i, j] for j in markets) <= supply[i], f"supply[{i}]")
dem = {j: s.Add(sum(x[i, j] for i in plants) >= demand[j], f"demand[{j}]") for j in markets}

status = s.Solve()
print(status == pywraplp.Solver.OPTIMAL, round(s.Objective().Value(), 3))
print({f"{i}->{j}": v.solution_value() for (i, j), v in x.items() if v.solution_value() > 1e-6})
print("dual of demand[new-york]:", dem["new-york"].dual_value())
# True 153.675
# {'seattle->chicago': 300.0, 'san-diego->new-york': 324.99999999999994, 'san-diego->topeka': 275.0}
# dual of demand[new-york]: 0.225

All three report the optimal cost of 153.675 and the same shadow price for New York: 0.225 thousand dollars, so each extra case demanded there costs 225 dollars, exactly the San Diego to New York freight rate. CBC returns the GAMS tutorial plan (50 cases from Seattle to New York), while HiGHS and GLOP return the other optimal vertex. GLOP's 324.99999999999994 is floating-point arithmetic at work, and without the > 1e-6 filter GLOP also prints a phantom Seattle to New York shipment of 6.1e-14 cases. The model code is almost identical across libraries. The real differences lie in scope (PuLP stops at MILP), in available solvers, and in specialized engines such as CP-SAT and the routing library.

Key takeaways

• Python libraries build models through operator overloading and hand them to interchangeable solvers.
• Choose PuLP for simple LP/MILP work, Pyomo for nonlinear or structured models, and OR-Tools for scheduling, routing and combinatorial problems through CP-SAT.
• Compare solver values against tolerances, look constraints up by their sanitized names, and pin library versions.
• Different solvers may return different optimal solutions with the same objective. Test the objective, not specific variable values.

Practice

1. (Beginner) Rewrite the PuLP model so that the data come from two CSV files, plants.csv (name, supply) and routes.csv (plant, market, distance), loaded with pandas. Which lines of the model change?
2. (Intermediate) Add a fixed charge of 10 thousand dollars for every route that is used: a binary $u_{ij}$ with $x_{ij} \le \min(a_i, b_j),u_{ij}$. Solve it with PuLP and CBC, and with OR-Tools using CreateSolver("SCIP"). How many routes are used now, and what is the new cost?
3. (Interview) "We prototyped our pricing model in PuLP. The next version needs a nonlinear demand curve and must run on two different solvers. What do you recommend?" Hint: consider which libraries support NLP and MINLP, how solver plugins work in Pyomo, whether the nonlinearity can be linearized piecewise (Section 4.1), and the cost of migrating.

Commercial Solvers (Gurobi, CPLEX)

Why it matters. On hard mixed-integer programs, the gap between the best commercial solvers and most alternatives is often one or two orders of magnitude in run time, and sometimes the difference between solving and not solving. In the benchmark of the next topic, Gurobi is 6 to 10 times faster than the best open-source solver on the same models, and CBC cannot solve two of the three instances within 5 minutes. Speed has direct business value: a plan that takes 5 minutes instead of 5 hours can be re-optimized during the day. Commercial solvers also provide robust numerics, support contracts and tools for tuning, infeasibility analysis and distributed solving.

Theory. A modern MIP solver is a pipeline of the techniques from Chapters 2 and 4, refined by decades of engineering:

  presolve ──> root LP ──> cutting planes ──> primal heuristics ──> branch and bound
  (remove,     (dual       (Gomory, MIR,      (rounding, RINS,      (branching rules,
   tighten,     simplex,    cover, flow,       local search)         node selection,
   detect       barrier,    implied bound,                           conflict analysis,
   structure)   concurrent) clique, ...)                             symmetry, parallel tree)

Each component matters: presolve often removes a large share of rows and columns, cuts close much of the root gap, and heuristics find good incumbents long before the tree search proves optimality. The leading products are:

Solver	Origin	Problem classes	Notes
Gurobi	founded 2008 by Gu, Rothberg and Bixby	LP, QP, SOCP, MILP, MIQP, nonconvex quadratic, general nonlinear in recent versions	often the fastest MIP solver; Python API gurobipy
IBM CPLEX	Bixby, 1988; ILOG, then IBM	LP, QP, SOCP, MILP, MIQP	mature and very robust; includes CP Optimizer for scheduling
FICO Xpress	Dash Optimization, then FICO	LP, MILP, MINLP	strong in finance and supply chain
MOSEK	MOSEK ApS	LP, conic (SOCP, SDP), MILP	the reference for conic optimization
COPT, Hexaly	Cardinal Operations; Hexaly (formerly LocalSolver)	LP/MILP/conic; heuristic global search	newer competitors

Licences are usually sold per machine, per user or per cloud hour. All vendors offer free academic licences, and some publish size-limited editions on PyPI: pip install gurobipy gives a licence limited to 2,000 variables and 2,000 constraints, and pip install cplex docplex a Community Edition limited to 1,000 of each. That is enough for learning and for the examples in this chapter.

Day to day, the parameters below matter most. Their names differ, but their meaning is shared:

Purpose	Gurobi (gurobipy)	CPLEX (docplex)
Stop after a time limit	Params.TimeLimit	parameters.timelimit
Relative optimality gap	Params.MIPGap (default $10^{-4}$)	parameters.mip.tolerances.mipgap (default $10^{-4}$)
Threads	Params.Threads	parameters.threads
Emphasis: feasibility or bound	Params.MIPFocus	parameters.emphasis.mip
Random seed	Params.Seed	parameters.randomseed
LP algorithm	Params.Method	parameters.lpmethod
Infeasibility diagnosis	computeIIS(), feasRelaxS()	conflict refiner, feasopt
Automatic tuning	model.tune(), grbtune	tuning tool

Example: capacitated facility location in Gurobi and CPLEX. The instance generator below is used throughout the rest of this chapter. It creates random single-source capacitated facility location problems (Section 4.3): open sites with fixed costs, and assign each customer to exactly one open site without exceeding capacities.

# cfl.py
import numpy as np

def make_cfl(n_fac, n_cust, seed=0):
    """Random capacitated facility location instance (Cornuejols-style)."""
    rng = np.random.default_rng(seed)
    fac, cust = rng.random((n_fac, 2)), rng.random((n_cust, 2))           # points in the unit square
    d = rng.integers(5, 36, n_cust)                                          # customer demands
    cap = rng.integers(10, 161, n_fac).astype(float)
    cap *= 2.0 * d.sum() / cap.sum()                                         # total capacity = 2 x demand
    cap = np.round(cap)
    fixed = np.round(rng.uniform(0, 90, n_fac) + rng.uniform(100, 110, n_fac) * np.sqrt(cap))
    c = np.round(10 * d[None, :] * np.linalg.norm(fac[:, None] - cust[None], axis=2))   # cost of serving j from i
    return fixed, cap, d, c

The Gurobi model:

import gurobipy as gp
from gurobipy import GRB
from cfl import make_cfl

fixed, cap, d, c = make_cfl(20, 40, seed=1)          # 20 sites, 40 customers: 820 binaries
I, J = range(len(fixed)), range(len(d))

m = gp.Model("cfl")
y = m.addVars(I, vtype=GRB.BINARY, name="open")
x = m.addVars(I, J, vtype=GRB.BINARY, name="assign")  # single sourcing: each customer uses one site
m.setObjective(gp.quicksum(fixed[i] * y[i] for i in I)
               + gp.quicksum(c[i, j] * x[i, j] for i in I for j in J), GRB.MINIMIZE)
m.addConstrs((x.sum("*", j) == 1 for j in J), name="serve")
m.addConstrs((gp.quicksum(d[j] * x[i, j] for j in J) <= cap[i] * y[i] for i in I), name="cap")
m.addConstrs((x[i, j] <= y[i] for i in I for j in J), name="link")   # redundant but tightens the LP

m.Params.OutputFlag = 0       # silence the log (Section 13.2 reads it)
m.Params.TimeLimit = 60       # seconds
m.Params.MIPGap = 1e-4        # stop when (incumbent - bound) / |incumbent| <= 0.01% (the default)
m.Params.Threads = 4
m.optimize()

print(f"status={m.Status} (2 = OPTIMAL)  obj={m.ObjVal:.0f}  bound={m.ObjBound:.1f}  "
      f"gap={m.MIPGap:.2e}  nodes={m.NodeCount:.0f}  time={m.Runtime:.2f}s")
print("open sites:", [i for i in I if y[i].X > 0.5])
# status=2 (2 = OPTIMAL)  obj=10040  bound=10040.0  gap=0.00e+00  nodes=1  time=0.34s
# open sites: [3, 4, 8, 12, 14, 15, 16]

The same model in CPLEX through docplex:

from docplex.mp.model import Model
from cfl import make_cfl

fixed, cap, d, c = make_cfl(20, 40, seed=1)
I, J = range(len(fixed)), range(len(d))

m = Model("cfl")
y = m.binary_var_dict(I, name="open")
x = m.binary_var_matrix(I, J, name="assign")
m.minimize(m.sum(fixed[i] * y[i] for i in I) + m.sum(c[i, j] * x[i, j] for i in I for j in J))
m.add_constraints(m.sum(x[i, j] for i in I) == 1 for j in J)
m.add_constraints(m.sum(d[j] * x[i, j] for j in J) <= cap[i] * y[i] for i in I)
m.add_constraints(x[i, j] <= y[i] for i in I for j in J)

m.parameters.timelimit = 60
m.parameters.mip.tolerances.mipgap = 1e-4
m.parameters.threads = 4
sol = m.solve()

det = m.solve_details
print(f"status={det.status}  obj={sol.objective_value:.0f}  bound={det.best_bound:.1f}  "
      f"gap={det.mip_relative_gap:.2e}  nodes={det.nb_nodes_processed}  time={det.time:.2f}s")
print("open sites:", [i for i in I if y[i].solution_value > 0.5])
# status=integer optimal solution  obj=10040  bound=10040.0  gap=0.00e+00  nodes=0  time=0.51s
# open sites: [3, 4, 8, 12, 14, 15, 16]

Both solvers prove the same optimum, 10040 with seven open sites, in about half a second and without branching. Presolve, cuts and heuristics close the gap at the root node; Gurobi counts the root as one node and CPLEX as zero. The two programs differ only in API names, which is typical: moving between commercial solvers is easy, and a modeling layer such as Pyomo or PuLP makes it a one-line change. The interesting differences appear on harder instances, in the next topic.

Key takeaways

• Commercial MIP solvers combine presolve, LP algorithms, cutting planes, heuristics and parallel branch-and-bound. Their engineering makes them the fastest option on hard MIPs.
• Gurobi and CPLEX share the same concepts and default tolerances; only API and parameter names differ.
• Free academic licences and size-limited pip editions make them accessible for learning. Production use needs a paid licence.
• The parameters that matter most are the time limit, the optimality gap, threads and the random seed.

Practice

1. (Beginner) With MIPGap = 1e-4 and an incumbent of 10040, how far above the true optimum can the reported solution be in the worst case? Why does a relative gap make more sense than an absolute one for most models?
2. (Intermediate) Solve make_cfl(30, 50, seed=3) with Gurobi using Threads = 1, 2 and 4, each with seeds 0 to 4. Report the median time per thread count and the parallel speedup. Why does doubling the threads rarely halve the time of a MIP?
3. (Interview) "A Gurobi licence costs tens of thousands per year. Our nightly planning model solves in 40 minutes with HiGHS. Should we buy it?" Hint: weigh the time window, model growth, re-optimization during the day, the value of a smaller gap in money, engineering time spent tuning open-source solvers, and the option to test on an evaluation licence.

Open-Source Solvers (HiGHS, CBC, SCIP)

Why it matters. Open-source solvers are free to deploy on any number of machines and containers, auditable, embeddable in other products, and free of licence servers. That matters for startups, for research that others must reproduce, for teaching, and for services that run thousands of small solves. Their quality has improved sharply in recent years. HiGHS now powers SciPy's linprog and milp, and is the usual open-source choice in JuMP, Pyomo and PuLP.

Theory. The main open-source options:

Solver	Developer and licence	Problem classes	Strengths	Weaknesses
HiGHS	University of Edinburgh, MIT	LP (dual and primal simplex, interior point, PDLP), MILP, convex QP	fast LP, clean API, inside SciPy, PuLP, Pyomo, JuMP	MIP behind the commercial leaders
CBC	COIN-OR, EPL 2.0	MILP	mature, embedded in many tools (PuLP, OpenSolver)	slow on harder MIPs, little active development
SCIP	Zuse Institute Berlin, Apache 2.0	MILP, nonconvex MINLP (global), constraint integer programs	research framework for branch-cut-and-price, plugins	slower LP engine
CP-SAT (OR-Tools)	Google, Apache 2.0	integer models, scheduling, assignment	very strong on combinatorial and scheduling problems	integer data only, no continuous LP
GLPK	GNU, GPL	LP, MILP	simple, educational	slow
Ipopt, Bonmin, Couenne	COIN-OR, EPL	NLP (local), convex MINLP, nonconvex MINLP	standard for nonlinear models	not for large MIPs

Benchmarking solvers requires care, because MIP run times are highly variable:

• Use many instances. The standard test set is MIPLIB 2017, and Hans Mittelmann's public benchmarks (Arizona State University) compare solvers on it. Several commercial vendors have withdrawn from those public comparisons, so they cover mainly open-source and newer solvers.
• Fix the hardware, the threads and a time limit, and count unsolved instances at the time limit.
• Summarize with the shifted geometric mean $\big(\prod_{k=1}^{n} (t_k + s)\big)^{1/n} - s$, usually with $s = 10$ seconds. It is not dominated by a single slow instance, and the shift stops very fast instances from distorting the ratio.
• Check the answers: all solvers must agree on the optimal objective of solved instances.

Example: swapping solvers behind one model. PuLP builds the facility location model once and hands it to four solvers. All runs use one thread and a 300-second limit, on a laptop with an Intel Core i7-1255U:

import time
import pulp
from cfl import make_cfl

def build(fixed, cap, d, c):
    """Single-source CFL written once in PuLP; any solver can be plugged in afterwards."""
    I, J = range(len(fixed)), range(len(d))
    m = pulp.LpProblem("cfl", pulp.LpMinimize)
    y = m.add_variable_dict("y_%s", I, cat="Binary")
    x = m.add_variable_dict("x_%s_%s", [(i, j) for i in I for j in J], cat="Binary")
    m += pulp.lpSum(fixed[i] * y[i] for i in I) + pulp.lpSum(c[i, j] * x[i, j] for i in I for j in J)
    for j in J:
        m += pulp.lpSum(x[i, j] for i in I) == 1
    for i in I:
        m += pulp.lpSum(d[j] * x[i, j] for j in J) <= cap[i] * y[i]
        for j in J:
            m += x[i, j] <= y[i]
    return m

solvers = {"CBC":    lambda: pulp.COIN_CMD(msg=False, timeLimit=300, threads=1),
           "HiGHS":  lambda: pulp.HiGHS(msg=False, timeLimit=300, threads=1),
           "SCIP":   lambda: pulp.SCIP_PY(msg=False, timeLimit=300, threads=1),
           "Gurobi": lambda: pulp.GUROBI(msg=False, timeLimit=300, threads=1)}   # reference

print("instance      solver   status     objective  time (s)")
for size, seed in [((20, 40), 1), ((25, 60), 2), ((30, 50), 3)]:
    data = make_cfl(*size, seed)
    for name, make in solvers.items():
        m = build(*data)
        t = time.perf_counter()
        st = m.solve(make())
        print(f"{size[0]}x{size[1]} #{seed}   {name:7s}  {st.status_str:9s}  {pulp.value(m.objective):9.0f}  "
              f"{time.perf_counter() - t:8.2f}")
# instance      solver   status     objective  time (s)
# 20x40 #1   CBC      Optimal        10040     24.05
# 20x40 #1   HiGHS    Optimal        10040      2.45
# 20x40 #1   SCIP     Optimal        10040     11.12
# 20x40 #1   Gurobi   Optimal        10040      0.39
# 25x60 #2   CBC      TimeLimit      12578    300.34
# 25x60 #2   HiGHS    Optimal        12507     51.81
# 25x60 #2   SCIP     Optimal        12507     51.70
# 25x60 #2   Gurobi   Optimal        12507      5.25
# 30x50 #3   CBC      TimeLimit      12697    300.31
# 30x50 #3   HiGHS    Optimal        12595     93.55
# 30x50 #3   SCIP     Optimal        12595    246.39
# 30x50 #3   Gurobi   Optimal        12595     10.71

The solvers that finish agree on the optimal objective of every instance, so their times can be compared. The shifted geometric means ($s = 10$ s, time-outs counted as 300 s) are 4.9 s for Gurobi, 33 s for HiGHS, 59 s for SCIP and at least 139 s for CBC. HiGHS is the fastest open-source solver overall, although SCIP matches it on the 25 x 60 instance; Gurobi is 6 to 10 times faster than HiGHS on each instance. CBC solves only the smallest instance. On the other two it stops at the time limit with plans 0.6% and 0.8% above the optimum, still usable but without a proof. Three instances on one laptop illustrate the method rather than rank the solvers: on other model families, and on pure LPs, the differences can be much smaller or ordered differently. Switching solvers took one line, which is the practical argument for a solver-independent modeling layer.

Key takeaways

• HiGHS is the strongest general open-source LP/MIP solver, SCIP the most flexible (MINLP, branch-cut-and-price), CBC the most widely embedded, and CP-SAT the best for scheduling and combinatorial models.
• A modeling layer makes switching solvers a one-line change, so benchmark on your own instances before choosing.
• Benchmark with several instances and seeds, fixed threads and time limits, shifted geometric means, and a check that the answers agree.
• Open source is often good enough for LPs and moderate MIPs. On hard MIPs, commercial solvers can be an order of magnitude faster.

Practice

1. (Beginner) Why does the benchmark fix threads=1 and a time limit, and why must the objectives of solved instances agree before the times are compared?
2. (Intermediate) Extend the benchmark to five instance sizes with three seeds each, and report the shifted geometric mean time ($s = 10$, time-outs counted as 300 s) and the number of solved instances per solver. Does the ranking change compared with the three instances above?
3. (Interview) "A vendor's slide says their solver is 5 times faster than open-source alternatives. What questions do you ask before believing it?" Hint: which instances and who chose them, threads and hardware, time limits and how time-outs were counted, mean or geometric mean, solution quality at the end, and whether the comparison matches your own models.

Spreadsheet Optimization

Why it matters. Spreadsheets are the most widely used modeling tool in business. Every copy of Excel includes Solver, managers can see and change every number, and many OR projects begin, and some end, as a spreadsheet. Knowing how to build a clean spreadsheet model, read its reports and recognize its limits makes this tool an asset, not a source of hidden errors.

Theory. A spreadsheet model maps directly onto the parts of an optimization model:

Model element	Spreadsheet element
Parameters	data cells (inputs, often colored)
Decision variables	changing cells (variable cells)
Objective	objective cell containing a formula
Constraint left-hand sides	formula cells, usually SUMPRODUCT of data and changing cells
Constraints	relations between cells entered in the Solver dialog

Excel Solver offers three engines:

Engine	For	Guarantees
Simplex LP	linear models, with integer constraints int, bin, dif	global optimum (within tolerance), Sensitivity Report for LPs
GRG Nonlinear	smooth nonlinear models	local optimum only; the Multistart option helps
Evolutionary	non-smooth models (IF, VLOOKUP, ABS)	none; a genetic-algorithm-style search (Section 11.2)

Good practice keeps the model correct and auditable:

• Separate data, decisions, calculations and constraints into visible blocks, and color the changing cells.
• Use SUMPRODUCT with absolute references ($B$3:$C$3) so the linearity of each constraint is obvious and the formulas can be copied down.
• Keep IF, MAX, ROUND and lookups away from the changing cells, or Simplex LP stops with "the linearity conditions required by this LP Solver are not satisfied". Reformulate with binary variables (Section 4.1) instead.
• Tick "Make Unconstrained Variables Non-Negative", and keep the data in sensible units (Section 13.2).

The limits are real. The standard Excel Solver accepts 200 decision variables and 100 constraints, apart from bounds and integrality. The free OpenSolver add-in uses CBC with no size limit, and commercial add-ins (Analytic Solver) go further. LibreOffice Calc and Google Sheets have their own solvers. Beyond a few thousand variables, or when a model must run unattended, it is time to move to a modeling language or a Python library.

Example: the bakery LP in a spreadsheet. The bakery of Chapter 1 makes bread ($x_1$) and cake ($x_2$) batches to maximize $3x_1 + 5x_2$ subject to the bread oven ($x_1 \le 4$), the cake oven ($2x_2 \le 12$) and the mixer ($3x_1 + 2x_2 \le 18$). The sheet after solving:

       A                    B          C          D                               E      F
  1    Bakery plan          Bread      Cake
  2    Profit per batch     3          5
  3    Batches              2          6          <- changing cells B3:C3
  4    Total profit         36                    <- objective cell B4: =SUMPRODUCT(B2:C2, B3:C3)
  5
  6    Resource (hours)     Bread      Cake       Used                            <=     Available
  7    Bread oven           1          0          2    =SUMPRODUCT(B7:C7, $B$3:$C$3)      4
  8    Cake oven            0          2          12   =SUMPRODUCT(B8:C8, $B$3:$C$3)      12
  9    Mixer                3          2          18   =SUMPRODUCT(B9:C9, $B$3:$C$3)      18

  Solver Parameters
  -----------------------------------------------------------------
  Set Objective:          $B$4              To:  (o) Max  ( ) Min
  By Changing Cells:      $B$3:$C$3
  Subject to:             $D$7:$D$9 <= $F$7:$F$9
  [x] Make Unconstrained Variables Non-Negative
  Select a Solving Method: Simplex LP

Solver finds $x = (2, 6)$ with profit 36 and offers an Answer Report (final values, binding constraints, slacks) and a Sensitivity Report:

Variable cell	Final value	Reduced cost	Objective coefficient	Allowable increase	Allowable decrease
$B$3 Batches Bread	2	0	3	4.5	3
$C$3 Batches Cake	6	0	5	1E+30	3

Constraint	Final value	Shadow price	Constraint R.H. side	Allowable increase	Allowable decrease
$D$7 Bread oven Used	2	0	4	1E+30	2
$D$8 Cake oven Used	12	1.5	12	6	6
$D$9 Mixer Used	18	1	18	6	6

1E+30 is Solver's way of writing infinity. Each extra mixer hour is worth 1 (hundred euros), valid for up to 6 extra hours. The plan stays optimal as long as bread's profit stays between 0 and 7.5 (Section 2.5). The same numbers can be checked outside the spreadsheet, which is good practice before a report goes to a manager:

import numpy as np
from scipy.optimize import linprog

# The spreadsheet's numbers: B2:C2 profit, B7:C9 hours per batch, F7:F9 hours available
profit = np.array([3, 5])
hours = np.array([[1, 0], [0, 2], [3, 2]])
avail = np.array([4, 12, 18])
names = ["bread oven", "cake oven", "mixer"]

res = linprog(-profit, A_ub=hours, b_ub=avail, method="highs")   # the "Simplex LP" engine
used = hours @ res.x                                            # column D: =SUMPRODUCT(B7:C7,$B$3:$C$3)
print(f"changing cells B3:C3 = {res.x}, objective B4 = {-res.fun:.0f}")
print("constraint   final value  shadow price  slack")
for n, u, y, a in zip(names, used, -res.ineqlin.marginals, avail):  # sign flip: we minimized -profit
    print(f"{n:11s}  {u:11.0f}  {y:12.2f}  {a - u:5.0f}")

# Check the mixer's shadow price and allowable increase: change F9 and re-solve
for extra in (1, 6, 7):
    r = linprog(-profit, A_ub=hours, b_ub=avail + np.array([0, 0, extra]), method="highs")
    print(f"mixer hours = {18 + extra}: profit = {-r.fun:.0f}  (+{-r.fun - 36:.0f})")
# changing cells B3:C3 = [2. 6.], objective B4 = 36
# constraint   final value  shadow price  slack
# bread oven             2          0.00      2
# cake oven             12          1.50      0
# mixer                 18          1.00      0
# mixer hours = 19: profit = 37  (+1)
# mixer hours = 24: profit = 42  (+6)
# mixer hours = 25: profit = 42  (+6)

The shadow prices match the report, and so does the allowable increase. Each mixer hour adds 1 up to 24 hours (+6), and the seventh extra hour adds nothing, because bread production has then dropped to zero and the basis changes.

Key takeaways

• In a spreadsheet model, changing cells are the variables, a formula cell is the objective, and SUMPRODUCT formulas form the constraints.
• Use Simplex LP for linear models (global optimum and Sensitivity Report), GRG Nonlinear for smooth models (local optimum) and Evolutionary for non-smooth ones (no guarantee).
• Keep non-linear functions away from the decision cells, lay the model out in clear blocks, and verify important results independently.
• Standard Solver is limited to 200 variables. OpenSolver or a real modeling tool takes over beyond that.

Practice

1. (Beginner) Build the bakery sheet, then add a market limit of at most 5 cake batches. Solve it, and explain from the new Sensitivity Report which constraints are now binding.
2. (Intermediate) The mixer can run overtime at a cost of 1.5 per hour for up to 10 hours. A manager raised the mixer limit in F9 to 28 and wrote the profit as =SUMPRODUCT(B2:C2,B3:C3) - IF(D9>18, 1.5*(D9-18), 0). Show that Simplex LP rejects this model, then reformulate it linearly with an overtime changing cell. Use the original shadow price of the mixer to predict whether overtime will be used.
3. (Interview) "Finance built a 5,000-variable budget allocation in Excel and Solver refuses to run it. What are the options?" Hint: OpenSolver or commercial add-ins, decomposing the model, moving to Python with the spreadsheet kept as the data source and report, and what each option means for the people who maintain it.

<a id="132-or-in-practice"></a>

13.2 OR in Practice

A correct model and a good solver are not yet a working application. Solvers compute in floating point and can be misled by badly scaled data, their output must be read correctly, models are often infeasible on the first run, run times grow with real data, and the results must survive contact with the business. This subchapter covers that craft and ends with three case studies. The examples reuse the facility location generator cfl.py of Section 13.1, plus small models built to show one problem at a time.

Numerical Stability and Scaling

Why it matters. Solvers work in double-precision floating point, about 16 significant digits, and accept a solution when its errors are below fixed tolerances. A model whose coefficients span many orders of magnitude, or which contains a huge "big-M", can make a solver declare a wrong answer optimal, report a feasible model as infeasible, or slow down dramatically. These failures are silent: the log says "optimal" and nothing looks wrong until the plan is executed.

Theory. Three facts explain most numerical trouble.

1. Rounding is everywhere. In floating point, 0.1 + 0.2 is 0.30000000000000004, and every simplex pivot adds a small error. When solving with a basis matrix $B$, the relative error grows roughly like $\kappa(B) \cdot \varepsilon$, where $\kappa$ is the condition number and $\varepsilon \approx 2.2 \times 10^{-16}$. With $\kappa = 10^{10}$, only about 6 digits can be trusted, which is exactly the scale of the tolerances below.
2. Tolerances define "feasible" and "integer". They are mostly absolute, so their meaning depends on the units of the model:

Tolerance	Meaning	Gurobi	CPLEX	HiGHS
Primal feasibility	allowed violation of a constraint	$10^{-6}$	$10^{-6}$	$10^{-7}$
Optimality (dual feasibility)	allowed wrong-signed reduced cost	$10^{-6}$	$10^{-6}$	$10^{-7}$
Integrality	distance from an integer that still counts as integer	$10^{-5}$	$10^{-5}$	$10^{-6}$
Relative MIP gap	when branch-and-bound may stop	$10^{-4}$	$10^{-4}$	$10^{-4}$

3. Big-M constraints multiply tolerances. In $x \le M y$ with binary $y$, the solver may accept $y = 10^{-6}$ as "zero" or as "integer", which allows $x$ up to $M \cdot 10^{-6}$. With $M = 10^9$ that is 1,000 units of flow through a facility that is supposedly closed, a trickle flow (Section 4.1).

Practical rules:

• Choose units that keep coefficients roughly between $10^{-3}$ and $10^{6}$: thousands of euros instead of euros, tonnes instead of grams. The solver log prints the ranges ("Matrix range", "Objective range", "RHS range"); read them first.
• Make every big-M as small as the data allow: the largest flow that can really occur, not "a large number". Better still, use indicator constraints or variable bounds.
• Remove noise coefficients such as $10^{-12}$ that come from data errors, and write exact fractions instead of rounded ones such as 0.333.
• Solvers scale automatically, multiplying rows and columns by powers of 2, but scaling cannot fix a model whose coefficients genuinely span 15 orders of magnitude.
• When in doubt, use the solver's numerical settings (Gurobi's NumericFocus, CPLEX's emphasis.numerical), check the condition number (Gurobi's KappaExact attribute), and verify the solution independently after rounding the integer variables.

Example: a trickle flow. Three plants can supply a customer who needs $D = 50$ units. Opening a plant costs 1000, 800 or 1200, and shipping costs 5, 7 or 4 per unit. The cheapest real plan opens plant 2 for a cost of $800 + 7 \cdot 50 = 1150$. The model is solved twice, with a tight $M = D$ and with a "safe" $M = 10^9$:

import numpy as np
from scipy.optimize import milp, LinearConstraint, Bounds

# Fixed-charge supply: 3 plants, one customer needing D = 50 units.
# Variables [x1, x2, x3, y1, y2, y3]: x = units shipped, y = 1 if the plant is opened.
fixed = np.array([1000.0, 800.0, 1200.0])
unit = np.array([5.0, 7.0, 4.0])
D = 50

def solve(M):
    A = np.zeros((4, 6))
    A[0, :3] = 1                                          # x1 + x2 + x3 >= D
    for i in range(3):
        A[1 + i, i], A[1 + i, 3 + i] = 1, -M              # x_i - M y_i <= 0  (big-M link)
    nz = np.abs(A[A != 0])
    print(f"M = {M:.0e}: matrix range [{nz.min():.0e}, {nz.max():.0e}], ratio {nz.max() / nz.min():.0e}")
    res = milp(np.concatenate([unit, fixed]),
               constraints=LinearConstraint(A, [D, -np.inf, -np.inf, -np.inf], [np.inf, 0, 0, 0]),
               integrality=[0, 0, 0, 1, 1, 1], bounds=Bounds(0, [np.inf] * 3 + [1] * 3))
    x, y = res.x[:3], res.x[3:]
    print(f"   solver: '{res.message}'  cost = {res.fun:.2f}")
    print(f"   x = {np.round(x, 6)}  y = {y}")
    y_round = np.round(y)                                 # what the business will actually do
    honest = unit @ x + fixed @ y_round
    viol = np.maximum(x - M * y_round, 0).max()           # re-check the links with integer y
    print(f"   after rounding y: cost = {honest:.2f}, max link violation = {viol:.1f}")

for M in (50, 1e9):                                       # tight M (= D) versus "safe" huge M
    solve(M)
# M = 5e+01: matrix range [1e+00, 5e+01], ratio 5e+01
#    solver: 'Optimization terminated successfully. (HiGHS Status 7: Optimal)'  cost = 1150.00
#    x = [ 0. 50.  0.]  y = [-0.  1. -0.]
#    after rounding y: cost = 1150.00, max link violation = 0.0
# M = 1e+09: matrix range [1e+00, 1e+09], ratio 1e+09
#    solver: 'Optimization terminated successfully. (HiGHS Status 7: Optimal)'  cost = 200.00
#    x = [ 0.  0. 50.]  y = [-0.e+00 -0.e+00  5.e-08]
#    after rounding y: cost = 200.00, max link violation = 50.0

With $M = 50$ the answer is correct. With $M = 10^9$, HiGHS reports "Optimal" with a cost of 200: it ships all 50 units from plant 3 while "opening" it at $y_3 = 5 \times 10^{-8}$, well inside the integrality tolerance of $10^{-6}$, and pays only $1200 \cdot 5 \times 10^{-8}$ of its fixed cost. Gurobi returns the same wrong answer when its presolve is switched off. With presolve on, it happens to repair this tiny model, but that cannot be relied on in a large one. Nothing in the status warns of the problem; only the coefficient range of $10^9$ and the independent check after rounding reveal it. The fix is not a tighter tolerance but a better model: $M = D$, the largest flow that can actually occur.

Key takeaways

• Solvers accept solutions within absolute tolerances (about $10^{-6}$), so the meaning of "feasible" depends on the units of the model.
• Keep coefficients roughly within $[10^{-3}, 10^{6}]$ by choosing units, and read the coefficient ranges in the log.
• Huge big-M values combine with the integrality tolerance into trickle flows that are reported as optimal. Derive $M$ from the data, or use indicator constraints.
• Always verify important solutions independently, with integer variables rounded.

Practice

1. (Beginner) A model has costs from 0.01 to 5,000,000 euros, capacities in grams up to $2 \times 10^9$, and demands in kilograms. Propose units that bring all coefficients into a safe range, and say which constraints change.
2. (Intermediate) Run the example for $M = 10^5, 10^6, 10^7, 10^8$. At which value does HiGHS first return the wrong plan? Predict the threshold from the integrality tolerance $10^{-6}$ and $D = 50$, and explain why.
3. (Interview) "The solver says optimal, but when operations executes the plan, a warehouse ships 0.4% more than its capacity. Is the solver buggy?" Hint: think about absolute tolerances, the units of the capacity row, scaling, and what the solution checker should do before a plan is published.

Interpreting Solver Output

Why it matters. A solver returns much more than $x^*$: a status, an objective bound, a gap, duals, reduced costs and statistics. Misreading them leads to wrong decisions. Typical mistakes are to report a time-limited solution as optimal, to read "optimal" in a MIP as "exactly optimal", to quote MIP duals as shadow prices, and to treat one of many optimal plans as the answer.

Theory.

1. Status. Always check the status before reading values:

Status	Meaning	What to do
Optimal	proven optimal within the tolerances and the MIP gap	read the solution and the bound
Infeasible	no solution satisfies all constraints	diagnose with an IIS (next topic)
Unbounded	the objective can improve without limit	look for a missing constraint or a wrong sign
Infeasible or unbounded	presolve proved one of the two	re-solve without presolve or dual reductions to find out which
Time, node or iteration limit	stopped early	check whether a solution exists (SolCount) and report the gap
Numerical trouble	the algorithm lost precision	see the previous topic

2. Bounds and gap in a MIP. For a minimization, the incumbent $z_{\text{inc}}$ is the best solution found (an upper bound on the optimum) and the best bound $z_{\text{bd}}$ is the lowest LP bound among open nodes (a lower bound):

$$ z_{\text{bd}} \le z^* \le z_{\text{inc}}, \qquad \text{gap} = \frac{|z_{\text{inc}} - z_{\text{bd}}|}{|z_{\text{inc}}|}. $$

"Optimal" means gap $\le$ MIPGap, $10^{-4}$ by default. A reported optimum of 1,000,000 may therefore be up to 100 above the true optimum. A time-limited run with gap 2% is a feasible plan that is at most 2% worse than optimal, a perfectly usable statement for managers.

3. LP sensitivity information. For LPs, the duals (shadow prices), reduced costs and slacks of Section 2.5 come with the solution. For MIPs they are not defined. Some solvers report the duals of the LP with all integer variables fixed, which describes only small changes that leave the integer decisions unchanged.

4. Alternative optima. In Section 13.1, CBC and HiGHS returned different shipping plans with the same cost 153.675. Solver versions, parameters, thread counts and even the order of constraints can change which optimal vertex is returned. Do not build tests or reports on specific variable values. Test the objective, or add a secondary objective (lexicographic optimization, Section 7.1) when the choice among optima matters.

5. Anatomy of a log. A MIP log follows the solver's pipeline: header (version, machine, threads), problem statistics and coefficient ranges, presolve, root relaxation, node log, cutting-plane summary and final result. The node log columns in Gurobi are:

Column	Meaning
Expl, Unexpl	nodes explored and nodes still open
Obj, Depth, IntInf	LP value of the current node, its depth, number of fractional integer variables
Incumbent, BestBd, Gap	best solution, best bound and their relative gap
It/Node, Time	simplex iterations per node, elapsed time
H or * in the first column	a new incumbent found by a heuristic (H) or at a node of the tree (*)

Example: reading a log. The facility location instance make_cfl(30, 50, seed=3) (1,530 binaries) was solved with Gurobi on one thread. The log is shortened; [...] marks omitted lines and <-- marks added comments:

Gurobi Optimizer version 13.0.3 build v13.0.3rc0 (linux64 - "Ubuntu 24.04.4 LTS")
CPU model: 12th Gen Intel(R) Core(TM) i7-1255U, instruction set [SSE2|AVX|AVX2]
Thread count: 6 physical cores, 12 logical processors, using up to 1 threads

Optimize a model with 1580 rows, 1530 columns and 6030 nonzeros (Min)
Variable types: 0 continuous, 1530 integer (1530 binary)
Coefficient statistics:                                  <-- first numerical check: all ranges are narrow
  Matrix range     [1e+00, 1e+02]
  Objective range  [2e+00, 1e+03]
  Bounds range     [1e+00, 1e+00]
  RHS range        [1e+00, 1e+00]

Found heuristic solution: objective 26374.000000         <-- a first plan before any LP is solved
Presolve removed 322 rows and 101 columns                <-- presolve shrinks the model
Presolved: 1258 rows, 1429 columns, 7379 nonzeros

Root relaxation: objective 1.236827e+04, 455 iterations, 0.01 seconds (0.01 work units)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 12368.2714    0   41 26374.0000 12368.2714  53.1%     -    0s
H    0     0                    15465.000000 12368.2714  20.0%     -    0s   <-- heuristics improve the incumbent
[... 13 more heuristic improvements ...]
H    0     0                    12906.000000 12368.2714  4.17%     -    0s
     0     0 12410.3162    0   53 12906.0000 12410.3162  3.84%     -    0s   <-- cuts raise the bound
[... 70 lines of root cut rounds and heuristics ...]
H    0     0                    12597.000000 12495.1333  0.81%     -    2s
     0     2 12495.1363    0  130 12597.0000 12495.1363  0.81%     -    2s   <-- branching starts
H    2     2                    12595.000000 12495.1363  0.79%  36.0    2s   <-- final incumbent found here
  1087   229 12564.5147   11   62 12595.0000 12552.5088  0.34%  28.6    5s   <-- the rest is proving it

Cutting planes:
  Gomory: 1
  Lift-and-project: 8
  Cover: 88
  MIR: 62
  StrongCG: 36
  GUB cover: 10
  Inf proof: 9
  Zero half: 2
  RLT: 1

Explored 2871 nodes (68090 simplex iterations) in 7.77 seconds (3.76 work units)
Thread count was 1 (of 12 available processors)

Solution count 10: 12595 12597 12601 ... 12728

Optimal solution found (tolerance 1.00e-04)
Best objective 1.259500000000e+04, best bound 1.259400000000e+04, gap 0.0079%

The story is typical of MIP solving. The root LP bound is 12368, and the root cut rounds raise it to 12495. Heuristics drive the incumbent from 26374 to 12597 within 2 seconds, so a plan within 0.8% of optimal exists after a quarter of the run. The final incumbent, 12595, appears at node 2. The remaining 5.8 seconds and 2,869 nodes only prove that nothing better exists. The final gap of 0.0079% (bound 12594) is below the $10^{-4}$ tolerance, so the status is "optimal" even though a solution of cost 12594 has not been formally excluded. The cut summary shows which families helped. Cover and MIR cuts dominate, as expected for knapsack-like capacity rows.

In code, the end state should be turned into an honest statement, whatever happened:

import gurobipy as gp
from gurobipy import GRB
from cfl import make_cfl

def build(fixed, cap, d, c):
    I, J = range(len(fixed)), range(len(d))
    m = gp.Model("cfl")
    m.Params.OutputFlag = 0
    y = m.addVars(I, vtype=GRB.BINARY, name="open")
    x = m.addVars(I, J, vtype=GRB.BINARY, name="assign")
    m.setObjective(y.prod(dict(enumerate(fixed))) + gp.quicksum(c[i, j] * x[i, j] for i in I for j in J))
    m.addConstrs((x.sum("*", j) == 1 for j in J), name="serve")
    m.addConstrs((gp.quicksum(d[j] * x[i, j] for j in J) <= cap[i] * y[i] for i in I), name="cap")
    m.addConstrs((x[i, j] <= y[i] for i in I for j in J), name="link")
    return m

def report(m):
    """Turn a solver's end state into a statement a decision maker can rely on."""
    name = {GRB.OPTIMAL: "OPTIMAL", GRB.TIME_LIMIT: "TIME_LIMIT", GRB.INFEASIBLE: "INFEASIBLE",
            GRB.INF_OR_UNBD: "INF_OR_UNBD", GRB.UNBOUNDED: "UNBOUNDED"}.get(m.Status, str(m.Status))
    if m.SolCount == 0:                                   # no incumbent: there is nothing to read
        return f"{name}: no solution available"
    gap = abs(m.ObjVal - m.ObjBound) / abs(m.ObjVal)
    verdict = "proven optimal" if m.Status == GRB.OPTIMAL else f"feasible, at most {gap:.2%} above optimum"
    return (f"{name}: cost {m.ObjVal:.0f}, bound {m.ObjBound:.1f} -> {verdict}; "
            f"max violation {m.MaxVio:.1e}")

data = make_cfl(30, 50, seed=3)
m = build(*data); m.optimize(); print("default      ", report(m))
m = build(*data); m.Params.TimeLimit = 0.5; m.optimize(); print("0.5 s limit  ", report(m))
fixed, cap, d, c = data
m = build(fixed, cap * 0.45, d, c); m.optimize(); print("capacity -55%", report(m))
# default       OPTIMAL: cost 12595, bound 12595.0 -> proven optimal; max violation 0.0e+00
# 0.5 s limit   TIME_LIMIT: cost 12675, bound 12441.0 -> feasible, at most 1.85% above optimum; max violation 0.0e+00
# capacity -55% INFEASIBLE: no solution available

With half a second instead of the full run, the plan costs 12675, and the bound guarantees it is at most 1.85% above optimal. In fact it is 0.6% above, because the bound is conservative. Accessing m.ObjVal after the infeasible solve would raise an error, which is why report checks SolCount first. The last field, MaxVio, is the largest constraint or integrality violation of the returned solution, a quick numerical sanity check.

Key takeaways

• Check the status and the number of solutions before reading any value.
• In a MIP, the incumbent and the best bound bracket the optimum. "Optimal" means within the gap tolerance, and a time-limited solution comes with a quality guarantee.
• Duals and reduced costs are meaningful for LPs. MIP duals only describe the LP with the integer decisions fixed.
• Optimal solutions are often not unique. Test objectives, not variable values.
• A log shows where the time goes: finding good solutions (heuristics) or proving them (bound and tree).

Practice

1. (Beginner) Using the log, compute the gap of the line H 0 0 12906.000000 12368.2714 and check it against the printed 4.17%. At the end, why is the status "optimal" when the best bound is 12594 and the incumbent 12595? What parameter would force a proof that 12595 is exactly optimal?
2. (Intermediate) For the transportation LP of Section 13.1, print all duals and the reduced costs of the unused routes with HiGHS. Interpret the reduced cost of Seattle to Topeka, and verify it by forcing one case onto that route.
3. (Interview) "Our nightly MIP reported 'optimal' after 2 hours. A colleague re-ran it with another seed and got a different plan with the same cost. Which one is right, and what should we change in the pipeline?" Hint: alternative optima, the gap tolerance, determinism and seeds, and secondary objectives that make the choice explicit.

Model Debugging and Infeasibility Analysis

Why it matters. Real models are often infeasible on their first run because of data errors, a sign mistake, or business rules that genuinely conflict. Others come back unbounded because a constraint is missing, or "optimal" with an absurd plan. A model with 50,000 constraints reports only "infeasible". Finding the cause requires a systematic workflow and tools that isolate the conflicting constraints.

Theory.

A debugging workflow.

1. Validate the data before modeling. Check units, missing values and duplicate keys, and run aggregate checks such as total capacity $\ge$ total demand. Many "infeasible models" are infeasible data.
2. Start small. Build and test on an instance small enough to solve by hand, then scale up.
3. Name everything. Constraint names like capacity[plant_7,2026-10] turn a list of row numbers into a diagnosis.
4. Read the model. Write it in LP format (model.write("m.lp"), prob.writeLP("m.lp")) and read a few constraints. Many bugs are visible at once: a wrong coefficient, a missing term, an unexpected bound.
5. Check statistics. Look at the numbers of variables and constraints, the coefficient ranges and the variable types: a binary should have bounds $[0, 1]$.
6. Verify the solution independently. Plug it into the original constraints with your own code, with integer variables rounded, and recompute the objective.
7. Compare with a simple plan. A manual or heuristic plan gives a sanity check on the order of magnitude of the cost.

Libraries have bugs too. In PuLP 4.0.0, prob.add_variable_dicts("x", (I, J), cat="Binary") with two index sets created integer variables without the $[0, 1]$ bounds, and the facility location model of Section 13.1 then returned negative costs. Steps 4 to 6 catch exactly this kind of error. The flat form add_variable_dict used in Section 13.1 is not affected.

Irreducible infeasible subsystems. An IIS is a subset of the constraints and bounds that is infeasible, while removing any one of its members makes it feasible. It is the smallest self-contained explanation of the infeasibility, often a handful of rows out of thousands. A model can have several IISs: fixing one may reveal another. IISs are found by:

• the deletion filter: remove each constraint in turn. If the rest is still infeasible, leave it out; otherwise, put it back. This takes one feasibility check per constraint.
• the additive and elastic filters, which are faster on large models;
• built-in tools: Gurobi's computeIIS() and CPLEX's conflict refiner. The underlying certificate is a Farkas ray (Section 2.2).

Feasibility relaxation (elastic programming). Add a violation variable $v_k \ge 0$ to each constraint, $a_k^\top x \le b_k + v_k$, and minimize $\sum_k w_k v_k$, or the original objective plus this penalty. The result is the cheapest way to repair the model. If the weights $w_k$ are real business costs (overtime per hour, penalties per late unit), the relaxation is often the plan management actually wants. Gurobi's feasRelax and CPLEX's feasopt automate it.

Unbounded models almost always miss a constraint or have a wrong sign or objective sense. Put a large bound (say $10^6$) on every variable, re-solve, and see which variables hit it.

Example: an impossible production week. A furniture plant makes chairs, tables and desks under machine, labor, wood and storage limits, with minimum quantities promised in contracts. The model is infeasible. A deletion filter written in a few lines finds the IIS, and an elastic relaxation prices the repair:

import numpy as np
from scipy.optimize import linprog

# Weekly plan for chairs, tables, desks. Every row is written as  a @ x <= b  and has a name.
rows = {                       #  chairs tables desks    rhs
    "machine_hours":         ([   2,     4,     3],    400),
    "labor_hours":           ([   3,     2,     4],    450),
    "wood_m3":               ([   1,     3,     2],    300),
    "storage":               ([   1,     1,     1],    200),
    "contract_chairs":       ([  -1,     0,     0],    -60),   # chairs >= 60
    "contract_tables":       ([   0,    -1,     0],    -50),   # tables >= 50
    "contract_desks":        ([   0,     0,    -1],    -40),   # desks  >= 40
    "desk_table_mix":        ([   0,   0.5,    -1],      0),   # desks >= 0.5 tables
}
profit = np.array([30, 70, 60])

def feasible(names):
    """Phase-1 check: is the subsystem made of these rows (plus x >= 0) feasible?"""
    if not names:
        return True
    A = np.array([rows[n][0] for n in names]); b = np.array([rows[n][1] for n in names])
    return linprog(np.zeros(3), A_ub=A, b_ub=b, method="highs").status == 0

print("full model feasible?", feasible(list(rows)))

# Deletion filter: drop each row in turn; if the rest is still infeasible, the row is not needed.
iis = list(rows)
for n in list(rows):
    trial = [r for r in iis if r != n]
    if not feasible(trial):
        iis = trial                                        # still infeasible without n: discard n
print("IIS:", iis)

# Elastic relaxation: allow each row to be violated by v_k >= 0 at a business-meaningful price.
penalty = {"machine_hours": 25,                            # overtime: 25 EUR per extra machine hour
           "contract_chairs": 80, "contract_tables": 200, "contract_desks": 150}   # late-delivery fees
names = list(rows)
A = np.array([rows[n][0] for n in names], float); b = np.array([rows[n][1] for n in names], float)
k = len(names)
A_el = np.hstack([A, -np.eye(k)])                          # a @ x - v_k <= b
w = np.array([penalty.get(n, 1e4) for n in names])        # physical limits: (almost) never violate
res = linprog(np.concatenate([-profit, w]), A_ub=A_el, b_ub=b, method="highs")
x, v = res.x[:3], res.x[3:]
print("relaxed plan (chairs, tables, desks):", np.round(x, 1), " profit:", round(profit @ x))
for n, vk in zip(names, v):
    if vk > 1e-6:
        print(f"  violate {n} by {vk:.1f}  (cost {penalty.get(n, 1e4) * vk:.0f} EUR)")
# full model feasible? False
# IIS: ['machine_hours', 'contract_chairs', 'contract_tables', 'contract_desks']
# relaxed plan (chairs, tables, desks): [60. 50. 40.]  profit: 7700
#   violate machine_hours by 40.0  (cost 1000 EUR)

Gurobi finds the same IIS with one call:

import gurobipy as gp
from ex_furniture import rows          # the rows dictionary above, saved as a module

m = gp.Model(); m.Params.OutputFlag = 0
x = m.addVars(3, name=["chairs", "tables", "desks"])
for n, (a, b) in rows.items():
    m.addConstr(gp.quicksum(a[k] * x[k] for k in range(3)) <= b, name=n)
m.optimize()
print("status:", m.Status, "(3 = INFEASIBLE)")
m.computeIIS()
print("Gurobi IIS rows:", [c.ConstrName for c in m.getConstrs() if c.IISConstr])
print("IIS bounds:", [v.VarName for v in m.getVars() if v.IISLB])
# status: 3 (3 = INFEASIBLE)
# Gurobi IIS rows: ['machine_hours', 'contract_chairs', 'contract_tables', 'contract_desks']
# IIS bounds: []

Four of the eight constraints form the conflict, and they translate into one sentence a plant manager understands: the three contracts need $2 \cdot 60 + 4 \cdot 50 + 3 \cdot 40 = 440$ machine hours, but only 400 are available. Labor, wood, storage and the product mix are irrelevant to the conflict, which saves a long search. The elastic model then chooses the cheapest repair. Forty hours of overtime cost 1,000 euros, while delaying any contract would cost more. It does not add further overtime to produce beyond the contracts, because the best profit per machine hour (desks, $60/3 = 20$ euros) is below the overtime cost of 25.

Key takeaways

• Debug systematically: validate the data, start small, name everything, read the LP file, check statistics, and verify solutions independently.
• An IIS is a minimal set of conflicting constraints. It turns "infeasible" into an explanation that can be acted on.
• A feasibility relaxation with business-meaningful penalties finds the cheapest repair, and is often the plan actually needed.
• Unboundedness usually signals a missing constraint or a wrong sign. Temporary large bounds reveal the culprit.

Practice

1. (Beginner) Explain the IIS of the example to the plant manager in two sentences, and list three different ways to make the plan feasible.
2. (Intermediate) Raise the overtime penalty to 100 euros per hour and re-solve the elastic model. Which constraints are violated now? Explain the choice by comparing, for each product, the late fee plus the lost profit per machine hour saved.
3. (Interview) "Our scheduling model runs at 2 a.m. and sometimes comes back infeasible, and nobody is awake to fix it. How should the system behave?" Hint: an elastic fallback model with penalties, logging the IIS for the morning, alerting, falling back to the last valid plan, and data validation before the solve.

Large-Scale Model Performance Tuning

Why it matters. Models that solve in seconds on test data can take hours on real data, and real deadlines are fixed: the plan must be ready before the trucks leave. Speed comes from several levers with very different effects, and some of them sit outside the solver. Knowing which lever to try first, and how to measure honestly, is the difference between systematic tuning and random parameter changes.

Theory. The time of an optimization run is spent in data preparation, model building, the solver (presolve, root LP, cuts, tree search) and solution extraction. The levers, roughly in order of their usual impact:

Lever	What to do	Notes
Formulation	tighter constraints (Section 4.1), symmetry breaking, fewer big-Ms, good bounds	biggest effect when the solver cannot derive it itself
Accept a gap	stop at 0.5% to 1% instead of 0.01%	most of the run often only proves optimality
Decomposition	Benders, column generation (Section 2.6), rolling horizon, aggregation	for models too large to solve in one piece
Warm starts	a MIP start from a heuristic or yesterday's plan; an LP basis for re-solves	helps when the solver's own heuristics struggle
LP algorithm	dual simplex, barrier (interior point, Section 2.6), concurrent, PDLP	barrier for large sparse LPs, simplex for re-solves
Parameters	MIPFocus, cut and heuristic levels, presolve; automatic tuning tools	defaults are strong; tune on a set of instances, not one
Model building	vectorized APIs, no nested loops over all pairs	in Python, building can take longer than solving
Hardware	more threads, more memory, a faster clock	MIP scales sub-linearly with threads

Measuring requires discipline because of performance variability: changing the random seed, the order of the constraints or the machine changes the path of branch-and-bound, and run times can vary by a factor of 2 or more with no change to the model. An "improvement" observed on one instance with one seed may be noise. Compare configurations on several instances and seeds, summarize with shifted geometric means (Section 13.1), and keep a hold-out set of instances, as in algorithm configuration (Section 11.3).

Example: what actually helps on a facility location model. The first experiment solves make_cfl(30, 50, seed=3) with Gurobi on one thread under different configurations: without the redundant links $x_{ij} \le y_i$, with two parameter changes, with a greedy MIP start, with a 1% gap, and with six different seeds:

import numpy as np
import gurobipy as gp
from gurobipy import GRB
from cfl import make_cfl

fixed, cap, d, c = make_cfl(30, 50, seed=3)
I, J = range(len(fixed)), range(len(d))

def run(strong=True, start=None, seed=0, **params):
    m = gp.Model(); m.Params.OutputFlag = 0; m.Params.Threads = 1; m.Params.Seed = seed
    m.Params.TimeLimit = 300
    for k, v in params.items():
        m.setParam(k, v)
    y = m.addVars(I, vtype=GRB.BINARY); x = m.addVars(I, J, vtype=GRB.BINARY)
    m.setObjective(gp.quicksum(fixed[i] * y[i] for i in I) + gp.quicksum(c[i, j] * x[i, j] for i in I for j in J))
    m.addConstrs(x.sum("*", j) == 1 for j in J)
    m.addConstrs(gp.quicksum(d[j] * x[i, j] for j in J) <= cap[i] * y[i] for i in I)
    if strong:
        m.addConstrs(x[i, j] <= y[i] for i in I for j in J)      # disaggregated links (Section 4.1)
    if start is not None:
        for i in I:
            y[i].Start = start[i]                                # partial MIP start: only the open sites
    m.optimize()
    return m.Runtime, m.NodeCount, m.ObjVal

def greedy_open():
    """Cheap heuristic: open sites by fixed cost per unit of capacity until 1.3 x demand is covered."""
    order = np.argsort(fixed / cap); y, tot = np.zeros(len(I)), 0
    for i in order:
        if tot >= 1.3 * d.sum(): break
        y[i], tot = 1, tot + cap[i]
    return y

print("configuration                 time (s)   nodes   objective")
for label, kw in [("strong formulation", dict()),
                  ("weak formulation", dict(strong=False)),
                  ("strong + MIPFocus=2", dict(MIPFocus=2)),
                  ("strong + Cuts=2", dict(Cuts=2)),
                  ("strong + greedy start", dict(start=greedy_open())),
                  ("strong + MIPGap=1%", dict(MIPGap=0.01))]:
    t, n, z = run(**kw)
    print(f"{label:28s}  {t:7.2f}  {n:7.0f}   {z:9.0f}")
times = [run(seed=s)[0] for s in range(6)]
print("strong, seeds 0-5:", " ".join(f"{t:.2f}" for t in times), f"-> max/min = {max(times) / min(times):.1f}")
# configuration                 time (s)   nodes   objective
# strong formulation               7.66     2871       12595
# weak formulation                 5.10     2329       12595
# strong + MIPFocus=2             13.43     3573       12595
# strong + Cuts=2                 14.17     3176       12595
# strong + greedy start            7.73     2871       12595
# strong + MIPGap=1%               1.67        1       12611
# strong, seeds 0-5: 7.68 6.80 9.23 9.99 10.63 7.17 -> max/min = 1.6

The table holds five lessons:

• Formulation advice must be tested. Without the 1,500 link rows the model solved faster (5.1 s against 7.7 s). Gurobi's presolve and implied-bound cuts generate the links it actually needs, and the full set only adds LP work. With HiGHS the two versions take about the same time. The strong formulation is still the safer default, because a solver that cannot derive the links pays dearly without them (Section 4.1). The lesson is to measure rather than assume.
• Default parameters are hard to beat. MIPFocus=2 (focus on the bound) and aggressive cuts almost doubled the run time. Vendors tune their defaults on thousands of models.
• A MIP start helps only when the solver struggles to find solutions. Here its heuristics found better plans than the greedy start in the first second, so nothing changed.
• Accepting a 1% gap gives the largest gain: 1.67 s instead of 7.66 s, solved at the root node. The plan is only 0.13% more expensive than the optimum, because the bound, not the plan, was the slow part.
• Seeds alone change the time by a factor of 1.6 (6.8 to 10.6 s) on an identical model. Any single-run comparison within that range is noise. Even the weak formulation's advantage should be confirmed with several seeds and instances before it is adopted.

The second experiment builds and solves a transportation LP with 400 warehouses and 400 stores (160,000 variables) in three ways:

import time
import numpy as np
import scipy.sparse as sp
import pulp
import pyomo.environ as pyo
from scipy.optimize import linprog

# A large transportation LP: 400 warehouses x 400 stores = 160,000 variables
rng = np.random.default_rng(0)
n = 400
supply = rng.integers(50, 150, n); demand = rng.integers(40, 120, n)
demand = demand * supply.sum() // demand.sum()             # keep total supply >= total demand
cost = rng.uniform(1, 100, (n, n))
N = range(n)

# (1) Term-by-term construction in two modeling layers (Python loops over 160,000 terms)
t = time.perf_counter()
m = pulp.LpProblem("big", pulp.LpMinimize)
x = m.add_variable_dict("x_%s_%s", [(i, j) for i in N for j in N], lowBound=0)
m += pulp.lpSum(cost[i, j] * x[i, j] for i in N for j in N)
for i in N:
    m += pulp.lpSum(x[i, j] for j in N) <= int(supply[i])
for j in N:
    m += pulp.lpSum(x[i, j] for i in N) >= int(demand[j])
print(f"PuLP build (loops):        {time.perf_counter() - t:6.2f} s")

t = time.perf_counter()
p = pyo.ConcreteModel()
p.x = pyo.Var(N, N, domain=pyo.NonNegativeReals)
p.obj = pyo.Objective(expr=sum(cost[i, j] * p.x[i, j] for i in N for j in N))
p.sup = pyo.Constraint(N, rule=lambda p, i: sum(p.x[i, j] for j in N) <= int(supply[i]))
p.dem = pyo.Constraint(N, rule=lambda p, j: sum(p.x[i, j] for i in N) >= int(demand[j]))
print(f"Pyomo build (loops):       {time.perf_counter() - t:6.2f} s")

# (2) Vectorized construction: build the constraint matrix directly with sparse Kronecker products
t = time.perf_counter()
A_sup = sp.kron(sp.eye(n), np.ones((1, n)))               # row i sums x[i, :]
A_dem = sp.kron(np.ones((1, n)), sp.eye(n))               # row j sums x[:, j]
A = sp.vstack([A_sup, -A_dem]).tocsr()
b = np.concatenate([supply, -demand])
print(f"sparse build (vectorized): {time.perf_counter() - t:6.2f} s")

# (3) Same LP, two algorithms
for method in ("highs-ds", "highs-ipm"):                  # dual simplex vs interior point
    t = time.perf_counter()
    res = linprog(cost.ravel(), A_ub=A, b_ub=b, method=method)
    print(f"{method:9s} solve: {time.perf_counter() - t:6.2f} s   cost = {res.fun:,.1f}   iterations = {res.nit}")
# PuLP build (loops):          1.66 s
# Pyomo build (loops):         2.54 s
# sparse build (vectorized):   0.02 s
# highs-ds  solve:   1.79 s   cost = 58,622.2   iterations = 1728
# highs-ipm solve:   2.76 s   cost = 58,622.2   iterations = 20

Building the model term by term takes as long as solving it, and about 100 times longer than vectorized construction. Model generation scales linearly, so at 10 million variables the loops take minutes on every run, and code with hidden quadratic loops (for example sum(x[i, j] for (i, j) in arcs if i == k) inside a loop over k) takes hours. Vectorized interfaces such as the Gurobi matrix API, highspy with NumPy arrays, linopy, or the sparse construction above remove that cost. On this LP the dual simplex beats the interior-point method, which needs only 20 iterations but much more work per iteration. On other large LPs the ranking reverses, which is why commercial solvers offer a concurrent mode that runs both and keeps the first answer.

Key takeaways

• Try the levers in order of impact: formulation, an acceptable gap, decomposition, then warm starts, the LP algorithm and parameters.
• Default parameters are strong. Tune on a representative set of instances and seeds, never on one run.
• Performance variability can reach a factor of 2 from the seed alone, so measure with several seeds and shifted geometric means.
• Formulation advice is a hypothesis to test: modern presolve and cuts recover some tightenings automatically.
• In Python, model building can dominate the run time. Use vectorized APIs for large models.

Practice

1. (Beginner) Classify each intervention by lever: (a) dropping the gap from 0.01% to 0.5%; (b) adding $y_1 \ge y_2 \ge \dots$ for identical trucks; (c) switching from barrier to dual simplex for a daily re-solve; (d) replacing a Python loop over all pairs by a sparse matrix product.
2. (Intermediate) Compare the default configuration with MIPFocus=1 on 5 instances from make_cfl(30, 50, seed) with 3 Gurobi seeds each. Report the shifted geometric mean time ($s = 1$ s) of each configuration. Is the difference larger than the variability between seeds?
3. (Interview) "Our nightly MIP takes 5 hours and must finish in 1. Walk me through your approach." Hint: profile where the time goes (build or solve, finding or proving), read the log, consider the acceptable gap, the formulation and symmetry, warm starts from yesterday's plan, decomposition or a rolling horizon, then tuning on a set of instances with seeds, and finally hardware.

Case Studies in Logistics, Finance and Cloud Computing

Why it matters. OR creates value only when its decisions are deployed. Logistics, finance and cloud computing are three domains where optimization runs every day at large scale. Each one illustrates a different model class, a different tool and different practical difficulties. This topic presents one compact but realistic model per domain, together with the lessons that apply when such a model goes into production.

Theory. Each case follows the same structure: the decision, the model, the tool, and the lessons.

	Logistics	Finance	Cloud computing
Decision	which van visits which customers, in what order	how to split capital among assets	which host runs each virtual machine
Model	capacitated vehicle routing (Section 4.3)	mean-CVaR portfolio LP (Sections 9.4 and 12.2)	multi-dimensional bin packing with side constraints (Section 4.3)
Tool	OR-Tools routing (construction plus local search)	LP solver (HiGHS)	CP-SAT, plus fast heuristics online
Frequency	daily, with re-optimization during the day	monthly or quarterly rebalancing	milliseconds online; periodic consolidation
Main difficulty	data quality and driver acceptance	estimation error in the inputs	scale, latency and side constraints

Case 1: parcel delivery routing. UPS's ORION system, which optimizes the routes of tens of thousands of drivers, won the 2016 INFORMS Edelman Award. UPS reports savings of about 100 million miles and 300 to 400 million dollars per year. Its hardest problems were not algorithmic: maps had to be corrected, business rules collected from experienced drivers, and drivers persuaded to trust routes that looked strange. The model below is the core of such a system on a small scale: 16 customers, 4 vans with a capacity of 40 parcels, total distance minimized. OR-Tools builds a first solution with a cheapest-arc rule and improves it with guided local search, a metaheuristic that penalizes the arcs of the current local optimum so that the search moves elsewhere (Section 11.2).

import numpy as np
from ortools.constraint_solver import pywrapcp, routing_enums_pb2

# Depot (node 0) and 16 customers on a 100 x 100 km grid; vans carry 40 parcels each.
rng = np.random.default_rng(4)
pts = np.vstack([[50, 50], rng.integers(0, 101, (16, 2))])
demand = [0] + [int(q) for q in rng.integers(3, 12, 16)]
dist = np.round(np.linalg.norm(pts[:, None] - pts[None], axis=2)).astype(int)   # integer km
n_vans, capacity = 4, 40

mgr = pywrapcp.RoutingIndexManager(len(pts), n_vans, 0)          # nodes, vehicles, depot
routing = pywrapcp.RoutingModel(mgr)
arc = routing.RegisterTransitCallback(lambda a, b: int(dist[mgr.IndexToNode(a), mgr.IndexToNode(b)]))
routing.SetArcCostEvaluatorOfAllVehicles(arc)                    # objective: total km
load = routing.RegisterUnaryTransitCallback(lambda a: demand[mgr.IndexToNode(a)])
routing.AddDimensionWithVehicleCapacity(load, 0, [capacity] * n_vans, True, "Load")

def solve(local_search):
    p = pywrapcp.DefaultRoutingSearchParameters()
    p.first_solution_strategy = routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC
    if local_search:
        p.local_search_metaheuristic = routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH
        p.time_limit.FromSeconds(5)                              # metaheuristics need a stopping rule
    return routing.SolveWithParameters(p)

print("construction only:    ", solve(False).ObjectiveValue(), "km")
sol = solve(True)
print("+ guided local search:", sol.ObjectiveValue(), "km")
for v in range(n_vans):
    i, route = routing.Start(v), []
    while not routing.IsEnd(i):
        route.append(mgr.IndexToNode(i)); i = sol.Value(routing.NextVar(i))
    print(f"  van {v}: {' -> '.join(map(str, route + [0]))}  ({sum(demand[k] for k in route)} parcels)")
# construction only:     564 km
# + guided local search: 463 km
#   van 0: 0 -> 0  (0 parcels)
#   van 1: 0 -> 3 -> 15 -> 1 -> 9 -> 7 -> 5 -> 0  (40 parcels)
#   van 2: 0 -> 11 -> 14 -> 12 -> 10 -> 6 -> 0  (38 parcels)
#   van 3: 0 -> 8 -> 4 -> 16 -> 13 -> 2 -> 0  (38 parcels)

Five seconds of local search cut the distance by 18%, and the solution needs only three of the four vans, since the 116 parcels fit into $3 \times 40$. In production the same model gains time windows, driver shifts, travel times from a map service that depend on the time of day, and re-optimization when new orders arrive. Each of these is one more "dimension" or constraint in the routing library.

Case 2: portfolio allocation with downside risk. A fund splits its capital among 8 asset classes. It requires an expected return of 7% per year, caps each class at 40%, and minimizes the 95% conditional value-at-risk (CVaR), the average loss in the worst 5% of months. Rockafellar and Uryasev showed in 2000 that minimizing CVaR over scenarios is a linear program:

$$ \min_{w, \alpha, u} \ \alpha + \frac{1}{(1 - \beta) S} \sum_{s=1}^{S} u_s \quad \text{s.t.} \quad u_s \ge -r_s^\top w - \alpha, \ \ u_s \ge 0, \ \ \mu^\top w \ge R, \ \ \textstyle\sum_k w_k = 1, \ \ 0 \le w \le 0.4, $$

where $r_s$ is the vector of asset returns in scenario $s$ and $\alpha$ becomes the value-at-risk at the optimum. The scenarios are drawn from a fat-tailed distribution. To isolate the effect of estimating risk, the expected returns $\mu$ are assumed known, and the portfolio optimized on $S$ scenarios is evaluated on 200,000 fresh ones:

import numpy as np
import scipy.sparse as sp
from scipy.optimize import linprog

# 8 asset classes; monthly returns with fat tails (Student-t, 4 degrees of freedom).
rng = np.random.default_rng(7)
n, beta, target = 8, 0.95, 0.07 / 12                                    # 95% CVaR, 7% per year
mu = np.array([0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.10, 0.12]) / 12
vol = np.array([0.03, 0.04, 0.05, 0.07, 0.08, 0.10, 0.12, 0.16]) / np.sqrt(12)
L = np.linalg.cholesky((0.3 + 0.7 * np.eye(n)) * np.outer(vol, vol))  # correlation 0.3
def scenarios(k):
    return mu + (rng.standard_t(4, (k, n)) / np.sqrt(2)) @ L.T          # t(4) scaled to unit variance

def min_cvar(R, w_max=0.4):
    """Rockafellar-Uryasev LP. Variables [w (n), alpha, u (S)]: CVaR = alpha + sum(u) / ((1 - beta) S)."""
    S = len(R)
    c = np.concatenate([np.zeros(n), [1], np.full(S, 1 / ((1 - beta) * S))])
    A_tail = sp.hstack([sp.csr_matrix(-R), -np.ones((S, 1)), -sp.eye(S)])   # loss_s - alpha - u_s <= 0
    A_ret = np.concatenate([-mu, [0], np.zeros(S)])[None]                  # expected return >= target
    res = linprog(c, A_ub=sp.vstack([A_tail, A_ret]), b_ub=np.concatenate([np.zeros(S), [-target]]),
                  A_eq=np.concatenate([np.ones(n), [0], np.zeros(S)])[None], b_eq=[1],
                  bounds=[(0, w_max)] * n + [(None, None)] + [(0, None)] * S, method="highs")
    return res.x[:n], res.fun

def cvar(r):
    """Empirical CVaR: average loss in the worst (1 - beta) share of scenarios."""
    loss = np.sort(-r)[::-1]
    return loss[: int(np.ceil((1 - beta) * len(r)))].mean()

R_test = scenarios(200_000)                                            # stands in for the future
print("scenarios   weights (%)                    CVaR95 in-sample   out-of-sample")
for S in (250, 1000, 5000):
    R = scenarios(S)
    w, cv = min_cvar(R)
    print(f"{S:7d}     {np.round(100 * w).astype(int)}     {cv:6.2%}             {cvar(R_test @ w):6.2%}")
# scenarios   weights (%)                    CVaR95 in-sample   out-of-sample
#     250     [ 9  0 35 22 23  4  5  3]      1.81%              2.25%
#    1000     [ 1 18 23 27 15  5  9  2]      2.05%              2.20%
#    5000     [10 16 22 13 13 12  9  4]      2.20%              2.15%

With 250 scenarios, roughly 20 years of monthly data, the optimizer reports a monthly CVaR of 1.81%, but the true value of that portfolio is 2.25%: it is 24% riskier than claimed. The LP has exploited the sampling noise and loaded up on assets that happened to behave well in the sample, the same optimizer's curse as in simulation-based optimization (Section 12.2). With more scenarios, the claim becomes honest (2.20% against 2.15%), the true risk falls, and the weights become more diversified. In practice the expected returns are even harder to estimate than the risk. This is why production portfolio models add robust estimators, limits on turnover and transaction costs, cardinality limits (which make the model a MILP), and above all out-of-sample backtests.

Case 3: virtual machine placement. A cloud provider must place 60 virtual machines (12 services with 5 replicas each, in 7 instance shapes) on hosts with 32 vCPU and 128 GB of RAM, using as few hosts as possible. For fault tolerance, an anti-affinity rule forbids two replicas of the same service on one host. Online, placement decisions are made in milliseconds by heuristics. Google's Borg cluster manager, for example, scores candidate machines with a hybrid of best-fit and worst-fit rules. Periodically, an exact model consolidates the fleet. The following compares first-fit decreasing (FFD), the classic heuristic from Chapter 1, with CP-SAT:

import math
import numpy as np
from ortools.sat.python import cp_model

# 60 VMs from 12 services (5 replicas each); hosts have 32 vCPU and 128 GB of RAM.
rng = np.random.default_rng(3)
shapes = [(2, 4), (2, 16), (4, 8), (4, 32), (8, 16), (8, 64), (16, 32)]       # (vCPU, GB) offerings
svc_shape = [shapes[k] for k in rng.integers(0, len(shapes), 12)]
vms = [(s, *svc_shape[s]) for s in range(12) for _ in range(5)]              # (service, cpu, ram)
CPU, RAM, H = 32, 128, 30                                                     # H = hosts available

lb = max(math.ceil(sum(v[1] for v in vms) / CPU), math.ceil(sum(v[2] for v in vms) / RAM))
print(f"VMs: {len(vms)}, lower bound on hosts: {lb}")

def first_fit_decreasing(anti_affinity):
    """Classic heuristic: largest VMs first (by normalized size), each on the first host where it fits."""
    hosts = []                                                               # [cpu used, ram used, services]
    for s, c, r in sorted(vms, key=lambda v: -(v[1] / CPU + v[2] / RAM)):
        for h in hosts:
            if h[0] + c <= CPU and h[1] + r <= RAM and not (anti_affinity and s in h[2]):
                h[0] += c; h[1] += r; h[2].add(s); break
        else:
            hosts.append([c, r, {s}])
    return len(hosts)

def cp_sat(anti_affinity, seconds=20):
    m = cp_model.CpModel()
    x = {(v, h): m.NewBoolVar(f"x{v}_{h}") for v in range(len(vms)) for h in range(H)}
    used = [m.NewBoolVar(f"used{h}") for h in range(H)]
    for v in range(len(vms)):
        m.AddExactlyOne(x[v, h] for h in range(H))                           # every VM placed once
    for h in range(H):
        m.Add(sum(vms[v][1] * x[v, h] for v in range(len(vms))) <= CPU * used[h])
        m.Add(sum(vms[v][2] * x[v, h] for v in range(len(vms))) <= RAM * used[h])
        if anti_affinity:                                                    # replicas on different hosts
            for s in range(12):
                m.AddAtMostOne(x[v, h] for v in range(len(vms)) if vms[v][0] == s)
    for h in range(H - 1):
        m.Add(used[h] >= used[h + 1])                                        # symmetry breaking
    m.Minimize(sum(used))
    solver = cp_model.CpSolver()
    solver.parameters.max_time_in_seconds = seconds
    solver.parameters.num_workers = 8
    st = solver.Solve(m)
    return int(solver.ObjectiveValue()), int(solver.BestObjectiveBound()), solver.StatusName(st)

for aa in (False, True):
    ffd, (best, bound, status) = first_fit_decreasing(aa), cp_sat(aa)
    print(f"anti-affinity={aa!s:5s}  FFD: {ffd} hosts   CP-SAT: {best} hosts (bound {bound}, {status})")
# VMs: 60, lower bound on hosts: 11
# anti-affinity=False  FFD: 13 hosts   CP-SAT: 11 hosts (bound 11, OPTIMAL)
# anti-affinity=True   FFD: 15 hosts   CP-SAT: 11 hosts (bound 11, OPTIMAL)

CP-SAT reaches the volume lower bound of 11 hosts and proves it optimal, with or without anti-affinity. FFD needs 13 hosts, and 15 once the anti-affinity rule blocks its first-fit choices: 18% and 36% more hardware than necessary. The pattern generalizes. Simple heuristics do well on pure bin packing but degrade as side constraints accumulate, while exact models absorb new rules as a few extra lines. The symmetry-breaking constraint used[h] >= used[h + 1] removes the many equivalent numberings of the hosts, a typical formulation lever (see the previous topic). Real systems combine both approaches: a fast heuristic for each arrival, and a periodic optimization that migrates VMs to free entire hosts, with a limit on the number of migrations.

Key takeaways

• Successful applications pair a model with the right tool: a routing library for vehicle routing, an LP for scenario-based risk, CP-SAT for packing with side constraints.
• In logistics, data quality and user acceptance matter as much as the algorithm.
• In finance, estimation error dominates: optimized portfolios look better in-sample than they are, so validate out of sample.
• In cloud computing, fast heuristics make online decisions and exact models consolidate periodically. Side constraints widen the gap between the two.
• Always compare against a lower bound or a simple baseline, so that the value of the optimization can be quantified.

Practice

1. (Beginner) For each case, name the decision variables, the objective, two constraints, and how often the model would run in production.
2. (Intermediate) Add time windows to the routing case: give each customer a random two-hour window within an 8-hour day, use a time dimension with service times of 10 minutes, and convert distances to minutes at 50 km/h. How much longer does the total route become, and how many vans are now needed?
3. (Interview) "Design the placement service for a cloud region: VM requests arrive continuously and must be placed within 50 ms, yet the fleet should stay consolidated." Hint: separate the online heuristic from periodic re-optimization, define the objective (hosts, fragmentation, migrations), consider anti-affinity and failure domains, set limits on migrations per hour, and decide how to measure the gap to a lower bound.