Temperature T controls whether simulated annealing accepts transitions to higher-energy states. Lower-energy proposals are accepted; higher-energy proposals may also be accepted, helping the search escape local optima. section 9.4 section 9.8
As the search proceeds:
- A high initial temperature permits many random transitions. section 9.4
- A cooling schedule lowers temperature, making energy differences more consequential and unfavorable transitions less likely. 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()
The implementation records temperatures and uses a schedule that multiplies the current temperature by 0.9. section 9.4
