Where SA Fits in the SO Landscape

General SO problem:

\[\min_{x \in X} \; f(x)\]

where \(f\) is a black-box objective observable only through a stochastic simulation oracle.

SA targets global search using a single solution at a time.

Key properties that make SA attractive for simulation optimization:

• Single-solution method — only one oracle evaluation per iteration; no population to maintain.
• Low memory footprint — critical when the simulation environment is large and cannot be duplicated.
• Noise-tolerant — the probabilistic acceptance step naturally absorbs stochastic noise in the objective.
• Simple to implement — the core algorithm requires only a neighbor generator and an acceptance rule.

Why Local Search Is Insufficient

Local search pseudocode (minimization):

1. Draw an initial solution i
2. Generate a solution j from the neighborhood N(i)
3. If f(j) < f(i) then j becomes the current solution
4. If f(j) ≥ f(i) for all j ∈ N(i) then STOP
5. Go to step 2

The fundamental problem: once the current solution falls into a locally convex subdomain, the algorithm is trapped.

• A local optimum \(i^*\) satisfies \(f(i^*) \leq f(j)\) for all \(j \in N(i^*)\).
• Unless the neighborhood is exact (implying full enumeration), local search cannot guarantee global optimality.
• In simulation optimization, local optima abound — the noisy objective surface creates many spurious basins.

SA’s solution: allow transitions that worsen the objective with a controlled probability, enabling escape from local optima.

SA Pseudocode

SA Algorithm (minimization):

1. Initialize:  i := i_start,  k := 0,  c_k := c_0,  L_k := L_0
2. Repeat:
   For l = 1 to L_k do:                    ← inner loop at temperature c_k
     a. Generate j from N(i)               ← neighbor generation
     b. Δf := f(j) - f(i)                  ← cost difference
     c. If Δf < 0:  i := j                 ← always accept improvement
        Else: accept j with probability exp(-Δf / c_k)
   k := k + 1                              ← advance to next temperature level
   Compute (L_k, c_k)                      ← update temperature and chain length
3. Until c_k ≈ 0

Two nested loops:

• Outer loop: decrements the temperature \(c_k\) according to the cooling schedule.
• Inner loop (Markov chain): generates \(L_k\) transitions at fixed temperature \(c_k\).

Note on KSL implementation: the inner loop is absorbed into the outer loop — each call to mainIteration() generates and evaluates one neighbor, then advances the temperature. The Markov chain length concept is managed implicitly through the iteration budget.

The Acceptance Function in KSL

fun acceptanceProbability(costDifference: Double, temperature: Double): Double {
    require(temperature > 0.0) { "The temperature must be positive" }
    // Improvements (Δf ≤ 0) are always accepted
    return if (costDifference <= 0.0) {
        1.0
    } else {
        // Worsening moves: Metropolis criterion  exp(-Δf / c)
        exp(-costDifference / temperature)
    }
}

Algorithm mapping:

Diagnostic properties exposed by SimulatedAnnealing:

• lastAcceptanceProbability — the probability computed at the last iteration.
• costDifference — the \(\Delta f\) computed at the last iteration.
• currentTemperature — the current value of \(c_k\).

All three are private set — readable for diagnostics, not writable externally.

Statistical Equilibrium

Theorem (Boltzmann stationary distribution): after a sufficient number of transitions at fixed \(c\), the SA algorithm visits solution \(i\) with probability:

\[q_i(c) = \frac{1}{N_0(c)} e^{-f(i)/c}, \qquad N_0(c) = \sum_{j \in S} e^{-f(j)/c}\]

Key limiting behaviors:

Practical meaning:

• At high \(c\): SA is a random walk — all solutions visited uniformly.
• At low \(c\): SA concentrates probability mass on global optima.
• The entropy \(H_c\) measures the degree of optimality achieved — decreasing \(H_c\) tracks convergence.

Asymptotic Convergence

Theorem (Theorem 2, Delahaye et al.): Let the transition probabilities satisfy the irreducibility condition — every state \(j\) is reachable from every state \(i\) via a finite sequence of generating steps. Then the Markov chain has a stationary distribution \(q(c)\) and:

\[\lim_{c \to 0} \lim_{k \to \infty} \Pr\{X^c_k \in S_\text{opt}\} = 1\]

In words: SA converges to a global optimum with probability 1, provided:

