Metaheuristics for a Flow Shop Scheduling Problem with Urgent Jobs and Limited Waiting Times

: This study considers a scheduling problem for a ﬂow shop with urgent jobs and limited waiting times. The urgent jobs and limited waiting times are major considerations for scheduling in semiconductor manufacturing systems. The objective function is to minimize a weighted sum of total tardiness of urgent jobs and the makespan of normal jobs. This problem is formulated in mixed integer programming (MIP). By using a commercial optimization solver, the MIP can be used to ﬁnd an optimal solution. However, because this problem is proved to be NP-hard, solving to optimality requires a signiﬁcantly long computation time for a practical size problem. Therefore, this study adopts metaheuristic algorithms to obtain a good solution quickly. To complete this, two metaheuristic algorithms (an iterated greedy algorithm and a simulated annealing algorithm) are proposed, and a series of computational experiments were performed to examine the effectiveness and efﬁciency of the proposed algorithms.


Introduction
This study deals with a scheduling problem in a flow shop with urgent jobs and limited waiting times. Urgent jobs refer to jobs that must be processed faster due to higher urgency than normal jobs. Limited waiting time is to limit the time that a job completed on the first machine waits to start the next process on the second machine. These two scheduling requirements are very common in semiconductor manufacturing processes.
In general, semiconductor manufacturing systems are operated in high-variety smallvolume production and make-to-order policies. Nowadays, because the semiconductor foundry market is growing rapidly, customer order management and fast demand response capability are considered as very important competitive factors. To cope with this environment, jobs are classified into two (or more) groups depending on the urgencies of customer orders, i.e., a normal job and an urgent job. In this case, for normal jobs, a schedule is established with the objective of minimizing the maximum completion time (makespan) to raise throughput rate. However, urgent jobs are required to start processing as early as possible when they arrive at workstations. To accomplish this, each urgent job has a calculated due date, assuming that it is processed through all steps without waiting immediately upon arrival. That is, the due date on each step is equal to the earliest possible completion time. Based on this assumption, the objective of minimizing total tardiness can be considered as the scheduling measure for the urgent jobs. If normal and urgent jobs are required to be scheduled together on the same machine or shop, the scheduling measure is treated as a multi-objective because the makespan for normal jobs and total tardiness for urgent jobs are combined. Especially in semiconductor manufacturing, the urgent job

