Algorithms and Heuristics for Combinatorial Optimization Problems: Case Study of Academic Course

Scheduling

by Saleh Elmohamed

Abstract of Dissertation November, 1999

Abstract

In this dissertation we try to answer the following question: How do we go about solving optimization problems, in general, and academic course scheduling (or timetabling) at the university level, in particular? we give a detailed answer to the latter part, and shed some light on the case of tackling optimization problems in general.

To summarize, we have investigated a variety of optimization algorithms, in particular, random search, hill climbing, knowledge-based systems (or expert systems), tabu search and simulated annealing-based methods to solve the course scheduling problem. The simulated annealing-based approach includes mean-field annealing, simulated annealing with three different cooling schedules, and the use of a rule-based preprocessor to provide a good initial solution for annealing.

by Saleh Elmohamed M.S., Syracuse University

Dissertation

November, 1999

By Saleh Elmohamed

M.S., Syracuse University

DISSERTATION

submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in computer science and engineering in the

Graduate School of Syracuse University November, 1999

Approved Dr. Geoffrey Fox

Date

\begin{matrix}c \\ ◯ \end{matrix} Copyright 2000 Saleh Elmohamed

Abstract

In this dissertation we try to answer the following question: How do we go about solving optimization problems, in general, and academic course scheduling (or timetabling) at the university level, in particular? we give a detailed answer to the latter part, and shed some light on the case of tackling optimization problems in general.

To summarize, we have investigated a variety of optimization algorithms, in particular, random search, hill climbing, knowledge-based systems (or expert systems), tabu search and simulated annealing-based methods to solve the course scheduling problem. The simulated annealing-based approach includes mean-field annealing, simulated annealing with three different cooling schedules, and the use of a rule-based preprocessor to provide a good initial solution for annealing.

Acknowledgements

iv

Contents

Acknowledgements iv

List of Tables ix

List of Figures x

1 Introduction 1 1.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Brief Notes on NP-Completeness . . . . . . . . . . . . . . . . . . . . 3

1.3.1 Characteristic of NP-Complete Problems . . . . . . . . . . . . 3 1.3.2 The Classes P and NP . . . . . . . . . . . . . . . . . . . . . . 4 1.3.3 The Meaning of NP-Completeness . . . . . . . . . . . . . . . . 5

1.4 Approximation Algorithms . . . . . . . . . . . . . . . . . . . . . . . . 6 1.5 Type of problems addressed . . . . . . . . . . . . . . . . . . . . . . . 7 1.6 Combinatorial Optimization . . . . . . . . . . . . . . . . . . . . . . . 8 1.7 Techniques for Constraint-Directed Scheduling . . . . . . . . . . . . . 9

1.7.1 Heuristic Commitment Techniques . . . . . . . . . . . . . . . 9 1.7.2 Intuitions on Heuristic Performance . . . . . . . . . . . . . . . 10

1.8 Local Search . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 1.8.1 The Reasoning Behind Local Search . . . . . . . . . . . . . . . 12 1.8.2 Examples of Neighborhoods . . . . . . . . . . . . . . . . . . . 13 1.8.3 Considerations to keep in mind . . . . . . . . . . . . . . . . . 15 1.8.4 Local Search Variations . . . . . . . . . . . . . . . . . . . . . 15

2 More Details on Heuristics 22

3 The Timetabling Problem 23 3.1 As a Generalized Assignment-Type Problem . . . . . . . . . . . . . . 23 3.2 Timetabling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

3.2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2.2 Problem representation as GATP . . . . . . . . . . . . . . . . 26 3.2.3 The Problem Constraints . . . . . . . . . . . . . . . . . . . . . 29 3.2.4 Cost Function . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

4 The Rule-Based (Expert) Systems 34 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.2 Examples of Rule-based Constraints . . . . . . . . . . . . . . . . . . . 36

4.2.1 Time and Space . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.2.2 Type of Room . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.2.3 Room Priority . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 4.2.4 The Rationale of Using a Preprocessor . . . . . . . . . . . . . 39

5 Simulated Annealing 40 5.1 Timetabling Using the Annealing Algorithm . . . . . . . . . . . . . . 42 5.2 Interpretation of the Solution Space . . . . . . . . . . . . . . . . . . . 43 5.3 The Annealing Schedules . . . . . . . . . . . . . . . . . . . . . . . . . 44

5.3.1 Cooling Parameters . . . . . . . . . . . . . . . . . . . . . . . . 45 5.3.2 Dependence on Initial Temperature . . . . . . . . . . . . . . . 46

5.4 Our Experimental Schedules . . . . . . . . . . . . . . . . . . . . . . . 46

6 Multi-Phase Approach to Scheduling 50 6.1 Our Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 6.2 Experimental Results For Class Scheduling . . . . . . . . . . . . . . . 50

6.2.1 The Choice of Moves . . . . . . . . . . . . . . . . . . . . . . . 50 6.3 More Analysis of the Simulation . . . . . . . . . . . . . . . . . . . . . 52 6.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

7 Optimization Networks 60

7.1 The Hopfield Network . . . . . . . . . . . . . . . . . . . . . . . . . . 60 7.1.1 Some Background . . . . . . . . . . . . . . . . . . . . . . . . . 60 7.1.2 Structure and Dynamics of the Net . . . . . . . . . . . . . . . 62 7.1.3 Operation of the Network . . . . . . . . . . . . . . . . . . . . 62 7.1.4 Updating of the States of the Neurons . . . . . . . . . . . . . 64 7.1.5 Stability of the Net . . . . . . . . . . . . . . . . . . . . . . . . 64 7.1.6 Updating the States of the Net . . . . . . . . . . . . . . . . . 64 7.1.7 Mapping the Problem onto the Network . . . . . . . . . . . . 65 7.1.8 Shortcomings of the Hopfield Model . . . . . . . . . . . . . . . 65

7.2 Mean Field Annealing for Timetabling . . . . . . . . . . . . . . . . . 65 7.2.1 Introduction and Motivation . . . . . . . . . . . . . . . . . . . 65 7.2.2 Brief Overview . . . . . . . . . . . . . . . . . . . . . . . . . . 68 7.2.3 The Mean-Field Annealing Algorithm . . . . . . . . . . . . . . 69 7.2.4 Observations on the Algorithm . . . . . . . . . . . . . . . . . 70 7.2.5 Definitions and Constraints . . . . . . . . . . . . . . . . . . . 71 7.2.6 More Definitions and Constraints . . . . . . . . . . . . . . . . 80 7.2.7 General Comments on the Algorithm . . . . . . . . . . . . . . 86 7.2.8 Some Implementation Details . . . . . . . . . . . . . . . . . . 86 7.2.9 Performance Issues of MFA on Timetabling . . . . . . . . . . 88 7.2.10 Mapping Timetabling onto the Potts Neural Nets . . . . . . . 88

8 Scheduling as a Graph Coloring Problem 92 8.1 introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 8.2 Mapping Scheduling onto GC . . . . . . . . . . . . . . . . . . . . . . 93 8.3 Graph Coloring Algorithm . . . . . . . . . . . . . . . . . . . . . . . . 95 8.4 Some Experimental Results . . . . . . . . . . . . . . . . . . . . . . . 98 9 Method Comparisons and Complexity Measures 99

9.1 Complexity of Course Scheduling . . . . . . . . . . . . . . . . . . . . 99 9.2 Airline Crew Scheduling . . . . . . . . . . . . . . . . . . . . . . . . . 99

9.2.1 Crew Scheduling vs. Course Scheduling . . . . . . . . . . . . . 101

10 Glossary and Definitions 103

Bibliography 105

List of Tables

6.1 Size of the data set for each of the three semesters. . . . . . . . . . . 56 6.2 The sparseness ratios of the problem for the data sets for each of the

three semesters. Lower values indicate a harder problem. . . . . . . . 57 6.3 Percentage of classes scheduled using the different methods. The aver-

ages and highest and lowest values were obtained using 10 independent runs for simulated annealing (SA) and mean-field annealing (MFA). The expert system (ES) is deterministic so the results are from a single run. No preprocessor was used with the three methods. . . . . . . . 58

6.4 Percentage of scheduled classes, averaged over 10 runs of the same initial temperature and other parameters, for three terms using simulated annealing with an expert system as preprocessor. . . . . . . . . . . . . 59

List of Figures

1.1 A very nice state space. . . . . . . . . . . . . . . . . . . . . . . . . . 13 1.2 A state space with many local optima; it will be hard to find the best

solution. The global optimum (x^{∗}) is at the dark circle, while open circles represent local optima. . . . . . . . . . . . . . . . . . . . . . . 14

1.3 Hill Climbing Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . 16 1.4 Tabu search paradigm . . . . . . . . . . . . . . . . . . . . . . . . . . 18 1.5 Genetic Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 1.6 Random Search Algorithm . . . . . . . . . . . . . . . . . . . . . . . . 21 5.1 The Simulated Annealing Algorithm . . . . . . . . . . . . . . . . . . 41 7.1 A schematic view of the Hopfield network. . . . . . . . . . . . . . . . 63 7.2 The Generic Mean-Field Annealing Algorithm . . . . . . . . . . . . . 70 7.3 A one-way arrow indicates that the “bottom” entity depends on the

“top” one. A two-way arrow indicates a mutual dependency between the two entities. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

8.1 A proper graph coloring using three colors. . . . . . . . . . . . . . . . 94 8.2 A graph representing temporal hard constraints. . . . . . . . . . . . . 95 8.3 DSatur Graph Coloring Algorithm – version 1 . . . . . . . . . . . . . 97 8.4 DSatur Graph Coloring Algorithm – version 2 . . . . . . . . . . . . . 97

Chapter 1

Introduction

1.1 Summary

In this dissertation we try to answer the following question: How do we go about solving optimization problems, in general, and academic course scheduling (or timetabling) at the university level, in particular? we give a detailed answer to the latter part, and shed some light on the case of tackling optimization problems in general.

To summarize, we have investigated a variety of optimization algorithms, in particular, random search, hill climbing, knowledge-based systems (or expert systems), and simulated annealing-based methods to solve the course scheduling problem. The simulated annealing-based approach includes mean-field annealing, simulated annealing with three different cooling schedules, and the use of a rule-based preprocessor to provide a good initial solution for annealing.

The best results were obtained using simulated annealing with adaptive cooling and reheating as a function of cost, and a rule-based preprocessor. This approach enabled us to obtain valid schedules for the timetabling problem for a large university, using a complex cost function that includes student preferences. None of the other methods were able to provide a complete valid schedule.

1.2 Objectives

The primary objective of this dissertation is to shed light on how we go about solving optimization problems in general, and to present detailed analysis on deriving an approximate solution to the problem of university course scheduling in particular.

In order to solve an optimization problem such as the TSP or scheduling , we design a series of steps, or a procedure, which we can tackle a problem instance to which it claims to adhere and, ideally, be sure of getting a good answer in a given amount of time^{∗}. Such a procedure is known as an ‘optimization algorithm’.

Sometimes we may be happy with (or forced to accept) a less than optimal solution, in which case it may be easier to find a suitable algorithm. We would like our algorithms to return as good an answer as possible in as little time as possible, and these will be important comparative judges of the merit of different algorithms. Again, this is more subtle than it seems in that the time will vary from problem instance to another, and we have to decide whether we want to minimize the best case, the worst case or the average case. The reason NP-complete problems are considered difficult is that all known algorithms which guarantee the optimal solution will operate in exponential time in the size of the problem in the worst case scenario. For reasonably large problem instances this is impractical and so we reduce our expectations of the solution from ’optimal’ to ’good’, and seek algorithms which can still search effectively in these difficult search spaces. As well as average versus best/worst case in time-requirements, we also want to understand average versus best/worst case in quality of solution found.

The optimization problem we have closely tackled is university course scheduling, which can be summarized as follows:

Given data sets of classes and their days, enrollments, and instructors; rooms and their capacities, types, and locations; distances between buildings; priorities of each building for different departments; and students and their class preferences; the problem is to construct a feasible class schedule satisfying all the hard constraints and minimizing the medium and soft constraints. Hard constraints are space and time constraints that must be satisfied, such as scheduling only one class at a time for any teacher, student, or classroom. Medium and soft constraints are student and teacher preferences that should be satisfied if possible.

∗Usually we can avoid complications due to the Halting problem, as we address problems of a simpler nature in which we are really fighting against NP-completeness and are sure that a solution does exist.

1.3 Brief Notes on NP-Completeness

1.3.1 Characteristic of NP-Complete Problems

If we pick three different and known NP-complete problems, let say Satisfiability, Knapsack, and Clique, we would notice the following characteristics about them.

• Each is solvable, and a relatively simple approach solves it (although the approach may be time-consuming). For each of them, we can simply enumerate all the possibilities: all ways of assigning the logical values of n variables, all subsets of the set S, all subsets of n vertices in G. If there is a solution, it will appear in the enumeration of all possibilities; if there is no solution, testing all possibilities will demonstrate that.

• There are 2^{n} cases to consider if we use the approach of enumerating all possibilities (where n depends on the problem). Each possibility can be tested in a relatively small amount of time, so the time to test all possibilities and answer yes or no is proportional to 2^{n}.

• They come from different fields; logic, number theory, and graph theory, respectively.

• If it were possible to guess perfectly, we could solve each problem in relatively little time. For example, if someone could guess the correct assignment or the correct subset, we could simply verify that the formula had been satisfied or a correct sum had been determined, or a clique had been identified, etc. The verification process could be done in time bounded by a polynomial of the size of the problem.

1.3.2 The Classes P and NP

Let P be the collection of all problems for which there is a solution that runs in time bounded by a polynomial function of the size of the problem. For example, you can determine whether an item is in a list in time proportional to the size of the list (simply by examining each element in the list to determine whether it is the correct one), and you can sort all items in a list into ascending order in time bounded by the square of the number of elements in the list (using, for example, the well-known bubble sort algorithm.) There may also be faster solutions; that is not important here. Both the searching problem and the sorting problem are in P because they can be solved in time n and n^{2}, respectively.

For most problems, polynomial time algorithms are about the limit of feasible complexity. Any problem that could be solved in time n^{1,000,000,000} would be in P, even though for large values of n, the time to perform such an algorithm might be prohibitive. Notice also that we do not have to know an explicit algorithm, we just have to be able to say that such an algorithm exists.

By contrast, let NP be the set of all problems that can be solved in time bounded by a polynomial function of the size of the problem, assuming the ability to guess perfectly. (in the literature, this “guess function” is called an oracle machine or a nondeterministic Turing machine.) This guessing is called nondeterminism.

Of course, no one can guess perfectly. Guessing is simulated by cloning an algorithm and applying one version of it to each possible outcome of the guess. Essentially, the idea is equivalent to a computer programming language in which IF statements could be replaced by GUESS statements: instead of testing a known condition and branching depending on the outcome of the test, the GUESS statements would cause the program to fork, following two or more paths concurrently.

The ability to guess can be useful. For example, instead of deciding whether to assign the value TRUE or FALSE to variable v_{1}, the nondeterministic algorithm can proceed in two directions: one assuming TRUE had been assigned to v_{1}, and the other assuming FALSE. As the number of variables increases so does the number of possible paths to be pursued concurrently.

Certainly every problem in P is also in NP because the guess function does not have to be invoked. There is also a class EXP, which consists of problems for which a deterministic solution exists in exponential time, c^{n} for some constant c. As noted earlier, every NP-complete problem has such a solution. Every problem in NP is also in EXP so P ⊆ NP ⊆ EXP.

1.3.3 The Meaning of NP-Completeness

The notion of NP-completeness tries to capture the “hardest” problems in the class NP. In other words, the way the class of NP-completeness is defined, if ever we were able to find a polynomial time algorithm for Q, an NP-complete problem, we can then find a polynomial time algorithm for any problem X that is in the class NP (and hence for every NP-complete problem). This follows by the definition of NP-completeness.

Stephen Cook (give reference[cook71]) showed that the satisfiability problem is NP-complete, meaning that it can represent the entire class NP. His important conclusion was that if there is a deterministic, polynomial time algorithm (one without guesses) for the satisfiability problem, then there is a deterministic, polynomial time algorithm for every problem in NP; that is, P = NP.

Richard Karp (give reference[karp72]) extended Cook’s result by identifying several other problems, all of which shared the property that if any one of them could be solved in a deterministic manner in polynomial time, then all of them could. The knapsack and clique problems were identified by Karp. The results of Cook and Karp included the converse: if for even one of these problems (or any NP-complete problem) it could be shown that there was no deterministic algorithm that ran in polynomial time, then no deterministic algorithm could exist for any of them.

Be careful to distinguish between a problem and an instance of a problem. An instance is a specific case: one formula, one specific graph, or one particular set S. Certain simple graphs or simple formulas may have solutions that are very easy and fast to identify. A problem is more general; it is the description of all instances of a given type. For example, the formal statements of the satisfiability, knapsack, and clique sections are statements of problems because they tell what each specific instance of that problem must look like. Solving a problem requires finding one general algorithm that will solve every instance of that problem.

