Genetic local search for multi-objective flowshop scheduling problems

European Journal of Operational Research 167 (2005) 717–738 www.elsevier.com/locate/dsw

Genetic local search for multi-objective flowshop scheduling problems Jose´ Elias Claudio Arroyo, Vinı´cius Amaral Armentano

*

Departamento de Engenharia de Sistemas Faculdade de Engenharia Ele´trica e de Computac¸a˜o Universidade Estadual de Campinas Caixa Postal 6101, CEP 13083-970, Campinas-SP, Brazil Available online 9 September 2004

Abstract This paper addresses flowshop scheduling problems with multiple performance criteria in such a way as to provide the decision maker with approximate Pareto optimal solutions. Genetic algorithms have attracted the attention of researchers in the nineties as a promising technique for solving multi-objective combinatorial optimization problems. We propose a genetic local search algorithm with features such as preservation of dispersion in the population, elitism, and use of a parallel multi-objective local search so as intensify the search in distinct regions. The concept of Pareto dominance is used to assign fitness to the solutions and in the local search procedure. The algorithm is applied to the flowshop scheduling problem for the following two pairs of objectives: (i) makespan and maximum tardiness; (ii) makespan and total tardiness. For instances involving two machines, the algorithm is compared with Branchand-Bound algorithms proposed in the literature. For such instances and larger ones, involving up to 80 jobs and 20 machines, the performance of the algorithm is compared with two multi-objective genetic local search algorithms proposed in the literature. Computational results show that the proposed algorithm yields a reasonable approximation of the Pareto optimal set.
2004 Elsevier B.V. All rights reserved. Keywords: Multi-objective combinatorial optimization; Metaheuristics; Genetic local search; Flowshop scheduling

1. Introduction It is well known that in a large number of real problems, production scheduling managers face multiple, often conflicting, decision criteria. Up to the 1980s, scheduling research was mainly concentrated on optimizing single performance measures such as makespan (Cmax), total flowtime (F), maximum tardiness (Tmax), *

Corresponding author. E-mail addresses: [email protected] (J.E.C. Arroyo), [email protected] (V.A. Armentano).

0377-2217/$ - see front matter
2004 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2004.07.017

718

J.E.C. Arroyo, V.A. Armentano / European Journal of Operational Research 167 (2005) 717–738

total tardiness (T) and number of tardy jobs (nT). Cmax and F are related to maximizing system utilization and minimizing work-in-process inventories, respectively, while the remaining measures are related to job due dates. Dileepan and Sen (1988) and Fry et al. (1989) survey the scarce research, at that time, on single-machine multi-objective scheduling problems. The survey by Nagar et al. (1995a) includes a few multiple machine problems. T
kindt and Billaut (2001) present a comprehensive state-of-art survey on multi-criteria machine scheduling, including single-machine, parallel machines and flowshop. They list 105 papers, which shows the increasing interest to such an important area. Multi-objective optimization was originally conceived with finding Pareto optimal solutions (Pareto, 1981), also called efficient solutions. Such solutions are nondominated, i.e., no other solution is superior to them when all objectives are taken into account. From a decision maker perspective, the choice of a solution from all presented efficient solutions, is called a posteriori approach. In this paper we adopt a posteriori approach for multi-criteria flowshop scheduling problems to obtain approximate efficient solutions for the following pairs of objectives: (Cmax, Tmax) and (Cmax, T). The sequence of jobs processed in every machine is the same, characterizing a permutational flowshop. The literature on this type of problem is rather scarce. Branch-and-Bound algorithms for two-machine problems have been proposed. Daniels and Chambers (1990) address the pair of objectives (Cmax, Tmax), Liao et al. (1997) deal with the pairs (Cmax, T) and (Cmax, nT), while Sayin and Karabati (1999) consider (Cmax, F). These papers report results for instances involving at most 30 jobs. From the heuristic viewpoint, Daniels and Chambers propose the only multi-objective constructive heuristic we are aware of, for two and more than two machines. Murata et al. (1996) develop a genetic algorithm to minimize (Cmax, T, F). Ishibuchi and Murata (1998) extend the above approach and include a local search procedure in the genetic algorithm to minimize (Cmax, Tmax, F) and (Cmax, Tmax). There is also a literature on multi-objective flowshop scheduling which can be classified as a priori and interactive approaches. In the a priori approach, the decision maker expresses his preference relative to the objectives in one of two ways. The first one consists of attaching weights to each objective and combining them in a linear function (Rajendran, 1995; Nagar et al., 1995b; S ß erifog˘lu and Ulusoy, 1998). In the second approach, objectives are ranked in decreasing order of importance; the problem is solved for the first objective, and then the second problem is solved for the second objective under the constraint that the optimal solution of the first objective does not change. This single-objective problem process is continued until the last objective (Rajendran, 1992; Gupta et al., 1999, 2001, 2002; T
kindt et al., 2001). In the interactive approach, the decision maker intervenes in the optimization process to guide it to most suitable solutions (Bernardo and Lin, 1994). This paper suggests a new genetic local search algorithm to obtain approximate efficient solutions for multi-objective combinatorial optimization (see Ehrgott and Gandibleux, 2000 for an annotated bibliography). The algorithm is then applied to the flowshop problems mentioned above. Metaheuristics such as genetic algorithms (Holland, 1975), tabu search (Glover and McMillan, 1986) and simulated annealing (Kirkpatrick et al., 1983) were originally conceived for single-objective combinatorial optimization and the success achieved in their application to a very large number of problems has stimulated researchers to extend them to multi-objective combinatorial optimization problems. Jones et al. (2002) list 115 articles published in the period of 1991–1999 and note that almost 80% of the articles are dedicated to real problems, especially in the discipline of engineering. This number reflects not only the increasing awareness of problems with multiple objectives, but also that metaheuristics are effective techniques to cope with such problems. Moreover, 70% of the articles use genetic algorithms as the primary metaheuristic, and the common argument for such a preference is that such algorithms generate at each iteration a population of solutions and in multi-objective optimization the goal is to obtain the set of efficient solutions. Nevertheless, local search based metaheuristics, also provide at each iteration a ‘‘population’’ comprised of the neighboring solutions. Therefore suitable implementation designs of tabu search and simulated annealing can lead to algorithms that are very competitive to genetic algorithms, as in single-objective optimization. Coello

J.E.C. Arroyo, V.A. Armentano / European Journal of Operational Research 167 (2005) 717–738

719

(1999) and Van Veldhuizen and Lamont (2000a) present critical reviews and classification of the most important approaches to genetic algorithms for multi-objective optimization. The new genetic local search multi-objective algorithm proposed here has the following features: elitism, use of the concept of Pareto dominance to assign suitable fitness values to the solutions in order to preserve dispersion, i.e., uniformly distribution of the nondominated points, and a parallel local search for search intensification on distinct regions. The algorithm is tested on instances of flowshop problems for the pairs of objectives (Cmax, Tmax) and (Cmax, T). For instances involving two machines and at most 25 jobs, the efficient solutions are obtained by the Branch-and-Bound algorithms proposed by Daniels and Chambers (1990) and Liao et al. (1997). For larger instances involving up to 80 jobs and 20 machines, the algorithm performance is compared to the genetic local search algorithms proposed by Ishibuchi and Murata (1998) and Jaszkiewicz (2002). The paper is organized as follows. Section 2 contains a description of the problem and Section 3 presents the proposed algorithm. Section 4 describes the application to flowshop problems and computational results are reported on Section 5. Conclusions are presented in Section 6.

