A critical-path based iterated local search for the green permutation flowshop problem

Computers & Industrial Engineering 169 (2022) 108276 Available online 28 May 2022 0360-8352/© 2022 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/). A critical-path based iterated local search for the green permutation flowshop problem Victor Fernandez-Viagas a , * , Bruno de Athayde Prata b , Jose M. Framinan a , c a Industrial Management, School of Engineering, University of Seville, Camino de los Descubrimientos s/n, 41092 Seville, Spain b Department of Industrial Engineering, Federal University of Ceara, Ceara, Brazil c Laboratory of Engineering for Environmental Sustainability, University of Seville, Spain ARTICLE INFO Keywords: Green scheduling Bi-objective Sustainable manufacturing Flow shop Carbon emissions reduction ABSTRACT The permutation flowshop scheduling problem is a widely studied combinatorial optimization problem with several real-world applications. In this paper we address a green variant of the problem with controllable processing times and two objective functions: one related to the service level of the factory (makespan) and another one related to the total cost or the total energy/carbon consumption. For this problem we propose a novel Critical-Path based Iterated Local Search. This metaheuristic incorporates several theoretical results to accelerate the search of solutions in the intensification phase. The proposed algorithm has been compared on an extensive benchmark with the most promising algorithms in the literature. The computational results show the excellent performance of the proposal. 1. Introduction The flowshop scheduling problem is one of the most studied problems in the optimization literature (see e.g. the reviews by Fernandez- Viagas et al., 2017; Framinan et al., 2004; Reza and Saghafian, 2005; Ruiz and Maroto, 2005). This problem consists of obtaining the best sequence of jobs to be processed on a series of machines in order to optimize one or more criteria, where each job follows the same route of machines. Traditionally, this problem is addressed in the literature considering the same order of jobs in all machines and it is denoted as the Permutation Flowshop Scheduling Problem (PFSP). In the PFSP as well as in most scheduling literature, the processing times required for the jobs to be processed on the machines has been commonly considered as known and constant. However, in many real-world scenarios, these processing times may change if additional resources are assigned (Lu et al., 2014; Nowicki and Zdrzalka, 1990), a case that is denoted as Controllable Processing Times (CPT). Thereby, in a scheduling problem assuming CPT, the processing time of a job on a machine depends on one or several factors, such as the amount of resources assigned to the machine (which in this case usually models a manual operation), the skills of these resources, or the machine configuration applied. Reviews of manufacturing scheduling with CPT have been carried out by Nowicki and Zdrzalka (1990), Shabtay and Steiner (2007), Shioura et al. (2018), Fernandez-Viagas and Framinan (2015c). Given the global pressure to reduce carbon emissions (Wang et al., 2011; Zhu et al., 2014), in the last years, a growing attention is being paid in the literature (see e.g. Amiri and Behnamian, 2020; Cota et al., 2021; ¨ Oztop et al., 2020) to a special case of scheduling with CPT –denoted as green scheduling– where the CPT are aimed to reduce the energy consumed by the machines during their processing. In this paper, we focus on green scheduling in a PFSP context with two criteria: one related to the service level of the factory (makespan) and another one related to the total cost or the total energy/carbon consumption. This problem can be denoted as Fmprmu,pij(uij)#(Cmax,Cost)according to Graham et al. (1979), T’Kindt and Billaut (2006) and it is NP-hard since the single-objective decision problem (i.e. the problem without CPT, Fm|prmu|Cmax) is already NP-hard for more than two machines (Rinnooy Kan, 1976). In the literature, most researchers have addressed this problem by developing approximated procedures able to provide good solutions within a reasonable computational effort (see e.g. Mokhtari et al., 2011; ¨ Oztop et al., 2020). Despite these contributions, we are not aware that problem-specific properties have been studied to improve the effectiveness of the approximated procedures. In this regard, recent advances from related scheduling problems (see e.g. Benavides and Ritt, 2018; Fernandez-Viagas et al., 2020; Fernandez-Viagas and Framinan, 2015b; Fernandez-Viagas and Framinan, 2015a) show that much more * Corresponding author. E-mail address: [email protected] (V. Fernandez-Viagas). Contents lists available at ScienceDirect Computers & Industrial Engineering journal homepage: www.elsevier.com/locate/caie https://doi.org/10.1016/j.cie.2022.108276 Received 4 November 2021; Received in revised form 1 April 2022; Accepted 23 May 2022 Computers & Industrial Engineering 169 (2022) 108276 2 efficient solutions can be obtained by incorporating problem-specific properties in the algorithms rather than by using different metaheuristics or more complex local search procedures. Therefore, we believe that there is room for improving the state-of-the-art of approximate solution procedures for the problem and present a novel iterated local search which is based on several theoretical results. To present this contribution, the rest of the paper is organised as follows: the problem is formally described and its related literature reviewed in Section 2; next, in Section 3 we develop some theoretical results which are incorporated in the proposed iterated local search detailed in Section 4. The efficiency of the proposal is analysed by comparison against the most effective metaheuristics in the literature in Section 5. Finally, the conclusions are discussed in Section 6. 2. Problem description and background The problem under consideration can be defined as follows: a set 𝒥 of n jobs (𝒥:= {1,2,…,n}) has to be scheduled in a set ℳof m machines. Each machine i can perform the processing of the jobs using K different speed levels or speed modes. Let sik denote the speed of machine i if speed mode k (k∈ {1,…,K}) is used. Then, the processing time of job j on machine i when speed mode k is employed is pijk =pij sik , where pij is a known reference processing time of job j on machine i. Each speed mode is associated to a time unit processing cost, i.e. processing each job on machine i at speed mode k incurs a cost per time unit of pcik, therefore the cost of processing job j on machine j at speed k would be pcik⋅pijk. Furthermore, we assume that there is a time-unit cost ici due to having machine i idle (i.e. in stand-by mode). These two sources of costs would be discussed later in the context of energy consumption. Next, we briefly detail the rest of commonly assumed constraints in the problem under consideration: •Each machine can process at most a job at the same time. •Each job can be processed by at most one machine at the same time. •To be processed, all jobs has to followed the same route of machines. •Release times are set to 0. •Data are assumed to be deterministic. •No preemption is assumed. •Setup times are considered as insignificant or non-anticipatory and sequence-independent (and therefore they could be added to the processing times). •Transportation times can be considered either insignificant or constant for all jobs. Given the problem above, since each job j can be processed on each machine using a different speed, let us uij denote the speed mode at which job j is to be processed on machine i. Furthermore, since the permutation constraint is considered in our study, 𝒮:= {Π,𝒰} a feasible solution for the problem is univocally defined by giving a permutation Π= ( π 1,…, π n)where π j indicates the job to be processed in order j in all the machines, and by giving 𝒰 = (u11,…,uij,…,umn)the set of speed modes to be employed for each machine when processing each job. For the problem under study, we minimize a bi-objective function of the makespan (Cmax) and total cost (TC). The makespan is the maximum completion time of any job on the last machine (i.e. Cmax =maxjCm,j), where Cij the completion time of job j on machine i is obtained using Eq. (1). Ci, π j=max{Ci−1, π j,Ci, π j−1}+pi, π j si,ui, π j ,i={1,2,…,m},j={1,2,…,n} (1) and Ci, π 0=C0, π j=0. Regarding the total cost TC, it can be computed using Eq. (2). TC =∑ m i=1∑ n j=1 pci,uij ⋅pij,uij +∑ m i=1 ici⋅ITi(2) where ITi is the total idle time on machine i, which can be defined according to Eq. (3). ITi=Ci, π n−∑ n j=1 pij,uij (3) The costs involved in the problem under consideration (processing costs and idle costs) can be linked to sustainable aspects, particularly if we take into account that, in many scenarios, most of these costs are due to state-dependent machine energy consumption (i.e. the energy consumption is different depending on if the machine is in stand-by, or processing a job using a certain speed mode). In addition, if such energy is obtained from the combustion of fossil fuels, then its consumption should be reduced as it is related to the emission of carbon dioxide into the atmosphere. Therefore, several authors have addressed the problem considering an energy-efficient approach. It is the case of Mansouri et al. (2016), Zhang et al. (2017), Zhang et al. (2019), Ramezanian et al. (2019), ¨ Oztop et al. (2020), Amiri and Behnamian (2020), which consider the processing cost as: pcik =λk⋅δi 60 (4) where λk and δi represent a conversion factor of the energy required to work at speed k and the power consumed by machine i, respectively. Similarly, for the specific problem of carbon emissions minimisation, Foumani and Smith-Miles (2019) add a conversion factor to CO 2 (denoted as γ) and apply: pcik =λk⋅γ⋅δi 60 (5) Regarding the idle cost, all previous authors assume for the green problem: ici=φi⋅δi 60 (6) where φi represent a conversion factor of the energy incurred by machine i when is in stand-by (or idle). The problem under consideration is a multi-objective optimization problem. In this kind of problem, the goal is to obtain the Pareto-optimal solution set (denoted as ℘) which is composed by all non-dominated solution (see e.g. Fathollahi-Fard et al., 2018). Note that, given two solutions 𝒮1 and 𝒮2, we say that solution 𝒮1 dominates solution 𝒮2 (or equivalently 𝒮1≺ 𝒮2) when all objectives obtained by 𝒮1 are equal or better than the ones obtained by 𝒮2 and at least one is better (Minella et al., 2008; Zitzler et al., 2003). As the decision problem associated with at least one of our multi-objective functions corresponds to an NP-hard problem, the multi-objective optimization problem is NP-hard. Thus, heuristic algorithms are required to find good approximations of the Pareto front, taking into account proximity to the front as well as its spreading across the front. In the literature, several researchers have addressed either the problem under consideration or related variants. Mixed integer linear programming models have been proposed by Fang et al. (2013), Mansouri et al. (2016), Foumani and Smith-Miles (2019). Using a single objective function, Fang et al. (2013) investigate a related flowshop scheduling problem taking into consideration the energy consumption as a constraint (i.e. using a peak-power consumption constraint) instead of as objective function. For this case, a maximal power consumption constraint limits the utilization of the resources in a given production period. Two mixed-integer linear programming formulations are proposed for the makespan minimisation. Furthermore, two fast heuristics V. Fernandez-Viagas et al. Computers & Industrial Engineering 169 (2022) 108276 3 are introduced to find feasible schedules as well as to evaluate the tradeoff between makespan and peak-power consumption. Secondly, Mansouri et al. (2016) address the scheduling of a two-machine flowshop environment for the minimisation of makespan and energy consumption. A mixed-integer linear programming formulation is proposed to find the Pareto front to both objectives. A lower bound for both objectives are proposed, as well as a fast heuristic that finds good approximations for the Pareto front. Finally, Foumani and Smith-Miles (2019) address the flowshop scheduling problem with the bi-objective of makespan and total carbon emission minimisation. A mixed-integer linear programming formulation is proposed, and a weight aggregation approach converted the original problem into a single objective one. These authors employed random data and real-based data in computational experiments. Regarding approximated algorithms proposed in the literature, the most related approaches are the proposals by Mokhtari et al. (2011), ¨ Oztop et al. (2020). The former considers a flowshop with CPT where a trade-off between makespan and the required amount of resources is sought to minimize. The authors decompose the problem into two subproblems (a sequencing problem and a resource allocation problem) and propose a solution procedure combining a discrete differential evolution algorithm (DDE) with a variable neighbourhood search (VNS). Regarding the latter, ¨ Oztop et al. (2020) address an energy-efficient flowshop scheduling problem for the minimisation of total flow time and total energy consumption. These authors present a bi-objective mixed-integer programming model formulation with a speed-scaling framework. An improved version of the well-known NEH algorithm (Nawaz et al., 1983) is proposed as initial solution. Furthermore, two variants of the iterated greedy metaheuristic and variable blockinsertion algorithm are proposed. The iterated greedy algorithms presented better results in the computational experiments carried-out with randomly generated test instances. Finally, note that several interesting researches have also been proposed for related scheduling problem considering either no-wait constraints (see e.g. Gao et al., 2018; Gomes et al., 2020), either different setup times approaches (see e.g. Ramezanian et al., 2019), or hybrid flowshop layouts (see e.g. Behnamian and Ghomi, 2011, 2015, 2019, 2002). Summarising, the present literature review highlights that several authors have addressed some related scheduling problems and therefore, although there is no optimization procedure for the problem under consideration, we found some metaheuristics dealing with related ones. These metaheuristics use novel mechanisms to conduct the local searches, yet they do not use problem-specific properties in their phases, a fact that can be exploited to substantially improve the performance of these methods. To cover this gap, in the following section we introduce some theoretical results, which will be incorporated in our proposal in Section 4. 3. Theoretical results In this section, we detail some theoretical results necessary to develop our proposal. Firstly, we apply some common definitions of the scheduling literature to the problem under consideration: critical path, forward and backward decoding procedures, and slack time. Using them, we can prove four theorems that are applied in the proposed iterated local search-based metaheuristic (see Section 4) to improve its efficiency. Let us start by defining two different decoding procedures to obtain a schedule from a solution of the problem 𝒮 = (Π,𝒰). Firstly the Forward Decoding obtains a semi-active schedule (see definition in Pinedo (2012)) inserting jobs from the left to the right of Π, while the Backward Decoding considers the reverse sequence executed in the reverse order of machines (see e.g. Ribas et al., 2010; Ribas et al., 2013). A formal definition of both procedures, applied to the problem under consideration, is detailed next: Definition 3.1.(Forward Decoding) Let 𝒮 = (Π,𝒰), with Π:= ( π 1,…, π n), be a solution for the ℐinstance of the Fmprmu,pij(uij)#(Cmax,Cost) problem. Then, the forward decoding constructs a solution by processing each job π j on each machine i (with j∈ {1,…,n}and i∈ {1,…,m}) according to Eq. (1), with CFD i, π 0=C0, π j=0 and CFD ij =Cij (∀j∈ {1,…,n},∀i∈ {1,…,m}). Definition 3.2.(Backward Decoding) Let 𝒮 = (Π,𝒰), with Π:= ( π 1, …, π n)and 𝒰 = (u11,…,uij,…,umn), be a solution for the ℐinstance of the Fmprmu,pij(uij)#(Cmax,Cost)problem. Then, the backward decoding constructs a solution by processing each job π j on each machine i following the reverse sequence of jobs and machines (j∈ {n,n−1,…,1} and i∈ {m,m−1,…,1}), i.e. starting by processing job π n on machine m. The completion time of operation Oi π j (denoted by CBD i, π j) is defined by Eq. (7). Note that Cij is the reverse completion time (i.e. starting from the final operation) defined in Eq. (8), where Ci, π n+1=Cm+1, π j=0 (∀j∈ {n, n−1,…,1},∀i∈ {m,m−1,…,1}). CBD i, π j=Cmax −Ci, π j+pi, π j,ui π j,∀j∈{n,n−1,…,1},∀i∈{m,m−1,…,1}(7) Ci, π j=max{Ci+1, π j,Ci, π j+1}+pi, π j,ui π j,∀j∈{n,n−1,…,1},∀i ∈{m,m−1,…,1}(8) An example of forward and backward codifications is shown in Fig. 1 for five jobs and four machines. At the top of the figure, a schedule is constructed from the sequence (1,2,3,4,5) using the forward decodification, while alternatively at the bottom of the figure its backward decoding is shown. We can observe that the makespan obtained by a solution 𝒮 = (Π,𝒰) using forward or backward decoding procedure is the same, which is a known property of the PFSP, denoted as Reversibility of the PFSP (see e.g. Ribas et al., 2010 for more details). Note that, when the speed modes (𝒰) are fixed, the problem under consideration, Fmprmu,pij(uij)#(Cmax, Cost), is obviously equivalent to the well-known Fm|prmu|Cmax problem, and in consequence this property applies for our problem. The reversibility in the problem is based on the existence of the Critical Path of a Fig. 1. Example for theoretical results. V. Fernandez-Viagas et al. Computers & Industrial Engineering 169 (2022) 108276 4 solution 𝒮, denoted by ℛ(𝒮), whose definition in the traditional PFSP can be found in Fernandez-Viagas et al. (2020). Nevertheless, for the sake of completeness, we give the definition for the problem under consideration: Definition 3.3.(Critical Path) Let 𝒮 = (Π,𝒰), with Π:= ( π 1,…, π n), be a solution of ℐan instance of the Fmprmu,pij(uij)(Cmax,Cost)problem. Then, considering semi-active schedules, we define the critical path of 𝒮,𝒞(𝒮), as a path composed of operations whose sum of processing times is equal to the makespan and it can go from the first operation O1, π 1 to the last one (Om, π n) moving always in direction of machines (from i to i+1) or jobs (from π j to π j+1). In Fig. 1, using the same example, we can observe the critical path for the sequence (1,2,3,4,5) in grey at the top of the figure, formed by Operations (O11,O21,O22,O23,O24,O34,O35,O45). These definitions are the basis of the accelerations proposed by Taillard (1990) (denoted as Taillard’s accelerations) for the Fm|prmu|Cmax problem (see Corollay 3.3 in Fernandez-Viagas et al. (2020) for more details), which explains an efficient mechanism to perform insertion-based methods in the Fm|prmu|Cmax problem. These accelerations can be applied to the problem under consideration when the modes 𝒰to execute the operations are fixed (Theorem 3.1). Theorem 3.1.Let 𝒮 = (Π,𝒰), with Π:= ( π 1,…, π n−1), be a solution of ℐ an instance of the Fmprmu,pij(uij)(Cmax,Cost)problem. Furthermore, let Ci, π j denote the reverse completion times of 𝒮and C σ i,j the completion time of job σ on machine i if this job σ is inserted in position j of Π using the set of modes ui σ (∀i∈ {1,…,m}). Then, the best makespan that can be obtained by inserting σ in any position of Π, maintaining the same modes 𝒰can be computed as Cmax =min j=1,…,n{max i=1,…,m{C σ ij +Ci, π j}} (9) The proof of this theorem is trivial in view of Corollary 3.3 by Fernandez-Viagas et al. (2020). Next, let us define the final term Slack Time (stij) as: Definition 3.4.(Slack Time) Let 𝒮 = (Π,𝒰), with Π:= ( π 1,…, π n), be a solution of ℐan instance of the Fmprmu,pij(uij)#(Cmax,Cost)problem. Then, we define sti, π j the operation slack time Oi, π j as the difference between the completion time of π j obtained using forward and backward procedures, i.e. CFD i, π j, and CBD i, π j respectively. More specifically: sti, π j=CBD i, π j−CFD i, π j(10) In Fig. 1, the slack times of operations O13 and O42 in the previous example are shown. Using this definition, we can present Theorems 3.2 and 3.4 that discuss which operations can be compressed or delayed in a schedule without increasing the makespan. Note that these theorems give conditions to reduce the speed modes of certain operations (and, therefore, reducing the total costs in TC) without worsening the makespan objective and will be used in the solution procedure presented in Section 4. Theorem 3.2.Let ℱbe the forward schedule obtained from solution 𝒮 = (Π,𝒰), with Π:= ( π 1,…, π n), using instance ℐ. Then, without altering solution 𝒮, the operation Oi, π j can be delayed for a maximum of sti, π j time units without increasing the makespan. Proof. According to the reversibility property, the makespan obtained using the backward decoding is the same as using the forward schedule. Then, by Definition 3.4, the backward schedule represents a feasible schedule where operation Oi, π j finishes sti, π j time units later than in the forward schedule, but with the same final makespan. Therefore, operation Oi, π j can be delayed for a maximum of sti, π j time units without increasing the makespan. Now, let us analyze the possibility to delay operation Oi, π j more time units than sti, π j without increasing the makespan. According to Definition 3.2, the backward decoding generates a semi-active schedule from the right to the left (see Eq. 7). Hence, any delay of an operation in the backward schedule (without altering either the orders of jobs in the machines or the modes applied) must necessary delay that time the final makespan. As a consequence, an operation Oi, π j can be delayed for at maximum sti, π j time units without increasing the makespan. square Theorem 3.3.Let ℱbe the forward schedule obtained from solution 𝒮 = (Π,𝒰), with Π:= ( π 1,…, π n), using instance ℐwith critical path ℛ(𝒮). All operations in ℛ(𝒮) do not have slack time, i.e. stOij =0,∀Oij ∈ ℛ. Proof. By the definition of the critical path, there can be no idle or waiting times between operations belonging to the critical path, i.e. they cannot be contracted or expanded without altering solution 𝒮. As a consequence, the completion times of these operations must remain unaltered if either forward or backward decoding is applied and, consequently their slack time is 0. square Corollary 3.1. Let ℱbe a forward schedule for instance ℐobtained from solution 𝒮and ℛ(𝒮) be the corresponding critical path. Then, a delay of δ time units in any operation in ℛ(𝒮) increases the makespan by exactly δ units. The proof of the corollary is trivial according to Theorem 3.3 and Definition 3.3. Theorem 3.4.Let ℱbe the forward schedule obtained from solution 𝒮 = (Π,𝒰), with Π:= ( π 1,…, π n), using instance ℐwith critical path ℛ(𝒮). Then, applying ℱ, any reduction in the completion time of an operation Oi, π j not belonging the critical path cannot reduce the final makespan. The proof of this theorem is trivial by Definition 3.3 and using the same reasoning as in 3.3. Theorem 3.5.Let ℱbe the forward schedule obtained from solution 𝒮 = (Π,𝒰), with Π:= ( π 1,…, π n), using instance ℐwith critical path ℛ(𝒮). Then, any reduction in the completion time of an operation Oi, π j belonging the critical path could potentially reduce the final makespan Cmax. Proof. By Definition 3.3, if some operation in the critical path is reduced, there are two possibilities: 1. Some idle time is inserted between operations in the critical path. In this case, there exists a new critical path with makespan less or equal than Cmax; 2. Otherwise, the critical path remains the same, and then, the makespan is reduced such amount of time. square The previous theorems and corollaries establish a guideline to be followed if some operation is advanced or delayed in the schedule. Note that in the problem under consideration, this fact is very relevant as the processing times of the operations (and consequently their completion times) depend on the modes employed to execute them. As a summary of the previous theoretical results, we can conclude that: •Given a solution 𝒮 = (Π,𝒰), Taillard’s accelerations represent an efficient method for reducing the complexity of computing the makespan if job σ is inserted in the best position of Π without changing the set of modes 𝒰(Theorem 3.1). •Regarding the operations in the critical path: - A delay in the completion time in any of these operations causes an analogous increase in the final makespan. V. Fernandez-Viagas et al. Computers & Industrial Engineering 169 (2022) 108276 5 - A reduction of the completion time in any of these operations (due, e.g., to reducing its processing time) can eventually reduce the final makespan. •For any operation Oij outside the critical path: - A delay δ (with δ<stij) in the completion time of such an operation does not increase the final makespan. - Reducing the processing time of such an operation does not reduce the final makespan. These theoretical results are embedded in a local search-based metaheuristic presented in the next section. 4. Critical-Path based Iterated Local Search, CP_ILS In this section, we detail the proposed critical-path based iterated local search, denoted as CP_ILS, to obtain a Pareto set of solutions, denoted as ℘. Starting from an initial solution, the traditional iterated local search iteratively perturbs a solution and searches in its neighbourhood for the best local solution. Despite being a basic procedure, this type of metaheuristic has found excellent results in solving permutation-based scheduling problems when a very efficient local search is applied, being state-of-the-art in several related scheduling problems (see e.g. Dong et al., 2013; M’Hallah, 2014; Subramanian et al., 2014). In each iteration of our proposal, after applying a simple perturbation mechanism (see Section 4.2 for more details), we carry out an intensive search using two different types of local search methods (types 𝒯S and 𝒯U) that do not modify either the sequence of jobs or the speed modes. More specifically, the specific search of the proposal is composed of the following four different local search methods (one of type 𝒯S and three of type 𝒯U) as follows: •Type 𝒯S: An improvement in the current solution 𝒮:= (Π,𝒰) is sought by altering thesequence Π while maintaining 𝒰. - The first proposed local search method performs insertion-based movements in Π without changing 𝒰. Under these conditions, the requirements of Theorem 3.1 are fulfilled. In order to take advantage of it, we search for the best solution after inserting each job in each position. This local search is denoted LS_S. See Section 4.3 for more details. •Type 𝒯U: An improvement in the current solution 𝒮:= (Π,𝒰) is sought by altering the modes configuration 𝒰while maintaining Π. For this type, the following local search methods are proposed: - LS_M: In this local search, we look for new solutions by changing the current speed mode of a job and assuming that this mode is used for all machines. The procedure is explained in detail in Section 4.4. Fig. 2. Critical-path based iterated local search. V. Fernandez-Viagas et al. Computers & Industrial Engineering 169 (2022) 108276 6 - LS_NCP: In this critical-path based local search, we use Theorem 3.2 and decrease the speed mode of the operations outside the critical path. The procedure is explained in detail in Section 4.5. - LS_CP: The second critical-path based local search applies Theorem 3.5 and increases the speed mode of the operations in the critical path. The procedure is explained in detail in Section 4.6. A global pseudo-code of CP_ILS is shown in Fig. 2, and in the next sections we describe its main parts. 4.1. Initial solution To generate an initial solution, we first generate the mode of each operation Oij. For simplicity, we assign a random mode for every job, and this mode is used for all machines, i.e. uij =uj,∀i,j. Once the set of nodes has been assigned, a simple NEH algorithm is applied to construct the initial sequence. This well-known algorithm (originally proposed by Nawaz et al. (1983) for the Fm|prmu|Cmax problem) obtains a sequence by sorting the jobs in Longest Processing Times order (denoted by ΠLPT = ( π LPT 1,…, π LPT n)). The initial sequence, denoted by ΠNEH, is then constructed by iteratively inserting job π LPT j (j∈ {2,…,n}) in ΠNEH, in the position that yields the minimal (partial) makespan. In the first iteration, ΠNEH is composed of a single job, i.e. ΠNEH = ( π LPT 1). The pseudo-code to construct the initial solution is shown in Fig. 3. Note that, as stated in Theorem 3.1, Taillard’s accelerations can be applied to compute such a minimal makespan. Hence, a calculation of the total cost is not performed in this phase. After a sequence is obtained using the NEH algorithm, the local search method LS_S is applied and the so-obtained solution is appended to ℘the Pareto set of solutions found by the algorithm, initially empty. Then, the local search method LS_M is applied to improve the solution. These two methods are discussed in detail in Sections 4.3 and 4.4, respectively. Once the two methods are applied, the algorithm performs a series of iterations until the stopping criterion (i.e. running time) is met. The first step in the iterations consists of the so-called perturbation phase, which is described in the next section. 4.2. Perturbation phase In the perturbation phase, a different element from the space of solutions is intensively explored using advanced local search mechanisms. To ensure the algorithm performance, this new element should satisfy the following issues: firstly, it should have some similarity to a previous element to avoid a completely new restart in each iteration, which would make it difficult to find a high-quality local optimum; secondly, the new element should be different enough from the previous one so that the new local optimum could potentially be different. In this regard, Fernandez-Viagas et al. (2018) explored this trade-off by comparing eight different perturbation mechanisms. Among them, the best result was obtained by a simple perturbation composed of several random adjacent swaps. This perturbation is also considered in the proposed metaheuristic. More specifically, in this phase, a job is selected randomly and interchanged with the next job of the sequence. This procedure is repeated d times, where d is a parameter of the algorithm to be calibrated. 4.3. Local search: LS_S In this section, we detail the local search method denoted by LS_S. This search starts with a sequence of jobs (Π) and a speed mode for each job that is maintained for all machines (uij =uj,∀i). Then, a local search is carried out by changing the sequence of jobs, without altering the set of speed modes (type 𝒯S of local search methods). Thereby, each job of the sequence is removed from the sequence and reinserted in the position that minimises the makespan. Note that, as the modes configuration is the same for all sequences tested, Theorem 3.1 can be applied to accelerate the computation. This procedure is repeated until no improvement is found. A detailed pseudo-code is shown in Fig. 4. 4.4. Local search: LS_M The LS M method is our first proposed local search of type 𝒯U. This simple search receives a sequence of jobs, Π, and a mode configuration Fig. 3. NEH algorithm. Fig. 4. LS_S function. V. Fernandez-Viagas et al. Computers & Industrial Engineering 169 (2022) 108276 7 for each job j,uj (with j∈ {1,…,n}), which is used to obtain the speed to be set in every machine. Iteratively, the procedure selects a job j (with j∈ {1,…,n}) and increases its mode one unit. If the resulting mode is greater than K, it is set to the first mode. Then, the procedure firstly checks if the new solution belongs to the set of non-dominated solutions and, second, it replaces the current solution if it has been dominated. A detailed pseudo-code is shown in Fig. 5. Note that this local search is not very intensive. The idea here is to have a good initial configuration of modes to start the next local method. Hence, by considering both the same mode in all machines and a unique change in each mode, we can accelerate the search to set up a good initial configuration of modes. An excess of intensification in this phase could be counterproductive, since in the following phases (i.e. LS_NCP and LS_CP) the mode in each machine is going to be altered. 4.5. Critical-path based local search: LS_NCP The second local search of type 𝒯U is based on the concept of critical path (see Definition 3.3) and and more specifically, on Theorem 3.2. Note that this theorem states that an operation outside the critical path can be delayed by stij time units without increasing the makespan of the solution. Hence, the proposed method first calculates the slack time stij of every operation Oij by performing the forward and backward decoding of the initial solution (see Eq. (10)). Then, these slack times are tried to be reduced by increasing the processing times of the corresponding operations. To do so, in each iteration, the speed of job j on machine is tried to be increased. If the increase in the processing time is lower than the corresponding slack time, the new speed replaces the last one and the values of the objective functions are computed. Finally, if the so-obtained solution is non-dominated, it is introduced in ℘, see the pseudo-code in Fig. 6 for more details. Note that the slack times (stij) are not recomputed in each iteration in order to reduce the computational effort of the search. 4.6. Critical-path based local search: LS_CP The last local search of type 𝒯U is based on the operations in the critical path (see Definition 3.3). In this case, according to Theorem 3.5 only these operations could potentially decrease the makespan. Therefore, as in the previous local search method, the slack time of each operation is first calculated using Eq. (10). Once these are computed, the speed of job j on machine i is iteratively decreased and the corresponding values of the objective functions are computed. If the soobtained solution is non-dominated, it is introduced in ℘, and is set as the current solution. Then, the whole procedure is repeated while there are some improvements, as illustrated in Fig. 7. As in the previous case, the slack time (stij) is not recomputed once a new solution is found to reduce the computational effort of the method. 4.7. Simulated annealing procedure Finally, in the last step of the metaheuristic, a simple simulated annealing procedure is applied to decide the current solution (𝒮 = {Π, uj}) 1 . To deal with that, a simulated annealing phase with a constant Fig. 5. LS_M function. Fig. 6. LS_NCP function. 1 Note that the current solution includes only one mode for each job j, and is consequently assuming that the same mode is applied to all machines V. Fernandez-Viagas et al. Computers & Industrial Engineering 169 (2022) 108276 8 temperature is applied (see e.g. Fernandez-Viagas and Framinan, 2015b; Pan et al., 2008; Ruiz and Stützle, 2007). In this phase, a new solution 𝒮′, is kept in the next iteration if the makespan in the last iteration has been outperformed: random⩽exp{ − (Cmax(𝒮′) − Cmax(𝒮))/Temperature}(11) where Temperature can be computed using parameter T as: Temperature =T⋅∑ ∀i∑ ∀j pij n⋅m⋅10 (12) 5. Computational evaluation In this section, we check the effectiveness of the proposed metaheuristic by comparing it against the most relevant metaheuristics proposed in the related literature for the same problem. This comparison is performed using three sets of medium-large instances. Both sets are detailed in Section 5.1. The indicators chosen to compare the algorithms are shown in Section 5.2, while the algorithms implemented are detailed in Section 5.3. Finally, the computational results are presented in Section 5.4. 5.1. Benchmark According to Section 2, the following parameters completely define an instance of the problem: n,m,S,pijk,sik,pcik, and ici. In this paper, these parameters have been chosen to cover a wide scope of the problem under consideration as well as its special case of energy consumption minimisation (green approach). Next, we explain the generation of the parameters: •Number of machines: Analogously, we consider equidistant values (m∈ {5,10,15,20}) based on Ramezanian et al. (2019). •Processing times: They are generated from a uniform distribution between 1 and 99 (see e.g. Foumani and Smith-Miles, 2019; Mansouri et al., 2016; ¨ Oztop et al., 2020). •Number of modes: Three modes (K=3) are generated to represent the different types of speeds of each machine, following the most common approach used in the literature (see e.g. Amiri and Behnamian, 2020; Foumani and Smith-Miles, 2019; Mansouri et al., 2016; ¨ Oztop et al., 2020; Ramezanian et al., 2019; Zhang et al., 2019) •Type of benchmark: For these parameters, three types of instances are generated (ℬ1,ℬ2, and ℬ3) representing different cases of the problem under consideration. The rest of the parameters (i.e. n,sik, pcik, and ipi) are defined according to the following: –Set of instances ℬ1. These instances are designed to test the algorithms on a generic random benchmark. For the number of jobs, we select equidistant values using bounds similar to Zhang et al. (2019), Mansouri et al. (2016), i.e., we select n∈ {40,80,120, 160}. Then, 30 iterations are generated for each combination of n and m, resulting in a total of 480 instances. The speed parameter (sik) is generated using a uniform distribution U[0.8,1.2](the same distribution as ¨ Oztop et al. (2020) but with bounds taken from Mansouri et al. (2016)). Three values are generated for each machine i and the highest value is assigned to the mode k=1, the medium value to the mode k=2, and the lowest value to the mode k=3. Regarding the costs, three processing costs (pcik) are generated using a uniform distribution U[0.5,2.0]for each machine i. Analogously, the highest, medium and lowest values are assigned to modes 1, 2 and 3, respectively. Finally, the idle costs are generated using U[0.025,0.100]. Fig. 7. LS_CP function. V. Fernandez-Viagas et al. Computers & Industrial Engineering 169 (2022) 108276 9 – Set of instances ℬ2. This set is designed to tackle the green special case of the problem and is composed of 480 instances varying the previous combinations of n and m (n∈ {40,80,120,160}and m∈ {5,10,15,20}). To generate this set, we use the same approach as ¨ Oztop et al. (2020), Amiri and Behnamian (2020). More specifically, we use λk∈ {1.5,1.0,0.6},φ=0.05, and δi=60. Replacing these values in Eqs. 4 and 6, we have ici=0.05 and pcik =pck= {1.5,1.0,0.6}for modes 1, 2, and 3, respectively. Regarding the speed parameter, sik, it is generated for three modes, regardless of the machine, as sik =sk∈ {1.2,1.0,0.8}. – Set of instances ℬ3. This set is designed to test the algorithms in big-size instances (n∈ {300,400,500}). The rest of the parameters are generated using the same distributions as ℬ1. This benchmark is composed of 30 instances for each combination of n and m, resulting in a total of 360 instances. 5.2. Indicators The literature is abundant in the use of indicators to analyse the performance of multi-objective approaches (for the interested reader, we refer to Zitzler et al. (2003)). In our case, in order to have a fair comparison with the algorithms implemented, we applied indicators similar to those of ¨ Oztop et al. (2020). Before introducing such indicators to evaluate the algorithms, let us assume that ℘A is the Pareto set of solutions found by Algorithm A. In addition, let ℘* be the reference Pareto set composed by the best non-dominated solutions found in an instance, formed by the solutions reached by any of the implemented algorithms. By doing so, based on the literature, we can measure the efficiency of the algorithms by using the following indicators to compare the Pareto set of each different algorithm. Note that these indicators try to measure the density, spread and quality of each set of solutions. These indicators are as follows: •Number of non-dominated solutions (see e.g. Fathollahi-Fard et al., 2018; Kacem and Dammak, 2019; Pan et al., 2009), NNDS. This indicator represents the number of global non-dominated solutions found by Algorithm A. •IR: Ratio of non-dominated solutions (see e.g. Pan et al., 2009). It represents the ratio of global non-dominated solutions found for Algorithm A with respect to the total solutions in ℘*, denoted as IR(see Eq. 13). IR(A)=NNDS(A) |℘*|(13) •IA: Average distance to ℘* (Pan et al., 2009; Zhang and Li, 2007). It is the mean of the normalized Euclidean distance between each element in ℘* and the closest element in the Pareto set of algorithm A, i.e. ℘A. This normalized distance for each element y∈ ℘*, denoted dy(℘A), can be defined by Eq. (14), where Cmaxmax (Cmaxmin) is the maximum (minimum) makespan value found in ℘* and similarly, Costmax (Costmin) is the maximum (minimum) cost in ℘*. Once dy(℘A) is known, IA is its mean value (IA=d(℘A)) according to Eq. (15). A lower value of this indicator is desirable. dy(℘A)=min x∈℘A{(Cmax(x) − Cmax(y) Cmaxmax −Cmaxmin )2 +(Cost(x) − Cost(y) Costmax −Costmin )2} (14) IA=d⎛ ⎜ ⎝℘A⎞ ⎟ ⎠=∑ y∈℘* dy(℘A) |℘*|(15) Table 1 NNDS of each metaheuristic grouped by n and m in benchmark ℬ1. In bold, we represent the best found value for each instance and stopping criterion. n ×m t =60 t=90 t=120 HDDE IG IGALL VBIH CP_ILS Total HDDE IG IGALL VBIH CP_ILS Total HDDE IG IGALL VBIH CP_ILS Total 40 ×5 13.58 2.17 2.46 0.69 179.75 198.65 17.28 2.45 2.41 0.61 177.07 199.81 15.77 2.59 2.80 0.68 178.63 200.46 40 ×10 10.08 1.73 1.95 0.59 255.25 269.60 8.24 1.95 2.50 0.54 250.18 263.41 11.10 2.03 2.65 0.56 256.24 272.58 40 ×15 1.52 1.77 1.95 0.45 403.08 408.77 4.01 1.99 2.10 0.45 400.40 408.95 4.91 2.43 2.41 0.48 398.79 409.01 40 ×20 1.13 1.93 1.79 0.55 476.93 482.32 1.82 1.97 2.22 0.51 480.65 487.18 3.26 2.23 2.46 0.46 473.49 481.90 80 ×5 24.79 1.82 1.93 0.87 275.81 305.22 30.14 2.13 2.22 0.84 288.14 323.47 35.79 2.36 2.45 0.74 280.70 322.03 80 ×10 2.16 1.27 1.32 0.45 517.72 522.91 2.32 1.35 1.78 0.51 520.51 526.47 5.54 1.65 1.88 0.47 522.13 531.67 80 ×15 0.19 1.12 1.36 0.51 696.96 700.15 0.43 1.57 1.54 0.51 691.22 695.26 1.04 1.61 1.66 0.48 703.53 708.33 80 ×20 0.41 1.25 1.30 0.53 949.59 953.07 0.46 1.29 1.45 0.49 948.30 951.99 0.41 1.35 1.75 0.45 958.71 962.66 120 ×5 8.12 1.65 1.77 0.86 454.79 467.18 19.18 2.06 1.88 0.89 447.92 471.93 23.69 2.09 2.03 0.97 444.23 473.00 120 ×10 1.37 1.56 1.57 0.75 741.49 746.74 1.00 2.01 2.19 0.66 746.77 752.62 2.35 2.19 1.95 0.73 749.01 756.23 120 ×15 0.04 0.79 1.00 0.43 1075.63 1077.90 0.48 1.11 1.45 0.44 1070.26 1073.74 0.31 1.13 1.16 0.47 1076.43 1079.51 120 ×20 0.00 0.93 1.10 0.51 1387.41 1389.95 0.00 1.05 1.18 0.44 1390.77 1393.45 0.10 1.45 1.22 0.48 1391.91 1395.17 160 ×5 0.00 2.25 2.27 1.19 675.84 681.55 8.47 2.17 2.09 1.14 671.61 685.49 18.21 2.17 2.18 1.17 651.23 674.96 160 ×10 0.00 1.39 1.59 0.70 1040.63 1044.32 1.01 1.43 1.49 0.71 1031.02 1035.66 0.75 1.57 1.49 0.74 1041.47 1046.02 160 ×15 0.40 0.96 1.24 0.59 1288.07 1291.27 1.16 1.11 1.19 0.49 1307.73 1311.68 1.37 0.97 1.24 0.44 1306.15 1310.18 160 ×20 0.25 1.23 1.27 0.63 1626.67 1630.04 0.34 1.13 1.25 0.55 1612.29 1615.55 0.20 0.95 1.29 0.63 1640.31 1643.38 Average 4.00 1.49 1.62 0.64 752.85 760.60 6.02 1.67 1.81 0.61 752.18 762.29 7.80 1.80 1.91 0.62 754.56 766.69 V. Fernandez-Viagas et al.