- Email: [email protected]
A prediction based iterative decomposition algorithm for scheduling large-scale job shops
Mathematical and Computer Modelling 47 (2008) 411–421 www.elsevier.com/locate/mcm
A prediction based iterative decomposition algorithm for scheduling large-scale job shops Min Liu ∗ , Jing-Hua Hao, Cheng Wu Department of Automation, Tsinghua University, Beijing 100084, PR China Received 28 February 2007; accepted 22 March 2007
Abstract In this paper, we present a prediction based iterative decomposition algorithm for solving large-scale job shop scheduling problems using the rolling horizon scheme and the prediction mechanism, in which the original large-scale scheduling problem is iteratively decomposed into several sub-problems. In the proposed algorithm, based on the job-clustering method, we construct the Global Scheduling Characteristics Prediction Model (GSCPM) to obtain the scheduling characteristics values, including the information of the bottleneck jobs and the predicted value of the global scheduling objective. Then, we adopt the above scheduling characteristics values to guide and coordinate the process of the problem decomposition and the sub-problem solving. Furthermore, we propose an adaptive genetic algorithm to solve each sub-problem. Numerical computational results show that the proposed algorithm is effective for large-scale scheduling problems. c 2007 Elsevier Ltd. All rights reserved.
Keywords: Large-scale; Job shop; Scheduling; Prediction; Decomposition
1. Introduction The job shop scheduling problem is one of the most complicated combinatorial optimization problems [1–3], and it has long been a significant research field since Akers and Jackson [4,5] proposed some efficient methods for 2×m and n × 2 job shop scheduling problems. From then on, many effective exact methods, including mathematical programming methods [6–8], etc., have been proposed to solve job shop scheduling problems with smaller scale. Since the job shop scheduling problem is a kind of typical NP-hard problem, the above exact methods are not suitable for solving job shop scheduling problems with larger scale. The dispatching rules [9,10] and the near-optimal methods, such as the genetic algorithm [11–13], tabu search [14–16] and simulated annealing [17–19], have become very popular and have gained much success in solving job shop scheduling problems with larger scale. However, in practical manufacturing environments the scale of job shop scheduling problems could be much larger; for example, in some big textile factories, the number of jobs may be up to 1000. For large-scale scheduling problems, the performances of the above scheduling methods, including exact methods, dispatching rules and near-optimal methods, are not very satisfactory. Decomposition based scheduling methods have long been effective and popular in handling larger-scale scheduling problems [20–27]. One of the key points in applying the above decomposition methods to solve large-scale scheduling ∗ Corresponding author.
E-mail address: [email protected] (M. Liu). c 2007 Elsevier Ltd. All rights reserved. 0895-7177/$ - see front matter doi:10.1016/j.mcm.2007.03.032
412
M. Liu et al. / Mathematical and Computer Modelling 47 (2008) 411–421
problems is how to evaluate globally a given schedule of each sub-problem. Traditional decomposition based methods in the literature, such as shifting bottleneck procedures [21], rolling horizon based methods [27] and graph decomposition based methods [24], evaluate a given schedule of the sub-problem using either some local scheduling indices or the lower/upper bound of the global scheduling objective. However, for large-scale job shop scheduling problems, the performance of decomposition algorithms based on the above two evaluation methods is usually not satisfactory because the local scheduling indices or the bounds of the global scheduling objective cannot reflect the global characteristics effectively in most cases. In this paper, aiming at solving large-scale job shop scheduling problems with the objective of minimizing the total tardiness, we propose a prediction based iterative decomposition algorithm (PIDA). In the PIDA, the original scheduling problem is iteratively decomposed into several sub-problems, each of which corresponds to a time window respectively, and each sub-problem is iteratively defined and solved with the rolling horizon scheme and the prediction mechanism. In the course of defining and solving each sub-problem, we first construct the Global Scheduling Characteristics Prediction Model (GSCPM) to obtain the scheduling characteristics values for guiding and coordinating the process of the problem decomposition and the sub-problem solving in order that the global scheduling performance can be improved. Then, an adaptive genetic algorithm based on the GSCPM is proposed to solve the current sub-problem, and the partial scheduling policy corresponding to the current sub-problem solved is fixed. After that, the next sub-problem will be defined and solved according to the above procedure until all the subproblems are solved. Finally, an overall schedule which is composed of all the fixed partial scheduling policies of the sub-problems can be obtained. The main difference between traditional rolling horizon based decomposition algorithms (take the RHA (Rolling Horizon Algorithm) proposed by Marcos Singer [27] as an example) and the proposed PIDA can be illustrated in three aspects: (1) in the RHA, all the sub-problems have been defined before the first sub-problem is solved, while in the PIDA, we define the current sub-problem only after the former sub-problem is already solved; (2) in the RHA, the schedule of each sub-problem is evaluated by a lower bound of the global scheduling objective value, while in the PIDA, we evaluate the schedule of each sub-problem by the predicted value of the global scheduling objective; (3) in the RHA, each sub-problem is solved by a modified shifting bottleneck heuristic, which has an exponential computational complexity, while in the PIDA, we solve each sub-problem with an adaptive genetic algorithm based on the GSCPM. Numerical computational results on large-scale job shop scheduling problem instances up to 1000 jobs and 20 machines indicate that the PIDA is effective and outperforms the RHA and some efficient dispatching rules (EDD/SLK/MDD). Also, the results of this paper are promising for large-scale job shop scheduling problems with due-date-related objectives in practical manufacturing environments. This paper is organized as follows. Problem formulation is given in Section 2. Section 3 describes the detailed PIDA. Numerical computations are made in Section 4. Finally, conclusions are drawn in Section 5. 2. Problem formulation The discussed job shop scheduling problem is a typical resource optimization problem with the following assumptions under a multistage manufacturing environment. There is no machine breakdown and no preemption of a processing unit. All jobs are released at zero time and each job has its own due date. The transportation time of jobs between different machines and the setup time of jobs are neglected. Each machine can process only one job at a time and each job can be processed by only one machine at a time. Based on the above assumptions, the job shop scheduling problem discussed in this paper can be formulated in detail as follows. The number of jobs is n and the number of machines is m. Job Ji (i = 1, 2, . . . , n) has its own due date di . Job Ji (i = 1, 2, . . . , n) consisting of m operations Oi,1 , Oi,2 , . . . , Oi,k , . . . , Oi,m must be processed on all the m machines in a predetermined order, where operation Oi,k is to be processed on machine Mk and its processing time is pi,k . Let Ci be the completion time of job Ji , the tardy time of job Ji can be defined as Ti = max(Ci − di , 0). The scheduling objective considered Pn in this paper is to determine the processing sequence of operations on each machine such that the total tardiness ( i=1 Ti ) is minimized without violating the precedence and capacity constraints: min
n X i=1
Ti = min
n X i=1
max (Ci − di , 0) .
413
M. Liu et al. / Mathematical and Computer Modelling 47 (2008) 411–421
The above job shop scheduling problem can also be described as follows. min
n X
max (Ci − di , 0)
i=1
s.t. Oi,k ∈ O sti,k ≥ 0 sti,k2 − sti,k1 ≥ pi,k1 Oi,k1 , Oi,k2 ∈ Ai , k1 = 6 k2 sti1 ,k − sti2 ,k ≥ pi2 ,k ∨ sti2 ,k − sti1 ,k ≥ pi1 ,k Oi1 ,k , Oi2 ,k ∈ E k , Mk ∈ M, i 1 6= i 2
(a) (b)
where sti,k : start time of operation Oi,k ; O: set of all the operations; J : set of all the jobs; M: set of all the machines; Ai : set of all pairs of adjacent operations constrained by the technological constraint for job Ji , i = 1, 2, . . . , n; (Oi,k1 , Oi,k2 ) ∈ Ai indicates that operation Oi,k1 is the immediate predecessor of operation Oi,k2 in the technological route of job Ji ; E k : set of all pairs of the operations to be processed on machine Mk , k = 1, 2, . . . , m. Constraint (a) represents that the processing sequence of each job’s operations follows the corresponding technological constraint, and constraint (b) guarantees that each machine can process only one job at a time. 3. Algorithm 3.1. Solution framework In the proposed prediction based iterative decomposition algorithm, the original scheduling problem is iteratively decomposed into several sub-problems with the rolling horizon scheme and the prediction mechanism. Each time window corresponds to a sub-problem, and the number of time windows is given in advance; that is, the number of operations in each sub-problem is predetermined, but which operations belong to each sub-problem is not predetermined. Each sub-problem is iteratively defined and solved. The solution framework mainly consists of two parts: (1) Constructing the Global Scheduling Characteristics Prediction Model (GSCPM) In the PIDA, the GSCPM, which is used for obtaining the scheduling characteristics values including the information of the bottleneck jobs and the predicted value of the global scheduling objective, needs to be constructed before each sub-problem is defined and solved. The above predicted scheduling characteristics values are used to guide and coordinate the process of the problem decomposition and the sub-problem solving in order that the global scheduling performance can be improved. In the course of constructing the GSCPM, all the unfinished jobs at the start time of the current time window are first clustered into a number of job groups according to their characteristics attributes, and then the RC rule is defined. Based on the above job-clustering results and the RC rule, the calculation formula for predicting each unfinished job’s completion time corresponding to the current sub-problem is constructed. Given a schedule of the current subproblem, each unfinished job’s completion time can be predicted based on the GSCPM. Thus, we can obtain the predicted value of the global scheduling objective and identify the bottleneck jobs. In the PIDA, such a class of jobs with the following characteristics is defined as the bottleneck jobs: their processing priorities have greater influence on the global scheduling objective value (the total tardiness) than those of the others. Based on the above definition and the predicted completion time of each unfinished job, a heuristic (called IBH) is proposed to identify the bottleneck jobs. (2) Defining and solving the sub-problem In the PIDA, based on the above constructed GSCPM, we iteratively define and solve each sub-problem with the rolling horizon scheme and the adaptive genetic algorithm. Considering that the number of sub-problems is predetermined, we need to determine only which operations are to be scheduled in the current time window (the current sub-problem) and optimize their scheduling policy; this is conducted by an adaptive genetic algorithm based on the GSCPM.
414
M. Liu et al. / Mathematical and Computer Modelling 47 (2008) 411–421
3.2. Detailed algorithm 3.2.1. Constructing the Global Scheduling Characteristics Prediction Model (GSCPM) The construction of GSCPM includes three modules: (a) clustering jobs; (b) predicting the completion time of each unfinished job; and (c) identifying the bottleneck jobs; these are detailed as follows. (a) Clustering jobs Since the scheduling objective is to minimize the total tardiness in the discussed job shop scheduling problem, the total remaining processing time and the slack time are two important scheduling characteristics; we adopt the above two scheduling characteristics to form the characteristics vector Vi (t) of job Ji (Ji ∈ Ju ): SLK(Ji , t) max(di − RPT(Ji , t) − t, 0) Vi (t) = = (1) w × RPT(Ji , t) w × RPT(Ji , t) where Ju is the set of unfinished jobs at time t, Ju ⊆ J , |Ju | = n, ¯ t is the start time of the current time window, n is the number of unfinished jobs at time t, RPT(Ji , t) and SLK(Ji , t) are the total remaining processing time and the modified slack time of job Ji at time t, respectively, and w is a weighting parameter whose value depends upon problem characteristics. Vi (t) can reflect the processing priority of job Ji to some extent. Based on the above characteristics vector Vi (t), we define the processing priority distance between two jobs Ji , J j as
d Vi (t), V j (t) = Vi (t) − V j (t) 2 q 2 2 SLK(Ji , t) − SLK J j , t + w × RPT(Ji , t) − w × RPT(J j , t) . (2) = Considering that d Vi (t), V j (t) describes the processing priority difference between two jobs Ji , J j , by means of job clustering based on the above defined priority distance, we can classify the unfinished jobs into job groups with different processing priorities. The job clustering is conducted through a fuzzy c-means algorithm (FCM), which is detailed as follows: 1. Initialize membership matrix U = {u iP j }n×c , where u i j is the membership value of job Ji belonging to job cluster j, Ji ∈ Ju , j = 1, 2, . . . , c, such that cj=1 u i j = 1. 2. Calculate the centroid vector C V j (t) of each job cluster by P (u i j )b Vi (t) C V j (t) =
Ji ∈Ju
P
(u i j )b
j = 1, 2, . . . , c,
(3)
Ji ∈Ju
where the fuzziness index b is a real number greater than 1. 3. Update the membership value u i j by 1/(b−1) 1/d 2 Vi (t), C V j (t) ui j = c Ji ∈ Ju ; j = 1, 2, . . . , c. 1/(b−1) P 2 1/d Vi (t), C V j (t)
(4)
j=1
4. Calculate the clustering error J f : Jf =
c X X
(u i j )b d 2 (Vi (t), C V j (t))
Ji ∈ Ju ; j = 1, 2, . . . , c.
(5)
j=1 Ji ∈Ju
Steps 2–4 are repeated until J f does not decrease in a given number of iterations. By the above procedure, an optimized membership matrix U = {u i j }n×c can be obtained. Then, each unfinished job is assigned to the job group with the maximal membership value and all the unfinished jobs can be classified into c job groups with different processing priorities. In the FCM, the number of job groups c is predetermined according to problem characteristics.
M. Liu et al. / Mathematical and Computer Modelling 47 (2008) 411–421
415
(b) Predicting the completion time of each unfinished job First, a dispatching rule called the RC rule is defined. Then, the jobs’ operations whose scheduling policies have been fixed before the current sub-problem is solved are excluded from the original scheduling problem, and the corresponding remaining technological route of each unfinished job is kept. Based on the above steps, the job shop scheduling problem for prediction (called PSP) is formed and the RC rule is adopted to solve the PSP in order that the optimized value of each unfinished job’s completion time corresponding to the q-th time window can be obtained, where q = 1, 2, . . . , Q and Q is the total number of time windows. The PSP corresponding to the q-th time window is formulated as follows. X X min max(CiP − di , 0) + max(Ci − di , 0) q
q
Ji ∈J /Ju
Ji ∈Ju
s.t. q Oi,k ∈ O q sti,k ≥ t q sti,k2 − sti,k1 ≥ pi,k1 Oi,k1 , Oi,k2 ∈ Ai , k1 = 6 k2 q sti1 ,k − sti2 ,k ≥ pi2 ,k ∨ sti2 ,k − sti1 ,k ≥ pi1 ,k Oi1 ,k , Oi2 ,k ∈ E k , Mk ∈ M, i 1 6= i 2 where t q : start time of the q-th time window; q CiP : the completion time of unfinished job Ji in the q-th PSP, Ji ∈ Ju ; O q : set of all the unscheduled operations at time t q ; n q : number of unfinished jobs at time t q ; q q q Ju : set of all the unfinished jobs at time t q , Ju ⊆ J , |Ju | = n q ; q q Ai : set of all pairs of unfinished adjacent operations belonging to job Ji at time t q , Ji ∈ Ju ; q
E k : set of all pairs of unfinished operations on machine Mk at time t q , k = 1, 2, . . . , m; P q q Ji ∈J /Ju max(Ci − di , 0): total tardiness of finished jobs at time t . p We construct the RC rule based on the 2-norm kVi (t)k2 of characteristics vector Vi (t) in which kVi (t)k2 = SLK2 (Ji , t) + (w × RPT(Ji , t))2 , t is the scheduling point, RPT(Ji , t) and SLK(Ji , t) are the total remaining processing time and the modified slack time of job Ji at time t, respectively. Definition 3.1. At each scheduling point, the operation with the minimal kVi (t)k2 has the highest priority; the operation with the second minimal kVi (t)k2 has the second highest priority, and so on. Also, those operations with the same kVi (t)k2 have the same priority, and will be selected randomly to be processed. Based on the above RC rule and job clustering results, the calculation formula of each unfinished job’s predicted completion time corresponding to the current time window (the q-th time window) can be formed through the following steps: 1. Obtain CiP∗ by means of solving the PSP with the RC rule, where CiP∗ is the optimized value of unfinished job’s completion time CiP corresponding to the current time window; 2. Compute the delay coefficient δ j of job group j by P CiP∗ − t Ji ∈Z j
δj = P
RPT(Ji , t)
j = 1, 2, . . . , c,
(6)
Ji ∈Z j
where Z j denotes job group j. 3. Based on the delay coefficient δ j of job group j ( j = 1, 2, . . . , c), each unfinished job’s predicted completion time C i can be calculated by C i = t + δki × RPT(Ji , t)
q
Ji ∈ Ju ,
where ki is the job group containing job Ji , δki is the delay coefficient of job group ki and ki = 1, 2, . . . , c.
(7)
416
M. Liu et al. / Mathematical and Computer Modelling 47 (2008) 411–421
Since each unfinished job’s predicted completion time corresponding to the current time window can be calculated by (7), the corresponding global scheduling objective value can also be predicted. (c) Identifying the bottleneck jobs After each unfinished job’s predicted completion time corresponding to the current time window is obtained, a heuristic (IBH) is proposed to identify the bottleneck jobs from the unfinished jobs at the start time of the current time window (the q-th time window). The procedure of the IBH is as follows: q
1. Calculate C i of each unfinished job Ji (Ji ∈ Ju ) by (7), where C i is the predicted completion time of unfinished job Ji and n¯ is the total number of unfinished jobs. 2. Calculate L i of each unfinished job Ji by L i = C i − di
q
Ji ∈ Ju ,
(8)
where L i is the predicted lateness of unfinished job Ji . q 3. Sequence all the jobs in a non-increasing order of L i (Ji ∈ Ju ) to obtain a job queue. Given a predefined bottleneck q ratio τ ∈ (0, 1), we define the first bn ×τ c jobs in the job queue as bottleneck jobs, the number of which is denoted q as n b , where τ depends on the problem characteristics such as the values of n q and di (Ji ∈ Ju ), etc. 3.2.2. Defining and solving each sub-problem In the PIDA, based on the number n sub (n sub = n × m/Q) of operations contained in each sub-problem, the length of each time window is determined, where the number Q of sub-problems is given in advance. The q-th sub-problem to be defined and solved can be formulated as follows. X X max(Ci − di , 0) + max(Ci − di , 0) min q
Ji ∈Ju
q
Ji ∈J /Ju
s.t. sti,k ≥ t q Oi,k ∈ O q q Oi,k1 , Oi,k2 ∈ Ai , k1 = 6 k2 , Oi,k ∈ O q sti,k2 − sti,k1 ≥ pi,k1 q sti1 ,k − sti2 ,k ≥ pi2 ,k ∨ sti2 ,k − sti1 ,k ≥ pi1 ,k Oi1 ,k , Oi2 ,k ∈ E k , Mk ∈ M, i 1 6= i 2 , Oi,k ∈ O q q |O | = n sub q Ci = f i t q , sti,k Oi,k ∈ O q where t q : start time of the q-th time window; q Ju : set of all the unfinished jobs at time t q ; q O : set of all the operations whose scheduling policies are to be fixed in the q-th sub-problem; |O q |: number of the elements in the set O q ; q q Ai : set of all pairs of unfinished adjacent operations belonging to job Ji in the q-th sub-problem, Ji ∈ Ju ; q E k : set of all pairs of unfinished operations on Mk in the q-th sub-problem, k = 1, 2, . . . , m; q q f i (t q , sti,k ): predicted completion time of job Ji (Ji ∈ Ju ) by the GSCPM corresponding to the q-th sub-problem. Two key points need to be highlighted in the above formulation of the q-th sub-problem: (1) it is not fixed which operation belong to the q-th sub-problem before the q-th sub-problem is solved; that is, O q may be different for different sub-problems corresponding to the q-th time window; (2) each unfinished job’s completion time, which is used to evaluate the schedule of the q-th sub-problem, is predicted by the GSCPM. Based on the GSCPM, we propose a rule based adaptive genetic algorithm (RBAGA) to define and solve the current sub-problem [28]. Based on the proposed RBAGA, the operations contained in the current sub-problem and their scheduling policies can be determined simultaneously. The proposed RBAGA is detailed as follows. • Coding In RBAGA, we adopt a priority-list/rule mixed coding method to form the chromosome, which can be formulated as
M. Liu et al. / Mathematical and Computer Modelling 47 (2008) 411–421
417
s , s , . . . , s , rule 11 12 1n b 1 .. . Ch = sk1 , sk2 , . . . , skn b , rulek . .. . sm1 , sm2 , . . . , smn b , rulem
•
•
•
•
•
The k-th row (k = 1, 2, . . . , m) represents the scheduling policies of the jobs to be processed on machine Mk in the above chromosome. The chromosome is composed of two sub-chromosomes: the sub-chromosome sk1 , sk2 , . . . , skn b (k = 1, 2, . . . , m) represents the priority list for sequencing the bottleneck jobs that are to be processed on machine Mk (k = 1, 2, . . . , m), and the sub-chromosome rulek represents the dispatching rule for sequencing the other jobs that are to be processed on machine Mk (k = 1, 2, . . . , m). The operations in the current sub-problem corresponding to the chromosome are obtained by the decoding process, and the bottleneck jobs are determined by the corresponding GSCPM. Decoding The decoding process of each chromosome constitutes of two phases: (1) simulation of the sub-problem in the current time window, and (2) prediction of the global scheduling objective value. In phase 1, we adopt the scheduling policies based on the priority lists and the rules in the chromosome to execute the simulation of the current sub-problem, in which the following principles are adopted to compare the processing priority between two jobs: for each machine Mk (k = 1, 2, . . . , m), the processing priorities of the bottleneck jobs are always higher than those of the other jobs, while the processing priorities among the bottleneck jobs are decided by their priority indices in the corresponding priority list in the chromosome, and the processing priorities among the other jobs are decided by the corresponding dispatching rule in the chromosome. Also, it is not permitted for any machine to be idle if there exists any job that can be processed on it at each scheduling point. The simulation of the current sub-problem starts at the beginning of the current time window and it stops when the number of scheduled operations reaches n sub . In phase 2, based on the simulation results in phase 1, the formula (7) for calculating each unfinished job’s predicted completion time and the completion time of each finished job, the predicted value of the global scheduling objective (the total tardiness) can be obtained. The current sub-problems corresponding to the chromosomes have the same number of operations, but the operations in the current sub-problem may be different; these are determined by the above decoding process. Population initialization We use some good dispatching rules to generate the initial priority list of each machine in the chromosomes, and each rulek in the chromosome is selected randomly from the rule set (EDD/MDD/SLK [10]). Evaluation and selection Through the above decoding process, we can obtain the predicted value of the global scheduling objective, which is used to evaluate the corresponding chromosome. Also, we adopt the roulette-wheel method to execute the selection operation. Crossover and mutation We adopt the PMX method [29] to execute the crossover operation between two matched chromosomes, but only the crossover operation between the two sub-chromosomes corresponding to the priority lists is executed, while the crossover operation between the two sub-chromosomes corresponding to the dispatching rules is not executed. For each chromosome, we use a trivial multi-point mutation method to execute the mutation operation of each sub-chromosome, respectively. Adaptive adjustment of bottleneck ratio τ An adaptive adjusting procedure for bottleneck ratio τ is proposed to improve the optimization efficiency of the current sub-problem (the q-th sub-problem), which is detailed as follows: q
(1) Calculate the total number of predicted tardy jobs n T by X q nT = Ui q
Ji ∈Ju
where U i =
n
1 if Ci − di > 0 0 else
unfinished job Ji (Ji ∈
q Ju )
denotes whether job Ji is tardy or not, and Ci is the predicted completion time of in the decoding process.
418
M. Liu et al. / Mathematical and Computer Modelling 47 (2008) 411–421
(2) Calculate 1τ by 1 q min τmax , n T /n q − τ k 1τ = P where 1τ is the increment of τ between two neighboring iterations, τmax is the maximal value of bottleneck ratio τ , τ k is the bottleneck ratio used in the k-th iteration of the genetic algorithm, and P (P > 1) is the parameter for adjusting the bottleneck ratio τ . (3) Update τ by τ k+1 = τ k + 1τ. By the above procedure, the bottleneck ratio τ can be adjusted gradually to approximate the ratio of jobs that are probably to be tardy in order that the scheduling performance can be improved. 4. Numerical computational results and analysis 4.1. Testing problem generation and PIDA parameters We randomly generate the job shop scheduling problem instances with different scale to test the performance of the PIDA. In the generated problem instances, the technological route of each job is a random sequence of m machines, the processing time of each operation ranges uniformly between [1,10], and the due date of job Ji is calculated by Pm di = r (1, 1 + n/20) × β × p , where k=1 i,k Pn is the total number of jobs, r (1, 1 + n/20) is a uniformly random integer number ranges between [1, 1 + n/20], m k=1 pi,k is the sum of processing times of job Ji , and β is the due date tightness. Also, we set β = 1.5 in the generated problem instances. Parameters for constructing the GSCPM and the adaptive genetic algorithm include: (1) fuzziness index b : 2; (2) maximal value of bottleneck ratio τmax : 0.2; (3) Number of generation: 20; (4) Population: 20; (5) Crossover probability: 0.2; (6) Mutation probability: 0.05. Other parameters (including the number of sub-problems and the number of job clusters) are set according to problem scale. Since the above two parameters have remarkable influence on the performance of the PIDA, we have carried out some special experiments to analyze the influence of the above algorithm parameters; these are shown in Section 4.2. 4.2. Numerical computations and analysis We compare the proposed PIDA with some of the most efficient dispatching rules (EDD/SLK/MDD) and the RHA [27] on randomly generated job shop scheduling problem instances with different scale. The results are shown in Table 1, in which we list the average results (the total tardiness and the running time) over five independent runs of each scale. From Table 1, we can see that for job shop scheduling problem instances whose scales are relatively small (the number of jobs is less than 50 and the number of machines is less than 10), the performance of the PIDA is comparative with that of the RHA while the running time of the PIDA is much less than that of the RHA. Besides that, since the subproblem in the RHA is solved by a partial enumeration heuristic whose computational complexity grows exponentially with the number of operations assigned to any of the machines in a given sub-problem, the RHA costs much more time than the PIDA does when the scale of the problem instance is larger than 100 × 20. We compare the PIDA with only dispatching rules on these problem instances with larger scale in this paper. From Table 1 we also can see that the PIDA outperforms the above dispatching rules within acceptable time. Also, in Fig. 1 (take the 300 × 20 instance as an example), we list the best global scheduling objective value predicted by the GSCPM corresponding to each sub-problem in order to explain the prediction effectiveness of the global scheduling objective value. The algorithm parameters are: (1) the number of job clusters: c = 18; (2) the number of sub-problems: Q = 15. We can see from Fig. 1 that the predicted value of the global scheduling objective (the total tardiness) is comparatively stable and decreases gradually with the increase of the sub-problem index, which shows that the proposed GSCPM has higher accuracy in predicting the global scheduling objective value and can reflect the global problem characteristics efficiently. Therefore, the optimization process of each sub-problem can be effectively guided by the predicted value of the global scheduling objective.
419
M. Liu et al. / Mathematical and Computer Modelling 47 (2008) 411–421 Table 1 Performance of the PIDA on problem instances with different scale n×m
10 × 10 30 × 10 50 × 10 100 × 10 100 × 20 300 × 20 500 × 20 800 × 20 1000 × 20
EDD/SLK/MDDa P Ti Time (s)
RHAb P Ti
Time (s)
PIDA P Ti
Time (s)
(Q, c)c
86.4 749.7 2 690.8 9 295.7 5 690.3 7 761.4 9 793.2 16 982.5 25 692.1
67.9 654.0 2292.7 8207.3 4577.5 M M M M
14.0 82.3 492.9 1384.2 3530.0 M M M M
74.1 650.2 2 308.4 7 927.2 4 287.9 6 892.5 8 877.0 14 327.9 21 892.4
2.0 5.2 11.1 38.3 214.4 467.5 951.8 1901.6 3210.4
(1, 3) (2, 4) (4, 7) (6, 11) (9, 12) (15, 18) (16, 26) (21, 37) (25, 51)
0.1 0.4 0.6 4.0 5.1 8.7 18.6 47 97.1
M: Since the RHA requires too much running time (more than 1 h) when the problem scale is larger than 100 × 20, we have not made experiments on these problemPinstances. a The column Ti is the minimal total tardiness obtained by EDD/SLK/MDD, respectively, and the column “Time” is the sum of running time of the above the dispatching rules. b The number of time windows and the overlapping factor used in the RHA are set according to Macros Singer [27]. c The column (Q, c) lists the best parameter combinations of the number of sub-problems and the number of clusters for solving each problem instance. Table 2 Performance of the PIDA with different Q on problem instances with smaller scale n×m
10 × 10 20 × 10 30 × 10 50 × 10
Q=1 P Ti 76a 175a 662 2451
Time (s)
Q=2 P Ti
Time (s)
Q=3 P Ti
Time (s)
Q=4 P Ti
Time (s)
Q=5 P Ti
Time (s)
2 3 3 4
81 176 630a 2408
3 4 5 6
84 189 645 2385
3 6 7 9
92 190 677 2378a
4 6 7 12
103 218 711 2463
5 8 9 14
a Minimal total tardiness of PIDA with different Q on the given problem instance.
Fig. 1. Predicted total tardiness for the 300 × 20 instance.
In order to illustrate the influence of different parameters (Q and c) on the performance of the PIDA, we have made a lot of numerical computations. Tables 2–4 show the performance of the PIDA with different Q and c on job shop problem instances, respectively. Tables 2 and 3 show that the optimal number of sub-problems grows with the problem scale; that is, for the problem instances with smaller scale, two or three sub-problems are suitable, while for the problem instances with larger scale, better global scheduling performance might be obtained when the number of sub-problems is larger. From Table 4 we can see that the obtained total tardiness does not decrease monotonically with the increase of the number of job clusters, which indicates that the constructed GSCPM cannot perform well when the number of job clusters is either too large or too small. Therefore, we must set c in a suitable range in order to obtain better scheduling performance.
420
M. Liu et al. / Mathematical and Computer Modelling 47 (2008) 411–421
Table 3 Performance of the PIDA with different Q on problem instances with larger scale n×m
100 × 10 150 × 10 100 × 20 200 × 10 200 × 20
Q=3 P Ti 8 741 9 723 5 319 11 892 5 793
Time (s) 22 27 56 62 92
Q=5 P Ti 8 403a 9 535 5 110 11 694 5 428
Time (s)
Q=8 P Ti
Time (s)
Q = 10 P Ti
Time (s)
Q = 13 P Ti
Time (s)
31 41 91 89 149
8 592 9 127a 4 911a 10 865 5 461
48 77 159 164 196
8 609 9 154 5 021 10 060a 5 187a
60 112 178 199 251
9 126 9 262 5 220 11 215 5 304
72 136 213 224 332
a Minimal total tardiness of PIDA with different Q on the given problem instance.
Table 4 Performance of the PIDA with different c on problem instances with different scale n×m
100 × 10 150 × 10 200 × 10 300 × 10
c=5 P Ti 8 438 9 601 11 092 15 493
Time (s) 30 75 179 253
c = 10 P Ti 8 258a 9 487 10 847 15 482
Time (s)
c = 15 P Ti
Time (s)
c = 20 P Ti
Time (s)
c = 25 P Ti
Time (s)
32 77 179 255
8 316 9 139a 10 608 14 783
34 78 180 258
8 409 9 267 10 202a 13 392a
35 79 181 258
8 493 9 368 10 512 14 582
35 80 184 257
a Minimal total tardiness of PIDA with different c on the given problem instance.
5. Conclusions In this paper, a prediction based iterative decomposition algorithm for solving large-scale job shop scheduling problem is proposed. In this algorithm, the original large-scale problem is decomposed into several sub-problems, and each sub-problem is iteratively defined and solved with the rolling horizon scheme and the prediction mechanism. Numerical computational results show that PIDA is effective for large-scale scheduling problems. The results of this paper are promising for optimizing large-scale due-date-related objective based scheduling problems in practical manufacturing environments. Acknowledgments This paper is supported by the National 973 Program of China (No. 2002CB312202), the National 863 High-Tech Program of China (No. 2006AA04Z163), the National Natural Science Foundation of China (No. 60004010, No. 60274045 and No. 60443009), the Beijing Science and Technology Planning Program of China (No. D0305005040321) and Program for New Century Excellent Talents in University. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10]
K.R. Baker, Introduction to Sequencing and Scheduling, Wiley, New York, 1974. J.K. Lenstra, K.A. Rinnooy, P. Brucker, Complexity of machine scheduling problems, Annals of Discrete Mathematics 1 (1977) 343–362. M.R. Garey, D.S. Johnson, Computers and Intractability, Freeman, San Francisco, CA, 1979. S.B. Akers, A graphical approach to production scheduling problems, Operations Research 4 (2) (1956) 244–245. J.R. Jackson, An extension of Johnson’s result on job lot scheduling, Naval Research Logistics Quarterly 3 (3) (1956) 201–203. J. Carlier, E. Pinson, An algorithm for solving the job-shop problem, Management Science 35 (2) (1989) 164–176. J.R. Barker, G.B. McMahon, Scheduling the general job-shop, Management Science 31 (5) (1985) 594–598. H.M. Wagner, An integer programming model for machine scheduling, Naval Research Logistics Quarterly 6 (1959) 131–140. S.S. Panwalker, W.A. Iskander, A survey of scheduling rules, Operations Research 25 (1) (1977) 45–61. Ari P.J. Vepsalainen, Thomas E. Morton, Priority rules for job shops with total weighted tardiness costs, Management Science 33 (8) (1987) 1035–1047. [11] R. Cheng, M. Gen, A tutorial survey of job-shop scheduling problems using genetic algorithms part I: Representation, Computers & Industrial Engineering 34 (4) (1996) 983–997. [12] L. Davis, Job-shop scheduling with genetic algorithm, in: J.J. Grefenstette (Ed.), Proceedings of the 1st International Conference on Genetic Algorithms and their Applications, Lawrence Erlbaum, Pittsburgh, PA, USA, 1985, pp. 136–140.
M. Liu et al. / Mathematical and Computer Modelling 47 (2008) 411–421 [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29]
421
F. Della Croce, R. Tadei, G. Volta, A genetic algorithm for the job shop problem, Computers & Operations Research 22 (1) (1995) 15–24. J.W. Barnes, J.B. Chambers, Solving the job shop scheduling problem with tabu search, IIE Transactions 27 (22) (1995) 257–263. M. Dell’Amico, M. Trubian, Applying Tabu search to the job-shop scheduling problem, Annals of Operations Research 41 (1993) 231–252. E. Nowicki, C. Smutnicki, A fast taboo search algorithm for the job-shop problem, Management Science 42 (6) (1996) 797–813. P. Van Laarhoven, E. Aarts, J. Lenstra, Job shop scheduling by simulated annealing, Operations Research 40 (1) (1992) 113–125. T. Yamada, R. Nakano, Job-shop scheduling by simulated annealing combined with deterministic local search, in: Meta-heuristics: Theory and Applications, Kluwer Academic Publishers, MA, USA, 1996, pp. 237–248. H. Matsuo, C.J. Suh, R.S. Sullivan, A controlled search simulated annealing method for the general job shop scheduling problem, Working Paper, 1988, 03-04-88, Graduate School of Business, University of Texas at Austin, Austin, TX, USA. C. Chu, M.C. Portman, J.M. Proth, A splitting-up approach to simplify job shop scheduling problem, International Journal of Production Research 30 (4) (1992) 859–870. J. Adams, E. Balas, D. Zawack, The shifting bottleneck procedure for job shop scheduling, Management Science 34 (3) (1988) 391–401. B.L. Peter, J.H. Debra, Scheduling of manufacturing systems using the lagrangian relaxation technique, IEEE Transactions on Automatic Control 38 (7) (1993) 1066–1080. D. Sun, R. Batta, Scheduling large job shops: A decomposition approach, International Journal of Production Research 34 (7) (1996) 2019–2033. B. Eui-Seok, S.D. Wu, H.S. Robert, Decomposition heuristics for robust job-shop scheduling, IEEE Transactions on Robotics and Automation 14 (2) (1998) 303–313. D.K. Jeffrey, L. John, Mann flowsheet decomposition heuristic for scheduling: A relax-and-fix method, Computers and Chemical Engineering 28 (11) (2004) 2193–2200. G.B. Anis, H. Mohamed, An approximate decomposition algorithm for scheduling on parallel machines with heads and tails, Computers & Operations Research 34 (3) (2007) 868–883. S. Macros, Decomposition methods for large job shops, Computers & Operations Research 28 (3) (2001) 193–207. M. Liu, C. Wu, Genetic algorithm using sequence rule chain for multi-objective optimization in re-entrant micro-electronic production line, Robotics and Computer-Integrated Manufacturing 20 (3) (2004) 225–236. D.E. Goldberg, R. Lingle, Alleles, loci and the traveling salesman problem, in: Proceedings of the International Conference on Genetic Algorithms and their Applications, Lawrence Erlbaum Associates, 1985, pp. 154–159.
A prediction based iterative decomposition algorithm for scheduling large-scale job shops Min Liu ∗ , Jing-Hua Hao, Cheng Wu Department of Automation, Tsinghua University, Beijing 100084, PR China Received 28 February 2007; accepted 22 March 2007
Abstract In this paper, we present a prediction based iterative decomposition algorithm for solving large-scale job shop scheduling problems using the rolling horizon scheme and the prediction mechanism, in which the original large-scale scheduling problem is iteratively decomposed into several sub-problems. In the proposed algorithm, based on the job-clustering method, we construct the Global Scheduling Characteristics Prediction Model (GSCPM) to obtain the scheduling characteristics values, including the information of the bottleneck jobs and the predicted value of the global scheduling objective. Then, we adopt the above scheduling characteristics values to guide and coordinate the process of the problem decomposition and the sub-problem solving. Furthermore, we propose an adaptive genetic algorithm to solve each sub-problem. Numerical computational results show that the proposed algorithm is effective for large-scale scheduling problems. c 2007 Elsevier Ltd. All rights reserved.
Keywords: Large-scale; Job shop; Scheduling; Prediction; Decomposition
1. Introduction The job shop scheduling problem is one of the most complicated combinatorial optimization problems [1–3], and it has long been a significant research field since Akers and Jackson [4,5] proposed some efficient methods for 2×m and n × 2 job shop scheduling problems. From then on, many effective exact methods, including mathematical programming methods [6–8], etc., have been proposed to solve job shop scheduling problems with smaller scale. Since the job shop scheduling problem is a kind of typical NP-hard problem, the above exact methods are not suitable for solving job shop scheduling problems with larger scale. The dispatching rules [9,10] and the near-optimal methods, such as the genetic algorithm [11–13], tabu search [14–16] and simulated annealing [17–19], have become very popular and have gained much success in solving job shop scheduling problems with larger scale. However, in practical manufacturing environments the scale of job shop scheduling problems could be much larger; for example, in some big textile factories, the number of jobs may be up to 1000. For large-scale scheduling problems, the performances of the above scheduling methods, including exact methods, dispatching rules and near-optimal methods, are not very satisfactory. Decomposition based scheduling methods have long been effective and popular in handling larger-scale scheduling problems [20–27]. One of the key points in applying the above decomposition methods to solve large-scale scheduling ∗ Corresponding author.
E-mail address: [email protected] (M. Liu). c 2007 Elsevier Ltd. All rights reserved. 0895-7177/$ - see front matter doi:10.1016/j.mcm.2007.03.032
412
M. Liu et al. / Mathematical and Computer Modelling 47 (2008) 411–421
problems is how to evaluate globally a given schedule of each sub-problem. Traditional decomposition based methods in the literature, such as shifting bottleneck procedures [21], rolling horizon based methods [27] and graph decomposition based methods [24], evaluate a given schedule of the sub-problem using either some local scheduling indices or the lower/upper bound of the global scheduling objective. However, for large-scale job shop scheduling problems, the performance of decomposition algorithms based on the above two evaluation methods is usually not satisfactory because the local scheduling indices or the bounds of the global scheduling objective cannot reflect the global characteristics effectively in most cases. In this paper, aiming at solving large-scale job shop scheduling problems with the objective of minimizing the total tardiness, we propose a prediction based iterative decomposition algorithm (PIDA). In the PIDA, the original scheduling problem is iteratively decomposed into several sub-problems, each of which corresponds to a time window respectively, and each sub-problem is iteratively defined and solved with the rolling horizon scheme and the prediction mechanism. In the course of defining and solving each sub-problem, we first construct the Global Scheduling Characteristics Prediction Model (GSCPM) to obtain the scheduling characteristics values for guiding and coordinating the process of the problem decomposition and the sub-problem solving in order that the global scheduling performance can be improved. Then, an adaptive genetic algorithm based on the GSCPM is proposed to solve the current sub-problem, and the partial scheduling policy corresponding to the current sub-problem solved is fixed. After that, the next sub-problem will be defined and solved according to the above procedure until all the subproblems are solved. Finally, an overall schedule which is composed of all the fixed partial scheduling policies of the sub-problems can be obtained. The main difference between traditional rolling horizon based decomposition algorithms (take the RHA (Rolling Horizon Algorithm) proposed by Marcos Singer [27] as an example) and the proposed PIDA can be illustrated in three aspects: (1) in the RHA, all the sub-problems have been defined before the first sub-problem is solved, while in the PIDA, we define the current sub-problem only after the former sub-problem is already solved; (2) in the RHA, the schedule of each sub-problem is evaluated by a lower bound of the global scheduling objective value, while in the PIDA, we evaluate the schedule of each sub-problem by the predicted value of the global scheduling objective; (3) in the RHA, each sub-problem is solved by a modified shifting bottleneck heuristic, which has an exponential computational complexity, while in the PIDA, we solve each sub-problem with an adaptive genetic algorithm based on the GSCPM. Numerical computational results on large-scale job shop scheduling problem instances up to 1000 jobs and 20 machines indicate that the PIDA is effective and outperforms the RHA and some efficient dispatching rules (EDD/SLK/MDD). Also, the results of this paper are promising for large-scale job shop scheduling problems with due-date-related objectives in practical manufacturing environments. This paper is organized as follows. Problem formulation is given in Section 2. Section 3 describes the detailed PIDA. Numerical computations are made in Section 4. Finally, conclusions are drawn in Section 5. 2. Problem formulation The discussed job shop scheduling problem is a typical resource optimization problem with the following assumptions under a multistage manufacturing environment. There is no machine breakdown and no preemption of a processing unit. All jobs are released at zero time and each job has its own due date. The transportation time of jobs between different machines and the setup time of jobs are neglected. Each machine can process only one job at a time and each job can be processed by only one machine at a time. Based on the above assumptions, the job shop scheduling problem discussed in this paper can be formulated in detail as follows. The number of jobs is n and the number of machines is m. Job Ji (i = 1, 2, . . . , n) has its own due date di . Job Ji (i = 1, 2, . . . , n) consisting of m operations Oi,1 , Oi,2 , . . . , Oi,k , . . . , Oi,m must be processed on all the m machines in a predetermined order, where operation Oi,k is to be processed on machine Mk and its processing time is pi,k . Let Ci be the completion time of job Ji , the tardy time of job Ji can be defined as Ti = max(Ci − di , 0). The scheduling objective considered Pn in this paper is to determine the processing sequence of operations on each machine such that the total tardiness ( i=1 Ti ) is minimized without violating the precedence and capacity constraints: min
n X i=1
Ti = min
n X i=1
max (Ci − di , 0) .
413
M. Liu et al. / Mathematical and Computer Modelling 47 (2008) 411–421
The above job shop scheduling problem can also be described as follows. min
n X
max (Ci − di , 0)
i=1
s.t. Oi,k ∈ O sti,k ≥ 0 sti,k2 − sti,k1 ≥ pi,k1 Oi,k1 , Oi,k2 ∈ Ai , k1 = 6 k2 sti1 ,k − sti2 ,k ≥ pi2 ,k ∨ sti2 ,k − sti1 ,k ≥ pi1 ,k Oi1 ,k , Oi2 ,k ∈ E k , Mk ∈ M, i 1 6= i 2
(a) (b)
where sti,k : start time of operation Oi,k ; O: set of all the operations; J : set of all the jobs; M: set of all the machines; Ai : set of all pairs of adjacent operations constrained by the technological constraint for job Ji , i = 1, 2, . . . , n; (Oi,k1 , Oi,k2 ) ∈ Ai indicates that operation Oi,k1 is the immediate predecessor of operation Oi,k2 in the technological route of job Ji ; E k : set of all pairs of the operations to be processed on machine Mk , k = 1, 2, . . . , m. Constraint (a) represents that the processing sequence of each job’s operations follows the corresponding technological constraint, and constraint (b) guarantees that each machine can process only one job at a time. 3. Algorithm 3.1. Solution framework In the proposed prediction based iterative decomposition algorithm, the original scheduling problem is iteratively decomposed into several sub-problems with the rolling horizon scheme and the prediction mechanism. Each time window corresponds to a sub-problem, and the number of time windows is given in advance; that is, the number of operations in each sub-problem is predetermined, but which operations belong to each sub-problem is not predetermined. Each sub-problem is iteratively defined and solved. The solution framework mainly consists of two parts: (1) Constructing the Global Scheduling Characteristics Prediction Model (GSCPM) In the PIDA, the GSCPM, which is used for obtaining the scheduling characteristics values including the information of the bottleneck jobs and the predicted value of the global scheduling objective, needs to be constructed before each sub-problem is defined and solved. The above predicted scheduling characteristics values are used to guide and coordinate the process of the problem decomposition and the sub-problem solving in order that the global scheduling performance can be improved. In the course of constructing the GSCPM, all the unfinished jobs at the start time of the current time window are first clustered into a number of job groups according to their characteristics attributes, and then the RC rule is defined. Based on the above job-clustering results and the RC rule, the calculation formula for predicting each unfinished job’s completion time corresponding to the current sub-problem is constructed. Given a schedule of the current subproblem, each unfinished job’s completion time can be predicted based on the GSCPM. Thus, we can obtain the predicted value of the global scheduling objective and identify the bottleneck jobs. In the PIDA, such a class of jobs with the following characteristics is defined as the bottleneck jobs: their processing priorities have greater influence on the global scheduling objective value (the total tardiness) than those of the others. Based on the above definition and the predicted completion time of each unfinished job, a heuristic (called IBH) is proposed to identify the bottleneck jobs. (2) Defining and solving the sub-problem In the PIDA, based on the above constructed GSCPM, we iteratively define and solve each sub-problem with the rolling horizon scheme and the adaptive genetic algorithm. Considering that the number of sub-problems is predetermined, we need to determine only which operations are to be scheduled in the current time window (the current sub-problem) and optimize their scheduling policy; this is conducted by an adaptive genetic algorithm based on the GSCPM.
414
M. Liu et al. / Mathematical and Computer Modelling 47 (2008) 411–421
3.2. Detailed algorithm 3.2.1. Constructing the Global Scheduling Characteristics Prediction Model (GSCPM) The construction of GSCPM includes three modules: (a) clustering jobs; (b) predicting the completion time of each unfinished job; and (c) identifying the bottleneck jobs; these are detailed as follows. (a) Clustering jobs Since the scheduling objective is to minimize the total tardiness in the discussed job shop scheduling problem, the total remaining processing time and the slack time are two important scheduling characteristics; we adopt the above two scheduling characteristics to form the characteristics vector Vi (t) of job Ji (Ji ∈ Ju ): SLK(Ji , t) max(di − RPT(Ji , t) − t, 0) Vi (t) = = (1) w × RPT(Ji , t) w × RPT(Ji , t) where Ju is the set of unfinished jobs at time t, Ju ⊆ J , |Ju | = n, ¯ t is the start time of the current time window, n is the number of unfinished jobs at time t, RPT(Ji , t) and SLK(Ji , t) are the total remaining processing time and the modified slack time of job Ji at time t, respectively, and w is a weighting parameter whose value depends upon problem characteristics. Vi (t) can reflect the processing priority of job Ji to some extent. Based on the above characteristics vector Vi (t), we define the processing priority distance between two jobs Ji , J j as
d Vi (t), V j (t) = Vi (t) − V j (t) 2 q 2 2 SLK(Ji , t) − SLK J j , t + w × RPT(Ji , t) − w × RPT(J j , t) . (2) = Considering that d Vi (t), V j (t) describes the processing priority difference between two jobs Ji , J j , by means of job clustering based on the above defined priority distance, we can classify the unfinished jobs into job groups with different processing priorities. The job clustering is conducted through a fuzzy c-means algorithm (FCM), which is detailed as follows: 1. Initialize membership matrix U = {u iP j }n×c , where u i j is the membership value of job Ji belonging to job cluster j, Ji ∈ Ju , j = 1, 2, . . . , c, such that cj=1 u i j = 1. 2. Calculate the centroid vector C V j (t) of each job cluster by P (u i j )b Vi (t) C V j (t) =
Ji ∈Ju
P
(u i j )b
j = 1, 2, . . . , c,
(3)
Ji ∈Ju
where the fuzziness index b is a real number greater than 1. 3. Update the membership value u i j by 1/(b−1) 1/d 2 Vi (t), C V j (t) ui j = c Ji ∈ Ju ; j = 1, 2, . . . , c. 1/(b−1) P 2 1/d Vi (t), C V j (t)
(4)
j=1
4. Calculate the clustering error J f : Jf =
c X X
(u i j )b d 2 (Vi (t), C V j (t))
Ji ∈ Ju ; j = 1, 2, . . . , c.
(5)
j=1 Ji ∈Ju
Steps 2–4 are repeated until J f does not decrease in a given number of iterations. By the above procedure, an optimized membership matrix U = {u i j }n×c can be obtained. Then, each unfinished job is assigned to the job group with the maximal membership value and all the unfinished jobs can be classified into c job groups with different processing priorities. In the FCM, the number of job groups c is predetermined according to problem characteristics.
M. Liu et al. / Mathematical and Computer Modelling 47 (2008) 411–421
415
(b) Predicting the completion time of each unfinished job First, a dispatching rule called the RC rule is defined. Then, the jobs’ operations whose scheduling policies have been fixed before the current sub-problem is solved are excluded from the original scheduling problem, and the corresponding remaining technological route of each unfinished job is kept. Based on the above steps, the job shop scheduling problem for prediction (called PSP) is formed and the RC rule is adopted to solve the PSP in order that the optimized value of each unfinished job’s completion time corresponding to the q-th time window can be obtained, where q = 1, 2, . . . , Q and Q is the total number of time windows. The PSP corresponding to the q-th time window is formulated as follows. X X min max(CiP − di , 0) + max(Ci − di , 0) q
q
Ji ∈J /Ju
Ji ∈Ju
s.t. q Oi,k ∈ O q sti,k ≥ t q sti,k2 − sti,k1 ≥ pi,k1 Oi,k1 , Oi,k2 ∈ Ai , k1 = 6 k2 q sti1 ,k − sti2 ,k ≥ pi2 ,k ∨ sti2 ,k − sti1 ,k ≥ pi1 ,k Oi1 ,k , Oi2 ,k ∈ E k , Mk ∈ M, i 1 6= i 2 where t q : start time of the q-th time window; q CiP : the completion time of unfinished job Ji in the q-th PSP, Ji ∈ Ju ; O q : set of all the unscheduled operations at time t q ; n q : number of unfinished jobs at time t q ; q q q Ju : set of all the unfinished jobs at time t q , Ju ⊆ J , |Ju | = n q ; q q Ai : set of all pairs of unfinished adjacent operations belonging to job Ji at time t q , Ji ∈ Ju ; q
E k : set of all pairs of unfinished operations on machine Mk at time t q , k = 1, 2, . . . , m; P q q Ji ∈J /Ju max(Ci − di , 0): total tardiness of finished jobs at time t . p We construct the RC rule based on the 2-norm kVi (t)k2 of characteristics vector Vi (t) in which kVi (t)k2 = SLK2 (Ji , t) + (w × RPT(Ji , t))2 , t is the scheduling point, RPT(Ji , t) and SLK(Ji , t) are the total remaining processing time and the modified slack time of job Ji at time t, respectively. Definition 3.1. At each scheduling point, the operation with the minimal kVi (t)k2 has the highest priority; the operation with the second minimal kVi (t)k2 has the second highest priority, and so on. Also, those operations with the same kVi (t)k2 have the same priority, and will be selected randomly to be processed. Based on the above RC rule and job clustering results, the calculation formula of each unfinished job’s predicted completion time corresponding to the current time window (the q-th time window) can be formed through the following steps: 1. Obtain CiP∗ by means of solving the PSP with the RC rule, where CiP∗ is the optimized value of unfinished job’s completion time CiP corresponding to the current time window; 2. Compute the delay coefficient δ j of job group j by P CiP∗ − t Ji ∈Z j
δj = P
RPT(Ji , t)
j = 1, 2, . . . , c,
(6)
Ji ∈Z j
where Z j denotes job group j. 3. Based on the delay coefficient δ j of job group j ( j = 1, 2, . . . , c), each unfinished job’s predicted completion time C i can be calculated by C i = t + δki × RPT(Ji , t)
q
Ji ∈ Ju ,
where ki is the job group containing job Ji , δki is the delay coefficient of job group ki and ki = 1, 2, . . . , c.
(7)
416
M. Liu et al. / Mathematical and Computer Modelling 47 (2008) 411–421
Since each unfinished job’s predicted completion time corresponding to the current time window can be calculated by (7), the corresponding global scheduling objective value can also be predicted. (c) Identifying the bottleneck jobs After each unfinished job’s predicted completion time corresponding to the current time window is obtained, a heuristic (IBH) is proposed to identify the bottleneck jobs from the unfinished jobs at the start time of the current time window (the q-th time window). The procedure of the IBH is as follows: q
1. Calculate C i of each unfinished job Ji (Ji ∈ Ju ) by (7), where C i is the predicted completion time of unfinished job Ji and n¯ is the total number of unfinished jobs. 2. Calculate L i of each unfinished job Ji by L i = C i − di
q
Ji ∈ Ju ,
(8)
where L i is the predicted lateness of unfinished job Ji . q 3. Sequence all the jobs in a non-increasing order of L i (Ji ∈ Ju ) to obtain a job queue. Given a predefined bottleneck q ratio τ ∈ (0, 1), we define the first bn ×τ c jobs in the job queue as bottleneck jobs, the number of which is denoted q as n b , where τ depends on the problem characteristics such as the values of n q and di (Ji ∈ Ju ), etc. 3.2.2. Defining and solving each sub-problem In the PIDA, based on the number n sub (n sub = n × m/Q) of operations contained in each sub-problem, the length of each time window is determined, where the number Q of sub-problems is given in advance. The q-th sub-problem to be defined and solved can be formulated as follows. X X max(Ci − di , 0) + max(Ci − di , 0) min q
Ji ∈Ju
q
Ji ∈J /Ju
s.t. sti,k ≥ t q Oi,k ∈ O q q Oi,k1 , Oi,k2 ∈ Ai , k1 = 6 k2 , Oi,k ∈ O q sti,k2 − sti,k1 ≥ pi,k1 q sti1 ,k − sti2 ,k ≥ pi2 ,k ∨ sti2 ,k − sti1 ,k ≥ pi1 ,k Oi1 ,k , Oi2 ,k ∈ E k , Mk ∈ M, i 1 6= i 2 , Oi,k ∈ O q q |O | = n sub q Ci = f i t q , sti,k Oi,k ∈ O q where t q : start time of the q-th time window; q Ju : set of all the unfinished jobs at time t q ; q O : set of all the operations whose scheduling policies are to be fixed in the q-th sub-problem; |O q |: number of the elements in the set O q ; q q Ai : set of all pairs of unfinished adjacent operations belonging to job Ji in the q-th sub-problem, Ji ∈ Ju ; q E k : set of all pairs of unfinished operations on Mk in the q-th sub-problem, k = 1, 2, . . . , m; q q f i (t q , sti,k ): predicted completion time of job Ji (Ji ∈ Ju ) by the GSCPM corresponding to the q-th sub-problem. Two key points need to be highlighted in the above formulation of the q-th sub-problem: (1) it is not fixed which operation belong to the q-th sub-problem before the q-th sub-problem is solved; that is, O q may be different for different sub-problems corresponding to the q-th time window; (2) each unfinished job’s completion time, which is used to evaluate the schedule of the q-th sub-problem, is predicted by the GSCPM. Based on the GSCPM, we propose a rule based adaptive genetic algorithm (RBAGA) to define and solve the current sub-problem [28]. Based on the proposed RBAGA, the operations contained in the current sub-problem and their scheduling policies can be determined simultaneously. The proposed RBAGA is detailed as follows. • Coding In RBAGA, we adopt a priority-list/rule mixed coding method to form the chromosome, which can be formulated as
M. Liu et al. / Mathematical and Computer Modelling 47 (2008) 411–421
417
s , s , . . . , s , rule 11 12 1n b 1 .. . Ch = sk1 , sk2 , . . . , skn b , rulek . .. . sm1 , sm2 , . . . , smn b , rulem
•
•
•
•
•
The k-th row (k = 1, 2, . . . , m) represents the scheduling policies of the jobs to be processed on machine Mk in the above chromosome. The chromosome is composed of two sub-chromosomes: the sub-chromosome sk1 , sk2 , . . . , skn b (k = 1, 2, . . . , m) represents the priority list for sequencing the bottleneck jobs that are to be processed on machine Mk (k = 1, 2, . . . , m), and the sub-chromosome rulek represents the dispatching rule for sequencing the other jobs that are to be processed on machine Mk (k = 1, 2, . . . , m). The operations in the current sub-problem corresponding to the chromosome are obtained by the decoding process, and the bottleneck jobs are determined by the corresponding GSCPM. Decoding The decoding process of each chromosome constitutes of two phases: (1) simulation of the sub-problem in the current time window, and (2) prediction of the global scheduling objective value. In phase 1, we adopt the scheduling policies based on the priority lists and the rules in the chromosome to execute the simulation of the current sub-problem, in which the following principles are adopted to compare the processing priority between two jobs: for each machine Mk (k = 1, 2, . . . , m), the processing priorities of the bottleneck jobs are always higher than those of the other jobs, while the processing priorities among the bottleneck jobs are decided by their priority indices in the corresponding priority list in the chromosome, and the processing priorities among the other jobs are decided by the corresponding dispatching rule in the chromosome. Also, it is not permitted for any machine to be idle if there exists any job that can be processed on it at each scheduling point. The simulation of the current sub-problem starts at the beginning of the current time window and it stops when the number of scheduled operations reaches n sub . In phase 2, based on the simulation results in phase 1, the formula (7) for calculating each unfinished job’s predicted completion time and the completion time of each finished job, the predicted value of the global scheduling objective (the total tardiness) can be obtained. The current sub-problems corresponding to the chromosomes have the same number of operations, but the operations in the current sub-problem may be different; these are determined by the above decoding process. Population initialization We use some good dispatching rules to generate the initial priority list of each machine in the chromosomes, and each rulek in the chromosome is selected randomly from the rule set (EDD/MDD/SLK [10]). Evaluation and selection Through the above decoding process, we can obtain the predicted value of the global scheduling objective, which is used to evaluate the corresponding chromosome. Also, we adopt the roulette-wheel method to execute the selection operation. Crossover and mutation We adopt the PMX method [29] to execute the crossover operation between two matched chromosomes, but only the crossover operation between the two sub-chromosomes corresponding to the priority lists is executed, while the crossover operation between the two sub-chromosomes corresponding to the dispatching rules is not executed. For each chromosome, we use a trivial multi-point mutation method to execute the mutation operation of each sub-chromosome, respectively. Adaptive adjustment of bottleneck ratio τ An adaptive adjusting procedure for bottleneck ratio τ is proposed to improve the optimization efficiency of the current sub-problem (the q-th sub-problem), which is detailed as follows: q
(1) Calculate the total number of predicted tardy jobs n T by X q nT = Ui q
Ji ∈Ju
where U i =
n
1 if Ci − di > 0 0 else
unfinished job Ji (Ji ∈
q Ju )
denotes whether job Ji is tardy or not, and Ci is the predicted completion time of in the decoding process.
418
M. Liu et al. / Mathematical and Computer Modelling 47 (2008) 411–421
(2) Calculate 1τ by 1 q min τmax , n T /n q − τ k 1τ = P where 1τ is the increment of τ between two neighboring iterations, τmax is the maximal value of bottleneck ratio τ , τ k is the bottleneck ratio used in the k-th iteration of the genetic algorithm, and P (P > 1) is the parameter for adjusting the bottleneck ratio τ . (3) Update τ by τ k+1 = τ k + 1τ. By the above procedure, the bottleneck ratio τ can be adjusted gradually to approximate the ratio of jobs that are probably to be tardy in order that the scheduling performance can be improved. 4. Numerical computational results and analysis 4.1. Testing problem generation and PIDA parameters We randomly generate the job shop scheduling problem instances with different scale to test the performance of the PIDA. In the generated problem instances, the technological route of each job is a random sequence of m machines, the processing time of each operation ranges uniformly between [1,10], and the due date of job Ji is calculated by Pm di = r (1, 1 + n/20) × β × p , where k=1 i,k Pn is the total number of jobs, r (1, 1 + n/20) is a uniformly random integer number ranges between [1, 1 + n/20], m k=1 pi,k is the sum of processing times of job Ji , and β is the due date tightness. Also, we set β = 1.5 in the generated problem instances. Parameters for constructing the GSCPM and the adaptive genetic algorithm include: (1) fuzziness index b : 2; (2) maximal value of bottleneck ratio τmax : 0.2; (3) Number of generation: 20; (4) Population: 20; (5) Crossover probability: 0.2; (6) Mutation probability: 0.05. Other parameters (including the number of sub-problems and the number of job clusters) are set according to problem scale. Since the above two parameters have remarkable influence on the performance of the PIDA, we have carried out some special experiments to analyze the influence of the above algorithm parameters; these are shown in Section 4.2. 4.2. Numerical computations and analysis We compare the proposed PIDA with some of the most efficient dispatching rules (EDD/SLK/MDD) and the RHA [27] on randomly generated job shop scheduling problem instances with different scale. The results are shown in Table 1, in which we list the average results (the total tardiness and the running time) over five independent runs of each scale. From Table 1, we can see that for job shop scheduling problem instances whose scales are relatively small (the number of jobs is less than 50 and the number of machines is less than 10), the performance of the PIDA is comparative with that of the RHA while the running time of the PIDA is much less than that of the RHA. Besides that, since the subproblem in the RHA is solved by a partial enumeration heuristic whose computational complexity grows exponentially with the number of operations assigned to any of the machines in a given sub-problem, the RHA costs much more time than the PIDA does when the scale of the problem instance is larger than 100 × 20. We compare the PIDA with only dispatching rules on these problem instances with larger scale in this paper. From Table 1 we also can see that the PIDA outperforms the above dispatching rules within acceptable time. Also, in Fig. 1 (take the 300 × 20 instance as an example), we list the best global scheduling objective value predicted by the GSCPM corresponding to each sub-problem in order to explain the prediction effectiveness of the global scheduling objective value. The algorithm parameters are: (1) the number of job clusters: c = 18; (2) the number of sub-problems: Q = 15. We can see from Fig. 1 that the predicted value of the global scheduling objective (the total tardiness) is comparatively stable and decreases gradually with the increase of the sub-problem index, which shows that the proposed GSCPM has higher accuracy in predicting the global scheduling objective value and can reflect the global problem characteristics efficiently. Therefore, the optimization process of each sub-problem can be effectively guided by the predicted value of the global scheduling objective.
419
M. Liu et al. / Mathematical and Computer Modelling 47 (2008) 411–421 Table 1 Performance of the PIDA on problem instances with different scale n×m
10 × 10 30 × 10 50 × 10 100 × 10 100 × 20 300 × 20 500 × 20 800 × 20 1000 × 20
EDD/SLK/MDDa P Ti Time (s)
RHAb P Ti
Time (s)
PIDA P Ti
Time (s)
(Q, c)c
86.4 749.7 2 690.8 9 295.7 5 690.3 7 761.4 9 793.2 16 982.5 25 692.1
67.9 654.0 2292.7 8207.3 4577.5 M M M M
14.0 82.3 492.9 1384.2 3530.0 M M M M
74.1 650.2 2 308.4 7 927.2 4 287.9 6 892.5 8 877.0 14 327.9 21 892.4
2.0 5.2 11.1 38.3 214.4 467.5 951.8 1901.6 3210.4
(1, 3) (2, 4) (4, 7) (6, 11) (9, 12) (15, 18) (16, 26) (21, 37) (25, 51)
0.1 0.4 0.6 4.0 5.1 8.7 18.6 47 97.1
M: Since the RHA requires too much running time (more than 1 h) when the problem scale is larger than 100 × 20, we have not made experiments on these problemPinstances. a The column Ti is the minimal total tardiness obtained by EDD/SLK/MDD, respectively, and the column “Time” is the sum of running time of the above the dispatching rules. b The number of time windows and the overlapping factor used in the RHA are set according to Macros Singer [27]. c The column (Q, c) lists the best parameter combinations of the number of sub-problems and the number of clusters for solving each problem instance. Table 2 Performance of the PIDA with different Q on problem instances with smaller scale n×m
10 × 10 20 × 10 30 × 10 50 × 10
Q=1 P Ti 76a 175a 662 2451
Time (s)
Q=2 P Ti
Time (s)
Q=3 P Ti
Time (s)
Q=4 P Ti
Time (s)
Q=5 P Ti
Time (s)
2 3 3 4
81 176 630a 2408
3 4 5 6
84 189 645 2385
3 6 7 9
92 190 677 2378a
4 6 7 12
103 218 711 2463
5 8 9 14
a Minimal total tardiness of PIDA with different Q on the given problem instance.
Fig. 1. Predicted total tardiness for the 300 × 20 instance.
In order to illustrate the influence of different parameters (Q and c) on the performance of the PIDA, we have made a lot of numerical computations. Tables 2–4 show the performance of the PIDA with different Q and c on job shop problem instances, respectively. Tables 2 and 3 show that the optimal number of sub-problems grows with the problem scale; that is, for the problem instances with smaller scale, two or three sub-problems are suitable, while for the problem instances with larger scale, better global scheduling performance might be obtained when the number of sub-problems is larger. From Table 4 we can see that the obtained total tardiness does not decrease monotonically with the increase of the number of job clusters, which indicates that the constructed GSCPM cannot perform well when the number of job clusters is either too large or too small. Therefore, we must set c in a suitable range in order to obtain better scheduling performance.
420
M. Liu et al. / Mathematical and Computer Modelling 47 (2008) 411–421
Table 3 Performance of the PIDA with different Q on problem instances with larger scale n×m
100 × 10 150 × 10 100 × 20 200 × 10 200 × 20
Q=3 P Ti 8 741 9 723 5 319 11 892 5 793
Time (s) 22 27 56 62 92
Q=5 P Ti 8 403a 9 535 5 110 11 694 5 428
Time (s)
Q=8 P Ti
Time (s)
Q = 10 P Ti
Time (s)
Q = 13 P Ti
Time (s)
31 41 91 89 149
8 592 9 127a 4 911a 10 865 5 461
48 77 159 164 196
8 609 9 154 5 021 10 060a 5 187a
60 112 178 199 251
9 126 9 262 5 220 11 215 5 304
72 136 213 224 332
a Minimal total tardiness of PIDA with different Q on the given problem instance.
Table 4 Performance of the PIDA with different c on problem instances with different scale n×m
100 × 10 150 × 10 200 × 10 300 × 10
c=5 P Ti 8 438 9 601 11 092 15 493
Time (s) 30 75 179 253
c = 10 P Ti 8 258a 9 487 10 847 15 482
Time (s)
c = 15 P Ti
Time (s)
c = 20 P Ti
Time (s)
c = 25 P Ti
Time (s)
32 77 179 255
8 316 9 139a 10 608 14 783
34 78 180 258
8 409 9 267 10 202a 13 392a
35 79 181 258
8 493 9 368 10 512 14 582
35 80 184 257
a Minimal total tardiness of PIDA with different c on the given problem instance.
5. Conclusions In this paper, a prediction based iterative decomposition algorithm for solving large-scale job shop scheduling problem is proposed. In this algorithm, the original large-scale problem is decomposed into several sub-problems, and each sub-problem is iteratively defined and solved with the rolling horizon scheme and the prediction mechanism. Numerical computational results show that PIDA is effective for large-scale scheduling problems. The results of this paper are promising for optimizing large-scale due-date-related objective based scheduling problems in practical manufacturing environments. Acknowledgments This paper is supported by the National 973 Program of China (No. 2002CB312202), the National 863 High-Tech Program of China (No. 2006AA04Z163), the National Natural Science Foundation of China (No. 60004010, No. 60274045 and No. 60443009), the Beijing Science and Technology Planning Program of China (No. D0305005040321) and Program for New Century Excellent Talents in University. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10]
K.R. Baker, Introduction to Sequencing and Scheduling, Wiley, New York, 1974. J.K. Lenstra, K.A. Rinnooy, P. Brucker, Complexity of machine scheduling problems, Annals of Discrete Mathematics 1 (1977) 343–362. M.R. Garey, D.S. Johnson, Computers and Intractability, Freeman, San Francisco, CA, 1979. S.B. Akers, A graphical approach to production scheduling problems, Operations Research 4 (2) (1956) 244–245. J.R. Jackson, An extension of Johnson’s result on job lot scheduling, Naval Research Logistics Quarterly 3 (3) (1956) 201–203. J. Carlier, E. Pinson, An algorithm for solving the job-shop problem, Management Science 35 (2) (1989) 164–176. J.R. Barker, G.B. McMahon, Scheduling the general job-shop, Management Science 31 (5) (1985) 594–598. H.M. Wagner, An integer programming model for machine scheduling, Naval Research Logistics Quarterly 6 (1959) 131–140. S.S. Panwalker, W.A. Iskander, A survey of scheduling rules, Operations Research 25 (1) (1977) 45–61. Ari P.J. Vepsalainen, Thomas E. Morton, Priority rules for job shops with total weighted tardiness costs, Management Science 33 (8) (1987) 1035–1047. [11] R. Cheng, M. Gen, A tutorial survey of job-shop scheduling problems using genetic algorithms part I: Representation, Computers & Industrial Engineering 34 (4) (1996) 983–997. [12] L. Davis, Job-shop scheduling with genetic algorithm, in: J.J. Grefenstette (Ed.), Proceedings of the 1st International Conference on Genetic Algorithms and their Applications, Lawrence Erlbaum, Pittsburgh, PA, USA, 1985, pp. 136–140.
M. Liu et al. / Mathematical and Computer Modelling 47 (2008) 411–421 [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29]
421
F. Della Croce, R. Tadei, G. Volta, A genetic algorithm for the job shop problem, Computers & Operations Research 22 (1) (1995) 15–24. J.W. Barnes, J.B. Chambers, Solving the job shop scheduling problem with tabu search, IIE Transactions 27 (22) (1995) 257–263. M. Dell’Amico, M. Trubian, Applying Tabu search to the job-shop scheduling problem, Annals of Operations Research 41 (1993) 231–252. E. Nowicki, C. Smutnicki, A fast taboo search algorithm for the job-shop problem, Management Science 42 (6) (1996) 797–813. P. Van Laarhoven, E. Aarts, J. Lenstra, Job shop scheduling by simulated annealing, Operations Research 40 (1) (1992) 113–125. T. Yamada, R. Nakano, Job-shop scheduling by simulated annealing combined with deterministic local search, in: Meta-heuristics: Theory and Applications, Kluwer Academic Publishers, MA, USA, 1996, pp. 237–248. H. Matsuo, C.J. Suh, R.S. Sullivan, A controlled search simulated annealing method for the general job shop scheduling problem, Working Paper, 1988, 03-04-88, Graduate School of Business, University of Texas at Austin, Austin, TX, USA. C. Chu, M.C. Portman, J.M. Proth, A splitting-up approach to simplify job shop scheduling problem, International Journal of Production Research 30 (4) (1992) 859–870. J. Adams, E. Balas, D. Zawack, The shifting bottleneck procedure for job shop scheduling, Management Science 34 (3) (1988) 391–401. B.L. Peter, J.H. Debra, Scheduling of manufacturing systems using the lagrangian relaxation technique, IEEE Transactions on Automatic Control 38 (7) (1993) 1066–1080. D. Sun, R. Batta, Scheduling large job shops: A decomposition approach, International Journal of Production Research 34 (7) (1996) 2019–2033. B. Eui-Seok, S.D. Wu, H.S. Robert, Decomposition heuristics for robust job-shop scheduling, IEEE Transactions on Robotics and Automation 14 (2) (1998) 303–313. D.K. Jeffrey, L. John, Mann flowsheet decomposition heuristic for scheduling: A relax-and-fix method, Computers and Chemical Engineering 28 (11) (2004) 2193–2200. G.B. Anis, H. Mohamed, An approximate decomposition algorithm for scheduling on parallel machines with heads and tails, Computers & Operations Research 34 (3) (2007) 868–883. S. Macros, Decomposition methods for large job shops, Computers & Operations Research 28 (3) (2001) 193–207. M. Liu, C. Wu, Genetic algorithm using sequence rule chain for multi-objective optimization in re-entrant micro-electronic production line, Robotics and Computer-Integrated Manufacturing 20 (3) (2004) 225–236. D.E. Goldberg, R. Lingle, Alleles, loci and the traveling salesman problem, in: Proceedings of the International Conference on Genetic Algorithms and their Applications, Lawrence Erlbaum Associates, 1985, pp. 154–159.