4. Integer and Combinatorial Optimization
Many real decisions are indivisible: open a warehouse or not, assign a crew to a flight or not, put a VM on a host or not. Integer programming extends LP with variables that must take whole-number values. This one change lets models express logic, fixed costs and combinatorial structure, and it also makes them NP-hard in general. This chapter covers how to model with integer variables, how exact solvers handle them (branch and bound, cutting planes, decomposition), the classical combinatorial problems that appear across OR, and the complexity theory that explains what can and cannot be solved efficiently.
<a id="41-integer-programming-modeling"></a>
4.1 Integer Programming Modeling
This subchapter shows how integer and binary variables turn yes/no decisions, logical rules, fixed costs and disjunctions into linear constraints, and why two formulations of the same problem can behave very differently inside a solver.
Several examples reuse the small integer program from Section 1.3:
$$ \max\ 5x_1 + 8x_2 \quad \text{s.t.} \quad x_1 + x_2 \le 6, \quad 5x_1 + 9x_2 \le 45, \quad x_1, x_2 \ge 0 \text{ integer}. $$
Its LP relaxation has optimum $(2.25, 3.75)$ with value $41.25$. The integer optimum is $(0, 5)$ with value $40$.
Pure, Mixed and Binary Integer Programs
Why it matters. LP assumes divisibility: producing 3.75 cakes or hiring 2.4 nurses is acceptable. That assumption breaks for indivisible items (trucks, people, machines) and especially for decisions (build, assign, select). Knowing which variables really need integrality, and which kind, is the first modeling choice, and it determines how hard the model is to solve.
Theory. An integer program is an LP with additional integrality restrictions. The general mixed integer linear program (MILP) is
$$ \max\ c^\top x + h^\top y \quad \text{s.t.} \quad Ax + Gy \le b, \quad x \in \mathbb{Z}^n_{\ge 0}, \quad y \in \mathbb{R}^p_{\ge 0}. $$
| Class | Variables | Typical use |
|---|---|---|
| Pure integer program (IP / ILP) | all integer | number of trucks, staff per shift, items to produce in batches |
| Mixed integer program (MILP / MIP) | some integer, some continuous | open plants (integer) and ship goods between them (continuous) |
| Binary integer program (BIP / 0-1 IP) | all in ${0, 1}$ | selection, assignment, routing, scheduling decisions |
| Mixed binary program | binaries plus continuous | the most common class in practice: logic on top of an LP |
The feasible region of an LP is a convex polyhedron. The feasible region of an IP is the set of lattice points (integer points) inside that polyhedron, which is a scattered, non-convex set:
x2
5 | * * feasible integer point (25 in total)
| o <- LP optimum (2.25, 3.75), value 41.25
4 | * *
3 | * * * *
2 | * * * * *
1 | * * * * * *
0 | * * * * * * *
+------------------------ x1
0 1 2 3 4 5 6
Integer optimum: (0, 5), value 40 -- far from the LP optimum
Three consequences follow:
- Optimality is no longer local. The integer optimum need not be near the LP optimum, and no simple "move to a better neighbor" rule certifies optimality.
- Bounds instead of certificates. The LP over the same constraints is a relaxation. Its optimum bounds the integer optimum (Section 4.2) but rarely equals it.
- General integers can be written as binaries. If $0 \le x \le U$, then $x = \sum_{k=0}^{K} 2^k z_k$ with $z_k \in {0, 1}$ and $K = \lfloor \log_2 U \rfloor$. So binary programs are as expressive as general integer programs, at the cost of more variables.
When can integrality be dropped? If a variable takes large values (10,000 litres, 3,712 cars), rounding the LP solution changes the objective by a negligible fraction, and an LP is fine. If the constraint matrix is totally unimodular (Section 3.1), the LP returns integer solutions anyway. Integrality matters most for small integers and, above all, for 0/1 decisions.
import numpy as np
from scipy.optimize import milp, LinearConstraint, Bounds
# Running example (Section 1.3): max 5 x1 + 8 x2 s.t. x1 + x2 <= 6, 5 x1 + 9 x2 <= 45
c = np.array([-5, -8]) # milp minimizes, so negate
cons = LinearConstraint([[1, 1], [5, 9]], ub=[6, 45])
for name, integrality in [("LP (x1 real, x2 real)", [0, 0]),
("MILP (x1 real, x2 int) ", [0, 1]),
("MILP (x1 int, x2 real)", [1, 0]),
("ILP (x1 int, x2 int) ", [1, 1])]:
res = milp(c, constraints=cons, integrality=integrality, bounds=Bounds(0, np.inf))
print(f"{name} x = {np.round(res.x, 3) + 0.0} value = {-res.fun:.2f}")
# Binary expansion: x2 in {0..7} becomes z0 + 2 z1 + 4 z2 with z binary
# Variables: [x1, z0, z1, z2]
c_bin = np.array([-5, -8, -16, -32])
cons_bin = LinearConstraint([[1, 1, 2, 4], [5, 9, 18, 36]], ub=[6, 45])
res = milp(c_bin, constraints=cons_bin, integrality=[1, 1, 1, 1],
bounds=Bounds([0, 0, 0, 0], [np.inf, 1, 1, 1])) # z_k in {0, 1}
x1, z = res.x[0], res.x[1:]
print(f"binary expansion: x1 = {abs(x1):.0f}, z = {z.round().astype(int)}, "
f"x2 = {z @ [1, 2, 4]:.0f}, value = {-res.fun:.0f}")
# LP (x1 real, x2 real) x = [2.25 3.75] value = 41.25
# MILP (x1 real, x2 int) x = [1.8 4. ] value = 41.00
# MILP (x1 int, x2 real) x = [2. 3.889] value = 41.11
# ILP (x1 int, x2 int) x = [0. 5.] value = 40.00 <- every restriction lowers the value
# binary expansion: x1 = 0, z = [1 0 1], x2 = 5, value = 40
Each added integrality restriction can only shrink the feasible set, so the optimal value can only get worse: $41.25 \ge 41.11 \ge 40$ and $41.25 \ge 41.00 \ge 40$.
Key takeaways
- ILP: all variables integer. MILP: a mix. BIP: all binary. Mixed binary programs are the most common in practice.
- The integer feasible set is the lattice points inside a polyhedron: non-convex and without local optimality certificates.
- Any bounded integer variable can be rewritten with $\lfloor \log_2 U \rfloor + 1$ binaries.
- Impose integrality only where it matters: small integers and yes/no decisions. Large quantities can usually be rounded.
Practice
- (Beginner) List all integer points of the running example with $x_2 = 3$ and compute their objective values. Which one is best, and how does it compare with $40$?
- (Intermediate) A model has a variable $0 \le x \le 1000$ that must be integer. How many binaries does its binary expansion need? Write the expansion, and explain why solvers usually prefer to keep $x$ as a general integer.
- (Interview) "A colleague solved an LP for a fleet plan and rounded 7.4 trucks down to 7. When is this acceptable and when is it dangerous?" Hint: consider the size of the values relative to the rounding error, and whether rounding can break feasibility (Section 1.3).
Logical Constraints with Binary Variables
Why it matters. Business rules are rarely pure arithmetic: "project B requires project A", "choose at most two suppliers", "open the plant only if both licences are granted". Binary variables turn propositional logic into linear constraints, so a single MILP can combine logic with the quantitative parts of the model.
Theory. Think of each binary variable $x_a \in {0, 1}$ as a light switch for the statement $a$: on means true, off means false. Logical operations become linear constraints on the switches:
| Logical rule | Linear constraint(s) |
|---|---|
| NOT $a$ | use $1 - x_a$ |
| at most one of $a, b, c$ | $x_a + x_b + x_c \le 1$ |
| exactly $k$ of the set $S$ | $\sum_{j \in S} x_j = k$ |
| at least one of $S$ | $\sum_{j \in S} x_j \ge 1$ |
| $a \Rightarrow b$ (if $a$ then $b$) | $x_a \le x_b$ |
| $a \Leftrightarrow b$ | $x_a = x_b$ |
| $a$ and $b$ mutually exclusive | $x_a + x_b \le 1$ |
| $z = a \wedge b$ (AND) | $z \le x_a,\ z \le x_b,\ z \ge x_a + x_b - 1$ |
| $z = a \vee b$ (OR) | $z \ge x_a,\ z \ge x_b,\ z \le x_a + x_b$ |
| $z = a \oplus b$ (XOR) | $z \ge x_a - x_b,\ z \ge x_b - x_a,\ z \le x_a + x_b,\ z \le 2 - x_a - x_b$ |
A general recipe. Write the rule in conjunctive normal form (an AND of OR-clauses). Each clause such as $(a \vee \neg b \vee c)$ becomes $x_a + (1 - x_b) + x_c \ge 1$: at least one literal in the clause must be true.
Products of variables. The AND row also linearizes the product $z = x_a x_b$ of two binaries. For a binary $x$ times a continuous $0 \le y \le U$, the product $z = xy$ is linearized by
$$ z \le Ux, \qquad z \le y, \qquad z \ge y - U(1 - x), \qquad z \ge 0 . $$
When $x = 0$ the first row forces $z = 0$. When $x = 1$ the second and third rows force $z = y$.
import numpy as np
from itertools import product
from scipy.optimize import milp, LinearConstraint, Bounds
# Linearizations of z = x AND y and z = x OR y, as rows a . (x, y, z) <= b
AND = [([-1, 0, 1], 0), ([0, -1, 1], 0), ([1, 1, -1], 1)] # z<=x, z<=y, z>=x+y-1
OR = [([1, 0, -1], 0), ([0, 1, -1], 0), ([-1, -1, 1], 0)] # z>=x, z>=y, z<=x+y
for name, rows in [("AND", AND), ("OR", OR)]:
for x, y in product([0, 1], repeat=2):
# which z in {0,1} satisfy every row? exactly one should survive
ok = [z for z in (0, 1) if all(np.dot(a, (x, y, z)) <= b for a, b in rows)]
print(f"{name}: x={x} y={y} -> feasible z = {ok}")
# AND: x=0 y=0 -> feasible z = [0]
# AND: x=0 y=1 -> feasible z = [0]
# AND: x=1 y=0 -> feasible z = [0]
# AND: x=1 y=1 -> feasible z = [1]
# OR: x=0 y=0 -> feasible z = [0]
# OR: x=0 y=1 -> feasible z = [1]
# OR: x=1 y=0 -> feasible z = [1]
# OR: x=1 y=1 -> feasible z = [1]
# Project selection with logic: values and costs of 5 projects, budget 10
value = np.array([7, 5, 8, 3, 8])
cost = np.array([4, 3, 5, 2, 6])
A = [cost, # budget
[0, 0, 1, 0, 1], # projects 3 and 5 are mutually exclusive: x3 + x5 <= 1
[-1, 0, 0, 1, 0], # project 4 requires project 1: x4 <= x1
[-1, -1, -1, 0, 0]] # at least one of projects 1, 2, 3: x1 + x2 + x3 >= 1
res = milp(-value, constraints=LinearConstraint(A, ub=[10, 1, 0, -1]),
integrality=np.ones(5), bounds=Bounds(0, 1))
print("selected projects:", [int(j) + 1 for j in np.flatnonzero(res.x > 0.5)], "value", -res.fun)
# selected projects: [1, 3] value 15.0
Without the rules, the best portfolio was projects 2, 3 and 4 with value 16 (Section 1.1). The rule "4 requires 1" makes it infeasible, and the optimum drops to 15.
Key takeaways
- A binary variable is a switch for a true/false statement. Logical connectives become linear inequalities on the switches.
- Implication $a \Rightarrow b$ is $x_a \le x_b$. Mutual exclusion is $x_a + x_b \le 1$.
- Any clause in conjunctive normal form becomes "sum of literals $\ge 1$", with $1 - x$ for negated literals.
- Products of binaries, or of a binary and a bounded continuous variable, can be linearized exactly.
Practice
- (Beginner) Write linear constraints for: "if project 2 is chosen, then projects 3 and 4 must both be chosen."
- (Intermediate) Linearize "at least two of the three plants $a, b, c$ must be open, and if plant $c$ is open then plant $a$ is closed." Verify your constraints by enumerating all 8 combinations, as in the code.
- (Interview) "How would you model $z = \max(x_1, x_2, x_3)$ for binary $x$, and how would it change if the $x_i$ were continuous in $[0, 10]$?" Hint: the binary case is an OR. The continuous case needs one extra binary per term to select which $x_i$ attains the maximum.
Big-M Formulations
Why it matters. Many constraints should hold only sometimes: "if the machine is set up, production is at most its capacity", "if job A goes first, B starts after A ends". The big-M technique switches a constraint on or off with a binary variable. It is the most widely used modeling trick in MILP, and also the most common source of slow, numerically fragile models.
Theory. To enforce $a^\top x \le b$ only when $y = 1$, write
$$ a^\top x \le b + M(1 - y). $$
- If $y = 1$, the constraint is $a^\top x \le b$, as intended.
- If $y = 0$, it becomes $a^\top x \le b + M$, which must be redundant. So $M$ must be at least $\max{a^\top x - b}$ over all solutions that are otherwise feasible.
Think of $M$ as the rating of a safety valve. Rated too low (M too small), it cuts off solutions that are actually valid, and the model returns a wrong answer. Rated absurdly high (M huge), it is useless in the LP relaxation. The solver can set $y$ to a tiny fraction and still get almost all of the effect for free.
Choosing M. Derive it from bounds you know: capacities, demands, horizon lengths. If $x \le u$ is implied by the rest of the model, then $x \le u,y$ is the tightest valid big-M linking constraint.
Why huge M hurts.
- Weak relaxation. In the LP, $x \le My$ is satisfied by $y = x/M$. For large $M$, this $y$ is nearly zero, so the fixed cost attached to $y$ almost disappears from the bound.
- Numerical trouble. Solvers accept $y$ as integer if it is within a tolerance such as $10^{-6}$ of an integer. With $M = 10^9$, the value $y = 10^{-7}$ counts as $0$ yet allows $x$ up to $100$. The solver then reports a solution that "uses the machine without setting it up".
Alternatives. Modern solvers support indicator constraints ($y = 1 \Rightarrow a^\top x \le b$) that are handled by branching rather than by an explicit $M$. When your solver supports them, prefer them over a guessed $M$.
import numpy as np
from scipy.optimize import milp, linprog, LinearConstraint, Bounds
# Meet demand 40 by producing in-house (setup 200, unit cost 2) or outsourcing (unit 9).
# Variables [y, x, w]: y = setup (binary), x = produced, w = outsourced.
# min 200 y + 2 x + 9 w s.t. x + w = 40, x - M y <= 0
cost = np.array([200, 2, 9])
for M in [40, 100, 1_000, 1_000_000]:
A = [[0, 1, 1], [-M, 1, 0]]
lb, ub = [40, -np.inf], [40, 0]
lp = linprog(cost, A_ub=[A[1]], b_ub=[0], A_eq=[A[0]], b_eq=[40],
bounds=[(0, 1), (0, None), (0, None)])
ip = milp(cost, constraints=LinearConstraint(A, lb, ub), integrality=[1, 0, 0],
bounds=Bounds([0, 0, 0], [1, np.inf, np.inf]))
print(f"M = {M:>9,}: LP bound = {lp.fun:7.2f} (y = {lp.x[0]:.5f}) IP optimum = {ip.fun:.0f}")
# M = 40: LP bound = 280.00 (y = 1.00000) IP optimum = 280 <- tight M: LP is exact
# M = 100: LP bound = 160.00 (y = 0.40000) IP optimum = 280
# M = 1,000: LP bound = 88.00 (y = 0.04000) IP optimum = 280
# M = 1,000,000: LP bound = 80.01 (y = 0.00004) IP optimum = 280 <- setup cost vanishes
All four models have the same integer optimum. Only the tight one ($M = 40$, the demand) gives an LP relaxation that sees the setup cost. Branch and bound (Section 4.2) relies on these LP bounds, so the loose models need far more work on larger instances.
Key takeaways
- $a^\top x \le b + M(1 - y)$ enforces the constraint when $y = 1$ and relaxes it when $y = 0$.
- $M$ must be large enough to make the relaxed constraint redundant, and no larger.
- Large $M$ values weaken the LP relaxation and cause numerical errors near integrality tolerances.
- Derive $M$ from known bounds, or use the solver's indicator constraints.
Practice
- (Beginner) In the example, what goes wrong if you set $M = 30$? Re-run the code and explain the result.
- (Intermediate) Variables satisfy $0 \le x_1 \le 8$ and $0 \le x_2 \le 5$. Find the smallest valid $M$ for "if $y = 1$ then $2x_1 + 3x_2 \le 12$" written as $2x_1 + 3x_2 \le 12 + M(1 - y)$.
- (Interview) "Your MILP reports an optimal solution where a binary is 0.0000004 and a linked flow is 380. What happened, and how do you fix it?" Hint: compute $M \cdot y$ and compare it with the flow. Think about integrality tolerances, tighter $M$ values and indicator constraints.
Fixed-Charge Problems
Why it matters. Many costs have a fixed part paid only if an activity happens at all: setting up a production line, opening a depot, chartering a truck, renting a server. Ignoring fixed costs, or pretending they are proportional, leads to plans that spread activity over too many options. Fixed-charge models are the core of lot sizing, facility location and network design.
Theory. The cost of an activity at level $x \ge 0$ is
$$ f(x) = \begin{cases} 0 & \text{if } x = 0 \ K + cx & \text{if } x > 0 \end{cases} $$
cost
| /
| / slope c (variable cost)
| /
| /
K +-----------o jump of size K at x > 0
|
|
0 *------------------------- x
0 (the point x = 0 costs nothing)
This function is discontinuous at $0$ and concave on $x \ge 0$. Minimizing a concave cost is outside what LP can represent. Introduce a binary $y$ (1 = the activity happens) and an upper bound $u$ on $x$:
$$ \min\ Ky + cx \quad \text{s.t.} \quad x \le u,y, \quad y \in {0, 1}. $$
If $x > 0$ then $y$ must be 1, and the fixed cost is paid. If $x = 0$, minimization sets $y = 0$ by itself, so no extra constraint is needed. In a maximization of profit minus costs, the same reasoning applies. The tightest valid choice of $u$ is the smallest upper bound on $x$ implied by the model: market size, capacity or demand, as discussed for big-M.
| Application | Activity $x$ | Fixed charge $K$ |
|---|---|---|
| Lot sizing | production in a period | machine setup |
| Facility location | flow through a depot | opening and running the depot |
| Network design | flow on a link | building or leasing the link |
| Cloud capacity | load on a server | starting the server |
import numpy as np
from scipy.optimize import milp, LinearConstraint, Bounds
# Three products. Profit per unit, setup cost, hours per unit, market limit.
profit = np.array([12, 9, 7])
setup = np.array([400, 250, 150])
hours = np.array([3, 2, 1])
market = np.array([80, 90, 120])
H = 360 # machine hours available
# Variables [x1, x2, x3, y1, y2, y3]
c = -np.concatenate([profit, -setup]) # maximize profit*x - setup*y
A = [np.concatenate([hours, [0, 0, 0]])] # hours constraint
for j in range(3): # x_j <= market_j * y_j (tightest valid M)
row = np.zeros(6); row[j] = 1; row[3 + j] = -market[j]
A.append(row)
cons = LinearConstraint(A, ub=[H, 0, 0, 0])
res = milp(c, constraints=cons, integrality=[0, 0, 0, 1, 1, 1],
bounds=Bounds(0, [np.inf] * 3 + [1] * 3))
x, y = res.x[:3], res.x[3:]
print("produce:", x.round(1) + 0.0, " setups:", y.round().astype(int), " net profit:", -res.fun)
# Naive plan: ignore setup costs, solve the LP, then pay setups for every product made
naive = milp(-profit, constraints=LinearConstraint([hours], ub=[H]), bounds=Bounds(0, market))
xn = naive.x
print("naive :", xn.round(1) + 0.0, " net profit:", profit @ xn - setup @ (xn > 1e-9))
# produce: [ 0. 90. 120.] setups: [0 1 1] net profit: 1250.0
# naive : [ 20. 90. 120.] net profit: 1090.0 <- 20 units of product 1 earn 240, cost 400
The naive plan uses leftover hours for 20 units of product 1. That earns 240 but triggers a 400 setup. The fixed-charge model sees the trade-off and leaves the hours unused.
Key takeaways
- A fixed charge is a cost $K$ paid only when $x > 0$. It is discontinuous, so LP cannot model it.
- Model it with $Ky + cx$ and $x \le u,y$, where $u$ is the tightest valid upper bound on $x$.
- Minimization sets $y = 0$ automatically when $x = 0$. No reverse constraint is needed.
- Fixed charges discourage small, scattered activities. This is the key difference from LP plans.
Practice
- (Beginner) At what setup cost for product 1 would the optimal plan start producing it? Reason with the 60 spare hours, then confirm with the code.
- (Intermediate) Add a minimum batch size: if a product is made, at least 30 units must be produced. Write the constraint with the existing $y_j$ and re-solve.
- (Interview) "A cloud team pays a flat monthly fee per active region plus traffic costs. How would you choose which regions to activate?" Hint: this is a fixed-charge (facility location) model. Name the binaries, the continuous flows and the linking constraints.
Either-Or and If-Then Constraints
Why it matters. Some requirements are disjunctions: "job A runs before job B or B runs before A", "a product is either not made or made in a batch of at least 50", "use either the old process or the new one". The feasible region becomes a union of polyhedra, which is non-convex. Binary variables with big-M constraints are the standard way to model such unions.
Theory.
Either-or. At least one of $g_1(x) \le b_1$ and $g_2(x) \le b_2$ must hold:
$$ g_1(x) \le b_1 + M_1 y, \qquad g_2(x) \le b_2 + M_2 (1 - y), \qquad y \in {0, 1}. $$
With $y = 0$ the first constraint is enforced. With $y = 1$ the second is enforced.
$k$ out of $m$ constraints. Give each constraint a binary $y_i$ that relaxes it, and limit how many can be relaxed:
$$ g_i(x) \le b_i + M_i y_i \quad (i = 1, \dots, m), \qquad \sum_i y_i \le m - k . $$
If-then. "If $f(x) > 0$ then $g(x) \ge 0$" is written as
$$ -g(x) \le M y, \qquad f(x) \le M (1 - y). $$
If $f(x) > 0$, the second row forces $y = 0$, and the first row then gives $g(x) \ge 0$.
Semi-continuous variables. "$x = 0$ or $l \le x \le u$" (a minimum batch size) is $l,y \le x \le u,y$.
Disjunctive scheduling. Two jobs $i, j$ on one machine cannot overlap. With start times $s_i$, processing times $p_i$, and $y_{ij} = 1$ meaning "$i$ before $j$":
$$ s_i + p_i \le s_j + M(1 - y_{ij}), \qquad s_j + p_j \le s_i + M y_{ij}. $$
y_ij = 1: |--- job i ---|--- job j ---| (i before j)
y_ij = 0: |--- job j ---|--- job i ---| (j before i)
forbidden: |--- job i ---|
|--- job j ---| (overlap)
import numpy as np
from itertools import permutations
from scipy.optimize import milp, LinearConstraint, Bounds
# One machine, four jobs: processing time, release date, weight. Minimize sum w_j C_j.
p = np.array([4, 2, 6, 3]); r = np.array([0, 1, 3, 5]); w = np.array([2, 1, 3, 2])
n, M = 4, int(r.max() + p.sum()) # M = latest possible completion time
pairs = [(i, j) for i in range(n) for j in range(i + 1, n)]
# Variables: start times s_0..s_3, then one binary y_ij per pair (1 = i before j)
nv = n + len(pairs)
A, lb, ub = [], [], []
for k, (i, j) in enumerate(pairs):
# either s_i + p_i <= s_j (y = 1) or s_j + p_j <= s_i (y = 0)
row = np.zeros(nv); row[i], row[j], row[n + k] = 1, -1, M
A.append(row); lb.append(-np.inf); ub.append(M - p[i]) # s_i - s_j + M y <= M - p_i
row = np.zeros(nv); row[j], row[i], row[n + k] = 1, -1, -M
A.append(row); lb.append(-np.inf); ub.append(-p[j]) # s_j - s_i - M y <= -p_j
c = np.concatenate([w, np.zeros(len(pairs))]) # sum w_j s_j (+ constant)
res = milp(c, constraints=LinearConstraint(A, lb, ub),
integrality=[0] * n + [1] * len(pairs),
bounds=Bounds(np.concatenate([r, np.zeros(len(pairs))]), # s_j >= r_j
np.concatenate([[M] * n, np.ones(len(pairs))])))
s = res.x[:n]
print("MILP order:", [int(j) + 1 for j in np.argsort(s)], "starts:", s.round().astype(int),
"objective:", round(res.fun + w @ p))
def weighted_completion(order): # brute-force check over all 4! orders
t, total = 0, 0
for j in order:
t = max(t, r[j]) + p[j]
total += w[j] * t
return total
best = min(permutations(range(n)), key=weighted_completion)
print("brute force:", [j + 1 for j in best], "objective:", weighted_completion(best))
# MILP order: [1, 2, 4, 3] starts: [0 4 9 6] objective: 77
# brute force: [1, 2, 4, 3] objective: 77
The objective uses $C_j = s_j + p_j$, so $\sum w_j C_j = \sum w_j s_j + \sum w_j p_j$, and the constant is added back after solving. Job 3 has the largest weight but is released at time 3 and is long, so the optimum runs it last.
Key takeaways
- Either-or: one binary chooses which of two constraints is enforced, and big-M relaxes the other.
- "$k$ of $m$" constraints: one binary per constraint, with at most $m - k$ relaxed.
- If-then constraints and semi-continuous variables follow the same pattern.
- Disjunctive scheduling uses one binary per job pair for the order, which is simple to write but has a weak LP relaxation.
Practice
- (Beginner) Write constraints for "either $x_1 + x_2 \le 10$ or $x_1 - x_2 \ge 4$", given $0 \le x_1, x_2 \le 20$. Choose the smallest valid $M$ values.
- (Intermediate) Add a fifth job ($p = 2$, $r = 2$, $w = 4$) to the scheduling code. How many binaries does the model have now? Confirm the MILP optimum with brute force.
- (Interview) "Disjunctive big-M scheduling models are known to be slow for 50+ jobs. Why, and what would you try instead?" Hint: look at the LP relaxation when all $y_{ij} = 1/2$. Then consider time-indexed formulations and constraint programming (Chapter 11).
Strength of Formulations
Why it matters. The same integer problem can be written in many valid ways, and solver times can differ by orders of magnitude between them. The difference is the LP relaxation: a formulation whose relaxation is close to the integer hull gives tight bounds, prunes branch-and-bound trees early and sometimes solves at the root. Learning to recognize strong formulations is the most valuable skill in integer programming modeling.
Theory. Let $X$ be the set of integer feasible solutions. A formulation of $X$ is any polyhedron $P$ with $P \cap \mathbb{Z}^n = X$. Many polyhedra contain the same lattice points:
- $P_1$ is stronger (tighter) than $P_2$ if $P_1 \subseteq P_2$. Then its LP bound is at least as good.
- The ideal formulation is the convex hull $\operatorname{conv}(X)$. Every vertex is integer, so the LP solves the integer problem directly. The TU matrices of Chapter 3 are ideal formulations.
+---------------------------------------------+ P_weak
| +-------------------------------+ |
| | +-------------------+ | | P_strong
| | | * * * | | |
| | | * * | | | conv(X): ideal formulation,
| | | * * * | | | every vertex is an integer point
| | +-------------------+ | |
| +-------------------------------+ |
+---------------------------------------------+
Same integer points (*), but LP optima of the outer regions lie at fractional corners.
A classic example: facility location. Depots $j$ with opening cost $f_j$, customers $i$ with service cost $c_{ij}$, $x_{ij} = 1$ if $i$ is served by $j$, $y_j = 1$ if $j$ is open. Customers can only use open depots, which can be written in two ways:
| Aggregated | Disaggregated | |
|---|---|---|
| Linking constraints | $\sum_i x_{ij} \le m, y_j$ for each depot | $x_{ij} \le y_j$ for each pair |
| Number of rows | $n$ | $mn$ |
| LP behavior | $y_j = (\text{customers served}) / m$ can be tiny | $y_j \ge \max_i x_{ij}$, so a used depot is nearly fully paid |
| Strength | weak | strong (the disaggregated rows imply the aggregated ones) |
More rows, but a much better formulation. Common ways to strengthen a model:
- Use the smallest valid big-M values.
- Disaggregate linking constraints, as above.
- Add valid inequalities that cut off fractional points (Section 4.2).
- Use extended formulations with extra variables whose projection is tighter.
- Break symmetry: identical machines, bins or vehicles create many equivalent solutions that the solver has to rule out one by one.
import numpy as np
from scipy.optimize import milp, LinearConstraint, Bounds
# Uncapacitated facility location: 4 candidate depots, 8 customers.
f = np.array([40, 30, 35, 38]) # opening cost of depot j
C = np.array([[ 9, 22, 18, 30], # C[i, j] = cost of serving customer i from j
[12, 15, 25, 24],
[20, 8, 30, 19],
[25, 10, 14, 27],
[30, 20, 7, 16],
[28, 26, 11, 9],
[15, 30, 21, 12],
[18, 24, 26, 6]])
m, n = C.shape # m customers, n depots
nv = n + m * n # variables: y_j, then x_ij (row-major)
x = lambda i, j: n + i * n + j
c = np.concatenate([f, C.ravel()])
rows, ub, lb = [], [], []
for i in range(m): # each customer served exactly once
r = np.zeros(nv); r[[x(i, j) for j in range(n)]] = 1
rows.append(r); lb.append(1); ub.append(1)
assign = (rows, lb, ub)
def link_aggregated(): # sum_i x_ij <= m * y_j (n rows)
out = []
for j in range(n):
r = np.zeros(nv); r[[x(i, j) for i in range(m)]] = 1; r[j] = -m
out.append(r)
return out
def link_disaggregated(): # x_ij <= y_j (m*n rows)
out = []
for i in range(m):
for j in range(n):
r = np.zeros(nv); r[x(i, j)] = 1; r[j] = -1
out.append(r)
return out
for name, link in [("aggregated", link_aggregated()), ("disaggregated", link_disaggregated())]:
A = assign[0] + link
cons = LinearConstraint(A, assign[1] + [-np.inf] * len(link), assign[2] + [0] * len(link))
lp = milp(c, constraints=cons, bounds=Bounds(0, 1)) # LP relaxation
ip = milp(c, constraints=cons, bounds=Bounds(0, 1), integrality=np.ones(nv))
print(f"{name:14s} rows = {len(A):2d} LP bound = {lp.fun:6.2f} "
f"IP = {ip.fun:.0f} LP/IP = {lp.fun / ip.fun:.1%} "
f"open = {[j + 1 for j in range(n) if ip.x[j] > 0.5]}")
# aggregated rows = 12 LP bound = 109.12 IP = 166 LP/IP = 65.7% open = [2, 4]
# disaggregated rows = 40 LP bound = 166.00 IP = 166 LP/IP = 100.0% open = [2, 4]
Both formulations have the same integer optimum, but the aggregated LP is off by a third. On this instance the disaggregated LP is already integral, so no branching is needed at all.
Key takeaways
- A formulation is any polyhedron with the right integer points. Smaller polyhedra give better LP bounds.
- The ideal formulation is the convex hull of the integer solutions, where LP alone solves the problem.
- More constraints can mean a much stronger model: compare disaggregated and aggregated facility location.
- Strengthen models with tight big-M values, disaggregation, valid inequalities, extended formulations and symmetry breaking.
Practice
- (Beginner) Show that the constraints $x_{ij} \le y_j$ for all $i$ imply $\sum_i x_{ij} \le m, y_j$. Why does the reverse not hold?
- (Intermediate) For $X = {x \in {0,1}^2 : x_1 + x_2 \le 1}$, compare the formulations ${2x_1 + 2x_2 \le 3,\ 0 \le x \le 1}$ and ${x_1 + x_2 \le 1,\ 0 \le x \le 1}$. Which is ideal? Maximize $x_1 + x_2$ over each.
- (Interview) "Your MILP has 200 identical trucks and solves very slowly. What is going on and what can you do?" Hint: count how many equivalent solutions exist when trucks can be permuted. Then think about ordering constraints such as $y_1 \ge y_2 \ge \dots$, or aggregating identical trucks into one integer variable.
<a id="42-exact-solution-methods"></a>
4.2 Exact Solution Methods
This subchapter explains how solvers find integer optima and prove that they are optimal. Every method combines two ingredients: feasible solutions, which bound the optimum from one side, and relaxations, which bound it from the other. The work consists of closing the gap between them.
LP Relaxation and Bounds
Why it matters. An integer program cannot be checked for optimality the way an LP can: there is no reduced-cost certificate. Instead, solvers bracket the unknown optimum between a feasible solution and a bound. Every exact method in this subchapter, and every gap a commercial solver prints, is built on this idea.
Theory. A relaxation of a maximization problem $\max{f(x) : x \in X}$ is a problem $\max{g(x) : x \in R}$ with $X \subseteq R$ and $g(x) \ge f(x)$ on $X$. Its optimum is an upper bound on the original optimum. The most common one is the LP relaxation: drop integrality and keep everything else.
For a maximization problem:
| Quantity | Source | Role |
|---|---|---|
| Primal bound $\underline{z}$ | value of the best feasible integer solution found (the incumbent) | lower bound on $z^*$ |
| Dual bound $\overline{z}$ | optimum of a relaxation (LP, Lagrangian, ...) | upper bound on $z^*$ |
| Optimality gap | $(\overline{z} - \underline{z}) / \lvert \underline{z} \rvert$ | how far the incumbent can be from optimal |
| Integrality gap | $z_{LP} / z_{IP}$ (worst case over instances) | how good a formulation's relaxation can be |
For minimization problems the roles flip: feasible solutions give upper bounds and relaxations give lower bounds.
value
16.8 ---- LP relaxation dual bound: no solution can beat this
16 ---- floor(16.8) integer objective coefficients -> optimum is an integer
16 .... true optimum (unknown until the gap closes)
15 ---- greedy solution primal bound: we already own a solution this good
Two useful refinements:
- If all objective coefficients are integers and all variables are integer, the optimum is an integer, so the bound can be rounded down to $\lfloor \overline{z} \rfloor$ (for maximization).
- When the gap reaches zero, the incumbent is proven optimal. Stopping at a gap of 1% is standard practice for large models: the solution is guaranteed to be within 1% of optimal.
Other relaxations drop or soften different constraints: the Lagrangian relaxation (later in this subchapter), or combinatorial relaxations such as the assignment relaxation of the TSP (Section 4.3).
import numpy as np
from scipy.optimize import linprog
# Project selection from Section 1.1: (cost, value), budget 10. Optimum (by enumeration) = 16.
cost = np.array([4, 3, 5, 2, 6])
value = np.array([7, 5, 8, 3, 8])
budget = 10
# Dual (upper) bound: LP relaxation, 0 <= x_j <= 1
lp = linprog(-value, A_ub=[cost], b_ub=[budget], bounds=[(0, 1)] * 5)
print("LP relaxation x =", lp.x.round(2) + 0.0, " upper bound =", -lp.fun)
# Primal (lower) bound: greedy by value/cost ratio, a feasible integer solution
x, room = np.zeros(5, dtype=int), budget
for j in np.argsort(-value / cost): # best ratio first
if cost[j] <= room:
x[j], room = 1, room - cost[j]
incumbent = value @ x
print("greedy x =", x, " incumbent =", incumbent)
gap = (-lp.fun - incumbent) / incumbent
print(f"optimality gap <= {gap:.1%} -> the optimum lies in [{incumbent}, {-lp.fun}]")
print(f"integer bound: values are integers, so optimum <= floor({-lp.fun}) = {int(np.floor(-lp.fun))}")
# LP relaxation x = [1. 1. 0.6 0. 0. ] upper bound = 16.8
# greedy x = [1 1 0 1 0] incumbent = 15
# optimality gap <= 12.0% -> the optimum lies in [15, 16.8]
# integer bound: values are integers, so optimum <= floor(16.8) = 16
Without knowing the optimum, we already know the greedy solution is within 12% of it. Finding any solution with value 16 would prove optimality immediately.
Key takeaways
- Feasible solutions give primal bounds. Relaxations give dual bounds. The optimum lies between them.
- The LP relaxation is the standard dual bound. Its quality depends on the formulation (Section 4.1).
- The optimality gap measures how far the incumbent may be from optimal. A zero gap proves optimality.
- With integer objective coefficients, a dual bound can be rounded to the nearest integer in the right direction.
Practice
- (Beginner) In a minimization problem, a heuristic finds a solution of cost 230 and the LP relaxation gives 212. What is the optimality gap, and what can you say about the optimum?
- (Intermediate) Show that the LP relaxation of a 0-1 knapsack is solved by the greedy ratio order with one fractional item (Dantzig's bound). Use it to compute the bound for the code's instance by hand.
- (Interview) "A solver has run for an hour and reports a 0.8% gap. Your manager asks whether the solution is good enough. What do you answer?" Hint: the gap is a guarantee, not an estimate. Relate 0.8% to the business value and to the uncertainty in the input data.
Branch and Bound
Why it matters. Branch and bound (Land and Doig, 1960) is the backbone of every MILP solver. It explores the exponential space of integer solutions while proving that most of it cannot contain the optimum, so it rarely needs to enumerate more than a tiny fraction.
Theory. Branch and bound is divide and conquer with bounds. Picture searching a huge library for the single best book: you split the shelves into sections, and if the catalog says that even the best book in a section is worse than one you already hold, you skip the whole section.
- Bound. Solve the LP relaxation of the current subproblem (node).
- Prune the node if:
- by infeasibility: the LP has no solution,
- by bound: the LP value is no better than the incumbent,
- by integrality: the LP solution is integer. It becomes the new incumbent if it is better.
- Branch otherwise: pick a variable $x_j$ with fractional value $x_j^$ and create two children with $x_j \le \lfloor x_j^ \rfloor$ and $x_j \ge \lceil x_j^* \rceil$. Every integer solution of the node survives in exactly one child, and $x^*$ is excluded from both.
- Repeat until no open nodes remain. The incumbent is then optimal.
The global dual bound at any moment is the best LP value among the open nodes, so the gap can be reported throughout.
| Design choice | Options | Trade-off |
|---|---|---|
| Node selection | depth-first | finds incumbents fast, low memory, but can explore poor regions |
| best-bound first | improves the global bound fastest, needs more memory | |
| hybrid (dive, then best bound) | what modern solvers do | |
| Branching variable | most fractional | simple, but not much better than random in practice |
| pseudocost | learns from past branchings how much each variable moves the bound | |
| strong branching | tentatively solves child LPs for candidates, giving small trees at high cost per node | |
| reliability branching | strong branching until pseudocosts are reliable, then pseudocosts |
The tree for the running example, with depth-first search exploring the "up" branch first:
root: z = 41.25 at (2.25, 3.75)
/ \
x1 >= 3 x1 <= 2
z = 39 at (3, 3), integral z = 41.11 at (2, 3.889)
-> incumbent 39 / \
x2 >= 4 x2 <= 3
z = 41 at (1.8, 4) z = 34 <= 40
/ \ -> prune (bound)
x1 >= 2 x1 <= 1
infeasible z = 40.56 at (1, 4.444)
-> prune / \
x2 >= 5 x2 <= 4
z = 40 at (0, 5) z = 37 <= 40
-> incumbent 40 -> prune (bound)
import math
import numpy as np
from scipy.optimize import linprog
# max 5 x1 + 8 x2 s.t. x1 + x2 <= 6, 5 x1 + 9 x2 <= 45, x >= 0 integer
c = np.array([5, 8]); A = [[1, 1], [5, 9]]; b = [6, 45]
best = {"value": -math.inf, "x": None}
nodes = 0
def branch_and_bound(bounds, depth=0, label="root"):
global nodes
nodes += 1
pad = " " * depth
lp = linprog(-c, A_ub=A, b_ub=b, bounds=bounds)
if lp.status != 0: # prune: infeasible
print(f"{pad}{label}: infeasible -> prune"); return
z, x = -lp.fun, lp.x
if z <= best["value"] + 1e-9: # prune: bound
print(f"{pad}{label}: z = {z:.2f} <= incumbent {best['value']:.0f} -> prune"); return
frac = [j for j in range(len(x)) if abs(x[j] - round(x[j])) > 1e-6]
if not frac: # prune: integral, new incumbent
best["value"], best["x"] = z, x.round()
print(f"{pad}{label}: z = {z:.2f} at {x.round() + 0.0} integral -> new incumbent"); return
j = min(frac, key=lambda k: abs(x[k] - math.floor(x[k]) - 0.5)) # most fractional
print(f"{pad}{label}: z = {z:.2f} at {x.round(3) + 0.0} -> branch on x{j + 1}")
lo, hi = bounds[j]
down = list(bounds); down[j] = (lo, math.floor(x[j]))
up = list(bounds); up[j] = (math.ceil(x[j]), hi)
branch_and_bound(up, depth + 1, f"x{j + 1} >= {math.ceil(x[j])}")
branch_and_bound(down, depth + 1, f"x{j + 1} <= {math.floor(x[j])}")
branch_and_bound([(0, None), (0, None)])
print(f"optimum x = {best['x'] + 0.0}, value = {best['value']:.0f}, nodes explored = {nodes}")
# root: z = 41.25 at [2.25 3.75] -> branch on x1
# x1 >= 3: z = 39.00 at [3. 3.] integral -> new incumbent
# x1 <= 2: z = 41.11 at [2. 3.889] -> branch on x2
# x2 >= 4: z = 41.00 at [1.8 4. ] -> branch on x1
# x1 >= 2: infeasible -> prune
# x1 <= 1: z = 40.56 at [1. 4.444] -> branch on x2
# x2 >= 5: z = 40.00 at [0. 5.] integral -> new incumbent
# x2 <= 4: z = 37.00 <= incumbent 40 -> prune
# x2 <= 3: z = 34.00 <= incumbent 40 -> prune
# optimum x = [0. 5.], value = 40, nodes explored = 9
Real solvers re-optimize each child from its parent's basis with the dual simplex method (Section 2.4): only one bound changed, so a few pivots usually suffice. Rounding the bound can also prune earlier. The objective coefficients are integers, so a node with $z = 40.56$ cannot contain anything better than 40, and once an incumbent of 40 is known, such nodes are pruned without branching.
Key takeaways
- Branch and bound splits the problem on a fractional variable and solves LP relaxations to bound each part.
- Nodes are pruned by infeasibility, by bound (no better than the incumbent) or by integrality.
- Node selection trades fast incumbents (depth-first) against fast bound improvement (best-bound).
- Good branching rules (pseudocost, strong, reliability branching) can shrink trees by orders of magnitude.
Practice
- (Beginner) Redraw the tree when the "down" branch $x_j \le \lfloor x_j^* \rfloor$ is explored first. Does the number of nodes change?
- (Intermediate) Modify the code to use best-bound node selection with a priority queue (
heapq) and to prune with $\lfloor z \rfloor$. Compare node counts on the running example and on the 5-project knapsack.- (Interview) "Branch and bound is exponential in the worst case. Why do solvers still solve problems with millions of binaries?" Hint: think about pruning by bound, strong formulations, presolve, cuts and primal heuristics, and how far typical instances are from the worst case.
Cutting Plane Methods
Why it matters. Branching splits the problem. Cutting planes instead tighten it: they add inequalities that every integer solution satisfies but the current fractional LP solution violates. Cuts improve the LP bound without growing a tree. Since the 1990s they have been the single largest contributor to the speed of MILP solvers.
Theory. An inequality $\pi^\top x \le \pi_0$ is valid for the integer set $X$ if every $x \in X$ satisfies it. A cut is a valid inequality violated by the current LP optimum $x^*$. Think of a sculptor: the LP polyhedron is the rough block, and each cut chips off a corner that contains no integer points, moving the block toward the shape of the integer hull.
Cutting plane loop
solve LP relaxation -> x*
|
v
x* integer? --- yes ---> optimal, stop
|
no
v
separation: find a valid inequality violated by x*
|
found? --- no ---> stop cutting (branch, Section "Branch and Cut")
|
yes
v
add the cut to the LP, re-solve with dual simplex --> back to the top
Finding a violated inequality from a given family is the separation problem. Its difficulty depends on the family:
| Cut family | Derived from | Typical use |
|---|---|---|
| Gomory / Chvátal-Gomory | any row of the simplex tableau | general, every MILP |
| Mixed integer rounding (MIR) | combinations of rows with continuous variables | general MILP, lot sizing |
| Knapsack cover | a single knapsack row $\sum a_j x_j \le b$ | capacity and budget constraints |
| Clique | conflict graph of binaries ($x_i + x_j \le 1$ pairs) | scheduling, set packing |
| Flow cover | fixed-charge network flow rows | network design |
| Subtour elimination | connectivity of tours | TSP, vehicle routing |
Knapsack cover inequalities. For $\sum_j w_j x_j \le W$ with binary $x$, a set $C$ is a cover if $\sum_{j \in C} w_j > W$: all items of $C$ cannot fit together. Then at most $\lvert C \rvert - 1$ of them can be chosen:
$$ \sum_{j \in C} x_j \le \lvert C \rvert - 1 . $$
import numpy as np
from itertools import combinations
from scipy.optimize import linprog, milp, LinearConstraint, Bounds
# 0-1 knapsack: max v x s.t. w x <= 20, x binary
v = np.array([20, 16, 13, 11, 8, 5])
w = np.array([12, 9, 8, 7, 5, 4])
cap, n = 20, len(v)
cuts_A, cuts_b = [], []
def separate(x):
"""Most violated cover inequality sum_{j in C} x_j <= |C| - 1 (exact, by enumeration)."""
best, best_viol = None, 1e-6
for k in range(2, n + 1):
for C in combinations(range(n), k):
if w[list(C)].sum() > cap: # C is a cover: cannot all fit
viol = x[list(C)].sum() - (k - 1)
if viol > best_viol:
best, best_viol = C, viol
return best
for it in range(1, 10):
A = [w] + cuts_A; b = [cap] + cuts_b
lp = linprog(-v, A_ub=A, b_ub=b, bounds=[(0, 1)] * n)
C = separate(lp.x)
print(f"round {it}: LP bound = {-lp.fun:6.2f} x = {lp.x.round(2) + 0.0}")
if C is None:
break
row = np.zeros(n); row[list(C)] = 1
cuts_A.append(row); cuts_b.append(len(C) - 1)
print(f" add cover cut " + " + ".join(f"x{j + 1}" for j in C) + f" <= {len(C) - 1}")
ip = milp(-v, constraints=LinearConstraint([w], ub=[cap]), integrality=np.ones(n), bounds=Bounds(0, 1))
print("integer optimum =", -ip.fun, "with x =", ip.x.round().astype(int))
# round 1: LP bound = 34.33 x = [0.92 1. 0. 0. 0. 0. ]
# add cover cut x1 + x2 <= 1
# round 2: LP bound = 33.80 x = [0. 1. 1. 0. 0.6 0. ]
# add cover cut x2 + x3 + x5 <= 2
# round 3: LP bound = 33.71 x = [0. 1. 1. 0.43 0. 0. ]
# add cover cut x2 + x3 + x4 <= 2
# round 4: LP bound = 33.60 x = [0.2 0.8 1. 0.2 0.2 0. ]
# integer optimum = 33.0 with x = [1 0 1 0 0 0]
The cuts lower the bound from 34.33 to 33.60. Then no violated cover exists, although $x^*$ is still fractional. This is typical: one cut family alone rarely reaches the integer hull, and the bound improvement per round tails off. Solvers therefore combine several families and stop cutting when progress slows, then branch.
Key takeaways
- A cut is valid for all integer solutions and violated by the current LP optimum.
- The separation problem finds a violated cut from a given family, or reports that none exists.
- Cover inequalities say "not all items of an overfull set fit". They strengthen every knapsack-type row.
- Pure cutting plane methods tail off. In practice cuts are combined with branching.
Practice
- (Beginner) For the constraint $5x_1 + 7x_2 + 4x_3 + 3x_4 \le 14$ with binary $x$, list all minimal covers (covers that stop being covers if any item is removed).
- (Intermediate) A cover inequality can be lifted: items outside $C$ get positive coefficients while staying valid. For the cover ${2, 3, 5}$ of the example, find the largest $\alpha$ such that $\alpha x_1 + x_2 + x_3 + x_5 \le 2$ is valid, and add it to the code.
- (Interview) "If cuts improve the bound without branching, why not just add cuts until the LP is integral?" Hint: think about tailing off, the size and density of the LP after many cuts, numerical accuracy, and the cost of separation.
Gomory Cuts
Why it matters. Gomory's fractional cut (1958) was the first general-purpose cutting plane: it can be read off any row of the optimal simplex tableau, for any pure integer program, with no problem-specific knowledge. Its mixed integer variant (GMI) is still among the most effective cuts in every modern solver.
Theory. Assume all variables (including slacks) must be integer, which holds when $A$ and $b$ have integer entries. Take a row of the optimal tableau whose basic variable has a fractional value:
$$ x_{B(i)} + \sum_{j \in N} \bar a_{ij} x_j = \bar b_i . $$
Write each number as integer part plus fractional part: $\bar a_{ij} = \lfloor \bar a_{ij} \rfloor + f_{ij}$ and $\bar b_i = \lfloor \bar b_i \rfloor + f_i$, with $0 \le f < 1$.
- Since $x \ge 0$, rounding coefficients down gives $x_{B(i)} + \sum_j \lfloor \bar a_{ij} \rfloor x_j \le \bar b_i$.
- The left side is an integer for integer $x$, so it is at most $\lfloor \bar b_i \rfloor$. This is the Chvátal-Gomory cut.
- Subtracting it from the tableau row gives the Gomory fractional cut:
$$ \sum_{j \in N} f_{ij} x_j \ge f_i . $$
At the current LP vertex every nonbasic $x_j$ is 0, so the left side is 0 while $f_i > 0$: the vertex is cut off. Every integer solution satisfies the cut.
Worked example. The running example has the optimal tableau row $x_2 - \tfrac{5}{4}s_1 + \tfrac{1}{4}s_2 = \tfrac{15}{4}$. The fractional parts are $f(-\tfrac{5}{4}) = \tfrac{3}{4}$, $f(\tfrac{1}{4}) = \tfrac{1}{4}$ and $f(\tfrac{15}{4}) = \tfrac{3}{4}$, so the cut is
$$ \tfrac{3}{4}s_1 + \tfrac{1}{4}s_2 \ge \tfrac{3}{4}. $$
Substituting $s_1 = 6 - x_1 - x_2$ and $s_2 = 45 - 5x_1 - 9x_2$ gives $2x_1 + 3x_2 \le 15$ in the original variables. The line passes through the integer points $(0, 5)$ and $(3, 3)$ and cuts off $(2.25, 3.75)$:
x2
5 | *. cut 2 x1 + 3 x2 <= 15 (dots)
| .. o <- (2.25, 3.75), now infeasible
4 | * * ..
3 | * * * *.
2 | * * * * *.
1 | * * * * * *
0 | * * * * * * *
+------------------------ x1
0 1 2 3 4 5 6
from fractions import Fraction as F
import math
# Gomory fractional cutting planes on: max 5 x1 + 8 x2
# s.t. x1 + x2 + s1 = 6, 5 x1 + 9 x2 + s2 = 45, all variables >= 0 and integer
# Tableau rows: [coefficients..., rhs]; obj row stores reduced costs (max problem).
names = ["x1", "x2", "s1", "s2"]
T = [[F(1), F(1), F(1), F(0), F(6)],
[F(5), F(9), F(0), F(1), F(45)]]
obj = [F(-5), F(-8), F(0), F(0), F(0)] # z - 5x1 - 8x2 = 0
basis = [2, 3]
def pivot(r, k):
piv = T[r][k]
T[r] = [a / piv for a in T[r]]
for i in range(len(T)):
if i != r and T[i][k] != 0:
T[i] = [a - T[i][k] * b for a, b in zip(T[i], T[r])]
f = obj[k]
obj[:] = [a - f * b for a, b in zip(obj, T[r])]
basis[r] = k
def primal_simplex():
while min(obj[:-1]) < 0:
k = min(range(len(obj) - 1), key=lambda j: obj[j]) # most negative reduced cost
r = min((i for i in range(len(T)) if T[i][k] > 0), key=lambda i: T[i][-1] / T[i][k])
pivot(r, k)
def dual_simplex():
while any(row[-1] < 0 for row in T):
r = min(range(len(T)), key=lambda i: T[i][-1]) # most negative rhs leaves
k = min((j for j in range(len(obj) - 1) if T[r][j] < 0),
key=lambda j: obj[j] / -T[r][j]) # dual ratio test
pivot(r, k)
def frac(a):
return a - math.floor(a)
def show():
sol = {names[basis[i]]: T[i][-1] for i in range(len(T))}
print(f" LP optimum z = {float(obj[-1]):.4f}, x1 = {sol.get('x1', 0)}, x2 = {sol.get('x2', 0)}")
primal_simplex(); show()
for it in range(1, 10):
rows = [i for i in range(len(T)) if frac(T[i][-1]) != 0]
if not rows:
print(" all basic variables integral -> integer optimum found"); break
r = max(rows, key=lambda i: frac(T[i][-1])) # most fractional source row
cut = [-frac(a) for a in T[r][:-1]] + [F(1), -frac(T[r][-1])] # -sum f_j x_j + t = -f_0
terms = " + ".join(f"{frac(a)}*{names[j]}" for j, a in enumerate(T[r][:-1]) if frac(a) != 0)
print(f"cut {it} from row of {names[basis[r]]}: {terms} >= {frac(T[r][-1])}")
for row in T:
row.insert(-1, F(0)) # new slack column
obj.insert(-1, F(0))
T.append(cut); names.append(f"t{it}"); basis.append(len(names) - 1)
dual_simplex(); show()
# LP optimum z = 41.2500, x1 = 9/4, x2 = 15/4
# cut 1 from row of x2: 3/4*s1 + 1/4*s2 >= 3/4
# LP optimum z = 40.0000, x1 = 0, x2 = 5
# all basic variables integral -> integer optimum found
Exact fractions matter here. With floating-point arithmetic, a coefficient such as $0.9999999$ produces a nearly meaningless cut, which is why early implementations were unstable. Here one cut is enough. In general, Gomory proved that a careful (lexicographic) version terminates after finitely many cuts, but the number can be huge. Solvers add a limited number of GMI cuts at the root and then branch.
Key takeaways
- A Gomory cut comes from any fractional row of the optimal tableau: $\sum_j f_{ij} x_j \ge f_i$.
- It is valid for all integer solutions and always cuts off the current LP vertex.
- In original variables, the first cut of the running example is $2x_1 + 3x_2 \le 15$, and it makes the LP integral.
- Gomory cuts need exact or carefully safeguarded arithmetic. GMI cuts are standard in modern solvers.
Practice
- (Beginner) From the tableau row $x_1 + \tfrac{1}{3}s_1 - \tfrac{2}{3}s_2 = \tfrac{7}{3}$, derive the Gomory fractional cut.
- (Intermediate) Change the right-hand side 45 to 44 in the code and run it again. How many cuts are needed? Translate each cut into the original variables $x_1, x_2$.
- (Interview) "Gomory cuts were invented in 1958 but considered impractical for 30 years. What changed?" Hint: think about numerical safeguards, adding cuts in rounds rather than one at a time, and combining them with branch and bound instead of using them alone.
Branch and Cut
Why it matters. Branch and cut is what modern MILP solvers (Gurobi, CPLEX, HiGHS, SCIP) actually run. It combines the strengths of both methods: cuts tighten the bound at each node, and branching guarantees progress when cuts tail off. It is also the only practical way to handle models with exponentially many constraints, such as the TSP.
Theory. Branch and cut is branch and bound in which each node runs a cutting plane loop before branching:
presolve (tighten bounds, remove redundancy, detect structure)
|
v
+--> select open node --> solve LP --> prune? ---- yes --> next node
| ^ |
| | no
| add cuts | v
| +--- separate cuts (rounds, while bound improves enough)
| |
| primal heuristics (rounding, diving, RINS ...)
| |
+------------------------------- branch: create children
Key design elements:
| Element | Idea |
|---|---|
| Global vs local cuts | Cuts valid for the whole problem go into a cut pool shared by all nodes. Cuts that rely on branching bounds are only valid in the subtree. |
| Cut management | Too many cuts make LPs slow and dense. Solvers age out inactive cuts and keep only effective ones. |
| Root emphasis | Most cutting happens at the root, where a better bound helps the whole tree. |
| Lazy constraints | Constraints too numerous to list (subtour elimination) are checked only when an integer candidate appears, and added if violated. |
| Primal heuristics | Rounding, diving and local search find incumbents early, so pruning by bound starts sooner. |
Reading a solver log. Every MILP solver prints the same core columns: nodes explored, open nodes, incumbent (primal bound), best bound (dual bound), and gap. The run is over when the gap closes or a limit is hit.
The code below adds cover cuts (with a global cut pool) to a simple branch and bound on a 14-item, 2-resource knapsack and compares tree sizes on three random instances:
import math
import numpy as np
from scipy.optimize import linprog
n = 14 # items; 2 resources (weight, volume)
def instance(seed):
rng = np.random.default_rng(seed)
W = rng.integers(10, 40, size=(2, n)) # resource use per item
cap = (W.sum(axis=1) * 0.45).astype(int) # about 45% of the total
v = (W.sum(axis=0) * rng.uniform(0.8, 1.2, n)).round()
return W, cap, v
def cover_cuts(x):
"""Greedy cover separation for each resource row: one cut per row if violated."""
cuts = []
for k in range(2):
order = sorted(range(n), key=lambda j: (1 - x[j]) / W[k, j]) # cheap-to-include first
C, load = [], 0
for j in order:
C.append(j); load += W[k, j]
if load > cap[k]:
break
if load > cap[k] and x[C].sum() > len(C) - 1 + 1e-6:
row = np.zeros(n); row[C] = 1
cuts.append((row, len(C) - 1))
return cuts
def solve(use_cuts):
best, nodes, pool = [-math.inf], [0], [] # pool: global cut pool (A row, rhs)
def node(bounds):
nodes[0] += 1
for _ in range(20 if use_cuts else 1): # cut loop at this node
A = list(W) + [r for r, _ in pool]; b = list(cap) + [r for _, r in pool]
lp = linprog(-v, A_ub=A, b_ub=b, bounds=bounds)
if lp.status != 0 or -lp.fun <= best[0] + 1e-9:
return # prune by infeasibility or bound
new = cover_cuts(lp.x) if use_cuts else []
if not new:
break
pool.extend(new) # cuts are globally valid: keep them
x = lp.x
frac = [j for j in range(n) if 1e-6 < x[j] < 1 - 1e-6]
if not frac:
best[0] = -lp.fun; return # integral: new incumbent
j = min(frac, key=lambda k: abs(x[k] - 0.5))
for val in (1, 0): # branch x_j = 1 first, then x_j = 0
nb = list(bounds); nb[j] = (val, val)
node(nb)
node([(0, 1)] * n)
return best[0], nodes[0], len(pool)
for seed in (1, 2, 7):
W, cap, v = instance(seed)
(z0, n0, _), (z1, n1, k1) = solve(False), solve(True)
print(f"seed {seed}: optimum {z0:.0f} = {z1:.0f} B&B nodes = {n0:3d} "
f"B&C nodes = {n1:3d} ({k1} cover cuts)")
# seed 1: optimum 346 = 346 B&B nodes = 327 B&C nodes = 157 (148 cover cuts)
# seed 2: optimum 324 = 324 B&B nodes = 303 B&C nodes = 385 (348 cover cuts)
# seed 7: optimum 343 = 343 B&B nodes = 573 B&C nodes = 289 (248 cover cuts)
Cuts halve the tree on two instances and enlarge it on the third. Even with fewer nodes, every LP is larger. This is why solvers tune cut aggressiveness carefully and measure success in total time, not in node counts.
Key takeaways
- Branch and cut runs a cutting plane loop at each node of branch and bound.
- Globally valid cuts are stored in a pool and reused across the tree. Local cuts apply only to a subtree.
- Lazy constraints let a model have exponentially many rows, checked only when needed.
- Cuts usually shrink trees but cost LP time. Their net effect is measured empirically.
Practice
- (Beginner) Explain why a cover inequality found at a deep node of the tree is still valid at the root, while a Gomory cut derived at that node in general is not.
- (Intermediate) Change the code so cuts are separated only at the root node. Compare node counts and the number of cuts with the version above on seeds 1, 2 and 7.
- (Interview) "How would you decide whether a family of problem-specific cuts is worth implementing in a production solver?" Hint: measure on a benchmark set: time to optimality, root gap closed and separation time, not just node count.
Branch and Price
Why it matters. Some of the strongest formulations in OR have one variable per pattern: a cutting pattern, a crew pairing, a vehicle route. These formulations have exponentially many variables, so they can only be solved with column generation (Section 2.6). Branch and price embeds column generation in branch and bound, and it is the method behind state-of-the-art solvers for crew scheduling, vehicle routing and cutting stock.
Theory. At every node of the branch-and-bound tree, the LP relaxation is solved by column generation: a restricted master problem with few columns, plus a pricing problem that generates columns with negative reduced cost (Section 2.6). The node's bound is valid only once pricing finds no improving column.
The branching difficulty. Branching on a master variable, such as $x_p \le \lfloor x_p^* \rfloor$, is a bad idea:
- In the down branch the pricing problem would keep regenerating the forbidden pattern $p$. Excluding it changes the pricing problem from a clean knapsack or shortest path into "best pattern except these", which is much harder.
- The resulting tree is very unbalanced: the down branch removes only one of exponentially many patterns.
Instead, branch on decisions the pricing problem can respect:
| Problem | Branching rule | Effect on pricing |
|---|---|---|
| Set partitioning (crew scheduling) | Ryan-Foster: rows $r, s$ covered by the same column, or by different columns | forbid or force the pair inside generated columns |
| Vehicle routing | branch on the original arc variable $x_{ij} = \sum_p a_{ij}^p \lambda_p$ | remove arc $(i, j)$ or force it in the pricing graph |
| Cutting stock | branch on the number of pieces of a width cut in some pattern position | add bounds to the knapsack |
Branch and price
node --> column generation until no column prices out --> LP bound for the node
|
prune? (bound, infeasible, integral) -- yes --> next node
|
no
v
branch on ORIGINAL decisions (Ryan-Foster, arcs)
and pass the restriction down to the pricing problem
A cheaper alternative is price and branch: run column generation only at the root, then solve the master as an integer program over the columns generated so far. It is a heuristic, because the optimal integer solution may need columns that were never generated. It becomes exact when its value matches the rounded root LP bound. The cutting stock instance of Section 2.6 shows this:
import math
import numpy as np
from scipy.optimize import linprog, milp, LinearConstraint, Bounds
# Cutting stock instance of Section 2.6: rolls of width 100
W = 100
widths = np.array([45, 36, 31, 14])
demand = np.array([97, 610, 395, 211])
def best_pattern(y):
"""Pricing: unbounded integer knapsack max y.a s.t. widths.a <= W (DP over capacity)."""
best = np.zeros(W + 1); take = [-1] * (W + 1)
for c in range(1, W + 1):
best[c], take[c] = best[c - 1], -1
for i, w in enumerate(widths):
if w <= c and best[c - w] + y[i] > best[c] + 1e-12:
best[c], take[c] = best[c - w] + y[i], i
a, c = np.zeros(len(widths), dtype=int), W
while c > 0:
if take[c] == -1: c -= 1
else: a[take[c]] += 1; c -= widths[take[c]]
return best[W], a
def column_generation(patterns):
while True:
P = np.array(patterns).T
rmp = linprog(np.ones(P.shape[1]), A_ub=-P, b_ub=-demand)
value, a = best_pattern(-rmp.ineqlin.marginals)
if 1 - value >= -1e-9: # no column prices out: LP optimal
return rmp, P
patterns.append(a)
patterns = [np.eye(4, dtype=int)[i] * (W // widths[i]) for i in range(4)]
rmp, P = column_generation(patterns)
print(f"root LP bound = {rmp.fun:.2f} -> any plan needs >= {math.ceil(rmp.fun - 1e-9)} rolls")
# Price-and-branch heuristic: solve the integer master over the columns generated so far
ip = milp(np.ones(P.shape[1]), constraints=LinearConstraint(P, lb=demand),
integrality=np.ones(P.shape[1]), bounds=Bounds(0, np.inf))
print(f"integer master on {P.shape[1]} columns: {ip.fun:.0f} rolls")
for k in np.flatnonzero(ip.x > 0.5):
print(f" {ip.x[k]:4.0f} rolls cut as {P[:, k]}")
print("proven optimal" if ip.fun <= math.ceil(rmp.fun - 1e-9) else "gap remains -> branch")
# root LP bound = 452.25 -> any plan needs >= 453 rolls
# integer master on 6 columns: 453 rolls
# 49 rolls cut as [2 0 0 0]
# 206 rolls cut as [0 2 0 2]
# 198 rolls cut as [0 1 2 0]
# proven optimal
The integer master reaches 453 rolls, which equals the rounded-up root bound, so it is optimal and no branching is needed. Cutting stock instances very often have this property (the "integer round-up property"). For crew scheduling and routing, the gap is usually larger, and full branch and price is required.
Key takeaways
- Branch and price solves the LP at every node by column generation. A node bound is valid only after pricing finds no improving column.
- Branch on original decisions (Ryan-Foster pairs, arcs), not on master variables, so the pricing problem keeps its structure.
- Price and branch (root columns only, then an integer master) is a fast heuristic, and exact when it matches the rounded root bound.
- Applications: crew scheduling, vehicle routing with time windows, cutting stock, generalized assignment.
Practice
- (Beginner) Check that the three patterns in the output fit in a roll of width 100 and that the plan meets all four demands.
- (Intermediate) Change the demands to $(44, 3, 48, 2)$ and re-run. Does price and branch still prove optimality? Describe what a full branch-and-price code would have to do if the integer master had returned one roll more than the rounded bound.
- (Interview) "In a branch-and-price code for crew scheduling, why can't you just fix a fractional pairing variable to 0 in the down branch?" Hint: think about what the pricing problem (a shortest path) would return next, and about how many pairings exist.
Lagrangian Relaxation
Why it matters. Many hard problems are easy problems with a few complicating constraints added: assignment plus capacities, shortest path plus a resource limit, independent subproblems plus a shared budget. Lagrangian relaxation removes the complicating constraints, charges a price for violating them, and solves what is left very quickly. It gives bounds that can be stronger than the LP bound, and it underlies many large-scale heuristics.
Theory. Consider a minimization problem with easy constraints $x \in X$ and complicating constraints $Ax \le b$:
$$ z_{IP} = \min { c^\top x : Ax \le b,\ x \in X }. $$
Move the complicating constraints into the objective with multipliers $\lambda \ge 0$:
$$ L(\lambda) = \min_{x \in X} \ c^\top x + \lambda^\top (Ax - b). $$
For any feasible $x$, the term $\lambda^\top(Ax - b)$ is $\le 0$, so $L(\lambda) \le z_{IP}$: every $\lambda$ gives a lower bound. Think of $\lambda_i$ as a fee per unit of resource $i$. Instead of a hard limit, each unit used costs money, and the fees are adjusted until the decisions respect the limits, like prices that clear a market.
The best bound is the Lagrangian dual $z_{LD} = \max_{\lambda \ge 0} L(\lambda)$. For equality constraints the multipliers are free in sign. A key theorem (Geoffrion, 1974) compares it with the LP bound:
$$ z_{LP} \le z_{LD} \le z_{IP}, $$
with $z_{LD} = z_{LP}$ when the subproblem over $X$ has the integrality property (its LP relaxation already has integer optima). So relax the constraints whose removal leaves a subproblem that is easy but not trivially integral. That is where Lagrangian relaxation beats LP.
Subgradient optimization. $L(\lambda)$ is concave and piecewise linear, so it is maximized with subgradients. If $x^\lambda$ solves the subproblem, then $g = Ax^\lambda - b$ is a subgradient, and the update is
$$ \lambda \leftarrow \max{0,\ \lambda + t, g}, \qquad t = \theta , \frac{UB - L(\lambda)}{\lVert g \rVert^2} \ \ \text{(Polyak step, with } UB \text{ from a heuristic)}. $$
Constraints that are violated ($g_i > 0$) get a higher price, and slack ones get a lower price.
Example: generalized assignment. Six jobs are assigned to three machines with capacities. Two relaxations are possible:
| Relaxed constraints | Remaining subproblem | Integrality property? | Expected bound |
|---|---|---|---|
| capacities ($\lambda \ge 0$) | each job picks its cheapest machine | yes | $= z_{LP}$ |
| assignments ($\mu$ free) | one 0-1 knapsack per machine | no | can be $> z_{LP}$ |
import numpy as np
from itertools import product
from scipy.optimize import milp, LinearConstraint, Bounds
# Generalized assignment: 6 jobs, 3 machines. Every job goes to one machine, loads <= capacity.
cost = np.array([[9, 2, 7], [6, 4, 3], [5, 8, 1], [7, 6, 9], [3, 10, 6], [8, 4, 5]]) # [job, machine]
load = np.array([[5, 8, 6], [7, 4, 5], [4, 6, 9], [6, 5, 7], [9, 7, 4], [3, 6, 8]])
cap = np.array([11, 12, 10])
J, I = cost.shape
# Reference values: LP relaxation and integer optimum (variables x[j, i] row-major)
A_assign = np.kron(np.eye(J), np.ones(I)) # sum_i x_ji = 1
A_cap = np.hstack([np.diag(load[j]) for j in range(J)]) # sum_j load_ji x_ji <= cap_i
cons = [LinearConstraint(A_assign, 1, 1), LinearConstraint(A_cap, ub=cap)]
lp = milp(cost.ravel(), constraints=cons, bounds=Bounds(0, 1))
ip = milp(cost.ravel(), constraints=cons, bounds=Bounds(0, 1), integrality=np.ones(J * I))
print(f"LP bound = {lp.fun:.2f} IP optimum = {ip.fun:.0f}")
def relax_capacity(lam): # lam_i >= 0; each job picks its cheapest machine
mod = cost + lam * load
pick = mod.argmin(axis=1)
used = np.array([load[pick == i, i].sum() for i in range(I)])
return mod.min(axis=1).sum() - lam @ cap, used - cap # L(lam), subgradient
SUBSETS = np.array(list(product([0, 1], repeat=J))) # 64 job subsets
def relax_assignment(mu): # mu_j free; each machine solves a 0-1 knapsack
L, count = mu.sum(), np.zeros(J)
for i in range(I):
ok = SUBSETS[SUBSETS @ load[:, i] <= cap[i]] # subsets that fit machine i
vals = ok @ (cost[:, i] - mu)
L += vals.min(); count += ok[vals.argmin()]
return L, 1 - count # L(mu), subgradient
def subgradient(oracle, dim, sign, ub=ip.fun, iters=300):
"""Maximize L with Polyak steps t = theta (UB - L) / ||g||^2; sign=+1 keeps mult >= 0."""
u, best, theta = np.zeros(dim), -np.inf, 1.0
for k in range(iters):
L, g = oracle(u)
if L > best + 1e-9: best = L
elif k % 20 == 19: theta /= 2 # shrink step if stuck
if not g.any(): break
u = u + theta * (ub - L) / (g @ g) * g
if sign > 0: u = np.maximum(u, 0)
return best
print(f"Lagrangian bound, capacities relaxed = {subgradient(relax_capacity, I, +1):.2f}")
print(f"Lagrangian bound, assignments relaxed = {subgradient(relax_assignment, J, 0):.2f}")
# LP bound = 27.59 IP optimum = 33
# Lagrangian bound, capacities relaxed = 27.59 <- integrality property: same as LP
# Lagrangian bound, assignments relaxed = 31.50 <- knapsack subproblems: beats the LP
The output confirms the theorem: relaxing capacities reproduces the LP bound exactly, while relaxing assignments closes more than half of the remaining gap. For simplicity the code uses the known optimum as $UB$. In practice $UB$ comes from a heuristic, often a Lagrangian heuristic that repairs the relaxed solution $x^\lambda$ into a feasible one.
Key takeaways
- Lagrangian relaxation prices complicating constraints into the objective. Every multiplier vector gives a bound.
- The Lagrangian dual satisfies $z_{LP} \le z_{LD} \le z_{IP}$ (minimization), with equality to the LP bound when the subproblem has the integrality property.
- Subgradient optimization updates multipliers in the direction of constraint violation.
- Choose which constraints to relax so the subproblem is easy but not trivially integral.
Practice
- (Beginner) For $\lambda = (1, 1, 1)$, compute $L(\lambda)$ by hand for the capacity relaxation. Which jobs go to which machine, and which capacities are violated?
- (Intermediate) Write a Lagrangian heuristic for the capacity relaxation: start from the jobs' choices, then move jobs from overloaded machines to the cheapest machine with room. Report the best feasible cost found during the subgradient run.
- (Interview) "When would you use Lagrangian relaxation instead of just calling a MILP solver?" Hint: think about very large models with decomposable structure, the need for fast bounds or heuristics, and problems where the subproblem has a specialized algorithm.
<a id="43-classical-combinatorial-problems"></a>
4.3 Classical Combinatorial Problems
This subchapter presents the problem families that appear again and again inside real applications: knapsack, set covering, facility location, traveling salesman, vehicle routing, bin packing and graph coloring. For each one it gives the standard formulation, a specialized algorithm or heuristic, and the bounds used to judge solution quality.
Knapsack Problem
Why it matters. The knapsack problem is the simplest integer program with a single constraint: pick items with maximum total value without exceeding a capacity. It models budget allocation, cargo loading, portfolio selection and ad placement. It is also a subproblem inside larger methods, such as column generation pricing (Section 2.6) and cover cuts (Section 4.2). It is the textbook example of a problem that is NP-hard yet solvable in practice.
Theory. Items $j = 1, \dots, n$ have weight $w_j$ and value $v_j$. The capacity is $W$.
$$ \max \sum_j v_j x_j \quad \text{s.t.} \quad \sum_j w_j x_j \le W, \quad x_j \in {0, 1}. $$
| Variant | Change | Example |
|---|---|---|
| 0-1 knapsack | each item at most once | project selection |
| Bounded / unbounded | item $j$ up to $u_j$ (or any number of) times | cutting patterns (Section 2.6) |
| Multiple knapsack | several knapsacks, each item in at most one | loading several trucks |
| Multidimensional | several resource constraints | budget and staff limits together |
| Multiple-choice | pick exactly one item from each class | choosing one configuration per server |
Dynamic programming. Let $B_j(c)$ be the best value using items $1..j$ with capacity $c$. Item $j$ is either skipped or taken:
$$ B_j(c) = \max\big{ B_{j-1}(c),\ B_{j-1}(c - w_j) + v_j \big} \quad (\text{second option only if } w_j \le c), \qquad B_0(c) = 0 . $$
The table has $n \times (W + 1)$ cells, so the running time is $O(nW)$. This is polynomial in the value of $W$ but not in its number of digits, a distinction explained in Section 4.4.
Greedy and bounds. Sorting by value density $v_j / w_j$ and adding items while they fit is fast but can be arbitrarily bad. With items (weight 1, value 2) and (weight $W$, value $W$), greedy takes the first item and gets 2 instead of $W$. Returning the better of the greedy solution and the single most valuable item guarantees at least half the optimum. Adding the first item that does not fit fractionally gives the LP relaxation optimum, Dantzig's upper bound.
DP table B_j(c) for 3 items: (w, v) = (2, 3), (3, 4), (4, 5), W = 6
c: 0 1 2 3 4 5 6
j = 0: 0 0 0 0 0 0 0
j = 1: 0 0 3 3 3 3 3 item 1 fits from c = 2
j = 2: 0 0 3 4 4 7 7 c = 5: max(B_1(5) = 3, B_1(2) + 4 = 7)
j = 3: 0 0 3 4 5 7 8 c = 6: max(B_2(6) = 7, B_2(2) + 5 = 8)
answer: B_3(6) = 8 (items 1 and 3)
import numpy as np
# Sprint planning: features with (effort in days, business value). Team capacity: 15 days.
features = {"search": (5, 10), "export": (4, 7), "sso": (6, 12), "dark mode": (2, 3),
"audit log": (3, 6), "webhooks": (7, 13), "i18n": (4, 6)}
names = list(features)
w = [features[k][0] for k in names]; v = [features[k][1] for k in names]
n, W = len(names), 15
# DP table: best[j][c] = max value using the first j items with capacity c
best = np.zeros((n + 1, W + 1), dtype=int)
for j in range(1, n + 1):
for c in range(W + 1):
best[j][c] = best[j - 1][c] # skip item j
if w[j - 1] <= c:
best[j][c] = max(best[j][c], best[j - 1][c - w[j - 1]] + v[j - 1]) # take item j
chosen, c = [], W # walk back through the table
for j in range(n, 0, -1):
if best[j][c] != best[j - 1][c]: # value changed -> item j taken
chosen.append(names[j - 1]); c -= w[j - 1]
print("DP optimum:", best[n][W], "with", sorted(chosen), f"({W - c} days)")
# Greedy by value density, for comparison
room, greedy = W, []
for k in sorted(names, key=lambda k: -features[k][1] / features[k][0]):
if features[k][0] <= room:
greedy.append(k); room -= features[k][0]
print("greedy :", sum(features[k][1] for k in greedy), "with", sorted(greedy))
# Dantzig upper bound: greedy order, last item taken fractionally
room, bound = W, 0.0
for k in sorted(names, key=lambda k: -features[k][1] / features[k][0]):
take = min(1, room / features[k][0]); bound += take * features[k][1]; room -= take * features[k][0]
if room == 0: break
print(f"LP (Dantzig) bound: {bound:.2f}")
# DP optimum: 29 with ['export', 'search', 'sso'] (15 days)
# greedy : 28 with ['audit log', 'search', 'sso']
# LP (Dantzig) bound: 29.86
Greedy fills 14 of 15 days with the three densest features and misses the better combination. The DP checks every capacity level and finds it. The Dantzig bound of 29.86 shows that 29 is optimal, since values are integers.
Key takeaways
- The 0-1 knapsack maximizes value under one capacity constraint. Many variants add bounds, knapsacks or resources.
- DP solves it in $O(nW)$ time by deciding "take or skip" for each item at each capacity.
- Greedy by density is fast but can be arbitrarily bad. Its fractional version gives the LP (Dantzig) bound.
- Knapsack problems appear as subproblems in column generation, cutting planes and Lagrangian relaxation.
Practice
- (Beginner) Extend the ASCII DP table with a fourth item of weight 1 and value 2. What is the new optimum for $W = 6$?
- (Intermediate) Rewrite the DP to use a single one-dimensional array of size $W + 1$. In which direction must the capacity loop run, and why does the other direction solve the unbounded knapsack instead?
- (Interview) "You have 10,000 candidate ads, each with a cost and expected revenue, and a daily budget of 2 million cents. Would you use the DP?" Hint: compute the table size $n \times W$. Then consider the Dantzig bound with greedy, a MILP solver, or scaling the costs (Section 4.4).
Set Covering, Packing and Partitioning
Why it matters. Three closely related models capture a large share of practical integer programs. Covering asks that everything is served at least once (emergency services, sensor placement). Packing asks that nothing is used twice (non-conflicting projects, combinatorial auctions). Partitioning asks for exactly once (airline crew scheduling, where each flight gets exactly one crew). Recognizing them gives you formulations, specialized cuts and approximation guarantees for free.
Theory. Let $A$ be a 0/1 matrix with $a_{ij} = 1$ if set (column) $j$ contains element (row) $i$, and let $x_j = 1$ if set $j$ is chosen.
| Model | Formulation | Meaning | Applications |
|---|---|---|---|
| Set covering | $\min c^\top x$, $Ax \ge \mathbf{1}$ | every element covered at least once | ambulance siting, test-suite minimization, sensor placement |
| Set packing | $\max p^\top x$, $Ax \le \mathbf{1}$ | no element covered twice | combinatorial auctions, conflict-free scheduling |
| Set partitioning | $\min c^\top x$, $Ax = \mathbf{1}$ | every element covered exactly once | crew scheduling, vehicle routing (columns = routes) |
stations (columns)
A B C D E F G
neighborhood 1 1 . . 1 . . 1 covering: each row has >= 1 chosen column
neighborhood 2 1 1 . . . . 1 packing: each row has <= 1 chosen column
neighborhood 3 . . 1 . . . . partitioning: each row has exactly 1
neighborhood 4 . 1 . 1 . 1 .
...
Facts worth knowing:
- Partitioning can be infeasible even when covering is easy. Exact covers are rare.
- Set packing with pairs of conflicting items is the independent set problem on a conflict graph. Clique cuts are its most important valid inequalities.
- Greedy covering (pick the set with the lowest cost per newly covered element) is an $H(m) \approx \ln m$ approximation, where $m$ is the size of the largest set. Under standard assumptions no polynomial algorithm does much better (Section 4.4).
- Set partitioning is the master problem of branch and price for crew scheduling and routing. Each column is a whole pairing or route.
import numpy as np
from scipy.optimize import milp, LinearConstraint, Bounds
# Ambulance siting: candidate stations cover neighborhoods within 8 minutes.
covers = {"A": {1, 2}, "B": {2, 4, 5}, "C": {3, 6, 7}, "D": {1, 4, 8},
"E": {5, 6, 8}, "F": {4, 5, 6, 7, 8}, "G": {1, 2, 7}}
cost = {"A": 4, "B": 3, "C": 4, "D": 3, "E": 3, "F": 4, "G": 4}
S = list(covers); U = sorted(set().union(*covers.values()))
A = np.array([[1 if u in covers[s] else 0 for s in S] for u in U]) # rows: neighborhoods
c = np.array([cost[s] for s in S])
for kind, lb, ub in [("covering (Ax >= 1)", 1, np.inf), ("partitioning (Ax = 1)", 1, 1)]:
res = milp(c, constraints=LinearConstraint(A, lb, ub), integrality=np.ones(len(S)),
bounds=Bounds(0, 1))
pick = [s for s, x in zip(S, res.x) if x > 0.5] if res.success else None
print(f"{kind}: stations {pick}, cost {res.fun if res.success else 'infeasible'}")
# Greedy covering: repeatedly pick the station with the lowest cost per newly covered node
uncovered, greedy = set(U), []
while uncovered:
s = min((s for s in S if covers[s] & uncovered),
key=lambda s: cost[s] / len(covers[s] & uncovered))
greedy.append(s); uncovered -= covers[s]
print(f"greedy covering: stations {greedy}, cost {sum(cost[s] for s in greedy)}")
# Packing: weight = population served, stations must not overlap (no shared neighborhoods)
pop = np.array([5, 4, 6, 4, 5, 9, 4])
res = milp(-pop, constraints=LinearConstraint(A, ub=1), integrality=np.ones(len(S)), bounds=Bounds(0, 1))
print(f"packing (Ax <= 1): stations {[s for s, x in zip(S, res.x) if x > 0.5]}, population {-res.fun:.0f}")
# covering (Ax >= 1): stations ['B', 'C', 'D'], cost 10.0
# partitioning (Ax = 1): stations None, cost infeasible
# greedy covering: stations ['F', 'A', 'C'], cost 12
# packing (Ax <= 1): stations ['A', 'F'], population 14
Greedy is attracted by station F, which covers five neighborhoods cheaply, and then needs two more stations to finish. The optimal cover avoids F entirely. No exact partition exists for this data, which is common: partitioning models usually need many candidate columns to be feasible.
Key takeaways
- Covering ($Ax \ge 1$), packing ($Ax \le 1$) and partitioning ($Ax = 1$) share the same 0/1 matrix and differ only in the constraint sense.
- Partitioning is the most restrictive and may be infeasible. Covering is the most forgiving.
- Greedy covering gives an $O(\ln m)$ approximation. Packing relates to independent sets and clique cuts.
- Crew scheduling and routing are solved as set partitioning over generated columns (branch and price).
Practice
- (Beginner) Verify by hand that stations B, C and D cover all 8 neighborhoods, and that no two stations do.
- (Intermediate) Formulate maximal covering: with a budget for at most 2 stations, maximize the population covered. Add variables $z_i$ (neighborhood $i$ is covered) and solve it with
milp.- (Interview) "How would you choose the smallest set of regression tests that still exercises every line of code changed in a PR?" Hint: identify the elements, the sets and the costs. Then discuss greedy versus exact, depending on suite size.
Facility Location Problems
Why it matters. Where to place warehouses, data centers, ambulances, cell towers or charging stations are long-term decisions with large fixed costs. Facility location models balance the cost of opening facilities against the cost of serving customers from them. They are among the most widely applied MILPs in supply chain design.
Theory. Candidate sites $j$ with opening cost $f_j$ and capacity $u_j$; customers $i$ with demand $d_i$; $c_{ij}$ = cost of serving all of customer $i$'s demand from $j$.
Uncapacitated facility location (UFL):
$$ \min \sum_j f_j y_j + \sum_{i, j} c_{ij} x_{ij} \quad \text{s.t.} \quad \sum_j x_{ij} = 1 \ \forall i, \qquad x_{ij} \le y_j \ \forall i, j, \qquad y \in {0, 1}^n,\ x \ge 0 . $$
With the disaggregated linking constraints (Section 4.1), the LP relaxation is often integral or nearly so.
Capacitated facility location (CFL) adds $\sum_i d_i x_{ij} \le u_j y_j$. Two versions exist:
- split deliveries: $x_{ij} \in [0, 1]$, so a customer may be served by several sites;
- single sourcing: $x_{ij} \in {0, 1}$, so each customer gets exactly one site. This is common in practice (one supplier per store) and much harder.
| Model | Objective | Typical use |
|---|---|---|
| UFL / CFL | fixed + service cost | warehouses, plants |
| p-median | min total (demand-weighted) distance with exactly $p$ facilities | public services with a fixed budget |
| p-center | min the maximum distance with $p$ facilities | emergency services (worst case matters) |
| Covering location | min facilities so that everyone is within distance $D$ | fire stations (set covering, previous topic) |
| Maximal covering | max population within distance $D$ with $p$ facilities | budget-limited coverage |
The code reuses the 4 depots and 8 customers from "Strength of Formulations", adding demands and capacities:
import numpy as np
from scipy.optimize import milp, LinearConstraint, Bounds
# Same depots and customers as in "Strength of Formulations", now with demands and capacities.
f = np.array([40, 30, 35, 38])
C = np.array([[9, 22, 18, 30], [12, 15, 25, 24], [20, 8, 30, 19], [25, 10, 14, 27],
[30, 20, 7, 16], [28, 26, 11, 9], [15, 30, 21, 12], [18, 24, 26, 6]])
d = np.array([3, 2, 4, 3, 2, 3, 2, 4]) # demand of each customer
m, n = C.shape
def solve(u=None, single_source=True, p=None):
"""u = capacities (None = uncapacitated), p = exact number of depots (p-median style)."""
nv = n + m * n
c = np.concatenate([f if p is None else np.zeros(n), C.ravel()])
A, lb, ub = [], [], []
for i in range(m): # serve every customer
r = np.zeros(nv); r[n + i * n: n + (i + 1) * n] = 1; A.append(r); lb.append(1); ub.append(1)
for i in range(m): # disaggregated linking x_ij <= y_j
for j in range(n):
r = np.zeros(nv); r[n + i * n + j] = 1; r[j] = -1; A.append(r); lb.append(-np.inf); ub.append(0)
if u is not None: # sum_i d_i x_ij <= u_j y_j
for j in range(n):
r = np.zeros(nv); r[[n + i * n + j for i in range(m)]] = d; r[j] = -u[j]
A.append(r); lb.append(-np.inf); ub.append(0)
if p is not None: # open exactly p depots
r = np.zeros(nv); r[:n] = 1; A.append(r); lb.append(p); ub.append(p)
integ = np.r_[np.ones(n), np.ones(m * n) if single_source else np.zeros(m * n)]
res = milp(c, constraints=LinearConstraint(A, lb, ub), integrality=integ, bounds=Bounds(0, 1))
y = res.x[:n] > 0.5
x = res.x[n:].reshape(m, n)
loads = (d[:, None] * x).sum(axis=0).round(1) + 0.0
return res.fun, [j + 1 for j in range(n) if y[j]], loads
print("UFL : cost %.1f, open %s, loads %s" % solve())
print("CFL, split deliveries : cost %.1f, open %s, loads %s" % solve(u=np.array([12, 13, 8, 10]), single_source=False))
print("CFL, single sourcing : cost %.1f, open %s, loads %s" % solve(u=np.array([12, 13, 8, 10])))
print("p-median (p = 3) : cost %.1f, open %s, loads %s" % solve(p=3))
# UFL : cost 166.0, open [2, 4], loads [ 0. 12. 0. 11.]
# CFL, split deliveries : cost 168.0, open [2, 4], loads [ 0. 13. 0. 10.]
# CFL, single sourcing : cost 188.0, open [2, 3, 4], loads [0. 9. 5. 9.]
# p-median (p = 3) : cost 82.0, open [1, 2, 4], loads [ 5. 7. 0. 11.]
The uncapacitated plan loads depot 4 with 11 units, one more than its capacity of 10. With split deliveries, one unit is shifted to depot 2 at a cost of 2. With single sourcing, customers cannot be split, and the cheapest feasible plan opens a third depot. The cost jumps by 22. The p-median line counts service cost only (no opening costs), with exactly 3 depots.
Key takeaways
- UFL trades opening costs against service costs. Its disaggregated formulation has a very strong LP relaxation.
- Capacities make the problem harder. Single sourcing makes it harder still, because it adds a bin-packing-like structure.
- p-median fixes the number of facilities and minimizes total distance. p-center minimizes the worst distance.
- Covering-type location models are set covering problems with distance-based coverage sets.
Practice
- (Beginner) In the UFL solution, list which customers are served by depot 2 and which by depot 4, using the cost matrix (each customer goes to its cheapest open depot).
- (Intermediate) Formulate the p-center problem for the same data with $p = 2$: minimize the largest $c_{ij}$ used. Add a variable $t$ with $c_{ij} x_{ij} \le t$ and solve it.
- (Interview) "A retailer wants to choose 5 new distribution centers out of 60 candidates for 2,000 stores. How would you approach it?" Hint: p-median or CFL with single sourcing. Discuss model size, aggregating stores into zones, a strong formulation, and Lagrangian or heuristic methods if the MILP is too large.
Traveling Salesman Problem
Why it matters. The traveling salesman problem (TSP) asks for the shortest tour visiting every city exactly once. It models delivery routes, drilling holes in circuit boards, sequencing jobs with setup times, and DNA fragment assembly. It is also the testing ground where most ideas of integer programming (cutting planes, branch and cut, heuristics) were developed. The Concorde solver has proven optimality for instances with 85,900 cities.
Theory. For a symmetric TSP on a complete graph with edge costs $c_e$, let $x_e = 1$ if edge $e$ is in the tour. The Dantzig-Fulkerson-Johnson (DFJ) formulation is
$$ \min \sum_e c_e x_e \quad \text{s.t.} \quad \sum_{e \in \delta(v)} x_e = 2 \ \ \forall v, \qquad \sum_{e \in E(S)} x_e \le \lvert S \rvert - 1 \ \ \forall S \subset V,\ 2 \le \lvert S \rvert \le n - 1, \qquad x_e \in {0, 1}. $$
- Degree constraints: every city has exactly two tour edges.
- Subtour elimination: a subset $S$ cannot contain a closed cycle of its own. Without these, the solution can be several disjoint loops.
There are about $2^n$ subtour constraints, so they are added lazily: solve with the degree constraints only, find the connected components of the solution, add the violated subtour constraints for them, and repeat. This is the TSP's branch-and-cut idea in its simplest form.
Iteration 1: degree constraints only -> three subtours (length 27.1)
8 9 Cluster C: 6-7-8-9 (a closed loop of its own)
6 7
add: sum of edges inside {6,7,8,9} <= 3
1 5 Cluster B: 3-4-5
2 3 4 add: sum of edges inside {3,4,5} <= 2
0
Cluster A: 0-1-2, same cut
Iteration 2: one tour through all 10 cities (length 39.5) -> optimal
MTZ formulation. The Miller-Tucker-Zemlin formulation for the directed TSP uses order variables $u_i$ and the constraints $u_i - u_j + n x_{ij} \le n - 1$ for $i, j \ne 1$. It has only polynomially many constraints, but its LP relaxation is much weaker than DFJ, so it only suits small instances.
| Method | Type | Notes |
|---|---|---|
| DFJ + branch and cut (Concorde) | exact | subtour, comb and other cuts. State of the art |
| Held-Karp DP | exact | $O(n^2 2^n)$, up to about 25 cities (Section 4.4) |
| Nearest neighbor, insertion | construction heuristic | fast, often 15-25% above optimum (Chapter 11) |
| 2-opt, Lin-Kernighan | local search | LKH is often within 1-2% of optimal (Chapter 11) |
| Double tree, Christofides | approximation | 2 and 1.5 times optimum for metric costs (Section 4.4) |
import numpy as np
from itertools import combinations
from scipy.optimize import milp, LinearConstraint, Bounds
# 10 cities in three clusters; symmetric Euclidean distances, rounded
pts = np.array([[0, 0], [1, 3], [3, 1], [10, 0], [12, 2], [11, 4], [5, 9], [7, 10], [4, 11], [6, 12]])
n = len(pts)
E = list(combinations(range(n), 2)) # one binary x_e per edge {i, j}
dist = np.array([round(np.hypot(*(pts[i] - pts[j])), 1) for i, j in E])
deg = np.array([[1 if k in e else 0 for e in E] for k in range(n)]) # degree rows
cuts = [] # subtour elimination rows
def components(x):
"""Connected components of the graph formed by edges with x_e = 1."""
adj = {k: set() for k in range(n)}
for (i, j), v in zip(E, x):
if v > 0.5: adj[i].add(j); adj[j].add(i)
seen, comps = set(), []
for s in range(n):
if s in seen: continue
stack, comp = [s], set()
while stack:
u = stack.pop()
if u not in comp: comp.add(u); stack.extend(adj[u] - comp)
seen |= comp; comps.append(sorted(comp))
return comps
for it in range(1, 20):
cons = [LinearConstraint(deg, 2, 2)] # every city has degree 2
if cuts:
cons.append(LinearConstraint([r for r, _ in cuts], ub=[b for _, b in cuts]))
res = milp(dist, constraints=cons, integrality=np.ones(len(E)), bounds=Bounds(0, 1))
comps = components(res.x)
print(f"iteration {it}: length {res.fun:5.1f}, components {comps}")
if len(comps) == 1:
break
for S in comps: # sum_{e inside S} x_e <= |S| - 1
cuts.append((np.array([1 if i in S and j in S else 0 for i, j in E]), len(S) - 1))
tour, prev = [0], None # walk the tour from city 0
while len(tour) < n:
nxt = [j if i == tour[-1] else i for (i, j), v in zip(E, res.x)
if v > 0.5 and tour[-1] in (i, j) and (j if i == tour[-1] else i) not in tour]
tour.append(nxt[0])
print("optimal tour:", " -> ".join(map(str, tour + [0])), f"length {res.fun:.1f}")
# iteration 1: length 27.1, components [[0, 1, 2], [3, 4, 5], [6, 7, 8, 9]]
# iteration 2: length 39.5, components [[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]]
# optimal tour: 0 -> 1 -> 6 -> 8 -> 9 -> 7 -> 5 -> 4 -> 3 -> 2 -> 0 length 39.5
Each iteration here solves an integer program, and the loop stops at the first connected solution. Because it only ever removes subtours, that solution is optimal. Real branch-and-cut codes separate subtour cuts on fractional LP solutions too, using minimum cut computations (Section 3.4).
Key takeaways
- The DFJ formulation uses degree-2 constraints plus exponentially many subtour elimination constraints.
- Subtour constraints are added lazily: solve, find disconnected components, cut them off, repeat.
- MTZ is compact but weak. DFJ with branch and cut is the exact method of choice.
- Large TSPs are solved by branch and cut (exactly) or by Lin-Kernighan-style local search (near-optimally).
Practice
- (Beginner) How many distinct tours does a symmetric TSP with 10 cities have? How many subtour elimination constraints does DFJ have for $n = 10$?
- (Intermediate) Implement the MTZ formulation for the same 10 cities (directed arcs, $u_i$ variables) and compare its LP relaxation value with the DFJ LP relaxation obtained after adding the two rounds of subtour cuts.
- (Interview) "A warehouse robot must pick 40 items per order, and orders arrive every few seconds. Would you solve each picking route to optimality?" Hint: compare the time budget with the solution quality of a fast heuristic plus 2-opt. Consider the warehouse's aisle structure, which allows special-purpose algorithms.
Vehicle Routing Problem
Why it matters. The vehicle routing problem (VRP) generalizes the TSP to a fleet: several vehicles leave a depot, serve customers with known demands, and return, without exceeding vehicle capacity. It is the core optimization problem of parcel delivery, grocery distribution, waste collection and field service. Small percentage savings translate into large fuel and labor costs.
Theory. The capacitated VRP (CVRP) has a depot 0, customers $1..n$ with demands $q_i$, identical vehicles of capacity $Q$, and travel costs $c_{ij}$. Find routes, each starting and ending at the depot and carrying at most $Q$, that visit every customer once at minimum total cost.
| Variant | Extra feature |
|---|---|
| CVRP | capacity only |
| VRPTW | customer time windows $[a_i, b_i]$ |
| Pickup and delivery | paired pickup and drop-off locations (ride sharing, couriers) |
| Multi-depot / heterogeneous fleet | several depots, vehicles of different sizes and costs |
| Periodic VRP | visits spread over several days |
Formulations.
- Two-index flow formulation: arc variables $x_{ij}$, degree constraints, and capacity cuts $\sum_{(i,j) \in \delta^+(S)} x_{ij} \ge \lceil q(S) / Q \rceil$. Any customer set $S$ needs at least that many vehicles entering and leaving it. These also eliminate subtours.
- MTZ-style load variables: $u_i$ is the load on the vehicle after visiting $i$, with $q_i \le u_i \le Q$ and $u_j \ge u_i + q_j - Q(1 - x_{ij})$. Compact but weak. Fine for small instances.
- Set partitioning over routes: one column per feasible route, solved by branch and price with a resource-constrained shortest path as pricing. This is the basis of the best exact algorithms.
Clarke-Wright savings heuristic. Start with one route per customer ($0 \to i \to 0$). Merging the routes of $i$ and $j$ via the edge $(i, j)$ saves
$$ s_{ij} = c_{0i} + c_{0j} - c_{ij}. $$
Process pairs in decreasing order of savings and merge whenever $i$ and $j$ are at the ends of different routes and the merged load fits.
before: 0 -> i -> 0 and 0 -> j -> 0 cost c0i + ci0 + c0j + cj0
after: 0 -> i -> j -> 0 cost c0i + cij + cj0
saving: c0i + c0j - cij
import numpy as np
from scipy.optimize import milp, LinearConstraint, Bounds
# Depot 0 and 7 customers; truck capacity Q = 10
pts = np.array([[0, 0], [-4, 3], [-2, 6], [3, 5], [5, 1], [4, -4], [-3, -5], [-6, -1]])
q = np.array([0, 4, 3, 4, 3, 5, 3, 4]) # demands (depot has none)
Q, n = 10, len(pts)
D = np.round(np.hypot(*(pts[:, None] - pts[None]).transpose(2, 0, 1)), 2)
def route_len(r):
stops = [0] + r + [0]
return sum(D[a, b] for a, b in zip(stops, stops[1:]))
# Clarke-Wright savings: start with one route per customer, merge by largest saving
routes = {i: [i] for i in range(1, n)} # route id -> customer list
where = {i: i for i in range(1, n)} # customer -> route id
savings = sorted(((D[0, i] + D[0, j] - D[i, j], i, j)
for i in range(1, n) for j in range(i + 1, n)), reverse=True)
for s, i, j in savings:
ri, rj = where[i], where[j]
if ri == rj or sum(q[routes[ri]]) + sum(q[routes[rj]]) > Q:
continue
a, b = routes[ri], routes[rj]
# merge only if i and j are route ends, so the depot -> ... -> depot shape is kept
if a[-1] == i and b[0] == j: merged = a + b
elif a[0] == i and b[-1] == j: merged = b + a
elif a[-1] == i and b[-1] == j: merged = a + b[::-1]
elif a[0] == i and b[0] == j: merged = a[::-1] + b
else: continue
routes[ri] = merged; del routes[rj]
for k in merged: where[k] = ri
cw = list(routes.values())
print("Clarke-Wright:", cw, f"total {sum(map(route_len, cw)):.2f}")
# Exact: directed arcs x_ij, load variables u_i (MTZ-style capacity/subtour elimination)
arcs = [(i, j) for i in range(n) for j in range(n) if i != j]
na, nv = len(arcs), len(arcs) + n
A, lb, ub = [], [], []
for k in range(1, n): # one arc in, one arc out per customer
for side in (0, 1):
r = np.zeros(nv); r[[a for a, e in enumerate(arcs) if e[side] == k]] = 1
A.append(r); lb.append(1); ub.append(1)
for i, j in arcs: # u_j >= u_i + q_j - Q (1 - x_ij)
if i and j:
r = np.zeros(nv); r[na + i] = 1; r[na + j] = -1; r[arcs.index((i, j))] = Q
A.append(r); lb.append(-np.inf); ub.append(Q - q[j])
c = np.r_[[D[i, j] for i, j in arcs], np.zeros(n)]
res = milp(c, constraints=LinearConstraint(A, lb, ub),
integrality=np.r_[np.ones(na), np.zeros(n)],
bounds=Bounds(np.r_[np.zeros(na), q], np.r_[np.ones(na), [0] + [Q] * (n - 1)]))
succ = {i: j for (i, j), v in zip(arcs, res.x) if v > 0.5}
exact = []
for start in [j for (i, j), v in zip(arcs, res.x) if v > 0.5 and i == 0]:
r, k = [], start
while k != 0: r.append(k); k = succ[k]
exact.append(r)
print("exact MILP :", exact, f"total {res.fun:.2f}")
print("loads :", [int(sum(q[r])) for r in exact], "capacity", Q)
# Clarke-Wright: [[1, 2], [3, 4], [5], [6, 7]] total 58.56
# exact MILP : [[1, 7], [4, 3, 2], [6, 5]] total 55.10
# loads : [8, 10, 8] capacity 10
The savings heuristic locks in the merges $1$-$2$, $6$-$7$ and $3$-$4$ (largest savings first). Each of these routes then carries 7 units, so customer 5 (demand 5) fits nowhere and gets a route of its own. The exact solution needs one truck fewer and is about 6% shorter. Savings plus local search (Chapter 11) usually closes most of that gap at a fraction of the cost.
Key takeaways
- The CVRP routes a fleet from a depot under capacity limits. Time windows, pickups and multiple depots are common extensions.
- Capacity cuts $\sum_{\delta^+(S)} x \ge \lceil q(S)/Q \rceil$ eliminate subtours and capacity violations in the two-index formulation.
- The best exact methods use set partitioning over routes with branch and price.
- Clarke-Wright savings is the classic fast construction heuristic, usually followed by local search.
Practice
- (Beginner) Compute the savings $s_{12}$ and $s_{23}$ from the coordinates. Why is $s_{23}$, the second-largest saving, rejected by the algorithm?
- (Intermediate) Add a time limit: each route may be at most 25 distance units long. Add this rule to the Clarke-Wright merge check and compare the result.
- (Interview) "A grocery delivery company has 300 orders per evening, 2-hour delivery windows and 25 vans. Outline your solution approach." Hint: VRPTW. Discuss a construction heuristic, local search or a metaheuristic (Chapter 11), runtime limits, and how you would measure solution quality without an exact optimum.
Bin Packing Problem
Why it matters. Bin packing asks how to put items of given sizes into the fewest bins of fixed capacity. It models VM placement on servers (Section 1.1), loading trucks or containers, cutting material, and scheduling jobs on identical machines with a deadline. Its simple heuristics have proven performance guarantees, which makes it a reference problem for approximation algorithms.
Theory. Items $i = 1..n$ with sizes $s_i$, bins of capacity $C$. With $y_j = 1$ if bin $j$ is used and $x_{ij} = 1$ if item $i$ goes into bin $j$ (for $j = 1..m$, where $m$ is any upper bound on the number of bins):
$$ \min \sum_j y_j \quad \text{s.t.} \quad \sum_j x_{ij} = 1 \ \forall i, \qquad \sum_i s_i x_{ij} \le C, y_j \ \forall j, \qquad x, y \in {0, 1}. $$
Weaknesses of this formulation.
- Symmetry: bins are identical, so every solution appears $m!$ times with bins relabeled. Adding $y_1 \ge y_2 \ge \dots \ge y_m$ removes part of it.
- Weak LP bound: the LP relaxation value equals $\sum_i s_i / C$, the trivial bound.
The pattern formulation, with one variable per feasible bin filling, is the cutting stock model of Section 2.6. It has a much stronger bound and is solved by column generation and branch and price.
Bounds and heuristics.
| Method | Guarantee |
|---|---|
| $L_1 = \lceil \sum_i s_i / C \rceil$ | lower bound (volume) |
| Next Fit | $\le 2 \cdot \text{OPT}$ |
| First Fit | $\le 1.7 \cdot \text{OPT}$ (asymptotically, and tight) |
| First Fit Decreasing (FFD) | $\le \tfrac{11}{9}\text{OPT} + \tfrac{6}{9}$ |
| Any polynomial algorithm | cannot guarantee better than $\tfrac{3}{2}\text{OPT}$ unless P = NP (Section 4.4) |
import math
import numpy as np
from scipy.optimize import milp, LinearConstraint, Bounds
# VM memory requests (GB) to place on hosts with 100 GB each
s = [54, 53, 36, 34, 24, 24, 24, 23, 20]
C, n = 100, len(s)
def first_fit_decreasing(sizes, cap):
bins = []
for x in sorted(sizes, reverse=True):
for b in bins:
if sum(b) + x <= cap: b.append(x); break
else: bins.append([x])
return bins
ffd = first_fit_decreasing(s, C)
print("FFD :", len(ffd), "bins", ffd)
print("L1 lower bound = ceil(sum / C) =", math.ceil(sum(s) / C))
# Assignment formulation: y_j = bin j used, x_ij = item i in bin j; m = FFD bins as upper bound
m = len(ffd); nv = m + n * m
x = lambda i, j: m + i * m + j
A, lb, ub = [], [], []
for i in range(n): # each item in exactly one bin
r = np.zeros(nv); r[[x(i, j) for j in range(m)]] = 1; A.append(r); lb.append(1); ub.append(1)
for j in range(m): # capacity, only if bin is used
r = np.zeros(nv); r[[x(i, j) for i in range(n)]] = s; r[j] = -C; A.append(r); lb.append(-np.inf); ub.append(0)
for j in range(m - 1): # symmetry breaking: y_j >= y_{j+1}
r = np.zeros(nv); r[j], r[j + 1] = -1, 1; A.append(r); lb.append(-np.inf); ub.append(0)
res = milp(np.r_[np.ones(m), np.zeros(n * m)], constraints=LinearConstraint(A, lb, ub),
integrality=np.ones(nv), bounds=Bounds(0, 1))
bins = [[s[i] for i in range(n) if res.x[x(i, j)] > 0.5] for j in range(m)]
print("exact :", round(res.fun), "bins", [b for b in bins if b])
# FFD : 4 bins [[54, 36], [53, 34], [24, 24, 24, 23], [20]]
# L1 lower bound = ceil(sum / C) = 3
# exact : 3 bins [[36, 34, 24], [54, 24, 20], [53, 24, 23]]
FFD pairs each large item with the next-largest that fits (54 with 36, 53 with 34), leaving 10 and 13 GB unusable. The optimum spreads the 24 GB items so that every bin is almost full. Unlike the VM example of Section 1.1, here the volume bound is reached only by the exact method.
Key takeaways
- Bin packing minimizes the number of fixed-capacity bins. It models VM placement, loading and cutting.
- The assignment formulation is symmetric and has a weak LP bound. The pattern (cutting stock) formulation is much stronger.
- FFD is simple and guaranteed within about 22% of optimal. The volume bound $L_1$ certifies quality.
- No polynomial algorithm can guarantee better than $3/2$ times optimal unless P = NP.
Practice
- (Beginner) Run Next Fit (only the last bin stays open) by hand on the code's instance. How many bins does it use?
- (Intermediate) Implement Best Fit Decreasing and compare FFD, BFD and the $L_1$ bound on 100 random instances with sizes uniform in $[20, 60]$ and $C = 100$. How often does each heuristic reach $L_1$?
- (Interview) "Your Kubernetes scheduler places pods using first fit. Utilization is 62%. How would you improve it?" Hint: consider sorting (decreasing order), multiple resources (CPU and memory make it a vector bin packing problem), periodic re-packing with an exact or column-generation model, and the cost of migrations.
Graph Coloring
Why it matters. Coloring a graph assigns a color to every node so that adjacent nodes get different colors, using as few colors as possible. Colors stand for time slots, frequencies, registers or machines, and edges stand for conflicts. Exam timetabling, compiler register allocation, radio frequency assignment and scheduling with incompatibilities are all coloring problems.
Theory. The minimum number of colors is the chromatic number $\chi(G)$. A set of pairwise adjacent nodes (a clique) needs as many colors as it has nodes, so the clique number $\omega(G)$ is a lower bound:
$$ \omega(G) \le \chi(G) \le \Delta(G) + 1, $$
where $\Delta(G)$ is the maximum degree. The upper bound is achieved by any greedy coloring.
ILP formulation with $K$ available colors, $x_{vk} = 1$ if node $v$ gets color $k$, $w_k = 1$ if color $k$ is used:
$$ \min \sum_k w_k \quad \text{s.t.} \quad \sum_k x_{vk} = 1 \ \forall v, \qquad x_{uk} + x_{vk} \le w_k \ \forall {u, v} \in E,\ \forall k . $$
Colors are interchangeable, so this formulation is highly symmetric. Symmetry breaking ($w_k \ge w_{k+1}$) and alternative formulations (set partitioning over independent sets, solved by branch and price) are standard improvements.
Heuristics.
- Greedy: color nodes in some order, each with the smallest color not used by its neighbors. The result depends heavily on the order.
- DSatur (Brélaz, 1979): always color next the uncolored node with the most distinct colors among its neighbors (highest saturation), breaking ties by degree. It colors the hardest nodes first and is exact on bipartite graphs.
| Application | Nodes | Edge means | Colors are |
|---|---|---|---|
| Exam timetabling | exams | a student takes both | time slots |
| Register allocation | program variables | live at the same time | CPU registers |
| Frequency assignment | transmitters | would interfere | frequencies |
| Map coloring | regions | share a border | colors (4 suffice for planar maps) |
import numpy as np
from itertools import combinations
from scipy.optimize import milp, LinearConstraint, Bounds
# Exam timetabling: an edge = two exams share a student, so they need different time slots
exams = ["Algebra", "Biology", "Chemistry", "Databases", "Economics",
"French", "Geology", "History", "Italian"]
edges = [(0, 3), (0, 5), (1, 3), (1, 6), (1, 8), (2, 4), (2, 6), (2, 7),
(3, 4), (4, 6), (4, 7), (5, 7)]
n = len(exams)
adj = {v: set() for v in range(n)}
for a, b in edges: adj[a].add(b); adj[b].add(a)
def first_free(v, col):
used = {col[u] for u in adj[v] if u in col}
return next(c for c in range(n) if c not in used)
col = {} # greedy in alphabetical order
for v in range(n): col[v] = first_free(v, col)
print("greedy (alphabetical):", max(col.values()) + 1, "slots")
col = {} # DSatur: most-constrained exam first
while len(col) < n:
v = max((v for v in range(n) if v not in col),
key=lambda v: (len({col[u] for u in adj[v] if u in col}), len(adj[v])))
col[v] = first_free(v, col)
print("DSatur :", max(col.values()) + 1, "slots",
{exams[v]: col[v] for v in range(n)})
# Lower bound: largest clique (exams that pairwise conflict)
clique = max((S for k in range(1, n + 1) for S in combinations(range(n), k)
if all(b in adj[a] for a, b in combinations(S, 2))), key=len)
print("clique lower bound :", len(clique), [exams[v] for v in clique])
# Exact ILP: x[v, k] = exam v in slot k, w[k] = slot k used; K = 4 slots available
K = 4; nv = K + n * K
x = lambda v, k: K + v * K + k
A, lb, ub = [], [], []
for v in range(n): # each exam gets one slot
r = np.zeros(nv); r[[x(v, k) for k in range(K)]] = 1; A.append(r); lb.append(1); ub.append(1)
for a, b in edges: # x[a,k] + x[b,k] <= w[k]
for k in range(K):
r = np.zeros(nv); r[[x(a, k), x(b, k)]] = 1; r[k] = -1; A.append(r); lb.append(-np.inf); ub.append(0)
for k in range(K - 1): # symmetry breaking: w[k] >= w[k+1]
r = np.zeros(nv); r[k], r[k + 1] = -1, 1; A.append(r); lb.append(-np.inf); ub.append(0)
res = milp(np.r_[np.ones(K), np.zeros(n * K)], constraints=LinearConstraint(A, lb, ub),
integrality=np.ones(nv), bounds=Bounds(0, 1))
print("exact ILP :", round(res.fun), "slots")
# greedy (alphabetical): 4 slots
# DSatur : 3 slots {'Algebra': 0, 'Biology': 0, 'Chemistry': 1, 'Databases': 1, 'Economics': 0, 'French': 1, 'Geology': 2, 'History': 2, 'Italian': 1}
# clique lower bound : 3 ['Chemistry', 'Economics', 'Geology']
# exact ILP : 3 slots
Greedy in alphabetical order wastes a fourth slot. DSatur finds 3 slots, and the clique {Chemistry, Economics, Geology} proves that 3 is optimal without running the ILP. On large graphs the clique bound can be far below $\chi(G)$, and that is where exact methods are needed.
Key takeaways
- Graph coloring assigns conflict-free colors (slots, registers, frequencies) using as few as possible.
- $\omega(G) \le \chi(G) \le \Delta(G) + 1$: cliques give lower bounds, greedy gives upper bounds.
- The assignment ILP is symmetric. Break symmetry or use independent-set columns (branch and price).
- DSatur colors the most constrained nodes first and usually beats plain greedy.
Practice
- (Beginner) Find an order of the exams for which plain greedy uses only 3 slots. Why does such an order always exist?
- (Intermediate) Each slot can hold at most 3 exams (room limit). Add this constraint to the ILP. Does the optimum change?
- (Interview) "How does a compiler decide which variables share a CPU register, and what happens when the graph needs more colors than there are registers?" Hint: build the interference graph from liveness. When coloring fails, the compiler spills variables to memory, choosing spill candidates by cost.
<a id="44-computational-complexity-in-or"></a>
4.4 Computational Complexity in OR
This subchapter explains why some problems in this book have fast algorithms (LP, shortest paths, flows) while others seem to resist them (TSP, bin packing, integer programming). It also covers what to do about the hard ones: exploit pseudo-polynomial structure, or settle for solutions with a proven quality guarantee.
P, NP and NP-Completeness
Why it matters. Before spending weeks looking for a fast exact algorithm, you want to know whether one is likely to exist. Complexity theory answers this. If your problem is NP-hard, the right goal changes: good formulations and branch and cut for moderate sizes, approximation algorithms or heuristics for large ones. It also explains why "the solver is slow" is sometimes a property of the problem, not of the code.
Theory. Complexity theory classifies decision problems (yes/no questions). An optimization problem has a decision version: "is there a tour of length $\le K$?" instead of "find the shortest tour". Running time is measured in the input size: the number of bits needed to write the instance.
| Class | Definition | Intuition |
|---|---|---|
| P | solvable in time polynomial in the input size | efficiently solvable |
| NP | a "yes" answer has a certificate checkable in polynomial time | efficiently verifiable |
| NP-hard | every problem in NP reduces to it in polynomial time | at least as hard as anything in NP |
| NP-complete | in NP and NP-hard | the hardest problems in NP |
NP stands for nondeterministic polynomial, not "not polynomial". Every problem in P is also in NP, because a problem you can solve quickly you can also verify quickly. Whether P = NP is the most famous open problem in computer science. Almost everyone believes P $\neq$ NP. If so, no NP-complete problem has a polynomial-time algorithm.
The first NP-complete problem was SAT (Cook and Levin, 1971). Karp (1972) then showed that 21 classic combinatorial problems are NP-complete, including several that OR had been struggling with for decades.
+-------------------------------------------+
| NP-hard |
| (optimization TSP, general MILP, ...) |
+-------------------+---------------------+ |
| NP | NP-complete | |
| | SAT, decision TSP, | |
| | knapsack, coloring,| |
| | bin packing | |
| +-------------+ +---------------------+---------------------+
| | P | |
| | LP, MST, | |
| | max flow, | |
| | matching | |
| +-------------+ |
+-----------------------------------------+
(drawn assuming P != NP)
| Polynomial (P) | NP-hard |
|---|---|
| Linear programming (ellipsoid, interior point; Section 2.6) | Integer linear programming |
| Shortest paths, minimum spanning trees (Chapter 3) | Traveling salesman, vehicle routing |
| Maximum flow, minimum cost flow | Fixed-charge network design, Steiner tree |
| Assignment, bipartite and general matching | Generalized assignment, 3-dimensional matching |
| 2-SAT | SAT, 3-SAT |
| Single-machine scheduling: SPT, EDD rules (Chapter 10) | Job shop scheduling, weighted tardiness |
| Fractional knapsack | 0-1 knapsack, bin packing, set cover, graph coloring |
What NP-hard does and does not mean. It is a statement about the worst case as instances grow. Many NP-hard problems are solved routinely at practical sizes: knapsacks with thousands of items, TSPs with thousands of cities, MILPs with millions of variables. Conversely, a problem in P can still be too slow if the polynomial has a high degree.
Certificates are asymmetric. For "is there a tour of length $\le K$?", a tour is a short proof of "yes". For "no", there is in general no short proof known: problems whose "no" answers have short certificates form the class co-NP. In practice, solvers prove "no better solution exists" with bounds, using LP relaxations and branch-and-bound trees, which can be exponentially large.
import math
# Decision version of TSP: "is there a tour of length <= K?" (cities from the TSP example)
pts = [(0, 0), (1, 3), (3, 1), (10, 0), (12, 2), (11, 4), (5, 9), (7, 10), (4, 11), (6, 12)]
d = lambda a, b: math.dist(pts[a], pts[b])
def verify(tour, K):
"""Polynomial-time verifier: O(n) checks on a proposed certificate (a tour)."""
n = len(pts)
if sorted(tour) != list(range(n)): # visits every city exactly once
return False
length = sum(d(tour[k], tour[(k + 1) % n]) for k in range(n))
return length <= K
certificate = [0, 1, 6, 8, 9, 7, 5, 4, 3, 2] # the optimal tour found earlier
print("certificate valid for K = 40?", verify(certificate, 40))
print("certificate valid for K = 39?", verify(certificate, 39)) # says nothing about other tours
# Finding a tour by brute force means checking (n-1)!/2 tours
print(f"{'n':>4} {'n^3':>10} {'2^n':>22} {'(n-1)!/2 tours':>28}")
for n in (10, 20, 30, 50):
print(f"{n:>4} {n**3:>10,} {2**n:>22,} {math.factorial(n - 1) // 2:>28.3e}")
# certificate valid for K = 40? True
# certificate valid for K = 39? False
# n n^3 2^n (n-1)!/2 tours
# 10 1,000 1,024 1.814e+05
# 20 8,000 1,048,576 6.082e+16
# 30 27,000 1,073,741,824 4.421e+30
# 50 125,000 1,125,899,906,842,624 3.041e+62
Checking a certificate takes a handful of operations. Finding one by enumeration does not scale beyond tiny $n$. For $K = 39$, the verifier rejects this tour, but that alone does not prove that no shorter tour exists. That proof required the full branch-and-cut run of Section 4.3.
Key takeaways
- P: solvable in polynomial time. NP: "yes" answers verifiable in polynomial time. NP-complete: the hardest problems in NP.
- If P $\neq$ NP (widely believed), NP-hard problems have no polynomial-time algorithms.
- Most problems of Chapter 3 are in P. Integer programming and the combinatorial problems of Section 4.3 are NP-hard.
- NP-hardness is a worst-case statement. Practical instances are often solved well by branch and cut or heuristics.
Practice
- (Beginner) For each problem, state the decision version and a certificate: (a) graph coloring with $k$ colors, (b) set cover with cost $\le K$, (c) bin packing with $\le B$ bins.
- (Intermediate) The TSP is NP-hard but the assignment problem is in P, although both choose $n$ arcs with one in and one out per node. Explain the difference using the assignment relaxation of the TSP and its subtours.
- (Interview) "A product manager asks you for an algorithm that always finds the optimal delivery routes for 5,000 stops in under a second. How do you respond?" Hint: explain NP-hardness without jargon, then propose what is achievable: near-optimal heuristics with quality bounds, and exact solutions for smaller subproblems.
Polynomial Reductions
Why it matters. Reductions are how hardness is proven and how solution methods are reused. To show that a new problem is hard, you transform a known hard problem into it. To solve a new problem, you transform it into one you already know how to solve, as bipartite matching was turned into max flow in Section 3.4.
Theory. Problem $A$ reduces to problem $B$ ($A \le_p B$) if there is a polynomial-time transformation $T$ that maps every instance $I$ of $A$ to an instance $T(I)$ of $B$ with the same answer:
$$ I \text{ is a yes-instance of } A \iff T(I) \text{ is a yes-instance of } B . $$
Then an algorithm for $B$ also solves $A$: transform, then solve. Read $A \le_p B$ as "$B$ is at least as hard as $A$".
Direction matters. To prove that $B$ is NP-hard, reduce a known NP-hard problem $A$ to $B$. Reducing $B$ to a hard problem proves nothing about $B$: a problem can always be made harder by encoding it badly.
A classic tree of reductions (arrow A -> B means A <=_p B)
SAT
|
3-SAT
/ | \
Independent Set Hamiltonian Subset Sum
/ \ Cycle / \
Vertex Clique | Partition Knapsack
Cover TSP |
| Bin Packing
Set Cover
Some reductions from this tree:
| Reduction | Transformation |
|---|---|
| Hamiltonian cycle $\to$ TSP | edges cost 1, non-edges cost 2. A Hamiltonian cycle exists iff the optimal tour costs exactly $n$ |
| Independent set $\to$ vertex cover | $S$ is independent iff $V \setminus S$ is a cover, so $\alpha(G) \ge k \iff$ a cover of size $\le n - k$ exists |
| Independent set $\to$ clique | an independent set in $G$ is a clique in the complement graph $\bar G$ |
| Vertex cover $\to$ set cover | elements = edges, sets = for each vertex, its incident edges |
| Partition $\to$ bin packing | items $a_i$, two bins of capacity $\tfrac{1}{2}\sum a_i$: they fit iff the numbers split into two equal halves |
| Subset sum $\to$ knapsack | $v_i = w_i = a_i$, $W = K$. Is the optimum $\ge K$? |
The code implements the first reduction and uses an exact TSP algorithm (Held-Karp dynamic programming) to decide Hamiltonicity of two graphs. The 3-dimensional cube has a Hamiltonian cycle. The Petersen graph famously does not.
from itertools import combinations
def held_karp(D):
"""Exact TSP by dynamic programming over subsets, O(n^2 2^n)."""
n = len(D)
best = {(1, 0): 0} # (visited set as bitmask, last city)
for size in range(2, n + 1):
for S in combinations(range(1, n), size - 1):
mask = 1 | sum(1 << k for k in S)
for j in S:
prev = mask ^ (1 << j)
best[(mask, j)] = min(best[(prev, i)] + D[i][j]
for i in [0, *S] if (prev, i) in best and i != j)
full = (1 << n) - 1
return min(best[(full, j)] + D[j][0] for j in range(1, n))
def ham_cycle_to_tsp(n, edges):
"""Reduction: weight 1 on graph edges, 2 on non-edges. G is Hamiltonian <=> OPT = n."""
E = {frozenset(e) for e in edges}
return [[0 if i == j else (1 if frozenset((i, j)) in E else 2) for j in range(n)]
for i in range(n)]
cube = [(0, 1), (1, 3), (3, 2), (2, 0), (4, 5), (5, 7), (7, 6), (6, 4), (0, 4), (1, 5), (2, 6), (3, 7)]
petersen = ([(i, (i + 1) % 5) for i in range(5)] + [(i, i + 5) for i in range(5)]
+ [(5 + i, 5 + (i + 2) % 5) for i in range(5)])
for name, n, E in [("cube (8 nodes)", 8, cube), ("Petersen (10 nodes)", 10, petersen)]:
opt = held_karp(ham_cycle_to_tsp(n, E))
print(f"{name:20s}: optimal tour = {opt:2d}, n = {n} -> "
f"{'Hamiltonian' if opt == n else 'not Hamiltonian'}")
# cube (8 nodes) : optimal tour = 8, n = 8 -> Hamiltonian
# Petersen (10 nodes) : optimal tour = 11, n = 10 -> not Hamiltonian
The transformation takes $O(n^2)$ time. So if the TSP had a polynomial algorithm, Hamiltonian cycle would too. Since Hamiltonian cycle is NP-complete, the TSP is NP-hard. Notice that the instances produced are metric (costs 1 and 2 satisfy the triangle inequality), so even the metric TSP is NP-hard.
Key takeaways
- $A \le_p B$: a polynomial transformation maps instances of $A$ to instances of $B$ with the same answer.
- To prove $B$ hard, reduce a known hard problem to $B$. The direction is the most common mistake.
- Reductions also transfer algorithms: matching via max flow, vertex cover via set cover.
- Hamiltonian cycle reduces to the TSP with costs 1 and 2, so even the metric TSP is NP-hard.
Practice
- (Beginner) Using the independent set to vertex cover reduction, find a maximum independent set of a 5-cycle from its minimum vertex cover.
- (Intermediate) Implement the partition to bin packing reduction and use
first_fit_decreasing(or an exact method) to test the numbers ${3, 1, 1, 2, 2, 1}$ and ${5, 3, 3, 2}$.- (Interview) "You found a polynomial algorithm that reduces your scheduling problem to a max-flow problem. A colleague says your problem is therefore easy. Another found a reduction of your problem to 3-SAT and says it is therefore hard. Who is right?" Hint: check the direction of each reduction and what each one proves.
Pseudo-Polynomial Algorithms
Why it matters. The knapsack DP of Section 4.3 runs in $O(nW)$ time and solves large instances quickly, yet knapsack is NP-hard. These facts are consistent because of how input size is measured. The distinction between weakly and strongly NP-hard problems tells you whether a DP over numeric values can help. It also explains why some problems allow arbitrarily good approximations and others do not.
Theory. An instance's size counts the bits of its numbers: the capacity $W$ takes about $\log_2 W$ bits. An algorithm is pseudo-polynomial if its running time is polynomial in $n$ and in the magnitude of the largest number, which can be exponential in its number of bits. Adding one digit to $W$ multiplies the work of the knapsack DP by 10, although the input grew by only about 3.3 bits.
| Category | Definition | Examples |
|---|---|---|
| Weakly NP-hard | NP-hard, but solvable in pseudo-polynomial time | knapsack, subset sum, partition, scheduling 2 parallel machines |
| Strongly NP-hard | stays NP-hard even when all numbers are bounded by a polynomial in $n$ | bin packing, 3-partition, TSP, set cover, coloring, job shop |
If a problem is strongly NP-hard, it has no pseudo-polynomial algorithm unless P = NP. Its difficulty is combinatorial, not numerical.
Practical consequences.
- Weakly NP-hard problems with moderate numbers (budgets in thousands, capacities in hundreds) are easy in practice with DP.
- Scaling and rounding the numbers trades accuracy for speed. This is the basis of the knapsack FPTAS (next topic): round values to multiples of $\varepsilon v_{\max} / n$, then run a DP over values.
- For strongly NP-hard problems, no such trick exists, and approximation guarantees have limits (Section "Approximation Ratios and Inapproximability").
import numpy as np
def knapsack_dp(w, v, W):
"""0-1 knapsack DP over capacities 0..W: O(n W) time, one row of W + 1 cells."""
best = np.zeros(W + 1, dtype=np.int64)
for wi, vi in zip(w, v):
best[wi:] = np.maximum(best[wi:], best[:-wi] + vi) # take item i (vectorized over c)
return best[W]
w = np.array([5, 4, 6, 2, 3, 7, 4]); v = np.array([10, 7, 12, 3, 6, 13, 6]); W = 15
print(f"{'scale':>9} {'capacity W':>11} {'input bits':>11} {'DP cells n*W':>14} {'optimum':>8}")
for k in range(6):
s = 10 ** k # same instance, weights written with more digits
bits = int(sum(np.ceil(np.log2(w * s + 1))) + np.ceil(np.log2(W * s + 1)))
print(f"{s:>9,} {W * s:>11,} {bits:>11} {len(w) * (W * s + 1):>14,} {knapsack_dp(w * s, v, W * s):>8}")
# scale capacity W input bits DP cells n*W optimum
# 1 15 23 112 29
# 10 150 49 1,057 29
# 100 1,500 75 10,507 29
# 1,000 15,000 100 105,007 29
# 10,000 150,000 129 1,050,007 29
# 100,000 1,500,000 155 10,500,007 29
This is the sprint instance from Section 4.3. Multiplying all weights by $10^k$ changes nothing about the problem: the optimum stays 29. The input grows linearly (about 26 more bits per factor of 10), but the DP work grows by a factor of 10 each time, which is exponential in the input size.
Key takeaways
- Pseudo-polynomial means polynomial in the numeric values, which is exponential in their bit length.
- Weakly NP-hard problems (knapsack, subset sum, partition) have pseudo-polynomial DPs and are easy when numbers are small.
- Strongly NP-hard problems (bin packing, TSP, 3-partition) stay hard even with small numbers.
- Scaling and rounding numbers turns pseudo-polynomial DPs into approximation schemes.
Practice
- (Beginner) Write the DP for subset sum ("does a subset of ${a_i}$ sum to exactly $K$?") with a boolean table. What is its running time?
- (Intermediate) Implement the knapsack DP over values instead of weights: $\text{minweight}[j][p]$ = least weight to reach value $p$ with the first $j$ items. When is it faster than the weight-based DP?
- (Interview) "Your knapsack DP worked well in testing but times out in production, although the item count is the same. What changed?" Hint: look at the units. A budget in cents instead of euros, or capacities in bytes instead of MB, multiplies $W$.
Approximation Algorithms
Why it matters. For large NP-hard problems, waiting for a proven optimum is often not an option. Heuristics are fast but give no guarantee. Approximation algorithms are fast and come with a mathematical promise: for every instance, the solution is within a known factor of optimal. They also provide design ideas, such as LP rounding and primal-dual methods, that good practical heuristics build on.
Theory. For a minimization problem, a polynomial-time algorithm is a $\rho$-approximation ($\rho \ge 1$) if, for every instance $I$,
$$ \text{ALG}(I) \le \rho \cdot \text{OPT}(I). $$
For maximization, $\text{ALG}(I) \ge \text{OPT}(I) / \rho$. The proof never needs OPT itself. It compares ALG with a lower bound on OPT (for minimization), just as branch and bound does.
Example: vertex cover (choose a minimum set of nodes touching every edge).
- Matching algorithm. Repeatedly pick an uncovered edge and add both endpoints. The picked edges form a matching $M$, and the cover has $2\lvert M \rvert$ nodes. Matching edges share no endpoints, so any cover needs a different node for each of them: $\text{OPT} \ge \lvert M \rvert$. Hence ALG $\le 2,\text{OPT}$.
- LP rounding. Solve the LP relaxation ($x_i + x_j \ge 1$, $0 \le x \le 1$) and take every node with $x_v \ge \tfrac{1}{2}$. Each edge has an endpoint with $x \ge \tfrac{1}{2}$, so this is a cover. Its size is at most $2 \sum_v x_v = 2, z_{LP} \le 2,\text{OPT}$.
| Problem | Algorithm | Ratio |
|---|---|---|
| Vertex cover | maximal matching, or LP rounding | 2 |
| Set cover | greedy (Section 4.3) | $H(m) \approx \ln m$ |
| Metric TSP | double tree: MST, double edges, Euler tour, shortcut | 2 |
| Metric TSP | Christofides: MST + min-weight perfect matching on odd-degree nodes, Euler tour, shortcut | 1.5 |
| Steiner tree | MST on the metric closure of the terminals (Section 3.3) | 2 |
| Bin packing | First Fit Decreasing | $\tfrac{11}{9}\text{OPT} + \tfrac{6}{9}$ |
| Knapsack | FPTAS: round values, DP over values | $1 + \varepsilon$, time poly$(n, 1/\varepsilon)$ |
Approximation schemes. A PTAS gives a $(1 + \varepsilon)$-approximation for every fixed $\varepsilon > 0$, in time polynomial in $n$ (but possibly like $n^{1/\varepsilon}$). An FPTAS is also polynomial in $1/\varepsilon$. Knapsack has an FPTAS. Euclidean TSP has a PTAS. Strongly NP-hard problems have no FPTAS unless P = NP.
Christofides for metric TSP
1. minimum spanning tree T cost(T) <= OPT
2. O = nodes with odd degree in T (always an even number)
3. minimum-weight perfect matching M on O cost(M) <= OPT / 2
4. T + M has all degrees even -> Euler tour cost <= 1.5 OPT
5. shortcut repeated nodes (triangle inequality) tour cost <= 1.5 OPT
import numpy as np
from scipy.optimize import milp, LinearConstraint, Bounds
def vertex_cover_instances(seed, n=30, p=0.12):
rng = np.random.default_rng(seed)
return n, [(i, j) for i in range(n) for j in range(i + 1, n) if rng.random() < p]
def matching_cover(n, edges):
"""2-approximation: take BOTH ends of every edge in a maximal matching."""
cover = set()
for i, j in edges:
if i not in cover and j not in cover: # edge still uncovered -> add to matching
cover |= {i, j}
return cover
def lp_rounding_cover(n, edges):
"""2-approximation: solve the LP relaxation, keep every vertex with x_v >= 1/2."""
A = np.zeros((len(edges), n))
for k, (i, j) in enumerate(edges): A[k, i] = A[k, j] = 1 # x_i + x_j >= 1
lp = milp(np.ones(n), constraints=LinearConstraint(A, lb=1), bounds=Bounds(0, 1))
return {v for v in range(n) if lp.x[v] >= 0.5 - 1e-9}, lp.fun
def exact_cover(n, edges):
A = np.zeros((len(edges), n))
for k, (i, j) in enumerate(edges): A[k, i] = A[k, j] = 1
return milp(np.ones(n), constraints=LinearConstraint(A, lb=1), bounds=Bounds(0, 1),
integrality=np.ones(n)).fun
print(f"{'seed':>4} {'edges':>5} {'LP':>6} {'OPT':>4} {'matching':>9} {'LP-round':>9} {'ratios':>13}")
for seed in range(5):
n, E = vertex_cover_instances(seed)
m_cov = matching_cover(n, E)
r_cov, lp_val = lp_rounding_cover(n, E)
opt = exact_cover(n, E)
assert all(i in m_cov or j in m_cov for i, j in E) # both are valid covers
print(f"{seed:>4} {len(E):>5} {lp_val:>6.1f} {opt:>4.0f} {len(m_cov):>9} {len(r_cov):>9}"
f" {len(m_cov) / opt:.2f} {len(r_cov) / opt:.2f}")
# seed edges LP OPT matching LP-round ratios
# 0 49 14.5 16 24 22 1.50 1.38
# 1 51 14.0 14 22 14 1.57 1.00
# 2 53 14.0 14 20 14 1.43 1.00
# 3 46 15.0 15 24 15 1.60 1.00
# 4 43 14.5 15 26 25 1.73 1.67
Both algorithms stay below the guaranteed factor of 2, as the proofs require. On seeds 1-3 the LP bound equals OPT, and LP rounding returns an optimal cover. When the LP uses many halves (seeds 0 and 4), rounding them all up costs more. The LP column also gives each instance its own a-posteriori guarantee: on seed 4, the matching cover is at most $26 / 14.5 \approx 1.79$ times the optimum, even without knowing OPT.
Key takeaways
- A $\rho$-approximation is polynomial and always within factor $\rho$ of optimal. The proof compares it with a bound on OPT.
- Standard techniques: greedy, LP rounding, primal-dual, and DP on rounded data.
- Classic results: vertex cover 2, set cover $\ln m$, metric TSP 1.5 (Christofides), knapsack FPTAS.
- PTAS and FPTAS reach $1 + \varepsilon$ for any $\varepsilon$. An FPTAS is only possible for weakly NP-hard problems (unless P = NP).
Practice
- (Beginner) Show that the matching algorithm can be exactly a factor 2 off: find a small graph where it returns twice the optimal cover.
- (Intermediate) Implement the double-tree 2-approximation for metric TSP (MST with Prim's algorithm, a depth-first preorder walk, then return to the start) and compare it with the optimal tour of the 10-city TSP instance in Section 4.3.
- (Interview) "Why would anyone use a 2-approximation algorithm when local search heuristics usually get within 5% of optimal?" Hint: think about worst-case guarantees, predictable behavior on adversarial inputs, approximations as starting points for local search, and lower bounds as a by-product.
Approximation Ratios and Inapproximability
Why it matters. Not every NP-hard problem can be approximated equally well. Knapsack can be approximated to any precision, vertex cover to a factor of 2, set cover only to about $\ln n$, and general TSP to no constant factor at all. Knowing these limits tells you when to stop searching for a better guaranteed algorithm and switch to heuristics, problem-specific structure, or exact methods for small instances.
Theory. For a minimization problem, the approximation ratio of an algorithm is the worst case of $\text{ALG}(I) / \text{OPT}(I)$ over all instances. A ratio is tight if instances exist where it is (nearly) reached. Two kinds of guarantee are used in practice:
| A-priori ratio | A-posteriori gap | |
|---|---|---|
| Source | worst-case proof over all instances | bound computed for this instance (LP, Lagrangian) |
| Example | Christofides is $\le 1.5 \cdot$ OPT | "this tour is within 2.3% of the LP bound" |
| Typical value | loose (e.g. 1.5) | usually much tighter |
Inapproximability results show that, unless P = NP, no polynomial algorithm achieves a ratio better than some threshold. The standard tool is a gap reduction: a reduction from an NP-complete problem that produces instances whose optimum is either small or far larger, with nothing in between. An algorithm with a good enough ratio would tell the two cases apart, and so would solve the NP-complete problem.
General TSP has no constant-factor approximation. Given a graph $G$ with $n$ nodes and any constant $\rho$, set edge costs to 1 and non-edge costs to $\rho n + 1$.
- If $G$ has a Hamiltonian cycle, $\text{OPT} = n$, so a $\rho$-approximation returns a tour of cost $\le \rho n$.
- If not, every tour uses a non-edge, so every tour costs $> \rho n$.
The approximate tour's cost would decide Hamiltonian cycle. The construction breaks the triangle inequality, which is why the metric TSP still allows Christofides' 1.5.
from itertools import combinations
def held_karp(D): # exact TSP, as in "Polynomial Reductions"
n = len(D); best = {(1, 0): 0}
for size in range(2, n + 1):
for S in combinations(range(1, n), size - 1):
mask = 1 | sum(1 << k for k in S)
for j in S:
prev = mask ^ (1 << j)
best[(mask, j)] = min(best[(prev, i)] + D[i][j]
for i in [0, *S] if (prev, i) in best and i != j)
return min(best[((1 << n) - 1, j)] + D[j][0] for j in range(1, n))
def gap_instance(n, edges, rho):
"""Edges cost 1, non-edges cost rho*n + 1: a rho-approximation would decide Ham. cycle."""
E = {frozenset(e) for e in edges}
big = rho * n + 1
return [[0 if i == j else (1 if frozenset((i, j)) in E else big) for j in range(n)]
for i in range(n)]
cube = [(0, 1), (1, 3), (3, 2), (2, 0), (4, 5), (5, 7), (7, 6), (6, 4), (0, 4), (1, 5), (2, 6), (3, 7)]
petersen = ([(i, (i + 1) % 5) for i in range(5)] + [(i, i + 5) for i in range(5)]
+ [(5 + i, 5 + (i + 2) % 5) for i in range(5)])
rho = 3
for name, n, E in [("cube", 8, cube), ("Petersen", 10, petersen)]:
opt = held_karp(gap_instance(n, E, rho))
if opt == n:
verdict = f"Hamiltonian: any {rho}-approximation returns a tour of cost <= {rho * n}"
else:
verdict = f"not Hamiltonian: every tour costs > {rho * n}"
print(f"{name:9s}: OPT = {opt:3d} -> {verdict}")
# cube : OPT = 8 -> Hamiltonian: any 3-approximation returns a tour of cost <= 24
# Petersen : OPT = 40 -> not Hamiltonian: every tour costs > 30
There is no instance with an optimum between 24 and 30. That empty band is the "gap" a 3-approximation would expose.
Known thresholds (assuming P $\neq$ NP unless noted):
| Problem | Best known ratio | Lower bound on achievable ratio |
|---|---|---|
| Knapsack | $1 + \varepsilon$ (FPTAS) | none: arbitrarily close to 1 |
| Euclidean TSP | $1 + \varepsilon$ (PTAS) | none |
| Bin packing | $\tfrac{11}{9}$ asymptotically (FFD), better schemes exist | $\tfrac{3}{2}$ in absolute terms (from partition: 2 bins or 3) |
| Metric TSP | 1.5 (Christofides; marginally below since 2020) | $\tfrac{123}{122}$ |
| Vertex cover | 2 | $\sqrt{2} - \varepsilon \approx 1.41$, and $2 - \varepsilon$ under the Unique Games Conjecture |
| Set cover | $\ln n$ (greedy) | $(1 - \varepsilon) \ln n$ |
| General TSP | none | no constant ratio at all |
| Max clique, graph coloring | about $n / \text{polylog}(n)$ | $n^{1 - \varepsilon}$ |
These results organize problems into classes: FPTAS $\subset$ PTAS $\subset$ APX (constant ratio) $\subset$ log-APX $\subset$ poly-APX. Knapsack sits at the top, set cover in log-APX, coloring and clique near the bottom.
Key takeaways
- The approximation ratio is a worst case over all instances. A-posteriori gaps from bounds are instance-specific and usually much tighter.
- Gap reductions prove inapproximability: a good enough approximation would solve an NP-complete problem.
- General TSP has no constant-factor approximation. Metric TSP has 1.5. Euclidean TSP has a PTAS.
- Problems range from arbitrarily approximable (knapsack) to nearly inapproximable (clique, coloring). Knowing where yours sits guides the method choice.
Practice
- (Beginner) Use the partition reduction to explain why no polynomial algorithm for bin packing can guarantee fewer than $\tfrac{3}{2}\text{OPT}$ bins (unless P = NP).
- (Intermediate) Modify the gap construction so that non-edges cost 2 instead of $\rho n + 1$. Which approximation ratio can this version rule out? Why does that match the fact that this instance is metric?
- (Interview) "Your graph coloring heuristic for exam scheduling uses 14 slots. Is it worth funding research into an algorithm with a better worst-case guarantee?" Hint: compare what the inapproximability results say about worst cases with the a-posteriori bound (largest clique) on your instances.
<a id="5-nonlinear-optimization"></a>