Local Search and Simulated Annealing
Build local search and simulated annealing from scratch on an eight-customer routing instance, then compare them across 20 seeds against the enumerated optimum.
Machine Learning from ScratchAll chapters Local Search and Simulated AnnealingFull chapter · Sections in this chapter1 / 10 ↓
Reading tools Definitions, links & notes
Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.
Article links
Make inline references clickable
Local Search and Simulated Annealing
Many of the optimization problems we care about do not have a differentiable objective we can hand to a gradient method. A delivery planner wants to order eight stops into two delivery routes; a compiler wants to schedule instructions to minimize register pressure; a logistics team wants to re-sequence a pick list to shave a few kilometers off the day. In each case the decision variable is a discrete arrangement (a permutation, a schedule, a subset), and the objective is defined on that arrangement. Gradient descent cannot move "a little bit" in a discrete space, so we need a different strategy: propose a small structural change, evaluate the result, and keep it if it helps.
Local search applies that idea to a chosen neighborhood. We start from some feasible solution, examine a curated set of nearby solutions obtained by a small edit called a neighborhood move, and step to a better one. If we only ever accept improvements, the procedure is called hill climbing (or descent, when minimizing). Its strict-improvement rule is easy to trace on a small routing example. Its weakness is that it can get trapped: it stops improving even though some distant solution is far better.
Simulated annealing addresses that trap by occasionally accepting a worse solution, with a probability that shrinks over time as a virtual "temperature" cools. The name draws on physical annealing: temperature controls thermal motion, and cooling affects how a material settles into lower-energy states. In optimization, the random uphill moves early on give the search a chance to escape a basin, and the gradual cooling lets it settle back down at the end.
In this chapter we build both algorithms from primitives on a concrete eight-customer routing instance, watch them in action, and compare them, under an equal evaluation budget and across twenty random seeds, against a deterministic construction heuristic and a random baseline. The instance is new and deliberately small so that we can enumerate every one of the 40,320 permutations and know the true optimum. That reference lets us measure the distance from each reported objective to the enumerated minimum. Reference enumeration, figure diagnostics, and post-run validation are outside the search budgets. Every method also builds the same nearest-neighbor start using 36 candidate-arc lookups; its complete objective score counts as evaluation one.
Advantages
For the finite combinatorial setting studied here, local search needs a representation, an objective evaluator, and a neighborhood generator. There is no requirement that the objective be differentiable, convex, or even continuous, and no matrix of derivatives to compute or invert. These interfaces can be defined for routing, scheduling, packing, and assignment. Continuous relaxations and specialized exact methods remain alternatives.
For a symmetric Euclidean tour, a within-route segment reversal replaces two boundary edges and preserves the lengths of internal edges. Selected reversals can remove crossings and shorten a route; arbitrary reversals need not improve it. This is the geometric idea behind 2-opt. Restarts provide additional starting points, at the cost of their construction and neighborhood evaluations.
Simulated annealing adds a controlled escape mechanism without much extra code. At a fixed positive temperature, symmetric Metropolis proposals leave the finite-state distribution proportional to invariant. This is a stationary-distribution statement, not a claim that a short run already samples from equilibrium. The search state uses current and incumbent permutations, while our experiment also stores the distance matrix, score histories, and repeated-run traces.
Finally, both methods are trivially parallelizable across seeds. Because each run is independent, running twenty seeds is embarrassingly parallel, and the spread of their final objectives describes seed variation for this fixed instance and algorithm configuration.
Disadvantages
Hill climbing's central failure is the local optimum. Because it only accepts strict improvements, it halts the moment every neighbor is at least as bad as the current point, even if a better solution sits one move further away. A local minimum can occur before the global optimum is reached. The fix, restarts, helps, but it spends the budget on independent attempts rather than systematically exploring, and there is no guarantee that the next restart will find a better basin. Perturbations and larger neighborhoods offer other ways to leave a local minimum.
Annealing permits uphill moves instead of stopping solely because no improving neighbor is available; a finite evaluation budget still terminates the run. The choice of initial temperature, cooling schedule, and number of moves at each temperature strongly affects quality. Too hot and the search wanders; too cold and it behaves like a slow hill climber. Classical asymptotic convergence results include sufficiently slow logarithmic schedules under finite-state and energy-barrier conditions. They do not certify a result after our 300-evaluation run; an independent bound or exhaustive reference is needed for that. We must therefore judge quality empirically, not by appeal to the asymptotic theory.
Both methods are also sensitive to the representation and the neighborhood. If the neighborhood is disconnected, a local walk from a fixed start stays in its component. Full-permutation random restarts can enter other components, but each epoch still refines only the component it starts in. If it is too wide, a full scan of the neighborhood becomes expensive. If it is too narrow, progress stalls. A move that changes the split between two vehicles in a routing problem has a very different effect from a move that only re-orders stops within one vehicle, and picking the wrong move set silently restricts local refinement between restarts. Finally, an equal number of objective evaluations does not imply equal runtime: a nearest-neighbor scan examines arcs cheaply, while each objective call here checks ten travel arcs for feasibility and traverses ten arcs again for distance, so two methods with the same nominal budget can take very different wall-clock time.
Formula
There are four objects to define, and they fit together in a natural order:
- A representation says what a candidate solution is.
- An objective assigns a number to each candidate.
- A neighborhood says which small edits produce the candidates the search is allowed to see next.
- An acceptance rule decides whether to move to a proposed neighbor.
The order matters. The representation fixes the universe of possible answers; the objective ranks them; the neighborhood defines the local graph the search can traverse; and the acceptance rule determines whether the walk is greedy or stochastic. Each choice below is made with the same routing instance in mind, so that the formal definitions can be checked immediately against concrete numbers.
Solution representations and feasible neighborhoods
We encode a solution as a permutation of the customer identifiers . The decoder splits the permutation into two routes of four customers each, because the total demand is eight and each of the two vehicles has capacity four:
where node is the depot. Because each vehicle carries exactly four customers, every permutation of the eight customers automatically produces two capacity-feasible routes, so capacity alone never rules out a permutation. Time-window and depot-return feasibility are checked separately, and on this instance they are deliberately loose. Each customer has demand one, service takes one time unit, and all customer and depot windows are . Travel time equals Euclidean distance in the chosen units. The coordinates fit in an 8-by-6 rectangle, so every arc is at most 10 units. Each four-customer route uses five arcs and four service periods, giving a return-time bound of . With window openings at zero, no waiting is needed; every valid permutation therefore satisfies these capacity and time limits. The loose windows isolate search behavior; this example does not test restrictive time-window handling.
The objective sums the Euclidean arc lengths over both closed tours:
where:
- : a permutation of the eight customer identifiers
- : the customer visited at position in the permutation
- : the Euclidean distance between two coordinates
- the first line traverses the first route's internal arcs (positions to ), and the two depot arcs and close that tour
- the second line traverses the second route's internal arcs (positions to ), and the two depot arcs and close that tour
There is no lateness penalty and no service-time term in the objective; service times only constrain the schedule validation.
A representation by itself is not enough: the search needs a notion of "nearby." That is the neighborhood. Here the neighborhood is the set of all inclusive segment reversals. A move selects indices and replaces the segment with its reverse. This move set is not arbitrary. A within-route reversal preserves the undirected internal-edge lengths and rewires two boundary arc positions. Reversing an entire route can leave its aggregate undirected edge multiset unchanged. A reversal crossing the split changes route assignments and route-labelled depot attachments, even when the aggregate depot-edge multiset is unchanged, so it must be scored through the full decoder here. When the reversed segment crosses the position- boundary, the move changes how customers are assigned to vehicles. Thus the neighborhood can reassign customers between routes.
This neighborhood has two properties we rely on. First, adjacent reversals (with ) are simple swaps, and swaps alone connect the entire permutation space; the wider reversals only add shortcuts. Consequently the state graph is connected and every permutation is reachable from every other permutation. Second, the wider reversals are geometrically meaningful: reversing a contiguous piece of a route is the 2-opt move, which can uncross edges and shorten a tour. A neighborhood that only reversed within each vehicle would be disconnected across the route boundary and would forbid exactly the cross-split moves a routing search needs.
PythonIn [3]Expand
import numpy as np
COORDS = np.array(
[[0, 0], [1, 3], [2, 6], [4, 5], [7, 6], [8, 3], [6, 1], [3, 1], [5, 3]],
dtype=float,
)
DEMANDS = np.array([0, 1, 1, 1, 1, 1, 1, 1, 1])
CAPACITY = 4
WINDOWS = [(0.0, 80.0)] * 9
SERVICE = np.array([0.0, 1, 1, 1, 1, 1, 1, 1, 1])
diff = COORDS[:, None, :] - COORDS[None, :, :]
D = np.sqrt((diff**2).sum(axis=-1))PythonIn [4]Expand
def tour_distance(perm):
total = 0.0
for route in (perm[:4], perm[4:]):
prev = 0
for c in route:
total += D[prev, c]
prev = c
total += D[prev, 0]
return total
def is_feasible(perm):
values = np.asarray(perm)
if (
values.shape != (8,)
or not np.issubdtype(values.dtype, np.integer)
or set(values.tolist()) != set(range(1, 9))
):
return False
for route in (perm[:4], perm[4:]):
load = 0
t = 0.0
prev = 0
for c in route:
load += DEMANDS[c]
if load > CAPACITY:
return False
arrival = t + D[prev, c]
start = max(arrival, WINDOWS[c][0])
if start > WINDOWS[c][1]:
return False
t = start + SERVICE[c]
prev = c
if t + D[prev, 0] > WINDOWS[0][1]:
return False
return True
def objective(perm):
return tour_distance(perm) if is_feasible(perm) else np.inf
def nearest_neighbor():
remaining = list(range(1, 9))
perm = []
for _ in range(2):
prev = 0
for _ in range(4):
best_c = min(remaining, key=lambda c: (D[prev, c], c))
perm.append(best_c)
remaining.remove(best_c)
prev = best_c
return tuple(perm)PythonIn [5]Expand
nn_perm = nearest_neighbor()
nn_length = objective(nn_perm)ConsoleOut [6]Expand
Nearest-neighbor permutation: (1, 7, 8, 3, 6, 5, 4, 2) Nearest-neighbor tour length: 40.856347 Nearest-neighbor permutation distance summary: mean arc = 4.085635
The nearest-neighbor heuristic, run through the decoder and objective above, returns a warm-start tour of length described by the printed output above. Because the customer count is only eight, we can enumerate every permutation and obtain the true minimum, which we will use as an exact reference later.
PythonIn [7]Expand
from itertools import permutations
best_exact, best_exact_perm = np.inf, None
for p in permutations(range(1, 9)):
d = objective(p)
if d < best_exact - 1e-15:
best_exact, best_exact_perm = d, pConsoleOut [8]Expand
Exact optimum: 37.35597706914 Achieved by: (6, 5, 4, 8, 1, 2, 3, 7) Nearest neighbor gap: 3.500370
Enumerating all permutations confirms the optimum is printed above, attained by the permutation printed above. This exhaustive reference search is outside the three budgeted heuristic runs, and we exclude it from every method's evaluation budget.
Hill climbing, local optima, and restarts
With the representation and neighborhood fixed, we can define an ideal steepest-descent rule: from the current permutation, choose a minimum-score neighbor and step only if it improves enough. The ideal strict-improvement rule is
where:
- : the current permutation at iteration
- : the set of all segment reversals of
- : the objective (total tour length) defined above
- : a small tolerance (here ) to guard against floating-point noise. The tolerance is not decoration: two mathematically equal tour lengths may differ in their last bits, and without a margin round-off can cause unproductive moves between tours that are equal for all practical purposes. Deterministic strict descent still cannot cycle, because every accepted step lowers the fixed computed score. If no neighbor improves by more than , the step is undefined: no improvement exceeding that tolerance is available in . This is a tolerance-based stopping point, which need not be an exact mathematical local minimum.
The code below uses an ordered, tolerance-aware scan rather than an exact argmin. The first scored neighbor becomes the scan candidate; a later neighbor replaces it only when its score is lower by more than . After a completed scan, the candidate is installed only if it improves the current score by more than . Differences within retain the earlier neighbor. The chosen candidate can therefore differ from the exact minimum; the final guard can even reject a scan when a later neighbor would have passed that guard. A budget-truncated scan returns the global best already scored, without installing a local step.
A local optimum is not a global one: there may exist a with that is several reversals away. Restarts handle local stopping by abandoning the current epoch when the tolerance-aware scan refuses to step and beginning a fresh epoch from a new permutation. Within an epoch, every accepted step strictly lowers the current computed score. The scan stopping rule is not an exact local-optimality certificate. The restart initialization is usually a fresh random permutation, and the overall best seen across all epochs is what we report.
Two subtleties matter because they change reported comparisons. First, when we scan the full neighborhood to select a neighbor, every neighbor evaluation consumes budget, including the neighbors we reject and including equal-cost duplicates. When the budget runs out mid-scan, the honest thing to do is return the best permutation already scored, not to spend one extra evaluation to finish the scan. Second, the fixed budget must include the evaluations spent starting the epochs, not only the improvement steps. If restart-initialization evaluations are excluded, the reported count understates the true scoring cost and can make the method appear more efficient than it is.
Annealing acceptance and temperature schedules
Simulated annealing uses the same representation, objective, and neighborhood. It proposes one random reversal instead of scanning for the best neighbor, then applies a stochastic acceptance rule. From the current state with objective , we propose a neighbor by sampling uniformly from the reversal neighborhood. We define the change in objective
and accept the proposal with probability
where:
- : the change in objective,
- : the temperature after evaluations (a positive scale that controls how freely uphill moves are accepted)
An improving move is always accepted; a worsening move is accepted with probability that falls as shrinks. When is large the search behaves like a random walk that pays little attention to the objective; as approaches zero, the probability of accepting any fixed positive increase approaches zero. Equal-cost moves still have acceptance probability one, so this limit differs from strict-improvement hill climbing.
The exponential acceptance rule can be read directly as a family of curves: for a fixed worsening move, higher temperature means more exploration, while lower temperature rapidly suppresses uphill moves.
PythonIn [9]Expand
T_GRID = np.geomspace(0.01, 5.0, 240)
DELTAS_TO_PLOT = np.array([0.1, 0.5, 1.0, 3.0])
ACCEPT_PROBS = np.exp(-DELTAS_TO_PLOT[:, None] / T_GRID[None, :])VisualizationOut [10]Expand

