An iterated greedy algorithm for solving the total tardiness parallel blocking flow shop scheduling problem

Accepted Manuscript

An Iterated Greedy Algorithm for Solving the Total Tardiness Parallel Blocking Flow Shop Scheduling Problem Imma Ribas , Ramon Companys , Xavier Tort-Martorell PII: DOI: Reference:

S0957-4174(18)30807-8 https://doi.org/10.1016/j.eswa.2018.12.039 ESWA 12380

To appear in:

Expert Systems With Applications

Received date: Revised date: Accepted date:

21 April 2018 30 October 2018 20 December 2018

Please cite this article as: Imma Ribas , Ramon Companys , Xavier Tort-Martorell , An Iterated Greedy Algorithm for Solving the Total Tardiness Parallel Blocking Flow Shop Scheduling Problem, Expert Systems With Applications (2018), doi: https://doi.org/10.1016/j.eswa.2018.12.039

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Highlights An iterated greedy algorithm for the total tardiness parallel blocking flow shop problem Constructive heuristics which considers the characteristics of the problem Design of experiments to calibrate and select the best variant Identification of the main characteristics and parts of the algorithm to increase its performance

AC

CE

PT

ED

M

AN US

CR IP T

   

1

ACCEPTED MANUSCRIPT

An Iterated Greedy Algorithm for Solving the Total Tardiness Parallel Blocking Flow Shop Scheduling Problem 1

Imma Ribas ,a, Ramon Companysb, Xavier Tort-Martorellc a, Departament d’Organització d’Empreses, DOE – ETSEIB - Universitat Politècnica de Catalunya. BarcelonaTech, Avda. Diagonal, 647, 7th Floor, 08028 Barcelona, Spain. E-mail: [email protected] b, CDE - EPSEB - Universitat Politècnica de Catalunya. BarcelonaTech, Gregorio Marañón 44-50, 3rd Floor, 08028 Barcelona, Spain. E-mail: [email protected]

CR IP T

c, Departament de Estadística e investigación Operativa- ETSEIB - Universitat Politècnica de Catalunya. BarcelonaTech, Avda. Diagonal, 647, 6th Floor, 08028 Barcelona, Spain. E-mail: [email protected]

Abstract

This paper proposes an iterated greedy algorithm for scheduling jobs in F parallel flow shops (lines), each consisting of a series of m machines without storage capacity between machines.

AN US

This constraint can provoke the blockage of machines if a job has finished its operation and the next machine is not available. The criterion considered is the minimization of the sum of tardiness of all the jobs to schedule, i.e., minimization of the total tardiness of jobs. Notice that the proposed method is also valid for solving the Distributed Permutation Blocking Flow Shop Scheduling Problem (DBFSP), which allows modelling the scheduling process in companies

M

with more than one factory when each factory has an identical flow shop configuration. Firstly, several constructive procedures have been implemented and tested to provide an efficient solution in terms of quality and CPU time. This initial solution is later improved upon with an

ED

iterated greedy algorithm that includes a variable neighbourhood search for interchanging or reassigning jobs from the critical line to other lines. Next, two strategies have been tested for

PT

selecting the critical line; the one with a higher total tardiness of jobs and the one with a job that has the highest tardiness. The experimental design chooses the best combination of initial

CE

solution and critical line selection. Finally, we compare the performance of the proposed algorithm against other benchmark algorithms proposed for the DPFSP, which have been adapted to the problem being considered here since, to the best of our knowledge, this is the first

AC

attempt to solve either the Parallel Blocking Flow Shop problem or the Distributed Blocking Flow Shop problem with the goal of minimizing total tardiness. This comparison has allowed us to confirm the good performance of the proposed method. Keywords: Parallel flow shop, Distribution flow shop, Blocking, Scheduling, Total Tardiness

1

Corresponding author. E-mail address: [email protected] Phone number +34 93 401 65 87

2

ACCEPTED MANUSCRIPT 1. Introduction Even though the Parallel Flow Shop Scheduling Problem has not receive proper attention in the scheduling literature, one can found it in different manufacturing plants. For instance, He, Kusiak, and Artiba (1996) studied the problem in a glass industry with F parallel flow shops having two

machines in each. In this productive configuration, deciding which jobs has to produce in each flow shop and in which sequence if very relevant. Not to mention when the company wants to

CR IP T

offer a high service levels, i.e., to serve their products according to the promised delivery dates. It is worth noticing that the parallel flow shop scheduling problem (PFSP) and the Distributed Permutation Flow Shop Scheduling Problem (DPFSP) are two types of scheduling problems that are conceptually similar. The former considers the scheduling of jobs in a plant where there are F identical parallel flow shops, and the latter considers the scheduling of jobs in F factories with an identical flow shop configuration in each one. Both problems can be decomposed into two sub-

AN US

problems: first, assigning each job to one of the F flow shops (or factories); then, scheduling the jobs in each flow shop (or factory) in order to optimize a given criterion. In this paper, the criterion considered is the total tardiness of jobs. The tardiness criterion is a realistic criterion to consider when scheduling jobs in any productive configuration. This is particularly true when there are parallel flow shops, since it helps assign and schedule jobs into lines in order to honour

M

their due dates, which is essential for the survival of companies in the current competitive market.

ED

There are few papers that consider the PFSP, and most of them consider parallel flow shops with only two machines, which is the simplest case. He, Kusiak, and Artiba (1996) proposed using mixed integer programming and an efficient heuristic to deal with this simplest case. In a case

PT

study, Vairaktarakis and Elhafsi (2000) analysed the deterioration in makespan performance when a two-machine hybrid flow shop problem is assimilated with a parallel two-machine flow

CE

shop problem. They justified the design of a parallel flow shop (i.e., independent manufacturing cells) instead of a hybrid flow shop configuration, since the deterioration of the makespan

AC

performance was less than 3%, and the former is easier to manage. They proposed an O(nP3)time dynamic programming algorithm for optimally solving 2 flow shops in parallel with two machines. Cao and Chen (2003) developed a mathematical model and a Tabu Search algorithm for the PFSP with two machines. Al-Salem (2004) proposed a polynomial-time algorithm to minimize the makespan in two-machine parallel flow shops with proportional processing time. Zhang and Van De Velde (2012) developed approximation algorithms with worst-case performance guarantees for scheduling jobs in 2 and 3 flow shops that are in parallel with 2 machines. Finally, Ribas et al., (2017) proposed an iterated greedy algorithm to schedule jobs in identical parallel flow shops without buffers between machines. In the productive configurations

3

ACCEPTED MANUSCRIPT without storage space between machines, a machine can be blocked by the job it has processed if the next machine is not available. Hence, accurate scheduling is necessary for minimizing machine blocking and idle time in order to increase productivity levels. This situation can be found, for example, in the iron and steel industry (Gong et al., 2010); in the treatment of industrial waste and the manufacture of metallic parts (Martinez et al., 2006); or in the use of robotic cells, where a job may block a machine while waiting for the robot to pick it up and move it to the next stage (Sethi et al., 1992).

CR IP T

Despite only one paper having been found for the PFSP with blocking, several procedures have been proposed for the case with only one flow shop in the productive system. The exact methods are: Bounded Dynamic Programming by Bautista et al. (2012) for makespan minimization, Branch and Bound algorithms by Ronconi (2005) and Companys and Mateo (2007) to minimize makespan, by Takano and Nagano(2017) to deal also with sequence dependent setup times, and by Ronconi and Armentano (2001) for the total tardiness minimization; constructive heuristics by

AN US

Companys et al. (2010), Pan and Wang (2012) and Ronconi (2004) to minimize the makespan and by Takano and Nagano (2019) to additionally consider setup times, by Ribas and Companys (2015) to minimize the total flowtime and by Ronconi and Henriques (2009) to minimize the total tardiness; and, finally, several Metaheuristics by Caraffa et al. (2001), TavakkoliMoghaddam et al. (2009), Companys et al. (2010), Gong et al. (2010), Wang et al (2010), Liang

M

et al. (2011), Wang et al (2011), Ribas et al.(2011), Han et al. (2012), Wang et al. (2012), Pan et al. (2013), Ribas et al. (2013), and Companys and Ribas (2014) for makespan minimization, by

ED

Han et al., (2012) and Ribas et al. (2015) to minimize the flowtime and by Ronconi and Henriques (2009), Pan et al., (2013), and Nagano et al. (2017) minimization.

for the total tardiness

PT

