1.. index::
2   single: simulated annealing
3   single: combinatorial optimization
4   single: optimization, combinatorial
5   single: energy function
6   single: cost function
7
8*******************
9Simulated Annealing
10*******************
11
12Stochastic search techniques are used when the structure of a space is
13not well understood or is not smooth, so that techniques like Newton's
14method (which requires calculating Jacobian derivative matrices) cannot
15be used. In particular, these techniques are frequently used to solve
16combinatorial optimization problems, such as the traveling salesman
17problem.
18
19The goal is to find a point in the space at which a real valued
20*energy function* (or *cost function*) is minimized.  Simulated
21annealing is a minimization technique which has given good results in
22avoiding local minima; it is based on the idea of taking a random walk
23through the space at successively lower temperatures, where the
24probability of taking a step is given by a Boltzmann distribution.
25
26The functions described in this chapter are declared in the header file
27:file:`gsl_siman.h`.
28
29Simulated Annealing algorithm
30=============================
31
32The simulated annealing algorithm takes random walks through the problem
33space, looking for points with low energies; in these random walks, the
34probability of taking a step is determined by the Boltzmann distribution,
35
36.. math:: p = e^{-(E_{i+1} - E_i)/(kT)}
37
38if
39:math:`E_{i+1} > E_i`, and
40:math:`p = 1` when
41:math:`E_{i+1} \le E_i`.
42
43In other words, a step will occur if the new energy is lower.  If
44the new energy is higher, the transition can still occur, and its
45likelihood is proportional to the temperature :math:`T` and inversely
46proportional to the energy difference
47:math:`E_{i+1} - E_i`.
48
49.. index::
50   single: cooling schedule
51   single: schedule, cooling
52
53The temperature :math:`T` is initially set to a high value, and a random
54walk is carried out at that temperature.  Then the temperature is
55lowered very slightly according to a *cooling schedule*, for
56example: :math:`T \rightarrow T/\mu_T`
57where :math:`\mu_T` is slightly greater than 1.
58
59The slight probability of taking a step that gives higher energy is what
60allows simulated annealing to frequently get out of local minima.
61
62Simulated Annealing functions
63=============================
64
65.. function:: void gsl_siman_solve (const gsl_rng * r, void * x0_p, gsl_siman_Efunc_t Ef, gsl_siman_step_t take_step, gsl_siman_metric_t distance, gsl_siman_print_t print_position, gsl_siman_copy_t copyfunc, gsl_siman_copy_construct_t copy_constructor, gsl_siman_destroy_t destructor, size_t element_size, gsl_siman_params_t params)
66
67   This function performs a simulated annealing search through a given
68   space.  The space is specified by providing the functions :data:`Ef` and
69   :data:`distance`.  The simulated annealing steps are generated using the
70   random number generator :data:`r` and the function :data:`take_step`.
71
72   The starting configuration of the system should be given by :data:`x0_p`.
73   The routine offers two modes for updating configurations, a fixed-size
74   mode and a variable-size mode.  In the fixed-size mode the configuration
75   is stored as a single block of memory of size :data:`element_size`.
76   Copies of this configuration are created, copied and destroyed
77   internally using the standard library functions :func:`malloc`,
78   :func:`memcpy` and :func:`free`.  The function pointers :data:`copyfunc`,
79   :data:`copy_constructor` and :data:`destructor` should be null pointers in
80   fixed-size mode.  In the variable-size mode the functions
81   :data:`copyfunc`, :data:`copy_constructor` and :data:`destructor` are used to
82   create, copy and destroy configurations internally.  The variable
83   :data:`element_size` should be zero in the variable-size mode.
84
85   The :data:`params` structure (described below) controls the run by
86   providing the temperature schedule and other tunable parameters to the
87   algorithm.
88
89   On exit the best result achieved during the search is placed in
90   :data:`x0_p`.  If the annealing process has been successful this
91   should be a good approximation to the optimal point in the space.
92
93   If the function pointer :data:`print_position` is not null, a debugging
94   log will be printed to :code:`stdout` with the following columns::
95
96      #-iter  #-evals  temperature  position  energy  best_energy
97
98   and the output of the function :data:`print_position` itself.  If
99   :data:`print_position` is null then no information is printed.
100
101The simulated annealing routines require several user-specified
102functions to define the configuration space and energy function.  The
103prototypes for these functions are given below.
104
105.. type:: gsl_siman_Efunc_t
106
107   This function type should return the energy of a configuration :data:`xp`::
108
109      double (*gsl_siman_Efunc_t) (void *xp)
110
111.. type:: gsl_siman_step_t
112
113   This function type should modify the configuration :data:`xp` using a random step
114   taken from the generator :data:`r`, up to a maximum distance of
115   :data:`step_size`::
116
117      void (*gsl_siman_step_t) (const gsl_rng *r, void *xp,
118                                double step_size)
119
120.. type:: gsl_siman_metric_t
121
122   This function type should return the distance between two configurations
123   :data:`xp` and :data:`yp`::
124
125      double (*gsl_siman_metric_t) (void *xp, void *yp)
126
127.. type:: gsl_siman_print_t
128
129   This function type should print the contents of the configuration :data:`xp`::
130
131      void (*gsl_siman_print_t) (void *xp)
132
133.. type:: gsl_siman_copy_t
134
135   This function type should copy the configuration :data:`source` into :data:`dest`::
136
137      void (*gsl_siman_copy_t) (void *source, void *dest)
138
139.. type:: gsl_siman_copy_construct_t
140
141   This function type should create a new copy of the configuration :data:`xp`::
142
143      void * (*gsl_siman_copy_construct_t) (void *xp)
144
145.. type:: gsl_siman_destroy_t
146
147   This function type should destroy the configuration :data:`xp`, freeing its
148   memory::
149
150      void (*gsl_siman_destroy_t) (void *xp)
151
152.. type:: gsl_siman_params_t
153
154   These are the parameters that control a run of :func:`gsl_siman_solve`.
155   This structure contains all the information needed to control the
156   search, beyond the energy function, the step function and the initial
157   guess.
158
159   ========================================= ============================================================
160   :code:`int n_tries`                       The number of points to try for each step.
161   :code:`int iters_fixed_T`                 The number of iterations at each temperature.
162   :code:`double step_size`                  The maximum step size in the random walk.
163   :code:`double k, t_initial, mu_t, t_min`  The parameters of the Boltzmann distribution and cooling
164                                             schedule.
165   ========================================= ============================================================
166
167Examples
168========
169
170The simulated annealing package is clumsy, and it has to be because it
171is written in C, for C callers, and tries to be polymorphic at the same
172time.  But here we provide some examples which can be pasted into your
173application with little change and should make things easier.
174
175Trivial example
176---------------
177
178The first example, in one dimensional Cartesian space, sets up an energy
179function which is a damped sine wave; this has many local minima, but
180only one global minimum, somewhere between 1.0 and 1.5.  The initial
181guess given is 15.5, which is several local minima away from the global
182minimum.
183
184.. include:: examples/siman.c
185   :code:
186
187:numref:`fig_siman-test` is generated by running
188:code:`siman_test` in the following way::
189
190  $ ./siman_test | awk '!/^#/ {print $1, $4}'
191   | graph -y 1.34 1.4 -W0 -X generation -Y position
192   | plot -Tps > siman-test.eps
193
194:numref:`fig_siman-energy` is generated by running
195:code:`siman_test` in the following way::
196
197  $ ./siman_test | awk '!/^#/ {print $1, $5}'
198   | graph -y -0.88 -0.83 -W0 -X generation -Y energy
199   | plot -Tps > siman-energy.eps
200
201.. _fig_siman-test:
202
203.. figure:: /images/siman-test.png
204   :scale: 60%
205
206   Example of a simulated annealing run: at higher temperatures (early in
207   the plot) you see that the solution can fluctuate, but at lower
208   temperatures it converges.
209
210.. _fig_siman-energy:
211
212.. figure:: /images/siman-energy.png
213   :scale: 60%
214
215   Simulated annealing energy vs generation
216
217.. index::
218   single: TSP
219   single: traveling salesman problem
220
221Traveling Salesman Problem
222--------------------------
223
224The TSP (*Traveling Salesman Problem*) is the classic combinatorial
225optimization problem.  I have provided a very simple version of it,
226based on the coordinates of twelve cities in the southwestern United
227States.  This should maybe be called the *Flying Salesman Problem*,
228since I am using the great-circle distance between cities, rather than
229the driving distance.  Also: I assume the earth is a sphere, so I don't
230use geoid distances.
231
232The :func:`gsl_siman_solve` routine finds a route which is 3490.62
233Kilometers long; this is confirmed by an exhaustive search of all
234possible routes with the same initial city.
235
236The full code is given below.
237
238.. include:: examples/siman_tsp.c
239   :code:
240
241Below are some plots generated in the following way::
242
243  $ ./siman_tsp > tsp.output
244  $ grep -v "^#" tsp.output
245   | awk '{print $1, $NF}'
246   | graph -y 3300 6500 -W0 -X generation -Y distance
247      -L "TSP - 12 southwest cities"
248   | plot -Tps > 12-cities.eps
249  $ grep initial_city_coord tsp.output
250    | awk '{print $2, $3}'
251    | graph -X "longitude (- means west)" -Y "latitude"
252       -L "TSP - initial-order" -f 0.03 -S 1 0.1
253    | plot -Tps > initial-route.eps
254  $ grep final_city_coord tsp.output
255    | awk '{print $2, $3}'
256    | graph -X "longitude (- means west)" -Y "latitude"
257       -L "TSP - final-order" -f 0.03 -S 1 0.1
258    | plot -Tps > final-route.eps
259
260This is the output showing the initial order of the cities; longitude is
261negative, since it is west and I want the plot to look like a map::
262
263  # initial coordinates of cities (longitude and latitude)
264  ###initial_city_coord: -105.95 35.68 Santa Fe
265  ###initial_city_coord: -112.07 33.54 Phoenix
266  ###initial_city_coord: -106.62 35.12 Albuquerque
267  ###initial_city_coord: -103.2 34.41 Clovis
268  ###initial_city_coord: -107.87 37.29 Durango
269  ###initial_city_coord: -96.77 32.79 Dallas
270  ###initial_city_coord: -105.92 35.77 Tesuque
271  ###initial_city_coord: -107.84 35.15 Grants
272  ###initial_city_coord: -106.28 35.89 Los Alamos
273  ###initial_city_coord: -106.76 32.34 Las Cruces
274  ###initial_city_coord: -108.58 37.35 Cortez
275  ###initial_city_coord: -108.74 35.52 Gallup
276  ###initial_city_coord: -105.95 35.68 Santa Fe
277
278The optimal route turns out to be::
279
280  # final coordinates of cities (longitude and latitude)
281  ###final_city_coord: -105.95 35.68 Santa Fe
282  ###final_city_coord: -103.2 34.41 Clovis
283  ###final_city_coord: -96.77 32.79 Dallas
284  ###final_city_coord: -106.76 32.34 Las Cruces
285  ###final_city_coord: -112.07 33.54 Phoenix
286  ###final_city_coord: -108.74 35.52 Gallup
287  ###final_city_coord: -108.58 37.35 Cortez
288  ###final_city_coord: -107.87 37.29 Durango
289  ###final_city_coord: -107.84 35.15 Grants
290  ###final_city_coord: -106.62 35.12 Albuquerque
291  ###final_city_coord: -106.28 35.89 Los Alamos
292  ###final_city_coord: -105.92 35.77 Tesuque
293  ###final_city_coord: -105.95 35.68 Santa Fe
294
295.. figure:: /images/siman-initial-route.png
296   :scale: 60%
297
298   Initial route for the 12 southwestern cities Flying Salesman Problem.
299
300.. figure:: /images/siman-final-route.png
301   :scale: 60%
302
303   Final (optimal) route for the 12 southwestern cities Flying Salesman Problem.
304
305Here's a plot of the cost function (energy) versus generation (point in
306the calculation at which a new temperature is set) for this problem:
307
308.. figure:: /images/siman-12-cities.png
309   :scale: 60%
310
311   Example of a simulated annealing run for the 12 southwestern cities
312   Flying Salesman Problem.
313
314References and Further Reading
315==============================
316
317Further information is available in the following book,
318
319* *Modern Heuristic Techniques for Combinatorial Problems*, Colin R. Reeves
320  (ed.), McGraw-Hill, 1995 (ISBN 0-07-709239-2).
321