# A PMBGA to Optimize the Selection of Rules for Job Shop Scheduling Based on the Giffler-Thompson Algorithm

Authors: Rui Zhang; Cheng Wu
Journal: Journal of Applied Mathematics
 (2012)

---

## Abstract

Most existing research on the job shop scheduling problem has been focused on the minimization of makespan (i.e., the completion time of the last job). However, in the fiercely competitive market nowadays, delivery punctuality is more important for maintaining a high service reputation. So in this paper, we aim at solving job shop scheduling problems with the total weighted tardiness objective. Several dispatching rules are adopted in the Giffler-Thompson algorithm for constructing active schedules. It is noticeable that the rule selections for scheduling consecutive operations are not mutually independent but actually interrelated. Under such circumstances, a probabilistic model-building genetic algorithm (PMBGA) is proposed to optimize the sequence of selected rules. First, we use Bayesian networks to model the distribution characteristics of high-quality solutions in the population. Then, the new generation of individuals is produced by sampling the established Bayesian network. Finally, some elitist individuals are further improved by a special local search module based on parameter perturbation. The superiority of the proposed approach is verified by extensive computational experiments and comparisons.


---

## Body

## 1. Introduction

The job shop scheduling problem (JSSP) has been known as an extremely difficult combinatorial optimization problem ever since the 1950s. In terms of computational complexity, JSSP is strongly𝒩𝒫-hard [1]. Because of its relevance to contemporary manufacturing systems, extensive research has been conducted on the problem [2]. In recent years, the metaheuristics—such as simulated annealing (SA) [3], genetic algorithm (GA) [4–6], tabu search (TS) [7, 8], particle swarm optimization (PSO) [9, 10], and ant colony optimization (ACO) [11, 12]—have clearly become the most popular methods.However, most existing research is still based on the standard JSSP model, so it is inconvenient to apply these algorithms directly to real-world scheduling scenarios. Especially, most algorithms are designed for the makespan criterion (i.e., minimizeCmax). Actually, in the make-to-order production environment nowadays, due date related performances are becoming increasingly significant, because the in-time delivery of goods is vital for maintaining a high service reputation. Therefore, the algorithms that aim at minimizing tardiness in JSSP deserve more investigations.In an attempt to contribute to the scheduling community, in this paper we study the JSSP with the objective of minimizing total weighted tardiness (TWT). TWT is the weighted sum of each job’s tardiness against its due date. In some respects, the TWT measure captures the critical factors related with the profitability of a firm more accurately than the makespan. Meanwhile, from the theoretical perspective, the complexity of solving JSSP with TWT objective (abbreviated as TWT-JSSP hereinafter) is much greater than that of solving JSSP with makespan objective [13]. The value of TWT is affected by all the tardy jobs and thus is more sensitive to the changes in the schedule. According to the three-field notation, the studied problem can be described as Jm||∑wjTj.Relatively few publications have discussed the minimization of TWT in JSSP. The only exact method is the branch-and-bound algorithm proposed by Singer and Pinedo [14], while the rest belong to the heuristic category. The first attempt of adapting the shifting bottleneck algorithm (which was very successful for makespan minimization) to TWT-JSSP is reported in [15]. In [16, 17], the authors propose modified shifting bottleneck heuristics for complex job shops (characterized by parallel batching machines, sequence-dependent setup times, and reentrant process flows) with TWT criterion. In these algorithms, the single/parallel machine subproblems are solved basically by the ATC rule or the extended BATCS rule (for the case of batching and setups). The shifting bottleneck heuristic proposed in [18], however, uses genetic algorithms to solve the subproblems. The computations show that using GA as subproblem solution procedures leads to improved results compared to dispatching-based subproblem solution procedures. Besides shifting bottleneck, a large step random walk (LSRW) algorithm is designed in [19], which uses different neighborhood sizes depending on whether a small step or a large step is performed. A small step consists of iterative improvement, while a large step consists of a Metropolis algorithm. In [20, 21], hybrid genetic algorithms are presented for solving TWT-JSSP. The former combines GA with different dispatching rules, while the latter combines GA with an iterated local search procedure that uses the longest path approach on a disjunctive graph model. In [22], a tabu search algorithm is presented for the generalized TWT-JSSP with release times and precedence constraints. Recently, electromagnetic algorithms (EM) have been successfully applied to TWT-JSSP [23]. In fact, most of these algorithms have relatively high computational complexity, and thus they are incapable of solving large-scale TWT-JSSP.Rule-based scheduling approaches have been investigated for job shop scheduling since the late 1990s [24, 25]. A remarkable advantage of using dispatching rules is that it helps to save computational time and generally produces satisfactory solutions. Recently, modern optimization strategies are applied to find high-quality combinations of dispatching rules. For example, a neural network is designed in [26] to determine the dispatching rule to use on each machine in a job shop (it is necessary to train the neural network beforehand with optimal rule combinations for some known instances). A memetic algorithm with tailored-encoding/decoding schemes, and a local search procedure is proposed in [27] for minimizing the number of tardy jobs in JSSP. The memetic algorithm defeats a multistart hill-climbing approach and a simulated annealing approach. A dispatching rule-based genetic algorithms (DRGA) is proposed in [28], which searches simultaneously for the best sequence of dispatching rules and the number of operations to be handled by each dispatching rule. The DRGA obtains better results than a GA using the conventional dispatching rule representation and a GA that uses the operation permutation representation. A genetic programming-based data mining approach is proposed in [29] to select dispatching rules under a given set of shop parameters (e.g., interarrival times). The results obtained from simulation show that the selected dispatching rules are appropriate according to the current shop status.The rest of the paper is organized as follows. The discussed problem is mathematically formulated in Section2. Section 3 makes a brief introduction to the principles of PMBGA. Section 4 proposes a rule-based PMBGA for solving TWT-JSSP. Section 5 presents the computational results and analysis. Finally, Section 6 concludes the paper.


