SA Lab Walkthrough

Unit 4: Optimization and Constraint Satisfaction — SA Lab Walkthrough

This walkthrough is an optional companion to the Simulated Annealing Lab. Review it after completing and submitting your own solution.

Overview

This walkthrough presents complete implementations for all four lab parts, with line-by-line explanations and analysis of results. The key insight throughout: SA’s ability to escape local optima comes entirely from the Metropolis acceptance criterion, which allows controlled acceptance of worse solutions.

Part 1 Solutions: Basic SA Implementation

Function 1: Random Neighbor

def random_neighbor(current, step_size=0.1):
    """Generate a random neighboring state."""
    change = np.random.uniform(-step_size, step_size)
    return current + change

Why this works:

np.random.uniform(-step_size, step_size) generates a random value uniformly distributed between −step_size and +step_size. Adding this to current creates a neighbor within a ball of radius step_size around the current state.

The step_size parameter controls the exploration radius. Too small: the algorithm moves slowly and may stay trapped near its starting point. Too large: moves are so wild that the algorithm rarely improves.

A step size of 0.1 to 0.5 works well for most 1D optimization problems. For higher-dimensional problems, each dimension typically gets its own step size, or a multivariate normal distribution is used.

Function 2: Acceptance Probability

def acceptance_probability(current_energy, neighbor_energy, temp):
    """Calculate probability of accepting a move (Metropolis criterion)."""
    if neighbor_energy < current_energy:
        return 1.0                                    # always accept better
    else:
        delta_e = neighbor_energy - current_energy   # positive: neighbor is worse
        return math.exp(-delta_e / temp)

The Metropolis criterion explained:

For better neighbors (lower energy): always accept. This is standard hill climbing behavior — we never refuse a move that improves things.

For worse neighbors (higher energy): accept with probability e^(−ΔE/T).

The exponential function ensures:

When T is high: −ΔE/T is close to zero, so e^(−ΔE/T) ≈ 1 → accept almost anything
When T is low: −ΔE/T is very negative, so e^(−ΔE/T) ≈ 0 → rarely accept worse moves
When ΔE is small: close to 0, so more likely to accept
When ΔE is large: more negative, so less likely to accept

Note: we use energy (minimize) rather than fitness (maximize) here. For maximization problems, flip the sign: delta_e = current_value - neighbor_value.

Function 3: Main SA Algorithm

def simulated_annealing(objective_function, initial_state,
                        initial_temp, cooling_rate, max_iterations):
    """Run simulated annealing to minimize objective_function."""
    current = initial_state
    current_energy = objective_function(current)

    best = current
    best_energy = current_energy
    history = [current_energy]

    temp = initial_temp

    for iteration in range(max_iterations):
        # Generate neighbor and evaluate
        neighbor = random_neighbor(current, step_size=0.5)
        neighbor_energy = objective_function(neighbor)

        # Metropolis acceptance
        prob = acceptance_probability(current_energy, neighbor_energy, temp)
        if np.random.random() < prob:
            current = neighbor
            current_energy = neighbor_energy

        # Track best solution found
        if current_energy < best_energy:
            best = current
            best_energy = current_energy

        history.append(current_energy)

        # Cool temperature with minimum floor
        temp *= cooling_rate
        temp = max(temp, 0.01)

    return best, best_energy, history

Key implementation decisions:

Two-pointer tracking (current vs. best): current is where the algorithm is now — it may have moved to a worse state temporarily. best records the lowest energy state seen across all iterations. We always return best, not current, because SA’s final position may not be its best.

Temperature floor: temp = max(temp, 0.01) prevents temperature from reaching exactly zero, which would cause a division-by-zero in acceptance_probability. It also ensures the algorithm maintains a tiny amount of randomness throughout.

Step size: Using step_size=0.5 (larger than the default 0.1) allows SA to explore broadly during the high-temperature phase. You can make step size a parameter and schedule it to decrease along with temperature for more sophisticated behavior.

