\chapter{Optimization Networks}\label{ch:ON}

\section{Introduction}

A typical human brain is superior to a digital computer in performing 
many tasks, e.g. learning, decision making, verbal behavior, visual recognition
and motor functions. It is also small, robust and energy efficient. Artificial
neural networks (ANNs) stem from simple models of how the neural networks in the 
brain is thought to work. The main reason for most people (including myself) for 
using and developing ANNs is not to build an artificial brain but to formulate 
models and algorithms with some attractive features of the brain. 

The brain's impressive performance still activates a lot of research on 
general properties of neural networks and will probably continue
to do so for a long time. The very simple model of the brain making the 
foundation of ANNs is based on the traditional viewpoint of the brain as
a computer with neurons working as calculators, adding and subtracting signal
spikes until a threshold for firing is reached. In 1943 McCulloch and Pitts made
a mathematical model of a neuron, the so called McCulloch-Pitts 
neuron~\cite{MP:logical}, see equation~\ref{MP-neuron}.

\be{MP-neuron}
s_i \rightarrow g(\sum_j w_{ij} s_j - \theta_i)
\ee

With the transfer function, $g(x)$, chosen to be the Heaviside step function 
($\Theta(x)$ = 1 if $x > 0$ otherwise 0), as McCulloch and Pitts did, a neuron
either has the value 0 or 1. Other often used transfer functions are 
$g(x) = sign(x)$ and
$g(x) = tanh(x)$. Even though this neuron is simple, a connected network
of such neurons is capable of performing any computations possible on an ordinary
digital computer~\cite{MP:logical}. This artificial neuron is the basic building 
block of most types of feed-forward networks, which are used mainly for 
classification and function approximation tasks.

A multitude of different directions within the artificial neural network field 
have emerged. One particular direction from the basis of model development 
is the Hopfield feed-back model, which will be disscused in the 
section~\ref{sec:HN}.

\section{Neural Methods as Optimization Paradigms}

While there already exist many effective algorithms for approximately
solving combinatorial optimization problems, this chapter is concerned 
with just one family of techniques, which we term {\it optimization
networks}. The origins of optimization networks can be traced back to
Hopfield and Tank's pioneering work in 1985~\cite{HT:neural}, though
it would perhaps be more appropriate to start with a review of general 
Artificial Neural Networks.

It has long been noted that while conventional computers are very good 
at performing numerically intensive tasks such as complex arithmetic,
they are less well adapted to the tasks which humans find
straightforward, like speech and vision recognition. For this reason
a new type of computer, modeled loosely on the architecture of the
human brain, was proposed. These Artificial Neural Networks (ANNs)
comprise a large number of simple processing elements (corresponding
to single neurons in the human brain) connected together in massive,
parallel arrays. An ANN can be trained to perform tasks such as speech 
and vision recognition and , like the human brain, has the ability to
learn from experience. Moreover, with the highly parallel structure of 
the ANN properly exploited in electronic hardware, extremely high
information processing speeds are possible, giving ANNs a huge
advantage over conventional computers.

\section{The Hopfield Network}\label{sec:HN}

One type of ANN is the Hopfield network~\cite{Hop:neural}, originally
proposed as a form of content addressable memory (a device which
allows stored patterns to be recalled by presentation of noisy,
corrupted versions of the same patterns). The Hopfield network is an
example of a feedback neural network, where the outputs of the
individual processing units are fed back to the inputs via a dense
array of interconnections, producing a nonlinear, continuous dynamic
system. 

Moreover, in~\cite{Hop:neural}, Hopfield introduced the idea of an energy
function into neural network theory (for the associative memory
problem). The term energy function originates from an analogy with
magnetic spin systems, in other contexts names like Lyapunov function, 
Hamiltonian, fitness function or objective function would have been
used.

The network described by Hopfield in his seminal 1982
paper~\cite{Hop:neural} is in fact a special case of the {\it additive
model} developed by Grossberg in the
1960's~\cite{Gro:neural}. Hopfield reviewed the network's application
as a content addressable memory with both binary~\cite{Hop:neural} and
continuous valued~\cite{Hop:network} outputs. However, 
after detailed investigations into the network's performance by a
number of authors, it is now clear that other types of content
addressable memory are far more efficient~\cite{TTR:neural}.

The Hopfield network was first proposed for combinatorial optimization 
applications in~\cite{HT:neural, HT:circuit}, where a penalty function 
mapping for the traveling salesman problem was given, more or less 
marking the start of using artificial neural networks techniques for 
optimization problems. In~\cite{HT:neural} the TSP was approached 
by mapping it onto an energy function of the type,

\be{TSP-energy-func}
E = cost + constraints.
\ee

This energy is expressed in terms of neurons in such a way that the
minimum corresponds to a neuron configuration representing a solution
to the problem. Hopfield and Tank's mapping was of a rather {\it ad
hoc} nature, with separately weighted penalty functions for each of
the hard constraints. The weights were set by trial and error,
resulting, not surprisingly, in poor performance~\cite{WP:hop}. 
Nevertheless, some researchers preserved with this
approach, sometimes proposing modified networks to correct what is
essentially a sloppy mapping.

Another optimization problem that can be mapped into such an energy
function as that of equation~\ref{TSP-energy-func} is the graph
bisection problem, see Figure~\ref{gb-problem}. The task is to divide
a given graph (nodes and links) into two equal parts with a minimal
cut size (number of links) between them. For each node $i$ a binary
neuron, $s_i$, is assigned $+1$ or $-1$ depending on whether the
neuron is in the left box or the right box, respectively. For each
pair of nodes, let $w_{ij} = 1$ if a link between them exists and
otherwise $0$. A suitable energy function is

\begin{figure}[tbp]
\begin{center}
\begin{picture}(195,75)
\psfig{figure=gb-problem.ps,height=195pt,width=75pt,angle=-90}
\end{picture}
\end{center}
\caption{A graph bisection problem.}
\label{gb-problem}
\end{figure}


\be{gb-energy-function}
E = - {1 \over 2}\sum_{ij} w_{ij} s_i s_j + {\alpha \over 2} {\left(
\sum_i s_i \right)}^2\;.
\ee

The first term expresses the cut size cost whereas the second term is
a constraint term penalizing situations where the nodes are not
equally partitioned, and $\alpha$ sets the relative strength between
the two terms. The cost and the constraint have to ``compete'',
i.e. the problem is frustrated, this typically leads to the appearance 
of many local minima. Once the energy is written down, the next step
would be to find the minimum. This could be done using methods such as 
mean field annealing, discussed in the next section.

Real progress in the field of optimization networks can be identified
with three main movements: rigorous investigation of network
convergence, the development of more sensible mapping techniques, and
the emergence of powerful annealing procedures. 
\begin{itemize}
\item {\bf Convergence}: Generalized feedback neural networks can
exhibit oscillatory \cite{CG:stability} or even
chaotic~\cite{LG:dynamics} behavior. It is important to determine
under what conditions a network will admit a Liapunov function for
optimization applications. For simulated networks, convergence
properties depend on whether the nodes are updated in serial or
parallel mode: such issues are discussed (for discrete-valued
networks) in~\cite{Br:convergence}.

\item {\bf Mapping}: The development of rigorous mappings was delayed
by the popularity of the original, {\it ad hoc} mapping of the TSP
problem. Researchers concentrated on finding better weights for the
penalty functions, often by way of an eigenvector analysis of
Liapunov function~\cite{ANF:theory}: this, though correct, is quite
complicated. Simpler, more reliable mappings began to appear for
specific problems in~\cite{AF:network}; some also dealt with inequality
constraints~\cite{TCP:neural}. One particular rigorous mapping for the
most general 0-1 quadratic programming problems onto standard Hopfield 
network dynamics was presented in~\cite{GP:function}

\item {\bf Annealing}: In parallel with the development of rigorous
mappings, modified dynamics were proposed as a way of improving the
quality of any valid solutions found by the network. The most
effective of these modifications typically represented some sort of
annealing procedure. Mean field annealing was developed as an
approximation to simulated annealing, and will be discussed from that
viewpoint in the next section. More direct annealing procedures
include hysteretic annealing~\cite{EDKT:neural}, convex
relaxation~\cite{OB:neural}, and matrix graduated
non-convexity~\cite{AF:network}.
\end{itemize}

\subsection{Structure and Dynamics of the Net}

\begin{itemize}
\item The network is built by connecting a large number of neurons to
each other to form a dense, parallel array. In general, any $i^{th}$
neuron in the net is described by two variables: its current state
(input) $u_i$,  and its output $v_i$.

The output is related to the input by a monotonically increasing
transfer function: $v_i = g_{\lambda}(u_i)$.

\item The net's dynamics are governed by the first order differential
equation:
\be{FODE.1}
{d{\bf u} \over {dt}} = 
- \eta {\bf u} + {\bf W v} + {\bf I}^b\;.
\ee

\item $g_{\lambda}(u_i) = {1 \over {1 + e^{-2\lambda u_i}}}\;$, for
large $\lambda$, this function approximates a step function,
forcing $v_i$ to $[0,1]\;\Longrightarrow\;$ (the output of the net
lies on the unit hypercube).

\item $W_{ij}$ is the connection strength between neuron $i$ and neuron
$j$. 

\item $I_i^b$ is an input bias to neuron $i$.

\item Neurons may be updated individually at random (asynchronously)
or all together (synchronously).
\end{itemize}

\subsection{Operation of the Network}

\vspace{0.1cm}
\begin{figure}[tbp]
\begin{center}
\begin{picture}(300,300)
\psfig{figure=hopfield_net.eps,height=300pt,width=300pt}
\end{picture}
\end{center}
\caption{A schematic view of the Hopfield network.}
\label{Hop-model}
\end{figure}

The network (see Figure~\ref{Hop-model}) operates as follows:

\begin{itemize}
\item The network's dynamics are governed by a set of nonlinear
coupled differential equations of the general form shown in 
equation~\ref{FODE.1},
where {\bf u}, {\bf W}, {\bf v}, and {\bf I} are vectors.

\item The $I_{i}^b$ corresponds to an input bias to neuron $i$.
\item The $\eta$ term introduces an element of decay into the
network's dynamics, which can be used to ensure that ${\bf u}$
variables remain bounded in magnitude.
\end{itemize}

\subsection{Updating of the States of the Neurons}