## 2. Problem Formulation

In a JSSP instance, a set ofn jobs {Jj}j=1n are to be processed on a set of m machines {Mk}k=1m under the following basic assumptions: (i) there is no machine breakdown; (ii) no preemption of operations is allowed; (iii) all jobs are released at time 0; (iv) the transportation times and the setup times are all neglected; (v) each machine can process only one job at a time; (vi) each job may be processed by only one machine at a time.Each job has a fixed processing route which traverses the relevant machines (in the standard JSSP benchmark instances, each job is required to visit all the machines. But actually, the number of operations for each job (mj) can be less than m) in a predetermined order. The manufacturing process of job j on machine k is noted as operation Ojk, with a duration of pjk. Besides, a preset due date dj (describing the level of urgency) and a preset weight wj (reflecting the importance of the order) are given for each job j. The objective function is defined as TWT=∑j=1nwjTj, where Tj=max{Cj-dj,0} defines the tardiness of job j.JSSP can be described by a disjunctive graphG(O,A,E) [30]. O={Ojk∣j=1,…,n,k=1,…,m} is the set of nodes. A is the set of conjunctive arcs which connect successive operations of the same job, so A describes the technological constraints in the JSSP instance. E=⋃k=1mEk is the set of disjunctive arcs, where Ek denotes the disjunctive arcs corresponding to the operations on machine k. Each arc in Ek connects a pair of operations to be processed by machine k and ensures that the two operations should not be processed simultaneously. Initially, the disjunctive arcs do not have fixed directions.Under the disjunctive graph representation, finding a feasible schedule for the JSSP is equivalent to orienting all the disjunctive arcs so that no directed cycles exist in the resulting graph. In this paper, we useσ to denote the set of directed disjunctive arcs which are transformed from the original E. Thus, if A∪σ is acyclic, the schedule corresponding to σ is feasible (in the rest of the paper, we do not distinguish between σ and the schedule. For the convenience of expression, we will write σ as a matrix. The kth row of σ represents the processing order of the operations on machine k).Based on the disjunctive graph model, the discussed TWT-JSSP can be mathematically formulated as follows:(2.1)min⁡TWT=∑j=1nwj(tjkj+pjkj-dj)+,(a)s.t.tjk+pjk⩽tjk′,∀(Ojk,Ojk′)∈A,(b)tjk+pjk⩽tj′k∨tj′k+pj′k⩽tjk,∀(Ojk,Oj′k)∈Ek,(c)tjk⩾0,j=1,…,n,k=1,…,m.In this formulation, (x)+=max{x,0}. tjk represents the starting time of operation Ojk. kj denotes the index of the machine that processes the last operation of job j, so the completion time of job j is tjkj+pjkj. The set of constraints (a) ensure that the processing order of the operations from each job is consistent with the technological routes. The set of constraints (b) ensure that any two operations on the same machine cannot be processed simultaneously.


## 3. A Brief Introduction to PMBGA

Recently, there has been a growing interest in the evolutionary algorithms that explore the search space by building and utilizing probabilistic models of high-quality solutions. Indeed, these algorithms use the following two steps to replace the conventional crossover and mutation operators in GA:(1)
Build a probabilistic model of the selected promising solutions;(2)
Sample the built model to produce a new generation of candidate solutions.The evolutionary algorithms based on such a principle are referred to as estimation of distribution algorithms (EDAs) or probabilistic model-building genetic algorithms (PMBGAs). The major steps of a PMBGA implementation are listed as follows, whereGN is the maximum number of generations.Step 1.
Set the generation indexg=0. Initialize the population of the first generation, that is, P(0).Step 2.
Select a subsetS of promising individuals from P(g).Step 3.
Establish a probabilistic modelM which somehow describes the distribution characteristics of S.Step 4.
Generate a setN of new individuals by sampling M.Step 5.
Select the best|P(g)| individuals from P(g)∪N and assign them to the next generation population P(g+1).Step 6.
Letg←g+1. If g<GN, return to Step 2. Otherwise, output the best solution in P(g).The PMBGA is especially useful for the complex optimization problems in which the decision variables are correlated. For such problems, the realized value of a certain decision variable can produce an impact on the optimal value for another decision variable. Therefore, if these variables are optimized in a separate way (or one by one), traditional GA will be very likely to converge to local optimum. PMBGA improves traditional GA by modeling the relationship between decision variables. Clearly, the most crucial element in the design of PMBGA is the type of the adopted probabilistic model (in Step3), which directly affects the algorithm’s capability of producing high-quality offspring solutions. In the artificial intelligence community, the commonly adopted graphical model for characterizing the relationship between a set of discrete variables is the Bayesian network [31]. A Bayesian network is a directed acyclic graph. A node in the Bayesian network indicates a variable under investigation (each variable actually corresponds to a coding gene for the solutions in PMBGA), and an arc indicates the probabilistic causal relationship between the two nodes connected by it. The direction of the arc implies that the variable corresponding to the head node of the arc is conditioned by the variable corresponding to the tail node. In general, the joint probabilistic distribution of an n-variate random vector X=(X1,…,Xn) described by a Bayesian network can be calculated as(3.1)P(x)=∏i=1nP(xi∣pa(xi)).
In this formulation, x=(x1,…,xn) is a vector of realized values for X; pa(xi) is a set of realized values for the parents (in a Bayesian network, if there exists a directed arc pointing from node Xj to Xi, then Xj is called a parent of Xi) of the random variable Xi.A detailed introduction to the PMBGA can be found in [32]. Some important advances and interesting applications of EDA are covered in [33, 34]. Successes in utilizing EDA to solve scheduling problems have been reported in [35, 36] (for flow shop scheduling).


