Simulated Annealing

Simulated annealing (SA) is a probabilistic global optimization metaheuristic in computation and algorithms designed to locate the global minimum of a cost or energy function E(x) over a large, rugged discrete or continuous search space Ω Containing numerous local minima. Introduced independently by Scott Kirkpatrick, C. Daniel Gelatt Jr., and Mario P. Vecchi at IBM Thomas J. Watson Research Center in Science in May 1983, and by Slovak physicist Vlado Černý in 1985, the algorithm adapts the statistical-mechanics Metropolis–Hastings algorithm by replacing its fixed temperature parameter with a monotonically decreasing cooling schedule T₀ ≥ T₁ ≥ T₂ ≥ ··· → 0.

Whereas a deterministic greedy algorithm or steepest-descent local search accepts only cost-reducing transitions (ΔE ≤ 0) and becomes permanently trapped in the first local minimum basin it encounters, simulated annealing accepts cost-increasing uphill moves (ΔE > 0) with Boltzmann probability exp(−ΔE / T_k). At high initial temperatures T_k ≫ ΔE, the walker traverses energy barriers freely across the entire state space; as T_k is gradually lowered, uphill transitions become progressively rarer until the trajectory settles into deep, low-energy basins of attraction.

The name and mathematical structure of the algorithm derive from the physical metallurgy of annealing. To grow a defect-free single crystal with minimum lattice potential energy, a solid is first heated above its recrystallization temperature so that atoms undergo random thermal motion across all configurational microstates. At thermal equilibrium at thermodynamic temperature T, the probability that the system occupies microstate x with energy E(x) is governed by the canonical Boltzmann distribution:

π_T(x) = (1 / Z(T)) · exp(−E(x) / (k_BT)), where Z(T) = ∑_y∈Ω exp(−E(y) / (k_BT))

As temperature approaches zero (T → 0⁺), the Boltzmann distribution π_T(x) concentrates its probability mass uniformly onto the set of global minimum-energy ground states Ω_min = { x* ∈ Ω : E(x*) = E_min }, because for any suboptimal state x with E(x) > E_min, the Boltzmann ratio exp(−(E(x) − E_min) / (k_BT)) decays to zero.

• Rapid Quenching (Greedy Local Search, T → 0 Abruptly): If a molten metal is plunged directly into a cold bath, atoms lose thermal mobility before finding long-range crystalline order, freezing into an amorphous glass or polycrystalline solid riddled with high-energy grain-boundary defects—the physical analog of a local optimization algorithm trapped in a poor local minimum.
• Controlled Annealing (Quasi-Equilibrium Cooling): If temperature is reduced slowly enough that the system re-establishes quasi-thermal equilibrium at each intermediate temperature stage, thermal fluctuations allow atoms to hop out of shallow metastable traps and lock into the global ground-state lattice.

To optimize an arbitrary mathematical objective function E : Ω → R without physical dimensions, the Boltzmann constant k_B is absorbed into the dimensionless control parameter computational temperature T > 0. Simulated annealing requires three design components:

1. State Space (Ω) and Objective Function (E(x)): A representation of candidate solutions x ∈ Ω and a scalar cost E(x) to be minimized.
2. Neighborhood Generator (q(x′ | x)): A stochastic perturbation rule that proposes a neighboring state x′ ∈ N(x) from the current state x (such as swapping two cities in a tour, flipping a spin in an Ising lattice, or perturbing a torsion angle in a molecule).
3. Cooling Schedule (T_k): A sequence of temperatures specifying the initial temperature T₀, the temperature decay law T_k, the number of Monte Carlo method steps L_k executed at each temperature plateau, and a stopping criterion.

Step-by-Step Transition Rule

At iteration k with current state x_k = x and temperature T_k:

1. Sample a candidate neighbor x′ ~ q(· | x) and evaluate the energy difference ΔE = E(x′) − E(x).
2. Compute the Metropolis acceptance probability:

P_accept(x → x′; T_k) = 1 if ΔE ≤ 0 P_accept(x → x′; T_k) = exp(−ΔE / T_k) if ΔE > 0

1. If ΔE ≤ 0, accept the downhill move unconditionally (x_k+1 ← x′). If ΔE > 0, draw a uniform random number u ~ Uniform(0, 1); if u < exp(−ΔE / T_k), accept the uphill move (x_k+1 ← x′); otherwise reject the move and retain x_k+1 ← x.

At any fixed temperature T_k = T with a symmetric proposal distribution q(x′ | x) = q(x | x′), this transition kernel defines a homogeneous Markov chain satisfying detailed balance π_T(x)P(x → x′) = π_T(x′)P(x′ → x) with stationary distribution π_T(x) ∝ exp(−E(x) / T). When T_k changes with k, the process becomes an inhomogeneous Markov chain.

Hajek's Necessary and Sufficient Condition (1988)

A central theoretical milestone in computational complexity theory and stochastic optimization is Bruce Hajek's 1988 theorem, which established the cooling rate required for simulated annealing to converge to the global optimum with probability 1 on any finite state space Ω.

