Studieson
Metaheuristics for
Jobshop and Flowshop Scheduling Problems

Takeshi YAMADA

Studieson
Metaheuristics for
Jobshop and Flowshop Scheduling Problems

Takeshi YAMADA
Submitted in partial fulfillment of
the requirement for the degree of
DOCTOR OF INFORMATICS(Applied Mathematics and Physics)

KYOTO UNIVERSITY
KYOTO, JAPANNovember, 2003

Preface
Scheduling has been a subject of a significant amount of literature in the operations research field
since the early 1950s. The main objective of scheduling is an efficient allocation of shared resources over time to competing activities. Emphasis has been on investigating machine scheduling problems where jobs represent activities and machines represent resources. The problem is
not only NP-hard, but also has a well-earned reputation of being one of the most computationally difficult combinatorial optimization problems considered to date. This intractability is oneof the reasons why the problem has been so widely studied. The problem was initially tackled by“exact methods” such as the branch and bound method (BAB), which is based on the exhaustive
enumeration of a restricted region of solutions containing exact optimal solutions. Exact methods are theoretically important and have been successfully applied to benchmark problems, butsometimes they are quite time consuming even for moderate-scale problems.
With a rapid progress in computer technology, it has become even more important to findpractically acceptable solutions by “approximation methods” especially for large-scale problems within a limited amount of time. Stochastic local search methods are such approximation
methods for combinatorial optimization. They provide robust approaches to obtain high-quality
solutions to problems of realistic sizes in reasonable amount of time. Some of stochastic local
search methods are proposed in analogies with the processes in nature, such as statistical physics
and biological evolution, and others are proposed in the artificial intelligence contexts. They
often work as an iterative master process that guides and modifies the operations of subordinate
heuristics; thus they are also called metaheuristics. Metaheuristics have been applied to wide
variety of combinatorial optimization problems with great successes.
The primary focus of this thesis is applications of metaheuristics, especially Genetic Algorithms (GAs), Simulated Annealing (SA) and Tabu Search (TS), to the jobshop schedulingproblem (and the flowshop scheduling problem as its special case) which is among the hardest
combinatorial optimization problems. The author hopes that the research in this dissertation will
help advance in the understanding of this significant field.
November, 2003
Takeshi YamadaiAcknowledgements
I wish to express my sincere gratitude to Professor Toshihide Ibaraki of Kyoto University for his
supervising this thesis. He read the manuscript very carefully and made many valuable suggestions and comments, which improved the accuracy and quality of this thesis. I also thank Professor Masao Fukushima, Professor Yutaka Takahashi, Professor Hiroyuki Kawano, Professor
Mutsunori Yagiura and Professor Koji Nonobe of Kyoto University for their useful comments.
The research reported in this thesis was supported by Nippon Telegraph and Telephone Corporation (NTT). I am grateful to Professor Seishi Nishikawa, Professor Tsukasa Kawaoka, Dr.
Kohichi Matsuda, Professor Yoh’ichi Tohkura, Professor Kenichiro Ishii, former directors of
NTT Communication Science Laboratories, Dr. Noboru Sugamura and Dr. Shigeru Katagiri,
director and vice director of NTT Communication Science Laboratories, for their warm encouragement and for providing me the opportunity to study these interesting subjects.
I am deeply indebted to Professor Ryohei Nakano of Nagoya Institute of Technology. He had
been my supervisor for more than ten years since I first started my research career at NTT fifteen
years ago. This thesis would not have been possible without his support and encouragement. I
am also indebted to Professor Colin Reeves of Coventry University. Some of the work have been
done while I was working with him as a visiting researcher at the university in 1996. I wish to
express my many thanks to Professor Bruce Rosen of University of California. The collaboration
with him while he was visiting NTT in 1994 is very important especially for the early stage of
the work.
I am also grateful to Dr.Ueda and Dr.Saito of NTT Communication Science Laboratories fortheir encouragement and long standing friendship.
Finally, I thank my parents for their endless support and encouragement, and my wife Kazumi
to whom I dedicate this work, for everything else.

iiContents
1

Introduction
1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Outline of the Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1
1
3

2

The Job Shop Scheduling Problem
2.1 The Problem Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Active Schedules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Disjunctive Graph Representation . . . . . . . . . . . . . . . . . . . . . . . . .
2.4 DG Distance and Binary Representation . . . . . . . . . . . . . . . . . . . . . .
2.5 Block Property . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6 The Shifting Bottleneck Heuristic . . . . . . . . . . . . . . . . . . . . . . . . .
2.7 The One-machine Scheduling Problem . . . . . . . . . . . . . . . . . . . . . . .
2.8 The Well-known Benchmark Problems . . . . . . . . . . . . . . . . . . . . . . .

5
5
7
13
15
16
18
20
25

3

Genetic Algorithms
3.1 Basic Concepts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 A Simple Genetic Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 The Procedure of a Simple Genetic Algorithm . . . . . . . . . . . . . . . . . . .

30
30
31
33

4

A Simple Genetic Algorithm for the Jobshop Scheduling Problem
4.1 Genetic Encoding of a Solution Schedule . . . . . . . . . . . . . . . . . . . . .
4.2 Local harmonization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Global harmonization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4 Forcing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5 Simple GA for the JSP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.6 The Limitation of the Simple Approach . . . . . . . . . . . . . . . . . . . . . .

35
35
36
38
38
40
40

5

GT-GA: A Genetic Algorithm based on the GT Algorithm
5.1 Subsequence Exchange Crossover . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Precedence Preservative Crossover . . . . . . . . . . . . . . . . . . . . . . . . .
5.3 GT Crossover . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4 GT-GA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.5 Computational Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42
43
43
45
48
48

iii

CONTENTS

iv
5.6

Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

6Neighborhood Search
6.1 The Concept of the Neighborhood Search . . . . . . . . . . . . . . . . . . . . .
6.2 Avoiding Local Optima . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6.3 The Neighborhood Structure for the Jobshop Scheduling Problem . . . . . . . .

52
52
54
54

7

Critical Block Simulated Annealing for the Jobshop Scheduling Problem
7.1 Simulated Annealing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.2 Critical block Simulated Annealing . . . . . . . . . . . . . . . . . . . . . . . . .
7.3 Reintensification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.4 Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.5 Methodology and Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.5.1 Random Search . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7.5.2 Low Temperature Greedy Search . . . . . . . . . . . . . . . . . . . . .
7.6 Performance on Benchmarks Problems . . . . . . . . . . . . . . . . . . . . . . .
7.7 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

57
57
59
61
61
62
64
65
67
70

8

Critical Block Simulated Annealing with Shifting Bottleneck Heuristics
8.1 Active Critical Block Simulated Annealing . . . . . . . . . . . . . . . . . . . .
8.2 Active CBSA Enhanced by Shifting Bottleneck . . . . . . . . . . . . . . . . . .
8.3 Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.3.1 Muth and Thompson’s Benchmark . . . . . . . . . . . . . . . . . . . . .
8.3.2 Other Benchmarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
8.4 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

71
71
72
76
76
76
78

9

Scheduling by Genetic Local Search with Multi-Step Crossover Fusion
9.1 Multi-step crossover fusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9.2 Scheduling in the reversed order . . . . . . . . . . . . . . . . . . . . . . . . . .
9.3 MSXF-GA for Job-shop scheduling . . . . . . . . . . . . . . . . . . . . . . . .
9.4 Benchmark Problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9.4.1 Muth and Thompson benchmark . . . . . . . . . . . . . . . . . . . . . .
9.4.2 The Ten Tough Benchmark Problems . . . . . . . . . . . . . . . . . . .

79
79
82
83
83
85
85

10 Permutation Flowshop Scheduling by Genetic Local Search
10.1 The Neighborhood Structure of the FSP . . . . . . . . . . . . . . . . . . . . . .
10.2 Representative Neighborhood . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10.3 Distance Measures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10.4 Landscape analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10.5 MSXF-GA for PFSP . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10.6 Experimental results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10.7 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

89
89
91
92
92
95
96
98

CONTENTS

v

11 C sum Permutation Flowshop Scheduling by Genetic Local Search99
11.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
11.2 Representative Neighborhood . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
11.3 Tabu List Style Adaptive Memory . . . . . . . . . . . . . . . . . . . . . . . . . 100
11.4 Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
11.5 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
12 Tabu Search with a Pruning Pattern List for the Flowshop Scheduling Problem
103
12.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
12.2 Tabu Search . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
12.3 Pruning Pattern . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
12.4 Pruning Pattern List Approach . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
12.5 Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
12.6 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
13 Conclusions

110

A List of Author’s Work118

List of Figures
The job sequence matrix {T jk } and the processing time matrix {p jk } for the 3 ×
3 problem given in Table2.1. T jk = r means that k-th operation for job J j isprocessed on machine Mr for p jk time units. . . . . . . . . . . . . . . . . . . . .
2.2 A Gantt chart representation of a solution for the 3 × 3 problem given in Table 2.1. Operation O31 can be shifted as early as at 5 time unit, as indicated bydotted lines, without altering the order of operations on any machine, and the
new schedule becomes semi-active. . . . . . . . . . . . . . . . . . . . . . . . .
2.3 The solution matrix S rk for the solution given in Figure2.2. S rk = j means thatthe k-th operation on machine Mr is job J j . . . . . . . . . . . . . . . . . . . . . .
2.4 An example of a permissible left shift, where in the upper picture, O12 can beshifted to the front of O32 without delaying any other operations resulted in amuch improved schedule given in the lower picture. . . . . . . . . . . . . . . . .
2.5 A snapshot in the middle of scheduling using Giffler and Thompson’s algorithm,
where O11 , O22 , O31 , O43 , O51 and O64 are schedulable. O11 is the earliest completable operation, and O11 and O31 are in conflict with O11 . O31 is selected forthe next operation to be scheduled and then O11 and O51 must be shifted forward
to avoid overlapping. O31 , which is the next operation to O31 in the technologicalsequence, now becomes schedulable. . . . . . . . . . . . . . . . . . . . . . . .
2.6 A disjunctive graph representation of a 3 × 3 problem . . . . . . . . . . . . . . .
2.7 The DG distance between two schedules: the distance = 2. . . . . . . . . . . . .
2.8 Labeling disjunctive arcs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.9 An example in which there are three critical blocks is illustrated. The blocks B1 ,B2 and B3 are on the critical path P(S ) and are corresponding to the different machines Mr1 , Mr2 and Mr3 respectively. The adjacent critical blocks are connectedby a conjunctive arc. The last operation on machine Mr3 is not a critical block,which must contain at least two operations. . . . . . . . . . . . . . . . . . . . .
2.10 The condition that no operation in a block in S is processed before the first or
after the last operation of the block in S 0 implies that all the operations in Bi areprocessed prior to those in B j both in S and S 0 if i < j, because each adjacent
blocks are connected by a conjunctive arc that cannot be altered. Hence, there
is a path P(S 0 ) in S 0 that contains all the operations on P(S ) and the length ofP(S 0 ) is greater than the length of P(S ). . . . . . . . . . . . . . . . . . . . . . .
2.11 A Schrage schedule represented by a conjunctive graph . . . . . . . . . . . . . .
2.1

vi

6

8
8

9

12
13
15
15

17

18
22

LIST OF FIGURES

vii

3.1

An example of the roulette wheel selection, where the roulette wheel is createdaccording to the fitness value of each individual shown in the upper left picture . 33

4.1
4.2

An optimal schedule for the mt06 (6 × 6) problem (makespan = 55) . . . . . . .
A binary representation of a solution schedule using the job-based ordering corresponding to the solution given in Figure 4.1. The first line corresponds to theprecedence relation between J1 and J2 . The first three digits of the bit-string on
the first line are 110. This corresponds to the fact that J1 is processed prior to J2
on J1 ’s first and second machines M3 and M1 , but is not prior to J2 on J1 ’s thirdmachine M2 and so on. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
An example of the local harmonization resolving cycles within six operations
O11 , O21 , . . . , O61 on the same machine M1 where the arcs O31 → O61 , O21 → O41
and O11 → O61 are reversed in this order and a consistent ordering O41 → O61 →O51 → O11 → O21 → O31 is eventually obtained. . . . . . . . . . . . . . . . . .
An example of global harmonization where a cycle O23 → O22 → O32 → O31 →O33 → O23 is resolved by reversing an arc O22 →O32 . . . . . . . . . . . . . . . .

4.3

4.4
5.1

5.2

5.3
5.4
5.5
5.6
5.7
6.1
6.2

The solution given in Figure 2.3 is converted to an m-partitioned permutation form = 3, where the permutation in the k-th partition corresponds to the processingorder of jobs on machine Mk . . . . . . . . . . . . . . . . . . . . . . . . . . . .
An example of subsequence exchange crossover (SXX), where each underlined
subsequence pair one from p0 and the other from p1 on each machine is identified
as exchangeable and interchanged to generate k0 and k1 . . . . . . . . . . . . . .
A job sequence (permutation with repetition) for a 3 × 3 problem defined in
Figure 2.1 is decoded to a schedule, which is equivalent to the one in Figure 2.3. .
An example of the precedence preservative crossover (PPX), where k is generated from p0 and p1 using h . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
GT crossover . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
The histgram of the best makespans obtained by the GT-GA after 200 generations
among 600 trials for the mt10 problem . . . . . . . . . . . . . . . . . . . . . . .
relationship between various crossover operators . . . . . . . . . . . . . . . . .

37

37

39
39

43

44
44
45
47
50
51

AE(S ), adjacent exchange neighborhood of S , consists of schedules obtained
from S by exchanging a pair of adjacent operations within a same critical block. . 55
CB(S ), critical block neighborhood of S , consists of schedules obtained from S
by moving an operation in a critical block to the front or the rear of the block. . . 55

7.1
7.2

Generated Makespans of 10,000 Greedy (mt10) Schedules. . . . . . . . . . . . . 66
Successive makespan differences between the current and optimal solution of the
mt10 problem, without reintensification (R=0) and with reintensification (R =3, 000). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

8.1

The time evolution of CBSA+SB trial for the la27 problem . . . . . . . . . . . . 77

LIST OF FIGURES

viii

A simple 2 × 2 problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Schedule reversal and activation . . . . . . . . . . . . . . . . . . . . . . . . . .
Distribution of solutions generated by an application of (a) MSXF and (b) ashort-term stochastic local search . . . . . . . . . . . . . . . . . . . . . . . . . .
Performance comparison using the la38 15 × 15 problem . . . . . . . . . . . . .

82
82

10.1 A grid graph representation of a solution to a problem of 8 jobs and 6 machines. .
10.2 The best move to the next/previous block is selected as a representative. . . . . .
10.3 1841 distinct local optima obtained from 2500 short term local search for theta011 (20 × 10) problem and 2313 distinct local optima for the ta021 (20 × 20)
problem are plotted in terms of (a) average distance from other local optima and(b) distance from global optima (x-axis), against their relative objective function
values (y-axis). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10.4 The correlation between the precedence-based distance (PREC) and the approximate number of steps (STEPS) . . . . . . . . . . . . . . . . . . . . . . . . . . .
10.5 The correlation between the precedence-based distance (PREC) and the positionbased distance (POSN) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10.6 Navigated local search by MSXF-GA: A new search is started from one of theparents and while no other good solutions are found, the search ‘navigates’ towards the other parent. In the middle of the search, good solutions would beeventually found somewhere between the parents. That direction is then pursued
to the top of a hill (or a bottom of the valley, if it is a minimization problem) —a new local optimum. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

90
91

9.1
9.2
9.3
9.4

84
87

93
94
95

97

