htswanalysis/MACS/lib/gsl/gsl-1.11/doc/siman.texi

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