A new iterated greedy algorithm for no-idle permutation flowshop scheduling with the total tardiness criterion

A new iterated greedy algorithm for no-idle permutation flowshop scheduling with the total tardiness criterion

Journal Pre-proof

A new iterated greedy algorithm for no-idle permutation flowshop scheduling with the total tardiness criterion Vahid Riahi, Raymond Chiong, Yuli Zhang PII: DOI: Reference:

S0305-0548(19)30281-3 https://doi.org/10.1016/j.cor.2019.104839 CAOR 104839

To appear in:

Computers and Operations Research

Received date: Revised date: Accepted date:

15 May 2019 15 September 2019 4 November 2019

Please cite this article as: Vahid Riahi, Raymond Chiong, Yuli Zhang, A new iterated greedy algorithm for no-idle permutation flowshop scheduling with the total tardiness criterion, Computers and Operations Research (2019), doi: https://doi.org/10.1016/j.cor.2019.104839

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.

Highlights • We propose a new iterated greedy algorithm for no-idle flowshop scheduling. • It is one of the few studies to have considered the total tardiness criterion. • Evaluation is carried out using 120 benchmark instances x 3 due date scenarios. • Our algorithm is shown to be superior than other state-of-the-art algorithms. • Of the 360 benchmark instances tested, > 50% of the best solutions are updated.

1

A new iterated greedy algorithm for no-idle permutation flowshop scheduling with the total tardiness criterion Vahid Riahia,, Raymond Chiongb,∗, Yuli Zhangc, a

Institute for Integrated and Intelligent Systems (IIIS), Griffith University, Nathan, QLD 4111, Australia b School of Electrical Engineering and Computing, The University of Newcastle, Callaghan, NSW 2308, Australia c School of Management and Economics, Beijing Institute of Technology, Beijing, 100081, P.R. China; Sustainable Development Research Institute for Economy and Society of Beijing, Beijing, China

Abstract With the no-idle constraint, a machine has to process a job after finishing the previous one without any interruption. The start time of the first job on each machine must thus be delayed to meet this condition. In this paper, a new Iterated Greedy Algorithm (IGA) is presented for no-idle flowshop scheduling with the objective of minimizing the total tardiness. For the initialization phase, a variant of the NEH procedure is developed. Then, we propose a new variable local search based on an insert move with two different job selection mechanisms. A tardiness-guided job selection procedure, a job-dependent parameter and an insert-swap based method are further introduced in the destruction-construction phases. While most of the related studies have used a fixed probability for accepting new or non-improving solutions, we propose a time-dependent probability that allows our algorithm to focus on exploration in early iterations and exploitation in later iterations. Comprehensive computational experiments show that the proposed IGA is superior in terms of solution quality than state-of-the-art algorithms for the problem at hand. As a result, more than 50% of the existing best solutions for the benchmark instances tested have been updated. ∗

Corresponding Author at: University Drive, Callaghan, NSW 2308, Australia Email addresses: [email protected] (Vahid Riahi), [email protected] (Raymond Chiong ), [email protected] (Yuli Zhang )

Preprint submitted to xxx

November 5, 2019

Keywords: Flowshop Scheduling, No-idle Permutation, NEH, IGA 1. Introduction Flowshop scheduling has its applicability and practicality in many industrial problems. Developing efficient and effective scheduling techniques along this line of research is therefore highly important. For flowshop scheduling, a number of jobs must be processed on a certain number of machines in the same processing order. The aim is to find permutations of jobs such that a given objective is minimized. Among the different variants of flowshop problems, relatively fewer studies have considered the no-idle case. For the no-idle permutation flowshop scheduling problem (NPFSP), the machines must process all jobs without any interruption. In other words, no idle time is allowed for any machine after the first job in the sequence is started. This constraint is commonly seen in real-world environments such as fiber glass processing (Kalczynski and Kamburowski, 2005) and foundry production (Saadani et al., 2003). The NPFSP was first studied by Adiri and Pohoryles (1982). Given that exact methods such as branch and bound are able to tackle only smallscale instances of the problem (e.g., see Vachajitpan (1982) and Baptiste and Hguny (1997)), heuristic approaches have been seen as alternatives to solve problem instances of larger sizes. The first heuristic method for the NPFSP can be found in the work of Woollam (1986). By modeling the NPFSP as a traveling salesman problem, a heuristic approach based on the well-known nearest insertion rule was later proposed by Saadani et al. (2005). Another heuristic method was developed by Kalczynski and Kamburowski (2005), and they showed that their algorithm is better than that of Saadani et al. (2005). A greedy algorithm based on an insert move was studied by Baraz and Mosheiov (2008) and shown to be better than the approach of Saadani et al. (2005) too. For a comprehensive review of heuristic algorithms used for the NPFSP, see Goncharov and Sevastyanov (2009). The ability of metaheuristics to find near-optimal solutions in polynomial time for large-scale problem instances has prompted Pan and Wang to propose two hybrid algorithms based on discrete Particle Swarm Optimization (Pan and Wang, 2008a) and discrete Differential Evolution (DE) (Pan and Wang, 2008b). In their studies, a speed-up method was developed for the insert move. As a result, both their algorithms are greatly dependent on the 3

insert move. After that, an Iterated Greedy Algorithm (IGA) was presented by Ruiz et al. (2009). They tested the proposed algorithm on newly generated benchmark instances and showed that it is better than the algorithms of Pan and Wang. Two improved algorithms based on DE were later proposed by Deng and Gu (2012) and Tasgetiren et al. (2013b) for the NPFSP. Another metaheuristic algorithm – known as Invasive Weed Optimization (IWO) – was presented by Zhou et al. (2014), and their results showed that IWO is able to outperform the algorithm of Pan and Wang (2008a). More recently, Shao et al. (2017) developed a Memetic Algorithm and showed that their algorithm is significantly better than the algorithms of Pan and Wang (2008a), Ruiz et al. (2009), and Tasgetiren et al. (2013b). They updated 89 of the 250 best known solutions. Nagano et al. (2019) studied no-idle flowshops with total flowtime and proposed a highly efficient but simple constructive heuristic as well as a local search algorithm. Cheng et al. (2019) considered mixed no-idle flowshops with the makespan objective in a distributed environment, and proposed an IGA-based algorithm for their problem. It is worth pointing out that all of the aforementioned studies were focusing on the makespan or flowtime objective. In our work, we focus on minimizing the total tardiness, which is of great importance in manufacturing systems because completing a job after its due date in such an environment would lead to loss of contracts, customers and trust. Tasgetiren et al. considered the total tardiness objective and presented a DE (Tasgetiren et al., 2011) as well as a Discrete Artificial Bee Colony (DABC) algorithm (Tasgetiren et al., 2013a) to address it. For both algorithms, a speed-up method was proposed for the insert move. They compared their DE with other variants of DE (Tasgetiren et al., 2011) and the DABC algorithm to a Genetic Algorithm (GA) (Tasgetiren et al., 2013a). Shen et al. (2015) later proposed a Bi-population Estimation of Distribution Algorithm (BEDA) for the total tardiness NPFSP. Their results indicated that the BEDA is more efficient and effective than both the GA and DABC of Tasgetiren et al. (2013a). Recently, Shao et al. (2018) proposed a hybrid discrete teaching-learning based metaheuristic (HDTLM) for the problem, and they demonstrated the effectiveness of their algorithm using some well-known benchmark instances. In this paper, we present a new IGA to solve the NPFSP considering the total tardiness criterion. The IGA is a simple and powerful algorithm with just a few parameters to tune (Ruiz and St¨ utzle, 2007). It has been successfully applied to different variants of the flowshop scheduling problems, such as permutation flowshop scheduling (Ruiz and St¨ utzle, 2007; Karabulut, 4

2016; Dubois-Lacoste et al., 2017), blocking flowshop scheduling (Ribas et al., 2011; Ding et al., 2016; Tasgetiren et al., 2017; Riahi et al., 2019), no-idle flowshop scheduling (Tasgetiren et al., 2013b; Pan and Ruiz, 2014), distributed no-idle flowshop scheduling (Ying et al., 2017), no-wait flowshop scheduling (Pan and Wang, 2008a; Ding et al., 2015), and distributed flowshop scheduling (Fernandez-Viagas and Framinan, 2015; Ruiz et al., 2019). Motivated by its simplicity and efficiency, we incorporate the following novel/improved elements to our proposed IGA: I) Initial solutions: A new heuristic based on the NEH algorithm (Nawaz et al., 1983) is designed for the total tardiness NPFSP. Here, a new initial sequence is generated by taking both the total tardiness objective and no-idle constraint into account. II) Local search: A variable local search based on an iterative insert move with two possible job selection mechanisms is developed. The two selection mechanisms are random job selection and greedy job selection. For random job selection, jobs are selected randomly and inserted in any possible positions of the sequence. For the greedy selection, jobs are selected based on their positions in a reference permutation. The iterative nature of our local search means the list of jobs will be updated when an improvement is observed. III) Destruction-construction phases: Unlike the original IGA’s destructionconstruction phases (Ruiz and St¨ utzle, 2007), where jobs are randomly selected and removed from the current solution, and then reinserted in the best positions of a partial sequence, we introduce a new job-dependent parameter such that its value differs depending on the problem size. Additionally, instead of random selection a tardiness-guided procedure is proposed to select jobs. The idea behind this procedure is that jobs with higher tardiness values would have priority over jobs with lower tardiness values (i.e., to fix the more problematic parts of a solution). We also apply a swap move alongside the insert move during the construction phase to increase solution diversity. IV) Acceptance probability: To accept new or non-improving solutions, most existing studies used a Simulating Annealing (SA)-based acceptance criterion as per Ruiz and St¨ utzle (2007), where the probability remains the same at all time as the algorithm progresses. We propose a new time-dependent acceptance criterion that gives higher probabilities to non-improving solutions in early iterations. Then, when final iterations approach, this acceptance probability becomes close to zero. The basic idea is to explore the search space more in early iterations, and focus on exploitation in the last itera5

