\chapter{Simulated Annealing and Energy Landscapes}\label{ch:SA}

Simulated annealing (SA) -- also known as Monte Carlo annealing,
statistical cooling, probabilistic hill-climbing, stochastic
relaxation, and probabilistic exchange ``algorithm'' -- has been
widely used for tackling different combinatorial optimization
problems, including academic scheduling~\cite{TD:cooling,
A:timetable, V:applied}. The basic algorithm is outlined in
Figure~\ref{SA-algorithm}.

As with any search algorithm, simulated annealing 
requires the answers for the following problem-specific questions:

\begin{itemize}
\item What is a solution?
\item What are the neighbors of a solution?
\item What is the cost of a solution?
\item How do we determine the initial solution?
\end{itemize}

These answers yield the structure of the search space together with
the definition of a neighborhood, the evaluation function, and the
initial starting point. Note, however, that simulated annealing also
requires answers for additional questions:

\begin{itemize}
\item How do we determine the initial temperature T?
\item How do we determine the cooling ratio of T (or the cooling schedule)?
\item How do we determine the length of each homogeneous Markov chain?
\item How do we determine the halting criterion?
\end{itemize}

Before talking about possible answers to these questions, first we
would like to briefly outline our intuition about annealing. It is
logically clear to us that the success of annealing in many ways
depends on both the cost function to which it is applied and the 
move set used. First, imagine the move set used is that of the complete
graph. Then by randomly generating moves, annealing is randomly
generating states. At the least, it could keep the best random state
found so far; even that would be a terrible algorithm. But if the
temperature is nonzero it does still worse, since it may relinquish an 
optimal state. On the other hand, if the move set is terribly sparse,
it may simply be impossible to get from a bad state to a good one in a 
reasonable number of moves. Finally, imagine that a ``good'' move set
has been chosen, but energies have been randomly assigned to the
states. In that case no information is known about a state before it
is visited, and the only way to find a particularly good state is to
visit many states. 

Note that this is not the case for, say, a smooth energy function. For 
a smooth energy function, bounds or expectations of all nearby states
are known from the current state alone, and it might be possible to
find a good state by visiting only a small fraction of the set of all
states. Of course, ``smooth'' in this context means that the energy
changes relatively little from move to move. So, for annealing to work 
well, the energies and the move set must be well matched to one
another. Therefore, the key question to ask is that, intuitively, what 
properties must the landscape have if annealing is to be efficient?

\section{Energy Landscapes}
In the statistical physics context, energy landscapes are 
visualizations of the control functions that determine the dynamics of 
multi-body (or many-body) systems. At one level, they describe true
energies in a high-dimensional space in which each ``grid-point''
specifies a complete micro-state and the micro-dynamics is
stochastically directed through the changes in energy. Alternatively,
the landscapes can describe free energies as a function of
macro-states. In other manifestations they can also be conceptual
guides to dynamics with detailed balance. Both microscopic and
macroscopic constraints can however lead to motion on the landscapes
which is more complex than that suggested by gradient descent. The
landscape structures are determined by global as well as local
interactions between the microscopic constituents with external
agencies. Energy landscapes come in different shapes; there are
smooth landscapes and there are rough, with many hills, valleys,
saddles and other non-trivial topology, describing correspondingly
non-trivial dynamics, flows and attractors. This is usually referred
to as `complex systems'. 

Such models typify systems with quenched disorder and
frustration in the microscopic rules of engagement, and also include
ones in which disorder is spontaneously self-induced through mutual
interaction. Examples occur in spin glasses, regular glasses,
proteins and many other condensed matter systems, where energy has its
usual Hamiltonian meaning and the extent of the stochasticity is
determined by the temperature. Conceptually and mathematically,
however, the description can also be employed fruitfully to analyze
multi-body systems outside the normal realm of physics, for example,
neural networks, hard optimization problems, evolution and many
problems involving large number of interacting agents and
constraints. In these cases `energy' and `temperature' are defined by
analogy. In systems with such rugged landscapes, it is necessary to
reconsider much of conventional statistical physics.