On the other hand, as we mentioned before, the PFSP and the DPFSP are two similar problems and because of this the procedures proposed for the DPFSP can be used for the PFSP.

CE

The DFSP was presented by Naderi and Ruiz (2010) and considers the problem of minimizing the makespan by scheduling jobs in several plants, each with an identical flow shop. After this first paper, several authors proposed different heuristics to solve this problem (Fernandez-Viagas

AC

and Framinan, 2014a; Gao et al., 2012; Gao andChen, 2012; Gao, Chen, and Deng, 2013; Gao, Chen, and Liu, 2012; Lin, Ying, and Huang, 2013; Liu and Gao, 2010; Bahman Naderi and Ruiz, 2014; Wang et al., 2013; Xu et al., 2013; Ruiz, Pan, and Naderi, 2018). However, none of the papers from the PFSP or DPFSP have considered the minimization of tardiness. Hence, the motivation of this paper is two-fold. First, our goal is to propose an efficient heuristic for scheduling jobs in a plant with parallel flow shops under the blocking constraint in order to minimize the total tardiness of jobs, since to the best of our knowledge there are not any papers dealing with this same problem. Second, we intend to provide a benchmark and the best obtained solutions that can help other researchers to test and compare other types of algorithms. 4

ACCEPTED MANUSCRIPT The method proposed here is an iterated greedy algorithm that comprises a constructive heuristic and a variable neighbourhood search. Firstly, several constructive heuristics for creating an initial solution were designed and tested. Next, several variants of the proposed algorithm were tested in order to know which parts and characteristics of the algorithm have greater influence on its performance. Finally, a comparison of the proposed algorithm with other adapted algorithms from the literature allowed us to demonstrate its better performance. After this introduction, we define the problem to solve. Next in Section 3, some constructive

CR IP T

procedures are presented. Section 4 details the different parts of the proposed algorithm, Section 5 shows the computational evaluation of the heuristics and the main findings, and finally, Section 6 concludes.

2. Problem description

In this problem n jobs have to be scheduled in one of the F identical flow shops with m

AN US

machines. Each flow shop (factory) is able to process all jobs. The jobs assigned to a flow shop have to be processed by all machines in the same order, from machine 1 to machine m. Each job i, i ϵ {1,2,. . .,n} requires a fixed non-negative processing time pj,i on every machine j, j ϵ {1,2,. . .,m}. Setup times are considered to be included in the processing time. The objective is to schedule the jobs to the different flow shops such that the total tardiness of jobs (TT) is minimized. We denote σf as the sequence of the nf jobs assigned to flow shop f. Therefore, a

M

solution  is formed by the sequence of jobs in each flow shop (=( σ1, σ2, …, σf)).

ED

Next, we introduce additional notation in order to define the mathematical program associated with this problem: let cj,k,f and Tk,f be the departure time and tardiness of a job that occupies position k in machine j at flow shop f, respectively. Notice that, with the blocking constraint, a

PT

job cannot leave the machine until the next machine is free, even if it has finished its operation.

CE

Therefore, according to this notation, the mathematical program can be formalized as follows:

Min TT = [MIN]  Tk,f F

n

1

AC

 xk ,i, f

(1)

i  1, 2,  n

(2)

f 1 k 1 n

 x k ,i , f

1

f  1, 2, F

k  1, 2, n

(3)

i 1

n

c 1,k , f  c1,k 1, f   x k ,i , f  p1,i

k  1, 2,. n

f  1, 2,  F

(4)

i 1 n

c j ,k , f  c j 1,k , f   x k ,i , f  p j ,i

j  2, 3,  , m k  2, 3,. n

f  1, 2,  F

(5)

i 1

5

ACCEPTED MANUSCRIPT c j ,k , f  c j 1,k 1, f

j  2, 3,, m  1 k  2, 3,.n

c j ,0, f  0

f  1, 2,F

j  1, 2, m f  1, 2, , F

c0,k , f  0

k  1, 2,, n F

n

T k , f  c m,k , f   x k ,i , f  d i

(7)

f  1, 2, F

k  2, 3,. n

(6)

(8)

f  1, 2,  F

(9)

f 1 k 1

k  2, 3,.n

f  1, 2,F

(10)

CR IP T

T k, f  0 xk ,i, f 0, 1

k  1, 2,, n i  1, 2,n

f  1, 2, F

(11)

c j ,k , f  0

k  1, 2,, n i  1, 2,n

f  1, 2, F

(12)

AN US

The decision variables are the binary variables xk,i,f ,which take value 1 if job i occupies position k in the sequence of flow shop f. Other variables are the continuous variable cj,k,f and Tk,f .

The objective function is set in equation (1). Constraint set (2) ensures that every job must be exactly at one position and only at one factory. Constraint set (3) ensures that only one job at

M

most can be allocated to each position at a factory. Constraint set (4) defines the departure time of the job that occupies position k in the first machine at factory f. Constraint set (5) specifies

ED

the relationship between the departure times of each job in two successive machines at the assigned factory. Constraint set (6) calculates the departure time of a job under the blocking conditions by considering that the next machine has to be available. Constraint sets (7) and (8)

PT

are the initial conditions. Constraint set (9) defines the tardiness of job that occupies position k in the first machine at factory f. Finally, constraint sets (10), (11) and (12) define the domain of

CE

the decision variables.

Since the problem considered is NP-hard, exact procedures are able to solve only small

AC

instances. Therefore, the next sections propose a heuristics procedure for solving large problems.

3. Constructive procedures Constructive procedures are simple heuristics that allow constructing an acceptable solution with little CPU time. These methods are normally used to obtain an initial solution for heuristics and metaheuristics. To the best of our knowledge, no constructive heuristic has been proposed for the PBFSP or DBFSP to minimize total tardiness; therefore, we present an adaptation of other methods that have been proposed for related problems (e.g., the permutation blocking flow shop 6

ACCEPTED MANUSCRIPT problem to minimize total tardiness and the distributed flow shop problem to minimize the makespan with and without blocking constraints). As stated before, both the PBFSP and the DBFSP need to deal with two related decisions: the allocation of jobs to flow shops (factories) and the sequence of jobs assigned to each line (plant). Hence, there are three types of constructive procedures: (1) those that first create a sequence of all jobs, which is then used to assign the job to each line (plant) according to an allocation rule; (2) those that first select the line (plant) according to a given criterion and then find the best job to be assigned; and, finally, (3)

CR IP T

those methods that, at each iteration, select the best combination of job and line (plant). In this section, we present methods for all three groups: 2 for the first group, 3 methods for the second group and 1 for the third group.

3.1 Methods for the first group

AN US

Naderi and Ruiz (2010) presented the NEH1 and NEH2 heuristic for the DPFSP. Both methods start by ordering the jobs according to the LPT rule. In the next step, NEH1 selects at each iteration the plant with a minimum makespan to allocate the job, which is then assigned to the best position of this plant. On the other hand, at each iteration in NEH2, a job is evaluated at each position of each plant and is then assigned to the best position that leads to the minimum

M

makespan at the plant. Here, to adapt these methods to the total tardiness, jobs are first sequenced according to the Early Due Date (EDD) rule. Next, during the insertion procedure, the best position of a job is selected according to the total tardiness instead of the makespan. To break

ED

ties, the minimum makespan was used. These methods are named NEH1(TT) and NEH2(TT),

PT

respectively.

3.2 Methods for the second group

CE

The methods of this group start by selecting the first job of each line according to index ind1(i). Next, an iterative process starts in order to assign and schedule the remaining jobs to the

AC

lines according to the ind2(i,k,f) index. At each iteration, the line which has the last machine available sooner is selected. The methods presented here are: RC3_m, which is an adaptation of the RC1_m proposed by Ribas, Companys, and Tort-Martorell (2017) for minimizing the PBSP makespan; FPD_m, which is an adaptation of the method proposed by (Ronconi and Henriques, 2009) for minimizing the total tardiness of the BFSP; and A3B2_m, which is an adaptation of the A3B2_m method proposed by Companys and Ribas (2014), also for the total tardiness of the BFSP. The three methods tested have similar structures, which consist of the two steps described as follows:

7

ACCEPTED MANUSCRIPT 

Step 1: For f=1 to F, select the unscheduled job with minimum Ind1(i) and put it in the first position in sequence σf. Set t=F.



Step 2: While t


Step 3: For f=1 to F, use the NEH-based insertion procedure to improve the sequence of

CR IP T

