@cindex simulated annealing
@cindex combinatorial optimization
@cindex optimization, combinatorial
@cindex energy function
@cindex cost function
Stochastic search techniques are used when the structure of a space is
not well understood or is not smooth, so that techniques like Newton's
method (which requires calculating Jacobian derivative matrices) cannot
be used. In particular, these techniques are frequently used to solve
combinatorial optimization problems, such as the traveling salesman
problem.

The goal is to find a point in the space at which a real valued
@dfn{energy function} (or @dfn{cost function}) is minimized.  Simulated
annealing is a minimization technique which has given good results in
avoiding local minima; it is based on the idea of taking a random walk
through the space at successively lower temperatures, where the
probability of taking a step is given by a Boltzmann distribution.

The functions described in this chapter are declared in the header file
@file{gsl_siman.h}.

@menu
* Simulated Annealing algorithm::  
* Simulated Annealing functions::  
* Examples with Simulated Annealing::  
* Simulated Annealing References and Further Reading::  
@end menu

@node Simulated Annealing algorithm
@section Simulated Annealing algorithm

The simulated annealing algorithm takes random walks through the problem
space, looking for points with low energies; in these random walks, the
probability of taking a step is determined by the Boltzmann distribution,
@tex
\beforedisplay
$$
p = e^{-(E_{i+1} - E_i)/(kT)}
$$
\afterdisplay
@end tex
@ifinfo

@example
p = e^@{-(E_@{i+1@} - E_i)/(kT)@}
@end example

@end ifinfo
@noindent
if 
@c{$E_{i+1} > E_i$}
@math{E_@{i+1@} > E_i}, and 
@math{p = 1} when 
@c{$E_{i+1} \le E_i$}
@math{E_@{i+1@} <= E_i}.

In other words, a step will occur if the new energy is lower.  If
the new energy is higher, the transition can still occur, and its
likelihood is proportional to the temperature @math{T} and inversely
proportional to the energy difference 
@c{$E_{i+1} - E_i$}
@math{E_@{i+1@} - E_i}.

The temperature @math{T} is initially set to a high value, and a random
walk is carried out at that temperature.  Then the temperature is
lowered very slightly according to a @dfn{cooling schedule}, for
example: @c{$T \rightarrow T/\mu_T$}
@math{T -> T/mu_T}
where @math{\mu_T} is slightly greater than 1. 
@cindex cooling schedule
@cindex schedule, cooling

The slight probability of taking a step that gives higher energy is what
allows simulated annealing to frequently get out of local minima.

@node Simulated Annealing functions
@section Simulated Annealing functions

@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})

This function performs a simulated annealing search through a given
space.  The space is specified by providing the functions @var{Ef} and
@var{distance}.  The simulated annealing steps are generated using the
random number generator @var{r} and the function @var{take_step}.

The starting configuration of the system should be given by @var{x0_p}.
The routine offers two modes for updating configurations, a fixed-size
mode and a variable-size mode.  In the fixed-size mode the configuration
is stored as a single block of memory of size @var{element_size}.
Copies of this configuration are created, copied and destroyed
internally using the standard library functions @code{malloc},
@code{memcpy} and @code{free}.  The function pointers @var{copyfunc},
@var{copy_constructor} and @var{destructor} should be null pointers in
fixed-size mode.  In the variable-size mode the functions
@var{copyfunc}, @var{copy_constructor} and @var{destructor} are used to
create, copy and destroy configurations internally.  The variable
@var{element_size} should be zero in the variable-size mode.

The @var{params} structure (described below) controls the run by
providing the temperature schedule and other tunable parameters to the
algorithm.

On exit the best result achieved during the search is placed in
@code{*@var{x0_p}}.  If the annealing process has been successful this
should be a good approximation to the optimal point in the space.

If the function pointer @var{print_position} is not null, a debugging
log will be printed to @code{stdout} with the following columns:

@example
number_of_iterations temperature x x-(*x0_p) Ef(x)
@end example

and the output of the function @var{print_position} itself.  If
@var{print_position} is null then no information is printed.
@end deftypefun