On the issue of the abovementioned dynamics, they need not be
restricted to motions on the landscapes. The landscapes themselves can 
evolve, either via changes in global controls or through a dynamical
evolution of the interactions between the particles of the system. The 
latter occurs, for example, in the modification of synapses in neural
learning (where the motion on the landscape represents the faster
process of neural association or recall), and in the slow evolutionary 
selection of biological and ecological species (whose reproduction is
determined by the faster motion of their micro-entities on their
energy or fitness landscape). These considerations then lead to
further considerations of co-evolving systems of agents and rules,
with many potentials for complex behavior.

\section{Good and Bad Energy Landscapes}

\begin{figure}[tbp]
\begin{center}
\begin{picture}(370,250)
\psfig{figure=bad-landscape.eps,height=370pt,width=250pt,angle=-90}
\end{picture}
\end{center}
\caption{A landscape that is bad for annealing.}
\label{bad-landscape}
\end{figure}

\begin{figure}[tbp]
\begin{center}
\begin{picture}(370,225)
\psfig{figure=good-landscape.eps,height=225pt,width=400pt}
\end{picture}
\end{center}
\caption{A landscape that is good for annealing.}
\label{good-landscape}
\end{figure}

Suppose that the energy landscape is as sketched in the ``bad''
landscape of Figure~\ref{bad-landscape}. States are represented 
along the x axis, with adjacent states being neighbors. On this
landscape, the energy differences between the low-energy states (the
valley bottoms) are fairly small, while the energy barriers separating
them (the mountains) are large. It is well known that the time
required to cross a barrier of height $h$ at temperature $T$ is
exponential in $h/T$, while the probability of a state of energy $f$
is exponential in $-f/T$. So crossing the high barriers in reasonable
time demands $T$ to be large, while favoring the better valleys over the
less good ones requires $T$ to be small, implying that annealing
cannot work both well and quickly on this space.

On the other hand, annealing should work well on a function like that
of the ``good'' landscape of Figure~\ref{good-landscape}. In this
case, annealing can work in a hierarchical fashion. Initially, the
goal can be just to choose the better (left) of the two valleys
separated by the tallest barrier $B_1$. While $B_1$ is large, the
overall energy difference $\Delta f_1$ of the two valleys is also
fairly large, so with a suitable value of $T_1$ we can be in
the lower-energy valley with high probability in short time.

Next we can aim to settle in the better (right) of the two valleys
separated by the smaller barrier $B_2$. While the energy difference
$\Delta f_2$ between these valleys is smaller than $\Delta f_1$ was,
$B_2$ is also smaller than $B_1$ and so using $T_2$ similarly smaller
than $T_1$ again allows us to be in the lower-energy valley with high 
probability in reasonable time. This process is repeated for smaller
energy scales. 

The success of annealing relies on the overall energy difference of
collection of states being large compared with the barriers dividing
these collections. So we would conjecture that either all the energy
barriers and energy differences are of the same scale (which really
does not seem to be the case in practice), or else in smaller and
smaller areas of the landscape the energy barriers must scale down
along with the energy differences, giving a general
``self-similarity'' or ``fractalness'' like that of the ``good''
landscape of Figure~\ref{good-landscape}.

\section{On the Application of Annealing}
Depending on the problem tackled, we observed that how the cooling
schedule is chosen plays a major part in shaping the quality of the
results obtained for academic scheduling. Initially, we used the most
commonly known and used schedule, which is the geometric cooling, but
later tried adaptive cooling, as well as the method of geometric
reheating based on cost \cite{KGV:sim, K:quant} and \cite{ADK:cooling}.