each line and minimize the total tardiness Therefore, the main differences between these methods are how ind1(i) and ind2(i,k,f) are

calculated. Next, the equations are detailed for each of these two indexes and each procedure.

RC3_m

AN US

The first index of this procedure, ind1(i), weights the front delay generated by the job (first term of equation) and the due date (second term). m

ind1(i)    2   (m  j )  p j ,i (m  1)  (1   )  di j 1

(13)

On the other hand, index ind2(i,k,f) weights the timeout (tmi) and the due date of each job.

M

Therefore, once the line is selected, choosing a job that minimizes this index is adequate for minimizing the total tardiness under the blocking constraint. (14)

ED

ind 2(i, k , f )    tmi  (1   )  d i

where, m

tmi   (c j ,k 1, f ( f * i )  c j ,k , f ( f )  p j ,i )

CE

FPD_m

(15)

PT

j 1

The first index of this method is a slight variation of the one used in the original FPD procedure, because in this situation the processing time of jobs in the first machine and the due date are

AC

weighted and the value of the parameter will be calibrated. It considers the processing time of jobs in the first machine and the due date.

ind1(i)    p1,i  (1   )  d i

(16)

The second index is similar to the one proposed for the FPD, but at each iteration it must consider the partial sequence, σf, in each line. m 1

ind 2(i, k , f )   

| b j 1

j

 p j ,i |  min | b j  p j ,i |

max | b j  p j ,i |  min | b j  p j ,i |

 (1   ) 

dynsl i  min dynsl i max dynsl i  min dynsl i

(17)

8

ACCEPTED MANUSCRIPT where, b j  (c j 1 ( f )  c j ( f ))

(18)

is the gap which can be filled by the new operation. Its fit is indicated by the difference between this value and the processing time of the job to be scheduled. Additionally, m