## 4. The Proposed PMBGA for Solving TWT-JSSP

### 4.1. Encoding

The proposed PMBGA relies on dispatching rules to record the scheduling policies. Eight scheduling rules are involved in this study. In the following expressions for the priority index, operationi belongs to job j, and wj and dj are the corresponding job’s weight and due date. JS(i) and JP(i), respectively, denote the set of job successors of operation i and the set of job predecessors of operation i. Z is used to indicate that the operation with the smallest index will be chosen from the conflict set, while Z′ is used to indicate that the operation with the largest index will be chosen from the conflict set:(1)
ATC (apparent tardiness cost):Zi′(t)=(wj/pi)×exp{-(dj-ri-pi - ∑i′∈JS(i)(Ŵi′+pi′))+/(K×p¯)} (in this expression, ri denotes the earliest starting time of operation i, and p¯ denotes the average processing time of the current operations in the conflict set. K is a scaling parameter (or called “look-ahead” parameter). Ŵi is the estimated lead time of operation i);(2)
SPT (shortest processing time):Zi=pi;(3)
LPT (longest processing time):Zi′=pi;(4)
WSPT (weighted shortest processing time):Zi′=wj/pi;(5)
SRPT (shortest remaining processing time):Zi=pi+∑i′∈JS(i)pi′;(6)
LRPT (longest remaining processing time):Zi′=pi+∑i′∈JS(i)pi′;(7)
EDD (earliest due date):Zi=dj;(8)
ODD (operation due date):Zi(t)=[dj/∑i′∈JS(i)∪{i}pi′]×∑i′∈JP(i)∪{i}pi′.In PMBGA, each encoding digit (i.e., gene) is expressed by the serial number (1~8) of the selected dispatching rule. A solution is represented by a rule sequence {Rik:i=1,…,n,k=1,…,m}, where Rik indicates the rule to use when scheduling the ith operation on machine k. Therefore, the encoding length for each solution is l=n×m.


### 4.2. Decoding

In order to evaluate the fitness of a solution, the Giffler-Thompson algorithm is applied to construct an active schedule based on the specified dispatching rules. The implementation of the Giffler-Thompson algorithm is detailed below.Input:  A sequence of rules {Rik:i=1,…,n,k=1,…,m}.Step 1.
LetQ(1)=O={1,…,nm} (the set of all operations), R(1)=F(O)={f1,…,fn} (the set of first operations of each job). Set t=1 and πk=1 (k=1,…,m).Step 2.
Find the operationi*=argmini∈R(t){ri+pi}, and let m* be the index of the machine on which this operation should be processed. Use B(t) to denote all the operations from R(t) which should be processed on machine m*.Step 3.
Delete fromB(t) the operations that satisfy ri⩾ri*+pi*.Step 4.
Use ruleRπm*m* to identify an operation ô from B(t) if currently there are more than one candidates. Schedule operation ô on machine m* at the earliest possible time. Let πm*←πm*+1.Step 5.
LetQ(t+1)=Q(t)∖{ô}, R(t+1)=R(t)∖{ô}∪{suc(ô)}, where suc(ô) denotes the immediate job successor of operation ô (if any).Step 6.
IfQ(t+1)≠∅, set t←t+1 and go to Step 2. Otherwise, the decoding procedure is terminated.In the above description, the release timeri equals the earliest possible starting time of operation i (determined from the already scheduled operations). So, (ri+pi) is the earliest possible completion time of operation i. Q(t) represents the set of operations yet to be scheduled at iteration t, while R(t) represents the set of ready operations (whose job predecessors have all been scheduled) at iteration t. In Step 4, the operation set B(t) is also called a conflict set.


### 4.3. Producing Offspring Individuals