11.1 Representative neighborhood . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
11.2 The framework of the proposed method . . . . . . . . . . . . . . . . . . . . . . 102
12.1 When v = (2, 7) is applied to π, (π(2), π(3)) = (x, y) is stored in T as tabu. Later,v0 = (1, 6) is not allowed to apply to β because it will restore the previously
banned precedence relation between (x, y). . . . . . . . . . . . . . . . . . . . . . 104
12.2 The time evolutions of makespans for the ta041 (50 jobs and 10 machines) problem averaged over 30 tabu search runs with and without the pruning pattern list
(left). The time evolutions for the ta051 (50 jobs and 20 machines) problem averaged over 10 tabu search runs and the computationally equivalent MSXF-GA
runs for comparison (right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
12.3 The time evolutions for the other nine Taillard problems of 50 jobs and 20 machines (ta052 – ta060) averaged over 10 tabu search runs with (labeled TS+PL)
and without (labeled TS) the pruning pattern list. . . . . . . . . . . . . . . . . . 109

List of Tables
2.1

2.2
2.3
2.4
2.5
2.6

An example of the jobshop scheduling problem with 3 jobs and 3 machines. Eachcolumn represents the technological sequence of machines for each job with theprocessing times in parentheses. . . . . . . . . . . . . . . . . . . . . . . . . . .
An example of the one-machine scheduling problem with 7 jobs . . . . . . . . .
Muth and Thompson’s 6 × 6 problem (mt06) . . . . . . . . . . . . . . . . . . . .
Muth and Thompson’s 10 × 10 problem (mt10) . . . . . . . . . . . . . . . . . .
Muth and Thompson’s 20 × 5 problem . . . . . . . . . . . . . . . . . . . . . . .
The ten tough benchmark problems (status reported by [1] in 1991) . . . . . . . .

4.1

Experimental results of the simple GA for mt benchmark problems . . . . . . . . 40

5.1

Experimental results of the GT-GA for mt benchmark problems . . . . . . . . . 49

7.1
7.2
7.3
7.4
7.5
7.6

Ten Trials using the Simulated Annealing Method (R = 3, 000). . . . . . . . . . .
Initial and Last Temperatures. Last temperature is the temperature when an optimal makespan was found, or the temperature after 1,000,000 iterations. . . . . .
Ten High Temperature Random Trials. . . . . . . . . . . . . . . . . . . . . . . .
Comparisons between CBSA and AESA. . . . . . . . . . . . . . . . . . . . . .
Ten difficult Benchmark Job Shop Scheduling Problems. . . . . . . . . . . . . .
Performances of the 40 Benchmark Job Shop Scheduling Problems. . . . . . . .

8.1
8.2
8.3

Comparisons between CBSA and CBSA+SB using MT benchmarks . . . . . . . 76
Results of 10 tough JSPs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
An optimal solution of la27 problem . . . . . . . . . . . . . . . . . . . . . . . . 78

9.1
9.2
9.3

Performance comparison using the MT benchmark problems . . . . . . . . . . . 85
Results of the 10 tough problems . . . . . . . . . . . . . . . . . . . . . . . . . . 86
Performance comparisons with various heuristic methods on the 10 tough problems 87

6
22
27
27
28
29

63
64
65
67
68
69

10.1 Results of the Taillard benchmark problems . . . . . . . . . . . . . . . . . . . . 97
11.1 Taillard’s benchmark results (ta001 – ta030) . . . . . . . . . . . . . . . . . . . . 101
11.2 Taillard’s benchmark results (ta031 – ta040) . . . . . . . . . . . . . . . . . . . . 102

ix

Chapter 1
Introduction
1.1 BackgroundIn late 1950s, B.Giffler and G.L.Thompson first showed in their paper titled “Algorithms forsolving production scheduling problems” [2] that it is not necessary to search for an optimalschedule over all possible schedules, but only over a subset of the feasible schedules, called
active schedules. They then proposed the GT algorithm, which is described later in Section 2.2,
to iteratively enumerate all active schedules for a given problem. H.Fisher and G.L.Thompson
proposed three well-known benchmark problems known as mt06 (6 jobs 6 machines), mt10 (10jobs 10 machines) and mt20 (20 jobs 5 machines) [3] in a book titled “Industrial Scheduling”
edited by J.F.Muth and G.L.Thompson in 1963 [4]. The “notorious” mt10 problem has been
unsolved for over 20 years. Their paper is also important in the sense that they first applied astochastic approximation method based on priority dispatching rules and active schedules.
As a pioneering work in the jobshop scheduling research, G.H.Brooks and C.R.White first
proposed a branch and bound method, a tree-based exact method to solve the problem optimally,
based on the GT algorithm [5]. E.Balas first pointed out the fact that the adjacent pairwise
exchange of operations on a certain part of the schedule, called “critical path”, of a feasible
schedule will always results in another new feasible schedule [6]. This fact will later play animportant role in metaheuristics context.
A great deal of efforts by Barker and McMahon [7] and then Carlier [8] among others havecontributed the progress of the exact approaches, which are mainly based on the branch andbound method. They have commonly used the mt benchmark problems (especially the mt10
problem) as a computational challenge to demonstrate the effectiveness of their algorithms, andthe best known solution for the mt10 problem has been improved. Finally in 1985, Carlier
and Pinson succeeded in solving the mt10 problem optimally by a branch and bound algorithm[8]. Since then, Brucker [9], Martin and Shmoys [10], and Carlier again [11] have improved theperformance of exact approaches. However, the NP-hardness of the problem barriers the efficient
application of exact methods to larger-scale problems.
In addtion to those exact methods, many approximation methods have been developed. Simulated annealing (SA) is one of the well-known stochastic local search methods, based on an1

2

Chapter 1. Introduction

analogy with the physical process of annealingl; heating up a solid in a heat bath until it melts,
then cooling it down slowly until it solidifies into a low-energy state results in a pure lattice
structure. Laarhoven et al. proposed a simulated annealing algorithm for the jobshop schedulingproblem, using the pairwise exchange of operations on the critical path proposed by Balas, asa transition operator [12]. However, very similar idea had already been proposed by Matsuo etal. [13].
Adams et al. proposed a very powerful method to find reasonably efficient schedules known
as shifting bottleneck heuristic in 1988 [14]. This method, as its name suggests, iteratively
identifies a bottleneck machine and optimize its job sequence. The details of the algorithm aredescribed in Section 2.6. In 1991, Applegate have combined the “shifting bottleneck heuristic”
and a branch and bound method to develop a powerful algorithm and demonstrated that the mt10
problem is no more a computational challenge. Instead, they proposed a new set of benchmark
problems known as the “ten tough benchmarks”, which contained the ten difficult problems including seven open problems that were not solved even by their approach [1].
In 1990s, Tabu Search (TS), proposed by Fred Glover [15, 16], has been used by many
researchers including Taillard [17], Dell’Amico, Trubian [18], Nowicki and Smutnicki [19, 20].
TS adopts a deterministic local search, whereby a ‘memory’ is implemented by the recording ofpreviously-seen solutions. Instead of storing solutions explicitly, this record is often an implicit
one in a sense that it stores the moves, or the modifications of the solution, that have been
made in the recent past of the search, and which are ‘tabu’ or forbidden for a certain number of
iterations. This prevents cycling, and also helps to promote a diversified coverage of the search
space. Taillard also proposed a new benchmark set consisting of 80 jobshop and 120 flowshop
problems known as “Taillard benchmark” [21].
Genetic Algorithms (GAs) model biological processes to optimize highly complex objectivefunctions. They allow a population composed of many individuals to evolve under specified
selection rules to a state that maximizes the “fitness”. The method was developed by John Holland over the course of the 1960s and 1970s [22], and popularized by one of his students, David
Goldberg who successfully applied a GA to the control of gas-pipeline transmission. He is also
well-known for a book titled “Genetic Algorithms in Search, Optimization, and Machine Learning” [23].
GAs have been used successfully in various fields of computer science, including machinelearning, control theory and combinatorial optimization. GAs can be uniquely characterized bytheir population-based search strategies and their operators: mutation, selection and crossover.
Nakano and Yamada were among the first who applied a conventional GA that uses binary representation of solutions, to the jobshop scheduling problem [24]. Yamada and Nakano [25]proposed a GA that uses problem-specific representation of solutions with crossover and mutation, which are based on the GT algorithm. The details of these approaches are described later in
Chapters 4 and 5.
Ulder and others first proposed Genetic Local Search (GLS), which is a hybridization of
GAs and local search [26]. In this framework, each individual, or search agent, performs local
search independently, while crossover is performed occasionally to the solutions of two selectedindividuals and a new solution is produced, which is then used as an initial solution for thesubsequent local search performed by an offspring. In this context, the embedded local search is a

1.2. Outline of the Thesis

3

main search engine to effectively improve solutions and crossover provides information exchange
between individuals who are performing independent local search in parallel.

1.2 Outline of the Thesis
This thesis is devoted to jobshop and flowshop scheduling by metaheuristics, especially by Genetic Algorithms, Simulated Annealing and Tabu Search. In Chapter 2, besic concepts and notations of the jobshop scheduling problem are described such as the semi-active and active schedules, the disjunctive graph representation and critical path and blocks. The main focus throughout
this thesis is the minimum-makespan problem, in which makespan, maximum completion timeof all the operations, is used as an objective function to be minimized. This is denoted as Cmax .
The sum of completion times of all operations, denoted as C sum , is also considered as an alternative objective function of the flowshop scheduling problem. The GT (Giffler & Thompson’s)
algorithm for generating active schedules and the well-known shifting bottleneck heuristic thatgenerates moderately good schedules by repeatedly solving one-machine scheduling problems
are also reviewed as well as some well-known benchmark problems.
In Chapter 3, Genetic Algorithms are reviewed with a major emphasis on conventional binarymodels for combinatorial optimization, in which a solution is encoded into a binary string of fixed
length and binary genetic operators, such as one-point, two-point and uniform crossover and bitflip mutation, are used. In Chapter 4, a conventional GA using a binary representation is applied
to the jobshop scheduling problem. By converting a solution of the problem into a bit-string,
conventional GAs can be applied without major modification.
In Chapter 5, a more advanced GA approach is described as the GT-GA method, which involves a non-binary representation of a solution schedule and domain knowledge, namely, active
schedules and the GT algorithm. GT-GA method consists of GT crossover and GT mutationthat are defined as simple modifications of the GT algorithm. One of the advantages of the GA
is its robustness over a wide range of problems with no requirement of domain specific adaptations. From this point of view, the conventional GA approach with binary encoding and binary
crossover that is obviously domain independent is reasonable. However, it is often more desireble to directly incorporate problem specific knowledge such as the GT algorithm into GA,
resulting the GT-GA method, from the performance point of view.
In Chapter 6, the concept of neighborhood search is described as a widely used local searchtechnique to solve combinatorial optimization problems and is extended to include metaheuristics. Especially, it is shown that Simulated Annealing (SA) and Tabu Search (TS) can be considered as advanced meta strategies for neighorhood search to avoid local optima. An efficientneighborhood for the jobshop scheduling problem, called Critical Block (CB) neighborhood, thatis defined based on the critical path and blocks, is also described.
In Chapter 7, a SA method for the jobshop scheduling problem that utilizes the CB neighborhood is described, and then in Chapter 8, it is further extended by incorporating the shifting bottleneck heuristic, which can be regarded as a problem specific deterministic local searchmethod.
In Chapter 9, it is shown that Genetic Algorithms (GAs) can be reagarded as a variant of

4

Chapter 1. Introduction

neighorhood search, that is called Genetic Local Search (GLS), and an approach called MultiStep Crossover Fusion (MSXF) method is proposed. In the MSXF method, one of the parent
solutions is used as an initial point of the new local search, while the other is used to define anorientation for the search. In other words, it is a local search that traces out a path from one
solution to another. The MSXF is applied to the jobshop scheduling problem in Chapter 9 andapplied to the Cmax and C sum flowshop scheduling problems in Chapter 10 and in Chapter 11
respectively.
In Chapter 12, a TS method with a data structure called the Pruning Pattern List (PPL) for theCmax flowshop scheduling problem is described. A pruning pattern is constructed from a solutionof the flowshop scheduling problem represented by a permutation of jobs numbers by replacingsome of its job numbers by a wild card. A list of pruning patterns generated from good schedules
collected in the course of a search process is used to inhibit the search to visit already searchedand no longer interesting region again and again and it is embedded into a TS method.
Finally, in Chapter 13, the study in this thesis is summerized and the future directions aresuggested.

Chapter 2
The Job Shop Scheduling Problem
Scheduling is the allocation of shared resources over time to competing activities. It is convenient to adopt manufacturing terminology, where jobs represent activities and machines represent
resources, while the range of application areas for scheduling theory is not limited to computers
and manufacturing but includes transportation, services, etc. In this chapter, we restrict our attention to deterministic jobshop scheduling, where all the data that define a problem instance are
known in advance.

2.1 The Problem Description
The n×m minimum-makespan general jobshop scheduling problem, designated by the symbols
n/m/G/Cmax and hereafter referred to as the JSP, can be described by a set of n jobs {Ji }1≤ j≤nwhich is to be processed on a set of m machines {Mr }1≤r≤m . The problem can be characterized asfollows:
1. Each job must be processed on each machine in the order given in a pre-defined technological sequence of machines.
2. Each machine can process only one job at a time.
3. The processing of job J j on machine Mr is called the operation O jr .
4. Operation O jr requires the exclusive use of Mr for an uninterrupted duration p jr , its processing time; the preemption is not allowed.
5. The starting time and the completion time of an operation O jr is denoted as s jr and c jr
respectively. A schedule is a set of completion times for each operation {c jr }1≤ j≤n,1≤r≤m thatsatisfies above constraints.
6. The time required to complete all the jobs is called the makespan, which is denoted as
Cmax . By definition, Cmax = max1≤ j≤n,1≤r≤m c jr .
5

Chapter 2. The Job Shop Scheduling Problem

6

The problem is “general”, hence the symbol G is used, in the sense that the technological sequence of machines can be different for each job as implied in the first condition and that theorder of jobs to be processed on a machine can be also different for each machine. The predefined technological sequence of each job can be given collectively as a matrix {T jk } in which
T jk = r corresponds to the k-th operation O jr of job Ji on machine Mr . The objective of optimizing the problem is to find a schedule that minimizes Cmax .
An example of a 3 × 3 JSP is given in Table 2.1. The data include the technological sequenceof machines for each job with the processing times in parentheses. In the table, operations
for job 1, for example, are processed in the order of O11 → O12 → O13 ; i.e., job 1 is firstprocessed on machine 1 with its processing time 3, and then processed on machine 2 with its
processing time 3, and then processed on machine 3 with its processing time 3. The problem
is equivalently represented by the job sequence matrix {T jk } and processing time matrix {p jk } asgiven in Figure 2.1.
Table 2.1: An example of the jobshop scheduling problem with 3 jobs and 3 machines. Each column represents the technological sequence of machines for eachjob with the processing times in parentheses.
job
1
2
3

1 (3)
1 (2)
2 (3)

machine (processing time)
2 (3)
3 (3)
3 (3)
2 (4)
1 (2)
3 (1)

1 2 3
{T jk } = 1 3 2 ,
2 1 3

3 3 3
{p jk } = 2 3 4
3 2 1

Figure 2.1: The job sequence matrix {T jk } and the processing time matrix {p jk } forthe 3 × 3 problem given in Table2.1. T jk = r means that k-th operation for job J j isprocessed on machine Mr for p jk time units.

Instead of minimizing the makespan Cmax , other objectives such as minimizing the sum ofP
the completion times of all the operations C sum (C sum = 1≤ j≤n,1≤r≤m c jr ) can also be considered
and is designated by the symbols n/m/G/C sum .
The Gantt chart is a convenient way of visually representing a solution of the JSP. The Gantt
chart shows time units at the abscissa and machine numbers at the axis of ordinate. An example
of a solution for the 3 × 3 problem in Table 2.1 is given in Figure 2.2. In the figure, each squarebox represents an operation O jr with its left edge placed at s jr as its x coordinate and with its
horizontal length representing processing time p jr . The makespan of this schedule is Cmax = 19
time unit.

2.2. Active Schedules

7

A schedule is called semi-active when no operation can be started earlier without altering theoperation sequences on any machine. Operation O31 in Figure 2.2, for example, can be started as
early as at 5 time unit, as indicated by dotted lines, without altering the order of operations on any
machine, and the new schedule is semi-active. By definition, a semi-active schedule is uniquely
obtained by specifying the operation sequences for all machines. In other words, an semi-active
schedule is represented by a m × n matrix S = {S rk } where S rk = j corresponds that the k-th
operation on Mr is job J j . Figure 2.3 shows the matrix representation of a solution for the 3 × 3
problem in Table 2.1. In the figure, operations on machine 2, for example, are processed in theorder of O23 → O22 → O12 . Given a matrix S = {S rk }, it is straightforward to obtain an associatedsemi-active schedule. Each operation O has two predecessor operations: job predecessor andmachine predecessor. The job predecessor of O, denoted by PJ(O), is the direct predecessor of
O in the technological sequence. The machine predecessor of O, denoted by PM(O), is the direct
predecessor of O in the solution matrix. For example, if we have a problem as given in Table 2.1,
we have PJ(O12 ) = O11 and if the solution is given as in Figure 2.3, then PM(O12 ) = O22 . An
operation O is scheduled at time unit 0 if it has no job and no machine predecessors. If only oneof job and machine predecessors exists, then O is scheduled immediately when the predecessor
is completed. Otherwise it is scheduled when both job and machine predecessors are completed
such that: s(O) = max{c(PJ(O)), c(PM(O))}, where s(O) and c(O) are the starting time andcompletion time of operation O respectively. The solution given in Figure 2.3 corresponds to theGannt chart representation given in Figure 2.2 except that O31 is started at 5 time unit.
Throughout this thesis, we assume that a schedule is always semi-active if not stated otherwise. By this formulation, jobshop scheduling can be interpreted as defining the ordering
between all operations that must be processed on the same machine, i.e., to fix precedences between these operations. In short, the jobshop scheduling problem is formulated as an ordering
problem.
As a special case, when the technological sequence of machines is the same for all jobs
and the order in which each machine processes the jobs is also same for all machines, then aschedule is uniquely represented by a permutation of jobs. This simplified problem is called the
permutation flowshop scheduling problem and it is designated by the symbols n/m/P/Cmax (orn/m/P/C sum when the objective is minimizing C sum ).

2.2 Active Schedules
The makespan of a semi-active schedule may often be reduced by shifting an operation to the left
without delaying other jobs. Consider a semi-active schedule S and two operations O jr and Okr
in S that use the same machine Mr . If Okr is processed prior to O jr and the machine Mr has anidle period longer than p jr before processing Okr , then reassignment is possible so that operation
O jr is processed prior to Okr without delaying any other operations. Such reassignment is calleda permissible left shift and a schedule with no more permissible left shifts is called an active
schedule. Figure 2.4 shows an example of permissible left shift. The schedule in the upper pictureof Figure 2.4 is identical to the one given in Figure 2.3 but O31 which is started at 5 time unit andits makespan = 19. On machine M2 in the schedule, O12 , the operation of job J1 , can be shifted toChapter 2. The Job Shop Scheduling Problem

8

M1

O11

M2

O32

O21

O31
O22

M3

O12

O23
0

2

4

6

O13
8

10

12

14

16

O33
18

time

Figure 2.2: A Gantt chart representation of a solution for the 3 × 3 problem given
in Table 2.1. Operation O31 can be shifted as early as at 5 time unit, as indicated bydotted lines, without altering the order of operations on any machine, and the newschedule becomes semi-active.

1 2 3{S rk } = 3 2 1
2 1 3
Figure 2.3: The solution matrix S rk for the solution given in Figure2.2. S rk = jmeans that the k-th operation on machine Mr is job J j .

2.2. Active Schedules

9

the front of O32 , the operation of J3 , without delaying any other operations. Operation O13 is thenshifted to the point immediately after the completion of its job and machine predecessors O12 andO23 . O33 on machine 3 is also shifted. The resulting schedule with improved makespan = 12 isshown in the lower picture of Figure 2.4. Because there always exists an optimal schedule that is
active, it should be safe and efficient to restrict the search space to the set of all active schedules.

M1

O11

M2

O32

O21

O31
O22

M3

O12

O23
0

2

4

6

M1

O11

O21

O31

M2

O32

O12

M3

O13
8

2

4

12

14

16

18

time

O22
O23

0

10

O33

6

O13
8

10

O33
12 timeFigure 2.4: An example of a permissible left shift, where in the upper picture, O12
can be shifted to the front of O32 without delaying any other operations resulted in amuch improved schedule given in the lower picture.

An active schedule can be generated by using the GT algorithm proposed by Giffler andThompson [2]. The algorithm is described in Algorithm 2.2.1. In the algorithm, the following
notations are used:• As in the previous section, an operation O has job and machine predecessors. The job
predecessor denoted by PJ(O) is the direct predecessor of O in the technological sequence.
The definition of the machine predecessor PM(O) is slightly modified and defined as thelast scheduled operation on the machine, that is, an operation with the largest completiontime among already scheduled operations on the same machine as O.
• An operation O that is not yet scheduled is called schedulable when both its job predecessor
PJ(O) and machine predecessor PM(O), if they exist, have already been scheduled. The
set of all schedulable operations is denoted as G.
• The earliest starting time ES (O) of O is defined as the maximum of the completion times
of PJ(O) and PM(O): ES (O) := max{c(PJ(O)), c(PM(O))}, where c(O) is the completion

Chapter 2. The Job Shop Scheduling Problem

10time of O. The earliest completion time EC(O) is defined as ES (O) plus its processing
time p(O).
• The earliest completable operation O?r in D with its machine Mr , is an operation whoseearliest completion time EC(O?r ) is the smallest in D (break ties arbitrarily) :O?r = arg min{EC(O) | O ∈ D}.

(2.1)• Given an earliest completable operation O?r and if there are i − 1 operations that havealready been scheduled on Mr , a conflict set C[Mr , i] is a set of candidate operations forthe next (i.e. i-th) processing on Mr defined as:
C[Mr , i] = {Okr ∈ D | ES (Okr ) < EC(O?r )}.

(2.2)

Note that O?r ∈ C[Mr , i].
The essence of GT algorithm is scheduling operations while avoiding an idle period thatis long enough to allow a permissible left shift. For this purpose, one has to carefully selectthe next operation that does not introduce such a long idle period among the set of schedulable
operations D. Hence the conflict set C[Mr , i] ⊂ D is maintained. As long as the next operation is
selected from the conflict set, an idle period is kept sufficiently short and the resulting scheduleis guaranteed to be active.
An active schedule is obtained by repeating Algorithm 2.2.1 until all operations are scheduled. In Step 4, instead of choosing one operation randomly, if all possible choices are considered, all active schedules will be generated, but the total number will still be very large. Practically, random choice is replaced by the application of so-called priority dispatching rules [27],which are the most popular and simplest heuristics for solving the JSP. For example, a rule called
SOT (shortest operation time) or SPT (shortest processing time) selects an operation with theshortest processing time from the conflict set, a rule called MWKR (most work remaining) selects an operation associated with the job with the longest total processing time remaining, a rulecalled FCFS (first come first serve rule) selects the first available operation among operations on
the same machine. Dorndorf and Pesch [28] proposed a priority rule-based GA for the JSP using
twelve such priority dispatching rules.
In the notations above, if the definition of the conflict set C[Mr , i] is simplified as C[Mr , i] =Gr , then the generated schedule is an semi-active schedule. Otherwise, if the definition of
C[Mr , i] is altered as C[Mr , i] = {Okr ∈ Gr | ES (Okr ) < ES (O?r )}, then the generated schedule is called a non-delay schedule. Unlike active schedules, an optimal schedule is not always anon-delay schedule.
Figure 2.5 shows how the GT algorithm works. The figure shows a snapshot in the middleof scheduling where operations O11 , O22 , O31 , O43 , O51 and O64 are schedulable and constitute G.
The earliest completable operation is identified as O11 , which results in G1 = {O11 , O31 , O51 }.
In G1 , only O11 and O31 satisfy the inequality in (2.2), therefore C[M1 , i] = {O11 , O31 }. If O31
is selected from C[M1 , i], then O11 and O51 are shifted forward according to Step 6 in Algorithm 2.2.1.

2.2. Active Schedules

11

Algorithm 2.2.1 (GT algorithm)
A scheduling problem represented by {T jk }, the technological sequence matrix, and {p jk }, theprocessing time matrix is given as an input.
1. Initialize G as a set of operations that are first in the technological sequence; i.e., G ={O1T11 , O2T21 , . . . , O2Tn1 }. For each operation O ∈ G, set ES (O) := 0 and EC(O) := p(O).
2. Find the earliest completable operation O∗r ∈ G by (2.1) with machine Mr . A subset of G
that consists of operations processed on machine Mr is denoted as Gr .
3. Calculate the conflict set C[Mr , i] ⊂ Gr by (2.2), where i−1 is the number of operations
already scheduled on Mr .
4. Select one of operations in C[Mr , i] randomly. Let the selected operation be Okr .
5. Schedule Okr as the i-th operation on Mr ; i.e. S ri := k, with its starting and completion
times equal to ES (Okr ) and EC(Okr ) respectively: s(Okr ) = ES (Okr ), c(Okr ) = E(COkr ).
6. For all O jr ∈ Gr \ {Okr }, Update ES (O jr ) as
ES (O jr ) := max{ES (O jr ), EC(Okr )} and EC(O jr ) as EC(Okr ) := ES (Okr ) + p(Okr ).
7. Remove Okr from G (and therefore from Gr ), and add an operation Oks that is the next to
Okr in the technological sequence to G if such Oks exits; i.e., if r = T ki and i < m, thens := T ki+1 and G := (G \ {Okr }) ∪ {Oks }.
Calculate ES (Oks ) and EC(Oks ) as:
ES (Oks ) := max{EC(Okr ), EC(PM(Oks ))} and EC(Oks ) := ES (Oks ) + p(Oks ) respectively.
8. Repeat from Step 1 to Step 7 until all operations are scheduled.
9. Output the solution matrix {S rk } as the active schedule obtained with the set of starting andcompletion times {s(O jr )} and {c(O jr )} respectively where j = S rk .

Chapter 2. The Job Shop Scheduling Problem

12G
C[M1,i]

O11
O22

O11
O22

O31

O31

O43

O32

O43
O51
O64

selected

O51
O64

G1
Figure 2.5: A snapshot in the middle of scheduling using Giffler and Thompson’salgorithm, where O11 , O22 , O31 , O43 , O51 and O64 are schedulable. O11 is the earliest
completable operation, and O11 and O31 are in conflict with O11 . O31 is selected forthe next operation to be scheduled and then O11 and O51 must be shifted forward
to avoid overlapping. O31 , which is the next operation to O31 in the technologicalsequence, now becomes schedulable.

2.3. Disjunctive Graph Representation

13

2.3 Disjunctive Graph Representation
The Gantt chart representation and the matrix representation described in the previous section
are simple and straightforward to identify a schedule. However it is not obvious to see whether
the resulting schedule is feasible or not: i.e., whether the job sequence on each machine does notcontradict with the pre-defined technological sequence of machines. A more informative problem
formulation based on a graph representation is first introduced by Roy and Sussman [29]. In thissection, we review a graph representation for the JSP using disjunctive graph formulation. The
following descriptions and notations are due to Adams et. al. [14].
The JSP can be described by a disjunctive graph G = (V, C ∪ D), where• V is a set of nodes representing operations of the jobs together with two special nodes, asource (0) and a sink ?, representing the beginning and end of the schedule, respectively.
• C is a set of conjunctive arcs representing technological sequences of machines for eachjob.
S
• D = mr=1 Dr , where Dr is a set of disjunctive arcs representing pairs of operations thatmust be performed on the same machine Mr .
The processing time for each operation is the weighted value pv attached to the corresponding
node v, and for the special nodes, p0 = p∗ = 0. Figure 2.6 shows a disjunctive graph representation of the problem given in Table 2.1.
conjunctive arc (technological sequences)disjunctive arc (pair of operations on the same machine)

sourcep11= 3O11
p21 = 4

0

O21

O32
p32= 3p12= 3O12
p23= 3O23

O31
p31= 2

p13 = 3O13
p22 = 2
O22

sink

*

O33
p33= 1

Oij : an operation of job i on machine jp ij : processing time of Oij

Figure 2.6: A disjunctive graph representation of a 3 × 3 problem
Let sv be the starting time of an operation corresponding to node v. By using the disjunctive
graph notation, the jobshop scheduling problem can be formulated as a mathematical programming model as follows:
minimize: s∗
subject to:sw − sv ≥ p v ,(v, w) ∈ C,
sv ≥ 0,v ∈ V,
sw − sv ≥ pv ∨ sv − sw ≥ pw , (v, w) ∈ Dr , 1 ≤ r ≤ m.

(2.3)

14

Chapter 2. The Job Shop Scheduling Problem

The formula A ∨ B means that either A or B is to be satisfied (but not both), thus it is called
disjunctive constraint. Note that ∗ is the dummy sink node that has a zero processing time. Thismeans that s∗ is equal to the completion time of the very last operation of the schedule, whichis therefore equal to Cmax . The first inequality in (2.3) ensures that when there is a conjunctive
arc from a node v to a node w, w must wait at least pv time period after v is started, thus the predefined technological sequence of machines for each job is not violated. The second condition isequivalent to s0 ≥ 0. According to the third constraints, when there is a disjunctive arc between anode v and a node w, one has to select either v to be processed prior to w (and w waits for at least
pv time period) or the other way around. This avoids any pair of operations on the same machine
to overlap in time. In the disjunctive graph, the selection corresponds to fixing the undirected(disjunctive) arc into a directed one.
To summarize, jobshop scheduling is to define the ordering between all operations that must
be processed on the same machine, as described in the previous section. This corresponds to thethird constraints in (2.3), and this is done by fixing all undirected (disjunctive) arcs into directed
ones: thus the disjunctive graph is turned into a conjunctive graph. A selection is defined as a setof directed arcs selected from the set of disjunctive arcs D. By definition, a selection is complete
if all the disjunctions in D are selected. It is consistent if the resulting directed graph is acyclic.
When a complete selection is consistent, one can define a unique and consistent ordering ofoperations on the same machine, namely a solution matrix {S rk } and this matrix corresponds to afeasible (semi-active) schedule. Hence a consistent complete selection, the obtained conjunctive
graph, and the corresponding (semi-active) schedule are all identified and represented by thesame symbol S without confusion. Given a selection, a path starting from a node v to any
destination node w is defined by following directed arcs from v to w (if they exist), and the lengthof the path is defined as the sum of the weights of all the nodes on the path including v and w.
It is clear that the makespan Cmax of a schedule S is given by the length of the longest weighted
path from source to sink in the graph of the corresponding complete selection. This path P (notnecessarily unique) is called a critical path and is composed of a sequence of critical operations.
By using the disjunctive graph model, we can easily show a well-known fact known as feasibility property for adjacent exchanges of operations on a critical path: any exchange of two
adjacent operations on a critical path will never lead to an infeasible schedule. Based on thisfact, Laarhoven et al. have proposed a simulated annealing algorithm using the pairwise exchange of operations on the critical path as a transition operator [12]. Taillard has proposed a
Tabu Search method by using the same transition operator [17].
Theorem 1 (feasibility for adjacent exchange) Let S be a consistent complete selection andP(S ) be a critical path in S . Consider a pair of adjacent critical operations (u, v) on a samemachine on P(S ), i.e.,there is an arc selected from u to v. Then a complete selection S uv obtained from S by reversing the direction of the arc between u and v is always acyclic (thus thecorresponding schedule is always feasible).
Proof: Assume the contrary, then the exchange introduces a cycle in S uv . This means that there
is a path Pu,v from u to v in S uv , and this Pu,v also exists in S . P(S ) can be represented asP(S ) = (0, . . . , t, u, v, w, . . . , ∗) and (0, . . . , t, Pu,v , w, . . . , ∗) is also a path from source to sink in2.4. DG Distance and Binary Representation

15

S but clearly longer than P(S ). This contradicts the assumption of the theorem that P(S ) is a
critical path of S .


2.4 DG Distance and Binary Representation
The distance between two schedules S and T can be measured by the number of differences in
the processing order of operations on each machine [24]. In other words, it can be calculated
by counting the directed (originally disjunctive) arcs whose directions are different between S
and T . We call this distance the disjunctive graph (DG) distance. Figure 2.7 shows the DG
distance between two schedules. The two directed arcs marked by thick lines in schedule T havedirections that differ from those of schedule S , and therefore the DG distance between S and T
is 2.
O11

S
0

O13

O12

O22

O23

O21

O32

O13

O12

O22

O23

O21

0

*

O33

O31

O11

T

O32

*

O33

O31

Figure 2.7: The DG distance between two schedules: the distance = 2.
As described in the previous section, a (semi-active) schedule is obtained by turning all undirected disjunctive arcs into directed ones. Therefore, by labeling each directed (originally) disjunctive arc of a schedule as 0 or 1 according to its direction, and rearrange them as a one dimensional vector, a schedule can be represented by a binary string of length mn(n − 1)/2. Figure 2.8
shows a labeling example, where an arc connecting Oi j and Ok j (i < k) is labeled as 1 if the arc is
directed from Oi j to Ok j (so Oi j is processed prior to Ok j ) or 0, otherwise. It should be noted thatthe DG distance between schedules and the Hamming distance between the corresponding binary strings can be identified through this binary mapping. Nakano and Yamada have proposed a
simple Genetic Algorithm based on this binary coding and using standard genetic operators [24].
O11

O12

1 1
0

O21

1

O23

1
0
O32

0

O31

0

O13

1
O22

1

*

O33

Figure 2.8: Labeling disjunctive arcs
By using the notion of DG distance, the following so called connectivity property for adjacentexchanges on the critical path can be derived easily from Theorem 1 as follows:Chapter 2. The Job Shop Scheduling Problem

16

Theorem 2 (connectivity for adjacent exchange) Let S and P(S ) be an arbitrarily scheduleand its critical path, then it is possible to construct a finite sequence of adjacent exchange on thecritical path that will lead to an optimal schedule.
Proof: Let S ∗ be an optimal schedule. Because S is not optimal, Cmax (S ) > Cmax (S ∗ ), and the
DG distance, denoted by d, between S and S ∗ is d(S , S ∗ ) > 0. Moreover, there is at least one pairof consecutive critical operations (u, v) in S such that (u, v) is processed in S in this order but inS ∗ , v is processed prior to u. This is true because if such pair does not exist, then the critical path
P(S ) in S exists also in S ∗ as a path and this means that Cmax (S ) < Cmax (S ∗ ) which contradicts
the assumption that S ∗ is optimal. Reverse the direction of (u, v) and obtain S uv . Theorem 1
guarantees that S uv is feasible. It is clear that d(S uv , S ∗ ) = d(S , S ∗ ) − 1, i.e., S uv is closer to S ∗than S by one step. Replace S by S uv and repeat this process at most d(S , S ∗ ) times until there isno such (u, v), then S becomes identical with S ∗ , or at least the critical paths of S and S ∗ become
identical, therefore Cmax (S ) = Cmax (S ∗ ) and so S is optimal.


2.5 Block Property
As described in the previous section, an operation on a critical path is called a critical operation.
A sequence of more than one consecutive critical operations on the same machine is called acritical block. More formally, let S be a feasible schedule associated with a disjunctive graph
G(S ) = G(V, C ∪ D) with all the disjunctive arcs in D being directed. Let P(S ) be a critical pathin G(S ) and L(S ) be the length of this critical path, which is equal to the makespan. A sequenceof successive nodes in P is called a critical block or just block if the following two properties are
satisfied:
1. The sequence contains at least two nodes,
2. The sequence represents a maximal number of operations to be processed on the samemachine.
The j-th block on the critical path is denoted by B j . Figure 2.9 shows an example of critical
blocks on a critical path. The following so-called block property gives us crucial informationin improving a current schedule by simple modifications and thus forms a basis for many of thejobshop scheduling solvers [30].
Theorem 3 (Block property) Let S , P(S ), L(S ) be a complete selection, its critical path, andthe length of the critical path respectively. If there exists another complete selection S 0 such thatL(S 0 ) < L(S ), then at least one operation of some block B of G(S) has to be processed in S 0
before the first or after the last operation of B.
Proof: Assume the contrary that L(S 0 ) < L(S ) and that there is no operation that satisfies theconclusion of the theorem; i.e., there is no operation of any block of S that is processed before

2.5. Block Property

17

Mr

Mr

2

Mr

B1

B2

B3

1

3

Mr

4

S
0

*

Figure 2.9: An example in which there are three critical blocks is illustrated. The
blocks B1 , B2 and B3 are on the critical path P(S ) and are corresponding to thedifferent machines Mr1 , Mr2 and Mr3 respectively. The adjacent critical blocks areconnected by a conjunctive arc. The last operation on machine Mr3 is not a critical
block, which must contain at least two operations.

the first or after the last operation of the corresponding block. Then, there is a path P(S 0 ) fromsource to sink in S 0 that contains all the operations on P(S ):{P(S 0 )} ⊃ {P(S )}(2.4)

where {P(S )} denotes a set of all nodes on P(S ) (See Figure 2.10). Let l(P) be the length of P,
then from ( 2.4) we have:l(P(S 0 )) ≥ L(S )(2.5)

The definition of the critical path indicates:
L(S 0 ) ≥ l(P(S 0 ))(2.6)

L(S 0 ) ≥ L(S )(2.7)

From (2.5) and (2.6), we have:This contradicts the assumption of the theorem.
The theorem gives us an important criterion about how to improve a current schedule. Namely,
if we wish to improve a schedule S , then either one of the following two situations must happen:
1. At least one operation in one block B, that is not the first one in B, has to be processed
before all the other operations in B.
2. At least one operation in one block B, that is not the last one in B, has to be processed after
all the other operations in B.
Let S be a complete selection, A critical path P(S ) in G(S ) defines critical blocks B1 , . . . , Bk .
Roughly speaking, the before-candidates BBj is defined as a set of all but the first operations in B j
and after-candidates BAj a set of all but the last operations in B j . More precisely, the first and theChapter 2. The Job Shop Scheduling Problem

18

S'

B1

B2

B3

0

*

Figure 2.10: The condition that no operation in a block in S is processed before thefirst or after the last operation of the block in S 0 implies that all the operations in Bi
are processed prior to those in B j both in S and S 0 if i < j, because each adjacent
blocks are connected by a conjunctive arc that cannot be altered. Hence, there is apath P(S 0 ) in S 0 that contains all the operations on P(S ) and the length of P(S 0 ) is
greater than the length of P(S ).

last blocks B1 and Bk need special treatment. If the first operation u11 of the first block B1 is also
the very first operation of the critical path P(S ), then we set B1B as empty. Likewise, if the lastoperation ukmk of the last block Bk is also the very last operation of the critical path P(S ), then we
set BkA as empty.
(∅
if j = 1 and u11 is the first in P(S )
B
(2.8)
Bj =
B j \ {u1j } otherwise.
(∅
if j = k and ukmk is the last in P(S )
A
(2.9)
Bj =
j
B j \ {um j } otherwise.
where u1j and umj j are the first and the last operations in B j . Note that B j contains at least two
operations, and so BBj ∪ BAj is always non-empty.
Brucker et al. have proposed an efficient branch and bound method based on the blockproperty [9]. It is natural to consider generating a new schedule by moving an operation in BBj(or BAj ) to the front (or rear) of B j aiming for possible improvements. The adjacent exchange ofcritical operations discussed in the end of Section 2.3 is a special case of this kind of transition.
Many metaheuristics approaches have been proposed based on this transition operators [18, 19].
Note that unlike the simpler adjacent exchange of critical operations, in which the feasibility is
guaranteed by Theorem 1, applying this transition may result in an infeasible schedule.

2.6 The Shifting Bottleneck Heuristic
The shifting bottleneck heuristic (SB) proposed by Adams, Balas and Zawack [14] is one of
the most powerful heuristic methods for the jobshop scheduling problem. Assume we have apartial schedule or corresponding selection S p , in which only some of the machines are alreadyscheduled and the ordering of operations on those machines has been determined and fully fixed.

2.6. The Shifting Bottleneck Heuristic

19

For this partial schedule, a critical path is defined exactly the same way as in the complete case;the longest path from source to sink in the graph of the corresponding selection S p . Hence themakespan L(S p ) is also defined as the length of the critical path.
Given a partial schedule S p , then we focus on a machine Mk not yet scheduled. When weschedule operations on machine Mk , we have to take into accounts the constraints imposed by
the already scheduled operations on other machines. For example, assume that we have a problem given in Figure 2.6 as a disjunctive graph, and that operations on machine M1 and M3 arescheduled and their starting and completion times are fixed. Then, operation O12 for job J1 on
machine M2 , which is not yet scheduled, has to be started at least after the completion of O11
and ended before the start of O13 . In short we have to start and complete the processing of O12
within a given limited time interval. Similar constraints are imposed to O22 and O32 as well, andwe have to determine the order of O12 , O22 and O32 on machine M2 such that these constraints
are not violated.
In general, we define head and tail for each unscheduled operation. Let i be an operation
on machine Mk not yet scheduled. Operation i can be identified as a node i in the graph corresponding to selection S p . Let ri be the length of the longest path from source to the node i: i.e.,ri = L(0, i), where L(i, j) is the length of the longest path from node i to node j. ri is called therelease time or the head of operation i. In a similar fashion, let qi be the length of the longest
path from i to the sink minus processing time of i, i.e., qi = L(i, ∗) − pi , where pi is the processing
time of i. qi is called the tail of operation i. The due date di is defined as di = L(0, ∗) − qi . Notethat L(0, ∗) is the makespan of S p .
When the head ri , tail qi and the processing time pi , summarized as {ri , pi , qi } are given foreach operation i on a machine Mk , and let C as a set of operations on machine Mk , we have a one
machine scheduling problem formulated as a mathematical programming model as follows:
minimize the makespan:maxi∈C {si + pi + qi }subject to:si ≥ riand the disjunctive constraints: s j − si ≥ pi ∨ si − s j ≥ p j(i ∈ C)(i, j ∈ C)

(2.10)Starting from a partial schedule S p , which is initially set as empty, we solve a one-machinescheduling problem for each machine not yet scheduled to optimality, and find a bottleneckmachine: a machine with the longest one-machine makespan. The algorithm to solve the onemachine scheduling problem will be described in the next section. The bottleneck machine isregarded as scheduled and S p is updated using the job sequence on the bottleneck machine obtained above. Every time a new machine has been scheduled, the job sequence on each previously
scheduled machine is subject to reoptimization. The original SB consists of two subroutines: thefirst one (SBI) repeatedly solves one-machine scheduling problems; the second one (SBII) builds
a partial enumeration tree where each path from the root to a leaf is similar to an application of
SBI.
The rough outline of SBI can be summarized as follows:
1. For each of unscheduled machines, solve the one-machine scheduling problem and obtain
the best makespan and corresponding job sequence.

Chapter 2. The Job Shop Scheduling Problem

20

2. Identify the most bottleneck machine which has the longest one-machine makespan obtained above.
3. Make the most bottleneck machine scheduled using the job sequence obtained above.
4. Reoptimize all the scheduled machines.
5. Repeat the above until all the machines are scheduled.
The more complete SBI algorithm is given in Algorithm 2.6.1.
Instead of considering the most bottleneck machine in Step 3 of Algorithm 2.6.1, if we consider the n-th highest bottleneck machines and apply the remaining steps of SBI for each bottleneck candidate, we have SBII.

2.7 The One-machine Scheduling Problem
In the SBI heuristic, we have to repeatedly solve the one-machine scheduling problem. Although
the problem is NP-hard, Carlier has developed an efficient branch and bound method [31]. In thissection, we focus on the one-machine problem and describe the algorithm proposed by Carlier.
In the one-machine case, each job has only one operation, so an operation and corresponding job
is identified. Furthermore, a simplified notation is used to identify job number and corresponding
job; i.e., we just say “job i” instead of saying “operation Oi of job Ji ”. The disjunctive graph ofthe problem is defined just as a special case of the jobshop scheduling problem. A schedule isobtained by determining the starting (or completion ) times of all jobs, or equivalently, turning
all undirected disjunctive arcs into directed ones, resulting in a conjunctive graph, or simply, a
job sequence.
Consider a one-machine scheduling problem P with n jobs I = {1, 2, . . . , n} characterized by{ri , pi , qi }, where ri , pi and qi are the head, processing time and tail of each job i respectively.
Hereafter, we omit a set of n jobs I when it is obvious and just say “P is defined by {ri , pi , qi }”.
The formal definition of the one-machine scheduling problem was already given by (2.10) and we
do not repeat it here. Table 2.2 shows an example of the one-machine scheduling problem with
7 jobs and Figure 2.11 shows an example of a schedule represented by a conjunctive graph forthis problem. Note that the conjunctive graph is simplified such that only n−1 arcs are presented
between adjacent job nodes instead of drawing n(n − 1)/2 arcs between all job node pairs, which
are apparently redundant. The number on each arc from source to a job node i is the head ri
and the numbers on each arc from a job node i to sink is the processing time pi “+” the tail qi .
The schedule presented in the figure is called a Schrage schedule, the definition of which will bepresented shortly.
A lower bound of the makespan for a one-machine scheduling problem defined by I ={J1 , . . . , Jn } and {ri , pi , qi } is calculated as follows:
Theorem 4 Let I1 be a subset of I, thenh(I1 ) = Mini∈I1 ri +X
i∈I1

pi + Mini∈I1 qi(2.11)

2.7. The One-machine Scheduling Problem

21

Algorithm 2.6.1 (SBI heuristic)
A scheduling problem given as a set of technological sequences of machines for each job withprocessing times (preferably represented by a disjunctive graph) are given as inputs. Let M be
the set of all machines: M = {M1 , . . . , Mr }.
1. A partial schedule S p and a set of already scheduled machines M0 in S p are initialized as
S p := {} and M0 := {} respectively.
2. For each Mk ∈ M \ M0 , do the following(a) Compute head ri and tail qi for each operation i on Mk given S p . Let pi be theprocessing time of i.
(b) Solve one-machine scheduling problem {ri , qi , pi } to optimality for machine Mk andobtain the best one-machine makespan v(Mk ) and corresponding job sequence (moreprecisely, corresponding selection) S Mk .
3. Let M ∗ be the bottleneck machine such that v(M ∗ ) ≥ v(M) for all M ∈ M \ M0
4. Set M0 := M0 ∪ {M ∗ } and S p := S p ∪ S M∗ .
5. Order the elements of M0 as {M1 , M2 , . . . , Ml } (l = |M0 |) in the order of its inclusion toM0 above.
6. [local optimization of S p with already scheduled machines M0 = {M1 , M2 , . . . , Ml }]do(a) For m = 1, 2, . . . , l, do the following:i. Reset the sequence of operations on Mm , i.e., let S p 0 = S p \ S Mm .
ii. Compute head ri and tail qi for each operation i on Mm of S p 0 .
iii. Solve one-machine scheduling problem {ri , qi , pi } and obtain new v(Mm ) andS Mm .
iv. Let S p 00 = S p 0 ∪ S Mm and compute makespan L(S p 00 ) for partial schedule S p 00 . If
L(S p 00 ) < L(S p ), then set S p := S p 00 , otherwise keep the original S p .
(b) Reorder {M1 , M2 , . . . , Ml } according to decreasing order of the new v(Mm )s.
while any improvements are made in Step 6a.
7. Repeat from Step 2 to Step 5 until M0 = M.
8. Output S p as a obtained complete selection.

Chapter 2. The Job Shop Scheduling Problem

22

Table 2.2: An example of the one-machine scheduling problem with 7 jobs
i (= job no.)
1
2
3
4
5
6
7
ri (head)
10
13
11
20
30
0
30
pi (processing time)
5
6
7
4
3
6
2
qi (tail)
7
26
24
21
8
17
0

is a lower bound on the optimal makespan.
Proof: In the conjunctive graph associated with the optimal schedule, there is a path PI1 fromsource to sink, passing through every job in I1 . It is clear that the length of the path PI1 is longer
than or equal to h(I1 ). Whereas, by definition of the critical path, the length of the path PI1 isshorter than or equal to the length of the critical path and which is equal to the makespan of the
optimal schedule. Thus, h(I1 ) is a lower bound.

There is a method to generate a reasonably good schedule called the longest tail heuristic
in which a job with the longest tail qi , therefore with the earliest due date, among ready jobs,
is regarded as the most urgent and scheduled first. The resulting schedule is called a Schrageschedule. An algorithm to generate a Schrage schedule is described in Algorithm 2.7.1. The
name “longest tail” heuristic stems from (2.13).
Figure 2.11 shows a Schrage schedule obtained by the longest tail heuristic to the problem
given in Table 2.2. The starting times of each job is: s1 = 10, s2 = 15, s3 = 21, s4 = 28, s5 = 32,
s0 = s6 = 0, s7 = 35, s? = 53. The critical path is 0, 1, 2, 3, 4, ? and the makespan is equal to
53. This example is taken from [31].
J6
6

J1

0
10
13

5

J2

17

7

6+2
6

7 + 24

J3
7

20

J4
30
30

5+

6

11

0

6+

4

4+2

1

3+

8
3+

J5

*

0

3

J7

Figure 2.11: A Schrage schedule represented by a conjunctive graph2.7. The One-machine Scheduling Problem

23

Algorithm 2.7.1 (longest tail heuristic)
A one-machine scheduling problem defined by {ri , pi , qi } for each job i ∈ I = {1, 2, . . . , n} is given
as input
1. Initialize a set of scheduled jobs U := {}, and the most recently scheduled job k := 0.
The starting times and processing times of the two dummy nodes 0 and ? are defined as
s0 := s? := 0 and p0 := p? := 0 respectively.
2. Identify R, a set of ready jobs as R := {i ∈ I \ U | ri ≤ sk + pk }.
3. If R is empty, then k0 , the next job to be scheduled, is selected as
k0 := arg Mini∈I\U ri(2.12)and schedule k0 as sk0 = rk0 . Otherwise k0 is selected as
k0 := arg Maxi∈R qi
and schedule k0 as sk0 = sk + pk . Update k := k0 .
4. Update U := U ∪ {k} and update s? as s? := max{s? , sk + pk + qk }.
5. Repeat Step 2 to Step 4 until U = I.
6. Output the set of obtained starting times {s0 , s1 , . . . , sn , s? }.

(2.13)

Chapter 2. The Job Shop Scheduling Problem

24

The following theorem implies that the Schrage schedule is in a sense “close” to the optimalschedule, and to improve it, a particular job c has to be rescheduled.
Theorem 5 Let L be the makespan of the Schrage schedule.
(a) If this schedule is not optimal, there is a critical job c and a critical set J such that:
X
h(J) = Mini∈J ri +pi + Mini∈J qi > L − rci∈J(b) If this schedule is optimal, there exists J such that h(J) = L.
Proof: Let G be the directed graph associated with the Schrage schedule and consider a critical
path of G that contains maximal number of jobs. By renumbering jobs if necessary, the critical
path P is denoted as P = (0, 1, 2, . . . , l, ?) without loss of generality.
The length of the critical path L is:X
L = s1 +pi + ql .
(2.14)
i=1,...,lIf job 1 was scheduled immediately when its predecessor k, if exists, was just finished; i.e.,
if sk + pk equals to s1 , then k must be on the critical path, too. This is in contradiction with the
maximality of P. Hence sk + pk < s1 , and this means R was empty and 1 is selected by (2.12) in
Algorithm 2.7.1. Therefore,
s1 = r1 = Mini=1,...,l ri .

(2.15)

If ql = Mini=1,...,l qi , then L = h(J) with J = {1, . . . , l} from (2.14) and (2.15). Because h(J)is a lower bound from Theorem 4, the schedule is optimal.
Otherwise, there exists i ∈ {1, . . . , l} such that qi < ql . Let c be the greatest number among such iand let J = {c + 1, . . . , l}, then:qc < qk for all k ∈ J(2.16)

ql = Mink∈J qk .

(2.17)

andFrom (2.16), c has the “shortest” tail among {c} ∪ J: the least urgent job in the “longest” tailheuristic. This means that when c was scheduled, c was not selected by (2.13) but selected by(2.12), meaning:sc = rc < rk for all k ∈ J.

(2.18)

Because c is on the critical path and s1 = r1 from (2.15),
sc = s1 + p1 + . . . pc−1 = r1 + p1 + . . . pc−1 .

(2.19)

2.8. The Well-known Benchmark Problems

25

From (2.18) and (2.19), we have
r1 + p1 + . . . pc−1 < rk for all k ∈ J.

(2.20)

Together with (2.17) and (2.20), we have:h(J) = Mink∈J rk +X

pk + Mink∈J qk > r1 + p1 + . . . pc−1 + pc +k∈JX

pk + ql − pc .

(2.21)k∈J

The right hand side of (2.21) is equal to L − pcAs a corollary of this theorem, it can be seen that in an optimal schedule, either c is processed
before all the jobs in J or after all the jobs in J. By using this fact, we consider two new onemachine scheduling problems S L and S R by modifying head or tail of the original problem S .
Let P be a one-machine scheduling problem corresponding to {ri , pi , qi }. Apply the longest tailheuristic to P and obtain a Schrage schedule S and the critical job c and critical set J. PL requires
job c scheduled before all jobs in J. Thus PL is obtained from P using the same {ri , pi , qi } but
only modifying qc as follows:X
pr + ql ).
(2.22)
qc := Max(qc ,k∈J

Likewise PR requires job c scheduled after all jobs in J. Therefore, PR is obtained from P bymodifying rc as follows:X
pr ).
(2.23)
rc := Max(rc , Mink∈J rr +k∈J

Now we are ready to describe the branch and bound algorithm. In the branch and bound algorithm, we consider a search tree in which each node is associated with a one-machine schedulingproblem P defined by {ri , pi , qi }. The root node P0 corresponds to the original problem to besolved. The active node Pa is initialized as P0 . The upper bound µ is initialized as µ := ∞.
We apply the longest tail heuristic to the active node Pa and obtain the Schrage scheduleand its makespan L s (Pa ), the critical operation c and the critical set J. The upper bound µ is
updated as µ := min{µ, L s (Pa )}. Then PL and PR are generated by using (2.22) and (2.23).The
lower bound λ of the two new nodes will be λ(PL ) := max{L s (Pa ), hL (J), hL (J ∪ {c})} andλ(PR ) := max{L s (Pa ), hR (J), hR (J ∪ {c})} respectively, where hL and hR correspond to h in (2.11)
but calculated with modified {ri , pi , qi } by (2.22) and (2.23) respectively. A new node will be
added to the tree only if its lower bound is less than the upper bound µ. The next active node is
identified as a node with the lowest bound.

2.8 The Well-known Benchmark Problems
The three well-known benchmark problems with sizes of 6 × 6, 10 × 10 and 20 × 5 (known as
mt06, mt10 and mt20) formulated by Muth and Thompson [4] are commonly used as test beds toChapter 2. The Job Shop Scheduling Problem

26

Algorithm 2.7.2 (Carlier’s branch and bound algorithm)
A one-machine scheduling problem P0 defined by {ri , pi , qi } is given as input
1. Initialize Pa = P0 and µ := ∞.
2. Apply the longest tail heuristic to Pa and obtain the Schrage schedule and its makespan
L s (Pa ), the critical operation c and the critical set J.
3. If L s (Pa ) < µ, then store Pa as the best schedule obtained so far:Pb := Pa , and updateµ := L s (Pa ). Make Pa visited.
4. Generate PL and PR defined by {riL , piL , qiL } and {riR , pRi , qRi } respectively using (2.22) and(2.23).
5. Calculate λ(PL ), the lower bound for PL , as λ(PL ) := max{L s (Pa ), hL (J), hL (J ∪ {c})},where hL is calculated by (2.11) with hL = h and {ri , pi , qi } = {riL , piL , qiL }. Calculate λ(PR )in the same way as λ(PL ).
6. Add the new node PL to the search tree as a child node of Pa if λ(PL ) < µ and add PR ifλ(PR ) < µ.
7. Update Pa as a node with the lowest bound among nodes not yet visited.
8. Repeat Step 2 to Step 7 until there is no node that is not yet visited.
9. Output Pb as the optimal schedule for P0 .

2.8. The Well-known Benchmark Problems

27

measure the effectiveness of a certain method. Figure 2.3 shows the mt06 problem, which is theeasiest in size and structure. Indeed, it employs only 6 jobs and 6 machines, and the technologicalsequence of jobs on each machine is similar to each other.

job
1
2
3
4
5
6

Table 2.3: Muth and Thompson’s 6 × 6 problem (mt06)
machine (processing time)
3 (1)
1 (3)
2 (6)
4 (7)
6 (3)
2 (8)
3 (5)
5 (10)
6 (10)
1 (10)
3 (5)
4 (4)
6 (8)
1 (9)
2 (1)
2 (5)
1 (5)
3 (5)
4 (3)
5 (8)
3 (9)
2 (3)
5 (5)
6 (4)
1 (3)
2 (3)
4 (3)
6 (9)
1 (10)
5 (4)

5 (6)
4 (4)
5 (7)
6 (9)
4 (1)
3 (1)

Table 2.4: Muth and Thompson’s 10 × 10 problem (mt10)
job
1
2
3
4
5
6
7
8
9
10

1 (29)
1 (43)
2 (91)
2 (81)
3 (14)
3 (84)
2 (46)
3 (31)
1 (76)
2 (85)

2 (78)
3 (90)
1 (85)
3 (95)1 (6)
2 (2)
1 (37)
1 (86)
2 (69)
1 (13)

machine (processing time)
3 (9) 4 (36) 5 (49) 6 (11) 7 (62) 8 (56) 9 (44) 10 (21)
5 (75) 10 (11) 4 (69) 2 (28) 7 (46) 6 (46) 8 (72) 9 (30)
4 (39) 3 (74) 9 (90) 6 (10) 8 (12) 7 (89) 10 (45) 5 (33)
1 (71) 5 (99) 7 (9) 9 (52) 8 (85) 4 (98) 10 (22) 6 (43)
2 (22) 6 (61) 4 (26) 5 (69) 9 (21) 8 (49) 10 (72) 7 (53)
6 (52) 4 (95) 9 (48) 10 (72) 1 (47) 7 (65)5 (6)8 (25)
4 (61) 3 (13) 7 (32) 6 (21) 10 (32) 9 (89) 8 (30) 5 (55)
2 (46) 6 (74) 5 (32) 7 (88) 9 (19) 10 (48) 8 (36) 4 (79)
4 (76) 6 (51) 3 (85) 10 (11) 7 (40) 8 (89) 5 (26) 9 (74)
3 (61) 7 (7) 9 (64) 10 (76) 6 (47) 4 (52) 5 (90) 8 (45)The mt10 and mt20 problems are like brothers. They are processing the same set of operations
and technological sequences are similar, but in the mt20 problem, the number of machines available is reduced to half of that of the mt10 problem. For example, the first operation of each job
in mt10 is exactly same as the first operation of each of the first 10 jobs in mt20 and the second
operation of each job in mt10 is exactly same as the first operation of each of the second 10 jobs
in mt20.
The mt10 and mt20 problems had been a good computational challenges for a long time.
Indeed, the mt10 problem has been referred as “notorious”, because it remained unsolved forover 20 years. The mt20 problem has also been considered as quite difficult. However they are
no longer a computational challenge.
Applegate and Cook proposed a set of benchmark problems called the “ten tough problems”
as a more difficult computational challenge than the mt10 problem, by collecting difficult problems from literature, some of which still remain unsolved [1]. The ten tough problems consist of

28

Chapter 2. The Job Shop Scheduling Problem

Table 2.5: Muth and Thompson’s 20 × 5 problem
job
machine (processing time)
1 1 (29) 2 (9) 3 (49) 4 (62) 5 (44)
2 1 (43) 2 (75) 4 (69) 3 (46) 5 (72)
3 2 (91) 1 (39) 3 (90) 5 (12) 4 (45)
4 2 (81) 1 (71) 5 (9) 3 (85) 4 (22)
5 3 (14) 2 (22) 1 (26) 4 (21) 5 (72)
6 3 (84) 2 (52) 5 (48) 1 (47) 4 (6)7 2 (46) 1 (61) 3 (32) 4 (32) 5 (30)
8 3 (31) 2 (46) 1 (32) 4 (19) 5 (36)
9 1 (76) 4 (76) 3 (85) 2 (40) 5 (26)
10 2 (85) 3 (61) 1 (64) 4 (47) 5 (90)
11 2 (78) 4 (36) 1 (11) 5 (56) 3 (21)
12 3 (90) 1 (11) 2 (28) 4 (46) 5 (30)
13 1 (85) 3 (74) 2 (10) 4 (89) 5 (33)
14 3 (95) 1 (99) 2 (52) 4 (98) 5 (43)
15 1 (6) 2 (61) 5 (69) 3 (49) 4 (53)
16 2 (2) 1 (95) 4 (72) 5 (65) 3 (25)
17 1 (37) 3 (13) 2 (21) 4 (89) 5 (55)
18 1 (86) 2 (74) 5 (88) 3 (48) 4 (79)
19 2 (69) 3 (51) 1 (11) 4 (89) 5 (74)
20 1 (13) 2 (7) 3 (76) 4 (52) 5 (45)

2.8. The Well-known Benchmark Problems

29

abz7, abz8, abz9, and la21, la24, la25, la27, la29, la38, la40. The abz problems are proposed byAdams in [14]. la problems are parts of 40 problems la01-la40 originally from [32]. Table 2.6,
which is taken from [1], shows for each of the ten tough benchmark problems, the problem size,
best solution reported in [1] and best lower bound. Those gaps between the best solution andthe best lower bound suggest the difficulties of the problems and no gap means the problem issolved. The more recent status of the best solutions will be reported in the later sections.
Table 2.6: The ten tough benchmark problems (status reported by [1] in 1991)
Prob
size (n×m)
Best Solution
Best Lower Bound
abz7
20×15
668
654
abz8
20×15
687
635
abz9
20×15
707
656
la21
15×10
1053
1040
la24
15×10
935
935
la25
20×10
977
977
la27
20×10
1269
1235
la29
20×10
1195
1120
la38
15×15
1209
1184
la40
15×15
1222
1222

Taillard proposed a set of 80 JSP and 120 FSP benchmark problems. They cover various
range of sizes and difficulties. They are randomly generated by a simple algorithm. It has
become more common to use the ten tough problems and/or Taillard benchmark than to use themt10 and mt20 problems as benchmark problems.

Chapter 3
Genetic Algorithms
Genetic Algorithms (GAs) are search strategies designed after the mechanics of natural selection and natural genetics to optimize highly complex objective functions. GAs have been quite
successfully applied to optimization problems including scheduling. In this chapter, the basic
concepts of GAs are reviewed.

3.1 Basic Concepts
Genetic Algorithms use a vocabulary borrowed from natural genetics. We have a set of individuals called population. An individual has two representations called phenotype and genotype.
The phenotype represents a potential solution to the problem to be optimized in a straightforward
way used in the original formulation of the problem. The genotype, on the other hand, gives anencoded representation of a potential solution by the form of a chromosome. A chromosome ismade of genes arranged in linear succession and every gene controls the inheritance of one orseveral characters or features. For example, a chromosome consists of a sequence of 0 or 1 (i.e. a
bit string), and the value at a certain position corresponds to on (the value = 1) or off (the value =0) of a certain feature. More complicated forms such as a sequence of symbols and a permutationof alphabets are chosen for chromosomes depending of the target problem.
Each individual has its fitness, which measures how suitable is the individual for the local
environment. The Darwinian theory tells us that among individuals in a population, the one thatis the most suitable for the local environment is most likely to survive to have greater numbersof offspring. This is called a rule of “survival of the fittest.”
The objective function f of the target optimization problem plays the role of an environment,
therefore, the fitness of an individual F measures how “good” is the corresponding potential
solution in terms of the original optimization criteria. When the target optimization is the maximization of the objective function, then fitness may be identical to the objective function value:
F(x) = f (x)(3.1)

where x is an individual in the current population P. When the target is the minimization, then theobjective function must be converted so that an individual with a small objective function value30

3.2. A Simple Genetic Algorithm

31

has a high fitness. The most obvious way to deal with it is to define the fitness as the maximumof objective function over the current population minus its own objective function value:
F(x) = maxy∈P { f (y)} − f (x).

(3.2)

Another method is known as ranking. In the ranking method, each individual in the current
population P is sorted in the descending order of its objective function value so that the worst
individual is numbered as x1 and the best as xn , where n is the size of P. Then the fitness F of anindividual xi , the i-th worst individual, is defined as
F(xi ) = i.

(3.3)

3.2 A Simple Genetic Algorithm
Meanwhile, let us consider a simple case in which the genotype is a bit string of length n. A
simple genetic algorithm is composed of the following three operators:
1. Crossover
2. Mutation
3. Reproduction
Crossover and Mutation are genetic recombination operators. Each individual in a population iscoupled to pairs which is called parents at random. Each pair of individuals undergoes crossover
described as follows.
Crossover operates on genotype (i.e. chromosomes) of two individuals called parents. It
generates new (usually two) individuals called offspring whose genes are inherited from either
parents. This can be done by splitting each of the two chromosomes into fragments and recombining them again to form new chromosomes.
Now we assume that the genotype is a bit string of length n. The 1-point crossover sets one
crossover point on a string at random and takes a section before the point from one parent andtakes another section after the point from the other parent and recombines two sections to forma new bit string. For example, consider A1 and A2 being bit strings of length n = 5 as parents asfollows:
A1 = 0000 | 0
A2 = 1111 | 1.

(3.4)

The symbol | indicates a crossover point, and in this case it is set after the fourth bit. The resulting
1-point crossover yields two new individuals A01 and A02 as follows:
A01 = 0000 | 1
A02 = 1111 | 0.

(3.5)

Chapter 3. Genetic Algorithms

32

The two-point crossover sets two crossover points at random, and takes a section between the
points from one parent and other sections outside the points from the other parent and recombinesthem. In the following example, the two crossover points are set after the first and fourth bits
respectively.
A1 = 0 | 000 | 0
A2 = 1 | 111 | 1.

(3.6)

The resulting two-point crossover yields the following two individuals:
A01 = 0 | 111 | 0
A02 = 1 | 000 | 1.

(3.7)

The uniform crossover is a generalization of the two above. A random mask bit vector oflength n is given, and the positions where the mask bit is zero are taken from one parent and theother positions where the mask bit is one are taken from the other. In the following example, themask bit vector M 0 is given as M = 01010.
M = 01010
A1 = 00000
A2 = 11111.

(3.8)

The resulting uniform crossover yields the following two individuals:
A01 = 01010
A02 = 10101.

(3.9)

Mutation operates on genotype of single individual. It corresponds to an error occurred
when chromosome is copied and duplicated. When exact copies are always guaranteed, then themutation rate is zero. However in real life, copy error can happen under some circumstancessuch as the presence of noise. Mutation changes values of certain genes with small probability.
An example of a typical bit-flip mutation is shown in (3.2), where the third gene from the left in
A is selected with a small probability and its bit is flipped resulting in A0 :A = 00000
A0 = 00100.

(3.10)

Reproduction is a process in which individuals in a population are copied according to theirfitness values to form a new population. The individuals evolve through successive iterationsof reproduction, called generations. Each individual makes number of its copies proportional toits fitness, therefore, an individual with higher fitness makes more copies than that with lower
fitness.This is an artificial version of natural selection; a Darwinian survival of the fittest among
string creatures.
A simple reproduction operator is called a roulette wheel selection where each individual
in a population has a roulette wheel slot sized in proportion to its fitness. To reproduce, we3.3. The Procedure of a Simple Genetic Algorithm

33

simply spin the weighted roulette wheel and obtain a reproduction candidate with probability
proportional to its fitness. Each time we require another offspring, a simple spin of the weighted
roulette wheel yields the reproduction candidate. An example of the roulette wheel selection isshown in Figure 3.1. In the upper left picture in the figure, there are total seven individuals, ID 1
to ID 7, with fitness value assigned. A roulette wheel is created and shown in the right picture.
The number in each wheel slot corresponds to the individual ID. The first individual with ID 1,
for example, has fitness value 9, which is the highest and therefore the largest slot is assigned.
We spin the roulette wheel seven times to select seven individuals. Individual ID 1 and ID 2 arelikely to be selected more than once but individual ID 4 and ID 6 are not likely to be selected
resulting in the seven individuals with some duplicates shown in the lower picture.

ID
1
2
3
4
5
6
7

Individuals
fitness
10101101001110
901100010100011
810011001100101
510001001100011
111011100010010
400110011101111
110111011010111
5
1
2
3
1
5
2
7

10101101001110
01100010100011
10011001100101
10101101001110
11011100010010
01100010100011
10111001010111

roulette wheel

7
6

1

5
4
3

2

Figure 3.1: An example of the roulette wheel selection, where the roulette wheel
is created according to the fitness value of each individual shown in the upper left
picture3.3 The Procedure of a Simple Genetic Algorithm
The general procedure of a simple GA can be summarized as in Algorithm 3.3.1. In the algorithm, we start from a random initial population P(0). P(t) is a population at generation t with N
individuals. Rc × N members are randomly selected from P(t) and crossover is applied to generate new Rc × N individuals that join into a new population P0 (t) in Step 2, where Rc < 1 is calledChapter 3. Genetic Algorithms

34
Algorithm 3.3.1(Simple Genetic Algorithm)

1. Initialize P(t) as a random population P(t = 0)
2. Recombine P(t) to yield P0 (t) by crossover and mutation
3. Evaluate P0 (t)4. Reproduce P(t + 1) from P0 (t) by selection
5. Set t ← t + 1
6. repeat from 2 to 5 until some termination condition is met.
7. Output the best individual in P(t).

a crossover ratio. The rest of P(t) is just copied to P0 (t). Rm × N members are then randomly
selected from P0 (t) and mutation is applied to generate new individuals that replace the original,
where Rm < 1 is call a mutation ratio. When the best individual in P(t) is preserved and copied
to P0 (t) without modification, it is called elitist strategy.
P0 (t) is evaluated in Step 3 and the new population P(t + 1) is obtained after the reproduction
using, for example, the roulette wheel selection in step4. The termination condition is usually
given as: when t is sufficiently large, when the best or average fitness in P(t) exceeds certain
value, or when the variation of the fitness in P(t) is small.
While the process described above is repeated for a sufficient number of generations, therecombination operators keep producing possibly new individuals with new fitness where someof them are possibly better than those of ever existing ones. The reproduction phase focuseson such good individuals and replicate them as occurred in the natural evolution. Eventually an
individual with a high fitness value is expected to emerge in a population. The natural evolution
process requires enormous amount of time. However, this simulated evolution process on acomputer runs much faster.

Chapter 4
A Simple Genetic Algorithm for the
Jobshop Scheduling Problem
In this chapter, the simple GA described in the previous chapter is applied to the jobshop scheduling problem. The approach described in this chapter was proposed by Nakano and Yamada [24].
An advantage of this approach is that conventional genetic operators, such as one-point, twopoint and uniform crossovers can be applied without any modification. However, one drawback
is that a new individual generated by crossover may not represent a feasible schedule. In otherwords, such genotype is called fatal or illegal. There are two approaches to solve this situation:
one is to repair a fatal genotype to a normal one, and the other is to impose a penalty for thefatality and to lower the fitness. One example of the former approach will be elaborated in thischapter.

4.1 Genetic Encoding of a Solution Schedule
We have a n×m jobshop scheduling problem (JSP). As described in Chapter 2, a solution of a
JSP can be represented as a directed graph. Therefore, by labeling each directed arc as 0 or 1
according to its direction, it can be represented as a bit string of length mn(n − 1)/2. For example,
consider a 3 × 3 problem given in Table 2.1 and a solution given in Figure 2.2 in Chapter 2.
According to Figure 2.8, each arc of the graph has a labelling 0 or 1. The only thing we need to
do is to specify the order of arcs. Note that each arc represents the precedence relation between
two jobs Ji and J j on the same machine Mr ; hence an arc is specified by a triplet (i, j, r). An
intuitive ordering between two arcs (i, j, r) and (i0 , j0 , r0 ) is a machine-based ordering defined as:(i, j, r) < (i0 , j0 , r0 ) ⇐⇒ (r < r0 or (r = r0 and (i, j) < (i0 , j0 )))(4.1)(i, j) < (i0 , j0 ) ⇐⇒ (i < i0 or (i = i0 and j < j0 ))(4.2)

where,The solution schedule given in Figure 2.2 and Figure 2.8 is encoded as follows:
111 | 100 | 011
35

(4.3)

36

Chapter 4. A Simple Genetic Algorithm for the Jobshop Scheduling Problem

The encoding from a schedule to a bit string based on the machine-based ordering of arcs in the
disjunctive graph corresponding to the schedule is called machine-based encoding and decoding.
Another ordering is called a job-based ordering defined as follows:(i, j, r) < (i0 , j0 , r0 )⇐⇒ ((i, j) < (i0 , j0 ) or ((i, j) = (i0 , j0 ) and o j (r) < o j (r0 )))(4.4)

where o j (r) < o j (r0 ) indicates that J j is processed on Mr prior to Mr0 . In other words, two arcs
corresponding to the same job pair (Ji , J j ) are ordered according to the processing order of thefirst job Ji . If we use the job-based ordering, then solution schedule given in Figure 2.2 andFigure 2.8 is encoded as follows:
100 | 101 | 101(4.5)

where the vertical bars are inserted as job-pair delimiters; the partitions correspond to job pair
(1, 2), (1, 3) and (2, 3) from left to right respectively. The encoding from a schedule to a bit string
based on the job-based ordering of arcs in the disjunctive graph corresponding to the schedule iscalled job-based encoding and decoding.
As another example, Figure 4.1 shows a simplified Gantt chart representation of an optimalschedule for the mt06 problem defined in Table 2.3. In the figure, the number indicates job
number and consecutive sequence of the same number represents an operation for the job. The
repetitions of the same number represents the processing time. For example, the left most sequence 111 on machine M1 represents an operation for job J1 on machine M1 with processing of
3 time units and it starts at time unit 7. Figure 4.2 shows a binary representation of the optimal
solution given in Figure 4.1 using the job-based ordering. For the ease of understanding, one
long bit string is partitioned and divided per each job pair. For example, the first bit substring
represents that J1 is processed prior to J2 on M3 and also on M1 but J2 is prior to J1 on M2 . Notethat in the optimal schedule, like the one in this example, there is a tendency that the same bit
continues in each substring of the same job pair. This confirms a heuristic that the processing
priority for each job pair tends to be unchanged. This is especially true for easy problems inwhich each technological sequence of jobs on each machine is similar to each other.
As described in the previous chapter, the simple one-point or two-point crossover exchange
chunks of bit sequences between parents. By using the job-based ordering, the consecutive same
bits are likely to be exchanged together and thus this tendency, once acquired, is not destroyed
easily.

4.2 Local harmonization
In the previous section, we have seen a couple of encoding methods to convert a solution schedule
into a bit string. In those methods, different solution schedules are mapped into different bit
strings. However, an arbitrary bit string generated by hand or crossover or mutation may notnecessarily mapped back into a feasible solution schedule. In fact, the directed graph obtained
from any bit string by selecting each arc’s direction according to zero or one of corresponding bit

4.2. Local harmonization

37

M1 :
111
44444333333333 66666666662222222222555
M2 : 2222222244444666111111555 3
M3 : 333331 2222255555555544444
6
M4 :
3333
666
4441111111
22225
M5 :
2222222222 555553333333444444446666111111
M6 :33333333 66666666622222222225555111444444444
Figure 4.1: An optimal schedule for the mt06 (6 × 6) problem (makespan = 55)(J1 , J2 ) : 110100(J1 , J3 ) : 011000(J1 , J4 ) : 110010(J1 , J5 ) : 111100(J1 , J6 ) : 110000(J2 , J3 ) : 101000(J2 , J4 ) : 111100(J2 , J5 ) : 111111(J2 , J6 ) : 111000(J3 , J4 ) : 111001(J3 , J5 ) : 111100(J3 , J6 ) : 111101(J4 , J5 ) : 110100(J4 , J6 ) : 111010(J5 , J6 ) : 101000
Figure 4.2: A binary representation of a solution schedule using the job-based ordering corresponding to the solution given in Figure 4.1. The first line corresponds to the precedence relation
between J1 and J2 . The first three digits of the bit-string on the first line are 110. This corresponds to the fact that J1 is processed prior to J2 on J1 ’s first and second machines M3 and M1 ,but is not prior to J2 on J1 ’s third machine M2 and so on.

38

Chapter 4. A Simple Genetic Algorithm for the Jobshop Scheduling Problem

value may contain cycles. In such cases, the bit string, i.e. genotype, is called fatal or illegal, anda bit string is called feasible when it corresponds to an executable schedule with corresponding
directed graph being acyclic.
A repairing procedure that generates a feasible bit string, as similar to an illegal one as possible, is called the harmonization algorithm [24]. The Hamming distance is used to assess thesimilarity between two bit strings. The harmonization algorithm goes through two phases: local harmonization and global harmonization. The former removes the ordering inconsistencies
within each machine, while the latter removes the ordering inconsistencies between machines.
This section explains the former and the next section will explain the latter harmonization.
The local harmonization works separately for each machine and resolves the cycle withineach machine by changing directions of arcs. Assume that we are resolving cycles on machine
Mr . A set of nodes Gr of machine Mr is initialized as Gr = {O1r , O2r , . . . , Onr } (operations andnodes are identified here). First the least priority node Olr ∈ Gr is identified as a node that has
the highest number of incoming arcs from Gr \ Olr (break ties arbitrarily). If this node has any
outgoing arc to any node in Gr \ Olr , then the direction of the arc is reversed so that Olr has only
incoming nodes from Gr \ Olr , and Gr is updated as Gr := Gr \ Olr . This process is repeated until
Gr becomes empty, and as a result, the local inconsistensy is completely removed. Figure 4.3
shows an example of the local harmonization for machine M1 . In the figure, O31 is first identified
as the least priority node, and the arc O31 → O61 is reversed such that O31 becomes the lastoperation on machine M1 . O21 is then identified as the second least priority node, and the arc O21→ O41 is reversed such that O21 becomes the second last operation on machine M1 and so on.
The obtained consistent ordering is O41 → O61 → O51 → O11 → O21 → O31 .

4.3 Global harmonization
The global harmonization removes ordering inconsistencies between machines. Even after the
local harmonization, there may exist cycles in the graph. In Figure 4.4, for example, there is acycle connecting O23 → O22 → O32 → O31 → O33 , and again, → O23 . The global harmonization
changes the directions of minimum number of arcs so that there exists no cycles. It is not always
guaranteed that the above harmonization will generate a feasible bit string closest to the original
illegal one, but the resulting one will be reasonably close and the harmonization algorithms arequite efficient.

4.4 Forcing
An illegal bit string produced by genetic operations can be considered as a genotype, and afeasible schedule generated by any repairing method can be regarded as a phenotype. Then theformer is an inherited character and the latter is an acquired one. Note that the repairing statedabove is only used for the fitness evaluation of the original bit string; that is, the repairing does
not mean the replacement of bit strings.
Forcing means the replacement of the original string with a feasible one. Hence forcing can4.4. Forcing

39

O11

O11

O61

O21

O61

O21

O51

O31

O51

O31

O41

O41

O11

O11

O61

O21

O61

O21

O51

O31

O51

O31

O41

O41

Figure 4.3: An example of the local harmonization resolving cycles within six operations O11 , O21 , . . . , O61 on the same machine M1 where the arcs O31 → O61 , O21 →O41 and O11 → O61 are reversed in this order and a consistent ordering O41 → O61→ O51 → O11 → O21 → O31 is eventually obtained.

O11

0

O12

O13

O22

O23

O21

*

10
O32

O31

O33

Figure 4.4: An example of global harmonization where a cycle O23 → O22 → O32 →O31 → O33 → O23 is resolved by reversing an arc O22 →O32 .

40

Chapter 4. A Simple Genetic Algorithm for the Jobshop Scheduling Problem

be considered as the inheritance of an acquired character, although it is not widely believed thatsuch inheritance occurs in nature. Since frequent forcing may destroy whatever potential anddiversity of the population, it is limited to a small number of elites, usually the best 5% in the
population. Such limited forcing brings about at least two merits: a significant improvement in
the convergence speed and the solution quality.

4.5 Simple GA for the JSP
Using the job-based encoding, standard crossover/mutation, global/local harmonizations andforcing, a simple GA for the JSP can be constructed. Because the JSP is a minimization problem, the fitness is defined by the ranking method (3.1) and standard roulette wheel selection
is utilized. The outline of the simple GA for the jobshop scheduling problem is described in
Algorithm 4.5.1.

4.6 The Limitation of the Simple Approach
The simple GA approach described in this chapter can be applied to small problems such as 6 × 6
problem given in Table 2.3. Table 4.1 summarizes the experimental results for the mt benchmark
problems. The column labeled SGA shows the best makespans obtained by the SGA and thecolumn labeled Optimal shows the known optimal makespans. In fact, the optimal scheduleshown in Figure 4.1 is obtained by the GA. However, larger problems such as 10 × 10 and 20 × 5
are not tractable by this simple approach.
Table 4.1: Experimental results of the simple GA for mt benchmark problems
Prob.
mt06
mt10
(size)(6×6) (10×10)
SGA
55
965
Optimal
55
930

mt20
(20×5)
1215
1165

4.6. The Limitation of the Simple Approach

41

Algorithm 4.5.1 (Simple GA for the JSP)
A n×m scheduling problem is given as an input. The GA parameters: N, the population size, Rc ,crossover ratio, Rm , mutation ratio are also given.
1. Initialize P(t) as a random population P(t = 0) of size N, where each random individual isa bit string of length m × n × (n − 1)/2.
2. Modify P(t) by applying one-point (3.2), two-point (3.2) or uniform (3.2) crossover to therandomly selected Rc × N members of P(t) and obtain P0 (t).
3. Modify P0 (t) by apply bit-flip mutation (3.2) to the randomly selected Rm × N members of
P0 (t) and obtain P00 (t).
4. Evaluate P00 (t) by the following steps;(a) Decode each individual p in P00 (t) by using the job-based decoding based on (4.1)
into S , with the local and global harmonization methods to repair illegal bit strings.
(b) Calculate the objective function f of p as f (p) = Cmax (S )(c) Calculate fitness F of p by using the ranking method shown in (3.1).
(d) Apply forcing to retain the phenotype of small number of elitest individuals to thenext generation.
5. Reproduce P(t + 1) from P00 (t) by the roulette wheel selection
6. Set t ← t + 1
7. repeat from 2 to 5 until some termination condition is met.
8. Output the best individual in P(t).

Chapter 5
GT-GA: A Genetic Algorithm based on the
GT Algorithm
As seen in the previous section, conventional GAs can be applied to the jobshop schedulingproblem in a rather straightforward way without major difficulties; a solution is represented as a
bit string and conventional genetic operators such as 1-point, 2-point and uniform crossover andbit-flip mutation are applied. Because of the complicated constraints of the problem, however,an individual generated by such genetic operators is often infeasible; its phenotype does not
represent an executable solution, and requires several steps of repairing process such as local andglobal harmonizations.
Obviously one of the advantages of the GA is its robustness over a wide range of problems
with no requirement of domain specific adaptations. Hence genetic operators deal with genotype, which is domain independent, and are separated from domain specific decoding processfrom genotype to phenotype. However from the performance viewpoint, it is often more efficient to directly incorporate domain specific features into the genetic operators and skip wastefulintermediate decoding steps. Thus the GT crossover and the genetic algorithm based on GT
crossover, denoted as GT-GA, has been proposed by Yamada and Nakano [25] and has the following properties.
• The GT crossover is a problem dependent crossover operator that utilizes the GT algorithm.
• The crossover operates directly on phenotype.
• In the crossover, parents cooperatively give a series of decisions; as which operation should
be processed next, to build a new schedule. These decisions are made based on their own
scheduling orders.
• The offspring represent active schedules, so there is no repairing process required.
Before describing the GT crossover, let us review some other crossover operators based on
non-binary encodings for comparisons and discuss their advantages and disadvantages in thefollowing sections.
42

5.1. Subsequence Exchange Crossover

43

5.1 Subsequence Exchange Crossover
As shown in Figure 2.3 in Section 2.1, a schedule of the JSP can be represented by a solution
matrix, in other words, the set of permutations of jobs on each machine. When the matrix is
expanded in one dimensional array as shown in Figure 5.1, it is called an m-partitioned permutation, where the permutation in the k-th partition (from the left) corresponds to the processingorder of jobs on machine Mk . A solution represented by a m-partitioned permutation is regarded
as genotype to which a crossover operator is applied.
M1
1 2 3

M2
3 1 2

M3
2 1 3

Figure 5.1: The solution given in Figure 2.3 is converted to an m-partitioned permutation for m = 3, where the permutation in the k-th partition corresponds to theprocessing order of jobs on machine Mk

The Subsequence Exchange Crossover (SXX) was proposed by Kobayashi, Ono and Yamamura [33]. The SXX is a natural extension of the subtour exchange crossover for TSPs presented
by the same authors [34]. Let two m-partitioned permutations be p0 and p1 , which correspond to
two feasible solution schedules. A pair of subsequences, one from p0 and the other from p1 on
the same machine, is called exchangeable if and only if they consist of the same set of job numbers. The SXX first identifies exchangeable subsequence pairs in p0 and p1 on each machine andinterchanges each pair to produce new m-partitioned permutations k0 and k1 . Figure 5.2 shows
an example of the SXX for a 6 × 3 problem. In the figure, each underlined subsequence pair is
identified as exchangeable and interchanged.
The SXX ensures that k0 and k1 are always valid m-partitioned permutations and therefore,
there are no inconsistencies within each machine to be resolved by the local harmonization described in the previous section. However, there may exist inconsistencies between machines thatmust be resolved by the global harmonization.

5.2 Precedence Preservative Crossover
Another representation that uses an unpartitioned permutation of n job numbers with m-repetitions
has been proposed by Bierwirth [35]. In this representation, we consider a permutation of n job
numbers but each identical job number occurs m times. When such a permutation with repetitions is given, it is decoded into a feasible schedule by scanning the permutation from left toright and referring the k-th occurrence of a job number to the k-th operation in the technologicalsequence of this job. Figure 5.3 shows an example of this decoding process. In the figure, a permutation of three job numbers with three repetitions is given in the top and there are three rows
labeled M1 , M2 and M3 in the bottom. We consider decoding this permutation into a solution
schedule of 3 × 3 problem given in Table 2.1 and Figure 2.1. From the job sequence matrix {T jk }given in Figure 2.1, we see that job J1 , for example, is first processed on M1 , therefore, the first44

Chapter 5. GT-GA: A Genetic Algorithm based on the GT Algorithm
M1

M2

M3

p0 123456 321564 235614
p1 621345 326451 635421
k0 213456 325164 263514
k1 612345 326415 356421
Figure 5.2: An example of subsequence exchange crossover (SXX), where eachunderlined subsequence pair one from p0 and the other from p1 on each machine is
identified as exchangeable and interchanged to generate k0 and k1

occurrence of 1 in the permutation should be moved straight down to the row of M1 . Likewise,
J1 is then processed on M2 and M3 in this order, so the second and third occurrences of 1 should
be moved down to the rows of M2 and M3 respectively. By moving all the job numbers down to
one of M1 , M2 , or M3 rows, a solution schedule is obtained. In this case the schedule obtained is
identical to the one given in Figure 2.3.
The advantage of this representation is that an arbitrary permutation with repetitions can be
decoded into a feasible schedule. Therefore, no repairing processes such as local and global
harmonizations are required.
A job permutation
is decodedtoa schedule

1 3 2 1 3 2 2 1 3

M1 1
M2
M3

2
3

3
1

2
2

1 3

Figure 5.3: A job sequence (permutation with repetition) for a 3 × 3 problem definedin Figure 2.1 is decoded to a schedule, which is equivalent to the one in Figure 2.3.

A crossover operator called Precedence Preservative Crossover (PPX) is proposed for thisrepresentation in [36]. The PPX perfectly respects the absolute order of genes in parental chromosomes as follows. Assume we have two parents p0 and p1 encoded in permutation withrepetitions representation and consider generating a new individual k also represented in a permutation with repetitions. First, a template bit string h of length mn is given that determines fromwhich parent, p0 or p1 , should genes to be drawn to generate a new individual. The bit value iszero means that the corresponding gene should be copied from p0 and one from p1 . When agene is drawn from one parent and then appended to the offspring chromosome, it is deleted
from the parent and the corresponding gene is also deleted from the other parent. This step isrepeated until both parent chromosomes are empty and the offspring contains all genes involved.

5.3. GT Crossover

45

Figure 5.4 shows an example of this crossover. Starting from the top left picture, the first two bitsof h are both zero, so the first and second job numbers 3 and 2 in k are copied from p0 as shownin the round boxes. The leftmost occurrences of job numbers 3 and 2 are deleted from both p0
and p1 in the top right picture. The first two bits in h are deleted as well. Then the leftmost
non-deleted bits in h become four ones shown in the square box, which means that the next fourjob numbers should be copied from p1 . The leftmost non-deleted four job numbers in p1 are
1121 which are copied to k as shown in the square boxes. In the bottom picture, the leftmost
non-deleted occurrences of job numbers 1121 are then deleted from p0 as well as from p1 . The
four ones in the square box are also deleted from h. The remaining non-deleted bits in h are threezeros, which indicates that the remaining job numbers should be copied from p0 (however theremaining non-deleted permutation 233 is identical both in p0 and p1 in this example).
p0 3 2 2 2 3 1 1 1 3

p0 3 2 2 2 3 1 1 1 3hh

0 0 1 1 1 1 0 0 0

0 0 1 1 1 1 0 0 0

p1 1 1 3 2 2 1 2 3 3

p1 1 1 3 2 2 1 2 3 3kk

3 2 1 1 2 1 2 3 3

3 2 1 1 2 1 2 3 3

p0 3 2 2 2 3 1 1 1 3h

0 0 1 1 1 1 0 0 0

p1 1 1 3 2 2 1 2 3 3k

3 2 1 1 2 1 2 3 3Figure 5.4: An example of the precedence preservative crossover (PPX), where k is
generated from p0 and p1 using h

5.3 GT Crossover
Unlike other crossover operators described in the previous two sections, the GT crossover, GTX
in short, directly operates on the solution matrix representation of a schedule given in Figure 2.3.
In this sense, we have no distinction between genotype and phenotype here. Assume we have twoparents p0 and p1 both represented by solution matrices and consider generating a new individualk also represented in a solution matrix. Let H be a binary matrix of size m × n,where Hri = 0
means that the i-th operation on machine r should be determined by the first parent p0 and Hri = 1
by the second parent p1 [25, 37]. H is called a inheritance matrix. The role of Hri is similar tothat of h described in Section 5.2. In fact, the idea of the GT crossover including the use of Hriis first proposed and later adopted to the precedence preservative crossover.

46

Chapter 5. GT-GA: A Genetic Algorithm based on the GT Algorithm

Algorithm 5.3.1 (GT crossover)A scheduling problem represented by {T jk }, the technological sequence matrix, and {p jk }, theprocessing time matrix as well as two solution schedules p0 and p1 represented by solution
01
matrices S 0 = {S rk} and S 1 = {S rk} respectively, are given as inputs.
1. Initialize G as a set of operations that are first in the technological sequence; i.e., G ={O1T11 , O2T21 , . . . , O2Tn1 }. For each operation O ∈ G, set ES (O) := 0 and EC(O) := p(O).
2. Find the earliest completable operation O∗r ∈ G by (2.1) with machine Mr . A subset of G
that consists of operations processed on machine Mr is denoted as Gr .
3. Calculate the conflict set C[Mr , i] ⊂ Gr by (2.2), where i−1 is the number of operations
already scheduled on Mr .
4. Select one of the parents {p0 , p1 } as p according to the value of Hri , that is, p := pHriand S p := S Hri . For each O jr ∈ C[Mr , i] with job number j, there exists an index l suchthat S rl = j. Let lm be the smallest index number among them; i.e., lm := min{l | S rl =j and O jr ∈ C[Mr , i]} and let k := S rlm . This results in selecting an operation Okr ∈C[Mr , i] that has been scheduled in p earliest among the members of C[Mr , i].
5. Schedule Okr as the i-th operation on Mr ; i.e. S ri := k, with its starting and completion
times equal to ES (Okr ) and EC(Okr ) respectively: s(Okr ) = ES (Okr ), c(Okr ) = E(COkr ).
6. For all O jr ∈ Gr \ {Okr }, Update ES (O jr ) as
ES (O jr ) := max{ES (O jr ), EC(Okr )} and EC(O jr ) as EC(Okr ) := ES (Okr ) + p(Okr ).
7. Remove Okr from G (and therefore from Gr ), and add operation Oks that is the next to
Okr in the technological sequence to G if such Oks exits; i.e., if r = T ki and i < m, thens := T ki+1 and G := (G \ {Okr }) ∪ {Oks }.
Calculate ES (Oks ) and EC(Oks ) as:
ES (Oks ) := max{EC(Okr ), EC(PM(Oks ))} and EC(Oks ) := ES (Oks ) + p(Oks ) respectively.
8. Repeat from Step 1 to Step 7 until all operations are scheduled.
9. Output the solution matrix {S rk } as the active schedule obtained with the set of starting andcompletion times {s(O jr )} and {c(O jr )} respectively where j = S rk .

5.3. GT Crossover

47

The GT crossover can be defined by modifying Step 4 of Algorithm 2.2.1, where the choiceof the next operation from the conflict set C[Mr , i] was at random. In the GT crossover, thechoice is made by looking at the processing order in one of the parents specified by H and an
operation that has been scheduled in the parent earliest among the members of the conflict set
is selected. By doing so, it tries to reflect the processing order of the parent schedules to theiroffspring. Note that if the parents are identical to each other, the resulting new schedule is also
identical to those of the parents. In general the new schedule inherits partial job sequences ofboth parents in different ratios depending on the number of zeros and ones contained in H.
Algorithm 5.3.1 describes the GT crossover. For the purpose of self-completeness, it is presented as a complete form, but the differences between the GT algorithm and the GT crossover
are only the inputs and Step 4. The other steps are just the exact copies of the GT algorithm.
The GT crossover generates only one schedule at once. Another schedule is generated by using
the same H but changing the roles of p0 and p1 . Thus two new schedules are generated thatcomplement each other.
Figure 5.5 shows an example of the GT crossover applied to the two parents p0 and p1 represented by solution matrices with an inheritance matrix H and generating k as an offspring when aproblem is given as shown in Figure 2.3. For better understanding, correspondig solutions represented by the simplified gantt chart introduced in Figure 4.1 are also shown in the square boxes.
Each number with an arrow pointing to an operation indicates the order of the corresponding
operation selected in Step 4 of Algorithm 5.3.1. This order is dynamically assigned in the algorithm. For example, the first operation on machine M1 should be first determined in this case.
The corresponding bit in H is consulted and here, we see H11 = 1 which means the operation
should be determined from p1 . Because the first operation on machine M1 in p1 is O11 , O11 isalso scheduled as the first operation on machine M1 in k. The next operation to be determined isthe first operation on machine M2 and H21 = 1 indicates that this operation should be determinedagain from p1 , thus O32 is selected in k. In the similar way, the next operation to be determined isthe second operation on machine M1 and H12 indicates that this operation should be determined
this time from p0 , thus O21 is selected in k. Repeating this process, the complete schedule shownin the right is finally obtained as k. One can see, for example, that the whole job sequence on M3
is consequently copied from p0 because corresponding bits in H are all zeros.
2 3 1

33111
p0= 3 2 1 ~~ 22
333 2222111
2 3 1

2223

1

111

3 5
9

1112233
333111 2222
22231111 0 0

H= 1 1 1
0 0 0
2
1 3 2

6

1113322

p1= 3 1 2 ~~ 333111
3 1 2

4
7 8

2222

3111222

Figure 5.5: GT crossover

1 2 3

k= 3 1 2
2 3 1

Chapter 5. GT-GA: A Genetic Algorithm based on the GT Algorithm

48

While applying the GT crossover, simulated random copy error is incorporated as mutation
built into the GT crossover. More precisely, in Step 4 of Algorithm 5.3.1, the n-th (n > 1) smallest
index number ln m is selected instead of lm = l1 m and the corresponding operation in C[Mr∗ , i] thatis the n-th earliest scheduled operation in p among the members of C[Mr∗ , i] is selected with asmall probability Rm When the two parents p0 and p1 are identical, then the offspring k generated
by the GT crossover is also identical to p0 and p1 , however, with the mutation incorporated, kcan be slightly different from the parents.

5.4 GT-GA
A Genetic Algorithm based on the GT crossover is straightforward. The general procedure described in Section 3.3 is used without major modifications. The following points should be
mentioned.
1. In the GT-GA, each individual is always a feasible schedule represented by a solution matrix. In fact, each individual is not only feasible but also an active schedule. As describedearlier, we have no distinction between genotype and phenotype here.
2. Because the problem is a minimization problem, the rank-based roulette wheel selection
method is used.
3. An elitest strategy to preserve the best individual in the current population to the next
generation is used.

5.5 Computational Experiments
GT-GA is applied to the mt benchmark problems to explore its efficiencies and limitations. Table 5.1 shows the best solutions obtained by the GT-GA for each problem. For the mt06 problem
the optimal schedule with makespan 55 is immediately obtained even with small population. For
the mt10 and mt20 problems, 600 trials are performed with different random number seeds ineach trial.
For the GA parameters, the population size N < 100 is used for the mt06 problem, N = 1000
and N = 2000 are used for the mt10 and mt20 problems respectively with high crossover rate
Rc > 0.9 and low mutation rate Rm < 0.01. Each GT-GA run is terminated after 200 generations.
Figure 5.6 shows the histgram of the obtained best solutions of the mt10 problem for 600 trials.
For example, the optimal schedule of the mt10 problem with makespan 930 was obtained fourtimes among 600 trials.
From the experimental results, we can observe that notorious mt10 problem can be solved
to optimally even with this simple algorithm. The success is limited in a sense that the optimal
makespan for the mt10 problem is obtained only occasionally among many trials and for the
mt20 problem, the optimal makespan cannot be obtained. However, considering the simplicityof the algorithm, the results are still interesting.

5.5. Computational Experiments

49

Algorithm 5.4.1 (A Genetic Algorithm using the GT crossover)As always, we are given a jobshop scheduling problem represented by {T jk }, the technologicalsequence matrix, and {p jk }, the processing time matrix. Besides, the following GA parameters
are given: population size N, crossover rate Rc and mutation rate Rm .
1. A random initial population P(t = 0) of size N is constructed in which each individual isgenerated using the GT algorithm with randomly selecting operations in Step 4 of Algorithm 2.2.1. The makespan of each individual is automatically calculated as an output of
the GT algorithm.
2. Select randomly N × Rc individuals from P(t) and pair them randomly. Apply the GT
crossover (with built-in mutation of probability Rm ) to each pair and generate new N × Rc
individuals that are inserted into P0 (t). The rest of P(t) members are just copied to P0 (t). As
a result of the GT crossover, the makespan of each individual is automatically calculated.
3. If the best makespan in P0 (t) is not as good as that in P(t), then the worst individual in P0 (t)is replaced by the best individual in P(t) (elitest strategy).
4. Reproduce P(t + 1) from P0 (t) by using the rank-based roulette wheel selection, in which
each individual in P0 (t) is sorted in the descending order of its makespan so that the worst
individual is numbered as x1 and the best as xN . Then the roulette wheel selection is appliedwith the fitness f of an individul xi defined as f (xi ) = i to obtain P(t + 1).
5. Set t ← t + 1
6. Repeat from Step 2 to 5 until some termination condition is met.
7. Output the best individual in P(t) as the obtained best solution.

Table 5.1: Experimental results of the GT-GA for mt benchmark problems
Prob.
mt06
mt10
(size)(6×6) (10×10)
SGA
55
965
GTGA
55
930
Optimal
55
930

mt20
(20×5)
1215
1184
1165

Chapter 5. GT-GA: A Genetic Algorithm based on the GT Algorithm

50

Schedules
90

80

70

60

50

40

30

20

10

0
940

960

980

1000

1020
MakespanFigure 5.6: The histgram of the best makespans obtained by the GT-GA after 200
generations among 600 trials for the mt10 problem

5.6. Concluding Remarks

51

5.6 Concluding Remarks
Figure 5.7 summarizes the relationship between various crossover operators discussed in thisand previous chapter. The ordinate indicates the level of problem independence of each solution
representation. The SGA described in the previous chapter uses the binary codings and does notuse any domain knowledge of the scheduling problem, therefore, gives the most general representation. The SXX uses the property that the scheduling problem can be represented by the
permutation problem but does not incorporate the fact that permutations on each machine are notmutually independent. The PPX utilizes this fact but does not directly use the fact that in the
makespan minimizing scheduling problem, the optimal schedules are always active schedules,
therefore, we can restrict the search for the optimal schedule within the set of all active schedules. The GT crossover denoted as GTX in the figure incorporates all these domain knowledges
therefore the most efficient but still simple enough.

general(any problem)

SGA with one-point crossoverdomain knowledges(genotype is represented by a bit string)

permutation problem

SXX(genotype is represented by a mpartitioned permutaion)jobshop scheduling problem

PPX(genotype is represented by a
permutation with m repetitions)

makespan-minimizing jobshop
scheduling problem in which
active schedules make sense

GTX(genotype and phenotype are identicaland represented by a solution matrix)

problem specificFigure 5.7: relationship between various crossover operators

Chapter 6
Neighborhood Search
As is now universally appreciated, it is not really likely that optimal solutions to large combinatorial problems will be found reliably by any exact method, although it is possible to find classes of
instances where problem-specific methods can achieve good results. However, for problems thatare NP-hard [38], it is now customary to rely on the application of heuristic techniques [39, 40].
These techniques include what some call the ‘metaheuristics’—simulated annealing (SA) andtabu search (TS)—as well as genetic algorithms (GAs) which are already discussed in the earlier chapters. Central to most heuristic search techniques is the concept of neighborhood search(NS). In this chapter, the general concept of the neighborhood search is first reviewed and thewell-known instances of metaheuristics, SA, TS, and GAs are formulated in this context so thatthe differences and characteristics of those methods become clear.

6.1 The Concept of the Neighborhood Search
If we assume that a solution is specified by a vector x, that the set of all (feasible) solutions isdenoted by X (which we shall also call the search space), and the cost of solution x is denotedby f (x), then every solution x ∈ X has an associated set of neighbors, N(x) ⊂ X, called theneighborhood of x. Each solution x0 ∈ N(x) can be reached directly from x by an operation
called a move, a single perturbation of x. Many different types of move are possible in any
particular case, and we can view a move as being generated by the application of a transition
operator ω. For example, if X is the n-dimensional binary hypercube ZZn2 , a simple transition
operator is the bit flip φ(k)(φ(k) : ZZn2 → ZZn2

zk 7→ 1 − zk
zi 7→ zi
if k , i(6.1)

where z is a binary vector of length l.
As another example, we can take the forward shift operator for the case where X is Πn —the52

6.1. The Concept of the Neighborhood Search

53

space of permutations π of length n. The operator F SH(i, j) (where we assume i < j) is

F SH(i, j) : Πn → Πn
π 7→ πk−1 if i < k ≤ j k
πi 7→ π j
 π 7→ πotherwise
k
k

(6.2)

The permutation flowshop scheduling problem with n jobs and m machines and with any objective function, such as n/m/P/Cmax or n/m/P/C sum introduced in Chapter 2 can be considered asa typical example of this permutation space.
An analogous backward shift operator BSH(i, j) can similarly be described; the composite
of BSH and F SH is denoted by SH. Another alternative for such problems is an exchange
operator EX(i, j) which simply exchanges the elements in the ith and jth positions.

Algorithm 6.1.1 A general structure of Neighborhood Search
A cost function of x ∈ X is given as f (x) and neighborhood of x as N(x). Certain criteria are
given to select y ∈ N(x) based on the value f (y).
1. Select a starting point x0 ∈ X at random and set x = xbest = x0 .
2. do(a) a candidate y is chosen from N(x) and is accepted or rejected according to the given
criteria based on the value f (y). Set x = y if y is accepted, otherwise repeat this step
until some y is accepted .
(b) If f (x) < f (xbest ) then set xbest = x.
until termination conditions are satisfied.
3. Output xbest as the best solution obtained.

A typical neighborhood search (NS) heuristic procedure is shown in Figure 6.1.1. As shownin the figure, NS operates by generating neighbors in an iterative process where a move to a new
solution is made whenever certain criteria are fulfilled in Step 2a. There is a great variety of ways
in which candidate moves can be chosen for consideration, and in defining criteria for accepting
candidate moves. Perhaps the most common case is that of ascent, in which the only moves accepted are to neighbors that improve the current solution. Steepest ascent corresponds to the casewhere all neighbors are evaluated before a move is made—that move being the best available.
Next ascent is similar, but the next candidate (in some pre-defined sequence) that improves thecurrent solution is accepted, without necessarily examining the complete neighborhood. Normally, the search terminates when no moves can be accepted.

54

Chapter 6. Neighborhood Search

6.2 Avoiding Local Optima
The trouble with NS is that the solution it generates is usually only a local optimum—a point in
the search space none of whose neighbors offer an improved solution and NS does not guarantee
to find the global optimum, the very best solution in the entire search space. In recent years manytechniques have been suggested for the avoidance of local optima. At the most basic level, we
could use iterative restarts of NS from many different initial points, thus generating a collectionof local optima from which the best can be selected. There are more popular and intelligent
principles. For completeness, We refer here briefly to some of the most popular ones. Simulated
annealing uses a controlled randomization strategy—inferior moves are accepted probabilistically, the chance of such acceptance decreasing slowly over the course of a search. By relaxing
the acceptance criterion in this way, it becomes possible to move out of the basin of attractionof a local optimum. Tabu search adopts a deterministic approach, whereby a ‘memory’ is implemented by the recording of previously-seen solutions. This record could be explicit, but isoften an implicit one, making use of simple but effective data structures. These can be thoughtof as a ‘tabu list’ of moves which have been made in the recent past of the search, and whichare ‘tabu’ or forbidden for a certain number of iterations. This prevents cycling, and also helps
to promote a diversified coverage of the search space. Perturbation methods improve the restart
strategy: instead of retreating to an unrelated and randomly chosen initial solution, the current
local optimum is perturbed in some way and the heuristic restarted from the new solution. Perhaps the most widely-known of such techniques is the ‘iterated Lin-Kernighan’ (ILK) method
introduced by Johnson [41] for the travelling salesman problem. On reaching a local optimum,a set of links is randomly chosen for removal and re-connection, in such a way that a new search
can start relatively close to a new local optimum. Such techniques can perhaps best be describedas perturbation methods. Genetic algorithms differ in using a population of solutions rather thanmoving from one point to the next. Furthermore, new solutions are generated from two (or,rarely) more solutions by applying a ‘crossover’ operator. However, they can also be encompassed within an NS framework, as we shall discuss later in this thesis.

6.3 The Neighborhood Structure for the Jobshop Scheduling
Problem
As shown in Section 5.1 and 5.2, the jobshop scheduling problem with n jobs and m machines
can be considered itself as a permutation problem; namely we have a permutation of n jobs
on each machine, which results in m-partitioned n-job permutations. However, the simple SH
and EX operators, for example, are not efficient because of the two reasons: (1) the size of theneighborhood becomes too large and (2) the resulting new permutation does not always correspond to a feasible schedule. One way to resolve these problems is to construct a neighborhood
structure based on Theorem 1 in Chapter 2. Namely, given a schedule S , a transition operator that exchanges a pair of adjacent critical operations (i.e., operations on a critical path) on asame machine in S as shown in Figure 6.1 forms a neighborhood which we call the AE (adjacent exchange) neighborhood and denote AE(S ). Theorem 1 guarantees that AE(S ) members

6.3. The Neighborhood Structure for the Jobshop Scheduling Problem

55

are always feasible and Theorem 2 guarantees that an optimal schedule is reachable from any
initial schedule by applying finite number of transitions. The transition operator was originallyproposed by Balas in his branch and bound approach[6] and has been used as a neighborhood
structure for SA in [12] and for TS in [17].
AE neighborhood

critical block

Figure 6.1: AE(S ), adjacent exchange neighborhood of S , consists of schedules
obtained from S by exchanging a pair of adjacent operations within a same critical
block.

Another very powerful transition operator was proposed in [9] using the notions of before
candidate and after candidate introduced in Section 2.5 of Chapter 2. Let a schedule be S and
let its critical blocks be B1 , . . . , Bk , then before candidate BBj and after candidate BAj in a critical
block B j are defined by Equation 2.5. Let NBB j (S ) and NAB j (S ) be sets of (maybe infeasible)
schedules obtained by moving each operation in BBj (or BAj ) to the front (or rear) of B j respectively
as shown in Figure 6.2. Because we have Theorem 2.5, it is tempted to define the CB neighorhood
as a set of all the shedules obtained from before and after candidates as follows:[{NBB j (S ) ∪ NAB j (S )}.
(6.3)
CB0 (S ) =Bj

However unfortunately, there is no theorem similar to Theorem 1 that guarantees the feasibilityof CB0 members. In fact, CB0 may contain infeasible schedules. Therefore the CB neighborhood
is given as follows:
CB neighborhood

u1

um

critical block

Figure 6.2: CB(S ), critical block neighborhood of S , consists of schedules obtained
from S by moving an operation in a critical block to the front or the rear of the block.

CB(S ) = {S 0 ∈ CB0 (S ) | S 0 is a feasible schedule}.

(6.4)

56

Chapter 6. Neighborhood Search

In the next chapter, we will see that using the CB neighorhood, an efficient simulated annealing
algorithm can be constructed.

Chapter 7
Critical Block Simulated Annealing for the
Jobshop Scheduling Problem
7.1 Simulated Annealing
Consider that we have a nonlinear optimization function f (x) defined over a continuous variable x
in multi-dimensional Euclidean space X. Such a nonlinear optimization function may be likened
to a mountainous state space landscape, with the algorithm’s objective being to locate the lowest
valley. Simulated Annealing (SA) methods1 are analogous to searching a state space landscapeby bouncing a rubber ball around the terrain. The ball bounces around the landscape, in andout of different valleys, probabilistically sampling different locations. As the degree of bounce
is reduced it becomes more difficult for the ball to bounce out of low valleys into higher ones,
than vice versa. Finally, when there is no bounce left in the ball, the ball will settle in the lowest
valley. This is mathematically guaranteed given time and a proper bounce (annealing) reduction
schedule [42].
Thus, SA algorithms employ noise to choose new parameter values. They generate a newstate x0 in the neighborhood of x, probabilistically. When x is a continuous variable, there areinfinite number of candidate states in the neighborhood of x and thus, a new state x0 is generatedusing a given distribution function g() which will be described shortly. The algorithms calculate
the value of the function cost E = f (x0 ), and then probabilistically decides to accept or reject it.
If accepted, the new state becomes the current state. The new state may be accepted even if it
has a larger function cost than the current state. The criteria for acceptance is determined by anacceptance function h(), the temperature parameter T , and the difference in the function valuesof the the two states. Initially, T is large, and a new state is accepted quite frequently. As thealgorithm progresses, T is reduced, lowering the probability that the acceptance function will
accept a new state if it’s functional cost is greater than that of the current state.
The general SA procedure [43] is defined below.
1. Choose an initial (high) temperature T 0 and a random state x0 .
1

Annealing (as in metallurgical annealing) refers to the process involving the slow reduction of a temperature.

57

58

Chapter 7. Critical Block Simulated Annealing for the Jobshop Scheduling Problem
T k=0 ← T 0 , x ← x0
2. Calculate the cost function value of the starting state.
Ek=0 ← f (x0 )3. For each iteration k, k = 1 . . . k f do the following:(a) Choose a new state x0 , using a generating function.
x0 ← g(x)(b) Calculate the cost of x0 .
E 0 ← f (x0 )(c) Set x ← x0 and E ← E 0 with probability determined by the acceptance function h().
(d) Reduce the temperature T by annealing(e. g. T k+1 ← γT k , 0 < γ < 1).
(e) When T is lower than a sufficiently small value T f , exit the loop.
4. Return x and E as the (near) optimal state and function cost value.

Because the algorithm occasionally chooses states uphill from its current state (i.e. chooses stateswith higher functional values than the current states), it can escape from local minima and more
effectively search the function space to find the global minimum.
The Simulated Annealing method in general consists of a system state x and the following
functional relationships:
1. f (x): a cost function to minimize,
2. g(x): a generating probability density function of new states.
3. P(x): an acceptance function that decides if a new state should become the current state,
and
4. T (k): an annealing temperature (T ) schedule.
For numeric optimization problems, x is often defined as an integer or real parameter vector,
x = {xi ; i = 1 . . . D}, and Boltzmann Annealing is used to generate new states. Boltzmann
Annealing employs a Gaussian probability density function,g(x0 ) =10 2e−(x−x ) /(2T )D/2(2πT )

where g(x0 ) is the probability of generating x0 from the currently accepted state x, and where thetemperature T is a measure of the fluctuations of the Boltzmann distribution.

7.2. Critical block Simulated Annealing

59

A Gaussian probability density distribution is not applicable for the generation of new states
for combinatorial optimization problems including the job shop scheduling problem. Instead, a
uniform random distribution is often usedg(x0 ) = 1/n, x0 ∈ N(x)(7.1)

where n is the number of states that can be directly generated by the generating function, i.e. n
is the number of states in the neighborhood of x, n = |N(x)|.
The acceptance probability function P(x) is based on the chances of accepting a new state x0
relative to the current state x, i.e. the difference of their function values0e− f (x )/T
1
P(x ) = − f (x0 )/T
=
.
−
f
(x)/T
(
f
(x
e
+e
1 + e 0 )− f (x))/T
0

(7.2)

If lower cost states are always accepted, as in [44], the acceptance function above can be redefined as(1
if
f (x0 ) ≤ f (x)
0
0
(7.3)
P(x ) =e(− f (x )− f (x))/T
otherwise.
The practical annealing schedule, T k , most often used to find the global minimum is of the form
T k = T 0 e−ck

(7.4)

where c is a positive constant.

7.2 Critical block Simulated Annealing
For the JSP, a state x is defined by a particular schedule S , and the cost f (x) is defined by the
makespan Cmax (S ). A neighborhood N(S ) of a schedule S can be defined as the set of feasible
schedules that can be reached from S by exactly one transition (a single perturbation of S ).
We use the critical block neighborhood CB(S ) defined by (6.4) in the previous chapter as theneighborhood structure.
The algorithm begins by setting the annealing temperature to an initial value and generating a random schedule S . The makespan and critical path of S is then calculated. Next, a newschedule S 0 in the neighborhood of S is randomly generated. The new schedule S 0 is comparedwith the current schedule S , and probabilistically accepted according to the makespan difference
between the two schedules, and the annealing temperature. The temperature, initially quite high,is decreased according to a given annealing schedule. This process is repeated until 1) the temperature is sufficiently low, 2) a given number of iterations have occurred, or 3) a schedule having
a (near) optimal makespan is found. Finally, the best generated schedule and its makespan are
printed. The algorithm is described in 7.2.1.

Chapter 7. Critical Block Simulated Annealing for the Jobshop Scheduling Problem

60

Algorithm 7.2.1 (The Critical Block Simulated Annealing Algorithm)
We are given a jobshop scheduling problem to optimize, and the initial temperature T 0 .
1. Set S = S 0 a randomly generated initial schedule, iteration step number k = 0, and T k=0 =T 0 the initial temperature.
2. do(a) doi. Pick S 0 ∈ CB(S )ii. Accept S 0 probabilistically according to the Metropolis Criterion distribution,0i.e. choose S 0 with probability one if Cmax (S 0 ) ≤ Cmax (S ), and e−(Cmax (S )−Cmax (S ))/T
otherwise, i.e., the probability P to accept S 0 is defined as follows:(1
if Cmax (S 0 ) ≤ Cmax (S )
0
0
(7.5)
P(S ) =
e−(Cmax (S )−Cmax (S ))/T
otherwise.
until S 0 is accepted.
(b) Set S = S 0 , increase k, and decrease T according to the annealing schedule.
until termination (i.e. T is sufficiently small or k is sufficiently large or the minimal
makespan or lower bound is found).
3. Print the best schedule found and its makespan.

7.3. Reintensification

61

7.3 Reintensification
Often during the search process 1) the system state wanders far from states that are leading to
the global minima, or 2) the system state may become trapped in a deep local minima. Then
all schedules generated from the current schedule will have longer makespans. If the acceptancetemperature is low, it may be difficult for the system to escape from this local minimum andcontinue the search.
Occasionally a reintensification strategy may be applied to improve the search. This reintensification is similar to reannealing [45], [46] which is used in numerical function optimization
to occasionally reset the system temperature and the state of the system after a sufficiently long
period of time without finding a new global minima.
It is useful to reintensify the search when a large number of acceptances have occurred without improving the problem makespan, or when the recent acceptance to generated ratio (AGratio) becomes lower than a prescribed threshold, indicating that the system is caught in a minimum. The reintensification process replaces the current state (schedule) with the best state foundso far, removing the system state from a local minimum if it has become trapped in a basinof attraction. Reintensification also alters the annealing temperature to a more current and appropriate value. A new annealing temperature is calculated from the standard deviation of thefunctional cost of states in the best neighborhood. If the new resulting temperature is greater thanthe current temperature, then the current temperature is reset to the new temperature.

7.4 Parameters
Simulated annealing algorithms often require some parameter specific values be determined apriori. These annealing scheduling parameters include the initial and final temperatures (T 0 , T f ),and the number of annealing steps (k f ). Reasonable values for the reintensification frequency
and the AG threshold must be chosen as well. For the annealing schedule, appropriate choicesof both the initial and final (lower bound) temperatures, and the maximum number of annealing
steps, must be determined.
An appropriate temperature reduction function is also needed. Since both inverse logarithmic
and inverse linear annealing schedules are too slow for practical consideration, it is useful toapply the exponential annealing schedule given in equation 7.4, with constant c determined byc = − log(T f /T 0 )/k f .

(7.6)

For example, if the annealing schedule is defined so that T 0 = 1, T f = 10−20 and k f = 10000,
then c = −log (10−20 )/10000 = 0.0046.
Since scheduling problems have different characteristics, constraints, and differing degrees
of difficulty, different annealing schedules must be chosen to fit different problems. Because theinitial and final lower bound temperatures are problem dependent parameters, they are difficult
to determine a priori. By defining these temperatures in terms of the desired uphill AG ratios, thetemperatures can be determined adaptively from problem independent values. The initial uphill
AG ratio should be relatively large so that a large number of uphill transitions are accepted.

62

Chapter 7. Critical Block Simulated Annealing for the Jobshop Scheduling Problem

Later, at the end of the annealing schedule, the final uphill AG ratio should be small forcing mostschedules having inferior makespans to be rejected.
The adaptive determination of the initial or final temperature can be incorporated in a shortwarmup sequence at the beginning of a simulation. A desired uphill AG threshold is chosen apriori which is used to adaptively determine the annealing temperatures. The reverse annealing
process is implemented by starting at a sufficiently low temperature such that no uphill states are
accepted, and increasing the temperature by a small percentage (e.g. 5%) until the actual uphill
AG ratio is greater than or equal to the desired threshold. This approach is similar to that definedby Aarts [47] (p59), however only uphill generated and accepted states determine the AG ratio.
For the initial temperature, T 0 , an uphill AG threshold of 50% of the makespans generated
that were larger than the current schedule’s makespan was found to be appropriate. For the
final (lower bound) temperature, T f , an uphill AG threshold of 0.2% proved most effective. T fprovides a lower bound for the final temperature. This value would only be arrived at if all
generated states would be accepted. Realistically this never occurs.
The total number of annealing steps k f can be chosen empirically, but should be governedby the desires 1) to have reasonable small differences between successive temperatures, and 2)
to have non-excessive trial run times, i.e. at most one or four hour per trial, and 3) to generate
as many potential schedules as possible within the time limits. Concerning to 1), the acceptancetemperature is assumed to be lowered according to equation 7.4 sufficiently slowly such that adetailed balance is maintained, and that the resulting distribution of the inhomogeneous Markov
chain generated using this schedule between temperatures T + δ and T − δ, (δ  1) approximates
the stationary distribution of a finite length homogeneous Markov chain at temperature T .
If reintensification is to be used, two parameters must be specified to determine when it
must be performed. The first parameter, the reintensification frequency R, determines how often
reintensification should be performed. Reintensification can be applied after a set number of newschedules are accepted without finding an improved makespan. The second parameter is an AG
threshold limit. Reintensification is performed when the current AG ratio falls below this value.

