Permutation flowshops in group scheduling with ... - Semantic Scholar

European J. Industrial Engineering, Vol. 6, No. 2, 2012

Permutation flowshops in group scheduling with sequence-dependent setup times B. Naderi Young Researchers Club, Qazvin Branch, Islamic Azad University, Qazvin, Iran E-mail: [email protected]

Nasser Salmasi* Department of Industrial Engineering, Sharif University of Technology, P.O. Box 11365-8639, Tehran, Iran Fax: (+9821)-66022702 E-mail: [email protected] *Corresponding author Abstract: This paper focuses on the flow shop sequence dependent group scheduling (FSDGS) problem with minimisation of total completion time as the criterion (Fm|fmls, prmu, Splk|∑Cj). The research problem is formulated in form of two different mixed integer linear programming (MILP) models. Comparing with the latest MILP model for the proposed problem in the literature, the complexity size of the proposed models are significantly reduced. One of the proposed mathematical models is so effective that even medium-sized instances (problems up to 60 jobs in all groups) are solved to optimality in a reasonable amount of time. Moreover, a metaheuristic hybridising genetic and simulated annealing algorithm, called GSA, is proposed to solve the problems heuristically. All the results and analyses show the high performance of the proposed mathematical models as well as the proposed metaheuristic algorithm compared to the available ones in literature. [Received 11 February 2010; Revised 24 May 2010; Accepted 30 September 2010] Keywords: sequence-dependent setup time scheduling; mixed integer linear programming; MILP; hybrid metaheuristics; group scheduling; permutation flowshop. Reference to this paper should be made as follows: Naderi, B. and Salmasi, N. (2012) ‘Permutation flowshops in group scheduling with sequence-dependent setup times’, European J. Industrial Engineering, Vol. 6, No. 2, pp.177–198. Biographical notes: B. Naderi has recently received his PhD in the area of Industrial Engineering from Amir Kabir University of Technology (September 2010). His areas of research interests are the applications of operations research in production scheduling, supply chain management, and portfolio management. More specifically, he is working on mathematical modelling and solution methods.

Copyright © 2012 Inderscience Enterprises Ltd.

177

178

B. Naderi and N. Salmasi Nasser Salmasi is currently an Assistant Professor at the Department of Industrial Engineering at Sharif University of Technology (SUT), Tehran, Iran. He received his PhD from the Department of Industrial and Manufacturing Engineering at Oregon State University (OSU). His primary areas of research interests are applied operations research, scheduling, and information systems.

1

Introduction

In this research, flow shop sequence dependent group scheduling (FSDGS) problem with minimisation of total completion time as the criterion (Fm|fmls, prmu, Splk|∑Cj) is considered. Assume that g sets (groups) of jobs are assigned to a flow shop cell with m machines to be processed. Each group includes nk jobs (k =1, 2, …, g). A major setup is required to start processing jobs belong to each group on each machine. The setup time of each group depends on the previous processed group on each machine. All jobs should be processed in the same order on all machines. FSDGS problem has been studied for the first time by Schaller et al. (2000). They developed several heuristic algorithms to solve the problem heuristically with minimisation of makespan as the criterion. Franca et al. (2005) developed several metaheuristic algorithms based on memetic algorithm (MA) and genetic algorithm (GA) to solve the proposed problem with makespan criterion. They showed that their MA has a superior performance compared to GA as well as the heuristics earlier developed by Schaller et al. (2000). Logendran et al. (2006) developed several metaheuristic algorithms based on tabu search (TS) for the two-machine FSDGS problem with minimisation of makespan as the criterion. They also developed several lower bounding techniques to evaluate the performance of the proposed metaheuristics. They showed that one of their lower bounding techniques which was developed centred on the mathematical model of the problem has a superior performance than the other proposed lower bounding techniques in that research. Hendizadeh et al. (2008) developed several metaheuristic algorithms based on TS for multi-stage problems and showed that the performance of their algorithm based on the value of the objective function (i.e., minimisation of makespan) is the same as MA proposed by Franca et al. (2005). Salmasi et al. (2010) investigated the FSDGS problem with minimisation of total completion time as the criterion, for the first time. In this research they developed a mixed integer linear programming model for FSDGS problems for the first time. They also developed several metaheuristic algorithms based on TS as well as a metaheuristic algorithm based on ant colony optimisation (ACO) to solve the problem. They performed a detailed experimental design based on split plot design to show that the developed TS algorithm with intensifying long term memory has a better performance compared to the other developed TS algorithms. Then, they compared the results of the best TS algorithm with the results of the ACO algorithm as a paired-t-test comparison. The results of the experiment showed that the ACO algorithm has a superior performance compared to the TS algorithm. Furthermore, they developed a lower bounding method based on branch-and-price algorithm by applying the mathematical model to estimate the quality of the developed upper bounds. By using this method, they provided lower bounds for the test problems. The result of the comparison was not very satisfactory since the mathematical model could not be solved optimally and thus, the lower bounding method