(a) Selection
In each iteration, we first sort all the individuals in the current population according to their fitness. Then, we select the best1/4 of individuals to form the set S, which will subsequently be used to build the Bayesian network in order to produce a new generation of individuals.(b) The Adopted Network Structure
A Bayesian network has to be built to describe the probabilistic distribution of favorableRik settings based on the elite individuals in S. Each individual can be characterized by a directed acyclic network as shown in Figure 1. In this network, a node Ni,kr (r∈{1,2,…,8}) indicates the fact that the dispatching rule used to schedule the ith operation on machine k is selected as rule Number r. The directed arc from node Ni,kr to node Ni,k+1r′ (or Ni+1,1r′) represents the dependency between the two nodes, so it characterizes the possible influence of the rule selection for the ith operation on machine k on the rule selection for the ith operation on machine k+1 (or the (i+1)-th operation on machine 1). Therefore, a directed path from a certain node in the first row to a certain node in the (n×m)-th row can completely describe an individual in the population (because a directed path records an instantiation of all the rule selections).Figure 1
The structure of the Bayesian network used in PMBGA.(c) Calculation of the Probability Values
Since we adopt a fixed network structure in PMBGA, building the Bayesian network is equivalent to determining the values of all the conditional probabilities according to the selected solution setS. After that, new individuals will be produced by iteratively sampling these probabilistic distributions, expecting to obtain high-quality offsprings.Given a number of individuals (i.e., the training setS), an estimate of the conditional probabilities can be obtained simply by counting the frequencies of occurrence.Example 4.1.
Here, we provide a concrete example to illustrate the probability calculation process. For a PMBGA optimization instance withn=3,m=1 and r=3, let us further suppose the current S contains 40 individuals. In Figure 2, the statistics of these individuals are displayed on a network as previously defined. The weight of each arc (the number placed on the arc) indicates the occurring frequency of this arc (counted for all the individuals in S). For example, if there exists an individual coded as “[3,1,2]”, then the path “N1,13→N2,11→N3,12” in the Bayesian network is used to record this individual, and consequently, the weights (counted frequencies) of the arcs N1,13→N2,11 and N2,11→N3,12 will be increased by 1, respectively. Note that, in the final network, the sum of the weights of all the incoming arcs of a certain node should be equal to the sum of the weights of all the outgoing arcs of the same node. This is because each individual corresponds to a complete path from the first row to the last row.Figure 2
An example network in PMBGA.By using frequency as an approximation for probability, the relevant (conditional) probabilities should be calculated as follows:(4.1)P(N1,11)=8+240,P(N1,12)=7+3+440,P(N1,13)=5+5+640,P(N2,11∣N1,11)=88+2,P(N2,12∣N1,11)=0,P(N2,13∣N1,11)=28+2,P(N2,11∣N1,12)=77+3+4,P(N2,12∣N1,12)=37+3+4,P(N2,13∣N1,12)=47+3+4,P(N2,11∣N1,13)=55+5+6,P(N2,12∣N1,13)=55+5+6,P(N2,13∣N1,13)=65+5+6,P(N3,11∣N2,11)=22+3+15,P(N3,12∣N2,11)=32+3+15,P(N3,13∣N2,11)=152+3+15,P(N3,11∣N2,12)=22+6,P(N3,12∣N2,12)=62+6,P(N3,13∣N2,12)=0,P(N3,11∣N2,13)=0,P(N3,12∣N2,13)=77+5,P(N3,13∣N2,13)=57+5.According to the above calculation method, a connection can never be rediscovered in the PMBGA if the corresponding conditional probability is zero (e.g., fromN2,12 to N3,13). To overcome this drawback, we can set the minimum count to 1. Taking N2,12 as an example, the conditional probabilities for the outgoing arcs will then become(4.2)P(N3,11∣N2,12)=33+7+1,P(N3,12∣N2,12)=73+7+1,P(N3,13∣N2,12)=13+7+1.
Now it is possible to discover N3,13 from N2,12, though the probability is small.(d) The Sampling Process
The sampling process for generating a new individual begins from the root nodes of the Bayesian network. By selecting an outgoing arc at each node based on the calculated conditional probabilities, the whole network can be gradually instantiated.


### 4.4. The Embedded Local Search

The search mechanism of GA guarantees a good performance in the “exploration” of the solution space. However, it has been reported that GA alone cannot achieve satisfactory solution quality for complex optimization problems. Actually, a local search procedure is usually added within the framework of GA in order to provide reliable “exploitation” ability. In this paper, we design such a local optimizer, which attempts to improve the selected solutions. In each iteration of PMBGA, the local search is carried out for the beste% of solutions in the current population. Thus, e is an important parameter for adjusting the frequency of local search and achieving a balance between exploration and exploitation.In the following, we will describe how to perform the local search on a given solutions (a sequence of rules).Step 1.
Use the Giffler-Thompson algorithm to generate an active scheduleσ1 based on s. The objective value of the obtained schedule is TWTσ1.Step 2.
Setu=1.Step 3.
Exert random perturbations on the processing times: generate a new set of processing times{pjk(u)} from the normal distribution 𝒩(pjk,(0.2pjk)2).Step 4.
Use the Giffler-Thompson algorithm to generate an active scheduleσ2 based on s. Note that, in this process, the processing time of each operation takes its new value, that is, pjk(u).Step 5.
Evaluate the objective value ofσ2 under the original processing times, obtaining TWTσ2.Step 6.
IfTWTσ2<TWTσ1, exit the local search. Otherwise, continue with Step 7.Step 7.
Letu←u+1. If u>U, terminate the local search procedure with no improvement found. Otherwise, go back to Step 3.The local search procedure attempts to find a better schedule from an initial starting solution. When the processing time of each operation varies due to the random perturbations in Step3, different schedules may be obtained by applying the same set of dispatching rules. This is the fundamental idea of the local search procedure. The relative deviation is set as 0.2 times the mean, which can produce a moderate level of perturbation. If the variance is too large or too small, the local search will be inefficient. Finally, if it turns out that TWTσ2<TWTσ1, which means that the local search has found an improvement, then the original schedule should be accordingly revised. On the other hand, if no better schedule is found within U trials, the local search will quit, leaving the current solution unchanged.Example 4.2.
Here, we provide a concrete example to illustrate the local search procedure. In the3×3 TWT-JSSP instance, the processing time of each operation is marked beside the corresponding node in Figure 3 (the disjunctive graph). The due dates of each job are set as d1=103, d2=146, and d3=137. The weights of each job are set as w1=1, w2=7, and w3=4.Figure 3
The disjunctive graph of the example instance.Suppose the initial solution iss=[2,2,…], which indicates using the SPT rule at all times.First, the solution is decoded by applying the Giffler-Thompson algorithm under the original values of processing times. The following schedule is obtained:(4.3)σ1=[O11O21O31O12O32O22O13O23O33].
The objective value is TWTσ1=2086 (see Figure 4).Figure 4
The Gantt chart for the initial schedule.Next, we add random perturbations to the processing times and thus generate different new instances of the JSSP. In one of these instances, the processing time of each operation is as follows:(4.4)p11=29.0,p12=85.0,p13=58.2,p21=52.3,p23=81.9,p22=65.3,p32=78.8,p31=84.0,p33=50.3.
The solution s is decoded again by applying the Giffler-Thompson algorithm under the above new values of pjk. The following schedule can be obtained:(4.5)σ2=[O11O21O31O32O22O12O23O33O13],
when evaluated under the original processing times, the objective value of σ2 is TWTσ2=1264 (see Figure 5).Figure 5
The Gantt chart for the new schedule generated by random perturbation.Therefore, an improvement has been found on the initial solution (because1264<2086).