7.5 Methodology and Results
The performance of the CBSA algorithms was tested by running several simulation trials withand without reintensification. For the reintensification trials, two reintensification frequencies
R = 3, 000 and R = 10, 000 were tested, both with an AG ratio threshold of 10−3 . Although
reintensification violates the theoretical ergodicity of simulated annealing by resetting the stateof the system, performances were found to be improved when reintensification was incorporated
into the system.
Table 7.1 shows the minimum makespans of the first ten trials when the CBSA algorithm(with reintensification every 3,000 acceptances) was initially applied to the mt10 problem. Different random number seeds were used in each trial, resulting in each trial starting from a different randomly generated schedule. The CBSA algorithm was executed for a warm up period
to generate new schedules and to gather acceptance rate statistics. The statistics were used toadaptively determine appropriate values for T 0 and T f . T f given in equation 7.6 was used to7.5. Methodology and Results
Run
1
2
3
4
5
6
7
8
9
10

Min Evaluations Generations
930
481429548175?
930
510050537087?
930
507396579957?
930
344341331749?
930
366680403856?
930
459323472286?
930
371984405170?
930
649431711167?
930
316954352034938
459372
1000000?

63
Initial Temp Last Temp
Time51.837858
6.194249 38m 0s
44.833024 10.612196 41m 28s
47.036501
6.203178 40m 8s
49.404748 17.333848 28m 45s
44.760441
7.891267 28m 40s
47.085840 14.834713 37m 59s
38.693533
9.067814 29m 8s
38.756278
9.568809 51m 39s
36.828203
8.406492 25m 0s
54.385611
0.500000 36m 6s

