How to escape from a local optimum · Simulated annealing

Part 4 · Heuristic algorithms · Chapter 32 · Chapter 32: Metaheuristics · Lecture notes, volume 4 (PDF) · Slides (1/2) · Slides (2/2) · Exercises: Metaheuristics · App: TSP: simulated annealing and tabu search

32.3 How to escape from a local optimum

• A local search moves from a solution to a better neighbour, and stops at the first local optimum: a solution with no better neighbour. In the landscape of the objective function, it descends to the bottom of the valley where it starts, and stays there. This is pure intensification.

• Repeating the local search from many starting solutions (multi-start) explores many valleys, but each start ignores everything learnt before, and good valleys are usually found close to other good valleys.

• The table below lists the ways to escape from a local optimum; it is the map of the chapter.

how to escape	idea	metaheuristic
accept worse moves at random	a worse neighbour is accepted with a probability that decreases with the worsening and with time	simulated annealing (Section 32.4)
accept worse moves deterministically, with memory	move to the best neighbour even if it is worse, and forbid going back to the solutions just visited	tabu search (Section 32.5)
perturb the local optimum and restart	apply a random modification that the local search cannot undo, and descend again	iterated local search (Section 32.6)
change the neighbourhood	when the search is stuck, use a larger neighbourhood, of which the local optimum is no longer a local optimum	variable neighbourhood search (Section 32.7)
restart from a randomized construction	build a different good solution at every restart, and improve it	GRASP (Section 32.8)
destroy and repair	remove a large part of the solution and rebuild it with a heuristic or an exact algorithm	large neighbourhood search (Section 32.9)
work with a population	combine the parts of several good solutions into new ones	genetic and memetic algorithms (Section 32.10)

• The first two ways keep one current solution and change the acceptance rule of the local search: they walk out of the valley through worse solutions. The next four jump: they generate a new starting point, far enough to leave the valley but close enough to keep what is good. The last one keeps many solutions at the same time.

• Every metaheuristic must balance intensification (descend to the bottom of the current valley) and diversification (leave it for another one). For each method we point out which component does what.

• A metaheuristic is a scheme, not an algorithm: to obtain an algorithm one must choose the representation, the neighbourhoods, the local search and the parameters, and these choices matter at least as much as the scheme. We present every metaheuristic on the problems where its idea is most natural, not on all four.

32.4 Simulated annealing

Simulated annealing (denoted by \(SA\)) is a local search that chooses a random neighbour of the current solution and moves to it if it is better; if it is worse, it moves to it anyway with a probability that decreases with the worsening and with a parameter \(T\), the temperature. The temperature starts high, so that at the beginning almost every move is accepted, and is lowered gradually, so that at the end only improving moves are accepted. The best solution visited is returned.

• The pseudo-code of \(SA\), for a minimization problem with objective \(f\) and neighbourhood \(N\), reads as follows:

Algorithm: Simulated annealing

Input: starting solution \(S\), initial temperature \(T_0\), cooling factor \(\alpha\in(0,1)\), number \(L\) of moves per temperature, final temperature \(T_{\min}\)

Output: best solution found

1Function SimulatedAnnealing(\(S, T_0, \alpha, L, T_{\min}\))

2\(cur := S\); \(fcur := f(cur)\) \(best := cur\); \(fbest := fcur\) \(T := T_0\) ▹ current solution best solution found temperature

3while \(T > T_{\min}\) do

4for \(l := 1\) to \(L\) do