## 4.1. Encoding

The proposed PMBGA relies on dispatching rules to record the scheduling policies. Eight scheduling rules are involved in this study. In the following expressions for the priority index, operationi belongs to job j, and wj and dj are the corresponding job’s weight and due date. JS(i) and JP(i), respectively, denote the set of job successors of operation i and the set of job predecessors of operation i. Z is used to indicate that the operation with the smallest index will be chosen from the conflict set, while Z′ is used to indicate that the operation with the largest index will be chosen from the conflict set:(1)
ATC (apparent tardiness cost):Zi′(t)=(wj/pi)×exp{-(dj-ri-pi - ∑i′∈JS(i)(Ŵi′+pi′))+/(K×p¯)} (in this expression, ri denotes the earliest starting time of operation i, and p¯ denotes the average processing time of the current operations in the conflict set. K is a scaling parameter (or called “look-ahead” parameter). Ŵi is the estimated lead time of operation i);(2)
SPT (shortest processing time):Zi=pi;(3)
LPT (longest processing time):Zi′=pi;(4)
WSPT (weighted shortest processing time):Zi′=wj/pi;(5)
SRPT (shortest remaining processing time):Zi=pi+∑i′∈JS(i)pi′;(6)
LRPT (longest remaining processing time):Zi′=pi+∑i′∈JS(i)pi′;(7)
EDD (earliest due date):Zi=dj;(8)
ODD (operation due date):Zi(t)=[dj/∑i′∈JS(i)∪{i}pi′]×∑i′∈JP(i)∪{i}pi′.In PMBGA, each encoding digit (i.e., gene) is expressed by the serial number (1~8) of the selected dispatching rule. A solution is represented by a rule sequence {Rik:i=1,…,n,k=1,…,m}, where Rik indicates the rule to use when scheduling the ith operation on machine k. Therefore, the encoding length for each solution is l=n×m.


## 4.2. Decoding

In order to evaluate the fitness of a solution, the Giffler-Thompson algorithm is applied to construct an active schedule based on the specified dispatching rules. The implementation of the Giffler-Thompson algorithm is detailed below.Input:  A sequence of rules {Rik:i=1,…,n,k=1,…,m}.Step 1.
LetQ(1)=O={1,…,nm} (the set of all operations), R(1)=F(O)={f1,…,fn} (the set of first operations of each job). Set t=1 and πk=1 (k=1,…,m).Step 2.
Find the operationi*=argmini∈R(t){ri+pi}, and let m* be the index of the machine on which this operation should be processed. Use B(t) to denote all the operations from R(t) which should be processed on machine m*.Step 3.
Delete fromB(t) the operations that satisfy ri⩾ri*+pi*.Step 4.
Use ruleRπm*m* to identify an operation ô from B(t) if currently there are more than one candidates. Schedule operation ô on machine m* at the earliest possible time. Let πm*←πm*+1.Step 5.
LetQ(t+1)=Q(t)∖{ô}, R(t+1)=R(t)∖{ô}∪{suc(ô)}, where suc(ô) denotes the immediate job successor of operation ô (if any).Step 6.
IfQ(t+1)≠∅, set t←t+1 and go to Step 2. Otherwise, the decoding procedure is terminated.In the above description, the release timeri equals the earliest possible starting time of operation i (determined from the already scheduled operations). So, (ri+pi) is the earliest possible completion time of operation i. Q(t) represents the set of operations yet to be scheduled at iteration t, while R(t) represents the set of ready operations (whose job predecessors have all been scheduled) at iteration t. In Step 4, the operation set B(t) is also called a conflict set.


## 4.3. Producing Offspring Individuals