For any local minimum m ∉ Ω_min that is not a global minimum, define its basin depth d(m) as the smallest maximum energy elevation ΔE along any path of transitions connecting m to a state of lower energy E(y) < E(m). Let d* = max_m∉Ωmin d(m) denote the depth of the deepest non-global local minimum basin in Ω. Hajek proved that under a logarithmic cooling schedule:

T(k) = c / ln(1 + k)

the inhomogeneous Markov chain satisfies lim_k→∞ Pr(x_k ∈ Ω_min) = 1 if and only if c ≥ d* (equivalently, ∑_k=1^∞ exp(−d* / T(k)) = ∞). Intuitively, the harmonic divergence of ∑ (1 + k)^−d*/c for c ≥ d* ensures that the walker retains enough cumulative probability across infinite time to climb over the tallest local barrier d* infinitely often until it enters the global minimum basin.

Comparison of Theoretical and Practical Cooling Schedules

Because logarithmic cooling T(k) = c / ln(1 + k) decreases too slowly for finite computational budgets, practical engineering implementations employ faster finite-time schedules:

Traveling Salesperson Problem (TSP) and VLSI Electronic Design Automation

In Kirkpatrick, Gelatt, and Vecchi's 1983 paper, simulated annealing was applied to two NP-hard combinatorial problems:

1. Traveling Salesperson Problem (TSP): Given N cities with pairwise distances d(i, j), a state x is a permutation tour of {1, ..., N} with length E(x) = ∑ d(x_i, x_i+1). Using Croes's 2-opt move (reversing a contiguous tour segment between two indices a < b, which changes only 2 edges in O(1) time), simulated annealing escapes crossed-edge local minima and consistently finds tours within 0.5% to 1.5% of the Held–Karp lower bound on instances with tens of thousands of cities.
2. VLSI Chip Placement and Wire Routing: When placing thousands of logic macrocells on a silicon wafer (such as the IBM 3081 processor designed by Kirkpatrick's team and Carl Sechen's TimberWolf placement suite), the energy function E(x) combines total Manhattan wirelength with cell-overlap penalties and routing-channel congestion. Simulated annealing reduced wire congestion by 15–30% over human layout engineers and became a standard industrial engine in electronic design automation (EDA).

Crystallographic Refinement and Biomolecular Conformational Search

In structural biology and computational chemistry, Axel T. Brünger introduced simulated annealing refinement (implemented in X-PLOR and CNS) for X-ray crystallography and nuclear magnetic resonance (NMR) spectroscopy in 1987. By coupling a physical force field E_chem with an experimental diffraction residual w_X-rayE_X-ray and heating a protein atomic model to T = 3,000–4,000 K in molecular dynamics before slow cooling to 300 K, simulated annealing flips misoriented peptide side chains across 10–30 kJ/mol steric barriers without manual human rebuilding. Similarly, small-molecule molecular docking packages (such as AutoDock) and de novo protein design suites (Rosetta fragment assembly) use simulated annealing over backbone torsion angles (φ, ψ, χ) to sample low-free-energy binding poses and tertiary folds.

How should the initial temperature T₀ be chosen in simulated annealing?

A standard empirical calibration runs a short preliminary random walk of 100 to 500 trial moves, computes the average positive uphill energy change ⟨ΔE₊⟩, and sets T₀ = −⟨ΔE₊⟩ / ln(p₀) so that the initial uphill acceptance rate p₀ = exp(−⟨ΔE₊⟩ / T₀) is approximately 80% to 90% (p₀ ≈ 0.85), ensuring the chain is initially melted and unrestrained by local barriers.

Why should cooling slow down when the energy variance (specific heat) is large?

In statistical physics, the dimensionless specific heat at temperature T is proportional to the energy variance of the chain: C(T) = Var_T(E) / T². A sharp spike in C(T) signals a configurational phase transition where the system is freezing into its macroscopic structural arrangement (such as global tour topology in TSP or secondary/tertiary core packing in a protein); cooling too rapidly through this critical temperature window locks in large-scale topological defects that lower temperatures cannot undo.

1. Kirkpatrick, S., Gelatt, C. D., Jr., & Vecchi, M. P. (1983). Optimization by simulated annealing. Science, 220(4598), 671–680. https://doi.org/10.1126/science.220.4598.671
2. Černý, V. (1985). Thermodynamical approach to the traveling salesman problem: An efficient simulation algorithm. Journal of Optimization Theory and Applications, 45(1), 41–51. https://doi.org/10.1007/BF00940812
3. Hajek, B. (1988). Cooling schedules for optimal annealing. Mathematics of Operations Research, 13(2), 311–329. https://doi.org/10.1287/moor.13.2.311
4. Brünger, A. T., Kuriyan, J., & Karplus, M. (1987). Crystallographic R factor refinement by molecular dynamics. Science, 235(4787), 458–460. https://doi.org/10.1126/science.235.4787.458