2. Problem statement and a multi-objective genetic local search review Given a vector function of r components f = (f1, . . . ,fr) defined on a finite set X, consider the multi-objective combinatorial optimization problem: Minimize f ðxÞ ¼ ðf1 xÞ ¼ z1 ; . . . ; fr ðxÞ ¼ zr Þ subject to x 2 X: The image of a solution x 2 X is the point z = f(x) in the objective space. A point z dominates z 0 , if zj ¼ fj ðxÞ 6 z0j ¼ fj ðx0 Þ, "j and zj < z0j for at least one j. A solution x dominates x 0 if the image of x dominates the image of x 0 . A solution x* 2 X is Pareto optimal (or efficient) if there is no x 2 X such that z = f (x) dominates z* = f (x*). Weighted Tchebycheff functions play a major role in the characterization of efficient Pr points for the above problem, as shown in Steuer (1986, chapter 14). Let K ¼ fk 2 Rr : kj P 0; j¼1 kj ¼ 1g be a set of weights and z
j ¼ minx fzj ¼ fj ðxÞg; j ¼ 1; . . . ; r, which is called the ideal point. Consider the following weighted Tchebycheff function, s1 ¼ max fkj ðfj ðxÞ  z
j Þg j¼1;...;r

and let M be the set of points which minimizes s1. Then there exists
z 2 M such that
z is efficient. Consider now the augmented weighted Tchebycheff function s1aug ¼ max fkj ðfj ðxÞ  z
j Þg þ q j¼1;...;r

r X

kj ðfj ðxÞ  z
j Þ;

j¼1

where q is a sufficiently small number. Then
z minimizes s1aug if and only if
z is efficient. The weighted linear function r X kj fj ðxÞ s1 ¼ j¼1

and the s1, s1aug functions have been used as search directions to guide multi-objective metaheuristics (see for example, Gandibleux and Fre´ville, 2000; Jaszkiewicz, 2002).

720

J.E.C. Arroyo, V.A. Armentano / European Journal of Operational Research 167 (2005) 717–738

The hybridization of genetic algorithms and local search, called genetic local search, has been applied to a variety of single-objective combinatorial problems. The role of the local search is to enhance the intensification process in the genetic search. As an example, Murata et al. (1996) show that a genetic local search algorithm outperforms a ‘‘pure’’ genetic algorithm for minimizing the makespan in a flowshop. Reeves and Yamada (1998) propose a genetic local search algorithm which is competitive with the tabu search algorithm developed by Nowicki and Smutnicki (1996), which is probably the best algorithm for this problem (see also Watson et al., 2002). We are aware of few multi-objective genetic local search algorithms (MOGLS). The first one is due to Ishibuchi and Murata (1998), and denoted IM-MOGLS algorithm. An iteration of the IM-MOGLS algorithm starts with a population P with N solutions, and a set of random weights is selected to determine the s1 function. The operators of selection, recombination and mutation are applied to the elements of P until reaching a population P 0 with NNelite solutions. Then Nelite solutions are randomly selected from the current set of nondominated solutions and included in the set P 0 . A restricted local search is applied to each solution in P 0 . More specifically, a random neighbor x 0 of x 2 P 0 is generated and if its fitness is better than that of x, it replaces x, otherwise, the local search initiated from x terminates. The number of neighborhoods examined from each x 2 P 0 is limited, and a new population P is formed to start a new iteration. The algorithm was applied to two instances of the flowshop scheduling problem: one with 10 jobs and 5 machines for the criteria (Cmax, Tmax) and the other with 20 jobs and 10 machines for the criteria (Cmax, Tmax, F). They show that their results are better than those obtained by the genetic algorithm VEGA, proposed by Schaffer (1985). More recently, Murata et al. (2000) proposed an extension of above method, denoted cellular multiobjective genetic local search algorithm. In this approach, a fixed number of weights, or cells with weight vectors are uniformly (not randomly) generated in the objective space. Each solution is located in a cell and is assigned its own weight vector. Two parents for generating a new solution in a cell C are selected from its k nearest cells, including cell C. The fitness of each neighbor is recalculated based on the weight vector of C. If the solution resulting from the crossover and mutation is far from C, then this solution is assigned to the weight vector of its nearest cell in order to guide the local search process. The authors test the cellular procedure on 10 trials for each of 10 randomly generated instances of the flowshop scheduling problem, involving 20 jobs and 10 machines for the criteria (Cmax, T). The results are not very conclusive, because although the proposed extension outperforms the original one, it is very competitive with the extension without local search. Another MOGLS algorithm was proposed by Jaszkiewicz (2002), and is denoted J-MOGLS. An iteration starts by drawing random weights from a pre-specified set of weights to define a scalarizing function. The B best and distinct solutions according to such a scalarizing function are selected from a current set CS to form a temporary population TP. Two randomly selected solutions of TP are recombined and generate an offspring that is submitted to a local search. If the solution resulting from the local search is better than the worst solution in TP and distinct from all solutions in TP then such a solution, say x 0 , is included in CS and the set of the best nondominated solutions is updated. The J-MOGLS algorithm does not use mutation and recombination is the only component of genetic algorithms that is preserved. From this viewpoint it has a similar structure to scatter search (Glover, 1998). An interesting feature lies in the dynamic updating of CS, which in turn generates TP. More specifically, the solution x may replace a solution, say x 0 0 , from CS that could have been selected as a parent for recombination in the next iteration. This is an aggressive strategy compared to the IM-MOGLS, where selection, recombination and mutation are applied to the population P to form population P 0 and then a restricted local search is applied to all elements of P 0 and the Nelite solutions. Jaszkiewicz (2002) applies his algorithm to several instances of the symmetric traveling salesman problem with two and three objectives and concludes that the quality (measured by a utility function) of the nondominated points generated by the J-MOGLS algorithm is much better than those generated by the IM-MOGLS

J.E.C. Arroyo, V.A. Armentano / European Journal of Operational Research 167 (2005) 717–738

721

algorithm. He also reports that better results are obtained when using the s1 function, compared to the s1 function. Talbi et al. (2001) propose a simplification of genetic local search, where multi-objective local search is applied when the genetic algorithm stops. The proposed algorithm was applied to an instance with 100 jobs and 10 machines of the flowshop scheduling problem for the criteria (Cmax, T).

3. The MOGLS algorithm The proposed algorithm makes use of the following strategies: • Elitism: a subset of the current set of nondominated solutions is included in the population. • Use of the concept of Pareto dominance in order to classify the elements of the population and a mechanism to assign suitable fitness values in order to promote dispersion. • Use of a parallel multi-objective local search based on the concept of Pareto dominance.

