Efficient non-population-based algorithms for the permutation flowshop scheduling problem with makespan minimisation subject to a maximum tardiness

Ecient non-population-based algorithms for the permutation owshop scheduling problem with makespan minimisation subject to a maximum tardiness Victor Fernandez-Viagas 1∗ , Jose M. Framinan 1 † 1 Industrial Management, School of Engineering, University of Seville, Ave. Descubrimientos s/n, E41092 Seville, Spain, {vfernandezviagas,framinan}@us.es July 19, 2015 Abstract This paper focuses on the problem of scheduling jobs in a permutation owshop with the objective of makespan minimisation subject to a maximum allowed tardiness for the jobs, a problem that combines two desirable manufacturing objectives related to machine utilisation and to customer satisfaction. Although several approximate algorithms have been proposed for this NP-hard problem, none of them can use the excellent speed-up method by Taillard (1990) for makespan minimisation due to the special structure of the problem under consideration. In this paper, several properties of the problem are dened in order to be able to partly apply Taillard's acceleration. This mechanism, together with a novel feasible tabu local search method, allows us to further exploit the structure of solutions of the problem, and are incorporated in two proposed algorithms: a bounded-insertion-based constructive heuristic and an advanced non-population-based algorithm. These algorithms are compared with state-of-the-art algorithms under the same computer conditions. The results show that both algorithms improve existing ones and therefore, constitute the new state- of-art approximate solution procedures for the problem. Keywords: Scheduling, Flowshop, Heuristics, NEH, PFSP, maximum tardiness, makespan, bounded insertion, non-population algorithm ∗ Corresponding author. Tel.: +34-954487220. † This is an Accepted Manuscript of an article published by Elsevier in Computers & Operations Research, Volume 64, Dec 2015, Pages 86-96, available online: http://dx.doi.org/10.1016/j.cor.2015.05.006 1 1 Introduction A permutation ow shop is a manufacturing layout in which a set of machines are arranged to be visited by a number of jobs in the same order, assuming that the sequence of the jobs remains the same for all machines. Usual additional hypotheses include the simultaneous availability of all jobs and all machines and deterministic processing times, among others (see e.g. Framinan et al., 2014 for a complete list of assumptions). Several criteria can be established to measure the performance of the dierent schedules (see e.g. Sun et al., 2011). Among them, the maximum completion time of a sequence or makespan is related to resource usage (see e.g. Ruiz and Maroto, 2005 and Fernandez-Viagas and Framinan, 2014b), while tardiness refers to the delay of the completion time of a job with respect to its committed due date (see e.g. Fernandez-Viagas and Framinan, 2015a and Framinan and Leisten, 2008). Since these are key aspects in manufacturing companies' competitiveness, it seems appropriate to consider both objectives together. Regarding tardiness minimisation, customer due dates may be regarded as `hard' constraints (i.e. deadlines) in some manufacturing scenarios, while in others some exibility is allowed by the customer as long as the deviation from the completion times of the jobs is limited. In contrast, makespan is an intracompany criteria that is related to maximising machine utilisation, which in turns minimises xed unit costs. Therefore, one option to balance both objectives is to seek the minimisation of the makespan while allowing only a given deviation from the committed due dates, expressed as the maximum tardiness allowed. Note that this problem includes the special case where no deviation from the jobs' due dates is allowed, thus forcing the fullment of the committed due dates. According to the notation by T'Kindt and Billaut (2006), the problem described in the previous paragraph can be denoted as Fm|prmu|ϵ(Cmax/Tmax) . This problem belongs to the class of ϵ -constrained multi-criteria scheduling problems, and it has been the subject of several research contributions in the last decades. Since the minimisation of any of the individual criteria (either makespan or maximum tardiness) in a ow shop is NP-hard, the research eort has focused on approximate procedures providing good but not necessarily optimal solutions in a relative short period of time. In this regard, the works by Daniels and Chambers (1990), Chakravarthy and Rajendran (1999), Framinan and Leisten (2006), and Ruiz and Allahverdi (2009) 2 develop dierent heuristics either for the problem, or for general cases of the problem. In this paper, we propose a constructive heuristic and a metaheuristic that exploits the specic structure of solutions of the problem to reduce the search space and to accelerate the evaluation of solutions. Both algorithms improve existing ones by a larger degree and constitute therefore the new state-of-art approximate solution procedures for the problem. The remainder of the paper is structured as follows: Section 2 describes the problem under consideration and its state-of-the-art. In Section 3, some denitions and properties of the problem are dened. Section 4 is devoted to propose two algorithms (a constructive heuristic and a metaheuristic) which use the properties discussed previously. The algorithms are compared with the (up to now) state-of-the-art algorithms in Section 4 and, nally, conclusions are discussed in Section 5. 2 Problem Statement and State of the Art In the Fm|prmu|ϵ(Cmax/Tmax) problem under study, n jobs have to be scheduled in a owshop composed of m one-machine stage. The processing time of job l on machine i is dened as pil . Following this notation, given a sequence of jobs Π := (π1, . . . , πn) , the processing time of job in position j , πj , is denoted as piπj . Analogously, Ciπj(Π) denotes the completion time of job πj on machine i according to the schedule given by Π . The makespan of the sequence Π is given by the completion time of last job in last machine, Cm,πn(Π) , which is denoted as Cmax(Π) . Whenever it does not lead to confusion, the sequence Π is omitted in the notation of the completion times and makespan. In a similar manner, if dπj is the due date of job πj , its tardiness is dened as Tπj(Π) = max{Cmπj(Π)−dπj,0} and the maximum tardiness of the sequence Π as Tmax(Π) = maxj=1,...,n{Tπj(Π)} . As with the makespan, Π is omitted when it is clear from the context. The goal of the problem is to nd a schedule Π for which the makespan is minimum subject to Tmax(Π) ≤ϵ . As mentioned in Section 1, the problem is NP-hard since the minimisation of each individual criterion is already an NP-hard problem for the permutation owshop (see e.g. T'Kindt and Billaut, 2006 for a detailed proof). Consequently, the interest lies in nding ecient approximate methods 3 or heuristics . Given the clear connection between our problem and that of makespan minimisation, most of the algorithms to solve the problem are based on the best heuristic for makespan minimisation: i.e. the NEH heuristic by Nawaz et al. (1983). It is then useful to recall the main steps in the NEH heuristic, which can be described as follows: 1. Jobs are ordered according to non-increasing sum of processing times. 2. A partial sequence is constructed only with rst job of the previous phase. 3. Each remaining job of initial phase is iteratively inserted in all positions of the partial sequence. The makespan of all these sequences is evaluated, and the partial schedule for which the lowest makespan is reached is selected for the next iteration. 4. The procedure is repeated until no more jobs are available. The above steps make clear that the computational burden of the NEH lies on the evaluation of all possible insertions in Step 3. In Taillard (1990), a mechanism named in the following Taillard's acceleration  is proposed so the computational complexity of evaluating all insertions is equivalent to that of evaluating one sequence. In order to explain Taillard's acceleration, let us rst dene three variables (for a more detailed description of the variables, see Taillard, 1990): •ei,πj : Earliest completion time of job in position j in machine i . •qi,πj . Once a sequence of jobs has been dened (and therefore the makespan of this sequence is obtained), qi,πj is the dierence between the makespan and the latest starting time of job in position j in machine i . •fi,πj . Earliest completion time of the new job σ when it is inserted before job in position j in machine i . These are computed using ei,πj and the processing times of σ . By means of these variables, the partial makespan Cj max when introducing job σ before job in position j can be determined using the following expression: Cj max = max i(fi,πj+qi,πj) (1) 4 As the job σ is inserted in the position with minimum makespan, the makespan of the sequence is dened by: Cmax = min j(Cj max) (2) As it can be seen, although the cost of evaluating the insertion slot with lowest makespan is greatly reduced by Taillard's acceleration, the completion time of each job cannot be obtained using this mechanism, and therefore its tardiness cannot be computed. As a consequence, none of the heuristics proposed up-to-now in the literature for the problem under consideration use this mechanism. Among the contributions on the Fm|prmu|ϵ(Cmax/Tmax) problem, Daniels and Chambers (1990) were the rst in proposing a constructive heuristic. In their heuristic, assuming a partial sequence Π formed by already scheduled jobs, a (partial) sequence is constructed for each nonscheduled job uk by placing it as the rst job, and then scheduling the jobs in Π after uk according to the NEH algorithm. Out of these so-obtained sequences, the one with the lowest makespan is chosen for the next iterations (consequently, uk is removed from the non-scheduled jobs set for the next iteration). Chakravarthy and Rajendran (1999) propose a simulated annealing algorithm to solve the Fm|prmu|ϵ(Z/Tmax) where Z=λ·Cmax + (1 −λ)·Tmax , λ∈[0,1] . Clearly, our problem is a special case of their problem when λ= 1 . Their algorithm begins with the best sequence among the solutions found by the NEH heuristic, the earliest due date rule and the least slack rule (jobs ordered according to ascending order of dj−∑m i=1 pij ). The procedure iteratively samples neighbour solutions (using an adjacent pairwise interchange neighbourhood) until the stopping criterion is fullled. Framinan and Leisten (2006) propose a constructive heuristic, denoted in the following as FL, based on the NEH algorithm to solve the Fm|prmu|ϵ(Cmax/Tmax) problem. The heuristic tries to improve the makespan without worsening the tardiness by using a property of the problem. The heuristic is compared with those of Daniels and Chambers (1990) and Chakravarthy and Rajendran (1999) for small and big instances. The results show that the FL outperforms the other ones in 5 terms of both the quality of the solutions and the number of the feasible solutions obtained. Finally, Ruiz and Allahverdi (2009) propose an iterated optimization algorithm to solve the Fm|prmu|ϵ(Z/Tmax) problem. More specically, they proposed a high-performance Genetic Algorithm (GA in the following) where the selection procedure is based on n -tournament (see Ruiz and Allahverdi, 2007). The tness values of the individuals are calculated depending on whether all individuals are feasible; feasible and infeasible; or only infeasible. The algorithm outperforms the FL for the Fm|prmu|ϵ(Z/Tmax) problem in an extended benchmark. Nevertheless, GA and FL were not compared for the specic Fm|prmu|ϵ(Cmax/Tmax) problem. To summarise the state of the art regarding the problem under consideration, there are some ecient heuristics for the problem, but their performance is not completely clear, as the comparison between the most ecient contributions (i.e. GA and FL) has been only partially conducted. In addition, both mechanisms made extensive use of insertion neighbourhoods, so it could be extremely interesting to devise a mechanism similar to that by Taillard to reduce the computational burden. Finally, it is also to note that all existing procedures make little use (or no use at all) of the knowledge on the problem domain. 3 Problem Properties As mentioned in Section 2, Taillard's acceleration does not compute the completion times of each job and therefore, it cannot be used to compute the maximum tardiness of the sequence. Indeed, our problem is complicated by the fact that, when inserting a new job σ in position r of an existing partial sequence, an infeasible solution can be obtained due to either the increase in the completion times of the jobs after σ , or due to the completion time of job σ itself. In order to further classify these two possibilities, let us introduce the following denitions: Denition 3.1 (First Feasible Position) . Given a feasible (partial) sequence Π := (π1, . . . , πk) , and a non scheduled job σ , let Π′ r:= (π1, . . . , πr−1, σ, πr, . . . , πk) be the (partial) sequence obtained by the insertion of σ in position r of Π . Then, the First Feasible Position ( FFP ) is dened as follows: FFP(Π, σ) := arg min 1≤r≤k+1 {Tπj(Π′ r)≤ϵ∀j=r, . . . , k} 6 As it can be seen from the denition, FFP is the lowest position where a new job can be inserted in an existing sequence without causing infeasible due dates in any of the jobs resulting in positions later than the insertion point. It is clear that, in a given instance of the problem and a partial sequence Π , it is not possible to obtain feasible schedules by inserting a non scheduled job σ into a position j < FFP(Π, σ) . Note also that obtaining FFP for a tuple Π and σ does not guarantee that Π′ j is feasible for j≥FFP , since the computation of FFP does not take into account the potential infeasibility caused by job σ . Denition 3.2 (Last Feasible Position) . Given a feasible (partial) sequence Π := (π1, . . . , πk) , and a non scheduled job σ , let Π′ r:= (π1, . . . , πr−1, σ, πr, . . . , πk) be the (partial) sequence obtained by the insertion of σ in position r of Π . Then, the Last Feasible Position ( LFP ) is dened as follows: LFP(Π, σ) := arg max 1≤r≤k+1 {Tσ(Π′ r)≤ϵ} In this manner LFP is the highest position r where job πr can be inserted without making its completion time infeasible. Note that the feasibility of the jobs in positions r+ 1, r + 2, . . . is not considered when computing LFP . The calculation of both limits is of interest due to some straightforward observations which follow from both denitions: 1. If FFP(Π, σ)> LF P (Π, σ) for a given tuple Π and σ , then no feasible sequence can be obtained by inserting σ into Π . 2. If FFP(Π, σ)≤LFP (Π, σ) for a given tuple Π and σ , then the sequence obtained by inserting σ in position FFP of Π is feasible. 3. If FFP(Π, σ)≤LF P (Π, σ) for a given tuple Π and σ , then at least one feasible sequence can be obtained by inserting σ in positions between FFP and LFP , inclusive. 4. For a given tuple Π and σ with FFP(Π, σ)≤LF P (Π, σ) , the set of feasible sequences obtained by inserting σ in positions between FFP and LFP represent all feasible sequences 7 that can be obtained by inserting σ into Π . Once FFP and LFP are obtained, the sequence with the lowest makespan can be computed by using Taillard's acceleration between both bounds, i.e.: Cmax = min j(Cj max)j=FFP, . . . , LFP (3) where Cmax is obtained using Expressions (1) and (2). Note that the so-found sequence is not necessarily feasible. The advantage of this mechanism lies in speeding up the computations. In view of the above expressions, the challenge now is to compute both FFP and LFP in an ecient manner. To do so, we introduce the following properties: Property 3.1. Given an instance of the Fm|prmu|ϵ(Cmax/Tmax) problem, and given a tuple Π and σ of a partial sequence and a non-scheduled job respectively, then FFP(Π, σ) := 1 + arg max j {Cmπj+ min 1≤i≤mpiσ −dπj> ϵ} is a lower bound for FFP(Π, σ) . Furthermore, FFP can be computed in O(n·m) Proof. Recall that the completion times of the jobs in Π placed after the insertion of a job σ must increase at least mini(pi,σ), i ∈[1,· · · , m] (see also Fernandez-Viagas and Framinan, 2014a). Therefore, Cmπj+ min1≤i≤mpiσ is a lower bound of completion time of job in position j after the insertion of σ in position r < j . As a consequence, if Cmπj+ min1≤i≤mpiσ −dπj> ϵ , then the due date of job in position j is always infeasible, so FFP is a lower bound of FFP . FFP can be computed in two steps: First M= min1≤i≤mpiσ is computed in O(m) . Next, the expression Ej=Cm,πj+M−dπj is computed in O(n·m) for j= 1,2, . . . until it veries that Ej> ϵ . It is thus clear that the computation of FFP is O(n·m) . Property 3.2. Given an instance of the Fm|prmu|ϵ(Cmax/Tmax) problem, and given a tuple Π and σ of a partial sequence and a non-scheduled job respectively, then LFP can be computed in O(n·m) . Proof. First ei,πj the earliest completion times of the jobs before σ are computed in O(n·m) (see Section 2). Then, for a position k where σ can be inserted, the completion time of σ is computed 8 in O(m) by adding the processing times of σ to ei,πj , and the result is compared to ϵ . Since this comparison is performed for all positions prior to the candidate position where the new jobs is to be inserted, it is clear that LFP can be computed in O(n·m) . The detailed pseudo code is presented in Figure 2. Equipped with these problem properties, in Section 4 we propose ecient approximate procedures based on the insertion of jobs into existing partial schedules. More specically, in Section 4.1 we present a constructive heuristic for the problem whereas in Section 4.2 we present a nonpopulation based metaheuristic. 4 Proposed Algorithms In this section, we propose two new approximate algorithms to solve the Fm|prmu|ϵ(Cmax/Tmax) problem. The rst one is a constructive heuristic based on a so-called bounded insertion using some properties of the problem, while the second one is a generic non-population-based algorithm searching for feasible solutions in each iteration that uses the proposed constructive heuristic as starting solution. In this manner, we attempt to provide ecient approximate solutions for a wide range of decision time intervals. To speed up the computation times, both algorithms use the problem properties presented in Section 3. 4.1 Bounded-Insertion-Based Constructive Heuristic, BICH In this Section, we present a constructive heuristic based on a bounded insertion (BICH) that also repairs infeasibility by means of a tabu local search in each iteration (see pseudo code in Figure 1). More specically, the algorithm obtains a sequence Π := (π1, . . . , πn) in the following manner: Initially, jobs are sorted in non ascending order of the sum of their processing times, so a sorted sequence α:= (α1, . . . , αn) is obtained. The rst job in the sorted sequence is also the rst job in Π , i.e. π1=α1 . Then, the remaining jobs in α are inserted in Π one by one in the following manner: in iteration k ( k∈[2, n] ), job αk is removed from α and the following steps are carried out: 9 solution, Π2 , is worse than the local search optimum, Π , it is maintained only if: random ≤exp {−(Cmax(Π2)−Cmax(Π)) Temperature } where random is a random number between 0 and 1 and the Temperature is a function that depends on parameter T : Temperature =T·∑∀i∑∀jpi,j n·m·10 The temperature parameter has been generated following the suggestions by Osman and Potts (1989) (see e.g. Fernandez-Viagas and Framinan, 2014a, Pan et al., 2008 and Ruiz and Stützle, 2007 for similar approaches). 4.3 Experimental Parameter Tuning The proposed algorithms use three parameters: T , d , and x . Therefore, it is interesting to investigate the values of these parameters for which the algorithms reach the best performance. Parameter x is used in both BICH and ANP A , while the other two parameters are included only in ANP A . In order to simplify the experimentation, the three parameters have been tested only for ANP A , and the value obtained for parameter x was also chosen for the constructive heuristic BICH . The level of parameters tested are: •T∈[0.1,0.2,0.3,0.4] •d∈[4,5,6,7] •x∈[10,20,30] ANP A is tested following the same calibration test as in Vallada and Ruiz (2010), which in turn is based in Vallada et al. (2008). This calibration benchmark can be summarized as follows. The number of jobs and machines is set to n={50,150,250,350} , m={10,30,50} and due dates are generated according to the procedure employed by Gelders and Sambandam (1978) using an uniform distribution between P·(1−T−R/2) and P·(1−T+R/2) . Parameters T and P take the 16 following values in the test: T={0.2,0.4,0.6} and R={0.2,0.6,1.0} . Additionally, processing times are generated according to a uniform distribution between 1 and 99. For each combination of parameters n , m , T and R , two instances are generated summing 216 instances. The stopping criterion adopted is to halt the procedure when the CPU time in milliseconds reaches the value n·(m/2) ·20 . To establish statistically signicant dierences between parameters T, d and x , a non-parametric Kruskal-Wallis test is performed, since the normality and homoscedasticity assumptions required for an analysis of variance were not satised. As a result of the test, statistically signicant dierences between the levels of the parameters x and T were found, but not for d since the signicance values were 0.011, 0.000 and 0.870 respectively. The best combination of parameters was found for d= 5 , T= 0.4 and x= 30 , so these were used in the computational experience carried out in the next section. 5 Computational Results In this Section, the performance of the proposed algorithms BICH and ANP A is compared with the best algorithms so far for the problem, i.e. the GA by Ruiz and Allahverdi (2009) and the constructive heuristic FL by Framinan and Leisten (2006). Additionally, two ecient heuristics for makespan minimisation (i.e. the NEH heuristic, and the iterated greedy algorithm, IGRS _ LS ) are included in the comparison as they are two of the most ecient constructive heuristic and iterative improvement algorithm, respectively. The adaptations of these heuristic to our problem are denoted A _ NEH and A _ IGA , respectively. When adapting both methods to the proposed problem, the following assumptions are adopted: • The objective function remains the original of these algorithms, i.e. the minimisation of makespan. • In the insertion phases of the algorithms, only feasible sequences are considered, i.e. the jobs to be inserted are placed in the feasible position with the lowest makespan. • Taillard's acceleration is removed from both algorithms since the calculation of the tardiness 17 for each job does not allow its use. To be able to determine the most ecient algorithms for the problem, all methods have been coded under the same programming language, C#, and under the same computer (an Intel Core i7-3770 with 3.4 GHz and 16GB RAM). The algorithms are tested using the set of instances of the benchmark of Vallada et al. (2008), which includes n={50,150,250,350} jobs and m={10,30,50} machines. Processing times are generated following an uniform distribution [1,99] as well as the due dates are generated according to Gelders and Sambandam (1978) as explained in Section 4.3 with T={0.2,0.4,0.6} and R={0.2,0.6,1.0} . 5 instances are obtained for each combination of n , m , T and R forming a total of 540 instances. Furthermore, in order to increase the accuracy of the iterated improvement algorithms, ve runs have been performed per instance and the average values are recorded for the makespan and for the CPU times. The same stopping criteria as in Ruiz and Allahverdi (2009) are applied for the iterative improvement algorithms. These stopping criteria depend on the size of the instance (i.e. the number of jobs and machines) following the expression t·n·m·/2 milliseconds where the values 2, 5, 20 and 60 are tested for the parameter t . The FL and BICH constructive heuristics stop naturally when their nal sequences are constructed. Due to the fact that the problem under consideration is subject to maximum tardiness, the evaluation of the quality of both constructive and iterative algorithms is not trivial. Usually, the decision maker would rst look for the feasibility of the solutions (i.e. tardiness of each job lower than the maximum tardiness) and, once it is achieved, he/she would look for a low value in the makespan. Finally, the quality of the sequences obtained by each algorithm has to be balanced against the time interval required to obtain the sequences, as the high CPU time requirements posed by some of the algorithms may not be acceptable for some scenarios. Therefore, there is a trade-o among these goals that increases the diculty of a direct comparison of the algorithms. To make an exhaustive analysis of all these aspects, three indicators have been chosen to determine the quality of the solutions obtained by the algorithms, as well as the CPU time required to obtain the solutions. The indicators are: • Number of feasible solutions obtained by each procedure. 18 • Makespan value of the solution (in terms of Average Relative Percentage Deviation) obtained by each procedure. • Number of instances with the best solution obtained by each procedure. These indicators are discussed in the next subsections. 5.1 Number of Feasible Solutions The average number of feasible solutions obtained by the algorithms are shown in Table 1. For 494 instances out of the 540 instances in the benchmark it was possible to nd a feasible solution by one/several algorithms. Among them, 492 feasible solutions were found by ANP A for the stopping criterion t= 60 , being the best algorithm in terms of the number of feasible solutions obtained. Next is the ANP A heuristic with 491 feasible for the stopping criteria t= 2 , t= 5 and t= 20 , 490 for the BICH heuristic. 457 feasible solutions were found by GA with t= 60 followed by the FL heuristic with 448 feasible solutions. The worst results were obtained for A _ NEH and A _ IGA algorithms. It is worth to note that both BICH and ANP A found more feasible solutions within lesser CPU time than the rest of the procedures, a remarkable result specially as BICH had very small CPU requirements. As it can be seen in the Table 1 and in Figure 7, the ecient algorithms following this criterion would be: BICH , ANP A(t= 2) and ANP A(t= 60) . 5.2 Average Relative Percentage Deviation The makespan of the solutions obtained by each algorithm can be evaluated by means of the Relative Percentage Deviation ( RPD ), which can be dened as follows: RPDi=Cmax,i −Best Best ·100 where Cmax,i is the makespan of algorithm i while Best is the best known value of the makespan (the lowest among the implemented heuristics). It has to be noted that RPD is 19 Figure 7: Number of feasible solutions vs Average CPU times for each algorithm. X-Axis is shown in logarithmic scale. 20 Table 1: Values of the three quality indicators and the average CPU time of each algorithm Algorithm #Feasible Solutions # Best Instances ARPD Average CPU Time (s.) A _ NEH 277 9 2.83 1.27 BICH 490 14 3.24 0.43 FL 448 1 6.96 15.22 ANPA(t= 2) 491 69 0.76 6 A _ IGA(t= 2) 303 26 1.32 6 GA(t= 2) 435 18 2.73 6 ANPA(t= 5) 491 80 0.52 15 A _ IGA(t= 5) 314 28 1.13 15 GA(t= 5) 439 26 2.23 15 ANPA(t= 20) 491 112 0.20 60 A _ IGA(t= 20) 333 43 0.80 60 GA(t= 20) 446 39 1.63 60 ANPA(t= 60) 492 491 0.00 180 A _ IGA(t= 60) 342 68 0.56 180 GA(t= 60) 457 47 1.37 180 computed only if a feasible solution is found by the algorithm, otherwise the results would be greatly biased. The average RPD (denoted as ARPD ) is shown in Table 1. Since certain algorithms do not nd feasible solutions for some instances for which others do, the ARPD is calculated with dierent sample sizes depending on the algorithm, e.g. 333 instances for A _ IGA(t= 20) and instances 491 by ANP A(t= 20) . This fact might cause that algorithms with lesser feasible solutions than other ones could have less ARPD , as it is the case with A _ IGA(t= 20) and A _ IGA(t= 60) as compared to GA(t= 60) . Nevertheless, we also include it in the analysis since it is the usual way in which this analysis is carried out (see e.g. Framinan and Leisten, 2006; Ruiz and Allahverdi, 2009). As it can be seen in Figure 8, the ecient heuristics with ARPD as indicator would be A _ NEH , BICH , ANP A(t= 2) , ANP A(t= 5) , ANP A(t= 20) and ANP A(t= 60) . In order to statistically justify this statement, we use Holm's procedure (Holm, 1979) with the following hypotheses: •H1 : ANP A(t= 2) = GA(t= 2) •H2 : ANP A(t= 2) = A _ IGA(t= 2) 21 Figure 8: ARP D vs Average CPU times for each algorithm. X-Axis is shown in logarithmic scale. 22 Table 2: Holm's procedure for multiple hypotheses. R indicate that hypothesis is reject by Mann-Whitney and/or Holm's procedure i Hip -value Mann-Whitney α/(k−i+ 1) Holm's Procedure 1ANPA(t= 2) = GA(t= 2) 0.000 R 0.0056 R 2ANPA(t= 2) = A _ IGA(t= 2) 0.000 R 0.0063 R 3ANPA(t= 5) = GA(t= 5) 0.000 R 0.0071 R 4ANPA(t= 5) = A _ IGA(t= 5) 0.000 R 0.0083 R 5ANPA(t= 5) = FL 0.000 R 0.0100 R 6ANPA(t= 20) = GA(t= 20) 0.000 R 0.0125 R 7ANPA(t= 20) = A _ IGA(t= 20) 0.000 R 0.0167 R 8ANPA(t= 60) = GA(t= 60) 0.000 R 0.0250 R 9ANPA(t= 60) = A _ IGA(t= 60) 0.000 R 0.0500 R •H3 : ANP A(t= 5) = GA(t= 5) •H4 : ANP A(t= 5) = A _ IGA(t= 5) •H5 : ANP A(t= 5) = FL •H6 : ANP A(t= 20) = GA(t= 20) •H7 : ANP A(t= 20) = A _ IGA(t= 20) •H8 : ANP A(t= 60) = GA(t= 60) •H9 : ANP A(t= 60) = A _ IGA(t= 60) The p -value of each hypothesis is calculated using a non-parametric Mann-Whitney test (see Pan et al., 2008). Then, Holm's procedure orders the hypotheses according to these p -values in non-decreasing order. The procedure rejects hypothesis i if its p -value is lower than α/(k−i+ 1) where k is the number of hypotheses. The results of this statistical analysis are shown in Table 2. As the p -values are always lower than α/(k−i+ 1) , each hypothesis is rejected justifying the statement regarding the ecient algorithms in terms of their ARPD . In fact, each p -value of the non-parametric Mann-Whitney analysis is 0.000. 23 5.3 Number of instances with the best makespan The third indicator used in this paper is the number of instances where each algorithm nds the best solution. This indicator is related to both the feasibility and the makespan in each instance of the algorithm. Results are shown in Table 1 for each algorithm. On the one hand, regarding the constructive heuristics, the BICH algorithm nds the best solution in 14 instances as compared to the 9 and 1 of the A _ NEH and FL heuristics respectively. On the other hand, the ANP A algorithm is clearly the best with 69, 80, 112 and 491 instances for the stopping criteria t= 2,5,20 and 60 respectively. Regarding the other two iterative improvement algorithms ( A _ IGA and GA ), A _ IGA slightly improves GA for each stopping criterion. Taking into account this indicator, the ecient algorithms would be BICH , ANP A(t= 2) , ANP A(t= 5) , ANP A(t= 20) and ANP A(t= 60) , as shown in Figure 9. 5.4 Summary Although the determination of the best algorithms for the problem under study is not trivial due to the existence of infeasible solutions, the proposed algorithms BICH and ANP A have been found to be the most ecient algorithms for each one of the three indicators analysed. With respect to the rest of the algorithms, it is not clear whether A _ IGA outperforms GA or vice versa, since the latter A _ IGA is better for the last two indicators, but nds less feasible solutions. The same happens when comparing FL and A _ NEH . Regarding the computational time requirements, note that the average CPU time is both instance- and instance-size- dependent indicator (see Fernandez-Viagas and Framinan, 2015b). However, similar results are also found when using a dimensionless time indicator. Particularly, the average relative percentage computation time described in Fernandez-Viagas and Framinan (2015b) has been tested with similar results. In order not to excessively increase the extension of the paper, these ndings are not detailed. 24 Figure 9: Amount of best solutions vs Average CPU times for each algorithm. X-Axis and Y-Axis are shown in logarithmic scale. 25