The acceptance rule is not an arbitrary hack. It is the Metropolis-Hastings rule for the Boltzmann (Gibbs) target
Here assigns probability to a permutation at a fixed . The proportionality sign means dividing by , with ranging over the finite permutation space. Every objective is finite on this deliberately feasible instance, so this normalizing sum is positive and finite.
To see why, note that if we propose symmetrically from (so ), then the ratio
and the Metropolis rule accepts with probability , which is exactly the formula above. This is why sampling two distinct indices, sorting them, and reversing the inclusive segment yields a symmetric proposal over all reversals: any reversal is as likely as its reverse.
The temperature must fall over time according to a cooling schedule. Here we choose geometric cooling:
where:
- : the initial temperature (the hot start)
- : the final temperature (near-greedy for worsening moves large relative to ; ties and sufficiently small increases can still be accepted)
- : the total evaluation budget
- : the evaluation index, running from to
For and , this schedule has endpoints and . Index zero is the initial nearest-neighbor score, where there is no acceptance test; the first proposal uses index one. Classical finite-state convergence results use schedules such as with a sufficiently large constant relative to the critical escape depth and suitable proposal conditions. Our finite geometric schedule has no certificate of global optimality by itself.
The schedule above has a simple shape that we can inspect before running the algorithm: most of the cooling happens early, and the final evaluations run at a near-greedy temperature.
PythonIn [11]Expand
K_EVALS = 300
T0 = 2.0
T1 = 0.01
EVAL_INDEX = np.arange(K_EVALS)
COOLING_TEMPS = T0 * (T1 / T0) ** (EVAL_INDEX / (K_EVALS - 1))VisualizationOut [12]Expand

