# An Improved Discrete Artificial Bee Colony Algorithm for Flexible Flowshop Scheduling with Step Deteriorating Jobs and Sequence-Dependent Setup Times

Authors: Hua Xuan; Huixian Zhang; Bing Li
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

This paper studies a flexible flowshop scheduling problem with step-deteriorating jobs and sequence-dependent setup times (FFSP-SDJ&SDST) where there are multiple unrelated parallel machines at each stage. The actual processing time of each job is modeled as a step function of its starting time. An integer programming model is first formulated with the objective of minimizing the total weighted completion time. Since this problem is NP-complete, it becomes an interesting and challenging topic to develop effective approximation algorithms for solving it. The artificial bee colony (ABC) algorithm has been successfully applied to solve both continuous and combinatorial optimization problems with the advantages of fewer control parameters and ease of implementation. So, an improved discrete artificial bee colony algorithm is proposed. In this algorithm, a dynamic generation mechanism of initial solutions is designed based on job permutation encoding. A genetic algorithm and a modified variable neighborhood search are introduced, respectively, to obtain new solutions for the employed and onlooker bees. A greedy heuristic is proposed to generate the solutions of the scout bees. Finally, to verify the performance of the proposed algorithm, an orthogonal test is performed to optimize the parameter settings. Simulation results on different scale problems demonstrate that the proposed algorithm is more effective compared against several presented algorithms from the existing literatures.


---

## Body

## 1. Introduction

Scheduling problem is one of key issues in many manufacturing systems. Flexible flowshop scheduling problem (FFSP) widely arises in discrete industry and process industry such as semiconductors, electronics manufacturing, cosmetics, and pharmaceuticals [1]. In the classical FFSP, the processing times of jobs are always assumed to be known and fixed. However, it might not be true in realistic industrial processes such as steel production and equipment maintenance [2]. The actual processing time may become longer due to equipment wear, physical properties of jobs, fatigue of the workers, and so on. This phenomenon is known as scheduling with deteriorating jobs. In addition, it often happens that some tasks such as equipment cleaning and tool replacement are performed when a machine handles different jobs. That is, a setup time occurs between adjacent jobs on the same machine and it always depends on the job sequence [3]. Therefore, this paper considers FFSP with step-deteriorating jobs and sequence-dependent setup times where step deterioration means that if the starting time of a job is more than a prespecified deteriorating date, it will require extra penalty time to complete its tasks successfully. The scheduling goal is to determine the assignment of jobs to machines and the sequence of jobs on each machine.Scheduling with deterioration has been presented and defined first by Gupta and Gupta [4] and Browne and Yechiali [5]. The current research studies on deteriorating job scheduling problems mostly assume that the job processing times are linear functions of their starting times. For single-machine scheduling, Wang and Liu [6] considered a group scheduling problem with setup times. A modified Lawler’s algorithm was proposed to minimize the total weighted completion time and the maximum cost simultaneously. In terms of unrelated parallel machine scheduling, Wang and Li [7] considered modifying maintenance activities with the objective of minimizing a cost function of total completion (waiting) time, total absolute differences in completion (waiting) time, and total machine load. They proved that the problem can be solved in polynomial time. For flowshop scheduling, Wang et al. [8] presented a multiverse optimizer algorithm so that the total tardy time is minimized. There is a relatively little research on step deterioration. To solve identical parallel machine scheduling, Guo et al. [2] considered sequence-dependent setup times and applied a hybrid discrete cuckoo algorithm to minimize total tardiness; Cheng et al. [9] presented a modified weight-combination search algorithm and a variable neighborhood search to minimize total completion time. For two-stage FFSP with identical parallel machines at stage 1 and a single batching machine at stage 2, Gong and Tang [10] proposed a heuristic algorithm to minimize makespan plus the total penalty cost of batching-machine utilization ratio where the job processing time at stage 2 is a step function of its waiting time between the two stages.No setup is commonly considered in many shop scheduling problems [11, 12]. However, separated setup times have to be included in scheduling decision so as to describe realistic production process more accurately. Recently, the sequence-dependent setup times (SDST) have attracted increasing attention among researchers. For the flowshop with SDST, Wang et al. [13] presented a hybrid local search algorithm to minimize makespan; Wang and Li [14] applied a hybrid chaotic biogeography-based optimization to minimize the total weighted tardiness; Li and Li [15] developed multiobjective local search based decomposition to solve permutation flowshop with the objective of minimizing makespan and total flowtime. Regarding the multistage FFSP with SDST, Pedro et al. [16] presented an improved genetic algorithm to minimize makespan. Marichelvam et al. [17] used a discrete firefly algorithm to minimize total tardiness. Pan et al. [18] developed nine metaheuristics based on trajectory and population to minimize makespan. Santosa and Riyanto [19] applied hybrid differential evolution-variable neighborhood search to minimize makespan and maximum lateness simultaneously.In past decades, many efforts have been made on heuristic and metaheuristic algorithms to obtain near-optimal solutions with high quality in a reasonable computational time [20, 21]. The artificial bee colony (ABC) algorithm is one of the population-based evolutionary metaheuristics, proposed by Karaboga [22], which simulates the intelligent foraging behavior of honey bees in nature. It has been applied to solve some shop scheduling problems. Chow et al. [23] used this algorithm to solve the flowshop makespan problem. Li and Ma [24] presented a novel discrete ABC (DABC) algorithm for the permutation flowshop scheduling with sequence-dependent setup times in order to minimize makespan and total flowtime. Peng et al. [25] improved the ABC algorithm to solve FFSP with serial-batching machines on the last stage in which the objective is minimizing sojourn time, earliness of casting starting penalty, and tardiness of cast starting penalty. Cui and Gu [26] proposed an improved DABC algorithm for FFSP to minimize makespan. Li et al. [27] used an improved DABC algorithm for hybrid FFSP with dynamic operation skipping so that the average sojourn time, the penalties of tardiness and earliness, and the penalty of the skipping rate are minimized.To the best of our knowledge, most research studies on scheduling with deterioration concentrate on linear deteriorating jobs and consider the machine environments such as single machine, parallel machine, flowshop, and two-stage FFSP. Little work has been carried out on a multistage FFSP with step-deteriorating jobs. This problem is NP-complete even if there is only one single machine [28]. In this paper, the more complicated multistage FFSP with step-deteriorating jobs and sequence-dependent setup times (FFSP-SDJ&SDST) is studied to minimize the sum of weighted completion times. A new improved DABC (IDABC) algorithm is then developed to obtain near-optimal solutions.The rest of this paper is organized as follows. Section2 describes the problem and formulates a mathematical model. Section 3 develops an IDABC algorithm. In Section 4, computational experiments are performed and the algorithm comparisons are provided. Finally, the conclusions are summarized in Section 5.


## 2. Problem Formulation

### 2.1. Problem Description

In the FFSP-SDJ&SDST, a set ofn jobs {J1, J2, … , Jn} has to be processed on S stages, and each stage k has mk unrelated machines (mk ≥ 2 for at least one stage). These unrelated machines are continuously available from time zero. Each job should be started at its release time ri (ri > 0). Job j can be handled on any one of the mk machines at stage k, and the basic processing time is different on each machine. A setup time δfijk occurs when job f is handled immediately before job i on machine j at stage k. δ0ijk = 0 means that no setup time is required if job i is the first job processed on machine j at stage k. At stage k, job i has a basic constant processing time Pijk on machine j if its starting time Stijk is no more than its deteriorating date hijk; otherwise, it requires an extra penalty time bijk. Transfer times between two successive stages have to be considered because job delivery is often carried out by using conveyors, carts, and so on. The goal of the problem is to minimize the total weighted completion time.


### 2.2. Problem Definition and Notation

To describe the above problem, the following notations are introduced:(1)
Parametersn: number of jobsS: number of stagesmk: number of machines on stage kTk,k+1: delivery time of jobs from stage k to stage k + 1ri: release time of job iwi: weight of job iδfijk: setup time when job f is processed immediately before job i on machine j at stage khijk: deteriorating date of job i on machine j at stage kbijk: deteriorating penalty time of job i on machine j at stage kPijk: basic processing time of job i on machine j at stage k(2)
Decision variablesCik: completion time of job i at stage kStijk: starting time of job i on machine j at stage kPijk′: actual processing time of job i on machine j at stage kXijk: binary, set to 1 if job i is assigned to machine j at stage k; 0 otherwiseYfijk: binary, set to 1 if job f is handled immediately before job i on machine j at stage k; 0 otherwise


### 2.3. Model Formulation

Applying the above notations, the studied problem can be formulated as an integer programming model, as shown below:(1)MinO=∑i=1nwiCiS,subject to(2)∑j=1mkXijk=1,i=1,…,n,k=1,…,S,(3)∑i=1n∑j=1mkY0ijk=mk,k=1,…,S,(4)Pijk′=Pijk,Stijk≤hijk,Pijk+bijk,Stijk>hijk,(5)Cik≥Stijk+Pijk′,i=1,…,n,j=1,…,mk,k=1,…,S,(6)Ci1−∑j=1m1Xij1Pij1′≥ri,i=1,…,n,(7)∑f=1n∑i=1nYfijkStijk≥∑f=1n∑i=1nYfijkCfk+δfijk,f≠i,j=1,…,mk,k=1,…,S,(8)Ci,k+1−Ci,k≥Tk,k+1+∑j=1mk+1Xij,k+1Pij,k+1′,i=1,…,n,k=1,…,S,(9)Xijk,Yfijk∈0,1,f=0,…,n−1,i=1,…,n,j=1,…,mk,k=1,…,S.The objective function (1) indicates the minimum total weighted completion time. Constraints (2) ensure that each job has to be processed in all stages and at each stage, a job can select only one available machine. Constraints (3) guarantee that all mk machines are scheduled at each stage when n ≥ max {m1, … , mS}. Constraints (4) define the actual processing time of a job in consideration of step-deteriorating effect. Constraints (5) indicate the relation of the completion time and beginning time of a job. Constraints (6) make sure that the beginning time of a job at stage 1 must be larger than or equal to its release time. Constraints (7) ensure that, for two consecutive jobs on the same machine, the next one can be started only after the completion of the preceding one plus the setup time between the two jobs. Constraints (8) indicate that, for two adjacent operations of a job, the latter one can be processed only after the former one has been finished and transferred to the corresponding machine. Constraints (9) represent the value ranges of decision variables.


## 2.1. Problem Description

In the FFSP-SDJ&SDST, a set ofn jobs {J1, J2, … , Jn} has to be processed on S stages, and each stage k has mk unrelated machines (mk ≥ 2 for at least one stage). These unrelated machines are continuously available from time zero. Each job should be started at its release time ri (ri > 0). Job j can be handled on any one of the mk machines at stage k, and the basic processing time is different on each machine. A setup time δfijk occurs when job f is handled immediately before job i on machine j at stage k. δ0ijk = 0 means that no setup time is required if job i is the first job processed on machine j at stage k. At stage k, job i has a basic constant processing time Pijk on machine j if its starting time Stijk is no more than its deteriorating date hijk; otherwise, it requires an extra penalty time bijk. Transfer times between two successive stages have to be considered because job delivery is often carried out by using conveyors, carts, and so on. The goal of the problem is to minimize the total weighted completion time.