Previous Studies
This section introduces the previous studies related to the scheduling problem considered in this study.
First, the considered problem with urgent jobs can be classified as a multi-objective scheduling problem because two types of measures are combined, i.e., the makespan for normal jobs and total tardiness for urgent jobs. Recently, multi-objective scheduling problems are becoming popular as the manufacturing environments become more complex. Liang et al. [6] proposed a self-adaptive differential evolution algorithm for minimizing the maximum tardiness and makespan simultaneously in a flow shop with limited buffers. Anjana et al. [7] provided four metaheuristics for a flow shop that requires sequencedependent setup times, with the objective of minimizing makespan and mean tardiness together, while Rifai et al. [8] considered a distributed reentrant permutation flow shop with sequence-dependent setup times for minimizing makespan, tardiness and production costs. Recently, multi-objective flow shop problems derived from energy consumption issues were considered in [9][10][11]. They suggested metaheuristic algorithms to minimize total energy consumption and makespan simultaneously.
However, the combined objective of this problem is different from a typical multiobjective. In this study, the objective functions for the two classes are different because jobs are classified into two classes depending on their urgencies. On the contrary, in the typical multi-objective problems, two or more objectives are applied in common to all jobs. Therefore, this problem can be regarded as a multi-agent problem that is the special case of multi-objective problems [12]. In multi-agent problems, there are two or more agents, and each agent has a different objective function for its own jobs. If the different classes with different objectives are defined as different agents, this problem can be regarded as a multi-agent problem.
Most studies on multi-agent scheduling considered problems for a single machine. Baker et al. [13] developed algorithms for the multi-agent problem of a single system, considering the maximum delay, total weighted completion times and makespan for objective function. Agnetis et al. [12] proposed algorithms for getting constrained optimal and pareto optimal solutions for the multi-agent on single machine problems and proved that the problem is NP-hard. Ng et al. [14] studied the problem of two-agents on a single machine with the objective functions that the total completion time of the first agent is minimized and the number of tardy jobs of the second agent cannot exceed a predetermined number. Leung et al. [15] proved that the same problem is NP-hard even if the number of tardy jobs for the second agent is limited to zero. Cheng et al. [16] proposed an algorithm for a two-agent problem on a single machine with the objective of minimizing the total weighted number of tardy jobs of each agent, and they also studied a multi-agent single machine problem with the max-form objective functions in [17]. Liu et al. [18] studied a single machine with a two-agent problem in which each job has a linear deteriorated processing time. Additionally, polynomial time solution algorithms for a single-machine with two agents were suggested in [19,20].
For a flow shop with multi agents, there have been only a few studies. Agnetis et al. [12] studied a two-agent problem on a flowshop for minimizing makespan of jobs in each agent. Lee et al. [21] developed a metaheuristic algorithm and a branch and bound algorithm to solve a two-agent problem in a two-machine flowshop. In addition, Lee et al. [22] studied a flowshop problem with two agents where the objective of each agent is to minimize the number of tardy jobs and the total tardiness, respectively. Mor et al. [23] studied three different flowshop problems with two agents and developed polynomial time solution algorithms for each problem. Fan and Cheng [24] proposed a linear programmingbased approximation algorithm for a two-agent flowshop problem. Jeong and Shim [1] proposed a metaheuristic algorithm for a reentrant flowshop problem with two agents. Jeong et al. [25] presented a two-machine flowshop problem considering urgent jobs and developed metaheuristic algorithms and a branch and bound algorithm. Azerine et al. [26] considered a two-machine no-wait flow shop problem with two competing agents and proposed a branch and bound algorithm for small size problems and tabu search metaheuristics for large sized problems. Now, a survey on flow shop problems with limited waiting times is provided. Most studies on flow shops with limited waiting times considered two-machine problems. For the objective of minimizing the makespan, the two-machine flow shop problem was proved to be NP-hard in [5]. To find an optimal solution for the problem, several branch and bound algorithms were developed by [5,27,28]. In addition, the reversibility property for the scheduling problem was proved and a constructive heuristic algorithm based on an insertion mechanism was developed in [2]. On the other hand, there are also studies with different scheduling measures. Hamdi and Loukil [29] developed a heuristic algorithm and a Lagrangian relaxation-based lower bound strategy for minimizing total tardiness, while Dhouib et al. [30] proposed simulated annealing algorithms to hierarchically minimize the number of tardy jobs and the makespan for a flow shop.
In addition, there have been studies considering variant problems with waiting time limits. Kim and Lee [31] developed a branch and bound algorithm to minimize the makespan in a three-machine flow shop with overlapping waiting time constraints. Furthermore, An et al. [3] suggested a branch and bound algorithm and heuristic algorithms for a two-machine flow shop with the waiting time constraints and sequence-dependent setup times to minimize the makespan, while Lee [32] provided a genetic algorithm to minimize total tardiness for the same flow shop. In addition, a case in which several jobs can skip the first stage was considered in [33], and mathematical properties for the problem were proposed. In addition, for a problem with a batch machine followed by a discrete machine, a hybrid membrane computing metaheuristic algorithm was developed in [34] to minimize the makespan. Additionally, there are recent studies considering flow shop problems with various scheduling requirements as well as limited waiting times [35][36][37][38][39][40][41].

