Simulated annealing (SA) searches for a low-energy state when an energy landscape may contain local optima. Unlike a search that accepts only improvements, it sometimes moves to a higher-energy state. Those temporary setbacks can open a route to a lower minimum later. section 9.4
Let E(x) be the energy at state x and T the temperature. Propose a new state x′ and compare E(x′) with E(x). Accept a lower-energy proposal with probability 1; accept a higher-energy proposal with a probability determined by the energy difference and T. A high initial temperature permits more unfavorable moves. As a cooling schedule lowers T, such moves become less likely. section 9.4
A run follows this loop:
- Initialize the state and temperature, then evaluate the initial energy. section 9.4
- Sample a proposal, evaluate its energy, and check for convergence. section 9.4
- Apply the temperature-dependent acceptance rule, update the temperature, and repeat. Return the coordinates with the lowest energy found. section 9.4
Two-dimensional energy search
A Python implementation uses a target function E(x, y) and proposes coordinates from a multivariate normal distribution centered on the current point, with covariance 1.1I. Its cooling schedule multiplies T by 0.9 each iteration. The example also sets a 1,000-iteration limit and convergence parameters, making the proposal, cooling, and stopping choices explicit. section 9.4
Listing 9.5 Simulated annealing
import numpy as np
import matplotlib.pyplot as plt
np.random.seed(42)
class simulated_annealing():
def __init__(self):
self.max_iter = 1000
self.conv_thresh = 1e-4
self.conv_window = 10
self.samples = np.zeros((self.max_iter, 2))
self.energies = np.zeros(self.max_iter)
self.temperatures = np.zeros(self.max_iter)
def target(self, x, y): #1
z = 3*(1-x)**2 * np.exp(-x**2 - (y+1)**2) \
- 10*(x/5 -x**3 - y**5) * np.exp(-x**2 - y**2) \
- (1/3)*np.exp(-(x+1)**2 - y**2)
return z
def proposal(self, x, y):
mean = np.array([x, y])
cov = 1.1 * np.eye(2)
x_new, y_new = np.random.multivariate_normal(mean, cov)
return x_new, y_new
def temperature_schedule(self, T, iter):
return 0.9 * T
def run(self, x_init, y_init):
converged = False
T = 1
self.temperatures[0] = T
num_accepted = 0
x_old, y_old = x_init, y_init
energy_old = self.target(x_init, y_init)
iter = 1
while not converged:
print("iter: {:4d}, temp: {:.4f}, energy = {:.6f}".format(iter,
T, energy_old))
x_new, y_new = self.proposal(x_old, y_old)
energy_new = self.target(x_new, y_new)
if iter > 2*self.conv_window: #2
vals = self.energies[iter-self.conv_window : iter-1]
if (np.std(vals) < self.conv_thresh):
converged = True
#end if
#end if
alpha = np.exp((energy_old - energy_new)/T)
r = np.minimum(1, alpha)
u = np.random.uniform(0, 1)
if u < r:
x_old, y_old = x_new, y_new
num_accepted += 1
energy_old = energy_new
#end if
self.samples[iter, :] = np.array([x_old, y_old])
self.energies[iter] = energy_old
T = self.temperature_schedule(T, iter)
self.temperatures[iter] = T
iter = iter + 1
if (iter > self.max_iter): converged = True
#end while
niter = iter - 1
acceptance_rate = num_accepted / niter
print("acceptance rate: ", acceptance_rate)
x_opt, y_opt = x_old, y_old
return x_opt, y_opt, self.samples[:niter,:], self.energies[:niter],
self.temperatures[:niter]
if __name__ == "__main__":
SA = simulated_annealing()
nx, ny = (1000, 1000)
x = np.linspace(-2, 2, nx)
y = np.linspace(-2, 2, ny)
xv, yv = np.meshgrid(x, y)
z = SA.target(xv, yv)
plt.figure()
plt.contourf(x, y, z)
plt.title("energy landscape")
plt.show()
#find global minimum by exhaustive search
min_search = np.min(z)
argmin_search = np.argwhere(z == min_search)
xmin, ymin = argmin_search[0][0], argmin_search[0][1]
print("global minimum (exhaustive search): ", min_search)
print("located at (x, y): ", x[xmin], y[ymin])
#find global minimum by simulated annealing
x_init, y_init = 0, 0
x_opt, y_opt, samples, energies, temperatures = SA.run(x_init, y_init)
print("global minimum (simulated annealing): ", energies[-1])
print("located at (x, y): ", x_opt, y_opt)
plt.figure()
plt.plot(energies)
plt.title("SA sampled energies")
plt.show()
plt.figure()
plt.plot(temperatures)
plt.title("Temperature Schedule")
plt.show()
Inspect the target function, proposal distribution, cooling schedule, and convergence parameters in the implementation. section 9.4
For difficult searches, different starting points and several restarts can improve the optima found. The traveling salesman problem is one example of an NP-hard problem that can be approximated with SA. section 9.6
