Simulated annealing in Python¶
This small notebook implements, in Python 3, the simulated annealing algorithm for numerical optimization.
References¶
- The Wikipedia page: simulated annealing.
- It was implemented in
scipy.optimizebefore version 0.14:scipy.optimize.anneal. - This blog post.
- These Stack Overflow questions: 15853513 and 19757551.
See also¶
- For a real-world use of simulated annealing, this Python module seems useful: perrygeo/simanneal on GitHub.
About¶
- Date: 20/07/2017.
- Author: Lilian Besson, (C) 2017.
- Licence: MIT Licence.
In [1]:
from __future__ import print_function, division # Python 2 compatibility if needed
In [2]:
import numpy as np
import numpy.random as rn
import matplotlib.pyplot as plt # to plot
import matplotlib as mpl
from scipy import optimize # to compare
import seaborn as sns
sns.set(context="talk", style="darkgrid", palette="hls", font="sans-serif", font_scale=1.05)
FIGSIZE = (19, 8) #: Figure size, in inches!
mpl.rcParams['figure.figsize'] = FIGSIZE
Algorithm¶
The following pseudocode presents the simulated annealing heuristic.
- It starts from a state $s_0$ and continues to either a maximum of $k_{\max}$ steps or until a state with an energy of $e_{\min}$ or less is found.
- In the process, the call $\mathrm{neighbour}(s)$ should generate a randomly chosen neighbour of a given state $s$.
- The annealing schedule is defined by the call $\mathrm{temperature}(r)$, which should yield the temperature to use, given the fraction $r$ of the time budget that has been expended so far.
Simulated Annealing:
- Let $s$ = $s_0$
- For $k = 0$ through $k_{\max}$ (exclusive):
- $T := \mathrm{temperature}(k ∕ k_{\max})$
- Pick a random neighbour, $s_{\mathrm{new}} := \mathrm{neighbour}(s)$
- If $P(E(s), E(s_{\mathrm{new}}), T) \geq \mathrm{random}(0, 1)$:
- $s := s_{\mathrm{new}}$
- Output: the final state $s$
Basic but generic Python code¶
Let us start with a very generic implementation:
In [3]:
def annealing(random_start,
cost_function,
random_neighbour,
acceptance,
temperature,
maxsteps=1000,
debug=True):
""" Optimize the black-box function 'cost_function' with the simulated annealing algorithm."""
state = random_start()
cost = cost_function(state)
states, costs = [state], [cost]
for step in range(maxsteps):
fraction = step / float(maxsteps)
T = temperature(fraction)
new_state = random_neighbour(state, fraction)
new_cost = cost_function(new_state)
if debug: print("Step #{:>2}/{:>2} : T = {:>4.3g}, state = {:>4.3g}, cost = {:>4.3g}, new_state = {:>4.3g}, new_cost = {:>4.3g} ...".format(step, maxsteps, T, state, cost, new_state, new_cost))
if acceptance_probability(cost, new_cost, T) > rn.random():
state, cost = new_state, new_cost
states.append(state)
costs.append(cost)
# print(" ==> Accept it!")
# else:
# print(" ==> Reject it...")
return state, cost_function(state), states, costs
Basic example¶
We will use this to find the global minimum of the function $x \mapsto x^2$ on $[-10, 10]$.
In [4]:
interval = (-10, 10)
def f(x):
""" Function to minimize."""
return x ** 2
def clip(x):
""" Force x to be in the interval."""
a, b = interval
return max(min(x, b), a)
In [5]:
def random_start():
""" Random point in the interval."""
a, b = interval
return a + (b - a) * rn.random_sample()
In [6]:
def cost_function(x):
""" Cost of x = f(x)."""
return f(x)
In [7]:
def random_neighbour(x, fraction=1):
"""Move a little bit x, from the left or the right."""
amplitude = (max(interval) - min(interval)) * fraction / 10
delta = (-amplitude/2.) + amplitude * rn.random_sample()
return clip(x + delta)
In [8]:
def acceptance_probability(cost, new_cost, temperature):
if new_cost < cost:
# print(" - Acceptance probabilty = 1 as new_cost = {} < cost = {}...".format(new_cost, cost))
return 1
else:
p = np.exp(- (new_cost - cost) / temperature)
# print(" - Acceptance probabilty = {:.3g}...".format(p))
return p
In [9]:
def temperature(fraction):
""" Example of temperature dicreasing as the process goes on."""
return max(0.01, min(1, 1 - fraction))
Let's try!
In [10]:
annealing(random_start, cost_function, random_neighbour, acceptance_probability, temperature, maxsteps=30, debug=True);
Step # 0/30 : T = 1, state = -7.45, cost = 55.5, new_state = -7.45, new_cost = 55.5 ... Step # 1/30 : T = 0.967, state = -7.45, cost = 55.5, new_state = -7.44, new_cost = 55.4 ... Step # 2/30 : T = 0.933, state = -7.44, cost = 55.4, new_state = -7.5, new_cost = 56.2 ... Step # 3/30 : T = 0.9, state = -7.5, cost = 56.2, new_state = -7.59, new_cost = 57.6 ... Step # 4/30 : T = 0.867, state = -7.59, cost = 57.6, new_state = -7.64, new_cost = 58.3 ... Step # 5/30 : T = 0.833, state = -7.59, cost = 57.6, new_state = -7.51, new_cost = 56.4 ... Step # 6/30 : T = 0.8, state = -7.51, cost = 56.4, new_state = -7.53, new_cost = 56.6 ... Step # 7/30 : T = 0.767, state = -7.53, cost = 56.6, new_state = -7.58, new_cost = 57.5 ... Step # 8/30 : T = 0.733, state = -7.53, cost = 56.6, new_state = -7.6, new_cost = 57.8 ... Step # 9/30 : T = 0.7, state = -7.53, cost = 56.6, new_state = -7.51, new_cost = 56.4 ... Step #10/30 : T = 0.667, state = -7.51, cost = 56.4, new_state = -7.24, new_cost = 52.4 ... Step #11/30 : T = 0.633, state = -7.24, cost = 52.4, new_state = -6.98, new_cost = 48.7 ... Step #12/30 : T = 0.6, state = -6.98, cost = 48.7, new_state = -6.6, new_cost = 43.5 ... Step #13/30 : T = 0.567, state = -6.6, cost = 43.5, new_state = -6.69, new_cost = 44.8 ... Step #14/30 : T = 0.533, state = -6.6, cost = 43.5, new_state = -6.84, new_cost = 46.8 ... Step #15/30 : T = 0.5, state = -6.6, cost = 43.5, new_state = -6.45, new_cost = 41.6 ... Step #16/30 : T = 0.467, state = -6.45, cost = 41.6, new_state = -6.24, new_cost = 38.9 ... Step #17/30 : T = 0.433, state = -6.24, cost = 38.9, new_state = -6.52, new_cost = 42.5 ... Step #18/30 : T = 0.4, state = -6.24, cost = 38.9, new_state = -5.92, new_cost = 35.1 ... Step #19/30 : T = 0.367, state = -5.92, cost = 35.1, new_state = -6.35, new_cost = 40.4 ... Step #20/30 : T = 0.333, state = -5.92, cost = 35.1, new_state = -5.98, new_cost = 35.8 ... Step #21/30 : T = 0.3, state = -5.92, cost = 35.1, new_state = -5.35, new_cost = 28.6 ... Step #22/30 : T = 0.267, state = -5.35, cost = 28.6, new_state = -4.67, new_cost = 21.8 ... Step #23/30 : T = 0.233, state = -4.67, cost = 21.8, new_state = -4.44, new_cost = 19.7 ... Step #24/30 : T = 0.2, state = -4.44, cost = 19.7, new_state = -4.59, new_cost = 21.1 ... Step #25/30 : T = 0.167, state = -4.44, cost = 19.7, new_state = -4.04, new_cost = 16.3 ... Step #26/30 : T = 0.133, state = -4.04, cost = 16.3, new_state = -4.77, new_cost = 22.8 ... Step #27/30 : T = 0.1, state = -4.04, cost = 16.3, new_state = -4.7, new_cost = 22.1 ... Step #28/30 : T = 0.0667, state = -4.04, cost = 16.3, new_state = -3.44, new_cost = 11.8 ... Step #29/30 : T = 0.0333, state = -3.44, cost = 11.8, new_state = -2.6, new_cost = 6.78 ...
Now with more steps:
In [11]:
state, c, states, costs = annealing(random_start, cost_function, random_neighbour, acceptance_probability, temperature, maxsteps=1000, debug=False)
state
c
Out[11]:
0.08501260465044064
Out[11]:
0.007227142949452121
Visualizing the steps¶
In [12]:
def see_annealing(states, costs):
plt.figure()
plt.suptitle("Evolution of states and costs of the simulated annealing")
plt.subplot(121)
plt.plot(states, 'r')
plt.title("States")
plt.subplot(122)
plt.plot(costs, 'b')
plt.title("Costs")
plt.show()
In [13]:
see_annealing(states, costs)
More visualizations¶
In [14]:
def visualize_annealing(cost_function):
state, c, states, costs = annealing(random_start, cost_function, random_neighbour, acceptance_probability, temperature, maxsteps=1000, debug=False)
see_annealing(states, costs)
return state, c
In [15]:
visualize_annealing(lambda x: x**3)
Out[15]:
(-10, -1000)
In [16]:
visualize_annealing(lambda x: x**2)
Out[16]:
(-0.007545319790391913, 5.6931850739279866e-05)
In [17]:
visualize_annealing(np.abs)
Out[17]:
(-0.007160473894364527, 0.0071604738943645274)
In [18]:
visualize_annealing(np.cos)
Out[18]:
(9.491236787419837, -0.9977924248899227)
In [19]:
visualize_annealing(lambda x: np.sin(x) + np.cos(x))
Out[19]:
(3.935615512590335, -1.4141609642965898)
In all these examples, the simulated annealing converges to a global minimum. It can be non-unique, but it is found.