## 2.2. Problem Definition and Notation

To describe the above problem, the following notations are introduced:(1)
Parametersn: number of jobsS: number of stagesmk: number of machines on stage kTk,k+1: delivery time of jobs from stage k to stage k + 1ri: release time of job iwi: weight of job iδfijk: setup time when job f is processed immediately before job i on machine j at stage khijk: deteriorating date of job i on machine j at stage kbijk: deteriorating penalty time of job i on machine j at stage kPijk: basic processing time of job i on machine j at stage k(2)
Decision variablesCik: completion time of job i at stage kStijk: starting time of job i on machine j at stage kPijk′: actual processing time of job i on machine j at stage kXijk: binary, set to 1 if job i is assigned to machine j at stage k; 0 otherwiseYfijk: binary, set to 1 if job f is handled immediately before job i on machine j at stage k; 0 otherwise


## 2.3. Model Formulation

Applying the above notations, the studied problem can be formulated as an integer programming model, as shown below:(1)MinO=∑i=1nwiCiS,subject to(2)∑j=1mkXijk=1,i=1,…,n,k=1,…,S,(3)∑i=1n∑j=1mkY0ijk=mk,k=1,…,S,(4)Pijk′=Pijk,Stijk≤hijk,Pijk+bijk,Stijk>hijk,(5)Cik≥Stijk+Pijk′,i=1,…,n,j=1,…,mk,k=1,…,S,(6)Ci1−∑j=1m1Xij1Pij1′≥ri,i=1,…,n,(7)∑f=1n∑i=1nYfijkStijk≥∑f=1n∑i=1nYfijkCfk+δfijk,f≠i,j=1,…,mk,k=1,…,S,(8)Ci,k+1−Ci,k≥Tk,k+1+∑j=1mk+1Xij,k+1Pij,k+1′,i=1,…,n,k=1,…,S,(9)Xijk,Yfijk∈0,1,f=0,…,n−1,i=1,…,n,j=1,…,mk,k=1,…,S.The objective function (1) indicates the minimum total weighted completion time. Constraints (2) ensure that each job has to be processed in all stages and at each stage, a job can select only one available machine. Constraints (3) guarantee that all mk machines are scheduled at each stage when n ≥ max {m1, … , mS}. Constraints (4) define the actual processing time of a job in consideration of step-deteriorating effect. Constraints (5) indicate the relation of the completion time and beginning time of a job. Constraints (6) make sure that the beginning time of a job at stage 1 must be larger than or equal to its release time. Constraints (7) ensure that, for two consecutive jobs on the same machine, the next one can be started only after the completion of the preceding one plus the setup time between the two jobs. Constraints (8) indicate that, for two adjacent operations of a job, the latter one can be processed only after the former one has been finished and transferred to the corresponding machine. Constraints (9) represent the value ranges of decision variables.


## 3. Improved Discrete Artificial Bee Colony Algorithm

The ABC algorithm is originally used to solve continuous function optimization. However, as mentioned above, the considered FFSP-SDJ&SDST is a discrete combination problem. Therefore, it is necessary to make discretization improvements so as to solve FFSP-SDJ&SDST. So, an improved discrete artificial bee colony (IDABC) algorithm is presented to solve the NP-complete problem.


### 3.1. The Framework of the Proposed IDABC Algorithm

The ABC algorithm starts with the generation of initial solutions and then performs an iteration of employed bee phase, onlooker bee phase, and scout bee phase. Since the considered problem belongs to discrete function optimization, this paper utilizes job sequence as solution representation and presents an IDABC algorithm. To make the algorithm more effective, a genetic algorithm (GA) is adopted including single-point crossover and multiple insertion mutation for employed bees; a modified variable neighborhood search (MVNS) using perturbation and search procedures is designed for onlooker bees; a greedy heuristic (GH) with destruction and reconstruction procedures is proposed for scout bees. The above modification or improvement is performed in the novel IDABC algorithm to obtain better food sources. The process of the IDABC algorithm is shown in Figure1. The detailed steps are provided as follows:Step 1: parameter setting phase. Set the following system parameters: maximum iteration numberImax, population size Nump, iteration number limit related to food sources to be abandoned, crossover probability CP, and mutation probability MPStep 2: initialization of population. Randomly generateNump job sequences (food sources) as initial solutionsStep 3: solution evaluation. Compute fitness values of food sources and memorize the best solutionStep 4: if the termination condition is met, output the optimal solution; otherwise, perform Step 5∼Step 9Step 5: employed bee phase. Utilize GA to generate the new food sourcesStep 6: calculate the fitness value of each newly generated food sources and use the roulette wheel method to select food sources for onlooker beesStep 7: onlooker bee phase. Use MVNS to search round the selected food sources. Calculate the fitness values and find the better solutionsStep 8: if there is a food source to be abandoned, perform Step 9; otherwise, go to Step 4Step 9: scout bee phase. For the food sources without improvement during thelimit iterations, a new solution is generated by using GHFigure 1
Framework of the proposed IDABC algorithm.


### 3.2. Encoding and Decoding Procedures

The FFSP-SDJ&SDST considers the sequence-dependent setup times and unrelated parallel machines. The scheduling goal consists of assigning the jobs to the appropriate machines at each stage and ordering the jobs assigned to the same machine. Thus, a dynamic generation mechanism of initial solutions is proposed where job sequence-based representation for stage 1 is applied to encode the solutions and the decoding scheme containing two steps is applied. The initial solutions of the population are then obtained.


#### 3.2.1. Encoding Solutions

To obtainNump food sources, each job sequence at stage 1 is yielded from discrete domain. So, a random perturbation procedure based on job numbers is developed to construct the initial job sequence ω1 = {ω1(1), ω1(2), … , ω1(i), … , ω1(n)}, where ωk(i) denotes the job on the ith position at stage k. Figure 2 illustrates an encoding scheme for an instance with five jobs where ω1(1) = 4.Figure 2
An instance with five jobs.


#### 3.2.2. Decoding Food Sources

Because the goal is to minimize the total weighted completion time, the job with the earlier completion time at the immediately preceding stage is allowed to have a higher priority in order to improve on-time delivery between adjacent stages. So, the job sequence at stagek(k > 1) is determined according to their completion times at stage k −1.At stage 1, according to the above generated job sequence, the jobω1(1) is firstly assigned to the machine with the shortest actual processing time. To balance the utility of all the machines, the job ω1(2) chooses one machine with the shortest actual processing time from the unscheduled machines until the job ω1(m1) is arranged. For the unscheduled n-m1 jobs, compute the possible completion time of job ω1(i) (i > m1) equal to the sum of the completion time of the immediately preceding job on the same machine, the setup time and its actual processing time, and assign it to the machine with the shortest completion time.The job sequence at stagek(k > 1) is obtained in the increasing order of the completion time of all jobs at stage k−1. If some jobs have the same completion time, order these jobs randomly. Similar with stage 1, the first mk job is, respectively, arranged on the mk machines with the shortest actual processing time based on the sequence so that all machines are used. When it turns to ωk(i) (i > mk), the starting time of job ωk(i) is set to the maximum between its completion time at stage k − 1 plus the delivery time and the completion time of the immediately preceding job on the same machine plus the setup time. So, its completion time can be computed by its starting time plus the actual processing time. The assigned machine is decided with the shortest completion time. The process continues until all jobs are arranged in all stages, as shown in Figure 3.Figure 3
Decoding scheme using the rule of minimum processing/completion time of the free machine.


#### 3.2.3. Population Initialization