Problem Description and an MIP Model
This section describes the considered problem in more detail with assumptions made in this study and provides a mixed-integer programming (MIP) formulation. There are n independent jobs (i = 1, . . . , n) to be scheduled in a two-stage flow shop with limited waiting times between the stages. The n jobs belong to one of two classes, A or B, according to their urgencies. Two classes A and B represent each class of urgent jobs and normal jobs, respectively. Accordingly, let J A and J B be the sets of jobs in classes A and B, respectively. The objective function of this scheduling problem is the weighted sum of makespan for normal jobs and total tardiness for urgent jobs, which is to be minimized. For this problem, only permutation schedules are considered. In other words, jobs are processed on the two machines in the same order. Note that, although the permutation schedule does not guarantee optimality for the scheduling measures, semiconductor manufacturing systems process wafer lots in permutation schedules for ease of lot management and flexibility in material handling, and because of limited buffer spaces [3]. The following assumptions are also made this study: 1.
no job can be preempted; 2.
all normal jobs are available to be processed at time zero (beginning of the scheduling horizon); 4.
release times and arriving times of urgent jobs are positive and given in advance; 5.
the due date of an urgent job i is set to the earliest possible completion time assuming that the urgent job i is started immediately as arrived at the flow shop and processed without waiting, i.e., d i = r i + p i1 + p i2 . Table 1 represents the notation for parameters and decision variables used to describe the MIP and proposed heuristic algorithms throughout this paper.

Symbol Definition
i index for jobs (i = 1, 2, . . . , n) k index for machines (k = 1, 2) t index for positions in a sequence [t] index of the job in the tth position in a sequence p ik processing time of job i on machine k The following is the MIP for the scheduling problem considered in this study.
The objective function (1) is to minimize the weighted sum of total tardiness for urgent jobs and the makespan for normal jobs. Constraints (2) and (3) ensure that each job can be and should be assigned to only one position, and each position requires only one job. Constraint (4) ensures that jobs can be started after their release times. Constraints (5) amd (6) define the completion times of jobs assigned to tth position on machines 1 and 2, respectively. Constraints (7) and (8) ensure that jobs should keep their limited waiting times between machines 1 and 2. Constraints (9) and (10) define the com-pletion time of each job on machine 2. Constraints (11) and (12) define the makespan of normal jobs and the tardiness of urgent jobs. Constraints (13) and (14) define the domain of decision variables.

Heuristic Algorithms
The MIP can be used to find an optimal solution using a commercial optimization solver. However, because the considered scheduling problem is NP-hard, the solver may require a significant computation time to obtain an optimal solution for large or practicalsized problems. Therefore, this study focuses on proposing heuristic algorithms that can find a good solution close to optimal or acceptable. Two metaheuristic algorithms (an iterated greedy algorithm and a simulated annealing algorithm) are adopted in this study.
Before introducing the metaheuristic algorithms, we provide the following equations to calculate completion times of jobs in a sequence obtained in procedures of the proposed heuristics. Note that, because this study considers only permutation schedules, the completion times are calculated sequentially one by one from the first to the last. Equations (15) and (16) calculate the completion times of normal jobs and urgent jobs on the first machine, respectively, while Equation (17) calculates the completion times on the second machine. In addition, Equation (18) calculates the tardiness of urgent jobs. Here, it is assumed that , 0}, if the tth position job is an urgent job (18)

Initial Seed Sequence
In general, the metaheuristic algorithm starts with an initial seed sequence, and the performance of the metaheuristic depends on the seed sequence. Thus, the method to obtain the seed is introduced. In this study, a two-step method is used. A feasible permutation is generated by a list-scheduling in the first step. After that, the permutation is improved by a constructive heuristic, called NEH algorithm [42], in the second step.
List scheduling is a common method in practice because it is intuitive and very easy to implement. This method is triggered when a machine becomes available, and a job with the highest priority is assigned to the machine. Defining the priority for the jobs is the only requirement of this method. Here, a modified earliest due date rule (MEDD) is used to focus on the urgent jobs with the objective of minimizing their tardiness. For an unscheduled job i and a given partial schedule S, the modified due date (d i ) is defined as Because the due dates of normal jobs are set to an infinite number, normal jobs may be placed on rear positions in the MEDD schedule. To avoid this schedule, the modified due dates of normal jobs are set to the completion times assuming that they are scheduled immediately after S.
In step 2, the NEH begins with the MEDD schedule and makes schedules in a constructive way. In each iteration of the NEH procedure, a partial schedule is created by inserting the job at the front position of the MEDD schedule into the best position of the current partial schedule, and the inserted job is removed from the MEDD schedule. After generating a complete schedule, a pairwise interchange method is applied to the complete schedule for improving it. A detailed procedure is given in Figure 1.