The nearest-neighbor tour supplies the initial state. The initial temperature controls acceptance of worsening perturbations; it does not bound how far the current tour can deteriorate over several accepted moves. The separate incumbent retains the best tour already scored. A final temperature near zero favors polishing for worsening moves large relative to . The temperature scale must be checked against actual objective changes; these two endpoint values are choices for this experiment.
Mathematical Properties
Several properties are worth recording so we do not overclaim what these algorithms achieve. First, the permutation graph under reversals is connected: any permutation can be transformed into any other by a finite sequence of inclusive segment reversals, because swaps (adjacent reversals) alone generate all permutations. This guarantees that with a sufficiently long search simulated annealing can in principle visit any solution, which a disconnected neighborhood would not allow.
Second, for a fixed temperature , the Metropolis chain with symmetric proposals has as its stationary distribution. That gives simulated annealing a probabilistic meaning at a frozen temperature, where sufficiently equilibrated samples can approximate the Boltzmann distribution, but stationarity alone does not establish that a changing-temperature chain tracks those distributions. Decreasing temperature changes the transition rule and target distribution over time.
Third, under Hajek's conditions, a logarithmic cooling schedule that decreases slowly enough relative to the critical depth of the deepest nonglobal minimum makes the probability of being in the global-optimum set tend to one. This is an asymptotic statement under its hypotheses. It does not supply a finite-budget success probability for the schedule used here. In practice we measure quality empirically, across many seeds, and report the spread.
Finally, hill climbing with strict improvement terminates in finitely many steps: the objective strictly decreases by at least at every accepted step, and the solution space is finite, so the loop cannot run forever. Hypothetically, independent uniform restart states followed by the same deterministic complete descent produce identically distributed completed-epoch outcomes. Our first epoch has a deterministic nearest-neighbor start, and a fixed budget can truncate the last epoch, so those outcomes do not all share that distribution. The best objective observed cannot increase as more candidates are scored.
The formulas now define a complete search process: a permutation is a solution, the objective ranks it, the reversal neighborhood supplies local moves, hill climbing accepts only strict improvements, and annealing accepts improvements always but worsening moves with a temperature-controlled probability. The next section makes the neighborhood visible, so we can see why these rules behave differently on the same eight-customer instance.
Key Parameters
The three algorithms have small parameter surfaces. Hill climbing has just two: the budget and the seed (and an implicit restart rule). Simulated annealing adds the initial and final temperatures and , which set the geometric cooling curve. Random search adds nothing beyond a budget and a seed.
In our runs we hold the budget fixed at evaluations and run all three across the same seeds. Both hill climbing and annealing start from the nearest-neighbor solution, which costs one evaluation immediately. That common warm start controls initialization: different starts would confound the effect of the search with the choice of starting point.
For annealing, we choose and . The initial temperature is chosen so that a typical worsening move is accepted with probability early, but a worsening move (an illustrative increase in these distance units) is only accepted with probability . At a final temperature of , a increase has extremely small acceptance probability, although sufficiently small worsening moves can still be accepted. Finite positive temperature is not a strict-improvement rule.
Visualizing Local Search
We can see the mechanism directly on this instance. Before plotting, compute the nearest-neighbor tour and its best single-reversal neighbor, so the figure cells only render precomputed results.
PythonIn [13]Expand
nn_perm = nearest_neighbor()
best_nb, best_nb_score = None, np.inf
for i in range(8):
for j in range(i + 1, 8):
p = list(nn_perm)
p[i : j + 1] = p[i : j + 1][::-1]
p = tuple(p)
s = objective(p)
if s < best_nb_score - 1e-15:
best_nb, best_nb_score = p, s
base = objective(nn_perm)
deltas = []
for i in range(8):
for j in range(i + 1, 8):
p = list(nn_perm)
p[i : j + 1] = p[i : j + 1][::-1]
p = tuple(p)
deltas.append(objective(p) - base)
deltas.sort()The next two plots show the nearest-neighbor routes and the routes after the best single reversal, on matching axes.
PythonIn [14]Expand
def draw_routes(ax, perm):
from book_plot_style import PALETTE, polish_axes
for index, route in enumerate((perm[:4], perm[4:])):
pts = [0] + list(route) + [0]
ax.plot(
COORDS[pts, 0],
COORDS[pts, 1],
color=PALETTE[("primary", "teal")[index]],
lw=1.1,
marker="o",
markersize=4,
label=f"Route {index + 1}",
)
ax.scatter(
*COORDS[0],
s=65,
marker="s",
color=PALETTE["amber"],
zorder=4,
label="Depot",
)
for customer in range(1, 9):
ax.annotate(
str(customer),
COORDS[customer],
xytext=(5, 5),
textcoords="offset points",
)
ax.set_aspect("equal")
ax.set_xlim(-0.5, 8.8)
ax.set_ylim(-0.5, 6.8)
ax.set_xlabel("x (distance units)")
ax.set_ylabel("y (distance units)")
ax.legend(
loc="upper center", bbox_to_anchor=(0.5, -0.21), ncol=3, frameon=False
)
polish_axes(ax)VisualizationOut [15]Expand

