Linear optimization for a box plant: what the solvers are for, and a job shop that proves it
Three kinds of solver, when each one fits, and when none of them is the right tool. Then two real models for a corrugated plant: a product-mix LP that tells you what an hour of the corrugator is worth, and a CP-SAT job shop that schedules eight orders across four machines to proven optimality — in a hundredth of a second.
The algorithms posts so far have been about writing the search: knapsack by hand, rod cutting by hand, a truck-loading heuristic by hand. This one is about when to stop writing it and hand the problem to a solver. Every number below comes from running the code at the bottom.
What "linear" buys you
A linear program is: maximize (or minimize) a linear function of some decision variables, subject to linear inequalities. CLRS chapter 29 opens with the two-variable picture, and it's still the right way to see it: each constraint is a half-plane, the feasible region is the polygon where they all overlap, and the objective is a family of parallel lines sliding across it. The optimum is always at a corner.
That picture is the whole theory. Because the region is convex and the objective is linear, a local optimum is the global optimum, and you can find it by walking corner to corner (simplex) or cutting through the interior (barrier methods). Both are polynomial in practice, and modern solvers handle millions of variables. If your problem is linear and your variables are continuous, it's solved.
The three solvers and what breaks the picture
| Solver | Variables | What it's for | What breaks it |
|---|---|---|---|
| LP (GLOP, HiGHS) | continuous | Product mix, blending, flows, capacity planning. "How much of each?" Shadow prices come free. | Anything that has to be whole: you can't run 2.4 trucks. |
| MILP (SCIP, CBC, HiGHS, Gurobi) | continuous + integer | Fixed costs, yes/no decisions, cutting stock patterns, facility location, lot sizing. The LP with "whole" enforced. | Integrality destroys convexity — the solver branches on every integer, and the tree is where the exponential hides. Tight formulations and good bounds decide whether it's seconds or days. |
| CP-SAT (OR-Tools) | integer / boolean | Scheduling, sequencing, assignment, rostering — anything with "no two jobs on the same machine at once," precedence, or "exactly one of these." | Large continuous quantities. It's a constraint solver with an LP inside, not an LP solver with constraints. |
The rough rule: quantities → LP; quantities with some switches → MILP; sequences and assignments → CP-SAT. And when the model is right but still too slow, that's the signal to reach back for the heuristics from the earlier posts — a first-fit-decreasing pass to hand the solver a warm start, or a decomposition like column generation where the solver only sees a piece at a time.
When a solver is the wrong tool
- The objective isn't known. If nobody can say whether an hour of corrugator time or a day of lateness costs more, the solver will optimize whichever one you typed. Get the weights first.
- The data is worse than the model. A proven-optimal schedule built on run rates that are 20% off is a precise answer to a different question. Fix the inputs before the formulation.
- It's a one-off. Eight orders on four machines, a planner with ten years on the floor does fine with a whiteboard. Solvers earn their keep when the problem recurs every shift, the instance is big enough that intuition fails, or the constraints change faster than habits.
- The exact answer isn't worth the wait. A 2%-worse plan in 5 seconds usually beats the optimal plan in 2 hours. That's what time limits and optimality gaps are for.
Model 1 · Product mix: what is an hour of the corrugator worth?
Three board grades, one week. Each grade has a margin per MSF, uses corrugator hours, converting hours, and tons of linerboard, and has a cap on what the order book will take. Decide how many MSF of each to run.
maximize 38·B + 42·C + 61·BC ($ contribution)
subject to 0.030·B + 0.028·C + 0.055·BC ≤ 120 corrugator hours
0.040·B + 0.040·C + 0.045·BC ≤ 190 converting hours
0.085·B + 0.085·C + 0.130·BC ≤ 400 linerboard tons
B ≤ 1800, C ≤ 2600, BC ≤ 900, all ≥ 0 demand caps
GLOP solves it instantly:
| Grade | Run (MSF) | Reduced cost | Meaning |
|---|---|---|---|
| C | 2,600.0 | +6.53 | At its demand cap. Every extra MSF of C demand is worth $6.53. |
| B | 1,573.3 | 0.00 | Basic — fills what's left of the corrugator. |
| BC | 0.0 | −8.67 | Not worth a slot. BC needs $8.67/MSF more margin before it displaces single wall. |
| Constraint | Used | Shadow price | Meaning |
|---|---|---|---|
| corrugator hours | 120.0 / 120 | $1,266.67 / hr | Binding. One more corrugator hour this week is worth $1,267 of margin. That's the number to put in front of whoever decides on overtime or a speed upgrade. |
| converting hours | 166.9 / 190 | $0 | Slack. Adding converting capacity changes nothing. |
| linerboard tons | 354.7 / 400 | $0 | Slack. |
Objective: $168,987. The surprise for most people is that BC — the highest-margin grade — runs zero. It uses nearly twice the corrugator time per MSF, and the corrugator is the only thing that's scarce, so per scarce hour it's the worst grade on the list. The LP sees that; a margin report doesn't. Duality (CLRS §29.4) is what turns "the constraint is tight" into a dollar figure, and it's the single most useful output a solver gives a plant manager.
Model 2 · Job shop: eight orders, four machines, proven optimal
Now the sequencing problem. Every order flows corrugator → converting → (sometimes) die cutter. There's one corrugator, two printer-slotters that can both do converting, and one die cutter. Each operation has a run time; there's a 10-minute setup between jobs on any machine; each order has a due time. Minimize the makespan, and don't be late.
This is Kleinberg–Tardos §4.2 (scheduling to minimize lateness) stacked three machines deep with a choice of machine in the middle — and that's enough to make it NP-hard. The greedy exchange argument that proves earliest-deadline-first optimal for one machine doesn't survive the second stage. Abdul Bari's job-sequencing lecture (41) is the one-machine greedy version; the branch-and-bound version (69) is what CP-SAT is doing underneath, with far better bounds.
CP-SAT returned OPTIMAL in 0.01 seconds: makespan 455 minutes, total lateness 0. The whiteboard baseline — release orders in the list order, each operation to whichever eligible machine frees up first — gives 485 minutes with 125 minutes of lateness across three orders.
And 455 is provably the floor, which you can check by hand: the corrugator has 340 minutes of work plus seven setups between eight jobs, 410 minutes, and whichever job runs last still needs its converting operation afterwards — the shortest is 45. 410 + 45 = 455. The solver found a sequence that puts the 45-minute converting job last on the corrugator and leaves nothing idle on it. That's not luck; that's the bound and the search meeting.
Look at what it chose. Aurora (due 300, with a long die-cut) goes first. The two double-wall RSCs, which have the longest corrugator runs, go early so their converting can overlap with everything else. The three die-cut orders are spread so the die cutter never has two waiting. No planner would get all of that right on a whiteboard, and this is only eight orders.
The model, in OR-Tools
The whole thing is about forty lines. What makes CP-SAT the right tool here is NoOverlap and optional intervals: you tell it "these intervals can't overlap on this machine" and "exactly one of these alternatives is real," and the disjunctive reasoning — the hard part of scheduling — is handled inside the solver.
from ortools.sat.python import cp_model
m = cp_model.CpModel()
for j, (name, grade, ops, due) in enumerate(ORDERS):
prev_end = None
for k, (stage, dur) in enumerate(ops):
s = m.NewIntVar(0, horizon, f"s{j}_{k}")
e = m.NewIntVar(0, horizon, f"e{j}_{k}")
m.Add(e == s + dur)
if prev_end is not None:
m.Add(s >= prev_end) # precedence within the order
prev_end = e
alts = []
for mach in MACHINES[stage]: # flexible: pick a machine
p = m.NewBoolVar(f"p{j}_{k}_{mach}")
iv = m.NewOptionalIntervalVar(s, dur + SETUP, m.NewIntVar(0, horizon + SETUP, ""), p, "")
intervals[mach].append(iv)
alts.append(p)
m.AddExactlyOne(alts)
job_end.append(prev_end)
L = m.NewIntVar(0, horizon, f"late{j}")
m.Add(L >= prev_end - due) # lateness, floored at 0 by the domain
late.append(L)
for mach, ivs in intervals.items():
m.AddNoOverlap(ivs) # one job at a time per machine
makespan = m.NewIntVar(0, horizon, "makespan")
m.AddMaxEquality(makespan, job_end)
m.Minimize(makespan * 10 + sum(late)) # makespan first, lateness as tiebreak
solver = cp_model.CpSolver()
solver.parameters.max_time_in_seconds = 20
status = solver.Solve(m) # → OPTIMAL, 0.01 s
Things this model leaves out that a plant's version would add, each a few lines: sequence-dependent setups (grade changes on the corrugator cost more than width changes — a AddCircuit or transition-matrix formulation), release times (the paper isn't here until 10:00), machine availability windows (planned maintenance as fixed intervals), and a weighted objective where a late key account costs more than a late one-off. The structure doesn't change. That's the point of modeling it instead of coding it.
The LP, in OR-Tools
from ortools.linear_solver import pywraplp
s = pywraplp.Solver.CreateSolver("GLOP")
x = {k: s.NumVar(0, demand_max[k], k) for k in ["B", "C", "BC"]}
corr = s.Add(sum(corr_hr[k] * x[k] for k in x) <= 120, "corrugator_hours")
conv = s.Add(sum(conv_hr[k] * x[k] for k in x) <= 190, "converting_hours")
lin = s.Add(sum(liner_t[k] * x[k] for k in x) <= 400, "linerboard_tons")
s.Maximize(sum(margin[k] * x[k] for k in x))
s.Solve()
corr.dual_value() # → 1266.67 the shadow price
x["BC"].reduced_cost() # → -8.67 what BC's margin is short by