Iterated Greedy Algorithm
The iterated greedy algorithm (IG) is a metaheuristic based on a stochastic local search strategy, developed in [43]. The original IG was specially devised for the permutation flow shop scheduling problem. Because of its simplicity of the structure and excellent performance, IGs have been popularly used in much flow shop research. The procedure of IG iterates four phases of destruction, construction, local search, and acceptance, to find better solutions. In the destruction phase, the predetermined number (d) of jobs are removed from a current schedule. Then, in the construction phase, a neighborhood schedule is generated by inserting the removed jobs in a constructive way such as that of NEH.

1.
In the proposed IG, d jobs are removed randomly from a current solution (S) in the destruction phase. Then, two partial schedules are generated. One (S D ) is a schedule consisting of d jobs removed from S, the other is a partial schedule (S R ) with (n − d) jobs. After that, in the construction phase, a complete schedule is constructed by inserting the jobs of S D into the best positions of S R , such as the for-loop in Figure 1. After these two phases, the complete schedule is improved by a local search procedure. In this IG, the interchange method used in NEH is implemented for the improvement.
Completing the local search, the new schedule (S') is compared to the current schedule (S). If S' is better than S, i.e., obj(S') < obj(S), S' replaces S. Otherwise, borrowing the concept of an acceptance in a simulated annealing algorithm, the replacement is accepted with a probability exp(-∆/T), where ∆ = obj(S') -obj(S), and T is an adjustable parameter (called temperature) in the IG. However, unlike SA, T is constant in IG. The temperature value T is calculated as suggested in [44] for the permutation flow shop scheduling problem. In this study, m = 2. These four phases are repeated until a termination condition is satisfied. Herein, the maximum computation time is used for the termination condition. The whole procedure of IG is summarized in Figure 2. In the procedure, let random() denote a function to return a random number greater than or equal to zero and less than one.

Simulated Annealing Algorithm
Simulated annealing (SA) is one of the most popular metaheuristic algorithms for combinatorial optimization problems including operations scheduling. In addition, SA showed good performance in a recent paper related to a two-machine flowshop with urgent jobs [1]. This study also adopts SA for the solution methodology. Generally, in flow shop scheduling problems, SA attempts to improve a current schedule (S) by making a small change in S. To accomplish this, insertion, which moves a randomly selected job to a randomly selected position, is used to generate a neighborhood solution (S'). If S' is better than S, S' replaces S. On the other hand, if S' is worse than S, S is replaced with S' with a probability exp(−∆/T). The initial temperature is obtained by Equation (19), and the temperature is gradually decreased by multiplying a cooling ratio γ if S cannot be improved for E iterations (E is called an epoch length.). The algorithm is terminated when the maximum computation time is elapsed. The whole procedure is given in Figure 3.

Computational Experiments
This section reports and analyzes the performance of the proposed algorithms through computational experiments on randomly generated problem instances. All the tested algorithms were coded in Java programming language, and the experiments were conducted on a personal computer with a 2.6 GHz CPU.
To generate problem instances, (urgent and normal) jobs and limited waiting times were generated by the methods used in [25] and [28], respectively. Processing times were randomly generated from a discrete uniform distribution with a range of [1, 100], i.e., U [1,100]. Let h be the proportion of urgent jobs to all jobs. For urgent jobs, three levels (0.3, 0.5 and 0.7) were considered for h, and three levels (0.5, 0.7 and 0.9) were considered for the weight for total tardiness of urgent jobs (α). The release times of urgent jobs were generated randomly from U[0, LB], where LB is a lower bound on the makespan of normal jobs. The lower bound is obtained by scheduling only normal jobs with Johnson's rule [45], assuming that their waiting times are infinite. Note that if all release times are zero or greater than LB, the problem can be solved more easily. In addition, limited waiting times were randomly generated from U[0, 100].
Prior to evaluating the proposed algorithms, calibration for the suggested two metaheuristic algorithms were performed to achieve better performance, because performance of metaheuristic algorithms generally depends on their parameters. For the calibration experiments, three levels of n = (50, 200, and 400) were considered, and three instances for each combination of (n, h, and α) were generated. For the termination condition, the maximum computation time was set to 50n milliseconds. Because the IG and SA are stochastic algorithms, they can find a good or bad solution coincidentally. To avoid such a coincidental solution, both algorithms solved each instance three times independently, and the average objective value was used. In addition, to compare the objective values obtained by the algorithms, the relative deviation index (RDI) was used as a measure, defined as (obj x − obj min )/(obj max − obj min ) where obj x , obj min and obj max are the objective function value from algorithm x, the minimum and the maximum, respectively.
The proposed IG has only one parameter which is the destruction size, d. To find the most appropriate value of d for the considered problem, eight IGs with different d = (4, 8, 10, 12, 14, 16, 18 and 20) were tested, respectively. Note that the original IG in [43] used d = 4, which showed the best performance in the computational experiment. The results are summarized in Figure 4 which shows the average RDIs. As shown in the figure, IG with d = 14 showed the lowest RDI, which means the best performance. Unlike the original IG, the RDI of IG with d = 4 was the highest among the test results. In addition, the bath-curve was observed, meaning that both lower and higher values from d = 14 deteriorate the performance and hence 14 is validated for the appropriate destruction size. Consequently, d = 14 was used in the subsequent experiments.  The performance of the proposed MIP was evaluated using a commercial optimization solver, CPLEX 12.10. For experiments, three levels of n = (10, 20, and 30) were considered, and five instances were generated for each combination of (n, h, and α). The maximum time limit for CPLEX was set to 3600 s to avoid an excessive computation time. The results were summarized in Table 2. From the table, CPLEX required a longer CPU time for the instances when h was increased and α was decreased. The reason is that as urgent jobs become more important, they must be scheduled closer to their arrival time, which reduces the number of sequencing candidates to consider. In addition, most of the problems with more than 20 jobs were not terminated within the maximum time limit. Considering that problem size with 20 jobs is not large enough, these results demonstrate the necessity to develop heuristic algorithms that perform well. We compared the performance of proposed algorithms with that of CPLEX. For experiments, heuristic algorithms solved the same instances used to evaluate the MIP. Let IG(x) and SA(x) denote the IG and SA algorithm, respectively, with the makespan nx milliseconds as a termination condition. Note that the destruction size for instances with n = 10 was set to 5 which is a half of n because using the determined value (d = 14) was impossible. The results are summarized in Table 3, which shows the average percentage errors (APE) of solutions from the proposed heuristics to those of CPLEX and the number of instances (out of 45) for which the algorithm found solutions better than or equal to those from CPLEX. If CPLEX failed to solve the problem to optimality within the time limit, the best solution obtained within the time limit was compared. In the table, negative APE indicates that the solutions from the heuristic were better than those from CPLEX. As can be seen from the table, NEH could significantly improve the MEDD solutions. In addition, because CPLEX could not solve most problems with more than 20 problems, IG and SA provided better solutions than CPLEX in the 20-30 size problems. In addition, IG showed better performance than SA. That is, the greedy search is more effective than a random search in small sized problems.  (27) measured by the relative percentage deviation (RPD), defined as (obj x − obj min )/obj min . The results are summarized in Table 4, which shows the average RPD. Overall, the proposed IG and SA outperformed by far TS from [26], except for only one case of (n = 50, h = 0.5, α = 0.9). Hence, we focused on the comparison between IG and SA. When the problem size is small and medium, IG found better solutions than SA, However, as the problem size increased, SA showed better performance than IG. This is probably because IG needs to check all insertion positions and this repeated procedure can be a significant computational burden for large size problems. In contrast, because SA generates a neighborhood solution with a simple insertion, the time to generate a new solution is almost the same for large size problems. According to the results, it can be said that IG is more efficient for a small size problem, whereas SA works better for a large size problem.
Based on the results, we plotted the interaction between algorithms and problem parameters (n, h and α) to check the effect on the performance of the algorithms. The interval plots are given in Figure 6, the 95% confidence interval for the mean RDIs. For the number of jobs n, as we stated, IG was better than SA with up to n = 200. However, SA outperformed IG in problem instances over n = 300. For the proportion (h) of urgent jobs to all jobs, SA showed better performance when h = 0.3 and 0.5. However, when h = 0.7, IG was superior to SA, and TS also showed good performance. In other words, three algorithms can be used in complement to one another when there are more urgent jobs than normal jobs. For the weight for total tardiness of urgent jobs (α), SA and IG far outperformed TS, and SA performed a little better than IG.
Finally, to see statistically significant differences in the performance of IG and SA, Kruskal-Wallis (KW) tests, which are non-parametric methods, were performed with a commercial statistical analysis problem SPSS. For these tests, instances were divided into two groups according to the number of jobs. This is because IG worked well on small sized instances, whereas SA performed better on large sized instances. Thus, one was a small group (n = 10, 20, 30, 50 and 100), and the other was a large group (n = 200, 300, 400, and 500). Since objective values have different scales according to n, RPD values among IG(50), IG(100), SA(50) and SA(100) were set as the dependent variable for the tests. To perform the tests, a significance level was set to 0.05. In addition, to check the effectiveness of problem parameters (h and α) on the performance of the proposed algorithms, KW tests were performed for (algorithms × h values) and (algorithms × α values) as well as algorithms.
For the small group, all p-values in the three KW tests were close to zero, i.e., less than the significance level. Detailed results were summarized in Table 5, which shows the results of pairwise comparisons from KW tests. The results confirmed a statistically significant difference between the performances of IG and SA. Thus, it can be said that IG is superior to SA for the small group. In addition, the performance of SA(50) and SA(100) was not different for the small group.   As in the KW tests for the small group, all p-values were close to zero and less than the significance level for the large group. Therefore, the results confirmed that there were statistically significant differences between algorithms. Table 6 shows the detailed test results. As stated earlier, p-values from pairwise comparisons between algorithms were less than 0.05, and hence it can be said that SA significantly outperformed IG. As shown in Figure 6, there were no significant differences between IG and SA when urgent jobs are more than normal jobs, i.e., h = 0.7, and the weights for two classes are equal, i.e., α = 0.05. Therefore, IG can be an alternative for these cases.

Discussion and Conclusions
This study considered a scheduling problem for a two-machine flow shop with urgent jobs and limited waiting times. The objective of this problem is minimizing the weighted sum of total tardiness of urgent jobs and the makespan of normal jobs. This is the first study considering urgent jobs and limited waiting time constraints together. We proposed a mixed integer programming for this problem. However, because this problem is NP-hard, it is very difficult to solve actual size problems with the MIP. Therefore, we developed two metaheuristic algorithms (an iterated greedy algorithm and a simulated annealing algorithm). Through a series of computational experiments, we suggested the best parameters for IG and SA. In addition, the effectiveness of IG and SA were verified by comparison to MIP. In addition, as the results showed, IG was efficient for small size problems, but SA showed superiority in large size problems. This problem can be extended in several directions in the future. For example, it is necessary to develop the lower bound to compare performance on large size problems. In addition, studies can investigate an optimal solution algorithm that is more efficient than MIP such as a branch and bound algorithm. If the proposed heuristic algorithms are appropriately modified, it is expected that this will become a standard for comparison of algorithms in future studies.