There are problems known to be solvable deterministically in polynomial time (P), and there are problems known not to have a polynomial time solution (EXP and beyond), so that P ⊆ EXP and P ≠ EXP, meaning P ⊂ EXP. The class NP fits somewhere between P and EXP: P ⊆ NP ⊂ EXP. It may be that P = NP, or that P ≠ NP.

The significance of Cook’s result is that NP-complete problems have been studied for a long time by many different groups of people: logicians, operations research specialists, electrical engineers, number theorists, operating systems specialists, and communications engineers. If there were a practical (polynomial time) solution to any one of these problems, you would hope that someone would have found it by now. See Garey and Johnson’s N-completeness catalog.

1.4 Approximation Algorithms

When facing an NP-complete optimization, we may want to consider algorithms that do not produce optimum solutions, but solutions guaranteed to be close to the optimum. Suppose that we wish to obtain such solutions for an optimization problem, maximization or minimization. For each instance x of this problem, there is an optimum solution with value opt(x); let us assume that opt(x) is always a positive integer.

Suppose now that we have a polynomial algorithm A which, when presented with instance x of the optimization problem, returns some solution with value A(x). Since the problem is NP-complete and A is polynomial, we cannot realistically hope that A(x) is always the optimum value. But suppose that we know that the following inequality hold:

\frac{|opt(x) − A(x)|}{opt(x)} ≤ ϵ  (1.1)

where ϵ is some positive real number, hopefully very small, that bounds from above the worst-case relative error of algorithm A. (The absolute value in this inequality allows us to treat both minimization and maximization problems within the same framework.) If algorithm A satisfies this inequality for all instances x of the problem, then it is called an ϵ-approximation algorithm.

Once an optimization problem has been shown to be NP-complete, the following question becomes most important: Are there ϵ-approximation algorithms for this problem? And if so, how small can ϵ be? Let us observe at the outset that such questions are meaningful only if we assume that p ≠ NP, because, if P = NP, then the problem can be solved exactly, with ϵ = 0.

1.5 Type of problems addressed

A question to ask here is that: What problems can stochastic optimization algorithms such as simulated annealing address?

Most of the algorithms we have looked at adhere to a simple set of requirements upon the space of problems to which they apply. In fact we require no more of a problem than that:

1. we can generate a random solution,

2. a ’local perturbation’ to a solution is defined,

3. such perturbations over all relevant parts of the solution space (this is usually naturally satisfied).

All of the algorithms discussed in this thesis fall within this category. Also note that because we search directly in the problem space, all the algorithms we consider will return both the actual solution and its cost (or value of utility depending upon the context). One can also envisage algorithms which only return a bound upon the solution cost, and not the actual solution. Whilst of theoretical value, in most practical instances we would ordinarily require the solution as well. Finally, it is important to note that:

• individual problem instances may vary enormously in difficulty. The average versus best or worse case difficulty (or complexity) of search and problem solving may depend in quite particular fashions on the method of problem selection. These issues were addressed later in greater depth.

• these algorithms may not specify a natural termination condition, and in such cases heuristic principles must be used. Possibilities include: maximum running time, solution quality, time since last improvement; or combinations thereof.

In the descriptions which follow, we shall assume we are supplied with a function ‘FITNESS: ∑ → R’ which supplies us with the solution quality. Without loss of generality, we assume higher fitnesses are better.

1.6 Combinatorial Optimization

Combinatorial optimization refers to the problem of minimizing a function of a large but finite set of variables. The cost function or energy function E(x⃗) assigns a real number to each possible state x⃗ of the system. The configuration space X is given by the finite set of all possible system configurations {x⃗}. The problem is then that of finding the global minimum of E(x⃗) over all x⃗ in X.

The existence of a large number of local but not global minima of E(x⃗) in configuration space conspires against the success of heuristic optimization methods. Algorithms such as simulated annealing is an improvement upon downhill algorithms in that it provides a mechanism for getting out of such local minima. It should be noted that while the notion of global minima is an absolute one that depends only on the cost function E(x⃗), the concept of local minima is relative to the topology since it implies comparing the value of E(x⃗) with that of E(x′) for those x′ which are neighbors of x⃗.

The topology of configuration space is defined by the neighbor relation. A natural definition of neighborhood for optimization problems is given by the move set: a small rearrangement of the system that produces a new trial configuration from the present one. Distance d(x⃗, x′) between two states is then defined as the minimum number of moves needed to turn one into other; and the neighborhood of a given state x⃗ is the set of all states x′ which can be reached from x⃗ with one move, i.e. d(x⃗, x′) = 1.

1.7 Techniques for Constraint-Directed Schedul-

ing

Christopher Beck, in his PhD thesis (ref: Beck-99), discusses among other things, constraint-directed techniques that have been applied to scheduling. These are heuristic commitment techniques, propagators, and retraction techniques. In the discussion he uses the term commitment to refer to a decision made during search. It is referred to or conceptualized as a constraint that is added to or removed from the constraint graph. For example, assigning a start time to an activity can be modeled as the assertion of a unary equals commitment on the activity, while sequencing two activities can be represented as the assertion of a binary precedence constraint between the two activities.

1.7.1 Heuristic Commitment Techniques

A heuristic commitment technique is a procedure that finds new commitments to be added to the graph in order to (heuristically) move toward a solution. In traditional constraint-directed search (and, indeed, artificial intelligence search more generally) heuristic commitments are the primary component of a search for a solution. A critical component of constraint-directed scheduling techniques, therefore, is the technique used to generate and select the heuristic commitments to make. We will be looking at a number of heuristics that have been applied to constraint-directed scheduling.

1.7.2 Intuitions on Heuristic Performance

Here we present intuitions about the nature of heuristic commitment techniques. These intuitions are speculative and their investigation forms the basis of our future work on heuristics.

We believe that there are two characteristics that correlate highly with the success of a heuristic commitment technique.

1. The Short dead-end Proof Intuition: the ability of a heuristic to quickly discover that the search is at a dead-end. Imagine a search that progresses through a series of states to state S. Assume that there are no solutions in the subtree below S, but that we are unable to prove this at state S. Some heuristic commitments must be made to explore the subtree below S in order to prove that S is a deadend. A heuristic that is able, through its heuristic commitments and subsequent propagation, to prove a dead-end while exploring a small number of states in the subtree will be superior to on that must explore a larger subtree.

2. The Few Mistakes Intuition: a high likelihood of making a commitment that leads to a solution. A heuristic that has a higher likelihood of making a commitment that leads to a solution, that is one that does not result in a dead-end, is likely to be a higher quality heuristic than one with a lower probability of making the correct decision.

These two characteristics are independent from the perspective that if either characteristic were infallible, search performance would be independent of the other characteristic. If a heuristic never made a commitment that resulted in a dead-end state, the ability to detect such a state is irrelevant. Similarly, if a search could immediately detect that it was a dead-end, then the quality of the commitments is irrelevant: any mistake is immediately detected and repaired.

With real, fallible heuristics, these two characteristics interact: it may be the case that a heuristic with a high probability of making the right decision results in a relatively large effort in discovering when it has made a mistake. Conversely, a heuristic with a lower probability of making a right decision may be able to find a shorter proof of a dead-end (on average) and therefore recover from mistakes quickly. The best heuristic will be the one that minimizes the size of the overall search tree and this may well be achieved by some trade-off between the two characteristics.

1.8 Local Search

Local search is meant to represent a large class of similar techniques that can be used to find a good solution from a problem. The idea is to think of the solution space as being represented by an undirected graph. That is, each possible solution is a node in the graph. An edge in the graph represent a possible move we can make between solutions. For an optimization problem such as scheduling, even if we want to, we could never hope of writing down the graph of the problem search space.

Therefore, the simple idea of local search is that we never actually try to write the whole graph down; we just move from one possible solution to a “nearby” possible solution, either for as long as we like, or until we happen to find an optimal solution.

To set up a local search algorithm, we need to have the following:

1. A set of possible solutions, which will be the vertices in our local search graph.

2. A notion of what the neighbors of each vertex in the graph are. For each vertex x, we will call the set of adjacent vertices N(x). The neighbors must satisfy several properties: N(x) must be easy to compute from x (since if we try to move from x we will need to compute the neighbors), if y ∈ N(x) then x ∈ N(y) (so it makes sense to represent neighbors as undirected edges), and N(x) cannot be too big, or more than polynomial in the input size (so that the neighbors of a node are easy to search through).

3. A cost function, from possible solutions to the real numbers. For example, f(s^{∗}) : R → R.

The most basic local search algorithm (say to minimize the cost function) is easily described:

1. Pick a starting point x.

2. While there is a neighbor y of x with f(y) < f(x), move to it; that is, set x to y and continue.

3. Return the final solution.

1.8.1 The Reasoning Behind Local Search

The idea behind local search is clear; if keep getting better and better solutions, we should end up with a good one. Pictorially, if we “project” the state space down to a two dimensional curve, we are hoping that the picture has a sink, or a global optimum, and that we will quickly move toward it. See Figure 1.1.

There are two possible problems with this line of thinking. First, even if the space does look this way, we might not move quickly enough toward the right solution. Also, for any reasonably large scale optimization problem, such as scheduling, if we start with a bad solution (i.e. something with a very large cost), it will take a lot of moves to reach the minimum. Generally, however, this is not much of a problem, as long as the cost function is reasonably simple.

The more important problem is that the solution space might not look at all the way as shown in Figure 1.1. For example, our cost function might not change smoothly when we move from a state to its neighbor. Also, it may be that there are several local optima, in which case our local search algorithm will hone in a local optimum and get stuck. See Figure 1.2.

The second problem, that the solution space might not “look nice”, is crucial, and it underscore the importance of setting up the problem. When we choose the possible

Figure 1.1: A very nice state space.

moves between solutions — that is, when we construct the mapping that gives us the neighborhood of each node — we are setting up how local search will behave, including how the cost function will change between neighbors, and how many local optima there are. How well local search will work depends tremendously on how smart one is in setting up the right neighborhoods, so that the solution space really does look the way we would like it to.

1.8.2 Examples of Neighborhoods

• Course Scheduling: A possible neighborhood structure is two assignments of a set of classes to a set of rooms over a give set of time slots, if the two assignments differ only in one class assignments. That is in the first assignment, class (A) is assigned to room (B) and class (C) is assigned to room (D); while in the second Figure 1.2: A state space with many local optima; it will be hard to find the best solution. The global optimum (x^{∗}) is at the dark circle, while open circles represent local optima.

assignment, class (A) is assigned to room (D) and class (C) to room (B).

• MAX3SAT: A possible neighborhood structure is two truth assignments are neighbors if they differ in only one variable. A more extensive neighborhood could make two truth assignments neighbors if they differ in at most two variables; this trades increased flexibility for increase size in the neighborhood.

• Traveling Salesperson: The K-opt neighborhood of x is given by all tours that differ in at most k edges from x. In practice, using the 3-opt neighborhood seems to perform better than the 2-opt neighborhood, and using 4-opt or larger increases the neighborhood size to a point where it is inefficient.

1.8.3 Considerations to keep in mind

There are several aspects of local search algorithms that we can vary, and all can have an impact on performance. For example:

1. What are the neighborhoods N(x)?

2. How do we choose an initial starting point?

3. How do we choose a neighbor y to move to? (Do we take the first one we find, a random neighbor that improves f, the neighbor that improves f the most, or do we use other criteria?)

4. What if there are ties?

There are other practical things to consider. Can we re-run the algorithm several times? Can we try several of the algorithms on different machines? Issues like these can have a big impact on actual performance. However, perhaps the most important issue is to think of the right neighborhood structure to begin with; if this is right, then other issues are generally secondary, and if this is wrong, the likelihood of success is not that great.

1.8.4 Local Search Variations

There are several approaches that take advantage of search strategies in which costdeteriorating neighbors are accepted. In particular, we mention simulated annealing, tabu search, genetic algorithms, and neural nets (via mean field annealing).

In the descriptions of the local search methods which follow, we will assume we are supplied with a function ‘FITNESS: ∑ → R’ which supplies us with the solution quality. Without loss of generality, we assume higher fitnesses are better.

• Hill-climbing: This is the name for the basic variation, where one moves to a vertex of lower (or possibly equal) cost. See Figure 1.3.

Hill Climbing (Σ)

1. Pick a random solution σ_{best} ∈ Σ

2. repeat

3. Pick a random solution σ ∈ Σ

4. repeat

5. σ’ ← Perturbation(σ)

6. if FITNESS(σ’) > FITNESS(σ)

7. then σ ← σ’

8. until No fitter perturbations of σ exist

9. if FITNESS(σ) > FITNESS(σ_{best})

10. then σ_{best} ← σ

11. until No time remains

12. return σ_{best}

Figure 1.3: Hill Climbing Algorithm

This algorithm is called hill-climbing, because the inner loop (steps 4-8) move uphill (not necessarily in the steepest direction in the basic algorithm) until a summit is reached. The outer loop (steps 2-11) mean that the algorithm tries to find further hills if given enough time, and finally returns the highest peak it has visited when time runs out (this algorithm is sometimes called ‘iterated hill climbing’ because of the existence of the outer loop). To use this algorithm we obviously require that the perturbation function be defined. Whether hill-climbing works well as an optimization algorithm depends upon the local^{†} structure of the landscape — if high points in the landscape are surrounded by a smooth slopes leading away long distances in the space, then any points picked on those slopes will lead to the top of the hill. The set of all such points is known as the ‘attracting set’, or ‘capture basin’ of the peak of the hill. If, on the other hand, tall peaks are surrounded by great crevasses then hill climbing will perform poorly. This image in terms of fitness landscape is extremely important in trying to understand when and why a particular algorithm will perform.

†Local being defined by the perturbation operator, of course.

• Metropolis rule: Pick a random neighbor, and if the cost is lower, move there. If the cost is higher, move there with some probability (that is usually set to depend on the cost differential). The idea is that possibly moving to a worse state helps avoid getting trapped at local minima.

• Simulated annealing: (see other chapters for more details) This method is similar to the Metropolis rule, except that the probability of going to a higher cost neighbor varies with time and is gradually decreased in the course of the algorithm’s execution. Lowering the acceptance probability is controlled by a set of parameters whose values are determined by a cooling schedule. Like most local search algorithms, simulated annealing is normally used as an approximation algorithm, for which much faster convergence rates are acceptable (ref: Aarts and Korst-89, Aarts-97).

• Tabu Search: It combines the deterministic iterative improvement algorithm with a possibility to accept cost-increasing solutions. In this way the search is directed away from local minima, such that other parts of the search space can be explored. The next solution visited is always chosen to be a legal neighbor of the current solution with the best cost, even if the cost is worse than that of the current solution. The set of legal neighbors is restricted by a tabu list designed to prevent us going back to recently visited solutions. The tabu list is dynamically updated during the execution of the algorithm. The tabu list defines solutions that are not acceptable in the next few iterations. However, a solution on the tabu list may be accepted if its quality is in some sense good 1. begin TABU

2. select an initial solution i ∈ S;

3. initialize tabu list T and aspiration function A;

4. while the stopping condition is not met do

5. compute the best i′ in N(i)

6. subject to tabu list T and aspiration function A;

7. i ← i′;

8. update tabu list T;

9. update aspiration function A;

10. end while

11. end TABU

Figure 1.4: Tabu search paradigm

enough, in which case it is said to attain a certain aspiration leve. See Figure 1.4.

• Genetic algorithms: (see Figure 1.5) This search strategy is based on a theme borrowed from Holland’s genetic algorithms approach. (ref: Holland-75, Goldberg-89). Holland uses concepts from population genetics and evolution theory to construct algorithms that try to optimize the fitness of a population of elements through recombination and mutation of their genes. There are many variations known in the literature of algorithms that follow these concepts. One example is genetic local search. The general idea is as follow:

– Step 1, initialize. Construct an initial population of n solutions.

– Step 2, improve. Use local search to replace the n solutions in the population by n local optima.

– Step 3, recombine. Augment the population by adding m offspring solutions; the population size now equals n + m.

Genetic Algorithm (Σ)

1. Pick a population P = {σ_{1}, …, σ_{n}} drawn randomly from Σ

2. repeat

3. Form a selection:

4. S ← Selection(P, Fitness())

5. P ← ∅

6. repeat

7. Pick a pair σ_{a},σ_{b} from S

8. Form the offspring:

9. σ_{ab,1},σ_{ab,2} ← Crossover(σ_{a},σ_{b})

10. P ← P ∪ {σ_{ab,1},σ_{ab,2}}

11. until Population is full: |P | = n

12. until P ≈ converged to a single solution

13. return σ = arg max_{σ∈Σ} {Fitness(σ)}

Figure 1.5: Genetic Algorithm

– Step 4, improve. Use local search to replace the m offspring solutions by m local optima.

– Step 5, select. Reduce the population to its original size by selecting n solutions from the current population.

– Step 6, evolute. Repeat Step 3 through 5 until a stop criterion is satisfied.