Table 7.1: Ten Trials using the Simulated Annealing Method (R = 3, 000).

determine reasonable values for c.
The table also shows the number of actual schedules evaluated, the number of new schedules
generated, and the cpu time for each of the trials2 .
Since new schedules are often regenerated from the same current schedule, their makespansneed not be reevaluated. Hence the actual number of schedules evaluated is always less than orequal to the number of generations. Optimal schedules are indicated by a star ? .
Table 7.2 shows the initial and final temperatures of the ten trials. The last temperature of thesuccessful runs shown can vary considerably depending upon when the algorithm terminated. Intable 7.2 the algorithm was terminated when an optimal solution was found, or when k was equal
to k f . Since T f was determined from accepted, rather than generated states, any solution foundon or before the last generation k f will have an actual final temperature larger than the initiallydetermined T f .

2
One test simulation performed during the initial programming of the mt10 problem found an optimal schedulein 47 seconds on a Sparcstation 2, however this was not to be representative of other trials.

64

Chapter 7. Critical Block Simulated Annealing for the Jobshop Scheduling Problem
Run
1
2
3
4
5
6
7
8
9
10

Initial Temp Last Temp
51.837858
6.194249
44.833024
10.612196
47.036501
6.20317849.404748
17.333848
44.760441
7.891267
47.085840
14.834713
38.693533
9.067814
38.756278
9.568809
36.828203
8.406492
54.385611
0.5

