A discrete differential evolution algorithm for multi-objective permutation flowshop scheduling M. Baioletti, A. Milani, V. Santucci Dipartimento di Matematica e Informatica Università degli Studi di Perugia Via Vanvitelli, 1 Perugia, Italy Abstract. Real-world versions of the permutation flowshop schedul- ing problem (PFSP) have a variety of objective criteria to be optimized simultaneously. Multi-objective PFSP is also a relevant combinatorial multi-objective optimization problem. In this paper we propose a multi- objective evolutionary algorithm for PFSPs by extending the previously proposed discrete differential evolution scheme for single-objective PF- SPs. The novelties of this proposal reside on the management of the evolved Pareto front and on the selection operator. A preliminary ex- perimental evaluation has been conducted on three bi-objective PFSPs resulting from all the possible bi-objective combinations of the criteria makespan, total flowtime and total tardiness. Introduction The Permutation Flowshop Scheduling Problem (PFSP) is an important type of schedul- ing problem which has many applications in manufacturing and large scale product fabrication. In this problem there are n jobs J1 , . . . , Jn and m machines M1 , . . . , Mm . Each job Ji is composed by m operations Oi1 , . . . , Oim . The generic operation Oij can be executed only by the machine Mj and its given processing time is pij . Moreover, the execution of any operation cannot be interrupted (no pre-emption) and job passing is not allowed, i.e., the jobs must be executed using the same order in every machine. The goal of PFSP is to find the optimal job permutation π = hπ(1), . . . , π(n)i with respect to a given objective function. Three important criteria are to minimize the total flowtime (TFT), the makespan (MS) and the total tardiness (TT) defined as follows: n X T F T (π) = C(π(i), m) (1) i=1 M S(π) = max C(π(i), m) = C(π(n), m) (2) i=1,...,n n X T T (π) = max{C(π(i), m) − dπ(i) , 0} (3) i=1 where C(h, j) is the completion time of the operation Ohj and is computed by the following recursive equation C(π(i), j) = pπ(i),j + max{C(π(i − 1), j), C(π(i), j − 1)} for i, j ≥ 1, while the terminal cases are C(π(0), j) = C(i, 0) = 0. In equation (3), for each job h, also a given delivery date dh is considered. The minimization of each one of these criteria is computationally hard. Indeed, both the TFT and TT problems are NP-hard for m ≥ 2, while the MS minimization becomes NP-hard when m > 2. Many single-optimization algorithms exists, either exact or approximate, for in- stance: heuristic techniques, local searches or evolutionary algorithms [2]. Anyway, in this paper we investigate the PFSP problem as a multi-objective optimization problem, in which the goal is to find a set of job permutations which are good enough with respect to two or more contrasting criteria, i.e. a set of Pareto optimal solutions. Given k objective functions f1 , . . . , fk , a solution x dominates a solution x0 (de- noted by x ≺ x0 ) if fi (x) ≤ fi (x0 ) for i = 1, . . . , k, and there exists at least an index j ∈ {1, . . . , k} such that fj (x) < fj (x0 ). A solution x is Pareto optimal if there exist no other solution x0 such that x0 ≺ x. The Pareto optimal set is the set of all the Pareto optimal solution. If two solutions x and x0 are such that neither x ≺ x0 nor x0 ≺ x, then x and x0 are incomparable. Since the Pareto set is in general very large, the goal is to find an approximation of this set, i.e., a set composed by incomparable solutions which is as close as possible to the Pareto optimal set. One of the most promising approaches to solve multi-objective optimization problems is to use evolutionary algorithms [1]. In the context of multi-objective PFSP, many approaches have been proposed. The surveys [3, 7] describe and compare many algorithms for PFSP with all the three possi- ble combinations of two objectives among TFT, MS and TT. In this paper we describe an algorithm for multi-objective optimization which is based on Differential Evolution for Permutation (DEP) [6]. DEP is a discrete differential evolution algorithm which directly operates on the permutations space and hence is well suited for permutation optimization problems like PFSP. Indeed, in [6] and in [5], it was shown that DEP reaches state-of-the-art results with respect to total flowtime and makespan single objective optimization. Here, DEP has been extended in order to handle multi-objective problems. The preliminary experimental results show that its performances are comparable with state-of-the-art algorithms. The rest of the paper is organized as follows. The second section describes the classical Differential Evolution algorithm. Its extension to the permutations space and multi-objective PFSP is introduced in the third section. An experimental investigation of the proposed approach is provided in the fourth section, while conclusions are drawn at the end of the paper. The Differential Evolution algorithm In this section we provide a short introduction to Differential Evolution (DE) algorithm. For more detail see [4]. Differential Evolution (DE) is a powerful population-based evo- lutionary algorithm for optimizing non-linear and even non-differentiable real functions in Rn . The main peculiarity of DE is to exploit the distribution of the solutions’ differ- ences in order to probe the search space. DE initially generates a random population of N P candidate solutions x1 , . . . , xN P uniformly distributed in the solutions space. At each generation, DE performs mutation and crossover in order to produce a trial vector ui for each individual xi , called target vector, in the current population. Each target vector is then replaced in the next genera- tion by the associated trial vector if and only if the produced trial is fitter than the target. This process is iteratively repeated until a stop criterion is met (e.g., a given amount of fitness evaluations has been performed). The differential mutation is the core operator of DE and generates a mutant vector vi for each target individual xi . The most used mutation scheme is “rand/1” and it is defined as follows: vi = xr0 + F · (xr1 − xr2 ) (4) where r0 , r1 , r2 are three random integers in [1, N P ] mutually different among them. xr0 is called base vector, xr1 −xr2 is the difference vector, and F > 0 is the scale factor parameter.In [4] it is argued that the differential mutation confers to DE the ability to automatically adapt the mutation step size and orientation to the fitness landscape at hand. After the mutation, a crossover operator generates a population of N P trial vectors, i.e. ui , by recombining each pair composed by the generated mutant vi and its corre- sponding target xi . The most used crossover operator is the binomial one that builds the trial vector ui taking some components from xi and some other ones from vi according to the crossover probability CR ∈ [0, 1]. Finally, in the selection phase, the next generation population is selected by a one- to-one tournament among xi and ui for 1 ≤ i ≤ N P . Discrete Differential Evolution for Multi-Objective Optimization In this section we describe the proposed Multi-Objective Differential Evolution for Per- mutation (MODEP) which directly evolves a population of N P permutations π1 , . . . , πN P . With respect to the classical DE, important variations have been made to the genetic operators of mutation, crossover and selection. Moreover, an additional archive of so- lutions is introduced to maintain the evolved Pareto front. To simplify our description, let us restrict to the case of two objective functions f1 and f2 . A population of N P permutations π1 , . . . , πN P is randomly generated at the beginning. At each iteration, a secondary population of trial elements υ1 , . . . , υN P is generated by means of the mutation and crossover operators. Then, a selection operator selects, for i = 1, . . . , N P , which element among υi and πi should be part of the population for the next iteration. The pseudo-code of MODEP is depicted in Alg. 1. Differential Mutation The mutation operator used is the same of DEP [6]. It produces a mutant νi for each population element πi using some algebraic concepts related to the symmetric group of permutations. Here we briefly recall its structure: Algorithm 1 MODEP 1: Initialize Population 2: Update N D 3: while num fit eval ≤ max fit eval do 4: for i ← 1 to N P do 5: νi ← DifferentialMutation(i) (1) (2) 6: υi , υi ← Crossover(πi , νi ) 7: Update N D (1) (2) 8: υi ← SelectChild(υi , υi ) 9: end for 10: for i ← 1 to N P do 11: πi ← Selection (πi , υi ) 12: end for 13: end while 1 Find r0 , r1 , r2 different to i and to each other 2 δ ← πr−12 ◦ π r1 3 S ← RandBS(δ) (S is a sequence of adjacent swaps) 4 L ← Length(S) 5 k ← dF · Le 6 ν i ← π r0 7 for j = 1, . . . , k apply Sj to νi where ◦ is the ordinary permutation composition operator, ·−1 denotes the inverse of a permutation, and RandBS is the randomized bubble sort procedure which allows to decompose a permutation in a sequence of adjacent swaps (that are themselves simple permutations). For more details, see [6]. It is worth to notice that this operator works directly with permutations, simulating from an algebraic point of view, the expression of equation (4). Crossover The crossover operator for permutation representations is the same of DEP and pro- (1) (2) duces two children υi and υi from πi and νi . The details are described in [6]. (1) (2) The two permutations υi and υi are compared with respect to both f1 and f2 . If (1) (2) (1) (2) (1) υi dominates υi , then the trial υi is υi . Analogously, if υi dominates υi , then (2) (1) (2) the trial υi is υi . When υi and υi are incomparable, then one of them is randomly selected to become the trial υi . Selection The selection operator chooses the new population element πi0 between the old element πi and the trial νi . If πi ≺ νi , then πi0 becomes πi , i.e., πi remains in the population. Otherwise, if νi ≺ πi or it is equal to πi , then πi0 becomes νi , that is νi replaces πi in the next generation population. However, if πi and νi are incomparable, then we use a probabilistic method somehow similar to the α-selection described in [6]. Suppose first that f1 (νi ) < f1 (πi ) but f2 (νi ) ≥ f2 (πi ). Then, πi0 becomes νi with (2) probability max{0, α2 − ∆i }, otherwise it retains the old element πi , where (2) f2 (νi ) − f2 (πi ) ∆i = f2 (πi ) is the relative worsening of νi with respect to πi according to f2 . Analogously, if f1 (νi ) ≥ f1 (πi ) and f2 (νi ) < f2 (πi ). Then, πi0 becomes νi with (1) probability max(0, α1 − ∆i ), where (1) f1 (νi ) − f1 (πi ) ∆i = . f1 (πi ) The rationale behind this selection operator is that νi enters the population if it dominates or is equal to πi or, with a small probability, if it is not too worse than πi in one of the objective functions, while it is better than πi in the other objective func- tion. Moreover, note that the probability of accepting a slightly worsening population (h) (h) element linearly shades from αh , when ∆i = 0, to 0, when ∆i = αh , for h = 1, 2. Therefore, the parameters αh regulates how worse νi can be in order to be accepted in the new population: if α1 = α2 = 0 only better elements (in the Pareto sense) can replace old elements in the population. Pareto Front The algorithm keeps updated the approximated Pareto front N D, which contains all the non-dominated elements ever generated and evaluated. Initially N D contains all the non-dominated population elements created during the random initialization. Then, (1) (2) at each generation, all the couples of children υi and υi are used to update N D. A new element υ enters N D if it is not dominated by any element of N D. Moreover, all the elements of N D which are dominated by υ are removed. Experimental Results In this section we report some preliminary experimental results obtained with an imple- mentation of MODEP. The experiments have been performed by solving the well known Taillard’s in- stances with the additional due times given in [3]. These instances are divided in 11 groups of 10 instances with the same values of n and m. The values of n are in the set {20, 50, 100, 200}, while m lies in {5, 10, 20}. The combination (n = 200, m = 5) is not considered. The processing time pij of each instance are randomly generated in {1, . . . , 99}, while the due date of each job Ji are generated by multiplying the value P m j=1 pij for a random factor in [1, 4]. MODEP has been run 10 times for each in- stance and the adopted stopping criterion is the maximum number of evaluations, which has been set to 2000 · n · m. Three combinations of objectives have been considered: (M S, T F T ), (M S, T T ), and (T F T, T T ). For each execution the obtained Pareto front (corresponding to N D) has been analyzed by computing two performance indices: the hypervolume IH and the unary multiplicative epsilon I1 . IH is computed as the area delimited by the solutions of N D and a reference point. I1 compares N D with the best known Pareto front B and is computed as fj (y) I1 = max min max . x∈B y∈N D j=1,2 fj (x) The indices have been computed by averaging over the multiple executions and instances for every combination of n × m. The value for the parameter N P has been set to 100 after some preliminary experi- ments. The parameter F used in the mutation operator is, as in [6], self-adapted during the evolution. Instead, the values for the selection parameters α1 and α2 have been set after a calibration phase according to Table 1. Table 1. Calibration values for α1 and α2 Opt. α1 α2 (M S, T F T ) 0.025 0.015 (M S, T T ) 0.01 0.01 (T F T, T T ) 0.01 0.01 The results of the optimization of (M S, T F T ) are shown in Table 2. MODEP works well on this problem and the values of the second index I1 (whose optimal value is 1) are quite good, while the values for IH (whose optimal value is 1.44) are however good, compared to those reported in [3].It is worth to notice that, fixing n, IH seems to have a decreasing behavior as m increases (except when n = 20). Table 2. Results for (M S, T F T ) n m IH I1 20 5 1.089 1.015 20 10 1.185 1.014 20 20 1.188 1.013 50 5 1.248 1.045 50 10 1.149 1.050 50 20 1.119 1.042 100 5 1.238 1.065 100 10 1.140 1.072 100 20 1.067 1.057 200 10 1.143 1.083 200 20 1.058 1.073 The results of the optimization of (M S, T T ) are shown in Table 3 and are similar to those for (M S, T F T ), even if the decreasing behavior of IH with respect to m is not so apparent. Table 3. Results for (M S, T T ) n m IH I1 20 5 1.241 1.048 20 10 1.136 1.162 20 20 1.071 1.038 50 5 1.219 1.078 50 10 1.140 1.175 50 20 1.195 1.237 100 5 1.221 1.086 100 10 1.137 1.119 100 20 1.138 1.167 200 10 1.138 1.105 200 20 0.976 1.117 Finally, the results of the optimization of (T F T, T T ) are shown in Table 4. Here, while the performances as measured by I1 are still satisfactory, the results of IH are slightly worse than in the previous cases. Table 4. Results for (T F T, T T ) n m IH I1 20 5 1.081 1.026 20 10 1.088 1.193 20 20 1.032 1.038 50 5 0.6525 1.024 50 10 0.899 1.083 50 20 1.021 1.206 100 5 0.540 1.027 100 10 0.660 1.055 100 20 0.799 1.103 200 10 0.588 1.051 200 20 0.641 1.077 Conclusion and Future Work In this paper we have described an algorithm for optimization of multi-objective per- mutation flowshop scheduling problems. Some preliminary experimental results show that this approach is promising and reaches results which are comparable to the state- of-the-art algorithms. As a future line of research, we would like to add to our algorithm some method to enhance the diversity of the population, as done in other evolutionary multi-objective algorithms, like crowding distance or niching techniques. References 1. Carlos A. Coello Coello, Arturo Hernández Aguirre, and Eckart Zitzler. Evolutionary multi- objective optimization. European Journal of Operational Research, 181(3):1617–1619, 2007. 2. J. Gupta and J.E. Stafford. Flowshop scheduling research after five decades. European Journal of Operational Research, (169):699–711, 2006. 3. Gerardo Minella, Rubén Ruiz, and Michele Ciavotta. A review and evaluation of multiob- jective algorithms for the flowshop scheduling problem. INFORMS Journal on Computing, 20(3):451–471, 2008. 4. K.V. Price, R.M. Storn, and J.A. Lampinen. Differential Evolution: A Practical Approach to Global Optimization. Springer, Berlin, 2005. 5. Valentino Santucci, Marco Baioletti, and Alfredo Milani. Solving permutation flowshop scheduling problems with a discrete differential evolution algorithm. submitted to AI Commu- nication. 6. Valentino Santucci, Marco Baioletti, and Alfredo Milani. A differential evolution algorithm for the permutation flowshop scheduling problem with total flow time criterion. In Parallel Problem Solving from Nature - PPSN XIII - 13th International Conference, Ljubljana, Slove- nia, September 13-17, 2014. Proceedings, pages 161–170, 2014. 7. M.M. Yenisey and B. Yagmahan. Multi-objective permutation flow shop scheduling problem: Literature review, classification and current trends. Omega, (45):119–135, 2014.