1. The generating mechanism is irreducible (can reach any solution from any other).
2. The cooling is infinitely slow (\(c \to 0\) continuously, with infinite transitions at each level).

Practical implications:

• Convergence requires an infinitely long cooling schedule — not achievable in practice.
• In finite time, SA approximates the global optimum, not guarantees it.
• The algorithm moves the search into the right basin of attraction; a local refinement step can then find the local minimum within that basin.
• The transition threshold \(c_t \approx \frac{2\sigma^2_\infty}{\langle f \rangle_\infty - f_\text{min}}\) divides the search into an exploration phase (\(c > c_t\)) and an exploitation phase (\(c < c_t\)).

Geometric (Exponential) Cooling — ExponentialCoolingSchedule

Formula:

\[c_{k+1} = \alpha \, c_k \implies c_k = c_0 \cdot \alpha^k, \qquad 0 < \alpha < 1\]

Typical range: \(0.80 < \alpha < 0.99\).

class ExponentialCoolingSchedule(
    initialTemperature: Double,
    val coolingRate: Double = defaultCoolingRate   // α, default 0.95
) : CoolingSchedule(initialTemperature) {

    override fun nextTemperature(iteration: Int): Double {
        return initialTemperature * coolingRate.pow(iteration)
    }
}

Auto-calculated cooling rate — given target \(c_0\), \(c_f\), and maxIterations \(K\):

\[\alpha = \left(\frac{c_f}{c_0}\right)^{1/K}\]

// Factory: calculates α to hit c_f exactly at iteration K
ExponentialCoolingSchedule.createTargetedSchedule(
    initialTemperature = 500.0,
    stoppingTemperature = 0.001,
    maxIterations      = 2000
)

Advantages: simple, robust, widely used, no problem-specific tuning of the decay form. Disadvantage: does not adapt to local objective landscape — may spend too many iterations in unproductive regions.

Linear Cooling — LinearCoolingSchedule

Formula:

\[c_k = c_0 - k \cdot \frac{c_0 - c_f}{K}, \qquad k = 0, 1, \ldots, K\]

class LinearCoolingSchedule(
    initialTemperature: Double = defaultInitialTemperature,
    val stoppingTemperature: Double = defaultStoppingTemperature,
    val maxIterations: Int = defaultMaxNumberIterations
) : CoolingSchedule(initialTemperature) {

    val temperatureDecreasePerIteration
        get() = (initialTemperature - stoppingTemperature) / maxIterations

    override fun nextTemperature(iteration: Int): Double {
        if (iteration >= maxIterations) return stoppingTemperature
        return initialTemperature - (temperatureDecreasePerIteration * iteration)
    }
}

Characteristics: temperature decreases at a constant rate per iteration. Spends equal time at every temperature level — suitable when the objective landscape is relatively uniform in scale.

Why \(c_0\) Matters

Too low \(c_0\): the algorithm behaves like greedy local search from the start — essentially no exploration, high risk of trapping in the first local minimum encountered.

Too high \(c_0\): wastes many iterations on a random walk, accepting essentially everything; exploration is redundant until \(c\) drops to a useful range.

Goal: set \(c_0\) so that the initial acceptance rate \(\chi(c_0)\) equals a target value (typically 0.8), ensuring nearly all transitions are accepted at the start.

Polynomial-time initialization formula (Eq. 2, Delahaye et al.):

Let \(m_1\) = number of improving transitions and \(m_2\) = number of worsening transitions observed in a random walk sample. Let \(\bar{\Delta f}^{(+)}\) be the average cost increase over worsening moves. Then:

\[c_0 = \frac{\bar{\Delta f}^{(+)}}{\ln\!\left(\dfrac{m_2}{m_2 \, \chi(c_0) - m_1\,(1-\chi(c_0))}\right)}\]

Simplified form used in KSL (when \(m_1 \ll m_2\), which is typical for a random walk): \(c_0 = \frac{-\bar{\Delta f}^{(+)}}{\ln(\chi_0)}\) where \(\chi_0\) is the target initial acceptance probability.

calibrateTemperature() — Annotated Code

