Simulated annealing in Python

This small notebook implements, in Python 3, the simulated annealing algorithm for numerical optimization.

References

See also

About


This notebook should be compatible with both Python versions, 2 and 3.

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.


That's it for today, folks!

More notebooks can be found on my GitHub page.