Simulated Annealing: Theory and Python Implementation

Introduction

Simulated Annealing (SA) is a metaheuristic optimization technique inspired by the annealing process in metallurgy. Just as slow cooling allows a metal’s crystal structure to reach a minimum energy state, SA searches for the global optimum of an objective function.

The key feature of SA is probabilistic acceptance of worsening solutions, enabling escape from local optima. Early iterations perform bold exploration at high temperature, transitioning to fine-grained local search as temperature decreases.

Algorithm

Metropolis Criterion

The core of SA is the Metropolis criterion. Transition from current solution \(x\) to neighbor \(x'\) is accepted with probability:

where \(\Delta E = f(x') - f(x)\) and \(T\) is the current temperature.

Deriving the Metropolis Criterion as a Special Case of MH

The acceptance rule (1) is not an arbitrary formula handed down from nowhere. It is precisely the symmetric-proposal special case of the Metropolis-Hastings (MH) acceptance probability we derived in a companion article, Markov Chain Monte Carlo (MCMC): Metropolis-Hastings and Gibbs Sampling . Here we make that connection explicit in the optimization context of this article.

Fix a temperature \(T\) , treat the objective value \(E(x) = f(x)\) as an energy, and regard the following Boltzmann distribution as the “target distribution”:

The key MH property that the normalizing constant \(Z(T)\) (the partition function) never needs to be computed explicitly (see the MCMC article) carries over here as well. Applying the MH acceptance probability (Eq. (4) in the MCMC article) to target \(\pi_T\) with a symmetric neighbor-proposal distribution \(q(x' | x) = q(x | x')\) (as is the case for both a Gaussian perturbation and a 2-opt swap: the probability of proposing \(x \to x'\) equals that of proposing \(x' \to x\) ) gives

which, setting \(\Delta E = E(x') - E(x) = f(x') - f(x)\) , is exactly the Metropolis criterion (1). Simulated annealing is therefore best understood as “repeatedly applying the Metropolis method targeting the Boltzmann distribution \(\pi_T\) at the current temperature \(T\) , while slowly lowering \(T\) .” Because the partition function \(Z(T)\) cancels in the ratio, we can avoid ever evaluating the intractable integral \(Z(T)\) as long as \(E(x)\) itself can be computed — the same benefit MH provides.

If the proposal is asymmetric (e.g., a neighborhood generator whose forward and backward proposal probabilities differ), the Hastings correction term \(q(x | x') / q(x' | x)\) from Eq. (3) of the MCMC article must be added to (1). Both Python implementations in this article (Gaussian perturbation, 2-opt swap) use symmetric proposals, so this correction is unnecessary here — but that is a consequence of that specific design choice, not something that can be omitted in general.

Boltzmann Distribution Connection

With sufficient exploration at temperature \(T\) , the solution distribution converges to the Boltzmann distribution:

As \(T \to 0\) , the distribution concentrates on the optimal solution.

Proof via Detailed Balance

Convergence to the Boltzmann distribution in (2) follows directly from the general theorem proved in the MCMC article: a transition kernel satisfying detailed balance has the target distribution as its stationary distribution. At a fixed temperature \(T\) , the transition kernel \(T_{\text{SA}}(x' | x) = q(x' | x) \alpha(x' | x)\) (for \(x' \neq x\) ) induced by the Metropolis criterion (1) satisfies detailed balance

which follows by substituting \(\pi = \pi_T\) into the proof given in the MCMC article (the case analysis with ratio \(r = \pi_T(x') q(x | x') / (\pi_T(x) q(x' | x))\) and \(\alpha = \min(1, r)\) , \(\alpha = \min(1, 1/r)\) ); for symmetric proposals the \(q\) terms cancel and the proof simplifies further. Hence, if we fix \(T\) and run simulated annealing long enough, the distribution of solutions converges to the Boltzmann distribution \(\pi_T\) in (2).

One important caveat: this convergence argument only holds for a fixed temperature. Real simulated annealing changes the temperature at every step, so the transition kernel \(T_{\text{SA}}\) itself varies with time — making the process a time-inhomogeneous Markov chain. The ergodic theorem introduced in the MCMC article assumes a fixed (time-homogeneous) kernel and therefore does not directly apply. We take up this point in the next section.

SA and MCMC: A Time-Inhomogeneous Chain That Degenerates into Greedy Search

SA Is MCMC with a Continuously Decreasing Temperature

As shown above, each step of simulated annealing is exactly “one step of the Metropolis method targeting the Boltzmann distribution \(\pi_{T_k}\) at the current temperature \(T_k\) .” Repeating this while following a cooling schedule \(T_0 > T_1 > T_2 > \cdots\) makes the whole process an MCMC whose target distribution keeps changing over time. This non-stationarity is exactly why SA is used as an optimization tool rather than a sampling tool, unlike ordinary MCMC. If \(T\) were held fixed, we would eventually just get samples from \(\pi_T\) ; but because \(\pi_T\) itself collapses onto the optimum as \(T \to 0\) , the chain tracks that shrinking target and is pulled toward the optimal solution.

Deriving the Degeneration into Greedy Search as \(T \to 0\)

The Boltzmann distribution \(\pi_T(x) \propto \exp(-E(x)/T)\) concentrates entirely on the minimum-energy state as \(T \to 0\) . To see this, let \(x^*\) be the optimal solution (\(E(x^*) = E_{\min}\) ) and \(x\) be any other state (\(E(x) > E_{\min}\) ). Their relative probability is

(since \(E(x) - E_{\min} > 0\) is a fixed constant, the exponent diverges to \(-\infty\) as \(T \to 0\) ), so the relative probability of any non-optimal state collapses to zero.

Applying the same limit to the acceptance probability (1), for \(\Delta E > 0\) ,

so worsening moves are never accepted in the limit. For \(\Delta E \leq 0\) , (1) always accepts regardless of \(T\) . Hence the Metropolis criterion degenerates as \(T \to 0\) into

which is exactly greedy search (hill climbing, accepting only improving moves). In other words, simulated annealing smoothly interpolates, along its cooling schedule, between two limiting behaviors: an essentially uniform random walk at high temperature, and greedy hill climbing at low temperature.

Verification: Confirming SA Matches Greedy Search at Sufficiently Low Temperature

We verified this limiting behavior directly on TSP. Using the same initial tour and the same sequence of neighbor proposals (fixing the random seed and pre-generating the \((i, j)\) proposals together with the uniform random numbers \(u\) used for acceptance decisions, so both variants share identical randomness), we ran (a) SA following the Metropolis criterion at a fixed constant temperature \(T\) , and (b) greedy search (hill climbing) accepting only improving proposals, on the exact same sequence of 20,000 proposals, and compared whether the accept/reject decision at every step was identical.

import numpy as np

np.random.seed(42)
n_cities = 30
cities = np.random.rand(n_cities, 2) * 100


def tour_length(tour, cities):
    n = len(tour)
    d = 0.0
    for i in range(n):
        d += np.linalg.norm(cities[tour[i]] - cities[tour[(i + 1) % n]])
    return d


def two_opt_delta(tour, cities, i, j):
    """Distance change from a 2-opt swap (reversing tour[i:j+1]), via delta computation"""
    n = len(tour)
    a, b = tour[i - 1], tour[i]
    c, d = tour[j], tour[(j + 1) % n]
    old = np.linalg.norm(cities[a] - cities[b]) + np.linalg.norm(cities[c] - cities[d])
    new = np.linalg.norm(cities[a] - cities[c]) + np.linalg.norm(cities[b] - cities[d])
    return new - old


def two_opt_swap(tour, i, j):
    new_tour = tour.copy()
    new_tour[i : j + 1] = tour[i : j + 1][::-1]
    return new_tour


def run_with_rule(cities, init_tour, proposals, accept_rule):
    """proposals: shared (i, j, u, T) sequence; only accept_rule differs between runs"""
    tour = init_tour.copy()
    current = tour_length(tour, cities)
    best_dist = current
    decisions = []
    for i, j, u, T in proposals:
        delta = two_opt_delta(tour, cities, i, j)
        accept = accept_rule(delta, T, u)
        decisions.append(accept)
        if accept:
            tour = two_opt_swap(tour, i, j)
            current += delta
            best_dist = min(best_dist, current)
    return best_dist, decisions


# --- Pre-generate a shared proposal sequence ---
max_iter = 20000
rng = np.random.RandomState(7)
init_tour = list(rng.permutation(n_cities))
proposals_template = []
for _ in range(max_iter):
    i, j = sorted(rng.choice(n_cities, 2, replace=False))
    if i == 0:
        i = 1
    if j <= i:
        j = min(i + 1, n_cities - 1)
    u = rng.random()
    proposals_template.append((i, j, u))

# --- SA at fixed temperature T vs. greedy search (identical proposal sequence) ---
for T_const in [1.0, 1e-3, 1e-6, 1e-9]:
    proposals = [(i, j, u, T_const) for (i, j, u) in proposals_template]

    def sa_rule(delta, T, u):
        return delta <= 0 or u < np.exp(-delta / T)

    def greedy_rule(delta, T, u):
        return delta <= 0

    sa_dist, sa_decisions = run_with_rule(cities, init_tour, proposals, sa_rule)
    greedy_dist, greedy_decisions = run_with_rule(cities, init_tour, proposals, greedy_rule)
    mismatches = sum(a != b for a, b in zip(sa_decisions, greedy_decisions))

    print(
        f"T={T_const:.0e}  SA={sa_dist:.4f}  greedy={greedy_dist:.4f}  mismatches={mismatches}/{max_iter}"
    )

The results were as follows.

At \(T = 1.0\) , 120 (0.6%) worsening proposals were stochastically accepted, and SA reached a better solution (453.9023) than greedy search (466.6980) — a concrete example of accepting worsening moves helping escape a local optimum. At \(T = 10^{-3}\) and below, all 20,000 accept/reject decisions matched greedy search exactly, and the final distances agreed to four decimal places. Given the typical magnitude of \(\Delta E\) in this TSP instance (on the order of single digits to tens), \(\exp(-\Delta E / T)\) becomes numerically indistinguishable from zero once \(T\) drops to around \(10^{-3}\) , confirming that SA degenerates into greedy search exactly as the theory predicts.

Temperature Schedules

Geometric Cooling (Most Common)

\(\alpha = 0.95\) is a typical value.

Linear Cooling

Logarithmic Cooling (Theoretical Guarantee)

Theoretically guarantees convergence to the global optimum, but convergence is extremely slow and impractical.

Why Only Logarithmic Cooling Guarantees Convergence

The claim that logarithmic cooling (5) theoretically guarantees convergence to the global optimum is not something to accept on faith — it can be understood via the following intuitive probabilistic argument.

Escaping a non-global local optimum on the energy landscape requires accepting at least one worsening transition that climbs an energy barrier of height \(d\) . By the Metropolis criterion (1), the probability that such a transition is accepted at step \(k\) is roughly of order \(\exp(-d/T_k)\) , where \(d\) is the largest energy gap that must be climbed along the escape path. Substituting the logarithmic schedule \(T_k = c/\ln(k+1)\) (with constant \(c\) ) gives

a power law. The infinite series \(\sum_{k=1}^{\infty} (k+1)^{-d/c}\) diverges (harmonic-series-like) when the exponent satisfies \(d/c \leq 1\) (i.e., \(c \geq d\) ), and converges when \(d/c > 1\) (i.e., \(c < d\) ). By the intuition behind the Borel-Cantelli lemma (in its simplified, independence-assuming form), if the series corresponding to the expected number of occurrences of “an event with probability \((k+1)^{-d/c}\) at each time step” diverges, that event (accepting a worsening move that clears barrier \(d\) ) is guaranteed to occur infinitely often over time. In other words: no matter how deep the energy barrier \(d\) is, choosing \(c \geq d\) guarantees that an escape opportunity will eventually arrive.

This intuition was made rigorous by Geman & Geman (1984) in the context of stochastic image restoration. Their theorem states that if the constant \(c\) in the logarithmic schedule \(T_k = c/\ln(k+1)\) is chosen at least as large as the maximum energy barrier needed to escape any non-global local optimum, then as \(k \to \infty\) the distribution of solutions converges to the global optimum with probability 1. The minimal constant giving this sufficient condition (which equals the maximum “depth” among all non-global local optima) was later characterized more precisely by Hajek (1988), who showed the condition to be necessary and sufficient.

By contrast, for geometric cooling \(T_k = T_0 \alpha^k\) (\(0 < \alpha < 1\) ), the same calculation gives

Since \(\alpha^{-k}\) grows exponentially in \(k\) , this probability collapses to zero doubly-exponentially fast — far faster than under logarithmic cooling. The series \(\sum_k \exp(-d \alpha^{-k}/T_0)\) converges to a finite value for any \(d > 0\) , so the Borel-Cantelli-type argument runs in reverse: once trapped in a deep local optimum, the expected number of future escape attempts saturates quickly, and escape with probability 1 is no longer guaranteed. This is precisely why geometric and linear cooling cannot theoretically guarantee convergence to the global optimum.

This does not mean such schedules are practically useless, however. As the numerical experiment below shows, geometric cooling reaches far higher-quality solutions than logarithmic cooling within a limited iteration budget. Logarithmic cooling’s theoretical guarantee is only meaningful under the unrealistic assumption of unlimited iterations; within any finite compute budget, cooling faster and spending more iterations on fine-grained low-temperature local search is practically superior — a central design trade-off in SA.

Setting the Initial Temperature and Cooling Rate: Edge Cases

SA’s practical performance depends far more on the choice of initial temperature \(T_0\) and cooling rate than on its theoretical convergence guarantee. Here we examine typical failure modes and quantify their behavior experimentally.

Initial temperature too high: nearly all worsening proposals are accepted, so the early phase is essentially a random walk. Little use is made of the objective function’s information, and if the iteration budget runs out before cooling catches up, the final solution quality can end up worse.

Initial temperature too low: worsening proposals are almost never accepted, so the search effectively becomes greedy search (see previous section) from the very start. It converges to a local optimum immediately, with almost no opportunity to escape.

To verify this, we varied the initial temperature \(T_0\) from \(0.01\) to \(10^5\) on TSP (30 cities, geometric cooling with \(\alpha=0.9995\) , 20,000 iterations, averaged over 5 seeds), measuring the overall acceptance rate, the acceptance rate in the early (first 10%), middle (40-50%), and late (last 10%) phases of the run, and the final tour length.

import numpy as np

np.random.seed(42)
n_cities = 30
cities = np.random.rand(n_cities, 2) * 100


def tour_length(tour, cities):
    n = len(tour)
    d = 0.0
    for i in range(n):
        d += np.linalg.norm(cities[tour[i]] - cities[tour[(i + 1) % n]])
    return d


def two_opt_delta(tour, cities, i, j):
    n = len(tour)
    a, b = tour[i - 1], tour[i]
    c, d = tour[j], tour[(j + 1) % n]
    old = np.linalg.norm(cities[a] - cities[b]) + np.linalg.norm(cities[c] - cities[d])
    new = np.linalg.norm(cities[a] - cities[c]) + np.linalg.norm(cities[b] - cities[d])
    return new - old


def two_opt_swap(tour, i, j):
    new_tour = tour.copy()
    new_tour[i : j + 1] = tour[i : j + 1][::-1]
    return new_tour


def sa_tsp_track_acceptance(cities, T0, alpha, max_iter, seed, n_bins=10):
    rng = np.random.RandomState(seed)
    n = len(cities)
    tour = list(rng.permutation(n))
    current = tour_length(tour, cities)
    best_dist = current
    T = T0
    bin_size = max_iter // n_bins
    accept_counts = np.zeros(n_bins)
    propose_counts = np.zeros(n_bins)

    for k in range(max_iter):
        i, j = sorted(rng.choice(n, 2, replace=False))
        if i == 0:
            i = 1
        if j <= i:
            continue
        delta = two_opt_delta(tour, cities, i, j)
        b = min(k // bin_size, n_bins - 1)
        propose_counts[b] += 1
        if delta <= 0 or rng.random() < np.exp(-delta / T):
            tour = two_opt_swap(tour, i, j)
            current += delta
            accept_counts[b] += 1
            best_dist = min(best_dist, current)
        T *= alpha

    acc_rate_per_bin = accept_counts / np.maximum(propose_counts, 1)
    overall_acc = accept_counts.sum() / propose_counts.sum()
    return best_dist, overall_acc, acc_rate_per_bin


max_iter = 20000
alpha = 0.9995
for T0 in [0.01, 1, 10, 100, 1000, 10000, 100000]:
    dists, overalls, bins_all = [], [], []
    for seed in range(5):
        best_dist, overall_acc, acc_bins = sa_tsp_track_acceptance(cities, T0, alpha, max_iter, seed)
        dists.append(best_dist)
        overalls.append(overall_acc)
        bins_all.append(acc_bins)
    bins_mean = np.mean(bins_all, axis=0)
    print(
        f"T0={T0:g}  overall_acc={np.mean(overalls)*100:.2f}%  "
        f"final_dist={np.mean(dists):.2f}  "
        f"bins(0-10%,40-50%,90-100%)={bins_mean[0]*100:.1f}%,{bins_mean[4]*100:.1f}%,{bins_mean[9]*100:.1f}%"
    )

The results (averaged over 5 seeds) were as follows.

At \(T_0=0.01\) and \(1\) , the early acceptance rate is already down to only 3-4%; the search behaves essentially greedily from the start, and the final distances (471-481) are clearly worse than the mid-range \(T_0\) values (10-1,000, giving 454-457). Conversely, at \(T_0=10^5\) , the acceptance rate is 100% in the early phase (0-10%) and still 98.3% at the midpoint (40-50%) — even after consuming half of the 20,000 iterations, the search has barely escaped a nearly pure random walk, and the final distance (507.92) was the worst of all settings. The best result was obtained at \(T_0=10\) (453.98), though \(T_0=100\) -\(1{,}000\) (456-457) performed comparably well, confirming that a mid-range \(T_0\) that keeps the early acceptance rate around 70-97% within the iteration budget is a practically reasonable choice. Note that this “keep the initial acceptance rate reasonably high” heuristic is a different criterion from the optimal acceptance rate for Random-Walk Metropolis discussed in the MCMC article (about 44%, converging toward about 23% in high dimensions): that criterion maximizes estimation efficiency from a stationary distribution, whereas SA’s initial temperature is tuned for the distinct goal of exploring a sufficiently wide region early in the search.

Cooling rate trade-off: making the cooling rate (e.g., \(\alpha\) in geometric cooling, or \(c\) in logarithmic cooling) smaller (faster cooling) lets the search transition into a local-search phase with fewer iterations, but raises the risk of converging to a local optimum before global exploration is complete. Slower cooling improves solution quality, as confirmed in the numerical experiment below, but at the cost of fewer iterations left for fine-grained low-temperature local search within the same total budget. This “slow cooling = high quality, low speed” versus “fast cooling = low quality, high speed” trade-off is a central consideration in SA hyperparameter design.

Reheating strategies: periodically, or whenever no improvement has been observed for a fixed number of iterations, the temperature is deliberately raised again and the search restarted. Under a standard monotone cooling schedule, once the chain falls into a deep local optimum at low temperature, the escape probability \(\exp(-d/T)\) becomes essentially pinned at zero, and the search is effectively stuck in greedy mode thereafter. Reheating addresses this by deliberately raising the temperature to temporarily restore the escape probability, then cooling again — allowing the search to reach solutions unreachable within a single monotone cooling cycle. In practice this is implemented either by repeating geometric or logarithmic cooling over multiple cycles (resetting to \(T_0\) at the start of each cycle), or by multiplying the temperature by some factor whenever no improvement has been seen for a fixed number of iterations.

Python Implementation: Continuous Optimization

Optimizing the Rastrigin function (multimodal benchmark) with SA:

import numpy as np
import matplotlib.pyplot as plt

def rastrigin(x):
    """Rastrigin function"""
    return 10 * len(x) + np.sum(x**2 - 10 * np.cos(2 * np.pi * x))

def simulated_annealing(objective, dim, bounds=(-5.12, 5.12),
                        T_init=100.0, alpha=0.995, max_iter=10000):
    """Simulated Annealing for continuous optimization"""
    x = np.random.uniform(*bounds, size=dim)
    best_x = x.copy()
    best_f = objective(x)
    current_f = best_f

    T = T_init
    history = [best_f]

    for i in range(max_iter):
        sigma = 0.5 * (T / T_init)
        x_new = x + np.random.normal(0, sigma, size=dim)
        x_new = np.clip(x_new, *bounds)

        new_f = objective(x_new)
        delta = new_f - current_f

        if delta <= 0 or np.random.random() < np.exp(-delta / T):
            x = x_new
            current_f = new_f
            if current_f < best_f:
                best_x = x.copy()
                best_f = current_f

        T *= alpha
        history.append(best_f)

    return best_x, best_f, history

# --- Run ---
np.random.seed(42)
best_x, best_f, history = simulated_annealing(rastrigin, dim=10)
print(f"Best fitness: {best_f:.6f}")

# --- Convergence plot ---
plt.figure(figsize=(10, 5))
plt.plot(history)
plt.xlabel('Iteration')
plt.ylabel('Best Objective Value')
plt.title('Simulated Annealing on Rastrigin Function (10D)')
plt.yscale('log')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

Python Implementation: Traveling Salesman Problem (TSP)

SA is effective for combinatorial optimization. Here we solve TSP using 2-opt neighborhood:

import numpy as np
import matplotlib.pyplot as plt

def total_distance(tour, cities):
    """Compute total tour distance"""
    n = len(tour)
    dist = 0
    for i in range(n):
        dist += np.linalg.norm(cities[tour[i]] - cities[tour[(i + 1) % n]])
    return dist

def two_opt_swap(tour, i, j):
    """2-opt swap: reverse tour[i:j+1]"""
    new_tour = tour.copy()
    new_tour[i:j+1] = tour[i:j+1][::-1]
    return new_tour

def sa_tsp(cities, T_init=1000.0, alpha=0.9995, max_iter=100000):
    """Simulated Annealing for TSP"""
    n = len(cities)
    tour = list(range(n))
    np.random.shuffle(tour)

    current_dist = total_distance(tour, cities)
    best_tour = tour.copy()
    best_dist = current_dist
    T = T_init
    history = [best_dist]

    for _ in range(max_iter):
        i, j = sorted(np.random.choice(n, 2, replace=False))
        new_tour = two_opt_swap(tour, i, j)
        new_dist = total_distance(new_tour, cities)

        delta = new_dist - current_dist
        if delta <= 0 or np.random.random() < np.exp(-delta / T):
            tour = new_tour
            current_dist = new_dist
            if current_dist < best_dist:
                best_tour = tour.copy()
                best_dist = current_dist

        T *= alpha
        history.append(best_dist)

    return best_tour, best_dist, history

# --- Run ---
np.random.seed(42)
n_cities = 30
cities = np.random.rand(n_cities, 2) * 100

best_tour, best_dist, history = sa_tsp(cities)
print(f"Best distance: {best_dist:.2f}")

# --- Visualization ---
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

tour_cities = cities[best_tour + [best_tour[0]]]
axes[0].plot(tour_cities[:, 0], tour_cities[:, 1], 'b-o', markersize=5)
axes[0].set_title(f'Best Tour (distance={best_dist:.2f})')
axes[0].set_aspect('equal')
axes[0].grid(True, alpha=0.3)

axes[1].plot(history)
axes[1].set_xlabel('Iteration')
axes[1].set_ylabel('Best Distance')
axes[1].set_title('SA Convergence on TSP')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Numerical Experiment: Comparing Cooling Schedules (TSP)

To quantify the trade-off between theoretical convergence guarantees and practical performance, we compared the logarithmic, geometric, and linear cooling schedules introduced above on the same TSP instance (30 cities, random seed 42), the same number of iterations (100,000), and the same initial temperature (\(T_0=1000\) ). Logarithmic cooling was normalized to match \(T_0\) as \(T_k = \frac{T_0 \ln 2}{\ln(k+2)}\) (choosing the constant so that \(T_0\) is recovered at \(k=0\) ), and linear cooling used \(\delta = T_0 / 100000\) so that the temperature reaches approximately zero exactly at the end of the run. Each schedule was run with 5 random seeds, recording the mean and standard deviation of the final tour length, and the best tour length at the 10%, 50%, and 100% iteration checkpoints.

import numpy as np

np.random.seed(42)
n_cities = 30
cities = np.random.rand(n_cities, 2) * 100


def tour_length(tour, cities):
    n = len(tour)
    d = 0.0
    for i in range(n):
        d += np.linalg.norm(cities[tour[i]] - cities[tour[(i + 1) % n]])
    return d


def two_opt_delta(tour, cities, i, j):
    n = len(tour)
    a, b = tour[i - 1], tour[i]
    c, d = tour[j], tour[(j + 1) % n]
    old = np.linalg.norm(cities[a] - cities[b]) + np.linalg.norm(cities[c] - cities[d])
    new = np.linalg.norm(cities[a] - cities[c]) + np.linalg.norm(cities[b] - cities[d])
    return new - old


def two_opt_swap(tour, i, j):
    new_tour = tour.copy()
    new_tour[i : j + 1] = tour[i : j + 1][::-1]
    return new_tour


def sa_tsp_schedule(cities, T0, schedule, max_iter, seed, T_min=1e-6):
    rng = np.random.RandomState(seed)
    n = len(cities)
    tour = list(rng.permutation(n))
    current = tour_length(tour, cities)
    best_dist = current
    checkpoints = {}
    cp_iters = sorted(set(int(max_iter * f) for f in (0.1, 0.5, 1.0)))

    for k in range(max_iter):
        if schedule == "log":
            c = T0 * np.log(2)
            T = max(c / np.log(k + 2), T_min)
        elif schedule == "exp":
            T = max(T0 * 0.9995**k, T_min)
        elif schedule == "linear":
            T = max(T0 - (T0 / max_iter) * k, T_min)

        i, j = sorted(rng.choice(n, 2, replace=False))
        if i == 0:
            i = 1
        if j <= i:
            continue
        delta = two_opt_delta(tour, cities, i, j)
        if delta <= 0 or rng.random() < np.exp(-delta / T):
            tour = two_opt_swap(tour, i, j)
            current += delta
            best_dist = min(best_dist, current)

        if (k + 1) in cp_iters:
            checkpoints[k + 1] = best_dist

    return best_dist, checkpoints


T0 = 1000.0
max_iter = 100000
for schedule in ["log", "exp", "linear"]:
    finals, cps_list = [], []
    for seed in range(5):
        best_dist, cps = sa_tsp_schedule(cities, T0, schedule, max_iter, seed)
        finals.append(best_dist)
        cps_list.append(cps)
    finals = np.array(finals)
    iters = sorted(cps_list[0].keys())
    cp_means = {it: np.mean([cps_list[s][it] for s in range(5)]) for it in iters}
    print(f"{schedule}: final={finals.mean():.3f}±{finals.std():.3f}  checkpoints={cp_means}")

The results (averaged over 5 seeds) were as follows.

Under the same budget of 100,000 iterations, geometric cooling was overwhelmingly the best, at 454.84, and also the most stable (standard deviation of only 1.87). Logarithmic cooling’s temperature at the end of the 100,000 iterations was still \(T = \frac{1000 \ln 2}{\ln(100002)} \approx 60.2\) — far from having cooled down — and its improvement from the 50% mark (939.34) to the 100% mark (906.02) was minimal, indicating the search is still effectively mid-exploration. This matches the prediction from the previous section: logarithmic cooling is guaranteed to eventually reach the global optimum, but the number of iterations required is orders of magnitude larger than what is practical. Linear cooling showed even worse tour lengths than logarithmic cooling at the 10% and 50% marks (1128.75, 1103.00), but improved rapidly (532.78) as the temperature dropped linearly to near-zero between the 50% and 100% marks — better than logarithmic cooling, though still far behind geometric cooling. This is likely because linear cooling does not reach a sufficiently low temperature until late in the run, leaving fewer iterations available for fine-grained local search compared to geometric cooling. Overall, the well-known rule of thumb that geometric cooling — despite lacking a theoretical convergence guarantee — outperforms logarithmic cooling in practice under a limited compute budget was quantitatively reproduced in this setting.

Comparison with Other Methods

Related Articles

• Genetic Algorithms: Fundamentals and Python Implementation - Another metaheuristic based on evolutionary computation.
• Bayesian Optimization: Fundamentals and Python Implementation - Efficient optimization for expensive evaluations.
• Markov Chain Monte Carlo (MCMC): Metropolis-Hastings and Gibbs Sampling - MCMC sampling using the same Metropolis criterion as SA.
• MPPI (Model Predictive Path Integral): A Unified View with the Cross-Entropy Method - Sampling-based optimal control.
• Cross-Entropy Method: A Practical Monte Carlo Optimization Technique - Useful comparison with CEM’s distribution-update approach.
• Particle Swarm Optimization (PSO) Python Implementation - Swarm intelligence-based optimization for comparison.
• Particle Swarm Optimization (PSO): Theory and Python Implementation - Updated population-based optimizer; useful contrast with SA’s single-solution stochastic acceptance.
• Gaussian Process Regression: Theory and Python Implementation - GP surrogate behind Bayesian optimization; complements SA when each evaluation is expensive.
• From SGD to Adam: Evolution of Gradient-Based Optimization - Gradient-based counterpart to SA’s gradient-free stochastic search; useful contrast across problem regimes.
• Monte Carlo Optimization Methods Comparison (CEM/SA/GA/MPPI/PSO) - Hub article organizing five Monte Carlo optimizers including SA under a common sample/evaluate/update skeleton, clarifying when single-solution SA versus population-based methods are appropriate.

References

• Kirkpatrick, S., Gelatt, C. D., & Vecchi, M. P. (1983). “Optimization by Simulated Annealing”. Science, 220(4598), 671-680.
• Metropolis, N., et al. (1953). “Equation of State Calculations by Fast Computing Machines”. The Journal of Chemical Physics, 21(6), 1087-1092.
• Bertsimas, D., & Tsitsiklis, J. (1993). “Simulated annealing”. Statistical Science, 8(1), 10-15.
• Geman, S., & Geman, D. (1984). “Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images”. IEEE Transactions on Pattern Analysis and Machine Intelligence, 6(6), 721-741.
• Hajek, B. (1988). “Cooling Schedules for Optimal Annealing”. Mathematics of Operations Research, 13(2), 311-329.