Table 7.2: Initial and Last Temperatures. Last temperature is the temperature whenan optimal makespan was found, or the temperature after 1,000,000 iterations.

Although the last mt10 trial did not find the optimal solution, it did terminate with a makespanof 938. Trial ten’s initial temperature, which was adaptively determined, was the highest of all of
the trials. Hence it is likely that the system spent excessive amounts of annealing time performing random search at high temperatures. The cpu time and the number of function evaluations
performed during the execution of trial ten was quite comparable with those of the successful
trials. It is likely that the the system state became caught in a deep local minimum, allowing few
if any new states to be explored.

7.5.1 Random Search

The effects of randomly searching the schedule space was investigated to determine if the transition operations were solely responsible for the generation of the high quality schedules shownin table 7.1. Ten trials of the simulated annealing algorithm were performed by setting the initialand final temperatures to large values, i.e. T 0 = 100.0 and T f = 99.0. Results of this random
search are shown in table 7.3.

7.5. Methodology and Results

65

Run Min Generations
Time
Acceptances
1 994
1000000 1h 41m 16s
844798
2 998
1000000 1h 39m 34s
845644
3 997
1000000 1h 40m 9s
846223
4 999
1000000 1h 39m 13s
847081
5 993
1000000 1h 39m 15s
846005
6 998
1000000 1h 39m 13s
846835
7 992
1000000 1h 39m 30s
8453028 995
1000000 1h 40m 33s
8461639 989
1000000 1h 41m 19s
84568810 1010
1000000 1h 40m 36s
846533
Table 7.3: Ten High Temperature Random Trials.

Approximately 84% of all schedules generated were accepted. Schedules with shorter makespans
were always accepted, i.e. from equation 7.5, S 0 is always accepted when L(S 0 ) ≤ L(S ). Sincehigh temperatures, T 0 ≈ T f  1, result in a large numbers of inferior schedules being accepted,
the method essentially performs like a random search with the critical block transition operators
being used to generate the new schedules. Performances of these random searches was noticeably
poorer than those in table 7.1.

7.5.2 Low Temperature Greedy Search

Searching with a very low temperature for a small number of generations essentially implements
a greedy (downhill only) search. When the critical block transition operators are applied withthis greedy algorithm, many of the generated schedules of poor quality were observed. Ten
thousand schedules were generated by randomly generating initial schedules and applying the
CBSA algorithm to them for 1,000 iterations at very low temperatures, i.e. T 0 = 0.01, T f ≈0.001, and k f = 1000. Figure 7.1 shows the a histogram of the 10,000 makespans generated. It isclear from those performances that the low temperature greedy algorithm is not solely responsible
for the performances in table 7.1.

66

Chapter 7. Critical Block Simulated Annealing for the Jobshop Scheduling Problem
10000 Greedy Trials (MT10x10)
Trials100.00
95.00
90.00
85.00
80.00
75.00
70.00
65.00
60.00
55.00
50.00
45.00
40.00
35.00
30.00
25.00
20.00
15.00
10.00
5.00
0.00
Makespan x 103
1.00

1.10

1.20

1.30

1.40

Figure 7.1: Generated Makespans of 10,000 Greedy (mt10) Schedules.

Van Laarhoven et. al. [12] describe a similar method called iterative improvement that consists of repeated generation of random schedules using the same neighborhood structure. They
tested that method using the previously described Balas transition operator and the best makespan
they found over 5 macro runs (averaging 9441.2 trials each) was 1006. In contrast, the best costmakespan found during the ten thousand greedy CBSA trials was 944. The relative difference isindicative of the power of the respective transition operators.

Figure 7.2 shows the time evolution of the makespans of two typical trials with and without
reintensification. The abscissa shows the number of schedules generated, and the ordinate shows
the makespan differences between the current and optimal schedules. The highly oscillatory behavior of the reintensification trial is due to reintensification.

7.6. Performance on Benchmarks Problems
Annealing Method
CBSA min/max
CBSA mean/std.
AESA min/max
AESA mean/std.
GREEDY min/max
GREEDY mean/std.

67

R=0
930
945
937.80 4.19
938
972
951.60 10.20
971
1491
1171.45 66.15

R = 3, 000
930
938
930.80 2.40
930
951
939.50 5.12

R = 10, 000
930
938
933.10 3.81
934
970
944.40 10.24

Table 7.4: Comparisons between CBSA and AESA.
Time Evolution of a CBSA Trial (R = 3,000)

220

220

200

200

180

180

160

160
Makespan Difference

Makespan DifferenceTime Evolution of a CBSA Trial (R = 0)

140

140

120

120

100

100

80

80

60

60

40

40

20

20
0

0
0

200

400

600

800

1000
Generations x 10 3

0

200

400

600

800

1000
3

Generations x 10

Figure 7.2: Successive makespan differences between the current and optimal solution of the
mt10 problem, without reintensification (R=0) and with reintensification (R = 3, 000).
The power of a reintensification and Critical Block neighborhood structure is shown in table
7.5.2. We show comparative performances of 10 trials of the CBSA and the simulated annealing algorithms using the AE (adjacent exchange) neighborhood (AESA), which was describedin Section 6.3, proposed by van Laarhoven et. al. [12] . All performances conditions wereidentical, except for the reintensification frequencies, (R=0 (no reintensification), R=3,000, and
R=10,000), and the neighborhood structure. The performances are best when R=3,000 and the
CBSA is used. Not shown are the cpu times which were similar for all runs given.

7.6 Performance on Benchmarks Problems
A set of benchmark problems has been established to judge the effectiveness of algorithms on
the JSP. 3
3

We are grateful for the benchmark problem set given to us by D. Applegate [48].

68

Chapter 7. Critical Block Simulated Annealing for the Jobshop Scheduling Problem
PROB
abz7
abz8
abz9
la21
la24
la25
la27
la29
la38
la40

NxM
20x15
20x15
20x15
15x10
15x10
15x10
20x10
20x10
15x15
15x15

LB
655
638
656
1046
935
977
1235
1130
1196
1222

CBSA
685
679
701
1050
943
985
1262
1188
1209
1235

Appl
668
687
707
1053
?
935
?
977
1269
1195
1217
?
1222

Laar Mats

Adams Bruck
710
716
735
1063 1071
1084 1059
952
973
976
935
992
991
1017
977
1269 1274
1291 1270
1203 1196
1239 1202
1215 1231
1255 1232
1234 1235
1269 1238

Table 7.5: Ten difficult Benchmark Job Shop Scheduling Problems.

The problem set includes the problems mt06, mt10, and mt20 from [4], and problems car1car8 from [49]. Problems abz5-9 are those given in Adams [14]. Also included is a set of 40
job shop scheduling benchmark problems la01-la40 originally from [32], The first column inboth tables gives the problem name. The next column, NxM, indicates the size of the problem,
i.e. N jobs by M machines. The LB column indicates the lower bound of the problem if the
optimal makespan is unknown. The CBSA column indicates the best makespan found from 5
CBSA trials. (Each trial used a reintensification frequency R = 3, 000 and was run for one
million generations or until the known optimal minimum makespan or lower bound was found.)
The column headings Appl, Laar, Mats, Adams and Bruck indicate the best performances ofApplegate [1], Van Laarhoven [12], Matsuo [13], Adams [14], and Brucker [9] respectively. A
single star, ? , indicates the optimum or best known minimum. Table 7.6 shows the performancesof the CBSA and some of the best known job shop algorithms on the la problem test set.

7.6. Performance on Benchmarks Problems
PROB
la01
la02
la03
la04
la05
la06
la07
la08
la09
la10
la11
la12
la13
la14
la15
la16
la17
la18
la19
la20
la21
la22
la23
la24
la25
la26
la27
la28
la29
la30
la31
la32
la33
la34
la35
la36
la37
la38
la39
la40

NxM
10x5
10x5
10x5
10x5
10x5
15x5
15x5
15x5
15x5
15x5
20x5
20x5
20x5
20x5
20x5
10x10
10x10
10x10
10x10
10x10
15x10
15x10
15x10
15x10
15x10
20x10
20x10
20x10
20x10
20x10
30x10
30x10
30x10
30x10
30x10
15x15
15x15
15x15
15x15
15x15

LB
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
—
1046
—
—
935
977
—
1235
—
1130
—
—
—
—
—
—
—
—
1196
—
1222

CBSA
?
666
?
655
?
597
?
590
?
593
?
926
?
890
?
863
?
951
?
958
?
1222
?
1039
?
1150
?
1292
?
1207
?
945
?
784
?
848
?
842
907
1050
935
?
1032
943
985
?
1218
1262
?
1216
1188
?
1355
?
1784
?
1850
?
1719
?
1721
?
1888
1291
1420
1209
1243
1235

69
Appl
?
666
?
655
?
597
?
590
?
593
?
926
?
890
?
863
?
951
?
958
?
1222
?
1039
?
1150
?
1292
?
1207
?
945
?
784
?
848
?
842
?
902
1053
?
927
?
1032
?
935
?
977
?
1218
1269
?
1216
1195
?
1355
?
1784
?
1850
?
1719
?
1721
?
1888
?
1268
?
1397
1209
?
1233
?
1222

Laar
666
655
606
590
593
926
890
863
951
958
1222
1039
1150
1292
1207
956
784
861
848
902
1063
938
1032
952
992
1218
1269
1224
1203
1355
1784
1850
1719
1721
1888
1293
1433
1215
1248
1234

Mats Adams Bruck
—
666
666
655
669
655
597
605
597
590
593
590
—
593
593
—
926
926
—
890
890
863
863
863
—
951
951
—
958
958
—
1222 1222—
1239 1039—
1150 1150—
1292 1292—
1207 1207
959
978
945
784
787
784
848
859
848
842
860
842
907
914
902
1071
1084 1059
927
944
927
1032
1032 1032
973
976
935
991
1017
977
1218
1224 1218
1274
1291 1270
1216
1250 1276
1196
1239 1202
1355
1355 1355—
1784 1784—
1850 1850—
1719 1719—
1721 1721—
1888 1888
1292
1305 1268
1435
1423 1424
1231
1255 1232
1251
1273 1233
1235
1269 1238

Table 7.6: Performances of the 40 Benchmark Job Shop Scheduling Problems.

70

Chapter 7. Critical Block Simulated Annealing for the Jobshop Scheduling Problem

7.7 Concluding Remarks
According to Aarts [47], The success of an approximation algorithm depends on a number of
aspects including performances, ease of implementation, and applicability and flexibility. It isclear that the SA algorithm is very simple and easy to implement, taking only a few hundredlines of code to implement. Regarding flexibility, we found that our SA code could be used
to implement both Aarts SA approach and the CBSA method by changing only the call to theneighborhood generation procedure.
Matsuo [13] has indicated that Aarts SA method has very few adjacent pairs that improve themakespan by exactly one interchange (Aarts SA transition). Compared with Aarts neighborhood
structure, the CB neighborhood contains more (or at least the same number of) transitions that
can immediately improve the a schedules makespan by exactly one transition.
Undoubtedly new and more powerful approaches will be developed to solve the JSP, methods having better performance and convergence characteristics then the methods described here.
Unlike SA methods, other approximation approaches have no theoretical guarantee that they will
converge to to an optimal solution.

Chapter 8
Critical Block Simulated Annealing with
Shifting Bottleneck Heuristics
In this chapter, we consider to improve the CBSA described in the previous chapter in two
aspects. One is to employ a new neighborhood by incorporating the notion of active schedule
described in Section 2.2 and the other is to combine with a deterministic heuristic called “shiftingbottleneck” described in Section 2.6.

8.1 Active Critical Block Simulated Annealing
As explained in the previous chapter, a solution obtained from a before or an after candidate
is not necessarily executable. In the following, we propose a new neighborhood by modifying
the CB neighborhood. Each element in the new neighborhood is not only executable, but also
active and structurally close to the corresponding member in the original CB neighborhood. Let
S be an active schedule and Bα,βγ be a critical block of S on a machine Mγ , where the first andα,βthe last operations of Bγ are the α-th and the β-th operations on Mγ respectively. Let Oλγ bea “moving” operation that is the λ-th operation on Mγ such that α < λ < β (i.e., Oλγ is inside
λ,β
λ,α
λ
Bα,β
γ ). We consider to generate an active schedule S γ (or S γ ) by moving Oγ into position α (orposition β). If the resulting schedule is active, then we will use it. Otherwise, we try to find analternate position that is as close to α (or β) as possible and is inside Bα,βγ such that the resulting
schedule becomes active and use it instead. This can be done by adopting the GT algorithm
described in Algorithm 2.2.1 and modifying Step 4 of the algorithm, where the choice of thenext operation from the conflict set C[Mr , i] was at random. Here, the choice is made by looking
at the processing order in S . An operation that is in C[Mr , i] and that was scheduled earliest in S
is selected when the operation is located outside Bα,βγ (i.e., i < [α, β]). This way, the processingλ,βλ,αorder of operations in S γ or S γ is kept mostly unchanged from that in S outside Bα,βγ . Howeverα,βλ
inside Bγ instead, the “moving” operation Oγ must be chosen from the conflict set with the top
most priority when generating S γλ,α and with the bottom least priority when generating S γλ,β . The
details are described in Algorithm 8.1.1. In the algorithm, if we are generating S γλ,α (i.e., d = f )and when i = α, Oλγ should be chosen from C[Mr , i] as soon as i becomes equal to α. However71

72

Chapter 8. Critical Block Simulated Annealing with Shifting Bottleneck Heuristics

at that point, Oλγ may not yet exist in C[Mr , i]. This means that if we move Oλγ to the front
position α ignoring the fact that Oλγ < C[Mr , α], then the resulting schedule becomes non-active,
or may even become infeasible. Hence, we have to “pass” this time and wait until Oλγ appearesin C[Mr , i]. Because S is an active schedule, it is guaranteed that Oλγ appeares in C[Mr , i] at thelatest when i = λ, in which case resulting schedule becomes identical to S . If we are generating
S γλ,β instead (i.e., d = r), then Oλγ should be moved to the rear postion β, in other words, Oλγ
should be chosen from C[Mr , i] only when i becomes equal to β and not while i < β. In fact, Oλγ
is guaranteed to appear in C[Mr , i] when i = β at the latest because S is active. However Oλγ may
become the only element of C[Mr , i] and then it is unavoidable to choose Oλγ even when i < β.
This means that if we move Oλγ to the rear position β ignoring the fact that Oλγ < C[Mr , β], thenthe resulting schedule becomes non-active, or may be even infeasible. Hence we have to choose
Oλγ even when i < β. The new neighborhood ACB(S ) is now defined as a set of all S γλ,α and S γλ,βover all critical blocks:


