The Permutation Flow Shop Scheduling Problem with Human Resources: MILP Models, Decoding Procedures, NEH-Based Heuristics, and an Iterated Greedy Algorithm

  Citation: Fernandez-Viagas, V.; Sanchez-Mediano, L.; Angulo-Cortes, A.; Gomez-Medina, D; Molina-Pariente, J.M. The Permutation Flow Shop Scheduling Problem with Human Resources: MILP Models, Decoding Procedures, NEH-Based Heuristics, and an Iterated Greedy Algorithm. Mathematics 2022,10, 3446. https:// doi.org/10.3390/math10193446 Academic Editor: Hsien-Chung Wu Received: 27 July 2022 Accepted: 19 September 2022 Published: 22 September 2022 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). mathematics Article The Permutation Flow Shop Scheduling Problem with Human Resources: MILP Models, Decoding Procedures, NEH-Based Heuristics, and an Iterated Greedy Algorithm Victor Fernandez-Viagas 1,2,* , Luis Sanchez-Mediano 1, Alvaro Angulo-Cortes 1, David Gomez-Medina 1,2 and Jose Manuel Molina-Pariente 1,2 1Department of Industrial Organization and Business Management I, School of Engineering, University of Seville, 41091 Seville, Spain 2 Grupo de Informática de la Salud Computacional, Hospital Universitario Virgen del Rocío, 41013 Seville, Spain *Correspondence: [email protected] Abstract: In this paper, we address the permutation flow shop scheduling problem with sequencedependent and non-anticipatory setup times. These setups are performed or supervised by multiple servers, which are renewable secondary resources (typically human resources). Despite the real applications of this kind of human supervision and the growing attention paid in the scheduling literature, we are not aware of any previous study on the problem under consideration. To cover this gap, we start theoretically addressing the problem by: proposing three mixed-integer linear programming models to find optimal solutions in the problem; and proposing different decoding procedures to code solutions in approximated procedures. After that, the best decoding procedure is used to propose a new mechanism that generates 896 different dispatching rules, combining different measures, indicators, and sorting criteria. All these dispatching rules are embedded in the traditional NEH algorithm. Finally, an iterated greedy algorithm is proposed to find near-optimal solutions. By doing so, we provide academics and practitioners with efficient methods that can be used to obtain exact solutions of the problem; applied to quickly schedule jobs and react under changes; used for initialisation or embedded in more advanced algorithms; and/or easily updated and implemented in real manufacturing scenarios. Keywords: scheduling; flow shop; MILP; decoding procedure; makespan; flow shop; human resources; multiple servers; sequence-dependent setups; iterated greedy MSC: 90B35; 68M20; 90C59 1. Introduction Minimising the makepan in the Permutation Flow shop Scheduling Problem (PFSP) is one of the most addressed problems in the Operations Research literature [ 1 ]. In the PFSP, there are n jobs, each one with m operations, and a set of m machines. Each machine is exclusively responsible for processing an operation of each job and, therefore, each job must be processed on all machines following the same route. The goal of this scheduling problem is to find the best sequence (the same in all machines) to process the jobs according to a certain objective function, with the makespan being the most common one in the literature. In this problem, machines are typically classified as primary resources that remain busy throughout the processing of each job. In addition, in several situations, the processing of jobs may also require the use of other resources (such as, e.g., raw materials, human resources, or setup tools), denoted as ‘secondary’ ([ 2 ]). Among these types of secondary resources, the use of servers is very common in the real manufacturing industry [ 3 ]: examples can be found in flexible manufacturing systems [ 4 ]; in scheduling problems with versatile machines and assembly components [ 5 ]; in computer-controlled material Mathematics 2022,10, 3446. https://doi.org/10.3390/math10193446 https://www.mdpi.com/journal/mathematics Mathematics 2022,10, 3446 2 of 32 handling systems [ 6 ]; in biomass truck scheduling problems in the context of supply chain optimisation [ 7 ]; and in container terminals [ 8 ]. Servers are secondary resources that are in charge of carrying out setups between jobs or families of jobs. In this sense, these servers could represent a robot (see, e.g., [ 4 ]), a human (see, e.g., [ 9 ]), or an automatic guided vehicle ([ 10 ]), among others. In the literature, this type of resource has also been denoted as a setup operator (see, e.g., [ 11 , 12 ]). In addition to the previous real implications of the problem, this scheduling problem is especially appropriate for the U-shaped manufacturing layout (see [ 13 ] for a definition of this production line layout), where machines are distributed in a flow-shop layout with a U shape, and human resources are located in the centre in charge of carrying out the changeovers. Despite the practical applications of servers in real scenarios and the recent interest in the literature to these kind of scheduling problems with secondary resources responsible for setup times (see, e.g., [ 3 , 9 , 14 – 16 ]), literature addressing this topic is still very scarce and there are many open research questions. Most papers address the problem with a single server (see, e.g., [ 6 , 17 ]) under very simple layouts, mostly identical parallel machines (see, e.g., [ 18 – 20 ]). The case with multiple human resources has been tackled only in unrelated parallel machines, hybrid flow shop, open shop, and no-wait flow shop. Although these papers propose relevant algorithms for solving such scheduling problems, there is still a need to further improve the knowledge of multiple servers, especially solving and analysing its influence in the permutation flow shop scheduling problem, the most relevant and addressed layout in the literature. To cover this gap, we address the Permutation Flow shop Scheduling problem with Multiple Servers, denoted as PFSMS in the following. Regarding the setup times which have to carry out these servers, according to [ 21 ], they may be classified first as non-anticipatory and anticipatory: the former (also denoted as inseparable or attached) must be executed after the job arrival; the latter (also denoted as separable or detached) can be performed at the machine at any time after the completion of the previous job (i.e., they may be performed before the job arrival). Second, setup times can be classified as sequence-dependent or sequence-independent if this amount of time depends or not on the previous job executed on the same machine, respectively. Following the common approach used in the literature, in this paper, we address the PFSMS with non-anticipatory and sequence-dependent setup times. This scheduling problem, with the objective of minimising the makespan, is denoted by Fm , Sr|prmu|Cmax , according to [ 18 , 22 ]. This problem is clearly NP-hard, since the same problem without considering the servers is already NP-hard [23]. To deal with this PFSMS for the first time, the contributions of this paper can be stated as follows: • A new PFSP problem is defined with multiple servers or human resources in charge of carrying out setups. • We identify efficient formulations to solve the proposed problem, by developing three different Mixed Integer Linear Programming (MILP) models. • Another contribution is the development of decoding procedures. These procedures are important to obtain feasible schedules and to reduce the solution space of the problem. Both issues are essential to guarantee the efficiency of approximated algorithms developed for the problem. • A new procedure to generate static dispatching rules is proposed. Using this procedure, we propose and compare 896 different dispatching rules. • In addition, we propose different NEH-based constructive heuristics to efficiently solve the problem. Using the traditional NEH heuristic ([ 24 ]), we analyse the influence and efficiency of its first phase, embedding all previous dispatching rules. In this regard, note that heuristics have been traditionally proposed in the literature either to obtain high-quality solutions in short times (required when decisions should be made almost instantaneously or with high computational requirements) or to provide good initial solutions for more advanced algorithms (as metaheuristics). Recently, they are also relevant for quickly reacting to the changing environments of industry 4.0 ([ 25 ]). Mathematics 2022,10, 3446 3 of 32 • Finally, starting with the best constructive heuristic, we develop an iterated greedy metaheuristic to find near-optimal solutions in large-sized instances. Note that with these contributions, we try to introduce future researchers and practitioners with the problem under study, providing them with tools that can be directly applied, combined or adapted to other related scheduling problems. It is also noteworthy to mention that, as most proposals in the literature, we propose deterministic approaches to solve the problem under consideration. Despite the inherent uncertainty in the supply chain ([ 26 – 28 ]), in many scheduling situations, deterministic approaches are very robust to stochasticity and uncertainty (see in this regard [25,29]). The rest of the paper is organised as follows: the problem is described in Section 2. In this section, we also review the literature related to the problem under study. In Section 3, we propose several new formulations for the problem using MILP models. The decoding procedures are explained in Section 4. Regarding approximate algorithms, the proposed dispatching rules and constructive heuristics are detailed in Section 5, while the iterated greedy algorithm is explained in Section 6. The computational results of all previous methods are shown in Section 7. Finally, the conclusions and future research lines are discussed in Section 8. 2. Problem Description and Background In the PFSP, there is a set N of n jobs that must be processed on a set M of m machines, following the same route of machines for each job. Each machine is always available and has to process jobs (one by one) following a certain sequence, Π= (π1 , . . . , πn) , which is identical for all machines (permutation constraint). Let Oij denote the operation of job j processed on machine i . Whenever it does not cause confusion, let Oi[k] be the operation corresponding to the job in position k on machine i (i.e., Oiπk ). Each operation Oij has a processing time pij . In addition, each job j requires a setup time silj , when it is processed after job l on machine i . This setup is both non-anticipatory and sequence-dependent, and has to be carried out on that machine i by a worker w . Let R be the set of r identical workers that can perform all setups. Denoting by Cij the completion time of job j on machine i , the goal of the problem is to find the best schedule that minimises the maximum completion time, Cmax, where Cmax is equal to Cmπn. In Figure 1, we show an example of Fm , Sr|prmu , silj|Cmax with four jobs, four machines, and two workers. We can observe, for example, that the setup operation O3[1] cannot be performed before, as all the workers are unavailable. Similarly, some idle time is also forced on machines 3 and 4 by the setup of operations O3[2] and O4[2] , which have to wait until a worker is available. As mentioned in the previous section, due to its importance, the permutation flow shop scheduling problem is one of the most active problems in the operation literature with hundreds of contributions in the last years, addressing different variants and constraints in the classical problems. Recent examples solving the permutation flow shop can be found in [ 30 – 36 ], and addressing different setups configurations in [ 37 – 44 ]. For comprehensive reviews of the problem under different constraints, we refer the interested reader to [1,45–49] . Despite the extensive flow-shop-based literature, we are not aware of any previous reference to the problem under study (PFSMS) so far. Therefore, we focus this review of the literature on scheduling problems that deal with single or multiple servers. Regarding the single server literature, most focus on solving the parallel machine scheduling problem and applying it . In this regard, this problem has been solved for two machines, for example [ 3 , 4 , 50 ] considering sequence-independent setup times and using objective functions based on the makespan or total idle time. The problem with m machines has been solved by [ 6 ] sequence-dependent setup times and makespan minimisation and by [ 17 ] for sequence-dependent setup times to minimise total weighted earliness and tardiness. The complexity of this kind of problem is addressed by [10]. The problem with a single server has also been solved by [ 51 ] for the parallel dedicated machine scheduling problem with sequence-dependent setup times and makespan minimisation, while [ 20 ] have proposed Mathematics 2022,10, 3446 4 of 32 a MILP model and two approximate algorithms to solve the unrelated parallel scheduling problem with sequence-dependent setup times and machine eligibility restrictions. Regarding other layouts in the literature with a single server, the problem has also been solved for the flow shop scheduling problem with two machines and minimisation of the makespan by [ 52 ]. The authors have addressed the problem of sequence-independent and anticipatory setup times. The complexity of this problem is addressed by [ 53 ], considering unit processing times, and, more generally, by [ 54 ]. Using the same configuration of setup times and traditional processing times, [ 55 ] have solved the no-wait variant of the problem to minimise the total completion time. Anticipatory and non-anticipatory setup times are addressed by [ 56 ], minimising the makespan, and also considering dismounting times. The problem with m machines and a single server has been solved by [ 5 ] for total completion times’ minimisation, but in this case, the server is responsible for both processing and setup operations. The authors have addressed the problem with non-anticipatory and sequence-independent setup times. Figure 1. Example of the Fm,Sr|prmu,sij|Cmax problem. Regarding the related scheduling problem with multiple servers, [ 18 ] have addressed the identical parallel machine scheduling problem with sequence-independent setup times. For the unrelated parallel machine layout, [ 15 ] have proposed a Grasp algorithm to solve the problem with sequence-dependent setup times for makespan minimisation, while [ 16 ] have proposed an iterated greedy-based algorithm to solve a bi-objective variant of the problem. Both researchers have addressed the case where a setup requires more than one server at the same time. Based on the ceramic tile manufacturing sector, [ 19 ] have solved a related unrelated parallel machine layout, where the setup times depend on the assignment of resources. The hybrid flow shop scheduling problem with sequenceindependent setup times and multiple servers has recently been addressed by [ 9 , 14 ]. The former has proposed a backtracking search optimisation algorithm to solve the anticipatory variant, while the latter have proposed several constructive and composite heuristics to solve both the anticipatory and non-anticipatory cases. An example in the open shop layout with two machines can be found in [ 57 ] for non-anticipatory and sequence-independent setup times. Regarding the flow shop layout with multiple servers, the problem has been solved only by [ 58 , 59 ]. However, they addressed the no-wait variant using anticipatory and sequence-dependent setup times by proposing a genetic algorithm. A summary of the literature review is presented in Table 1. Mathematics 2022,10, 3446 5 of 32 Table 1. Summary of the literature review. Paper #Machines Layout #Workers Setup Other Objective Approach m=2m≥2 FS O r=1r≥1 SD SI A NA Constraint Cmax ∑CjOM MO MILP H MH O [3] X P X X - - X X [4] X P X X - - X X [5] X X X X X SPS X X X [6] X P X X - - X X X [9] X HF X X X X X X [14] X HF X X X X X X [15] X R X X - - X X X X [16] X R X X - - X X [50] X P X X - - X X X [17] X P X X - - X X [51] X DP X X - - X X X [20] X R X X - - E X X X [52] X X X X X CPT X X X [55] X X X X X X X X [56] X X X X X X DT X X [57] X O X X X X X [19] X R X X - - STAR X X X [58] X X X X X NW X X [59] X X X X X NW X X This X X X X X X X X X Notation: FS, Flow shop; P, Identical Parallel Machines; DP, Dedicated Machines; R, Unrelated Parallel Machines; O, Other; NW, no-wait; STAR, setup times depending on the number of assigned resources; SPS, servers must process the jobs and carry out the setup; DT, Dismounting Times; CPT, Constant Processing Times; E, Eligibility; SD, Sequence dependent setup times; SI, sequence independent setup times; A, Anticipatory setup times; NA, Non-anticipatory setup times; OM, Other Mono-objective; MO, Multiobjective; H, Heuristic; MH, Metaheuristic. Mathematics 2022,10, 3446 6 of 32 3. MILP Models In this section, we elaborate three different MILP models to exactly solve the problem under consideration. The parameters and variables used by the models are presented in Section 3.1. Then, a detailed description of each model is shown in Sections 3.2–3.4 for Models 1, 2, and 3, respectively. 3.1. Notation: Variables and Common Parameters The three proposed MILP models use the following common indices and parameters: •i(i∈ M): machine index. •j(j∈ N ): job index. •k(k=1, . . . , n): position index. •w(w∈ R): worker index. •sijj0 ( i∈ M , j , j0∈ N ): setup time of job j0 on machine i when it is processed immediately after job j . If the job j0 is the first job of the sequence, its setup time is defined by si,0,j0. •pij (i∈ M,j∈ N ): processing time of job jon machine i. The proposed models use some of the following variables: •Xijw ( i∈ M , j∈ N , w∈ R ): 1 if the setup associated with the operation Oij is carried out by the worker w. •αjkj0 ( j , j0 , k∈ N ): 1 if job j is assigned to position k and is processed before job j0 , 0 otherwise. Additionally, let α0,0,j0= 1 if job j0 is the first job in the sequence (0 otherwise). •βjj0(j,j0∈ {0} ∪ N ): 1 if job jis processed before job j0, 0 otherwise. •Zjk (j∈ N ,k∈ N ): 1 if job jis assigned to position k, 0 otherwise. •γikw ( i∈ M , k∈ N , w∈ R ): 1 is setup of the job in position k is performed by worker won machine i. •Cij (i∈ M,j∈ N ): Completion time of job jon machine i. •Cmax: Makespan. •ˆ Cik (∈ M,k∈ N ): Completion time of the job in position kon machine i. •Bij (i∈ M,j∈ N ): Starting time of job jon machine i. •ˆ Bik (i∈ M,k∈ N ): Starting time of the job in position kon machine i. •ηwiji0j0 ( w=∈ R , i , i0∈ M , j , j0∈ N ): 1 if worker w performs the setup of operation Oij immediately before operation Oi0j0 . Furthermore, let ηw00ij = 1 if operation Oij is the first operation in the worker w. •λiji0j0 ( i , i0∈ M , j , j0∈ N ): if the setup for the operations Oij and Oi0j0 are performed on the same machine, then it takes 1 if the setup of Oi0j0 is carried out after the setup of Oij, and 0 otherwise. •ˆ λiki0k0 ( i , i0∈ M , k , k0∈ N ): if the setup for the operations Oi[k] and Oi0[k0] is performed on the same machine, then it takes 1 if the setup of Oi0[k0] is carried out after setup of Oi[k], and 0 otherwise. 3.2. Model 1 Our first proposal, denoted Model 1, emerges from the Wilson model, introduced by [ 60 ] in the classical Fm|prmu|Cmax , which belongs to the Wagner family. This model used Zjk to define the sequence of jobs in the shop. Then, this auxiliary variable is used to calculate the starting time of the job in position k ( ˆ Bik ). In our proposal, we introduce the setup times and workers that are defined by the decision variables αjkj0 , γikw , and ˆ λiki0k0 . Using the previous parameters and variables, Model 1 can be formulated as follows: Mathematics 2022,10, 3446 7 of 32 Min ˆ Bmn +∑ j∈N pm,jZj,n subject to n ∑ k=1 Zjk =1j∈ N (1) ∑ j∈N Zjk =1k∈ N (2) ∑ j0∈N −{j} αjkj0=Zjk j∈ N ,k∈ N − {n}(3) ∑ j0∈N −{j} αj0,k−1,j=Zjk j∈ N ,k∈ N − {1}(4) α0,0,j=Zj,1 j∈ N (5) ˆ B1,k+∑ j∈N p1,jZj,k+∑ j∈N ∑ j0∈N −{j} (s1,j,j0·αjkj0)≤ˆ B1,k+1k∈ N − {n}(6) ˆ Bi,1 +∑ j∈N pi,jZj,1 +∑ j0∈N (si+1,0,j0·α0,0,j0)≤ˆ Bi+1,1 i∈ M − {m}(7) ˆ B1,1 =∑ j0∈N (s1,0,j0·α0,0,j0) (8) ˆ Bi,k+∑ j∈N pijZjk +∑ j∈N ∑ j0∈N −{j} (si+1,j,j0·αj,k−1,j0)≤ˆ Bi+1,ki∈ M − {m},k∈ N − {1}(9) ˆ Bi,k+∑ j∈N pijZjk +∑ j∈N ∑ j0∈N −{j} (sijj0·αjkj0)≤ˆ Bi,k+1i∈ M − {1},k∈ N − {n}(10) ∑ w∈R γikw =1i∈ M,k∈ N (11) ˆ Bi0k0−ˆ Bik +V(3−γikw −γi0k0w−ˆ λiki0k0)≥∑ j∈N ∑ l∈N −{j} (si0jj0αj,k0−1,j0),w∈ R,i0∈ M,i>i0,k∈ N ,k0>k,(12) ˆ Bik −ˆ Bi0k0+V(2−γikw −γi0k0w+ˆ λiki0k0)≥∑ j∈N ∑ j0∈N −{j} (sijj0αj,k−1,j0),w∈ R,i0∈ M,i>i0,k∈ N ,k0>k,(13) ˆ Bi0k0−ˆ Bik +V(2−γikw −γi0k0w)≥∑ j∈N ∑ l∈N −{j} (si0jj0αj,k0−1,j0),w∈ R,i∈ M,i0≥i,k∈ N ,k0>k,(14) ˆ Bi0k−ˆ Bik +V(2−γikw −γi0kw)≥∑ j∈N ∑ l∈N −{j} (si0jj0αj,k0−1,j0),w∈ R,i∈ M,i0>i,k∈ N ,(15) The set of constraints (1) ensures that every job is assigned to a unique position, while the constraints (2) achieve that exactly one job is assigned to each position. The set (3) establishes that each job has a successor, whereas the predecessors are defined in constraints (4) and (5). The starting time of every job on the first machine is defined in the set (6). Analogously, the starting time of the first job on every machine is calculated using constraints (7) and (8). Constraints (9) establish that a job must start after its previous operation is completed. The set (10) ensures that a job starts once its previous job in the same machine has been processed. Constraints (11) ensure that every job is assigned to a worker. The set of constraints (12) and (13) achieves that the same worker does not process two setups at the same time. Finally, on the one hand, constraints (14) establish that an operation Oi[k] has to start before an operation Oi0[k0] if k0>k and i0≥i . On the other hand, (15) enforces that operation Oi[k]must start before Oi0[k]if i0>i. Mathematics 2022,10, 3446 8 of 32 3.3. Model 2 Our second proposal is denoted by Model 2. This model belongs to the Manne family of models ([ 60 ]), also denoted as a disjunctive formulation. Some examples of this formulation in flow-shop layouts can be found in the SGST model proposed by [ 60 ], in the MILP model proposed by [ 9 ] or in the MIP1 model proposed by [ 5 ]. This family starts by defining the sequence of jobs by variable βjj0 . In our case, we incorporate variable λiji0j0 for the sequences in the workers. In doing so, Model 3 can be formulated as follows: Min Cmax subject to ∑ j∈{0}∪N −{j0} βjj0=1j0∈ N (16) ∑ j0∈N −{j} βjj0≤1j∈ {0} ∪ N (17) ∑ w∈R Xijw =1i∈ M,j∈ N (18) λiji0j0+λi0j0ij ≥Xijw +Xi0j0w−1i∈ M,i0∈ M − {i},j∈ N ,j0∈ N − {j},w∈ R (19) Bij0≥Cij −V·(1−βjj0)i∈ M,j∈ {0} ∪ N ,j0∈ N − {j}(20) Bi0j0≥Bij +∑ j00{0}∪N −{j} (sij00 j·βj00 j)−V·(1−λiji0j0)i∈ M,i0∈ M − {i},j∈ N ,j0∈ N − {j}(21) Cij ≥Bij +∑ j0{0}∪N −{j} (sij0j·βj0j) + pij i∈ M,j∈ N (22) Ci−1,j≤Bij i∈ M − {1},j∈ N (23) Cmax ≥Cm,jj∈ N (24) The set of constraints (16) ensures that each job has a predecessor. Note that a dummy job is introduced to represent the predecessor of the first job in the sequence. Furthermore, each job (including the dummy job) has at most one successor, which is bounded by constraints (17). The set (18) enforces that the setup of operation Oij is assigned to a unique worker. Constraints (19) define binary variables λiji0j0 . The sets (20) and (21) ensure that an operation has to start processing after its predecessor in the machine and after its predecessor in the worker, respectively. The set of constraints (22) defines the completion time of an operation and (23) ensures that an operation is not processed after its previous operation (in the previous stage) is completed. Finally, the set (24) defines the makespan of the sequence. 3.4. Model 3 Our last MILP model arises from Model 2, but modifies how workers are treated. In this case, we introduce variable ηwiji0j0 to define the sequence followed by each worker (the use of a related decision variable can be found in the first MILP model proposed in [ 61 ]). The formulation of this model is presented below: Mathematics 2022,10, 3446 9 of 32 Min Cmax subject to ∑ j∈{0}∪N −{j0} βjj0=1j0∈ N (25) ∑ j0∈N −{j} βjj0≤1j∈ {0} ∪ N (26) C1,j≥p1,j+∑ j0{0}∪N −{j} (s1j0j·βj0j)j∈ N (27) Cij −Ci−1,j≥pij +∑ j0{0}∪N −{j} (sij0j·βj0j)i∈ M − {1},j∈ N (28) Cij0−Cij +V·(1−βjj0)≥pij0+sijj0i∈ M,j∈ {0} ∪ N ,j0∈ N − {j}(29) ∑ w∈R ηw00i0j0+∑ i∈M ∑ j∈N ηwiji0j0!=1i0∈ M,i∈ M − {i0},j0N,j∈ N − {j0}(30) ∑ w∈R ∑ i0∈M−{i} ∑ j0∈N −{j} ηwiji0j0≤1i∈ M,j∈ N (31) ∑ i∈M ∑ j∈N ηw00ij =1w∈ R (32) ηw00ij +∑ i0∈M−{i} ∑ j0∈N −{j} ηwi0j0ij ≥∑ i0∈M−{i} ∑ j0∈N −{j} ηwiji0j0i∈ M,jN,w∈ R (33) Ci0j0−pi0j0+V˙ (1−ηwiji0j0)≥Cij −pij +∑ j00∈{0}∪N −j0 (si0j00 j0·βj00 j0)w∈ R,i∈ M,i0∈ M − {i},j∈ N ,j0∈ N − {j}(34) Cij −pij +V˙ (1−ηw00ij)≥∑ j0∈{0}∪N −j (sij0j·βj0j)w∈ R,i∈ M,j∈ N (35) Cmax ≥Cm,jj∈ N (36) As in the previous model, constraints (25) and (26) fully define variables βjj0 . The completion times of the first machine are established in the set (27), while the constraints (28) and (29) bound the completion times of each job in the other machines. To do so, constraints (28) limit these completion times considering its previous machine, while constraints (29) consider its previous job. Variable ηwi0j0ij is defined in Equations (30)–(33). More specifically, the set (30) ensures that the setup of each operation has either a predecessor in the same worker or a dummy setup. Constraints (31) achieve that each setup has at most a successor in a worker, while constraints (32) and (33) ensure that each worker starts with a setup. The completion times of the jobs and the precedence relationships between the operations are linked in constraints (34) and (35). Finally, the set of constraints (36) defines the makespan. 4. Proposed Decoding Procedures and Complete Enumeration Due to the limitation of MILP models to solve medium-large sized instances, we also propose different approximate algorithms to tackle the problem (described in Section 5). Each of the proposed algorithms uses a sequence of jobs to easily represent the solutions, since each machine must process the jobs in the same order (permutation constraint, see Section 2). This sequence can be formally defined as follows: • Sequence of jobs, Π= (π1 , . . . , πk , . . . , πn) : It represents the order that each machine has to follow to process the jobs, i.e., πk has to be processed before job πk+1 in any machine i(k∈ {1, . . . , n−1}). Obviously, for the problem under consideration, the same sequence of jobs could theoretically lead to different semi-active schedules depending on how and in which order the operations of the jobs are assigned to the workers. Therefore, in the proposed algorithms, it is necessary to include mechanisms to construct a unique schedule to represent the sequence. In this regard, the procedure to find a specific schedule (and consequently an objective function value) from a representation of solutions is denoted by decoding procedure . For the problem under consideration, a decoding procedure can be determined by establishing the following rules: Mathematics 2022,10, 3446 16 of 32 7. Computational Results In this section, we present the computational evaluations carried out in our study. A total of four experiments have been performed: first, we compare the performance of the proposed MILP models in Section 7.2; second, the proposed decoding procedures are compared in Section 7.3 by analysing all their complete solutions (complete enumeration procedure); thirdly, the efficiency of each proposed dispatching rule is tested in Section 7.4; fourthly, the constructive heuristics proposed are analysed in Section 7.6; then, the iterated greedy algorithm is calibrated and its performance is evaluated in Section 7.6. To perform all these experiments, we generate three sets of benchmarks explained in Section 7.1. All procedures tested have been compared using the same computer conditions: on the same Intel Core i7-3770 with 3.4 GHz and 16 GB RAM; under the same programming languages (C# and Gurobi 9.5.0); and by the same person and using the same common functions and libraries. Finally, a sensitivity analysis is performed in Section 7.7 to study the impact of each factor on the problem. 7.1. Benchmarks In this section, we describe three sets of instances generated to test the performance of the proposals. The first set, denoted by β1 , is specifically constructed to compare exact methods, i.e., the proposed MILP models (see Section 3) and the complete enumerations for each of the proposed decoding procedure (see Section 7.3). The second set, denoted by β2 is constructed to compare approximate methods, i.e., the proposed dispatching rules (see Section 5.1), the proposed NEHV (see Section 5.2), and to evaluated the performance of the proposed IGV (see Section 6). Finally, a last benchmark, denoted by β3 , is generated to calibrated IGV. •β1 for comparison of exact methods. As the complete enumeration procedure has a high computational requirement (due to the evaluation of all solutions), the set is composed of 540 small-sized instances with processing and setup times following uniform distributions from [ 1, 99 ] and [ 1, γ] , respectively ([ 14 , 91 , 92 ]). γ is a parameter of the testbed in the range γ∈ { 25, 50, 100 } ([ 14 ]). The other parameters are generated for the following values: n∈ { 4, 5, 6, 7, 8, 9 } , m∈ { 3, 4, 5 } , and r∈ { 1, 2 } . Finally, five instances are constructed for each combination of the parameters. •β2 for comparison of approximate methods: β2 is a set of medium and large instances. This benchmark is generated following a procedure similar to β1 . Thereby, pij and silj follow the uniform distributions [ 1, 99 ] and [ 1, γ] , respectively. Regarding the parameters of the testbed, five instances are generated for each combination of the following parameters: γ∈ { 25, 50, 100 } , n∈ { 50, 100, 150, 200 } , m∈ { 10, 20 } , and r=r0·m, with r0∈ {0.1, 0.2, 0.3, 0.4, 0.5}. •β3 for calibration of the iterated greedy algorithm. Similarly as in β2 , we generate two instances for each combination of the following parameters: γ∈ { 25, 50, 100 } , n∈ {50, 100, 150, 200},m∈ {10, 20}, and r=r0·m, with r0∈ {0.1, 0.2, 0.3, 0.4, 0.5}. 7.2. Computational Evaluation: Comparison of MILP Models In this section, we test the performance of the MILP models proposed in Section 3. The proposed models are solved using Gurobi 9.5.0 with 500 seconds as the stopping criterion in each instance of the β1 set. The models are then compared in terms of the quality of the solution by the Relative Percentage Deviation ( RPDpi ) indicator (see Equation (37)). RPDpi measures (in percentage) the difference between the objective function found by procedure p in instance i ( OFhi ) and the best solution found in this instance (denoted besti ). In addition, to compare the computational effort required by the MILP models, we use the average CPU time. Both indicators are summarised in Table 2, grouped by n , m , r , and γ . Finally, in Table 3, we also present the number of optimal, feasible, and no solutions found. RPDpi =OFhi −besti besti 100 (37) Mathematics 2022,10, 3446 17 of 32 According to these computational results, we can observe that the best performance is found in Model 1. This model clearly outperforms the other MILP models in terms of average RPD (denoted ARPD) and CPU times. In this regard, Model 1 has an ARPD value of 0.01 (in the whole benchmark), while Models 2 and 3 found 0.14 and 0.56, respectively. This difference is statistically significant with a p -value of 0.004, using a non-parametric Mann-Whitney test between Model 1 and 2. With respect to average CPU times, we can observe the same trend. Model 1 requires a lower average CPU (61.61 seconds) as compared to Model 2 (84.78) and Model 3 (214.48). Similarly, Models 1 and 2 clearly outperform Model 3 in terms of number of optimal solutions, both finding a total of 487 optimal solutions versus 332 by Model 3. Table 2. Computational results of the proposed MILP models: ARPD values and CPU times (in seconds). Parameter Value ARPD CPU Time (s) Model 1 Model 2 Model 3 Model 1 Model 2 Model 3 n 4 0.00 0.00 0.00 0.24 0.27 2.03 5 0.00 0.00 0.00 2.54 0.66 19.90 6 0.00 0.00 0.35 19.18 6.24 75.67 7 0.00 0.00 1.01 63.08 22.66 272.94 8 0.04 0.21 2.68 112.70 86.66 421.57 9 0.02 0.61 3.65 171.94 392.16 494.75 m 3 0.00 0.00 0.00 6.85 50.19 134.54 4 0.00 0.06 0.25 56.17 84.77 221.45 5 0.03 0.35 1.60 121.82 119.36 287.44 r1 0.01 0.26 0.55 91.13 104.58 213.69 2 0.01 0.01 0.57 32.10 64.97 215.26 γ 25 0.00 0.00 0.02 5.30 71.31 204.16 50 0.01 0.06 0.42 54.59 79.43 207.71 100 0.02 0.35 1.18 124.96 103.59 231.55 Average 0.01 0.14 0.56 61.61 84.78 214.48 Table 3. Computational results of the proposed MILP models: optimal solutions (#O), feasible solutions (#F), and no solution is found (#N). Parameter Value Model 1 Model 2 Model 3 #O #F #N #O #F #N #O #F #N n 4 90 0 0 90 0 0 90 0 0 5 90 0 0 90 0 0 88 2 0 6 87 3 0 90 0 0 85 5 0 7 82 8 0 88 2 0 47 15 28 8 73 17 0 84 6 0 19 5 66 9 65 25 0 45 45 0 3 10 77 m 3 179 1 0 176 4 0 138 4 38 4 164 16 0 164 16 0 108 9 63 5 144 36 0 147 33 0 86 24 70 r1 229 41 0 232 38 0 166 23 81 2 258 12 0 255 15 0 166 14 90 γ 25 180 0 0 165 15 0 114 2 64 50 165 15 0 166 14 0 114 10 56 100 142 38 0 156 24 0 104 25 51 All 487 53 0 487 53 0 332 37 171 Mathematics 2022,10, 3446 18 of 32 7.3. Computational Evaluation: Complete Enumeration Applying Different Decoding Procedures In this section, using the same set of instances β1 , we compare the decoding procedures proposed in Section 4. Each procedure is embedded in a complete enumeration procedure and the best solution found is compared using the RPD indicator (see Equation (37), where, again, besti is the best solution found in instance i , i.e., considering both MILP models and the complete enumeration procedures). The computational results are demonstrated in Table 4. The best ARPD (2.32) is found using the PRS SI procedure of the PRS family. In this family, the tie-breaking mechanisms have a high influence on reducing the ARPD value from 3.76 ( PRS SD procedure) to 2.32, while these mechanisms are not as relevant in the PRC family, oscillating the ARPD from 7.65 to 7.76. Despite the excellent performance found by PRS SI , it is not statistically significant with PRS M (the second best proposal), according to the non-parametric Mann-Whitney test ( p -value equals to 0.221). However, with the exception of PRS M and PRS SPD (whose ARPD is similar to 2.64 and 2.65, respectively), it is statistically significant with all other decoding procedures (in this regard, when compared with the PRS SPI procedure, the p-value found in the Mann-Whitney test is 0.000). Table 4. Computational results in small-sized instances: ARPD values of the proposed decoding procedures, the best NEHV (i.e., NEHV(PMS,ABS,D)), and IGV. Parameter Value PRS JPRS MPRS SPI PRS SPD PRS SI PRS SD PRC JPRC MPRC SPI PRC SPD PRC SI PRC SD NEHV IGV n 4 3.02 2.89 3.49 2.21 2.17 3.41 6.54 6.49 6.49 6.49 6.54 6.49 3.86 1.92 5 4.08 2.96 3.85 3.11 2.57 4.33 8.04 8.04 8.15 7.96 8.01 8.07 5.53 2.02 6 4.08 2.80 3.76 2.96 2.60 4.31 7.87 7.85 7.85 7.83 7.86 7.87 6.35 1.97 7 4.06 2.61 3.47 2.62 2.39 3.62 7.78 7.81 7.83 7.74 7.83 7.75 7.61 1.42 8 3.35 2.36 3.12 2.68 2.20 3.54 7.82 7.83 7.96 7.68 7.82 7.87 7.80 1.48 9 3.56 2.24 2.92 2.31 1.98 3.34 8.20 8.27 8.27 8.20 8.17 8.27 9.21 1.27 m 3 3.39 1.67 2.41 2.00 1.43 3.03 4.73 4.74 4.78 4.69 4.69 4.78 6.15 1.04 4 3.73 2.91 3.47 2.76 2.47 3.95 7.64 7.69 7.68 7.65 7.64 7.66 6.80 1.85 5 3.95 3.36 4.43 3.19 3.06 4.30 10.74 10.72 10.82 10.61 10.78 10.71 7.23 2.16 r1 5.94 4.15 5.33 4.26 3.33 6.26 12.43 12.47 12.51 12.37 12.41 12.49 8.13 2.52 2 1.44 1.14 1.54 1.04 1.31 1.26 2.99 2.96 3.01 2.93 3.00 2.95 5.32 0.85 γ 25 2.50 1.83 2.81 1.49 1.40 2.78 4.45 4.50 4.49 4.44 4.45 4.50 4.36 0.88 50 4.09 2.90 3.79 2.85 2.47 4.18 7.26 7.23 7.32 7.17 7.25 7.26 6.61 1.73 100 4.49 3.20 3.70 3.60 3.10 4.31 11.40 11.42 11.47 11.34 11.42 11.40 9.20 2.43 Average 3.69 2.64 3.43 2.65 2.32 3.76 7.71 7.72 7.76 7.65 7.71 7.72 6.73 1.68 7.4. Computational Evaluation: Dispatching Rules In this section, we compare the dispatching rules generated following the procedure described in Section 5.1. Using this procedure, we construct 896 (potentially) different sequences of jobs in each instance, based on the specific indicator, measure, and sorting criterion applied. To generate a schedule (and consequently to obtain an objective function value), we apply the best decoding procedure (i.e., the PRS SI procedure) to each proposed dispatching rule as mentioned previously. Then, each dispatching rule is tested in the set of instances β2 and its performance is evaluated by the RPD indicator (see Equation (37) ). In the same way, besti is the best solution found in instance i , considering both dispatching rules and NEHV variants. The computational results are shown in Tables 5–7. All ARPD values for the proposed dispatching rules are obtained between 23.97 and 28.75. Regarding the best results, they are found by dispatching rules using the WSUM indicator, the C sorting criterion, and any measure that considers the processing times of the operations (i.e., P, PS, PMS, PmS). Thereby, the ARPD values of this group of dispatching rules are 24.01, 23.97, 24.00, and 24.10 (for { P,WSUM,C } , { PS,WSUM,C } , { PMS,WSUM,C } , and { PmS,WSUM,C } , respectively). In fact, there is a statistically significant difference between this group of rules and { P,SUM,C } (the following best rule with an ARPD of 25.28), finding ap-value of 0.000 (using the non-parametric Mann-Whitney test). Mathematics 2022,10, 3446 19 of 32 Table 5. Computational results for the dispatching rules: ARPD values I. Measure Indicator Sorting C D H V HIH HIL LOH LOL P SUM 26.42 26.30 25.28 27.25 26.48 26.35 26.34 26.20 P WSUM 24.01 28.64 25.60 27.24 26.47 26.20 26.32 26.44 P ABS 26.40 26.23 26.23 26.57 26.34 26.28 26.37 26.27 P WABS 26.02 26.52 26.43 26.62 26.31 26.23 26.24 26.38 P SRA 26.46 26.40 26.42 26.44 26.27 26.21 26.45 26.37 P WSRA 25.98 26.58 26.30 26.40 26.29 26.28 26.27 26.34 P SRS 26.43 26.11 26.31 26.30 26.34 26.22 26.30 26.34 P WSRS 25.93 26.65 26.22 26.31 26.28 26.38 26.24 26.31 P SRN 26.22 27.25 26.62 26.42 26.49 26.61 26.40 26.50 P WSRN 28.13 26.13 26.55 26.92 26.13 26.26 26.50 26.47 P SRAN 26.15 27.17 26.33 26.50 26.45 26.47 26.27 26.33 P WSRAN 27.34 26.53 26.15 26.89 26.01 25.96 26.78 26.89 P SRSN 26.23 27.16 26.23 26.51 26.46 26.42 26.32 26.30 P WSRSN 27.31 26.55 26.46 27.12 25.92 26.06 26.90 26.92 P SRAN2 25.85 27.21 26.58 26.11 26.49 26.42 26.26 26.25 P WSRAN2 27.39 25.86 26.45 26.67 26.21 26.22 26.75 26.70 S SUM 26.33 26.34 26.25 26.23 26.16 26.21 26.32 26.34 S WSUM 26.08 26.33 26.36 26.38 26.25 26.12 26.21 26.27 S ABS 26.32 26.34 26.29 26.38 26.30 26.32 26.34 26.28 S WABS 26.40 26.14 26.28 26.40 26.18 26.22 26.21 26.23 S SRA 26.24 26.22 26.27 26.32 26.22 26.30 26.36 26.34 S WSRA 26.28 26.35 26.26 26.39 26.17 26.12 26.27 26.46 S SRS 26.24 26.35 26.31 26.19 26.20 26.18 26.41 26.34 S WSRS 26.34 26.45 26.30 26.33 26.29 26.21 26.27 26.51 S SRN 26.21 26.35 26.16 26.24 26.34 26.25 26.20 26.06 S WSRN 26.26 26.19 26.18 26.37 26.26 26.28 26.28 26.28 S SRAN 26.39 26.33 26.34 26.37 26.20 26.32 26.29 26.28 S WSRAN 26.38 26.13 26.20 26.36 26.43 26.29 26.38 26.48 S SRSN 26.37 26.18 26.33 26.51 26.30 26.33 26.45 26.35 S WSRSN 26.26 26.21 26.20 26.27 26.24 26.30 26.32 26.35 S SRAN2 26.26 26.20 26.15 26.23 26.14 26.21 26.37 26.23 S WSRAN2 26.23 26.33 26.22 26.24 26.32 26.28 26.31 26.28 MS SUM 26.30 26.18 26.32 26.37 26.33 26.24 26.21 26.34 MS WSUM 26.20 26.29 26.25 26.38 26.35 26.36 26.36 26.30 MS ABS 26.24 26.29 26.29 26.12 26.28 26.22 26.21 26.20 MS WABS 26.24 26.33 26.41 26.12 26.36 26.25 26.28 26.26 MS SRA 26.23 26.27 26.29 26.18 26.33 26.34 26.23 26.33 MS WSRA 26.33 26.44 26.29 26.10 26.18 26.23 26.13 26.15 MS SRS 26.30 26.34 26.34 26.06 26.27 26.43 26.16 26.16 MS WSRS 26.31 26.27 26.27 26.11 26.33 26.27 26.14 26.23 MS SRN 26.33 26.29 26.37 26.32 26.26 26.34 26.29 26.19 MS WSRN 26.28 26.31 26.30 26.14 26.25 26.19 26.28 26.28 MS SRAN 26.29 26.17 26.23 26.14 26.34 26.33 26.21 26.30 MS WSRAN 26.32 26.23 26.31 26.11 26.20 26.38 26.20 26.30 MS SRSN 26.29 26.35 26.21 26.19 26.31 26.25 26.27 26.27 MS WSRSN 26.18 26.30 26.36 26.12 26.20 26.08 26.20 26.25 MS SRAN2 26.23 26.28 26.31 26.19 26.27 26.36 26.23 26.20 MS WSRAN2 26.24 26.37 26.32 26.12 26.21 26.34 26.25 26.27 Mathematics 2022,10, 3446 20 of 32 Table 6. Computational results for the dispatching rules: ARPD values II. Measure Indicator Sorting C D H V HIH HIL LOH LOL mS SUM 26.46 26.51 26.27 26.39 26.37 26.38 26.22 26.18 mS WSUM 26.31 26.40 26.29 26.36 26.39 26.29 26.20 26.29 mS ABS 26.39 26.40 26.24 26.28 26.30 26.36 26.23 26.28 mS WABS 26.52 26.32 26.34 26.22 26.32 26.27 26.18 26.33 mS SRA 26.42 26.30 26.24 26.34 26.39 26.26 26.07 26.11 mS WSRA 26.44 26.39 26.23 26.29 26.22 26.17 26.18 26.30 mS SRS 26.37 26.45 26.14 26.25 26.34 26.29 26.20 26.30 mS WSRS 26.45 26.41 26.29 26.30 26.43 26.31 26.22 26.38 mS SRN 26.35 26.16 26.29 26.25 26.33 26.26 26.36 26.38 mS WSRN 26.33 26.26 26.29 26.30 26.18 26.25 26.30 26.27 mS SRAN 26.44 26.33 26.32 26.40 26.26 26.33 26.23 26.21 mS WSRAN 26.39 26.34 26.28 26.31 26.28 26.37 26.29 26.28 mS SRSN 26.32 26.34 26.22 26.33 26.24 26.35 26.30 26.22 mS WSRSN 26.40 26.36 26.28 26.23 26.34 26.37 26.17 26.23 mS SRAN2 26.38 26.38 26.27 26.35 26.43 26.33 26.19 26.27 mS WSRAN2 26.54 26.30 26.18 26.22 26.33 26.29 26.10 26.34 PS SUM 26.58 26.33 25.36 27.17 26.53 26.34 26.47 26.42 PS WSUM 23.97 28.75 25.46 27.20 26.63 26.25 26.43 26.44 PS ABS 26.38 26.31 26.30 26.59 26.16 26.40 26.31 26.30 PS WABS 26.15 26.67 26.25 26.45 26.23 26.33 26.21 26.34 PS SRA 26.38 26.34 26.37 26.38 26.24 26.33 26.26 26.34 PS WSRA 26.01 26.64 26.17 26.30 26.39 26.28 26.28 26.29 PS SRS 26.46 26.19 26.31 26.28 26.14 26.29 26.32 26.31 PS WSRS 25.93 26.58 26.34 26.42 26.39 26.41 26.21 26.23 PS SRN 26.24 27.41 26.70 26.43 26.42 26.49 26.32 26.47 PS WSRN 28.01 26.13 26.73 27.11 26.28 26.17 26.47 26.57 PS SRAN 26.07 27.26 26.53 26.61 26.44 26.54 26.19 26.29 PS WSRAN 27.26 26.49 26.22 27.01 25.88 26.16 26.90 26.82 PS SRSN 26.15 27.21 26.37 26.60 26.54 26.53 26.26 26.28 PS WSRSN 27.31 26.60 26.24 27.16 25.85 26.07 26.96 27.00 PS SRAN2 25.92 27.27 26.63 26.07 26.47 26.38 26.29 26.39 PS WSRAN2 27.34 25.87 26.59 26.91 26.13 26.18 26.58 26.71 Table 7. Computational results for the dispatching rules: ARPD values III. Measure Indicator Sorting C D H V HIH HIL LOH LOL PMS SUM 26.36 26.28 25.46 27.44 26.42 26.41 26.31 26.19 PMS WSUM 24.00 28.66 25.55 27.26 26.57 26.18 26.42 26.49 PMS ABS 26.42 26.29 26.16 26.43 26.26 26.37 26.40 26.38 PMS WABS 26.03 26.59 26.18 26.38 26.26 26.23 26.32 26.32 PMS SRA 26.55 26.31 26.50 26.46 26.16 26.24 26.28 26.24 PMS WSRA 26.02 26.71 26.25 26.31 26.25 26.28 26.20 26.31 PMS SRS 26.35 26.03 26.37 26.45 26.25 26.36 26.38 26.33 PMS WSRS 26.00 26.60 26.17 26.20 26.32 26.33 26.23 26.25 PMS SRN 26.21 27.42 26.74 26.40 26.51 26.58 26.48 26.49 PMS WSRN 28.14 26.14 26.58 26.96 26.17 26.21 26.52 26.42 PMS SRAN 26.14 27.18 26.41 26.63 26.36 26.51 26.28 26.18 PMS WSRAN 27.31 26.46 26.34 26.99 25.98 26.08 26.89 26.93 PMS SRSN 26.27 27.12 26.40 26.55 26.41 26.47 26.21 26.30 PMS WSRSN 27.39 26.53 26.36 27.22 25.98 26.06 26.81 26.88 PMS SRAN2 25.91 27.09 26.62 26.21 26.47 26.43 26.34 26.31 PMS WSRAN2 27.30 25.94 26.57 26.67 26.13 26.24 26.67 26.76 PmS SUM 26.47 26.30 25.39 27.33 26.51 26.43 26.34 26.27 Mathematics 2022,10, 3446 21 of 32 Table 7. Cont. Measure Indicator Sorting C D H V HIH HIL LOH LOL PmS WSUM 24.10 28.63 25.59 27.37 26.64 26.16 26.23 26.53 PmS ABS 26.40 26.24 26.21 26.59 26.29 26.43 26.20 26.36 PmS WABS 26.12 26.44 26.44 26.48 26.32 26.32 26.28 26.30 PmS SRA 26.50 26.24 26.35 26.26 26.29 26.30 26.25 26.23 PmS WSRA 26.03 26.76 26.27 26.39 26.27 26.26 26.34 26.29 PmS SRS 26.33 26.20 26.32 26.41 26.24 26.30 26.19 26.26 PmS WSRS 26.07 26.59 26.27 26.38 26.29 26.43 26.20 26.28 PmS SRN 26.15 27.26 26.78 26.37 26.43 26.50 26.46 26.48 PmS WSRN 28.12 26.05 26.55 26.89 26.08 26.27 26.49 26.46 PmS SRAN 26.10 27.06 26.51 26.63 26.39 26.38 26.27 26.31 PmS WSRAN 27.28 26.54 26.16 26.89 25.95 26.01 26.80 26.89 PmS SRSN 26.21 27.06 26.31 26.54 26.28 26.50 26.36 26.33 PmS WSRSN 27.35 26.61 26.37 27.08 25.97 26.07 26.79 26.83 PmS SRAN2 25.81 27.09 26.69 26.19 26.41 26.38 26.35 26.32 PmS WSRAN2 27.34 25.94 26.56 26.61 26.19 26.28 26.73 26.71 7.5. Computational Evaluation: NEHV Variants In this section, we analyse the performance of the proposed NEHV variants and compare it with the best proposals in the related literature. These 896 proposed heuristics are tested in the set of instances β2 . The computational results are shown in Tables 8–10. In addition, computational results comparing the best proposal with the most promising NEH variants in the literature are shown in Table 11 (RA, NM, KK1, KK2, AD, and PRSKE are the NEH heuristics that apply the initial order proposed by [ 74 – 79 ], respectively, while NEH represents the traditional NEH ordering the jobs in non-increasing sum of processing times). As in the previous section, the quality of the solution of the proposals is evaluated using the RPD indicator. In this case, the ARPDs of the proposals range between 2.58 and 3.30. Despite this difference, there are several variants whose ARPD is very similar, close to the best-obtained value. In this regard, the best value is found using the { PMS,ABS,D } rule as an initial sequence (with an ARPD equal to 2.58), while, e.g., the { PS,WABS,H } , { P,SRA,H } , { P,WSRS,H } , { P,WABS,HIH } , and { PS,ABS,D } found an ARPD value equal to 2.61. In this case, no statistically significant differences (with the Mann-Whitney test) were found when comparing { PMS, ABS, D } with any of the previous variants. However, the proposed NEHV( { PMS,ABS,D } clearly outperforms every previous proposal in the literature. In this regard, we obtain p values equal to 0.031 and 0.004 (using the non-parametric Mann-Whitney test) comparing our best proposal with RA and with the traditional NEH algorithm, whose ARPD values are 2.73 and 2.82, respectively.Finally, in Table 4we demonstrate the ARPD values obtained by NEHV( { PMS,ABS,D } ) in the small-sized instances, β1 . The global ARPD obtained is 6.73, which outperforms the best results obtained by several decoding procedures. Note that the best ARPD achievable using the PRS SI decoding procedure is 2.32. Mathematics 2022,10, 3446 22 of 32 Table 8. Computational results for the NEH algorithm: ARPD values I. Measure Indicator Sorting C D H V HIH HIL LOH LOL P SUM 2.80 2.82 2.74 2.79 2.83 2.79 2.73 2.68 P WSUM 2.79 3.19 2.77 2.88 2.70 2.66 2.70 2.80 P ABS 2.86 2.73 2.68 2.80 2.75 2.78 2.73 2.82 P WABS 2.84 2.65 2.73 2.70 2.61 2.71 2.84 2.67 P SRA 2.85 2.74 2.61 2.79 2.73 2.72 2.74 2.77 P WSRA 2.78 2.67 2.68 2.78 2.71 2.76 2.76 2.74 P SRS 2.86 2.67 2.68 2.73 2.71 2.64 2.74 2.75 P WSRS 2.80 2.72 2.61 2.81 2.73 2.76 2.75 2.70 P SRN 2.86 3.01 2.80 2.73 2.76 2.76 2.73 2.76 P WSRN 3.14 2.90 2.86 2.88 2.75 2.73 2.76 2.76 P SRAN 2.78 2.91 2.76 2.80 2.79 2.86 2.63 2.64 P WSRAN 3.01 2.71 2.77 2.81 2.72 2.67 2.83 2.80 P SRSN 2.77 2.88 2.71 2.79 2.81 2.78 2.65 2.61 P WSRSN 3.01 2.81 2.79 2.81 2.64 2.67 2.80 2.79 P SRAN2 2.73 2.89 2.76 2.83 2.74 2.80 2.72 2.68 P WSRAN2 3.02 2.75 2.85 2.86 2.69 2.76 2.75 2.77 S SUM 2.75 2.66 2.74 2.76 2.67 2.71 2.77 2.75 S WSUM 2.68 2.76 2.62 2.77 2.71 2.79 2.80 2.70 S ABS 2.83 2.67 2.64 2.69 2.72 2.68 2.71 2.71 S WABS 2.74 2.70 2.69 2.74 2.76 2.74 2.72 2.69 S SRA 2.74 2.69 2.75 2.73 2.75 2.73 2.73 2.78 S WSRA 2.72 2.72 2.73 2.75 2.75 2.73 2.74 2.69 S SRS 2.72 2.74 2.72 2.78 2.74 2.74 2.70 2.69 S WSRS 2.77 2.75 2.72 2.73 2.78 2.71 2.74 2.75 S SRN 2.70 2.78 2.76 2.78 2.70 2.74 2.68 2.71 S WSRN 2.76 2.73 2.68 2.77 2.72 2.70 2.75 2.72 S SRAN 2.72 2.70 2.72 2.66 2.74 2.77 2.72 2.71 S WSRAN 2.68 2.75 2.73 2.73 2.82 2.68 2.71 2.75 S SRSN 2.63 2.68 2.71 2.70 2.79 2.76 2.69 2.68 S WSRSN 2.68 2.71 2.69 2.71 2.74 2.75 2.68 2.66 S SRAN2 2.69 2.76 2.78 2.72 2.72 2.77 2.70 2.69 S WSRAN2 2.73 2.72 2.64 2.74 2.69 2.72 2.78 2.72 MS SUM 2.77 2.76 2.74 2.75 2.73 2.71 2.71 2.72 MS WSUM 2.70 2.70 2.69 2.77 2.75 2.69 2.67 2.80 MS ABS 2.72 2.79 2.70 2.71 2.70 2.74 2.70 2.73 MS WABS 2.72 2.70 2.78 2.70 2.71 2.71 2.71 2.72 MS SRA 2.78 2.73 2.69 2.78 2.71 2.72 2.71 2.65 MS WSRA 2.75 2.76 2.78 2.67 2.72 2.70 2.65 2.79 MS SRS 2.81 2.74 2.72 2.73 2.67 2.78 2.71 2.69 MS WSRS 2.72 2.77 2.79 2.67 2.73 2.79 2.72 2.68 MS SRN 2.71 2.68 2.73 2.76 2.76 2.71 2.76 2.68 MS WSRN 2.72 2.73 2.75 2.79 2.67 2.76 2.77 2.64 MS SRAN 2.76 2.70 2.72 2.70 2.73 2.70 2.69 2.73 MS WSRAN 2.77 2.73 2.74 2.72 2.77 2.70 2.74 2.71 MS SRSN 2.74 2.72 2.70 2.75 2.79 2.73 2.68 2.80 MS WSRSN 2.75 2.76 2.70 2.77 2.71 2.71 2.74 2.72 MS SRAN2 2.77 2.70 2.74 2.75 2.75 2.69 2.72 2.65 MS WSRAN2 2.72 2.77 2.74 2.63 2.71 2.72 2.77 2.70 Mathematics 2022,10, 3446 23 of 32 Table 9. Computational results for the NEH algorithm: ARPD values II. Measure Indicator Sorting C D H V HIH HIL LOH LOL mS SUM 2.73 2.65 2.72 2.75 2.76 2.79 2.80 2.73 mS WSUM 2.67 2.70 2.67 2.72 2.78 2.65 2.76 2.69 mS ABS 2.70 2.75 2.71 2.73 2.67 2.81 2.74 2.75 mS WABS 2.76 2.77 2.70 2.70 2.78 2.66 2.75 2.76 mS SRA 2.72 2.75 2.77 2.76 2.71 2.75 2.73 2.73 mS WSRA 2.66 2.73 2.70 2.72 2.76 2.74 2.70 2.70 mS SRS 2.71 2.69 2.77 2.77 2.77 2.78 2.71 2.75 mS WSRS 2.75 2.72 2.67 2.78 2.76 2.67 2.73 2.69 mS SRN 2.73 2.73 2.75 2.71 2.67 2.77 2.68 2.73 mS WSRN 2.73 2.66 2.75 2.69 2.65 2.68 2.68 2.76 mS SRAN 2.74 2.75 2.72 2.70 2.73 2.70 2.68 2.74 mS WSRAN 2.71 2.72 2.75 2.74 2.70 2.77 2.69 2.78 mS SRSN 2.67 2.75 2.80 2.76 2.76 2.76 2.71 2.74 mS WSRSN 2.70 2.77 2.72 2.75 2.75 2.75 2.77 2.71 mS SRAN2 2.67 2.75 2.78 2.67 2.73 2.77 2.76 2.77 mS WSRAN2 2.69 2.74 2.67 2.75 2.74 2.79 2.69 2.73 PS SUM 2.77 2.80 2.69 2.67 2.78 2.76 2.75 2.71 PS WSUM 2.79 3.24 2.80 2.89 2.65 2.66 2.79 2.76 PS ABS 2.85 2.61 2.67 2.92 2.73 2.75 2.78 2.70 PS WABS 2.90 2.63 2.61 2.84 2.76 2.71 2.74 2.70 PS SRA 2.88 2.74 2.62 2.78 2.68 2.70 2.76 2.70 PS WSRA 2.82 2.65 2.65 2.80 2.71 2.81 2.67 2.74 PS SRS 2.81 2.69 2.63 2.82 2.66 2.67 2.75 2.71 PS WSRS 2.83 2.70 2.68 2.80 2.77 2.75 2.73 2.68 PS SRN 2.90 2.95 2.72 2.77 2.76 2.81 2.68 2.66 PS WSRN 3.15 2.95 2.88 2.89 2.67 2.76 2.70 2.78 PS SRAN 2.87 2.92 2.74 2.78 2.80 2.79 2.72 2.75 PS WSRAN 2.97 2.75 2.81 2.87 2.76 2.70 2.78 2.76 PS SRSN 2.86 2.88 2.76 2.80 2.78 2.78 2.75 2.68 PS WSRSN 3.00 2.82 2.77 2.78 2.68 2.72 2.78 2.77 PS SRAN2 2.75 2.87 2.77 2.81 2.76 2.75 2.74 2.66 PS WSRAN2 3.03 2.74 2.67 2.79 2.69 2.74 2.76 2.86 Table 10. Computational results for the NEH algorithm: ARPD values III. Measure Indicator Sorting C D H V HIH HIL LOH LOL PMS SUM 2.84 2.76 2.70 2.85 2.80 2.80 2.71 2.71 PMS WSUM 2.84 3.30 2.80 2.91 2.70 2.70 2.78 2.85 PMS ABS 2.81 2.58 2.72 2.87 2.71 2.78 2.79 2.74 PMS WABS 2.89 2.64 2.72 2.79 2.66 2.75 2.83 2.73 PMS SRA 2.82 2.67 2.62 2.87 2.69 2.71 2.83 2.68 PMS WSRA 2.76 2.64 2.68 2.82 2.74 2.75 2.72 2.71 PMS SRS 2.73 2.66 2.70 2.79 2.68 2.72 2.76 2.66 PMS WSRS 2.81 2.70 2.64 2.78 2.74 2.72 2.79 2.74 PMS SRN 2.82 3.04 2.78 2.80 2.78 2.82 2.69 2.71 PMS WSRN 3.18 2.97 2.82 2.85 2.78 2.71 2.75 2.76 PMS SRAN 2.82 2.81 2.66 2.85 2.79 2.70 2.70 2.70 PMS WSRAN 3.04 2.75 2.83 2.78 2.67 2.67 2.76 2.74 PMS SRSN 2.77 2.93 2.75 2.83 2.71 2.77 2.77 2.69 PMS WSRSN 3.10 2.82 2.79 2.82 2.66 2.71 2.79 2.82 PMS SRAN2 2.81 2.89 2.72 2.82 2.80 2.82 2.75 2.71 Mathematics 2022,10, 3446 24 of 32 Table 10. Cont. Measure Indicator Sorting C D H V HIH HIL LOH LOL PMS WSRAN2 3.10 2.81 2.82 2.77 2.74 2.72 2.76 2.79 PmS SUM 2.84 2.78 2.72 2.81 2.78 2.73 2.81 2.71 PmS WSUM 2.80 3.20 2.80 2.90 2.64 2.67 2.72 2.76 PmS ABS 2.78 2.68 2.64 2.83 2.73 2.75 2.76 2.71 PmS WABS 2.85 2.63 2.68 2.78 2.66 2.68 2.78 2.66 PmS SRA 2.86 2.68 2.73 2.84 2.72 2.66 2.74 2.67 PmS WSRA 2.76 2.69 2.62 2.80 2.77 2.70 2.80 2.72 PmS SRS 2.79 2.65 2.70 2.81 2.71 2.70 2.79 2.72 PmS WSRS 2.87 2.66 2.62 2.81 2.73 2.71 2.82 2.77 PmS SRN 2.82 3.00 2.75 2.77 2.81 2.78 2.69 2.73 PmS WSRN 3.14 2.92 2.80 2.84 2.68 2.70 2.72 2.79 PmS SRAN 2.88 2.88 2.75 2.76 2.73 2.81 2.72 2.70 PmS WSRAN 2.98 2.74 2.81 2.81 2.71 2.67 2.80 2.82 PmS SRSN 2.78 2.92 2.76 2.84 2.79 2.83 2.77 2.68 PmS WSRSN 3.02 2.77 2.80 2.80 2.68 2.70 2.78 2.77 PmS SRAN2 2.80 2.85 2.77 2.77 2.76 2.78 2.69 2.75 PmS WSRAN2 3.11 2.72 2.82 2.77 2.66 2.68 2.83 2.77 Table 11. Computational results for the NEH algorithm: comparison with the most promising proposals in the literature. Parameter Value Our RA NM KK1KK2AD PRSKE NEH n 50 3.38 3.58 3.71 3.69 3.59 3.70 3.74 3.46 100 2.71 2.90 2.86 2.82 3.06 2.62 2.86 2.95 150 2.25 2.34 2.55 2.50 2.55 2.48 2.34 2.54 200 2.00 2.09 2.41 2.13 2.34 2.27 2.09 2.33 m10 2.96 3.10 3.28 3.27 3.25 3.15 3.18 3.20 20 2.21 2.36 2.49 2.30 2.52 2.39 2.34 2.44 r0 0.1 3.19 3.12 3.39 3.39 3.41 3.47 2.88 3.38 0.2 3.05 3.19 3.54 3.55 3.48 3.45 2.93 3.40 0.3 2.54 2.79 3.02 2.93 2.92 2.89 2.89 2.93 0.4 2.25 2.40 2.37 2.30 2.43 2.24 2.66 2.40 0.5 1.89 2.14 2.10 1.75 2.17 1.78 2.41 1.98 γ 25 2.07 2.18 2.27 2.12 2.28 2.06 2.50 2.08 50 2.54 2.82 2.90 2.81 2.90 2.73 2.80 2.86 100 3.15 3.17 3.48 3.42 3.47 3.51 2.97 3.51 Average 2.58 2.73 2.78 2.88 2.88 2.77 2.76 2.82 7.6. Computational Evaluation: IGV In this section, we analyse the performance of the proposed iterated greedy algorithm, IGV. To do so, we first compare it with both the MILP models and the complete enumeration procedures, on small-sized instances ( β1 ). Second, we compare our proposal with some of the most promising iterated greedy algorithms proposed in the related literature, on medium-large sized instances ( β2 ). To deal with this last issue, the following metaheuristics are re-implemented: • IGP: The iterated greedy algorithm proposed by [ 93 ] for the mixed no-idle permutation flow shop scheduling problem and makespan minimisation. • IGF: The iterated greedy algorithm proposed by [ 94 ] for the classical permutation flow shop scheduling problem with total tardiness minimisation. • IGL: The iterated greedy algorithm proposed by [ 95 ] for the distributed permutation flow shop with makespan minimisation. Mathematics 2022,10, 3446 25 of 32 • IGH: The iterated greedy algorithm proposed by [ 96 ] for the distributed assembly permutation flow shop scheduling problem with sequence-dependent setup times. • IGD: The iterated greedy algorithm proposed by [ 97 ] for the classical permutation flow shop scheduling problem with makespan minimisation. • IGV: The proposed iterated greedy algorithm. Regarding the parameters of the algorithm, we calibrate them based on [ 97 ]. More specifically, we vary parameter d between 3 and 8, setting T to 0.4. This parameter has not been found statistically significant according to the non-parametric Kruskal-Wallis ( p -value equal to 0.609). A similar result is found for parameter T when we vary it in { 0.2, 0.3, 0.4, 0.5, 0.6 } , setting d to 4. In this case, we obtain a p -value equal to 0.776. In Figure 4, we show the 95% confidence intervals obtained in both calibrations. Due to the fact that there are no statistically significant differences in d and T and the difference between the levels is very narrow, we use d= 4 and T= 0.4 as applied in the traditional algorithm by [84]. Figure 4. Calibration of dand T. To adapt these metaheuristics to the problem under consideration, we apply the best decoding procedure (i.e., PRS SI ), and initialise them with the best NEH algorithm (i.e., NEHV( { PMS,ABS,D } )). As stopping criterion, we use n×m×m× 10 / 1000, which depends of the size of the problem. The computational results obtained by IGV on the small-sized instances are shown in Table 4, while the comparison on medium-large sized instance is shown in Table 12. Regarding the former, the ARPD obtained by IGV is 1.68 compared to 0.01 found by the best Model 1, requiring only an average computational time of 1.08 seconds, instead of 61.61 seconds required by Model 1. Note that the proposed iterated greedy algorithm explores the space of solutions of PRS SI , whose optimum is 2.32, and incorporates a local search to try to capture the global optimum of the problem. Therefore, this intensive local search helps to efficiently escape from the local optimum. The computational results also demonstrate the efficiency of the proposed mechanism, as IGV clearly outperforms the best solutions obtained by any of the proposed complete enumeration procedures. Regarding the comparison with other metaheuristics from the related literature, we again compare the performance of the metaheuristics using the RPD indicator (see Equation (37) ), with besti as the best solution found in instance i by any of the iterated greedy algorithms. We also include in the comparison the best obtained NEH algorithm. The following conclusions can be derived from the results: • IGV obtains an ARPD of 0.49, clearly outperforming all other metaheuristics. This conclusion is also statistically confirmed by a non-parametric Mann-Whitney test comparing IGV and IGD (which obtains an ARPD of 0.61). A statistically significant difference has been found with the p-value equal to 0.003. • The constructive heuristic (NEHV( { PMS,ABS,D } )) improves its performance when n increases compared to the metaheuristics. This fact can be explained by the hardness of the problem. When n increases, the local search included in the metaheuristics requires much more computational time, and the number of global iterations decreases compared to small instances. Mathematics 2022,10, 3446 32 of 32 87. Zou, W.Q.; Pan, Q.K.; Tasgetiren, M. An effective iterated greedy algorithm for solving a multi-compartment AGV scheduling problem in a matrix manufacturing workshop. Appl. Soft Comput. 2021,99, 106945. [CrossRef] 88. Missaoui, A.; Ruiz, R. A parameter-Less iterated greedy method for the hybrid flowshop scheduling problem with setup times and due date windows. Eur. J. Oper. Res. 2022,303, 99–113. [CrossRef] 89. Lu, C.; Liu, Q.; Zhang, B.; Yin, L. A Pareto-based hybrid iterated greedy algorithm for energy-efficient scheduling of distributed hybrid flowshop. Expert Syst. Appl. 2022,204, 117555. [CrossRef] 90. Wang, Z.Y.; Pan, Q.K.; Gao, L.; Wang, Y.L. An effective two-stage iterated greedy algorithm to minimize total tardiness for the distributed flowshop group scheduling problem. Swarm Evol. Comput. 2022,74, 101143. [CrossRef] 91. Shao, Z.; Pi, D.; Shao, W. A novel discrete water wave optimization algorithm for blocking flow-shop scheduling problem with sequence-dependent setup times. Swarm Evol. Comput. 2018,40, 53–75. [CrossRef] 92. Vallada, E.; Ruiz, R. A genetic algorithm for the unrelated parallel machine scheduling problem with sequence dependent setup times. Eur. J. Oper. Res. 2011,211, 612–622. [CrossRef] 93. Pan, Q.K.; Ruiz, R. An effective iterated greedy algorithm for the mixed no-idle permutation flowshop scheduling problem. Omega 2014,44, 41–50. [CrossRef] 94. Fernandez-Viagas, V.; Valente, J.; Framinan, J. Iterated-greedy-based algorithms with beam search initialization for the permutation flowshop to minimise total tardiness. Expert Syst. Appl. 2018,94, 58–69. [CrossRef] 95. Lin, S.W.; Ying, K.C.; Huang, C.Y. Minimising makespan in distributed permutation flowshops using a modified iterated greedy algorithm. Int. J. Prod. Res. 2013,51, 5029–5038. [CrossRef] 96. Hatami, S.; Ruiz, R.; Andrés-Romano, C. Heuristics and metaheuristics for the distributed assembly permutation flowshop scheduling problem with sequence dependent setup times. Int. J. Prod. Econ. 2015,169, 76–88. [CrossRef] 97. Dubois-Lacoste, J.; Pagnozzi, F.; Stützle, T. An iterated greedy algorithm with optimization of partial solutions for the makespan permutation flowshop problem. Comput. Oper. Res. 2017,81, 160–166. [CrossRef]