As described above, an individual is denoted as a permutation of all jobs at stage 1. To generate a population, a simple random procedure is adopted whereNump individuals (food sources) are consistent with the size of employed bee or onlooker bee.To clearly illustrate the FFSP-SDJ&SDST and the decoding procedure, a simple example is introduced with 5 jobs and 3 stages where each stage has three machines. The release time and weight of jobs are given byri = [3, 1, 4, 2, 3] and wi = [5, 6, 9, 7, 3]; delivery time Tk,k+1 is, respectively, set to 3 and 2. For simplicity, it is assumed that the setup time δfijk is dependent on jobs and irrelevant to machines and denoted as δfi′. The basic processing time Pijk, deteriorating date hijk, deteriorating penalty time bijk, and setup time δfi′ are, respectively, given below:(10)The initial job sequence [4, 5, 1, 3, 2] at stage 1 is given and consistent with the instance illustrated in Figure2. The machine allocation and completion time of all jobs are computed as follows. The Gantt chart of a solution is drawn in Figure 4, where each operation is denoted by a rectangle labeled with a job number.Stage 1:ω1 = {4, 5, 1, 3, 2}C41 = r4 + min{P411, P421, P431} = 2 + min{5, 13, 11} = 7X411 = 1, X421 = X431 = 0C51 = r5 + min{P421, P431} = 3 + min{10, 19} = 13X521 = 1, X511 = X531 = 0C11 = r1 + P131 = 3 + 8 = 11X131 = 1, X111 = X121 = 0C31 = min{C41 + δ43′ + P311, C51 + δ53′ + P321, C11 + δ13′ + P331} = min{7 + 5 + 17, 13 + 1 + 6, 11 + 4 + 12} = 20 (C41 + δ43′ < h311, C51 + δ43′ < h321, C11 + δ13′ < h331)X321 = 1, X311 = X331 = 0C21 = min{C41 + δ42′ + P211, C31 + δ32′ + P221, C11 + δ12′ + P231} = {7 + 4 + 6, 20 + 2 + 4, 11 + 3 + 17} = 17 (C41 + δ42′ < h211, C31 + δ32′ < h221, C11 + δ12′ < h231)X211 = 1, X221 = X231 = 0Stage 2:ω2 = {4, 1, 5, 2, 3}C42 = C41 + T12 + min{P412, P422, P432} = 7 + 3 + min{14, 16, 12} = 22 (C41 + T12 < h412, C41 + T12 < h422, C41 + T12 < h432)X432 = 1, X412 = X422 = 0C12 = C11 + T12 + min{P112, P122} = 11 + 3 + min{16, 12} = 26 (C11 + T12 < h112, C11 + T12 < h122)X122 = 1, X112 = X132 = 0C52 = C51 + T12 + P512 = 13 + 3 + 12 = 28 (C51 + T12 < h512)X512 = 1, X522 = X532 = 0C22 = min{max{C21 + T12, C52 + δ52′} + P212, max{C21 + T12, C12 + δ12′} + P222, max{C21 + T12, C42 + δ42′} + P232} = min{max{17 + 3, 28 + 1} + 16, max{17 + 3, 26 + 3} + 10, max{17 + 3, 22 + 4} + 20 = 39 (max{C21 + T12, C52 + δ52′} < h212, max{C21 + T12, C12 + δ12′} < h222, max{C21 + T12, C42 + δ42′} < h232)X222 = 1, X212 = X232 = 1C32 = min{max{C31 + T12, C52 + δ53′} + P312 + b312, max{C31 + T12, C22 + δ23′} + P322 + b322, max{C31 + T12, C42 + δ43′} + P332 + b332} = min{max{20 + 3, 28 + 1} + 16 + 2, max{20 + 3, 39 + 2} + 11 + 6, max{20 + 3, 22 + 5} + 10 + 4} = 41 (max{C31 + T12, C52 + δ53′} > h312, max{C31 + T12, C22 + δ23′} > h322, max{C31 + T12, C42 + δ43′} > h332)X332 = 1, X312 = X322 = 0Stage 3:ω3 = {4, 1, 5, 2, 3}C43 = C42 + T23 + min{P413, P423, P433} = 22 + 2 + min{11, 14, 19} = 35 (C42 + T23 < h413, C42 + T23 < h423, C42 + T23 < h433)X413 = 1, X423 = X433 = 0C13 = C12 + T23 + min{P123, P133} = 26 + 2 + min{13, 10} = 38(C12 + T23 < h123, C12 + T23 < h133)X133 = 1, X113 = X123 = 0C53 = C52 + T23 + P523 + b523 = 28 + 2 + 15 + 5 = 50 (C52 + T23 > h523)X523 = 1, X513 = X533 = 0C23 = min{max{C22 + T23, C43 + δ42′} + P213 + b213, max{C22 + T23, C53 + δ52′} + P223 + b223, max{C22 + T23, C13 + δ12′} + P233 + b233} = min{max{39 + 2, 35 + 4} + 16 + 3, max{39 + 2, 50 + 1} + 11 + 1, max{39 + 2, 38 + 3} + 9 + 4} = 54(max{C22 + T23, C43 + δ42′} > h213, max{C22 + T23, C53 + δ52′} > h223, max{C22 + T23, C13 + δ12′} > h233)X233 = 1, X213 = X223 = 0C33 = min{max{C32 + T23, C43 + δ43′} + P313 + b313, max{C32 + T23, C53 + δ53′} + P323 + b323, max{C32 + T23, C23 + δ23′} + P333 + b333} = min{max{41 + 2, 35 + 5} + 7 + 2, max{41 + 2, 50 + 1} + 14 + 10, max{41 + 2, 54 + 2} + 14 + 10} = 52 (max{C32 + T23, C43 + δ43′} > h313, max{C32 + T23, C53 + δ53′} > h323, max{C32 + T23, C23 + δ23′} > h333)X313 = 1, X323 = X333 = 0Therefore, the total weighted completion time is:∑iwiCiS = 5 × 38 + 6 × 54 + 9 × 52 + 7 × 35 + 3 × 50 = 1377.Figure 4
Example Gantt chart.


### 3.3. Employed Bee Phase

In the employed bee phase, the new solutions are yielded by performing a genetic algorithm in order to find a neighboring food source from the current food source. The detailed GA is described below.


#### 3.3.1. Fitness Evaluation and Selection Operation

Because the objective of this model is to minimize total weighted completion time, the reciprocal of the objective function is used as the fitness function. Therefore, the fitness function of theyth food source is estimated by f(y) = 1/∑iwiCiS. The roulette wheel selection is applied to determine the offspring where the individual with the higher fitness value has the higher probability to be selected.


#### 3.3.2. Single-Point Crossover Operator

Step 1: randomly generate a real number Rc from the interval [0, 1] and compare it with the crossover probabilityCP. If Rc ≤ CP, two food sources A = {ωA(1), ωA(2), … , ωA(i), … , ωA(n)} and B = {ωB(1), ωB(2), … , ωB(i), … , ωB(n)} are picked out randomly from the population after the selection operation.
Step 2: a crossover positionθ (0 < θ<n) is generated randomly, and the single-point crossover is performed on A and B. As shown in Figure 5, the genes after position θ in A (B) are copied to A′ (B′) by keeping their position unchanged. The lacking genes in A′ (B′) are filled with the first θ genes in B (A) by maintaining their relative orders.
Step 3: find the job number (denoted asμ(θ1)) in {ωB(1), … , ωB(θ)} which emerges repeatedly in {ωA(θ + 1), … , ωA(n)} for A′. Similarly, find the job number v (θ2) from {ωA(1), … , ωA(θ)} with the same as {ωB(θ + 1), … , ωB(n) for B′.
Step 4: replaceμ(θ1) in the θ1th position with v (θ2) in A′ and v (θ2) in the θ2th position with μ(θ1) in B′. Therefore, two new food sources are formed.Figure 5
Crossover operator.


#### 3.3.3. Multi-Insertion Mutation Operator

A real numberRm is randomly produced from the interval [0, 1]. If Rm is not larger than mutation probability MP, two job positions (θ1 and θ2) are randomly chosen in a permutation Z = {ω(1), ω(2), … , ω(n)}. Suppose that IT is an insertion number. Insert job ω(θ1) into (θ2 + q − 1)th (q ∈ {1, … , IT}) position, as shown in Figure 6. If θ2 + q − 1 > n, let θ2 + q − 1 equal to 1. Thus, IT new food sources are obtained to increase the diversity of population.Figure 6
Mutation operator.


### 3.4. Onlooker Bee Phase

This phase uses the same roulette wheel selection applied in the employed bee phase for the food source returned by the employed bee. An onlooker bee chooses a food sourceω depending on its probability value FPω = f(ω)/∑yf(y).The modified variable neighborhood search (MVNS) is then applied for the above food sources. Firstly, the following neighborhood structures NSφ (φ = 1, 2, 3) are designed to help the onlooker bee explore possible promising neighbors.


#### 3.4.1. Exchange Mutation

Two randomly generated positions (θ1 and θ2) are chosen and the corresponding jobs are swapped, as shown in Figure 7.Figure 7
Exchange mutation.


#### 3.4.2. 2-Opt Invert Mutation

Two randomly selected jobs (ω(θ1) and ω(θ2)) are picked out and the subsequence between positions θ1 and θ2 are inverted, as illustrated in Figure 8.Figure 8
2-opt invert mutation.


#### 3.4.3. Fragment Insert Mutation

Randomly select two different jobsω(θ1) and ω(θ2) and remove the subsequence between positions θ1 and θ2. A randomly selected position τ from the remaining positions is chosen, and the subsequence between θ1 and θ2 is inserted into position τ, as shown in Figure 9.Figure 9
Fragment insert mutation.Next, the new neighboring solution obtained by applying these search rules is evaluated and compared to the primary food source. If the new solution obtained by MVNS is better than the primary one, the primary food source is replaced and the new food source becomes a new member in the population.The procedure of MVNS is listed in Algorithm1.Algorithm 1:Procedure of MVNS.
Initialization: define NSφφ=1,2,3
for every food source ω do
Perturbation:
x = 1;
whilex < 3;
executeNS1 to ω
x = x + 1;
end while
Search:
φ=1
while the stopping criterion is not met
σ = execute NSφ to ω
If fσ>fω
ω=σ;
φ=φ;
else
φ=φ+1;
end if
end while
end for


### 3.5. Scout Bee Phase

If a food source has no improvement duringlimit iterations, this employed bee will become a scout bee. A greedy heuristic is developed to produce a new food source including two steps.


#### 3.5.1. Destruction Procedure

Jobω(θ) is randomly chosen and removed from the food source Z = {ω(1), ω(2), … , ω(θ − 1), ω(θ), ω(θ + 1), … , ω(n)} so that a remaining food source Z − θ = {ω(1), ω(2), … , ω(θ−1), ω(θ + 1) … , ω(n)} is generated.


#### 3.5.2. Reconstruction Procedure

The selected jobω(θ) is sequentially inserted into another randomly selected position of Z-θ, and the corresponding fitness value is calculated. Compare these new generated neighboring solutions and choose the best one with the largest fitness value as the current solution.


## 3.1. The Framework of the Proposed IDABC Algorithm

The ABC algorithm starts with the generation of initial solutions and then performs an iteration of employed bee phase, onlooker bee phase, and scout bee phase. Since the considered problem belongs to discrete function optimization, this paper utilizes job sequence as solution representation and presents an IDABC algorithm. To make the algorithm more effective, a genetic algorithm (GA) is adopted including single-point crossover and multiple insertion mutation for employed bees; a modified variable neighborhood search (MVNS) using perturbation and search procedures is designed for onlooker bees; a greedy heuristic (GH) with destruction and reconstruction procedures is proposed for scout bees. The above modification or improvement is performed in the novel IDABC algorithm to obtain better food sources. The process of the IDABC algorithm is shown in Figure1. The detailed steps are provided as follows:Step 1: parameter setting phase. Set the following system parameters: maximum iteration numberImax, population size Nump, iteration number limit related to food sources to be abandoned, crossover probability CP, and mutation probability MPStep 2: initialization of population. Randomly generateNump job sequences (food sources) as initial solutionsStep 3: solution evaluation. Compute fitness values of food sources and memorize the best solutionStep 4: if the termination condition is met, output the optimal solution; otherwise, perform Step 5∼Step 9Step 5: employed bee phase. Utilize GA to generate the new food sourcesStep 6: calculate the fitness value of each newly generated food sources and use the roulette wheel method to select food sources for onlooker beesStep 7: onlooker bee phase. Use MVNS to search round the selected food sources. Calculate the fitness values and find the better solutionsStep 8: if there is a food source to be abandoned, perform Step 9; otherwise, go to Step 4Step 9: scout bee phase. For the food sources without improvement during thelimit iterations, a new solution is generated by using GHFigure 1
Framework of the proposed IDABC algorithm.


## 3.2. Encoding and Decoding Procedures

The FFSP-SDJ&SDST considers the sequence-dependent setup times and unrelated parallel machines. The scheduling goal consists of assigning the jobs to the appropriate machines at each stage and ordering the jobs assigned to the same machine. Thus, a dynamic generation mechanism of initial solutions is proposed where job sequence-based representation for stage 1 is applied to encode the solutions and the decoding scheme containing two steps is applied. The initial solutions of the population are then obtained.


### 3.2.1. Encoding Solutions

To obtainNump food sources, each job sequence at stage 1 is yielded from discrete domain. So, a random perturbation procedure based on job numbers is developed to construct the initial job sequence ω1 = {ω1(1), ω1(2), … , ω1(i), … , ω1(n)}, where ωk(i) denotes the job on the ith position at stage k. Figure 2 illustrates an encoding scheme for an instance with five jobs where ω1(1) = 4.Figure 2
An instance with five jobs.


### 3.2.2. Decoding Food Sources

Because the goal is to minimize the total weighted completion time, the job with the earlier completion time at the immediately preceding stage is allowed to have a higher priority in order to improve on-time delivery between adjacent stages. So, the job sequence at stagek(k > 1) is determined according to their completion times at stage k −1.At stage 1, according to the above generated job sequence, the jobω1(1) is firstly assigned to the machine with the shortest actual processing time. To balance the utility of all the machines, the job ω1(2) chooses one machine with the shortest actual processing time from the unscheduled machines until the job ω1(m1) is arranged. For the unscheduled n-m1 jobs, compute the possible completion time of job ω1(i) (i > m1) equal to the sum of the completion time of the immediately preceding job on the same machine, the setup time and its actual processing time, and assign it to the machine with the shortest completion time.The job sequence at stagek(k > 1) is obtained in the increasing order of the completion time of all jobs at stage k−1. If some jobs have the same completion time, order these jobs randomly. Similar with stage 1, the first mk job is, respectively, arranged on the mk machines with the shortest actual processing time based on the sequence so that all machines are used. When it turns to ωk(i) (i > mk), the starting time of job ωk(i) is set to the maximum between its completion time at stage k − 1 plus the delivery time and the completion time of the immediately preceding job on the same machine plus the setup time. So, its completion time can be computed by its starting time plus the actual processing time. The assigned machine is decided with the shortest completion time. The process continues until all jobs are arranged in all stages, as shown in Figure 3.Figure 3
Decoding scheme using the rule of minimum processing/completion time of the free machine.


### 3.2.3. Population Initialization

As described above, an individual is denoted as a permutation of all jobs at stage 1. To generate a population, a simple random procedure is adopted whereNump individuals (food sources) are consistent with the size of employed bee or onlooker bee.To clearly illustrate the FFSP-SDJ&SDST and the decoding procedure, a simple example is introduced with 5 jobs and 3 stages where each stage has three machines. The release time and weight of jobs are given byri = [3, 1, 4, 2, 3] and wi = [5, 6, 9, 7, 3]; delivery time Tk,k+1 is, respectively, set to 3 and 2. For simplicity, it is assumed that the setup time δfijk is dependent on jobs and irrelevant to machines and denoted as δfi′. The basic processing time Pijk, deteriorating date hijk, deteriorating penalty time bijk, and setup time δfi′ are, respectively, given below:(10)The initial job sequence [4, 5, 1, 3, 2] at stage 1 is given and consistent with the instance illustrated in Figure2. The machine allocation and completion time of all jobs are computed as follows. The Gantt chart of a solution is drawn in Figure 4, where each operation is denoted by a rectangle labeled with a job number.Stage 1:ω1 = {4, 5, 1, 3, 2}C41 = r4 + min{P411, P421, P431} = 2 + min{5, 13, 11} = 7X411 = 1, X421 = X431 = 0C51 = r5 + min{P421, P431} = 3 + min{10, 19} = 13X521 = 1, X511 = X531 = 0C11 = r1 + P131 = 3 + 8 = 11X131 = 1, X111 = X121 = 0C31 = min{C41 + δ43′ + P311, C51 + δ53′ + P321, C11 + δ13′ + P331} = min{7 + 5 + 17, 13 + 1 + 6, 11 + 4 + 12} = 20 (C41 + δ43′ < h311, C51 + δ43′ < h321, C11 + δ13′ < h331)X321 = 1, X311 = X331 = 0C21 = min{C41 + δ42′ + P211, C31 + δ32′ + P221, C11 + δ12′ + P231} = {7 + 4 + 6, 20 + 2 + 4, 11 + 3 + 17} = 17 (C41 + δ42′ < h211, C31 + δ32′ < h221, C11 + δ12′ < h231)X211 = 1, X221 = X231 = 0Stage 2:ω2 = {4, 1, 5, 2, 3}C42 = C41 + T12 + min{P412, P422, P432} = 7 + 3 + min{14, 16, 12} = 22 (C41 + T12 < h412, C41 + T12 < h422, C41 + T12 < h432)X432 = 1, X412 = X422 = 0C12 = C11 + T12 + min{P112, P122} = 11 + 3 + min{16, 12} = 26 (C11 + T12 < h112, C11 + T12 < h122)X122 = 1, X112 = X132 = 0C52 = C51 + T12 + P512 = 13 + 3 + 12 = 28 (C51 + T12 < h512)X512 = 1, X522 = X532 = 0C22 = min{max{C21 + T12, C52 + δ52′} + P212, max{C21 + T12, C12 + δ12′} + P222, max{C21 + T12, C42 + δ42′} + P232} = min{max{17 + 3, 28 + 1} + 16, max{17 + 3, 26 + 3} + 10, max{17 + 3, 22 + 4} + 20 = 39 (max{C21 + T12, C52 + δ52′} < h212, max{C21 + T12, C12 + δ12′} < h222, max{C21 + T12, C42 + δ42′} < h232)X222 = 1, X212 = X232 = 1C32 = min{max{C31 + T12, C52 + δ53′} + P312 + b312, max{C31 + T12, C22 + δ23′} + P322 + b322, max{C31 + T12, C42 + δ43′} + P332 + b332} = min{max{20 + 3, 28 + 1} + 16 + 2, max{20 + 3, 39 + 2} + 11 + 6, max{20 + 3, 22 + 5} + 10 + 4} = 41 (max{C31 + T12, C52 + δ53′} > h312, max{C31 + T12, C22 + δ23′} > h322, max{C31 + T12, C42 + δ43′} > h332)X332 = 1, X312 = X322 = 0Stage 3:ω3 = {4, 1, 5, 2, 3}C43 = C42 + T23 + min{P413, P423, P433} = 22 + 2 + min{11, 14, 19} = 35 (C42 + T23 < h413, C42 + T23 < h423, C42 + T23 < h433)X413 = 1, X423 = X433 = 0C13 = C12 + T23 + min{P123, P133} = 26 + 2 + min{13, 10} = 38(C12 + T23 < h123, C12 + T23 < h133)X133 = 1, X113 = X123 = 0C53 = C52 + T23 + P523 + b523 = 28 + 2 + 15 + 5 = 50 (C52 + T23 > h523)X523 = 1, X513 = X533 = 0C23 = min{max{C22 + T23, C43 + δ42′} + P213 + b213, max{C22 + T23, C53 + δ52′} + P223 + b223, max{C22 + T23, C13 + δ12′} + P233 + b233} = min{max{39 + 2, 35 + 4} + 16 + 3, max{39 + 2, 50 + 1} + 11 + 1, max{39 + 2, 38 + 3} + 9 + 4} = 54(max{C22 + T23, C43 + δ42′} > h213, max{C22 + T23, C53 + δ52′} > h223, max{C22 + T23, C13 + δ12′} > h233)X233 = 1, X213 = X223 = 0C33 = min{max{C32 + T23, C43 + δ43′} + P313 + b313, max{C32 + T23, C53 + δ53′} + P323 + b323, max{C32 + T23, C23 + δ23′} + P333 + b333} = min{max{41 + 2, 35 + 5} + 7 + 2, max{41 + 2, 50 + 1} + 14 + 10, max{41 + 2, 54 + 2} + 14 + 10} = 52 (max{C32 + T23, C43 + δ43′} > h313, max{C32 + T23, C53 + δ53′} > h323, max{C32 + T23, C23 + δ23′} > h333)X313 = 1, X323 = X333 = 0Therefore, the total weighted completion time is:∑iwiCiS = 5 × 38 + 6 × 54 + 9 × 52 + 7 × 35 + 3 × 50 = 1377.Figure 4
Example Gantt chart.


## 3.2.1. Encoding Solutions

To obtainNump food sources, each job sequence at stage 1 is yielded from discrete domain. So, a random perturbation procedure based on job numbers is developed to construct the initial job sequence ω1 = {ω1(1), ω1(2), … , ω1(i), … , ω1(n)}, where ωk(i) denotes the job on the ith position at stage k. Figure 2 illustrates an encoding scheme for an instance with five jobs where ω1(1) = 4.Figure 2
An instance with five jobs.


## 3.2.2. Decoding Food Sources

Because the goal is to minimize the total weighted completion time, the job with the earlier completion time at the immediately preceding stage is allowed to have a higher priority in order to improve on-time delivery between adjacent stages. So, the job sequence at stagek(k > 1) is determined according to their completion times at stage k −1.At stage 1, according to the above generated job sequence, the jobω1(1) is firstly assigned to the machine with the shortest actual processing time. To balance the utility of all the machines, the job ω1(2) chooses one machine with the shortest actual processing time from the unscheduled machines until the job ω1(m1) is arranged. For the unscheduled n-m1 jobs, compute the possible completion time of job ω1(i) (i > m1) equal to the sum of the completion time of the immediately preceding job on the same machine, the setup time and its actual processing time, and assign it to the machine with the shortest completion time.The job sequence at stagek(k > 1) is obtained in the increasing order of the completion time of all jobs at stage k−1. If some jobs have the same completion time, order these jobs randomly. Similar with stage 1, the first mk job is, respectively, arranged on the mk machines with the shortest actual processing time based on the sequence so that all machines are used. When it turns to ωk(i) (i > mk), the starting time of job ωk(i) is set to the maximum between its completion time at stage k − 1 plus the delivery time and the completion time of the immediately preceding job on the same machine plus the setup time. So, its completion time can be computed by its starting time plus the actual processing time. The assigned machine is decided with the shortest completion time. The process continues until all jobs are arranged in all stages, as shown in Figure 3.Figure 3
Decoding scheme using the rule of minimum processing/completion time of the free machine.


## 3.2.3. Population Initialization

As described above, an individual is denoted as a permutation of all jobs at stage 1. To generate a population, a simple random procedure is adopted whereNump individuals (food sources) are consistent with the size of employed bee or onlooker bee.To clearly illustrate the FFSP-SDJ&SDST and the decoding procedure, a simple example is introduced with 5 jobs and 3 stages where each stage has three machines. The release time and weight of jobs are given byri = [3, 1, 4, 2, 3] and wi = [5, 6, 9, 7, 3]; delivery time Tk,k+1 is, respectively, set to 3 and 2. For simplicity, it is assumed that the setup time δfijk is dependent on jobs and irrelevant to machines and denoted as δfi′. The basic processing time Pijk, deteriorating date hijk, deteriorating penalty time bijk, and setup time δfi′ are, respectively, given below:(10)The initial job sequence [4, 5, 1, 3, 2] at stage 1 is given and consistent with the instance illustrated in Figure2. The machine allocation and completion time of all jobs are computed as follows. The Gantt chart of a solution is drawn in Figure 4, where each operation is denoted by a rectangle labeled with a job number.Stage 1:ω1 = {4, 5, 1, 3, 2}C41 = r4 + min{P411, P421, P431} = 2 + min{5, 13, 11} = 7X411 = 1, X421 = X431 = 0C51 = r5 + min{P421, P431} = 3 + min{10, 19} = 13X521 = 1, X511 = X531 = 0C11 = r1 + P131 = 3 + 8 = 11X131 = 1, X111 = X121 = 0C31 = min{C41 + δ43′ + P311, C51 + δ53′ + P321, C11 + δ13′ + P331} = min{7 + 5 + 17, 13 + 1 + 6, 11 + 4 + 12} = 20 (C41 + δ43′ < h311, C51 + δ43′ < h321, C11 + δ13′ < h331)X321 = 1, X311 = X331 = 0C21 = min{C41 + δ42′ + P211, C31 + δ32′ + P221, C11 + δ12′ + P231} = {7 + 4 + 6, 20 + 2 + 4, 11 + 3 + 17} = 17 (C41 + δ42′ < h211, C31 + δ32′ < h221, C11 + δ12′ < h231)X211 = 1, X221 = X231 = 0Stage 2:ω2 = {4, 1, 5, 2, 3}C42 = C41 + T12 + min{P412, P422, P432} = 7 + 3 + min{14, 16, 12} = 22 (C41 + T12 < h412, C41 + T12 < h422, C41 + T12 < h432)X432 = 1, X412 = X422 = 0C12 = C11 + T12 + min{P112, P122} = 11 + 3 + min{16, 12} = 26 (C11 + T12 < h112, C11 + T12 < h122)X122 = 1, X112 = X132 = 0C52 = C51 + T12 + P512 = 13 + 3 + 12 = 28 (C51 + T12 < h512)X512 = 1, X522 = X532 = 0C22 = min{max{C21 + T12, C52 + δ52′} + P212, max{C21 + T12, C12 + δ12′} + P222, max{C21 + T12, C42 + δ42′} + P232} = min{max{17 + 3, 28 + 1} + 16, max{17 + 3, 26 + 3} + 10, max{17 + 3, 22 + 4} + 20 = 39 (max{C21 + T12, C52 + δ52′} < h212, max{C21 + T12, C12 + δ12′} < h222, max{C21 + T12, C42 + δ42′} < h232)X222 = 1, X212 = X232 = 1C32 = min{max{C31 + T12, C52 + δ53′} + P312 + b312, max{C31 + T12, C22 + δ23′} + P322 + b322, max{C31 + T12, C42 + δ43′} + P332 + b332} = min{max{20 + 3, 28 + 1} + 16 + 2, max{20 + 3, 39 + 2} + 11 + 6, max{20 + 3, 22 + 5} + 10 + 4} = 41 (max{C31 + T12, C52 + δ53′} > h312, max{C31 + T12, C22 + δ23′} > h322, max{C31 + T12, C42 + δ43′} > h332)X332 = 1, X312 = X322 = 0Stage 3:ω3 = {4, 1, 5, 2, 3}C43 = C42 + T23 + min{P413, P423, P433} = 22 + 2 + min{11, 14, 19} = 35 (C42 + T23 < h413, C42 + T23 < h423, C42 + T23 < h433)X413 = 1, X423 = X433 = 0C13 = C12 + T23 + min{P123, P133} = 26 + 2 + min{13, 10} = 38(C12 + T23 < h123, C12 + T23 < h133)X133 = 1, X113 = X123 = 0C53 = C52 + T23 + P523 + b523 = 28 + 2 + 15 + 5 = 50 (C52 + T23 > h523)X523 = 1, X513 = X533 = 0C23 = min{max{C22 + T23, C43 + δ42′} + P213 + b213, max{C22 + T23, C53 + δ52′} + P223 + b223, max{C22 + T23, C13 + δ12′} + P233 + b233} = min{max{39 + 2, 35 + 4} + 16 + 3, max{39 + 2, 50 + 1} + 11 + 1, max{39 + 2, 38 + 3} + 9 + 4} = 54(max{C22 + T23, C43 + δ42′} > h213, max{C22 + T23, C53 + δ52′} > h223, max{C22 + T23, C13 + δ12′} > h233)X233 = 1, X213 = X223 = 0C33 = min{max{C32 + T23, C43 + δ43′} + P313 + b313, max{C32 + T23, C53 + δ53′} + P323 + b323, max{C32 + T23, C23 + δ23′} + P333 + b333} = min{max{41 + 2, 35 + 5} + 7 + 2, max{41 + 2, 50 + 1} + 14 + 10, max{41 + 2, 54 + 2} + 14 + 10} = 52 (max{C32 + T23, C43 + δ43′} > h313, max{C32 + T23, C53 + δ53′} > h323, max{C32 + T23, C23 + δ23′} > h333)X313 = 1, X323 = X333 = 0Therefore, the total weighted completion time is:∑iwiCiS = 5 × 38 + 6 × 54 + 9 × 52 + 7 × 35 + 3 × 50 = 1377.Figure 4
Example Gantt chart.


## 3.3. Employed Bee Phase

In the employed bee phase, the new solutions are yielded by performing a genetic algorithm in order to find a neighboring food source from the current food source. The detailed GA is described below.


### 3.3.1. Fitness Evaluation and Selection Operation

Because the objective of this model is to minimize total weighted completion time, the reciprocal of the objective function is used as the fitness function. Therefore, the fitness function of theyth food source is estimated by f(y) = 1/∑iwiCiS. The roulette wheel selection is applied to determine the offspring where the individual with the higher fitness value has the higher probability to be selected.


### 3.3.2. Single-Point Crossover Operator

Step 1: randomly generate a real number Rc from the interval [0, 1] and compare it with the crossover probabilityCP. If Rc ≤ CP, two food sources A = {ωA(1), ωA(2), … , ωA(i), … , ωA(n)} and B = {ωB(1), ωB(2), … , ωB(i), … , ωB(n)} are picked out randomly from the population after the selection operation.
Step 2: a crossover positionθ (0 < θ<n) is generated randomly, and the single-point crossover is performed on A and B. As shown in Figure 5, the genes after position θ in A (B) are copied to A′ (B′) by keeping their position unchanged. The lacking genes in A′ (B′) are filled with the first θ genes in B (A) by maintaining their relative orders.
Step 3: find the job number (denoted asμ(θ1)) in {ωB(1), … , ωB(θ)} which emerges repeatedly in {ωA(θ + 1), … , ωA(n)} for A′. Similarly, find the job number v (θ2) from {ωA(1), … , ωA(θ)} with the same as {ωB(θ + 1), … , ωB(n) for B′.
Step 4: replaceμ(θ1) in the θ1th position with v (θ2) in A′ and v (θ2) in the θ2th position with μ(θ1) in B′. Therefore, two new food sources are formed.Figure 5
Crossover operator.


### 3.3.3. Multi-Insertion Mutation Operator

A real numberRm is randomly produced from the interval [0, 1]. If Rm is not larger than mutation probability MP, two job positions (θ1 and θ2) are randomly chosen in a permutation Z = {ω(1), ω(2), … , ω(n)}. Suppose that IT is an insertion number. Insert job ω(θ1) into (θ2 + q − 1)th (q ∈ {1, … , IT}) position, as shown in Figure 6. If θ2 + q − 1 > n, let θ2 + q − 1 equal to 1. Thus, IT new food sources are obtained to increase the diversity of population.Figure 6
Mutation operator.


## 3.3.1. Fitness Evaluation and Selection Operation

Because the objective of this model is to minimize total weighted completion time, the reciprocal of the objective function is used as the fitness function. Therefore, the fitness function of theyth food source is estimated by f(y) = 1/∑iwiCiS. The roulette wheel selection is applied to determine the offspring where the individual with the higher fitness value has the higher probability to be selected.


## 3.3.2. Single-Point Crossover Operator

Step 1: randomly generate a real number Rc from the interval [0, 1] and compare it with the crossover probabilityCP. If Rc ≤ CP, two food sources A = {ωA(1), ωA(2), … , ωA(i), … , ωA(n)} and B = {ωB(1), ωB(2), … , ωB(i), … , ωB(n)} are picked out randomly from the population after the selection operation.
Step 2: a crossover positionθ (0 < θ<n) is generated randomly, and the single-point crossover is performed on A and B. As shown in Figure 5, the genes after position θ in A (B) are copied to A′ (B′) by keeping their position unchanged. The lacking genes in A′ (B′) are filled with the first θ genes in B (A) by maintaining their relative orders.
Step 3: find the job number (denoted asμ(θ1)) in {ωB(1), … , ωB(θ)} which emerges repeatedly in {ωA(θ + 1), … , ωA(n)} for A′. Similarly, find the job number v (θ2) from {ωA(1), … , ωA(θ)} with the same as {ωB(θ + 1), … , ωB(n) for B′.
Step 4: replaceμ(θ1) in the θ1th position with v (θ2) in A′ and v (θ2) in the θ2th position with μ(θ1) in B′. Therefore, two new food sources are formed.Figure 5
Crossover operator.


## 3.3.3. Multi-Insertion Mutation Operator

A real numberRm is randomly produced from the interval [0, 1]. If Rm is not larger than mutation probability MP, two job positions (θ1 and θ2) are randomly chosen in a permutation Z = {ω(1), ω(2), … , ω(n)}. Suppose that IT is an insertion number. Insert job ω(θ1) into (θ2 + q − 1)th (q ∈ {1, … , IT}) position, as shown in Figure 6. If θ2 + q − 1 > n, let θ2 + q − 1 equal to 1. Thus, IT new food sources are obtained to increase the diversity of population.Figure 6
Mutation operator.


## 3.4. Onlooker Bee Phase

This phase uses the same roulette wheel selection applied in the employed bee phase for the food source returned by the employed bee. An onlooker bee chooses a food sourceω depending on its probability value FPω = f(ω)/∑yf(y).The modified variable neighborhood search (MVNS) is then applied for the above food sources. Firstly, the following neighborhood structures NSφ (φ = 1, 2, 3) are designed to help the onlooker bee explore possible promising neighbors.


### 3.4.1. Exchange Mutation

Two randomly generated positions (θ1 and θ2) are chosen and the corresponding jobs are swapped, as shown in Figure 7.Figure 7
Exchange mutation.


### 3.4.2. 2-Opt Invert Mutation

Two randomly selected jobs (ω(θ1) and ω(θ2)) are picked out and the subsequence between positions θ1 and θ2 are inverted, as illustrated in Figure 8.Figure 8
2-opt invert mutation.


### 3.4.3. Fragment Insert Mutation

Randomly select two different jobsω(θ1) and ω(θ2) and remove the subsequence between positions θ1 and θ2. A randomly selected position τ from the remaining positions is chosen, and the subsequence between θ1 and θ2 is inserted into position τ, as shown in Figure 9.Figure 9
Fragment insert mutation.Next, the new neighboring solution obtained by applying these search rules is evaluated and compared to the primary food source. If the new solution obtained by MVNS is better than the primary one, the primary food source is replaced and the new food source becomes a new member in the population.The procedure of MVNS is listed in Algorithm1.Algorithm 1:Procedure of MVNS.
Initialization: define NSφφ=1,2,3
for every food source ω do
Perturbation:
x = 1;
whilex < 3;
executeNS1 to ω
x = x + 1;
end while
Search:
φ=1
while the stopping criterion is not met
σ = execute NSφ to ω
If fσ>fω
ω=σ;
φ=φ;
else
φ=φ+1;
end if
end while
end for


## 3.4.1. Exchange Mutation

Two randomly generated positions (θ1 and θ2) are chosen and the corresponding jobs are swapped, as shown in Figure 7.Figure 7
Exchange mutation.


## 3.4.2. 2-Opt Invert Mutation

Two randomly selected jobs (ω(θ1) and ω(θ2)) are picked out and the subsequence between positions θ1 and θ2 are inverted, as illustrated in Figure 8.Figure 8
2-opt invert mutation.


## 3.4.3. Fragment Insert Mutation

Randomly select two different jobsω(θ1) and ω(θ2) and remove the subsequence between positions θ1 and θ2. A randomly selected position τ from the remaining positions is chosen, and the subsequence between θ1 and θ2 is inserted into position τ, as shown in Figure 9.Figure 9
Fragment insert mutation.Next, the new neighboring solution obtained by applying these search rules is evaluated and compared to the primary food source. If the new solution obtained by MVNS is better than the primary one, the primary food source is replaced and the new food source becomes a new member in the population.The procedure of MVNS is listed in Algorithm1.Algorithm 1:Procedure of MVNS.
Initialization: define NSφφ=1,2,3
for every food source ω do
Perturbation:
x = 1;
whilex < 3;
executeNS1 to ω
x = x + 1;
end while
Search:
φ=1
while the stopping criterion is not met
σ = execute NSφ to ω
If fσ>fω
ω=σ;
φ=φ;
else
φ=φ+1;
end if
end while
end for


## 3.5. Scout Bee Phase

If a food source has no improvement duringlimit iterations, this employed bee will become a scout bee. A greedy heuristic is developed to produce a new food source including two steps.


### 3.5.1. Destruction Procedure

Jobω(θ) is randomly chosen and removed from the food source Z = {ω(1), ω(2), … , ω(θ − 1), ω(θ), ω(θ + 1), … , ω(n)} so that a remaining food source Z − θ = {ω(1), ω(2), … , ω(θ−1), ω(θ + 1) … , ω(n)} is generated.


### 3.5.2. Reconstruction Procedure

The selected jobω(θ) is sequentially inserted into another randomly selected position of Z-θ, and the corresponding fitness value is calculated. Compare these new generated neighboring solutions and choose the best one with the largest fitness value as the current solution.


## 3.5.1. Destruction Procedure

Jobω(θ) is randomly chosen and removed from the food source Z = {ω(1), ω(2), … , ω(θ − 1), ω(θ), ω(θ + 1), … , ω(n)} so that a remaining food source Z − θ = {ω(1), ω(2), … , ω(θ−1), ω(θ + 1) … , ω(n)} is generated.


## 3.5.2. Reconstruction Procedure

The selected jobω(θ) is sequentially inserted into another randomly selected position of Z-θ, and the corresponding fitness value is calculated. Compare these new generated neighboring solutions and choose the best one with the largest fitness value as the current solution.


## 4. Experimental Testing and Analyses

In order to verify the performance of the presented algorithm, an extensive experimental comparison with GA, VNS, and GH are provided. All algorithms are coded in MATLAB R2014a and run on a Micro Core i5-4210U, 1.7 GHz PC with 4.00 GB. The orthogonal test method is used to optimize the parameters in the IDABC algorithm. The maximum iteration number is set to 100, and the longest running time is 800 s. Test instances are randomly generated as follows:(1)
Job numbern ∈ {10, 40, 70, 100, 120} and stage number S ∈ {2, 3, 5}. The setting of mk at stage k is divided into two cases: X and Y. In case X, mk = 3 for k = 1,… , S and in case Y, mk ∈ U[1, 5].(2)
The basic processing time, job release time, delivery time between stages, setup time, and job weight are generated from the following discrete uniform distributions:Pijk ∈ U[1, 20], ri ∈ U[1, 5], Tk,k+1 ∈ U[1, 6], δfijk ∈ U[1, 6], and wi ∈ U[1, 10].(3)
The penalty timebijk is produced according to the uniform distribution U[1, 20ψ] where ψ = 0.5.(4)
Parameter Δ is determined by equation (11), where P¯ik represents the average processing time of job i at stage k defined as equation (12). The deterioration dates are randomly generated from a uniform distribution on the intervals H1 ∈ U[1, 0.7Δ] and H2 ∈ U[0.3Δ, Δ].(11)Δ=∑i=1n∑k=1SPik¯,(12)Pik¯=∑j=1mkPijkmk.(5)
Based on the simulation experiments of different parameter values, the insertion numberIT applied for employed bees is set to IT = n × 20%. The number of cycles is used as the stopping criterion of MVNS for onlooker bees and it is fixed to 5.


### 4.1. Parameter Settings

In the presented IDABC algorithm, four main parameters are included, i.e.,Nump, limit, CP, and MP. They are optimized through orthogonal experiments where each parameter is viewed as a factor changing at three different levels (see Table 1), resulting in 32 = 9 experiments.Table 1
Factor settings of orthogonal experiments.
Factor level
Factors
Nump
limit
CP
MP
1
60
3
0.2
0.2
2
80
5
0.5
0.5
3
90
8
0.8
0.8An instance 30 × 3 × 3 (n = 30, S = 3, and mk = 3) is taken as an example, and every parameter combination is tested ten times. Table 2 gives the orthogonal parameter table L9(34) with nine groups of the parameters, and the mean objective values where Qy (y = 1, 2, 3) in column z (z = 1, 2, 3, 4) denotes the average of the results for three sets of factor z at level y. According to the experimental results, the parameters are set as follows: Nump = 80, limit = 8, CP = 0.8, and MP = 0.5.Table 2
Orthogonal array and average objective value.
Factor
Average objective value
Nump
limit
CP
MP
1
60
3
0.2
0.2
7328.7
2
60
5
0.5
0.5
7319.5
3
60
8
0.8
0.8
7325.3
4
80
3
0.5
0.8
7330.6
5
80
5
0.8
0.2
7312.3
6
80
8
0.2
0.5
7302.7
7
90
3
0.8
0.5
7317.0
8
90
5
0.2
0.8
7335.4
9
90
8
0.5
0.2
7308.5
Q
1
7324.5
7325.4
7322.3
7316.5
Q
2
7315.2
7322.4
7319.5
7313.1
Q
3
7320.3
7312.2
7318.3
7330.4


### 4.2. Simulation Results and Analyses

To make a fair comparison, all the algorithms utilize the same CPU time limit as a termination criterion. All of the cited parameters {n, S, mk, hijk} mentioned above result in a total of 5 × 3 × 2 × 2 = 60 different combinations. For each combination, ten random instances are generated. Thus, there are 360 small- and medium-sized instances and 240 large-sized instances to test these algorithms. The relative percentage increase (RPI) is defined as follows:(13)RPI=Oζ−ObestObest×100,where Oζ is the objective value obtained by different metaheuristic algorithm and Obest is the best objective value found by all algorithms.Therefore, average RPI (ARPI) can be computed and used as a performance measure to compare the performance of the algorithms. Tables3 and 4 list the testing results for different sized problems where columns Avg and Min represent the average and minimum of 10 runs of each algorithm, respectively.Table 3
Comparison of the four algorithms for small and medium instances.
n×S×mk×h
GA
GH
VNS
IDABC
CPU time
Avg
Min
ARPI
Avg
Min
ARPI
Avg
Min
ARPI
Avg
Min
ARPI
n×S×mk×H1
10×2×X×H1
1754.4
1710
3.32
1735.8
1704
2.23
1725.6
1704
1.63
1709.4
1698
0.07
25.52
10×2×Y×H1
778.2
765
2.74
787.8
765
2.98
776.7
765
1.53
765.3
765
0.04
27.19
10×3×X×H1
1081.8
1065
1.73
1081.2
1053
2.67
1130.4
1071
7.35
1056.9
1053
0.37
49.66
10×3×Y×H1
1592.8
1552
3.32
1594.4
1572
2.73
1577.2
1552
1.62
1564.8
1552
0.82
42.11
10×5×X×H1
4236.1
4135
2.62
4225.1
4119
3.05
4225.9
4101
3.07
4216.0
4100
2.83
78.03
10×5×Y×H1
4430.3
4354
1.75
4435.9
4368
1.88
4363.8
4354
0.22
4358.2
4354
0.10
81.33
40×2×X×H1
17268.0
16626
3.35
17857.8
17622
11.01
17208.6
16836
6.98
16736.4
16808
4.04
146.63
40×2×Y×H1
16813.6
16192
4.14
17879.2
17480
10.75
17276.0
16872
7.01
16409.6
16144
1.65
139.72
40×3×X×H1
39717.0
38140
4.55
40121.3
39170
5.61
40267.2
38860
5.99
38991.2
37990
2.63
262.19
40×3×Y×H1
16295.5
15995
3.39
16700.5
16345
5.97
16326.0
15875
3.59
15999.5
15760
1.52
189.74
40×5×X×H1
17268.0
16060
7.52
17357.8
16626
8.08
17208.6
16836
7.15
16736.4
16086
4.21
346.06
40×5×Y×H1
45633.6
44064
6.29
45929.7
44784
6.99
44965.8
44379
4.74
44474.4
42930
3.60
482.23
70×2×X×H1
7312.4
7115
8.70
7445.0
7204
10.67
7397.3
7247
9.96
7163.6
6727
6.49
294.42
70×2×Y×H1
42046.9
40474
8.64
43504.3
41874
12.41
43201.9
42287
11.62
40496.4
38703
4.63
281.67
70×3×X×H1
39585.9
38406
8.68
39309.6
38298
7.93
39376.8
37869
8.11
37216.2
36423
2.18
496.63
70×3×Y×H1
67027.4
64744
6.36
67655.0
65840
7.35
66655.0
64952
5.76
63835.7
63033
1.29
429.50
70×5×X×H1
51786.4
50873
7.44
51987.9
51003
7.84
51972.9
51030
7.81
48987.1
48201
1.61
734.33
70×5×Y×H1
81461.2
80544
7.00
81714.2
80664
7.33
81300.2
80210
6.79
77476.8
76134
1.76
789.32
Average
25338.3
24601
5.09
25629.0
25027
6.53
25386.4
24822
5.61
24344.1
23803
2.21
272.02
n×S×mk×H2
10×2×X×H2
1126.0
1108
1.62
1120.0
1108
1.08
1122.4
1108
1.30
1108.4
1108
0.04
25.52
10×2×Y×H2
1020.6
999
2.47
1029.6
1011
3.37
1020.9
1014
2.50
1001.7
996
0.58
27.19
10×3×X×H2
469.9
461
2.37
470.5
461
2.51
468.0
459
1.96
460.4
459
0.31
49.66
10×3×Y×H2
3003.0
2926
2.63
2974.0
2926
0.71
2946.3
2926
0.69
2926.4
2926
0.01
42.11
10×5×X×H2
4823.6
4802
2.02
4831.9
4801
2.20
4856.8
4802
2.72
4804.1
4728
1.61
78.03
10×5×Y×H2
4261.6
4200
2.05
4274.4
4200
2.36
4215.2
4176
0.94
4176.8
4176
0.02
81.33
40×2×X×H2
8677.2
8304
4.76
8999.1
8607
8.65
8837.4
8676
6.69
8530.8
8232
2.99
146.63
40×2×Y×H2
6795.0
6507
5.30
7229.7
7050
12.04
6890.4
6711
6.78
6614.4
6453
2.50
139.72
40×3×X×H2
22859.4
22218
7.41
25345.8
22872
10.64
22761.0
22110
6.95
21937.8
21282
3.08
262.19
40×3×Y×H2
3179.7
3090
3.30
3248.5
3195
5.54
3196.8
3127
3.86
3122.6
3078
1.45
189.74
40×5×X×H2
39642.2
39099
4.21
39505.9
38739
3.89
39667.5
39494
4.32
38536.0
38025
1.34
346.06
40×5×Y×H2
34203.4
33327
4.14
34980.4
34174
6.50
34365.8
33782
4.63
33717.6
32884
2.66
482.23
70×2×X×H2
52081.9
50935
5.76
53207.0
52353
8.05
52888.6
52086
7.40
50476.3
49245
2.50
294.42
70×2×Y×H2
68669.4
65230
7.97
68177.7
66990
7.20
68780.6
66240
8.14
66003.8
63600
3.78
281.67
70×3×X×H2
50135.6
48936
4.57
51162.1
49422
6.71
50552.6
49806
5.44
48794.7
47946
1.77
496.63
70×3×Y×H2
24507.3
24162
5.15
24948.9
24729
7.04
24859.2
24399
6.66
24277.6
23308
4.15
429.50
70×5×X×H2
61215.0
60342
7.73
61804.5
61083
8.77
61518.9
60497
8.26
57294.3
56823
0.83
734.33
70×5×Y×H2
54368.4
53136
5.93
55100.4
53332
7.35
54555.6
53724
6.29
52626.0
51328
2.53
789.32
Average
24502.2
23877
4.32
24911.7
24281
5.81
24639.1
24174
4.75
23689.4
23144
1.79
272.02Table 4
Comparison of the four algorithms for large instances.
n × S × mk × h
GA
GH
VNS
IDABC
CPU time
Avg
Min
ARPI
Avg
Min
ARPI
Avg
Min
ARPI
Avg
Min
ARPI
n×S×mk×H1
100×2×X×H1
85002.6
80394
8.89
86004.0
83364
10.18
86182.8
85224
10.41
81648.6
78060
4.60
563.30
100×2×Y×H1
164413.6
162340
8.69
165181.4
164110
9.20
164004.7
160820
8.42
154190.4
151270
1.93
570.34
100×3×X×H1
46886.7
45298
6.64
47590.7
45573
8.24
47375.1
45339
7.75
44838.0
43968
1.98
674.16
100×3×Y×H1
162497.4
162030
7.42
163215.7
162620
7.90
163404.8
162140
8.02
154372.7
151270
2.05
691.13
100×5×X×H1
111831.0
103973
9.28
113080.3
109763
10.50
111534.2
108734
8.99
105269.5
102336
2.87
800.00
100×5×Y×H1
27299.4
26998
8.45
27476.4
27194
9.15
27325.2
26392
8.55
25724.9
25173
2.19
800.00
120×2×X×H1
67987.9
64257
8.94
69349.3
68724
11.13
68148.8
65922
9.20
64419.3
62406
3.23
755.67
120×2×Y×H1
181648.0
174368
8.44
182875.2
175072
9.18
182338.4
176600
8.86
172412.8
167504
2.93
759.46
120×3×X×H1
250104.5
244510
11.00
254437.3
251670
12.93
249121.0
242280
10.57
232410.3
225310
3.15
800.00
120×3×Y×H1
156395.6
155288
8.31
158088.0
157040
9.48
157724.8
157128
9.23
147907.2
144400
2.43
800.00
120×5×X×H1
353503.5
349220
7.68
355892.7
351950
8.40
534840.2
351250
8.08
331862
328300
1.08
800.00
120×5×Y×H1
317292.8
303048
6.99
318821.6
312832
7.50
314614.4
311240
6.08
302820.8
296576
2.11
800.00
Average
160405.3
155977
8.39
161834.4
159159
9.48
175551.2
157756
8.68
151489.7
148048
2.55
734.51
n×S×mk×H2
100×2×X×H2
30774.4
30210
4.67
31132.6
30408
5.89
30855.8
30238
4.94
30072.0
29402
2.28
563.30
100×2×Y×H2
36300.3
34101
6.45
36957.0
36483
8.38
36704.4
36333
7.63
35929.5
34662
5.36
570.34
100×3×X×H2
46820.3
44035
9.09
48176.5
46020
12.25
47780.5
45895
11.32
42967.4
42920
0.11
674.16
100×3×Y×H2
96520.6
93430
5.07
97543.1
93552
6.18
97208.4
93372
5.82
92908.5
91886
1.13
691.13
100×5×X×H2
233298.4
229788
9.26
235390.2
222711
10.24
233203.8
227548
9.22
215706.6
213516
1.03
800.00
100×5×Y×H2
26102.1
25189
9.10
26158.0
25455
9.33
25978.5
25208
8.58
24693.1
23925
3.21
800.00
120×2×X×H2
146292.3
141300
3.81
150106.5
148113
6.52
148338.0
144729
5.26
143712.9
140992
1.98
755.67
120×2×Y×H2
110558.5
107555
9.29
109579.0
108495
8.32
108672.5
107310
7.42
105177.0
101165
3.97
759.46
120×3×X×H2
244468.6
238810
6.30
248267.6
237611
7.96
246137.8
240340
7.03
235697.3
229970
2.49
800.00
120×3×Y×H2
73899.5
70101
7.02
74268.4
72523
7.55
74104.5
70653
7.31
70228.2
69054
1.70
800.00
120×5×X×H2
282084.6
271104
8.07
283192.5
271120
8.50
279734.2
270824
7.18
266991.9
261000
2.30
800.00
120×5×Y×H2
258686.4
245560
7.55
267808.8
264159
11.34
264243.0
257390
9.86
247565.3
240536
2.92
800.00
Average
132150.5
127599
7.14
134048.4
129721
8.54
132746.8
129153
7.63
125970.8
123252
2.37
734.51It can be concluded from Tables3 and 4 that the proposed IDABC algorithm performs better than GA, GH, and VNS consistently. Precisely, the IDABC algorithm has the least ARPI 2.00% and 2.46%, respectively, for small-medium-sized problems and large-sized problems, which are smaller than the ARPI of the GA, GH, and VNS, 4.75% and 7.77%, 6.17% and 9.01%, and 5.18% and 8.16%. This indicates that the results obtained by the presented algorithm are closer to the optimal or near-optimal solutions. For all instances, the average objective values of IDABC outperform consistently the other algorithms, and the best optimal or near-optimal solutions are almost found by IDABC except two combinations (40 × 5 × X and 100 × 2 × Y) where the IDABC is slightly worse than the GA.When the scale of the problem increases, the gap of ARPI values become bigger compared with the other algorithms. For the two tested deterioration dates ofH1 and H2, it can be seen that the algorithms under H2 perform better.To explain the difference statistically among the algorithms, the analysis of variance (ANOVA) is carried out to check the results. The method of residual diagnosis is applied, and the results support the main hypotheses of the ANOVA. The means plot of ARPI and least significant deviation (LSD) intervals for all the algorithms are shown in Figure10 with LSD intervals at a 95% confidence level. It can be seen from the figure that the intervals obtained by the IDABC algorithm differ significantly from the other algorithms and it provides better stability.Figure 10
Means plot and LSD intervals with a 95% confidence level for different algorithms. (a) Small- and medium-scale problems. (b) Large-scale problems.
(a)
(b)


### 4.3. IDABC Effect Experiments

This section will further study the relationship between the algorithm performance and iteration number. Taken = {10, 70, 120}, S = 3 and mk = 3, and the job deterioration date is set to H1. Figure 11 shows the convergence trend of the three instances for the four algorithms. As seen in the figure, as the job size increases, the IDABC can find better objective values within a less iteration number and have a faster convergence.Figure 11
The convergence of different problem sizes. (a) 10 jobs. (b) 70 jobs. (c) 120 jobs.
(a)
(b)
(c)


## 4.1. Parameter Settings

In the presented IDABC algorithm, four main parameters are included, i.e.,Nump, limit, CP, and MP. They are optimized through orthogonal experiments where each parameter is viewed as a factor changing at three different levels (see Table 1), resulting in 32 = 9 experiments.Table 1
Factor settings of orthogonal experiments.
Factor level
Factors
Nump
limit
CP
MP
1
60
3
0.2
0.2
2
80
5
0.5
0.5
3
90
8
0.8
0.8An instance 30 × 3 × 3 (n = 30, S = 3, and mk = 3) is taken as an example, and every parameter combination is tested ten times. Table 2 gives the orthogonal parameter table L9(34) with nine groups of the parameters, and the mean objective values where Qy (y = 1, 2, 3) in column z (z = 1, 2, 3, 4) denotes the average of the results for three sets of factor z at level y. According to the experimental results, the parameters are set as follows: Nump = 80, limit = 8, CP = 0.8, and MP = 0.5.Table 2
Orthogonal array and average objective value.
Factor
Average objective value
Nump
limit
CP
MP
1
60
3
0.2
0.2
7328.7
2
60
5
0.5
0.5
7319.5
3
60
8
0.8
0.8
7325.3
4
80
3
0.5
0.8
7330.6
5
80
5
0.8
0.2
7312.3
6
80
8
0.2
0.5
7302.7
7
90
3
0.8
0.5
7317.0
8
90
5
0.2
0.8
7335.4
9
90
8
0.5
0.2
7308.5
Q
1
7324.5
7325.4
7322.3
7316.5
Q
2
7315.2
7322.4
7319.5
7313.1
Q
3
7320.3
7312.2
7318.3
7330.4


## 4.2. Simulation Results and Analyses

To make a fair comparison, all the algorithms utilize the same CPU time limit as a termination criterion. All of the cited parameters {n, S, mk, hijk} mentioned above result in a total of 5 × 3 × 2 × 2 = 60 different combinations. For each combination, ten random instances are generated. Thus, there are 360 small- and medium-sized instances and 240 large-sized instances to test these algorithms. The relative percentage increase (RPI) is defined as follows:(13)RPI=Oζ−ObestObest×100,where Oζ is the objective value obtained by different metaheuristic algorithm and Obest is the best objective value found by all algorithms.Therefore, average RPI (ARPI) can be computed and used as a performance measure to compare the performance of the algorithms. Tables3 and 4 list the testing results for different sized problems where columns Avg and Min represent the average and minimum of 10 runs of each algorithm, respectively.Table 3
Comparison of the four algorithms for small and medium instances.
n×S×mk×h
GA
GH
VNS
IDABC
CPU time
Avg
Min
ARPI
Avg
Min
ARPI
Avg
Min
ARPI
Avg
Min
ARPI
n×S×mk×H1
10×2×X×H1
1754.4
1710
3.32
1735.8
1704
2.23
1725.6
1704
1.63
1709.4
1698
0.07
25.52
10×2×Y×H1
778.2
765
2.74
787.8
765
2.98
776.7
765
1.53
765.3
765
0.04
27.19
10×3×X×H1
1081.8
1065
1.73
1081.2
1053
2.67
1130.4
1071
7.35
1056.9
1053
0.37
49.66
10×3×Y×H1
1592.8
1552
3.32
1594.4
1572
2.73
1577.2
1552
1.62
1564.8
1552
0.82
42.11
10×5×X×H1
4236.1
4135
2.62
4225.1
4119
3.05
4225.9
4101
3.07
4216.0
4100
2.83
78.03
10×5×Y×H1
4430.3
4354
1.75
4435.9
4368
1.88
4363.8
4354
0.22
4358.2
4354
0.10
81.33
40×2×X×H1
17268.0
16626
3.35
17857.8
17622
11.01
17208.6
16836
6.98
16736.4
16808
4.04
146.63
40×2×Y×H1
16813.6
16192
4.14
17879.2
17480
10.75
17276.0
16872
7.01
16409.6
16144
1.65
139.72
40×3×X×H1
39717.0
38140
4.55
40121.3
39170
5.61
40267.2
38860
5.99
38991.2
37990
2.63
262.19
40×3×Y×H1
16295.5
15995
3.39
16700.5
16345
5.97
16326.0
15875
3.59
15999.5
15760
1.52
189.74
40×5×X×H1
17268.0
16060
7.52
17357.8
16626
8.08
17208.6
16836
7.15
16736.4
16086
4.21
346.06
40×5×Y×H1
45633.6
44064
6.29
45929.7
44784
6.99
44965.8
44379
4.74
44474.4
42930
3.60
482.23
70×2×X×H1
7312.4
7115
8.70
7445.0
7204
10.67
7397.3
7247
9.96
7163.6
6727
6.49
294.42
70×2×Y×H1
42046.9
40474
8.64
43504.3
41874
12.41
43201.9
42287
11.62
40496.4
38703
4.63
281.67
70×3×X×H1
39585.9
38406
8.68
39309.6
38298
7.93
39376.8
37869
8.11
37216.2
36423
2.18
496.63
70×3×Y×H1
67027.4
64744
6.36
67655.0
65840
7.35
66655.0
64952
5.76
63835.7
63033
1.29
429.50
70×5×X×H1
51786.4
50873
7.44
51987.9
51003
7.84
51972.9
51030
7.81
48987.1
48201
1.61
734.33
70×5×Y×H1
81461.2
80544
7.00
81714.2
80664
7.33
81300.2
80210
6.79
77476.8
76134
1.76
789.32
Average
25338.3
24601
5.09
25629.0
25027
6.53
25386.4
24822
5.61
24344.1
23803
2.21
272.02
n×S×mk×H2
10×2×X×H2
1126.0
1108
1.62
1120.0
1108
1.08
1122.4
1108
1.30
1108.4
1108
0.04
25.52
10×2×Y×H2
1020.6
999
2.47
1029.6
1011
3.37
1020.9
1014
2.50
1001.7
996
0.58
27.19
10×3×X×H2
469.9
461
2.37
470.5
461
2.51
468.0
459
1.96
460.4
459
0.31
49.66
10×3×Y×H2
3003.0
2926
2.63
2974.0
2926
0.71
2946.3
2926
0.69
2926.4
2926
0.01
42.11
10×5×X×H2
4823.6
4802
2.02
4831.9
4801
2.20
4856.8
4802
2.72
4804.1
4728
1.61
78.03
10×5×Y×H2
4261.6
4200
2.05
4274.4
4200
2.36
4215.2
4176
0.94
4176.8
4176
0.02
81.33
40×2×X×H2
8677.2
8304
4.76
8999.1
8607
8.65
8837.4
8676
6.69
8530.8
8232
2.99
146.63
40×2×Y×H2
6795.0
6507
5.30
7229.7
7050
12.04
6890.4
6711
6.78
6614.4
6453
2.50
139.72
40×3×X×H2
22859.4
22218
7.41
25345.8
22872
10.64
22761.0
22110
6.95
21937.8
21282
3.08
262.19
40×3×Y×H2
3179.7
3090
3.30
3248.5
3195
5.54
3196.8
3127
3.86
3122.6
3078
1.45
189.74
40×5×X×H2
39642.2
39099
4.21
39505.9
38739
3.89
39667.5
39494
4.32
38536.0
38025
1.34
346.06
40×5×Y×H2
34203.4
33327
4.14
34980.4
34174
6.50
34365.8
33782
4.63
33717.6
32884
2.66
482.23
70×2×X×H2
52081.9
50935
5.76
53207.0
52353
8.05
52888.6
52086
7.40
50476.3
49245
2.50
294.42
70×2×Y×H2
68669.4
65230
7.97
68177.7
66990
7.20
68780.6
66240
8.14
66003.8
63600
3.78
281.67
70×3×X×H2
50135.6
48936
4.57
51162.1
49422
6.71
50552.6
49806
5.44
48794.7
47946
1.77
496.63
70×3×Y×H2
24507.3
24162
5.15
24948.9
24729
7.04
24859.2
24399
6.66
24277.6
23308
4.15
429.50
70×5×X×H2
61215.0
60342
7.73
61804.5
61083
8.77
61518.9
60497
8.26
57294.3
56823
0.83
734.33
70×5×Y×H2
54368.4
53136
5.93
55100.4
53332
7.35
54555.6
53724
6.29
52626.0
51328
2.53
789.32
Average
24502.2
23877
4.32
24911.7
24281
5.81
24639.1
24174
4.75
23689.4
23144
1.79
272.02Table 4
Comparison of the four algorithms for large instances.
n × S × mk × h
GA
GH
VNS
IDABC
CPU time
Avg
Min
ARPI
Avg
Min
ARPI
Avg
Min
ARPI
Avg
Min
ARPI
n×S×mk×H1
100×2×X×H1
85002.6
80394
8.89
86004.0
83364
10.18
86182.8
85224
10.41
81648.6
78060
4.60
563.30
100×2×Y×H1
164413.6
162340
8.69
165181.4
164110
9.20
164004.7
160820
8.42
154190.4
151270
1.93
570.34
100×3×X×H1
46886.7
45298
6.64
47590.7
45573
8.24
47375.1
45339
7.75
44838.0
43968
1.98
674.16
100×3×Y×H1
162497.4
162030
7.42
163215.7
162620
7.90
163404.8
162140
8.02
154372.7
151270
2.05
691.13
100×5×X×H1
111831.0
103973
9.28
113080.3
109763
10.50
111534.2
108734
8.99
105269.5
102336
2.87
800.00
100×5×Y×H1
27299.4
26998
8.45
27476.4
27194
9.15
27325.2
26392
8.55
25724.9
25173
2.19
800.00
120×2×X×H1
67987.9
64257
8.94
69349.3
68724
11.13
68148.8
65922
9.20
64419.3
62406
3.23
755.67
120×2×Y×H1
181648.0
174368
8.44
182875.2
175072
9.18
182338.4
176600
8.86
172412.8
167504
2.93
759.46
120×3×X×H1
250104.5
244510
11.00
254437.3
251670
12.93
249121.0
242280
10.57
232410.3
225310
3.15
800.00
120×3×Y×H1
156395.6
155288
8.31
158088.0
157040
9.48
157724.8
157128
9.23
147907.2
144400
2.43
800.00
120×5×X×H1
353503.5
349220
7.68
355892.7
351950
8.40
534840.2
351250
8.08
331862
328300
1.08
800.00
120×5×Y×H1
317292.8
303048
6.99
318821.6
312832
7.50
314614.4
311240
6.08
302820.8
296576
2.11
800.00
Average
160405.3
155977
8.39
161834.4
159159
9.48
175551.2
157756
8.68
151489.7
148048
2.55
734.51
n×S×mk×H2
100×2×X×H2
30774.4
30210
4.67
31132.6
30408
5.89
30855.8
30238
4.94
30072.0
29402
2.28
563.30
100×2×Y×H2
36300.3
34101
6.45
36957.0
36483
8.38
36704.4
36333
7.63
35929.5
34662
5.36
570.34
100×3×X×H2
46820.3
44035
9.09
48176.5
46020
12.25
47780.5
45895
11.32
42967.4
42920
0.11
674.16
100×3×Y×H2
96520.6
93430
5.07
97543.1
93552
6.18
97208.4
93372
5.82
92908.5
91886
1.13
691.13
100×5×X×H2
233298.4
229788
9.26
235390.2
222711
10.24
233203.8
227548
9.22
215706.6
213516
1.03
800.00
100×5×Y×H2
26102.1
25189
9.10
26158.0
25455
9.33
25978.5
25208
8.58
24693.1
23925
3.21
800.00
120×2×X×H2
146292.3
141300
3.81
150106.5
148113
6.52
148338.0
144729
5.26
143712.9
140992
1.98
755.67
120×2×Y×H2
110558.5
107555
9.29
109579.0
108495
8.32
108672.5
107310
7.42
105177.0
101165
3.97
759.46
120×3×X×H2
244468.6
238810
6.30
248267.6
237611
7.96
246137.8
240340
7.03
235697.3
229970
2.49
800.00
120×3×Y×H2
73899.5
70101
7.02
74268.4
72523
7.55
74104.5
70653
7.31
70228.2
69054
1.70
800.00
120×5×X×H2
282084.6
271104
8.07
283192.5
271120
8.50
279734.2
270824
7.18
266991.9
261000
2.30
800.00
120×5×Y×H2
258686.4
245560
7.55
267808.8
264159
11.34
264243.0
257390
9.86
247565.3
240536
2.92
800.00
Average
132150.5
127599
7.14
134048.4
129721
8.54
132746.8
129153
7.63
125970.8
123252
2.37
734.51It can be concluded from Tables3 and 4 that the proposed IDABC algorithm performs better than GA, GH, and VNS consistently. Precisely, the IDABC algorithm has the least ARPI 2.00% and 2.46%, respectively, for small-medium-sized problems and large-sized problems, which are smaller than the ARPI of the GA, GH, and VNS, 4.75% and 7.77%, 6.17% and 9.01%, and 5.18% and 8.16%. This indicates that the results obtained by the presented algorithm are closer to the optimal or near-optimal solutions. For all instances, the average objective values of IDABC outperform consistently the other algorithms, and the best optimal or near-optimal solutions are almost found by IDABC except two combinations (40 × 5 × X and 100 × 2 × Y) where the IDABC is slightly worse than the GA.When the scale of the problem increases, the gap of ARPI values become bigger compared with the other algorithms. For the two tested deterioration dates ofH1 and H2, it can be seen that the algorithms under H2 perform better.To explain the difference statistically among the algorithms, the analysis of variance (ANOVA) is carried out to check the results. The method of residual diagnosis is applied, and the results support the main hypotheses of the ANOVA. The means plot of ARPI and least significant deviation (LSD) intervals for all the algorithms are shown in Figure10 with LSD intervals at a 95% confidence level. It can be seen from the figure that the intervals obtained by the IDABC algorithm differ significantly from the other algorithms and it provides better stability.Figure 10
Means plot and LSD intervals with a 95% confidence level for different algorithms. (a) Small- and medium-scale problems. (b) Large-scale problems.
(a)
(b)


## 4.3. IDABC Effect Experiments

This section will further study the relationship between the algorithm performance and iteration number. Taken = {10, 70, 120}, S = 3 and mk = 3, and the job deterioration date is set to H1. Figure 11 shows the convergence trend of the three instances for the four algorithms. As seen in the figure, as the job size increases, the IDABC can find better objective values within a less iteration number and have a faster convergence.Figure 11
The convergence of different problem sizes. (a) 10 jobs. (b) 70 jobs. (c) 120 jobs.
(a)
(b)
(c)


## 5. Conclusions

In this paper, the FFSP with step-deteriorating jobs and sequence-dependent setup times is studied, and the integer programming model is established. An improved discrete artificial bee colony (IDABC) algorithm is then presented for minimum total weighted completion time. In the IDABC algorithm, a permutation-based encoding scheme is applied, decoding with the dynamic generation mechanism to obtain the initial solutions. A genetic algorithm is designed to new offspring solutions in the employed bee phase; a modified variable neighborhood search with perturbation and searching operators is then developed in onlooker bee phase to generate new neighborhood food sources; scouting bees adopt a greedy heuristic to further find better solutions for the nonimproved individual. The comparison with GA, GH, and MVNS verifies the effectivity and efficiency of the IDABC algorithm, especially for large-scale problems.Some future work may be done from the following aspects: (1) our study involves the single-objective problem of total weighted completion time; the application of the developed algorithm may be extended to solve other single or multiobjective scheduling problems; (2) it is interesting to present some heuristic algorithms to generate initial population instead of the random procedure used in the paper so as to further improve the algorithm performance.



---

*Source: 8520503-2019-12-12.xml*