3.1. Elitism The elitism scheme is similar to that used by Ishibuchi and Murata (1998). A set Et with Nelite solutions is selected from the set of the current nondominated solutions NDt at iteration t and is included in the population Pt+1, i.e., at iteration t+1 we have the combined population Pt+1 [ Et. 3.2. Classification of the population In order to assign fitness to the elements of the combined population Pt [ Et1 at iteration t, we first classify them. For this, we use the Nondominated Sorting technique suggested by Goldberg (1989) and implemented by Srinivas and Deb (1995) and Deb et al. (2002). This technique orders the population into K sub-populations or frontiers, F1, F2, . . . ,FK, according to the dominance relationship between the solutions. The set F1 corresponds to the set of nondominated solutions. The set Fi, i = 2, . . . ,K contains the nondominated solutions when the sets F1, . . . ,Fi1 are removed from the combined population. Fig. 1 shows an example for a bi-objective minimization problem in the objective space, where the population f(Pt [ Et1) is classified into K = 3 frontiers f(F1), f(F2) and f(F3). The nondominated points correspond to f(F1) and the points in f(F2) dominate the points in f(F3).

f ( Pt ∪ Et –1 )

f2 f ( F2 )

f ( F3 )

f ( F1 )

f1 Fig. 1. Classification of the combined population.

722

J.E.C. Arroyo, V.A. Armentano / European Journal of Operational Research 167 (2005) 717–738

We use the fast nondominated sort algorithm proposed by Deb et al. (2002) for the classification. This algorithm has computational complexity O(rÆN2), where r is the number of objectives and N is the size of the population. 3.3. Fitness assignment After classifying the population Pt [ Et1 into K sub-populations F1, F2, . . . ,FK, a fitness value should be assigned to the solution of each set. It is clear that the fitness values of the solutions in set Fk should be greater than those of the set Fj, j > k. The fitness values assigned to the solutions in Fk are based on a technique proposed by Deb et al. (2002) that tries to preserve the dispersion of the population Pt [ Et1. For each point z 2 f(Fk), an estimate of the points which are located around z is determined. This estimate is given by the crowding distance, which is computed by the following algorithm: Algorithm 3.1. (Crowding distance) Input: Fk,jFkj = np (sub-population of size np) Output: dist(z), "z 2 f(Fk) (crowding distance in the objective space) (1) For each point zi 2 f(Fk) do dist(zi) = 0. (2) For each objective j do (2.1) Order the points in f(Fk) according to increasing values of the jth objective obtaining z1, . . . ,znp. (2.2) dist(z1) = v, dist(znp) = v, where v is a large number. (2.3) For i = 2 to (np1) do dist(zi) = dist(zi) + (fj(zi+1)fj(zi1)), where fj(zi) is the value of the jth objective function at the point zi. Fig. 2 shows the cuboid and the crowding distance of the point zi, for the case of two objectives. In the genetic algorithm proposed Deb et al. (2002), binary tournament is used for selection in the following way. Randomly select solutions two points from Pt [ Et1 If both points belong to the frontier Fk, then the solution with greatest crowding distance is chosen. Otherwise, the point that lies in the lowest frontier Fk is selected. In order to assign fitness values to all points, we suggest a normalization of the crowding distances, involving the maximum and minimum distances. This is shown in the algorithm below: Algorithm 3.2. (Normalized crowding distance) Input: F1, F2, . . . ,FK, (K sub-populations) Output: dist(z), "z 2 f(Fk), k = 1, . . . ,K (normalized crowding distance). f2

f (S ) z1

cuboid of zi zi-1 d2

zi d1

zi+1

dist(zi ) = d 1 (zi ) + d 2 (zi )

znp

f1

d 1 (zi ) = f1 (zi+1 ) – f1 (zi-1 ), d 2 (zi ) = f2 (zi-1 ) – f2 (zi+1 ) Fig. 2. Crowding distance for two objectives.

J.E.C. Arroyo, V.A. Armentano / European Journal of Operational Research 167 (2005) 717–738

723

(1) For each set Fk do: (2) Compute the crowding distance (dist) using Algorithm 3.1 for all points in f(Fk). Note that the extreme points of f(Fk) have a crowding distance v arbitrarily large. (3) Determine the maximum distance Md < v and the minimum distance md 5 0 in the set f(Fk): If np = jFkj 6 2 then do Md = 1 and md = 1. Else Md = max{dist(zj) < v: "zj 2 f(Fk)} and md = min{dist(zj) 5 0: "z j 2 f(Fk)}. (4) For each point zi 2 f(Fk) do If dist(zi) P v then dist(zi) = (dist(zi)v) + (Md + md)/md. Else dist(zi) = dist(zi))/md. Fig. 3(a) shows an example of crowding distances for the entire population obtained by Algorithm 3.1 and Fig. 3(b) shows the normalized crowding distances computed by Algorithm 3.2. For the bi-objective case, the extreme points always have the same normalized crowding distance, namely (Md + md)/md. Note that the more isolated points have larger values of normalized crowding distance. By using the normalized crowding distances, the algorithm below assigns fitness values to all points in the population. Algorithm 3.3. (Fitness assignment) (1) Set maxfitnessK+1 = 0, (maximum fitness value in frontier K + 1). (2) For j = K, K1, . . . ,1 do (2.1) For each solution x 2 Fj do fitnessðxÞ ¼ maxfitnessjþ1 þ distðf ðxÞÞ: (2.2) maxfitnessj = max{fitness(y): y 2 Fj} (maximum fitness of Fj). Fig. 4 shows the fitness values computed by Algorithm 3.3. The fitness value of a solution x 2 Fj is the sum of its crowding distance and the maximum fitness of Fj+1, therefore the solutions in Fj are more likely to be reproduced than the solution in Fj+1. Note that for maximization and minimization problems the fitness is maximized. i.e., high fitness values correspond to high reproduction probabilities. Once the fitness values of the combined population Pt [ Et1 are computed, N solutions are selected by the roulette wheel method to form the mating pool MPt. The recombination and the mutation operators are then applied to MPt to obtain a new population composed of N solutions.

f2

f2 f ( F1 ) f ( F2 )

f ( F1 ) f ( F2 )

f ( F3 )

v 6

v

(12+ 6)/6

v 8

2/1

6/6

(14+ 5)/5

v

14

f ( F3 )

14/5

v

10

12/6 (12+ 6)/6

10/5

5

5/5

v

(14+ 5)/5

f1 (a)

2/1

8/6

12

f1

(b)

Fig. 3. (a) Example of crowding distances computed by Algorithm 3.1. (b) Normalized crowding distances.

724

J.E.C. Arroyo, V.A. Armentano / European Journal of Operational Research 167 (2005) 717–738

f2 f ( F1 ) f ( F2 )

f ( F3 )

5 3

8.8

2

3.3

2

7.8

4

5

7

6

8.8

f1 Fig. 4. Fitness assignment.

3.4. Multi-objective local search In this section we propose a multi-objective local search (MOLS) procedure with the objective of improving a set of nondominated solutions. The method is based on the concept of Pareto dominance and it explores several regions of the solution space in parallel. It starts with an initial set S of nondominated solutions and from each solution x 2 S a neighborhood N(x) is generated. We then construct a subset S 0 of [x 2 SN(x), where S 0 contains the set of neighboring solutions that are nondominated by the solutions of S. If S 0 5 Ø, the process is repeated until a stopping criterion. The Algorithm 3.4 below describes the steps of the (MOLS) method. Algorithm 3.4. (Multi-objective local search (MOLS)) Input: S (a set of dominant solutions). Nmax (maximum number of parallel paths to be explored). Niter (maximum number of iterations). Output: NDI (set of improved nondominated solutions). (0) Initialization: Set t = 1 (iteration counter) and NDI = S. (1) Reduction of the set S: If jSj > Nmax then Reduce the set S by selecting Nmax solutions. Let S
be the set of solutions of the set S that are not selected. (2) Construction of the set of nondominated neighboring solutions S 0 : Set S 0 = Ø For each xk 2 S do Generate the neighborhood of xk, N(xk). For each y 2 N(xk) do If y 62 NDI or y is not dominated by NDI then Set S 0 = nondominated solutions of (S 0 [ {y}). (3) Updating of the set NDI with the new solutions in S 0 . NDI = nondominated solutions of (NDI [ S 0 ). (4) Updating of the set S 0 : Set S 0 = nondominated solutions of (S 0 [ S). If (S 0 5 Ø and t 6 Niter) then Set S = S 0 , t = t + 1 and return to step 1. Else stop.

J.E.C. Arroyo, V.A. Armentano / European Journal of Operational Research 167 (2005) 717–738

725

Note that the set of solutions S
that were not selected in step 1 are taken into account in step 4, when defining the new set of local nondominated solutions S 0 . The reduction scheme used here is based on a clustering technique proposed by Morse (1980), which can be applied to any number of objectives. This technique partitions the set of nondominated solutions S in Nmax clusters, Nmax < jSj, according to the proximity of the solutions. For each cluster the centroid is selected to define a new search direction. It should be pointed out that multi-objective local search based on the concept of Pareto dominance was also suggested by Baykasoglu et al. (1999) and Talbi et al. (2001). In the first work, a set of locally nondominated solutions which are not dominated by the current nondominated approximation set is selected. A random solution is chosen to proceed with the search and the remaining locally nondominated solutions are stored in a list. The process continues until the set of locally nondominated solutions is empty. In this case the first solution of the list is chosen to proceed with the search. In the second work, local search is applied to every point of the nondominated approximation set generated by a genetic algorithm. This generates a new nondominated approximation set, and the process is repeated until no new nondominated point is generated. The following section shows how the multi-objective local search algorithm described above is integrated into the genetic algorithm. 3.5. Pseudocode of MOGLS Algorithm 3.5. (MOGLS) Input: N (size of the population). pR, pM (probability of recombination and mutation). MaxE (maximum number of elite solutions). Nmax (maximum number of parallel paths to be explored by the local search). Niter (maximum number of iterations of the local search). Output: ND (set of nondominated solutions). (1) Initialization: Generate an initial population P1 with N solutions. For each solution in P1 compute the value of the r objectives. Initialize the set of nondominated solutions: ND1 = Ø. Initialize the set of elite solutions: E0 = Ø. Let L be a list containing the solutions explored by the local search. Set L = Ø and t = 1. (2) Classification of the population: Classify the population Pt [ Et1 according to the dominance between solutions and obtain a list of K sets F1, F2, . . . ,FK. (3) Updating the nondominated set: Set NDt = nondominated solutions of (NDt [ F1). (4) Compute fitness assignment: For each x 2 Fi (i = 1,2, . . . ,K) compute (x). (5) Selection: Let MPt = Ø (mating pool). Use a selection method to choose N solutions from Pt [ Et1 and store them in MPt. (6) Updating of the set of elite solutions: If jNDtj > MaxE, then set Nelite = MaxE, otherwise set Nelite = jNDtj. Randomly select Nelite solutions in NDt and store them in Et. (7) Recombination and mutation: Let P 0t ¼ [ be a new population. Repeat N/2 times • Choose two solutions x1, x2 2 MPt and set MPt = MPt{x1,x2}. Recombine x1 and x2 by using a recombination operator and obtain the solutions y1 and y2. • Add y1 and y2 to P 0t with probability pR. Otherwise, add x1 and x2 to P 0t . For each solution x 2 P 0t apply a mutation operator with probability pM.

726

J.E.C. Arroyo, V.A. Armentano / European Journal of Operational Research 167 (2005) 717–738

(8) Evaluation: For each solution in P 0t , compute the value of the r objectives. (9) If t = 1 or t is a multiple of b, then go to step 10. Otherwise, set t = t + 1, P t ¼ P 0t1 and go to step 2. (10) Determining the nondominated solutions in P 0t to be explored by the local search: S = nondominated solutions of P 0t . For each solution x 2 S do If x 2 L (i.e., x was explored by the local search) then remove x from S. Update L: L = nondominated solutions of (L [ S). (11) Multiobjective local search: Apply MOLS (Algorithm 3.4) to the set S (use parameters Nmax, Niter). Let NDI be the set of nondominated solutions at the end of the local search. Update the population P 0t with NDI. (12) Stopping criterion: If the stopping criterion is satisfied, then ND = nondominated solutions of ðNDt [ P 0t Þ and stop. Otherwise, set t = t + 1, P t ¼ P 0t1 and return to step 2. Note that the local search is usually time consuming and for this reason, we apply it every b iterations of MOGLS. In addition, the points explored by local search are stored in a list L to avoid exploring them in future iterations. As a result of the local search we obtain the set NDI of improved solutions. We then update the population P 0t with the solutions of NDI in the following way: (i) jNDIj = jSj. Replace the solutions of S by the solutions of NDI in P 0t . (ii) jNDIj > jSj. Insert the set NDI in P 0t and remove from P 0t the set S and jNDIjjSj randomly chosen solutions. (iii) jNDIj < jSj. Insert the set NDI in P 0t , randomly choose jSjjNDIj solutions of S and remove them from P 0t . Fig. 5 illustrates the structure of the proposed MOGLS algorithm, which is identical to the IM-MOGLS algorithm. However, their components are quite different. 4. Application to flowshop scheduling problems This section describes the implementation of the proposed MOGLS algorithm to flowshop scheduling problems, where n jobs should be processed on m machines. All jobs are available at time zero and preemp-

Fig. 5. Steps and basic procedures of MOGLS.

J.E.C. Arroyo, V.A. Armentano / European Journal of Operational Research 167 (2005) 717–738

727

tion is not allowed. Each job i, i = 1, . . . ,n, has a processing time pij on machine j, j = 1, . . . ,m, and a due date di. 4.1. Representation of solutions and initial population The sequence of n jobs is represented by a vector with n elements. The initial population of size N was generated by a partial enumeration heuristic proposed by Arroyo and Armentano (2004) which constructs a set of Nh nondominated solutions. If Nh < N, then NNh solutions are constructed by the method of generation of diverse solutions for permutation problems proposed by Glover (1998). 4.2. Evaluation of solutions and fitness function Let f(x) = (f1(x), f2(x)) denote the vector objective function of a solution x. The completion time of the job i in machine j, Ci,j, is computed recursively in the following way: C i;j ¼ maxfC i;j1 ; C i1;j g þ pij ;

i ¼ 1; . . . ; n;

j ¼ 1; . . . ; m;

where C 0;j ¼ 0;

j ¼ 1; . . . ; m;

C i;0 ¼ 0;

i ¼ 1; . . . ; n:

A solution x for a bi-objective flowshop is evaluated the values of the objectives f1(x) = Cmax(x) = Pby n Cn,m(x) and f2(x) = Tmax(x) = maxi{Ti(x)} (or f2 ðxÞ ¼ i¼1 T i ðxÞ), where Ti(x) = max{diCi,m(x),0} is the tardiness of job i. The fitness of a solution x is computed by the Algorithm 3.3. 4.3. Selection, recombination and mutation Selection is performed as described in Section 3.3. 4.3.1. Recombination There exists a variety of crossover operators that are suitable for permutation problems. We tested the following operators: order crossover (OX), two-point crossover (2PX), described in (Goldberg, 1989), and the edge recombination crossover (ER), which is presented in (Michalewicz, 1996). 4.3.2. Mutation The mutation operator used here is the insertion operator, which randomly selects a job in the sequence and inserts it in a random position of the sequence. This operator is usually used in flowshop problems (Watson et al., 2002; Ishibuchi and Murata, 1998). 4.4. Local search and a strategy for neighborhood reduction Initially, we identify the nondominated solutions of the new population obtained by the application of the genetic operators. The neighborhood of a solution x is obtained by using insert moves. A job i, located in position k of the sequence is inserted in all positions of the sequence, generating a neighborhood of size (n1)2. The strategy of scanning the entire neighborhood usually consumes the largest part of an Algorithm 3.5. For this reason, we apply a strategy of neighborhood reduction, similar to that proposed by

728

J.E.C. Arroyo, V.A. Armentano / European Journal of Operational Research 167 (2005) 717–738

Armentano and Ronconi (1999). Computational tests disclosed that the best insert position of a job i in position k decreases with the number of local search iterations. Due to this behavior we defined the following expression for the maximum insertion distance, which showed good results. MaxDistance ¼ dðn  1Þ=ðiter  1 þ IÞe; where, iter is the current local search iteration, dae is the least integer P a,
1 if c ¼ 0 I¼ 0 otherwise and c = Niter module iter. For example, if Niter is fixed in 12, the maximum insert distance in the first iteration of the local search is n1, i.e., all insertions are allowed. At iterations 2 and 3 the maximum distance is reduced to d(n1)/2e and d(n1)/3e, respectively. For the iterations 4 and 5 the maximum distance is d(n1)/4e, and for the subsequent iterations it is d(n1)/5e. 4.5. Stopping criterion We adopted a criterion that is suitable to compare distinct algorithms, that consists of the number of evaluated solutions. Computational tests have shown that a suitable expression for this number is (800 + 10n)n2, where n is the number of jobs. This expression is motivated by the size of the neighborhood, which has quadratic order with the number of jobs. For example, for an instance with n = 20 jobs, 1000n2 solutions are evaluated and for an instance with n = 80 jobs, 1600n2 solutions are evaluated. Such an increase is needed because the number of nondominated solutions increases with the number of jobs. 4.6. Determination of parameters Table 1 shows the combination of parameters which yielded the best results. If the local search is activated at every one to three iterations (parameter b) the performance deteriorates. That is, the population should evolve for more iterations in order to obtain a reasonable degree of diversification before the next local search is applied. If the maximum number of parallel paths (parameter Nmax) is small (2 or 3), then we obtain clusters of nondominated points, which do not cover several regions. On the other hand, large values of Nmax (11 or 12) tend to generate a fairly homogeneously disperse, but poor frontier, since its points are in general, dominated by the points obtained with Nmax = 6. The mutation probability is very small for the ER operator, since it usually produces solutions which are very different from their parents even if the parents are the same. This fact does not occur for the other crossover operators, which results in a higher mutation probability. The performance of MOGLS was slightly better for the operator OX, which was then Table 1 Parameters N

100

Probabilities of recombination and mutation OX

2PX

ER

pR = 0.9 pM = 0.3

pR = 0.9 pM = 0.3

pR = 0.8 pM = 0.05

N = size of the population. MaxE = maximum number of elite solutions. b = activation of the local search (every b = 5 iterations). Nmax = maximum number of parallel paths in the local search. Niter = maximum number of iterations in the local search.

MaxE

b

Nmax

Niter

20

5

6

12

J.E.C. Arroyo, V.A. Armentano / European Journal of Operational Research 167 (2005) 717–738

729

selected. The size of the population (N), the degree of elitism (MaxE) and the number of local search iterations (Niter) were obtained experimentally. 5. Computational results All implemented algorithms were programmed in the language C++, and computational tests were executed in a workstation SUN Ultra 60. The algorithm was tested on two sets of instances. The first set consists of instances with two machines and at most 25 jobs, and the efficient solutions are obtained by the Branch-and-Bound (B & B) algorithms proposed by Daniels and Chambers (1990) for the pair (Cmax, Tmax), and Liao et al. (1997) for the pair (Cmax, T). The second set contains instances involving up to 80 jobs and 20 machines and in this case efficient solutions are not known. The performance of the proposed algorithm is compared to the multi-objective local search genetic algorithms suggested by Ishibuchi and Murata (1998) and Jaszkiewicz (2002). In the implementation of the algorithms MOGLS, J-MOGLS and IM-MOGLS, we used the same initial solutions, the same genetic operators and the same stopping criterion as in MOGLS. The following parameters were used for these algorithms: Number of initial solutions N = 100 (initial size of set CS in J-MOGLS); Probability of recombination pR = 0.8; Probability of mutation pM = 0.3. In the implementation of J-MOGLS algorithm, 30 different solutions were used to form the temporary population TP and the maximum size of the set CS is 500. This set is implemented as a queue and when the size of the queue is greater than 500, then the last solution is removed (Jaszkiewicz, 2002). At each iteration, we randomly select the weights and one of the scalarizing functions: the weighted linear utility function and the weighted Tchebycheff utility function. Regarding IM-MOGLS, we use the weighted linear utility function with weights randomly generated at each iteration. In this algorithm, we have to define the maximum number of elite solutions Nelite to be included in the population before the local search and the maximum number of solutions h evaluated at each neighborhood obtained from insert moves. Such parameters are not provided in (Ishibuchi and Murata, 1998), for generic instances. Best results were achieved for Nelite = min{20,jNDtj} and h = n/4, where NDt is the set of nondominated solutions at iteration t. 5.1. Generation of instances The scheme for instance generation is similar to that used by Armentano and Ronconi (1999). Job processing times are uniformly distributed between 1 and 99. We considered four scenarios for the due dates, which are uniformly distributed over the interval [P(1ab/2), P(1a + b/2)] (Potts and Van Wassenhove, 1982), where a and b represent the tardiness factor of jobs and dispersion range of due dates, respectively, and P is a lower bound of the makespan (Taillard, 1993) given by ( ( ) ) j1 n m m X X X X P ¼ max max pij þ min pil þ min pil ; max pij : 16j6m

i¼1

i

l¼1

i

l¼jþ1

i

j¼1

Each scenario is defined by a combination of values for a and b: Scenario Scenario Scenario Scenario

1: 2: 3: 4:

low tardiness factor (a = 0.2) and small due date range (b = 0.6); low tardiness factor (a = 0.2) and wide due date range (b = 1.2); high tardiness factor (a = 0.4) and small due date range (b = 0.6); high tardiness factor (a = 0.4) and wide due date range (b = 1.2).

730

J.E.C. Arroyo, V.A. Armentano / European Journal of Operational Research 167 (2005) 717–738

For each combination of number of jobs n, number of machines m and a matrix of processing times we generated four instances, according to the scenarios above defined. The performance of the genetic algorithms was tested on the following sets of instances. Set 1: n = 10, 15, 20, 25; m = 2; 10 processing times matrices; total of instances: 4 · 10 · 4 = 160. Set 2: n = 20, 50, 80; m = 5, 10, 20; 10 processing times matrices; total of instances: 9 · 10 · 4 = 360. 5.2. Evaluation of computational results Hansen and Jaszkiewicz (1998) discuss several measures and propose a general framework to compare and evaluate an approximation set A to the Pareto optimal set (the reference set R). A cardinal measure defines the proportion of the number of efficient solutions in the approximate set to the number of efficient solutions in R. Note that this proportion can be zero and yet A can be a good approximation. A distance measure aims at establishing an average (maximum) distance between the points of A and R. Van Veldhuizen and Lamont (2000b) propose an additional metric that measures the spread of the points in A. We measure the performance of the nondominated set ND generated by the genetic algorithms relative to the reference set R by using the cardinal measure and the following distance metrics proposed by Czyak and Jaskiewicz (1998) and Ulungu et al. (1998): 1 X Dav ¼ min dðz0 ; zÞ j R j z2R z0 2ND and Dmax ¼ maxf min dðz0 ; zÞg; 0 z2R

z 2ND

where d is defined by
 1 0 0 ðzj  zj Þ ; dðz ; zÞ ¼ max j¼1;...;r Dj

z0 ¼ ðz01 ; . . . ; z0r Þ 2 ND;

z ¼ ðz1 . . . ; zr Þ 2 R;

where Dj is the range of the objective fj among all reference and heuristic solutions. Note that Dav is the average distance from a solution z 2 R to its closest solution in ND while Dmax yields the maximum of the minimum distance from a solution z 2 R to any solution in ND. When the Pareto optimal set is not known and ND 0 is the set of nondominated solutions generated by another heuristic we define the reference set R as the nondominated solutions of ND [ ND 0 and use the same measures mentioned above to assess the quality of ND and ND 0 relative to R. 5.3. Computational results for f = (Cmax, Tmax) 5.3.1. Results for the set of instances 1 The 160 instances of this set were optimally solved (in the Pareto sense) by the B & B algorithm developed by Daniels and Chambers (1990), resulting in 507 efficient solutions. Table 2 shows, for each problem size, the number of (efficient) solutions found by B & B and the algorithms. The algorithms MOGLS, IMMOGLS and J-MOGLS identified the exact optimal Pareto set in 158, 153 and 155 instances, respectively. Table 3 shows the distance measure for the algorithms and Table 4 shows the mean computational times spent by the algorithms. These results show that finding efficient solutions for problems with two machines and few jobs is relatively easy, and that the algorithms MOGLS and J-MOGLS perform slightly better.

J.E.C. Arroyo, V.A. Armentano / European Journal of Operational Research 167 (2005) 717–738

731

Table 2 Performance of the algorithms: number of efficient solutions Problem n · m

B&B

MOGLS

IM-MOGLS

J-MOGLS

NES

NS

NES

NS

NES

NS

NES

10 · 2 15 · 2 20 · 2 25 · 2

122 140 144 101

122 140 142 101

122 140 142 100

121 139 138 100

120 138 133 100

122 140 139 100

122 140 135 98

Total

507

505

504

498

491

501

495

NS: number of solutions found by the algorithm. NES: number of efficient solutions found by the algorithm.

Table 3 Performance of the algorithms: distance measure Problem n · m

MOGLS

IM-MOGLS

J-MOGLS

Dav

Dmax

Dav

Dmax

Dav

Dmax

10 · 2 15 · 2 20 · 2 25 · 2

0.0 0.0 0.0003 0.00006

0.0 0.0 0.0022 0.0003

0.0009 0.0028 0.0078 0.0005

0.0045 0.005 0.0233 0.002

0.0 0.0 0.0029 0.0018

0.0 0.0 0.0106 0.0041

Mean

0.0001

0.0006

0.0030

0.0087

0.0012

0.0037

Table 4 Mean computational time (in seconds) Problem n · m

10 · 2 15 · 2 20 · 2 25 · 2

Branch and Bound Scenario 1

Scenario 2

Scenario 3

Scenario 4

0.014 0.124 0.436 0.268

0.026 0.229 51.957 720.605

0.001 0.183 0.845 0.414

0.012 0.474 499.171 3695.813

MOGLS

IM-MOGLS

J-MOGLS

0.918 1.854 3.626 6.254

0.766 1.624 3.524 6.052

0.927 1.905 3.783 7.122

5.3.2. Results for the set of instances 2 The algorithms MOGLS, IM-MOGLS and J-MOGLS were tested on 360 larger instances of this set, for which efficient solutions are not known. Therefore, the performance is evaluated with respect to a reference set R, which consists of the nondominated solutions from all solutions generated by the three algorithms. Table 5 shows the number of nondominated solutions in the set R and the number of nondominated solutions generated by each algorithm, and Fig. 6 shows the mean number of nondominated solutions for each size. It is clear that the number of such solutions increases as the ratio n/m decreases for n fixed. This conclusion cannot be inferred from the results of IM-MOGLS, and moreover this algorithm yields the least number of nondominated solutions for all problem sizes. Table 5 also shows that for larger instances involving 50 and 80 jobs, 10 and 20 machines, MOGLS generates more nondominated solutions, when compared to J-MOGLS. Fig. 7 shows that the average number of nondominated solutions generated by MOGLS and J-MOGLS is higher for scenarios 2 and 4, associated with wide due date range. In addition, the algorithms of MOGLS, IM-MOGLS and J-MOGLS obtained the least values of Cmax for 248, 151 and

732

J.E.C. Arroyo, V.A. Armentano / European Journal of Operational Research 167 (2005) 717–738

Table 5 Performance of the algorithms: number of nondominated solutions Problem n · m

jRj

NS

NNDS

NS

NS

NNDS

20 · 5 20 · 10 20 · 20 50 · 5 50 · 10 50 · 20 80 · 5 80 · 10 80 · 20

518 928 1097 538 1257 1804 398 1043 2205

491 867 1017 488 1137 1626 368 922 2096

398 592 661 347 786 1054 284 624 1324

447 674 841 446 887 1250 328 740 1518

305 249 339 228 140 112 186 158 208

480 862 1032 436 1120 1721 371 962 2060

387 595 652 362 694 949 288 531 1031

Total

9788

9012

6070

7131

1925

9044

5489

MOGLS

IM-MOGLS

J-MOGLS NNDS

Mean number of nondominated solutions

jRj: number of reference solutions. NS: number of solutions found by the algorithm. NNDS: number of nondominated solutions found by the algorithm.

60

R 50 40 30

MOGLS J-MOGLS IM-MOGLS

20 10 0 20x5

20x10

20x20

50x5

50x10

50x20

80x5

80x10

80x20

Mean number of nondominated solutions

Fig. 6. Mean number of nondominated solutions by problem size.

35 30

R MOGLS J-MOGLS IM-MOGLS

25 20 15 10 5 0 1

2

Scenario

3

4

Fig. 7. Mean number of nondominated solutions by scenario.

196 instances, and the least values of Tmax in 291, 188 and 226 instances, respectively. This simply measures the capability of generating the extreme points of the reference set R.

J.E.C. Arroyo, V.A. Armentano / European Journal of Operational Research 167 (2005) 717–738

733

The nonparametric Wilkoxon signed rank test (Mosteller and Rourke, 1973) was applied to compare the expected values E[NNDS] and E[Dav] of the random variables NNDS(X) and Dav(X), which result from the application of the heuristic X to flowshop problems. This test is very useful in comparing the performance of two heuristics X and Y (Golden and Stewart, 1985) and consists of computing the difference of the random variables values for each instance of every group of 40 instances associated with each problem size in Table 5. For a given group, ranks are assigned to each difference and we test the null hypothesis, say, EX[NNDS] = EY[NNDS] with the alternative hypotheses EX[NNDS] > EY[NNDS] or EX[NNDS] < EY[NNDS] at the significance level a = 0.05. This test was applied to compare MOGLS with J-MOGLS and MOGLS with IM-MOGLS. For the groups of sizes involving n = 20 jobs and the group 80 · 5, the null hypothesis EMOGLS[NNDS] = EJ-MOGLS[NNDS] and the null hypothesis EMOGLS[Dav] = EJ-MOGLS[Dav] are accepted, except for the size 20 · 20, where MOGLS performs better. For the two random variables and the group of size 50 · 5, J-MOGLS is superior, while for the groups of sizes 50 · 10, 50 · 20, 80 · 10, 80 · 20, MOGLS outperforms J-MOGLS. MOGLS is superior to IM-MOGLS for all groups of instances. Table 6 shows that MOGLS performs better than the other algorithms regarding the distance measure, again for larger instances involving 50 and 80 jobs, 10 and 20 machines. Table 7 shows that the mean computational times spent by the algorithms are very similar. Figs. 8(a)–(c) illustrate the initial set and the final sets of points obtained by each algorithm for an instance with 50 jobs and 20 machines. The initial set is the same for all algorithms. Fig. 8(d) shows the reference set composed of 53 nondominated solutions, from which 43 are generated by MOGLS, 11 by IM-MOGLS and 27 by J-MOGLS. Figs. 8(a)–(c) show that the three algorithms generate a set of nondominated points which are fairly disperse along the frontier. The reference set in Fig. 8(d) shows that MOGLS

Table 6 Performance of the algorithms: distance measure Problem n · m

MOGLS

IM-MOGLS

J-MOGLS

Dav

Dmax

Dav

Dmax

Dav

Dmax

20 · 5 20 · 10 20 · 20 50 · 5 50 · 10 50 · 20 80 · 5 80 · 10 80 · 20

0.009 0.008 0.008 0.041 0.010 0.014 0.004 0.011 0.012

0.044 0.048 0.041 0.064 0.051 0.063 0.018 0.053 0.047

0.017 0.027 0.022 0.029 0.051 0.053 0.063 0.038 0.052

0.065 0.091 0.086 0.053 0.125 0.152 0.082 0.101 0.137

0.011 0.008 0.009 0.037 0.022 0.019 0.003 0.019 0.024

0.058 0.042 0.046 0.055 0.058 0.078 0.013 0.078 0.072

Mean

0.012

0.043

0.037

0.092

0.016

0.049

Table 7 Mean computational time (in seconds) n:

20

50

80

m:

5

10

20

5

10

20

5

10

20

MOGLS IM-MOGLS J-MOGLS

5.75 5.62 5.70

8.16 8.36 8.47

12.25 12.34 12.88

80.19 86.05 84.16

152.42 170.22 169.74

247.26 264.56 267.32

351.03 354.33 352.65

625.82 643.62 638.72

1106.2 1143.3 1148.5

734

J.E.C. Arroyo, V.A. Armentano / European Journal of Operational Research 167 (2005) 717–738 (a) MOGLS

4750

Tmax

3750

2750

Initial population Final population Nondominated points

1750

750 3800

4000

4200

4400

C max

4600

4800

5000

(b) IM-MOGLS

4750

Tmax

3750

2750

Initial population Final population Nondominated points

1750

750 3800

4000

4200

4400

C max

4600

4800

5000

(c) J-MOGLS

4750

Tmax

3750

2750

Initial set Final set Nondominated points

1750

750 3800

4000

4200

4400

4600

4800

5000

C max (d) Nondominated points

2750

MOGLS J-MOGLS

2250

Tmax

IM-MOGLS

1750

1250

750 3800

3850

3900

3950

4000

4050

4100

4150

4200

C max

Fig. 8. Initial and final solution sets, and the set of nondominated solutions for an instance.

generates a larger number of nondominated points along the frontier and that the nondominated points from IM-MOGLS are located in very restricted regions of the frontier. 5.4. Computational results for f = ðCmax, TÞ For the sets of instances 1 we aggregate all results obtained for each size, and also omit the analysis of mean number of nondominated solutions for each size and scenario. The heuristic parameters used here are identical to those for f = (Cmax, Tmax). 5.4.1. Results for the set of instances 1 The efficient solutions for the instances of this set were generated by the B & B algorithm proposed by Liao et al. (1997). Instances involving 25 jobs could not be solved due to the excessive computational time required by the B & B, and therefore we tested 120 instances. Table 8 shows the number of solutions found

J.E.C. Arroyo, V.A. Armentano / European Journal of Operational Research 167 (2005) 717–738

735

Table 8 Number of solutions found by the algorithms

Total

B&B

MOGLS

IM-MOGLS

J-MOGLS

NES

NS

NES

NS

NES

NS

NES

474

473

461

468

423

471

452

by the algorithms. The algorithms MOGLS, IM-MOGLS and J-MOGLS found the exact optimal Pareto set in 114, 96 and 110 instances, respectively. Table 9 shows the distance measure for the algorithms. These results show again that finding efficient solutions for problems with two machines and few jobs is relatively easy, and that the algorithms MOGLS and J-MOGLS perform better. 5.4.2. Results for the set of instances 2 Table 10 shows the number of nondominated solutions in the set R and the number of (nondominated) solutions generated by each algorithm, while Table 11 shows the distance measure. From these tables, it follows again that MOGLS outperforms the other algorithms for larger problems, involving 50 and 80 jobs, 10 and 20 machines. Moreover, the algorithms of MOGLS, IM-MOGLS and J-MOGLS obtained the least values of Cmax for 262, 154 and 218 instances, and the least values of T in 271, 119 and 226 instances, respectively. The Wilkoson test was applied to the groups of instances of Table 10. For the random variable NNDS, EMOGLS[NNDS] = EJ-MOGLS[NNDS] for the groups of sizes 20 · 5 and 20 · 10. For the groups of sizes 50 · 5 and 80 · 5, J-MOGLS performs better, while MOGLS outperforms J-MOGLS for the remaining sizes. For the random variable Dav, J-MOGLS is superior for the groups of sizes 205, 50 · 5 and 80 · 5, while MOGLS performs better than J-MOGLS for the remaining sizes. MOGLS is superior to IMMOGLS for all groups of instances. Table 9 Performance of the algorithms: distance measure MOGLS

Mean

IM-MOGLS

J-MOGLS

Dav

Dmax

Dav

Dmax

Dav

Dmax

0.0017

0.0049

0.0184

0.0301

0.0032

0.0124

Table 10 Number of solutions found by the algorithms Problem n · m

j Rj

MOGLS NS

20 · 5 20 · 10 20 · 20 50 · 5 50 · 10 50 · 20 80 · 5 80 · 10 80 · 20 Total

IM-MOGLS NNDS

NS

J-MOGLS NNDS

NS

NNDS

813 1134 1226 811 2068 2435 634 2038 3374

782 1047 1088 771 1966 2260 633 1971 3250

578 822 937 403 1304 1495 492 1686 2051

685 817 852 678 1518 1730 547 1576 2542

346 266 292 155 218 118 121 182 262

786 1087 1122 795 1924 2416 646 1964 3361

581 803 924 415 986 1168 502 1358 1685

14533

13768

9768

10945

1960

14101

8422

736

J.E.C. Arroyo, V.A. Armentano / European Journal of Operational Research 167 (2005) 717–738

Table 11 Performance of the algorithms: Distance measure Problem n · m

MOGLS Dav

Dmax

Dav

Dmax

Dav

Dmax

20 · 5 20 · 10 20 · 20 50 · 5 50 · 10 50 · 20 80 · 5 80 · 10 80 · 20

0.0105 0.009 0.006 0.011 0.012 0.009 0.023 0.016 0.012

0.0501 0.045 0.042 0.064 0.054 0.056 0.029 0.063 0.079

0.0185 0.038 0.037 0.046 0.058 0.081 0.078 0.061 0.069

0.0661 0.119 0.127 0.092 0.120 0.148 0.152 0.217 0.168

0.0107 0.010 0.007 0.008 0.016 0.042 0.009 0.049 0.032

0.0554 0.052 0.052 0.041 0.071 0.118 0.019 0.120 0.152

Mean

0.011

0.048

0.052

0.127

0.019

0.069

IM-MOGLS

J-MOGLS

6. Conclusions This paper presented a new multi-objective genetic local search algorithm, which was applied to multi-objective flowshop problems in order to find an approximation of the Pareto optimal set. The algorithm was tested on several instances involving two pairs of objectives. Comparative results with two algorithms from the literature have shown that, for larger instances, the algorithm is able to find a larger number of approximate Pareto optimal solutions with higher quality in terms of the measures used. The proposed algorithm has few parameters and the same structure as the algorithm proposed by Ishibuchi and Murata (1998). In our opinion, the better performance achieved by our algorithm is due to three components: parallel multi-objective local search and use of the concept of Pareto dominance to compare neighbors and to select elements from the population. The proposed selection scheme uses the crowding distances proposed by Deb et al. (2002) in order to define fitness values, however the improvement achieved by this new information is marginal compared to the improvement generated by the multi-objective local search. Computational results have shown that the algorithm proposed by Jaszkiewicz (2002) is very competitive with ours. This raises the issue whether to use scalarizing functions (SF) or the concept of Pareto dominance (PD) in order to guide the search to a good set of nondominated solutions. This topic has been recently discussed by Jaszkiewicz (2003), who introduces an efficiency index to compare the running times needed by metaheuristics proposed in the literature that use SF or PD to generate solutions of the same quality. Computational tests for the multi-objective set covering problem have shown that both approaches are competitive. An interesting research topic would involve keeping the structure of the most promising metaheuristic implementations proposed in the literature and test the use of SF and PD for several combinatorial problems. For example, the temporary population of J-MOGLS could contain a set of scattered nondominated solutions, and the parallel local search in the proposed MOGLS would be guided by scalarizing functions. This would allow us to compare different algorithm structures and analyze when the use of SF or PD is more advantageous.

Acknowledgments This research was partially funded by the Conselho Nacional de Desenvolvimento Cientıfico e Tecnolo´gico (CNPq). Part of this work was developed while Vinıcius Amaral Armentano was on a sabbatical year at the College of Business and Administration of the University of Colorado at Boulder with a fellowship

J.E.C. Arroyo, V.A. Armentano / European Journal of Operational Research 167 (2005) 717–738

737

from the Fundac¸a˜o de Amparo a` Pesquisa do Estado de Sa˜o Paulo (FAPESP). The authors are also grateful to the valuable comments made by the referees.

References Armentano, V.A., Ronconi, D.P., 1999. Tabu search for total tardiness minimization in flowshop scheduling problems. Computers and Operations Research 26, 219–235. Arroyo, J.E.C., Armentano, V.A., 2004. A partial enumeration heuristic for multi-objective flowshop scheduling problems. The Journal of the Operational Research Society, accepted. Baykasoglu, A., Owen, S., Gindy, N., 1999. A Taboo search based approach to find the pareto optimal set in multiple objective optimization. Engineering Optimization 31, 731–748. Bernardo, J.J., Lin, K.-S., 1994. An interactive procedure for bi-criteria production scheduling. Computers and Operations Research 21, 677–688. Coello, C.A.C., 1999. A comprehensive survey of evolutionary-based multiobjective optimization. Knowledge and Information Systems 1, 269–308. Czyak, P., Jaskiewicz, A., 1998. Pareto simulated annealing––a metaheuristic technique for multiple objective combinatorial optimization. Journal of Multi-Criteria Decision Making 7, 34–47. Daniels, R.L., Chambers, R.J., 1990. Multiobjective flow-shop scheduling. Naval Research Logistics 37, 981–995. Deb, K., Pratap, A., Agarwal, S., Meyarivan, T., 2002. A Fast and Elitist Multi-Objective Genetic Algorithm-NSGA-II. IEEE Transactions on Evolutionary Computation 6, 182–197. Dileepan, P., Sen, T., 1988. Bicriterion static scheduling research for a single machine. OMEGA 16, 53–59. Ehrgott, M., Gandibleux, X., 2000. A survey and annotated bibliography of multicriteria combinatorial optimization. OR Spektrum 22, 425–460. Fry, T.D., Armstrong, R.D., Lewis, H., 1989. A framework for single machine multiple objective sequencing research. OMEGA 17, 595–607. Gandibleux, X., Fre´ville, A., 2000. Tabu search based procedure for solving the 0–1 multiobjective knapsack problem: the two objectives case. Journal of Heuristics 6, 361–383. Glover, F., McMillan, C., 1986. The general employee scheduling problem: An integration of management science and artificial intelligence. Computers and Operations Research 15, 563–593. Glover, F.W., 1998. A Template for Scatter Search and Path Relinking. In: Hao, J.-K., Lutton, E., Ronald, E., Schoenauer, M., Snyers, D. (Eds.), Artificial Evolution, Lecture Notes in Computer Science, 1363. Springer, pp. 13–54. Goldberg, D., 1989. Genetic Algorithms in Search, Optimization, and Machine Learning. Addison-Wesley. Golden, B.L., Stewart, W.R., 1985. Empirical Analysis of Heuristics. In: Lawler, E.L., Lenstra, J.K., Rinnoy Kan, A.H.G., Shmoys, D.B. (Eds.), The Traveling Salesman Problem. Wiley, pp. 207–250. Gupta, J.N.D., Palanimuthu, N., Chen, C.-L., 1999. Designing a tabu search algorithm for the two-stage flow shop problem with secondary criterion. Production Planning and Control 10, 251–265. Gupta, J.N.D., Neppalli, V.R., Werner, F., 2001. Minimizing total flow time in a two-machine flowshop problem with minimum makespan. International Journal of Production Economics 69, 323–338. Gupta, J.N.D., Hennig, K., Werner, F., 2002. Local search heuristics for two-stage flow shop problems with secondary criterion. Computers and Operations Research 29, 123–149. Hansen, M.P., Jaszkiewicz, A., 1998. Evaluating the quality of approximations to the nondominated set, Technical Report, Institute of Mathematical Modelling (1998-7), Technical University of Denmark. Holland, J.H., 1975. Adaptation in Natural and Artificial Systems. University of Michigan Press. Ishibuchi, H., Murata, T., 1998. A multi-objective genetic local search algorithm and its application to flowshop scheduling. IEEE Transactions on Systems, Man and Cybernetics–Part C: Applications and Reviews 28, 392–403. Jaszkiewicz, A., 2002. Genetic local search for multi-objective combinatorial optimization. European Journal of Operational Research 137, 50–71. Jaszkiewicz, A., 2003. Do multiple-objective metaheuristics deliver on their promises. A computational experiment on the set-covering problem. IEEE Transactions on Evolutionary Computation 7, 133–143. Jones, D.F., Mirrazavi, S.K., Tamiz, M., 2002. Multi-objective meta-heuristics: An overview of the current state-of-the-art. European Journal of Operational Research 137, 1–9. Kirkpatrick, S.C., Gellat Jr., D., Vecchi, M.P., 1983. Optimization by simulated annealing. Science 220, 671–680. Liao, C.J., Yu, W.C., Joe, C.B., 1997. Bicriterion scheduling in the two-machine flowshop. Journal of Operational Research Society 48, 929–935. Michalewicz, Z., 1996. Genetic Algorithms + Data Structures = Evolution Programs, third ed. Springer.

738

J.E.C. Arroyo, V.A. Armentano / European Journal of Operational Research 167 (2005) 717–738

Morse, J.N., 1980. Reducing the size of the nondominated set: Pruning by clustering. Computers and Operations Research 7, 55–66. Mosteller, F., Rourke, R.E.K., 1973. Sturdy Statistics. Addison-Wesley. Murata, T., Ishibuchi, H., Tanaka, H., 1996. Multi-objective genetic algorithm and its applications to flowshop scheduling. Computers and Industrial Engineering 30, 957–968. Murata, T., Ishibuchi, H., Gen, M. 2000. Cellular genetic local search for multi-objective optimization. In: Proceedings of Genetic and Evolutionary Computation Conference, Las Vegas, Nevada, USA, pp. 307–314. Nagar, A., Haddock, J., Heragu, S., 1995a. Multiple and bicriteria scheduling: A literature survey. European Journal of Operational Research 81, 88–104. Nagar, A., Heragu, S., Haddock, J., 1995b. A branch and bound approach for a 2-machine flowshop scheduling problem. Journal of the Operational Research Society 46, 721–734. Nowicki, E., Smutnicki, C., 1996. A fast tabu search algorithm for the permutation flow-shop problem. European Journal of Operational Research 91, 160–175. Pareto, V., 1981. Manuel d
E´conomie Politique, fifth ed., Geneva, Librairie Droz. Potts, C.N., Van Wassenhove, L.N., 1982. A decomposition algorithm for the single machine total tardiness problem. Operations Research Letters 1, 177–181. Rajendran, C., 1992. Two-stage flow shop scheduling problem with bicriteria. Journal of the Operational Research Society 43, 871– 884. Reeves, C., Yamada, T., 1998. Genetic algorithms, path relinking and the flowshop sequencing problem. Evolutionary Computation 6, 45–60. Rajendran, C., 1995. Heuristics for scheduling in flowshop with multiple objectives. European Journal of Operational Research 82, 540–555. Sayin, S., Karabati, S., 1999. A bicriteria approach to the two-machine flow shop scheduling problem. European Journal of Operational Research 113, 435–449. Schaffer, J.D., 1985. Multiple objective optimization with vector evaluated genetic algorithms. International Conference on Genetic Algorithms. Morgan Kaufmann, New York, pp. 93–100. S ß erifog˘lu, F.S., Ulusoy, G., 1998. A bicriteria two-machine permutation flowshop problem. European Journal of Operational Research 107, 414–430. Srinivas, N., Deb, K., 1995. Multi-objective optimization function optimization using non-dominated sorting genetic algorithms. Evolutionary Computation 2, 221–248. Steuer, R.E., 1986. Multiple Criteria Optimization––Theory Computation and Application. Wiley. Taillard, E., 1993. Benchmarks for basic scheduling problems. European Journal of Operational Research 64, 278–285. Talbi, El.-G., Rahoual, M., Mabed, M.H., Dhaenens, C., 2001. A hybrid evolutionary approach for multicriteria optimization problems: Application to flow shop. In: Zietler, E., Deb, K., Thiele, L, Coello Coello, C.A., Corne, D. (Eds.), Proceedings of the First International Conference on Evolutionary Multi-Criterion Optimization, Lecture Notes in Computer Science 1993. Springer, Berlin, pp. 416–428. T
kindt, V., Billaut, J.-C., 2001. Multicriteria scheduling problems: A survey. RAIRO Operations Research 35, 143–163. T
kindt, V., Billaut, J.-C., Proust, C., 2001. Solving bicriteria scheduling problem on unrelated parallel machines occurring in the glass bottle industry. European Journal of Operational Research 135, 42–49. Ulungu, E.L., Teghem, J., Ost, C., 1998. Efficiency of interactive multi-objective simulated annealing through a case study. Journal of the Operational Research Society 49, 1044–1050. Van Veldhuizen, D.A., Lamont, G.B., 2000a. Multiobjective evolutionary algorithms: Analysing the state-of art. Evolutionary Computation 8, 125–147. Van Veldhuizen, D.A., Lamont, G.B., 2000. On Measuring Multiobjective Evolutionary Algorithm Performance, 2000 Congress on Evolutionary Computation, vol. 1. Piscataway, New Jersey, pp. 204–211. Watson, J.-P., Barbulescu, L., Whitley, L.D., Howe, A.E., 2002. Contrasting structured and random permutation flow-shop scheduling problems: Search-space topology and algorithm performance. INFORMS Journal on Computing 14, 98–123.