(a) Selection
In each iteration, we first sort all the individuals in the current population according to their fitness. Then, we select the best1/4 of individuals to form the set S, which will subsequently be used to build the Bayesian network in order to produce a new generation of individuals.(b) The Adopted Network Structure
A Bayesian network has to be built to describe the probabilistic distribution of favorableRik settings based on the elite individuals in S. Each individual can be characterized by a directed acyclic network as shown in Figure 1. In this network, a node Ni,kr (r∈{1,2,…,8}) indicates the fact that the dispatching rule used to schedule the ith operation on machine k is selected as rule Number r. The directed arc from node Ni,kr to node Ni,k+1r′ (or Ni+1,1r′) represents the dependency between the two nodes, so it characterizes the possible influence of the rule selection for the ith operation on machine k on the rule selection for the ith operation on machine k+1 (or the (i+1)-th operation on machine 1). Therefore, a directed path from a certain node in the first row to a certain node in the (n×m)-th row can completely describe an individual in the population (because a directed path records an instantiation of all the rule selections).Figure 1
The structure of the Bayesian network used in PMBGA.(c) Calculation of the Probability Values
Since we adopt a fixed network structure in PMBGA, building the Bayesian network is equivalent to determining the values of all the conditional probabilities according to the selected solution setS. After that, new individuals will be produced by iteratively sampling these probabilistic distributions, expecting to obtain high-quality offsprings.Given a number of individuals (i.e., the training setS), an estimate of the conditional probabilities can be obtained simply by counting the frequencies of occurrence.Example 4.1.
Here, we provide a concrete example to illustrate the probability calculation process. For a PMBGA optimization instance withn=3,m=1 and r=3, let us further suppose the current S contains 40 individuals. In Figure 2, the statistics of these individuals are displayed on a network as previously defined. The weight of each arc (the number placed on the arc) indicates the occurring frequency of this arc (counted for all the individuals in S). For example, if there exists an individual coded as “[3,1,2]”, then the path “N1,13→N2,11→N3,12” in the Bayesian network is used to record this individual, and consequently, the weights (counted frequencies) of the arcs N1,13→N2,11 and N2,11→N3,12 will be increased by 1, respectively. Note that, in the final network, the sum of the weights of all the incoming arcs of a certain node should be equal to the sum of the weights of all the outgoing arcs of the same node. This is because each individual corresponds to a complete path from the first row to the last row.Figure 2
An example network in PMBGA.By using frequency as an approximation for probability, the relevant (conditional) probabilities should be calculated as follows:(4.1)P(N1,11)=8+240,P(N1,12)=7+3+440,P(N1,13)=5+5+640,P(N2,11∣N1,11)=88+2,P(N2,12∣N1,11)=0,P(N2,13∣N1,11)=28+2,P(N2,11∣N1,12)=77+3+4,P(N2,12∣N1,12)=37+3+4,P(N2,13∣N1,12)=47+3+4,P(N2,11∣N1,13)=55+5+6,P(N2,12∣N1,13)=55+5+6,P(N2,13∣N1,13)=65+5+6,P(N3,11∣N2,11)=22+3+15,P(N3,12∣N2,11)=32+3+15,P(N3,13∣N2,11)=152+3+15,P(N3,11∣N2,12)=22+6,P(N3,12∣N2,12)=62+6,P(N3,13∣N2,12)=0,P(N3,11∣N2,13)=0,P(N3,12∣N2,13)=77+5,P(N3,13∣N2,13)=57+5.According to the above calculation method, a connection can never be rediscovered in the PMBGA if the corresponding conditional probability is zero (e.g., fromN2,12 to N3,13). To overcome this drawback, we can set the minimum count to 1. Taking N2,12 as an example, the conditional probabilities for the outgoing arcs will then become(4.2)P(N3,11∣N2,12)=33+7+1,P(N3,12∣N2,12)=73+7+1,P(N3,13∣N2,12)=13+7+1.
Now it is possible to discover N3,13 from N2,12, though the probability is small.(d) The Sampling Process
The sampling process for generating a new individual begins from the root nodes of the Bayesian network. By selecting an outgoing arc at each node based on the calculated conditional probabilities, the whole network can be gradually instantiated.


## 4.4. The Embedded Local Search