Part 2 Solutions: Cooling Schedules

def linear_cooling(t, t0, max_t):
    """T(t) = T0 * (1 - t/max_t)"""
    return t0 * (1 - t / max_t)

def exponential_cooling(t, t0, alpha):
    """T(t) = T0 * alpha^t"""
    return t0 * (alpha ** t)

def logarithmic_cooling(t, t0):
    """T(t) = T0 / log(1 + t)"""
    return t0 / math.log(1 + t)

Schedule comparison:

Linear cooling is the simplest but often too aggressive — temperature reaches zero at step max_t, leaving the algorithm no time to refine later.

Exponential cooling (α = 0.95) is the standard choice. Temperature decays by 5% each step; after 100 steps, temperature is T₀ × 0.95^100 ≈ 0.006 × T₀. Increasing α slows cooling (more exploration); decreasing it speeds cooling (faster convergence but higher risk of getting stuck).

Logarithmic cooling is theoretically correct (it guarantees convergence to the global optimum as iterations → ∞) but impractically slow. After 1,000 steps with T₀ = 100: T ≈ 100 / log(1001) ≈ 14.5. After 1,000,000 steps: T ≈ 100 / log(1,000,001) ≈ 7.2. The temperature barely decreases.

Part 3 Results: SA vs. Hill Climbing

Expected Output

Results on Rastrigin function (20 runs each):
  Hill Climbing:  mean=4.21, std=2.13, best=0.98
  SA (T=100, α=0.95): mean=0.83, std=0.47, best=0.01

Results on Schwefel function (20 runs each):
  Hill Climbing:  mean=87.3, std=44.1, best=12.5
  SA (T=100, α=0.95): mean=11.2, std=7.9, best=0.3

Analysis

SA achieves roughly 80% better mean solution quality on the Rastrigin function and 87% better on the Schwefel function compared to hill climbing. The difference arises because both functions have many local optima: hill climbing gets trapped in the first one it finds, while SA’s temperature mechanism allows it to escape and search more broadly.

On the unimodal quadratic function, both algorithms perform nearly identically — there are no local optima to escape, so the temperature mechanism provides no advantage.

The standard deviation for SA is also much lower (more consistent results), which is valuable in practice: you can rely on SA to produce good results reliably, not just occasionally.

Part 4 Analysis: Parameter Effects

Temperature Analysis Results

Higher initial temperatures lead to better solutions up to a point. Below T₀ ≈ 10, SA behaves like hill climbing and gets trapped. Above T₀ ≈ 50, SA explores well and finds near-optimal solutions. Excessively high T₀ > 500 slightly hurts because too much time is spent on random exploration with little refinement.

The sweet spot for the Rastrigin function in 1D: T₀ ∈ [50, 200].

Cooling Rate Analysis Results

Faster cooling (lower α) produces worse solutions because the algorithm does not spend enough time at each temperature level. At α = 0.80, cooling is so rapid that SA behaves like a few hundred random-restart hill climbing runs. At α = 0.999, cooling is so slow that SA barely reaches a cool phase within 5,000 iterations.

Recommended: α ∈ [0.93, 0.97] for most problems with 5,000 iterations. Increase α (slow cooling) for harder, higher-dimensional problems; decrease it when function evaluations are expensive.

Key Takeaways

The simulated annealing algorithm has three essential components:

Random neighbor selection — ensures full exploration of the neighborhood
Metropolis acceptance criterion — allows controlled acceptance of worse moves, scaled by temperature
Cooling schedule — transitions the algorithm from broad exploration to focused refinement

Remove any one component and SA degrades: without random neighbors you get hill climbing, without Metropolis acceptance you get random search, without cooling you get a random walk. The power of SA comes from all three working together.

Next: Sudoku CSP Solver Lab →

This work is licensed under CC BY-SA 4.0.