Permutation flowshops in group scheduling

179

could not provide high quality results. This is our main motivation to develop an efficient mathematical model for the proposed research problem. Mathematical models are regarded as one of the best methods to characterise a scheduling problem. Moreover, due to the progress obtained in computer capacities and advent of specified software, it is not a surprise to indicate that they are recently treated as a solution strategy to tackle the small even medium sized problems. In this paper, two high performing mathematical models, in the form of mixed integer linear programming (MILP), are developed to formulate the problem configuration. To evaluate the proposed models, their specifications, such as the number of binary and continuous variables as well as the number of constraints are compared with the mathematical model proposed by Salmasi et al. (2010). The capability of all three models (two proposed in this paper and the one proposed by Salmasi et al. (2010) are assessed to solve a set of instances existing in the literature. Salmasi et al. (2010) showed that the FSDGS problem with minimisation of total completion time (TCT) as the criterion is an NP-hard problem. Thus, the exact methods are handicapped to solve the large-sized problems. In this regard, the FSDGS problem is dealt with by a metaheuristic approach, a hybridisation of genetic and simulated annealing algorithms (GSA). Both GA and SA seem to be sensitive towards great choice of their operators and parameters. Therefore, the performance of several operators and parameters in the proposed hybrid metaheuristics is examined. This parameter tuning is conducted by means of Taguchi method as a robust experimental design. After well calibrating GSA, its general performance against the optimal solution obtained by MILP models and an existing ACO which was proposed by Salmasi et al. (2010) is evaluated. The rest of the paper is organised as follows. Section 2 introduces the proposed MILP models. Section 3 presents the proposed hybrid metaheuristic. Section 4 evaluates the proposed models and algorithm. Finally, Section 5 gives the conclusions and interesting future research directions.

2

Mathematical models

For a long time, mathematical models were just utilised to specify all the characteristics of a problem. If the problem is NP-hard, solving the mathematical model is not regarded as an effective solution approach. Besides, mathematical models play the role of a key starting point in other solution strategies such as branch-and-price algorithm. More effective a model is constructed, more capable these strategies become. Considering all these together, researchers put more attempts on designing mathematical models to contribute a basis for future studies (Tseng et al., 2004; Stafford et al., 2005; Tseng and Stafford, 2008). Scheduling problems are commonly formulated through MILP models. Despite MILP models’ importance, in the literature for FSDGS, only one MILP model is proposed by Salmasi et al. (2010). By analysing the strength and weakness of Salmasi et al. (2010) mathematical model, in this research, two new mathematical models for the research problem are proposed. In these proposed models by reducing the size complexity of MILP model, the required computational time to optimally solve small- and medium-sized problems is remarkably decreased compared to the mathematical model proposed by Salmasi et al. (2010). The proposed two mathematical models are described

180

B. Naderi and N. Salmasi

in the following sections. The notations and parameters used in both models are as follows: no

the number of jobs

m

the number of machines

g

the number of groups

l, j

index for jobs, l, j ∈ {1, 2,…, no }

i

index for machines, i ∈ {1, 2,…, m}

k, t

index for groups, k, t ∈ {1, 2,… , g}

Pj,i

the run time of job j on machine i

at,k,i the setup time of group k if it is processed immediately after group t on machine i Gk

a set including the jobs belonging to group k ∈ {1, 2,…, g}

nk

the number of jobs in group k; that is nk = Gk

M

a large positive number.

2.1 Model 1 In order to evaluate the value of the objective function of a sequence for a FSDGS problem, the sequence of groups as well as the sequence of the jobs in each group should be determined. In Model 1, two different types of sequence-based variables are proposed. Since the setup times between groups are sequence dependent, the binary variables related to group sequence need to keep track of every two consecutive groups. Therefore, the immediate sequence-based variable is applied. In this case, it is assumed that a dummy group, called reference group or group zero, that is processed as the first group in the sequence is needed to be defined. This assumption is basically used to calculate the required setup time for the first group on each machine. On the other hand, there is no need to directly clarify the exact position of the jobs inside each group. Armed with this, the model exploits another type of binary variables that merely express whether a job precedes/succeeds another job. Therefore, it entails less variables in comparison with immediate sequence-based one. The following decision variables are used in this model. Decision variables:

X l, j

⎧1, if job j is processed after job l ⎪ =⎨ where ∀k , j ,l∈Gk , j >l ⎪ ⎩0, otherwise

⎧1, if group k is processed immediately after group t ⎪ U t ,k = ⎨ where t ≠ k ⎪0, otherwise ⎩

(1)

(2)


181

Fk ,i the finishing time of the last job of group k on machine i

(3)

Sk ,i the starting time of the first job of group k on machine i

(4)

C j ,i the completion time of job j on machine i

(5)

Model 1 formulates FSDGS problems with minimisation of TCT as the criterion as follows: Minimise

∑

no j =1

(6)

C j ,m

Subject to: C j ,i ≥ C j ,i −1 + p j ,i

∀ j ,i

(7)

C j ,i ≥ Cl ,i + p j ,i − 1 − X l , j × M

∀k , j ,l∈Gk ,i ,l < j

(8)

Cl ,i ≥ C j ,i + pl ,i − X l , j × M

∀k , j ,l∈Gk ,i ,l < j

(9)

C j ,i ≥ S k ,i + p j ,i

∀k , j∈Gk , i

(10)

Fk ,i ≥ C j ,i

∀k , j∈Gk , i

(11)

Sk ,i ≥ Ft ,i + at ,k ,i − (1 − U t ,k ) × M

∀t∈{0,1,…, g},k ≠t ,i

(12)

∑

g t = 0, t ≠ k

U t ,k = 1

∀k

(13)

∑

g k =1, t ≠ k

U t ,k ≤ 1

∀t

(14)

∑

g k =1

(

U o,k = 1

)

(15)

U t , k + U k ,t ≤ 1

(16)

C j ,i , S k ,i , Fk ,i ≥ 0

(17)

X l , j , U t , k ∈ {0,1}

(18)

where C j ,0 = F0,i = 1. As it is mentioned the objective function is based on calculation of the TCT. The objective function is presented in equation (6). It is clear that a job cannot be processed on two different machines at the same time. Thus, the completion time of a job on a machine should be greater than its completion time on the previous machine plus its runtime on the machine. Constraint set (7) is incorporated to the model for this reason. A machine cannot also process two different jobs simultaneously. In other words, the difference between the completion times of two different jobs on a machine should be greater than the runtime of the job is processed later. Constraint sets (8) and (9) are incorporated to the model in order to support this fact. Constraint sets (10) and (11) are incorporated to the model to ensure that every job must start and finish between starting

182


and finishing times of the group it belongs to. The process of each group on each machine can start after the process of the last job of the immediately previous group is finished. Constraint set (12) is incorporated to the model for this reason. The start time of processing a job belonging to a certain group (say group k) on a machine (say machine i) is greater than the finish time of the last group (say group t) process on the machine plus the setup of the group k (i.e., at,k,i). Each group should have exactly one group as the previous processed group. Constraint set (13) is incorporated to the model to support this fact. Constraint set (14) is incorporated to the model to make sure that each group has at most one successor group. Constraint set (15) is incorporated to the model in order to make sure that the reference group is assigned to the first slot. Each group (say group t) can be processed either before another group (say group k) or after that group. Constraint set (16) is incorporated to the model to support this fact. Constraint sets (17) and (18) define the decision variables.

2.2 Model 2 In this model, the position-based concept is used to find job sequence of each group. In order to find the sequence of groups, the immediate sequence-based is still used as the best possible choice. The idea of considering the reference group still exists in Model 2. Another alternative for group sequence decision could be the consideration of position-based variables. In this case, numerous additional binary variables should be used to find relative order of groups, as it is used in Salmasi et al. (2010) model. Due to huge number of required binary variables, it is very likely that this possible alternative appears to be ineffective. In this model the domain of indices and are changed as follows: j

job index

{1, 2,… , nk }

l

slot index

{1, 2,… , nk }

In addition to decision variables defined in (2), (3), and (4), this model consists of the following decision variables:

X j ,l

⎧1, if job j occupies the lth position of group k ⎪ =⎨ where ∀k , j∈Gk ,i ={1,…, nk } ⎪ ⎩0, otherwise

Ck ,l ,i Completion time of the job in the lth position of group k on machine i.

(19)

(20)

A FSDGS problem is characterised by Model 2 as the followings: Minimise

∑

g k =1

∑

nk l =1

Ck ,l , m

(21)

Subject to: Constraints (12), (13), (14), (15) and (16)

∑

nk l =1

X j ,l = 1

∀k , j∈Gk

(22)


∑

j∈Gk

X j ,l = 1

Ck ,l ,i ≥ Ck ,l −1,i +

∀k ,l∈{1,…,nk } j ,l

× p j ,i

∀k ,l∈{2,…, nk },i

(24)

∑X

j ,l

× p j ,i

∀k ,l∈{1,…, nk },i

(25)

∀k ,i

(26)

∀ k ,i

(27)

j∈Gk

Ck ,1,i ≥ Sk ,i +

(23)

∑X

j∈Gk

Ck ,l ,i ≥ Ck ,l ,i −1 +

183

∑X

j ,i

× p j ,i

j∈Gk

Fk ,i ≥ Ck , nk ,i Ck ,l ,i , Sk ,i Fk ,i ≥ 0

(28)

X j ,l , Uτ ,k ∈ {0, 1}

(29)

where Ck ,l ,0 = F0,i = 0. TCT value as the objective function is obtained by equation (21). Constraint set (22) specifies that every job occupies exactly one slot in the group it belongs to. Constraint set (23) is incorporated to the model to make sure that each slot of any groups is assigned to one of the jobs in that group. The process of a job on a machine can start if the process of the job in the previous slot is finished. Constraint set (24) is incorporated to the model for this reason. Constraint set (25) is incorporated to the model to ensure that the process of each job starts after the process of the job is concluded on the previous machine. Constraint set (26) is incorporated to the model to ensure that the process of the job in the first slot of each group starts after the starting point of the group. It is clear that the completion time of a group is greater than the completion time of the job in the last slot. Constraint set (27) is incorporated to the model for this reason.

2.3 Comparison of the models In this section by comparing the number of variables, the number of constraints and status of the big M inside constraints, the proposed models are compared with the one proposed by Salmasi et al. (2010). Table 1 summarises the number of binary variables (NBV), the number of continuous variables (NCV), and the number of constraints (NC) required by each model to formulate a same problem size with a total no of jobs and m machines. There are g groups with n1 jobs in the first group, n2 jobs in the second group until ng jobs in the gth group

(∑ (b

g k =1

max

)

nk = no . The maximum number of jobs in a group is considered as bmax

{

= max n1 , n2 ,… , ng

mathematical model.

}) .

This parameter is used in Salmasi et al. (2010) proposed

184


Table 1

The comparison of the models

Model

NBV

Salmasi et al. (2010)

g 2 ( g + 1) + gbmax ( bmax − 1)

∑

Model 1

g k =1

NCV

NC

m ( 2 g + 1) + gmbmax

2 g 2 ( g − 1) + gmbmax ( bmax + 2 )

nk ( nk − 1) 2

m ( 2 g + no )

n 2k + g 2

m ( 2 g + no )

− g ( bmax + m + 4 )

m⎡ ⎣

nk ( nk − 1) + 3n0 ⎤ ⎦ 1⎞ 3 2⎛ +g ⎜ m + ⎟ + g +1 2⎠ 2 ⎝

+g2

∑

Model 2

g k =1

∑

g

k =1

2n0 ( m + 1) 1 3⎞ ⎛ + g ⎜ gm + m + g + ⎟ + 1 2 2⎠ ⎝

2.3.1 The models’ comparison based on complexity size in binary variables Considering the NBV, the complexity size of Salmasi et al. (2010) model is calculated

(

)

based on g 3 + gb 2max . Therefore, the NBVs is cubic in g and quadratic in bmax. The

(∑

)

g NBVs of the Model 1 has the complexity size of 0 ⎡⎢ g 2 + n 2k ⎤⎥ which means k =1 ⎣ ⎦ quadratic in both g and nk. Interestingly, the complexity reduction is obtained in both dimensions (i.e., the number of groups and the number of jobs in each group) in Model 1 compared to Salmasi et al. (2010). The NBV of Salmasi et al. (2010) model soars cubically in the number of groups while for the first model it increases in quadratic manner (i.e., one degree complexity reduction). The most important weakness of Salmasi et al. (2010) model is that by increasing the maximum number of jobs in a group, the problem becomes more complex. That makes the model unable to solve any problem optimally if the maximum number of jobs in a group be greater than five. This could be a deadly drawback when the numbers of jobs in groups are not monotonous. This shortcoming is smoothed by Model 1. The NBVs of Model 2 is calculated by g 0 ⎡⎢ g 2 + n 2k ⎤⎥ . It is clear that the complexity of Model 2 is the same as the k =1 ⎣ ⎦ complexity of Model 1. We should notice that, in nk dimension, the NBV of Model 1 increases with half slope of the NBV of Model 2, although they have the same complexity size.

(∑

Table 2

)

NBV comparison of the models Problem size

Example

(

g n1 ,…, ng

The models

)

n0

m

1

20

5

3 (6, 9, 5)

2

50

5

4 (10, 15, 15, 10)

3

100

10

4

200

10


Model 1

Model 2

252

70

142

950

300

666

5 (20, 30, 10, 15, 25 )

4,500

1,075

2,275

6 (25, 40, 15, 50, 30, 40)

14,952

6,625

7,486


185

Table 2 includes several numerical examples about the number of BVs of the proposed models for different problem sizes. For instance, in a problem with no = 20, m = 5, g = 3, n1 = 6, n2 = 9 and n3 = 5, Salmasi et al. (2010) model has 252 BVs, while Models 1 and 2 consist of 70 and 142 BVs, respectively. These differences become more significant when the problem size grows up.

2.3.2 The models’ comparison based on the complexity size in continuous variables As noted in Table 1, the NCV of Salmasi et al. (2010) model is calculated based on 0(mgbmax), while Models 1 and 2 have the same complexity size of 0(gm + n0m). Depending on the values of m, g, bmax and n, each model could have less complexity size. Generally, since the complexity orders of the models in the above parameters are linear, overly the NCV is not an influential factor. Table 3 provides the NCV for each formulation for different problem sizes defined in Table 1. It is clear that in realistic sizes, Models 1 and 2 consist of less NCV than Salmasi et al. (2010) model has. Table 3

NCV comparison of the models Problem size

Example

(

The models

)


Model 1

Model 2

3 (6, 9, 5)

170

130

130

4 (10, 15, 15, 10)

345

290

290

10

5 (20, 30, 10, 15, 25 )

1,610

1,100

1,100

10

6 (25, 40, 15, 50, 30, 40)

3,130

2,120

2,120

n0

m

g n1 ,…, ng

1

20

5

2

50

5

3

100

4

200

2.3.3 The model’s comparison based on the complexity size in constraints The complexity size of Salmasi et al. (2010) model in NC is calculated based on

(

)

0 g 3 + gmb 2max . Therefore, its number of NCs is cubic in g and quadratic in bmax.

(∑

)

g Model 1 has 0 ⎡⎢ g 2 m + mn0 + m n 2k ⎤⎥ in NCs. Comparing these two models, the k =1 ⎣ ⎦ complexity is reduced to quadratic in both g and nk. In Model 1, the parameters n0 and m are linearly influencing the NC. The model developed by Salmasi et al. (2010) suffers from presence of parameter bmax that might make the model totally ineffective. The

(

)

complexity size of Model 2 in NCs outstandingly is calculated based on 0 g 2 m + mn0 . In this model, the required NC is remarkably reduced. Its complexity is linear in both m and n0, and quadratic in g. Table 4 shows the numerically comparison of the models based on NC. Model 2 impressively has fewer variables than the two other models. In the largest size (i.e., n0 = 200, m = 10 and g = 6), the NC goes up to only 4,908 constraints while the Model 1 and Salmasi et al. (2010) model consist of 78,888 and 155,976 constraints, respectively.

186


Table 4

NC comparison of the models Problem size

Example 1

The models

(

n0

m

g n1 ,…, ng

20

5

3 (6, 9, 5)

)


Model 1

Model 2

1,467

965

445

2

50

5

4 (10, 15, 15, 10)

5,100

3,845

735

3

100

10

5 (20, 30, 10, 15, 25 )

47,980

24,796

2,571

4

200

10

6 (25, 40, 15, 50, 30, 40)

155,976

78,888

4,908

2.3.4 Status of big M In MILP models, big Ms are used to linearly relax dichotomous constraints. It has been concluded that the presence of big M would result in poor model performance by increasing the size of the problem. With this in mind, it is also interesting to compare the status of big M in the proposed models. As earlier mentioned, in FSDGS problems, there are two decision dimensions, i.e., finding the sequence of groups and finding the sequence of jobs belonging to each group. Consequently, a fraction of constraints in each model are to follow each decision. Now, we consider each model to see whether there are dichotomous constraints for each decision or not. The results are presented in Table 5. Group sequencing

the concept, used in Models 1 and 2, consist of dichotomous constraints. In Salmasi et al. (2010) model, the authors developed the model based on another concept in which there is no dichotomous constraint, but it suffers from excessive NBV that make it completely ineffective. Therefore, both concepts are sensitive towards an increase in the number of groups i.e., g.

Job sequencing

Salmasi et al. (2010) model and Model 1 have dichotomous constraints for this decision whereas Model 2 is released from presence of them. Considering very less NCs along with low NBVs, Model 2 is expected to be very effective and could solve problems with larger number of jobs.

Table 5

Comparison of the models via big M

The models

Big M status Job sequencing

Group sequencing


Yes

No

Model 1

Yes

Yes

Model 2

No

Yes

3

Solution methods

As mentioned, the FSDGS has two decision dimensions: group sequencing and job sequencing inside each group. In this research, a hybrid metaheuristic algorithm, called GSA, is utilised including two algorithms GA, and simulated annealing (SA). Among the


187

successful applications of GA and SA in scheduling problems, one can point out to papers by Damodaran et al. (2009) and Kim et al. (2002), respectively. In the proposed GSA, the GA searches for better group sequences while the SA is used to search for better job sequence inside each group for a given group sequence. GSA starts from a population of encoded solutions (popsize). Then, at each iteration of GSA, the population evolves by a set of operators so long as a stopping criterion is met. The GSA procedure could be seen in three clusters: Steps linked to: 1

population evaluation

2

group sequencing evolution

3

job sequencing evolution.

Each of these clusters consists of the steps shown by Table 6. The used mechanisms are detailed in the following subsections. Table 6

The steps of the clusters of GSA ## Steps linked to population evaluation ##

Step 1

Initialising the population by generating initial solutions.

Step 2

Assess each solution in current population based on its TCT and assign a value called fitness.

Step 3

If the stopping criterion is not met, go to Steps 4, …, 11; otherwise, stop.

Step 4

Control the trend of improvement by an operator, called restart phase. ## Steps linked to group sequencing evolution ##

Step 5 Step 6

Directly copy two of the best solutions to the next population (elite operator) and make i = 1. If i ≤

popsize − 2 , then go to the Step 7; otherwise, go to the Step 10. 2

Step 7

Generate two new solutions in this way: a selection mechanism selects two solutions of the current population so as to give higher chance of being selected to solutions with better fitness values. The group sequences of selected solutions (parents) are combined and produce two new offspring by an operator called crossover.

Step 8

Perturb each offspring by another mechanism called mutation. This mechanism still focuses on the group sequence.

Step 9

Increase i by one unit and go to Step 6. ## Steps linked to job sequencing evolution ##

Step 10

Apply a local search engine to each produced offspring for a better job sequence. This engine is based on SA.

Step 11

Construct the new population by the offspring, and then go to Step 2.

3.1 Solution encoding, initialisation, fitness evaluation and selection mechanism The encoding scheme is used to represent a solution to metaheuristics. The widely used encoding scheme in the literature of the permutation flow shop is the permutation list. The permutation list merely expresses the relative order of all objects in which they are

188


processed. In FSDGS, a solution can be encoded by g + 1 permutation lists, one permutation list to represent group sequence and g permutation lists to represent job sequence inside per group. The initial solutions are randomly generated from feasible space. Since the higher values of fitness are desirable and the objective is minimisation of TCT, the fitness value of each solution equals to 1/TCT. There are two commonly used selection mechanisms, Roulette wheel and tournament selection (Goldberg, 1989). In Roulette wheel selection, each solution, namely solution π, is given a normalised probability obtained by the following formula: 1 TCT (π )

∑

popsize y =1

1 TCT ( y )

(30)

In tournament selection, two solutions are randomly picked and the better is selected.

3.2 Restart phase In order to prevent to be trapped in a local optimal solution, an extra operator, the restart phase, is used in the algorithm (Stafford et al., 2005). The procedure is as follows: If no improvement is obtained in one generation, a counter increases by one unit. Clearly, when any improvement is found, the counter restarts from 0. If counter shows a number more than a specified number (MAX), 50% of worse solutions in the population are replaced by randomly regenerated solutions.

3.3 Elite, crossover and mutation operators Two of the best solutions in current generation are copied into the next generation by elite operator. By crossover operators, the remaining solutions of the next generation are produced. In fact, the crossover operator combines two selected solutions to build two new solutions. Among the possible options, we intend to examine the performances of three crossovers shown high performances in scheduling problems. These crossovers are applicable to the case of permutation list encoding scheme. The first crossover is ‘similar block order crossover’ or SBOX. This crossover has been shown to be very high performing in permutation flowshops against many other crossover operators (Ruiz et al., 2006). The second crossover is ‘uniform parameterised crossover’ or UPC used by Kurz and Askin (2004) for flexible flow line scheduling. It is necessary to state that UPC works through random keys (RK). In other words, g random numbers ranging (0, 1) are generated. Then, each RK based on ascending order is assigned to a single group in the permutation from left to right. Then, for each group a random number is generated. If the value is less than 0.7, then the RK of the corresponding group from the first parent is copied to the first child; otherwise, the RK from the second parent is selected. Finally, groups are sorted according to the ascending order of RKs. The same procedure is applied to produce the second child. The third crossover is in the literature of the open shop (Senthilkumar and Shahabudeen, 2006). This crossover is a “modification of two-point crossover” or MTPC so as to consider the characteristics of permutation list. The performances of three above-mentioned crossovers in calibration section are compared.


189

Mutation operators are to further change new crossed solutions. There are three different types of mutation operators widely used in the literature. The first one is shift mutation by which the position of one randomly selected group is randomly changed. The second is swap mutation by which the positions of two randomly selected groups are swapped. The third is inversion mutation by which the positions of groups between two randomly selected cut points are reversed. The probability of mutation (Pm) is defined as such: Pm is the probability by which each crossed solution undergoes the mutation operator. The type of mutation and Pm value is chosen in the calibration section.

3.4 Local search engine (SA) Since group sequencing decision of GPFS is improved by the GA’s operator, a local search engine is required to search for better job sequence. In this case, a simple and fast, yet effective SA, is applied. The applied SA has the following characteristics: it starts from an initial solution, each solution of current generation, and makes a series of moves until a stopping criterion is met. The basic idea of the proposed SA is to generate a new solution, namely θ, by shift operator from the neighbourhood of the current solution, namely π. This new solution is accepted or rejected by another random rule. A parameter t, called the temperature, controls the acceptance rule. The difference between the objective function values of two candidate solutions are computed based on ΔC = TCT (θ ) − TCT (π ) . If ΔC ≤ 0, solution θ is accepted. Otherwise, solution θ is accepted with probability equal to exp ( ΔC / t ) . That is, a random number in the range of (0, 1) is generated, then if it is less than exp ( ΔC / t ) , solution θ is accepted. Table 7

The steps of the proposed local search engine

Step 1

Take the solution π, and k = l.

Step 2

If k ≤ g , then t = T0 initialise counter = 1 and go to Step 3; otherwise, stop.

Step 3

If the counter ≤ CT, then make i = 1 and go to Steps 4, …, 9; otherwise, increase k by one unit and go to Step 2.

Step 4

Generate a new neighbour θ from current solution π by concentrating on job sequence of k-th group and calculate ∆C = TCT (θ) – TCT (π).

Step 5

If ∆C ≤ 0, then replace solution π with solution θ and go to Step 7; otherwise, go to Step 6.

Step 6

Generate a random number between (0, 1). If the number is less than exp (∆C/t) then replace solution π with solution θ and go to Step 7.

Step 7

Increase i by one unit. If i ≤ FN, then go to Step 4; otherwise, go to Step 8.

Step 8

If the best solution is improved, then restart the counter from zero; otherwise, increase the counter by one unit.

Step 9

Decrease the temperature by multiplying it to 0.95 (i.e., t = 0.95t) and go to Step 3.

The algorithm proceeds by trying a fixed number (FN) of neighbourhood moves at each temperature t, while the temperature is gradually decreased. The geometric cooling schedule is used by which in each iteration the temperature decreases as follows: ti = α ⋅ ti −1 where α is a constant and ti is i-th temperature. The SA proceeds until no improvement is made in a predetermined number of consecutive temperatures (CT). The SA has three parameters, i.e., FN, CT, and initial temperature (t0), that are tuned in

190


calibration section. The steps of the proposed local search engine are described in Table 7.

4

Experimental evaluation

In this section, the general performance of the proposed models and the metaheuristic algorithm are evaluated. The operators and the parameters of the proposed GSA are tuned at first. Then, the proposed mathematical models’ capability of solving the FSDGS instances against Salmasi et al. (2010) model is evaluated. Then, the performance of the proposed GSA metaheuristic algorithm is compared to the TS as well as the ACO algorithms proposed by Salmasi et al. (2010). The evaluation experiment is conducted based on the same benchmark previously used in Salmasi et al. (2010) and Salmasi (2005). This benchmark includes 270 instances in two, three, and six machine problems. The number of groups is between two to 16 groups and the number of jobs in a group is between two to ten jobs. The run time and the setup times are randomly generated from uniform distributions between (1, 20) and (1, 400), respectively. More details about the benchmark could be found in Salmasi (2005). The algorithms are implemented in Borland C++ and run on a PC with 2.0 GHz Intel Core 2 Duo and 2 GB of RAM memory. The relative percentage deviation (RPD) is used as a common performance measure to compare the methods. RPD is obtained from the following formula: RPD =

Alg sol − Minsol ⋅100 Minsol

(31)

where Minsol is the best TCT obtained for each problem instance. In other words, it is the optimal solution for those instances solved by MILP models to optimality, or it is the best TCT achieved by any of the algorithms under consideration. Algsol is the TCT obtained for a given algorithm. The stopping criterion used to test algorithms is set at a limit CPU time fixed to nm seconds. This stopping criterion allows for more time as the number of jobs or machines increases. After each 1/3nom seconds, the algorithm refresh all the population except the best one.

4.1 Parameter tuning In this section, the GSA’s performance considering the operators and parameters mentioned above is studied. To design an experiment, there exist several statistical approaches. Although the most widely used approach is a full factorial design, it is not always efficient because it becomes costly and difficult to perform when the number of factors is significantly large. To reduce the number of required trials, in Taguchi method (Phadke, 1989), a set of orthogonal arrays is developed. They eventually reduce the number of experiments, but still provide adequate information. Taguchi separates the factors into two main groups: controllable and noise factors. Noise factors are those over which we have no direct control. Since elimination of the noise factors is impractical and often impossible, the Taguchi method seeks to minimise the effect of noise and to determine optimal level of important controllable factors based


191

on the concept of robustness (Phadke, 1989). Besides determining the optimal levels, Taguchi establishes the relative significance of individual factors in terms of their main effects on the objective function. Taguchi has created a transformation of the repetition data to another value which is the measure of variation. The transformation is the signal-to-noise (S/N) ratio which explains why this type of parameter design is called robust design (Phadke, 1989). Here, the term ‘signal’ denotes the desirable value (mean response variable) and ‘noise’ denotes the undesirable value (standard deviation). So the S/N ratio indicates the amount of variation presents in the response variable. The aim is to maximise the signal-to-noise ratio. Since the objective function is minimisation of TCT, its corresponding S/N ratio is: S /N ratio = − 10log10 ( RPD)2

(32)

As mentioned, the GSA controllable factors are: popsize (A), selection mechanism (B), crossover (C), mutation (D), Pm (E), MAX (F), FN (G), CT (H), and T0 (L). Different levels of these factors are shown in Table 8. The associated degree of freedom for these nine factors is 17. Therefore, the selected orthogonal array should have a minimum of 13 rows and nine columns to assess the five factors. From standard table of orthogonal arrays, the L27 is selected as the fittest orthogonal array design which fulfils our all minimum requirements. But this orthogonal array still entails some modifications to adapt itself to our experimental design (see more details about how to modify in Phadke, 1989). Table 8

Factors and their levels

Factor

Symbol

Level

Type

Popsize

A

3

A(1) – 15, A(2) – 30, A(3)–45

Selection

B

2

B(1) – Roulette wheel, B(2) – tournament

Crossover

C

3

C(1) – SBOX, C(2) – UPC, C(3) – MTPC

Mutation

D

3

D(1) – Swap, D(2) – shift, D(3) – inversion

Pm

E

3

E(1) – 0.1, E(2) – 0.2, E(3) – 0.3

MAX

F

3

F(1) – 10, F(2) – 20, F(3) – 30

FN

G

3

G(1) – 10, G(2) – 20, G(3) – 30

CT

H

3

H(1) – 2, H(2) – 5, H(3) – 10

T0

L

3

L(1) – 10, L(2) – 20, L(3) – 30

To conduct the experiment, we generate a set of 45 instances in different sizes. After obtaining the results of Taguchi experiment, RPDs are transformed into S/N ratio. Figure 1 shows the average S/N ratio obtained at different levels of each factor. As clear in Figure 1, the results demonstrate that factors A, D, E, F, and L have almost insignificant impact on the response variable. Therefore, the best level of these factors is chosen [i.e., A(2), D(2), E(3), F(3) and L(2)], they, then, are taken out of the further experiment. In this case, 4 factors of B, D, G and H are remaining. The new degree of freedom becomes 7. Considering all the situations, the experiment is designed based on orthogonal array . The results show that B(2), C(2), G(1) and H(1) are the best levels for factors B, C, G and H, respectively.

192


Figure 1

The mean S/N ratio plot for each level of the factors (see online version for colours)

To explore the relative significance of individual factors in terms of their main effects on the objective function, the analysis of variance (ANOVA) is conducted. The results of analysis are presented in Table 9. Crossover operator has the greatest effect on the quality of the algorithm with relative importance of 45.5%. The next important factors are FN and selection mechanism with relative importance of 20.3% and 17.161%, respectively. Factor CT, among the four factors, has the least impact on performance of our GSA with 7.9%. Generally, whole variation existing in the experiment is defined by these four factors, and only 8.7% of the variation is defined by noise factors. Table 9

ANOVA table for S/N ratio

Factor

df

SS

MS

F

Relative importance

Cumulative

p-value

Selection

1

11.021

11.021

17.161

17.6

17.6