private fun calibrateTemperature(config: TemperatureConfiguration.AutoCalibrate): Double {
    var totalWorseningCost = 0.0
    var worseningMovesCount = 0
    var previousWalkSolution = currentSolution   // baseline from initializeIterations()
    var bestWalkSolution     = currentSolution
    for (i in 0 until config.numRandomWalkSteps) {    // random walk of numRandomWalkSteps steps
        val nextPoint = generateNeighbor(previousWalkSolution.inputMap, rnStream)
        // Evaluate the neighbor — uses numRepsPerStep if > 0, else model default
        val nextSolution = if (config.numRepsPerStep > 0)
            requestEvaluation(nextPoint, config.numRepsPerStep)
        else
            requestEvaluation(nextPoint)
        // Track best solution seen during walk (warm-start option)
        if (nextSolution.penalizedObjFncValue < bestWalkSolution.penalizedObjFncValue)
            bestWalkSolution = nextSolution
        val costDiff = nextSolution.penalizedObjFncValue - previousWalkSolution.penalizedObjFncValue
        if (costDiff > 0.0) {             // worsening move: accumulate Δf⁺
            totalWorseningCost  += costDiff
            worseningMovesCount += 1
        }
        previousWalkSolution = nextSolution
    }
    if (worseningMovesCount == 0) return defaultInitialTemperature   // fallback
    // Apply starting-point strategy
    if (config.useBestStepForStartingPoint) currentSolution = bestWalkSolution
    val avgWorseningCost = totalWorseningCost / worseningMovesCount
    return -avgWorseningCost / ln(config.targetProbability)   // c_0 = -Δf̄⁺ / ln(χ_0)
}

Algorithm mapping:

Stopping Temperature, Iterations, and Replication Noise

Stopping temperature \(c_f\):

• Sets the exploitation floor — below \(c_f\), virtually no worsening moves are accepted.
• Default 0.001 is suitable when \(f\) is scaled to order of magnitude \(\sim 1\)–\(100\).
• If \(f\) values are very large (e.g., \(\sim 10^6\)), \(c_f\) must be scaled proportionally — otherwise the algorithm effectively never stops cooling.

Max iterations maxIterations:

• Acts as a hard safety budget, independent of temperature.
• Should be set large enough that the temperature criterion (\(c < c_f\)) triggers first under normal operation.

Replications per evaluation replicationsPerEvaluation:

• Higher replications \(\Rightarrow\) lower noise in \(\hat{f}(x)\) \(\Rightarrow\) more reliable acceptance decisions.
• More costly: each iteration runs multiple simulation replications.
• Low replications are acceptable at high \(c\) (noise is dominated by randomness anyway); consider increasing them as \(c\) decreases if the budget allows.

No-Improve Threshold and Stagnation Stopping

isStoppingCriteriaSatisfied():

override fun isStoppingCriteriaSatisfied(): Boolean {
    return solutionQualityEvaluator?.isStoppingCriteriaReached(this)
        ?: checkTemperature() || solutionChecker.checkSolutions()
}

// Criterion 1: temperature has dropped below stopping threshold
private fun checkTemperature() = currentTemperature < stoppingTemperature

// Criterion 2: no improvement in best solution for d iterations
// solutionChecker.checkSolutions() returns true when stagnation detected

Two independent stopping criteria — either triggers termination:

Key Takeaways

1. SA escapes local optima through probabilistic acceptance. The Metropolis criterion allows worsening moves with probability \(e^{-\Delta f / c}\), which decreases as \(c\) cools — balancing exploration and exploitation.

2. Asymptotic convergence is guaranteed but impractical. The theoretical convergence proof requires an infinitely slow schedule. In practice, SA finds good approximate solutions within a finite budget.

3. The cooling schedule is the primary design choice. The geometric schedule with \(\alpha \in (0.80, 0.99)\) is the robust default; createTargetedSchedule() automates \(\alpha\) selection given \(c_0\), \(c_f\), and \(K\).

4. Initial temperature strongly influences solution quality. Auto-calibration via a random walk is the recommended approach when no domain knowledge is available — it sets \(c_0\) to achieve a target acceptance rate of \(\approx 0.8\).

5. SA is ideally suited for memory-constrained simulation optimization. A single solution is maintained at all times — no population to evaluate, no simulation environment duplication required.

6. KSL separates temperature strategy from the SA loop. Swapping TemperatureConfiguration or CoolingScheduleIfc implementations changes behavior without modifying SimulatedAnnealing itself.