A comprehensive discussion of the theoretical and practical details of
SA is given in~\cite{AKvL:anneal, OG:annealing, vLA:anneal}. It
suffices here to say that the elementary operation in the Metropolis
method for a combinatorial problem such as scheduling is the
generation of some new candidate configuration, which is then
automatically accepted if it lowers the cost ($C$), or accepted with
probability $\exp(-\Delta C /T)$, where $T$ is the temperature, if it
would increase the cost by $\Delta (C)$. Also, in
Figure~\ref{SA-algorithm}, $s$ is the current schedule and $s^{'}$ is
a neighboring schedule obtained from the current neighborhood  space
${\cal N}_{s}$ by swapping two classes in time and/or space. 

Thus the technique is essentially a generalization of the local optimization
strategy, where, at non-zero temperatures, thermal excitations can
facilitate escape from local minima. 

\begin{figure}
\begin{center}
\begin{enumerate}
\item Generate an initial schedule $s$.
\item Set the initial best schedule $s^{*} = s$.
\item Compute cost of $s$ :  $C(s)$.
\item Compute initial temperature $T_0$.
\item Set the temperature $T\; =\; T_0$.
\item While {\em stop criterion} is not satisfied do:
\begin{enumerate}
\item Repeat {\em Markov chain length} (M) times:
\begin{enumerate}
\item Select a random neighbor $s^{'}$ to the current schedule, 
$( s^{'} \subset {\cal N}_{s} )\;$.
\item Set $\Delta(C) = C(s^{'}) - C(s)\;$.
\item If ($\Delta(C) \leq 0$ $\{$downhill move$\}$):
\begin{itemize}
\item Set $s = s^{'}\;$.
\item If $C(s) < C(s^{*})$ then set $s^{*} = s$.
\end{itemize}
\item If ($\Delta(C) > 0$ $\{$uphill move$\}$):
\begin{itemize}
\item Choose a random number $r$ uniformly from $[0,1]$.
\item If $r < e^{-\Delta(C) /T}\;$ then set $s = s^{'}$.
\end{itemize}
\end{enumerate}
\item Reduce (or update) temperature $T$.
\end{enumerate}
\item Return the schedule  $s^{*}$.
\end{enumerate}
\end{center}
\caption{The Simulated Annealing Algorithm}
\label{SA-algorithm}
\end{figure}

The SA algorithm has advantages and disadvantages compared to other
global optimization techniques. 
Among its advantages are the relative ease of implementation,
the applicability to almost any combinatorial optimization problem, 
the ability to provide reasonably good solutions for most problems
(depending on the cooling schedule and update moves used), 
and the ease with which it can be combined with other heuristics, such as 
expert systems, forming quite useful hybrid methods for tackling a range of 
complex problems.
SA is a robust technique, however, it does have some drawbacks.
To obtain good results the update moves and the various tunable parameters 
used (such as the cooling rate) need to be carefully chosen, 
the runs often require a great deal of computer time, and many runs may be
required.

Using the search space (or the landscape) metaphor, an intuitive
picture of how it works is as follows: the shape of the landscape is
randomly smoothed or blurred at high temperatures, so that local peaks 
don't persist for too long. The solution can easily move to the fitter 
parts of the landscape. As the temperature is gradually reduced, the
peaks begin to stick together, but the solution has escaped from their 
clutches and converges gradually to the top of the highest (or a high) 
peak.


\begin{figure}[tbp]
\begin{center}
\begin{picture}(320,250)
\psfig{figure=state-energy.eps,height=220pt,width=350pt}
\end{picture}
\end{center}
\caption{Illustrating a problem for which annealing won't work very
well. The state here is a one-dimensional, with the energy of a state 
as shown. At low temperature an energy barrier confines the simulation 
for long periods of time to either the small region on the left or the 
larger region on the right. The small region contains the global
minimum, but annealing will likely direct the simulation to the larger 
region, since the tiny region in the vicinity of the global minimum has 
negligible influence on the distribution at high temperatures where 
the energy barrier can easily be traversed. (See also Figure~\ref{bad-landscape})}
\label{state-energy}
\end{figure}

As shown in Figure~\ref{state-energy}, the distribution at high
temperature is dominated by the size of a two-energy regions, not by
the depths of their lowest points. Since the global minimum happens to
be in the smaller region, annealing will in this case direct the
simulation away from the region where the equilibrium distribution
will be concentrated at very low temperatures. To avoid such a
situation we would need to be very careful about the kind of annealing
schedule we choose.

In general, depending on the problem to which it is applied, SA
has been shown to be quite competitive with many of the best
heuristics, as shown in the work of Johnson {\it et
al.}~\cite{JM:traveling}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Timetabling Using the Annealing Algorithm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The most obvious mapping of the timetabling problem into the SA algorithm
involves the following constructs:

\begin{enumerate}
\item a {\bf state} is a timetable containing the following sets:
\begin{itemize}
\item $P$: a set of professors.
\item $C$: a set of classes.
\item $S$: a set of students.
\item $R$: a set of classrooms.
\item $I$: a set of time intervals.
\end{itemize}
\item a {\bf cost} or ``energy'' E($P, C, S, R, I$) such that:
\begin{itemize}

\item E($P$): is the cost of assigning more than maximum number of 
allowed classes $M_p$ to the same professor, 
plus scheduling one or more classes that cause a conflict
in the professor's schedule.

\item E($C$): is the cost of scheduling certain classes at/within the
same time period in violation of the exclusion constraint, for
example.

\item E($S$): is the cost of having two or more classes conflict in
time; plus the cost of having in the schedule one or more classes that
really don't meet the student's major, class requested, or class
requirements; plus the cost of not having the classes evenly spread
out over the week, etc. 

\item E($R$): is the cost resulting from assigning room(s) of the wrong
size and/or type to a certain class.

\item E($I$): is the cost of having more or less time periods than
required, plus the cost of imbalanced class assignments (a certain
period will have more classes assigned than others, etc.).
\end{itemize}

\item A {\bf swap} (or a move) is the exchange of one or more of the
following:
class $c_i$ with class $c_j$ in the set $C$ with
respect to time periods $I_i$ and $I_j$, and/or with
respect to classroom $R_i$ and $R_j$, respectively. Generally, this
step is referred to as class swapping.
\end{enumerate}

Along with all of the necessary constraints, the simulated annealing
algorithm also takes as input data the following:  the preprocessor
output in the form of lists of scheduled and non-scheduled classes and
their associated professors and room types, a list of rooms provided
by the registrar's office, a department to building distance matrix, a
list of students and their class preferences, a list of all
departments and their building preferences, and a list of classes
that are not allowed to be scheduled simultaneously.


\section{Interpretation of the Solution Space}

After setting up or mapping the problem we need to determine the
following: 
\begin{itemize}
\item 
A {\it feasible solution} is defined as a schedule that
satisfies a subset ($C^{'}$) of constraints ($C$) : $C^{'} \subseteq C\;$;
provided that all hard constraints are part of $C^{'}\;$.
\item 
A {\it neighboring solution} is obtained by modifying the current solution
space (i.e. swapping two or more classes in time or place).
\item
Here are two guidelines in determining whether a given constraint to be
included in $C^{'}$ or not:
\begin{enumerate}
\item A class should never be assigned to a time period if such an
assignment induces only class schedules that do not satisfy all the
constraints (specifically the hard constraints).
\item A conflict between two or more classes should not be a
sufficient condition for hindering an assignment. For example,
let us consider two courses $c_1$ and $c_2$ involving common
students. If $c_1$ cannot be scheduled at period $p$ for some
reason (maybe the teacher is not available at that time), we will
never visit a solution where $c_1$ is given at time $p$. However, if this
constraint does not exist and $c_2$ is scheduled at time $p$, we are
allowed to schedule $c_1$ at time $p$, even if this violates the
constraint of no overlapping of courses involving common students. 
The reason is that $c_1$ is possibly given at time $p$ in an optimal
schedule and the violation of the constraint can be canceled by
moving course $c_2$ to another period.
%% \item
%% In the context of graph coloring (see chapter~\ref{ch:GC}), if all courses 
%% last one time period, this is equivalent to changing the color of a vertex, 
%% for example. 
\end{enumerate}
\end{itemize}


To use simulated annealing effectively, it is crucial to use a good 
cooling schedule, and a good method for choosing new trial schedules, 
in order to efficiently sample the search space. We have experimented
with both these areas, which are discussed in the following sections.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{The Annealing Schedules}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The search for efficient cooling schedules may be thought of in two
quite different ways, which could be called ``off-line'' versus
``on-line'', or ``oracular'' versus ``constructive''. In all cases,
of course, what is desired is a cooling schedule which is as efficient
as possible: one which within a specified run time produces the
minimum possible expected energy, or one which produces a specified
expected energy in the minimum possible run time.

In an off-line or oracular construction, the computational cost of
finding the schedule itself is neglected. The extreme example of this
approach is the work of Strenski and Kirkpatrick~\cite{SK:annealing},
where (via a complete expansion of the transition matrices as a
function of temperature) the expected energy is computed as a function
of the temperature sequence $T_1$,$\ldots$,$T_t$ (one move at each
temperature).

The on-line or constructive approach is the one more commonly
taken. Here, all quantities of interest are estimated during an
annealing run, or with a limited amount of pre-computation. For
example, it is common to try to estimate the mean and variance of the 
energy in equilibrium at temperature $T$ from the energies of the
states produced during annealing at temperature $T$.

In practice it is the second approach (on-line) that is usually most
relevant. On the other hand, the first approach shows what is
possible, providing a yardstick against which on-line approaches may
be compared. Furthermore, on-line construction of a good cooling
schedule cannot proceed blindly. Off-line studies of optimal cooling 
schedules, whether theoretical or empirical, can indicate the
qualities that a good schedule should have; it is then the job of the
on-line algorithm to produce as good as possible an approximation to
this ideal.

In short, let $K$ be the number of ``generations'' used for annealing, 
$t_k$ is the number of attempted moves made (or ``time spent'') in
generation $k$, and $T_k$ is the temperature used in that
generation. Then the sequence $\{(T_k, t_k)\}^{K}_{k=1}$ and the value 
$K$, together, are referred to as the annealing schedule.


\subsection{Cooling Parameters}

The basic set of parameters required by each cooling is:
\begin{enumerate}
\item the temperature reduction factor,
\item time at each temperature,
\item initial temperature, and
\item final temperature,
\end{enumerate}

all of which must be determined.

%% In practice, among these parameters, the one with the least
%% significance is the temperature reduction factor. (The experiments
%% showed that if it is replaced with, say, its square root, and the
%% number of moves at the initial temperature is doubled, the outcome
%% will be the same.)

Since we desire a comparison of the cooling schedules over a range of
run times and solution qualities, the time spent at each temperature
was left as an independent variable.

The choice of initial temperature poses more of a challenge for
non-geometric schedules as we will see in the next sections. Also, for
non-geometric schedules the stopping criterion would vary, greatly
depending on the structure of the schedule. For geometric cooling, the
standard stopping criterion is to terminate the run when no moves had
been accepted for some predetermined number of generations. There are
also other stopping criteria.

\subsection{Dependence on Initial Temperature}

The initial temperature, $T_1$, is generally chosen to be high enough
that the canonical distribution is close to uniform. When the
Metropolis algorithm is used with a fixed distribution for candidate
moves, a related criterion is that the temperature be high enough that 
the rejection rate is very low. On other words, we only want to choose
an initial temperature to be high enough to give low rejection rates if
we are using a random start, whereas we certainly do not want low
rejection rates if we are using a preprocessor to start with a good
solution, otherwise our good solution will soon turn into a random
solution!

For $T_1$, we have simply used the upper bound suggested by
White~\cite{Wh:scale}, namely the standard deviation of the energies
of randomly-selected states in the space. As indicated by
Sorkin~\cite{Sorkin:landscape}, this can be found with sufficient
accuracy with a short run at infinite temperature. For further
discussion of the run lengths needed and the accuracy of the estimates
produced, the reader is referred to~\cite{OG:annealing}.

\subsubsection{Interpretation, and identification of a good initial
temperature}

Based on the work of White~\cite{Wh:scale},
Sorkin~\cite{Sorkin:landscape} adopted an approach in which 
``cooling curves'' (graphs of the mean energy at 
each temperature during annealing, plotted against the temperature) 
``quench out'' at different temperatures depending on how many moves
are made at each temperature. 

%% The interpretation of what is happening there is that, at high
%% temperatures both processes are in equilibrium, and so the energy
%% difference between two states at a given instant is a random number
%% symmetric about 0. The difference at a later instant is also random,
%% but if the times are close together, the two differences will be
%% correlated. The time over which the correlation is significant, the
%% ``correlation length'', is the mixing time. For two energies to be
%% effectively independent, they must be sampled at times more widely
%% separated than the correlation length.

In regard to the issues of the `set-up' runs to compute the initial
temperature: in order for an accurate  estimate it would be best to
make the set-up runs with as nearly as possible the same number of
moves $n$ (per temperature) as the actual runs planned. However, for
the sake of efficiency, pre-computations would need to be done over
shorter runs. This question of ``efficiency'' is significant here
since the pre-computation runs tend to range over the full temperature 
scale, while the final run will be made from a lower initial
temperature, and will be that much faster.

Overall, it is quite difficult to avoid having the pre-computations
consume more time than the run itself, and it is even worse with short 
run times. With a very long run, the optimal temperatures change very
slowly, so it is quite reasonable to determine the initial temperature 
for a very long run from a much shorter one; but that is not always true 
for short runs.

%% \subsubsection{Few Other Comments on the Initial Temperature}
%%
%% \begin{itemize}
%%
%% \item Based on the work of White~\cite{Wh:scale}, in principle we wish to
%% choose the initial temperature as the lowest temperature where the
%% energy timeseries is stationary. But the difficulty with this is that
%% the energy timeseries of a single annealing run cannot be stationary
%% if the temperature is being changed. 
%%
%% \item We note that if it is true that the relevant temperature range
%% increases with size of the problem instance, this would mean that the
%% effect of a high estimate for the initial temperature is relatively
%% small, thus, making the task that much easier.
%%
%% \item Although the simulated annealing literature already includes
%% numerous suggestions for stationarity tests (for example, see article
%% of Romeo and Sangiovanni-Vincentelli~\cite{RSV:framework}, book of van
%% Laarhoven and Aarts~\cite{vLA:anneal}, and book of Otten and van
%% Ginneken~\cite{OG:annealing}), none is appropriate for our problem
%% here. Since almost all of them claim ``convergence'' (in whatever
%% sense) down to the lowest temperatures used during annealing, this
%% would imply that it is below the temperature at which global
%% equilibrium is maintained (?!!) 
%%
%% \end{itemize}

\section{Our Experimental Schedules}
Three annealing schedules have been used in our experiments to update 
the temperature of the SA algorithm in Figure~\ref{SA-algorithm}:
geometric cooling, adaptive cooling, and adaptive reheating as a
function of cost.

The first schedule we have used is {\bf geometric cooling}, where the 
new temperature ($T^{'}$) of the SA algorithm is computed using
\be{geometric}
T^{'} = \alpha T\; ,
\ee
where $\alpha$ ($0 < \alpha < 1$) denotes the cooling factor. 
Typically the value of $\alpha$ is chosen in the range $0.90$ to
$0.99$. This cooling schedule has the advantage of being well
understood, having a solid theoretical foundation, and being the most
widely used annealing schedule. Its main disadvantage is that it
doesn't take account of the state of the system in any way, and thus
cannot adapt the intensity of the search depending on the difficulty
of the problem. 

Our results obtained from using this standard cooling schedule
will be used as a baseline for comparison with those using the other
two schedules, which allow the rate of cooling to be varied. 

It was noted by Sorkin~\cite{Sorkin:landscape} that if the same amount
of time was spent at each temperature, then subproblems of all energy
scales would be solved equally well. But subproblems of larger energy
scale will have a larger effect on the final cost. Therefore, it is
sensible to spend somewhat more time solving these large-energy
problems well, spending less time in the less important problems of
smaller energy scales. We can deduce from this statement that an
optimal cooling schedule should not be geometric, spending equal time
at each temperature, but adaptive in nature.

The second annealing schedule we used is the method of 
{\bf reheating as a function of cost} (RFC), 
which was used for timetabling by Abramson {\it et
al.}~\cite{ADK:cooling}, but the idea behind it is due to Kirkpatrick
{\it et al.}~\cite{KGV:sim, K:quant} and White~\cite{Wh:scale}.   
Before introducing this schedule we first summarize a few relevant
points on the concept of specific heat ($C_H$). Specific heat is a
measure of the variance of the energy (or cost) values of states at a
given temperature. The higher the variance, the longer it presumably
takes to reach equilibrium, and so the longer one should spend at the
temperature, or alternatively, the slower one should lower the
temperature.

Generally, in combinatorial optimization problems, phase 
transitions~\cite{HHW:phase, M:phase} can be observed as sub-parts of
the problem are resolved.
In some of the work dealing with the traveling salesman problem using 
annealing~\cite{L:fractal}, the authors often observe that the resolution of
the overall structure of the solution occurs at high temperatures, and
at low temperatures the fine details of the solution are resolved.
As reported in~\cite{ADK:cooling}, applying a reheating type procedure,
depending on the phase, would allow the algorithm to spend more time
in the low temperature phases, thus reducing the total amount of time
required to solve a given problem.

In order to calculate the temperature at which a phase transition
occurs, it is necessary to compute the specific heat of the
system. A phase transition occurs at a temperature $T(C_H^{max})$ when
the specific heat is maximal ($C_H^{max}$), and this triggers the
change in the state ordering. If the best solution found to date has a
high energy or cost then the super-structure may require
re-arrangement. This can be done by raising the temperature to a level
which is higher than the phase transition temperature $T(C_H^{max})$.
Generally, the higher the current best cost, the higher the
temperature which is required to escape the local minimum.
To compute the aforementioned maximum specific heat, we employ the
following steps~\cite{ADK:cooling, vLA:anneal, OG:annealing}.

At each temperature $T$, the annealing algorithm generates a set of 
configurations ${\cal C}(T)$. Let $C_i$ denote the cost of 
configuration $i$, $C(T)$ is the average cost at temperature $T$, 
and $\sigma(T)$ is the standard deviation of the cost at $T$. 

At temperature $T$, the probability distribution for configurations is: 
\be{T.1}
P_i(T) = {e^{-C_i \over {kT}} \over \sum_j{e^{-C_j \over {kT}}}} \,.
\ee

The average cost is computed as: 
\be{T.2}
<C(T)> = \sum_{i \in {\cal C}} C_i P_i(T) \,.
\ee

Therefore, the average square cost is: 
\be{T.3}
<C^{2}(T)> = \sum_{i \in {\cal C}} C_i^2 P_i(T) \,.
\ee

The variance of the cost is:
\be{T.4}
\sigma^{2}(T) = <C^{2}(T)> - <C(T)>^{2} \,.
\ee

Now, the specific heat is defined as: 
\be{T.5}
C_H(T) = {{\sigma^{2}(T)} \over {T^{2}}} \,.
\ee
The temperature $T(C_H^{max})$ at which the maximum specific heat
occurs, or at which the system undergoes a phase transition, can thus
be found. 

Reheating sets the new temperature to be
\be{T.6}
T = K \cdot C_b + T(C_H^{max}) \, ,
\ee
where $K$ is a tunable parameter and $C_b$ is the current best cost. 
Reheating is done when the temperature  drops below the phase
transition (the point of maximum specific heat) 
and there has been no decrease in cost for a specified number of
iterations, i.e. the system gets stuck in a local minimum. Reheating
increases the temperature above the phase transition (see equation~\ref{T.6}), 
in order to produce enough of a change in the
configuration to allow it to explore other minima when the temperature
is reduced again.

The third cooling schedule we have tried is {\bf adaptive cooling}. 
In this case, a new temperature is computed based on the specific
heat, i.e. the standard deviation of all costs obtained at the current
$T$. The idea here is to keep the system close to equilibrium, by
cooling slower close to the phase transition, where the specific heat
is large. There are many different ways of implementing this idea,
we have chosen the approach taken by Huang {\it et
al.}~\cite{HRS:efficient}, which was shown to yield an efficient
cooling schedule. Let $T_j$ denote the current temperature, at step
$j$ of the annealing schedule. After calculating $\sigma(T_j)$ from
equation \ref{T.4}, the new temperature $T_{j+1}$ is computed as
follows:

\be{T.6.1}
T_{j+1} = T_j \cdot e^{-{{a T_j} \over \bar\sigma(T_j)}} \, ,
\ee

where $a$ is a tunable parameter.
Following suggestions by Otten and van Ginneken~\cite{OG:annealing}
and Diekmann {\it et al.}~\cite{DLS:distributed}, $\sigma(T_j)$ is
smoothed out in order to avoid any dependencies of the temperature
decrement on large changes in the standard deviation $\sigma$. 
We used the following standard method to provide a smoothed standard
deviation $\bar\sigma$:

\be{T.6.0}
\bar\sigma(T_{j+1}) = (1 - \omega) \sigma(T_{j+1}) + \omega \sigma(T_j)
{T_{j+1} \over T_j}
\ee

and set $\omega$ to $0.95$.
This smoothing function is used because it follows (from the form of
the Boltzmann distribution, see~\cite{Sorkin:landscape, Wh:scale})
that it preserves the key relationship:
\be{T.7}
{d \over {dT}}C(T) = {\bar\sigma^{2}(T) \over T^{2}} = C_H
\ee

Note that reheating can be used in conjunction with any cooling
schedule. We have used it with adaptive cooling.

\subsection{On Studies of Cooling Schedules}
It is unfortunate that in most published studies of automatic cooling
schedules we are aware of, the actual form of the schedule is not
presented. That is, while an algorithm for setting the temperatures
and times is given, the resulting empirical dependence of temperature
on time is not, and it is impossible to reconstruct this without
re-implementing the algorithm. There is much to be learned about these 
schedules. For start, answers to these questions could be used
to improve the scheduling algorithm: Is the schedule geometric? Could
it start at a lower temperature? Should it stop at a higher
temperature? How does the form change as parameters are altered?

Studies of cooling schedules also generally compare two or three
algorithms -- or rather, algorithm implementations -- and it is hard
to judge any of them in an absolute sense. Since geometric schedules
were the {\it de facto} standard before good adaptive schedules were
developed; since a number of studies, and reasoning going back to the
work of White~\cite{Wh:scale}, gives intuitive support for their
efficiency, and since some adaptive schedules may wind up being
roughly geometric anyway; comparisons with ``ideal'' geometric
schedules, found by exhaustive search or any other means, would be
something quite useful to have. 

\section{Annealing in Comparison with Other Algorithms}

While it is valuable to know about the time versus quality tradeoff of 
annealing itself, knowing whether annealing is a good algorithm to use 
depends on its relative performance compared with other
algorithms. The most significant set of experimental comparisons is
that of~\cite{JAMS:optimization-I, JAMS:optimization-II,
JM:traveling}, where for each of a number of problems annealing was
compared against the best algorithm known for that problem. Generally
speaking the specialized algorithms tailored to the particular problem 
at hand ran far faster than annealing, but annealing's results were
often as good or better. Overall, if annealing was not always the best 
algorithm, it was at least competitive.
