11. Heuristics and Metaheuristics
Chapter 4 showed that many combinatorial problems are NP-hard. Exact methods can still prove optimality, but their running time can explode, and many real decisions must be made in seconds. This chapter covers the methods that give up the proof in exchange for speed. Constructive heuristics build a solution in one pass, and local search improves it step by step. Metaheuristics steer local search away from poor local optima, and hybrid methods combine heuristics with integer programming, constraint programming and machine learning. Most production systems for routing, scheduling and assignment run some combination of these ideas, so knowing how to design, compare and tune them is a core OR skill.
<a id="111-constructive-and-local-search-heuristics"></a>
11.1 Constructive and Local Search Heuristics
A heuristic finds good solutions quickly, using rules that work well in practice but prove nothing about optimality. This subchapter covers the two building blocks of almost every heuristic method. Constructive heuristics build a complete solution from scratch, one decision at a time. Local search starts from a complete solution and keeps applying small changes that improve it. Most examples share one running instance: 80 cities placed uniformly at random in a 100 x 100 square (seed 42). Its optimal tour has length 685.27. The subtour-cut loop of Section 4.3 proves this in about two seconds, so every heuristic below can be measured by its exact gap to the optimum.
Greedy Heuristics
Why it matters. Greedy is the first heuristic most people write, and often the one that ships. It builds a solution by always taking the choice that looks best right now. It is used to produce starting solutions for local search, as a baseline that every smarter method must beat, and as the decision rule itself when there are only milliseconds to decide (dispatching, ad allocation, task assignment). Some problems are solved exactly by greedy, while on others it fails badly. Knowing which situation you are in, and how to make greedy less short-sighted, is essential.
Theory. A greedy heuristic needs three ingredients: a set of candidate elements (items, edges, assignments), a greedy function $g$ that scores how attractive a candidate is right now, and a feasibility test. The template is:
S = empty partial solution
while S is not complete:
C = candidates c such that S + c can still be completed feasibly
c* = the candidate in C with the best greedy score g(c, S)
S = S + c* # never undone: no backtracking
Greedy is fast (often $O(n \log n)$ to $O(n^2)$) and simple. Its weakness is that decisions are never revisited.
When greedy is optimal. For some problems a locally best choice can always be extended to a global optimum. The proof is usually an exchange argument: take any optimal solution that differs from the greedy one, and swap in the greedy choice without making it worse. Examples are Kruskal's and Prim's algorithms for minimum spanning trees (Section 3.3), Dijkstra's algorithm (Section 3.2) and the SPT rule for $1 ,||, \sum C_j$ (Section 10.3). The general theory is that greedy is optimal for every weight function exactly when the feasible sets form a matroid.
When greedy fails. On most NP-hard problems, greedy has three typical failure modes:
- Myopia. An early choice that looks cheap uses up a resource that a later element needed more. For the knapsack problem, greedy can be arbitrarily bad (Section 4.3).
- Dead ends. On tightly constrained problems, greedy can run into a state where no candidate is feasible, even though feasible solutions exist. Greedy has no way to back out.
- Order dependence. The result can depend heavily on the processing order or on tie-breaking, as with greedy graph coloring (Section 4.3).
Making greedy smarter. Several standard refinements keep the speed but reduce the myopia:
| Refinement | Idea | Example |
|---|---|---|
| Better greedy function | normalize the benefit by the resource it consumes | value density $v_j / w_j$ for knapsack |
| Dynamic ordering | recompute scores after every decision | DSatur for graph coloring |
| Regret | first decide the element that loses most if it does not get its best option | regret insertion in vehicle routing |
| Look-ahead (pilot method) | score each candidate by completing the solution greedily from it | expensive but much stronger |
| Randomization | choose at random among the best few candidates | GRASP (Section 11.2) |
The regret of an element $j$ is the difference between the cost of its second-best and its best feasible option:
$$ \text{regret}_j = c_j^{(2)} - c_j^{(1)}. $$
A large regret means that postponing $j$ is dangerous: if its best option is taken by someone else, the cost jumps. An analogy is a sports draft. A good team does not simply pick the best player available. It picks the player whose position will be hardest to fill later.
Example: assigning jobs to servers. The generalized assignment problem (GAP) assigns each of $n$ tasks to one of $m$ agents (servers, technicians, vehicles). Giving task $j$ to agent $i$ costs $c_{ij}$ and uses $a_{ij}$ of the agent's capacity $b_i$:
$$ \min \sum_{i,j} c_{ij} x_{ij} \quad \text{s.t.} \quad \sum_i x_{ij} = 1 \ \ \forall j, \qquad \sum_j a_{ij} x_{ij} \le b_i \ \ \forall i, \qquad x_{ij} \in {0, 1}. $$
Even deciding whether a feasible assignment exists is NP-hard. The code compares three greedy variants with the exact optimum on six random instances with 5 agents and 40 tasks. The capacities are tight: each agent can take only its fair share of the work.
import numpy as np
from scipy.optimize import milp, LinearConstraint, Bounds
def instance(seed, m=5, n=40):
r = np.random.default_rng(seed)
a = r.integers(5, 26, size=(m, n)) # a[i, j]: capacity task j uses on agent i
c = r.integers(10, 51, size=(m, n)) # c[i, j]: cost of giving task j to agent i
b = (a.sum(axis=1) / m).astype(int) # tight: each agent can take a fair share only
return a, c, b
def greedy(a, c, b, order):
"""Assign tasks in the given order, each to its cheapest agent with room left."""
room, cost = b.copy(), 0
for j in order:
ok = [i for i in range(len(b)) if a[i, j] <= room[i]]
if not ok: return None # dead end: no agent can take task j
i = min(ok, key=lambda i: c[i, j])
room[i] -= a[i, j]; cost += c[i, j]
return cost
def regret_greedy(a, c, b):
"""Assign next the task that would lose most if it missed its best agent."""
room, cost, left = b.copy(), 0, set(range(a.shape[1]))
while left:
best = None
for j in left:
ok = sorted(c[i, j] for i in range(len(b)) if a[i, j] <= room[i])
if not ok: return None
regret = ok[1] - ok[0] if len(ok) > 1 else np.inf # one option left: assign now
if best is None or regret > best[0]: best = (regret, j)
j = best[1]
i = min((i for i in range(len(b)) if a[i, j] <= room[i]), key=lambda i: c[i, j])
room[i] -= a[i, j]; cost += c[i, j]; left.remove(j)
return cost
def optimum(a, c, b):
m, n = a.shape # x[i, j] flattened row by row
assign = np.kron(np.ones(m), np.eye(n)) # each task to exactly one agent
cap = np.kron(np.eye(m), np.ones(n)) * a.ravel()
res = milp(c.ravel(), integrality=np.ones(m * n), bounds=Bounds(0, 1),
constraints=[LinearConstraint(assign, 1, 1), LinearConstraint(cap, ub=b)])
return res.fun
fmt = lambda v, opt: " infeasible" if v is None else f"{v:5d} ({v / opt - 1:+5.1%})"
print("seed index order cheapest first regret optimum")
for seed in range(6):
a, c, b = instance(seed)
opt = optimum(a, c, b)
by_cost = np.argsort(c.min(axis=0)) # tasks with the cheapest option first
print(f"{seed:4d} {fmt(greedy(a, c, b, range(a.shape[1])), opt):17s} "
f"{fmt(greedy(a, c, b, by_cost), opt):17s} {fmt(regret_greedy(a, c, b), opt):17s} {opt:5.0f}")
# seed index order cheapest first regret optimum
# 0 756 (+13.9%) infeasible 678 (+2.1%) 664
# 1 infeasible infeasible infeasible 692
# 2 682 (+10.9%) infeasible 619 (+0.7%) 615
# 3 infeasible infeasible 667 (+0.9%) 661
# 4 749 (+11.1%) 720 (+6.8%) 675 (+0.1%) 674
# 5 678 (+5.9%) 680 (+6.2%) 646 (+0.9%) 640
Plain greedy in index order runs into a dead end on two instances. Taking the tasks with the cheapest option first is even worse: it fails on four of the six, because tasks with cheap options grab capacity early and leave the expensive-to-place tasks with nowhere to go. Regret greedy finds a feasible solution on five instances and stays within 2.1% of the optimum. Seed 1 shows the limit of any one-pass method: a feasible assignment exists (the MILP finds cost 692), but even regret greedy commits early to choices it cannot undo.
Key takeaways
- Greedy builds a solution by repeatedly adding the best-scoring feasible candidate and never revisits a decision.
- Greedy is optimal on matroid-like problems (spanning trees, SPT), which an exchange argument proves. On NP-hard problems it is a heuristic.
- Its failure modes are myopia, dead ends on tight constraints, and sensitivity to order and ties.
- Regret, dynamic re-scoring, look-ahead and randomization make greedy much stronger at a modest cost.
Practice
- (Beginner) Coins have values 1, 3 and 4. Show that the greedy rule "always take the largest coin that fits" does not minimize the number of coins for an amount of 6. Then explain why greedy is optimal for the coin values 1, 5, 10 and 25.
- (Intermediate) Extend the GAP example with a repair step: when a task has no feasible agent, move one already assigned task to another agent with room, then retry. On how many of the six instances does the plain greedy then become feasible, and at what cost?
- (Interview) "Our ride-hailing app assigns each request, the moment it arrives, to the nearest free driver. When is this far from optimal?" Hint: think about two requests arriving a second apart that both want the same driver. Consider batching requests for a few seconds and solving an assignment problem (Section 3.6), and what regret would say about a request in a remote area.
Nearest Neighbor and Insertion Heuristics
Why it matters. Tour construction heuristics for the TSP are the classic examples of constructive heuristics. They produce the starting tours for local search, the routes inside vehicle routing heuristics, and quick estimates of route length for planning. They take milliseconds, and a good choice among them can halve the gap that the local search must close later.
Theory. Nearest neighbor (NN) starts at some city and repeatedly moves to the closest unvisited city, then returns to the start. It runs in $O(n^2)$ time. Its flaw shows at the end: cities that were skipped early must be collected at the end with long detours, and the final edge back to the start can be very long. Its worst-case ratio grows like $\log n$ (Rosenkrantz, Stearns and Lewis, 1977).
Insertion heuristics keep a closed subtour and insert the remaining cities one at a time. Each step makes two decisions:
- Which city to insert (the selection rule): the city nearest to the tour, the farthest one, a random one, or, in cheapest insertion, the city whose best insertion is cheapest overall.
- Where to insert it: between the consecutive tour cities $a$ and $b$ that minimize the added length
$$ \Delta(a, k, b) = d(a, k) + d(k, b) - d(a, b). $$
inserting city k into edge (a, b)
a --------------- b a b
-> \ /
k \ /
---- k -----
added length: d(a,k) + d(k,b) - d(a,b)
Nearest and cheapest insertion are 2-approximations for metric instances. Farthest insertion has no comparable worst-case guarantee, yet it is usually the best of the family in practice. Inserting the farthest cities first sketches the overall shape of the tour, much like the convex hull, and the remaining cities are then local corrections. Nearest insertion grows a compact blob instead, and the far cities at the end force expensive detours.
| Method | Time | Worst-case ratio (metric) | Typical behavior on random instances |
|---|---|---|---|
| Nearest neighbor | $O(n^2)$ | $\Theta(\log n)$ | about 25% above optimum, long closing edges |
| Nearest / cheapest insertion | $O(n^2)$ / $O(n^2 \log n)$ | 2 | 15-25% |
| Farthest / random insertion | $O(n^2)$ | $O(\log n)$ | often 10-15% |
| Greedy edge (Kruskal-like) | $O(n^2 \log n)$ | $O(\log n)$ | often 15-20% |
| Savings (Clarke-Wright, Section 4.3) | $O(n^2 \log n)$ | $O(\log n)$ | often around 10-12% |
| Christofides (Section 4.4) | $O(n^3)$ | 1.5 | about 10%, slower |
Example: construction heuristics on the running instance. The code also defines the instance, D, n, OPT, length() and nearest_neighbor() that later examples reuse.
import numpy as np
rng = np.random.default_rng(42) # running instance of this chapter
pts = rng.uniform(0, 100, size=(80, 2)) # 80 cities in a 100 x 100 square
D = np.linalg.norm(pts[:, None] - pts[None, :], axis=2)
n, OPT = len(pts), 685.2665 # optimum, from the cutting loop of Section 4.3
def length(tour):
return sum(D[tour[k - 1], tour[k]] for k in range(len(tour)))
def nearest_neighbor(start):
tour, left = [start], set(range(n)) - {start}
while left:
nxt = min(left, key=lambda j: D[tour[-1], j]) # closest unvisited city
tour.append(nxt); left.remove(nxt)
return tour
def insertion(rule, seed=0):
"""Generic insertion: pick a city by `rule`, insert it where it adds least length."""
r = np.random.default_rng(seed)
i, j = np.unravel_index(D.argmax(), D.shape) if rule == "farthest" else (0, int(D[0].argmin()))
tour, left = [int(i), int(j)], set(range(n)) - {int(i), int(j)}
d_tour = np.minimum(D[i], D[j]) # distance from each city to the tour
while left:
cand = list(left)
if rule == "nearest": k = min(cand, key=lambda c: d_tour[c])
elif rule == "farthest": k = max(cand, key=lambda c: d_tour[c])
elif rule == "random": k = cand[r.integers(len(cand))]
else: k = None # cheapest: decided jointly with the position
best = (np.inf, None, None)
for c in ([k] if k is not None else cand):
for p in range(len(tour)): # insert c between tour[p-1] and tour[p]
a, b = tour[p - 1], tour[p]
delta = D[a, c] + D[c, b] - D[a, b]
if delta < best[0]: best = (delta, c, p)
_, c, p = best
tour.insert(p, c); left.remove(c)
d_tour = np.minimum(d_tour, D[c])
return tour
results = {"nearest neighbor (start 0)": nearest_neighbor(0)}
results["nearest neighbor (best start)"] = min((nearest_neighbor(s) for s in range(n)), key=length)
for rule in ("nearest", "cheapest", "random", "farthest"):
results[f"{rule} insertion"] = insertion(rule)
for name, tour in results.items():
L = length(tour)
print(f"{name:30s} {L:7.2f} gap {L / OPT - 1:6.1%}")
# nearest neighbor (start 0) 840.67 gap 22.7%
# nearest neighbor (best start) 771.11 gap 12.5%
# nearest insertion 800.95 gap 16.9%
# cheapest insertion 847.40 gap 23.7%
# random insertion 785.69 gap 14.7%
# farthest insertion 733.50 gap 7.0%
Nearest neighbor from city 0 is 22.7% above the optimum. Trying all 80 start cities and keeping the best brings it to 12.5%, at 80 times the cost. Farthest insertion gives the best tour, only 7.0% above the optimum, while cheapest insertion is the worst here despite its guarantee of a factor 2. The exact ranking changes from instance to instance, but the pattern that far-first and random insertion beat near-first methods on uniform instances is robust.
Key takeaways
- Nearest neighbor is simple and fast, but it leaves the far cities for the end and pays for them with long closing edges.
- Insertion heuristics combine a selection rule with cheapest-position insertion, $\Delta = d(a, k) + d(k, b) - d(a, b)$.
- Farthest and random insertion usually beat nearest and cheapest insertion, even though only the latter have a factor-2 guarantee.
- Construction quality matters mainly as a starting point. Local search (next topics) closes most of the remaining gap.
Practice
- (Beginner) Five cities lie on a line at positions 0, 1, -2, 4 and -8. Run nearest neighbor from position 0 and compute the tour length. Compare it with the optimal tour, whose length is twice the distance between the extreme points.
- (Intermediate) Implement the greedy edge heuristic: sort all edges by length and accept an edge if both endpoints still have degree below 2 and it does not close a cycle too early (use a union-find structure). Compare it with the results of the example.
- (Interview) "Drivers complain that our route app sends them back across town at the end of every shift. What is the likely cause, and what is the cheapest fix?" Hint: recall how nearest neighbor ends. Then consider insertion heuristics and a 2-opt pass (two topics ahead).
Neighborhood Structures
Why it matters. Local search changes a solution step by step, and the neighborhood defines which steps are allowed. It is the most important design decision in local search. It determines which solutions are local optima, how good they are, and how much each step costs. Two local search methods that differ only in their neighborhood can give very different results.
Theory. A neighborhood assigns to every solution $x$ a set $N(x)$ of solutions reachable by one move: flipping a bit, swapping two jobs, replacing two edges of a tour. A solution $x$ is a local optimum with respect to $N$ (for minimization) if
$$ f(x) \le f(y) \quad \text{for all } y \in N(x). $$
Every global optimum is a local optimum for every neighborhood, but not the other way around. A neighborhood is exact if every local optimum is also global. The simplex method works because adjacent vertices form an exact neighborhood for linear programs (convexity), and edge exchange is exact for spanning trees. For the TSP, Papadimitriou and Steiglitz showed that no exact neighborhood can be searched in polynomial time unless P = NP. Practical local search for NP-hard problems always works with inexact neighborhoods.
An analogy is hiking in fog with a goal of reaching the lowest point. You can only see the spots within one step of you (your neighborhood). Short sight means you stop in the first hollow. Longer sight means fewer false bottoms, but each look around takes longer.
The main trade-offs are:
- Size versus power. Larger neighborhoods have fewer and better local optima, but each step evaluates more neighbors.
- Evaluation cost. Good neighborhoods allow delta evaluation: computing $f(y) - f(x)$ from the few changed parts in $O(1)$ instead of re-evaluating $f$ from scratch.
- Connectivity. Every solution, and in particular the optimum, should be reachable by a sequence of moves.
- Feasibility. If moves often break constraints, the search either needs compound moves or must accept infeasible solutions with penalties.
| Solution type | Neighborhood | Move | Size |
|---|---|---|---|
| Binary vector | $k$-flip | flip $k$ bits | $\binom{n}{k}$ |
| Assignment (GAP) | shift / swap | move one task to another agent / exchange two tasks | $n(m-1)$ / $O(n^2)$ |
| Sequence (scheduling) | adjacent swap | exchange two neighbors | $n - 1$ |
| Sequence | swap | exchange any two jobs | $n(n-1)/2$ |
| Sequence | insertion | remove a job, reinsert it elsewhere | $(n-1)^2$ |
| TSP tour | 2-opt / 3-opt | replace 2 / 3 edges | $O(n^2)$ / $O(n^3)$ |
| Vehicle routes | relocate, exchange, 2-opt* | move customers or route tails between routes | $O(n^2)$ |
| Job shop | critical-path moves | swap adjacent operations on a critical block | small |
Neighborhoods that are exponentially large but can still be searched efficiently, for example by dynamic programming or shortest paths, are called very large-scale neighborhoods (VLSN). Large Neighborhood Search (Section 11.2) is the most widely used example.
Example: counting local optima. For 8 jobs on one machine (minimize total weighted tardiness), all $8! = 40{,}320$ sequences can be enumerated. That makes it possible to count the local optima of each neighborhood exactly.
import numpy as np
from itertools import permutations
r = np.random.default_rng(3) # 8 jobs on one machine
p = r.integers(1, 11, size=8) # processing times
w = r.integers(1, 6, size=8) # weights
d = r.integers(5, int(p.sum() * 0.7), size=8) # due dates
def twt(seq): # total weighted tardiness
t, total = 0, 0
for j in seq:
t += p[j]; total += w[j] * max(0, t - d[j])
return total
f = {s: twt(s) for s in permutations(range(8))} # objective of all 8! = 40,320 sequences
best = min(f.values())
def adjacent_swap(s):
for i in range(len(s) - 1):
yield s[:i] + (s[i + 1], s[i]) + s[i + 2:]
def any_swap(s):
for i in range(len(s)):
for j in range(i + 1, len(s)):
t = list(s); t[i], t[j] = t[j], t[i]; yield tuple(t)
def insertion(s): # remove one job, reinsert it elsewhere
for i in range(len(s)):
rest = s[:i] + s[i + 1:]
for k in range(len(s)):
if k != i: yield rest[:k] + (s[i],) + rest[k:]
def swap_or_insert(s):
yield from any_swap(s); yield from insertion(s)
print(f"global optimum TWT = {best}")
for name, N in [("adjacent swap", adjacent_swap), ("any swap", any_swap),
("insertion", insertion), ("swap + insertion", swap_or_insert)]:
size = len(set(N(next(iter(f)))))
loc = [s for s in f if all(f[t] >= f[s] for t in N(s))] # no strictly better neighbor
worst = max(f[s] for s in loc)
print(f"{name:17s} |N| = {size:2d} local optima = {len(loc):5d} "
f"worst local optimum = {worst:4d} ({worst / best - 1:+.0%})")
# global optimum TWT = 46
# adjacent swap |N| = 7 local optima = 82 worst local optimum = 184 (+300%)
# any swap |N| = 28 local optima = 6 worst local optimum = 51 (+11%)
# insertion |N| = 49 local optima = 2 worst local optimum = 46 (+0%)
# swap + insertion |N| = 70 local optima = 2 worst local optimum = 46 (+0%)
With adjacent swaps (7 neighbors), the landscape has 82 local optima, and the worst is 300% above the optimum. Allowing any swap (28 neighbors) leaves only 6 local optima, all within 11%. Insertion (49 neighbors) is exact on this instance: its two local optima are both global (two sequences tie for the optimum). Adding swaps to insertion does not help further and makes each step more expensive. This is typical for scheduling: insertion moves fit the structure of the problem, because moving one job shifts a whole block of others.
Key takeaways
- The neighborhood defines the moves and therefore the local optima. Changing it changes the landscape.
- Larger neighborhoods give fewer and better local optima but cost more per step. Delta evaluation keeps steps cheap.
- A good neighborhood matches the structure of the problem: insertion for sequencing, 2-opt for tours, relocate and exchange for routes.
- Neighborhoods must be connected, and constraints may require compound moves or penalties.
Practice
- (Beginner) For $n = 10$ jobs, compute the sizes of the adjacent swap, swap and insertion neighborhoods. Show that every adjacent swap is also a swap and also an insertion move.
- (Intermediate) Modify the example to run steepest descent from every one of the 40,320 sequences, and compute for each neighborhood the fraction of starts that end at a global optimum (the size of its basin of attraction).
- (Interview) "Our nurse-rostering local search uses the move 'change one shift of one nurse' and gets stuck quickly. What would you try?" Hint: think about coverage constraints that every single change violates. Consider swaps between nurses, compound moves, penalty functions for infeasibility, and ejection chains.
Hill Climbing
Why it matters. Hill climbing (called descent or iterative improvement for minimization) is the simplest local search: keep moving to a better neighbor until none exists. It is the core inside almost every metaheuristic, and it already turns poor solutions into decent ones. Its weakness, stopping at the first local optimum, is exactly what the metaheuristics of Section 11.2 address.
Theory.
x = initial solution
repeat:
find y in N(x) with f(y) < f(x) # pivoting rule decides which one
if no such y exists: return x # x is a local optimum
x = y
The pivoting rule decides which improving neighbor is chosen:
- First improvement: scan the neighborhood and take the first improving neighbor found. Steps are cheap, and the search takes many small steps.
- Best improvement (steepest descent): evaluate the whole neighborhood and take the best neighbor. Each step costs a full scan, but fewer steps are needed.
Neither rule dominates. Best improvement is greedy at the level of moves, and greedy can lead into a narrow valley. First improvement, especially with a randomized scan order, behaves more like a random walk downhill and often ends in slightly better local optima.
Other practical points:
- Restarts. Running hill climbing from many random starts and keeping the best (multi-start) is the simplest way to escape a single poor local optimum.
- Plateaus. If many neighbors have equal cost, allowing a limited number of sideways moves helps cross flat regions.
- Running time. The number of steps is usually small in practice, but it can be exponential in the worst case. 2-opt for the TSP has instances that need exponentially many improving steps (Englert, Röglin and Vöcking, 2007).
Example: swap and relocate moves for the TSP. Two simple tour neighborhoods are used: swap exchanges the positions of two cities, and relocate takes one city out and reinserts it elsewhere. Each runs from the same 10 random tours, with both pivoting rules.
import numpy as np
# Reuses D, n, OPT and length() from the Nearest Neighbor and Insertion topic.
def swap_moves(t):
"""Yield (delta, apply) for every exchange of two cities' positions."""
for i in range(n):
for j in range(i + 1, n):
u = t.copy(); u[i], u[j] = u[j], u[i]
changed = {(i - 1) % n, i, (j - 1) % n, j} # edges (t[k], t[k+1]) that change
delta = sum(D[u[k], u[(k + 1) % n]] - D[t[k], t[(k + 1) % n]] for k in changed)
yield delta, lambda u=u: u
def relocate_moves(t):
"""Yield (delta, apply) for every 'take one city out and reinsert it elsewhere'."""
for i in range(n):
a, c, b = t[i - 1], t[i], t[(i + 1) % n]
gain = D[a, c] + D[c, b] - D[a, b] # length saved by removing c
rest = t[:i] + t[i + 1:]
for k in range(n - 1): # insert c between rest[k-1] and rest[k]
u, v = rest[k - 1], rest[k]
if (u, v) == (a, b): continue # that is where it came from
yield D[u, c] + D[c, v] - D[u, v] - gain, lambda r=rest, k=k, c=c: r[:k] + [c] + r[k:]
def hill_climb(t, neighborhood, first):
evals = 0
while True:
best_delta, best_apply = -1e-9, None
for delta, apply in neighborhood(t):
evals += 1
if delta < best_delta:
best_delta, best_apply = delta, apply
if first: break # first improvement: move at once
if best_apply is None:
return t, evals # no improving neighbor: local optimum
t = best_apply()
r = np.random.default_rng(0)
starts = [list(r.permutation(n)) for _ in range(10)]
print(f"random tours: mean gap {np.mean([length(s) for s in starts]) / OPT - 1:.0%}")
for name, N in [("swap", swap_moves), ("relocate", relocate_moves)]:
for first in (True, False):
runs = [hill_climb(s, N, first) for s in starts]
L = [length(t) for t, _ in runs]
print(f"{name:8s} {'first' if first else 'best':5s}: mean gap {np.mean(L) / OPT - 1:5.1%}, "
f"best gap {min(L) / OPT - 1:5.1%}, evaluations per run {np.mean([e for _, e in runs]):9,.0f}")
# random tours: mean gap 489%
# swap first: mean gap 61.2%, best gap 43.8%, evaluations per run 180,507
# swap best : mean gap 70.9%, best gap 54.6%, evaluations per run 273,024
# relocate first: mean gap 26.4%, best gap 15.9%, evaluations per run 412,984
# relocate best : mean gap 30.7%, best gap 24.8%, evaluations per run 658,320
Random tours are almost 6 times the optimal length. Hill climbing with swaps stops 61-71% above the optimum: swapping two cities changes four edges at once, which rarely helps in a tour that is already reasonable. Relocate fits the TSP better and stops at 26-31%. The neighborhood matters much more than the pivoting rule. First improvement was both better and cheaper here than best improvement. Even the best of the 10 relocate runs (15.9%) is worse than farthest insertion (7.0%), which shows that a good neighborhood (next topic) and a good start are both needed.
Key takeaways
- Hill climbing moves to improving neighbors until it reaches a local optimum with respect to its neighborhood.
- First improvement takes cheap steps. Best improvement takes expensive but larger steps. Neither is always better.
- The choice of neighborhood dominates the outcome. Restarts and sideways moves are the simplest remedies for poor local optima.
- Plain hill climbing stops at the first local optimum. Metaheuristics exist to go beyond it.
Practice
- (Beginner) On the integers 0 to 10, with neighbors $x - 1$ and $x + 1$, the function to minimize takes the values $f = (5, 3, 4, 6, 2, 1, 3, 7, 0, 4, 6)$. List the local minima. From which starting points does best-improvement hill climbing reach the global minimum? Which local minima disappear if the neighborhood becomes $x \pm 1$ and $x \pm 2$?
- (Intermediate) Run relocate hill climbing (first improvement) from 100 random tours and tabulate the distribution of final gaps. How many restarts would you need for a 90% chance of a final gap below 20%?
- (Interview) "My teammate says best improvement is always better than first improvement, because it always takes the best move. How do you respond?" Hint: compare the cost per step with the number of steps, consider what greediness does to the path through the landscape, and look at the results above.
2-Opt and k-Opt Moves
Why it matters. The 2-opt neighborhood is the workhorse of TSP and vehicle routing local search. Its extensions, Or-opt, 3-opt and the variable-depth Lin-Kernighan method, are the heart of the best heuristic TSP solvers. LKH, an implementation of Lin-Kernighan, finds tours within 1-2% of the optimum, and often the optimum itself, for instances with many thousands of cities (Section 4.3). The ideas carry over to any problem whose solutions are sequences.
Theory. A 2-opt move removes two edges $(a, b)$ and $(c, d)$ from the tour and reconnects the two paths the only other way that gives a single tour, with edges $(a, c)$ and $(b, d)$. In an array, this reverses the segment from $b$ to $c$.
before: ... a -> b -> ... -> c -> d ... after: ... a -> c -> ... -> b -> d ...
(the path b ... c is traversed in reverse)
a d a d
\ / | |
\ / edges cross | | no crossing
\ / | |
/ \ | |
/ \ | |
c b c b
change in length: delta = d(a,c) + d(b,d) - d(a,b) - d(c,d) (computed in O(1))
For Euclidean instances, a tour with two crossing edges can always be shortened by the 2-opt move that uncrosses them. Every 2-optimal tour is therefore free of crossings. The neighborhood has $n(n-3)/2$ moves.
k-opt. A $k$-opt move removes $k$ edges and reconnects the pieces differently. For 3-opt there are 7 ways to reconnect three paths: 3 are 2-opt moves and 4 are "pure" 3-opt moves. A tour is $k$-optimal if no $k$-opt move improves it. Larger $k$ gives better local optima, but the neighborhood grows like $O(n^k)$. Two restricted variants are popular:
- Or-opt moves a segment of 1 to 3 consecutive cities to another position, possibly reversed. It is a subset of 3-opt with only $O(n^2)$ moves, and it fixes exactly the defects 2-opt cannot: a city that sits on the wrong side of the tour.
- Lin-Kernighan (1973) builds a $k$-opt move step by step, extending a chain of exchanges as long as the cumulative gain stays positive, with $k$ chosen adaptively. This variable-depth idea is what makes it so strong.
Making it fast. Three standard techniques make 2-opt run in near-linear time per pass:
- Neighbor lists. In an improving move, at least one new edge must be shorter than the edge it replaces at its endpoint. So for each city $a$ it is enough to try partners $c$ among $a$'s $K$ nearest cities (typically $K$ = 5 to 10), stopping once $d(a, c) \ge d(a, b)$.
- Don't-look bits. A city whose neighborhood produced no improvement is skipped until one of its tour edges changes.
- Tour representation. Reversing a segment costs $O(n)$ in an array. Very large instances use doubly linked lists or two-level lists.
On random uniform instances, 2-opt, 3-opt and Lin-Kernighan typically end about 5%, 3% and 2% above the optimum (Johnson and McGeoch).
Example: 2-opt, neighbor lists and Or-opt. The code defines two_opt(), or_opt_move() and two_opt_or_opt(), which later examples reuse. Each move is counted as one evaluation.
import numpy as np
# Reuses D, n, OPT, length() and nearest_neighbor() from the Nearest Neighbor and Insertion topic.
NEIGH = np.argsort(D, axis=1)[:, 1:11] # 10 nearest cities of each city
def length(t):
return sum(D[t[k - 1], t[k]] for k in range(n))
def gain_2opt(t, e1, e2):
"""Change in length when edges e1 = (t[e1], t[e1+1]) and e2 are replaced by
(t[e1], t[e2]) and (t[e1+1], t[e2+1])."""
a, b, c, d = t[e1], t[(e1 + 1) % n], t[e2], t[(e2 + 1) % n]
return D[a, c] + D[b, d] - D[a, b] - D[c, d]
def apply_2opt(t, e1, e2):
lo, hi = sorted((e1, e2))
t[lo + 1:hi + 1] = t[lo + 1:hi + 1][::-1] # reverse the path between the two edges
def two_opt(t, neighbor_lists=False):
t, evals, improved = list(t), 0, True
while improved:
improved = False
pos = {city: k for k, city in enumerate(t)}
for i in range(n):
if neighbor_lists: # new edge (t[i], c) with c among the 10 nearest
pairs = [(i, pos[c]) for c in NEIGH[t[i]]] + \
[((i - 1) % n, (pos[c] - 1) % n) for c in NEIGH[t[i]]] # predecessor side
else:
pairs = [(i, j) for j in range(i + 2, n)]
for e1, e2 in pairs:
if abs(e1 - e2) < 2 or {e1, e2} == {0, n - 1}: continue # edges share a city
evals += 1
if gain_2opt(t, e1, e2) < -1e-9:
apply_2opt(t, e1, e2)
pos = {city: k for k, city in enumerate(t)}
improved = True
break # positions changed: move on to the next city
return t, evals
def or_opt_move(t):
"""First improving move of a 1-3 city segment to another place (either orientation)."""
for L in (1, 2, 3):
for i in range(n):
seg = [t[(i + k) % n] for k in range(L)]
p, q = t[i - 1], t[(i + L) % n] # cities on both sides of the segment
saved = D[p, seg[0]] + D[seg[-1], q] - D[p, q]
rest = [c for c in t if c not in seg]
for k in range(len(rest)): # try between rest[k-1] and rest[k]
u, v = rest[k - 1], rest[k]
for s in (seg, seg[::-1]):
if D[u, s[0]] + D[s[-1], v] - D[u, v] < saved - 1e-9:
return rest[:k] + s + rest[k:]
return None
def two_opt_or_opt(t):
"""Alternate the two neighborhoods until the tour is a local optimum of both."""
while True:
t, _ = two_opt(t)
u = or_opt_move(t)
if u is None: return t
t = u
starts = {"nearest neighbor": nearest_neighbor(0)}
starts.update({f"random tour {s}": list(np.random.default_rng(s).permutation(n)) for s in range(3)})
print(f"{'start':16s} {'start gap':>9s} {'2-opt':>7s} {'evals':>7s} {'2-opt NL':>9s} {'evals':>6s} {'+ Or-opt':>9s}")
for name, s in starts.items():
full, ev_full = two_opt(s)
fast, ev_fast = two_opt(s, neighbor_lists=True)
both = two_opt_or_opt(s)
print(f"{name:16s} {length(s) / OPT - 1:9.1%} {length(full) / OPT - 1:7.1%} {ev_full:7,d} "
f"{length(fast) / OPT - 1:9.1%} {ev_fast:6,d} {length(both) / OPT - 1:9.1%}")
# start start gap 2-opt evals 2-opt NL evals + Or-opt
# nearest neighbor 22.7% 8.5% 17,236 6.5% 3,584 3.2%
# random tour 0 475.5% 7.9% 19,665 9.9% 4,430 1.2%
# random tour 1 420.8% 7.1% 30,138 11.6% 4,695 0.8%
# random tour 2 471.4% 14.0% 28,237 12.6% 4,762 8.5%
From any start, 2-opt removes almost all of the initial gap: random tours more than 5 times too long end 7-14% above the optimum. The neighbor-list version reaches similar quality (6.5-12.6%) with 4 to 6 times fewer evaluations, a factor that grows with $n$ because it scans $O(K)$ instead of $O(n)$ partners per city. Adding Or-opt is the biggest single improvement, bringing three of the four tours to 0.8-3.2%. Random tour 2 still ends 8.5% above the optimum: a local optimum of both neighborhoods can still be poor, which motivates the next topic.
Key takeaways
- A 2-opt move replaces two edges and reverses the path between them. Its gain is computed in $O(1)$, and 2-optimal Euclidean tours have no crossings.
- $k$-opt generalizes this at cost $O(n^k)$. Or-opt and Lin-Kernighan capture most of the benefit of larger $k$ cheaply.
- Neighbor lists, don't-look bits and good data structures make these searches fast enough for very large instances.
- Combining neighborhoods (2-opt plus Or-opt) gives much better local optima than either alone.
Practice
- (Beginner) Cities $A = (0, 0)$, $B = (1, 0)$, $C = (1, 1)$ and $D = (0, 1)$ are visited in the order A, C, B, D. Compute the tour length, identify the crossing edges, apply the 2-opt move that removes them, and compute the new length.
- (Intermediate) Add don't-look bits to
two_opt(..., neighbor_lists=True): keep a queue of "active" cities, and reactivate the four endpoints of every applied move. Measure the number of evaluations for 80 and 300 random cities.- (Interview) "Why does 2-opt need changes for an asymmetric TSP, such as routing with one-way streets?" Hint: a 2-opt move reverses a segment. Ask what that does to the cost of the segment when $d(i, j) \ne d(j, i)$, and why Or-opt and orientation-preserving 3-opt moves are preferred.
Variable Neighborhood Search
Why it matters. A local optimum for one neighborhood is often not a local optimum for another, and the global optimum is a local optimum for all of them. Variable Neighborhood Search (VNS), introduced by Mladenović and Hansen in 1997, turns these observations into a simple and effective metaheuristic. It has very few parameters and is easy to add on top of an existing local search.
Theory. VNS rests on three observations:
- A local optimum with respect to one neighborhood need not be one with respect to another.
- A global optimum is a local optimum with respect to every neighborhood.
- For many problems, local optima for different neighborhoods are close to each other, so good solutions are found near other good solutions.
Its variants combine a family of neighborhoods $N_1, N_2, \dots, N_{k_{\max}}$ in different ways:
- Variable Neighborhood Descent (VND): run local search in $N_1$; when stuck, try $N_2$; after any improvement, return to $N_1$. The result is a local optimum for all the neighborhoods.
two_opt_or_opt()from the previous topic is a VND with two neighborhoods. - Basic VNS: shake the current best solution by jumping to a random point in $N_k$, then apply local search. If the result is better, move there and reset $k = 1$. Otherwise try the next, larger neighborhood $k + 1$.
- General VNS: basic VNS with VND as the local search.
x = local_search(initial solution); k = 1
while budget left:
x1 = random point in N_k(x) # shaking: bigger k, bigger jump
x2 = local_search(x1)
if f(x2) < f(x): x, k = x2, 1 # success: recenter, small jumps again
else: k = k % k_max + 1 # failure: jump farther next time
The neighborhoods for shaking are typically nested by distance: $N_k$ = "solutions $k$ random moves away". For the TSP, a good shaking move is the double bridge, which cuts the tour into four segments A B C D and reconnects them as A C B D. This 4-opt move cannot be undone easily by 2-opt or 3-opt, so the local search does not just fall back into the same local optimum. Iterated local search (ILS) is the same idea with a fixed perturbation strength: kick, descend, and accept or reject.
Example: VNS against multi-start. Both methods get 60 local searches (2-opt plus Or-opt). Multi-start begins each search from a new random tour. VNS begins each one from a shaken copy of its best tour, with $k = 1$ to 3 double bridges.
import numpy as np
# Reuses D, n, OPT, length() and nearest_neighbor() from the Nearest Neighbor and Insertion topic,
# and two_opt_or_opt() from the 2-Opt and k-Opt Moves topic.
def shake(t, k, r):
"""A random point of N_k: k random double-bridge moves (segments A B C D -> A C B D)."""
t = list(t)
for _ in range(k):
i, j, l = sorted(r.choice(np.arange(1, n), size=3, replace=False))
t = t[:i] + t[j:l] + t[i:j] + t[l:]
return t
def vns(t, k_max, budget, seed):
r = np.random.default_rng(seed)
best, searches, k, log = two_opt_or_opt(t), 1, 1, []
while searches < budget:
cand = two_opt_or_opt(shake(best, k, r)) # shake in N_k, then descend with VND
searches += 1
if length(cand) < length(best) - 1e-9:
best, k = cand, 1 # improvement: move there, restart with N_1
log.append(f"{length(best) / OPT - 1:.1%}@{searches}")
else:
k = k % k_max + 1 # no luck: try a larger neighborhood next
return best, log
BUDGET = 60 # local searches allowed per method
r = np.random.default_rng(1)
runs = [length(two_opt_or_opt(list(r.permutation(n)))) for _ in range(BUDGET)]
print(f"multi-start: gaps of the first 10 local optima "
f"{' '.join(f'{L / OPT - 1:.1%}' for L in runs[:10])}")
print(f"multi-start: best gap {min(runs) / OPT - 1:.1%}, first reached at search {np.argmin(runs) + 1}")
for seed in range(3):
best, log = vns(nearest_neighbor(0), k_max=3, budget=BUDGET, seed=seed)
print(f"VNS seed {seed}: gap {length(best) / OPT - 1:.1%} improvements: {' '.join(log)}")
# multi-start: gaps of the first 10 local optima 0.8% 1.2% 4.1% 2.3% 5.0% 3.6% 1.5% 2.5% 5.3% 3.8%
# multi-start: best gap 0.0%, first reached at search 20
# VNS seed 0: gap 0.0% improvements: 2.3%@3 1.9%@6 0.7%@7 0.3%@10 0.0%@17
# VNS seed 1: gap 0.0% improvements: 3.0%@4 1.2%@5 0.0%@6
# VNS seed 2: gap 0.0% improvements: 3.0%@2 1.9%@3 1.5%@10 0.3%@13 0.0%@23
Local optima from random starts vary between 0.8% and 5.3% above the optimum, and multi-start first reaches the optimum at search 20. VNS reaches it after 6, 17 and 23 searches depending on the seed, each improvement coming from a nearby jump. On this small instance the two are close, partly because 2-opt plus Or-opt is already strong. On large instances the difference is dramatic. Local optima from random starts concentrate around the same average gap, so restarts rarely find a better one, while VNS and ILS keep the good structure of the best tour and re-optimize only where it was perturbed. With neighbor lists and don't-look bits, a search after a kick costs a tiny fraction of a search from scratch. Here it was already roughly 15-40% cheaper.
Key takeaways
- VNS changes neighborhoods systematically: deterministically in VND, and by random shaking of increasing size in basic VNS.
- After an improvement it returns to the smallest neighborhood. After a failure it jumps farther.
- Shaking must be big enough to leave the current basin but small enough to keep good structure. The double bridge is the classic TSP kick.
- VNS and iterated local search beat random restarts on large instances because good local optima tend to be close to each other.
Practice
- (Beginner) In the VNS log, why does $k$ return to 1 after each improvement? What does VNS turn into if $k_{\max} = 1$?
- (Intermediate) Turn the code into iterated local search: always one double bridge, and accept the new local optimum if it is at most 1% worse than the current one (not only if it is better). Compare the number of searches needed to reach the optimum over 5 seeds.
- (Interview) "Our local search keeps finding the same few good solutions. Should we add more random restarts or switch to VNS?" Hint: consider how different the local optima are from each other, how expensive a restart is compared with a perturbation, and the size of the instances.
<a id="112-metaheuristics"></a>
11.2 Metaheuristics
A metaheuristic is a general strategy that guides simpler heuristics (construction and local search) beyond the first local optimum. Every metaheuristic balances two forces. Intensification searches thoroughly near good solutions, and diversification moves the search to unexplored regions. Trajectory methods (simulated annealing, tabu search, GRASP, large neighborhood search) improve one solution at a time. Population methods (genetic algorithms, ant colonies, particle swarms) evolve many solutions together. All the TSP examples use the running instance of Section 11.1, so their gaps can be compared directly.
| Method | Escapes local optima by | Memory | Typical strength |
|---|---|---|---|
| Simulated annealing | accepting some worse moves at random | none | simple, any neighborhood |
| Tabu search | taking the best non-tabu move, even uphill | tabu list, long-term statistics | scheduling, routing, coloring |
| Genetic algorithm | recombining solutions | population | black-box objectives, many objectives |
| Ant colony | sampling from a learned model | pheromone matrix | constructive problems, dynamic data |
| Particle swarm | particles attracted by good points | personal and swarm bests | continuous nonconvex problems |
| GRASP | randomized greedy restarts | none (or an elite pool) | easy parallelism |
| Large neighborhood search | destroying and rebuilding large parts | operator weights (ALNS) | routing and scheduling with side constraints |
Simulated Annealing
Why it matters. Simulated annealing (SA) is probably the easiest way to turn hill climbing into a method that escapes local optima. It needs only a neighborhood, a way to evaluate moves, and a temperature schedule. Since its introduction by Kirkpatrick, Gelatt and Vecchi (1983) and Černý (1985), it has been used for chip placement, timetabling, scheduling and many problems where better-tailored methods were not worth building.
Theory. In metallurgy, annealing heats a metal and cools it slowly, so the atoms settle into a low-energy crystal structure. Quenching (fast cooling) freezes defects instead. SA mimics this. At each step it picks a random neighbor $y$ of the current solution $x$ and computes $\Delta = f(y) - f(x)$. It then applies the Metropolis rule:
$$ P(\text{accept } y) = \begin{cases} 1 & \text{if } \Delta \le 0, \ e^{-\Delta / T} & \text{if } \Delta > 0. \end{cases} $$
At high temperature $T$, almost every move is accepted and the search wanders freely. As $T$ falls, uphill moves become rare, and at $T \to 0$ SA becomes randomized hill climbing. Small uphill steps are accepted much more often than large ones:
| $\Delta$ | $T = 41$ | $T = 5$ | $T = 1$ |
|---|---|---|---|
| 1 | 0.976 | 0.819 | 0.368 |
| 10 | 0.784 | 0.135 | 0.00005 |
Why it works. At a fixed temperature, SA is a Markov chain (Section 9.2) whose stationary distribution is the Boltzmann distribution $\pi_T(x) \propto e^{-f(x)/T}$. As $T \to 0$, this distribution concentrates on the global optima. Hajek (1988) proved that a logarithmic schedule $T_k = c / \log(1 + k)$, with $c$ at least the depth of the deepest non-global local minimum, reaches a global optimum with probability 1. That schedule is far too slow for practice, so SA in practice is a heuristic.
Practical schedule.
- Initial temperature. Choose $T_0$ so that a given share of uphill moves is accepted at the start, for example 50%: $T_0 = -\bar{\Delta}^+ / \ln 0.5$, where $\bar{\Delta}^+$ is the mean of sampled uphill deltas.
- Cooling. The most common schedule is geometric cooling $T \leftarrow \alpha T$ with $\alpha \in [0.8, 0.99]$ and $L$ moves per temperature, often proportional to the neighborhood size.
- Stopping. Stop when the acceptance rate has been near zero for a while, or when the time budget ends. Always return the best solution seen, not the last one.
- Variants. Reheating after stagnation, adaptive schedules that keep a target acceptance rate, and simpler threshold rules (accept if $\Delta < $ threshold, or "record-to-record travel").
Example: SA with 2-opt moves. The trace shows one run. The comparison then varies the cooling rate with 2,000 moves per temperature, from the same random start and over 5 seeds.
import numpy as np
# Reuses D, n, OPT and length() from the Nearest Neighbor and Insertion topic.
def simulated_annealing(t, alpha, steps_per_T, seed, verbose=False):
r = np.random.default_rng(seed)
t = list(t); cur = length(t); best, best_t = cur, t[:]
def random_move(): # random 2-opt move: edges i and j, i < j
i, j = sorted(r.choice(n, size=2, replace=False))
a, b, c, d = t[i], t[i + 1], t[j], t[(j + 1) % n]
return i, j, D[a, c] + D[b, d] - D[a, b] - D[c, d]
# Initial temperature: accept an average uphill move with probability 0.5
ups = [dl for dl in (random_move()[2] for _ in range(1000)) if dl > 0]
T = -np.mean(ups) / np.log(0.5)
stage, moves = 0, 0
while T > 0.1: # frozen: almost no uphill move is accepted
accepted = 0
moves += steps_per_T
for _ in range(steps_per_T):
i, j, delta = random_move()
if j - i < 2 or (i == 0 and j == n - 1): continue # degenerate move
if delta < 0 or r.random() < np.exp(-delta / T): # Metropolis rule
t[i + 1:j + 1] = t[i + 1:j + 1][::-1]
cur += delta; accepted += 1
if cur < best - 1e-9: best, best_t = cur, t[:]
if verbose and stage % 12 == 0:
print(f" T = {T:7.2f} accepted {accepted / steps_per_T:5.1%} "
f"current gap {cur / OPT - 1:6.1%} best gap {best / OPT - 1:6.1%}")
T *= alpha # geometric cooling
stage += 1
return best_t, moves
r0 = np.random.default_rng(0)
start = list(r0.permutation(n))
print("trace of one run (alpha = 0.95, 2,000 moves per temperature):")
simulated_annealing(start, 0.95, 2000, seed=0, verbose=True)
for alpha, steps in [(0.80, 2000), (0.95, 2000), (0.99, 2000)]:
runs = [simulated_annealing(start, alpha, steps, seed=s) for s in range(5)]
gaps = [length(t) / OPT - 1 for t, _ in runs]
print(f"alpha = {alpha:.2f}: {runs[0][1]:9,d} moves per run, mean gap {np.mean(gaps):5.1%}, "
f"best {min(gaps):5.1%}")
# trace of one run (alpha = 0.95, 2,000 moves per temperature):
# T = 41.00 accepted 60.8% current gap 365.8% best gap 285.1%
# T = 22.15 accepted 38.3% current gap 221.0% best gap 205.6%
# T = 11.97 accepted 18.9% current gap 143.2% best gap 129.6%
# T = 6.47 accepted 5.8% current gap 65.7% best gap 53.3%
# T = 3.50 accepted 2.8% current gap 31.8% best gap 24.9%
# T = 1.89 accepted 0.8% current gap 9.6% best gap 6.6%
# T = 1.02 accepted 0.4% current gap 6.6% best gap 5.7%
# T = 0.55 accepted 0.0% current gap 4.6% best gap 4.6%
# T = 0.30 accepted 0.1% current gap 4.4% best gap 4.3%
# T = 0.16 accepted 0.0% current gap 4.2% best gap 4.2%
# alpha = 0.80: 54,000 moves per run, mean gap 7.7%, best 3.3%
# alpha = 0.95: 236,000 moves per run, mean gap 3.7%, best 1.4%
# alpha = 0.99: 1,198,000 moves per run, mean gap 2.3%, best 0.3%
At $T = 41$, about 60% of moves are accepted and the tour is still almost random. Most of the progress happens between $T \approx 12$ and $T \approx 2$, where the acceptance rate falls from about 20% to 1%. In this range the search can still leave bad configurations but increasingly prefers good ones. Below $T \approx 0.5$ the tour is frozen. Slower cooling gives better tours (mean gaps of 7.7%, 3.7% and 2.3%) at proportionally more moves. Note also that SA with plain 2-opt moves needs over a million moves to reach about 2%, while 2-opt plus Or-opt local search reaches 1-3% from most starts in a fraction of a second. On the TSP, SA is beaten by tailored local search. Its strength is generality.
Key takeaways
- SA accepts a worse neighbor with probability $e^{-\Delta/T}$, so it can leave local optima. As $T$ falls it turns into hill climbing.
- At fixed $T$ it samples the Boltzmann distribution. Logarithmic cooling guarantees the optimum in the limit, but practice uses geometric cooling.
- Calibrate $T_0$ from a target acceptance rate. Most progress happens in a middle temperature range, and slower cooling trades time for quality.
- SA is simple and general, but tailored local search usually beats it on well-studied problems.
Practice
- (Beginner) An uphill move has $\Delta = 5$. Compute its acceptance probability at $T = 20$, $5$ and $1$. At what temperature is it accepted 10% of the time?
- (Intermediate) With a fixed budget of about 240,000 moves, compare $\alpha = 0.95$ with 2,000 moves per temperature, $\alpha = 0.99$ with 400 moves per temperature, and $\alpha = 0.8$ with 9,000 moves per temperature. Which matters more for a fixed budget, the number of temperatures or the moves per temperature?
- (Interview) "Evaluating one move in our problem requires a 2-second simulation. Would you choose simulated annealing or tabu search?" Hint: count the evaluations per iteration of each method. Then think about surrogate models and parallel evaluation.
Tabu Search
Why it matters. Tabu search, introduced by Glover in 1986, is one of the most successful metaheuristics for scheduling, vehicle routing, graph coloring and assignment problems. Unlike SA, it is mostly deterministic and makes systematic use of memory. Some of the best known methods for the job shop (the TSAB algorithm of Nowicki and Smutnicki) and for vehicle routing (the unified tabu search of Cordeau, Laporte and Mercier) are tabu searches.
Theory. At every iteration, tabu search moves to the best admissible neighbor, even if it is worse than the current solution. That lets it walk out of a local optimum. The obvious danger is cycling: from a local optimum, the best move goes uphill, and the best move after that goes straight back down. Tabu search prevents this with memory:
- Short-term memory (the tabu list). Recently used attributes of moves are forbidden for a number of iterations called the tenure. For a 2-opt move that drops the edges $(a, b)$ and $(c, d)$, the attribute can be "these edges may not be re-added for $t$ iterations". Storing attributes instead of whole solutions is cheap, but it also forbids some solutions never visited before.
- Aspiration criterion. A tabu move is allowed anyway if it produces a new best solution. The tabu rule is a heuristic and must not block a clear improvement.
- Long-term memory. Frequency counts of attributes drive diversification (penalize moves that were used very often) and intensification (restart from elite solutions and explore their common features).
- Candidate lists. For large neighborhoods, only a promising subset of moves is evaluated, for example with neighbor lists.
x = best = initial solution
repeat until budget is spent:
choose the best y in N(x) such that
the move x -> y is not tabu, OR f(y) < f(best) # aspiration
make the attributes of the move tabu for `tenure` iterations
x = y # even if f(y) > f(x)
if f(x) < f(best): best = x
The tenure is the key parameter. If it is too short, the search cycles. If it is too long, too many moves are forbidden, and the search cannot reach the good regions. Common remedies are random tenures drawn from a range, and reactive tabu search (Battiti and Tecchiolli), which increases the tenure when solutions repeat and decreases it otherwise.
Example: the effect of the tenure. The code evaluates all 3,000 2-opt moves of the running instance per iteration in one vectorized NumPy computation. A dropped edge becomes tabu for tenure iterations, with aspiration. As a cycling detector, it counts the distinct tour lengths visited in the last 500 iterations.
import numpy as np
# Reuses D, n, OPT, length() and nearest_neighbor() from the Nearest Neighbor and Insertion topic.
def tabu_search(t, iters, tenure, seed=0, aspiration=True):
r = np.random.default_rng(seed)
t = np.array(t); cur = length(t); best, best_t = cur, t.copy()
tabu_until = np.zeros((n, n), dtype=int) # tabu_until[a, b]: edge a-b may not be re-added before
I, J = np.triu_indices(n, k=2) # all 2-opt moves: edges i and j, i < j - 1
keep = ~((I == 0) & (J == n - 1)); I, J = I[keep], J[keep]
visited, best_it = [], 0
for it in range(iters):
a, b, c, d = t[I], t[I + 1], t[J], t[(J + 1) % n]
delta = D[a, c] + D[b, d] - D[a, b] - D[c, d] # all moves at once
is_tabu = (tabu_until[a, c] > it) | (tabu_until[b, d] > it) # would re-add a dropped edge
allowed = ~is_tabu
if aspiration:
allowed |= cur + delta < best - 1e-9 # tabu, but leads to a new best: allow it
delta = np.where(allowed, delta, np.inf)
k = np.flatnonzero(delta == delta.min())
k = k[r.integers(len(k))] # best admissible move, even if uphill
i, j = I[k], J[k]
for x, y in ((t[i], t[i + 1]), (t[j], t[(j + 1) % n])): # dropped edges become tabu
tabu_until[x, y] = tabu_until[y, x] = it + tenure
t[i + 1:j + 1] = t[i + 1:j + 1][::-1]
cur += delta[k]
visited.append(round(cur, 6))
if cur < best - 1e-9: best, best_t, best_it = cur, t.copy(), it + 1
return best_t, best_it, visited
start = nearest_neighbor(0)
print(f"start (nearest neighbor): gap {length(start) / OPT - 1:.1%}")
for tenure in (0, 5, 15, 25, 40, 60):
best_t, best_it, visited = tabu_search(start, iters=3000, tenure=tenure)
print(f"tenure {tenure:2d}: best gap {length(best_t) / OPT - 1:5.1%} found at iteration {best_it:4d}, "
f"distinct tour lengths in the last 500 iterations: {len(set(visited[-500:])):3d}")
# start (nearest neighbor): gap 22.7%
# tenure 0: best gap 4.0% found at iteration 16, distinct tour lengths in the last 500 iterations: 2
# tenure 5: best gap 4.0% found at iteration 16, distinct tour lengths in the last 500 iterations: 10
# tenure 15: best gap 2.5% found at iteration 42, distinct tour lengths in the last 500 iterations: 78
# tenure 25: best gap 1.0% found at iteration 927, distinct tour lengths in the last 500 iterations: 500
# tenure 40: best gap 0.8% found at iteration 2642, distinct tour lengths in the last 500 iterations: 500
# tenure 60: best gap 4.0% found at iteration 16, distinct tour lengths in the last 500 iterations: 500
With tenure 0, the search reaches a 2-opt local optimum (4.0%) at iteration 16 and then oscillates between two tours forever: the best move goes uphill and the next one undoes it. Tenure 5 produces a longer cycle of 10 tours, and tenure 15 a longer cycle (about 78 distinct tours) that still finds 2.5% on the way. With tenures 25 and 40 the search keeps visiting new tours and reaches 1.0% and 0.8%. Tenure 60 is too restrictive: about 120 recently dropped edges are forbidden at any time. In an 80-city tour, each city's good edges go to its few nearest neighbors, so the search is pushed into poor regions and never improves on its first local optimum.
Key takeaways
- Tabu search always moves to the best admissible neighbor, even uphill, and uses short-term memory to avoid cycling.
- Tabu attributes are cheap to store. The aspiration criterion overrides them when a move gives a new best solution.
- The tenure must be long enough to prevent cycling and short enough to keep good moves available. Random or reactive tenures reduce the tuning effort.
- Long-term memory adds intensification and diversification. Candidate lists keep large neighborhoods tractable.
Practice
- (Beginner) Explain, using the output, why tenure 0 cycles between exactly two tours once it reaches a local optimum.
- (Intermediate) Give each tabu entry a random tenure drawn uniformly from $[15, 45]$. Compare it with fixed tenures 25 and 40 over 5 different start tours (nearest neighbor from cities 0 to 4).
- (Interview) "How would you adapt tabu search to a vehicle routing problem with time windows, where most moves are infeasible?" Hint: consider accepting infeasible solutions with penalty weights that adapt during the search, attributes such as "customer $i$ in route $r$", and restricting moves to short edges (granular neighborhoods).
Genetic Algorithms
Why it matters. Genetic algorithms (GAs), introduced by Holland in the 1970s, evolve a population of solutions by selection, recombination and mutation, inspired by natural evolution. They need little knowledge of the problem, parallelize naturally, and handle several objectives at once (NSGA-II is the standard method for Pareto fronts, Section 7.1). Their combination with local search, the memetic algorithm, is among the best methods known for vehicle routing (the hybrid genetic search of Vidal and coauthors).
Theory. A GA needs a representation (chromosome) of solutions and five operators:
- Fitness: the objective, possibly with penalties for infeasibility.
- Selection of parents, biased towards fitness: tournament (the best of $k$ random individuals), roulette (probability proportional to fitness) or rank-based.
- Crossover (recombination): build a child from two parents, ideally passing on good features of both. For bit strings, use one-point or uniform crossover. For permutations, use special operators that keep the child a valid permutation: order crossover (OX), partially mapped crossover (PMX), and edge recombination or edge assembly crossover (EAX) for the TSP.
- Mutation: a small random change (flip a bit, swap or invert a segment) to keep diversity.
- Replacement: generational (children replace parents) or steady-state (one child at a time), usually with elitism (the best solution always survives).
Order crossover (OX): keep a slice of parent 1, fill the rest in parent 2's order
parent 1: 1 2 | 3 4 5 | 6 7 8 slice positions 3-5: 3 4 5
parent 2: 8 6 | 4 2 7 | 5 3 1 parent 2 read from after the slice: 5 3 1 8 6 4 2 7
skip 3, 4, 5: 1 8 6 2 7
child: 2 7 | 3 4 5 | 1 8 6 fill positions 6, 7, 8, then 1, 2
Encodings and constraints. If a direct representation produces many infeasible children, use a decoder: the chromosome is a priority list or a vector of random keys, and a constructive heuristic turns it into a feasible solution. The serial schedule generation scheme for project scheduling (Section 10.2) is the classic decoder: a GA searches over activity lists, and the scheme turns each list into a feasible schedule.
The main risk is premature convergence. The population becomes nearly identical before it finds good solutions, and crossover of identical parents produces nothing new. Diversity management (penalizing similar individuals, restarts, islands) is the remedy.
Example: pure GA versus memetic GA. The pure GA uses OX, inversion mutation, tournament selection and elitism, with 100 individuals for 500 generations, so 50,000 children. The memetic GA applies 2-opt with neighbor lists to every individual, with 20 individuals for 5 generations.
import numpy as np
# Reuses D, n, OPT and length() from the Nearest Neighbor and Insertion topic,
# and two_opt() from the 2-Opt and k-Opt Moves topic.
def order_crossover(p1, p2, r):
"""OX: copy a random slice of parent 1, fill the rest in parent 2's order."""
i, j = sorted(r.choice(n, size=2, replace=False))
child = [None] * n
child[i:j] = p1[i:j]
used = set(p1[i:j])
fill = [c for c in p2[j:] + p2[:j] if c not in used] # parent 2's order, starting after the slice
for k, c in zip(list(range(j, n)) + list(range(i)), fill):
child[k] = c
return child
def mutate(t, r, rate):
t = list(t)
if r.random() < rate: # inversion mutation: reverse a random segment
i, j = sorted(r.choice(n, size=2, replace=False))
t[i:j] = t[i:j][::-1]
return t
def genetic_algorithm(pop_size, generations, seed, memetic=False, rate=0.3):
r = np.random.default_rng(seed)
improve = (lambda t: two_opt(t, neighbor_lists=True)[0]) if memetic else (lambda t: t)
pop = [improve(list(r.permutation(n))) for _ in range(pop_size)]
fit = [length(t) for t in pop]
history = [min(fit) / OPT - 1] # generation 0: the initial population
for g in range(generations):
def tournament(): # best of 3 random individuals
k = r.choice(pop_size, size=3, replace=False)
return pop[min(k, key=lambda i: fit[i])]
children = [improve(mutate(order_crossover(tournament(), tournament(), r), r, rate))
for _ in range(pop_size)]
elite = int(np.argmin(fit)) # elitism: the best individual always survives
pop = [pop[elite]] + children[:-1]
fit = [length(t) for t in pop]
history.append(min(fit) / OPT - 1)
return min(fit) / OPT - 1, history
_, h = genetic_algorithm(pop_size=100, generations=500, seed=0)
print("pure GA, 100 x 500: best gap at generation 0, 100, ..., 500:",
" ".join(f"{x:.1%}" for x in h[::100]))
_, h = genetic_algorithm(pop_size=20, generations=5, seed=0, memetic=True)
print("memetic GA, 20 x 5 : best gap at generation 0, 1, ..., 5: ",
" ".join(f"{x:.1%}" for x in h))
# pure GA, 100 x 500: best gap at generation 0, 100, ..., 500: 401.3% 172.0% 125.1% 111.8% 99.7% 90.2%
# memetic GA, 20 x 5 : best gap at generation 0, 1, ..., 5: 4.0% 2.8% 0.8% 0.8% 0.8% 0.0%
The pure GA improves steadily but slowly. After 50,000 children it is still 90% above the optimum, which is worse than nearest neighbor. OX preserves the relative order of the cities, but in a TSP what matters is which cities are adjacent. Most children inherit few good edges, and blind mutation rarely repairs them. The memetic GA starts from 2-opt local optima (4.0%) and reaches the optimum in the fifth generation, after only 120 local searches. On combinatorial problems, the GA's role is to recombine good local optima, while local search does the fine work. Pure GAs are rarely competitive there.
Key takeaways
- A GA evolves a population with selection, crossover, mutation and replacement. Elitism keeps the best solution.
- Crossover must fit the representation. For permutations use OX or PMX, and for the TSP edge-based operators such as EAX.
- Decoders (priority lists, random keys) handle constraints by turning any chromosome into a feasible solution.
- Pure GAs are weak on combinatorial problems. Memetic algorithms (GA plus local search) are among the strongest methods.
Practice
- (Beginner) Apply OX to parents $(1, 2, 3, 4, 5, 6)$ and $(3, 6, 1, 5, 2, 4)$, keeping positions 2 to 4 of parent 1. Check that the child is a permutation.
- (Intermediate) Replace OX in the pure GA by an edge-preserving crossover: build the child from the union of the parents' edges, always moving to a neighbor (in either parent) that has the fewest remaining neighbors. How much does the pure GA improve?
- (Interview) "When would you pick a genetic algorithm over tabu search?" Hint: consider objectives that come from a black-box simulator, several objectives at once, cheap parallel evaluation, and whether a natural recombination of two solutions exists.
Ant Colony Optimization
Why it matters. Ant colony optimization (ACO), introduced by Dorigo in 1992, builds solutions with many simple agents that learn from each other through a shared memory. It works best on problems that are naturally solved by construction, such as routing, sequencing and assignment. Its pheromone memory adapts gradually, which also makes it useful on problems whose data changes over time, such as routing in communication networks.
Theory. Real ants lay pheromone as they walk and prefer trails with more pheromone. If two paths lead to food, ants on the shorter path return sooner, so pheromone accumulates faster there and more ants choose it. The colony "finds" the shortest path without any ant knowing it.
ACO for the TSP copies this. Each artificial ant builds a tour city by city. From city $i$ it moves to an unvisited city $j$ with probability
$$ p_{ij} = \frac{\tau_{ij}^{\alpha}, \eta_{ij}^{\beta}}{\sum_{l \text{ unvisited}} \tau_{il}^{\alpha}, \eta_{il}^{\beta}}, $$
where $\tau_{ij}$ is the pheromone (learned desirability of edge $ij$), $\eta_{ij} = 1 / d_{ij}$ is the heuristic information, and $\alpha$ and $\beta$ weigh the two. After all ants finish:
- Evaporation: $\tau_{ij} \leftarrow (1 - \rho), \tau_{ij}$ forgets old information, at rate $\rho$.
- Deposit: edges of good tours receive $\Delta\tau = 1 / L$, where $L$ is the tour length.
The main variants differ in who deposits and how:
| Variant | Deposit | Special feature |
|---|---|---|
| Ant System (1992) | every ant | the original version, weak on large instances |
| Ant Colony System | best-so-far ant | greedy choice with probability $q_0$, local pheromone update while building |
| MAX-MIN Ant System (Stützle and Hoos, 2000) | iteration-best or best-so-far ant | pheromone kept in $[\tau_{\min}, \tau_{\max}]$ to avoid stagnation, initialized at $\tau_{\max}$ |
The pheromone matrix is a probability model of good solutions that is learned from samples. In this sense ACO belongs to the model-based search methods, together with the cross-entropy method and estimation-of-distribution algorithms. In practice, the best ACO methods apply local search to the ants' tours before depositing.
Example: MAX-MIN Ant System and the role of $\beta$. 20 ants per iteration, $\alpha = 1$, $\rho = 0.1$, and no local search.
import numpy as np
# Reuses D, n, OPT, length() and nearest_neighbor() from the Nearest Neighbor and Insertion topic.
def max_min_ant_system(ants, iters, alpha, beta, rho, seed):
r = np.random.default_rng(seed)
eta = 1.0 / (D + np.eye(n)) # heuristic desirability: inverse distance
tau_max = 1.0 / (rho * length(nearest_neighbor(0)))
tau_min = tau_max / (2 * n)
tau = np.full((n, n), tau_max) # MMAS starts at the upper bound: maximal exploration
best_t, best_L, log = None, np.inf, []
for it in range(iters):
for _ in range(ants): # each ant builds a tour city by city
t = [int(r.integers(n))]
free = np.ones(n, dtype=bool); free[t[0]] = False
for _ in range(n - 1):
w = (tau[t[-1]] ** alpha) * (eta[t[-1]] ** beta) * free # attractiveness of each free city
nxt = int(r.choice(n, p=w / w.sum()))
t.append(nxt); free[nxt] = False
L = length(t)
if L < best_L: best_t, best_L = t, L
tau *= 1 - rho # evaporation forgets old information
for k in range(n): # only the best-so-far tour deposits pheromone
a, b = best_t[k - 1], best_t[k]
tau[a, b] += 1.0 / best_L; tau[b, a] = tau[a, b]
tau_max = 1.0 / (rho * best_L)
np.clip(tau, tau_min, tau_max, out=tau) # the MAX-MIN bounds prevent stagnation
if (it + 1) % 25 == 0: log.append(best_L / OPT - 1)
return best_t, log
for beta in (0, 2, 5):
_, log = max_min_ant_system(ants=20, iters=150, alpha=1, beta=beta, rho=0.1, seed=0)
print(f"beta = {beta}: best gap after 25, 50, ..., 150 iterations: "
+ " ".join(f"{g:5.1%}" for g in log))
# beta = 0: best gap after 25, 50, ..., 150 iterations: 368.9% 352.0% 352.0% 352.0% 352.0% 345.0%
# beta = 2: best gap after 25, 50, ..., 150 iterations: 73.2% 13.6% 8.7% 4.7% 3.6% 2.9%
# beta = 5: best gap after 25, 50, ..., 150 iterations: 13.7% 5.7% 3.8% 1.0% 1.0% 1.0%
Without heuristic information ($\beta = 0$), ants choose almost at random. The first poor tours deposit pheromone, and the colony reinforces them: after 3,000 tours it is still 345% above the optimum. With $\beta = 2$, the distance guides the construction and the pheromone learns which short edges belong together: 2.9% after 150 iterations. With $\beta = 5$, the colony reaches 1.0% after 100 iterations, using only 2,000 constructed tours. That beats the average SA run with over a million moves. A strong heuristic bias makes the learned pheromone useful much sooner.
Key takeaways
- Ants build solutions with probabilities proportional to $\tau^{\alpha} \eta^{\beta}$, combining learned pheromone with greedy heuristic information.
- Evaporation forgets, and deposits on good solutions reinforce. Only the best ants should deposit (MMAS, ACS).
- Pheromone bounds (MMAS) prevent premature stagnation. The heuristic weight $\beta$ is critical.
- ACO is a learned sampling model for constructive problems. Strong implementations add local search.
Practice
- (Beginner) An ant at city $i$ can go to three cities with pheromone $\tau = (1, 2, 1)$ and distances $d = (10, 20, 5)$. Compute the transition probabilities for $\alpha = 1$, $\beta = 2$, and again for $\beta = 0$.
- (Intermediate) Add 2-opt with neighbor lists (Section 11.1) to the iteration-best tour before the pheromone update. How many iterations does MMAS with $\beta = 2$ then need to reach 1%?
- (Interview) "ACO is said to be good for routing packets in networks where traffic changes constantly. Why?" Hint: think about evaporation as a forgetting mechanism, about continuous small updates, and about many agents sampling routes in parallel.
Particle Swarm Optimization
Why it matters. Many optimization problems have continuous variables but an objective that is nonconvex, noisy or available only as a black box: calibrating a simulation model, tuning controller gains, or fitting a model with many local optima. Gradient methods (Chapter 5) find only a local optimum, and the gradient may not even exist. Particle swarm optimization (PSO), proposed by Kennedy and Eberhart in 1995, is a simple population method for these problems, inspired by birds flocking.
Theory. A swarm of particles moves through $\mathbb{R}^d$. Particle $i$ has a position $x_i$ (a candidate solution), a velocity $v_i$, and remembers its personal best position $p_i$. It also knows the best position $g_i$ found by its neighborhood in the swarm. In each iteration:
$$ v_i \leftarrow w, v_i + c_1 r_1 \odot (p_i - x_i) + c_2 r_2 \odot (g_i - x_i), \qquad x_i \leftarrow x_i + v_i, $$
where $r_1, r_2$ are vectors of uniform random numbers in $[0, 1]$ and $\odot$ multiplies element by element. The three terms are:
- inertia $w, v_i$: keep moving in the same direction,
- cognitive term: return towards your own best memory,
- social term: move towards the best point your neighbors found.
The standard constriction parameters $w = 0.729$ and $c_1 = c_2 = 1.494$ (Clerc and Kennedy, 2002) keep the swarm stable, so that it converges instead of exploding.
The topology decides who counts as a neighbor:
- Global best (gbest): every particle follows the best of the whole swarm. Information spreads fast, but the swarm may converge prematurely to one basin.
- Ring (lbest): each particle sees only itself and two neighbors. Good positions spread slowly, so several basins are explored in parallel, which makes it more robust on multimodal functions.
PSO handles bounds by clipping and constraints by penalties. Discrete variants exist (binary PSO, random-key encodings for permutations), but PSO is mainly a continuous method.
Example: the Rastrigin function. The function
$$ f(x) = 10 d + \sum_{i=1}^{d} \left( x_i^2 - 10 \cos(2 \pi x_i) \right), \qquad x \in [-5.12, 5.12]^d, $$
has its global minimum $f = 0$ at the origin and a local minimum near every integer point, more than $10^{10}$ of them for $d = 10$. Every method gets 40,000 function evaluations, 10 runs each.
import numpy as np
from scipy.optimize import minimize
def rastrigin(X): # global minimum 0 at x = 0, about 10^d local minima
return 10 * X.shape[-1] + np.sum(X ** 2 - 10 * np.cos(2 * np.pi * X), axis=-1)
DIM, LO, HI, EVALS = 10, -5.12, 5.12, 40_000 # same evaluation budget for every method
def pso(seed, topology, swarm=40, w=0.729, c1=1.49445, c2=1.49445):
r = np.random.default_rng(seed)
X = r.uniform(LO, HI, size=(swarm, DIM)) # positions
V = r.uniform(-1, 1, size=(swarm, DIM)) # velocities
P, fP = X.copy(), rastrigin(X) # personal bests
for _ in range(EVALS // swarm - 1):
if topology == "global": # every particle follows the swarm's best
G = np.repeat(P[[fP.argmin()]], swarm, axis=0)
else: # ring: best of itself and its two neighbors
idx = np.stack([np.arange(swarm) - 1, np.arange(swarm), (np.arange(swarm) + 1) % swarm])
G = P[idx[fP[idx].argmin(axis=0), np.arange(swarm)]]
r1, r2 = r.random((2, swarm, DIM))
V = w * V + c1 * r1 * (P - X) + c2 * r2 * (G - X) # inertia + memory + social pull
X = np.clip(X + V, LO, HI)
f = rastrigin(X)
better = f < fP
P[better], fP[better] = X[better], f[better]
return fP.min()
def random_search(seed):
return rastrigin(np.random.default_rng(seed).uniform(LO, HI, size=(EVALS, DIM))).min()
def multistart_lbfgs(seed):
r, best, used = np.random.default_rng(seed), np.inf, 0
while used < EVALS: # restart from random points until the budget is spent
res = minimize(lambda x: rastrigin(x), r.uniform(LO, HI, DIM), method="L-BFGS-B",
bounds=[(LO, HI)] * DIM)
used += res.nfev; best = min(best, res.fun)
return best
for name, method in [("random search", random_search), ("multi-start L-BFGS", multistart_lbfgs),
("PSO, global best", lambda s: pso(s, "global")),
("PSO, ring", lambda s: pso(s, "ring"))]:
vals = [method(s) for s in range(10)]
print(f"{name:19s} best f over 10 runs: median {np.median(vals):6.2f}, "
f"min {min(vals):6.2f}, max {max(vals):6.2f}")
# random search best f over 10 runs: median 58.93, min 51.06, max 74.04
# multi-start L-BFGS best f over 10 runs: median 21.39, min 11.94, max 28.85
# PSO, global best best f over 10 runs: median 8.46, min 3.98, max 13.93
# PSO, ring best f over 10 runs: median 5.97, min 3.98, max 8.96
Random search gets nowhere in 10 dimensions. Multi-start L-BFGS finds a local minimum quickly from each start, but each start lands in a random basin, and estimating gradients by finite differences costs 11 evaluations each. PSO uses the swarm's shared information to concentrate on good regions. The ring topology is better and more consistent (median 5.97, worst 8.96) than the global-best topology, which converges faster but sometimes to a worse basin. A value of 3.98 means 4 coordinates sit near $\pm 1$, a local minimum next to the global one, and no method found the global optimum. On multimodal landscapes, PSO is a good tool for a good solution, not for a certified optimum.
Key takeaways
- Each particle combines inertia, attraction to its own best position and attraction to its neighbors' best position.
- Constriction parameters ($w = 0.729$, $c_1 = c_2 = 1.494$) keep the swarm stable. The topology trades speed (gbest) for robustness (ring).
- PSO suits continuous, nonconvex, black-box problems. Use gradient methods for smooth convex problems and LP or MIP solvers for linear models.
- Like all metaheuristics, PSO gives no proof of optimality.
Practice
- (Beginner) A particle in 2D is at $x = (1, 2)$ with velocity $v = (0.5, -0.5)$, personal best $p = (0, 2)$ and neighborhood best $g = (2, 0)$. With $w = 0.7$, $c_1 = c_2 = 1.5$, and random numbers $r_1 = (0.2, 0.6)$, $r_2 = (0.5, 0.1)$, compute the new velocity and position.
- (Intermediate) Run the comparison on the Rosenbrock function, which is unimodal but has a narrow curved valley. Which method wins now, and why?
- (Interview) "Would you use PSO to solve a linear program with 10,000 variables?" Hint: compare what an LP solver guarantees (optimality, duals, speed) with what PSO offers. Then name problems where PSO is the right tool.
GRASP
Why it matters. The Greedy Randomized Adaptive Search Procedure (GRASP), introduced by Feo and Resende, is a disciplined form of multi-start. Every iteration builds a new solution with a randomized greedy heuristic and improves it with local search. It is easy to implement on top of an existing greedy algorithm and local search, its iterations are independent and run in parallel, and it has been applied to hundreds of problems, from set covering to network design.
Theory. Each GRASP iteration has two phases, and the best solution over all iterations is returned:
- Greedy randomized construction. At each step, score the candidates with the greedy function, put the best ones in a restricted candidate list (RCL), and pick one of them at random.
- Local search from the constructed solution.
The RCL can be built in two ways:
- Cardinality-based: the $k$ best candidates. $k = 1$ is pure greedy.
- Value-based: every candidate with $c \le c_{\min} + \alpha,(c_{\max} - c_{\min})$. $\alpha = 0$ is pure greedy and $\alpha = 1$ is pure random.
The parameter controls a trade-off. Pure greedy builds good starting points, but always the same ones, so local search always ends in the same local optimum. Pure random gives maximum diversity, but poor starts that need long local searches and often end in worse local optima. Good values are in between, and reactive GRASP adapts $\alpha$ by sampling values with probabilities proportional to their past success.
GRASP iterations have no memory of each other. Two standard extensions add some:
- Path relinking: walk from one elite solution towards another, one move at a time, and apply local search to the best intermediate solution. Good solutions often share structure, and the path between them explores it.
- Elite pools and bias functions: keep the best solutions, and bias the construction towards their features.
Example: the size of the RCL. The randomized nearest neighbor always starts at city 0 and picks at random among the $k$ nearest unvisited cities. Each setting runs 30 iterations with 2-opt plus Or-opt.
import numpy as np
# Reuses D, n, OPT and length() from the Nearest Neighbor and Insertion topic,
# and two_opt_or_opt() from the 2-Opt and k-Opt Moves topic.
def randomized_nearest_neighbor(k, r):
"""Greedy randomized construction: pick at random among the k nearest unvisited cities."""
tour, left = [0], set(range(1, n))
while left:
cand = np.array(sorted(left))
rcl = cand[np.argsort(D[tour[-1], cand])[:k]] # restricted candidate list
nxt = int(r.choice(rcl))
tour.append(nxt); left.remove(nxt)
return tour
def grasp(k, iters, seed):
r = np.random.default_rng(seed)
built, improved = [], []
for _ in range(iters):
t = randomized_nearest_neighbor(k, r) # phase 1: greedy randomized construction
built.append(length(t))
improved.append(length(two_opt_or_opt(t))) # phase 2: local search
return np.array(built), np.array(improved)
print(" RCL size k construction after LS (mean) after LS (best) local optima hits of OPT")
for k in (1, 2, 3, 5, 10, n):
built, improved = grasp(k, iters=30, seed=0)
print(f"{k:11d} {built.mean() / OPT - 1:14.1%} {improved.mean() / OPT - 1:17.1%} "
f"{improved.min() / OPT - 1:17.1%} {len(set(np.round(improved, 4))):14d} "
f"{np.sum(improved < OPT + 1e-3):13d}")
# RCL size k construction after LS (mean) after LS (best) local optima hits of OPT
# 1 22.7% 3.2% 3.2% 1 0
# 2 63.3% 4.2% 0.3% 29 0
# 3 100.3% 3.6% 0.0% 30 1
# 5 149.9% 3.5% 0.0% 28 2
# 10 234.7% 3.6% 0.0% 30 1
# 80 486.5% 4.2% 0.0% 28 1
With $k = 1$ the construction is deterministic, so all 30 iterations produce the same local optimum, 3.2% above the optimum: no amount of repetition helps. With $k = 2$, the 30 iterations give 29 different local optima and the best is 0.3%. From $k = 3$ on, GRASP finds the optimum. The construction quality gets much worse with $k$ (486% for a random start), but the average after local search barely moves on this small instance. What the randomization adds is diversity. On larger instances, poor starts cost more local-search time and give worse local optima, which is why small values ($k$ = 2 to 5, or $\alpha$ = 0.1 to 0.3) are common.
Key takeaways
- GRASP repeats two phases: a greedy construction randomized by a restricted candidate list, then local search. It keeps the best result.
- The RCL size ($k$ or $\alpha$) trades the quality of starting points against diversity. Pure greedy produces one local optimum only.
- Reactive GRASP adapts $\alpha$. Path relinking and elite pools add memory between iterations.
- Iterations are independent, so GRASP parallelizes almost perfectly.
Practice
- (Beginner) Five candidates have greedy costs 4, 5, 7, 12 and 20. List the value-based RCL for $\alpha = 0$, $0.25$ and $1$.
- (Intermediate) Add path relinking to the GRASP example: for the two best tours found, apply 2-opt moves that each add an edge of the target tour, evaluate every intermediate tour, and run local search from the best one. Does it find better tours than more GRASP iterations would?
- (Interview) "GRASP is embarrassingly parallel. What do we lose by running 1,000 independent iterations on 1,000 cores?" Hint: think about what later iterations could learn from earlier ones (reactive $\alpha$, elite pools, path relinking) and how to share it cheaply.
Large Neighborhood Search
Why it matters. Real routing and scheduling problems come with many side constraints: time windows, driver breaks, vehicle types, precedences. Classic small moves such as 2-opt often turn a feasible solution into an infeasible one, and designing a move for every constraint is impractical. Large Neighborhood Search (LNS), introduced by Shaw in 1998, avoids this. It repeatedly destroys part of the solution and repairs it, and the repair step handles the constraints. Its adaptive version (ALNS) by Ropke and Pisinger (2006) is one of the most widely used methods for vehicle routing in industry.
Theory. Each iteration:
- Destroy: remove $q$ elements (customers, jobs, assignments) from the current solution. The size $q$ is typically 10-40% of the solution.
- Repair: reinsert them with a heuristic (greedy or regret insertion), or with an exact method such as constraint programming or a small MIP.
- Accept or reject the new solution: only if better, with an SA-like rule, or within a threshold of the best.
The neighborhood (all solutions reachable by destroying $q$ elements and repairing them) is exponentially large, but the repair heuristic searches it efficiently. Diversity comes from the choice of destroy operator:
| Destroy operator | Removes | Purpose |
|---|---|---|
| Random removal | $q$ random customers | diversification |
| Worst removal | customers whose removal saves the most (randomized) | fix badly placed customers |
| Related (Shaw) removal | a seed customer and the customers most similar to it (close in space or time) | let similar customers swap places |
| Route or cluster removal | a whole route or a geographic cluster | reduce the number of vehicles, rebuild regions |
Repair operators are greedy insertion (each customer at its cheapest feasible position) and regret-$k$ insertion, which first inserts the customer with the largest difference between its best and $k$-th best insertion cost (the regret idea of the Greedy Heuristics topic).
Adaptive LNS (ALNS) keeps a weight $w$ for every operator and chooses operators with probability proportional to their weights. After each iteration, the operators used get a score $\sigma$, for example 10 for a new best solution, 5 for an improvement of the current one and 1 for an accepted one. Their weights are then updated as
$$ w \leftarrow (1 - r), w + r, \sigma, $$
with a reaction factor $r$ such as 0.1. ALNS thus learns online which operators work on the instance at hand.
Example: capacitated vehicle routing. The 80 cities of the running instance become customers with demands from 1 to 9. A depot is at the center, and each vehicle carries 60 units. Each method runs 1,500 iterations with SA-like acceptance, removing 5 to 15 customers per iteration. The starting solution inserts the customers greedily in random order. The alns() function can also record every route it produces, a feature the Matheuristics topic uses.
import numpy as np
# Reuses pts from the Nearest Neighbor and Insertion topic: the 80 cities become customers.
r0 = np.random.default_rng(1)
P = np.vstack([[50, 50], pts]) # node 0 is a central depot
C = np.linalg.norm(P[:, None] - P[None, :], axis=2)
q = np.concatenate([[0], r0.integers(1, 10, size=80)]) # demands
Q = 60 # vehicle capacity
def cost(routes):
return sum(C[0, r[0]] + sum(C[a, b] for a, b in zip(r, r[1:])) + C[r[-1], 0] for r in routes if r)
def insertion_options(routes, c):
"""(extra cost, route index, position) of every feasible place for customer c."""
opts = [(C[0, c] * 2, len(routes), 0)] # a new route of its own
for k, r in enumerate(routes):
if sum(q[r]) + q[c] > Q: continue
path = [0] + r + [0]
for p in range(len(path) - 1):
opts.append((C[path[p], c] + C[c, path[p + 1]] - C[path[p], path[p + 1]], k, p))
return sorted(opts)
def place(routes, c, k, p):
if k == len(routes): routes.append([c])
else: routes[k].insert(p, c)
# ---- destroy operators: remove `m` customers, return the partial solution and the removed list
def random_removal(routes, m, r):
out = set(r.choice(np.arange(1, 81), size=m, replace=False))
return [[c for c in rt if c not in out] for rt in routes], list(out)
def worst_removal(routes, m, r):
saving = {} # how much each customer costs where it is
for rt in routes:
path = [0] + rt + [0]
for i in range(1, len(path) - 1):
saving[path[i]] = C[path[i - 1], path[i]] + C[path[i], path[i + 1]] - C[path[i - 1], path[i + 1]]
ranked = sorted(saving, key=lambda c: -saving[c])
out = set()
while len(out) < m: # randomized: favor, but do not force, the worst
out.add(ranked[int(len(ranked) * r.random() ** 4)])
return [[c for c in rt if c not in out] for rt in routes], list(out)
def related_removal(routes, m, r):
seed = int(r.integers(1, 81)) # a customer and its nearest customers
near = [c for c in np.argsort(C[seed]) if c != 0][:m]
out = set(near)
return [[c for c in rt if c not in out] for rt in routes], list(out)
# ---- repair operators
def greedy_repair(routes, removed, r):
for c in r.permutation(removed): # each customer at its cheapest feasible place
_, k, p = insertion_options(routes, c)[0]
place(routes, c, k, p)
return [rt for rt in routes if rt]
def regret_repair(routes, removed, r):
removed = list(removed)
while removed: # insert first the customer with most to lose
best = None
for c in removed:
o = insertion_options(routes, c)
regret = o[1][0] - o[0][0] if len(o) > 1 else np.inf
if best is None or regret > best[0]: best = (regret, c, o[0])
_, c, (_, k, p) = best
place(routes, c, k, p); removed.remove(c)
return [rt for rt in routes if rt]
def alns(destroys, repairs, iters, seed, adaptive, pool=None):
r = np.random.default_rng(seed)
cur = greedy_repair([], list(range(1, 81)), r) # initial solution: insert everyone greedily
best, start = cur, cost(cur)
w_d, w_r = np.ones(len(destroys)), np.ones(len(repairs))
T = 0.02 * start # SA-style acceptance, cooled linearly
for it in range(iters):
i = r.choice(len(destroys), p=w_d / w_d.sum()) if adaptive else r.integers(len(destroys))
j = r.choice(len(repairs), p=w_r / w_r.sum()) if adaptive else r.integers(len(repairs))
m = int(r.integers(5, 16)) # remove 5-15 customers
partial, removed = destroys[i]([list(rt) for rt in cur], m, r)
cand = repairs[j](partial, removed, r)
if pool is not None: # remember every route seen (Matheuristics topic)
for rt in cand:
key = frozenset(rt)
if key not in pool or cost([rt]) < cost([pool[key]]): pool[key] = rt
temp = T * (1 - it / iters)
score = 0
if cost(cand) < cost(best) - 1e-9: best, score = cand, 10 # new best solution
if cost(cand) < cost(cur) - 1e-9: score = max(score, 5) # improves the current one
if cost(cand) < cost(cur) - 1e-9 or r.random() < np.exp(-(cost(cand) - cost(cur)) / max(temp, 1e-9)):
cur, score = cand, max(score, 1)
if adaptive: # reward operators that produced good solutions
w_d[i] = 0.9 * w_d[i] + 0.1 * score
w_r[j] = 0.9 * w_r[j] + 0.1 * score
w_d, w_r = np.maximum(w_d, 0.1), np.maximum(w_r, 0.1)
return start, best, w_d / w_d.sum(), w_r / w_r.sum()
print(f"total demand {q.sum()}, capacity {Q}: at least {int(np.ceil(q.sum() / Q))} vehicles")
for name, D_ops, R_ops, adaptive in [
("LNS (random removal, greedy repair)", [random_removal], [greedy_repair], False),
("ALNS (3 destroy x 2 repair, adaptive)", [random_removal, worst_removal, related_removal],
[greedy_repair, regret_repair], True)]:
for seed in range(2):
start, best, wd, wr = alns(D_ops, R_ops, iters=1500, seed=seed, adaptive=adaptive)
print(f"{name} seed {seed}: initial {start:6.1f} -> best {cost(best):6.1f} "
f"with {len(best)} routes" + (f" destroy weights {np.round(wd, 2)}, repair {np.round(wr, 2)}"
if adaptive else ""))
# total demand 392, capacity 60: at least 7 vehicles
# LNS (random removal, greedy repair) seed 0: initial 2178.8 -> best 979.0 with 7 routes
# LNS (random removal, greedy repair) seed 1: initial 2038.9 -> best 966.9 with 7 routes
# ALNS (3 destroy x 2 repair, adaptive) seed 0: initial 2178.8 -> best 964.4 with 7 routes destroy weights [0.43 0.4 0.18], repair [0.42 0.58]
# ALNS (3 destroy x 2 repair, adaptive) seed 1: initial 2038.9 -> best 954.6 with 7 routes destroy weights [0.42 0.41 0.17], repair [0.53 0.47]
Both methods turn a poor starting solution (over 2,000) into routes of about 955-980 with 7 vehicles, the minimum the total demand allows. ALNS with three destroy and two repair operators beats plain LNS on both seeds. Its learned weights are informative: random and worst removal get about 40% each, while related removal gets less than 20% on this instance. A different instance, for example with time windows, could produce quite different weights. That is exactly why the weights are learned rather than fixed.
Key takeaways
- LNS alternates destroy and repair. Each step explores an exponentially large neighborhood with a heuristic or exact repair method.
- Destroy operators (random, worst, related, route) give diversity. Repair operators (greedy, regret, exact) handle the side constraints.
- ALNS chooses operators adaptively, with weights updated from their recent success.
- LNS is the standard approach for rich vehicle routing and scheduling problems with many constraints.
Practice
- (Beginner) Customer $c$ sits between $a$ and $b$ with $d(a, c) = 4$, $d(c, b) = 5$ and $d(a, b) = 6$. Compute its removal saving. Another customer has insertion costs 3, 4 and 9 in three routes, and a third has 2, 8 and 9. Which one does regret-2 insert first?
- (Intermediate) Add a route removal destroy operator (remove all customers of a random route) and set the capacity to 70. Does ALNS find solutions with 6 vehicles, and which operators gain weight?
- (Interview) "Our dispatchers add constraints every month (new vehicle types, driver breaks, customer time windows). Which solution method keeps up best?" Hint: compare how much code changes in a local search with specialized moves, in an LNS whose repair step checks feasibility, and in a MIP or CP model.
<a id="113-hybrid-and-learning-based-methods"></a>
11.3 Hybrid and Learning-Based Methods
The best modern solvers rarely use one technique in isolation. This subchapter covers four kinds of combination. Matheuristics put LP and MIP solvers inside heuristics. Constraint programming brings propagation-based search for scheduling and logical constraints. Machine learning extracts regularities from many similar instances. Algorithm configuration turns the choice of parameters into an optimization problem of its own. Each topic builds on the methods of Sections 11.1 and 11.2.
Matheuristics
Why it matters. Heuristics are fast but blind to bounds, while MIP solvers (Chapter 4) are rigorous but struggle with very large models. Matheuristics combine the two. A heuristic reduces the problem to a size that a MIP can handle, or a MIP solves the hard part of a heuristic exactly. Many record-breaking methods for vehicle routing, lot sizing and network design are matheuristics, and commercial MIP solvers themselves rely on MIP-based heuristics such as RINS and local branching to find good incumbents (Section 4.2).
Theory. The main families are:
| Family | Idea | Examples |
|---|---|---|
| MIP-based neighborhoods | add constraints that keep the MIP close to an incumbent $\bar{x}$ and solve the restricted MIP | local branching, RINS, proximity search, polishing |
| Decomposition over time or variables | solve a sequence of smaller MIPs, fixing or relaxing parts of the model | relax-and-fix, fix-and-optimize (lot sizing), kernel search |
| Recombination by set partitioning | a heuristic generates many good columns (routes, schedules), then a MIP picks the best combination | route-pool set partitioning, the petal method |
| Relaxation-based | turn LP or Lagrangian relaxations into feasible solutions | LP rounding and diving, Lagrangian heuristics (Section 4.2) |
| Exact repair inside LNS | the repair step of LNS solves a small MIP or CP model | MIP-based LNS |
Local branching (Fischetti and Lodi, 2003) restricts a binary MIP to the solutions that differ from the incumbent $\bar{x}$ in at most $k$ variables:
$$ \sum_{j:\ \bar{x}j = 1} (1 - x_j) + \sum{j:\ \bar{x}_j = 0} x_j \le k. $$
The left side is the Hamming distance to $\bar{x}$. For small $k$ (10 to 20), the restricted MIP is a neighborhood that the solver searches exactly, using its LP bounds to prune. RINS (relaxation induced neighborhood search) fixes every variable on which the incumbent and the LP optimum agree, and solves a sub-MIP over the rest.
Set partitioning over a pool of columns. For routing, every route $r$ produced by a heuristic is a column with cost $c_r$, covering a set of customers. The recombination MIP chooses the cheapest set of routes that serves every customer exactly once:
$$ \min \sum_{r} c_r y_r \quad \text{s.t.} \quad \sum_{r \ni i} y_r = 1 \ \ \forall \text{ customers } i, \qquad y_r \in {0, 1}. $$
The heuristic supplies diversity, and the MIP supplies optimal recombination. Different solutions often contain excellent routes that never appeared together in one solution.
ALNS iterations ----> route pool -----> set partitioning MIP -----> best combination
(every route seen) (thousands of (one row per customer, (can mix routes from
columns) one column per route) different solutions)
Example: recombining the routes of ALNS. ALNS from the previous topic records every route it produces, keeping the cheapest sequence for each set of customers. The set partitioning MIP then selects the best combination of pooled routes. Its LP relaxation is a lower bound for the pool only, not for the vehicle routing problem.
import numpy as np
from scipy.optimize import milp, LinearConstraint, Bounds
# Reuses C, cost(), alns() and the destroy/repair operators from the Large Neighborhood Search topic.
def set_partitioning(routes):
"""Choose routes so that every customer is served exactly once, at minimum total cost."""
A = np.zeros((80, len(routes)))
for k, rt in enumerate(routes):
A[np.array(rt) - 1, k] = 1 # column k covers the customers of route k
c = np.array([cost([rt]) for rt in routes])
lp = milp(c, constraints=LinearConstraint(A, 1, 1), bounds=Bounds(0, 1)) # LP relaxation
ip = milp(c, constraints=LinearConstraint(A, 1, 1), bounds=Bounds(0, 1),
integrality=np.ones(len(routes)))
chosen = [routes[k] for k in np.flatnonzero(ip.x > 0.5)]
return chosen, lp.fun
destroys = [random_removal, worst_removal, related_removal]
repairs = [greedy_repair, regret_repair]
for iters in (300, 1500):
for seed in range(3):
pool = {} # every route ALNS produced, rejected or not
_, best, _, _ = alns(destroys, repairs, iters=iters, seed=seed, adaptive=True, pool=pool)
sp, lp = set_partitioning(list(pool.values()))
print(f"{iters:4d} iterations, seed {seed}: ALNS best {cost(best):6.1f} | pool of {len(pool):4d} routes: "
f"LP bound {lp:6.1f}, set partitioning {cost(sp):6.1f} (saves {round(cost(best), 1) - round(cost(sp), 1):4.1f})")
# 300 iterations, seed 0: ALNS best 966.9 | pool of 445 routes: LP bound 966.9, set partitioning 966.9 (saves 0.0)
# 300 iterations, seed 1: ALNS best 1013.6 | pool of 474 routes: LP bound 994.2, set partitioning 994.2 (saves 19.4)
# 300 iterations, seed 2: ALNS best 1100.2 | pool of 541 routes: LP bound 1089.9, set partitioning 1100.2 (saves 0.0)
# 1500 iterations, seed 0: ALNS best 964.4 | pool of 1308 routes: LP bound 951.8, set partitioning 954.6 (saves 9.8)
# 1500 iterations, seed 1: ALNS best 954.6 | pool of 1265 routes: LP bound 954.6, set partitioning 954.6 (saves 0.0)
# 1500 iterations, seed 2: ALNS best 954.6 | pool of 1334 routes: LP bound 954.4, set partitioning 954.6 (saves 0.0)
The MIP over pools of 445 to 1,334 routes solves in a fraction of a second. Every route of the ALNS solution is in the pool, so the recombination can never be worse than ALNS, and sometimes it is clearly better. After 300 iterations it improves seed 1 from 1,013.6 to 994.2 (1.9%). After 1,500 iterations it lifts seed 0 from 964.4 to 954.6, the best value any run found. When the LP bound equals the integer value, as in several rows, the LP alone proves that the pool contains no better combination. In a real implementation this step runs periodically: ALNS continues from the recombined solution and keeps feeding the pool.
Key takeaways
- Matheuristics embed LP and MIP solvers in heuristics, or heuristics in MIP solvers.
- Local branching and RINS define neighborhoods by adding constraints to the MIP. Relax-and-fix and kernel search solve sequences of smaller MIPs.
- Set partitioning over a pool of heuristic columns recombines good parts of many solutions optimally and never loses the incumbent.
- LP bounds inside a matheuristic are valid only for the restricted problem, unless the restriction is provably harmless.
Practice
- (Beginner) Write the local branching constraint for the incumbent $\bar{x} = (1, 0, 1, 1, 0)$ with $k = 2$. Is $x = (1, 1, 1, 0, 0)$ inside the neighborhood?
- (Intermediate) Make the example iterative: run ALNS for 300 iterations, solve the set partitioning, restart ALNS from its solution with the same pool, and repeat five times. Compare with a single ALNS run of 1,500 iterations.
- (Interview) "A colleague suggests replacing our ALNS with 'just give the whole model to Gurobi'. How do you decide?" Hint: consider instance size, the time limit, the MIP gap reached within that time, and hybrids: warm-starting the MIP with the heuristic solution, or using the MIP inside the heuristic.
Constraint Programming Integration
Why it matters. Some problems are dominated by logical and combinatorial constraints rather than by a linear objective: scheduling with machine conflicts and precedences, timetabling, rostering with complex labor rules, configuration. MIP formulations of these constraints are often weak. Section 10.3 noted that big-M scheduling models have poor LP relaxations. Constraint programming (CP) reasons about these constraints directly and is often the best tool for them. Combining CP with MIP and with heuristics such as LNS gives some of the strongest scheduling solvers available.
Theory. A CP model has variables with finite domains, and constraints among them. A CP solver alternates two steps:
- Propagation: each constraint removes from the domains the values that cannot be part of any solution, and the removals trigger further propagation until nothing changes.
- Search: when propagation stops, the solver branches (for example, fixes a variable or orders two tasks) and backtracks on failure.
The power of CP comes from global constraints, which capture a whole substructure with a dedicated propagation algorithm:
| Global constraint | Meaning | Propagation technique |
|---|---|---|
AllDifferent(x1..xn) | all variables take distinct values | bipartite matching |
NoOverlap / disjunctive | tasks on one machine do not overlap | edge finding, not-first/not-last |
Cumulative | resource usage stays within capacity at all times | time tabling, energetic reasoning |
Circuit | successor variables form one Hamiltonian cycle | subtour elimination on the graph |
Element, Table | indexing, allowed combinations | domain consistency |
A small example of disjunctive propagation, with two tasks on one machine:
task A: duration 4, may start in [0, 2], must finish by 6
task B: duration 3, may start in [2, 4], must finish by 7
try "B before A": B ends >= 5, so A starts >= 5 > 2 -> impossible
hence A before B: B starts >= A's end >= 4 -> B starts exactly at 4
A must end by B's start 4 -> A starts exactly at 0
propagation alone fixes the whole schedule, with no search
Interval variables (start, duration, end, and optional presence) make scheduling models compact. Modern solvers such as OR-Tools CP-SAT combine CP propagation with SAT-style conflict learning (lazy clause generation), LP relaxations for bounds, and portfolios of LNS workers.
| MIP | CP | |
|---|---|---|
| Strength | LP bounds, cutting planes, optimality gaps | propagation of logical and combinatorial constraints |
| Natural objectives | linear costs and profits | makespan, lateness, feasibility |
| Natural constraints | linear inequalities, flows | disjunctions, sequences, all-different, resources |
| Weak spots | big-M disjunctions, logical rules | weak bounds for sum-type objectives |
Integration patterns.
- Logic-based Benders decomposition (Hooker and Ottosson): a MIP master problem assigns tasks to machines and minimizes cost, and CP subproblems check or optimize each machine's schedule. Infeasible or costly assignments return cuts ("not these tasks together on machine $k$") to the master.
- CP pricing in column generation: a column generation master (a MIP) chooses among shifts or crew pairings, and CP generates new feasible columns under complex rules.
- CP-based LNS: release part of a schedule, and let CP re-optimize it under all constraints (Section 11.2).
- LP inside CP: the LP relaxation acts as a global constraint for bounding, as in CP-SAT.
Example: a 10 x 10 job shop. The job shop of Section 10.3, with 10 jobs and 10 machines, random routes and processing times from 1 to 99, in the style of Taillard's benchmarks. The CP model uses interval variables and one NoOverlap constraint per machine. The MIP uses start times and one binary per pair of operations on the same machine, with big-M constraints. Each solver gets one thread and 30 seconds (times vary by machine).
import time
import numpy as np
from ortools.sat.python import cp_model
from scipy.optimize import milp, LinearConstraint, Bounds
# Random 10 x 10 job shop in the style of Taillard's benchmarks: each job visits every machine once,
# in a random order, with processing times from 1 to 99.
r = np.random.default_rng(0)
J = M = 10
jobs = [[(int(k), int(r.integers(1, 100))) for k in r.permutation(M)] for _ in range(J)]
H = sum(t for job in jobs for _, t in job) # horizon: all work done one after another
# ---- Constraint programming: interval variables + one global NoOverlap per machine
m = cp_model.CpModel()
start, intervals = {}, {k: [] for k in range(M)}
for j, job in enumerate(jobs):
for o, (k, t) in enumerate(job):
s = m.new_int_var(0, H, f"s{j}{o}")
start[j, o] = s
intervals[k].append(m.new_interval_var(s, t, s + t, f"op{j}{o}"))
if o > 0: # route order inside the job
m.add(s >= start[j, o - 1] + job[o - 1][1])
for k in range(M):
m.add_no_overlap(intervals[k]) # one global constraint replaces all big-M pairs
makespan = m.new_int_var(0, H, "Cmax")
for j, job in enumerate(jobs):
m.add(makespan >= start[j, M - 1] + job[M - 1][1])
m.minimize(makespan)
solver = cp_model.CpSolver()
solver.parameters.num_workers = 1 # single thread, like HiGHS
solver.parameters.max_time_in_seconds = 30
t0 = time.perf_counter()
status = solver.solve(m)
print(f"CP-SAT : makespan {solver.objective_value:.0f}, bound {solver.best_objective_bound:.0f}, "
f"{solver.status_name(status)}, {solver.num_branches:,} branches, {time.perf_counter() - t0:.2f} s")
# ---- MILP: start times + one binary per pair of operations on the same machine (big-M)
ops = [(j, o) for j in range(J) for o in range(M)]
idx = {op: i for i, op in enumerate(ops)}
pairs = [(a, b) for a in ops for b in ops if a < b and jobs[a[0]][a[1]][0] == jobs[b[0]][b[1]][0]]
nv = len(ops) + len(pairs) + 1 # start times, order binaries, makespan
rows, lo, hi = [], [], []
def row(coefs, lb, ub):
v = np.zeros(nv)
for i, c in coefs: v[i] += c
rows.append(v); lo.append(lb); hi.append(ub)
for j, job in enumerate(jobs):
for o in range(1, M): # route order
row([(idx[j, o], 1), (idx[j, o - 1], -1)], job[o - 1][1], np.inf)
row([(nv - 1, 1), (idx[j, M - 1], -1)], job[M - 1][1], np.inf) # Cmax >= completion of the job
for p, (a, b) in enumerate(pairs): # y = 1: a before b, y = 0: b before a
y, pa, pb = len(ops) + p, jobs[a[0]][a[1]][1], jobs[b[0]][b[1]][1]
row([(idx[b], 1), (idx[a], -1), (y, -H)], pa - H, np.inf) # s_b >= s_a + p_a - H (1 - y)
row([(idx[a], 1), (idx[b], -1), (y, H)], pb, np.inf) # s_a >= s_b + p_b - H y
c = np.zeros(nv); c[-1] = 1
integ = np.zeros(nv); integ[len(ops):len(ops) + len(pairs)] = 1
ub = np.full(nv, H); ub[len(ops):len(ops) + len(pairs)] = 1
cons = LinearConstraint(np.array(rows), lo, hi)
lp = milp(c, constraints=cons, bounds=Bounds(0, ub)) # LP relaxation
t0 = time.perf_counter()
res = milp(c, constraints=cons, bounds=Bounds(0, ub), integrality=integ, options={"time_limit": 30})
print(f"MILP : makespan {res.fun:.0f}, bound {res.mip_dual_bound:.0f}, {res.message.split('.')[0].lower()}, "
f"{len(pairs)} binaries, LP relaxation {lp.fun:.0f}, {time.perf_counter() - t0:.2f} s")
# CP-SAT : makespan 910, bound 910, OPTIMAL, 25,656 branches, 4.04 s
# MILP : makespan 963, bound 751, time limit reached, 450 binaries, LP relaxation 656, 30.02 s
CP-SAT proves the optimal makespan of 910 in a few seconds. The MIP stops at the time limit with a schedule of 963 (6% worse) and a lower bound of 751, a remaining gap of 22%. Its LP relaxation is only 656, which is just the length of the longest job: with fractional ordering variables, the big-M constraints effectively disappear. Even the busiest machine's total load (699) is a better bound, and a CP solver obtains that kind of bound from propagation alone. For disjunctive scheduling, CP is the natural tool. For cost-driven decisions with linear structure, MIP is. Real systems often combine the two.
Key takeaways
- CP solves problems by propagating constraints on variable domains and searching with backtracking. Global constraints carry most of its power.
- CP excels at scheduling, sequencing and logical constraints, where MIP formulations with big-M are weak.
- Integration patterns include logic-based Benders decomposition, CP pricing in column generation, CP-based LNS, and LP relaxations inside CP.
- Modern hybrid solvers such as CP-SAT combine propagation, conflict learning, LP bounds and LNS in one engine.
Practice
- (Beginner) Two tasks share one machine. Task A has duration 3, may start in $[0, 4]$ and must end by 7. Task B has duration 4, may start in $[1, 3]$ and must end by 7. Use disjunctive reasoning to narrow both start-time domains.
- (Intermediate) Add the constraints $C_{\max} \ge \sum_{o \text{ on } k} p_o$ for every machine $k$ to the MIP and rerun it. How much do the LP relaxation and the final bound improve? Then give CP-SAT 8 workers and compare its time.
- (Interview) "We must roster 200 nurses under dozens of labor rules, minimizing overtime cost. MIP or CP?" Hint: consider which rules are sequence-like (no more than 5 night shifts in a row) and which parts of the objective are linear. Discuss column generation with CP or dynamic programming for the pricing step.
Machine Learning for Combinatorial Optimization
Why it matters. In practice, the same optimization problem is solved again and again on instances from one distribution: the same city's deliveries every morning, or the same network's schedule every week. Classical algorithms treat each instance as brand new. Machine learning can exploit what the instances have in common: which edges tend to appear in good tours, which branching decisions pay off, which heuristic works best on which instance. It is one of the most active research areas in OR. It is also an area where careful comparison with strong classical baselines is essential.
Theory. Bengio, Lodi and Prouvost (2021) distinguish three ways to combine learning and optimization:
| Approach | What is learned | Examples |
|---|---|---|
| End-to-end | a model that outputs solutions directly | pointer networks (Vinyals et al., 2015), attention models trained by reinforcement learning (Kool et al., 2019), edge "heatmaps" from graph neural networks |
| Learning to configure | which algorithm or parameters to use for an instance | algorithm selection, parameter prediction, warm starts |
| ML inside the algorithm | a policy for a decision the solver makes repeatedly | learned branching that imitates strong branching (Gasse et al., 2019), node and cut selection, learned candidate lists, learned LNS destroy operators |
Learning can be supervised, imitating good solutions or expensive expert decisions, or reinforcement learning, which optimizes the solution cost directly by trial and error. A related idea is predict-then-optimize, where ML forecasts uncertain inputs (Section 10.4). Decision-focused learning trains the forecast for the quality of the resulting decision rather than for accuracy.
Practical cautions:
- Feasibility. Model outputs can violate constraints. The safest designs let the model guide a solver (ordering, pruning, warm starts) rather than replace it.
- Generalization. Models trained on 50-node instances often degrade on 500 nodes or on a different spatial distribution. Scale-free features and size-invariant architectures help.
- Baselines. Compare against the best classical methods (LKH, HGS, Gurobi with tuned parameters) under equal time limits, and count the training time.
- Guarantees. Learned components inside exact solvers (branching, cuts) keep optimality proofs. End-to-end models do not.
Example: learning candidate edges for the TSP. Heuristics such as LKH, and exact solvers too, restrict the search to a sparse candidate graph. If an optimal edge is missing from it, it cannot be found, and if the graph is too dense, the search is slow. The code trains a logistic regression on 40 random 50-city instances, labeled with their exact optimal tours. Each edge between near neighbors has five features: its relative length, the neighbor rank at both ends, membership in the minimum spanning tree (Section 3.3), and its value in the LP relaxation with degree-2 constraints only. The test measures how many edges per city each scoring rule must keep to retain 95% and 100% of the optimal edges, on unseen instances with 50 cities and with 100 cities, a size never seen in training.
import numpy as np
from itertools import combinations
from scipy.optimize import milp, LinearConstraint, Bounds
from scipy.sparse import csr_matrix
from scipy.sparse.csgraph import connected_components, minimum_spanning_tree
def optimal_tour_edges(D):
"""Exact TSP by the subtour-cut loop of Section 4.3; returns the set of tour edges."""
n = len(D); E = list(combinations(range(n), 2)); c = np.array([D[e] for e in E])
rows = [[k for k, e in enumerate(E) if v in e] for v in range(n)]
A = [np.isin(np.arange(len(E)), r).astype(float) for r in rows]; lb = [2] * n; ub = [2] * n
while True:
x = milp(c, constraints=LinearConstraint(np.array(A), lb, ub),
integrality=np.ones(len(E)), bounds=Bounds(0, 1)).x
used = [e for e, v in zip(E, x) if v > 0.5]
g = csr_matrix((np.ones(len(used)), ([i for i, _ in used], [j for _, j in used])), shape=(n, n))
k, lab = connected_components(g, directed=False)
if k == 1: return set(used)
for comp in range(k): # cut every subtour: sum inside S <= |S| - 1
S = set(np.flatnonzero(lab == comp))
A.append(np.array([float(i in S and j in S) for i, j in E])); lb.append(0); ub.append(len(S) - 1)
def edge_features(D):
"""Features of every edge (i, j), i < j, among each city's 10 nearest."""
n = len(D)
rank = np.argsort(np.argsort(D, axis=1), axis=1) # rank[i, j]: j is i's rank-th nearest
nn1 = np.sort(D, axis=1)[:, 1] # distance to the nearest city
E = sorted({tuple(sorted((i, j))) for i in range(n) for j in np.argsort(D[i])[1:11]})
mst = minimum_spanning_tree(D).toarray() # cheap structure: the MST (Section 3.3)
A = np.zeros((n, len(E))) # ... and the 2-matching LP:
for k, (i, j) in enumerate(E): A[i, k] = A[j, k] = 1 # degree 2, no subtour cuts
lp = milp([D[e] for e in E], constraints=LinearConstraint(A, 2, 2), bounds=Bounds(0, 1)).x
X = np.array([[np.log(D[i, j] / np.mean(nn1)), # length relative to the instance
np.log(min(rank[i, j], rank[j, i])), # mutual nearness, best side
np.log(max(rank[i, j], rank[j, i])), # mutual nearness, worst side
float(mst[i, j] + mst[j, i] > 0), # edge is in the MST
lp[k]] # LP relaxation value of the edge
for k, (i, j) in enumerate(E)])
return E, X
def instance(seed, n):
p = np.random.default_rng(seed).uniform(0, 100, size=(n, 2))
return np.linalg.norm(p[:, None] - p[None, :], axis=2)
# ---- training data: 40 random 50-city instances, labels from exact solutions
Xs, ys = [], []
for s in range(40):
D = instance(1000 + s, 50)
E, X = edge_features(D); opt = optimal_tour_edges(D)
Xs.append(X); ys.append(np.array([e in opt for e in E], dtype=float))
X, y = np.vstack(Xs), np.concatenate(ys)
mu, sd = X.mean(axis=0), X.std(axis=0)
Z = np.hstack([(X - mu) / sd, np.ones((len(X), 1))])
w = np.zeros(Z.shape[1])
for _ in range(3000): # logistic regression by gradient descent
p = 1 / (1 + np.exp(-Z @ w))
w -= 0.5 * Z.T @ (p - y) / len(y)
print(f"training edges {len(y):,}, of which in an optimal tour {y.mean():.1%}; weights {np.round(w, 2)}")
def edges_needed(D, opt, score):
"""Keep the best-scored edges first: how many are needed to keep 95% / all optimal edges?"""
E, X = edge_features(D)
hit = np.array([E[k] in opt for k in np.argsort(-score(E, X))])
found = np.cumsum(hit) # (an optimal edge outside the 10 nearest is lost anyway)
return np.searchsorted(found, 0.95 * found[-1]) + 1, np.searchsorted(found, found[-1]) + 1
scorers = {"shortest first": lambda E, X: -X[:, 0],
"neighbor rank": lambda E, X: -X[:, 1] - 0.01 * X[:, 0], # k-nearest-neighbor lists
"learned model": lambda E, X: np.hstack([(X - mu) / sd, np.ones((len(X), 1))]) @ w}
for n in (50, 100): # unseen instances, and a size never trained on
tests = [instance(2000 + s, n) for s in range(10)]
opts = [optimal_tour_edges(D) for D in tests]
for name, f in scorers.items():
need = np.array([edges_needed(D, o, f) for D, o in zip(tests, opts)]) / n
print(f"n = {n:3d}, {name:14s}: edges per city to keep 95% of optimal edges {need[:, 0].mean():.2f}, "
f"for all {need[:, 1].mean():.2f}")
# training edges 11,944, of which in an optimal tour 16.7%; weights [ 0.25 -0.87 -0.68 0.39 0.73 -2.78]
# n = 50, shortest first: edges per city to keep 95% of optimal edges 3.43, for all 5.06
# n = 50, neighbor rank : edges per city to keep 95% of optimal edges 2.63, for all 3.80
# n = 50, learned model : edges per city to keep 95% of optimal edges 2.52, for all 3.72
# n = 100, shortest first: edges per city to keep 95% of optimal edges 3.45, for all 5.09
# n = 100, neighbor rank : edges per city to keep 95% of optimal edges 2.37, for all 4.20
# n = 100, learned model : edges per city to keep 95% of optimal edges 2.18, for all 3.78
Only 16.7% of the candidate edges belong to an optimal tour. The learned weights make sense: a high LP value and MST membership raise the score, while a high neighbor rank lowers it. Keeping the globally shortest edges is a poor rule, because short edges cluster in dense regions. Rank-based k-nearest-neighbor lists are much better. The model is better still: to keep every optimal edge on the 100-city instances it needs 3.78 edges per city, against 4.20 for rank and 5.09 for shortest-first. A sparser candidate graph speeds up every later search. Because the features are scale-free, the model transfers to the larger instances without retraining. Most of the predictive power, however, comes from classical OR quantities (neighbor ranks, MST, LP values). Here the learned model refines good hand-made features rather than replacing them, which is typical of successful applications.
Key takeaways
- ML for combinatorial optimization exploits regularities across instances from one distribution.
- Three approaches: end-to-end solution models, learning to configure algorithms, and learned policies inside solvers.
- Guiding a solver (pruning, ordering, warm starts) keeps feasibility and often keeps guarantees. Replacing the solver risks both.
- Always evaluate against strong classical baselines, on unseen and larger instances, and include training costs.
Practice
- (Beginner) Why do the features use logarithms of relative lengths and ranks instead of raw distances? What would happen to a model trained on raw distances if the test cities were spread over a 1,000 x 1,000 square?
- (Intermediate) Use the model's scores to build tours with the greedy edge heuristic (accept edges in order of score, keeping degrees at most 2 and no early cycles), then apply 2-opt. Compare the tour lengths with greedy edge by distance.
- (Interview) "A startup claims their neural network solves our vehicle routing problems better than our current solver. What do you ask before believing it?" Hint: instance sizes and distribution, time limits, which baseline and with which settings, feasibility of every output, training data and cost, and performance on your own instances.
Algorithm Configuration and Tuning
Why it matters. Every heuristic in this chapter has parameters: tabu tenure, cooling rate, population size, RCL size, operator choices. Section 11.2 showed that a poor value can make the same algorithm look useless (tabu tenure 0 or 60, $\beta = 0$ in ACO). MIP solvers also have hundreds of parameters. Choosing them by hand is slow and biased. Algorithm configuration turns the choice into an optimization problem that is solved automatically, with proper experimental methodology.
Theory. Let $\theta \in \Theta$ be a parameter configuration of algorithm $A$, and $\mathcal{D}$ a distribution of instances. The configuration problem is
$$ \theta^* = \arg\min_{\theta \in \Theta} \ \mathbb{E}{I \sim \mathcal{D}} \big[ m(A\theta, I) \big], $$
where $m$ is a performance measure such as the solution cost after a fixed time, or the time to reach a target. The expectation is estimated on training instances, and the result is validated on separate test instances. This is the same discipline as out-of-sample testing in forecasting (Section 10.4). What makes the problem hard:
- each evaluation is an expensive run of a heuristic,
- results are noisy: randomized algorithms and different instances give different values,
- parameters are mixed (continuous, integer, categorical) and sometimes conditional (a cooling rate only matters if SA is chosen),
- configurations can overtune to the training instances.
| Method | Idea | Tools |
|---|---|---|
| Grid or manual search | try a few values per parameter | simple, but scales badly |
| Random search | sample configurations at random | strong baseline in high dimension |
| Racing | evaluate many configurations on a growing set of instances, drop the clearly worse ones early | F-Race, irace |
| Model-based (Bayesian) optimization | fit a surrogate model of performance and sample where it promises most | SMAC, Optuna, Hyperopt |
| Local search in configuration space | hill climbing with restarts over configurations | ParamILS |
Racing rests on a simple observation: bad configurations are recognizable after a few instances, so there is no need to run them on all of them. Successive halving is the simplest version. Start with many configurations and a few instances, keep the better half, double the number of instances, and repeat. F-Race and irace use statistical tests (Friedman's test) to decide which configurations to drop, and irace iterates by sampling new configurations around the survivors.
Related distinctions:
- Offline tuning (before deployment) versus online parameter control (during the run, such as ALNS weights and reactive tabu search).
- Configuration (one setting for the whole distribution) versus per-instance algorithm selection (choose the method from instance features).
Example: racing the parameters of simulated annealing. SA on random 60-city TSPs with a fixed budget of 10,000 moves has two parameters: p0, the initial acceptance probability of a typical uphill move, and t_ratio, the final temperature as a fraction of the initial one. 32 random configurations are raced by successive halving on 32 training instances. The winner, the default and all sampled configurations are then evaluated on 40 unseen instances with 2 seeds each.
import numpy as np
def instance(seed, n=60):
p = np.random.default_rng(seed).uniform(0, 100, size=(n, 2))
return np.linalg.norm(p[:, None] - p[None, :], axis=2)
def sa(D, p0, t_ratio, moves=10_000, seed=0):
"""Simulated annealing with 2-opt moves and a FIXED move budget. Parameters:
p0 - initial probability of accepting an average uphill move (sets T0)
t_ratio - final temperature as a fraction of T0 (sets the cooling rate)"""
n, r = len(D), np.random.default_rng(seed)
t = list(r.permutation(n)); cur = sum(D[t[k - 1], t[k]] for k in range(n)); best = cur
I, J = np.sort(r.integers(0, n, size=(2, moves)), axis=0)
U = r.random(moves)
scale = 0.5 * D[np.triu_indices(n, 1)].mean() # rough size of an uphill move
T = -scale / np.log(p0) # accept such a move with probability p0
cool = t_ratio ** (1 / moves)
for i, j, u in zip(I, J, U):
if j - i >= 2 and not (i == 0 and j == n - 1):
a, b, c, d = t[i], t[i + 1], t[j], t[(j + 1) % n]
delta = D[a, c] + D[b, d] - D[a, b] - D[c, d]
if delta < 0 or u < np.exp(-delta / T):
t[i + 1:j + 1] = t[i + 1:j + 1][::-1]; cur += delta
best = min(best, cur)
T *= cool
return best
# ---- configuration by racing (successive halving) on TRAINING instances
r = np.random.default_rng(7)
configs = [(float(10 ** r.uniform(-2, -0.05)), float(10 ** r.uniform(-5, -1))) for _ in range(32)]
train = [instance(100 + s) for s in range(32)]
cache = {}
def score(cfg, k): # tour length of config cfg on training instance k
if (cfg, k) not in cache: cache[cfg, k] = sa(train[k], *cfg, seed=k)
return cache[cfg, k]
alive, used = configs, 2
while len(alive) > 1:
L = np.array([[score(c, k) for k in range(used)] for c in alive])
ranks = L.argsort(axis=0).argsort(axis=0).mean(axis=1) # mean rank per config, instance by instance
alive = [alive[i] for i in np.argsort(ranks)[:len(alive) // 2]]
print(f"{used:2d} instances: {len(L):2d} configs raced, keep {len(alive):2d}; leader p0 = {alive[0][0]:.3f}, "
f"t_ratio = {alive[0][1]:.1e}")
used = min(2 * used, len(train))
print(f"runs spent on tuning: {len(cache)} (a full grid of 32 configs x 32 instances would be 1024)")
# ---- evaluation on UNSEEN test instances, against the default and every sampled config
test = [instance(900 + s) for s in range(40)] # many instances x 2 seeds: test noise matters
default = (0.5, 1e-3)
perf = {c: np.mean([sa(D, *c, seed=s) for D in test for s in (1, 2)]) for c in configs + [default]}
ranked = sorted(perf, key=perf.get)
for name, c in [("best sampled config (hindsight)", ranked[0]), ("tuned by racing", alive[0]),
("default p0 = 0.5, t_ratio = 1e-3", default),
("median sampled config", ranked[len(ranked) // 2]), ("worst sampled config", ranked[-1])]:
print(f"{name:33s}: mean test length {perf[c]:6.1f} (rank {ranked.index(c) + 1:2d} of {len(ranked)})")
# 2 instances: 32 configs raced, keep 16; leader p0 = 0.038, t_ratio = 3.1e-02
# 4 instances: 16 configs raced, keep 8; leader p0 = 0.038, t_ratio = 3.1e-02
# 8 instances: 8 configs raced, keep 4; leader p0 = 0.038, t_ratio = 3.1e-02
# 16 instances: 4 configs raced, keep 2; leader p0 = 0.038, t_ratio = 3.1e-02
# 32 instances: 2 configs raced, keep 1; leader p0 = 0.098, t_ratio = 3.1e-02
# runs spent on tuning: 192 (a full grid of 32 configs x 32 instances would be 1024)
# best sampled config (hindsight) : mean test length 658.4 (rank 1 of 33)
# tuned by racing : mean test length 660.0 (rank 4 of 33)
# default p0 = 0.5, t_ratio = 1e-3 : mean test length 669.7 (rank 24 of 33)
# median sampled config : mean test length 668.5 (rank 17 of 33)
# worst sampled config : mean test length 817.7 (rank 33 of 33)
Racing spends 192 runs instead of the 1,024 of a full grid over the same instances, because most configurations are dropped after 2 or 4 instances. On the test set, the tuned configuration ranks 4th of 33, within 0.3% of the best configuration in hindsight. The default ranks 24th, 1.5% worse than the tuned configuration. The worst sampled configuration ($p_0 = 0.87$) is 24% worse than the best: it starts so hot that its 10,000 moves are wasted on random wandering. The parameters matter, and the main gain of tuning is often avoiding such bad regions. The two finalists of the race differ by only 0.1% on the test set, so the final decision between them hardly matters.
The test itself must also be large enough. In an earlier version of this experiment, evaluating on only 12 test instances with one seed each, the same tuned configuration ranked 13th. That ranking was noise.
Key takeaways
- Algorithm configuration minimizes expected performance over an instance distribution. Tune on training instances and report results on separate test instances.
- Racing (successive halving, F-Race, irace) and model-based methods (SMAC) save most of the cost of grid search by dropping bad configurations early.
- Results are noisy. Use enough instances and seeds, fixed budgets, and statistical comparisons.
- Offline tuning sets parameters before deployment. Online control (ALNS weights, reactive tabu) adapts them during the run.
Practice
- (Beginner) Why must the tuning instances and the test instances be different? Relate your answer to overfitting in forecasting (Section 10.4).
- (Intermediate) Replace successive halving by random search with the same budget: evaluate 6 random configurations on all 32 training instances (192 runs). Compare the test rank of its winner with the racing result. Repeat with three different random seeds.
- (Interview) "Our routing heuristic runs on 1,000 new instances every day. How would you set its parameters and keep them good over the next year?" Hint: build a representative training set from history, tune offline with racing or SMAC, keep a hold-out set, monitor solution quality for drift, schedule retuning, and consider per-instance algorithm selection.
<a id="12-simulation"></a>