VisualizationOut [16]Expand

For the nearest-neighbor tour there are distinct reversals. We score all of them and plot the change in objective sorted from most improving to most worsening. Any bar below the zero line is a downhill move that hill climbing could accept immediately. The bars above the line are moves that hill climbing must reject but annealing may accept with some probability. This neighborhood contains one improving move and three changes that are zero within the comparison tolerance. After moving, hill climbing must scan the new neighborhood before declaring a local minimum; the budgeted implementation then restarts rather than ending the experiment.
VisualizationOut [17]Expand

These figures together summarize why each method behaves as it does. Within each epoch, hill climbing accepts only a scan-selected move that lowers the current score by more than . It can stop with tiny negative bars still present, and the evaluation budget can end before a complete scan establishes even that stopping condition. Annealing occasionally climbs a small positive bar early on, which can drop it into a different basin with a lower floor. Neither strategy guarantees success. Both depend on how many evaluations they get and how good their starting point is.
Example
Let us work one reversal entirely by hand to see how the bookkeeping fits together. Start from the nearest-neighbor permutation
where the vertical bar separates the two vehicles. The current objective is the sum of the arc lengths on the two closed tours:
Substituting the coordinates (with and given in the array above) gives:
so the first route contributes . For the second route:
Carrying unrounded distances through the sums, the second route contributes and the total is . Summing the displayed rounded arc values instead can change the final reported decimal.
Now consider the reversal with and . Reversing positions of the permutation (that is, the segment ) produces
The reversal leaves the first customer, 1, fixed but changes the last customer of the first route from 3 to 7. Since distances are symmetric, reversing the internal edges keeps their lengths unchanged. The two changed boundary edges are replaced by , and replaced by . Therefore
Here is the new total distance minus the old total distance, and is the Euclidean distance between nodes and , including depot node . The negative sign means the reversal improves the tour. The two edges strictly inside the reversed segment retain their total length. The changed entering edge lengthens, while the changed return edge shortens by more, producing the net improvement. Checking every changed boundary edge is essential when a reversed segment ends at a route boundary.
PythonIn [18]Expand
def reverse_at(perm, i, j):
p = list(perm)
p[i : j + 1] = p[i : j + 1][::-1]
return tuple(p)
base = objective(nn_perm)
reverse_deltas = []
for i, j in [(1, 3), (2, 3), (0, 3), (3, 4)]:
p = reverse_at(nn_perm, i, j)
reverse_deltas.append((i, j, objective(p) - base, p))ConsoleOut [19]Expand
reversal (1,3): delta = -2.4637 perm = (1, 3, 8, 7, 6, 5, 4, 2) reversal (2,3): delta = +0.7225 perm = (1, 7, 3, 8, 6, 5, 4, 2) reversal (0,3): delta = +0.0000 perm = (3, 8, 7, 1, 6, 5, 4, 2) reversal (3,4): delta = +1.6437 perm = (1, 7, 8, 6, 3, 5, 4, 2)
ConsoleOut [20]Expand
Best single reversal so far: 38.3926 at (1, 3, 8, 7, 6, 5, 4, 2) Full tour improvement from NN to exact optimum: 3.5004
The output confirms that reversing to lowers the length by about 2.4637. Other reversals can worsen it, as the printed positive delta illustrates. Every complete candidate score costs one objective evaluation. A full scan contains 28 such calls; initialization and restarts also consume the 300-call limit. Annealing spends that same limit on its initial score and proposals, while NN-seeded random search uses one nearest-neighbor score followed by 299 random permutations.
For the acceptance rule, take a hypothetical worsening move with . At temperature the acceptance probability is , so even a small worsening is accepted most of the time. At temperature it is , so this worsening move is rarely accepted. And at temperature it is , effectively zero. That exponential controls acceptance of worsening proposals; annealing also requires a proposal rule, a temperature schedule, state updates, and a stopping budget. In the mathematical acceptance rule, a larger temperature accepts a fixed worsening move more often, while any finite positive temperature leaves a positive probability. In floating-point code, sufficiently small exponentials can underflow to zero.
Implementation in NumPy
The design center is a tiny Budget object that counts every objective evaluation and remembers the best permutation it has scored. This centralizes the accounting, so that hill climbing, annealing, and random search all spend budget by the same rule.
We now package the pieces into runnable algorithms. The search implementation uses NumPy and standard Python; figures use Matplotlib and the shared book style. We implement the acceptance and budgeting rules directly so their behavior is inspectable. The budget accounting, the fact that a scoring call is the currency of comparison, is what we want to control ourselves. It also makes the comparison honest: a method cannot quietly spend extra evaluations on rejected proposals or on restart initialization without those calls showing up in the trace.
PythonIn [21]Expand
class Budget:
def __init__(self, limit):
if (
isinstance(limit, (bool, np.bool_))
or not isinstance(limit, (int, np.integer))
or limit < 1
):
raise ValueError("limit must be a positive integer")
self.limit = int(limit)
self.evals = 0
self.scores = []
self.best_score = np.inf
self.best_perm = None
def score(self, perm):
if self.evals >= self.limit:
raise RuntimeError("budget exceeded")
s = objective(perm)
self.evals += 1
self.scores.append(s)
if s < self.best_score:
self.best_score = s
self.best_perm = perm
return s
def trace(self):
return np.minimum.accumulate(np.array(self.scores))The trace method returns the running best-so-far length after every evaluation, which is the incumbent curve we will plot. When we call it at the end of a run, it captures the entire trajectory rather than just the final number. That distinction matters: two methods can end at the same objective while one found it early and the other found it only at the last evaluation.
Key Methods
Hill climbing scans the reversals in the order shown, retaining a candidate through the tolerance-aware update described above. After a completed scan, it steps to that candidate only if the candidate beats the current length by more than ; otherwise it restarts from a fresh random permutation if budget remains. Equal-cost neighbors still consume budget and the scan continues; they do not themselves authorize an installed step. If the budget interrupts a partial scan, the function returns the best permutation scored so far without spending another evaluation.
Simulated annealing proposes a uniformly random reversal, accepts it by the Metropolis rule, and updates the running best. Nothing about the acceptance uses the budget structures differently from hill climbing: every accepted and every rejected proposal costs one evaluation. Random search uses the same budget object but proposes uniformly random permutations, which makes it a baseline for how much structure the local moves add.
The implementation also rejects a zero or subnormal ratio in the quotient's floating-point type. This range check leaves equal positive temperatures valid, even when their absolute values are subnormal. The executable schedule is a floating-point approximation to the ideal geometric sequence.
PythonIn [22]Expand
def reverse_segment(perm, i, j):
p = list(perm)
p[i : j + 1] = p[i : j + 1][::-1]
return tuple(p)
def hill_climb(limit, seed):
rng = np.random.default_rng(seed)
b = Budget(limit)
current = nearest_neighbor()
current_score = b.score(current)
while b.evals < limit:
best_local, best_local_score = None, np.inf
for i in range(8):
for j in range(i + 1, 8):
if b.evals >= limit:
return b.best_perm, b.best_score, b.trace()
n = reverse_segment(current, i, j)
s = b.score(n)
if best_local is None or s < best_local_score - 1e-12:
best_local, best_local_score = n, s
if best_local_score < current_score - 1e-12:
current, current_score = best_local, best_local_score
else:
if b.evals >= limit:
break
current = tuple(rng.permutation(np.arange(1, 9)))
current_score = b.score(current)
return b.best_perm, b.best_score, b.trace()
def simulated_annealing(limit, seed, T0=2.0, T1=0.01):
if not (np.isfinite(T0) and np.isfinite(T1) and 0 < T1 <= T0):
raise ValueError("temperatures must be finite with 0 < T1 <= T0")
temperature_ratio = T1 / T0
if temperature_ratio < np.finfo(np.asarray(temperature_ratio).dtype).tiny:
raise ValueError(
"temperature ratio must be normal in its floating-point dtype"
)
rng = np.random.default_rng(seed)
b = Budget(limit)
current = nearest_neighbor()
current_score = b.score(current)
while b.evals < limit:
progress = b.evals / (limit - 1)
temp = T0 * (T1 / T0) ** progress
i, j = sorted(rng.choice(8, size=2, replace=False))
n = reverse_segment(current, i, j)
s = b.score(n)
delta = s - current_score
if delta < 0 or rng.random() < np.exp(-delta / temp):
current, current_score = n, s
return b.best_perm, b.best_score, b.trace()
def random_search(limit, seed):
rng = np.random.default_rng(seed)
b = Budget(limit)
b.score(nearest_neighbor())
while b.evals < limit:
p = tuple(rng.permutation(np.arange(1, 9)))
b.score(p)
return b.best_perm, b.best_score, b.trace()Budgeted Comparisons and Repeated-Run Uncertainty
We now run each method on each of seeds with the same -evaluation budget, and record both the final incumbent and the full incumbent trace. We also check every reported final solution: its identifiers must form the required permutation, it must pass the feasibility checks, and its recomputed objective must match the reported best score within tolerance. These assertions are part of the experiment, not decoration; they catch accounting mistakes before those mistakes become a figure.
PythonIn [23]Expand
SEEDS = list(range(20))
BUDGET = 300
METHODS = {
"Hill climbing": hill_climb,
"Simulated annealing": simulated_annealing,
"NN-seeded random search": random_search,
}
results = {name: {"final": [], "traces": []} for name in METHODS}
for seed in SEEDS:
for name, fn in METHODS.items():
perm, best, trace = fn(BUDGET, seed)
assert is_feasible(perm), (
f"infeasible final perm for {name}, seed {seed}"
)
assert np.isclose(objective(perm), best, rtol=0, atol=1e-12)
assert len(trace) == BUDGET
assert np.all(np.isfinite(trace))
assert np.all(np.diff(trace) <= 0)
assert np.isclose(trace[-1], best, rtol=0, atol=1e-12)
results[name]["final"].append(best)
results[name]["traces"].append(trace)ConsoleOut [24]Expand
Method mean std min max Hill climbing 37.8378 0.4567 37.3560 38.3926 Simulated annealing 37.8173 0.6258 37.3560 38.9898 NN-seeded random search 39.2541 0.8956 38.0002 40.7200 Exact optimum: 37.3560
The recorded traces show how quickly each run finds its incumbent. Compare the full distributions and repeated-run summaries rather than claiming universal superiority from this tiny instance. NN-seeded random search reports the best of its starting score and 299 random draws, so its final values do not have the distribution of one uniformly random permutation. The fixed instance, neighborhood, initialization, and temperature schedule all limit how these results generalize.
PythonIn [25]Expand
trajectory_data = {}
for name in METHODS:
traces = np.asarray(results[name]["traces"])
trajectory_data[name] = (
traces.mean(axis=0),
*np.quantile(traces, [0.1, 0.9], axis=0),
)
plot_positions = np.arange(1, BUDGET + 1)
rng_jitter = np.random.default_rng(0)
final_plot_data = {
name: (
np.asarray(results[name]["final"]),
rng_jitter.uniform(-0.08, 0.08, len(SEEDS)),
)
for name in METHODS
}VisualizationOut [26]Expand