To make sure that stability of the of the net is guaranteed,
$u_i$ is updated continuously by evaluating ${d \over dt}u_i$ as a
function of the total input to the $i^{th}$ neuron. See equation~\ref{FODE.1}, 
where 

\be{State.2}
v = g(u)\;.
\ee

Equations \ref{FODE.1} and \ref{State.2} govern the dynamics of the
network.

\subsection{Stability of the Net}

For symmetric {\bf W}, Hopfield proved stability for the
network by showing that with the above dynamic equations
(\ref{FODE.1} and \ref{State.2}) the network output {\bf v} evolves
so as to minimize the {\it Liapunov function}:

\be{Liapunov}
E^{liap}(v) = 
\left( -{1 \over 2} \sum_{ij} W_{ij} v_i v_j -\sum_i I_{i}^b v_i \right) 
+ \eta{\displaystyle \sum_i}\, {\displaystyle \int_{0}^{v_i}}\, g^{-1}(V)\,dV\;.
\ee

\begin{itemize}
\item $v = (v_1 \ldots v_n)$ is denoting the net state, therefore we
have an energy function defined throughout $v \in [0,1]^n$ and
dynamics which guarantee convergence to minima.

\item For a carefully chosen set of connection weights $W$, the minima of 
the first term lie at the hypercube vertices. The second term is minimized at 
the hypercube center but has negligible impact at low $\eta$.
\end{itemize}

\subsection{Updating the States of the Net}
Updating of the states (or the output of the neurons) is carried
out either asynchronously or synchronously. From a simulation point of
view, we found that synchronous updating a bit more efficient, but
from a theoretical point of view, it does not make much of a 
difference on the operation of the network as to which updating mode
is used.

\subsection{Mapping the Problem onto the Network}
A representation of the problem is needed in which the feasible
solutions lie at hypercube vertices. Each hypercube vertex is regarded
as a configuration of the problem, with an associated energy (or cost)
function. Next, the partition function is then formed by summing the
Boltzmann factors of the admissible configurations. A quadratic energy
function (to be mapped onto the first term of \ref{Liapunov}) for
which the minima correspond to solutions and the depth of each minimum
reflects the solution quality. Aside from its intellectual interest,
this approach also offers potentially great benefits, for example,
equation \ref{State.1} can also describe the behavior of a network of
interconnected electrical amplifiers, and hence hardware circuits
might be capable of solving hard optimization problems in real
time. In practice however, even at a software level implementation,
considerable difficulties have been encountered carrying this program
out for the TSP problem, which was used by Hopfield and
Tank~\cite{HT:neural} to test their neural model.

In summary, here are the necessary four steps to be taken to map a
combinatorial optimization problem such as the TSP or scheduling onto
the Hopfield-Tank model:
\begin{enumerate}
\item Choose a representation scheme which allows the activation
levels of the units to be decoded into a solution of the problem.
\item Design an energy function whose minimum corresponds to the best
solution of the problem.
\item Derive the connectivity of the network from the energy function.
\item Set up the initial activation levels of the units.
\end{enumerate}

These ideas can easily be applied to the design of a Hopfield-tank
network in a course scheduling context:

\begin{itemize}
\item First, a suitable representation of the problem must be
chosen. For example, one representation is an $N$ x $M$ matrix of
units, where each row correspond to a particular class and each
column to a particular time slot that the class is given at. Another 
representation would be to have each row to correspond to a particular
teacher and each column to a particular class taught by that teacher
at a particular time slot. Therefore, on the network itself if the
activation level of a given unit, for example $V_{nm}$, is close to 1, 
it is then assumed that class $n$ is assigned to time slot $m$ in the 
scheduling table. In this way, the final configuration of the network
can be interpreted as a ``solution'' to the problem. Since we are
dealing with a large scale scheduling problem, a very large number of
units are needed to encode into the network any worthwhile
``solution'' to our problem. 

\item Second, the energy function must be defined. It can be expressed
as a summation of various terms, and each term represents a constraint
or a set of constraints. The cost of each term is controlled by a
coefficient (i.e. A, B, C, etc.) multiplied by each term. The terms
collectively define a feasible schedule and each term's contribution
to the final cost is specified by its corresponding coefficient.

\item Third, the bias and connection weights are derived. To do so,
the energy function derived in the previous step is compared to the
generic energy function defined by the problem from the given set of
hard constraints. It is quite essential that any worthwhile solution
the network gives would need to satisfy this class of constraints. 

\item The last step is to set the initial activation value of each unit
of the network to a small constant, perhaps proportional to the size
of the net or proportional to the number of courses, plus or minus a
small random perturbation (this way the sum of the initial activations
is approximately equal to the number of courses $N$).

\end{itemize}


\subsection{Shortcomings of the Hopfield Model}

The main shortcomings of the original Hopfield-Tank
model~\cite{HT:neural}, as pointed out by~\cite{WP:hop} are as follow: 

\begin{itemize}
\item In the case of dealing with an N-city TSP, the network would
need $O(N^2)$ units and $O(N^4)$ connections. Thus, for a reasonably
large $N$ we would end up with a very large number of connections.

\item The net's partition function sums over a vast number of
configurations (e.g. for the TSP, it is not solved in space
$O(N!)$, but in a space ${2^{N}}^2$). 
These configurations, again, for the TSP, are nothing like valid
tours, and even though these offenders have small Boltzmann weights
their large number inevitably affects the thermal average quantities. 

\item The model performs a gradient descent of the energy function in
the configuration space, and is thus plagued with the limitations of
``hill-climbing'' approaches, where a local optimum is found. As a
consequence, the performance of the model is very sensitive to the
initial starting configuration.

\item The model does not guarantee feasibility. In other words, many
local minima of the energy function correspond to infeasible
solutions. This is related to the fact that the constraints of the
problem, in case of TSP namely that each city must be visited exactly
once, are not strictly enforced but rather introduced into the energy
function as penalty terms.

\item Setting the values of the coefficients (i.e. A, B, C, etc.) in
the energy function is much more of an art than a science and requires 
a long ``trial-and-error'' process. Setting the coefficients to small
values usually leads to short but infeasible tours in the case of the
TSP. Alternatively, setting the penalty parameters or coefficients to
large values forces the network to converge to any feasible solution
regardless of the total length. Moreover, again in the TSP, it seems
to be increasingly difficult to find ``good'' coefficient settings as
the number of cities increases.

\item It usually takes a large number of iterations (in the thousands) 
before the network converges to a solution. Moreover, the network can
``freeze'' far from a corner of the hypercube in the configuration
space, where it is not possible to interpret the configuration as a
solution to the problem. This phenomenon can be explained by the shape 
of the sigmoidal activation function which is very flat for large
positive and large negative input. Consequently, if the activation
level of a given unit is close to zero or one, even large
modifications to the input will produce only slight modifications to
the activation level. If a large number of units are in this
state, the network will evolve very slowly, a phenomenon referred
to as ``network paralysis''. Paralysis far from a corner of the
hypercube can occur if the slope of the activation function is not
very steep. In that case, the flat regions of the sigmoidal function
extend further and affect a larger number of units (even those with
activation levels far from zero and one).

\item The network is not adaptive, because the weights of the network
are fixed and derived from the problem's data, rather than taught from 
it.
\end{itemize}

\section{Mean Field Annealing for Timetabling}

\subsection{Introduction and Motivation}

One of the potential drawbacks of using simulated annealing for hard 
optimization problems is that finding a good solution can often take an 
unacceptably long time. 
Mean-field annealing (MFA) attempts to avoid this problem by using a
deterministic approximation to simulated annealing, by attempting to
average over the statistics of the annealing process. Essentially, the
{\it mean field approximation} involves averaging statistics of two
coupled  spins, and replacing a complex function with its truncated
Taylor expansion around a {\it saddle point}. It is these
approximations which differentiate mean field annealing from simulated
annealing. 

The result is improved execution speed at the expense of solution
quality. Although not strictly a continuous descent technique, MFA is
closely related to the Hopfield neural network, on which we would like
to state few more points. We clearly know, and outlined in the
previous section, that the original Hopfield 
algorithm had difficulty constraining the network into valid
solutions. Therefore, one possible way forward is to reduce the
computational burden being placed on the net. For example, in the case
of the traveling salesman problem, this can be done by constraining
each city to be ``on'' only once, by enforcing $\sum_a V_{ia} = 1\;
\forall i$, ($V_{ia} = 1/0$, which means that the city $i$ is (is not)
the $a^{th}$ city visited in the  tour), rather than relying on an
energy penalty term to try to explicitly do this, as is the case in
the Hopfield model.  

This idea was first utilized and analyzed for the TSP 
case by, among others, Peterson and
S$\ddot{o}$derberg~\cite{PS:mapping}. It greatly enhances the degree 
of convergence to valid solutions. Later, Peterson and colleagues
somewhat extended the technique and applied it to a case of high
school timetabling~\cite{GSP:complex}. 
%% add to the above paragraph: Other similar work that roughly 
%% appeared around the same time as~\cite{PS:mapping} is that
%% of Simic, Fox et al (give citations) on deterministic annealing.  

In~\cite{PS:mapping} the normalization (see definition~\ref{norm-dfn})
occurs in the context of a {\it mean field annealing} approach. This
generally yields the same final network equations as the ``neuronal
circuit'' approach of Hopfield, but the exposition is a little bit
clearer, as it is laid out  in a statistical mechanics
framework\footnote{This equivalence was also noted in
Hopfield and Tank's 1985 paper~\cite{HT:neural}. Martin Simmen (in a
personal communication) brought this point to our attention.}. The
idea is to regard each hypercube vertex as a configuration, with an
associated energy (or cost) ${\it E}$. The partition function is then
formed by summing the Boltzmann factors of the  admissible
configurations. For example, in the Hopfield algorithm, for the TSP,
all of the $2^{N^2}$ configurations are admissible, whereas if each
city is restricted to being visited only once, then only $N^N$
vertices are admissible. After taking a mean field approximation, {\it
saddle point equations} are derived (see section~\ref{sec:MFA}), the solutions
of which pick out the dominant states of the network at the current
temperature $T$.