dynsl i  (d i   p j ,i  c1 ( f )

(19)

j 1

A3B2_m

CR IP T

is an approximation to the dynamic slack of job i, as defined in Ronconi and Henriques (2009).

In this method, ind1(1) weights the normalized total processing time of the job with its normalized slack.

Pi  Pmin sl  sl min  (1   )  i Pmax  Pmin sl max  sl min

AN US

ind1(i)   

(20)

where Pi is the sum of the processing time of job i, and Pmin and Pmax are the minimum and maximum total processing times of jobs; sli is the slack of job i; and slmin and slmax are the minimum and maximum job slacks, respectively. The second index is calculated as in (21):

tmi  tmmin dynsl 'i dynsl ' min  (1   )  tmmax  tmmin dysl ' max dynsl ' min

M

ind 2(i)   

(21)

ED

where tmi is the timeout generated when job i is scheduled at the end of the partial sequence σf at position k+1; tmi is the timeout calculated as in (15); and dynsl’i is the real dynamic slack of job i when it is scheduled at the end of σ, and it is calculated as in (22)

PT

dynsl'i  d i  cm,i

(22)

CE

3.3 Method for the third group The method for this group is named RC4.

AC



Step 1: For f=1 to F, assign the first job according to ind1(i), which is calculated by equation (13). Set t=F.



Step 2: While t < n, select the job and line simultaneously in order to minimize the second index, which is calculated by equation (23). t=t+1.

ind 2(i, k )    d i  (1   )  (cm,k 1, f ( f * i)  c0 )

(23)

where c0 is the minimum makespan among all lines before the schedule of the new job i. 

Step 3: For f=1 to F, use the NEH-based insertion procedure to improve the sequence of each line in order to minimize the total tardiness. 9

ACCEPTED MANUSCRIPT Notice that all these methods use two parameters that will be calibrated before the computational evaluation.

4. Improvement heuristics An iterated greedy algorithm is the type of algorithm chosen to design an efficient procedure for solving the problem at hand. This type of algorithm has proven to be very effective for scheduling jobs under several productive configurations and constraints. In particular, Fernandez-

CR IP T

Viagas and Framinan (2014) proposed this type of algorithm for dealing with the DFSP, and Ribas et al., (2017) did so for the PBFSP. Although the performance criterion considered in both cases was minimization of makespan, some ideas can be taken from both algorithms. However, as the tardiness criterion affects the structure of the problem, some considerations have to be

AN US

taken into account.

The algorithm starts with an initial solution, which is improved upon by applying, to the sequence of each line, a variable neighbourhood search (VNS) based on swap and insert neighbourhood structures. Next, the Reassignment and Permutation procedures are applied in order to try and improve the current solution. The Reassignment procedure reassigns a job from the critical line to another randomly selected line. On the other hand, the Permutation procedure

M

tries to swap two jobs, one from the critical line and another from another randomly selected line. Then the solution goes to the main part of the algorithm, where three procedures are

ED

iteratively used, the solution is improved upon, and the CPU time does not reach its time limit. These procedures are: Perturbation, which modifies the current solution ; and the

AC

CE

PT

IGA Algorithm : = Generate_Initial_solution := VNS()  := Reassignment and Permutation () *:=  TT*:=∑Ti while CPUtime < timelimit do if TT*=0 then exit do  := Perturbation () if ∑Ti < TT* then *:=  TT*=∑Ti end := Reassignment and Permutation () if ∑Ti > TT* and random>0.5 then := * end wend end

10

ACCEPTED MANUSCRIPT Reassignment and Permutation functions, which intensify the search to find better solutions in the insert and swap neighbourhoods. Finally, the current solution is compared with the best solution found so far. If it is worse, the former is kept with a probability of 50%. Figure 1 shows the general outline of the proposed algorithm. Notice that if the total tardiness (TT) is 0 (i.e., an optimal solution is obtained), the algorithm finishes with this solution. Figure 1. Outline of the IGA algorithm

CR IP T

The outline of the Reassignment and Permutation modules, which manage the insertion and swap movements between lines, is shown in Figure 2. Notice that the Reassignment and Permutation procedures are applied while the solution is improved. Index nml is used to access the Permutation procedure during the first iteration if the Reassignment procedure does not improve the solution.

AN US

In the following sections, the different parts of the algorithm are explained.

ED

M

Reassignment and Permutation () nml=0 do if TT*=0 then exit do nml=nml+1 ’:=Reassignment () if TT(’) = TT() and nml>1 then exit do  = ’ :=Permutation() if TT(’) = TT() then exit do  = ’ loop end

PT

Figure 2. Outline of the Reassignment and Permutation procedures

CE

4.1 Variable Neighbourhood Search procedure This procedure is similar to the one proposed in Ribas et al. (2017) but has been adapted to

AC

the total tardiness criterion. The VNS uses two local searches sequentially (named LS1 and LS2), which explore the swap and insert neighbourhoods of the current sequence at each flow shop (factory). The neighbourhood to be explored first is selected randomly. In LS1, for each job in the sequence, neighbours are generated by swapping a job with all jobs that follow it in the sequence. If the best neighbour (σf′) is better than the current solution (σf), it becomes the new current solution σf, and the process continues until all jobs have been considered. To prevent the neighbourhoods from always being explored in the same order, the jobs are selected randomly.

11

ACCEPTED MANUSCRIPT In contrast, in LS2, for each job in the sequence, neighbours are generated by removing the job from its position and inserting it into all other possible positions. If the best neighbour (σf ′) is better than the current solution (σf), it becomes the new current solution σf, and the process continues until all jobs have been considered. As in LS1, jobs are selected randomly. After exploring the neighbouring solutions of the current solution σf, the local optimum σf ′ is compared with σf. If the solution has been improved, σf ′ replaces σf and the search continues in the other neighbourhoods. This process continues until the current solution is no longer improved of the solution. If TT(σf ′) is less than TT(σf *), σf′ replaces σf *.

4.2 Perturbation mechanism

CR IP T

upon. Next, the local optimum σf′ is compared with the best solution σf * in terms of the quality

This mechanism helps the algorithm escape from local minima by perturbing the current

AN US

solution, which thus allows the algorithm to jump and explore other solution spaces. The implemented procedure consists of randomly removing d jobs from their lines and selecting the best line and position in which to schedule them in order to minimize the total tardiness of jobs. Therefore, each position is evaluated for each line and the job is scheduled at the position that leads to the minimum increment in tardiness. Finally, the selected line and position in which to

M

schedule each job is the one that leads to a minimum increment in the total tardiness.

ED

The outline of this procedure can be seen in Figure 3.

AC

CE

PT

Perturbation Mechanism_IG () 1:= remove d jobs randomly without repetition from their line for k=1 to d i:= job of position k in 1 ΔTT0=infinite for f=1 to F insert i in the best position of f Calculate ΔTTf if ΔTTf
Figure 3. Perturbation mechanism of IGA

4.3 Reassignment procedure 12

ACCEPTED MANUSCRIPT

The reassignment procedure starts by selecting nlimr jobs, without repetition, from the critical line (fc). The jobs are kept in 1. The critical line can be defined in two different ways: (1) it is the line where the total tardiness of jobs assigned to this line is maximum; or (2) it is the line that has the job with maximum tardiness. These two strategies define two variants of the algorithm. Next, the first job in 1 is removed from the critical (fc) line and the algorithm tries to insert it into the best position of another line (f), which is selected randomly. The change is evaluated by

CR IP T

comparing the diminishment in the tardiness of the fc line with the increase in tardiness of the line where the job has been inserted. If the balance is favourable, the change is kept and the process re-starts by detecting the new critical line, since it may have changed; if the movement is not favourable, the job is not reassigned and the process continues by trying to insert the next selected job in σ1. The process finishes when all the nlimr jobs from fc have been tested and the

AC

CE

PT

ED

M

AN US

solution has not been improved. The outline of this procedure is shown in Figure 4.

13

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

AN US

CR IP T

Reassignment() ’ =  do fc=Detect_ fcrit(’) 1 =select(’,fc,nlimr) k=0 do move=0 k=k+1 if k>1 then exit do i=1(k) Tfc = tardiness variation in fc by extracting job i f=0 do f=f+1 if fF and ffc insert i in the best position of f Tf=tardiness variation in f due to insertion of i if Tfc + Tf < 0 then move=1 if move=1 then exit do endif loop if move=1 then exit do loop if move=0 then exit do ’:= extract i from fc and insert in the best position of f if TT(’)
Figure 4. Outline of the reassignment procedure

4.4 Permutation procedure The permutation procedure consists of swapping two jobs, one from the critical line (fc) and another from another randomly selected line (f). Here again, the two strategies for selecting the

14

ACCEPTED MANUSCRIPT critical line have been implemented: (1) the line with the maximum total tardiness of jobs assigned; (2) the line which has the job with maximum tardiness.

AC

CE

PT

ED

M

AN US

CR IP T

Procedure Permutation() ’ =  do fc=Detect_ fcrit(’) 1 =select(’,fc,nlimp) k1=0 do move=0 k1=k1+1 if k1>1 then exit do i=1(k) f=0 do f=f+1 if fF and ffc 2 =select(’,f,nlimp) k2=0 do move=0 k2=k2+1 if k2>2 then exit do h=2(k2) Tfc = tardiness variation in fc by extracting job i and inserting job h in the best position Tf = tardiness variation in f by extracting job h and inserting job i in the best position if Tfc + Tf < 0 then move=1 if move=1 then exit do loop endif if move=1 then exit do loop if move=1 then exit do loop if move=0 then exit do extract i from fc and insert h in the best position extract h from f and insert i in the best position the new solution is ’ if TT(’)
Figure 5. Outline of permutation procedure.

15

ACCEPTED MANUSCRIPT The process starts by detecting the critical line and selecting nlimp jobs, without repetition, which are kept in σ1. Next, nlimp jobs are selected from another line and kept in σ2. Both the lines and jobs are randomly selected. Then, the first job i from σ1 and the first job h from σ2 are removed from their lines. It is then attempted to insert the job from σ1 into the best position of line f, i.e., the position which leads to minimized tardiness variation after removing job i and inserting job h in this line. Next, it is attempted to insert h into the best position of line fc. Notice that the tardiness variation is the difference between the original tardiness of the line and the

CR IP T

tardiness after swapping the jobs. Each line can have a positive, negative or null variation. Then, if the sum of tardiness variation of each line is positive, it means that the solution has not improved. Therefore, the tentative swap of jobs is not done and the process is restarted by trying to swap the job from σ1 with the next job from σ2. This process continues for all jobs in σ2 for as long as the solution is not improved upon. Once all jobs in σ2 have been tested without improvement, the process restarts by selecting the next job from σ1, which the process again tries

AN US

to swap with the jobs from σ2 until either an improvement is found or the algorithm has tried to swap all the jobs from σ1 with all the jobs from σ2.

If the solution improves, i.e., the sum of tardiness variations is negative, the change is kept and the process restarts from the beginning by detecting the critical line and selecting a set of

M

nlimp jobs from it.

The outline of this procedure is shown in Figure 5.

ED

5. Computational Evaluation

This section presents the computational evaluation of the constructive procedures and

PT

heuristics methods. All the algorithms were coded in the same language (QB64) and tested on

AC

CE

the same computer, a 3 GHz Intel Core 2 Duo E8400 CPU with 2 GB of RAM.

16

ACCEPTED MANUSCRIPT After considering the number of variables related to the initial solution (λ, μ,); the algorithm we want to optimize for the proposed method; the levels of n, m and F that we consider necessary to explore; and the time needed to run each experiment, we chose to use a sequential experimentation strategy, as recommended by Box et al. (2009) and Montgomery (2013). Firstly, the constructive procedures were calibrated and evaluated in order to select the best methods as well as to find a good compromise between the quality of the solution and the CPU time required to obtain the initial solution of the heuristics. Next, we built, calibrated and

CR IP T

evaluated some variants of the presented heuristic, not only to evaluate the influence of the quality of the initial solution but also the strategies for selecting the critical line. Finally, we compared the best variant against two benchmark algorithms proposed for the DPFSP that were adapted to the total tardiness parallel blocking flow shop problem.

It is worth noting that the benchmark sets used to adjust the constructive procedure and for

AN US

calibration are different from those used to compare the heuristics. To compare the heuristics, we used the adaptation of Taillard’s benchmark (Taillard, 1993) proposed by Naderi and Ruiz (2010) to which due dates were added whereas we generated, in the same way, a specific benchmark set for the calibration.

5.1 Calibration of Constructive Procedures

M

As stated before, RC3_m, FPD_m, A3B2_m and RC4 all have two parameters that must be calibrated. Parameters λ and µ of each heuristic were selected by measuring the performance of

ED

the algorithm, which was done by combining several λ and µ values. The values of λ and µ ranged from 0 to 1 in increments of 0.05. Therefore, we tested 21 values for each parameter.

PT

For this test, we generated 100 instances, grouped into 20 sets of size n x m (5 instances per set), where n= {25, 50, 100, 200, 400} and m = {5, 10, 15, 20}. The processing times of jobs were uniformly distributed between [1, 99]. All these 100 instances were considered with a different

CE

number of factories, F= {2,3,4,5,6,7}, which gave us 600 instances in total. Since there are no previous studies for scheduling jobs in parallel lines (or distributed flow

AC

shops) to minimize the total tardiness, we developed a method for generating the due dates of each job in these contexts. The due dates were generated according to a uniform distribution between [α·LB, β·LB], being α=1-T-R/2 and β=-T+R/2 as defined in Potts and Van Wassenhove (1982) and where T and R are, respectively, the tardiness factor of jobs and the dispersion range of due dates, while LB is the lower bound for the makespan. By giving different values to T and R, we can obtain different scenarios. In this case, as in Potts and Van Wassenhove (1982), we define 4 scenarios: (1) T=0.2, R=0.6; (2) T=0.2, R=1.2; (3) T=0.4, R=0.6; (4) T=0.4, R=1.2. The lower bound (LB) has been calculated by taking the one proposed in (Companys and Mateo, 2007) for the blocking flow shop problem and multiplying it by (n/F+m-1)/(n+m-1), because the 17

ACCEPTED MANUSCRIPT critical path that defines the makespan is composed of n+m-1 operations. Notice that when m=1, this factor is 1/F, which means the lower bound is divided by the F lines. Hence, by considering the 600 generated instances for the 4 due date scenarios, we obtain a total of 2400 instances for calibrating the heuristics. The performance was measured by the relative deviation index (RDI), calculated as (24) for each procedure: Heurh , s  Best s

(24)

CR IP T

RDI =

Worst s  Best s

where Heurh,s is the average of tardiness values obtained by heuristic h in instance s, and Bests and Worsts are the minimum and maximum values of tardiness obtained for this instance from among all combinations of parameters.

AN US

For each procedure, we conducted an Analysis of Variance, including in the model the terms: F, μ, λ, F*μ, F*λ and μ* λ. This allowed checking the effects of the two parameters as well as determining via their interactions with F whether or not their best values were dependent on the number of lines. In all four cases, the result was a dominating, very strong and highly significant effect of μ and significant but very weak effects of λ, F*μ and μ* λ. As an example, we show the results for the A3B2 procedure. Figure 6 shows the dominating μ main effect. Figure 7 shows an

M

enlarged plot of the best region for the interaction F*μ, and the best μ value, which in this case

AC

CE Figure 6. Main effect of μ

Contour Plot of RDI vs F; MU 7

0.41

6

0.39

0.39

0.41

5

F

PT

ED

is 0.65.

4

3

0.42

0.40

0.42 0.40

2 0.40

0.45

0.50

0.55

0.60

μ

0.65

0.70

0.75

0.80

Figure 7. Blown up of the F*μ interaction for 0.4>μ>0.8

Figure 8 shows the F*λ interaction. The best value for λ is 0.65 and the interaction is very weak, especially when λ takes values above 0.5.

18

ACCEPTED MANUSCRIPT Figure 9 shows the μ* λ interaction. Although it is also weak, it reaches a plateau at RDI values

CR IP T

of 0.35

Figure 9. μ* λ interaction

Figure 8. F*λ interaction.

AN US

Although the pattern has different λ and μ values, it is basically the same for each procedure, especially regarding the small influence of the interactions. Thus, we have been able to find

M

(Table 1) the best values of λ and μ for each procedure.

λ

μ

A3B2

0.65

0.65

RC3_m

0.75

0.65

FPD_m

0.8

0.25

RC4

0.8

0.35

AC

CE

PT

ED

Constructive Procedures

Table 1. Value of parameters for each procedure

5.2 Computational Evaluation of Constructive Heuristics An adaptation of the well-known Taillard benchmark (Taillard, 1993) was used to compare the DPFSP, as proposed by ( Naderi and Ruiz, 2010). The due dates for these test beds were calculated as explained in Section 5.1, and they can be consulted at the following link: http://hdl.handle.net/2117/123116.

19

ACCEPTED MANUSCRIPT The benchmark comprises 12 sets of 10 instances, ranging from 20 jobs and 5 machines to 500 jobs and 20 machines, and where n ϵ {20, 50, 100, 200, 500} and m ϵ {5, 10, 20}, although not all combinations of n and m are available. In particular, sets 200x5, 500x5 and 500x10 are missing. These 120 instances were augmented with six values of F ϵ {2, 3, 4, 5, 6, 7} and four due dates, one per each scenario, giving us a total of 2880 instances. In this test, the RDI index is calculated by considering as Heurh,s the tardiness value obtained by constructive procedure h in instance s, and as Bests and Worsts the respective minimum and

CR IP T

maximum values of tardiness obtained for this instance from among all compared procedures. Figure 10 shows the 95% confidence interval for the mean RDI value of each heuristic. It is very clear that there are two groups: The three heuristics FPD+, NEH1 and RC4+ are much worse than the other three heuristics A3B2+, NEH2 and RC3_m+. Two other results are also clear: NEH2 is the one that provides, by far, the best RDI; and RC3_m+ is not significantly better than

AN US

Interval Plot of RDI 95% CI for the Mean

0,8 0,7

0,5 0,4

M

RDI

0,6

0,3 0,2

ED

0,1

A3B2+

FPD+

NEH1

NEH2

RC3_m+

RC4+

Heuristic

Individual standard deviations are used to calculate the intervals.

PT

A3B2.

Figure 10. 95% confidence interval for the mean RDI value

CE

However, it is worth mentioning that NEH2 has the inconvenience of being much slower than

AC

A3B2+ and RC3_m+, as can be seen in Table 2. F

RC3_m+

A3B2+

RC4+

FPD+

NEH1

NEH2

2

181.961

186.801

192.242

186.340

171.422

339.223

3

87.777

93.258

107.040

92.398

77.391

224.320

4

55.090

60.203

82.238

59.871

44.597

168.508

5

39.402

45.148

75.250

44.661

29.492

135.121

6

32.450

36.852

75.019

36.359

21.371

112.980

7

26.851

31.918

77.832

31.473

16.528

97.110

Average

70.589

75.697

101.604

75.184

60.134

179.544

Table 2. Average CPU time per heuristic and number of lines, in seconds

20

ACCEPTED MANUSCRIPT We observed, as Ruiz et al. (2018) also remark, that in IGA the effect of the initial solution is quickly neutralized by the local search methods. Hence, it is sometimes more effective to reserve time for doing more iterations than to spend it generating the initial solution. On the other hand, in some cases the structure of the initial solution leads to the algorithm becoming trapped in a local optimum from which it is difficult to escape. Hence, we decided to test the performance of the algorithm with different initial solution procedures. In particular, we tested three variants: one

5.3 Experimental Parameter Adjustment of IGA

CR IP T

with RC3_m+, another with A3B2_m+ and finally the third with NEH2.

The next step was to calibrate the IGA’s parameters, i.e.: d, the number of jobs extracted during the perturbation phase; and nlim, the number of iterations for the reassignment and permutation movements of jobs in the improvement phase. We tested four values of d (2, 3, 4, 5) and three of nlim (10, 15, 20). The process was done for the two variants implemented for choosing the

AN US

critical line: (1) the line where the total tardiness of jobs assigned to this line is maximum or (2) the line which has the job with maximum tardiness. Here, as an example, we show the process for the first variant. The set of instances used in this test were the same as in the calibration of constructive procedures, i.e., the 2400 instances that resulted from combining n, m and F values with the due dates of the four tardiness scenarios. The algorithms were stopped after 20·n2·m

M

milliseconds. The performance of the algorithms was measured with the RDI index calculated by (24). In this case, Heurh,s was the average of tardiness values obtained in 5 runs by the heuristic h in instance s, and Bests and Worsts represent, respectively, the minimum and

ED

maximum values of tardiness obtained for this instance from among all algorithms compared in any run.

PT

Given some previous experiences, we were assuming that nlimp=nlimpr and, further, that this factor would not interact with the remaining factors. Thus, we designed an experiment in which

CE

the numbers of iterations in the reassignment and permutation procedures were the same (here named nlim in order to avoid misunderstandings). To our surprise the results of the first experiment showed that all factors and their pairwise interactions were significant. The huge

AC

difference in contribution to total variability between the main effects and their interactions (see the Mean Square Error in Table 3) suggests that the main effects dominate the interactions and the results can be interpreted by omitting them, as Figures 11 and 12 confirm.

Source F n m d nlim

DF 5 3 3 3 2

SS 17.67537 14.55136 22.07832 1.09691 0.54040

MS 3.53507 4.85045 7.35944 0.36564 0.27020

F 166.19 228.02 345.97 17.19 12.70

P 0.000 0.000 0.000 0.000 0.000

21

ACCEPTED MANUSCRIPT

CR IP T

15 18.90656 1.26044 59.25 0.000 F*n 15 1.10924 0.07395 3.48 0.000 F*m 15 2.44949 0.16330 7.68 0.000 F*d 10 1.89568 0.18957 8.91 0.000 F*nlim 9 6.26954 0.69662 32.75 0.000 n*m 9 1.60863 0.17874 8.40 0.000 n*d 6 1.06537 0.17756 8.35 0.000 n*nlim 9 0.62069 0.06897 3.24 0.001 m*d 6 0.31022 0.05170 2.43 0.024 m*nlim 6 0.70532 0.11755 5.53 0.000 d*nlim 22923 487.60874 0.02127 Error 23039 578.49183 Total Table 3. ANOVA table for F, n, m, d and nlim

The 95% confidence interval plot for the d main effects (Figure 11) shows that the differences between d=3 regarding d=2 and d=5 are highly significant. The difference from d=4 is also significant, although by a small margin. In any case, it is clear that d=2 is better than the rest.

AN US

Furthermore, let us look at the two bigger interactions (although, as we have seen, they are much smaller than the main effects) involving d, the d*F and d*n interactions. It is clear that, for all values of F and n, d=3 always gives either the best RDI value or one that is not significantly different from the best (Figure 12).

Interval Plot of rdi 95% CI for the Mean

M

0.460

0.455

0.445

PT

0.440

ED

rdi

0.450

0.435

CE

2

3

d

4

5

Individual standard deviations are used to calculate the intervals.

AC

Figure 11. Confidence interval for RDI by d values

22

ACCEPTED MANUSCRIPT

Interval Plot of rdi

Interval Plot of rdi

95% CI for the Mean

95% CI for the Mean

0.52

0.50

0.50

0.48

0.48

0.46

0.44

rdi

rdi

0.46

0.44

0.42 0.42

0.40

0.40

0.38 0.36 F d

2 3 4 5 6 7 2

2 3 4 5 6 7 3

2 3 4 5 6 7 4

2 3 4 5 6 7 5

Individual standard deviations are used to calculate the intervals.

n d

25 50 100 200 2

25 50 100 200 3

25 50 100 200 4

25 50 100 200 5

CR IP T

Individual standard deviations are used to calculate the intervals.

Figure 12. Confidence interval for RDI grouped by F and d values and by n and d, respectively.

AN US

It is also clear that nlim=15 is overall the best value, independently of F, n and m (Figure 13).

Interval Plot of rdi 95% CI for the Mean

0.455

0.445

M

rdi

0.450

0.440

0.435

15

20

nlim

ED

10

Individual standard deviations are used to calculate the intervals.

Only slightly significant interactions appear when nlim=10.

PT

Figure 13. Confidence interval for RDI by nlim values

Interval Plot of rdi

Interval Plot of rdi

95% CI for the Mean

95% CI for the Mean

0.49

CE

0.500

0.48

0.475

0.47 0.46

AC

0.45

rdi

rdi

0.450

0.425

0.43 nlim 10 15 20

0.400

F

nlim

2

3

0.44

4

5 10

6

7

2

3

4

5

6

15

Individual standard deviations are used to calculate the intervals.

7

2

3

4

5 20

6

7

0.42

nlim 10 15 20

0.41 0.40 n nlim

25

50

100 10

200

25

50

100

200

25

50

15

100

200

20

Individual standard deviations are used to calculate the intervals.

Figure 14. Confidence interval for RDI showing interactions of parameters

These interactions made us think that the behaviour of the reassignment and permutation movements could be different, or even perhaps that the best value of this parameter could be different in both procedures. Therefore, we performed another experiment with the 23

ACCEPTED MANUSCRIPT consideration that the values of nlimr and nlimp could be different. The results showed that the improvements when nlim =15 were due to nlimp, i.e., the permutation movements. This can be seen in Table 4, where aside from n, m and F, the only significant factors are nlimp and its interactions with n and m. Of these, the nlimp*m interaction is very weak, and the nlimp*n interaction is very strong, as indicated by their MSE. SS 8.020 11.517 11.904 0.075 1.437 0.034 0.086 0.060 7.046 0.581 0.346 0.069 406.893 448.067

MSE 1.60399 3.83906 3.96786 0.03745 0.71844 0.00566 0.01430 0.00596 1.17436 0.09676 0.03461 0.01715 0.02363

F 67.87 162.43 167.88 1.58 30.40 0.24 0.60 0.25 49.69 4.09 1.46 0.73

P 0.000 0.000 0.000 0.205 0.000 0.964 0.727 0.991 0.000 0.000 0.146 0.574

CR IP T

DF 5 3 3 2 2 6 6 10 6 6 10 4 17216 17279

AN US

Source F n m nlimr nlimp n*nlimr m*nlimr F*nlimr n*nlimp m*nlimp F*nlimp nlimr*nlimp Error Total

Table 4. Analysis of Variance for RDI

Figure 15 shows that nlimr does not affect RDI whereas nlimp =15 is significantly better

M

than the others. Therefore, we decided to set nlimr=10 (reassignment movement) in order to

ED

reserve CPU time for more improvement steps.

Interval Plot of rdi 95% CI for the Mean

PT

0.47

CE

0.46

AC

rdi

0.45

0.44

0.43

nlimr nlimp

10

15 10

20

10

15 15

20

10

15

20

20

Individual standard deviations are used to calculate the intervals.

Figure 15. Confidence interval for RDI grouped by nlimr and nlimp

24

ACCEPTED MANUSCRIPT 5.4 Computational Evaluation of the Heuristics Once the two variants of IGA were calibrated, we performed the computational evaluation. In this test, we evaluated the combination of the two variants with the three selected methods to obtain the initial solution: RC3_m, A3B2 and NEH2. To distinguish between each variant, we gave the algorithms the name of the initial solution plus “_1” to indicate that the critical line is the one with the maximum total tardiness of jobs assigned to it, or plus “_2” to indicate that the critical line is the one that has the job with the maximum tardiness.

CR IP T

A comparison was made again of the 2280 instances that were used to compare the constructive procedure, i.e., the Taillard instances that had been adapted to the problem. Once again, 5 runs per instance were conducted. The algorithms were stopped after a CPU time of 30·n2·m milliseconds. The performance of the algorithms was measured with the RDI index, calculated by (24). In this case, Heurh,s is the average of tardiness values obtained by the heuristic h in instance s after 5 runs; and Bests and Worsts are the minimum and maximum values of tardiness

AN US

obtained for this instance, in any run and from among all compared procedures.

As can be seen in Figure 16, A3B2_1 is the best-performing method, with a huge difference over the others. This confirms that, for this problem, using the best initial solution on hand does not always lead to obtaining the best final solution and, subsequently, it is probably necessary to make a compromise between quality and the required CPU time. Remember that the best

M

performing constructive procedure was NEH2, but it was also the procedure that needed higher CPU time to obtain a solution. Hence, it is recommended to select an initial solution procedure

ED

that provides a quality solution but without huge computational effort. Additionally, we recommend setting as the critical line the one with higher tardiness, since our results show that the procedures with this setting always lead to similar or even better results than the procedures

PT

that consider the critical line to be the one that has the job with the highest tardiness.

CE

Interval Plot of RDI 95% CI for the Mean

0.43

0.41

0.40 0.39

RDI

AC

0.42

0.38 0.37 0.36 0.35 0.34 A3B2_2

A3B2_1

NEH2_2

NEH2_1

RC3_m_2

RC3_m_1

Heuristic-Method Individual standard deviations are used to calculate the intervals.

Figure 16. Confidence interval for RDI by variants of IGA’s

25

ACCEPTED MANUSCRIPT Finally, it is interesting to note that A3B2_1 is best or tied to best overall, no matter the values of n, m or F – with one exception: when n=20, it loses to NEH2_2, as shown in Figure 17.

Interval Plot of RDI 95% CI for the Mean

0.6

CR IP T

0.5

RDI

0.4

0.3

0.2

0.1 n

20 50 100 200 500

Heuristic-Method

_2 B2 A3

20 50 100 200 500

N

20 50 100 200 500

2 2_ EH

N

2 0 50 100 200 500

1 "_ EH

_2 m 3_ RC

20 50 100 200 500

RC

_1 m 3_

AN US

A

20 50 100 2 00 5 00

1 2_ 3B

Individual standard deviations are used to calculate the intervals.

Figure 17. Confidence interval for RDI by procedure and n

Finally, to the best of our knowledge, this is the first paper dealing with the parallel blocking flow shop problem with the goal of minimizing total tardiness. To this end, we took two

M

algorithms that had been previously proposed for solving the distributed flow shop problem for minimizing makespan and adapted them to the problem on hand. Specifically, these were the

ED

Bounded–Search iterated greedy algorithm proposed in Fernandez-Viagas and Framinan (2014) and the Scatter Search of Naderi and Ruiz (2014). We named these modified algorithms B_IGA and SS. The first adaptation of the B_IGA consisted of changing the ordering criterion of NEH2

PT

from the Largest Processing Time rule (LPT) to the Early Due Date (EDD) rule. Next, in the improvement phase, the critical line selected was the one with the highest tardiness. Finally, the

CE

evaluation criterion was changed to the total tardiness. As the authors proposed, we generated 25 solutions for the SS, and from them we kept 10 to serve as the initial population. These 25 solutions were generated with NEH2 method but with changing the initial sequence. One used

AC

EDD as the ordering rule while all the others considered one generated randomly. Interval Plot of A3B2; IGAS_l; SS_l 95% CI for the Mean

0.8 0.7 0.6

rdi

0.5 0.4 0.3 0.2 0.1 Lines

26 2

3

4

5

A3B2_l

6

7

2

3

4

5

B_IGA

Individual standard deviations are used to calculate the intervals.

6

7

2

3

4

5 SS

6

7

ACCEPTED MANUSCRIPT Figure 18. Confidence interval for RDI by algorithm

The comparison of these two algorithms against A3B2_l was done using the adapted Taillard’s instances to our problem. Figure 18 shows the RDI confidence interval for the three algorithms. From this figure it is clear that the proposed A3B2_l outperforms the other two algorithms. Moreover, this algorithm is very stable, since its behaviour is not dependent on the number of lines, as can be seen in Figure 19.

95% CI for the Mean

0.7

0.6

rdi

0.5

0.4

AN US

0.3

CR IP T

Interval Plot of A3B2_l; IGAS; SS

0.2 A3B2_l

B_IGA

SS

Individual standard deviations are used to calculate the intervals.

Figure 19. Confidence interval for RDI by algorithm and number of lines

To make clear which features of the proposed algorithm have greater influence on the good

M

behaviour, we compare and describe the main differences among them.

The SS uses only the insertion neighbourhood to improve the solution. According to our

ED

experiments, this neighbourhood, which consists of reassigning a job from the critical line to the best position of another line (Reassignment procedure), is less effective than the Permutation

PT

procedure. Furthermore, SS consumes a lot of time to generate the initial population (25 solutions generated in order to keep only the 10 best) with a procedure that requires more CPU

CE

time than the method used in A3B2_l. Therefore, A3B2_l has more CPU time for iterating and intensifying the search to find better solutions.

AC

On the other hand, the initial solution procedure in B_IGA is NEH2, which we have already demonstrated in the computational evaluation section is more time-consuming than A3B2, although it obtains better solutions. Moreover, the reassignment movements of B_IGA are exhaustive, which means that jobs from the critical line are tested in all positions of the other lines and the best solution is kept. However, in A3B2_1, the first best solution found is kept and the process restarts by selecting another job from the new critical line – if it has changed. Finally, in B_IGA, the permutation movements (which are more effective than the reassignment movements) are made only if n < F·20, as this avoids expending too much CPU time, which is critical because this process is also performed exhaustively and thus implies a huge amount of CPU time. Conversely, in A3B2_1, a greater number of iterations is set for the Permutation 27

ACCEPTED MANUSCRIPT process than for the Reassignment one in order to afford more opportunities for improving the solution. The process restarts when a better solution is found. Hence, we recommend using initial solution procedures that allow attaining solutions with certain quality while maintaining moderate CPU time, thus giving the algorithm more time to iterate and find better solutions. In addition, we recommend using a non-exhaustive variable neighbourhood search that performs permutation and reassigns jobs between lines in order to

are more effective.

CR IP T

improve the solution while offering more opportunities for permutation movements, since they

Finally, the best solutions obtained for the Taillar’s benchmark to the Parallel blocking flow shop problem can be consulted in : http://hdl.handle.net/2117/123116. These solutions can help other researchers to compare other algorithms for this problem.

AN US

6. Conclusions

This paper proposes an Iterated Greedy Algorithm to solve both the parallel blocking flow shop problem and the distributed blocking flow shop problem by minimizing the total tardiness of jobs. To the best of our knowledge, there is no paper that proposes a heuristic for scheduling jobs in these two environments in a way that minimizes total tardiness; therefore, we began our

M

investigation by designing constructive procedures that considered the characteristics of the problem and allowed us to obtain a good initial solution for the IGA.

ED

Our preliminary tests showed us that the quality of the initial solution does not determine the quality of the final solution reached by the IGA. This is because the effect of the initial solution

PT

in these types of algorithms is quickly neutralized by the local search methods. Hence, it is sometimes more effective to reserve time for more iterations than to spend it on generating the initial solution. Therefore, we selected three constructive procedures by considering not only

CE

the solution quality but also the CPU time. We confirmed the initial findings: the best performing algorithm uses the A3B2 procedure, which was especially designed for this problem

AC

and which is less time consuming than NEH2. However, A3B2 was outperformed by NEH2 when testing the constructive procedures. On the other hand, we have shown that interchanging jobs between the critical line and the others during the improvement phase (i.e., the permutation procedure) is more effective than reassigning jobs from the critical line to another, even though it is recommended that both should give higher weight to the permutation movements. Moreover, we recommend selecting the line with the highest job tardiness as the critical line.

28

ACCEPTED MANUSCRIPT Finally, a set of instances with due date are provided for the parallel flow shop and distributed flow shop. What is more, we present the best solution found for the total tardiness parallel flow shop blocking problem. Both the instances and the best solutions can help other researchers who wish to contrast new procedures. Future lines of research can consider setup times or the possibility of using a uniform or unrelated parallel flow shop, which is often the case when companies buy new lines to increase

most industrial environments. Acknowledgement

CR IP T

capacity and they outperform the older lines. Both considerations link the problem closely to

This work was partially supported by the Spanish Ministry of Science and Innovation under the project SCHEYARD with reference DPI2015-65895-R and by the PROREC program of the

AN US

business administration department of UPC.

CRediT author statement

CE

References

PT

ED

M

Imma Ribas: Conceptualization, Methodology, Investigation, Formal Analysis, Writing- Original draft preparation. Ramon Companys: Conceptualization, Investigation, Software, Data curation. Xavier Tort-Martorell: Formal Analysis

AC

Al-Salem, A., 2004. A Heuristic to Minimize Makespan in Proportional Parallel Flow Shops Ameer Al-Salem Pages 98-107 A Heuristic to Minimize Makespan in Proportional Parallel Flow Shops. Int. J. Comput. Inf. Sci. 2. Bautista, J., Cano, A., Companys, R., Ribas, I., 2012. Solving the Fm|block|C max problem using Bounded Dynamic Programming. Eng. Appl. Artif. Intell. 25. https://doi.org/10.1016/j.engappai.2011.09.001 Box, G., Hunter, J.S., Hunter, W., 2009. Statistics Experimenters. Book 1–655. Cao, D., Chen, M., 2003. Parallel flowshop scheduling using Tabu search. Int. J. Prod. Res. 41, 3059–3073. https://doi.org/10.1080/0020754031000106443 Caraffa, V., Ianes, S., Bagchi, T., Sriskandarajah, C., 2001. Minimizing makespan in a blocking flowshop using genetic algorithms. Int. J. Prod. Econ. 70, 101–115. https://doi.org/DOI: 10.1016/S0925-5273(99)00104-8

29

ACCEPTED MANUSCRIPT Companys Pascual, R., Ribas Vila, I., 2014. A New Constructive Heuristic for the Fm|block|ΣT. Springer, Cham, pp. 285–291. https://doi.org/10.1007/978-3-319-04705-8_33 Companys, R., Mateo, M., 2007. Different behaviour of a double branch-and-bound algorithm on Fm|prmu|Cmax and Fm|block|Cmax problems. Comput. Oper. Res. 34, 938–953. Companys, R., Ribas, I., Mateo, M., 2010. Improvement tools for NEH based heuristics on permutation and blocking flow shop scheduling problems, IFIP Advances in Information and Communication Technology. https://doi.org/10.1007/978-3-642-16358-6_5

CR IP T

Fernandez-Viagas, V., Framinan, J.M., 2014. A bounded-search iterated greedy algorithm for the distributed permutation flowshop scheduling problem. Int. J. Prod. Res. 53, 1111– 1123. https://doi.org/10.1080/00207543.2014.948578 Gao, J., Chen, R., 2012. A hybrid genetic algorithm for the distributed permutation flowshop scheduling problem. Gao, J., Chen, R., Deng, W., 2013. An efficient tabu search algorithm for the distributed permutation flowshop scheduling problem. Int. J. Prod. Res. 51, 641–651. https://doi.org/10.1080/00207543.2011.644819

AN US

Gao, J., Chen, R., Deng, W., Liu, Y., 2012a. Solving multi-factory flowshop problems with a novel variable neighbourhood descent algorithm. J. Comput. Inf. Syst. 8, 2025–2032.

Gao, J., Chen, R., Liu, Y., 2012b. A Knowledge-based Genetic Algorithm for Permutation Flowshop Scheduling Problems with Multiple Factories. Int. J. Adv. Comput. Technol. 4, 121–129. https://doi.org/10.4156/ijact.vol4.issue7.13

M

Gong, H., Tang, L., Duin, C.W., 2010. A two-stage flow shop scheduling problem on a batching machine and a discrete machine with blocking and shared setup times. Disrupt. Manag. 37, 960–969. https://doi.org/DOI: 10.1016/j.cor.2009.08.001

ED

Han, Y.-Y., Liang, J.J., Pan, Q.-K., Li, J.-Q., Sang, H.-Y., Cao, N.N., 2012. Effective hybrid discrete artificial bee colony algorithms for the total flowtime minimization in the blocking flowshop problem. Int. J. Adv. Manuf. Technol. 67, 397–414. https://doi.org/10.1007/s00170-012-4493-5

PT

He, D.W., Kusiak, A., Artiba, A., 1996. A scheduling problem in glass manufacturing. IIE Trans. 28, 129–139. https://doi.org/10.1080/07408179608966258

CE

Liang, J.J., Pan, Q.-K., Tiejun, C., Wang, L., 2011. Solving the blocking flow shop scheduling problem by a dynamic multi-swarm particle swarm optimizer. Int. J. Adv. Manuf. Technol. 55, 755–762. https://doi.org/10.1007/s00170-010-3111-7

AC

Lin, S.-W., Ying, K.-C., Huang, C.-Y., 2013. Minimising makespan in distributed permutation flowshops using a modified iterated greedy algorithm. Int. J. Prod. Res. 51, 5029–5038. https://doi.org/10.1080/00207543.2013.790571 Liu, H., Gao, L., 2010. A discrete electromagnetism-like mechanism algorithm for solving distributed permutation flowshop scheduling problem. Proc. - 2010 Int. Conf. Manuf. Autom. ICMA 2010 156–163. https://doi.org/10.1109/ICMA.2010.17 Martinez, S., Dauzère-Pérès, S., Guéret, C., Mati, Y., Sauer, N., 2006. Complexity of flowshop scheduling problems with a new blocking constraint. Eur. J. Oper. Res. 169, 855–864. https://doi.org/DOI: 10.1016/j.ejor.2004.08.046 Montgomery, D.C., 2013. Design and Analysis of Experiments. Design. https://doi.org/10.1198/tech.2006.s372 30

ACCEPTED MANUSCRIPT Naderi, B., Ruiz, R., 2014. A scatter search algorithm for the distributed permutation flowshop scheduling problem. Eur. J. Oper. Res. 239, 323–334. https://doi.org/10.1016/j.ejor.2014.05.024 Naderi, B., Ruiz, R., 2010. The distributed permutation flowshop scheduling problem. Comput. Oper. Res. 37, 754–768. https://doi.org/10.1016/j.cor.2009.06.019 Nagano, M.S., Komesu, A.S., Miyata, H.H., 2017. An evolutionary clustering search for the total tardiness blocking flow shop problem. J. Intell. Manuf. https://doi.org/10.1007/s10845-017-1358-7

CR IP T

Pan, Q.-K., Wang, L., 2012. Effective heuristics for the blocking flowshop scheduling problem with makespan minimization. Omega 40, 218–229. https://doi.org/DOI: 10.1016/j.omega.2011.06.002 Pan, Q., Wang, L., Sang, H., Li, J., Liu, M., 2013. A High Performing Memetic Algorithm for the Flowshop Scheduling Problem With Blocking. IEEE Trans. Autom. Sci. Eng. 10, 741– 756. https://doi.org/10.1109/TASE.2012.2219860

AN US

Potts, C.N., Van Wassenhove, L.N., 1982. A decomposition algorithm for the single machine total tardiness problem. Oper. Res. Lett. 1, 177–181. Ribas, I., Companys, R., 2015. Efficient heuristic algorithms for the blocking flow shop scheduling problem with total flow time minimization. Comput. Ind. Eng. 87. https://doi.org/10.1016/j.cie.2015.04.013

M

Ribas, I., Companys, R., Tort-Martorell, X., 2017. Efficient heuristics for the parallel blocking flow shop scheduling problem. Expert Syst. Appl. 74, 41–54. https://doi.org/10.1016/J.ESWA.2017.01.006

ED

Ribas, I., Companys, R., Tort-Martorell, X., 2015. An efficient Discrete Artificial Bee Colony algorithm for the blocking flow shop problem with total flowtime minimization. Expert Syst. Appl. 42, 6155–6167. https://doi.org/10.1016/j.eswa.2015.03.026 Ribas, I., Companys, R., Tort-Martorell, X., 2013. A competitive variable neighbourhood search algorithm for the blocking flowshop problem. Eur. J. Ind. Eng. 7, 729–754.

CE

PT

Ribas, I., Companys, R., Tort-Martorell, X., 2011. An iterated greedy algorithm for the flowshop scheduling problem with blocking. Omega 39, 293–301. https://doi.org/10.1016/j.omega.2010.07.007 Ronconi, D.P., 2005. A Branch-and-Bound Algorithm to Minimize the Makespan in a Flowshop with Blocking. Ann. Oper. Res. 138, 53–65. https://doi.org/10.1007/s10479-005-2444-3

AC

Ronconi, D.P., 2004. A note on constructive heuristics for the flowshop problem with blocking. Int. J. Prod. Econ. 87, 39–48. Ronconi, D.P., Armentano, V.A., 2001. Lower bounding schemes for flowshops with blocking in-process. J. Oper. Res. Soc. 52, 1289–1297. https://doi.org/10.1057/palgrave.jors.2601220 Ronconi, D.P., Henriques, L.R.S., 2009. Some heuristic algorithms for total tardiness minimization in a flowshop with blocking. Omega 37, 272–281. https://doi.org/DOI: 10.1016/j.omega.2007.01.003 Ruiz, R., Pan, Q.K., Naderi, B., 2018. Iterated Greedy methods for the distributed permutation flowshop scheduling problem. Omega, In press. ttps://doi.org/10.1016/j.omega.2018.03.004 31

ACCEPTED MANUSCRIPT Sethi, S.P., Sriskandarajah, C., Sorger, G., Blazewicz, J., Kubiak, W., 1992. Sequencing of parts and robot moves in a robotic cell. Int. J. Flex. Manuf. Syst. 4, 331–358. https://doi.org/DOI: 10.1007/BF01324886 Taillard, E., 1993. Benchmarks for basic scheduling problems. Eur. J. Oper. Res. 64, 278–285. Takano, M.I., Nagano, M.S., 2019. Evaluating the performance of constructive heuhristics for the blocking flow shop scheduling problem with setup times. Int. J. Ind. Eng. Comput. 37– 50. https://doi.org/10.5267/j.ijiec.2018.5.002 Takano, M.I., Nagano, M.S., 2017. A branch-and-bound method to minimize the makespan in a permutation flow shop with blocking and setup times. Cogent Engineering 4: 1389638,1-16. https://doi.org/10.1080/23311916.2017.1389638

CR IP T



Tavakkoli-Moghaddam, R., Safaei, N., Sassani, F., 2009. A memetic algorithm for the flexible flow line scheduling problem with processor blocking. Comput. Oper. Res. 36, 402–414. Vairaktarakis, G., Elhafsi, M., 2000. The use of flowlines to simplify routing complexity in twostage flowshops. IIE Trans. 32, 687–699. https://doi.org/10.1023/A:1007652626510

AN US

Wang, C., Song, S., Gupta, J.N.D., Wu, C., 2012. A three-phase algorithm for flowshop scheduling with blocking to minimize makespan. Comput. Oper. Res. 39, 2880–2887. https://doi.org/10.1016/j.cor.2012.02.020 Wang, L., Pan, Q.-K., Suganthan, P.N., Wang, W.-H., Wang, Y.-M., 2010. A novel hybrid discrete differential evolution algorithm for blocking flow shop scheduling problems. Hybrid Metaheuristics 37, 509–520. https://doi.org/DOI: 10.1016/j.cor.2008.12.004

M

Wang, L., Pan, Q.-K., Tasgetiren, M.F., 2011. A hybrid harmony search algorithm for the blocking permutation flow shop scheduling problem. Comput. Ind. Eng. https://doi.org/10.1016/j.cie.2011.02.013

ED

Wang, S., Wang, L., Liu, M., Xu, Y., 2013. An effective estimation of distribution algorithm for solving the distributed permutation flow-shop scheduling problem. Int. J. Prod. Econ. 145, 387–396. https://doi.org/10.1016/j.ijpe.2013.05.004

CE

PT

Xu, Y., Wang, L., Wang, S., Liu, M., 2013. An effective hybrid immune algorithm for solving the distributed permutation flow-shop scheduling problem. Eng. Optim. 46, 1269–1283. https://doi.org/10.1080/0305215X.2013.827673

AC

Zhang, X., Van De Velde, S., 2012. Approximation algorithms for the parallel flow shop problem. Eur. J. Oper. Res. https://doi.org/10.1016/j.ejor.2011.08.007

32