The search mechanism of GA guarantees a good performance in the “exploration” of the solution space. However, it has been reported that GA alone cannot achieve satisfactory solution quality for complex optimization problems. Actually, a local search procedure is usually added within the framework of GA in order to provide reliable “exploitation” ability. In this paper, we design such a local optimizer, which attempts to improve the selected solutions. In each iteration of PMBGA, the local search is carried out for the beste% of solutions in the current population. Thus, e is an important parameter for adjusting the frequency of local search and achieving a balance between exploration and exploitation.In the following, we will describe how to perform the local search on a given solutions (a sequence of rules).Step 1.
Use the Giffler-Thompson algorithm to generate an active scheduleσ1 based on s. The objective value of the obtained schedule is TWTσ1.Step 2.
Setu=1.Step 3.
Exert random perturbations on the processing times: generate a new set of processing times{pjk(u)} from the normal distribution 𝒩(pjk,(0.2pjk)2).Step 4.
Use the Giffler-Thompson algorithm to generate an active scheduleσ2 based on s. Note that, in this process, the processing time of each operation takes its new value, that is, pjk(u).Step 5.
Evaluate the objective value ofσ2 under the original processing times, obtaining TWTσ2.Step 6.
IfTWTσ2<TWTσ1, exit the local search. Otherwise, continue with Step 7.Step 7.
Letu←u+1. If u>U, terminate the local search procedure with no improvement found. Otherwise, go back to Step 3.The local search procedure attempts to find a better schedule from an initial starting solution. When the processing time of each operation varies due to the random perturbations in Step3, different schedules may be obtained by applying the same set of dispatching rules. This is the fundamental idea of the local search procedure. The relative deviation is set as 0.2 times the mean, which can produce a moderate level of perturbation. If the variance is too large or too small, the local search will be inefficient. Finally, if it turns out that TWTσ2<TWTσ1, which means that the local search has found an improvement, then the original schedule should be accordingly revised. On the other hand, if no better schedule is found within U trials, the local search will quit, leaving the current solution unchanged.Example 4.2.
Here, we provide a concrete example to illustrate the local search procedure. In the3×3 TWT-JSSP instance, the processing time of each operation is marked beside the corresponding node in Figure 3 (the disjunctive graph). The due dates of each job are set as d1=103, d2=146, and d3=137. The weights of each job are set as w1=1, w2=7, and w3=4.Figure 3
The disjunctive graph of the example instance.Suppose the initial solution iss=[2,2,…], which indicates using the SPT rule at all times.First, the solution is decoded by applying the Giffler-Thompson algorithm under the original values of processing times. The following schedule is obtained:(4.3)σ1=[O11O21O31O12O32O22O13O23O33].
The objective value is TWTσ1=2086 (see Figure 4).Figure 4
The Gantt chart for the initial schedule.Next, we add random perturbations to the processing times and thus generate different new instances of the JSSP. In one of these instances, the processing time of each operation is as follows:(4.4)p11=29.0,p12=85.0,p13=58.2,p21=52.3,p23=81.9,p22=65.3,p32=78.8,p31=84.0,p33=50.3.
The solution s is decoded again by applying the Giffler-Thompson algorithm under the above new values of pjk. The following schedule can be obtained:(4.5)σ2=[O11O21O31O32O22O12O23O33O13],
when evaluated under the original processing times, the objective value of σ2 is TWTσ2=1264 (see Figure 5).Figure 5
The Gantt chart for the new schedule generated by random perturbation.Therefore, an improvement has been found on the initial solution (because1264<2086).


## 5. Computational Results