tions. This time-dependent acceptance function is not only computationally simpler, its value would also reduce directly as the algorithm progresses. To the best of our knowledge, our work represents one of the few studies in this area. Systematic computational experiments based on more than 100 benchmark problem instances from the literature with three different due date scenarios (> 300 instances in total) clearly show that our proposed algorithm is able to outperform the best-performing algorithms in the literature. As a result, over 50% of the existing best solutions for the benchmark instances used have been updated. The rest of this paper is structured as follows. In Section 2, we present the problem description. Details of the proposed algorithm are then given in Section 3. In Section 4, we discuss the setup of our computational evaluation as well as results obtained. Finally, we conclude the paper in Section 5, with possible future research directions highlighted. 2. Problem Description The NPFSP consists of n jobs (j = 1, 2, . . . , n) that should be processed on m machines (i = 1, 2, . . . , m) with the same sequence. p(j, i) denotes the time of processing job j on machine i (including the setup time), while dj denotes the due date of job j. At any time, each machine can process only one job, and each job can be processed on only one machine. To satisfy the no-idle constraint, machines must process jobs without any interruption from the start of processing the first job to the completion of the last job. Our formulation of the total tardiness NPFSP is based on that of Tasgetiren et al. (2013a). The purpose is to find a sequence of jobs to minimize the total tardiness criterion. Suppose π n = {π1 , π2 , . . . , πn } is a scheduling solution while π k = {π1 , π2 , . . . , πk } is a partial sequence, where πk refers to the job at position k. In addition, D(π k , i, i + 1) denotes the minimum difference between the completion of processing the last job of π k on machines i + 1 and i, which can be calculated as follows: D(π 1 , i, i + 1) = p(π1 , i + 1) ∀ i = 1, 2, . . . , m − 1,  D(π k , i, i + 1) = max D(π k−1 , i, i + 1) − p(πk , i), 0 + p(πk , i + 1) ∀ k = 2, 3, . . . , n, i = 1, 2, . . . , m − 1. 6

(1)

(2)

Then, the makespan of job πn on machine m, denoted as Cmax (πn ), can be calculated as follows: Cmax (πn ) =

m−1 X

n

D(π , i, i + 1) +

i=1

n X

p(πk , 1).

(3)

k=1

As there is no idle time between machines, after obtaining the makespan of the last job (πn ) on machine m, the makespan of other jobs can be calculated as follows: Cmax (πk ) = Cmax (πk+1 ) − p(πk+1 , m) ∀ k = n − 1, . . . 1.

(4)

Finally, the total tardiness, denoted as T T , can be computed as follows: TT =

n X k=1

(max(Cmax (πk ) − dπk , 0)) .

(5)

3. The Proposed IGA The IGA is a simple and effective algorithm proposed by Ruiz and St¨ utzle (2007) for flowshop scheduling. It has a number of main steps. At first, the algorithm begins with an initial solution, generated either randomly or through a constructive heuristic, and then it goes through a loop until the stopping criterion is reached. The loop consists of two main phases: destruction and construction. In the destruction phase, some jobs are removed from the current sequence. In the construction phase, the removed jobs are reinserted into the remaining partial sequence by a greedy constructive method to create a new and complete sequence. After that, an acceptance criterion is applied to decide if the new solution should replace the incumbent. Although local search is an optional component for the algorithm, it has been demonstrated that a local search-based IGA can produce better results than one without local search (Ruiz and St¨ utzle, 2007). In this work, we propose an IGA incorporated with some novel/improved elements. For simplicity’s sake, we called it the Enhanced IGA (EIGA). Our EIGA uses a new NEH heuristic to generate the initial solution, and then applies three steps – a tardiness-guided destruction procedure with jobdependent selection and a swap-insert move during the construction phase (applied as perturbations), a variable insertion-based local search before and after destruction-construction (to enhance the quality of solutions), and a 7

new time-dependent acceptance probability during the decision-making stage (to decide if a solution is kept) – iteratively. Details of these are described in the following subsections. 3.1. Initial Solutions A new heuristic based on the NEH procedure (Nawaz et al., 1983), denoted as NNEH, is designed to generate the initial sequence. With the original NEH algorithm, the initial order of jobs is obtained by sorting jobs in descending order according to the total processing time. After that, the best two-job partial sequence is chosen by selecting two of the first jobs from the initial order. Then, other unscheduled jobs from this initial order are picked to construct a complete sequence by inserting them at all possible positions of the current best partial sequence. The basic idea of NEH is to give priority to jobs with higher total processing times than jobs with less total processing times. For the NPFSP, however, early processing of jobs with smaller total processing times is preferable for no-idle machines. So, variants of NEH such as Rajendran’s algorithm (Rajendran, 1993), known as Raj in short, are more efficient for addressing the NPFSP. With the Raj algorithm, jobs with smaller processing times on earlier machines are given more priority than jobs with larger processing times on earlier machines. This is an appropriate assumption based on the Weighted Total Processing Times (WTPT) characteristic of the NPFSP. On the other hand, given the total tardiness objective, an initial order based on the Earliest Due Date (EDD) (Kim, 1993) is desirable. With the EDD, the initial sequence is obtained by arranging jobs in increasing order of their due dates. We decide that a new initial order should be extended to include both Raj’s initial order as well as the EDD’s. The jobs will therefore be sorted based on non-decreasing of αW T P Tj + (1 − α)dj , where m X W T P Tj = (m − i + 1) × p(j, i)

(6)

i=1

Here, dj is the due date of job j, and α is a parameter with values between 0 and 1. Details of the NNEH heuristic can be found in Algorithm 1. 3.2. Local Search Our local search is based on the insert move: first a job is removed from its position in the permutation process, and then reinserted in another ran8

Algorithm 1: NNEH Arrange all n jobs in non-decreasing order of αW T P Tj + (1 − α)dj ; th 2 For each k ∈ [2, n], insert the k job of the initial order one by one to all positions of the best partial sequence. Then, find the best sequence for next iteration;

1

dom position. The insert move is selected because of the proposed speed-up method by Tasgetiren et al. (2013a), which reduces computational complexity of the whole insertion neighborhood in a sequence from O(n3 m) to O(n2 m). Specifically, we use an iterative local search (iLS), which means when an improved solution is found the procedure is repeated again. Here, the job selection process can be either random or greedy, i.e., random iterative local search (riLS) or greedy iterative local search giLS. For riLS, at each iteration one of the jobs is selected randomly with no repetition, and it is inserted in any possible positions of current solution π. If the insertion of this job does not result in a better solution, a different job is selected. Details can be seen in Algorithm 2. Algorithm 2: riLS Let π be the current solution and π R the reference permutation that is constructed from a random permutation. 2 for k = 1 to n do 3 Find job πkR in π and remove it from the sequence of π. 4 Find the best permutation, π 0 , by inserting job πkR in any other positions of π. 5 If T T (π 0 ) < T T (π), π = π 0 , and go to Step 1. 6 return π as the solution after local search. 1

For giLS, instead of random job selection, jobs are extracted based on an order given by a reference permutation, π R , which is the best solution found so far in the search (Pan et al., 2008). As an example, let π R = {3, 4, 2, 1} and π = {2, 3, 1, 4} be the reference and current permutations, respectively. First, job 3 is selected from the current solution, and then jobs 4, 2, and 1 in the given order. We use both riLS and giLS in our proposed algorithm. At each iteration, the algorithm draws at random a number between 0 and 1. If this number 9