Evidently, recombination is an important step, since here one must try to take advantage of the fact that more than local optimum is available, i.e., one most exploit the structure present in the available local optima.

• Artificial Neural networks: Consist of networks of elementary nodes (neurons) that are linked through weighted connections. The nodes represent computational units, which are capable of performing a simple computation, consisting of a summation of the weighted inputs, followed by the addition of a constant called the threshold or bias, and the application of a nonlinear response function. The result of the computation of a unit constitutes its output. This output is used as an input for the nodes to which it is linked through an outgoing connection. The overall task of the network is to achieve a certain network configuration, for instance a required input-output relation, by means of the collective computation of the nodes. This process is often called self-organization. Over the years a large variety of neural network models have been proposed, which differ in the net architecture and the computational model used for selforganization. In many models, the self-organization of the network tries to find a stable network configuration. For models such as the Hopfield network or the Boltzmann machine, the self-organization can be modeled in terms of local search. Here the issue of escaping from local optima is again of prime importance. In a a number of cases it is resolved by applying probabilistic methods, as in simulated annealing. This takes to the mean field annealing model outlined in another chapter.

• Random Search: If we don’t have enough time to search all possible solutions, we could just generate solutions randomly for the length of time available, and retain the best one as our result. We formalize this as follows(see Figure 1.6): It is quite clear that unless good solutions are pretty dense in the search space, a random search will almost always not perform well. However, another use of random search is as a benchmark — a standard level from which to compare other algorithms and with which to produce a very simple measure of how difficult a problem is (such comparisons will never be ideal, however, for at least the reason that different computer architectures and compilers may have different levels of performance across different types of computer code). Choices for termination condition in the random search include time or solution quality. However since we cannot expect solutions of very good quality, the only practical condition is on the allowable search time.

Random Search (Σ)

1. Pick a random solution σ ∈ Σ

2. repeat

3. Pick a random solution σ’ ∈ Σ

4. If FITNESS(σ’) > FITNESS(σ)

5. then σ ← σ’

6. until No time remains

7. return σ

Figure 1.6: Random Search Algorithm

It is not at all clear on which problem classes these algorithms are expected to perform well. If we are presented with a new problem demanding a solution we would like to be able to know a priori which of these algorithms is expected to perform best. To do this, we need a formal interpretation of their dynamics, motion in search space and convergence properties. From these we can derive a set of tests or summary of inherent advantages and disadvantages which will let to pick the best algorithm for a given purpose. Given that our search problem is difficult, it is unlikely that a first guess approach will yield sufficiently good answers (for whatever purpose we had in mind), so we cannot be happy with that status quo. Somehow the structure of the search space (the almost infinite set of possible perturbation operators and the generated fitness landscapes) is is vital in dictating how the space should be searched efficiently.

In more direct terms, we must think about a process for generating or transforming solutions in a manner which we intuitively feel takes into account the deep structure of the problem. In short, without an understanding of the manner in which optimization algorithms operate, for every new case we need to invent a new heuristic, or rule of thumb, which seems particularly relevant. Clearly, this poses great difficulties for the unfortunate algorithm designer.

Chapter 2 More Details on Heuristics The Timetabling Problem

3.1 As a Generalized Assignment-Type Problem

The notion of the assignment-type problem and its generalized version appear to be powerful modeling tools. In particular, such models are very appropriate for formulating timetabling and scheduling problems where time periods have to be determined for activities according to several specific constraints.

The Assignment-Type Problem (ATP) can be summarized as follows:

Given n items and m resources, the problem is to determine an assignment of each item to a resource in order to optimize an objective function and to satisfy K additional side constraints.

The mathematical model associated with an ATP is formulated as follows:

(ATP)

Min F (x)  (3.1)

subject to

G_{k}(x) ≤ 0 1 ≤ k ≤ K  (3.2)

and

x ∈ X(1)  (3.3)

where

X(1) = {x : ∑_{j∈J_{i}} x_{ij} = 1, 1 ≤ i ≤ n; xij = 0 or 1, 1 ≤ i ≤ b, j ∈ J_{i}}  (3.4)

and J_{i} ⊂ {1, 2, …, m} is the set of admissible resources for item i, 1 ≤ i ≤ n. The decision variables x_{ij} are such that:

x_{ij} = \{ \begin{matrix}1 & if item i assigned to resource j \\ 0 & otherwise. \end{matrix}  (3.5)

The objective function F and the side constraints G_{k}, 1 ≤ k ≤ K, are only required to be calculable.

A Generalized Assignment-Type Problem (GATP) is an ATP where each item i has to be assigned to a_{i} resources, a_{i} ≥ 1. The mathematical model associated with a GATP is:

(GATP)

Min F (x)  (3.6)

subject to

G_{k}(x) ≤ 0 1 ≤ k ≤ K  (3.7)

and

x ∈ X(a)  (3.8)

where

X(a) = {x : ∑_{j∈J_{i}} x_{ij} = a_{i}, 1 ≤ i ≤ n; xij = 0 or 1, 1 ≤ i ≤ n, j ∈ J_{i}}  (3.9)

and

a = [a_{1}, a_{2}, …, a_{n}]^{T}  (3.10)

and a_{1} ≥ 1 and integer 1 ≤ i ≤ n. It is easy to see that an ATP is a GATP where a = [1, 1, …, 1]^{T}.

The constraints x ∈ X(a) underly several timetabling and scheduling problems. Furthermore, to take advantage of the simple structure of the constraints x ∈ X(a), a penalty approach is used to deal with GATP. In this approach, let V_{k}(x) denote the violation of constraint G_{k}(x) ≤ 0; i.e.

V_{k}(x) = Max {0, G_{k}(x)}.  (3.11)

Hence, for each x ∈ X(a), let V (x) denote the total weighted violation:

V (x) = ∑_{k=1}^{K} β_{k}V_{k}(x)  (3.12)

where β_{k} ≥ 0 is proportional to the relative importance of the kth constraint. Then, instead of solving GATP, the solution techniques already introduced, are applied to the following penalized version of GATP:

(PGATP)

Min_{x∈X(a)}P (x) = αF (x) + V (x)  (3.13)

where α ≥ 0. The values of α and β_{k}, 1 ≤ k ≤ K, are parameters that can be adjusted according to the specific application. For instance, larger values for β_{k}, 1 ≤ k ≤ K, allow putting more emphasis on feasibility.

3.2 Timetabling

3.2.1 Introduction

Timetabling is the assignment of time slots to a set of events, subject to constraints on these assignments. The NP-complete professors and classes timetabling problem is a constraint satisfaction problem that can be briefly stated as follows:

For a certain school with N_{p} professors, N_{q} classes, N_{x} classrooms and lecture halls, and N_{s} students, it is required to schedule N_{l} professor-class pairs within a time limit of N_{t} time slots producing a legal schedule. A legal schedule needs to be found such that no professor, class, or student is in more than one place at a time, and no room is expected to accommodate more than one lesson at a time or more students than its capacity.

The constraints for this problem can be hard, medium or soft. The medium and soft constraints have an associated cost (or penalty), and if they are not satisfied, the goal is to minimize this cost. Soft constraints have a lower priority (and thus lower cost) than medium constraints. The hard constraints must be satisfied, so their associated cost must be reduced to zero. A feasible schedule is one that satisfies all the hard constraints.

3.2.2 Problem representation as GATP

Consider the problem of establishing the schedule of lectures accounting for individual student registrations, lecturers, and classroom availabilities. Note that this formulation allows for lectures of different length. Hence, to formulate this problem as a GATP, the lectures are the items and the starting times allowed are the resources; then for each lecture i, 1 ≤ i ≤ n, and starting time j ∈ J_{i} (the set of admissible starting times for lecture i) such that J_{i} = {1, 2, …, m}:

x_{ij} = \{ \begin{matrix}1 & if lecture i starts at j \\ 0 & otherwise. \end{matrix}  (3.14)

The objective function accounts for the lecturer preferences:

F (x) = ∑_{i=1}^{n} ∑_{j∈J_{i}} c_{ij}x_{ij}  (3.15)

where c_{ij} is the cost of starting lecture i at j specified in terms of lecturer preferences (c_{ij} = 0 if j is the most preferred starting time for lecture i).

A first set of side constraints is specified to eliminate the conflicting situations where students or lecturers are involved in lectures taking place simultaneously. Define:

Γ_{ij} = {(k, l): lectures i and k have students or lecturers in common, and they overlap in time if lecture i starts at j and lecture k at l }.

The conflicting constraints are specified as follows:

x_{ij}x_{kl} ≤ 0 (k, l) ∈ Γ_{ij} , 1 ≤ i ≤ n , j ∈ J_{i} .  (3.16)

Other side constraints are introduced to account for classroom availability. Partition the set of classrooms available into several subsets including classrooms of the same type. Define:

B = number of classroom types U_{b} = number of classrooms of type b , 1 ≤ b ≤ B. For each lecture, a unique classroom type is specified. Define: K_{bt} = {(i, j): lecture i requires a classroom of type b and it takes place during

teaching hour t of the week whenever its starting time is j }.

The classroom availability constraints are as follow:

∑_{(i,j)∈K_{bt}} x_{ij} − U_{b} ≤ 0 1 ≤ b ≤ B , 1 ≤ t ≤ T  (3.17)

where T is the total number of teaching hours of the week. The course timetabling problem can be summarized as follows:

Min ∑_{i=1}^{n} ∑_{j∈J_{i}} c_{ij}x_{ij} .  (3.18)

Subject to:

x_{ij}x_{kl} ≤ 0 (k, l) ∈ Γ_{ij}, 1 ≤ i ≤ n , j ∈ J_{i}  (3.19)

Next, we need to derive the associated penalized version of this problem (i.e. the penalized GATP), with weights specified to indicate the relative importance of violating the corresponding constraints. For the conflicting constraints, the weights account for the length of the conflict and the number of individuals involved. Define:

• L_{ijkl} = length (in number of hours) of the conflict (k, l) ∈ Γ_{ij} , 1 ≤ i ≤ n , j ∈ J_{i} .

• δ_{ik} = number of students taking both lectures i and k, 1 ≤ i ≤ n , 1 ≤ k ≤ n .

• γ_{ik} = \{ \begin{matrix}M + δ_{ik} & if lectures i and k have lecturers in common \\ δ_{ik} & otherwise \end{matrix}

For 1 ≤ i ≤ n, 1 ≤ k ≤ n .

where M is a parameter to indicate the relative importance of conflicts due to lecturers. Let V_{kl}^{ij}(x) denote the violation of the conflicting constraints:

V_{kl}^{ij}(x) = Max{0, x_{ij}x_{kl}} = x_{ij}x_{kl} .  (3.20)

Let V_{bt}(x) denote the violation of the classroom availability constraints:

V_{bt}(x) = Max{0, ∑_{(i,j)∈K_{bt}} x_{ij} − U_{b}} .  (3.21)

Furthermore, let the weight β_{kl}^{ij} associated with V_{kl}^{ij}(x) be specified as follows:

β_{kl}^{ij} = γ_{ik}L_{ijkl} .  (3.22)

Also, let the weight β_{bt} associated with V_{bt}(x) be equal to a scalar ρ for all 1 ≤ b ≤ B , 1 ≤ t ≤ T . Then, the penalized GATP for this problem is as follows:

Min_{x∈X(1)} ∑_{i=1}^{n} ∑_{j∈J_{i}} [αc_{ij}x_{ij} + ∑_{(k,l)∈Γ_{ij}} γ_{ik}L_{ijkl}x_{ij}x_{kl}] + ρ ∑_{b=1}^{B} ∑_{t=1}^{T} Max {0, ∑_{(i,j)∈K_{bt}} x_{ij} − U_{b}} .

(3.23)

The weight parameters α, ρ, and M are selected according to the specific context. For instance, if the priority is to reduce the conflicts, then α and ρ take smaller values.

3.2.3 The Problem Constraints

Hard constraints are usually constraints that physically cannot be violated. This includes events that must not overlap in time, such as:

• classes taught by the same professor,

• classes held in the same room,

• a class and a recitation or a lab of the same class.

Another examples are space or room constraints:

• Room Type: Some classes, such as laboratories, require a certain type of room.

• Room Size: A class cannot be assigned to a particular room unless the capacity of the room (denoted by R_{c}) is greater than or equal to the class enrollment (denoted by C_{e}). The smaller the difference between these two sizes, the better the assignment.

There is a cost associated with any violation of these constraints. For the Room Size constraint our cost function C_{rs} computes the magnitude of the difference between the room capacity R_{c} and the class enrollment C_{e} as follows:

C_{rs} = \{ \begin{matrix}(C_{e} − R_{c})^{2} & if C_{e} > R_{c} , \\ tanh((1.0 − \frac{C_{e}}{R_{c}})βγ)(R_{c} − C_{e}) & otherwise. \end{matrix}  (3.24)

In our experiments, the parameters β and γ were set to 0.4 and 0.8, respectively. When the class enrollment is larger than the room size, this function will always yield a much higher cost than the cost incurred from the assignment of a class with an enrollment that is smaller than the room size. The use of tanh function is to nonlinearize the cost and not to make its increase directly proportional to the difference between room size and class enrollment.

There is also a cost assigned to any violation to the constraint that certain classes cannot be scheduled at the same time or have any overlap in the scheduled times. This is called exclusion cost C_{x} and derived as follows:

C_{x} = \{ \begin{matrix}(end(c_{i}) − s(c_{j}))Rlen(c_{i}, c_{j}) & if (s(c_{j}) > s(c_{i})) and (end(c_{i}) ≥ s(c_{j})) , \\ (end(c_{j}) − s(c_{i}))Rlen(c_{i}, c_{j}) & if (s(c_{j}) < s(c_{i})) and (end(c_{j}) ≥ s(c_{i})) , \\ (end(c_{j}) − s(c_{i}))Rlen(c_{i}, c_{j}) & if (s(c_{j}) ≥ s(c_{i})) and (end(c_{i}) ≥ end(c_{j})) , \\ (end(c_{i}) − s(c_{j}))Rlen(c_{i}, c_{j}) & if (s(c_{i}) ≥ s(c_{j})) and (end(c_{j}) ≥ end(c_{i})) . \end{matrix}

(3.25)

where the relation between the two class durations is defined as follows:

Rlen(c_{i}, c_{j}) = \{ \begin{matrix}len(c_{i})/len(c_{j}) & if len(c_{i}) ≥ len(c_{j}) , \\ len(c_{j})/len(c_{i}) & otherwise. \end{matrix}  (3.26)

The first and second choices of C_{x} take care of any partial overlaps between c_{i} and c_{j}, while the third and fourth choices take care of any complete overlaps between the two classes. Medium constraints are usually considered to be those constraints that fall into the gray area between the hard and soft constraints [?] (ref: Eiselt-Laporte-87). In our implementation, we define medium constraints to be constraints such as time and space conflicts which, like hard constraints, cannot physically be violated (for example, it is not possible for one person to be in two different classes at the same time). However we consider these constraints to be medium rather than hard if they can be avoided by making adjustments to the specification of the problem. The primary example is student preferences. We cannot expect to be able to satisfy all student class preferences, in some cases, certain students will have to adjust their preferences since certain classes will clash, or will be oversubscribed.

Medium constraints have a high penalty attached to them, although not as high as that associated with the hard constraints. In the final schedule the penalty of these constraints should be minimized and preferably reduced to zero. Some examples of medium constraints are:

• Avoid time conflicts for classes with students in common.

• Eligibility criteria for the class must be met.

• Do not enroll athletes in classes that conflict with their sport practice time (of course, depending on the sport).

• The predecessor-successor relation between the topics is satisfied.

• Two courses of a same topic are scheduled on two different working days. Soft constraints are preferences that do not deal with time conflicts, and have a lower penalty (or cost) associated with them. We aim to minimize the cost, but do not expect to be able to reduce it to zero. Some examples are:

• For each student, balance the three-day (Mon, Wed, Fri) as well as the two-day (Tue, Thu) schedules.

• Balance or spread out the lectures over the week.

• Classes may request contiguous time slots.

• Balance enrollment in multi-section classes.

• Lunch and other break times may be specified.

• Professors may request periods in which their classes are not taught.

• Professors may have preferences for specific rooms or types of rooms.

• Once a class enrollment exceeds the upper limit, the class is split into multisections.

• Minimize the distance between the room where the class is assigned and the building housing its home department.

Some soft constraints may have higher priority (and thus higher cost) than others. For example, preferences involving professors will have higher priority than the preferences of students.

We have dealt with the distance minimization constraint in the following way. Given the various building preferences of the departments, we have constructed a matrix M_{IJ} between all the academic departments I and all the buildings involved in scheduling J. Using M_{IJ} in conjunction with the (appropriately scaled) distance between all buildings involved in the process, a final distance matrix Q is derived and directly used in the scheduling process. Also, unless otherwise stated, a department’s home building is always the first preference for the department classes to be assigned.

Let B denote the set of all k buildings, D denote the set of all n departments. Also for department d_{i} where 1 ≤ i ≤ n the distance to all buildings B is a vector denoted by < d_{i}, b_{j} > where b_{j} ∈ B and 1 ≤ j ≤ k. At each step of the scheduling process and according to d_{i} space preferences, d_{i} courses would be scheduled into buildings B′ where |B′| ≤ |B| and the distance is denoted by the vector ≪ d_{i}, b_{j} ≫ where b_{j} ∈ B′. In addition, let |B′| = k′.

Now the cost associated with the distance of department d_{i} is computed as follows:

C_{d_{i}}^{dist} = \frac{∑_{j=1}^{k′} ≪ d_{i}, b_{j} ≫}{∑_{j=1}^{k} < d_{i}, b_{j} >} k′  (3.27)

The cost function is the distance ratio of department(i)’s buildings of preference to all available buildings multiplied by the number of department(i)’s preferred buildings. Notice that for a particular department the denominator of equation 3.27 stays fixed the scheduling process.

Other constraints that may need to be taken into account during the scheduling process:

• It is not the case, in general, that all classrooms are available for all subjects. In reality, many subjects require special purpose rooms.

• Some rooms may be declared unavailable during any time period.

• Daily limits must be specified to restrict the assignment of more than a specified number of classes per day to each professor.

. . . . . .

3.2.4 Cost Function

The cost function measures the quality of the current schedule and generally involves the weighted sum of penalties associated with different types of constraint violations. The aim of the optimization technique is to minimize the cost function. We consider the following as part of the cost function:

• Class cost is the sum of costs incurred from scheduling a class in a particular time period, including having a certain class and its lab scheduled in the same time period; spreading classes over the week and within each of the weekdays; and gluing to fill in holes in the timetable and allow contiguous time periods for a class that requests them.

• Professor cost is the sum of costs associated with each professor, including teaching two or more classes that overlap in time; being assigned more classes than the maximum number allowed to be taught; and not having a long enough time interval between classes taught by the same professor.

• Student cost is the sum of costs associated with each student, including clashes in class times; not having the classes evenly spread out over the week; having to schedule a class from the list of alternates rather than the list of first preferences; or classes that do not satisfy eligibility criteria.

• Room cost is the sum of costs associated with each classroom, including the distance cost obtained from the distance constraint; not being able to schedule a room of the specified type; and the room size cost obtained from the room constraint.

The Rule-Based (Expert) Systems

4.1 Introduction

We have implemented a fairly complex rule-based expert system for solving the timetabling problem, for three reasons. Firstly, it gives us a benchmark as to how well other methods do in comparison to this standard technique. Secondly, a simplified version of the rule-based system is used to provide sensible choices for moves in the simulated annealing algorithm, rather than choosing swaps completely at random, and this greatly improves the proportion of moves that are accepted. Thirdly, we have used this system as a preprocessor for simulated annealing, in order to provide a good initial solution.

Simulated annealing is a very time-consuming, computationally intensive procedure. Using an expert system as a preprocessor is a way of quickly providing a good starting point for the annealing algorithm, which reduces the time taken in the annealing procedure, and improves the quality of the result. Our results clearly support this rationale for the case of academic scheduling.

The rule-based expert system consists of a number of rules (or heuristics) and conventional recursion to assist in carrying out class assignments. We have developed this system specifically for the problem of academic scheduling. The basic data structures or components of the system are:

1. Distance matrix of values between each academic department and every other building under use for scheduling.

2. Class data structure of each class scheduled anywhere in campus. These structures are capable of linking with each other.

3. Room data structure of each room (regardless of type) involved in the scheduling process. Like classes, room structures are also linked with each other.

4. Data structures for time periods to keep track of which hour or time slot was occupied and which was not.

5. Department inclusion data structure giving department inclusion within other larger departments or colleges.

6. Students structures indicating classes of various degree of requirements and preferences for each student.

The basic function of the system is as follows: given data files of classes, rooms and buildings, department-to-building distance matrix, students data, and the inclusion data, using the abovementioned data structures, the system builds an internal database which in turn is used in carrying out the scheduling process. This process involves a number of essential sub-processes such as checking the distances between buildings, checking building, room type and hours occupied, checking and comparing time slots for any conflicts, checking rooms for any space conflict, and keeping track of and updating the hours already scheduled.

The rule-based system uses an iterative approach. The basic procedure for each iteration is as follows. The scheduling of classes is done by department, so each iteration consists of a loop over all departments. The departments are chosen in order of size, with those having the most classes being scheduled first. The system first loops over all the currently unscheduled classes, and attempts to assign them to the first unoccupied room and timeslot that satisfies all the rules governing the constraints. Since constraints involving capacity of rooms are very difficult to satisfy, larger classes are scheduled first, to try to avoid not having large enough rooms later for those class sections with large enrollments.

In some cases the only rooms and timeslots that satisfy all the rules will already be occupied by previously scheduled classes. In that case, the system attempts to move one of these classes into a free room and timeslot, to allow the unscheduled class to be scheduled.

Next, the system searches through all the scheduled classes, and selects those that have a high cost, by checking the medium and soft constraints such as how closely the room size matches the class size, how many students have time conflicts, whether the class is in a preferred time period or a preferred building, and so on. Selecting threshold values for defining what is considered a “high” cost in each case is a subjective procedure, but it is straightforward to choose reasonable values. When a poorly scheduled class is identified, the system searches for a class to swap it with, so that the hard constraints are still satisfied, but the overall cost of the medium and soft constraints is reduced.

This process of swapping rooms continues provided all the rules are satisfied and no “cycling” (swapping of the same classes) occurs. Once all the departments have been considered, this completes one iteration. The system continues to follow this iterative procedure until a complete iteration produces no changes to the schedule.

4.2 Examples of Rule-based Constraints

There are many rules dealing with space and hours, type of room, and priority of room. Many are quite complex, but some of the basic rules, such as those implementing the hard constraints, can be quite straightforward.

4.2.1 Time and Space

The following is the basic rule for dealing with time and space conflicts for a room:

IF [room(capacity) > class(space-requested)] and [no time conflict in this room] THEN assign the room to the class.

[If more than a room having an identical availability time are present in the database, then choose the one with the minimal magnitude of the difference between room(capacity) and course(space-requested)].

\underline{Input:} CPS-615 requests space for 75 students and meeting hours of [MWF 10:40 - 11:35 am].

In the database we have, for example, the following rooms:

• room: 1-218 CST , capacity: 45 students, and pre-occupied hours: none.

• room: 207 HL, capacity 120 students, and pre-occupied hours: [MW 09:45 - 11:10 am].

• room: Kitt-Aud HBC, capacity: 127 students, and pre-occupied hours: [TTH 10:40 - 12:00 am].

\underline{Output:} Room Kitt-Aud is assigned to CPS-615.

4.2.2 Type of Room

IF [room(capacity) > course(space-requested)] and [no time conflict in assigning this room] and [type(room) matches type(requested)] THEN assign the room to the course: [If identical rooms in type are present in the database, then choose the one with the minimal magnitude of the difference between room(capacity) and course(spacerequested)].

\underline{Input:} BIO-615 requests space (not biology lab) for 35 students and meeting hours in [MW 10:40 - 11:35].

In the database we have, for example, the following rooms:

• room: 1-218 CST, capacity: 45 students, type: classroom, and pre-occupied hours: none.

• room: 124-SIMS, capacity: 40 students, type: laboratory, and pre-occupied hours: none.

\underline{Output:} Room 124-SIMS is assigned to BIO-615.

4.2.3 Room Priority

Rooms are classified as:

1. General Purpose: used by all users.

2. Priority Rooms: are used based on certain priority of the class or the department.

3. Special Purpose: are only assigned to specific users.

• IF [room(priority) == general] and [room(capacity) > course(space-requested)] and [no time conflict in assigning this room] and [type(room) matches type(requested)] THEN assign any user to this room.

• IF [room has a priority user] and [the user is already assigned] and [room(capacity) > course(space-requested)] and [no time conflict in assigning this room] and [type(room) matches type(requested)] THEN assign other users.

• IF [room(priority) == special purpose] and [user belongs to the special users of this room] and [room(capacity) > course(space-requested)] and [no time conflict in assigning this room] and [type(room) matches type(requested)] THEN assign the room to the user.

• IF [room has a priority user] and [the user is already assigned] and [room(capacity) > course(space-requested)] and [no time conflict in assigning this room] and [type(room) matches type(requested)] THEN assign other users.

• IF [room(priority) == special purpose] and [user belongs to the special users of this room] and [room(capacity) > course(space-requested)] and [no time conflict in assigning this room] and [type(room) matches type(requested)] THEN assign the room to the user.

4.2.4 The Rationale of Using a Preprocessor

• It is a good practice when using heuristics such as simulated annealing, to start at a good point in the search space or with a partial solution of the problem; regardless of the difficulty of the problem.

• A good starting point can only have a positive effect on the overall solution, and that is improving its quality.

• A good starting point or a partial solution will also have a positive effect on the convergence of the system; faster convergence to the final solution.

When the rule-based system is used as a preprocessor, it produces a partial schedule as an output, since it is usually unable to assign all of the given classes to rooms and times slots. The output is divided into two parts: the first consists of classes, with their associated professors and students, assigned to various rooms; and the second is a list of classes that could not be assigned due to constraint conflicts.

Simulated Annealing

Simulated annealing (SA) has been widely used for tackling different combinatorial optimization problems, particularly academic scheduling [?, ?, ?]. The basic algorithm is described in Figure 5.1. 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 [?].

A comprehensive discussion of the theoretical and practical details of SA is given in [?, ?, ?, ?]. 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(−ΔC/T), where T is the temperature, if it would increase the cost by Δ(C). Also, in Figure 5.1, s is the current schedule and s′ is a neighboring schedule obtained from the current neighborhood space (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.

The SA algorithm has advantages and disadvantages compared to other global 1. Generate an initial schedule s.

2. Set the initial best schedule s^{∗} = s.

3. Compute cost of s : C(s).

4. Compute initial temperature T_{0}.

5. Set the temperature T = T_{0}.

6. While stop criterion is not satisfied do:

(a) Repeat Markov chain length (M) times: i. Select a random neighbor s′ to the current schedule, (s′ ⊂ N_{s}) . ii. Set Δ(C) = C(s′) − C(s) . iii. If (Δ(C) ≤ 0 {downhill move}):

• Set s = s′ .

• If C(s) < C(s^{∗}) then set s^{∗} = s. iv. If (Δ(C) > 0 {uphill move}):

• Choose a random number r uniformly from [0, 1].

• If r < e^{−Δ(C)/T} then set s = s′.

(b) Reduce (or update) temperature T.

7. Return the schedule s^{∗}.

Figure 5.1: The Simulated Annealing Algorithm

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 et al. [?].

5.1 Timetabling Using the Annealing Algorithm

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

1. a state is a timetable containing the following sets:

• P: a set of professors.

• C: a set of classes.

• S: a set of students.

• R: a set of classrooms.

• I: a set of time intervals.

2. a cost or “energy” E(P, C, S, R, I) such that:

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

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

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

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

• 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.).

3. A 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.

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.

5.2 Interpretation of the Solution Space

After setting up or mapping the problem we need to determine the following:

• A feasible solution is defined as a schedule that satisfies a subset (C′) of constraints (C) : C′ ⊆ C ; provided that all hard constraints are part of C′ .

• A neighboring solution is obtained by modifying the current solution space (i.e. swapping two or more classes in time).

• Here are two guidelines in determining whether a given constraint to be included in C′ or not:

1. 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).

2. A conflict between two or more classes should not be a sufficient condition for hindering an assignment.

3. If all courses last one time period, this is equivalent to changing the color of a vertex, for example.

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.

5.3 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” 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 [see Sorkin -ref 62], 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},…,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.

5.3.1 Cooling Parameters

Each cooling takes a set of parameters:

1. the temperature reduction factor,

2. time at each temperature,

3. initial temperature, and

4. final temperature,

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.

5.3.2 Dependence on Initial Temperature

…

5.4 Our Experimental Schedules

Three annealing schedules have been used in our experiments to update the temperature of the SA algorithm in Figure 5.1: geometric cooling, adaptive cooling, and adaptive reheating as a function of cost.

The first schedule we have used is geometric cooling, where the new temperature (T′) of the SA algorithm is computed using

T′ = αT ,  (5.1)

where α (0 < α < 1) denotes the cooling factor. Typically the value of α 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 reheating as a function of cost (RFC), which was used for timetabling by Abramson et al. [?], but the ideas behind it are due to Kirkpatrick et al. [?, ?] and White [?]. 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 [?, ?] can be observed as sub-parts of the problem are resolved. In some of the work dealing with the traveling salesman problem using annealing [?], 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 [?], 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 [?, ?, ?].

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

At temperature T, the probability distribution for configurations is:

P_{i}(T) = \frac{e^{\frac{−C_{i}}{kT}}}{∑_{j} e^{\frac{−C_{j}}{kT}}} .  (5.2)

The average cost is computed as:

< C(T) >= ∑_{i∈C} C_{i}P_{i}(T) .  (5.3)

Therefore, the average square cost is:

< C^{2}(T) >= ∑_{i∈C} C_{i}^{2}P_{i}(T) .  (5.4)

The variance of the cost is:

σ^{2}(T) =< C^{2}(T) > − < C(T) >^{2} .  (5.5)

Now, the specific heat is defined as:

C_{H}(T) = \frac{σ^{2}(T)}{T^{2}} .  (5.6)

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

T = K · C_{b} + T (C_{H}^{max}) ,  (5.7)

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 5.7), 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 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 et al. [?], 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 σ(T_{j}) from equation 5.5, the new temperature T_{j+1} is computed as follows:

T_{j+1} = T_{j} · e^{− \frac{aT_{j}}{σ̄(T_{j})}} ,  (5.8)

where a is a tunable parameter. Following suggestions by Otten and van Ginneken [?] and Diekmann et al. [?], σ(T_{j}) is smoothed out in order to avoid any dependencies of the temperature decrement on large changes in the standard deviation σ. We used the following standard method to provide a smoothed standard deviation σ̄:

σ̄(T_{j+1}) = (1 − ω)σ(T_{j+1}) + ωσ(T_{j})\frac{T_{j+1}}{T_{j}}  (5.9)

and set ω to 0.95. This smoothing function is used because it follows (from the form of the Boltzmann distribution, see [?, ?]) that it preserves the key relationship:

\frac{d}{dT} C(T) = \frac{σ̄^{2}(T)}{T^{2}} = C_{H}  (5.10)

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

Multi-Phase Approach to

Scheduling

6.1 Our Method

6.2 Experimental Results For Class Scheduling

6.2.1 The Choice of Moves

The performance of any application of simulated annealing is highly dependent on the method used to select a new trial configuration of the system for the Metropolis update. In order for the annealing algorithm to work well, it must be able to effectively sample the parameter space, which can only be done with efficient moves.

The simplest method for choosing a move is to swap the rooms or timeslots of two randomly selected classes. However this is extremely inefficient, since most of the time random swapping of classes will increase the overall cost, especially if we are already close to obtaining a valid solution (i.e. at low temperature), and will likely be rejected in the Metropolis procedure. This low acceptance of the moves means this simple method is very inefficient, since a lot of computation is required to compute the change in cost and do the Metropolis step, only to reject the move.

What is needed is a strategy for choosing moves that are more likely to be accepted. A simple example is in the choice of room. If we randomly choose a new room from the list of all rooms, it will most likely be rejected, since it may be too small for the class, or an auditorium when, for example, a laboratory is needed. One possibility is to create a subset of all the rooms which fulfill the hard constraints on the room for that particular class, such as the size and type of room. Now we just make a random selection for a room for that class only from this subset of feasible rooms, with an acceptance probability that is sure to be much higher. In addition, each class in our data set comes with a “type-of-space-needed” tag which is used along with other information to assign the class to the right room. This effectively separates the updates into independent sets based on room type, so for example, laboratories are scheduled separately from lectures. In our method we carry out the scheduling of lectures first, followed by scheduling of laboratories making sure that during the course of this process no lecture and its associated laboratory are scheduled in the same time period.

In effect, we have embedded a simple expert system into the annealing algorithm in order to improve the choice of moves, as well as using a more complex expert system as a preprocessor for the annealing step. When used to choose the moves for annealing, the main function of the rule-based system is to ensure that all the trial moves satisfy the hard constraints. Many of the rules dealing with the medium and soft constraints are softened or eliminated, since reducing the cost of these constraints is done using the Metropolis update in the annealing algorithm.

Another of the modifications to the rule-based system is that while the version used in the preprocessor is completely deterministic, the version used in choosing the moves for annealing selects at random from multiple possibilities that satisfy the rules equally well. This extra freedom in choosing new schedules, plus the extra degree of randomness inherent in the annealing update, helps prevent the system from getting trapped in a local minimum before it can reach a valid schedule, which is the problem with the standard deterministic rule-based system.

To improve further on the move strategy, we can take the subset of possible move choices that we have created for each class, and choose from them probabilistically rather than randomly. There may be certain kinds of moves that are more likely to be effective, so our move strategy is to select these moves with a higher probability. For example, swapping a higher level class (e.g. graduate) with a lower level class (e.g. a first or a second year type) generally has a higher acceptance, since there is little overlap between students taking these classes. Furthermore, we have experimented with two kinds of swaps, those that only involve classes offered by the same department or college and the second, swaps between classes of different departments and colleges.

Generally, the swap methods we have taken here can be considered as heuristics for pruning the neighborhood or narrowing the search space, which provides much more efficient moves and in turn an overall improvement in the results.

6.3 More Analysis of the Simulation

Our computations were done with a number of goals in mind. The main objective was to provide a schedule which satisfied all hard constraints and minimized the cost of medium and soft constraints, using real-life data sets for a large university. We also aimed to find an acceptable set of annealing parameters and move strategies for general class scheduling problems of this kind, and to study the effect of using a preprocessor to provide the annealing program with a good starting point. Finally, we wanted to make a comparison of the performance of the three different cooling schedules, geometric cooling, adaptive cooling, and reheating based on cost.

We spent quite some time finding optimal values for the various parameters for the annealing schedule, such as the initial temperature, the parameters controlling the rate of cooling (α for geometric cooling, a for adaptive cooling) and reheating (K), and the number of iterations at each temperature. Johnson et al. [?] noted in their SA implementation for the traveling salesman problem (TSP) that the number of steps at each temperature (or the size of the Markov chain) needed to be at least proportional to the “neighborhood” size in order to maintain a high-quality result.

From our experiments we found the same to be true for the scheduling problem, even though it is very different from the TSP. Furthermore, in a few tests for one semester we fixed the number of classes and professors but varied the number of rooms and time slots, and found that the final result improves as the number of iterations in the Markov chain becomes proportional to a combination of the number of classes, rooms and time slots. We also observed the same behavior when we fixed the number of rooms and time slots but varied number of classes.

Our study case involved real scheduling data covering three semesters at Syracuse University. The size and type of the three-semester data is shown in Table 6.1. Nine types of rooms were used: auditoriums, classrooms, computer clusters, conference rooms, seminar rooms, studios, laboratories, theaters, and unspecified types. Staff and teaching assistants are considered part of the set of professors. Third semester (summer) data was much smaller than other semesters, however, there were additional space and time constraints and fewer available rooms. Our data was quite large in comparison to data used by other researchers. For example, high school data used by Peterson and colleagues [?, ?] consists of approximately 1000 students, 20 different possible majors, and an overall periodic school schedule (over weeks). In the case of Abramson et al. [?], their data set was created randomly and was relatively small, and they stated that problems involving more than 300 tuples were very difficult to solve.

Table 6.1 lists all major components of the data we have used. Timetabling problems can be characterized by their sparseness. After the required number of classes N_{l} have been scheduled, there will be N_{sp} = (N_{x}N_{t} − N_{l}) spare space-time slots, hence, the sparseness ratio of the problem is defined as the ratio N_{sp}/(N_{x}N_{t}). The denser the problem, the lower the sparseness ratio, and the harder the problem is to solve. Also, for dense problems, there is an additional correlation involving the problem size. Table 6.2 shows the sparseness of the three-semester data. For university scheduling, the sparseness ratio generally decreases as the data size (particularly the number of classes) increases, so the problem becomes harder to solve. Including student preferences makes the problem much harder, but these are viewed as medium constraints and thus are not necessarily satisfied in a valid solution.

Our overall results are shown in Tables 6.3 and 6.4. These tables show the percentage of classes that could be scheduled in accordance with the hard constraints. In each case (apart from the expert system, which is purely deterministic), we have done 10 runs (with the same parameters, just different random numbers), and the tables show the average of the 10 runs, as well as the best and worst results. The MFA results are different only due to having different initial conditions. Each simulated annealing run takes about 10 to 20 hours on a Unix workstation, while a single MFA run takes approximately an hour and an expert system run takes close to two hours.

As expected, each of the methods did much better for the third (summer) semester data, which has a higher sparseness ratio. Our results also confirm what we expected for the different cooling schedules for simulated annealing, in that adaptive cooling performs better than geometric cooling, and reheating improves the result even further.

When a random initial configuration is used, simulated annealing performs very poorly, even worse than the expert system (ES). However, there is a dramatic improvement in performance when a preprocessor is used to provide a good starting point for the annealing. In that case, using the best cooling schedule of adaptive cooling with reheating as a function of cost, we are able to find a valid class schedule every time.

In the case of mean-field annealing, the overall results are generally below those of SA and ES. In addition, we have found in the implementation of this method that the results were quite sensitive to the size of the data as well the type of constraints involved. If we confine ourselves to the set of hard constraints, the results are as good as or even better than the other methods. However if we take into account the medium and soft constraints, that is, the overall cost function, this method does not perform as well.

Student preferences are included only as medium constraints in our implementation, meaning that these do not have to be satisfied for a valid solution, but they have a high priority. For the valid schedules we have produced, approximately 75% of the student preferences were satisfied. This is reasonably good (particularly since other approaches do not deal with student preferences at all), but we are working to improve upon this result.

6.4 Conclusions

We have successfully applied simulated annealing to the difficult problem of academic scheduling for a large university. Feasible schedules were obtained for real data sets, including student preferences, without requiring enormous computational effort.

Mean-field annealing works well for small scheduling problems, but does not appear to scale well to large problems with many complex constraints. For this problem, both simulated annealing and the rule-based system were more effective than MFA. It is more difficult to tune the parameters for MFA than for simulated annealing, and because of the complexity and size of the Potts neural encoding, there seems to be no clear way of preserving the state of a good initial configuration provided by a preprocessor when using MFA.

Using a preprocessor to provide a good initial state greatly improved the quality of the results for simulated annealing. In theory, using a good initial state should not be necessary, and any initial state should give a good result, however in practice, we do not have an ideal cooling schedule for annealing, or an ideal method for choosing trial moves and efficiently exploring the search space, and there are restrictions on how long the simulation can take. In general, for very hard problems with large parameter spaces that can be difficult to search efficiently, and for which very slow cooling would be much too time-consuming, we might expect that a good initial solution would be helpful. We used a fairly complex rule-based expert system for the preprocessor, however the type of preprocessor may not be crucial. Other fast heuristics could possibly be used, for example a graph coloring approach [?], or it may be possible to just utilize the schedule from the same semester for the previous year. A modified version of the rule-based system was used to choose the trial moves Table 6.1: Size of the data set for each of the three semesters.

First Semester Second Semester Third Semester
Rooms  509  509  120
Classes  3839  3590  687
Professors  1190  1200  334
Students  13653  13653  2600
Buildings  43  43  11
Schools and/or Colleges  20  21  17
Departments or  143  141  108

Course Prefixes

Areas of Study (majors)  200  200  200

for the simulated annealing, and the high acceptance rate provided by this system was crucial to obtaining good results.

As expected, for the simulated annealing, adaptive cooling performed better than geometric cooling, and using reheating improved the results even further. The best results were obtained using simulated annealing with adaptive cooling and reheating as a function of cost, and with a rule-based preprocessor to provide a good initial solution. Using this method, and with careful selection of parameters and update steps, we were able to generate solutions to the class scheduling problem using real data for a large university. None of the other methods were able to provide a complete solution.

Our main conclusion from this work is that simulated annealing, with a good cooling schedule, optimized parameters, carefully selected update moves, and a good initial solution provided by a preprocessor, can be used to solve the academic scheduling problem at a large university, including student preferences. Similar approaches should prove fruitful for other difficult scheduling problems.

CHAPTER 6. MULTI-PHASE APPROACH TO SCHEDULING 57 Table 6.2: The sparseness ratios of the problem for the data sets for each of the three semesters. Lower values indicate a harder problem.

Academic  Sparseness
Time Period  ratio
First Semester  0.50
Second Semester  0.53
Third Semester  0.62

Table 6.3: Percentage of classes scheduled using the different methods. The averages and highest and lowest values were obtained using 10 independent runs for simulated annealing (SA) and mean-field annealing (MFA). The expert system (ES) is deterministic so the results are from a single run. No preprocessor was used with the three methods.

Academic  Algorithm  Scheduled Highest  Lowest
Time Period  (average) Scheduled Scheduled
%  %  %
First Semester SA (geometric)  65.00  67.50  56.80
SA (adaptive)  67.80  70.15  61.20
SA (cost-based)  70.20  72.28  68.80
ES  76.65  76.65  76.65
MFA  65.60  71.00  61.00
Second Semester SA (geometric)  65.65  68.00  57.10
SA (adaptive)  68.50  70.10  60.77
SA (cost-based)  75.14  77.68  70.82
ES  79.00  79.00  79.00
MFA  67.20  75.00  65.00
Third Semester SA (geometric)  83.10  86.44  68.50
SA (adaptive)  85.80  89.00  70.75
SA (cost-based)  91.20  95.18  85.00
ES  96.80  96.80  96.80
MFA  88.00  95.00  82.00

Table 6.4: Percentage of scheduled classes, averaged over 10 runs of the same initial temperature and other parameters, for three terms using simulated annealing with an expert system as preprocessor.

Academic  Algorithm  Scheduled Highest  Lowest
Time Period  (average) Scheduled Scheduled
%  %  %
First Semester SA (geometric)  93.90  95.12  85.20
SA (adaptive)  98.80  99.20  95.00
SA (cost-based)  100.0  100.0  100.0
Second Semester SA (geometric)  95.00  98.95  89.40
SA (adaptive)  99.00  99.50  98.50
SA (cost-based)  100.0  100.0  100.0
Third Semester SA (geometric)  97.60  98.88  90.90
SA (adaptive)  100.0  100.0  100.0
SA (cost-based)  100.0  100.0  100.0

Optimization Networks

As a prelude to mean field annealing (MFA), a summary of the Hopfield model is given below.

7.1 The Hopfield Network

7.1.1 Some Background

The network described by Hopfield in his seminal 1982 paper (Hopfield-82) is in fact a special case of the additive model developed by Grossberg in the 1960’s. Hopfield reviewed the network’s application as a content addressable memory with both binary (Hopfield-82) and continuous valued (Hopfield-84) outputs. However, after detailed investigations into the network’s performance (see for example (Psaltis-89), it is now clear that other types of content addressable memory are far more efficient (Tarassenko-91).

The Hopfield network was first proposed for combinatorial optimization applications in (Hopfield-85 and 86), where a penalty function mapping for the TSP was given. The mapping was of a rather 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 (Wilson-88). Nevertheless, some researchers preserved with this approach, sometimes proposing modified networks to correct what is essentially a sloppy mapping.

Real progress in the field 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.

• Convergence: Generalized feedback neural networks can exhibit oscillatory (Cohen-83) or even chaotic (Lewis-92) 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 (Bruck-90).

• Mapping: The development of rigorous mappings was delayed by the popularity of the original, ad hoc mapping of the TSP problem. Researchers concentrated on finding better weights for the penalty functions, often by way of an eigen vector analysis of Liapunov function ((Aiyer-90): this, though correct, is quite complicated. Simpler, more reliable mappings began to appear for specific problems in (Aiyer-91); some also dealt with inequality constraints (Tagliarini-91). One particular rigorous mapping for the most general 0-1 quadratic programming problems onto standard Hopfield network dynamics was presented in (Gee-92).

• 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 hysteretc annealing (Eberhart-91), convex relaxation (Ogier-90), and matrix graduated non-convexity (Aiyer-91).

7.1.2 Structure and Dynamics of the Net

• 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_{λ}(u_{i}).

• The net’s dynamics are governed by the first order differential equation:

\frac{du}{dt} = −ηu + Wv + I^{b} .  (7.1)

• g_{λ}(u_{i}) = \frac{1}{1+e^{−2λu_{i}}} , for large λ, this function approximates a step function, forcing v_{i} to [0, 1] ⟹ (the output of the net lies on the unit hypercube).

• W_{ij} is the connection strength between neuron i and neuron j.

• I_{i}^{b} is an input bias to neuron i.

• Neurons may be updated individually at random (asynchronously) or all together (synchronously).

7.1.3 Operation of the Network

The network (see Figure 7.1) operates as follows:

• The network’s dynamics are governed by a set of nonlinear coupled differential equations of the form:

\frac{du_{i}}{dt} = −u_{i} + ∑_{j=1}^{n} W_{ij}v_{j} + I_{i}^{b} .  (7.2)

or the general form of it:

Figure 7.1: A schematic view of the Hopfield network.

\frac{du}{dt} = −ηu + Wv + I^{b} .  (7.3)

where u, W, v, and I are vectors.

• The I_{i}^{b} corresponds to an input bias to neuron i.

• The η term introduces an element of decay into the network’s dynamics, which can be used to ensure that u variables remain bounded in magnitude.

7.1.4 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 \frac{d}{dt}u_{i} as a function of the total input to the i^{th} neuron.

u̇ = −ηu + Wv + I^{b} ,  (7.4)

where

v = g(u) .  (7.5)

Equations 7.4 and 7.5 govern the dynamics of the network.

7.1.5 Stability of the Net

Given that W is symmetric, Hopfield proved stability for the network by showing that with the above dynamic equations (7.4 and 7.5) the network output v evolves so as to minimize the Liapunov function:

E^{liap}(v) = (−\frac{1}{2} ∑_{ij} W_{ij}v_{i}v_{j} − ∑_{i} I_{i}^{b}v_{i}) + η∑_{i} ∫_{0}^{v_{i}} g^{−1}(V) dV .  (7.6)

• v = (v_{1} … v_{n}) is denoting the net state, therefore we have an energy function defined throughout v ∈ [0, 1]^{n} and dynamics which guarantee convergence to minima.

• For a suitable 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 η.

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

7.1.7 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 7.6) for which the minima correspond to solutions and the depth of each minimum reflects the solution quality.

Hopfield and Tank used the TSP problem as an example, which they have mapped onto the net using a permutation matrix representation of the TSP tours.

7.1.8 Shortcomings of the Hopfield Model

The net’s partition function sums over a vast number of configurations (e.g. for the TSP, all of the 2^{N}^{2} configurations are admissible). 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.

7.2 Mean Field Annealing for Timetabling

7.2.1 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 mean field approximation involves averaging statistics of two coupled spins, and replacing a complex function with its truncated Taylor expansion around a 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 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 TSP, this can be done by constraining each city to be “on” only once, by enforcing ∑_{a} V_{ia} = 1 ∀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 Hopfield model.

This idea was first utilized and analyzed for the TSP case by, among others, Peterson and Söderberg (ref: Peterson-89). 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 (ref: Peterson-91). In (ref: Peterson-89) the normalization (definition 7.2.1) occurs in the context of a 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^{∗}. The idea is to regard each hypercube vertex as a configuration, with an associated energy (or cost) 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, saddle point equations are derived (see next sections), the solutions of which pick out the dominant states of the network at the current temperature T.

∗This equivalence was also noted in Hopfield and Tanks’s 1985 paper (ref: Hop-85). Martin Simmen brought this point to our attention.

This statistical mechanics characterization gives a somewhat clearer picture of why the normalized model (which corresponds to a type of Potts model (ref: Wu-82) in physics, in that each city represented by a spin which can be in only one of N states) ought to perform better than Hopfield model. As indicated in (ref: Peterson-89), the partition function of Hopfield algorithm sums over a vast number of configurations which are nothing like valid tours in the TSP, and even though these offenders have small Boltzmann weights their large number inevitably affects the thermal average quantities. 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.

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 spacetime 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 ∑_{a} S_{ia} = 1.

For a full derivation of the mean-field annealing algorithm from its roots in statistical physics, see Krogh and Hertz (ref: Hertz-91).

Definition 7.2.1 In (Peterson-89), it was stated that the performance of MFA on certain problems can be improved using neuron normalization. This idea was originally inspired by an analogy with Potts glasses in statistical physics (Kanter-87). Neuron normalization modifies the MFA algorithm so that constraints of the form ∑_{a∈ States} ν_{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.

7.2.2 Brief Overview

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 saddle point equations: \frac{∂E′}{∂V_{ia}} = 0 and \frac{∂E′}{∂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:

U_{ia} = −\frac{1}{T} \frac{∂E}{∂V_{ia}}  (7.7)
V_{ia} = \frac{e^{U_{ia}}}{∑_{b} e^{U_{ib}}} .  (7.8)

The MFA algorithm involves solving equations 7.7 and 7.8 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 7.8 around the trivial fixed-point (ref: Peterson-89, Peterson-91) V_{ia}^{(0)} = \frac{1}{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_{pq;xt}^{(0)} = \frac{1}{N_{x}N_{t}}.

Equations 7.7 and 7.8 can be solved iteratively using either synchronous or serial updating. The iterative dynamics to evolve the mean field variables toward a selfconsistent solution is explained in detail by Peterson at al. (ref: Peterson-91). The solutions correspond to stable states of the Hopfield model (ref: Hopfield-85).

Observe from equation 7.8 that any solution to the MF equations respects a continuous version of the Potts condition

∑_{a} V_{ia} = 1 ∀ i.  (7.9)

7.2.3 The Mean-Field Annealing Algorithm

The generic MFA algorithm appears in Figure 7.2. At high temperatures T, the meanfield 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:

• The values of the coefficients of terms in the energy function.

• The types of dynamics used to find solutions of the MFT equations at each T.

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

Peterson et al. (ref: Peterson-92) introduced a quantity called saturation, Σ, defined as

Σ = \frac{1}{N_{i}} ∑_{ia} V_{ia}^{2} ,  (7.10)

1. Choose a problem and encode the constraints into weights {w_{ij}}.

2. Find the approximate phase transition temperature by linearizing equation (7.8).

3. Add a self-coupling β-term if necessary. In a neural net, this corresponds to a feedback connection from a neuron to itself.

4. Initialize the neurons V_{ia} to high temperature values \frac{1}{N_{a}} plus a small random term such as rand[−1, 1] × 0.001; and set T (0) = T_{c}.

5. Until (Σ ≥ 0.99) do:

• At each T (n), update all U_{ia} and V_{ia} by iterating to a solution of the mean field equations.

• T (n + 1) = αT (n), we chose α = 0.9

6. 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).

7. Perform greedy heuristics if needed to account for possible imbalances or rule violations.

Figure 7.2: The Generic Mean-Field Annealing Algorithm

where N_{i} is the number of events (in this case the number of professor-class pairs). This characterizes the degree of clustering^{†} of events in time and/or space, Σ_{min} = \frac{1}{N_{a}} corresponds to high temperature, whereas Σ_{max} = 1 means that all the V_{ia} have converged to 0 or 1 values, indicating that each event has been assigned to a spacetime slot.

7.2.4 Observations on the Algorithm

The first step of Figure 7.2 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 7.2) performs one update per neural variable (defined as one sweep) with sequential updating using equations 7.7 and 7.8. After reaching a saturation value close to 1 (we chose Σ = 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 et al. (ref: Peterson-92).

†Other measures could also do this, e.g. an entropy measure − ∑_{ia} ν_{ia}ln(ν_{ia}).

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.

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

Definition 7.2.2 Let S_{pq;xt} denote the Potts (or multi-state) neuron, where

S_{pq;xt} = \{ \begin{matrix}1 & if event (p, q) takes place in (x, t) , \\ 0 & otherwise. \end{matrix}  (7.11)

Definition 7.2.3 We have dealt with two types of constraints:

1. Those involving professors, rooms, subjects, and time periods.

2. Those dealing with students scheduling and preferences.

In addition, some of these constraints are usually classified to be hard, others are medium and the rest are soft.

Definition 7.2.4 Here, we define a number of hard constraints of the first type (see definition 7.2.3) .

• An event (p, q) should occupy precisely one space-time slot (x, t).

This constraint can be embedded in the neural net in terms of the Potts normalization condition:

∑_{x,t} S_{pq,xt} = 1  (7.12)

• Different events (p_{1}, q_{1}) and (p_{2}, q_{2}) should not occupy the same space-time slot (x, t).

E_{XT} = \frac{1}{2} ∑_{x,t} ∑_{p_{1},q_{1}} ∑_{p_{2},q_{2}} S_{p_{1}q_{1};xt}S_{p_{2}q_{2};xt} = \frac{1}{2} ∑_{x,t} (∑_{pq} S_{pq;xt})^{2}  (7.13)

• A teacher p should have at most one class at a time.

E_{P T} = \frac{1}{2} ∑_{p,t} ∑_{q_{1},x_{1}} ∑_{q_{2},x_{2}} S_{pq_{1};x_{1}t}S_{pq_{2};x_{2}t} = \frac{1}{2} ∑_{p,t} (∑_{qx} S_{pq;xt})^{2}  (7.14)

• A class q should have at most one teacher at a time.

E_{QT} = \frac{1}{2} ∑_{q,t} ∑_{p_{1},x_{1}} ∑_{p_{2},x_{2}} S_{p_{1}q;x_{1}t}S_{p_{2}q;x_{2}t} = \frac{1}{2} ∑_{q,t} (∑_{px} S_{pq;xt})^{2}  (7.15)

In (ref: Hop-84, Hop-85), 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 Temperature, the behavior of x_{i} can be stated as:

x_{i} = \frac{1}{2} + \frac{1}{2} · tanh ((∑_{i≠j} W_{ij}x_{j} − I_{i}) × \frac{1}{T} )  (7.16)

When T ⟶ 0, function (7.16)’s response takes on a step-wise form. These models of (ref: Hop-85) 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.

Definition 7.2.5 At a finite T, x_{i} (of eq. 7.16) can be replaced by an MFT variable υ_{ia}, then differentiating to obtain the Potts MFT equations:

υ_{ia} = −\frac{1}{T} · \frac{∂E[ν⃗]}{∂ν_{ia}}  (7.17)
ν_{ia} = \frac{e^{υ_{ia}}}{∑_{b} e^{υ_{ib}}}  (7.18)

Which enforces

∑_{a} ν_{ia} = 1 ∀i .  (7.19)

Where a denotes the vector components or all the elements of the group. Also, we could state that ν_{ia} ≈ Probability(neuron(i) ∈ state(a)). Equations (7.17, 7.18) are computed in an iterative manner.

Definition 7.2.6 To carry out the minimization of the energy functions of definition (7.2.2) we need to convert the neural variables {0, 1} of S_{i} to their corresponding mean field variables at T. These correspondences are derived (see (ref: Peterson-92)) to be:

V_{pq;xt} ≈< S_{pq;xt} >_{T} ;  (7.20)

which is the thermal average of the original binary spines, with the following multistate MF T equations (see (7.17, 7.18) for a general format):

U_{pq;xt} = −\frac{1}{T} · \frac{∂E}{∂V_{pq;xt}}  (7.21)
V_{pq;xt} = \frac{e^{U_{pq;xt}}}{∑_{x′t′} e^{U_{pq;x′ t′}}}  (7.22)

The MF T equations of υ and ν are solved 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 X-neurons S_{pq;x}^{(X)} and T ime-neurons S_{pq;t}^{(T ime)} :

Definition 7.2.7

S_{pq;xt} = S_{pq;x}^{(X)} · S_{pq;t}^{(T ime)}  (7.23)

With separate Potts conditions:

∑_{x} S_{pq;x}^{(X)} = 1  (7.24)
∑_{t} S_{pq;t}^{(T ime)} = 1  (7.25)

For X-neurons and T ime-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 exits in a look-up table which contains all the relevant information to process each event (i). Now, the above (7.24, 7.25) conditions can be re-stated as follows:

∑_{x} S_{i;x}^{(X)} = 1  (7.26)
∑_{t} S_{i;t}^{(T ime)} = 1  (7.27)

During annealing, as we go from a higher temperature to a lower one, a phase transition is passed at (T = T_{c}), where T_{c}^{‡} is the critical temperature. The value of T_{c} set the scale or the range of the value of T, and it is estimated by expanding 7.22 around some trivial fixed-point such as

V_{pq;xt}^{(0)} = \frac{1}{N_{x}N_{t}} , at T = 0 .  (7.28)

‡See pp 29-31, (ref: Hertz-91) for an approximate value in the case of Ising ferromagnets.

Also T_{c} depends on the self-coupling coefficient β, and on the method of the update. This parameter β acts as a Lagrangian multiplier for the given energy constraint(s). See section on the generic Potts algorithm and (ref: Peterson-89).

Definition 7.2.8 Relative Clamping: During the revision process of 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:

i ⟶ (j, k) ;  (7.29)

where (j) denotes the cluster and (k) an event within (j). In this situation, one has a common t-neuron (S_{j;t}^{(T ime)}) for cluster j but distinct x-neurons (S_{jk;x}^{(X)}) for the individual events.

Definition 7.2.9 To deal with events that take place over two consecutive hours, what (ref: Peterson-89) called effective time-neurons S̃_{i;t}^{(T ime)} are introduced:

S̃_{i;t}^{(T ime)} = ∑_{k=0}^{g_{i}−1} S_{i;t−k}^{(T)} ;  (7.30)

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.

Utilizing definition (7.2.9), constraints (7.13, 7.14, 7.15) can be put into the following form:

Definition 7.2.10

E_{XT} = \frac{1}{2} ∑_{x,t} ∑_{i≠i′} S_{i;x}^{(X)}S_{i′;x}^{(X)}S̃_{i;t}^{(T ime)}S̃_{i′;t}^{(T ime)}  (7.31)
E_{P T} = \frac{1}{2} ∑_{p,t} ∑_{i≠i′} δ_{p(i),p(i′)}S̃_{i;t}^{(T ime)}S̃_{i′;t}^{(T ime)}  (7.32)
E_{QT} = \frac{1}{2} ∑_{q,t} ∑_{i≠i′} δ_{q(i),q(i′)}S̃_{i;t}^{(T ime)}S̃_{i′;t}^{(T ime)}  (7.33)

Definition 7.2.11 The Kronecker δ_{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.

δ_{r,s} = \{ \begin{matrix}1 & if r = s \\ 0 & otherwise. \end{matrix}  (7.34)

Definition 7.2.12 One important soft constraint that need to be handled is 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 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 δ with a more structured overlap matrix Γ which is given by a look-up table.

δ_{q(i),q(i′)} ⟶ Γ_{q(i),q(i′)} .  (7.35)

Definition 7.2.13 From definitions (7.2.8, 7.2.12), we need to make sure that (∀k) events within a given (j) cluster p(j, k) ≠ 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 δ̃_{j,k} according to:

δ̃_{j,k} = δ_{j,k}(1 − δ_{j,0}δ_{k,0}) .  (7.36)

If j = 0 then ∀k , k = 0 .

Definition 7.2.14 Utilizing the concepts of the above definition (7.2.12, 7.2.13), the collision terms of definition (7.2.10) are extended to ensure that all the hard constraints are satisfied (or the solutions obtained are legal). The extension takes the following forms:

E_{XT} (form 1) = \frac{1}{2} ∑_{jk≠j′k′} ∑_{x} S_{jk;x}^{(X)}S_{j′k′;x}^{(X)} ∑_{t} S̃_{j;t}^{(T ime)}S̃_{j′;t}^{(T ime)}  (7.37)

E_{XT} (form 2) = \frac{1}{2} ∑_{j≠j′} ∑_{kk′} ∑_{x} S_{jk;x}^{(X)}S_{j′k′;x}^{(X)} ∑_{t} S̃_{j;t}^{(T ime)}S̃_{j′;t}^{(T ime)} + \frac{1}{2} ∑_{j} g_{j} ∑_{k≠k′} ∑_{x} S_{jk;x}^{(X)}S_{jk′;x}^{(X)}

(7.38)

E_{P T} (form 1) = \frac{1}{2} ∑_{p,t} [∑_{jk} δ̃_{p,p(j,k)}S̃_{j;t}^{(T ime)}] [ ∑_{j′≠j} ∑_{k′} δ̃_{p,p(j′,k′)}S̃_{j′;t}^{(T ime)}]  (7.39)

E_{P T} (form 2) = \frac{1}{2} ∑_{j≠j′} ∑_{kk′} δ̃_{p(j,k),p(j′,k′)} ∑_{t} S̃_{j;t}^{(T ime)}S̃_{j′;t}^{(T ime)}  (7.40)

E_{QT} (form 1) = \frac{1}{2} ∑_{qt} [∑_{jk} δ_{q,q(j,k)}S̃_{j;t}^{(T ime)}] [ ∑_{j′≠j} ∑_{k′} Γ_{q,q(j′,k′)}S̃_{j′;t}^{(T ime)}]  (7.41)

E_{QT} (form 2) = \frac{1}{2} ∑_{j≠j′} ∑_{kk′} Γ_{q(j,k),q(j′,k′)} ∑_{t} S̃_{j;t}^{(T ime)}S̃_{j′;t}^{(T ime)}  (7.42)

Note: Each “form 2” above is a derivation from its corresponding “form 1”.

The factorization idea of (dfn 7.2.14) can be extended to mean field variables V_{jk;x}^{(X)} and V_{j;t}^{(T ime)}, corresponding to neurons (S_{jk;x}^{(X)}) and (S_{j;t}^{(T ime)}), respectively. V_{jk;x}^{(X)} and V_{j;t}^{(T ime)} 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 (7.26, 7.27) .

Substituting V for S in (7.26, 7.27), the mean field equations at a finite T, are given in terms of the local fields U_{jk;x}^{(X)} and U_{j;t}^{(T ime)}. Now, the mean field dynamics proceed over the following local and mean field variables:

Definition 7.2.15

U_{jk;x}^{(X)} = −\frac{1}{T} · \frac{∂E}{∂V_{jk;x}^{(X)}}  (7.43)

U_{j;x}^{(T ime)} = −\frac{1}{T} · \frac{∂E}{∂V_{j;t}^{(T ime)}}  (7.44)

V_{jk;x}^{(X)} = \frac{e^{U_{jk;x}^{(X)}}}{∑_{x′} e^{U_{jk;x′}^{(X)}}}  (7.45)

V_{j;t}^{(T ime)} = \frac{e^{U_{j;t}^{(T ime)}}}{∑_{t′} e^{U_{j;t′}^{(T ime)}}}  (7.46)

The iterations of the MF dynamics over equations (7.43 - 7.46) can either be synchronous or serial (see (ref: Peterson-92)).

Definition 7.2.16 Synchronous Updating: A single sweep consists of first computing all the local fields, and then updating all the neurons. That means do equations (7.43,7.44) then perform equations (7.45,7.46).

Definition 7.2.17 Serial Updating: For each neuron S_{jk;x}^{(X)} and/or S_{j;t}^{(T ime)}, the local field is computed immediately before the corresponding neuron state is updated. That means either do equation (7.43) followed by (7.45) and/or do equation (7.44) followed by (7.46).

Next, we need to define 2 important soft constraints encountered when scheduling classes.

Definition 7.2.18

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

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

S_{i;d}^{(D)} = ∑_{h} S̃_{j;dh}^{(T ime)} = g_{j} · S_{j;dh}^{(T ime)}  (7.47)

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

E_{QSD} = \frac{1}{2} ∑_{qsd} [∑_{jk} ∑_{h} δ_{q,q(j,k)}δ_{s,s(j,k)}S̃_{jk;dh}^{(T ime)}]^{2}  (7.48)

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

E_{QSD} = \frac{1}{2} ∑_{d} ∑_{j≠ĵ} ∑_{kh} ∑_{k̂ĥ} δ_{q(j,k),q(ĵ,k̂)}δ_{s(j,k),s(ĵ,k̂)}Ŝ_{j;dh}^{(T ime)}Ŝ_{j;dĥ}^{(T ime)}  (7.49)

Now for some day (d) of the week:

E_{QSD} = \frac{1}{2} ∑_{j≠ĵ} ∑_{kk̂} δ_{q(j,k),q(ĵ,k̂)}δ_{s(j,k),s(ĵ,k̂)} ∑_{d} S_{j,d}^{(D)}S_{ĵ,d}^{(D)}  (7.50)

• Glueing: The schedule should have as few “holes” as possible; that is lessons should be glued together.

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

E_{QDH} = −α ∑_{q} ∑_{d,h} [∑_{j,k} δ_{q,q(j,k)}S̃_{j;dh}^{(T ime)}] [∑_{ĵ≠j} Γ_{q,q(ĵ,k̂)}S̃_{j;d(h−1)}^{(T ime)}]  (7.51)

The energy term of (7.51) can be summarized to say that: 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 (k̂) from a different cluster (ĵ) to occur on hour (h − 1) of that day.

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

E_{QDH} = −α ∑_{ĵ≠j} ∑_{k,k̂} Γ_{q(j,k),q(ĵ,k̂)} ∑_{d,h>1} S̃_{j;dh}^{(T ime)}S̃_{ĵ,d(h−1)}^{(T ime)}  (7.52)

The presence of the (−) sign in front of (7.51, 7.52) indicates that such energy functions are to be minimized. Also, the (α) in front of the the above forms governs the strength of this reward relative to the energies of equations (7.13,7.14,7.15). This (α) takes on a value between zero and 1.

Finding legal solutions with no collisions corresponds to minimizing the hard energy:

E_{hard} = combine(E_{XT} , E_{P T} , E_{QT})  (7.53)

to zero. The soft energy is given by:

E_{soft} = combine(E_{QSD}, E_{QDH}) .  (7.54)

The total energy E to be minimized is then:

E = combine(E_{hard}, E_{medium}, E_{soft}) .  (7.55)

7.2.6 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 ⊆ R.

Here, we need to introduce the R-neuron: S_{r(i);(d(t)d(x)}^{(R)} where k = 1, 2, . . .. Each student r(i) ∈ R should have a balanced schedule over the two-group of days (MW F) and (T uT h). 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, …, n stand for places and time slots (or events), respectively. Using such notation, an example of the Potts normalization will be: S_{p(i)q(i);x(j)}^{(X)} and S_{p(i)q(i);t(j)}^{(T ime)} respectively, where i = 1, 2, …, n is a section index.

Definition 7.2.19 The athletes case: Let r⃗ ⊆ 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)) 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 :

1. S_{r^{∗}b;d(x)}^{(r⃗,X)} , S_{r^{∗}b;d(t)}^{(r⃗,T ime)} : It takes a value of (1) if event (r^{∗}, b) takes place in the place and time slot, respectively, of a weekday d, and zero otherwise.

2. S_{r^{+}a;d(x)}^{(r⃗,X)} , S_{r^{+}a;d(t)}^{(r⃗,T ime)} : It takes a value of (1) if event (r^{+}, a) takes place in the place and time slot, respectively, of a weekday d, and zero otherwise.

3. S_{jk;d(x)}^{(r⃗,X)} : It takes a value of (1) if event (k) of cluster (j) that is associated with the x-neuron of r⃗ occurs in place d(x), and zero otherwise.

4. S_{ĵk̂;d(x)}^{(r⃗,X)} : It takes a value of (1) if event (k̂) of cluster (ĵ) that is associated with the x-neuron of r⃗ occurs in place d(x), and zero otherwise.

5. S̃_{j;d(t)}^{(r⃗,T ime)} : It takes a value of (1) if cluster (j) that is associated with the t-neuron of r⃗ occurs in time d(t), and zero otherwise.

6. S̃_{ĵ;d(t)}^{(r⃗,T ime)} : It takes a value of (1) if cluster (ĵ) that is associated with the t-neuron of r⃗ occurs in time d(t), and zero otherwise.

In the above list, the last 2 neurons act as what Peterson calls effective timeneurons. Not only will they handle periods of single hours but also variants of longer and shorter intervals.

Definition 7.2.20 Set (ϵ) and (c) to some desired values depending on the time period of the set of classes, now, the above effective t-neurons can be defined as follows:

S̃_{j;d(t)}^{(r⃗,T)} = ∑_{k=0}^{g_{j}−c} S_{j;d(t−k+ϵ)}^{(r⃗,T)}  (7.56)

S̃_{ĵ;d(t)}^{(r⃗,T ime)} 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 (ϵ = 0). If (c = 1) and (g_{i} = 1) then (ϵ = 0) and

S̃_{j;d(t)}^{(r⃗,T ime)} = S_{j;d(t)}^{(r⃗,T ime)}.  (7.57)

Definition 7.2.21 Let S_{q(i)r(i);b(x,t)}^{(R)} 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_{q(i)r(i);(d(x)[d(t)])}^{(X[T ime])}|EC ∈ {0, 1} is set to0 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 link concept, we assumed in definition (7.2.18) 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, 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(ĵ, k̂)).

Let x̂(i) be a single location (or a seat) in class x(i) such that i = 1, 2, …, num ; where (num) is the number of class-rooms.

Day-Class Balancing

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

Let D be the set {{MW F }, {T uT h}} = {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}}.

Definition 7.2.22 The Potts condition of balancing S_{q(m)r(m);b(x,t)}^{(R)} (take a value of 1 for balanced and 0, otherwise) is accomplished through one (or all) of these energies:

E_{QSD}^{r} = \frac{1}{2} ∑_{j≠ĵ} ∑_{kk̂} δ_{q(j,k),q(ĵ,k̂)}δ_{s(j,k),s(ĵ,k̂)} ∑_{d_{1},d_{2}} S_{j,d_{1}}^{(D)}S_{ĵ,d_{2}}^{(D)}  (7.58)
E_{QSD_{1}}^{r} = \frac{1}{2} ∑_{j_{1}≠ĵ_{1}} ∑_{k_{1}k̂_{1}} δ_{q(j_{1},k_{1}),q(ĵ_{1},k̂_{1})}δ_{s(j_{1},k_{1}),s(ĵ_{1},k̂_{1})} ∑_{d_{1}} ∑_{i≠î} S_{i,d_{1}}^{(D)}S_{î,d_{1}}^{(D)}  (7.59)
E_{QSD_{2}}^{r} = \frac{1}{2} ∑_{j_{2}≠ĵ_{2}} ∑_{k_{2}k̂_{2}} δ_{q(j_{2},k_{2}),q(ĵ_{2},k̂_{2})}δ_{s(j_{2},k_{2}),s(ĵ_{2},k̂_{2})} ∑_{d_{2}} ∑_{i≠î} S_{i,d_{2}}^{(D)}S_{î,d_{2}}^{(D)}  (7.60)

The Case of Student Athletes

The main constraint of this case very much analogous to the hard constraint stated somewhere above : different events should not occupy the same space-time slot.

Definition 7.2.23

Er⃗ ; d(XT) = \frac{1}{2} ∑_{d} ∑_{xt} ∑_{ṙb} ∑_{r̈a} S_{ṙb;d(x)[d(t)]}^{(r⃗,X[T ime])}S_{r̈a;d(x)[d(t)]}^{(r⃗,X[T ime])}  (7.61)

Where (ṙb) and (r̈a) are some two different events. After some further refinement, we obtain:

Er⃗ ; d(XT) = \frac{1}{2} ∑_{j≠ĵ} ∑_{kk̂} ∑_{d(x)} S_{jk;d(x)}^{(r⃗,X)}S_{ĵk̂;d(x)}^{(r⃗,X)} ∑_{d(t)} S̃_{j;d(t)}^{(r⃗,T ime)}S̃_{ĵ;d(t)}^{(r⃗,T ime)} + \frac{1}{2} ∑_{j} g_{j} ∑_{k≠k̂} ∑_{d(t)} S_{jk;x}^{(r⃗,X)}S_{jk̂;x}^{(r⃗,X)}

(7.62)

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:

Er⃗ ; d(T ime) = \frac{1}{2} ∑_{j≠ĵ} ∑_{d(t)} S̃_{j;d(t)}^{(r⃗,T ime)}S̃_{ĵ;d(t)}^{(r⃗,T ime)}  (7.63)

Eligibility for Enrollment

Definition 7.2.24 Before going through the process of balancing, each student (r(i) = r_{s}) eligibility for enrolling in a given class need 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 :

E_{QR} = \frac{1}{2} ∑_{qr} ∑_{jk} ∑_{i} δ_{q,q(j,k)}δ_{r,r(i)}S_{ij;x(i)}^{X}|(EC ∈ \{0 , 1 \})  (7.64)

∑_{i} O_{q(i)} = 1 ⇒ Test: {EC(R) + satisfy-constraints (Hard, Soft)}  (7.65)

Course Pre-requisite and Subsequent Assignment

Utilizing the link concept (lin(s, s′)) between two different subjects (see definition 7.2.21), a general energy assignment can be put together as follow:

Definition 7.2.25 Let (i) stands for the situation at hand (also can be viewed as a space-time event).

E_{link(1)} = −γ ∑_{is} [∑_{jk} δ_{s,s(j,k)}S̃_{ln(s,s′);i}^{(X[T ime])}]^{2}  (7.66)

E_{link(2)} = −γ ∑_{s} ∑_{j≠j′} ∑_{k≠k′} δ_{s(j,k),s(j′,k′)} ∑_{i} S̃_{ln(s(j,k),s(j′,k′));i}^{(X[T ime])}S̃_{ln(s(j,k),s(j,k′));i}^{(X[T ime])}  (7.67)

E_{link(3)} = −γ ∑_{j≠j′} ∑_{k≠k′} δ_{s(j,k),s(j′,k′)}S̃_{ln(s(j,k),s(j′,k′));i}^{(X[T ime])}S̃_{ln(s(j,k),s(j,k′));i}^{(X[T ime])}  (7.68)

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

Enroll in all parts of a course There are a number course 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 all associated sections of it.

Let us assume the existence of some subject attribute s(j, k) consisting of 2 components: ṡ denoting lectures, and s̈ stands for labs. Now the energy equation depicting this assignment can be outlined as follow:

E_{enroll(1)} = \frac{1}{2} ∑_{s} ∑_{jk} ∑_{i} δ_{s,s(j,k)}δ_{r,r(i)} ∑_{t} S_{jk;t}^{(T ime)}  (7.69)
E_{enroll(2)} = \frac{1}{2} ∑_{si} ∑_{t≠t′} [∑_{jk} δ_{s,ṡ(j,k)}δ_{r,r(i)}S_{jk;t}^{(T ime)}] [∑_{jk′} δ_{s,s̈(j,k′)}δ_{r,r(i)}S_{jk′;t′}^{(T ime)}]  (7.70)
E_{enroll(3)} = \frac{1}{2} ∑_{k≠k′} δ_{s,ṡ(j,k)}δ_{s,s̈(j,k′)} ∑_{t≠t′} S_{jk;t}^{(T ime)}S_{jk′;t′}^{(T)}  (7.71)

Lecture hall capacity In order to satisfy this constraint, the following inequality need to hold.

∀x ∑_{i} ||x(i)| − ∑_{j} ∑_{k} S_{jk;x(i)}^{(X)}| ≥ 0 .  (7.72)

Student status and priority mapping

E_{assign} = \frac{1}{2} ∑_{m,x̂,x} ∑_{i≠i′} S_{x̂;x}^{(X)} (S_{r(i);(x̂(m)∈x(m))}^{(R)}|(M = Req) + S_{r(i′);(x̂(m)∈x(m))}^{(R)}|(M = χ))

(7.73)

On the Energy Equation of the Potts Model The energy function of the Potts spin model is of the form:

E = E_{A} + E_{B} + E_{C} + …  (7.74)

In (ref: Peterson-89), in which the TSP is used as study case, parameters α and β are used to control the strength of E’s restraints. The same parameters can also be used when dealing with non-Euclidean type CO problems such as scheduling.

Observation 7.2.26 For timetabling, β is split into two parts (β^{(T ime)} and β^{(X)}), one with respect to the time-based neurons and the other with respect to the spacebased neurons. E_{β} becomes as follows:

−\frac{β^{(T ime)}}{2} ∑_{j} ∑_{t} [S_{j;t}^{(T)}]^{2} − \frac{β^{(X)}}{2} ∑_{jk} ∑_{x} [S_{jk;x}^{(X)}]^{2}  (7.75)

Here, the appearance of β^{(T ime)} and β^{(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.

7.2.7 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: ∀(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 (ref: Hop-85)).

7.2.8 Some Implementation Details

The general idea is to initialize the system close to the symmetric state at a temperature near T_{c}, then, allowing 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^{§}: (Martin Simmen (ref: Simmen-93) investigated in depth these and most of the other issues pertaining to the Potts model, using the TSP as a test case.)

• Initial State Another possible initial value of ν_{ia}^{n=0} is (1 + ξ_{ia})/N, with the ξ_{ia} being random variables drawn uniformly from [−ξ, ξ]. Clearly ξ should be ≪ 1.

• Convergence At each T the ν_{ia} variables are updated as above or by the following rule^{¶}:

ν_{ia}^{(n+1)} = ν_{ia}^{(n)} + γ(f_{ia} − ν_{ia}^{(n)}), 0 < γ ≤ 1 .  (7.76)

Setting γ = 1 recovers the dynamics of (ref: Peterson-89), and the γ → 0 limit gives the first order Euler integration method for the differential equation:

dν_{ia}/dt = f_{ia} − ν_{ia} .  (7.77)

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 some sort of a “tolerance” quantity of the form ∑_{ia} (ν_{ia} − f_{ia})^{2}.

• Annealing Schedule A simple exponential scheme is used, i.e. after reaching a fixed point at T, T ↦ T · T_{r}, where T_{r} ∈ (0, 1) controls the cooling rate. Roughly speaking, T_{r} should be set to 0.9.

• 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} ≪ T_{c}), without ever satisfying the Σ > Σ^{thresh} criterion. T_{min} was set to 10^{−3}.

§These follow the method used in (ref: Peterson-89) but with some additions/improvements done by the authors -see algorithm above.

¶Thanks to M. Simmen for pointing this out. Also see (ref: Peterson-89)

7.2.9 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, which empirically have an overall O(N^{2}) running time (ref: Johnson-90). Ideally, one would like to study the likely performance of the method through numerical simulations using a small step-size γ to approximate continuous time dynamics. However, this is computationally very expensive: not only do the simulation times grow as \frac{1}{γ}, but also small values of (γ) 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 γ runs tend to depress T_{c}^{exp} artificially and thereby degrade solution quality, and as we observed, result in a number of constraint violations. Many of these constraints are hard, but the majority are classified to be soft. This later 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.

7.2.10 Mapping Timetabling onto the Potts Neural Nets

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

Neural Mapping

Involve mapping the problem onto the Potts model.

Figure 7.3: A one-way arrow indicates that the “bottom” entity depends on the “top” one. A two-way arrow indicates a mutual dependency between the two entities.

• These multi-state neurons, forming the network, represent different activities of the timetabling problem.

• Construct the constraints of the problem using these neurons.

• Derive those MF equations controlling the dynamics of the net:

– Introduce the MF variables , for example,

V_{pq;xt} =< S_{pq;xt} >_{T} .  (7.78)

– Start off with the partition function:

Z = ∑_{[s⃗]}e^{−E[s⃗]/T} .  (7.79)

– Re-write Z in terms of a multidimensional integral over continuous vectors (u⃗_{i}, v⃗_{i}).

– The MFT approximation to < s⃗_{i} > is given by the value u⃗_{i} at the saddle point equations:

\frac{∂E}{∂v_{ia}} = 0  (7.80)

and

\frac{∂E}{∂u_{ia}} = 0 ,  (7.81)

generating a set of self-consistency equations in the MF variables:

u_{ia} = −\frac{1}{T} \frac{∂E[v⃗]}{∂v_{ia}} ,  (7.82)

and

v_{ia} = \frac{e^{u_{ia}}}{∑_{b}e^{u_{ib}}} .  (7.83)

Where a denotes the vector components. For events (pq, xt), we get the following MF variables:

u_{pq,xt} = −\frac{1}{T} \frac{∂E}{∂v_{pq,xt}} ,  (7.84)

and

v_{pq,xt} = \frac{e^{u_{pq,xt}}}{∑_{x′t′}e^{u_{pq,x′t′}}} .  (7.85)

From the later follows, that ∑_{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.

• 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}.

CHAPTER 7. OPTIMIZATION NETWORKS 91

• 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

v_{pq,xt}^{(0)} = \frac{1}{N_{x}N_{t}} .  (7.86)

Scheduling as a Graph Coloring

Problem

8.1 introduction

Many scheduling problems involve allowing for a number of pairwise restrictions on which jobs can be done simultaneously. For instance, in attempting to schedule classes at a university, two courses taught by the same faculty member cannot be scheduled for the same time slot. Similarly, two courses that are required by the same group of students also should not conflict. The problem of determining the minimum number of time slots needed subject to these restrictions is a graph coloring problem (GC).

GC has considerable application to a large variety of complex problems involving optimization. Among those problems is the one we have dealt with, which is a large scale academic class scheduling. Constraints for these kind of problems –class scheduling or timetabling– are usually expressible in the form of pairs of incompatible objects (e.g., pairs of classes that cannot be assigned to the same room at the same time period). Such incompatibilities are usefully embodied through the structure of a graph. Each object (e.g. class) is represented by a node and each incompatibility is represented by an edge joining the two nodes. A coloring of this graph is then simply a partitioning of the objects into blocks (or colors) such that no two incompatible objects end up in the same block. Thus, optimal solutions to such problems may be found by determining minimal coloring for the corresponding graphs. Unfortunately, this may not always be accomplished in a reasonable amount of time.

8.2 Mapping Scheduling onto GC

Again, the mapping between a simple case of class scheduling and a graph-based representation goes as follows: The events are represented by the vertices of the graph, and a pair of vertices are joined by an indirected edge if and only if the corresponding events cannot take place at the same time. Scheduling the events subject to the given constraints is therefor equivalent to coloring the corresponding graph such that no two adjacent vertices are of the same color. The determination of the minimum number of intervals of time needed for the schedule is therefore the same as finding the minimum number of colors required for the graph. This is known as the chromatic number of the graph, and its determination for an arbitrary graph is yet unsolved problem.

Let us look at a simple example to illustrate the mapping pictorially. Suppose we have a graph with vertex set V and edge set E, where the ordered pair (R,S) is in E if and only if an edge exists between the vertices R and S. Two vertices in a graph are said to be adjacent if an edge exists between them. Given these sets the overall GC problem is to partition the vertices into a minimum number of sets in such a way that no two adjacent vertices are placed in the same set. Then, a different color is assigned to each set of vertices, as shown in Figure 8.1. This is known as a proper graph coloring.

In regard to the scheduling problem, our objective is an attempt in finding a schedule satisfying all of the hard constraints and if possible, most of the soft constraints. One essential subset of the hard constraints set is the time-based constraints, and the way to map them onto Figure 8.1 is as follows:

Let each course be represented by a vertex. Two vertices are connected by an edge if and only if there is a reason that the courses they represent may not be offered at

Figure 8.1: A proper graph coloring using three colors.

the same time. Initially, there are two such reasons for vertices to be linked: either the courses they represent are taught by the same instructor, or required by the same set of students. Upon adding in the links between vertices, we assign colors to the graph. Each color represents a different time slot, so every vertex with the same color will be offered at the same time.

Now assume we have, for example, four physics courses: PHY 101, PHY 234, PHY 300, and PHY 300 LAB. The first three which are lectures are taught by the same professor, and the fourth which is a laboratory is taught by another professor. In addition, students enrolled in the lecture section of PHY 300 need also to be enrolled in the lab section of the same course, and vice versa. These time-based hard constraints are represented in the mapping as shown in Figure 8.2. An edge exists between each of the pairs of courses taught by the same professor, since these courses cannot be offered at the same time. Also, an edge exists between the nodes of PHY 300 and PHY 300 LAB, because the same set of students must enroll in both courses. Figure 8.2 also shows a proper coloring with each node of the graph representing a different time slot, so only three different time slots are needed to properly schedule these four courses. PHY 101 and PHY 300 LAB will be offered during the same time slot, PHY 234 will be in a second time slot, and PHY 300 in a third time slot. Note that this system only partitions the courses into groups which will be offered at the same time; it does not assign the actual time of day. This assignment of groups to actual hours of the day can be done in another phase of the scheduling process in which soft constraints are handled.

Figure 8.2: A graph representing temporal hard constraints.

8.3 Graph Coloring Algorithm

When dealing with this kind of mapping between course scheduling and graph coloring, the first algorithm that a researcher would investigate would be the simple saturation algorithm (DSatur, for short) of Brelaz (ref: Brelaz-79). As we mentioned in the previous section that the timetabling problem need to be simple in order for methods such as DSatur to be effective. Before outlining this algorithm, we would like to briefly state one of the characteristics of the problem we have dealt with to show our justification or reasoning for using this algorithm. In addition, we also believe that our observations apply to other similar instance of the problems.

It is often the case that after carrying out the mapping, the representative graphs tend to be loosely connected. One particular reason for this that any college professor cannot be assigned to teach more than few courses per semester, perhaps maximum of four or five. Therefore, graphs sometimes would show fully-connected clumps of vertices representing those courses. That is partly the reason of having fully (or almost fully) connected components within the overall graph. Also, these components are usually loosely connected among themselves, resulting in a somewhat sparse graphs.

Graph coloring algorithms such as DSatur have been shown to work particularly quite well on sparse graphs. Figure 8.3 and Figure 8.4 show two different versions of the same algorithm. At the time when we used this technique, in the summer of 95, we tweaked with it and made two modifications to accommodate our problem setup.

1. One of the main constraints of course scheduling is the fixed number of available time slots per week, and the preference to distribute the courses among each of these time slots so that students have many options when selecting courses. Therefore, we modified DSatur to select colors for vertices so that, among the previously used legal colors for a vertex in V, the one used the fewest number of times previously was selected. This small modification to the algorithm served to distribute the courses more evenly over the time slots. In a typical graph coloring setup, a user usually tries to minimize the number of colors used.

2. The second modification to DSatur is that, due to our use of two sets of timeslots (lecture and laboratory timeslots) and also due to our criterion of scheduling one set prior to the other, we had to modify the method to color the lecture-of-course set of vertices first, followed by the lab-of-course set of vertices. Without this modification there is no distinction between various sets of vertices. Also, during the scheduling process we needed to do “re-coloring” of the graph that was constructed to get a legal schedule.

It is interesting to notice that the DSatur method is an example of a heuristic procedure for general graphs which is derived from an exact method designed for coloring bipartite graphs (for these we have |F (x)| ≤ 1 at each step). There are other coloring methods that slightly differ from DSatur in their choice of the next node to be colored as well as the color which it will receive.

Let V denote the vertex set, and E denotes the set of edges of the graph.

• Repeat until all vertices are colored

1. Choose the next vertex v from V to be colored by selecting the one whose neighbors have already used the most colors. Break ties by choosing the one with the most uncolored neighbors.

2. Color vertex v with the first color in the set of all colors that have not been already used by one of the neighbors, x of V, such that v and x are not equal.

• End Repeat.

Figure 8.3: DSatur Graph Coloring Algorithm – version 1

Let V denote the vertex set, and E denotes the set of edges of the graph.

1. Initially take any node with the largest possible number of neighbors and color it with the smallest color.

2. In general, for each v of V, there is a set F(v) of forbidden colors (these may be colors which have already been assigned to neighbors of v).

3. At each step, we choose a node x of V for which the cardinality, —F(x)—, of F(x) is maximum (highest connectivity with the neighbors), and we color it with the smallest possible color.

Figure 8.4: DSatur Graph Coloring Algorithm – version 2 CHAPTER 8. SCHEDULING AS A GRAPH COLORING PROBLEM 98

8.4 Some Experimental Results

— More on the work done in 1995 using a small set of data (about 450 science and engineering courses from Fall 93 at SU) here ... —

Method Comparisons and

Complexity Measures

9.1 Complexity of Course Scheduling

9.2 Airline Crew Scheduling

Airline crew scheduling is classified as a resource allocation problem, where a given flight schedule is to be covered by a set of crew rotations. Each rotation consists of a a connected sequence of flights or legs, each starting and ending at a given home base or hub. Then the objective is to minimize the total waiting time of the crew subject to a number of constraints on the rotations. Its topological structure and the restrictions imposed on it is quite similar to a multi-task phone routing structure.

One particular approach in tackling this problem is first to convert it into a set covering problem, by

1. generating a large number of legal rotation templates, and

2. seeking a subset of these templates that precisely covers the entire flight schedule.

CHAPTER 9. METHOD COMPARISONS AND COMPLEXITY MEASURES 100

Then, solutions to the set covering problem are found with conventional methods such as linear programming or other methods such as feedback neural networks.

This will have a big disadvantage and that is, on one hand, from a computational point of view, an exhaustive generation of rotations is not feasible for large real world problems; on the other hand, a non-exhaustive generation will only cover a fraction of the solution space.

One other approach followed by Peterson et. al (Neural Computation-98), proceeds in two phases:

1. the full solution space is narrowed down by using a reduction method that removes a large part of the sub-optimal solutions; then

2. a mean field annealing approach is applied to handle the topology, leg-counting, etc.

It has been determined that the computational requirement for random artificial problems with resemblance to real-world situations grows in the order of N_{f}^{3} where N_{f} is the number of flights.

So, given a schedule in terms of a set of N_{f} flights (let say per week), with specified times and airport departure and arrival, a crew is to be assigned to each flight such that the total crew waiting time is minimized subject to the following restrictions:

• Each flight crew must follow a connected path, or rotation, starting and ending at the hub.

• The number of flight legs in a rotation must not exceed a given upper bound, i.e. L_{max}.

• The total duration (flight + waiting time) of a rotation is similarly bounded by a maximum time, i.e. T_{max}.

In a real-world case there are many constraints to be satisfied, but in general the above three are the most crucial and difficult ones. Furthermore, without them the problem would reduce to that of minimizing waiting times independently at each airport — the local problems; which can be solved exactly in polynomial time; e.g. by pairwise connecting arrival flights with subsequent departure flights. Therefore, without the global structural requirements, the crew scheduling problem is not much of a challenge.

9.2.1 Crew Scheduling vs. Course Scheduling

The most noticeable difference between academic course scheduling and crew scheduling is in the topological structures. Course scheduling problem has a non-Euclidean global structure while the structure (or the graphical representation of it) of crew scheduling is quite similar to a multi-task phone routing structure. Neural networks based methods such as the abovementioned mean field annealing approach has been shown by several authors to be unsuited for tackling non-Euclidean problems. (ref: Peterson90).

The other distinction is in the duration. Typically, a real-world flight schedule has a basic period of about one week; while the duration of a university course schedule is an academic semester or a quarter.

Some of Geoffrey’s suggestions to be taken into account in the thesis.

• A section (medium to large) on the relationship between airline crew scheduling and course scheduling.

• A section (medium to large) giving the complexity of the problem. This would include runtime of various algorithms used, number of variables, size of the overall problem and data, etc.

• A chapter explaining why some optimization problems can also be tackled by methods such as used in this thesis, particularly the annealing approaches... While others are not suitable. For example, Genetic Algorithms are known to be quite good in tackling certain problems but not so good in tackling others. It would be quite useful to shed some light on the reasons behind this. Also, it would be interesting to have a section on the reasons why the neural net approach does not do well on scheduling problems.

• A section on the comparison between the type of constraints used in university course scheduling and those typically used for airline or crew scheduling.

Glossary and Definitions

Definitions of some of the key terms appearing in this thesis.

• Fitness Function: A fitness function is a mathematical equation that describes the relative value of each member of a family of objects with respect to some criterion. The fitness function is maximized or minimized (depending upon the criterion) to find the best object in the family. For example, one could determine which set of operating characteristics maximized the (revenue – operating expenses) criterion.

• Landscape Analysis: Landscape analysis is a mathematical process used to search for the best solution to a multi-variable optimization problem. It utilized a fitness function. Within the search space (i.e., a plot of all the possible solutions), the fitness of a particular solution can be shown graphically as its height. The resulting “landscape” will have peaks in the regions of the search space that contain better solutions; further analysis of the peak regions can be used to refine the results and achieve even more precise solutions.

• Multi-Objective Optimization: Multi-objective optimization is the task of maximizing or minimizing a number of criteria simultaneously. Multi-objective optimization is more difficult than most kinds of optimization problems because the trade-offs between the multiple objectives may be unknown.

• Simulated Annealing: A mathematical technique for general optimization problems. The same comes from the physical process of annealing, during which a material is first heated and then slowly cooled. During annealing, the component atoms of a material are allowed to settle into a lower energy state so that a more stable arrangement of atoms is maintained throughout the cooling process. SA has been applied with success to a wide variety of optimization tasks.

• Spin Glass: Spin glasses are simplified models from physics describing dilute magnetic alloys. They define an energy function (or equivalently a fitness function) over tiny magnetic dipoles called spins. Spins may point either up (+1) or down (-1). For every configuration of spins (e.g., +1-1+1+1-1 for a configuration of 5 spins), the energy function assigns that configuration’s energy by summing up the interactions between spins. many aspects of spin glasses are utilized by physicists and they comprise a well-understood family of fitness landscapes. etc. etc.

Bibliography

[AD92] R. W. Ashford and R. C. Daniel. Some lessons in solving practical integer programs. Journal of the Operational Research Society, 43(5):425–433, 1992.

[AK89] Emile H.L. Aarts and Jan Korst. Simulated Annealing and Boltzmann Machines: A Stochastic Approach to Combinatorial Optimization and Neural Computing. Wiley, Chichester, 1989.

[AL87] Emile H. L. Aarts and Peter J. M. Van Laarhoven. Simulated Annealing: Theory and Applications. Kluwer Academic Publishers, 1987.

[AMO93] Ravindra K. Ahuja, Thomas L. Magnanti, and James B. Orlin. Network Flows : Theory, Algorithms, and Applications. Prentice Hall, Englewood Cliffs, NJ, 1993.

[Bab92] László Babai. Combinatorial optimization is hard. Focus, 12(5):3–18, September 1992.

[BBI92] Robert E. Bixby, E. Andrew Boyd, and Ronni R. Indovina. MIPLIB: A test set of mixed integer programming problems. SIAM News, 25(20):16, March 1992.

[Bea90] J. E. Beasley. OR-library: Distributing test problems by electronic mail. Journal of the Operational Research Society, 41(11):1069–1072, 1990.

[Ber91] Dimitri P. Bertsekas. Linear Network Optimization: Algorithms and Codes. MIT Press, Cambridge, MA, 1991.

[BJN^{+}93] Cynthia Barnhart, Ellis L. Johnson, George L. Nemhauser, Gabriele Sigismondi, and Pamela Vance. Formulating a mixed integer programming problem to improve solvability. Operations Research, 41(6):1013–1019, 1993.

[Bla93] Jacek Blazewicz. Scheduling in computer and manufacturing systems. Springer-Verlag, New York, 1993.

[CLR90] Thomas H. Cormen, Charles E. Leiserson, and Ronald L. Rivest. Introduction to Algorithms. MIT Press, Cambridge, MA, 1990.

[Cof76] E. G. Coffman. Computer and Job Shop Scheduling Theory. Wiley, New York, 1976.

[Dav87] Lawrence Davis. Genetic Algorithms and Simulated Annealing. Morgan Kaufmann, Los Altos, CA, 1987.

[Dav91] Lawrence Davis. Handbook of Genetic Algorithms. Van Nostrand Reinhold, New York, 1991.

[Fis81] Marshall L. Fisher. The lagrangian relaxation method for solving integer programming problems. Management Science, 27(1):1–18, January 1981.

[Fou84] L. R. Foulds. Combinatorial Optimization for Undergraduates. Springer Verlag, New York, 1984.

[GJ79] Michael R. Garey and David S. Johnson. Computers and Intractibility: A Guide to the Theory of NP-Completeness. W. H. Freeman and Company, New York, 1979.

[GKP92] Fred Glover, Darwin Klingman, and Nancy Phillips. Network Models in Optimization and Their Applications in Practice. Wiley, New York, 1992.

[Glo89] Fred Glover. Tabu search– part I. ORSA Journal on Computing, 1(3):190–206, 1989.

[Glo90a] Fred Glover. Tabu search– part II. ORSA Journal on Computing, 2(1):4–32, 1990.

[Glo90b] Fred Glover. Tabu search: A tutorial. Interfaces, 20(4):74–94, 1990.

[GLS88] Martin Grötschel, László Lovász, and Alexander Schrijver. Geometric Algorithms and Combinatorial Optimization. Springer-Verlag, New York, 1988.

[Gol89] David E. Goldberg. Genetic Algorithms in Search, Optimization, and Machine Learning. Addison–Wesley, Reading, MA, 1989.

[HS84] Ellis Horowitz and Sartaj Sahni. Fundamentals of Computer Algorithms. Computer Science Press, Rockville, MD, 1984.

[Hu82] T. C. Hu. Combinatorial Algorithms. Addison-Wesley, Reading, MA, 1982.

[Joh92] David S. Johnson. The NP-completeness column: An ongoing guide. Journal of Algorithms, 13(3):502–524, September 1992.

[Jun85] Michael Junger. Polyhedral Combinatorics and the Acyclic Subdigraph Problem. Helderman Verlag, Berlin, 1985.

[KLS91] Bernhard Korte, László Lovász, and Rainer Schrader. Greedoids. Springer-Verlag, New York, 1991.

[Law76] Eugene Lawler. Combinatorial Optimization: Networks and Matroids. Saunders College Publishing, Fort Worth, 1976.

[Law85] Eugene Lawler. The Traveling Salesman Problem: A Guided Tour of Combinatorial Optimization. Wiley, New York, 1985.

[MT90] Silvano Martello and Paolo Toth. Knapsack Problems: Algorithms and Computer Implementations. Wiley, New York, 1990.

[NW88] George L. Nemhauser and Laurence A. Wolsey. Integer and Combinatorial Optimization. John Wiley and Sons, New York, 1988.

[OLK85] M. O’hEigeartaigh, J. K. Lenstra, and A. H. G. Rinnooy Kan. Combinatorial Optimization: Annotated Bibliographies. Wiley, New York, 1985.

[Oxl92] J. G. Oxley. Matroid Theory. Oxford University Press, New York, 1992.

[Pap94] Christos H. Papadimtriou. Computational Complexity. Addison-Wesley, Reading, MA, 1994.

[PR88] R. Gary Parker and Ronald L. Rardin. Discrete Optimization. Academic Press, Boston, 1988.

[PR91] Manfred Padberg and Giovanni Rinaldi. A branch–and–cut algorithm for the resolution of large–scale symmetric travelling salesman problems. SIAM Review, 33(1):60–100, 1991.

[PS82] Christos H. Papadimtriou and Kenneth Steiglitz. Combinatorial Optimization: Algorithms and Complexity. Prentice-Hall, Englewood Cliffs, NJ, 1982.

[Ree93] Colin R. Reeves. Modern Heuristic Techniques for Combinatorial Problems. Halsted Press, New York, 1993.

[Rei91] Gerhard Reinelt. TSPLIB— a traveling salesman problem library. ORSA Journal on Computing, 3(4):376–384, 1991.

[Sch87] Alexander Schrijver. Theory of Linear and Integer Programming. John Wiley and Sons, New York, 1987.

[Sha79] Jeremy F. Shapiro. A survey of lagrangian techniques for discrete optimization. 5:113–138, 1979.

[SM89] Harvey M. Salkin and Kamlesh Mathur. Foundations of Integer Programming. North-Holland, New York, 1989.

[Sto92] Mechthild Stoer. Design of Survivable Networks, volume 1531 of Lecture Notes in Mathematics. Springer Verlag, Heidelberg, 1992.

[Tar83] Robert Endre Tarjan. Data Structures and Network Algorithms. SIAM, Philadelphia, PA, 1983.

[Wil85] H. P. Williams. Model Building in Mathematical Programming. Wiley, New York, second edition, 1985.