In order to test the performance of the proposed PMBGA, the same benchmark instances as in [15, 19] are used in our computational experiment. In these instances, the due date of each job is set as dj=f×Pj, where Pj denotes the total processing time of job j, and f∈{1.6,1.5,1.3} is a coefficient that controls the tightness level of the due date setting. The first 20% of jobs are assigned weighting 4 (very important), the next 60% are assigned weighting 2, and the remaining jobs are assigned weighting 1 (not important). That is, w1=w2=4, w3=⋯=w8=2, and w9=w10=1.Based on extensive computational tests, the algorithm parameters are set as follows. The population sizePS=50; the proportion of individuals selected for local search e%=30%; the maximum number of random perturbations U=100. In the implementation of the ATC rule, we set K=2 and Ŵi=0.4pi.We compare the performance of the proposed PMBGA with the hybrid optimization algorithm PSO-SA [37] and a rule-based genetic algorithm (RBGA) which uses the same dispatching rules for encoding but adopts classical crossover and mutation operators. In order to make the comparisons meaningful, we set a computational time limit for all the algorithms. In the following, the time limit for solving each instance is determined as 60 seconds. Each algorithm is run for 20 independent times on each TWT-JSSP instance. Tables 1, 2, and 3 report the average objective value obtained from the 20 runs. “# opt” indicates the number of times that the optimum (the optimum for each instance (listed in the second column in the tables) is first given by the branch-and-bound algorithm [14] and recently updated by the hybrid genetic algorithm [21]) has been reached during the 20 runs. The results demonstrate the superiority of the PMBGA.Table 1
The computational results and comparisons underf=1.6.
InstanceOptimalPMBGAPSO-SARBGAAverage# optAverage# optAverage# optABZ500.0200.0200.020ABZ600.0200.0200.020LA1600.0200.0200.318LA176565.02065.02065.020LA1800.0200.0200.020LA1900.0200.0200.818LA2000.0200.0200.020LA2100.0200.0200.020LA2200.0200.2190.020LA2300.0200.0200.020LA2400.0200.0200.020MT10141143.116145.914149.911ORB1566588.74602.01630.70ORB24444.51845.71345.612ORB3422425.210428.88452.18ORB46666.31766.71871.013ORB5163166.06168.94171.81ORB62828.51528.91029.112ORB700.0200.7150.713ORB8621633.42646.70661.40ORB96666.91668.71271.19ORB107683.2283.5084.20Table 2
The computational results and comparisons underf=1.5.
InstanceOptimalPMBGAPSO-SARBGAAverage# optAverage# optAverage# optABZ56969.02069.31869.517ABZ600.0200.0200.517LA16166166.020167.915166.917LA17260260.020260.719262.115LA183434.61834.51834.718LA192121.02021.11921.316LA2000.0200.1190.714LA2100.0200.0200.615LA22196196.020196.517197.115LA2322.0202.0202.020LA248282.11985.91287.58MT10394394.219407.77414.51ORB110981108.4131131.511163.00ORB2292292.020299.15310.92ORB3918922.910937.76989.00ORB4358360.117366.23381.21ORB5405405.318414.87433.82ORB6426427.115441.88455.54ORB75050.11652.3753.46ORB810231032.3141060.2101082.56ORB9297308.812313.89329.25ORB10346377.13388.70396.80Table 3
The computational results and comparisons underf=1.3.
InstanceOptimalPMBGAPSO-SARBGAAverage# optAverage# optAverage# optABZ514031415.1121458.911466.91ABZ6436436.020446.95456.32LA1611691170.2191212.331239.50LA17899900.417914.77942.53LA18929943.015947.62960.11LA19948950.618967.010989.96LA20805822.114829.38840.54LA21463466.815480.012490.57LA2210641095.7101103.981120.93LA23835837.512866.42867.80LA24835844.013855.64860.51MT1013631375.0171403.911436.40ORB125682613.352617.272718.62ORB214081440.541452.031494.71ORB321112132.982154.832205.60ORB416231637.1111656.151707.12ORB515931617.431638.731673.80ORB617901835.991837.371884.03ORB7590592.317602.210620.05ORB824292483.102509.602563.70ORB913161316.7191359.6121393.88ORB1016791683.9171709.991747.13According to the computational results, the proposed PMBGA systematically outperforms the comparative methods. In addition, the following comments can be made:(1)
The proposed PMBGA performs better than the PSO-SA which adopts operation permutation-based encoding scheme. The advantage of PMBGA is even stronger when the due dates in the TWT-JSSP instances are set tighter. This suggests that the rule-based optimization approach is more effective than sequence-based search when many jobs are prone to be tardy and thus are competing fiercely for the limited machine resources. Applying dispatching rules turns out to be a satisfactory and robust strategy in this situation;(2)
Also, PMBGA outperforms RBGA to a greater extent. This reveals the effectiveness of the proposed approach from two aspects:
(i)
Since RBGA does not involve a local search module, the results show that the specialized local search procedure can help to promote the overall performance of GA. In particular, the local optimizer in PMBGA is closely based on the specific characteristics of the considered JSSP instance: the magnitude of the random sampling is consistent with the processing time (cf. Step3 of the procedure), which ensures a reasonable size of the search scope.(ii)
The results show that using the estimation of distribution principle to optimize the rule combinations is more effective than the traditional crossover and mutation operators. Noticeably, the essential point in this process is to model the relationship and interactions between the rule selections for different jobs and different machines.Meanwhile, we also test the impact of the parametere on the final solution quality of PMBGA. A reasonable selection of e will result in an effective balance between exploration and exploitation. In the following experiment, instance ORB1 under f=1.5 is used and the time limit is set as 40 sec and 60 sec, respectively. The computational results are displayed in Figure 6, where the vertical axis gives the average objective value obtained from 20 independent executions of the proposed PMBGA under each e.Figure 6
The influence of parametere.According to the results, the setting ofe has a considerable impact on the solution quality, especially when the computational time is scarce (40 sec), which verifies that the proposed local search module is effective in accelerating the overall convergence of PMBGA. A small e means that only a few solutions in each generation can be improved by the local search, which has little effect on the entire population. A large e suggests that too much time is consumed on local search, which may impair the normal function of PMBGA because of the reduced generations.The best setting ofe under each constraint level is 60 (for tight time budget) and 30 (for loose time budget). When the exogenous restriction on computational time is tight, PMBGA has to rely on frequent local search to find good solutions. This is because, in the short term, local search is more efficient than PMBGA’s mechanism (Bayesian network modeling) in improving a solution. However, the price to pay is possibly a premature convergence of the whole optimization process. On the other hand, when the computational time is more sufficient, PMBGA will prefer a larger number of generations to conduct a systematic exploration of the solution space. In this case, the local search need not be used very frequently, otherwise the steady searching and learning process may be disturbed.Finally, we observe the impact of the parameters in the ATC rule. We write the estimated waiting timeŴi of operation i as proportional to its processing time: Ŵi=b×pi. Now, we use instance MT10 with f=1.3 to test the two parameters, K and b. The average objective values obtained by PMBGA under each (K,b) combination is shown in Figure 7.Figure 7
The influence of parameters(K,b).In fact, the influence of the ATC parameters is not so significant as the other parameters. But a clearly inferior setting isK=1, which could be eliminated. Based on additional experimental results which are not listed here, both K=2 and K=3 seem acceptable. Overall, setting K=2 and b=0.4 yield satisfactory solution quality for most scheduling instances involved in this study.


## 6. Conclusion

In this paper, we propose a new probabilistic model-building genetic algorithm for solving the job shop scheduling problem under the total weighted tardiness criterion. Since the TWT objective is systematically more difficult to optimize than the conventional makespan objective, we rely on some dispatching rules for schedule construction. PMBGA is used to search for good combinations of these rules, and a specific local search algorithm is embedded into the optimization framework. The computational experiments have shown the superiority of the proposed methodology.Future research can be conducted from the following aspects. Although the standard Bayesian network is capable of modeling the interactions between decision variables, it is necessary to improve the computational efficiency in the rule optimization process. Meanwhile, it is interesting to try other encoding schemes, which may be beneficial for the discovery of useful structural properties and may enhance the overall performance of the PMBGA.



---

*Source: 623230-2012-03-08.xml*