The final distribution is the honest summary of quality. It shows not only how good each method is on average, but how much its answer depends on the seed. On this instance the exact optimum is attained by some seeds and missed by others; the spread of a method's final values is the practical price of stochasticity, and it is what we should report when we compare search heuristics.
VisualizationOut [27]Expand

Key takeaways from the implementation:
- The methods receive the same number of complete objective evaluations; this does not equalize runtime, construction work, or memory.
- On these 20seeds, annealing has a slightly lower mean than hill climbing and a larger descriptive standard deviation. This is an observation about one instance and schedule.
- NN-seeded random search has a higher final mean here, while providing a baseline with the same nearest-neighbor initial score.
- Report seed variation and the experimental configuration alongside the mean.
Practical Implications
Local search and simulated annealing are appropriate when the decision variable is a discrete arrangement, the objective is cheap to evaluate or can be updated incrementally, and a good feasible solution is more useful than a proof of optimality. Routing, scheduling, packing, and assignment are natural fits because a small structural edit (reversing a segment, swapping two jobs, moving an item) produces a candidate that can be scored directly. If an exact method can handle the instance, include it as a reference, as we do by enumerating these 40,320 permutations. A disconnected neighborhood or a budget too small for a meaningful scan also limits the heuristic. In those settings the heuristic may still produce a useful answer, but the comparison should be framed as baseline improvement under a fixed budget, not as a claim of global optimality.
The representation, decoder, and neighborhood jointly determine which candidate changes the algorithm can explore. In this chapter the permutation, the fixed four-customer split, and inclusive reversals work together because adjacent reversals generate all permutations, so the state graph is connected, while wider reversals can cross the route boundary and reassign customers between vehicles. Restricting moves to within-route reversals would silently forbid that reassignment and could leave the search stuck on a solution that is poor only because a customer is on the wrong vehicle. The practical test is to ask whether the move set can express the changes the objective rewards, and whether repeated application can reach every feasible solution. Symmetric within-route reversals give routing a geometric neighborhood, but arbitrary reversals need not improve a route. Scheduling and packing require moves suited to their own decisions and constraints.
Data requirements are modest but not zero. The instance must supply the coordinates, demands, capacities, time windows, service times, and depot information needed by the decoder, and the distance matrix should be computed once and reused rather than recalculated at every proposal. A feasibility checker is part of the data pipeline, not an afterthought: it should be called on every reported final solution, and in this chapter the experiment checks the required permutation, feasibility, objective equality, and monotone incumbent traces. Unlike feature-based models, local search does not require feature scaling or missing-value imputation in the usual sense, but inconsistent units, asymmetric distances, or an infeasible warm start can dominate the outcome. When several objectives are combined, state their units, scaling, and policy weights before scalarizing; equal numerical scales alone do not determine appropriate priorities.
Check experimental comparisons at the level of scored candidates. Reporting only a successful seed hides run-to-run variation. Counting only accepted moves understates cost, because rejected proposals, restarts, and initialization all consume evaluations. Equal evaluation counts should be accompanied by runtime measurements if computational cost is the comparison target: scoring a permutation traverses every arc in both tours, while a nearest-neighbor scan examines candidate arcs more cheaply. Inspect uphill acceptance separately from ties and improvements. A high total acceptance rate can reflect equal-cost or improving proposals; it does not alone prove that temperature is too high. Very few accepted worsening moves may reflect low temperature or the proposed objective changes. Hill climbing should be paired with restarts or perturbation, and annealing should be monitored rather than tuned by a single final number.
Computational cost is governed by how often the objective is evaluated and how expensive each evaluation is. A full scan here examines 28 reversals, so the 300-evaluation budget supports roughly ten scans, enough to reveal plateaus but not enough to certify optimality. For a symmetric within-route 2-opt reversal, the distance delta uses two removed and two replacement boundary edges and can be computed in distance lookups. Cross-split moves require the correct decoder-specific depot accounting, and time-window feasibility can still need a route scan. This implementation deliberately uses full candidate scoring and copies permutations; its distance matrix uses storage, each budget trace uses storage, and the repeated-run experiment retains all traces. Independent seeds can be run separately. If the neighborhood grows quadratically with the number of customers, candidate lists or first-improvement scans may be necessary to keep each iteration affordable; the general principle is to preserve the evaluation accounting while making each evaluation as light as the problem structure allows.
Performance should be reported as a distribution, not a point. Across 20 seeds with the same 300-evaluation budget, report the mean, standard deviation, minimum, maximum, and incumbent traces for each method, and compare against both a construction heuristic and a random baseline. On small instances, full enumeration supplies an exact reference; on larger ones, use the best known solution or a lower bound. In deployment, log the acceptance rate, the incumbent improvement, and the feasibility rate to help diagnose poor schedule behavior. Inspect uphill acceptance separately from ties and improvements; these logs do not guarantee a unique diagnosis or detection before quality degrades. Re-run the search when the instance changes, keep the seed ensemble as part of the evaluation protocol, and treat the method as a budgeted improvement procedure rather than a source of guaranteed optima. That framing also sets up the next chapter, where population-based methods reuse this instance, decoder, and evaluation accounting to ask whether a different search dynamic can spend the same evaluations more effectively.
Summary
Local search and simulated annealing are the entry point to heuristic optimization of discrete, non-differentiable objectives. We defined a solution as a permutation of customers, decoded it into two depot tours, and measured it by total Euclidean arc length. We defined a neighborhood of all inclusive segment reversals, established that this state graph is connected, and used that connectivity to motivate cross-route reassignment. Hill climbing used ordered, tolerance-aware neighborhood scans. Every installed step strictly improved the current state within its epoch; the global incumbent retained the best candidate scored across all epochs. When a completed scan refused to step, a random restart spent the remaining budget on a fresh starting point.
Simulated annealing replaced strict improvement with a probabilistic acceptance rule for worsening moves, justified through the Metropolis-Hastings stationary distribution , and annealed the temperature geometrically from a hot start to a near-greedy end. Under a common -evaluation budget and seeds, hill climbing was low-variance and prone to plateaus, annealing had a slightly lower mean and larger descriptive standard deviation on this instance, and NN-seeded random search had a higher mean. Full enumeration of the permutations supplied the optimality yardstick and certified runs that matched it. The finite geometric schedule alone supplies no such certificate.
The larger lesson is about how to reason about heuristics: define the representation and neighborhood clearly, count every scored candidate, validate feasibility on every reported final solution, and report a distribution across seeds rather than a single lucky result. Those habits carry directly into the population-based methods that follow, which reuse this instance and budget to ask whether a fundamentally different search dynamic can spend the same evaluations more effectively.
Quiz
Ready to test your understanding? Take this quick quiz to reinforce what you've learned about local search and simulated annealing on the eight-customer routing instance.
Join the discussionAsk a question or share what you learned.
Keep exploring
Reference
Citation details
Cite or share this article.


Comments
Add to the conversation
A question, a correction, or an example from your own work.
No comments yet. Be the first to share your thoughts!