exceeds a threshold, denoted as P LS, riLS is adopted as local search. Otherwise, giLS is employed as local search. The best solution found so far is used as the reference solution for giLS. 3.3. Tardiness-guided Destruction-Construction The standard destruction-construction phases (Ruiz and St¨ utzle, 2007) work as follows. First, d jobs are randomly selected from the current solution (the destruction phase). Then, those jobs will be reinserted in the partial permutation, one by one, similar to the insertion procedure of NEH, to obtain a new complete sequence (the construction phase). This process has been adopted by almost all of the IGAs presented in the literature (e.g., see Ribas et al. (2011), Tasgetiren et al. (2013b), Pan and Ruiz (2014) and Pan and Wang (2008a)). We attempt to improve the standard destruction-construction phases described above by asking the following questions: Which jobs should be selected? How many of them should be selected? How selected jobs should be re-inserted in the partial sequence? One of the most important issues here is to decide which jobs should be selected. In the literature, jobs are randomly selected from the current solution, and then placed in the best possible spot of the partial sequence. We propose a strategy based on the tardiness objective, referred as Tardinessguided Destruction-Construction (TgDC ). Instead of random selection, we select jobs with the highest tardiness values. First, jobs are sorted decreasingly based on their tardiness values and saved in list ϕ. Next, the first d/2 jobs of list ϕ are chosen and inserted into π d . The remaining jobs are then selected at random to guarantee a different solution at each iteration. The basic idea behind our TgDC is that at each permutation, very likely some of the jobs are early while others are tardy. Among those tardy jobs, some of them are more crucial than others, especially those with higher tardiness values. These crucial jobs are given the priority to be selected. In other words, TgDC tries to fix the most problematic parts of the current solution rather than selecting jobs at random, which may lead to the possibility of destroying the proper parts of a solution. Another important issue is the number of jobs that should be removed in the destruction phase. Note that the choice of d here is vital. When the d value is small, a new solution is similar to the current solution and the algorithm may not be able to escape from local optima; whereas a large d value is no different from a randomized NEH procedure (Pan and Ruiz, 10

2014). In most of the related studies, the value of d remains the same for every problem instance despite the different job sizes, e.g., 6 in Ding et al. (2015), 7 in Tasgetiren et al. (2017), 8 in Pan and Ruiz (2012), and 10 in Pan and Ruiz (2014) for problem instances with sizes of 20, 50, 100, 200, and 500. We opt for a job-dependent d as follows: ( 3 if n ≤ 50, d= (7) n × β if n > 50. where n is the number of jobs and β is an algorithm parameter. In Equation 7, we divide the d value into two parts, i.e., when n is less than or equal to 50 (small and medium sizes) and when n is greater than 50 (large sizes). By doing this, the value of d directly depends on the number of jobs when it comes to large-scale instances. It is worth noting that the value of d=3 for small and medium-size problems was determined based on some preliminary analysis. Last but not least is the question of how to re-insert selected jobs in the construction phase? In addition to using the standard insert move, we propose a swap move alongside it. With the swap move, two jobs are selected in either a random or greedy manner from the current sequence and then, they exchange their positions with each other. Our insert and swap moves work as follows. First, a number between 0 and 1 is drawn randomly. If this number is less than or equal to 0.5, then the insert move is applied. Otherwise, the swap move is used. Our motivation for introducing the swap move is that most of the related studies on no-idle flowshop scheduling have devoted their attention to speed up the insert move (e.g., see Pan and Wang (2008b) and Tasgetiren et al. (2013a)). The insert move alone exposes the algorithm to a higher risk of being trapped at local optima. In our attempt to address this issue, we apply both insert and swap to increase diversity of the solutions. To better understand the details, consider the following example with 7 jobs and 5 machines. Their processing times and due dates are given in Table 1. Suppose that the initial sequence before destruction-construction is π = (2, 3, 6, 4, 5, 1, 7), with the objective value of 1435. The jobs will be sorted decreasingly based on their tardiness values. Tardiness values for jobs in the current permutation are 245, 51, 93, 296, 255, 260, and 235 for jobs 2, 3, 6, 4, 5, 1, and 7, respectively. So, ϕ = (4, 1, 5, 2, 7, 6, 3). Since the number 11

Table 1: An example of 7 jobs and 5 machines with their processing times and due dates.

Job j 1 2 3 4 5 6 7

1 1 48 18 2 98 38 26

2 56 35 86 10 45 53 66

Machine i 3 4 20 81 89 82 71 51 37 15 12 1 57 60 41 38

5 58 74 31 17 1 61 86

Job due dates 290 137 362 195 237 381 401

of jobs is less than 50, 3 jobs will be selected. The first two jobs of list ϕ are first selected, and then the third one is selected at random. Assume the randomly selected job is 5, we have π d = (4, 1, 5) for the construction phase. Next, three random numbers will be generated – assume that they are 0.37, 0.62, 0.51. These random numbers indicate that the first job, job 4, must be reinserted and the other two jobs, jobs 1 and 5, must be swapped in the remaining partial sequence. The first job will be reinserted in the partial sequence as per Table 2. Table 2: Current and candidate solutions with their total tardiness after job 4 is reinserted.

current solution (2,3,6,7,4)

candidate solutions (2,3,6,4,7) (2,3,4,6,7) (2,4,3,6,7) (4,2,3,6,7)

total tardiness 690 621 652 638 441

Based on Table 2, the new partial solution (4,2,3,6,7) will be selected for the next step, in which job 1 is swapped as per Table 3. From Table 3, the sequence (4,2,3,6,7,1) is selected for the last step. According to Table 4, after swapping job 5, the solution (4,2,3,6,7,5,1) is the new solution for the next phase. 12

Table 3: Current and candidate solutions with their total tardiness after job 1 is swapped.

current solution (4,2,3,6,7,1)

candidate solutions (4,2,3,6,1,7) (4,2,3,1,7,6) (4,2,1,6,7,3) (4,1,3,6,7,2) (1,2,3,6,7,4)

total tardiness 731 949 1073 1136 787 735

Table 4: Current and candidate solutions with their total tardiness after job 5 is swapped.

current solution (4,2,3,6,7,1,5)

candidate solutions (4,2,3,6,7,5,1) (4,2,3,6,5,1,7) (4,2,3,5,7,1,6) (4,2,5,6,7,1,3) (4,5,3,6,7,1,2) (5,2,3,6,7,1,4)

total tardiness 1245 1244 1418 1338 1937 1706 1758

3.4. Time-dependent Acceptance Probability Acceptance probability is one of the important considerations of local search-based algorithms, since it helps algorithms to escape from local optima by accepting non-improving solutions. With the IGA, at each iteration, it compares the quality of a new solution generated by the local search to the current one based on the acceptance criterion and decides whether to accept or reject the new solution for the next iteration. In the literature, two methods are commonly seen: an SA-based acceptance criterion and a random-based acceptance criterion. • SA-based acceptance criterion: Most previous studies have used the SA-based acceptance criterion presented in Ruiz and St¨ utzle (2007). According to this acceptance criterion, a new solution π 0 is accepted with a probability of 13

    T T (π) − T T (π 0 ) ,1 P = min exp T

(8)

where π is the current solution and the temperature, T , is determined as follows: T = where T0 is a parameter.

Pn

j=1

Pn

i=1

pij

10.m.n

× T0

(9)

This acceptance criterion has some drawbacks. With the original SA algorithm, the value of T is typically high during the initial iterations, and then it will be reduced such that most solutions will be accepted initially, but as the T value approaches zero, most non-improving solutions Pn Pn will be rejected. Let us consider a simple example, based on j=1 i=1 pij ≈ 5 and T0 = 0.5 (used in most of the published papers), 10.m.n T ≈ 2.5. With this value of T , when T T (π) − T T (π 0 ) = −20, the probability of accepting a new solution is only 0.03 %. For problem instances of medium and large sizes, however, a solution with a difference of 20 is likely to be a good solution and should have been given a higher probability to be accepted especially at the earlier iterations. • Random-based acceptance criterion: Another acceptance criterion commonly used in published IGAs is a random mechanism (Ribas et al., 2011). As the name suggests, non-improving solutions will be accepted with a probability of P = 0.5. Contrary to the SA-based acceptance criterion, this random mechanism accepts almost half of the non-improving solutions. This may help an algorithm to escape from local optima, however, having the same probability of 0.5 throughout the search process is not ideal. To overcome the above limitations, we propose a time-dependent acceptance criterion as follows: P =

(M SC − CSC) M SC × γ

(10)

where M SC is the maximum run time, CSC is the time for local search to complete, and γ is the algorithm parameter. This time-dependent procedure 14

allows our algorithm to focus more on exploration in earlier iterations and exploitation in later iterations. Additionally, it is computationally simpler than the SA-type acceptance criterion, because it only needs a random number generation, while the SA-based one requires an exponential function (Talbi, 2009). 3.5. Overview of the EIGA Taken as a whole, the EIGA is shown in Algorithm 3. It should be noted that π ∗ is the best solution found so far, and the stopping criterion is based on the maximum CPU time (see Section 4 for details). Algorithm 3: The EIGA Set the parameters, d and β, and best solution so far π ∗ = ∞. 2 Establish the initialized solution, π ← NNEH(π). 0 3 Apply local search to the initial solution, π ←LocalSearch(π). 0 ∗ 0 4 If T T (π ) < T T (π), π = π = π . 5 Apply the destruction-construction procedure to π, π 0 ← destruction − construction(π). 0 00 0 6 Apply local search to π , π ←LocalSearch(π ). 00 ∗ ∗ 00 7 If T T (π ) < T T (π ), then π = π = π . Otherwise, if 00 T T (π ) < T T (π) or rand[0, 1] ≤ P (equation 10), π = π 00 . 8 If the stopping criterion is met, return the best solution found so far; otherwise go to Step 5. 1

4. Computational Experiments and Results In this section, we report on comprehensive computational experiments conducted to validate the effectiveness of the proposed EIGA. Experimental settings are first given in Section 4.1, and parameter tuning is discussed in Section 4.2. A comparison between our NNEH and the original NEH (Nawaz et al., 1983) is then presented in Section 4.3. After that, we compare the EIGA with existing IGAs in the literature, including IGARS1 (Ruiz and St¨ utzle, 2007), IGARS2 (Ruiz and St¨ utzle, 2008), IGAP T L (Pan et al., 2008), and IGAP R (Pan and Ruiz, 2014), in Section 4.4. Finally, results obtained by the EIGA are compared to state-of-the-art algorithms designed for NPFSPs with the total tardiness criterion, namely the 15

GA and DABC of Tasgetiren et al. (2013a) and BEDA of Shen et al. (2015), in Section 4.5. It is worth noting that we have not included the HDTLM of Shao et al. (2018) in our comparison, as we failed to replicate their algorithm. 4.1. Experimental Settings All algorithms were implemented using the C programming language with the same random seeds and termination criterion, i.e., the same maximum elapsed CPU time limit of Tmax = ρn m2 milliseconds, where ρ was tested with three values: 20, 40 and 60. These algorithms were run on the same computing cluster named Gowonda at Griffith University, Australia. In addition, all statistical tests were carried out with minitab 18.1. We used the well-known Taillard’s benchmark set (Taillard, 1993) for our computational experiments. This benchmark set has 120 instances with 12 groups of different problem sizes ranging from 20 jobs with 5 machines to 500 jobs with 20 machines. Each group of instances contains 10 test problems. The due date of a job is determined by the total work content (TWC) rule Pm (Brah, P 1996), i.e., dj = τ × k=1 Pjk , where τ is a due date tightness factor and m k=1 Pjk is the total processing time of job j on all the machines. As per Tasgetiren et al. (2013a), we set τ to 1, 2 and 3 to have tight, medium and loose job due dates, respectively. To compare the performances of the algorithms, we have selected the relative percentage deviation (RPD) as a measure, calculated as follows: RP D =

(T T A − T T B ) × 100 TTB

(11)

where T T A is the total tardiness of the solution obtained by algorithm A, and T T B is the reference total tardiness value. A smaller RPD value indicates better algorithm performance. It should be noted that, for fairness’ sake, the used acceleration method for the insertion operator was employed in all algorithms being compared. The proposed constructive heuristic was also used in all the tested algorithms. For algorithms proposed for the problem at hand from the literature, their parameter values were taken directly from the respective papers. For the IGA variants tested, we performed parameter tuning using the same experimental design as per the parameter tuning of our proposed EIGA, to find the best values for their parameters. All the four implemented IGA variants have two parameters: T0 and d. As suggested by Ruiz and St¨ utzle (2007) and later 16

confirmed by Pan and Ruiz (2014) in their study on the NPFSP, T0 is very robust and basically any value except 0 would be acceptable. We therefore set the T0 value for each IGA variant based on their respective paper. As for d, after calibration the following values were used: 4 for IGARS1 (Ruiz and St¨ utzle, 2007), IGARS2 (Ruiz and St¨ utzle, 2008), and IGAP T L (Pan et al., 2008), and 6 for IGAP R (Pan and Ruiz, 2014). 4.2. Parameter Tuning and Calibration We calibrated three parameters of the proposed EIGA: β in the destructionconstruction phases, γ for the acceptance criterion, and P LS for the local search. We also tested the idea of using tardiness guided job selection in the destruction phase, the new time-dependent acceptance criterion employing both insert and swap in the construction phase, and the variable jobdependent d value in the destruction phase. To be specific, the following combinations of factors were tested: • β: 0.02, 0.03, 0.04, 0.05. • γ : 1, 2, 3, 4. • P LS: 0.0, 0.2, 0.5, 0.8, 1.0. • Job selection in the destruction phase (DCJobSelection in short): – TGDC (the proposed Tardiness Guided Job Selection in the destruction phase) – RBDC (common Random Based Job Selection in the destruction phase (Ruiz and St¨ utzle, 2007)) • Acceptance criteria (Acceptance in short): – RBAC (the proposed random-based acceptance criterion) – SABAC (the SA-based acceptance criterion (Ruiz and St¨ utzle, 2007)) • d values: – Cd (a constant d value for all problem sizes in the destruction phase, i.e., d = 4 (Ruiz and St¨ utzle, 2007)) – Vd (the proposed variable d value) 17

• Construction-destruction move strategies (CD move in short): – Ins (only insert as per the original (Ruiz and St¨ utzle, 2007)) – Swap (only the swap move instead of insert) – InsSwap (both swap and insert) These resulted in a total of 4 × 4 × 5 × 2 × 2 × 2 × 3 = 1920 different configurations. To avoid biased results, the benchmark problem set used for parameter calibration was different from that used for algorithm comparisons. Here, we randomly generated 30 basic instances (n ∈ {20, 50, 100, 200, 500}, m ∈ {10, 20}, and 3 instances for each n × m combination) and 3 due date scenarios with the TWC rule, resulting in a total of 90 instances. Additionally, to have more reliable results, each combination was tested 5 times with the specific instances. The maximum elapsed CPU time was set as Tmax = 20n m2 milliseconds. The average RPD was calculated based on Equation (11), where T T B was the minimum total tardiness obtained (only in this series of experiments). As shown in Table 5, we used Analysis of Variance (ANOVA) to analyze the experimental results. From the table, we see that all factors except P LS have high F-ratios and are statistically significant, with p-values close to zero. Job selection in the destruction phase (DCJobSelection) has the highest Fratio, which means it has a huge effect on the response variable, i.e., the proposed algorithm. In addition, we also show the 95% confidence interval plots of the factors in Figure 1. Here, non-overlapping confidence intervals of each pair level indicate significant differences. From Figure 1, we see that level 0.04 of β, the RBAC level of acceptance criterion, the Vd level of d value, and the TGDC level of DCJobSelection have statistically better results than the other tested levels. For factor γ, levels 2, 3, and 4 are statistically equivalent but level 2 has a lower average. For factor P LS, the levels are also statistically equivalent, with level 0.5 having a smaller average. For factor CD move, levels Ins and InsSwap yield statistically better results compared to Swap, but there is no significant difference between Ins and InsSwap. Nevertheless, InsSwap is the preferred choice because of a lower average. Based on these results, the following levels have been selected: β = 0.04, γ = 2, P LS = 0.5, DCJobSelection = TGDC, Acceptance = RBAC, d value = Vd, and CD move = InsSwap. It is worth noting that when an interaction between two factors in the ANOVA table is significant, the interaction plot needs to be checked. Due to 18

space constraints, however, we show just a few of them: β ∗ CD move, γ ∗ CD move, CD move∗γ, d value∗Acceptance, and d value∗DCJobSelection. They can be found in Figure 2. Table 5: ANOVA test results for parameter calibration. Source β γ P LS DCJobSelection Acceptance d value CD move β∗γ β ∗ P LS β ∗ DCJobSelection β ∗ d value β ∗ CD move γ ∗ P LS γ ∗ DCJobSelection γ ∗ Acceptance γ ∗ d value γ ∗ CD move P LS ∗ DCJobSelection P LS ∗ Acceptance P LS ∗ d value P LS ∗ CD move DCJobSelection ∗ Acceptance DCJobSelection ∗ d value DCJobSelection ∗ CD move Acceptance ∗ d value Acceptance ∗ CD move d value ∗ CD move Error T otal

Degree of freedom Sum of square Mean square F-ratio p-value 3 3 4 1 1 1 2 9 12 3 3 6 12 3 3 3 6 4 4 4 8 1 1 1 1 2 2 9492 9599

16.07 35.54 1.36 16.65 8.56 11.95 4.90 25.85 4.48 3.34 5.52 2.15 3.24 3.36 3.35 4.92 3.57 0.96 2.98 3.23 1.94 9.53 7.45 5.02 8.97 6.38 6.19 1282.67 1495.34

5.3574 11.8475 0.3389 16.6461 8.5644 11.9494 2.4516 2.8727 0.3730 1.1119 1.8389 0.3587 0.2702 1.1211 1.1150 1.640 0.5943 0.2390 0.7460 0.8067 0.2422 9.5287 7.4509 2.5111 8.9859 3.1902 3.0946 0.1351

39.65 87.67 2.51 123.18 63.38 88.43 18.14 21.26 2.76 8.23 13.61 2.65 2.00 8.25 8.25 12.14 4.40 1.77 5.52 5.97 1.79 70.51 55.14 18.58 66.50 23.61 22.90

0.0000 0.0000 0.061 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0000 0.0000 0.014 0.021 0.0000 0.0000 0.0000 0.0000 0.123 0.0000 0.0000 0.073 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

4.3. NNEH as a New Heuristic Next, we carried out evaluation on the proposed NNEH of our algorithm. We would like to present it as a new heuristic that other researchers can adopt 19

7 (a)

(b)

(c)

(d)

(e)

(f)

20 (g) Figure 1: 95% confidence interval plots of parameters.

(a)

(b)

(c)

(d)

(e) Figure 2: Interaction plots of CD move versus β, γ versus β, CD move versus γ, d value versus Acceptance, and d value versus DCJobSelection (confidence intervals omitted for clarity’s sake).

21

for the problem at hand, and therefore we compared it to NEH and other variants here. The evaluation was done using all the 120 problem instances from Taillard (1993) with 11 different values of α. We calculated the RPD of each α based on Equation (11), where T T B is the best total tardiness value found by NNEH with different α. The results are shown in Figure 3. Note that when α = 0 and α = 1, the algorithm basically resembles the EDD (Kim, 1993) and Raj (Rajendran, 1993) algorithms, respectively. As can be seen in Figure 3, NNEH with α = 0.2 has the best performance. From the figure, we also observe that the best results for the algorithm are obtained from α = 0.1 to α = 0.3. We decided to carry out additional analysis with this range, based on values from α = 0.11 to α = 0.29 (see the dots highlighted in red in Figure 3). The additional analysis shows that α = 0.23 leads to the best performance of NNEH. It is worth pointing out that the proposed heuristic has obtained much better results compared to α = 1, which is the EDD algorithm.

Figure 3: NNEH with different values of α

Results of a direct comparison between the performances of our proposed NNEH and the classical NEH algorithm on problem instances from Taillard (1993) can be found in Table 6. From Table 6, we note that NNEH has better performances than NEH for most of the problem instances. This is true not only with increasing sizes of the instances but also with looser due dates. Statistical analysis based on simple t-tests (instead of ANOVA, 22

Table 6: Comparison of NEH and NNEH. n×m 20 × 5 20 × 10 20 × 20 50 × 5 50 × 10 50 × 20 100 × 5 100 × 10 100 × 20 200 × 10 200 × 20 500 × 20 Average

τ =1 NEH NNEH 4.01 2.31 3.43 4.19 3.69 2.10 7.43 0.00 5.79 0.45 0.39 1.48 8.88 0.00 4.15 0.32 0.58 1.04 6.15 0.00 2.74 0.02 4.07 0.00 4.28 0.99

Average makespan τ =2 NEH NNEH 5.10 4.08 5.93 2.99 6.69 3.26 8.98 0.00 5.50 2.00 1.25 2.71 9.16 0.00 4.15 0.68 1.07 1.25 6.03 0.00 3.25 0.17 4.32 0.00 5.12 1.43

NEH 10.59 14.09 12.55 9.11 7.07 2.83 9.83 4.95 2.32 7.33 3.07 4.77 7.37

τ =3 NNEH 2.43 5.23 5.39 0.00 1.90 5.00 0.00 0.00 1.04 0.00 0.01 0.05 1.75

as there are only two algorithms in comparison) with a significance level of 0.05 confirmed that the differences between results obtained by NNEH and NEH are significant in all three due date scenarios (p = 0.00 in all three scenarios). In addition, the 95% Tukey’s HSD intervals, shown in Figure 4, further confirmed that differences between the results are significant in all cases. However, it is worth highlighting that, as this entire evaluation was done based on Taillard’s instances, the results here basically mean that NNEH is better than NEH only on these problem instances and the findings cannot be generalized. 4.4. Comparing the EIGA with Existing IGAs Subsequently, we carried out experiments to compare the proposed EIGA with existing IGAs. To this end, four IGAs were re-coded: IGARS1 (Ruiz and St¨ utzle, 2007), IGARS2 (Ruiz and St¨ utzle, 2008), IGAP T L (Pan et al., 2008), and IGAP R (Pan and Ruiz, 2014). A brief description about these algorithms are given as follows: • IGARS1 (Ruiz and St¨ utzle, 2007). This algorithm is the original IGA proposed for the permutation flowshop scheduling problem. First, the well-known NEH algorithm is employed to obtain an initial solution. In the destruction-construction phases, d jobs are randomly removed from the current solution, and then they are reinserted, one by one, 23

Figure 4: 95% Tukey’s HSD intervals for NNEH and NEH on all three due date scenarios

in the best position of the partial permutation. An SA-based function with a fixed T value is used as the acceptance criterion. Finally, an iLS with random job selection is applied as the local search method. • IGARS2 (Ruiz and St¨ utzle, 2008). This algorithm is developed for the sequence dependent setup times flowshop problem with both makespan and weighted tardiness objectives. This algorithm is similar to IGARS1 , with only one major difference. That is, N EHEDD is used to generate the initial solution when it comes to the tardiness objective. • IGAP T L (Pan et al., 2008). This algorithm is developed for the flowshop problem with total flowtime objective. This algorithm is also similar to IGARS1 but with one major difference. That is, a greedy local search (not iterative) with the NEH as the reference permutation is applied as the local search method. • IGAP R (Pan and Ruiz, 2014). This algorithm is proposed to deal with the mixed no-idle flowshop scheduling problem. To generate the initial n solution, a heuristic algorithm named LR( m ) proposed by Liu and Reeves (2001) is used. giLS is employed as the local search method. In the destruction-construction phases, a minor but useful modification is proposed, in which after reinserting one job, the jobs occupying the previous and posterior positions are also reinserted in all positions of the current partial sequence. The acceptance criterion is the same as the original version. 24

It is worth mentioning that these four algorithms have been selected because the first three, IGARS1 , IGARS2 , and IGAP T L , are the most cited ones to date, and the last one, IGAP R , was proposed for a variant of the NPFSP. As can be seen, most of these existing IGAs are very similar compared to the original version. Our algorithm differs quite substantially from the original IGA, as we have introduced a new heuristic component, a time-dependent acceptance criterion, and tardiness-guided destruction-construction with both the swap-insert move and job-dependent d value. We tested these algorithms using the Taillard’s problem instances (Taillard, 1993) with three different due date parameters. All algorithms used the same stopping criterion: Tmax = 20×n×( m2 ) milliseconds. For every problem instance, each algorithm was repeated for five times. The results, averaged over the five runs, are reported in Table 7, where T T B is the minimum total tardiness known for each instance. From this table, it is obvious that the proposed EIGA is superior compared to other existing IGAs, especially for large instances. Among the four existing algorithms, the performances of IGARS1 , IGARS2 and IGARS2 are almost identical owing to the fact that they are very similar to each other. The results of IGAP R are better than the other three. Table 7: Comparing the EIGA with existing IGAs.

20*5 20*10 20*20 50*5 50*10 50*20 100*5 100*10 100*20 200*10 200*20 500*20 Average

IGA RS1 0.06 0.02 0.00 0.62 0.98 0.96 0.76 1.09 1.24 1.18 0.86 1.14 0.74

IGA RS2 0.01 0.00 0.00 0.61 0.99 0.93 0.69 1.02 1.13 1.19 0.79 1.06 0.70

IGA PTL 0.00 0.01 0.00 0.59 0.93 0.98 0.78 1.05 1.21 1.23 0.89 1.03 0.73

IGA PR 0.01 0.00 0.00 0.53 0.85 0.75 0.63 0.82 0.91 0.95 0.70 0.96 0.59

EIGA 0.00 0.00 0.00 0.15 0.25 0.27 0.30 0.25 0.41 0.31 0.28 0.32 0.21

We also performed an ANOVA test to check the statistical differences of the results by considering each algorithm as a factor. The test results are given in Table 8. At a confidence level of α = 0.05, we see that there are significant differences among the algorithms compared. An interval plot is 25

shown in Figure 5. From this figure, it is clear that the proposed algorithm is able to significantly outperform the other algorithms. Table 8: ANOVA test results for the existing IGAs in comparison. Source

Degree of freedom

Sum of square

Mean square

F-statistic

p-value

Algorithms Error T otal

4 595 599

23.64 88.45 112.09

5.9095 0.1487

39.75

0.0000

Figure 5: An interval plot of the EIGA and other IGAs

4.5. Comparing the EIGA to State-of-the-art Algorithms Finally, we compared the proposed EIGA to four state-of-the-art algorithms designed for the same problem, namely the GA and DABC from Tasgetiren et al. (2013a), BEDA from Shen et al. (2015), and IGAP R from Pan and Ruiz (2014), under three different due date scenarios. We carefully re-coded all the algorithms and ran them with the same maximum elapsed CPU time limit: Tmax = ρn m2 milliseconds, where ρ = 20, 40 and 60. The average RPD values of each algorithm over five runs are reported in Tables 911 and the best solutions are given in Table 14. Here, T T B is the minimum total tardiness known for each instance reported in Table 14. 26

From Table 9, we can see that our proposed EIGA provides the lowest average RPD with 0.24%, 0.21%, and 0.19% for τ = 1, τ = 2, and τ = 3, respectively, when ρ = 20. The proposed algorithm also produces the best results for all problem instances. In all three due date scenarios, IGAP R remains the closest competitor. It is interesting to note that the EIGA has larger gaps in results compared to other algorithms when the due dates are smaller. From Tables 10-11, when ρ = 40 and ρ = 60, the same observation can be made, i.e., the proposed algorithm is superior compared to other algorithms. Again, IGAP R is the closest competitor. Interestingly, as the allowed CPU time increases, improvement in the results produced by the proposed algorithm is relatively small. This indicates that the proposed EIGA converges rapidly and additional CPU time does not improve the current solution much. Table 9: Computational results produced by algorithms in comparison with ρ = 20. Best values are in bold. Instance n×m 20*5 20*10 20*20 50*5 50*10 50*20 100*5 100*10 100*20 200*10 200*20 500*20 Average

GA 0.04 0.00 0.05 1.07 2.35 1.58 0.96 1.34 1.51 1.47 1.12 1.48 1.08

DABC 0.01 0.00 0.02 0.72 1.55 1.35 0.99 0.85 1.38 1.35 1.17 1.28 0.89

τ =1 BEDA 0.02 0.00 0.03 0.54 1.70 2.17 1.03 0.90 1.37 1.18 0.75 1.02 0.89

IGA PR 0.00 0.00 0.00 0.39 0.93 0.66 0.67 1.05 0.88 0.90 0.96 1.24 0.64

EIGA 0.00 0.00 0.00 0.15 0.23 0.25 0.29 0.36 0.44 0.34 0.37 0.45 0.24

GA 0.04 0.02 0.08 1.04 1.72 1.77 1.15 1.64 1.23 1.19 1.00 1.16 1.00

DABC 0.03 0.00 0.04 0.87 1.80 0.89 0.99 0.65 0.84 1.23 1.10 0.73 0.76

τ =2 BEDA IGA PR 0.02 0.01 0.00 0.00 0.00 0.00 0.77 0.62 1.25 0.98 1.22 1.02 1.02 0.76 1.07 0.76 1.46 1.12 1.37 0.86 1.00 0.65 0.83 0.82 0.83 0.63

EIGA 0.00 0.00 0.00 0.13 0.26 0.27 0.32 0.29 0.41 0.36 0.26 0.22 0.21

GA 0.04 0.03 0.04 0.91 0.73 1.62 1.06 0.83 1.35 1.36 1.04 1.89 0.91

DABC 0.04 0.00 0.02 0.76 0.64 1.35 0.65 0.67 0.95 1.34 0.73 1.50 0.72

To ascertain the statistical significance of the results, we performed statistical tests to check whether the differences observed in Tables 9-11 are indeed statistically significant. First, an ANOVA test with a confidence level of α = 0.05 was performed, with each algorithm as a factor. The test results are shown in Table 12. The p-value of < 0.05 shows that there is a significant difference between the algorithms compared. Then, we show the 95% confidence interval plot in Figure 6, to compare the behavior of the algorithms. Again, the proposed algorithm significantly outperforms the other algorithms. Subsequently, we would like to know more about the performance of the 27

τ =3 BEDA IGA PR 0.03 0.02 0.00 0.00 0.01 0.00 0.71 0.64 0.64 0.63 1.01 0.59 0.77 0.60 0.69 0.68 0.79 0.74 0.89 0.89 0.87 0.57 1.18 0.79 0.63 0.51

EIGA 0.00 0.00 0.00 0.18 0.21 0.32 0.23 0.20 0.42 0.34 0.19 0.25 0.19

Table 10: Computational results produced by algorithms in comparison with ρ = 40. Best values are in bold. Instance n×m 20*5 20*10 20*20 50*5 50*10 50*20 100*5 100*10 100*20 200*10 200*20 500*20 Average

GA 0.00 0.00 0.00 0.83 1.81 1.26 0.86 1.07 1.27 1.30 1.03 1.40 0.90

DABC 0.00 0.00 0.00 0.56 1.22 1.09 0.83 0.68 1.08 1.13 1.10 1.15 0.74

τ =1 BEDA 0.00 0.00 0.00 0.42 1.36 1.67 0.84 0.69 1.16 1.05 0.70 0.85 0.73

IGA PR 0.00 0.00 0.00 0.32 0.74 0.59 0.54 0.82 0.80 0.75 0.87 1.12 0.55

EIGA 0.00 0.00 0.00 0.12 0.20 0.21 0.26 0.31 0.34 0.29 0.32 0.41 0.20

GA 0.00 0.00 0.00 0.81 1.46 1.50 0.99 1.27 1.07 1.04 0.92 1.02 0.74

DABC 0.00 0.00 0.00 0.68 1.40 0.73 0.90 0.52 0.76 1.06 0.93 0.69 0.56

τ =2 BEDA 0.00 0.00 0.00 0.61 0.97 0.99 0.87 0.82 1.14 1.22 0.94 0.72 0.60

IGA PR 0.00 0.00 0.00 0.48 0.80 0.93 0.64 0.61 0.91 0.79 0.60 0.74 0.48

EIGA 0.00 0.00 0.00 0.10 0.24 0.21 0.29 0.22 0.36 0.33 0.23 0.21 0.16

GA 0.00 0.00 0.00 0.70 0.59 1.42 0.85 0.74 1.06 1.26 0.88 1.67 0.76

DABC 0.00 0.00 0.00 0.59 0.52 1.09 0.54 0.57 0.74 1.13 0.66 1.26 0.59

τ =3 BEDA 0.00 0.00 0.00 0.57 0.50 0.86 0.61 0.60 0.72 0.82 0.82 1.09 0.55

IGA PR 0.00 0.00 0.00 0.51 0.56 0.50 0.46 0.62 0.64 0.76 0.54 0.74 0.45

EIGA 0.00 0.00 0.00 0.15 0.16 0.27 0.19 0.18 0.36 0.29 0.17 0.24 0.17

τ =3 BEDA 0.00 0.00 0.00 0.48 0.45 0.77 0.54 0.51 0.59 0.81 0.73 1.07 0.50

IGA PR 0.00 0.00 0.00 0.43 0.48 0.44 0.41 0.50 0.52 0.73 0.52 0.66 0.39

EIGA 0.00 0.00 0.00 0.12 0.14 0.24 0.17 0.15 0.31 0.26 0.15 0.23 0.15

Table 11: Computational results produced by algorithms in comparison with ρ = 60. Best values are in bold. Instance n×m 20*5 20*10 20*20 50*5 50*10 50*20 100*5 100*10 100*20 200*10 200*20 500*20 Average

GA 0.00 0.00 0.00 0.69 1.55 1.05 0.76 0.91 1.08 1.23 0.93 1.24 0.79

DABC 0.00 0.00 0.00 0.49 1.06 0.94 0.74 0.57 0.97 1.04 1.03 1.12 0.66

τ =1 BEDA 0.00 0.00 0.00 0.37 1.17 1.36 0.73 0.60 0.94 0.93 0.71 0.78 0.63

IGA PR 0.00 0.00 0.00 0.26 0.65 0.52 0.45 0.66 0.67 0.65 0.76 0.98 0.47

EIGA 0.00 0.00 0.00 0.10 0.16 0.18 0.23 0.26 0.29 0.26 0.32 0.40 0.18

GA 0.00 0.00 0.00 0.71 1.20 1.32 0.87 1.09 0.96 0.93 0.88 0.94 0.74

DABC 0.00 0.00 0.00 0.62 1.20 0.66 0.74 0.43 0.63 1.01 0.81 0.65 0.56

τ =2 BEDA 0.00 0.00 0.00 0.52 0.79 0.89 0.75 0.75 0.92 1.12 0.87 0.64 0.60

IGA PR 0.00 0.00 0.00 0.43 0.71 0.76 0.55 0.49 0.72 0.75 0.55 0.74 0.48

EIGA 0.00 0.00 0.00 0.11 0.19 0.18 0.24 0.19 0.31 0.30 0.21 0.20 0.16

GA 0.00 0.00 0.00 0.58 0.53 1.17 0.72 0.66 0.92 1.09 0.84 1.56 0.67

DABC 0.00 0.00 0.00 0.51 0.47 0.91 0.48 0.48 0.67 0.99 0.65 1.12 0.52

Table 12: ANOVA test results for the algorithms in comparison. Source

Degree of freedom

Sum of square

Mean square

F-statistic

p-value

Algorithm Error T otal

4 26995 26999

1359.01 4821 6180

339.755 0.179

1902.32

0.0000

algorithms for each timeout and due date scenario. Therefore, for each of the three timeout and due date scenarios (9 in total), we separately performed an ANOVA test. The p-value of 0.000 for all the 9 ANOVA tests confirmed the significant differences between the factors (algorithms compared). In 28

Figure 6: An interval plot of the algorithms in comparison

Table 13, we provide the Tukey and Fisher’s post-hoc test results. In each comparison, different letters mean that there is a significant difference between the algorithms. For example, for scenarios τ = 1 and ρ = 1, there are significant differences between the algorithms except the DABC and BEDA that share the same letter C i.e., these two algorithms perform equivalently. We can see that the proposed EIGA is statistically superior compared to other algorithms in all the 9 scenarios. Table 13: Grouping information using the Tukey and Fisher’s Least Significant Difference (LSD) methods. Algorithms ρ=1 T* F** GA D D DABC C C BEDA C C IGA PR B B EIGA A A T = Tukey Method F = Fisher LSD Method

τ =1 ρ=2 T F D E C D C C B B A A

ρ=3 T F E E D D C C B B A A

ρ=1 T F D D C C C C B B A A

τ =2 ρ=2 T F D E C D C C B B A A

ρ=3 T F D E C D C C B B A A

ρ=1 T F D D C C C C B B A A

τ =3 ρ=2 T F D E C D C C B B A A

ρ=3 T F D D C C C C B B A A

In Table 14, the best known solutions for all instances are summarized. Here, sources 1, 2, 3, 4 and 5 represent the GA, DABC, BEDA, IGAP R , and EIGA, respectively. As we can see, best solutions for small instances can 29

be obtained by almost all the five algorithms. For medium and especially large-scale instances, however, our proposed EIGA acquires the majority of the best solutions. In fact, all the best solutions for problem instances with 500 jobs and 20 machines are obtained by the EIGA. 5. Conclusion In this paper, we proposed a new, improved IGA for no-idle permutation flowshop scheduling with the total tardiness criterion. No-idle flowshop scheduling is widely seen in practical applications related to processing technology, e.g., fiber glass processing where glass batches are reduced to molten glass (Kalczynski and Kamburowski, 2005), or casting/assembly lines in which maximizing throughputs is of utmost importance, e.g., foundry production (Saadani et al., 2003). Core production machines in these environments must work with no idle time for technical and economic reasons. While most of the previous studies related to no-idle flowshop scheduling considered the makespan criterion, this work is among the few that have examined the total tardiness criterion, which has great managerial implications in manufacturing systems as managers become concerned with costs incurred when incomplete jobs approach their predetermined due dates. Situations like this can be perfidious and would lead to activation of penalty clauses in contracts or losing the entente with customers (Sen and Gupta, 1984). Improved aspects of the proposed algorithm include using a new variant of NEH for generating initial solutions, tardiness-guided destructionconstruction with a variable d value as well as a swap move alongside the typical insert move, and a time-dependent probability for dealing with nonimproving solutions. Computational evaluation and statistical analysis showed that each of the new elements has a significant impact on the algorithm’s performance. In addition, direct comparisons with state-of-the-art algorithms from the literature confirmed that the proposed EIGA is superior, with 189 new upper bounds obtained out of the 360 problem instances tested. Future work will investigate the possibility of applying some of the ideas used in this work to solve other real-world combinatorial problems, or considering other objectives such the total flow time or number of tardy jobs for no-idle permutation flowshop scheduling. There is also room to improve the speed-up method.

30

Table 14: Best known solutions obtained by the algorithms in comparison. Instance tai1 tai2 tai3 tai4 tai5 tai6 tai7 tai8 tai9 tai10 tai11 tai12 tai13 tai14 tai15 tai6 tai17 tai18 tai19 tai20 tai21 tai22 tai23 tai24 tai25 tai26 tai27 tai28 tai29 tai30 tai31 tai32 tai33 tai34 tai35 tai36 tai37 tai38 tai39 tai40 tai41 tai42 tai43 tai44 tai45 tai46 tai47 tai48 tai49 tai50 tai51 tai52 tai53 tai54 tai55 tai56 tai57 tai58 tai59 tai60

τ BKS 11967 10752 12041 11163 12812 13264 9017 10564 11887 10057 22317 17657 18568 19259 14414 19029 16089 18196 17898 16051 34821 32635 34341 31682 35635 37514 33469 36358 34532 1 42470 74628 63898 60972 64026 75718 63172 68517 65011 62556 61244 81900 79129 79990 86283 72261 101486 87530 100118 87320 87030 148793 176527 131053 162768 132781 172487 165312 144215 145502 159778

=1 Source 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 5 4,5 5 3,4,5 4,5 5 3,4,5 5 3,4,5 5 5 4,5 5 5 5 4,5 3,4,5 4,5 2 5 5 4,5 5 5 4,5 5 5 5 2 5

τ BKS 6941 6042 7463 5784 7844 8209 4395 5542 6645 5341 11988 7321 8863 10329 5289 9889 6879 8294 8053 5881 14548 13863 14114 12042 15168 17813 13378 16601 14331 23132 62494 51117 48915 51193 62841 50329 56355 52776 50949 48438 57245 55386 56250 60904 47768 76280 61896 75706 63424 61396 97060 127285 82986 113747 85087 123848 115714 95810 94775 109373

=2 Source 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 2,4,5 2,3,4,5 5 5 5 5 5 2,4,5 5 5 5 3 5 5 4,5 4 3 5 5 5 5 5 1,4,5 5 4 5 5 2 5 5

τ BKS 2543 2423 3637 1901 3174 3592 1036 1657 2618 1670 2993 1026 1430 2754 44 2305 825 993 800 522 989 630 336 233 888 1978 119 920 642 4659 50426 37421 37315 38611 50306 37577 44099 40585 39429 36137 33054 30994 32680 36270 25118 51227 37392 51309 39683 35829 46188 77923 35849 65049 38528 74931 66617 47422 45378 58575

=3 Source 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 1,2,3,4,5 2,5 3,4,5 4,5 3,5 5 2,3 4,5 3,4 5 2 5 4 1,3 5 5 5 2 3,5 4 2,4 5 5 4 5 2 5 4 5 3,4 5

Instance tai61 tai62 tai63 tai64 tai65 tai66 tai67 tai68 tai69 tai70 tai71 tai72 tai73 tai74 tai75 tai76 tai77 tai78 tai79 tai80 tai81 tai82 tai83 tai84 tai85 tai86 tai87 tai88 tai89 tai90 tai91 tai92 tai93 tai94 tai95 tai96 tai97 tai98 tai99 tai100 tai101 tai102 tai103 tai104 tai105 tai106 tai107 tai108 tai109 tai110 tai111 tai112 tai113 tai114 tai115 tai116 tai117 tai118 tai119 tai120

31

τ =1 BKS Source 271743 4,5 229957 5 227661 5 233252 5 254522 1 228838 5 277095 2 219424 5 237540 5 237485 5 333920 5 290198 5 356305 4,5 310758 5 314698 5 326151 3,4,5 369802 5 359609 5 367992 5 353311 5 533716 5 465508 5 499366 5 567817 5 580377 5 454265 5 487221 5 458026 5 518107 5 469132 5 1169937 4 1336088 4 1296716 5 1062142 5 1180247 5 1243417 3 1052148 5 1266999 5 1048700 5 1242593 5 1737592 5 1661558 5 1615834 5 1524487 5 1738186 5 1696078 5 1758017 5 1717058 5 1584909 5 1585713 5 8070550 5 8268024 5 9084684 5 7919050 5 8909261 5 8230729 5 8862049 5 9181333 5 7795109 5 8938644 5

τ =2 BKS Source 246030 5 205654 3 203594 3 208630 5 229770 5 203719 5 251984 5 195123 5 211016 4 212186 5 281557 5 240787 3,5 306024 4 259897 5 269233 5 278123 2 320583 5 309354 2,4 317298 5 301495 5 436495 4 365732 5 400682 2 467800 5 480954 5 356945 5 387776 5 361439 2 418068 5 367354 5 1063466 4 1233766 5 1198548 5 959858 3 1079056 5 1143142 5 952286 5 1164421 5 947510 1 1147959 5 1536847 5 1456097 4 1406059 5 1324866 3 1541546 5 1500083 5 1555221 5 1515864 5 1388513 5 1379115 5 7610936 5 7694540 5 8572106 5 7399708 5 8465129 5 7755439 5 8344824 5 8634395 5 7403881 5 8449019 5

τ =3 BKS Source 219781 4 180093 5 178541 3 184419 5 204638 4 179705 5 227584 5 171466 2 186197 5 187206 5 230105 5 192313 5 255711 4 206724 4 218915 3 230864 5 271789 5 260076 5 260076 5 250865 5 338076 5 265785 5 301721 4 369704 2 381685 5 256551 5 287419 3 255342 3 318146 5 269867 4 967522 5 1135778 5 1092913 5 857786 5 975410 5 1046325 3 848157 5 1060855 4 850362 5 1046142 4 1332345 5 1252401 4 1212453 5 1126264 4 1327727 5 1299676 5 1360587 5 1327552 5 1193121 5 1176529 5 7119751 5 7247543 5 8018526 5 6921652 5 7895349 5 7192395 5 7869830 5 8193202 5 6816062 5 7955856 5

References References Adiri, I. and Pohoryles, D. (1982). Flowshop/no-idle or no-wait scheduling to minimize the sum of completion times. Naval Research Logistics Quarterly, 29(3):495–504. Baptiste, P. and Hguny, L. (1997). A branch and bound algorithm for the f/no- idle/cmax. In Proceedings of the international conference on industrial engineering and production management, IEPM, volume 97, pages 429–438. Baraz, D. and Mosheiov, G. (2008). A note on a greedy heuristic for flowshop makespan minimization with no machine idle-time. European Journal of Operational Research, 184(2):810–813. Brah, S. (1996). A comparative analysis of due date based job sequencing rules in a flow shop with multiple processors. Production Planning & Control, 7(4):362–373. Cheng, C.-Y., Ying, K.-C., Chen, H.-H., and Lu, H.-S. (2019). Minimising makespan in distributed mixed no-idle flowshops. International Journal of Production Research, 57(1):48–60. Deng, G. and Gu, X. (2012). A hybrid discrete differential evolution algorithm for the no-idle permutation flow shop scheduling problem with makespan criterion. Computers & Operations Research, 39(9):2152–2160. Ding, J.-Y., Song, S., Gupta, J. N., Wang, C., Zhang, R., and Wu, C. (2016). New block properties for flowshop scheduling with blocking and their application in an iterated greedy algorithm. International Journal of Production Research, 54(16):4759–4772. Ding, J.-Y., Song, S., Gupta, J. N., Zhang, R., Chiong, R., and Wu, C. (2015). An improved iterated greedy algorithm with a tabu-based reconstruction strategy for the no-wait flowshop scheduling problem. Applied Soft Computing, 30:604–613. Dubois-Lacoste, J., Pagnozzi, F., and St¨ utzle, T. (2017). An iterated greedy algorithm with optimization of partial solutions for the makespan permutation flowshop problem. Computers & Operations Research, 81:160–166. 32

Fernandez-Viagas, V. and Framinan, J. M. (2015). A bounded-search iterated greedy algorithm for the distributed permutation flowshop scheduling problem. International Journal of Production Research, 53(4):1111–1123. Goncharov, Y. and Sevastyanov, S. (2009). The flow shop problem with no-idle constraints: A review and approximation. European Journal of Operational Research, 196(2):450–456. Kalczynski, P. J. and Kamburowski, J. (2005). A heuristic for minimizing the makespan in no-idle permutation flow shops. Computers & Industrial Engineering, 49(1):146–154. Karabulut, K. (2016). A hybrid iterated greedy algorithm for total tardiness minimization in permutation flowshops. Computers & Industrial Engineering, 98:300–307. Kim, Y.-D. (1993). Heuristics for flowshop scheduling problems minimizing mean tardiness. Journal of the Operational Research Society, 44(1):19–28. Liu, J. and Reeves, C. R. (2001). Constructive and composite heuristic solutions to the p// ci scheduling problem. European Journal of Operational Research, 132(2):439–452. Nagano, M. S., Rossi, F. L., and Martarelli, N. J. (2019). High-performing heuristics to minimize flowtime in no-idle permutation flowshop. Engineering Optimization, 51(2):185–198. Nawaz, M., Enscore, E. E., and Ham, I. (1983). A heuristic algorithm for the m-machine, n-job flow-shop sequencing problem. Omega, 11(1):91–95. Pan, Q.-K. and Ruiz, R. (2012). Local search methods for the flowshop scheduling problem with flowtime minimization. European Journal of Operational Research, 222(1):31–43. Pan, Q.-K. and Ruiz, R. (2014). An effective iterated greedy algorithm for the mixed no-idle permutation flowshop scheduling problem. Omega, 44:41–50. Pan, Q.-K., Tasgetiren, M. F., and Liang, Y.-C. (2008). A discrete differential evolution algorithm for the permutation flowshop scheduling problem. Computers & Industrial Engineering, 55(4):795–816.

33

Pan, Q.-K. and Wang, L. (2008a). No-idle permutation flow shop scheduling based on a hybrid discrete particle swarm optimization algorithm. The International Journal of Advanced Manufacturing Technology, 39(7-8):796– 807. Pan, Q.-K. and Wang, L. (2008b). A novel differential evolution algorithm for no-idle permutation flow-shop scheduling problems. European Journal of Industrial Engineering, 2(3):279–297. Rajendran, C. (1993). Heuristic algorithm for scheduling in a flowshop to minimize total flowtime. International Journal of Production Economics, 29(1):65–73. Riahi, V., Newton, M. H., Su, K., and Sattar, A. (2019). Constraint guided accelerated search for mixed blocking permutation flowshop scheduling. Computers & Operations Research, 102:102–120. Ribas, I., Companys, R., and Tort-Martorell, X. (2011). An iterated greedy algorithm for the flowshop scheduling problem with blocking. Omega, 39(3):293–301. Ruiz, R., Pan, Q.-K., and Naderi, B. (2019). Iterated greedy methods for the distributed permutation flowshop scheduling problem. Omega, 83:213–222. Ruiz, R. and St¨ utzle, T. (2007). A simple and effective iterated greedy algorithm for the permutation flowshop scheduling problem. European Journal of Operational Research, 177(3):2033–2049. Ruiz, R. and St¨ utzle, T. (2008). An iterated greedy heuristic for the sequence dependent setup times flowshop problem with makespan and weighted tardiness objectives. European Journal of Operational Research, 187(3):1143– 1159. Ruiz, R., Vallada, E., and Fern´andez-Mart´ınez, C. (2009). Scheduling in flowshops with no-idle machines. In Computational Intelligence in Flow Shop and Job Shop Scheduling, pages 21–51. Springer. Saadani, N. E. H., Guinet, A., and Moalla, M. (2003). Three stage no-idle flow-shops. Computers & industrial engineering, 44(3):425–434.

34

Saadani, N. E. H., Guinet, A., and Moalla, M. (2005). A travelling salesman approach to solve the f/no-idle/c max problem. European Journal of Operational Research, 161(1):11–20. Sen, T. and Gupta, S. K. (1984). A state-of-art survey of static scheduling research involving due dates. Omega, 12(1):63–76. Shao, W., Pi, D., and Shao, Z. (2017). Memetic algorithm with node and edge histogram for no-idle flow shop scheduling problem to minimize the makespan criterion. Applied Soft Computing, 54:164–182. Shao, W., Pi, D., and Shao, Z. (2018). A hybrid discrete teaching-learning based meta-heuristic for solving no-idle flow shop scheduling problem with total tardiness criterion. Computers & Operations Research, 94:89–105. Shen, J.-n., Wang, L., and Wang, S.-y. (2015). A bi-population eda for solving the no-idle permutation flow-shop scheduling problem with the total tardiness criterion. Knowledge-Based Systems, 74:167–175. Taillard, E. (1993). Benchmarks for basic scheduling problems. european journal of operational research, 64(2):278–285. Talbi, E.-G. (2009). Metaheuristics: from design to implementation, volume 74. John Wiley & Sons. Tasgetiren, M. F., Kizilay, D., Pan, Q.-K., and Suganthan, P. N. (2017). Iterated greedy algorithms for the blocking flowshop scheduling problem with makespan criterion. Computers & Operations Research, 77:111–126. Tasgetiren, M. F., Pan, Q.-K., Suganthan, P., and Jin Chua, T. (2011). A differential evolution algorithm for the no-idle flowshop scheduling problem with total tardiness criterion. International Journal of Production Research, 49(16):5033–5050. Tasgetiren, M. F., Pan, Q.-K., Suganthan, P., and Oner, A. (2013a). A discrete artificial bee colony algorithm for the no-idle permutation flowshop scheduling problem with the total tardiness criterion. Applied Mathematical Modelling, 37(10):6758–6779.

35

Tasgetiren, M. F., Pan, Q.-K., Suganthan, P. N., and Buyukdagli, O. (2013b). A variable iterated greedy algorithm with differential evolution for the noidle permutation flowshop scheduling problem. Computers & Operations Research, 40(7):1729–1743. Vachajitpan, P. (1982). Job sequencing with continuous machine operation. Computers & Industrial Engineering, 6(3):255–259. Woollam, C. (1986). Flowshop with no idle machine time allowed. Computers & industrial engineering, 10(1):69–76. Ying, K.-C., Lin, S.-W., Cheng, C.-Y., and He, C.-D. (2017). Iterated reference greedy algorithm for solving distributed no-idle permutation flowshop scheduling problems. Computers & Industrial Engineering. Zhou, Y., Chen, H., and Zhou, G. (2014). Invasive weed optimization algorithm for optimization no-idle flow shop scheduling problem. Neurocomputing, 137:285–292.

36