@noindent
The simulated annealing routines require several user-specified
functions to define the configuration space and energy function.  The
prototypes for these functions are given below.

@deftp {Data Type} gsl_siman_Efunc_t
This function type should return the energy of a configuration @var{xp}.

@example
double (*gsl_siman_Efunc_t) (void *xp)
@end example
@end deftp

@deftp {Data Type} gsl_siman_step_t
This function type should modify the configuration @var{xp} using a random step
taken from the generator @var{r}, up to a maximum distance of
@var{step_size}.

@example
void (*gsl_siman_step_t) (const gsl_rng *r, void *xp, 
                          double step_size)
@end example
@end deftp

@deftp {Data Type} gsl_siman_metric_t
This function type should return the distance between two configurations
@var{xp} and @var{yp}.

@example
double (*gsl_siman_metric_t) (void *xp, void *yp)
@end example
@end deftp

@deftp {Data Type} gsl_siman_print_t
This function type should print the contents of the configuration @var{xp}.

@example
void (*gsl_siman_print_t) (void *xp)
@end example
@end deftp

@deftp {Data Type} gsl_siman_copy_t
This function type should copy the configuration @var{source} into @var{dest}.

@example
void (*gsl_siman_copy_t) (void *source, void *dest)
@end example
@end deftp

@deftp {Data Type} gsl_siman_copy_construct_t
This function type should create a new copy of the configuration @var{xp}.

@example
void * (*gsl_siman_copy_construct_t) (void *xp)
@end example
@end deftp

@deftp {Data Type} gsl_siman_destroy_t
This function type should destroy the configuration @var{xp}, freeing its
memory.

@example
void (*gsl_siman_destroy_t) (void *xp)
@end example
@end deftp

@deftp {Data Type} gsl_siman_params_t
These are the parameters that control a run of @code{gsl_siman_solve}.
This structure contains all the information needed to control the
search, beyond the energy function, the step function and the initial
guess.

@table @code
@item int n_tries          
The number of points to try for each step.

@item int iters_fixed_T    
The number of iterations at each temperature.

@item double step_size     
The maximum step size in the random walk.

@item double k, t_initial, mu_t, t_min
The parameters of the Boltzmann distribution and cooling schedule.
@end table
@end deftp


@node Examples with Simulated Annealing
@section Examples

The simulated annealing package is clumsy, and it has to be because it
is written in C, for C callers, and tries to be polymorphic at the same
time.  But here we provide some examples which can be pasted into your
application with little change and should make things easier.

@menu
* Trivial example::             
* Traveling Salesman Problem::  
@end menu

@node Trivial example
@subsection Trivial example

The first example, in one dimensional cartesian space, sets up an energy
function which is a damped sine wave; this has many local minima, but
only one global minimum, somewhere between 1.0 and 1.5.  The initial
guess given is 15.5, which is several local minima away from the global
minimum.

@smallexample
@verbatiminclude examples/siman.c
@end smallexample

@need 2000
Here are a couple of plots that are generated by running
@code{siman_test} in the following way:

@example
$ ./siman_test | awk '!/^#/ @{print $1, $4@}' 
 | graph -y 1.34 1.4 -W0 -X generation -Y position 
 | plot -Tps > siman-test.eps
$ ./siman_test | awk '!/^#/ @{print $1, $4@}' 
 | graph -y -0.88 -0.83 -W0 -X generation -Y energy 
 | plot -Tps > siman-energy.eps
@end example

@iftex
@sp 1
@center @image{siman-test,2.8in} 
@center @image{siman-energy,2.8in}

@quotation
Example of a simulated annealing run: at higher temperatures (early in
the plot) you see that the solution can fluctuate, but at lower
temperatures it converges.
@end quotation
@end iftex

@node Traveling Salesman Problem
@subsection Traveling Salesman Problem
@cindex TSP
@cindex traveling salesman problem