This statistical mechanics characterization gives a somewhat clearer
picture of why the normalized model (which corresponds to a type of
Potts model~\cite{Wu:potts} in physics, in that, again for the TSP,
each city is represented by a spin which can be in only one of $N$
states) ought to perform better than Hopfield model. The other
advantage that draws us to this derivation is that, by introducing a
temperature, the concept of annealing, i.e. reducing $T$ during a run,
can be justifiably employed. 

\subsection{Potts Neurons}

It is advantageous to use a multi-valued neurons (or Potts neurons)
instead of binary neurons for many types of optimization problems,
e. g. course scheduling, frequency assignment, airline crew scheduling, 
and routing problems. The use of Potts neurons when encoding a problem 
can reduce the number of terms needed in the energy function
(equation~\ref{TSP-energy-func} of Hopfield type). Fewer energy terms
in general result in a less complex energy landscape and hence the
global minimum would be much easier to locate.

\subsection{Mean-Field Approximation}\label{sec:MFA}

Given a problem mapped onto an energy function $E(s)$ (of Hopfield
type) the next step is to locate a global minimum. A gradient descent
method would lead to the nearest minimum which probably is not global
minimum. Alternatively, instead of using a stochastic method like
simulated annealing, SA equations could be replaced with deterministic
equations. The basic idea is that it is possible to approximate the
actual cost or energy function $E$, which is a function of discrete
neural variables $S_{ia}$, by an effective energy function $E^{'}$
that can be represented in terms of continuous variables $U_{ia}$ and
$V_{ia}$. These are known as mean field variables, since $V_{ia}$ is
an approximation to the average value of $S_{ia}$  at a given
temperature $T$.

This approach effectively smoothes out the energy function and makes it
easier to find the minimum value, which is obtained by solving 
the {\it saddle point equations}:
${\partial E^{'} \over \partial V_{ia}} = 0$ and  
${\partial E^{'} \over \partial U_{ia}} = 0$, 
which generate a set of self-consistent mean field theory (MFT)
equations in terms of the mean field variables $U$ and $V$:

\be{spe-1}
U_{ia} = - {1 \over T} {\partial E \over \partial V_{ia}}
\ee

\be{spe-2}
V_{ia} = {e^{U_{ia}} \over \sum_b e^{U_{ib}}} \;\;\; (\equiv f_{ia})
\ee

The MFA algorithm involves solving equations \ref{spe-1} and
\ref{spe-2} at a series of progressively lower temperatures $T$: this
process is known as temperature annealing. The critical temperature
$T_c$, which sets the scale of $T$, is estimated by expanding equation
\ref{spe-2} around the trivial fixed-point~\cite{PS:mapping,
GSP:complex} $V^{(0)}_{ia} = {1 \over N_a}$, where $N_a$ is the number
of possible states of each of the network neurons. For example, for
the events defined by professor-class pairs $(p,q)$ mapped onto
classroom-timeslots $(x,t)$, we have $N_p N_q$ neurons, each of which
has $N_x N_t$ possible states, in which case $V^{(0)}_{pq;xt} = {1
\over N_x N_t}$. 

Equations \ref{spe-1} and \ref{spe-2} can be solved iteratively using
either synchronous or serial updating. The iterative dynamics to
evolve the mean field variables toward a self-consistent solution 
is explained in detail by Peterson et al.~\cite{GSP:complex}
The solutions correspond to stable states of the Hopfield
model~\cite{HT:neural}.

Observe from equation \ref{spe-2} that any solution to the MF
equations respects a continuous version of the Potts condition 

\be{potts-cond}
\sum_a V_{ia} = 1  \;\;\;\forall\;i .
\ee

\bd{norm-dfn}
\ \ In~\cite{PS:mapping}, it was stated that the performance of MFA on
certain problems can be improved using {\it neuron
normalization}. This idea was originally inspired by an analogy with
Potts glasses in statistical physics~\cite{KS:graph}. {\it Neuron
normalization} modifies the MFA algorithm so that constraints of the
form  $\sum_{a \in \text{states}} \nu_{ia} = 1$ (to be consistent with 
our notation we instead use the form $\sum_a S_{ia} = 1$) are
explicitly enforced, without the need for a penalty term. It was found
by Peterson and other to be particularly useful for mapping the TSP
and graph partitioning (GP) type problems; though the technique does
not have general applicability.
\ed

\subsection{Mapping of Optimization Problems onto the Potts Model}

Th Mean-field annealing algorithm has been successfully applied to, as
mentioned above, high school class scheduling. For scheduling, it is
advantageous to use a  Potts neural encoding to specify discrete
neural variables (or neurons) for the problem. This is defined in its
simplest form as a mapping of events onto space-time slots, for
example an event $i$, in this case a professor-class pair ($p,q$), is
mapped onto a space-time slot $a$, in this case a classroom-timeslot
pair ($x,t$).  Now, the Potts neurons $S_{ia}$ are defined to be $1$
if event $i$ takes place in space-time slot $a$, and $0$ otherwise.
In this way, the constraints involved can be embedded in the neural
net in  terms of the weights $w_{i,j}$ of the network, which encode a 
Potts normalization condition such as $\sum_{a} S_{ia} = 1$.

Another mapping example is the general graph partition problem. Here
the task is to evenly partition $N$ nodes into $K$ boxes with a
minimal cut-size between them (see Figure~\ref{gp-problem}). If binary 
neurons were to be used (as was suitable for the graph-bisection
problem), then a neuron, $S_{ia}$, could encode whether node $i$ is in 
box $a$ or not. Since a neuron only can be in one box (in the final
result), an energy term enforcing

\be{gp-energy}
\sum_b S_{ib} = 1
\ee

must be used. To reduce the available solution space an encoding using 
the Potts neurons can be used. A Potts neuron has $K$ possible values
or states, $\vec{S} = (S_{i1}, S_{i2}, \ldots, S_{iK})$ where the
$j$th possible state is given by $S_{ij} = 1$ and $S_{ik} = 0$ for $k
\neq j$. These states are then mutually orthogonal and automatically
satisfy equation~\ref{gp-energy}. A Potts neuron could then encode to
which box a neuron is assigned, e. g. $\vec{S} = \{0\; 0\; 1\; 0\}$ means
that node $i$ is in box 3.

\begin{figure}[tbp]
\begin{center}
\begin{picture}(200,175)
\psfig{figure=gp-problem.ps,height=200pt,width=175pt,angle=-90}
\end{picture}
\end{center}
\caption{A graph partition problem.}
\label{gp-problem}
\end{figure}

\subsection{The Mean-Field Annealing Algorithm}\label{sec:MFAalg}

The generic MFA algorithm appears in Figure~\ref{MFA-algorithm}. 
At high temperatures $T$, the mean-field solutions will be states 
near the fixed-point symmetrical maximum entropy state $V_{ia} =
1/N_a$. At low temperatures, finding a mean-field solution will be
equivalent to using the Hopfield model, which is highly sensitive to
the initial conditions and known to be ineffective for hard
problems. MFA improves over the Hopfield model by using
annealing to slowly decrease the temperature in order to sidestep
these problems.
 
These characteristics are similar to those of simulated annealing,
which is no surprise since both it and the mean-field method compute
thermal averages over Gibbs distributions of discrete states, the
former stochastically and the latter through a deterministic
approximation. It is therefore natural to couple the mean-field method
with the concept of annealing from high to low temperatures. 

In addition to the structure of the energy function, there are three
major interdependent issues which arise in completely specifying 
a mean-field annealing algorithm for a timetabling problem:

\begin{itemize}
\item The values of the coefficients of terms in the energy function.
\item The types of dynamics used to find solutions of the MFT
 equations at each $T$.
\item The annealing schedule details, i.e. the initial temperature
 $T(0)$, the rules for deciding when to reduce $T$ and
 by how much, and the termination criteria.
\end{itemize}

Peterson {\it et al.}~\cite{GSP:complex} introduced a quantity 
called {\it saturation}, $\Sigma$, defined as
\be{saturation}
\Sigma = {1 \over N_i} \sum_{ia} V^{2}_{ia} \;\; ,
\ee

where $N_{i}$ is the number of events (in this case the number of 
professor-class pairs).
This characterizes the degree of clustering\footnote{Other measures could
also do this, e.g. an entropy measure $-\sum_{ia}\nu_{ia}ln(\nu_{ia})$.}
 of events in time and/or space, $\Sigma_{min} = {1 \over N_a}$
corresponds to high temperature, whereas $\Sigma_{max} = 1$ means that
all the $V_{ia}$ have converged to 0 or 1 values, indicating that each
event has been assigned to a space-time slot.


\begin{figure}[tbp]
% \vspace{0.25cm}
% \hrule
\medskip
\begin{center}
\begin{enumerate}
\item Choose a problem and encode the constraints into weights
$\{w_{ij}\}$. That is, map the problem onto a Hopfield type of energy
function using Potts neurons.
\item Find the approximate phase transition temperature by linearizing
  equation (\ref{spe-2}).
\item Add a self-coupling $\beta$-term if necessary. In a neural net, 
  this corresponds to a feedback connection from a neuron to itself.
\item Initialize the neurons $V_{ia}$ to high temperature values
  ${1 \over N_a}$ plus a small random term such as ${rand[-1,1] \times 0.001}$; 
  and set $T(0) = T_c$.
\item Until ($\Sigma \geq 0.99$) do:
\begin{itemize}
  \item At each $T(n)$, update all $U_{ia}$ and $V_{ia}$ by iterating to
     a solution of the mean field equations.
  \item $T(n+1) = \alpha T(n)$, we choose $\alpha = 0.9$
\end{itemize}
\item The discrete values $S_{ia}$ that specify the schedule are obtained 
   by rounding the mean field values $V_{ia}$ to the nearest integer 
   (0 or 1). 
\item Perform {\it greedy heuristics} if needed to account for possible 
   imbalances or rule violations.
\end{enumerate}
\end{center}
% \smallskip
% \hrule
% \smallskip 
\caption{The Generic Mean-Field  Annealing Algorithm}
\label{MFA-algorithm}
% \vspace{0.25cm}
\end{figure}

\subsection{Observations on the Algorithm}

The first step of Figure~\ref{MFA-algorithm} is to map the constraints
of the problem into the neural net connection weights.
In our implementation, at each $T_n$ the MFA algorithm
(Figure~\ref{MFA-algorithm}) performs one update per neural variable 
(defined as one sweep) with sequential updating using equations
\ref{spe-1} and \ref{spe-2}. If the saturation $\Sigma$, has changed
more than $10\%$ after all neurons have been updated once, the system
is re-initialized with $T(0) \rightarrow 2 T(0)$. After reaching a
saturation value close to 1 (we choose $\Sigma = 0.99$) we check
whether the obtained solutions are valid, i.e. $E_{hard} = 0$. If this
is not the case the network is re-initialized and is allowed to
resettle. We repeat this procedure a number of times until the best
solution is found. A similar procedure was carried out on high school 
scheduling by Peterson {\it et al.}~\cite{GSP:complex}. 

The MFA implementation was a little more complicated 
than the implementation of simulated annealing and the expert system,
since it had many more parameters to handle, and it was 
often more difficult to find optimal values for these parameters.  
For example, one complication is the computation of the critical
temperature $T_c$, which involved an iterative procedure of a
linearized dynamic system.  On the other hand, we observed that the 
convergence time was indeed much less than any of the convergence
times of the simulated annealing using the three annealing schedules
studied.

For a full derivation of the mean-field annealing algorithm from its
roots in statistical physics, see Krogh, Hertz and
Palmer~\cite{HKP:neural}. 

\subsection{Definitions and Constraints}
Here we define all of the terms and concepts used throughout this
chapter for the implementation of class scheduling on the Potts
model. 

Let ($N_p$) to be a number of professors that give ($N_q$) classes in 
($N_x$) rooms at ($N_t$) time slots. Event-wise, $(p,q)$ event
-- denoting a professor $p$ giving a class $q$ -- takes place in the
space-time slot $(x,t)$ -- denoting room $x$ at time $t$.

\bd{dfn-0}
\ \ Let $S_{pq;xt}$ denote the Potts (or multi-state) neuron, where
\be{form-1}
S_{pq;xt} =
\twoif {1} {\text{if event $(p,q)$ takes place in $(x,t)$ ,}}
       {0} {\text{otherwise.}}
\ee
\ed

\bd{const-dfn.1}
\ \ We have dealt with two types of constraints:
\begin{enumerate}
\item Those involving professors, rooms, subjects, and time
periods.
\item Those dealing with students scheduling and preferences.
\end{enumerate}
In addition, some of these constraints are usually classified to be
{\bf hard}, others are {\bf medium} and the rest are {\bf soft}.
\ed

\bd{const-dfn.2}
Here, we define a number of {\bf hard constraints} of the first
type (see definition~\ref{const-dfn.1}) .
\begin{itemize}
\item An event $(p,q)$ should occupy precisely one space-time slot
$(x,t)$.
\end{itemize}
This constraint can be embedded in the neural net in terms of the
Potts normalization condition:
\be{form-2}
\sum_{x,t} S_{pq,xt} = 1
\ee
\begin{itemize}
\item Different events $(p_1,q_1)$ and $(p_2,q_2)$ should not occupy
the same space-time slot $(x,t)$.
\end{itemize}
\be{form-3}
E_{XT} =
{1 \over 2}\sum_{x,t}\sum_{p_1,q_1}\sum_{p_2,q_2}S_{p_1q_1;xt}S_{p_2q_2;xt} =  
{1 \over 2}\sum_{x,t}{\left(\sum_{pq}S_{pq;xt}\right)}^2
\ee
\begin{itemize}
\item A teacher $p$ should have at most one class at a time.
\end{itemize}
\be{form-4}
E_{PT} =
{1 \over 2}\sum_{p,t}\sum_{q_1,x_1}\sum_{q_2,x_2}S_{pq_1;x_1t}S_{pq_2;x_2t} = 
{1 \over 2}\sum_{p,t}{\left(\sum_{qx}S_{pq;xt}\right)}^2
\ee
\begin{itemize}
\item A class $q$ should have at most one teacher at a time.
\end{itemize}
\be{form-5}
E_{QT} =
{1 \over 2}\sum_{q,t}\sum_{p_1,x_1}\sum_{p_2,x_2}S_{p_1q;x_1t}S_{p_2q;x_2t} = 
{1 \over 2}\sum_{q,t}{\left(\sum_{px}S_{pq;xt}\right)}^2
\ee
\ed

In~\cite{Hop:network, HT:neural}, Hopfield introduced a model with neurons
having a graded response. So given $I_i$ a fixed input bias to neuron
$x_i$, $W_{ij}$ the weight of the connection between neurons
$(x_i,x_j)$, and parameter $T$ to denote the ${\it Temperature}$, the
behavior of $x_i$  can be stated as:

\be{form-JH}
 x_i =
 {1 \over 2} + {1 \over 2} \cdot 
 tanh\left(({\sum_{i \neq j}W_{ij}x_j - I_i}) \times {1 \over T}\right)
\ee
When $T \longrightarrow 0$, function ({\ref{form-JH})'s
response takes on a step-wise form. These models of~\cite{HT:neural}
are known to be a bit sensitive to parameter tuning for real world
applications, especially those of the type that involves ``finding''
or ``searching for'' $1$ out-of $N$ type problems.  

\bd{don-1}
\ \ At a finite $T$, $x_i$ (of eq. \ref{form-JH})
can be replaced by an MFT variable $\upsilon_{ia}$, then
differentiating to obtain the Potts MFT equations:

\be{form-5.1}
\upsilon_{ia} =
- {1 \over T} \cdot {{\partial E}[{\vec\nu}] \over {\partial \nu_{ia}}}
\ee
\be{form-5.2}
\nu_{ia} =
{e^{\upsilon_{ia}} \over {\sum_b e^{\upsilon_{ib}}}}
\ee

Which enforces

\be{form-5.3}
\sum_{a} \nu_{ia} = 1\,\,\,\, \forall i\, .
\ee
Where $a$ denotes the vector components or all the elements of the
group. Also, we could state that $\nu_{ia}$ $\approx$
Probability(neuron(i) $\in$ state(a)). Equations (\ref{form-5.1},
\ref{form-5.2}) are computed in an iterative manner.
\ed

\bd{dfn-2}
\ \ To carry out the minimization of the energy functions of
definition (\ref{dfn-0}) we need to convert the neural variables
$\{0,1\}$ of $S_i$ to their corresponding {\it mean field variables} 
at $T$. These correspondences are derived (see~\cite{GSP:complex}) to
be:

\be{form-6.1}
V_{pq;xt} \approx <S_{pq;xt}>_T\,\,;
\ee
which is the thermal average of the original binary spines, with the
following multi-state $MFT$ equations (see~(\ref{form-5.1},
\ref{form-5.2}) for a general format):

\be{form-6}
U_{pq;xt} =
- {1 \over T} \cdot {{\partial E} \over {\partial V_{pq;xt}}}
\ee
\be{form-7}
V_{pq;xt} =
{e^{U_{pq;xt}} \over {\sum_{x^{'}t^{'}}e^{U_{pq;x^{'}t^{'}}}}}
\ee
\ed

The $MFT$ equations of $\upsilon$ and $\nu$ are 
solved {\it iteratively} either synchronously or asynchronously, under
$T$. Also, these equations give rise to 2 phase  transitions; one in
$x$ (space) and the other in $t$ (time). This leads to factorization
of the system by replacing $S_{pq;xt}$ by {\bf $X$-neurons}
$S^{(X)}_{pq;x}$ and {\bf $Time$-neurons} $S^{(Time)}_{pq;t}$ :
 
\bd{dfn-3} \ \
\be{form-8}
S_{pq;xt} =
S^{(X)}_{pq;x} \cdot S^{(Time)}_{pq;t}
\ee
With separate Potts conditions:

\be{form-9}
\sum_{x}S^{(X)}_{pq;x} = 1
\ee
\be{form-10}
\sum_{t}S^{(Time)}_{pq;t} = 1
\ee

For {\bf $X$-neurons} and {\bf $Time$-neurons}, respectively. 
Rather than dealing with $p$ and $q$ directly, an independent variable
($i$) -acting as an event index- to which $p$ and $q$ are attributes
$p(i)$ and $q(i)$. The value of these are kept in a look-up table which
contains all the relevant information to process each event ($i$).
Now, the above (\ref{form-9}, \ref{form-10}) conditions can be 
re-stated as follows:

\be{form-9.1}
\sum_{x}S^{(X)}_{i;x} = 1
\ee
\be{form-10.1}
\sum_{t}S^{(Time)}_{i;t} = 1
\ee

\ed

During annealing, as we go from a higher temperature to a lower one,
a {\it phase transition} is passed at ($T = T_c$), where 
$T_c$\footnote{See pp 29-31 of~\cite{HKP:neural} for an approximate
value in the case of Ising ferromagnets.} is the {\it critical
temperature}. The value of $T_c$ set the scale or the range of the
value of $T$, and it is estimated by expanding \ref{form-7} around
some trivial fixed-point such as

\be{trivial-form}
V^{(0)}_{pq;xt} = {1 \over {N_x N_t}} \;\text{, at $T$ = 0}.
\ee

Also $T_c$ depends on the self-coupling coefficient $\beta$, and on
the method of the update.  This parameter $\beta$ acts as a
Lagrangian multiplier for the given energy constraint(s). See 
section~\ref{sec:MFAalg} and~\cite{PS:mapping}.

\bd{dfn-4}\ \
{\em Relative Clamping}: During the process of revising an existing
schedule, sometimes it is desirable to have a certain cluster of
events to stick together in time (or within some specified period of
time). This amounts to making the following notational replacement: 
\be{form-12}
i \longrightarrow (j,k) \;\; ;
\ee
where $(j)$ denotes the cluster and $(k)$ an event within $(j)$. In
this situation, one has a common $t$-neuron $(S^{(Time)}_{j;t})$ for
cluster $j$ but distinct $x$-neurons $(S^{(X)}_{jk;x})$ for the
individual events.
\ed

\bd{dfn-3.9}
To deal with events that take place over two consecutive hours, what
~\cite{PS:mapping} called effective $time$-neurons
$\tilde{S}^{(Time)}_{i;t}$ are introduced: 
\be{form-11}
\tilde{S}^{(Time)}_{i;t} =
\sum_{k=0}^{g_i - 1} S^{(T)}_{i;t-k} \;\;{\text{;}}
\ee
where $g_i = 1$ for single hours and $= 2$ for double hours. For those
events that occur between one and two hours, $g_i$ can be adjusted
accordingly to handle these cases.
\ed

Utilizing definition (\ref{dfn-3.9}), constraints (\ref{form-3},
\ref{form-4}, \ref{form-5}) can be put into the following form:

\bd{dfn-3.9.1}

\be{form-11.1}
 E_{XT} =
 {1 \over 2}\sum_{x,t}\sum_{i \neq i^{'}}
 S^{(X)}_{i;x}S^{(X)}_{i^{'};x} \tilde{S}^{(Time)}_{i;t}
 \tilde{S}^{(Time)}_{i^{'};t}
\ee

\be{form-11.2}
 E_{PT} =
 {1 \over 2}\sum_{p,t} \sum_{i \neq i^{'}} \delta_{p(i),p(i^{'})}
 \tilde{S}^{(Time)}_{i;t} \tilde{S}^{(Time)}_{i^{'};t}
\ee

\be{form-11.3}
 E_{QT} =
 {1 \over 2}\sum_{q,t} \sum_{i \neq i^{'}} \delta_{q(i),q(i^{'})}
 \tilde{S}^{(Time)}_{i;t} \tilde{S}^{(Time)}_{i^{'};t}
\ee
\ed

\bd{dfn-6}\ \
The Kronecker $\delta_{q(i),q(i')}$ ensures that only
events $i$ and $i'$ with identical ``primordial'' classes are,  for
example, summed over. Here, $q(i)$ and $q(i')$ refer to classes of the
same type or 2 sections derived from the same original class. 

\be{Kronecker}
\delta_{r,s} =
\twoif {1} {\text{if} r = s}
       {0} {\text{otherwise.}}
\ee
\ed

\bd{dfn-3.9.2}\ \
One important soft constraint that needs to be handled is {\it group
formation}. Briefly, it can be stated as follows: for some optional
subjects the classes are broken up into option groups temporarily
forming new classes. In order to account for this, we need a formalism
that allows for the breaking up of these primordial classes and
subsequent recombinations into {\it quasi}-classes. This is a common
practice with those classes where the students have many choices. 
To carry out such an extension to quasi-classes one might
have contributing pairs of events where the primordial
classes are different. Hence one should replace $\delta$ with
a more structured overlap matrix $\Gamma$ which is given by a look-up 
table.
\be{form-11.4}
\delta_{q(i),q(i^{'})} \longrightarrow \Gamma_{q(i),q(i^{'})} \; .
\ee
\ed

\bd{dfn-3.9.3}\ \
From definitions ($\ref{dfn-4}$, $\ref{dfn-3.9.2}$), we need to make
sure that ($\forall k$) events within a given ($j$) cluster 
$p(j,k) \neq 0$. This says that there can be no ``uncovered'' class
period or a class event (that is, periods when no professor has been
assigned to an already scheduled class).
For this purpose we introduce $\tilde{\delta}_{j,k}$
according to:
\be{form-3.9.3}
\tilde{\delta}_{j,k} = \delta_{j,k}(1 - \delta_{j,0}\delta_{k,0}) \; .
\ee
If $j = 0$ then $\forall k$ , $k = 0$ .
\ed

\bd{dfn-3.9.4}\ \
Utilizing the concepts of the above definition ($\ref{dfn-3.9.2}$,
$\ref{dfn-3.9.3}$), the collision terms of definition
($\ref{dfn-3.9.1}$) are extended to ensure that all the hard
constraints are satisfied (or the solutions obtained are legal).
The extension takes the following forms:

\be{form-3.9.4.1.1}
 E_{XT}{\text{(form 1)}} =
 {1 \over 2}\sum_{jk \neq j^{'}k^{'}}\sum_{x}
 S^{(X)}_{jk;x}S^{(X)}_{j^{'}k^{'};x}
\sum_{t}\tilde{S}^{(Time)}_{j;t}\tilde{S}^{(Time)}_{j^{'};t}
\ee

\be{form-3.9.4.1.2}
 E_{XT}{\text{(form 2)}} =
 {1 \over 2}\sum_{j \neq j^{'}}\sum_{k k^{'}}\sum_{x}
 S^{(X)}_{jk;x}S^{(X)}_{j^{'}k^{'};x}
 \sum_{t}\tilde{S}^{(Time)}_{j;t}\tilde{S}^{(Time)}_{j^{'};t}
 + {1 \over 2}\sum_{j}g_j \sum_{k \neq k^{'}} \sum_{x} 
 S^{(X)}_{jk;x}S^{(X)}_{jk^{'};x}
\ee

\be{form-3.9.4.2.1}
 E_{PT}{\text{(form 1)}} =
 {1 \over 2}
\sum_{p,t} \left[ \sum_{jk} \tilde{\delta}_{p,p(j,k)}
 \tilde{S}^{(Time)}_{j;t} \right]
\left[ \sum_{j^{'}\neq j}\sum_{k^{'}}
\tilde{\delta}_{p, p(j^{'},k^{'})} 
 \tilde{S}^{(Time)}_{j^{'};t} \right]
\ee

\be{form-3.9.4.2.2}
 E_{PT}{\text{(form 2)}} =
        {1 \over 2}\sum_{j \neq j^{'}}\sum_{k k^{'}} 
        \tilde{\delta}_{p(j,k), p(j^{'},k^{'})}
        \sum_{t}\tilde{S}^{(Time)}_{j;t} \tilde{S}^{(Time)}_{j^{'};t}
\ee

\be{form-3.9.4.3.1}
 E_{QT}{\text{(form 1)}} =
 {1 \over 2} 
\sum_{qt} \left[ \sum_{jk} \delta_{q,q(j,k)} \tilde{S}^{(Time)}_{j;t}
\right] 
\left[ \sum_{j^{'} \neq j} \sum_{k^{'}} \Gamma_{q,q(j^{'}, k^{'})}
        \tilde{S}^{(Time)}_{j^{'};t} \right]
\ee

\be{form-3.9.4.3.2}
 E_{QT}{\text{(form 2)}} =
        {1 \over 2}
        \sum_{j \neq j^{'}}\sum_{kk^{'}} 
        \Gamma_{q(j,k), q(j^{'},k{'})} \sum_{t} \tilde{S}^{(Time)}_{j;t}
        \tilde{S}^{(Time)}_{j^{'};t}
\ee

Note: Each ``form 2'' above is a derivation from its corresponding
``form 1''. 
\ed

The factorization idea of (dfn \ref{dfn-3.9.4}) can be extended to
{\it mean field} variables
$V^{(X)}_{jk;x}$ and $V^{(Time)}_{j;t}$, corresponding to neurons
$(S^{(X)}_{jk;x})$ and $(S^{(Time)}_{j;t})$, respectively.
$V^{(X)}_{jk;x}$ and $V^{(Time)}_{j;t}$ can be interpreted as being 
the probabilities that the corresponding events occur
at given $x$-value and $t$-value with normalization given by the above
equations (\ref{form-9.1}, \ref{form-10.1}) .

Substituting $V$ for $S$ in (\ref{form-9.1}, \ref{form-10.1}), 
the {\it mean field} equations at a finite 
$T$, are given in terms of the {\it local fields} 
$U^{(X)}_{jk;x}$ and $U^{(Time)}_{j;t}$. Now, the {\it mean field} 
dynamics proceed over the following local and mean field variables:

\bd{dfn-3.99}
\be{form-4.1.3}
U^{(X)}_{jk;x} =
- {1 \over T} \cdot {{\partial E} \over {\partial V^{(X)}_{jk;x}}}
\ee

\be{form-4.1.4}
U^{(Time)}_{j;x} =
- {1 \over T} \cdot {{\partial E} \over {\partial V^{(Time)}_{j;t}}}
\ee

\be{form-4.1.5}
V^{(X)}_{jk;x} =
{e^{U^{(X)}_{jk;x}} \over {\sum_{x^{'}} e^{U^{(X)}_{jk;x^{'}}}}}
\ee

\be{form-4.1.6}
V^{(Time)}_{j;t} =
{e^{U^{(Time)}_{j;t}} \over {\sum_{t^{'}} e^{U^{(Time)}_{j;t^{'}}}}}
\ee

The iterations of the {\it MF} dynamics over equations 
({\ref{form-4.1.3}} - {\ref{form-4.1.6}})
can either be {\bf synchronous} or {\bf serial} (see~\cite{GSP:complex}).
\ed

\bd{dfn-4.2}
{\em Synchronous Updating}: \ \
A single sweep consists of first computing all the local fields, and
then updating all the neurons. That means do equations
({\ref{form-4.1.3}},{\ref{form-4.1.4}}) then perform equations
({\ref{form-4.1.5}},{\ref{form-4.1.6}}). 
\ed

\bd{dfn-4.3}
{\em Serial Updating}: \ \
For each neuron $S^{(X)}_{jk;x}$ and/or $S^{(Time)}_{j;t}$, the {\it
local field} is computed immediately before the corresponding {\it
neuron state} is updated. That means either do equation
({\ref{form-4.1.3}}) followed by ({\ref{form-4.1.5}}) and/or do
equation ({\ref{form-4.1.4}}) followed by ({\ref{form-4.1.6}}).
\ed

Next, we need to define 2 important {\bf soft constraints} encountered 
when scheduling classes.


\bd{dfn-5}\ \
\begin{itemize}
\item {\em Spreading}: different classes in a particular area
(i.e. computer science) should be spread over the week in a (MWF) and
(TuTh) type of spread. 
\end{itemize}

First, assume we have access to a subject attribute $s(j,k)$. Second,
we define an effective ``day-neuron'' as follows  
($\tilde{S}^{(Time)}_{j;dh}$ is reduced to a $D$-type neuron by
summing over  all hours of the day):

\be{day-neuron}
S^{(D)}_{i;d} = \sum_h \tilde{S}^{(Time)}_{j;dh} = g_j \cdot
S^{(Time)}_{j;dh}
\ee

So the penalty term $E_{QSD}$ that spreads classes $Q$ in a particular
subject $S$ over different week-days $D$ :

\be{spread-1}
E_{QSD} = 
{1 \over 2} \sum_{qsd}{\left[\sum_{jk}\sum_{h} \delta_{q,q(j,k)}
\delta_{s,s(j,k)} \tilde{S}^{(Time)}_{jk;dh}\right]}^2
\ee

Refining the above a bit more by considering 2 different clusters
$(j)$ and $(\hat{j})$ of classes over $D$ : ($q$ and $s$ are not used
in the summation, since no change occurs to them) 

\be{spread-2}
E_{QSD} =
{1 \over 2} \sum_{d}\sum_{j \neq \hat{j}}\sum_{kh}\sum_{\hat{k}\hat{h}}
\delta_{q(j,k),q(\hat{j},\hat{k})}\delta_{s(j,k), s(\hat{j},\hat{k})}
\hat{S}^{(Time)}_{j;dh}\hat{S}^{(Time)}_{j;d\hat{h}}
\ee

Now for some day $(d)$ of the week:

\be{spread-3}
E_{QSD} =
{1 \over 2} \sum_{j \neq \hat{j}} \sum_{k\hat{k}} 
\delta_{q(j,k),q(\hat{j},\hat{k})}\delta_{s(j,k), s(\hat{j},\hat{k})}
\sum_d S^{(D)}_{j,d} S^{(D)}_{\hat{j},d}
\ee

\begin{itemize}
\item {\em Glueing}: The schedule should have as few ``holes'' as
possible; that is lessons should be glued together. 
\end{itemize}

The penalty term that rewards situations where $(d,h)$-events are
glued to $(d,h-1)$-events can be defined as follows:

\be{glueing-1}
E_{QDH} = -\alpha \sum_{q}\sum_{d,h}
\left[\sum_{j,k}\delta_{q,q(j,k)}\tilde{S}^{(Time)}_{j;dh}\right]
\left[\sum_{\hat{j} \neq j}\Gamma_{q,q(\hat{j},\hat{k})}
\tilde{S}^{(Time)}_{j;d(h-1)}\right]
\ee

The energy term
of (\ref{glueing-1}) can be summarized to say that: {\it if cluster
($j$) with events each symbolized by ($k$) occurs on hour ($h$) of day
($d$) then it would be quite desirable to have another event
($\hat{k}$) from a different cluster ($\hat{j}$) to occur on hour
($h-1$) of that day.} 

Next, since ($q$) is fixed we could re-write the above equation as
follows:

\be{glueing-2}
E_{QDH} = -\alpha\sum_{\hat{j} \neq j}\sum_{k,\hat{k}}\Gamma_{q(j,k),
q(\hat{j},\hat{k})}
\sum_{d,h>1}\tilde{S}^{(Time)}_{j;dh}\tilde{S}^{(Time)}_{\hat{j},d(h-1)}
\ee

The presence of the ($-$) sign in front of (\ref{glueing-1}, 
\ref{glueing-2}) indicates that such energy functions are to be
minimized. Also, the ($\alpha$) in front of the the above forms
governs the strength of this reward relative to the energies of
equations  (\ref{form-3},\ref{form-4},\ref{form-5}). 
This ($\alpha$) takes on a value between zero and $1$.

Finding {\em legal} solutions with no collisions corresponds to
minimizing the {\em hard} energy:
\be{Hard}
E_{hard} = combine (E_{XT}, E_{PT}, E_{QT})
\ee
to zero. The {\em soft energy} is given by:
\be{Soft}
E_{soft} = combine (E_{QSD}, E_{QDH})\;.
\ee
The total energy $E$ to be minimized is then:
\be{total}
E = combine (E_{hard}, E_{medium}, E_{soft})\;.
\ee
\ed

\subsection{More Definitions and Constraints}

Let $q(i)$ represents a section of class $q$, $O_{q(i)}$ indicates 
that section $q(i)$ is open, and $C_{q(i)}$ is closed. Also, let 
$b(x,t)$ be a balanced interval for some time period $(t)$, $R$ is the
overall set of students, and an arbitrary given subset is  
$r \subseteq R$.

Here, we need to introduce the $R$-neuron: $S^{(R)}_{r(i);(d(t)d(x)}$
where $k= 1,2,\ldots$. Each student $r(i) \in R$ should have a
balanced schedule over the two-group of days $(MWF)$ and
$(TuTh)$. We will also use $r_s$ to denote a single student.

For $X$ and $T$ neurons, let $x(j)$ and $t(j)$ where $j= 1,2,\ldots, n$
stand for places and time slots (or events), respectively. 
Using such notation, an example of the Potts normalization will be:
$S^{(X)}_{p(i)q(i);x(j)}$ and $S^{(Time)}_{p(i)q(i);t(j)}$ respectively,
where $i=1,2,\ldots,n$ is a section index.

\bd{athletes}\ \
{\it The athletes case}: Let $\vec{r} \subseteq R$ denotes the group
of athletes  that are dealt with in this study. 
Let ($r^{*}$,$b$) be playing/training events scheduled during
the weekdays at (or during) times when academic (attending lectures)
events ($r^{+}$,$a$)  are also scheduled, or may possibly occur. Let
($d(x)$, $d(t)$) be a place/time during the week at which any of these
events can occur. 

In addition to the above stated Potts neurons, we also will be using :

\begin{enumerate}
\item $S^{(\vec{r},X)}_{r^{*}b;d(x)}$ ,
$S^{(\vec{r},Time)}_{r^{*}b;d(t)}$ : It takes a value of $1$ if
event ($r^*,b$) takes place in the place and time slot, respectively,
of a weekday $d$, and $0$ otherwise.

\item $S^{(\vec{r},X)}_{r^{+}a;d(x)}$ ,
$S^{(\vec{r},Time)}_{r^{+}a;d(t)}$ : It takes a value of $1$ if
event ($r^+,a$) takes place in the place and time slot, respectively,
of a weekday $d$, and $0$ otherwise.

\item $S^{(\vec{r},X)}_{jk;d(x)}$ : It takes a value of $1$ if event
($k$) of cluster ($j$) that is associated with the $x$-neuron of
$\vec{r}$ occurs in place $d(x)$, and $0$ otherwise. 

\item $S^{(\vec{r},X)}_{\hat{j}\hat{k};d(x)}$ : It takes a value of
$1$ if event ($\hat{k}$) of cluster ($\hat{j}$) that is associated
with the $x$-neuron of $\vec{r}$ occurs in place $d(x)$, and $0$
otherwise. 

\item $\tilde{S}^{(\vec{r},Time)}_{j;d(t)}$ : 
It takes a value of $1$ if cluster ($j$) that is associated with
the $t$-neuron of $\vec{r}$ occurs in time $d(t)$, and $0$
otherwise.

\item $\tilde{S}^{(\vec{r},Time)}_{\hat{j};d(t)}$ : 
It takes a value of $1$ if cluster ($\hat{j}$) that is associated
with the $t$-neuron of $\vec{r}$ occurs in time $d(t)$, and $0$
otherwise. 
\end{enumerate}

In the above list, the last two neurons act as what Peterson et al. calls {\it
effective} $time$-neurons. Not only will they handle periods of single
hours but also variants of longer and shorter intervals.
\ed

\bd{time-int.dfn}\ \
Set ($\epsilon$) and ($c$) to some desired values depending on the
time period of the set of classes, now, the above {\it effective} 
$t$-neurons can be defined as follows:

\be{time-int.1}
\tilde{S}^{(\vec{r},T)}_{j;d(t)} = 
\sum^{g_j - c}_{k=0} S^{(\vec{r},T)}_{j;d(t-k+\epsilon)}
\ee

$\tilde{S}^{(\vec{r},Time)}_{\hat{j};d(t)}$ can be defined in a
similar manner. If ($c = 1$) and ($g_i = 2$), then we are dealing with
double-hour intervals so we set ($\epsilon = 0$). If ($c = 1$) and
($g_i = 1$) then ($\epsilon = 0$) and

\be{time-int.2}
\tilde{S}^{(\vec{r},Time)}_{j;d(t)} = S^{(\vec{r},Time)}_{j;d(t)}.
\ee
\ed

\bd{consts-dfn}\ \
Let $S^{(R)}_{q(i)r(i);b(x,t)}$ for balancing the enrollment in a
multi-section courses. Section $q(i)$ of students $r(i)$ are balanced
over space-time $b(x,t)$. 

To enroll student $r(i)$ in a particular section $q(i)$ given the tag
$EC$ to indicate the status of the student,  
we have $S^{(X[Time])}_{q(i)r(i);(d(x)[d(t)])}|EC\in \{0,1\}$  is set
to$0$ if $EC$ is $0$, otherwise, it is set to $1$ and the action
$q(i)r(i)$ occurs during place $d(x)$ and/or time $d(t)$. 

Let $lin(c_i,c_j)$ representing the concept that course $c_i$ is
linked to some other course $c_j$ ; in this case, $c_i$ is a
pre-requisite to $c_j$. The degree of ``linkage'' varies between $0$
and $1$. It takes a value of $1$ if under all circumstances, course
$c_i$ must be taken before $c_j$. 

To deal with this linkage concept, we assumed in definition (\ref{dfn-5}) 
that we have an attribute $s(j,k)$ to represent a subject. For a
strong linkage, we will assume that we're dealing with subjects
that belong to the same cluster ($j$) (e.g. physics courses), so the 
link between various events of ($j$) is: $lin(s(j,k), s(j,\hat{k}))$.

Other degrees of linkage can involve subjects or courses from
different groups or clusters (e.g. physics and math courses, or 
engineering and math courses), and this connection is in the 
form: $lin(s(j,k), s(\hat{j},\hat{k}))$.

Let $\hat{x}(i)$ be a single location (or a seat) in class
$x(i)$ such that $i=1,2,\ldots, num$ ; where ($num$) is the 
number of class-rooms.
\ed

\subsubsection{Day-Class Balancing}

Each student of any particular section (i.e. $r(k)$) should have a
balanced schedule between ({\it MWF}) and ({\it TuTh}). The balancing
operation is accomplished by spreading as much as possible the courses
of the student over each of ({\it MWF}) and ({\it TuTh}).

Let $D$ be the set $\{\{MWF\}, \{TuTh\}\} = \{ D_1, D_2 \}$. Assume
the existence of the subject attribute $s(j,k)$. We introduce the
energy function $E_{QSD}$ that spreads the classes ($Q$) in a
particular field or subject ($S$) over the  set $\{D_1, D_2\}$.

\bd{balance-dfn.1}
The Potts condition of balancing $S^{(R)}_{q(m)r(m);b(x,t)}$
(take a value of $1$ for balanced and $0$, otherwise) is accomplished
through one (or all) of these energies:

\be{balance.1}
E^{r}_{QSD} =
{1 \over 2}\sum_{j \neq \hat{j}}\sum_{k \hat{k}} 
\delta_{q(j,k),q(\hat{j},\hat{k})}
\delta_{s(j,k),s(\hat{j},\hat{k})}
\sum_{d_1,d_2} S^{(D)}_{j,d_1} S^{(D)}_{\hat{j},d_2}
\ee

\be{balance.2}
E^{r}_{QSD_1} =
{1 \over 2}\sum_{j_1 \neq \hat{j_1}}\sum_{k_1 \hat{k_1}} 
\delta_{q(j_1,k_1),q(\hat{j_1},\hat{k_1})}
\delta_{s(j_1,k_1),s(\hat{j_1},\hat{k_1})}\sum_{d_1}
\sum_{i \neq \hat{i}} S^{(D)}_{i,d_1} 
S^{(D)}_{\hat{i},d_1}
\ee


\be{balance.3}
E^{r}_{QSD_2} =
{1 \over 2}\sum_{j_2 \neq \hat{j_2}}\sum_{k_2 \hat{k_2}} 
\delta_{q(j_2,k_2),q(\hat{j_2},\hat{k_2})}
\delta_{s(j_2,k_2),s(\hat{j_2},\hat{k_2})}\sum_{d_2}
\sum_{i \neq \hat{i}} S^{(D)}_{i,d_2} 
S^{(D)}_{\hat{i},d_2}
\ee
\ed

\subsubsection{The Case of Student Athletes}
The main constraint of this case very much analogous to the hard
constraint stated somewhere above : {\it different events should not
occupy the same space-time slot}.

\bd{athletic}

\be{athletic.1}
E^{\vec{r}}_{d(XT)} =
{1 \over 2}\sum_d\sum_{xt}\sum_{\dot{r}b}\sum_{\ddot{r}a} 
S^{(\vec{r},X[Time])}_{\dot{r}b;d(x)[d(t)]} 
S^{(\vec{r},X[Time])}_{\ddot{r}a;d(x)[d(t)]}
\ee

Where $(\dot{r}b)$ and $(\ddot{r}a)$ are some two different events.
After some further refinement, we obtain:

\be{athletic.2}
E^{\vec{r}}_{d(XT)} =
{1 \over 2}\sum_{j \neq \hat{j}}\sum_{k \hat{k}}\sum_{d(x)}
S^{(\vec{r},X)}_{jk;d(x)} S^{(\vec{r},X)}_{\hat{j}\hat{k};d(x)}
\sum_{d(t)} \tilde{S}^{(\vec{r},Time)}_{j;d(t)} 
\tilde{S}^{(\vec{r},Time)}_{\hat{j};d(t)} + 
{1 \over 2} \sum_j g_j \sum_{k \neq \hat{k}}\sum_{d(t)}
S^{(\vec{r},X)}_{jk;x} S^{(\vec{r},X)}_{j\hat{k};x} 
\ee

Since it is always the case that athletic training facilities 
are kept separate from the academic class-rooms and lecture halls, 
a more refined energy function is obtained:

\be{athletic.3}
E^{\vec{r}}_{d(Time)} =
{1 \over 2}\sum_{j \neq \hat{j}}\sum_{d(t)}
\tilde{S}^{(\vec{r},Time)}_{j;d(t)} 
\tilde{S}^{(\vec{r},Time)}_{\hat{j};d(t)}
\ee
\ed

\subsubsection{Eligibility for Enrollment}
\bd{The Eligibility Criteria}
Before going through the process of balancing, each student 
($r(i) = r_s$) eligibility for enrolling in a given class needs to be
determined.
Let the student's eligibility criteria be denoted by ($EC$), using 
the notation stated above, we could derive the
following energy function ---note: if $EC$ is set to ``No'' then it
takes a $0$ value in the energy function, otherwise, it is a $1$ : 

\be{eligibility-criteria.1}
E_{QR} = 
{1 \over 2}\sum_{qr}\sum_{jk}\sum_{i}
\delta_{q,q(j,k)}\delta_{r,r(i)} S^{X}_{ij;x(i)}|({\it EC\in\{0,1\}})
\ee

\be{eligibility-criteria.2}
\sum_{i} O_{q(i)}=1
\Rightarrow \text{Test:}\left\{ EC(R) + \text{satisfy-constraints}(Hard,Soft)\right\}
\ee
\ed

\subsubsection{Course Pre-requisite and Subsequent Assignment}

Utilizing the link concept $lin(s,s')$ between two different
subjects (see definition~\ref{consts-dfn}), 
a general energy assignment can be put together as follow: 

\bd{pre-sub.dfn}\ \
Let ($i$) stands for the situation at hand (also can be viewed as 
a space-time event).

\be{pre-requisite.1}
E_{link(1)} =
 -\gamma \sum_{is}{\left[\sum_{jk} \delta_{s,s(j,k)}
\tilde{S}^{(X[Time])}_{ln(s,s');i}\right]}^{2}
\ee

\be{pre-requisite.2}
E_{link(2)} =
 -\gamma \sum_{s}\sum_{j \neq j'}\sum_{k \neq k'} \delta_{s(j,k),s(j',k')}
\sum_{i} \tilde{S}^{(X[Time])}_{ln(s(j,k),s(j',k')); i}
\tilde{S}^{(X[Time])}_{ln(s(j,k),s(j,k')); i}
\ee

\be{pre-requisite.3}
E_{link(3)} =
 -\gamma \sum_{j \neq j'}\sum_{k \neq k'} \delta_{s(j,k),s(j',k')}
\tilde{S}^{(X[Time])}_{ln(s(j,k),s(j',k')); i}
\tilde{S}^{(X[Time])}_{ln(s(j,k),s(j,k')); i}
\ee

If ($\gamma > 0$) then function minimization is the target here,
otherwise maximization will take place.

\ed

\subsubsection{Enroll in all parts  of a course}
There are courses of different types that not only
have lecture periods but also associated sections that consist
of either labs or recitations or both. Signing up for such a course
entails enrolling in associated sections of it. 

Let us assume the existence of some subject attribute $s(j,k)$
consisting of two components: $\dot{s}$ denoting lectures, and
$\ddot{s}$ stands for labs. Now the energy equation depicting
this assignment can be outlined as follow:

\be{enroll-in-all-parts.1}
E_{enroll(1)} =
{1 \over 2}\sum_{s}\sum_{jk}\sum_{i} \delta_{s,s(j,k)}\delta_{r,r(i)} 
\sum_{t} S^{(Time)}_{jk;t}
\ee

\be{enroll-in-all-parts.2}
E_{enroll(2)} =
{1 \over 2}\sum_{si} \sum_{t\neq t'}\left[\sum_{jk} 
\delta_{s,\dot{s}(j,k)} \delta_{r,r(i)}
S^{(Time)}_{jk;t}\right]\left[\sum_{jk'} \delta_{s,\ddot{s}(j,k')}
\delta_{r,r(i)} S^{(Time)}_{jk';t'}\right]
\ee

\be{enroll-in-all-parts.3}
E_{enroll(3)} =
{1 \over 2}\sum_{k \neq k'}\delta_{s,\dot{s}(j,k)}
\delta_{s,\ddot{s}(j,k')}
\sum_{t \neq t'} S^{(Time)}_{jk;t} S^{(T)}_{jk';t'}
\ee

\subsubsection{Lecture hall capacity}
In order to satisfy this constraint, the following inequality
need to hold.
\be{part-1}
\forall x \sum_{i}||x(i)| - \sum_{j}\sum_{k} S^{(X)}_{jk;x(i)}| \geq 0
\; .
\ee

\subsubsection{Student status and priority mapping}
\be{status-and-priority-mapping}
E_{assign} =
{1 \over 2}\sum_{m, \hat{x}, x}\sum_{i\neq i^{'}}
S^{(X)}_{\hat{x};x} 
\left( S^{(R)}_{r(i);(\hat{x}(m) \in x(m))}|({\it M = Req})
+
S^{(R)}_{r(i^{'});(\hat{x}(m) \in x(m))}|({\it M = \chi})\right)
\ee

\subsubsection{On the Energy Equation of the Potts Model}
The energy function of the Potts spin model is of the form:
\be{enrg}
E = E_{part(1)} + E_{part(2)} + E_{part(3)} + \ldots + E_{part(n)}
\ee

where each component of this function stands for each of the energy
functions defined previously.

In~\cite{PS:mapping}, in which the TSP is used as study case,
parameters $\alpha$ and $\beta$ are used to control the strength of
E's restraints. The same parameters can also be used when dealing with
non-Euclidean type combinatorial optimization problems such as scheduling. 

\bo{ob-1}
\ \ For timetabling, $\beta$ is split into two parts ($\beta^{(Time)}$
and $\beta^{(X)}$), one with respect to the time-based neurons and the
other with respect to the space-based neurons. $E_{\beta}$ becomes as
follows: 
\be{diagonal-energy}
- {\beta^{(Time)} \over 2}\sum_j\sum_t {[S_{j;t}^{(T)}]}^2 - 
{\beta^{(X)} \over 2}\sum_{jk}\sum_{x} {[S_{jk;x}^{(X)}]}^2 
\ee
Here, the appearance of $\beta^{(Time)}$ and $\beta^{(X)}$ is only in
the diagonal terms. They play a crucial role in monitoring the
dynamics (phase transitions, etc.). These diagonal terms correspond to
self coupling interactions for the neurons of the network. In our 
simulations, each is set to $-0.1$.
\eo

\subsection{General Comments on the Algorithm}
Let ($N$) denotes the total number of the utilized neurons. At high 
$T$, the mean field solutions will tend to be states near to the
(symmetrical) maximum entropy state: $\forall (i, a) v_{ia} = {1/N}$.
Conversely at low $T$, finding a mean field solution will be
equivalent to using a local optimization method on the internal energy
--- a procedure highly sensitive to the initial conditions and known
to be ineffective (see p. 149 of~\cite{HT:neural}).

\subsection{Some Implementation Details}
The general idea is to initialize the system close to the symmetric  
state at a temperature near $T_c$, then, allow the system
to converge to a fixed point at each temperature, anneal (i.e. reduce 
$T$) until either the system reaches a stable high saturation
state or some other termination criterion is satisfied.
Specifically, the implementation used here has the following
ingredients\footnote{These follow the method used in
~\cite{PS:mapping} but with some additions/improvements done by the
authors -see algorithm above.}\footnote{Martin Simmen --in a private
communication-- investigated in depth these and most of the other
issues pertaining to the Potts model, using the TSP as a test case.}:

\begin{itemize}
\item {\bf Initial State}\ \
Another possible initial value of $V^{n=0}_{ia}$ is $(1 + \xi_{ia})/N$,
with the $\xi_{ia}$ being random variables drawn uniformly from 
${[-\xi, \xi]}$. Clearly $\xi$ should be $\ll 1$.

\item {\bf Convergence}\ \
At each $T$ the $V_{ia}$ variables are updated as above or by
the following rule\footnote{Thanks to M. Simmen for pointing this out.
Also see~\cite{PS:mapping}}:
\be{convergence-rule}
V^{(n+1)}_{ia} = V^{(n)}_{ia} + \gamma(f_{ia} - V^{(n)}_{ia}),
\;\; 0 < \gamma \leq 1\; .
\ee
Setting $\gamma = 1$ recovers the dynamics of~\cite{PS:mapping},
and the $\gamma \rightarrow 0$ limit gives the first order Euler
integration method for the differential equation:
\be{diff-equation}
d V_{ia} / dt = f_{ia} - V_{ia}\;.
\ee
The fixed points (also referred to as the free energy minima)
of this equation constitute the mean field solutions.
The degree of convergence is monitored through
a ``tolerance'' quantity $\Delta^{expected}$ of the form 
$\sum_{ia}{(V_{ia} - f_{ia})}^2$.
If $\Delta^{expected} < \Delta^{'}$, where $\Delta{'} \equiv \Delta (0.05/N)$, the system is deemed to have reached a fixed point and the
temperature is reduced. $\Delta$ will henceforth be called the
tolerance parameter.

\item {\bf Annealing Schedule} \ \
A simple exponential scheme is used, i.e. after reaching a fixed point
at $T$, $T \mapsto T \cdot T_r$, where 
$T_r \in (0,1)$ controls the cooling rate. Roughly speaking,
$T_r$ should be set to $0.9$.

\item {\bf Termination}\ \
One criterion is as outlined above in the algorithm. The other two
involve terminating runs which either fail to reach a fixed point
within a certain number of sweeps at some temperature, or which allow
$T$ to fall below some $T_{min}$, ($T_{min} \ll T_c$), without ever
satisfying the $\Sigma > \Sigma^{thresh}$ 
criterion. $T_{min}$ was set to $10^{-3}$.

\end{itemize}

\subsection{Performance Issues of MFA on Timetabling}

Given $M_x$ = average number of rooms for each category of classes,
$N_x$ = the overall total number of rooms, and $N_t$ = total number of
slots , putting the figures of the students aside, with a computational
complexity of O($M_x N_x N_t$) per sweep, software implementation of
the MFA are unlikely to be competitive with conventional serial
algorithms such as 3-opt and Lin-Kernighan's (for TSP like problems),
which empirically have an overall O($N^2$) running
time~\cite{JM:traveling}. Ideally, one would like to study the likely
performance of the method through numerical simulations using a small
step-size $\gamma$ to approximate continuous time dynamics. However,
this is computationally very expensive: not only do the simulation
times grow as $1 \over \gamma$, but also small values of ($\gamma$)
accentuate the difference between the expected value $T_c^{exp}$ of
the critical temperature $T_c$ and the theoretically computed value
$T_c^{theo}$. That is, low $\gamma$ runs tend to depress $T_c^{exp}$
artificially and thereby degrade solution quality, and as we observed,
result in a number of constraint violations for mainly soft
constraints and a number of hard ones as well. This latter point can
be compensated for by more stringent annealing parameters, but at the
cost of longer simulation times. No attempts have been made to carry
this out. Also, we observe that the overall quality of the solution
degrades as the problem size increases. This degradation moves from
soft constraint violation to major violations of hard constraints, and
hence, illegal timetables. 

\subsection{Mapping Timetabling onto the Potts Neural Nets}

\subsubsection{The Relationship Between the Key Timetabling Entities}

Our problem consist of the following sets: teachers (P), classes
(C), students (S), rooms (R), and time periods (I).
\vspace{0.1cm}
\begin{figure}[htb]
\begin{center}
\begin{picture}(200,200)
\psfig{figure=scheduling-net.eps,height=200pt,width=200pt}
\end{picture}
\end{center}
\caption{A one-way arrow indicates that the ``end-of-arrow'' entity
depends on the ``begin-of-arrow'' entity; while a two-way arrow
indicates a mutual dependency between the two entities.}
\label{relations}
\end{figure}

\subsubsection{Neural Mapping}

Involve mapping the course scheduling problem onto the Potts model.
\begin{itemize}
\item These multi-state neurons, forming the network, represent
different activities of the timetabling problem. 

\item Construct the constraints of the problem using these neurons.

\item Derive those MF equations controlling the dynamics of the net:

\begin{itemize}
\item Introduce the MF variables , for example, 

\be{rel.1}
V_{pq;xt} = <S_{pq;xt}>_T\;\;.
\ee

\item Start off with the partition function:
\be{part}
Z = {\displaystyle \sum_{[\vec{s}]}} e^{-E[\vec{s}] / T} \;.
\ee

\item Re-write $Z$ in terms of a multidimensional integral
over continuous vectors ($\vec{u}_i, \vec{v}_i$).

\item The MFT approximation to $<\vec{s}_i>$ is given by the value 
$\vec{u}_i$ at the saddle point equations:\\

\be{mft.1}
{{\partial E} \over {\partial v_{ia}}} = 0
\ee
and
\be{mft.2}
{{\partial E} \over {\partial u_{ia}}} = 0\;,
\ee
generating a set of self-consistency equations in the MF variables:

\be{mfv.1}
u_{ia} = 
- {1 \over T} {{\partial E[\vec{v}]} \over {\partial v_{ia}}}\;,
\ee
and
\be{mfv.2}
v_{ia} = {e^{u_{ia}} \over {{\displaystyle \sum_b} e^{u_{ib}}}}\;.
\ee

Where $a$ denotes the vector components.

For events (pq, xt), we get the following MF variables:

\be{mfv.3}
u_{pq,xt} = 
- {1 \over T} {{\partial E} \over {\partial v_{pq,xt}}} \;,
\ee
and
\be{mfv.4}
v_{pq,xt} = 
        {e^{u_{pq,xt}} \over {{\displaystyle \sum_{x' t'}}
e^{u_{pq,x't'}}}} \;.
\ee

From the latter follows, that ${\displaystyle \sum_{a}} v_{ia} =
1\;\;$. One can think of the mean field $v_{ia}$ as the probability
for the Potts neuron $i$ to be in state $a$.
\end{itemize}
\item Iterate $v_{pq, xt}$ equation with annealing, i.e. starting at a
high temperature, and successively lowering it in the course of the
process.  A phase transition is passed at $T = T_c$.

\item Estimate the critical temperature $T_c$, which sets the scale of
$T$, by expanding the above $v_{pq,xt}$ equation around the trivial
fix-point, at $T = 0$ :
\be{fixed}
v^{(0)}_{pq,xt} = {1 \over N_x N_t}\;.
\ee
\end{itemize}


%% \section{Experimental Results}
%%
%% \vspace{0.25cm}
%% \begin{figure}[tbp]
%% \centerline{\psfig{figure=scheduled-courses.ps,height=330pt,width=350pt,angle=-90}}
%% \vspace{0.75cm}
%% \caption{Percentage of scheduled classes versus number of spare
%% room/time slots}
%% \label{class-scheduling}
%% \end{figure}
%%
%% Figure~\ref{class-scheduling} shows three sets of problems, each has
%% $m$ number of pofessors, $n$ number of classes, $k$ number of
%% classrooms (including auditoriums, seminar rooms, etc.), and $l$
%% %% number of time slots. The values of $k$ and $l$ are the same (fixed)
%% for the three problems but each has different values of $m$ and
%% $n$. The largest has an $m = 1190$ lecturers (including visiting
%% professors and teaching assistants), and $n = 3839$ classes of the
%% first semester as shown in Table~\ref{table.01}. The medium size
%% problem has an $m = 400$ and $n = 1600$, while the smallest has and $m
%% = 150$ professors and $n = 450$ classes. At the start of scheduling
%% when all space-time slots are open, the three sets, as shown in the 
%% Figure~\ref{class-scheduling}, would be solvable. But as the available 
%% number of room/time slots reduces, the three sets would become harder
%% and harder to fully schedule. On other words as more space
%% and time resources get used the problem sets become harder to schedule
%% to satisfy the given constraints. Also, from our experiments we observed
%% that the less sparse the set is, the more difficult it is to
%% obtain a full schedule as space and time slots get reduced.


\section{Combinatorial applications of optimization networks}

In addition to the classic benchmarks like the TSP, optimization
networks have also been proposed for the solution of many more
practical problems. Although some of the papers cited below do not
describe proper mappings or effective annealing procedures, they
nonetheless provide a useful index of possible combinatorial
applications of optimization networks. These include decoding error
correcting codes~\cite{BB:neural}, image
segmentation~\cite{GPP:image}, stereo
correspondence~\cite{Tresp:neural}, load
balancing~\cite{FF:balancing}, classroom
scheduling~\cite{GSP:teachers}, high school course
scheduling~\cite{GSP:complex}, multiprocessor
scheduling~\cite{HK:asymmetric}, synthesis of digital
circuits~\cite{UDP:solving}, invariant pattern
recognition~\cite{BD:elastic}, analogue-to-digital
conversion~\cite{TH:neural}, data rearrangement~\cite{IN:data},
assignment~\cite{EDKT:neural, TCP:neural}, communication link
scheduling~\cite{OB:neural}, trackfinding~\cite{Fox:note}, navigation
problems~\cite{GF:solving}, clustering~\cite{RGF:constrained,
RGF:deterministic, RGF:statistical, Fox:optimization}, vector
quantization~\cite{RGF:vector}, and implementing the viterbi algorithm for 
Hidden Markov Models~\cite{AF:network}.













