\chapter{Simulated Annealing}

Simulated annealing (SA) has been
widely used for tackling different combinatorial optimization
problems, particularly academic scheduling~\cite{39, 3, 4}. 
The basic algorithm is described in Figure~\ref{SA-algorithm}.
The results obtained depend heavily on the cooling schedule used. We
initially 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{2}.

A comprehensive discussion of the theoretical and practical details of
SA is given in~\cite{57, 60, 80, 17}. 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.

Depending on the problem to which it is applied, 
SA appears competitive with many of the best 
heuristics, as shown in the work of Johnson {\it et al.}~\cite{51}. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\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 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 cost of an imbalanced class assignments (a certain
period will have more classes assigned to 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, 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).
\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.
\item
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 of efficient cooling schedules may be thought of in two
quite different ways, which could be called ``off-line'' versus
``on-line'', or ``oracular'' versues ``constructive''. In all cases,
of course, what is desired is a colling 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 extereme example of this
approach is the work of Strenski and Kirkpatrick [see Sorkin -ref 62],
where (via a complete expansion of the transition matrices as a
dunction 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 precomputation. 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 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.

\subsection{Cooling Parameters}

Each cooling takes a set of parameters:
\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 controlled for.

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
depended 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}
$\ldots$

\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. 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.

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{2}, 
but the ideas behind it are due to Kirkpatrick {\it et al.}~\cite{52, 53}
and White~\cite{41}. 
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 cost (or energy) 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{66, 8} 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{7}, 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{2}, 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{2, 17, 60}.

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{59}, 
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{60} and Diekmann
{\it et al.}~\cite{61}, $\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{80, 41}) 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.