The TSP (@dfn{Traveling Salesman Problem}) is the classic combinatorial
optimization problem.  I have provided a very simple version of it,
based on the coordinates of twelve cities in the southwestern United
States.  This should maybe be called the @dfn{Flying Salesman Problem},
since I am using the great-circle distance between cities, rather than
the driving distance.  Also: I assume the earth is a sphere, so I don't
use geoid distances.

The @code{gsl_siman_solve} routine finds a route which is 3490.62
Kilometers long; this is confirmed by an exhaustive search of all
possible routes with the same initial city.

The full code can be found in @file{siman/siman_tsp.c}, but I include
here some plots generated in the following way:

@smallexample
$ ./siman_tsp > tsp.output
$ grep -v "^#" tsp.output  
 | awk '@{print $1, $NF@}'
 | graph -y 3300 6500 -W0 -X generation -Y distance 
    -L "TSP - 12 southwest cities"
 | plot -Tps > 12-cities.eps
$ grep initial_city_coord tsp.output 
  | awk '@{print $2, $3@}' 
  | graph -X "longitude (- means west)" -Y "latitude" 
     -L "TSP - initial-order" -f 0.03 -S 1 0.1 
  | plot -Tps > initial-route.eps
$ grep final_city_coord tsp.output 
  | awk '@{print $2, $3@}' 
  | graph -X "longitude (- means west)" -Y "latitude" 
     -L "TSP - final-order" -f 0.03 -S 1 0.1 
  | plot -Tps > final-route.eps
@end smallexample

@noindent
This is the output showing the initial order of the cities; longitude is
negative, since it is west and I want the plot to look like a map.

@smallexample
# initial coordinates of cities (longitude and latitude)
###initial_city_coord: -105.95 35.68 Santa Fe
###initial_city_coord: -112.07 33.54 Phoenix
###initial_city_coord: -106.62 35.12 Albuquerque
###initial_city_coord: -103.2 34.41 Clovis
###initial_city_coord: -107.87 37.29 Durango
###initial_city_coord: -96.77 32.79 Dallas
###initial_city_coord: -105.92 35.77 Tesuque
###initial_city_coord: -107.84 35.15 Grants
###initial_city_coord: -106.28 35.89 Los Alamos
###initial_city_coord: -106.76 32.34 Las Cruces
###initial_city_coord: -108.58 37.35 Cortez
###initial_city_coord: -108.74 35.52 Gallup
###initial_city_coord: -105.95 35.68 Santa Fe
@end smallexample

The optimal route turns out to be:

@smallexample
# final coordinates of cities (longitude and latitude)
###final_city_coord: -105.95 35.68 Santa Fe
###final_city_coord: -106.28 35.89 Los Alamos
###final_city_coord: -106.62 35.12 Albuquerque
###final_city_coord: -107.84 35.15 Grants
###final_city_coord: -107.87 37.29 Durango
###final_city_coord: -108.58 37.35 Cortez
###final_city_coord: -108.74 35.52 Gallup
###final_city_coord: -112.07 33.54 Phoenix
###final_city_coord: -106.76 32.34 Las Cruces
###final_city_coord: -96.77 32.79 Dallas
###final_city_coord: -103.2 34.41 Clovis
###final_city_coord: -105.92 35.77 Tesuque
###final_city_coord: -105.95 35.68 Santa Fe
@end smallexample

@iftex
@sp 1
@center @image{initial-route,2.2in} 
@center @image{final-route,2.2in}

@quotation
Initial and final (optimal) route for the 12 southwestern cities Flying
Salesman Problem.
@end quotation
@end iftex

@noindent
Here's a plot of the cost function (energy) versus generation (point in
the calculation at which a new temperature is set) for this problem:

@iftex
@sp 1
@center @image{12-cities,2.8in}

@quotation
Example of a simulated annealing run for the 12 southwestern cities
Flying Salesman Problem.
@end quotation
@end iftex

@node Simulated Annealing References and Further Reading
@section References and Further Reading

Further information is available in the following book,

@itemize @asis
@item
@cite{Modern Heuristic Techniques for Combinatorial Problems}, Colin R. Reeves
(ed.), McGraw-Hill, 1995 (ISBN 0-07-709239-2).
@end itemize