[



λ,α
λ,β
ACB(S ) =  {S γ }α<λ<β ∪ {S γ }α<λ<β  \ {S }.
(8.1)
 α,β

Bγ

Once we have the neighborhood structure ACB(S ) defined, the basic framework of the simulated annealing algorithm described in Algorithm 7.2.1 can be applied without major modifications by using ACB(S ) in place of CB(S ) in Step 2(a)i of Algorithm 7.2.1.
In Algorithm 7.2.1, if the acceptance probabilities are low for all members in ACB(S ), thesystem will remain trapped in a local minimum S and it will take a long time to move to a newstate. The algorithm may stay in S even after all members are selected in Step 2(a)i and evaluatedin Step 2(a)ii. To avoid this, relative acceptance probability defined by (8.1)
P(S i ) = P

P(S i )
S j ∈ACB(S ) P(S j )for each S i ∈ ACB(S )(8.2)

will be introduced after all members in ACB(S ) are visited and evaluated without being chosen.
The member S i in ACB(S ) is then randomly selected in proportion to P(S i ), and the systemmoves unconditionally to S i . This modification is effective when, like ACB(S ), the neighborhood
size is limited. Therefore, Algorithm 7.2.1 is modified as in Algorithm 8.1.2.

8.2 Active CBSA Enhanced by Shifting Bottleneck
As described in Section 2.6, Shifting bottleneck (SB) proposed by [14] is a powerful heuristic
method for solving a JSP. Here we consider to incorporate the local optimization procedure usedin the SBI heuristic described in Algorithm 2.6.1 into the active CBSA described in the previous
section.
SBI is a constructive method that generates a complete schedule from scratch. Modifying
the method is necessary in order to refine a certain complete schedule for improvement. The
BottleRepair shown in Algorithm 8.2.1 describes an iterative version of the basic SB. The reoptimization process used here is the same as used in Algorithm 2.6.1. The basic idea of BottleRepair

8.2. Active CBSA Enhanced by Shifting Bottleneck73

Algorithm 8.1.1 (An algorithm to generate S γλ,α or S γλ,β from S )A scheduling problem is given as in Algorithm 2.2.1. An original active schedule S of the
problem and a moving operation Oλγ on a critical block Bα,βγ of S is given, where α < λ < β. Theλ,αdirection d = f (front, when generating S γ ) or d = r (rear, when generating S γλ,β ) is given.
1. Initialize G (as in Step 1 of Algorithm 2.2.1). Initialize S 0 as an empty schedule.
2. Find the earliest completable operation O∗r ∈ G (as in Step 2 of Algorithm 2.2.1).
3. Calculate the conflict set C[Mr , i] ⊂ Gr (as in Step 3 of Algorithm 2.2.1).
4. If r , λ, then select and schedule from C[Mr , i] the earliest operation Okr in S as the i-th
operation on Mr in S 0 , i.e., S ri0 := k, where k is the job number of Okr .
If r = λ,• If d = f (i.e., generating S γλ,α ),– If α ≤ i ≤ λ and Oλγ ∈ C[Mr , i], then select and schedule Oλγ as the i-th operation
on Mr in S .
– Otherwise, select and schedule from C[Mr , i] the earliest operation in S as thei-th operation on Mr in S 0 .
• If d = r (i.e., generating S γλ,β ),– If λ ≤ i ≤ β and C[Mr , i] contains any operation other than Oλγ , then select andschedule from C[Mr , i] \ {Oλγ } the earliest operation in S as the i-th operation on
Mr in S 0 .
– If C[Mr , i] = {Oλγ } or i = β and Oλγ ∈ C[Mr , i] , then select and schedule Oλγ as thei-th operation on Mr in S 0 .
– Otherwise, select and schedule from C[Mr , i] the earliest operation in S as thei-th operation on Mr in S 0 .
Let Okr be the operation selected above, where k is the job number of Okr , then
S ri0 := k.
5. Do Step 6 and Step 7 of Algorithm 2.2.1.
6. Repeat from Step 1 to Step 5 until all operations are scheduled.
7. Output the active schedule S 0 as S γλ,α if d = f , or as S γλ,β if d = r.

Chapter 8. Critical Block Simulated Annealing with Shifting Bottleneck Heuristics

74

Algorithm 8.1.2 (The Active Critical Block Simulated Annealing)
We are given a jobshop scheduling problem to optimize, and the initial temperature T 0 .
1. Set S = S 0 a random initial active schedule, generated by the GT algorithm in Algorithm 2.2.1. Set k = 0, and T k=0 = T 0 .
2. do(a) Set N be the size of ACB(S ) and set n = 0, the number of elements in ACB(S ) thatare already evaluated.
(b) doi. Pick S i from ACB(S ) randomly and if S i is first time to be picked and not yet
evaluated, then evaluate S i by calculating Cmax (S i ) and set n = n + 1.
ii. Accept S i with probability:(1
if Cmax (S i ) ≤ Cmax (S )P(S i ) =−(Cmax (S i )−Cmax (S ))/Teotherwise.
until S i is accepted or n = N.
(c)

• If S i is accepted, then set S = S i and n = 0.
• Otherwise (i.e., if n = N) select S 0 from ACB(S ) in proportion to the probability
P(S i ) = P

P(S i )
S j ∈ACB(S ) P(S j )and set S = S i and n = 0.
(d) Set k = k + 1 and decrease T k according to the annealing schedule.
until termination (i.e. T is sufficiently small or k is sufficiently large or the minimal
makespan or lower bound is found).
3. Print the best schedule found and its makespan.

8.2. Active CBSA Enhanced by Shifting Bottleneck

75

Algorithm 8.2.1 (BottleRepair: Iterative SB)
A complete selection S is given as an input. Let M be a set of all machines: M = {M1 , . . . , Mr }.
1. Let M0 = {M1 , M2 , . . . , Ml } be a set of all the critical machines (machines that containparts of critical path in S ).
2. Reset all sequences of all machines in M \ M0 and make the machines unsequenced. Let
S p be a partial selection obtained from S by this resetting.
3. Reoptimize all sequenced machines in M0 by applying Step 6 of Algorithm 2.6.1 with M0
and S p obtained above and obtain new S p .
4. Do Step 2 to Step 8 of Algorithm 2.6.1 with M0 and S p obtained above and obtain acomplete selection.
5. Output the complete selection as an obtained schedule.

Algorithm 8.2.2(SB incorporation for Active CBSA)

2(b)iii. If S i is rejected, apply BottleRepair to S i and obtain S i ∗ .
Accept S i ∗ and set S i = S i ∗if Cmax (S i ∗ ) < Cmax (S ).

comes from the original paper of SB [14] where the last α noncritical machines are temporarily
removed for the reoptimization.
As shown in Algorithm 8.1.2, S i is selected from ACB(S ) and is probabilistically accepted.
BottleRepair is applied to S i only when S i is rejected. The resulting schedule S i ∗ is accepted if
its makespan is shorter than that of S . To summarize, Step 2(b)iii as defined in Algorithm8.2.2
is added to Algorithm 8.1.2 just after Step 2(b)ii.
BottleRepair gives a systematic way to inspect the schedule’s critical path and permutes operations again and again by repeatedly solving one machine problems in a deterministic manner.
If it generates an improved schedule S 0 from S , the critical path of S 0 becomes different from
S and the difference is much greater than that between S and its active CBSA neighbor. On theother hand, active CBSA gives a stochastic more focused local search around the current critical
path. The proposed integration of active CBSA and SB is expected to have the synergistic effectas: SB gives a long jump to active CBSA so that it can omit many time-consuming inferior transitions and active CBSA adds stochastic perturbations to SB so that it can escape from the local
minima.

76

Chapter 8. Critical Block Simulated Annealing with Shifting Bottleneck Heuristics

8.3 Experimental Results
8.3.1 Muth and Thompson’s Benchmark
A 10 × 10 problem (mt10) and 20 × 5 problem (mt20) formulated by [4] (MT benchmarks) are
well known benchmark JSPs. CBSA with and without SB modification was evaluated using
these problems. Table 8.1 shows the results of 20 trials with different random number seeds on a
SUN SPARC station 10. All programs are written in the C language.
Table 8.1: Comparisons between CBSA and CBSA+SB using MT benchmarks
CBSA (R = 6, 000)
CBSA+SB
Prob n × m
bestmean std BT best
mean std BT
mt10 10 × 10 930 933.65 4.04 190 930 932.45 3.01 786
mt20 20 × 5 1178 1179.45 1.94 235 1165 1165.00 0.00 449
std: standard deviation
BT: average cpu time (sec.) to find the best solution

Results for the mt10 problem using CBSA without SB show that the optimal solutions of
L = 930 were found in 11 trials. The average cpu time was 3 min. 10 sec., and the fastest was 1
min. 21 sec. Although the solutions of L = 930 were found in only half of the trials, the cpu timein successful runs were satisfactorily short. If the temperature is more slowly lowered, though it
takes longer, the rate of finding optimal solutions will become higher as in [50]. CBSA without
SB could not find any optimal solution for mt20 problem. In most cases, solutions of L = 1178
were found instead of the optimal L = 1165.
The results for the mt10 problem using active CBSA with SB modification show that the
number of trials finding optimum solutions increased slightly, but the average cpu time increased
about four times. This fact indicates that CBSA without SB is powerful enough to solve mt10
problem. On the other hand, all 20 trials with SB modification for mt20 problem found the
optimal solutions of L = 1165 in an average cpu time of 7 min. 29 sec., and 1 min. 22 sec. wasthe best time. The effect of SB enhancement is obvious from this problem. Reintensification did
not work well because the optimal or near optimal solutions were obtained at an early stage ofthe search.

8.3.2 Other Benchmarks
Results in the previous section indicate that if CBSA without SB can solve a problem skillfully,
applying CBSA+SB has no advantage. However, if CBSA fails to work well, CBSA+SB may
improve solution quality and compensate for the extra cpu time needed with SB enhancement. A
set of benchmark problems has been established to evaluate different algorithms for JSPs. Table
8.2 shows the makespan performances of CBSA+SB and various other algorithms for the ten
difficult benchmark JSPs. All experiments of CBSA+SB runs were done on a HP 730 (HP 730
is about 1.5 times faster than SUN SPARC station 10).

8.3. Experimental Results

77

The LB column indicates the theoretical lower bound of the problem if the optimal makespanis unknown. The CBSA+SB column indicates the best makespans found from ten trials. Eachtrial used a reintensification frequency of R = 1, 000 and ran for three hours or until the known
optimal makespan or lower bound was found. The column headings Aart, Matt, Appl and Tailindicate the best performances of those in [51], [52], [1] and [17] respectively.
Table 8.2: Results of 10 tough JSPs
CBSA+SB
Prob n×m LB
Aart Matt Appl Tailbest mean std. BT
abz7 20×15 654 665 671.0 3.92 7814 668 672 668 665
abz8 20×15 635 675 680.0 3.13 8775 670 683 687 676
abz9 20×15 656 686 698.6 7.42 8749 691 703 707 691
la21 15×10 1040 1046 1049.3 3.32 361 1053 1053 1053 1047
la24 15×10 –935 939.2 1.99 6226 935 938 935
la25 20×10 –977 979.3 1.62 4117 983 977 977
la27 20×10 1235 1235 1242.4 6.15 7805 1249 1236 1269 1240
la29 20×10 1120 1154 1162.4 7.10 5434 1185 1184 1195 1170
la38 15×15 1184 1198 1206.8 4.53 3479 1208 1201 1209 1202
la40 15×15 – 1228 1230.2 2.32 3331 1225 1228 1222

Initial Schedule: L = 1665
1635

1585

1535

Makespan1485

T = 17.19

1435

Temperature (T)1385

Makespan (L)T = 6.24

1335

L = 1291

1285

1235
0

T = 0.02

L = 1269 4

6

L = 1235

83
Cpu time x 10

Figure 8.1: The time evolution of CBSA+SB trial for the la27 problem
Figure 8.1 shows the time evolution of the makespan (L) and temperature (T ) of the best
trial of CBSA+SB for the la27 problem. The abscissa shows the cpu time in seconds, and the

78

Chapter 8. Critical Block Simulated Annealing with Shifting Bottleneck Heuristics

Table 8.3: An optimal solution of la27 problemm
Job sequences on each machine
1 18 9 4 2 11 3 12 8 10 7 13 15 20 1 6 19 5 14 17 16
2 14 9 8 12 3 17 15 4 10 13 6 7 11 19 1 2 18 16 20 5
3 11 14 3 19 10 4 9 5 6 2 16 17 13 15 7 20 18 12 8 1
4 6 11 10 1 19 3 8 13 15 18 2 12 14 17 4 7 9 16 20 5
5 12 3 5 6 20 13 11 8 1 17 9 19 7 14 2 16 18 15 10 4
6 19 18 11 8 14 7 3 16 2 17 1 10 9 5 13 6 20 4 12 15
7 17 2 13 15 9 8 19 7 6 20 10 18 14 16 5 4 3 1 11 12
8 2 12 16 14 11 10 3 5 19 20 6 7 13 17 8 18 15 1 4 99 18 12 3 14 11 10 16 8 19 13 4 6 15 2 17 20 1 9 7 5
10 13 20 18 9 4 5 11 6 3 8 2 19 10 14 1 17 7 15 16 12
m: machinesolid and dotted lines show the makespan and the temperature respectively. Starting from the
makespan value L = 1665, it rapidly decreases to L = 1291 during the first warmup interval.
After twelve times of reintensification, it finally reached the optimal value L = 1235. In thisexperiment, 56059 schedules were generated and 8850 schedules accepted. About 25% of theaccepted schedules were accepted by BottleRepair to add CBSA long jumps and the rest served
as stochastic perturbations to BottleRepair. The oscillatory behavior is due to reintensification.
Although CBSA+SB outperformed other methods most of the cases in Table 8.2, the required
computational time is much longer. For example, [53] reported that Applegate’s method in theAppl colum found a schedule of L = 1269 in 604.2 sec. on SUN SPARC station ELC which isabout 10 times slower than HP 730. This is because each CBSA+SB experiment includes a lotof unsuccessful BottleRepair trials. But this gap can be filled to some extent by the fact that in
the same experiment, CBSA+SB passed a point L = 1269 in 164 sec.

8.4 Concluding Remarks
The proposed method CBSA+SB is an improved CBSA enhanced by integrating with a problem
specific method called shifting bottleneck. The performance of CBSA+SB was evaluated using
difficult benchmark problems. The results show that for eight problems of ten difficult benchmark
problems, CBSA+SB could find schedules better than or equal to the best schedules published so
far in the literature, when enough computational time is given. A new solution of L = 1235 was
found for the la27 problem; it is optimal because the value equals the theoretical lower bound.
Further research is necessary to reduce the computational time.

Chapter 9
Scheduling by Genetic Local Search with
Multi-Step Crossover Fusion
As we have seen in Chapter 5, Genetic Algorithms can be applied to the job-shop schedulingproblem with good success. However, the well-known fact that GAs are, in general, not well
suited for fine-tuning structures which are very close to optimal solutions also applies to thiscase and obstructs further improvements. The general remedy to this problem is to incorporate
local search methods, such as neighborhood search described in Chapter 6, into GAs. The resultof such incorporation is often called Genetic Local Search (GLS) [26]. In this framework, anoffspring obtained by a recombination operator, such as crossover, is not included in the next
generation directly but is used as a “seed” (initial solution) for the subsequent local search. The
local search moves the offspring from its current point to the nearest locally optimal point, whichis included in the next generation.
In solving combinatorial optimization problems, it is often difficult to construct an efficient
crossover operator, because a crossover operator that “recombines” solutions, allows to cause
global changes that alters the structure of a solution in large parts and therefore may violate theconstraints of the target problem, resulting in generating many infeasible solutions.
On the other hand, a neighborhood search operator which modifies a solution only locally israther easier to construct as we have seen in Chapter 7 and Chapter 8. A simple neighborhood
search operator exchanges a pair of consecutive jobs in a job sequence, another operator removes
a job from its original position and re-insert it in another position on the same job sequence. As
we have seen in the previous chapters, they should be improved by focusing on the jobs on thecritical path. Unfortunately the same method cannot be applied directly to construct an effective
crossover operator.

9.1 Multi-step crossover fusion
Reeves has been exploring the possibility of integrating local optimization directly into a Simple GA with bit string representations and has proposed the Neighborhood Search Crossover(NSX) [54]. Let any two individuals be x and z. An individual y is called intermediate between
79

Chapter 9. Scheduling by Genetic Local Search with Multi-Step Crossover Fusion

80

Algorithm 9.1.1 (Multi-Step Crossover Fusion — MSXF)
Parent individuals p0 , p1 are given as inputs. The distance measure d is defined. For a given
individual x, f (x) is the function to be minimized.
1. Initialize both current solution x and output q as p0 ; x = q = p0 .
2. do(a) For each member yi ∈ N(x), estimate d(yi , p1 ).
(b) Sort yi ∈ N(x) in ascending order of the d(yi , p1 ) estimation.
doi. Select yi from N(x) randomly, but with a bias in favor of yi with a small index i.
ii. Evaluate yi and obtain f (yi ) if yi has not yet been evaluated.
iii. Accept yi with probability one if f (yi ) ≤ f (x), and with Pc (yi ) otherwise.
iv. (optional) Change the index of yi from i to n, and the indexes of yk (k ∈ {i+1, i+
2, . . . , n}) from k to k − 1.
until yi is accepted.
(c) Set x = yi .
(d) If f (x) < f (q) then set q = x.
until some termination condition is satisfied.
3. Output q as an individual for the next generation.

9.1. Multi-step crossover fusion

81

x and z, written as x  y  z, if and only if d(x, z) = d(x, y) + d(y, z) holds, where x, y and z
are represented in binary strings and d(x, y) is the Hamming distance between x and y. Thenthe kth -order 2 neighborhood of x and z is defined as the set of all intermediate individuals at a
Hamming distance of k from either x or z. Formally,
Nk (x, z) = {y | x  y  z and (d(x, y) = k or d(y, z) = k)}.
Given two parent bit strings p0 and p1 , the neighborhood search crossover of order k (NSXk ) will
examine all individuals in Nk (p0 , p1 ), and pick the best as the new offspring.
In this chapter, we extend the idea of the NSX to make it applicable to more complicatedproblems such as job-shop scheduling and propose the Multi-Step Crossover Fusion (MSXF): a
new crossover operator with a built-in local search functionality [55, 56, 57]. The MSXF has thefollowing characteristics compared to the NSX.
• It can handle more general class of representations (i.e., it is not limited to the binary
representation) and neighborhood structures.
• It is based on a stochastic local search algorithm.
• Instead of restricting the neighborhood by the intermediateness, a biased stochastic replacement method is used.
A stochastic local search algorithm is used for the base algorithm of the MSXF. Although the
SA is a well-known stochastic method and has been successfully applied to many problems aswell as to the JSP, it would be unrealistic to apply SA repeatedly in a GA run that would consume
too much time. In such a case, a restricted method with a fixed temperature parameter T = c is agood alternative. Accordingly, the acceptance probability defined in (7.1) is modified as:(1
if
f (x0 ) ≤ f (x)0(9.1)
P(x ) =(− f (x0 )− f (x))/ceotherwise.
Let parent individuals be p0 and p1 , and let the distance between any two individuals x andy in any representation be d(x, y). If x and y are schedules, then d(x, y) is the DG distancedefined in Section 2.4. Let f be the function to be minimized. In the case of the schedulingproblem, f is the makespan Cmax . The basic idea to incorporate crossover functionality into theneighborhood search described in Algorithm 6.1.1 is to set the initial starting point as one ofthe parent: x0 = p0 and when choosing a candidate y from N(x), give a greater acceptance bias
for y with small d(y, p1 ), which we call biased stochastic replacement. This bias control in the
MSXF is achieved easily by sorting N(x) members in ascending order of d(yi , p1 ) so that yi witha smaller index i has a smaller distance d(yi , p1 ). The calculation of d(yi , p1 ) is not an expensive
task if d(x, p1 ) and the the nature of the move from x to yi are known; it is not necessary toactually generate and evaluate yi . Then yi is chosen from N(x) randomly, but with a bias for yi
with small index i. The outline of the MSXF is described in Algorithm 9.1.1.
In place of d(yi , p1 ), one can also use sign(d(yi , p1 ) − d(x, p1 )) + r to sort N(x) members in
Algorithm 9.1.1. Here sign(x) denotes the sign of x: sign(x) = 1 if x > 0, sign(x) = 0 if x = 0,

82

Chapter 9. Scheduling by Genetic Local Search with Multi-Step Crossover Fusion

sign(x) = −1 otherwise. A small random fraction r is added to randomize the order of memberswith the same sign. The termination condition can be given, for example, as the fixed number of
iterations in the outer loop.
The MSXF is not applicable if the distance between p0 and p1 is too small compared to thenumber of iterations. In such a case, a mutation operator called the Multi-Step Mutation Fusion(MSMF) is applied instead. The MSMF can be defined in the same manner as the MSXF is
except for one point: the bias is reversed, i.e. sort the N(x) members in descending order ofd(yi , p1 ) in Algorithm 9.1.1.

9.2 Scheduling in the reversed order
The GT algorithm in Algorithm 2.2.1 and all its variants determine the job sequences from left
to right in temporal order. This is because active schedules are defined to have no extra idle
periods of machines prior to their operations. However the idea described below enables thesame algorithms to determine the job sequences from right to left with only small modifications.
In general, a given problem of the JSP can be converted to another problem by reversing allof the technological sequences. The new problem is equivalent to the original one in the sensethat reversing the job sequences of any schedule for the original problem results in a schedule forthe reversed problem with the same critical path and makespan. It can be seen, however, that an
active schedule for the original problem may not necessarily be active in the reversed problem;the activeness is not necessarily preserved.
job
Routing
1 2(3) 1(4)2 2(2) 1(4)Figure 9.1: A simple 2 × 2 problem

(1)

(2)

M1

J2

J1

J2

L=13

M1

J2

J1

M2 J2

J1

M2

J1

(4)

M1
M2

J1

J2

L=13

J2

J1

(3)

M1
L=10

J1

J2

M2 L=10J1

J2

Figure 9.2: Schedule reversal and activation
For example, the simple 2 × 2 problem described in Table 9.1 is considered. Figure 9.2(1)shows a solution of this problem, which is active and no more left shifts can improve its makespan.

9.3. MSXF-GA for Job-shop scheduling

83

Figure 9.2(2), obtained by reversing Figure 9.2(1), is not active and can be improved by a left
shift that moves job 1 prior to job 2 on machine 2, resulting in Figure 9.2(3). Finally Figure 9.2(4)is obtained by reversing Figure 9.2(3) again, which is optimal. As things turn out, Figure 9.2(1)
is improved by moving job 1 posterior to job 2 on machine 2, resulting in Figure 9.2(4).
Although repairing a semi-active schedule to the active one improves the makespan, it canbe seen from the example above that there sometimes are obvious improvements that cannot
be attained only by left shifts. We call a schedule left active if it is an active schedule for theoriginal problem and right active if it is such for the reversed problem. It sometimes happensthat a reserved problem is easier to solve compared to the original. Searching only in the setof left (or right) active schedules may bias the search toward the wrong direction and result inpoor local minima. Therefore left active schedules as well as right active ones should be taken
into account together in the same algorithm. In most local search methods, many schedules are
generated in a single run; therefore it would be better to apply this reversing and repairing method
periodically to change the scheduling directions rather than to reverse and repair every scheduleeach time it is generated.

9.3 MSXF-GA for Job-shop scheduling
The MSXF is applied to the JSP by using the active CB neighborhood and the DG distancedefined in Section 2.4. Algorithm 9.3.1 describes the outline of the MSXF-GA routine for the
JSP using the steady state model proposed in [58, 59]. To avoid premature convergence evenunder a small-population condition, an individual whose fitness value is equal to someone’s in
the population is not inserted into the population in Step 2d.
A mechanism to search in the space of both the left and right active schedules is introduced
into the MSXF-GA as follows. First, there are equal numbers of left and right active schedules
in the initial population. The schedule q generated from p0 and p1 by the MSXF ought to beleft (or right) active if p0 is left (or right) active, and with some probability (0.1 for example) the
direction is reversed.
Figure 9.3 shows all of the solutions generated by an application of (a) the MSXF and (b) astochastic local search computationally equivalent to (a) for comparison. Both (a) and (b) started
from the same solution (the same parent p0 ), but in (a) transitions were biased toward the othersolution p1 . The x axis represents the number of disjunctive arcs whose directions are different
from those of p1 on machines with odd numbers, i.e. the DG distance was restricted to odd
machines. Similarly, the y axis representing the DG distance was restricted to even machines.

9.4 Benchmark Problems
The two well-known benchmark problems with sizes of 10 × 10 and 20 × 5 (known as mt10
and mt20) formulated by Muth and Thompson [4] are commonly used as test beds to measure
the effectiveness of a certain method. The mt10 problem used to be called a “notorious” problem, because it remained unsolved for over 20 years; however it is no longer a computational

Chapter 9. Scheduling by Genetic Local Search with Multi-Step Crossover Fusion

84

Algorithm 9.3.1(MSXF-GA for the JSP)

1. Initialize population: randomly generate a set of left and right active schedules in equal
number and apply the local search to each of them.
2. do(a) Randomly select two schedules p0 , p1 from the population with some bias dependingon their makespan values.
(b) Change the direction (left or right) of p1 by reversing the job sequences with probability Pr .
(c) Do step 2(c)i with probability Pc , or otherwise do Step 2(c)ii.
i. If the DG distance between p1 , p2 is shorter than some predefined small value
dmin , apply MSMF to p1 and generate q.
Otherwise, apply MSXF to p1 , p2 using the active CB neighborhood N(p1 ) andthe DG distance and generate a new schedule q.
ii. Apply Algorithm 6.1.1 with acceptance probability given by (9.1) and the active
CB neighborhood.
(d) If q’s makespan is shorter than the worst in the population, and no one in the population has the same makespan as q, replace the worst individual with q.
until some termination condition is satisfied.
3. Output the best schedule in the population.

(a) Multi-Step Crossover Fusion

50(b) Stochastic Local Search

50

Best (957)
40

40

30

30

20

P0
Start (961)10Target (951)
P1

0
0

10

20

P0
Start (961)10Best (935)

Target (951)
P1

0
20

30

40

50

60

0

10

20

30

40

50

60

Figure 9.3: Distribution of solutions generated by an application of (a) MSXF and (b) a shortterm stochastic local search9.4. Benchmark Problems

85

Table 9.1: Performance comparison using the MT benchmark problems
1963
1991
1992

1995
1995
1996

Muth-Thompson
Nakano/Yamada
Yamada/Nakano
Dorndorf/Pesch
Dorndorf/Pesch
Kobayashi/Ono/Yamamura
Bierwirth
Yamada/Nakano

Test problems
Simple GA
Giffler-Thompson GT-GA
Priority-Rule based P-GA
Shifting-Bottleneck SB-GA
Subsequence Exchange Crossover
SXX-GA
Generalized-Permutation GP-GA
Multi-step Crossover Fusion MSXF-GA

10 × 10
965
930
960
938
930

20 × 5
1215
1184
1249
1178
1178

936
930

1181
1165

challenge.
Applegate and Cook proposed a set of benchmark problems called the “ten tough problems”
as a more difficult computational challenge than the mt10 problem, by collecting difficult problems from literature, some of which still remain unsolved [1].

9.4.1 Muth and Thompson benchmark
Table 9.1 summarizes the makespan performance of the methods described in this chapter. The
Simple GA described in Chapter 4 has only limited success. It would be improved by beingcombined with the GT algorithm and/or the schedule reversal. The other results excluding the
MSXF-GA results are somewhat similar to each other, although the SXX-GA (a GA with Subsequence Exchange Crossover described in Section 5.1) is improved over the GT-GA describedin Section 5.4 in terms of speed and the number of times needed to find optimal solutions for the
mt10 problem. The SB-GA produces better results using the very efficient and tailored shiftingbottleneck procedure. The MSXF-GA which combines a GA and local search obtains the best
results.
For the MSXF-GA, the population size = 10, constant temperature c = 10, number of iterations for each MSXF = 1000, Pr = 0.1 and Pc = 0.5 are used. The MSXF-GA experiments were
performed on a DEC Alpha 600 5/226 which is about four times faster than a Sparcstation 10,
and the programs were written in the C language. The MSXF-GA finds the optimal solutions forthe mt10 and mt20 problems almost every time in less than five minutes on average.

9.4.2 The Ten Tough Benchmark Problems
Table 9.2 shows the makespan performance statistics of the MSXF-GA for the ten difficult benchmark problems proposed in [1]. The parameters used here were the same as those for the MT
benchmark except for the population size = 20. The algorithm was terminated when an optimal
solution was found or after 40 minutes of cpu time passed on the DEC Alpha 600 5/266. In thetable, the column named lb shows the known lower bound or known optimal value (for la40)

86

Chapter 9. Scheduling by Genetic Local Search with Multi-Step Crossover Fusionof the makespan, and the columns named bst, avg, var and wst show the best, average, variance and worst makespan values obtained, over 30 runs respectively. The columns named noptand topt show the number of runs in which the optimal schedules are obtained and their average
cpu times in seconds. The problem data and lower bounds are taken from the OR-library [60].
Optimal solutions were found for half of the ten problems, and four of them were found very
quickly. The small variances in the solution qualities indicate the stability of the MSXF-GA as
an approximation method.
Table 9.3 shows the performance comparisons with various heuristic methods for the 10
tough problems. The column headings Nowi and Dell indicate the best performances of TABU
search proposed in [19] and [18]. CBSA+SB indicates SA results described in Chapter 8. Aart,
Matt, Appl indicate SA results proposed in [51], GA results in [52] and results in [1] respectively.
It is interesting to observe from those results that the two approaches CBSA+SB are mutuallycomplementary. For the problem instances abz7, abz8, abz9 and la29, SA+SB outperforms
MSXF-GA. In fact, SA+SB performs the best for these problems. However, for the problem
instance la38 for which SA+SB fails to find the global optimum, MSXF-GA successfully andfrequently find the global optimum. MSXF-GA also shows reasonably good performances forthe problem instance la40 for which SA+SB performs relatively poorly.
Table 9.2: Results of the 10 tough problems
prob sizelb
bst avg var wst nopt
topt
abz7 20×15 655
678 692.5 0.94 703 ––abz8 20×15 638
686 703.1 1.54 724 ––abz9 20×15 656
697 719.6 1.53 732 ––?
la21 15×10 –1046 1049.9 0.57 1055 9 687.7?
la24 15×10 –935 938.8 0.34 941 4 864.1?
la25 20×10 –977 979.6 0.40 984 9 765.6?
la27 20×10 –1235 1253.6 1.56 1269 1 2364.75
la29 20×10 1130 1166 1181.9 1.31 1195 ––la38 15×15 – ? 1196 1198.4 0.71 1208 21 1051.3
la40 15×15 ? 1222 1224 1227.9 0.43 1233 ––Figure 9.4 shows a performance comparison of GLS with and without MSXF using the la38
problem. A total of 100 experiments (runs) were performed for each under the same conditions
used in Table 9.2 but with different random seeds. In the figure, the solid line gives the results of
MSXF-GA (in other words, GLS with MSXF) and the dotted line gives the equivalent results of
GLS without MSXF (i.e., in place of MSXF, CPU equivalent short-term stochastic local searchis used). Each of the 100 runs is numbered from No.1 to No.100 in ascending order of cpu
time at which each run is terminated. For example, run No.1 successfully found the optimalschedule and was terminated the most quickly. The cpu time value = 2400 means that the run
was terminated before it found the optimal schedule. In the figure, instead of standard time
evolution graph, the x axis represents run numbers and the y axis represents the cpu time. The9.4. Benchmark Problems

87

Table 9.3: Performance comparisons with various heuristic methods on the 10 tough problems
prob our bst Nowi Dell CBSA+SB Aarts Matt Appl
abz7
678
– 667
665 668 672 668
abz8
686
– 678
675 670 683 687
abz9
697
– 692
686 691 703 707?
la21 ? 1046 1047 1048
1046 1053 1053 1053?
la24 ? 935 939 941
935 ? 935 938 ? 935?
la25 ? 977 ? 977 979
977 983 ? 977 ? 977?
?
la27 1235 1236 1242
1235 1249 1236 1269
la29 1166 1160 1182
1154 1185 1184 1195
la38 ? 1196 ? 1196 1203
1198 1208 1201 1209
la40 1224 1229 1233
1228 1225 1228 ? 1222

2500

CPU time (sec.)

2000

without MSXF

1500
with MSXF
1000

500

0
0

10

20

30

40

50

60

70

80

90
100
No. of runs

Figure 9.4: Performance comparison using the la38 15 × 15 problem

88

Chapter 9. Scheduling by Genetic Local Search with Multi-Step Crossover Fusion

fact that the solid line increases slower and has shorter horizontal tail part than the dotted line
means that the experiments with MSXF outperforms those without MSXF both in terms of the
cpu time and in the number of successful runs.

Chapter 10
Permutation Flowshop Scheduling by
Genetic Local Search
So far, we have mainly considered the jobshop scheduling problem (JSP), but hereafter, we will
shift our focus to the permutation flowshop scheduling problem, abbreviated as PFSP or just
FSP, which is a special case of the JSP in a sense that the technological sequence of machines
is the same for all jobs and the order in which each machine processes the jobs is also same forall machines. A solution schedule is then represented by a permutation of n jobs, instead of m
permutations of n jobs. We do not have the concept of active schedule, but still we have critical
path and blocks, if the objective is Cmax but not if the objective is C sum , the sum of the completion
times of all the operations. In this chapter, the MSXF method described in the previous chapter
is applied to the n/m/P/Cmax and in the next chapter, the MSXF is applied to the n/m/P/C sum .

10.1 The Neighborhood Structure of the FSP
Because the FSP is a special case of the JSP, the concepts of critical path and blocks also apply
to the FSP in a similar but much simpler fashion. We briefly review those concepts again in the
FSP context for the better understanding.
The permutation flowshop scheduling problem designated by the symbols n/m/P/Cmax has
n jobs that have to be processed on m machines in the same order. Because of this simplicity
against the general JSP, we can assume that the machines are indexed in the processing orderof jobs, so that the M1 is the first machine to process jobs and Mm is the last. We are given theprocessing time p jr of each operation of job J j on machine Mr . A schedule can be represented
by a permutation of jobs, or to put it more simply, a permutation of job numbers π.
A critical path is a sequence of operations starting from the first operation on the first machine
M1 and ending with the last operation on the last machine Mm . The starting time of each operationon the path, except for the first one, is equal to the completion time of its preceding operation—that is, there is no idle time along the path. Thus, the length of the critical path is the sum of theprocessing times of all the operations on the path and equals to Cmax . There can be more than
one critical path on a schedule.
89

90

Chapter 10. Permutation Flowshop Scheduling by Genetic Local Searchjob 4
block machine
1
1

2

2

3

5

6

1

2

3

8

7

4
3

5

4

6

Figure 10.1: A grid graph representation of a solution to a problem of 8 jobs and 6 machines.

The operations on a critical path can be partitioned into subsequences, called critical blocks,
according to their associated machines. A critical block consists of maximal consecutive operations on the same machine, as in the JSP, but here we make it simple and define a critical block
as a subsequence of associated jobs, instead of operations.
Consider a schedule represented by a permutation π. Let B1 , . . . Bk be a set of all critical
blocks that contains more than one job and let ml be the index of the machine associated with Bl .
Let π(ul ) be the first job of Bl (and the last job of Bl−1 ). Then the ‘inside’ of Bl , denoted by B̂l , isdefined as follows:
B \ {π(ul+1 )}if l = 1 and ml = 1 l
Bl \ {π(ul )}if l = k and ml = n
B̂l = (10.1)

 B \ {π(u ), π(u )} otherwise.
l

l

l+1

Figure 10.1 shows an example of schedule π = 4, 5, 6, 1, 2, 3, 8, 7 for a problem with n = 8
jobs and m = 6 machines represented by so-called a grid graph, which is a simplification of
the disjunctive graph for the JSP. In the figure, the vertical axis corresponds to machines andthe horizontal axis to jobs. Each circle represents an operation, and arrows precedence relation
between operations. A critical path is marked by thick lines. In this example, there are four
critical blocks B1 , . . . , B4 that contain more than one job. B2 on machine 3, for example, consistsof four jobs 5, 6, 1 and 2, and B̂2 consists of jobs 6 and 1. Likewise B̂4 on machine 6 consists ofjobs 8 and 7. Note that the machines are indexed in the processing order of jobs.
As described in Chapter 6, a neighborhood N(x) of a point x in a search space can be defined
as a set of new points that can be reached from x by exactly one transition or move (a single
perturbation of x). One of the well-known transition operators for PFSP is the shift move thattakes a job from its current position and re-inserts it in another position. Let v = (a, b) be a pair of
positions in π. Here, v defines a move that removes the job π(a) from a position a and re-inserts
it in a position b. If a < b, the resulting schedule is represented by πv = π(1), . . . , π(a − 1)π(a +1), . . . , π(b), π(a), π(b + 1), . . . , π(n), and if a > b, πv = π(1), . . . , π(b), π(a), π(b + 1), . . . , π(a −1)π(a + 1), . . . , π(n). A neighborhood N(V, π) is defined as the set of all schedules obtained byshift moves in V = {(a, b) : b < {a − 1, a}, a, b ∈ {1, . . . , n}}.

10.2. Representative Neighborhood

91

s jb
b

Sj

the best move is selected
as a representative moves ja

the job moved to
the next blockj
the job moved tothe previous block

S ja

Figure 10.2: The best move to the next/previous block is selected as a representative.
Let Wl (π) be a set of moves restricted to the inside of Bl , namely
Wl (π) = {(a, b) ∈ V|a, b ∈ B̂l }andW(π) =(10.2)k[

Wl (π)(10.3)

l=1

then the block property described in Theorem 2.5 for the JSP is reformulated as follows:
Corollary 1 (Block property for the FSP) For any schedule β,
if β ∈ N(W(π), π) then Cmax (β) ≥ Cmax (π).
According to Corollary 1 above, no move in W(π) can directly improve schedule π. Therefore,
it is reasonable, for computational efficiency, to reduce the size of the neighborhood N(V, π) byeliminating moves in W(π), and to use a new neighborhood N(V \ W(π), π), which we call herea “critical block neighborhood”.

10.2 Representative Neighborhood
Nowicki and Smutnicki have proposed the representative neighborhood method[20], where a reduced neighborhood is generated from the original neighborhood by first patitioning the originalneighborhood members into clusters and then picking up the best move (representative) fromeach cluster as a representative. A new neighborhood is the set of all representative moves.
In this chapter, the representative neighorhood method is applied as follows. For each job
j in a critical block, let S aj be a set of moves that shift the job j to some position in the next
block; similarly S bj shifts j to the previous block. Evaluate schedules obtained from each movein S aj and denote the best one by saj . Similarly sbj is obtained from S bj . Then the representative
neighborhood is defined as a set of all schedules obtained by representative moves {saj , sbj } for alljobs j in all critical blocks (see Figure 10.2).

92

Chapter 10. Permutation Flowshop Scheduling by Genetic Local Search

10.3 Distance Measures
To measure the difference between two permutation schedules S and T , an appropriately defined distance is required. In the case of the JSP, we have used the DG distance introduced in
Section 2.4. For the FSP, two well-known distances are considered as follows:
precedence-based: This distance counts the number of job pairs {i, j} in which j is preceded byi in S but not in T .
position-based: This distance sums up the positional differences for each job in S and T .
The first precedence-based distance is equivalent to the DG distance if the FSP is viewed as a
special case of the JSP. The relationship between these two distances will be discussed in thelater section.

10.4 Landscape analysis
According to Höhn and Reeves [61], a landscape is defined by a triple of a search space, anobjective function and a distance measure. The link between landscape and search algorithm isgiven by the NS operators used in the algorithm. Because these operators generate new points inthe search space relative to a given point, they define a distance dN (s, t) on the search space givenby the minimum number of applications of operator N that will convert element t into elements.
One can understand the degree of difficulty of the given combinatorial optimization problem
by looking at its landscape: if the landscape is simple and has only one peak, it is very easy
to find the global optimum by using simple best ascent search. Unfortunately most NP-hardcombinatorial optimization problems, including PFSP, have very ‘rugged’ landscapes with many
false peaks under any NS operator.
Recently, Boese et al. [62] have shown that an appropriate choice of NS operator introducessome neat structure into the landscape. In this ‘big valley’ structure, local optima occur in clusters– good candidate solutions are usually to be found ‘fairly close’ to other good solutions. If
a landscape has this structure, it would support the idea of generating new starting points forsearch from a previous local optimum rather than from a random point in the search space.
Before we apply our GLS method to PFSP, we investigate whether there is a big valley
structure for the PFSP and the NS operator using the representative neighborhood and a stochastic
search using (9.1) (constant temperature). For the same PFSP but with simpler NS operators,
similar experiments reported in Reeves [63] found such a landscape did occur.
As discussed in [62, 64, 63], the existence of a big valley structure can be examined byfirst generating a set of random local optima and then observing the correlation between theirobjective function values and their distances to the nearest global optimum, and/or their average distances to other local optima. The distance used here should be dN for an operator N.
However this distance is difficult to compute, and precedence-based distance is used here as an
approximation.

10.4. Landscape analysis

93

ta011(a)

OBJFN
140
130
120
110
100
90
80
70
60
50
40
30
20
10
0
40

50

60

70

80

MEAND
90

ta021(a)

OBJFN

80

90

BESTD
0

20

40

100

MEAND
110

60

80

100

ta021(b)

OBJFN
210
200
190
180
170
160
150
140
130
120
110
100
90
80
70
60
50
40
30
20
10
0

200
190
180
170
160
150
140
130
120
110
100
90
80
70
60
50
40
30
20
10
0
70

ta011(b)

OBJFN
200
190
180
170
160
150
140
130
120
110
100
90
80
70
60
50
40
30
20
10
0

20

40

60

80

100

120

BESTD
140

Figure 10.3: 1841 distinct local optima obtained from 2500 short term local search for the ta011(20 × 10) problem and 2313 distinct local optima for the ta021 (20 × 20) problem are plottedin terms of (a) average distance from other local optima and (b) distance from global optima(x-axis), against their relative objective function values (y-axis).

94

Chapter 10. Permutation Flowshop Scheduling by Genetic Local Search

STEPS
28
26
24
22
20
18
16
14
12
10
8
6
4
2
0
0

20

40

60

80

PREC
100

Figure 10.4: The correlation between the precedence-based distance (PREC) and the approximate number of steps (STEPS)
Figure 10.3 shows a scatter plot of random local optima for problems ta011 and ta021, being
respectively the first of Taillard’s 20 × 10 and 20 × 20 groups of problems [21]. Each local
optimum is generated by running the neighborhood search described in Algorithm 6.1.1 with
L = 5000 based on the stochastic method with acceptance probability Pc , c = 5. Extensive
preliminary experiments found only two distinct global optima for the ta011 problem, very close
to each other in terms of the precedence-based distance (the distance is two) and only one global
optimum for ta021 problem; although one cannot rule out the possibility of finding other different
global optima by continuing the search. However, more than 2500 global optima were found forthe smaller ta001 (20 × 5) problem by spending the same amount of CPU time.
The x-axis in Figure 10.3 represents (a) the average precedence-based distance from otherlocal optima (MEAND), and (b) the precedence-based distance from one of the nearer global
optima (BESTD). The y-axis represents their objective function values relative to the global
optimum. These plots clearly show that there are good correlations between the distances andobjective function values. The calculated correlation coefficients for each plot are: ta011(a):0.74, ta011(b): 0.50, ta021(a): 0.62 and ta021(b): 0.44. These values are statistically significant
at the 0.1% level, on the basis of 1000 replications in a randomization test [63]. These highcorrelations suggest that the local optima are radially distributed in the problem space with theglobal optima as the centre, and the more distant are the local optima from the centre, the worse
are their objective function values. Hence, by tracing local optima step by step, moving from one
optimum to nearby slightly better one, without being trapped, one can eventually reach a near
global optimal solution.
In the analysis above, the precedence-based distance is used as a surrogate for dN , because
the minimum number of steps for the neighborhood operator to reach the global optimum isdifficult to compute. Although the precedence-based distance seems to be a good alternative, theapproximation still need to be justified. For this purpose, the approximate number of steps to10.5. MSXF-GA for PFSP

95

POSN
140
130
120
110
100
90
80
70
60
50
40
30
20
10
0
0

20

40

60

80

PREC
100

Figure 10.5: The correlation between the precedence-based distance (PREC) and the positionbased distance (POSN)
reach the global optimum from each local optimum was calculated by choosing the closest moveto the global optimum each time from the neighborhood. While this does not necessarily givethe best distance between two points, it seems likely to give a fairly close upper bound.
Figure 10.4 shows the correlation between the precedence-based distance (PREC) and theapproximate number of steps (STEPS) for the local optima shown in Figure 10.3 ta011(a) (correlation coefficient is 0.66). Figure 10.5 shows that there is a strong correlation between theprecedence-based distance (PREC) and the position-based distance (POSN) for the same local
optima (correlation coefficient is 0.91). Thus it does not matter which distance is used. The samekind of experiments were carried out for all Taillard’s 20×10 and 20×20 benchmarks, and similar
results were obtained in every case. Therefore, the use of the easily-computed precedence-based
distance appears to be justified, and the ‘big valley’ structure can be assumed to hold for thisneighborhood.

10.5 MSXF-GA for PFSPAs described in Section 9.1, the MSXF operator is designed to find a new local optimum basedon previous ones. MSXF-GA provides a framework for traversing local optima without being
trapped, by concentrating its attention on the area between the parent solutions and thus eventually finding a very good solution under the assumption of a ‘big valley’. MSXF-GA was applied
to PFSP using the representative neighborhood described in Section 10.1 and the precedencebased distance.
Algorithm 10.5.1 describes the outline of the MSXF-GA routine for the PFSP. In this model,
the population is ranked according to the makespan values, and the parents are selected fromthe population with a probability inversely proportional to their ranks. The newly generated

Chapter 10. Permutation Flowshop Scheduling by Genetic Local Search

96

solution q is inserted into the population only if its makespan is better than the worst in thecurrent population. To avoid premature convergence even under a small-population condition, if
an individual with the same makespan already exists in the population, then q is not inserted into
the population in Step 2c.
Algorithm 10.5.1 MSXF-GA for the PFSPThe population size P is given.
1. Initialize population: randomly generate a set of permutation schedules. Sort the population
members in descending order of their makespan values.
2. do(a) Select two schedules p1 , p2 from the population with a probability inversely proportional
to their ranks.
(b) Do Step 2(b)i with probability PX , or otherwise do Step 2(b)ii.
i. If the precedence-based distance between p1 , p2 is less than dmin , apply MSMF to p1
and generate q.
Otherwise, apply MSXF to p1 , p2 using the representative neighborhood and theprecedence-based distance and generate a new schedule q.
ii. Apply Algorithm 6.1.1 with acceptance probability Pc and the representative neighborhood.
(c) If q’s makespan is less than the worst in the population, and no member of the current
population has the same makespan as q, replace the worst individual with q.
until some termination condition is satisfied.
3. Output the best schedule in the population.

10.6 Experimental results
In Section 10.4, the existence of a big valley structure became clear for the relatively small-size
PFSP instances. An adaptive multi-start method (AMS) in which new local search is concentrated in a region between previously found local optima should be effective at least for these
problems. MSXF-GA for PFSP is especially designed as one of the AMS approaches for PFSP.
Preliminary experiments show that MSXF-GA is very effective for the problem instances discussed in Section 10.4, and the global optima are found very quickly. In this section we will
extend our investigations to larger-size problems and apply MSXF-GA to a subset of Taillard’s
benchmark problems.
Table 10.1 summarizes the performance statistics of MSXF-GA for a subset of Taillard’s
benchmark problems together with the results found by Nowicki and Smutnicki using their tabu10.6. Experimental results

97

Offspring

Parent2

Parent1

MSXF

Figure 10.6: Navigated local search by MSXF-GA: A new search is started from one of the parents and while no other good solutions are found, the search ‘navigates’ towards the other parent.
In the middle of the search, good solutions would be eventually found somewhere between theparents. That direction is then pursued to the top of a hill (or a bottom of the valley, if it is aminimization problem) — a new local optimum.

Table 10.1: Results of the Taillard benchmark problems
50 × 20
100 × 20
200 × 20
No. best avg. nowi lb – ub best avg. nowi lb – ub
best avg. nowi
lb – ub
1 3861 3880 3875 3771–3875 6242 6259 6286 6106–6228 11272 11316 11294 11152–11195
2 3709 3716 3715 3661–3715 6217 6234 6241 6183–6210 11299 11346 11420 11143–11223
3 3651 3668 3668 3591–3668 6299 6312 6329 6252–6271 11410 11458 11446 11281–11337
4 3726 3744 3752 3631–3752 6288 6303 6306 6254–6269 11347 11400 11347 11275–11299
5 3614 3636 3635 3551–3635 6329 6354 6377 6262–6319 11290 11320 11311 11259–11260
6 3690 3701 3698 3667–3687 6380 6417 6437 6302–6403 11250 11288 11282 11176–11189
7 3711 3723 3716 3672–3706 6302 6319 6346 6184–6292 11438 11455 11456 11337–11386
8 3699 3721 3709 3627–3700 6433 6466 6481 6315–6423 11395 11426 11415 11301–11334
9 3760 3769 3765 3645–3755 6297 6323 6358 6204–6275 11263 11306 11343 11145–11192
10 3767 3772 3777 3696–3767 6448 6471 6465 6404–6434 11335 11409 11422 11284–11313
best, avg.: our best and average makespan values
nowi: results of Nowicki and Smutnicki
lb, ub: theoretical lower bounds and best known makespans taken from OR-library

98

Chapter 10. Permutation Flowshop Scheduling by Genetic Local Searchsearch implementation[20] and the lower and upper bounds, taken from the OR-library [60].
(Upper bounds are the currently best-known makespans, most of them found by a branch andbound technique with computational time unknown). In all, 30 runs were completed for eachproblem under the same conditions but with different random number seeds. For each MSXFGA run, population size = 15, constant temperature c = 3, number of iterations for each MSXF
= 1000, dmin = n/2 and PX = 0.5 are used. Each run is terminated after 700 iterations, whichtakes about 12, 21 and 47 minutes of CPU time respectively for each 50 × 20, 100 × 20 and
200 × 20 problems on a DEC Alpha 600 5/226.
It can be seen that the results for 50 × 20 problems are remarkable: the solution qualities ofour best results are improved over those found in [20] for most of the problems, and some results(marked in bold letters) are even better than the existing best results reported in the OR-library.
The results for larger problems are not as impressive as those of 50 × 20 problems, but still goodenough to support our hypothesis. The degradation is probably due to the increasing complexityof the neighborhood calculation. In fact for problems where the ratio n/m > 3, Nowicki and
Smutnicki abandoned their representative neighborhood and used a simple one instead: justmoving a job to the beginning or the end of its critical block. They also implemented an efficient
way of evaluating all the members in the neighborhood in a specific order. This method is usefulfor the tabu search, but not directly applicable to our stochastic search.

10.7 Concluding Remarks
The landscape for the Permutation Flowshop Scheduling Problem with stochastic local searchand the representative neighborhood structure has been investigated. The experimental analysis
using 20 × 10 and 20 × 20 Taillard benchmark problems shows the existence of a ‘big valley’
structure for PFSP. This suggests a well-designed AMS method, such as MSXF-GA in which
new local search is concentrated in a region between previously found local optima should be
effective in finding near-optimal solutions. MSXF-GA for the PFSP is implemented using theneighborhood operator and applied to more challenging benchmark problems. Experimental
results demonstrates the effectiveness of the proposed method.

Chapter 11
C sum Permutation Flowshop Scheduling by
Genetic Local Search
11.1 Introduction
We have already proposed an efficient method based on the genetic local search with the MSXF
to solve the n/m/P/Cmax in the previous chapter. In this chapter we deal with the n/m/P/C sum .
Compared to the Cmax problem, the C sum problem is more difficult to optimize, mainly because thecalculation of the objective function is more time consuming, and problem specific knowledge
such as critical blocks is not applicable. However, this difficulty can be partly overcome byextending the idea of the representative neighborhood discussed in Section 10.2.

11.2 Representative Neighborhood
Let s be a job sequence of the current solution, and s[i, k] be a new job sequence obtained from s
by moving a job from the ith position in s and re-inserting it in the kth position. Nia (s) and Nib (s),subsets of N(s), are defined as follows:
Nia (s) = {s[i, k] | i < k ≤ n}, Nib (s) = {s[i, k] | 1 ≤ k < i}.
Thus the original neighborhood N(s) is divided into clusters consisting of Nia (s) and Nib (s). Let
Nia (s) and Nib (s) be one of the best members in Nia (s) and Nib (s) respectively. The representative
neighborhood N ∗ (s) can be denoted as:N ∗ (s) = {Nia (s) | 1 ≤ i < n} ∪ {Nib (s) | 1 < i ≤ n}.
In most local search algorithms, the NS operators to choose a new member from the neighborhood of the current solution can be categorized into two types according to their choice criteria
discussed in Section 6.1; one is best descent, and the other is first descent. The best descentmethod scans all the members in the neighborhood and choose the best as a new current solution. This is suitable when the neighborhood size is small and the cost of evaluating all the
99

100

Chapter 11. C sum Permutation Flowshop Scheduling by Genetic Local Search

S
ll
ea )
t
a
u S
val n Ni(ei
(2)

Ni (S)(3) choose the best as
a representative(1) choose Ni(S) at random
Figure 11.1: Representative neighborhood
members is negligible. Tabu search can be seen as an extension of this method. The first descentmethod selects one member (at random) and accepts it if it is sufficiently good, otherwise selects
another one. This can be used even when the neighborhood size is large. The stochastic sampling in stochastic local search including SA can be seen as an extension of this method. The
representative neighborhood fills the gap between these two criteria: a cluster Ni (s) is chosenrandomly by using first descent, then best decent is applied to evaluate all the members in Ni (s)of which the best is chosen as a representative. Figure 11.1 illustrates this process. As we will
see in the later section, this enables the TS and a stochastic local search method to integrate intoa single unified method.

11.3 Tabu List Style Adaptive Memory
As described in Chapter 6, Tabu Search (TS) adopts a deterministic local search approach witha ‘memory’ implemented as a ‘tabu list’ of moves which have been made in the recent past of
the search, and which are ‘tabu’ or forbidden for a certain number of iterations. The use of therepresentative neighborhood makes it easy to have ’tabu list’ style adaptive memory.
If a solution s[i, k] generated from the current solution s by moving a job j = s[i] to the kth
position is accepted, the pair ( j, i), i.e. the job and its original position, is stored on the top of a listof length l and recorded as tabu. The oldest element in the list is then deleted. In the subsequent
iterations, a solution generated by moving job j to the ith position should not be accepted as longas ( j, i) is on the list. In our representative neighborhood scheme, this is achieved by excluding
the tabu solution from the calculation of the representative best solution (Nia (s) and Nib (s)) sothat the representative neighborhood will contain no tabu solutions. Because the tabu solutions
have already been excluded from the representative neighborhood, there is no need to modify thestochastic local search procedure described in Section 10.5 by this modification.

11.4. Experimental Results

101

11.4 Experimental Results
We applied our method to some of Taillard’s benchmark problems (ta problems, in short) [21].
First it was applied to relatively easy problems from ta001 to ta030 (the number of jobs is 20
and the number of machines is 5, 10 and 20: denoted by 20x5, 20x10 and 20x20). Six runs
were carried out for each problem with different random seeds. The parameters used in theseexperiments are: P = 5, L1 = 1000, L2 = 700, PX = 0.5 and the length of the tabu list l = 7.
Here quite consistent results were obtained, i.e. almost all of the 6 runs converged to the samejob sequence in a short time (from a few seconds to a few minutes) before the limit of L2 = 700
was reached on a HP workstation. The best results (and they are also the average results in mostcases) are reported in Table 11.1 together with the results obtained by the constructive method(NSPD) due to J.Liu [65].
Table 11.1: Taillard’s benchmark results (ta001 – ta030)
prob
001
002
003
004
005
006
007
008
009
010

best
14033
15151
13301
15447
13529
13123
13548
13948
14295
12943

NSPD
14281
15599
14121
15925
13829
13420
13953
14235
14552
13054

prob
011
012
013
014
015
016
017
018
019
020

best
20911
22440
19833
18710
18641
19245
18363
20241
20330
21320

NSPD
21520
23094
20561
18867
19580
20010
19069
21048
21138
22212

prob
021
022
023
024
025
026
027
028
029
030

best
33623
31587
33920
31661
34557
32564
32922
32412
33600
32262

NSPD
34119
32706
35290
32717
35367
33153
33763
33234
34416
33045

Problems ta031 to ta050 (50x5 and 50x10 problems) are much more difficult and the best
results were different in each run. Ten runs were carried out for each problem with different
random seeds. The parameters used in these experiments were: P = 30, L1 = 10000, L2 =700, PX = 0.5. It takes 45 minutes per run for 50x5 problems (ta031 to ta040) and 90 minutes
for 50x10 problems (ta041 to ta050).
It is difficult to say how good these solutions are, in other words, how far they are from theglobal optima. Even for the easier problems in Table 11.1, there is no guarantee that the bestsolutions obtained so far are optimal, although we believe that they are at least very close tobeing so. For the problems in Table 11.2, it is almost certain that our best results are not optimal.
In fact we found one solution of C sum = 64803 for problem ta031 by an overnight run.

11.5 Concluding Remarks
A new genetic local search method is proposed to solve the C sum permutation flowshop scheduling problem. This method effectively integrates the stochastic sampling of Simulated Annealing,Chapter 11. C sum Permutation Flowshop Scheduling by Genetic Local Search

102

Table 11.2: Taillard’s benchmark results (ta031 – ta040)
prob
031
032
033
034
035
036
037
038
039
040

best
64860
68134
63304
68259
69491
67006
66311
64412
63156
68994

average64934.868247.2
63523.2
68502.769619.6
67127.666450.064550.1
63223.8
69137.4

NSPD
66590
68887
64943
70040
71911
68491
67892
66037
64764
69985

prob
041
042
043
044
045
046
047
048
049
050

best
87430
83157
79996
86725
86448
86651
89042
86924
85674
88215

average87561.4
83305.880303.4
86822.4
86703.7
86888.089220.787180.5
85924.3
88438.6

NSPD
90373
86926
83213
89527
89190
91113
93053
90614
91289
91622

the adaptive memory using the tabu list of Tabu Search and the population-based search of Genetic Algorithms into a single unified framework as summarized in Figure 11.2. The method isapplied to Taillard’s benchmark problems. Experimental results demonstrate the effectiveness ofthe proposed method.

Simulated Annealing(Stochastic local search)

Stochastic
sampling

representative
neighbourhood

first descent

AdaptivememoryTabu Searchbest descent+mult-step
crossover fusion

Genetic Local Search+

Population,
path relinking
crossover, mutation

Figure 11.2: The framework of the proposed method

Chapter 12
Tabu Search with a Pruning Pattern Listfor the Flowshop Scheduling Problem
12.1 Introduction
In this chapter, an approximation method based on Tabu Search with an additional memory
structure called “pruning pattern list” is described [66]. A pruning pattern is constructed froma solution, which is represented by a permutation of job numbers, by replacing some of its job
numbers by a “wild card” or a “don’t care” symbol. The job numbers to be replaced are determined by investigating the critical path of the schedule. A list of pruning patterns are generated
from “good” schedules that are obtained in the course of a search process, and maintained. The
list is used to inhibit the search to visit already searched and no longer interesting region again
and again.

12.2 Tabu Search
As described in Chapter 6, Tabu Search (TS) adopts a deterministic local search approach with a‘memory’ implemented as a ‘tabu list’ of moves which have been made in the recent past of thesearch, and which are ‘tabu’ or forbidden for a certain number of iterations. A tabu move may
be accepted (even if it is ’tabu’) if certain criteria are satisfied, such as the solution obtained by
the application of the move being better than the best solution obtained so far. Such criteria arecalled aspiration criteria.
Nowicki and Smutnicki [20] proposed a Tabu Search (TS) method for the Permutation Flowshop Scheduling Problem (PFSP). The tabu list used in their approach is described as follows byusing the notations used in Chapter 10.
Let π be a permutation of job numbers representing a schedule and v = (a, b) be a movethat removes the job π(a) from position a and re-inserts it into position b in π. When v = (a, b)with a < b is performed on π, the precedence relation between a pair of jobs (π(a), π(a + 1)) isreversed. To inhibit the future recovery, i.e., re-reversal, of the precedence relation, (π(a), π(a+1))
is stored as ‘tabu’ in the tabu list T of length tl. Then, a move v0 = (a0 , b0 ) with a0 < b0 cannot be103

104 Chapter 12. Tabu Search with a Pruning Pattern List for the Flowshop Scheduling Problem

v = (2,7)
π

x

y

z

tabu
v' = (1,6)
β

y

x

w

Figure 12.1: When v = (2, 7) is applied to π, (π(2), π(3)) = (x, y) is stored in T
as tabu. Later, v0 = (1, 6) is not allowed to apply to β because it will restore thepreviously banned precedence relation between (x, y).

performed on a new schedule β if v0 induces a recovery of any precedence relation in T , i.e., if{(β( j), β(a0 )) | j = a0 + 1, . . . , b0 } ∩ T , φ. Likewise, (π(a − 1), π(a)) is stored as ‘tabu’ if v = (a, b)with a > b is performed on π.
The neighborhood N(V \ W(π), π), which is a subset of N(V, π), is still too big for problemsof even moderate size. Thus, it is impracticable to evaluate all the neighbors in N(V \ W(π), π),so Nowicki and Smutnicki use a subset of N(V \ W(π), π) as a neighborhood structure. Alternatively, one can make the evaluation probabilistic: a member α from N(V \ W(π), π) is selected at
random and accepted if Cmax (α) < Cmax (π), otherwise probabilistically accepted according to theMetropolis probability Pc (α) = exp (−∆Cmax /c), where ∆Cmax = Cmax (α) − Cmax (π).

12.3 Pruning Pattern
A pruning pattern [π]l associated with π and its critical block Bl is derived from π by replacing
the jobs that belong to B̂l by ? as follows:(?
if j ∈ B̂l ,[π]l ( j) =π( j) otherwise.
For example, the pruning pattern corresponding to π = 4, 5, 6, 1, 2, 3, 8, 7 and B2 in Figure 10.1
is [π]2 = 4, 5, ?, ?, 2, 3, 8, 7. The makespan of [π]l is defined as the makespan of π, namely,
Cmax ([π]l ) := Cmax (π). The ? means a wild card, and a permutation β ‘matches’ [π]l if β( j) =[π]l ( j) at all but ? positions. For example 4, 5, 6, 1, 2, 3, 8, 7 and 4, 5, 1, 6, 2, 3, 8, 7 both match[π]2 . [π]l and a set of all permutations that match [π]l are identified. It is clear that N(Wl (π), π) ⊂[π]l . The block property described in Corollary 1 in Section 10.1 can be reformulated using [π]l12.4. Pruning Pattern List Approach105

as follows:
Corollary 2 For any schedule β and any pruning pattern [π]l , β ∈ [π]l =⇒ Cmax (β) ≥ Cmax ([π]l ).
This again can be reformulated as follows:
Corollary 3 For any two pruning patterns [π]i and [φ] j , [π]i ⊂ [φ] j =⇒ Cmax ([π]i ) ≥ Cmax ([φ] j ).
Corollary 2 suggests that when a new, possibly good, solution π is found during the search, [π]lidentifies a region where no solutions are better than π, and that excluding [π]l from the search
space can reduce the size of the search space without eliminating the global optima.

12.4 Pruning Pattern List Approach
The basic idea of the pruning pattern list approach is to reduce the size of the search space
effectively through storing the ‘important’ pruning patterns that correspond to a long critical
block of a good solution. Let [π] (without suffix) be the pruning pattern that corresponds to thelongest critical block of schedule π. The pruning pattern list PL with length pl is maintained
and updated as shown in Algorithm 12.4.1. Here, N(π) represents the neighborhood of π. N(V \W(π), π) or its subset is normally used as N(π). There is no good reason to keep patterns that are
rarely accessed in the list. Therefore such patterns are replaced by new patterns. According toCorollary 3, it is also not necessary to keep pruning patterns in the list that are subsets of otherpatterns. Thus they are removed.
Starting from an initial solution, the local search iteratively replaces the current solution with
one of its neighbors. In advanced local search strategies, such as Simulated Annealing (SA)and TS, cost-increasing neighbors can be accepted as well as cost-decreasing ones. Accepting
cost-increasing moves enables the search to escape from the local optima, but this may also
cause revisiting of previously evaluated points in the search space – something that is wastefulof computing resources in itself, and which also means the search is not adequately diversified.
One of the main aims of tabu search is to discourage such revisiting. This can be accomplishedby means of an explicit ‘tabu list’ of points previously visited, but normally it is easier, and more
efficient, to record specific attributes of such points, or of moves that would lead towards them.
Nevertheless, recording attributes of the points, and not the points themselves, can risk treatingeven moves that are better than any solution obtained so far as tabu. This is one of the main
reasons why aspiration criteria should be introduced into TS.
In the case of the pruning pattern list approach, Corollary 2 guarantees that a move thatmatches a pruning pattern is never better than the best solution obtained so far. This means thatan aspiration criterion is not necessary, and it also means that the pruning pattern list can serve asa longer-term memory to prevent the search getting stuck in an already searched and no longer
interesting region. Unlike the tabu list, an old pattern can remain in the pruning list as longas it is accessed frequently. If N(V \ W(π), π) or its subset is used as the neighborhood of thecurrent solution π, it does not contain any solution that matches [π]l , even without PL. However,
it is still possible that the region [π]l would be revisited in some later stages without any better106 Chapter 12. Tabu Search with a Pruning Pattern List for the Flowshop Scheduling Problem
Algorithm 12.4.1 The pruning pattern list approach to the PFSP.
A flowshop scheduling problem is given as an input. The size of PL is given as pl.
1. PL is initialized as a list of empty pl elements.
2. Generate a starting solution schedule π0 at random, and set π = π0 .
3. do(a) Calculate N(π) from π.
(b) N P (π) is initialized as N(π). For each α ∈ N P (π), if α matches any pattern [β] ∈ PL,then α is removed from N P (π), resulting N P (π) := N P (π) \ α, and the ‘access count’of [β] is incremented. N P (π) is used as the new neighborhood of π.
(c) doi. A candidate schedule φ is chosen from N P (π) and is accepted or rejected basedon the value Cmax (φ).
until φ is accepted.
(d) When the accepted φ is a cost-decreasing solution (i.e., Cmax (φ) < Cmax (π)), then its
longest pruning pattern [φ] is stored in PL.
(e) All the existing patterns [β] in PL such that [β] ⊂ [φ] are removed and substituted byempty patterns.
(f) If the number of non-empty patterns on PL exceeds pl, then the least accessed pattern[γ] is removed from the list.
(g) Set π = φ.
until termination conditions are satisfied.
4. Output the best solution obtained.

solutions than π being found. The size of PL is fixed to pl mainly for computational efficiency.
Theoretically, N P (π) may not satisfy the connectivity property that N(V \W(π), π) does. It should
be noted that the pruning pattern list is treated in a similar way that the population is treated in
Genetic Algorithms (GAs), especially in the steady state model [67]. The access count of thepruning pattern corresponds to the fitness function of the individual in GAs.

12.5 Experimental Results
The pruning pattern list approach described in the previous section can be embedded into any
local search method, including TS, SA, and even GAs [68, 69]. Here a simple probabilistic
version of the TS described in Section 12.2 is used as a test case. The programs are coded in C.

12.6. Concluding Remarks

107

The graph on the left side of Figure 12.2 shows the time evolution of the makespan averaged
over 30 runs of the TS with and without the pruning pattern list. This uses the ta041 problem,
which is one of Taillard’s benchmarks for size n = 50, m = 10 [21]. The parameters used arec = 3.0, tl = 7 and pl = 10. One run takes about 10 minutes on an HP 700 workstation. The
best and worst makespans obtained among 30 runs are 3000 and 3016 respectively when TS withthe pruning pattern list is used, and 3010 and 3025 when TS without the pruning pattern list isused. It can be seen that the pruning pattern list approach improves both the solution quality andthe speed. Without the pruning pattern list 12 runs out of 30 were trapped at a local minimumof makespan = 3025, which Nowicki and Smutnicki reported as the new reference makespanfor this problem. This suggests that their TS method can also be improved by incorporating thepruning pattern list approach.
The graph on the right side of Figure 12.2 shows the results of the ta051 problem (50 jobs
and 20 machines) based on 10 replications rather than 30. It also includes the computationallyequivalent MSXF-GA results reported in [69] for comparison. The pruning pattern list approachis also applied to other Taillard’s benchmarks of the same sizes, and the similar behaviors wereobserved for most of the cases. The results are summarized in Figure 12.3. The parameters used
are the same as in the ta041 and ta051 cases and the results are averaged over 10 replications.
The results of TS with the pruning pattern list generally outperform the results of TS without thepruning pattern list , except for the ta057 problem, where the difference is not so significant.

12.6 Concluding Remarks
The pruning pattern list approach to the makespan-minimizing permutation flowshop schedulingproblem is proposed and embedded into a Tabu Search method. The preliminary experimental results demonstrate the effectiveness of the proposed approach. Future research will aim to
implement an way to efficiently scan the pruning pattern list to find a match quickly, because
currently the scanning time becomes not negligible as the problem size increases. The proposed approach can also be embedded into GAs. However, the implementation details including
whether each individual should have its own pruning pattern list and should occasionally exchange it with others, or a single list should be shared by all the members in the population are
yet to be investigated. It is also expected to apply the pruning pattern list approach to the jobshop
scheduling problem.

108 Chapter 12. Tabu Search with a Pruning Pattern List for the Flowshop Scheduling Problem

ta041

makespanta051

makespan
3930

3052
3048

3920
3044
3910

3040
3036

3900
3032

TS without pruning pattern list
3028

TS without pruning pattern list

3890

3024

MSXF-GA
3880

3020

TS with pruning pattern list
TS with pruning pattern list

3016

3870
3012
3860

3008
0

100

200

300

400

500

600

700

iterations

0

100

200

300

400

500

600

700

iterations

Figure 12.2: The time evolutions of makespans for the ta041 (50 jobs and 10 machines) problem
averaged over 30 tabu search runs with and without the pruning pattern list (left). The time
evolutions for the ta051 (50 jobs and 20 machines) problem averaged over 10 tabu search runs
and the computationally equivalent MSXF-GA runs for comparison (right).

12.6. Concluding Remarks

Cmax

ta052

109

Cmax

ta053

3705

3785

3750

3695

3775

3740

3685

3765

3730

TS
TS+PL

3720
3710

Cmax

0

Cmax
3770
3760
3750
3740
3730
3720
3710
3700

3755

TS
TS+PL

3665

0

200 400 600
iterations

ta055

3735
0

Cmax

200 400 600
iterations

ta056

3770

3740

3760

3730

3750

3720
TS
3710
TS+PL
3700

3740

3690
0
200 400 600
iterations

Cmax

3825
3815
3805
3795
3785
TS
TS+PL 3775
3765
3755
200 400 600
0iterations

0

Cmax

3750

ta058

TS
TS+PL

3745

3655
0

3685
3675
3665
3655
3645
3635
3625
3615

3675

ta054

Cmax

3760

TS

3730

TS+PL

200 400 600
iterations

ta057

TS+PL
TS

3720
200 400 600 3710 0
iterations

ta059

Cmax

200 400 600
iterations

ta060

3815
3805
3795
TS

3785
TS+PL

3775

3765
0
200 400 600
iterations

TS
TS+PL

200 400 600
iterations

Figure 12.3: The time evolutions for the other nine Taillard problems of 50 jobs and 20 machines(ta052 – ta060) averaged over 10 tabu search runs with (labeled TS+PL) and without (labeledTS) the pruning pattern list.

Chapter 13
Conclusions
In this thesis, we have investigated various approaches to solve scheduling problems by metaheuristics, including Genetic Algorithms, Simulated Annealing and Tabu Search, and demonstrated the effectiveness of the proposed methods. We have basically started from a simple
problem-independent approach to more tailored problem-specific ones that involve more domain
specific knowledge of the scheduling problem.
The key features of the proposed methods include the effective use of the concept of active
schedule and the GT algorithm as well as the concept of critical path and blocks. The use of
these problem-specific knowledge greatly improves the performance of the proposed methods.
This is because the use of problem-specific knowledge enables to reduce the size of the search
space.
It appears to be always desirable to keep the size of the search space as small as possible.
For example, it is more efficient to search in the sub space of active schedules rather than thespace of all the semi-active schedules. The use of the critical block neighborhood reduces thesize of the neighorhood. TS also reduces the neighborhood size by temporally eliminating someof the members in the neighborhood as “tabu”. Likewise, the pruning pattern list introduced in
Chapter 12 also provides a mechanism to reduce the neighborhood size by eliminating alreadysearched and no longer interesting solutions from the neighborhood. In GAs, a redundant encoding of phenotype to genotype should be avoided because it increases the size of the search spacein which more than one genotypes correspond to the same phenotype. A time-consuming genedecoding process as well as allowing fatal genotypes to be generated and later to be repaired
must be avoided.
That said, a simple approach presented in Chapter 4 cannot be dismissed. Because it is
implemented easily and sometimes it is more robust, especially in the real-world situation theobjective function becomes more complicated and incorporating its domain specific knowledge
into the search structure may become more difficult.
Another key feature is the hybrid of GAs and other local search methods such as SA and
TS. It appears more promising to consider such hybrid rather than adhering to a single approach.
In fact, GAs are known to be unsuited for fine-tuning structures which are very close to optimal
solutions as opposed to SA and TS. On the other hand, SA and TS are inherently serial algorithms
and not straight-forward to be parallelized. Thus, the MSXF method proposed in Chapter 9 as a
110

111
new approach for Genetic Local Search would be promising. In this framework, each individual,or search agent in the population, performs local search using SA and/or TS as a main searchengine. Whereas crossover operates occasionally on the solutions of two selected individuals
in the population and produces a new solution, which is then used as an initial solution for thesubsequent local search. Here, the role of the crossover is to exchange informations between
the search agents which otherwise perform independent local search in parallel, rather than toperform search itself as in the conventional GAs. One promising future direction might be toconsider multiple parents in the MSXF.
Recently, the factory automation has been so advanced that each machine has self diagnostic sensors as well as a network interface and is monitored online even remotely. The flexibility
and the accuracy of controlling machines for reconfiguration and rescheduling are also improved.
Therefore, the need for generating an efficient and fine-tuned schedule in reasonable time is moreenvisaged. However unfortunately, the jobshop and flowshop scheduling problems investigatedin this thesis might be too simplistic compared to the real-world problems that have more complicated constraints, more flexible objective functions and more dynamic features. It is highly
expected to extend these approaches to incorporate more realistic settings. The author believes
that the ideas presented in this thesis, at least some of them, will be helpful for future research.

Bibliography[1] D. Applegate and W. Cook. A computational study of the job-shop scheduling problem.
ORSA Journal on Computing, Vol. 3, No. 2, pp. 149–156, 1991. (document), 1.1, 2.8, 2.6,
7.6, 8.3.2, 9.4, 9.4.2[2] B. Giffler and G.L. Thompson. Algorithms for solving production scheduling problems.
Operations Research, Vol. 8, pp. 487–503, 1960. 1.1, 2.2[3] C. Fischer and G.L. Thompson. Probabilistic learning combinations of local job-shop
scheduling rules. In Industrial Scheduling, pp. 225–251, Englewood Cliffs, N.J., 1963.
Prentice-Hall. 1.1[4] J.F. Muth and G.L. Thompson. Industrial Scheduling. Prentice-Hall, Englewood Cliffs,
N.J., 1963. 1.1, 2.8, 7.6, 8.3.1, 9.4[5] G.H. Brooks and C.R. White. An algorithm for finding optimal or near optimal solutions to
the production scheduling problem. The Journal of Industrial Engineering, Vol. 16, No. 1,
pp. 34–40, 1969. 1.1[6] E. Balas. Machine sequencing via disjunctive graphs: an implicit enumeration algorithm.
Operations Research, Vol. 17, pp. 941–957, 1969. 1.1, 6.3[7] J.R. Barker and G.B. McMahon. Scheduling the general job-shop. Management Science,Vol. 31, No. 5, pp. 594–598, 1985. 1.1[8] J. Carlier and E. Pinson. An algorithm for solving the job-shop problem. Management
Science, Vol. 35, No. 2, pp. 164–176, 1989. 1.1[9] P. Brucker, B. Jurisch, and B. Sievers. A branch & bound algorithm for the job-shop
scheduling problem. Discrete Applied Mathematics, Vol. 49, pp. 107–127, 1994. 1.1, 2.5,
6.3, 7.6[10] P. Martin and D.B. Shmoys. A new approach to computing optimal schedules for the jobshop scheduling problem. In the 5th International IPCO Conference, pp. 389–403, 1996.
1.1[11] J. Carlier and E. Pinson. Adjustment of heads and tails for the job-shop problem. European
Journal of Operational Research, Vol. 78, pp. 146–161, 1994. 1.1
112

Bibliography113[12] P.J.M. van Laarhoven, E.H.L. Aarts, and J.K. Lenstra. Job shop scheduling by simulated
annealing. Operations Research, Vol. 40, No. 1, pp. 113–125, 1992. 1.1, 2.3, 6.3, 7.5.2,
7.5.2, 7.6[13] H. Matsuo, C.J. Suh, and R.S. Sullivan. A controlled search simulated annealing method
for the general jobshop scheduling problem. Department of Management, The Universityof Texas at Austin, Vol. Working Paper, 03-04-88, , 1988. 1.1, 7.6, 7.7[14] J. Adams, E. Balas, and D. Zawack. The shifting bottleneck procedure for job shop scheduling. Management Science, Vol. 34, No. 3, pp. 391–401, 1988. 1.1, 2.3, 2.6, 2.8, 7.6, 8.2[15] F. Glover. Tabu search - part 1. ORSA Journal on Computing, Vol. 1, No. 3, pp. 190–206,
1989. 1.1
[16] F. Glover. Tabu search - part 2. ORSA Journal on Computing, Vol. 2, No. 1, pp. 4–32, 1990.
1.1[17] E.D. Taillard. Parallel taboo search techniques for the job-shop scheduling problem. ORSA
Journal on Computing, Vol. 6, No. 2, pp. 108–117, 1994. 1.1, 2.3, 6.3, 8.3.2[18] M. Dell’Amico and M. Trubian. Applying tabu search to the job-shop scheduling problem.
Annals of Operations Research, Vol. 41, pp. 231–252, 1993. 1.1, 2.5, 9.4.2[19] E. Nowicki and C. Smutnicki. A fast tabu search algorithm for the job-shop problem.
Management Science, Vol. 42, pp. 797–813, 1996. 1.1, 2.5, 9.4.2[20] E. Nowicki and C. Smutnicki. A fast tabu search algorithm for the permutation flow-shop
problem. European Jouranl of Operational Research, Vol. 91, pp. 160–175, 1996. 1.1,
10.2, 10.6, 12.2
[21] E. Taillard. Benchmarks for basic scheduling problems. European Journal of Operational
Research, Vol. 64, pp. 278–285, 1993. 1.1, 10.4, 11.4, 12.5
[22] J.H. Holland. Adaptation in Natural and Artificial Systems. Unuv. of Michigan Press, 1975.
1.1[23] D.E. Goldberg. Genetic Algorithms in Search, Optimization, and Machine Learning.
Addison-Wesley, Reading,Mass., 1986. 1.1[24] R. Nakano and T. Yamada. Conventional genetic algorithm for job shop problems. In
Proceedings of International Conference on Genetic Algorithms (ICGA ’91), pp. 474–479,
1991. 1.1, 2.4, 2.4, 4, 4.2[25] T. Yamada and R. Nakano. A genetic algorithm applicable to large-scale job-shop problems.
In Proceedings of The Second International Conference on Parallel Problem Solving fromNature PPSN ’92, pp. 281–290, 1992. 1.1, 5, 5.3

114

Bibliography[26] N.L.J. Ulder, E. Pesch, P.J.M. van Laarhoven, J. Bandelt, H, and E.H.L. Aarts. Genetic
local search algorithm for the traveling salesman problem. In Parallel Problem Solvingfrom Nature, 1, pp. 109–116, 1994. 1.1, 9[27] S. S. Panwalkar and Wafix Iskander. A survey of scheduling rules. Operations Research,
Vol. 25, No. 1, pp. 45–61, 1977. 2.2[28] U. Dorndorf and E. Pesch. Evolution based learning in a job shop scheduling environment.
Computers Ops Res, Vol. 22, pp. 25–40, 1995. 2.2[29] B. Roy and B. Sussmann. Les probl‘emes d’ordonnancement avec contraintes disjonctives.
Note DS no 9 bis, SEMA, Paris, 1964. 2.3[30] P. Brucker. Scheduling Algorithms. Springer-Verlag, Berlin, 1995. 2.5[31] J. Carlier. The one-machine sequencing problem. European Journal of Operational Research, Vol. 11, pp. 42–47, 1982. 2.7, 2.7[32] S. Lawrence. Resource constrained project scheduling: an experimental investigation ofheuristic scheduling techniques (supplement). Technical report, Graduate School of Industrial Administration, Carnegie Mellon University, 1984. 2.8, 7.6[33] S. Kobayashi, I. Ono, and M. Yamamura. An efficient genetic algorithm for job shop
scheduling problems. In Proceedings of 6th International Conference on Genetic Algorithms and their Applications (Pittsburgh, PA), pp. 506–511, 1995. 5.1[34] I. Yamamura, M. Ono and S. Kobayashi. Character-preserving genetic algorithms for traveling salesman problem (in japanese). Journal of Japanese Society for Artificial Intelligence, Vol. 7, pp. 1049–1059, 1992. 5.1[35] C. Bierwirth. A generalized permutation approach to job shop scheduling with genetic
algorithms. OR Spektrum, Vol. 17, pp. 87–92, 1995. 5.2[36] C. Bierwirth, D. Mattfeld, and H. Kopfer. On permutation representations for schedulingproblems. In Parallel Problem Solving from Nature, 4, pp. 310–318, 1996. 5.2[37] Y. Davidor, T. Yamada, and R. Nakano. The ecological framework II: Improving ga performance at virtually zero cost. In Proceedings of 5th International Conference on Genetic
Algorithms and their Applications (Urbana - Champaign), pp. 171–176, 1993. 5.3[38] M.R. Garey and D.S. Johnson. Computers and Intractability - A Guide to the Theory of
NP-Completeness. Freeman and Company, New York, 1979. 6[39] C.R. (ed.) Reeves. Modern Heuristic Techniques for Combinatorial Problems. Blackwell
Scientific Publications, Oxford, UK; re-issued by McGraw-Hill, London, UK., 1995. 6

Bibliography115[40] C.R. Reeves. Heuristic search methods: A review. In D.Johnson and F.O’Brien Operational
Research: Keynote Papers, pp. 122–149. Operational Research Society, Birmingham, UK,1996. 6[41] D.S. Johnson. Local optimization and the traveling salesman problem. In G.Goos and
J.Hartmanis (Eds.) Automata, Languages and Programming, Lecture Notes in Computer
Science 443, pp. 446–461. Springer-Verlag, Berlin, 1990. 6.2[42] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller. Equationof state calculations by fast computing machines. Journal of Chemical Physics, Vol. 21,
No. 6, pp. 1087–1092, 1953. 7.1[43] P. D. Wasserman. Neural Computing: Theory and Practice. Van Nostrand Reinhold, New
York, 1989. 7.1[44] S. Kirkpatrick, C.D.Jr. Gelatt, and M.P. Vecchi. Optimization by simulated annealing.
Science, Vol. 220, pp. 671–680, 1983. 7.1[45] L. Ingber and B E.Rosen. Genetic algorithms and very fast simulated reannealing: A comparison. Mathematical and Computer Modeling, Vol. 16, No. 11, pp. 87–100, 1992. 7.3[46] B. E. Rosen. Function optimization based on advanced simulated annealing. Workshop on
Physics and Computation, PhysComp 92, 1992. 7.3[47] E. H. L. Aarts and J. H. M. Korst. Simulated Annealing and Boltzmann machines. Wiley,
Chichester, 1989. 7.4, 7.7[48] D. Applegate. Jobshop benchmark problem set. Personal Communication, 1992. 3[49] J. Carlier. Ordonnancements à contraintes disjonctives. RAIRO, Vol. 12, pp. 333–351,
1978. 7.6[50] T. Yamada, B.E. Rosen, and R. Nakano. A simulated annealing approach to job shop
scheduling using critical block transition operators. In Proceedings of IEEE International
Conference on Neural Networks (Orlando, Florida.), pp. 4687–4692, 1994. 8.3.1[51] E.H.L. Aarts, P.J.M. van Laarhoven, J.K. Lenstra, and N.L.J. Ulder. A computational studyof local search algorithms for job shop scheduling. ORSA Journal on Computing, Vol. 6,
No. 2, pp. 118–125, 1994. 8.3.2, 9.4.2[52] H. Kopfer D.C. Mattfeld and C. Bierwirth. Control of parallel population dynamics bysocial-like behavior of ga-individuals. In Parallel Problem Solving from Nature, 3, 1994.
8.3.2, 9.4.2[53] R.J.M. Vaessens, E.H.L. Aarts, and J.K. Lenstra. Job shop scheduling by local search.
Technical report, Eindhoven University of Technology, Dpt. of Math. and CS, February
1994. 8.3.2

116

Bibliography[54] C. R. Reeves. Genetic algorithms and neighbourhood search. In Evolutionary Computing,
AISB Workshop (Leeds, U.K.), pp. 115–130, 1994. 9.1[55] T. Yamada and R. Nakano. A genetic algorithm with multi-step crossover for job-shop
scheduling problems. In Proceedings of the 1st IEE/IEEE International Conference on
Genetic ALgorithms in Engineering Systems (GALESIA ’95), pp. 146–151, 1995. 9.1[56] T. Yamada and R. Nakano. Scheduling by genetic local search with multi-step crossover.
In Proceedings of The Fourth International Conference on Parallel Problem Solving fromNature (PPSN ’96), pp. 960–969, 1996. 9.1[57] T. Yamada and R. Nakano. A fusion of crossover and local search. In Proceedings of IEEE
International Conference on Industrial Technology, pp. 426–430, 1996. 9.1[58] D. Whitley. The genitor algorithm and selection pressure: why rank-based allocation of
reproductive trials is best. In Proceedings of 3rd International Conference on Genetic
Algorithms and their Applications (Arlington,VA), pp. 116–121, 1989. 9.3[59] G. Syswerda. Uniform crossover in genetic algorithms. In Proceedings of 3rd International
Conference on Genetic Algorithms and their Applications (Arlington,VA), pp. 2–9, 1989.
9.3[60] J. E. Beasley. Or-library: distributing test problems by electronic mail. European Journalof Operational Research, Vol. 41, pp. 1069–1072, 1990. 9.4.2, 10.6
[61] C. R. Reeves and C. Höhn. Are long path problems hard for genetic algorithms? In 4thInternational Conference on Parallel Problem Solving from Nature, pp. 134–153, 1996.
10.4[62] K. D. Boese, A. B. Kahng, and S. Muddu. A new adaptive multi-start technique for combinatorial global optimization. Operations Research Letters, Vol. 16, pp. 101–113, 1994.
10.4[63] C. R. Reeves. Landscapes, operators and heuristic search. Annals of Operations Research,
Vol. to appear, , 1998. 10.4
[64] T. Jones and S. Forrest. Fitness distance correlation as a measure of problem difficulty
for gas. In Proceedings of 6th International Conference on Genetic Algorithms and their
Applications (Pittsburgh, PA), pp. 184–192, 1995. 10.4
[65] J. Liu. A new heuristic algorithm for csum flowshop scheduling problems. Personal Communication, 1997. 11.4
[66] T. Yamada. A pruning pattern list approach to the permutation flowshop scheduling problem. Kluwer academic publishers, MA, USA, 2002. 12.1

Bibliography

117[67] D. Whitley, T. Starkweather, and D. Fuquay. Scheduling problems and traveling salesman:
The genetic edge recombination operator. In Proceedings of 3rd International Conference
on Genetic Algorithms and their Applications (Arlington,VA), pp. 133–140, 1989. 12.4
[68] T. Yamada and C.R. Reeves. Permutation flowshop scheduling by genetic local search.
In Proceedings of the 2nd IEE/IEEE International Conference on Genetic ALgorithms in
Engineering Systems (GALESIA ’97), pp. 232–238, 1997. 12.5
[69] C.R. Reeves and T. Yamada. Genetic algorithms, path relinking and the flowshop sequencing problem. Evolutionary Computation journal (MIT press), Vol. 6, No. 1, pp. 45–60,
1998. 12.5

A List of Author’s Work
Journals
1. T.Yamada, B.E.Rosen and R.Nakano, “Critical Block Simulated Annealing for Job Shop
Scheduling (in Japanese),” The Transaction of The Institute of Electrical Engineers ofJapan, Vol.114-C, No.4, pp.476–482 (1994).
2. T.Yamada and R.Nakano, “Job-Shop Scheduling by Simulated Annealing Combined with
Deterministic Local Search (in Japanese),” Transactions of Information Processing Societyof Japan, Vol.37 No.4, pp. 597–604 (1996).
3. T.Yamada and R.Nakano, “Job-Shop Scheduling by Genetic Local Search (in Japanese),”Transactions of Information Processing Society of Japan, Vol.38 No.6, pp. 1126–1138(1997).
4. T.Yamada and C.R.Reeves, “Landscape Analysis of the Flowshop Scheduling Problem
and Genetic Local Search (in Japanese),” Transactions of Information Processing Societyof Japan, Vol.39 No.7, pp. 2112–2123 (1998).
5. C.R.Reeves and T.Yamada, “Genetic Algorithms, Path Relinking and the Flowshop Sequencing Problem,” Evolutionary Computation journal, MIT press, Vol.6 No.1, pp.45–60,
Spring (1998).

Books
1. T.Yamada and R.Nakano, “Job-Shop Scheduling by Simulated Annealing Combined with
Deterministic Local Search,” in Meta-Heuristics: Theory & Applications, Kluwer academic publishers, pp. 237–248 (1996).
2. T.Yamada and R.Nakano, “Chapter 7: Job Shop Scheduling,” in Genetic algorithms in
Engineering Systems, The Institution of Electrical Engineers, pp.134–160 (1997).
3. C.R.Reeves and T.Yamada, “Goal-Oriented Path Tracing Methods,” in New Ideas in Optimization, The MacGraw-Hill Companies, pp.341–355 (1999).
118

A List of Author’s Work

119

4. T.Yamada “A Pruning Pattern List Approach to the Permutation Flowshop Scheduling
Problem,” in Essays and Surveys in Metaheuristics, Kluwer academic publishers, pp. 641–
651 (2002).

International Conferences
1. R.Nakano and T.Yamada, “Conventional Genetic Algorithm for Job Shop Problems,” Proceedings of International Conference on Genetic Algorithms (ICGA ’91), pp.474–479(1991).
2. T.Yamada and R.Nakano, “A Genetic Algorithm Applicable to Large-Scale Job-Shop Problems,” Proceedings of The Second International Conference on Parallel Problem Solvingfrom Nature (PPSN ’92), pp.281–290 (1992).
3. Y.Davidor, T.Yamada and R.Nakano, “The ECOlogical Framework II: Improving GA Performance At Virtually Zero Cost,” Proceeding of International Conference on Genetic
Algorithms (ICGA ’93), pp.171–176 (1993).
4. T.Yamada, B.E.Rosen and R.Nakano, “A Simulated Annealing Approach to Job Shop
Scheduling using Critical Block Transition Operators,” Proceedings of IEEE International
Conference on Neural Networks (ICNN ’94), pp.4687–4692 (1994).
5. T.Yamada and R.Nakano, “A Genetic Algorithm with Multi-Step Crossover for Job-Shop
Scheduling Problems,” Proceedings of the 1st IEE/IEEE International Conference on Genetic ALgorithms in Engineering Systems (GALESIA ’95), pp.146–151 (1995).
6. T.Yamada and R.Nakano, “Scheduling by Genetic Local Search with Multi-Step Crossover,”
Proceedings of The Fourth International Conference on Parallel Problem Solving from Nature (PPSN ’96), pp.960–969 (1996).
7. T.Yamada and R.Nakano, “A Fusion of Crossover and Local Search,” Proceedings of IEEE
International Conference on Industrial Technology, pp.426–430 (1996).
8. T.Yamada and C.R.Reeves, “Permutation Flowshop Scheduling by Genetic Local Search,”
Proceedings of the 2nd IEE/IEEE International Conference on Genetic ALgorithms in Engineering Systems (GALESIA ’97), pp. 232–238, 1997.
9. T.Yamada and C.R.Reeves, “Solving the Csum Permutation Flowshop Scheduling Problem by Genetic Local Search”, Proceedings of 1998 IEEE International Conference on
Evolutionary Computation, pp.230–234 (1998).
10. C.R.Reeves and T.Yamada, “Implicit Tabu Search Methods for Flowshop Sequencing,”
Proceedings of IMACS International Conference on Computational Engineering in Systems Applications, pp.78–81 (1998).

120

A List of Author’s Work

11. T. Yamada, K. Yoshimura and R. Nakano, “Information Operator Scheduling by Genetic
Algorithms,” Simulated Evolution and Learning: Proceedings of the Second Asia-PacificConference on Simulated Evolution and Learning (SEAL’98), Lecture Notes in Computer
Science 1585, pp. 50–57 (1999).