5\(S' :=\) random neighbour in \(N(cur)\)

6\(\delta := f(S') - fcur\)

7\(v := \textsc{Rand}()\)

8if \(\delta \le 0\) or \(v < e^{-\delta/T}\) then

9\(cur := S'\); \(fcur := fcur + \delta\) ▹ move accepted

10if \(fcur < fbest\) then

11\(best := cur\); \(fbest := fcur\)

12end

13end

14end

15\(T := \alpha\cdot T\) ▹ cooling

16end

17return \(best\)

Run in Colab:PythonC++

• The name comes from metallurgy: annealing is the process of heating a metal and cooling it slowly, so that its atoms settle into a state of minimum energy; cooled too fast, the metal freezes into a disordered state of higher energy. Kirkpatrick, Gelatt and Vecchi (1983) and, independently, Černy (1985) proposed to simulate this process on an optimization problem, with the objective \(f\) in the role of the energy. The acceptance rule is the one used by Metropolis et al. (1953) to simulate the atoms.

Metropolis acceptance rule. At temperature \(T\), a neighbour \(S'\) of the current solution, with \(\delta=f(S')-f(cur)\), is accepted with probability

\[ P(\text{accept } S') \;=\; \begin{cases} 1 & \text{if } \delta\le 0,\\[1ex] e^{-\delta/T} & \text{if } \delta>0. \end{cases} \]

• The probability is implemented by drawing \(v=\textsc{Rand}()\), uniform in \((0,1)\), and accepting if \(v<e^{-\delta/T}\), which happens with probability \(e^{-\delta/T}\).

• A worsening move is accepted more easily if the worsening \(\delta\) is small, and if the temperature \(T\) is high. As \(T\to\infty\) every move is accepted, and \(SA\) becomes a random walk (pure diversification); as \(T\to 0\) only the moves with \(\delta\le 0\) are accepted, and \(SA\) becomes a local search with random choice of the neighbour (pure intensification). Lowering the temperature moves gradually from the first behaviour to the second.

• Moves with \(\delta=0\) are always accepted: \(SA\) walks freely on the plateaus of flat objectives.

• \(SA\) never stops at a local optimum: a local optimum is left as soon as a worsening move is drawn and accepted. For this reason the best solution visited must be stored separately: the current solution at the end is not necessarily the best one.

• Only one random neighbour is evaluated per iteration, not the whole neighbourhood: an iteration of \(SA\) is very cheap, and \(SA\) performs many more iterations than a local search. For the TSP, a random 2-opt move is drawn by numbering the \(M=\frac{n(n-3)}{2}\) moves \((i,j)\) in lexicographic order and choosing the one in position \(1+\lfloor \textsc{Rand}()\cdot M\rfloor\); its \(\delta\) is computed in \(O(1)\).

• The figure shows the acceptance probability as a function of the worsening, for four temperatures. At \(T=10\) a worsening of \(5\) is accepted with probability \(e^{-0.5}\approx 0.61\); at \(T=1\) with probability \(e^{-5}\approx 0.007\). The temperature has the same unit as the objective: \(T\) is “large” or “small” only compared with the typical values of \(\delta\).

32.4.1 The cooling schedule and the parameters

• The cooling schedule is the sequence of temperatures. The pseudo-code uses the geometric schedule \(T\leftarrow\alpha T\), the most common in practice, with \(\alpha\) between \(0.8\) and \(0.99\), and \(L\) moves at each temperature. Other schedules are the linear one, \(T_k=T_0-k\beta\), and the logarithmic one, \(T_k=\frac{\Gamma}{\ln(k+1)}\), which matters in the theory of \(SA\) but is far too slow in practice.

• Initial temperature \(T_0\). It must be high enough to accept most worsening moves at the beginning. A common rule fixes a target initial acceptance ratio \(\chi_0\) (for instance \(\chi_0=0.8\)), estimates the average worsening \(\bar\delta\) of a sample of random moves, and sets \(T_0=-\bar\delta/\ln\chi_0\), so that \(e^{-\bar\delta/T_0}=\chi_0\).

• Moves per temperature \(L\). The search must be given time to adapt to each temperature; \(L\) is often taken proportional to the size of the neighbourhood.

• Stopping rule. The pseudo-code stops when \(T\le T_{\min}\). With the geometric schedule the number of temperatures is the smallest integer \(K\) with \(T_0\alpha^K\le T_{\min}\), i.e., \(K=\big\lceil \ln(T_{\min}/T_0)/\ln\alpha\big\rceil\) when the ratio is not an integer, and the number of iterations is \(K\cdot L\). Other common rules stop after a number of temperatures without any accepted move, or at a time limit.

• Slow cooling (\(\alpha\) close to \(1\)) gives better solutions but more iterations: the parameters trade quality for time, as every heuristic does.

32.4.2 Examples

Example 32.3: \(SA\) with the 2-opt neighbourhood on the TSP reference instance

We start from the tour computed by the nearest neighbour heuristic from city \(1\) in the chapter «Constructive heuristics», \(\sigma=(1,3,2,4,5)\) of length \(25\), with \(T_0=5\), \(\alpha=0.8\), \(L=2\) and \(T_{\min}=3\): the temperatures are \(5\), \(4\) and \(3.2\) (then \(2.56\le 3\) stops the search), for a total of \(6\) iterations. The \(M=5\) moves of 2-opt on \(5\) cities are numbered

\[ 1:(1,3),\quad 2:(1,4),\quad 3:(2,4),\quad 4:(2,5),\quad 5:(3,5), \]

and the random numbers, with \(seed=1\), are those of Section 32.2.2: the first number of each iteration (\(u\)) chooses the move \(1+\lfloor 5u\rfloor\), the second (\(v\)) is used in the acceptance test.

it	\(T\)	\(u\)	move	\(\delta\)	\(v\)	\(e^{-\delta/T}\)	decision	current tour	\(L\)
\(1\)	\(5\)	\(0.0000\)	\(1:(1,3)\)	\(4\)	\(0.1315\)	\(0.4493\)	accept	\((1,2,3,4,5)\)	\(29\)
\(2\)	\(5\)	\(0.7556\)	\(4:(2,5)\)	\(-7\)	—	—	accept	\((1,2,5,4,3)\)	\(22\)
\(3\)	\(4\)	\(0.5328\)	\(3:(2,4)\)	\(2\)	\(0.2190\)	\(0.6065\)	accept	\((1,2,4,5,3)\)	\(24\)
\(4\)	\(4\)	\(0.0470\)	\(1:(1,3)\)	\(-6\)	—	—	accept	\((1,4,2,5,3)\)	\(18\)
\(5\)	\(3.2\)	\(0.6793\)	\(4:(2,5)\)	\(4\)	\(0.9347\)	\(0.2865\)	reject	\((1,4,2,5,3)\)	\(18\)
\(6\)	\(3.2\)	\(0.3835\)	\(2:(1,4)\)	\(4\)	\(0.5194\)	\(0.2865\)	reject	\((1,4,2,5,3)\)	\(18\)

At iterations \(1\) and \(3\) a worse tour is accepted (\(25\to 29\) and \(22\to 24\)), because \(v\) is smaller than the acceptance probability; at iterations \(5\) and \(6\) the temperature is lower, the worsening moves are less likely to be accepted, and both are rejected. The best tour, \((1,4,2,5,3)\), is the optimal tour \((1,3,5,2,4)\) traversed in the opposite direction: \(SA(I)=18=z(I)\).

• On this small instance the optimal tour is the only 2-opt local optimum among the \(12\) tours, so a 2-opt local search would reach it as well: the example shows the mechanics of \(SA\); its ability to escape from local optima shows on instances with bad local optima.

32.4.3 Infeasible solutions: simulated annealing with a penalty

• For the KP, the flip neighbourhood of a solution close to the capacity contains many infeasible solutions. Instead of discarding them, \(SA\) can visit them, with a penalty on the excess weight that makes them less attractive:

\[ f_\beta(x) \;=\; \sum_{j\in\mathscr{I}} p_j x_j \;-\; \beta\cdot \max\Big\{0,\ \sum_{j\in\mathscr{I}} w_j x_j - c\Big\}, \qquad \beta>0 . \]

\(SA\) maximizes \(f_\beta\) over all \(x\in\{0,1\}^n\); the worsening of a move from \(x\) to \(y\) is \(\delta=f_\beta(x)-f_\beta(y)\), and the best solution is updated with feasible solutions only.

• On a feasible solution \(f_\beta\) is the profit. Passing through an infeasible solution, \(SA\) can reach feasible solutions that are not reachable through feasible flips only: the infeasible region acts as a bridge between valleys.

• The penalty coefficient \(\beta\) is a parameter: if it is too small the search wanders in the infeasible region, where the profit is larger; if it is too large the infeasible solutions are never accepted, and the bridge is closed.

Example 32.4: \(SA\) with a penalty on the KP reference instance

We start from \(x=(1,0,1,0)\), the solution \(\{1,3\}\) of the greedy algorithm by efficiency, of profit \(16\). It is a local optimum for the flip neighbourhood (chapter «Heuristic algorithms: definitions»): its feasible neighbours have profit \(6\) and \(10\). We use \(\beta=3\), \(T_0=4\), \(\alpha=0.8\), \(L=2\), \(T_{\min}=2.5\) (temperatures \(4\), \(3.2\), \(2.56\): \(6\) iterations) and \(seed=3\). At each iteration \(u\) chooses the item \(k=1+\lfloor 4u\rfloor\) to flip, and \(v\) is used in the acceptance test.

it	\(T\)	\(u\)	\(k\)	\(y\)	\(\sum w_jy_j\)	\(f_3(y)\)	\(\delta\)	\(v\) vs \(e^{-\delta/T}\)	decision
\(1\)	\(4\)	\(0.0000\)	\(1\)	\((0,0,1,0)\)	\(3\)	\(6\)	\(10\)	\(0.3946\ge 0.0821\)	reject
\(2\)	\(4\)	\(0.2668\)	\(2\)	\((1,1,1,0)\)	\(12\)	\(14\)	\(2\)	\(0.3760<0.6065\)	accept
\(3\)	\(3.2\)	\(0.5983\)	\(3\)	\((1,1,0,0)\)	\(9\)	\(17\)	\(-3\)	—	accept
\(4\)	\(3.2\)	\(0.1411\)	\(1\)	\((0,1,0,0)\)	\(4\)	\(7\)	\(10\)	\(0.0366<0.0439\)	accept
\(5\)	\(2.56\)	\(0.0379\)	\(1\)	\((1,1,0,0)\)	\(9\)	\(17\)	\(-10\)	—	accept
\(6\)	\(2.56\)	\(0.1505\)	\(1\)	\((0,1,0,0)\)	\(4\)	\(7\)	\(10\)	\(0.5582\ge 0.0201\)	reject

At iteration \(2\) the search moves to the infeasible solution \((1,1,1,0)\), of weight \(12>9\): its profit is \(23\) and its penalized value \(23-3\cdot 3=14\), a worsening of \(2\) accepted with probability \(0.6065\). From there, removing item \(3\) gives the feasible solution \((1,1,0,0)\) of profit \(17\): the optimum, \(SA(I)=17=z(I)\). Iteration \(4\) shows that at this temperature even a worsening of \(10\) can be accepted, with a small probability; iteration \(5\) goes back to the optimum.

32.4.4 An experiment on the medium TSP instance

• We run \(SA\) with the 2-opt neighbourhood on TSP20 (Example 32.1), from the nearest neighbour tour of length \(326\), with \(L=100\) and \(T_{\min}=0.5\). The figure shows one run, with \(T_0=20\), \(\alpha=0.95\) and \(seed=1\): the length of the current tour (gray) and of the best tour (blue) along the \(7200\) iterations, every \(10\) iterations. The data are produced by the C++ code of the chapter, in its Colab notebook.

• At the beginning the temperature is high compared with the typical worsening of a 2-opt move, and the current tour gets much worse than the starting one (up to length \(765\)): the search is almost a random walk. As \(T\) decreases, fewer worsening moves are accepted and the current tour descends into a valley; at the end the search is a local search. The best tour, of length \(298\), is found at iteration \(4479\): it is better than the 2-opt local optimum (\(319\)), but it is not optimal.

• The table reports, for \(10\) runs with \(seed=1,\dots,10\) and several parameters, the number of iterations, the minimum, the average and the maximum length of the best tour, and the number of runs that reach the optimum \(291\).

\(\alpha\)	\(T_0\)	iterations	min	average	max	optimal runs
\(0.9\)	\(1\)	\(700\)	\(291\)	\(310.3\)	\(322\)	\(1\)
\(0.9\)	\(5\)	\(2200\)	\(291\)	\(297.3\)	\(305\)	\(4\)
\(0.9\)	\(20\)	\(3600\)	\(291\)	\(300.3\)	\(323\)	\(4\)
\(0.9\)	\(100\)	\(5100\)	\(291\)	\(297.8\)	\(311\)	\(3\)
\(0.95\)	\(1\)	\(1400\)	\(291\)	\(302.3\)	\(319\)	\(3\)
\(0.95\)	\(5\)	\(4500\)	\(291\)	\(294.3\)	\(298\)	\(5\)
\(0.95\)	\(20\)	\(7200\)	\(291\)	\(294.3\)	\(298\)	\(5\)
\(0.95\)	\(100\)	\(10400\)	\(291\)	\(293.2\)	\(299\)	\(7\)
\(0.99\)	\(1\)	\(6900\)	\(291\)	\(298.8\)	\(305\)	\(4\)
\(0.99\)	\(5\)	\(23000\)	\(291\)	\(291.7\)	\(298\)	\(9\)
\(0.99\)	\(20\)	\(36800\)	\(291\)	\(291.0\)	\(291\)	\(10\)
\(0.99\)	\(100\)	\(52800\)	\(291\)	\(291.0\)	\(291\)	\(10\)

• With \(T_0=1\) the temperature is too low from the beginning: few worsening moves are accepted, and \(SA\) behaves like a randomized local search. Slower cooling (\(\alpha=0.99\)) spends more iterations at every temperature range and reaches the optimum in almost every run, at the price of many more iterations. A single run says little: the comparison must look at the distribution over several seeds, as recalled in the chapter «Heuristic algorithms: definitions».

32.4.5 Properties

• At a fixed temperature \(T\), the sequence of current solutions of \(SA\) is a Markov chain: the next solution depends only on the current one. If the neighbourhood is symmetric (\(S'\in N(S)\) if and only if \(S\in N(S')\), all neighbourhoods having the same size) and connected (every solution can reach every other one through moves), the chain has a unique stationary distribution, the Boltzmann distribution:

\[ \pi_T(S) \;=\; \frac{e^{-f(S)/T}}{\sum_{S'} e^{-f(S')/T}} \qquad \text{for every solution } S. \]

• For two solutions \(S_1\) and \(S_2\), \(\frac{\pi_T(S_1)}{\pi_T(S_2)}=e^{(f(S_2)-f(S_1))/T}\): better solutions are more probable, and the more so the lower the temperature. As \(T\to 0\), \(\pi_T\) concentrates on the optimal solutions.

• Hajek (1988) proved that with the logarithmic schedule \(T_k=\frac{\Gamma}{\ln(k+1)}\), with \(\Gamma\) at least the depth of the deepest valley that does not contain an optimal solution, the probability that the current solution is optimal tends to \(1\). The result is a guarantee in the limit: the number of iterations needed can be larger than the number of solutions, and in practice the faster geometric schedule is used, without guarantees.

---

Next: 32.5 Tabu search