Applied Soft Computing 12 (2012) 2765–2780 Contents lists available at SciVerse ScienceDirect Applied Soft Computing journal homepage: www.elsevier.com/locate/asoc A decomposition-based hybrid multiobjective evolutionary algorithm with dynamic resource allocation Wali Khan Mashwani ∗ , Abdellah Salhi Department of Mathematical Sciences, University of Essex, Wivenhoe Park, Colchester CO4 3SQ, UK a r t i c l e i n f o a b s t r a c t Article history: Different crossover operators suit different problems. It is, therefore, potentially problematic to chose Received 17 November 2011 the ideal crossover operator in an evolutionary optimization scheme. Using multiple crossover operators Received in revised form 25 February 2012 could be an effective way to address this issue. This paper reports on the implementation of this idea, Accepted 8 March 2012 i.e. the use of two crossover operators in a decomposition-based multi-objective evolutionary algorithm, Available online 24 May 2012 but not simultaneously. After each cycle, the operator which has helped produce the better offspring is rewarded. This means that the overall algorithm uses a dynamic resource allocation to reward the better Keywords: of the crossover operators in the optimization process. The operators used are the Simplex Crossover oper- Multi-objective optimization Pareto optimality ator (SPX) and the Center of Mass Crossover operator (CMX). We report experimental results that show Decomposition that this innovative use of two crossover operators improves the algorithm performance on standard test Crossover problems. Results on the sensitivity of the suggested algorithm to key parameters such as population size, Resource allocation neighborhood size and maximum number of solutions to be altered for a given subproblem in the the decomposition process are also included. © 2012 Elsevier B.V. All rights reserved. 1. Introduction Definition 1. Let u = (u1 , u2 , . . ., um )T and v = (v1 , v2 , . . . , vm )T be two vectors in Rm ; u is said to dominate v, denoted as u ≺ v, if and 1.1. Problem definition and basic concepts only if A multiobjective optimization problem (MOP) can be stated as 1. ui ≤ vi for every i ∈ {1, 2, . . ., m} follows: 2. uj < vj for at least one index j ∈ {1, 2, . . ., m}. minimize F(x) = (f1 (x), . . . , fm (x))T (1) subject to x ∈ , (1) Remark 1. For any two given vectors, u and v, there are two pos- sibilities: where is the decision variable space, x = (x1 , x2 , . . ., xn )T is the the vector of decision variables, and F(x) : → Rm is a function that 1. Either u dominates v or v dominates u. associates with every decision vector in an m-vector in Rm with 2. Neither u dominates v nor v dominates u. each entry corresponding to the value of some objective function of the MOP. Rm is referred to as the objective space. When is a closed and connected region in Rn , and all the objectives are continuous Definition 2. A solution x* ∈ is said to be Pareto optimal for prob- in x, problem (1) is a continuous MOP. lem (1), if there is no other solution x ∈ such that F(x) dominates Often, the objectives of problem (1) are in conflict with one F(x* ). F(x* ) is called the Pareto optimal (objective) vector. another or are incommensurable. This means, a single solution Remark 2. Any Pareto optimal point improvement in one objec- in the search space that minimizes all the objectives func- tive must lead to deterioration in at least one other objective. tions simultaneously does not exit. Instead, one tries to find the best trade-offs among the different objectives. These tradeoffs are Definition 3. The set of all Pareto optimal solutions is called the defined in terms of Pareto optimality [10], a formal definition of Pareto set (PS), written as PS = {x ∈ | y ∈ , F(y) ≺ F(x)}. which can be as follows [5–7,18]. Definition 4. The image of the Pareto optimal set (PS) in the objec- tive space is called the Pareto front (PF), PF = {F(x)|x ∈ PS}. In other ∗ Corresponding author. words, the set of all Pareto optimal solutions forms a tradeoff sur- E-mail address:
[email protected](W. Khan Mashwani). face in the objective space. 1568-4946/$ – see front matter © 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.asoc.2012.03.067 2766 W. Khan Mashwani, A. Salhi / Applied Soft Computing 12 (2012) 2765–2780 Remark 3. In MOP, Pareto domination does not define a complete 3. It uses two crossover operators each of which deals with ordering among the solutions in the objective space. It, also, does its own set of subproblems. not measure how much one solution is better than another. 4. The number of subproblems in the two different sets is adjusted dynamically according to a parameter (proba- Over the last two decades, MOP has received growing attention, bility) of the used crossover operator. possibly because it has a wide application. Evolutionary Compu- tation (EC) has also seen significant developments in recent years. The performance of MOEA/D-DRA-CMX+SPX is compared with Multi-objective evolutionary algorithms (MOEAs) have received a that of MOEA/D-DRA-CMX and that of MOEA/D-DRA-SPX on CEC’09 lot of attention due to their ability to find multiple tradeoff solutions test instances [31]. Its sensitivity to key parameters is also empiri- in a single run. They are population-based and, therefore, work with cally investigated. To summarize, the paper addresses the following a population of candidate solutions. issues. Most existing MOEAs dealing with MOP are Pareto dominance based, i.e. they determine the fitness of candidate solutions in the population by Pareto dominance relations with other solu- 1. Would the use of multiple crossovers enhance MOEA/D- tions found in previous searches [4,9,13,20,32,34]. To promote type algorithms? diversity, they use various techniques such as fitness sharing, 2. How is resource allocation used to reward different niching, the Kernel approach, the nearest neighbor approach, the crossover operators? histogram technique, crowding/clustering, and a relaxed form of dominance and restricted mating [5]. Non-dominated sorting algo- The rest of the paper is organized as follows. Section 2 presents rithms (NSGA), such as NSGA-II [9], and SPEA2 [32], are examples the framework of MOEA/D-DRA, and explains how the use of two of Pareto dominance based MOEAs. genetic operators can be implemented by adjusting a probability Decomposition is a basic strategy in traditional mathematical parameter. Section 3 presents MOEA/D-DRA-CMX+SPX. Section 4 programming [11,18]. It has not been widely used in Evolu- describes the benchmark problems, the parameter settings for the tionary Multi-objective Optimization (EMO). Fairly recently, a empirical investigation, and the metric used to measure the per- multi-objective evolutionary algorithm based on decomposition formance of the suggested hybrid algorithm. Section 5 discusses (MOEA/D) has been developed [29]. It combines decomposition experimental results. Section 6 discusses the sensitivity of the techniques and evolutionary algorithms for dealing with MOP. It algorithm performance to key parameters, and finally, Section 7 converts a MOP into a number of scalar optimization subproblems concludes the paper and points out some worthwhile issues to which it then optimizes simultaneously [16,29]. The key features of follow up. MOEA/D are that it solves each subproblem using information from its neighborhood which consists of other subproblems. The neigh- 2. The MOEA/D-DRA framework borhood relationships among the subproblems are based on the distances between their corresponding weight vectors. The popu- MOEA/D-DRA builds up on MOEA/D [30], which, itself, builds up lation kept by MOEA/D includes the best solution found so far for on the MOEA framework which is perhaps epitomized by NSGA-II each subproblem. [9]. In the following we provide the different components required MOEA-type algorithms use mainly two genetic operators, by the suggested hybrid algorithm for a successful implementation. crossover and mutation, for creating offspring. Recall that mutation is applied to a single solution, unlike crossover which is generally 2.1. Decomposition: implementation applied to at least two [1]. A number of crossover operators can be found in the literature MOEA/D implements decomposition to transform the approxi- such as the blend crossover operator (BLX-˛) [12], the simulated mation of the PF of the overall problem by that of N single objective binary crossover (SBX) [15,19], the Simplex Crossover (SPX) [2,27], subproblems. Decomposition can be implemented in many ways. the Center of Mass Crossover (CMX) [25,26], the unimodal normally The approach adopted in this paper is Tchebycheff aggregation distributed crossover (UNDX) [21,22], the parent-centric crossover [18,29,30]. (PCX) [8], and many other real coded crossover operators for deal- ing with continuous problems [7,23]. Here, we consider the use of 2.1.1. The Tchebycheff aggregation function multiple crossover operators within the MOEA/D search algorithm In this approach, the scalar optimization problem is of the form [30]. These operators are SPX [2,27], and CMX [25,26]. The result minimize g te (x|, z) = max {j |fj (x) − zj∗ |}, (2) is a hybrid algorithm denoted MOEA/D-DRA-CMX+SPX in which 1≤j≤m the probability of the use of each crossover operator is automati- m cally adjusted based on its performance. A reward or credit of 1 is where x ∈ , = (1 , 2 , . . ., m )T , j ≥ 0, j = 1, . . ., m, and j=1 j = 1. assigned whenever the generated offspring solution of the current z∗ = ∗ )T (z1∗ , . . . , zmis the reference point, i.e., zj∗ = min{fj (x)|x ∈ } subproblem successfully replaces at least one solution of a neigh- for each j = 1, . . ., m. It is well known that, under mild conditions, boring subproblem. Otherwise, the reward is 0. Hence the concept for each Pareto optimal solution there exists a weight vector such of dynamic resource allocation (DRA). Note that the resource we that it is the optimal solution to (2) and each optimal solution of are concerned with is really “computing”. Computing facilities can (2) is a Pareto optimal solution to problem (1) [29]. Therefore, a be explicitly or implicitly allocated. Here, the facility is implicitly set of Pareto optimal solutions can be obtained to problem (1) by allocated by allowing an operator to be more likely used rather than optimizing a set of single optimization subproblems defined by the another, thus making it absorb computing power every time it is Tchebycheff approach with different weight vectors. used. The key features of the proposed hybrid algorithm, i.e. MOEA/D- 2.1.2. Weight vector selection in the decomposition process DRA-CMX+SPX, are as follows. A set of N weight vectors, W, is generated according to the fol- lowing criteria as were adopted in [30]. 1. It is adaptive. 2. It maintains two different sets of subproblems in its evo- 1. Uniformly randomly generate 5000 weight vectors for lutionary process. forming the set W1. Set W is initialized as the set W. Khan Mashwani, A. Salhi / Applied Soft Computing 12 (2012) 2765–2780 2767 containing all the weight vectors (1, 0, . . ., 0, 0), (0, 1, . . ., 2.2.2.2. The Simplex Crossover. Given three parents, x(1) , x(2) , x(3) , 0, 0), . . ., (0, 0, . . ., 0, 1). SPX generates three offspring as follows: 2. Find the weight vector in set W1 with the largest distance to set W; add it to set W and delete it from set W1 . y(k) = (1 + )(x(k) − o), k = 1, 2, 3. (7) 3. If the size of set W is N, stop and return set W. Otherwise, 3 go to 2. where o = (1/3) k=1 x(k) is the centroid of the parents and ≥ 0 is the scaling parameter which controls the expansion of the simplex. 2.2. Genetic operators In our implementation, we generate one offspring as in [25]: Two main operators are implemented in evolutionary algo- 3 rithms. They are mutation and crossover. The literature of EA also z= ˛i yi + o, refers to a third one, reproduction which is the process of copying i=1 individuals from the current population into the next one with- 3 out alteration. Here, however, reproduction refers to the concept of where ˛i > 0 are randomly generated which satisfy i=1 ˛i = 1. generating new offspring. Following are some key features of SPX [27]: 2.2.1. The polynomial mutation operator 1. SPX is simple and easy to implement. The polynomial mutation operator of [7] is used in the imple- 2. It has only one control parameter, . mentation of MOEA/D-DRA. Let yk be the decision variable to be 3. SPX is based on simplex. It can generate offsprings inside mutated. Then yk is computed as follows. the expanded region of the simplex. 4. SPX is independent of the coordinate system. yk + k (uk − lk ) with probability pm , yk = (3) yk with probability 1 − pm , 2.3. Resource allocation: a probabilistic approach where lk and uk are the lower and upper bounds of the kth decision variable, respectively. And Let p1t and p2t be the probabilities of the use of crossovers C1 and C2 at generation t. These probabilities are updated as follows. (2 × rand)(1/(+1)) − 1 if rand < 0.5 k = 1 − (2 − 2 × rand)(1/(+1)) otherwise 1. p11 and p21 are initialized to 0.5 in generation 1. 2. In Step 3 of the algorithm (the details are given in the where rand ∈ [0, 1], is a uniform random number. The mutation rate next section), in p1t × 100% of the subproblems, C1 will pm and the distribution index are two control parameters. be used for generating new solutions. The rest of the sub- problems will use C2 . 2.2.2. The crossover operators 3. Compute the total reward, r1 and r2 of each crossover As mentioned earlier, there are many crossover operators. Here operator at the current generation. Then set: we are concerned with two of them, namely the Center of Mass Crossover (CMX) [25,26], and Simplex Crossover (SPX) [2,27]. ri pit+1 = 0.5 × pit + 0.5 × r1 + r2 2.2.2.1. The Center of Mass Crossover. Given three parents, x(1) , x(2) , x(3) , the CMX operator works as follows: for i = 1, 2. • Compute the center of mass as A crossover is said to be successful if its new generated solution can replace at least one solution. Whenever a crossover operator is 1 3 o= xi (4) successful it will get a reward of 1. Otherwise, it will get a reward 3 of 0. This means, the more successful crossover operator will be i=1 applied to more subproblems than the less successful one. and then create a set of virtual mates by mirroring each parent Let 1 , . . ., N be a set of even spread weight vectors and z∗ be across the center of mass the reference point. The problem of approximating the PF of (1) can be decomposed into N scalar optimization subproblems and vi = 2o − xi . (5) the objective function of the ith subproblem is In our implementation, we randomly select a virtual mate vj and a parent xi and then generate one offspring as: g te (x|i , z ∗ ) = max {ij |fj (x) − zj∗ |}, 1≤j≤m y = (1 − ˛)xi + ˛vj , (6) where i = (i1 , . . . , im )T , and i = 1, . . ., N. where ˛ = 2r − 0.5, r is a random number. MOEA/D-DRA minimizes all of these N objective functions simultaneously, in a single run. In [17,29], all the subproblems are CMX is expected to have the following characteristics: treated equally; each one of them receives the same amount of computational time. However, it is reasonable to assign different 1. It is the natural generalization of the 2-parent crossover computational efforts to different problems. Indeed, MOEA/D-DRA, operator. as in [30], does just that for different subproblems. It defines and 2. It is simple and easy to implement. computes a utility i for each subproblem i. Computational efforts 3. It has one control parameter which changes dynamically. are distributed to these subproblems based on their utilities. Dur- 4. It tends to generate offspring uniformly around the par- ing the search, MOEA/D-DRA with the Tchebycheff approach (2) ents. maintains: 2768 W. Khan Mashwani, A. Salhi / Applied Soft Computing 12 (2012) 2765–2780 • a population of N points x1 , . . ., xN ∈ , where xi is the current Table 1 Details of the used benchmark functions. solution to the ith subproblem; • FV1 , . . ., FVN , where FVi is the F-value of xi , i.e. FVi = F(xi ) for each Function Search space range Characteristics of PF i = 1, . . ., N; UF1 [0, 1] × [−1, 1]n−1 Concave • z = (z1 , . . ., zm )T , where zi is the best (lowest) value found so far UF2 [0, 1] × [−1, 1]n−1 Concave for each objective fi ; UF3 [0, 1]n Concave • utilities 1 , . . ., N for subproblem i, i = 1, . . ., N. UF4 [0, 1] × [−2, 2]n−1 Convex UF5 [0, 1] × [−1, 1]n−1 21 point front • gen, the current generation number. UF6 [0, 1] × [−1, 1]n−1 One isolated point and two disconnected parts UF7 [0, 1] × [−1, 1]n−1 Continuous straight line 3. A new hybrid algorithm: MOEA/D-DRA-CMX+SPX UF8 [0, 1]2 × [−2, 2]n−2 Parabolic UF9 [0, 1]2 × [−2, 2]n−2 Planar UF10 [0, 1]2 × [−2, 2]n−2 Parabolic Here we describe the hybrid algorithm MOEA/D-DRA-CMX+SPX and how it is implemented. Note that the general framework, as already mentioned earlier, is that of MOEA/D which can be found Step 4 Stopping Criteria: If the stopping criteria are satisfied, in many sources such as [17,29]. The core of the contribution of this then stop and output {x(1) , . . ., x(N) } and {F(x(1) ), . . ., F(x(N) )}. paper is in Step 3.2 where new offspring are generated according Step 5 gen = gen + 1. to either CMX or SPX. Initially, they have the same probability of If gen is a multiple of 50, then being used. Thereafter, their use is dependent on how successful Compute i , the relative decrease in the objective value of each they are at generating good offspring that are used to replace others subproblem i during the last 50 generations, as previously generated. The MOEA/D framework is repeated here for completeness and the benefit of the novice in MOEA matters. old function value − new function value ; old function value Algorithm Update utility i using Input: ⎧ 1: MOP (1); ⎨1 if i > 0.001 2: a stopping criterion; i = i (8) 3: N, the number of subproblems to be considered in MOEA/D- ⎩ (0.95 + 0.05 )i otherwise. DRA; 0.001 4: a uniform spread of N weight vectors 1 , . . ., N ; Endif 5: T, the number of weight vectors in the neighborhood of each Go to Step 2. weight vector. Output: {x(1) , . . ., xN } and {F(x(1) ), . . ., F(xN )} In the 10-tournament selection of Step 2, the index with highest Step 1 Initialization i value from 10 uniformly randomly selected indices is chosen to Step 1.1 Compute the Euclidean distances between any two enter I. We should do this selection [N/5] − m times. weight vectors and then find the T closest weight vectors to each In Step 3.2, we used two crossover operators for generating weight vector. For i = 1, . . ., N, set B(i) = {i1 , . . ., iT } where i1 , . . ., new solutions to the subproblems. The probability of the use of iT are the T closest weight vectors to i . each crossover is updated according to their individual successes. Step 1.2 Generate an initial population x1 , . . ., xN by uniformly Each new individual solution undergoes the polynomial mutation randomly sampling in the search space. [7] operation which is formulated in Eq. (3). The probabilities of Step 1.3 Initialize z = (z1 , . . ., zm )T by setting zi = min{fi (x(1) ), fi (x(2) ), the crossover operators in MOEA/D-DRA-CMX+SPX are updated as . . ., fi (x(N) )}, for i = 1, . . ., m. explained in Section 2.3. Step 1.4 Set gen = 0 and i = 1 for all i = 1, . . ., N. Step 2 Selection of Subproblems for Search: the indices of the 4. Empirical investigation: benchmark, settings and subproblems whose objectives are MOP individual objectives fi performance measure are selected to form initial set I. By using 10-tournament selection based on i , select other [N/5] − m indices and add them to I. 4.1. The benchmark problems: the CEC’09 test set Step 3 For each i ∈ I, do: Step 3.1 Selection of Mating/Update Range: Uniformly randomly Three versions of the new hybrid algorithm based on MOEA/D a number rand from [0, 1]. Then set generate enhanced with resource allocation and one of two crossovers, i.e. B(i) if rand < ı, MOEA/D-DRA-CMX and MOEA/D-DRA-SPX, and MOAE/D enhanced P= {1, . . . , N} otherwise. with resource allocation and both crossovers, i.e. MOEA/D-DRA- Step 3.2 Reproduction: Set r1 = i and randomly select two indices CMX+SPX, are tested on CEC’09 benchmark problems [31]. These r2 and r3 from the current population P. Then apply SPX or CMX problems, their search space ranges and the characteristics of their to generate a new solution y using the selected parent solution, Pareto fronts are given in Table 1. xr1 , xr2 and xr3 . Then apply the polynomial mutation operator [7], to y with probability of mutation pm to get a new solution y. 4.2. Parameter settings Step 3.3 Repair: If the elements of y are not within the specified domain in the parameter space , reset them inside this domain. The following parameters settings were used in our experi- Step 3.4 Update z: For each j = 1, . . ., m, if zj > fj (y), then set zj = fj (y). ments. Step 3.5 Update Solutions: Set c = 0 and then do the following: (1) If c = nr or P is empty, go to Step 4. Otherwise, randomly pick • N = 600, the population size for 2-objective test instances an index j from P. (UF1–UF7), (2) If g(y|j , z) ≤ g(xj |j , z), then set xj = y, FVj = F(y) and c = c + 1. • N = 1000, the population size for 3-objective test instances (3) Delete j from P and go to (1). (UF8–UF10), W. Khan Mashwani, A. Salhi / Applied Soft Computing 12 (2012) 2765–2780 2769 Table 2 The IGD-metric values of the three algorithms over 30 independent runs on UF1, . . ., UF5. CEC’09 Min Median Mean St.dev. max Algorithm Ranks MOEA/D-DRA-CMX (a), MOEA/D-DRA-SPX (b), MOEA/D-DRA-CMX+SPX (c) UF1 0.004084 0.0046905 0.0051250 0.0034 0.0079980 a 2 0.005077 0.008203 0.008652 0.004079 0.026533 b 3 0.003985 0.004171 0.004292 0.000263 0.005129 c 1 UF2 0.0050340 0.0066460 0.0071868 0.00190157 0.012719 a 3 0.004952 0.005596 0.005717 0.000487 0.007131 b 2 0.005149 0.005472 0.005615 0.000412 0.006778 c 1 UF3 0.0049430 0.0432895 0.04153040 0.0239478 0.0855450 a 3 0.005455 0.022374 0.030286 0.03823 0.091197 b 2 0.004155 0.005313 0.011165 0.013093 0.068412 c 1 UF4 0.0508390 0.0632330 0.0628110 0.00476440 0.0734990 a 2 0.0502 0.0573 0.0573 0.0038 0.0657 b 1 0.055457 0.063524 0.064145 0.004241 0.075361 c 3 UF5 0.1676280 0.3791160 0.363792 0.0813385 0.5085780 a 1 0.2186 0.3889 0.3963 0.0742 0.5704 b 2 0.211058 0.379241 0.418508 0.135554 0.707093 c 3 • T = 0.1N, the number of neighbors for each subproblem, approximation. There are several different ways of computing and • nr = 0.01N, the maximal number of solutions to update success- averaging the distances in the GD performance metrics [9,28]. The fully by a new solution of subproblem, version of the IGD-metric used in this thesis inverts the ϒ version • ı = 0.9, the probability for the mating pool to select solutions from of GD as described in [9]. the neighborhood, • = 20, the distribution index used in real mutation, • pm = 1/n, the probability of mutation, where n = 30 denoted the 5. Experimental results √ of parametric variables, number • = n + 1, the control parameter of the simplex, used in SPX, The three versions of MOEA/D-DRA algorithm, namely, • The algorithm stops after 300,000 function evaluations; this is MOEA/D-DRA-CMX, MOEA/D-DRA-SPX, and MOEA/D-DRA- the stopping criterion for all tested versions of MOEA/D-DRA in CMX+SPX, are tested on UF1 through to UF10 [31]. All were run 30 our experimental work. times independently on each test instance. The same parameter settings are used by all algorithms for a fair comparison. 4.3. Performance metric 5.1. Discussion of IGD-metric statistical results for A performance metric or indicator is used to measure the quality MOEA/D-DRA-CMX+SPX of the approximated set of Pareto optimal solutions. In multi- objective optimization, the quality of the obtained set of solutions Tables 2 and 3 present the best (i.e., minimum), median, can be measured in several ways, for example, the closeness of the mean, standard deviation (std), and worst (i.e., maximum) of the solutions to true Pareto optimal front (i.e., convergence property), IGD-metric values of the final approximation over 30 indepen- the uniformity of the solutions (i.e., diversity preservation) and dent executions of MOEA/D-DRA-CMX, MOEA/D-DRA-SPX, and the dominance relationship between two sets of solutions. Lim- MOEA/D-DRA-CMX+SPX on each test instance. itations and suitability of various performance indicators/metrics The last columns of Tables 2 and 3 show the performance ranks have been discussed in the specialized literature such as [24,33,35]. of the individual algorithms; bold data highlight rank 1 (best), tele- The inverted generational distance (IGD)-metric [35] is used in type data, rank 2 and Italic data, rank 3 (i.e., worst). Based on these assessing the performance of the proposed versions of the MOEA/D- rankings, MOEA/D-DRA-CMX+SPX comes first on nine out of ten DRA [30], in this paper. It measures both the convergence and the test problems. This clearly shows that the incorporation of two spread of the obtained solutions; the smaller the value of the IGD- crossover operators is a good idea. It, indeed, improves the per- metric, the better is the obtained set of Pareto optimal solutions. formance of MOEA/D-DRA, on nearly all problems of the CEC’09 Let P∗ be a set of uniformly distributed points in the objective collection [31]. space along the PF. Let P be an approximate set to the PF, the average distance from P∗ to P is defined [29,31], as 5.2. Discussion of Pareto front graphical results v∈P ∗ d(v, P) D(P ∗ , P) = , Figs. 1–4 plot the final approximate solutions in the objective |P ∗ | space with the best (i.e., smallest) IGD-metric value over 30 inde- where d(v, P) is the minimum Euclidean distance between v and pendent runs of MOEA/D-DRA-SPX (1st panel), MOEA/D-DRA-CMX the points in P. If P∗ is large enough to represent the PF very well, (2nd panel) and MOEA/D-DRA-CMX+SPX (3rd panel), respectively, D(P∗ , P) could measure both the diversity and convergence of P in a on all test problems. These figures clearly show that MOEA/D-DRA- sense. To have a lower value of D(P∗ , P), P must be very close to the CMX+SPX is more effective based on the IGD-metric values found. PF and cannot miss any part of it. It is well known, [14], that no single metric can always provide In our experiments, P∗ is 1000 evenly distributed points in the the necessary information to judge the performance of MOEAs. For PF of 2-objective problems and it is 10,000 for problems with 3 this reason, the Pareto fronts generated by the algorithms are plot- objective functions. ted in the objective space to show their distribution range for a The IGD-metric is the inverted version of the widely used gen- given problem. Figs. 5 and 6 show such plots. These figures clearly erational distance (GD) performance metric [28]. GD represents indicate that MOEA/D-DRA-CMX+SPX has produced better PFs in the average distance of the points in an approximation to the each of 30 runs, in the objective space for all CEC’09 test instances true PF, which cannot effectively measure the diversity of the considered except for UF5 and UF6. These two problems are much 2770 W. Khan Mashwani, A. Salhi / Applied Soft Computing 12 (2012) 2765–2780 Table 3 The IGD-metric values of the 3 algorithms over 30 independent runs on UF6–UF10. CEC’09 Min Median Mean St.dev. Max Algorithm Ranks MOEA/D-DRA-CMX (a), MOEA/D-DRA-SPX (b), MOEA/D-DRA-CMX+SPX (c) UF6 0.0710280 0.44046750 0.3665337 0.12022382 0.60331490 a 2 0.445548 0.507600 0.532296 0.070171 0.670946 b 3 0.056972 0.248898 0.327356 0.185717 0.792910 c 1 UF7 0.0045910 0.006100 0.00762145 0.00526699 0.0317640 a 2 0.005776 0.010898 0.062154 0.152973 0.533897 b 3 0.003971 0.004745 0.006262 0.003307 0.014662 c 1 UF8 0.0559350 0.070072250 0.077114571 0.02153013 0.1310650 a 2 0.078084 0.142313 0.136907 0.053748 0.272108 b 3 0.051800 0.056872 0.057443 0.003366 0.065620 c 1 UF9 0.2075150 0.2875470 0.2807007 0.046701797 0.343310 a 3 0.064034 0.181923 0.159370 0.053077 0.216422 b 2 0.033314 0.144673 0.097693 0.054285 0.151719 c 1 UF10 0.4175320 0.4898737 0.4945660 0.0460605756 0.5991810 a 2 0.445548 0.507600 0.532296 0.070171 0.670946 b 3 0.391496 0.467715 0.462653 0.038698 0.533234 c 1 harder compared to the others; indeed, each algorithm has found that of the other crossover must be low. This follows from their a weak set of final solutions for both of them. respective performances. In other words, the crossover that per- formed well will get the higher probability of use. Plots of the 5.3. Discussion of the evolution of the IGD-metric values allocated probability values of CMX and SPX versus the number of iterations of MOEA/D-DRA-CMX+SPX on each of the CEC’09 test Fig. 7 plots the evolution of the average IGD-metric values found instances show just that. The important conclusion to draw is that over 30 independent runs against the number of generations spent where a crossover is not so good, the other one will get used instead. by MOEA/D-DRA-SPX+CMX on each of the CEC’09 test instances This means that overall, the search process uses the most suitable [31]. These indicates that our proposed hybrid algorithm effectively crossover, therefore, increasing the search effectiveness of the pro- handle all tested problems in terms of reducing the IGD-metric cess. It can also be seen as if the two crossovers work synergistically. value compared to one crossover based algorithm. Note that the combination of well performing crossover oper- Fig. 8 plots the evolution of the IGD-metric values of the best run ators will not necessarily lead to a better performance on all among 30 independent runs against the number of generations for problems. Indeed, MOEA/D-DRA-CMX+SPX is more effective than all tested algorithms, MOEA/D-DRA-CMX, MOEA/D-DRA-SPX and the versions of MOEA/D-DRA with individual crossovers overall, MOEA/D-DRA-CMX+SPX on each of the CEC’09 test instance. These but not on UF5 and UF6. This may be explained by the fact that figures visually show that MOEA/D-DRA-CMX+SPX performs signif- for some problems, a given crossover operator is more effective icantly better in terms of convergence towards the true the PF of all throughout the search. The combination may distract the over- CEC’09 test instances except on UF5 and UF6. These two problems all search from concentrating on the use of that one crossover. are quite hard for all versions of MOEA/D-DRA. There is also the cost of overheads involved in the combined version. 5.4. Discussion of the dynamic allocation of resources to CMX and Another way to explain this shortcoming is perhaps that the PF SPX of UF5 is discontinuous. Note that the combined version still per- formed as well as the SPX-based version and only slightly worse The resource, here computing, is in terms of the probability with than the CMX-based version. This may point, as said above, to the which a given crossover is called during the search process. It is fact that the CMX crossover operator is well adapted to that prob- easy to see that when the probability of using a crossover is high, lem. The same goes for UF6. UF1 Pareto Front UF1 Pareto Front UF1 Pareto Front 1.2 1.2 1.2 Real PF Real PF Real PF MOEAD/DRA−SPX MOEAD/DRA−CMX MOEAD/DRA−CMX+SPX 1.0 1.0 1.0 0.8 0.8 0.8 f2 f2 f2 0,6 0,6 0,6 0.4 0.4 0.4 0.2 0.2 0.2 0.0 0.0 0.0 0.0 0.2 0.4 0,6 0.8 1.0 1.2 0.0 0.2 0.4 0,6 0.8 1.0 1.2 0.0 0.2 0.4 0,6 0.8 1.0 1.2 f1 f1 f1 Fig. 1. Plots of the final non-dominated solutions with the lowest IGD-metric value over 30 independent runs in the objective space of UF1. The first figure is for MOEA/D- DRA-SPX, the second figure is for MOEA/D-DRA-CMX and the third figure is for MOEA/D-DRA-CMX+SPX. W. Khan Mashwani, A. Salhi / Applied Soft Computing 12 (2012) 2765–2780 2771 UF2 Pareto Front UF2 Pareto Front UF2 Pareto Front 1.2 1.2 1.2 Real PF Real PF Real PF MOEAD/DRA−SPX MOEAD/DRA−CMX MOEAD/DRA−CMX+SPX 1.0 1.0 1.0 0.8 0.8 0.8 f2 f2 f2 0,6 0,6 0,6 0.4 0.4 0.4 0.2 0.2 0.2 0.0 0.0 0.0 0.0 0.2 0.4 0,6 0.8 1.0 1.2 0.0 0.2 0.4 0,6 0.8 1.0 1.2 0.0 0.2 0.4 0,6 0.8 1.0 1.2 f1 f1 f1 UF3 Pareto Front UF3 Pareto Front UF3 Pareto Front 1.2 1.2 1.2 Real PF Real PF Real PF MOEAD/DRA−SPX MOEAD/DRA−CMX MOEAD/DRA−CMX+SPX 1.0 1.0 1.0 0.8 0.8 0.8 f2 f2 f2 0,6 0,6 0,6 0.4 0.4 0.4 0.2 0.2 0.2 0.0 0.0 0.0 0.0 0.2 0.4 0,6 0.8 1.0 1.2 0.0 0.2 0.4 0,6 0.8 1.0 1.2 0.0 0.2 0.4 0,6 0.8 1.0 1.2 f1 f1 f1 UF4 Pareto Front UF4 Pareto Front UF4 Pareto Front 1.2 1.2 1.2 Real PF Real PF Real PF MOEAD/DRA−SPX MOEAD/DRA−CMX MOEAD/DRA−CMX+SPX 1.0 1.0 1.0 0.8 0.8 0.8 f2 f2 f2 0,6 0,6 0,6 0.4 0.4 0.4 0.2 0.2 0.2 0.0 0.0 0.0 0.0 0.2 0.4 0,6 0.8 1.0 1.2 0.0 0.2 0.4 0,6 0.8 1.0 1.2 0.0 0.2 0.4 0,6 0.8 1.0 1.2 f1 f1 f1 Fig. 2. Plots of the final non-dominated solutions with the lowest IGD-metric value over 30 independent runs in the objective space of UF2, . . ., UF4. The 1st panel is for MOEA/D-DRA-SPX, the 2nd is for MOEA/D-DRA-CMX, and the 3rd is for MOEA/D-DRA-CMX+SPX. 6. Impact of parameter settings on MOEA/D-DRA-CMX+SPX 6.1. Performance of MOEA/D-DRA-CMX+SPX with different population sizes In the following, we report on our experiences with MOEA/D- DRA-CMX+SPX when solving standard test problems with MOEA/D-DRA-CMX+SPX was tested with population sizes different settings of key parameters such as the population N = 200, 300, 400, and 500 on 2-objective problems and with popu- size (N), the neighborhood size (T), and the maximum num- lation sizes N = 250, 500, 600, and 800 on 3-objective problems. The ber of subproblem solutions replaced (nr ) after each search results were compared with those obtained with MOEA/D-DRA- cycle. The impact of these parameter settings on the perfor- CMX and MOEA/D-DRA-SPX. The test problems are from the CEC’09 mance of the suggested hybrid algorithm is in terms of the test suit. calculated IGD-metric values from the corresponding search Tables 4 and 5 show the IGD-metric values including best, results. median, mean, standard deviation and maximum obtained for 2772 W. Khan Mashwani, A. Salhi / Applied Soft Computing 12 (2012) 2765–2780 UF5 Pareto Front UF5 Pareto Front UF5 Pareto Front 1.2 1.2 1.2 Real PF Real PF Real PF MOEAD/DRA−SPX MOEAD/DRA−CMX MOEAD/DRA−CMX+SPX 1.0 1.0 1.0 0.8 0.8 0.8 f2 f2 f2 0,6 0,6 0,6 0.4 0.4 0.4 0.2 0.2 0.2 0.0 0.0 0.0 0.0 0.2 0.4 0,6 0.8 1.0 1.2 0.0 0.2 0.4 0,6 0.8 1.0 1.2 0.0 0.2 0.4 0,6 0.8 1.0 1.2 f1 f1 f1 UF6 Pareto Front UF6 Pareto Front UF6 Pareto Front 1.2 1.2 1.2 Real PF Real PF Real PF MOEAD/DRA−SPX MOEAD/DRA−CMX MOEAD/DRA−CMX+SPX 1.0 1.0 1.0 0.8 0.8 0.8 f2 f2 f2 0,6 0,6 0,6 0.4 0.4 0.4 0.2 0.2 0.2 0.0 0.0 0.0 0.0 0.2 0.4 0,6 0.8 1.0 1.2 0.0 0.2 0.4 0,6 0.8 1.0 1.2 0.0 0.2 0.4 0,6 0.8 1.0 1.2 f1 f1 f1 UF7 Pareto Front UF7 Pareto Front UF7 Pareto Front 1.2 1.2 1.2 Real PF Real PF Real PF MOEAD/DRA−SPX MOEAD/DRA−CMX MOEAD/DRA−CMX+SPX 1.0 1.0 1.0 0.8 0.8 0.8 f2 f2 f2 0,6 0,6 0,6 0.4 0.4 0.4 0.2 0.2 0.2 0.0 0.0 0.0 0.0 0.2 0.4 0,6 0.8 1.0 1.2 0.0 0.2 0.4 0,6 0.8 1.0 1.2 0.0 0.2 0.4 0,6 0.8 1.0 1.2 f1 f1 f1 Fig. 3. Plots of the final non-dominated solutions with the lowest IGD-metric value over 30 independent runs in the objective space of UF5, . . ., UF7. The 1st panel is for MOEA/D-DRA-SPX, the 2nd is for MOEA/D-DRA-CMX, and the 3rd is for MOEA/D-DRA-CMX+SPX. MOEA/D-DRA-CMX+SPX with different population settings (given of the same collection. Note that, apart from N, the pop- in the last columns of these tables), when applied to each of the ulation size, the rest of parameters is the same as set in CEC’09 test instances [31]. Section 4.2. Table 7 presents the IGD-metric statistics (i.e., best, median, Tables 4 and 5 clearly indicate that MOEA/D-DRA-CMX+SPX mean, standard deviation and maximum) for MOEA/D-DRA-CMX is not very sensitive to population size compared to MOEA/D- with two different population sizes N = 400, 500, when applied to DRA-CMX and MOEA/D-DRA-SPX. Furthermore, Tables 4 and 5 2-objective problems, and N = 500, 600 when applied to 3-objective show that the performance of MOEA/D-DRA-CMX+SPX improves problems in the CEC’09 collection. Table 6 records the statis- with increasing population size unlike that of MOEA/D-DRA- tics of MOEA/D-DRA-SPX with two different population sizes, CMX and MOEA/D-DRA-SPX, in terms of reducing the average 400 and 500, when applied to 2-objective problems and IGD-metric values on seven (i.e., UF1–UF7) out of the ten test sizes 500, and 600, when applied to 3-objective problems instances. W. Khan Mashwani, A. Salhi / Applied Soft Computing 12 (2012) 2765–2780 2773 UF8 Pareto Front UF8 Pareto Front UF8 Pareto Front Real PF Real PF Real PF MOEAD/DRA−SPX MOEAD/DRA−CMX MOEAD/DRA−CMX+SPX 1.2 1.2 1.2 1.0 1.0 1.0 0.8 0.8 0.8 0,6 0,6 0,6 0.4 0.4 0.4 0.2 0.2 0.2 0.0 0.0 0.0 1.2 1.2 1.2 1.2 1.0 1.2 1.0 1.0 1.0 1.0 0.8 1.0 0.8 0.8 0.8 0.8 0,6 0.8 0,6 0,6 0,6 0,6 0 0.4 0 ,6 0.4 0.2 0.2 0.4 0.4 0.2 0.2 0.4 0.2 0.2 .4 0.0 0.0 0.0 0.0 0.0 0.0 f2 f1 f2 f1 f2 f1 UF9 Pareto Front UF9 Pareto Front UF9 Pareto Front Real PF Real PF Real PF MOEAD/DRA−SPX MOEAD/DRA−CMX MOEAD/DRA−CMX+SPX 1.2 1.2 1.0 1.0 1.2 0.8 0.8 1.0 0,6 0.8 0,6 0.4 0,6 0.4 0.4 0.2 0.2 0.2 0.0 0.0 0.0 1.2 1.2 1.2 1.2 1.2 1.2 1.0 1.0 1.0 1.0 1.0 1.0 0.8 0.8 0.8 0.8 0.8 0.8 0,6 0,6 0,6 0,6 0,6 0,6 0.4 0.4 0.4 0.4 0.4 0.4 0.2 0.2 0.2 0.2 0.2 0.2 0.0 0.0 0.0 0.0 0.0 0.0 f2 f1 f2 f1 f2 f1 UF10 Pareto Front UF10 Pareto Front UF10 Pareto Front Real PF Real PF Real PF MOEAD/DRA−SPX MOEAD/DRA−CMX MOEAD/DRA−CMX+SPX 1.2 1.0 0.8 1.2 1.0 0,6 0.8 1.2 1.0 0.4 0,6 0.8 0.4 0,6 0.2 0.4 0.2 0.2 0.0 0.0 0.0 1.2 1.2 1.2 1.2 1.2 1.2 1.0 1.0 1.0 1.0 1.0 1.0 0.8 0.8 0.8 0.8 0.8 0.8 0,6 0,6 0,6 0,6 0,6 0,6 0.4 0.4 0.4 0.4 0.4 0.4 0.2 0.2 0.2 0.2 0.2 0.2 0.0 0.0 0.0 0.0 0.0 0.0 f2 f1 f2 f1 f2 f1 Fig. 4. Plots of the final non-dominated solutions with the lowest IGD-metric value of the three algorithms over 30 independent runs in the objective space of UF8, . . ., UF10. The 1st panel is for MOEA/D-DRA-SPX, the 2nd is for MOEA/D-DRA-CMX, and the 3rd is for MOEA/D-DRA-CMX+SPX. 6.2. Impact of the neighborhood size on MOEA/D-DRA-CMX+SPX parameter in MOEA/D and its setting should be smaller than N as recommended in [29]. The neighborhood concept is one of the key features of the In this section, we study the sensitivity of MOEA/D-DRA- MOEA/D paradigm. Its size T is the most important control CMX+SPX to parameter T. All other parameters are kept the same as 2774 W. Khan Mashwani, A. Salhi / Applied Soft Computing 12 (2012) 2765–2780 Fig. 5. Plots of all 30 Pareto fronts in the objective space found by MOEA/D-DRA-SPX for UF1, . . ., UF4. The 1st panel is for MOEA/D-DRA-SPX, the 2nd is for MOEA/D-DRA-CMX, and the 3rd is for MOEA/D-DRA-CMX+SPX. W. Khan Mashwani, A. Salhi / Applied Soft Computing 12 (2012) 2765–2780 2775 UF7 Pareto Front UF7 Pareto Front UF7 Pareto Front 1.2 1.2 1.2 Real PF Real PF Real PF MOEA/D−DRA−SPX−30Runs MOEA/D−DRA−CMX−30Runs MOEA/D−DRA−CMX+SPX−30Runs 1.0 1.0 1.0 0.8 0.8 0.8 f2 f2 f2 0,6 0,6 0,6 0.4 0.4 0.4 0.2 0.2 0.2 0.0 0.0 0.0 0.0 0.2 0.4 0,6 0.8 1.0 1.2 0.0 0.2 0.4 0,6 0.8 1.0 1.2 0.0 0.2 0.4 0,6 0.8 1.0 1.2 f1 f1 f1 UF8 Pareto Front UF8 Pareto Front UF8 Pareto Front Real PF Real PF MOEA/D−DRA−CMX−30Runs MOEA/D−DRA−CMX+SPX−30Runs 1.2 1.2 1.2 1.0 1.0 1.0 0.8 0.8 0.8 0,6 0,6 0,6 0.4 0.4 0.4 0.2 0.2 0.2 0.0 0.0 0.0 1.2 1.2 1.2 1.2 1.0 1.0 1.0 1.0 1.2 1.2 0.8 0.8 0.8 0.8 1.0 1.0 0,6 0,6 0,6 0.8 0.8 0.4 0.4 0,6 0,6 0,6 0.2 0.2 0.4 0.4 0.4 0.4 0.0 0.0 0.2 0.2 0.2 0.2 0.0 0.0 0.0 0.0 f2 f1 f2 f1 f2 f1 Fig. 6. Plots of all 30 Pareto fronts in the objective space found by MOEA/D-DRA-SPX+CMX for UF7, UF8. The 1st panel is for MOEA/D-DRA-SPX, the 2nd is for MOEA/D-DRA- CMX, and the 3rd is for MOEA/D-DRA-CMX+SPX. set in Section 4.2, except T, of course. Fig. 9(center) plots the evolu- Note that we have only included an exemplar graph of the results tion of the average IGD-metric value against different settings of T obtained for UF1. in MOEA/D-DRA-CMX+SPX. This figure clearly shows that this algo- The IGD-metric values in Tables 4 and 5 are obtained with rithm performs very well with T in the range 50–100. It also works MOEA/D-DRA-CMX+SPX with different population sizes (i.e., as very well for all values of T except the very small ones and T = N. explained in Section 6) and with different neighborhood sizes, Thus, we can say that MOEA/D-DRA-CMX+SPX is not sensitive to T = 20, 30, 40, 50 for 2-objective and T = 25, 50, 60, 80 for 3-objective T on at least the CEC’09 problems that we have considered here. CEC’09 test instances [31], respectively. Table 4 The IGD-metric values of MOEA/D-DRA-CMX+SPX over 30 independent runs with different population sizes N on the 2-objective CEC’09 test instances [31]. CEC’09 Min Median Mean St.dev. Max N UF1 0.012823 0.026487 0.041614 0.038193 0.153849 200 0.004129 0.005037 0.006873 0.007852 0.046767 300 0.005343 0.009221 0.026405 0.045996 0.172712 400 0.004223 0.005579 0.006141 0.002063 0.013711 500 UF2 0.008142 0.010616 0.011191 0.002201 0.016848 200 0.004974 0.005656 0.005798 0.000633 0.007276 300 0.006223 0.007343 0.007427 0.000985 0.011199 400 0.005426 0.006348 0.006437 0.000676 0.008557 500 UF3 0.047194 0.112486 0.141920 0.081489 0.325030 200 0.004188 0.009407 0.013399 0.015863 0.086123 300 0.005075 0.043814 0.047143 0.030288 0.112878 400 0.004654 0.010830 0.011821 0.006879 0.034609 500 UF4 0.066841 0.080706 0.081409 0.008061 0.105816 200 0.004188 0.009407 0.013399 0.015863 0.086123 300 0.057589 0.070794 0.071627 0.005536 0.083088 400 0.050508 0.056913 0.057854 0.005591 0.071978 500 UF5 0.368941 0.606737 0.623653 0.138410 0.882580 200 0.214192 0.381650 0.399385 0.090454 0.707107 300 0.329791 0.417681 0.432884 0.083484 0.631771 400 0.223677 0.373221 0.421708 0.154778 0.707106 500 2776 W. Khan Mashwani, A. Salhi / Applied Soft Computing 12 (2012) 2765–2780 UF1 UF2 UF3 1.4 0.35 0.8 Average IGD−metric Values Variation in 30 Runs Average IGD−metric Values Variation in 30 Runs MOEA/D−DRA−CMX MOEA/D−DRA−CMX Average IGD−metric Values Variation in 30 Runs MOEA/D−DRA−CMX MOEA/D−DRA−SPX MOEA/D−DRA−SPX MOEA/D−DRA−SPX 1.2 0.7 MOEA/D−DRA−CMX+SPX 0.3 MOEA/D−DRA−CMX+SPX MOEA/D−DRA−CMX+SPX 0.6 1 0.25 0.5 0.8 0.2 0.4 0.6 0.15 0.3 0.4 0.1 0.2 0.2 0.05 0.1 0 0 0 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 Function Evaluations Function Evaluations Function Evaluations UF4 UF5 UF6 0.2 5 5 Average IGD−metric Values Variation in 30 Runs Average IGD−metric Values Variation in 30 Runs MOEA/D−DRA−CMX MOEA/D−DRA−CMX Average IGD−metric Values Variation in 30 Runs MOEA/D−DRA−CMX MOEA/D−DRA−SPX 4.5 MOEA/D−DRA−SPX 4.5 MOEA/D−DRA−SPX 0.18 MOEA/D−DRA−CMX+SPX MOEA/D−DRA−CMX+SPX MOEA/D−DRA−CMX+SPX 4 4 0.16 3.5 3.5 3 3 0.14 2.5 2.5 0.12 2 2 0.1 1.5 1.5 1 1 0.08 0.5 0.5 0.06 0 0 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 Function Evaluations Function Evaluations Function Evaluations UF7 UF8 UF9 1.4 1.4 1.4 Average IGD−metric Values Variation in 30 Runs Average IGD−metric Values Variation in 30 Runs Average IGD−metric Values Variation in 30 Runs MOEA/D−DRA−CMX MOEA/D−DRA−CMX MOEA/D−DRA−CMX MOEA/D−DRA−SPX MOEA/D−DRA−SPX MOEA/D−DRA−SPX 1.2 MOEA/D−DRA−CMX+SPX 1.2 MOEA/D−DRA−CMX+SPX 1.2 MOEA/D−DRA−CMX+SPX 1 1 1 0.8 0.8 0.8 0.6 0.6 0.6 0.4 0.4 0.4 0.2 0.2 0.2 0 0 0 0 500 1000 1500 2000 2500 0 500 1000 1500 0 500 1000 1500 Function Evaluations Function Evaluations Function Evaluations UF10 9 Average IGD−metric Values Variation in 30 Runs MOEA/D−DRA−CMX 8 MOEA/D−DRA−SPX MOEA/D−DRA−CMX+SPX 7 6 5 4 3 2 1 0 0 500 1000 1500 Function Evaluations Fig. 7. Evolution of the average IGD-metric values versus the number of generations spent by MOEA/D-DRA-CMX, MOEA/D-DRA-SPX, and MOEA/D-DRA-CMX+SPX on UF1, . . ., UF6. The left of Fig. 9 shows the evolution of the best run IGD-metric D-DRA-CMX+SPX is worst when T = 600 in terms of the increase in values over 30 runs on the CEC’09 instances. To address the issue the IGD-metric value (here displayed for instance UF1) compared of the impact on the performance of MOEA/D-DRA-CMX+SPX, we to other T values. These results indicate that T should be set less set T = N. The center of Fig. 9 shows that the performance of MOEA/ than the population size for better performance on the IGD-metric. W. Khan Mashwani, A. Salhi / Applied Soft Computing 12 (2012) 2765–2780 2777 UF1 UF2 UF3 Variation in IGD Values in best run 1.4 0.35 0.9 Variation in IGD Values in best run MOEA/D−DRA−CMX MOEA/D−DRA−CMX Variation in IGD Values in best run MOEA/D−DRA−CMX MOEA/D−DRA−SPX MOEA/D−DRA−SPX 0.8 MOEA/D−DRA−SPX 1.2 MOEA/D−DRA−CMX+SPX 0.3 MOEA/D−DRA−CMX+SPX MOEA/D−DRA−CMX+SPX 0.7 1 0.25 0.6 0.8 0.2 0.5 0.6 0.4 0.15 0.3 0.4 0.1 0.2 0.2 0.05 0.1 0 0 0 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 Function Evaluations Function Evaluations Function Evaluations UF5 UF6 UF4 6 7 Variation in IGD Values in best run 0.2 MOEA/D−DRA−CMX MOEA/D−DRA−CMX Variation in IGD Values in best run Variation in IGD Values in best run MOEA/D−DRA−CMX MOEA/D−DRA−SPX MOEA/D−DRA−SPX MOEA/D−DRA−SPX MOEA/D−DRA−CMX+SPX 6 MOEA/D−DRA−CMX+SPX 0.18 5 MOEA/D−DRA−CMX+SPX 0.16 5 4 0.14 4 3 0.12 3 0.1 2 2 0.08 1 1 0.06 0.04 0 0 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 Function Evaluations Function Evaluations Function Evaluations UF7 UF8 UF9 1.4 1.4 1.4 Variation in IGD Values in best run Variation in IGD Values in best run MOEA/D−DRA−CMX MOEA/D−DRA−CMX MOEA/D−DRA−CMX MOEA/D−DRA−SPX MOEA/D−DRA−SPX Variation in IGD Values in best run MOEA/D−DRA−SPX 1.2 MOEA/D−DRA−CMX+SPX 1.2 MOEA/D−DRA−CMX+SPX 1.2 MOEA/D−DRA−CMX+SPX 1 1 1 0.8 0.8 0.8 0.6 0.6 0.6 0.4 0.4 0.4 0.2 0.2 0.2 0 0 0 0 500 1000 1500 2000 2500 0 500 1000 1500 0 500 1000 1500 Function Evaluations Function Evaluations Function Evaluations UF10 9 Variation in IGD Values in best run MOEA/D−DRA−CMX 8 MOEA/D−DRA−SPX MOEA/D−DRA−CMX+SPX 7 6 5 4 3 2 1 0 0 500 1000 1500 Function Evaluations Fig. 8. Evolution of the IGD-metric values versus the number of generations of the best run over 30 independent runs of MOEA/D-DRA-CMX, MOEA/D-DRA-SPX, and MOEA/D-DRA-CMX+SPX on UF1, . . ., UF6. 2778 W. Khan Mashwani, A. Salhi / Applied Soft Computing 12 (2012) 2765–2780 Table 5 The IGD-metric values of MOEA/D-DRA-CMX+SPX over 30 independent runs with different population sizes on the 3-objective CEC’09 test instances [31]. CEC’09 Min Median Mean St.dev. Max N UF6 0.173843 0.443582 0.400633 0.134066 0.759079 200 0.164252 0.428657 0.374936 0.139382 0.722266 300 0.167471 0.444579 0.415866 0.130556 0.716430 400 0.143108 0.442482 0.0382341 0.153729 0.804836 500 UF7 0.009082 0.019022 0.146966 0.183822 0.556820 200 0.004466 0.005955 0.006501 0.002456 0.016846 300 0.005969 0.009837 0.083261 0.162711 0.542298 400 0.005045 0.006124 0.00690 0.001872 0.014013 500 UF8 0.112062 0.135966 0.148222 0.035861 0.244353 250 0.060292 0.098290 0.105947 0.035883 0.217259 500 0.051823 0.067814 0.070392 0.016783 0.110421 600 0.055175 0.0758446 0.0885577 0.029123 0.182638 800 UF9 0.114800 0.202512 0.202146 0.023488 0.260696 250 0.048576 0.162500 0.148970 0.039373 0.189787 500 0.059098 0.152641 0.111974 0.032441 0.118181 600 0.053043 0.079842 0.07828 0.020178 0.15234 800 UF10 0.428427 0.514462 0.518051 0.057482 0.611057 200 0.364067 0.457048 0.459822 0.045823 0.546396 500 0.406185 0.460197 0.465373 0.037577 0.530965 600 0.314962 0.446823 0.451374 0.045025 0.530874 800 Table 6 The IGD-metric values of MOEA/D-DRA-SPX over 30 independent runs with population sizes N used on CEC’09 test instances [31]. CEC’09 Best Median Mean St.dev. Worst N UF1 0.010524 0.024126 0.032884 0.030600 0.167038 300 0.007577 0.013552 0.025604 0.041755 0.210340 400 UF2 0.007600 0.008815 0.008911 0.000829 0.010604 300 0.006783 0.008209 0.008298 0.000699 0.010667 400 UF3 0.012746 0.058641 0.059622 0.030882 0.140573 300 0.007767 0.037966 0.041621 0.026521 0.108513 400 UF4 0.061380 0.067608 0.067832 0.004257 0.077741 300 0.057845 0.063680 0.063797 0.004520 0.074215 400 UF5 0.350164 0.448279 0.459738 0.079113 0.630960 300 0.266307 0.406658 0.421402 0.073858 0.577894 400 UF6 0.174771 0.419950 0.370526 0.096431 0.462408 300 0.185385 0.443906 0.383862 0.097282 0.475299 400 UF7 0.008347 0.191911 0.190371 0.162191 0.526301 300 0.007253 0.022904 0.173042 0.213450 0.574394 400 UF8 0.064486 0.093309 0.091978 0.018175 0.131712 500 0.053873 0.064945 0.071863 0.017748 0.116492 600 UF9 0.048333 0.151325 0.120079 0.047077 0.193034 500 0.040426 0.148211 0.131620 0.039884 0.158120 600 UF10 0.385331 0.493508 0.489788 0.050616 0.574868 500 0.245726 0.506193 0.493731 0.077428 0.635953 600 Table 7 The IGD-metric values of MOEA/D-DRA-CMX over 30 independent runs on CEC’09 test instances [31]. CEC’09 Best Median Mean St.dev. Worst N UF1 0.010524 0.024126 0.032884 0.030600 0.167038 300 0.086596 0.113851 0.115079 0.014754 0.162955 400 UF2 0.022165 0.027052 0.027051 0.001940 0.030585 300 0.020010 0.022568 0.022696 0.001549 0.025608 400 UF3 0.082230 0.119640 0.230061 0.195813 0.594848 300 0.065553 0.114643 0.269749 0.227045 0.595418 400 UF4 0.081325 0.097801 0.097550 0.010116 0.116335 300 0.077429 0.095759 0.094752 0.010616 0.112952 400 UF5 1.889568 2.794576 2.726526 0.406303 3.513249 300 1.249339 2.392959 2.364813 0.401596 3.123953 400 UF6 0.619284 1.364154 1.280091 0.363662 1.862737 300 0.590681 1.148391 1.165530 0.278422 1.799595 400 UF7 0.036360 0.050242 0.065372 0.088558 0.532118 300 0.024297 0.033658 0.047508 0.054138 0.320080 400 UF8 0.095381 0.195370 0.192308 0.045898 0.281929 500 0.103358 0.192468 0.182471 0.041294 0.242652 600 UF9 0.061337 0.169943 0.143623 0.049655 0.193465 500 0.055389 0.172894 0.153690 0.046937 0.199654 600 UF10 0.184626 0.196215 0.320123 0.215095 0.787695 500 0.178995 0.197305 0.258797 0.150136 0.778593 600 W. Khan Mashwani, A. Salhi / Applied Soft Computing 12 (2012) 2765–2780 2779 UF1 UF1 −3 x 10 UF1 0.03 0.04 7 MOEA/D−DRA−CMX+SPX MOEA/D−DRA−CMX+SPX MOEA/D−DRA−CMX+SPX The average IGD−metric Values 0.035 6 The Average IGD−metric Values 0.025 The IGD−metric Values 0.03 5 0.02 0.025 4 0.015 0.02 3 0.015 0.01 2 0.01 0.005 0.005 1 0 0 0 0 100 200 300 400 500 600 0 100 200 300 400 500 600 5 10 15 20 T: Neighbourhood size T: Neighbourhood size The maximal mumber of solutions to replace Fig. 9. (Left) Plot of IGD-metric values against T values in MOEA/D-DRA-CMX+SPX on UF1; (center) plot of average IGD-metric values against T values in MOEA/D-DRA- CMX+SPX on UF1; (right) plot of average IGD-metric values against values of nr in MOEA/D-DRA-CMX+SPX for UF1. 6.3. Impact of the maximal replacement number restriction on shapes [3]. Again, we are looking at carrying out such experiments MOEA/D-DRA-CMX+SPX in the future. Furthermore, it is perhaps worthwhile considering a generaliza- We have also investigated the sensitivity of MOEA/D-DRA- tion of the multi-crossover approach described here, to any number CMX+SPX to parameter nr , i.e. the maximal number of solutions of crossover operators when no obvious one can be chosen. to be replaced in the neighborhood of each subproblem, to main- tain diversity in the current population. Here we only report an Acknowledgments exemplar outcome of this investigation in the case problem UF1. Fig. 9(right) plots the average IGD-metric values of MOEA/D- We are grateful to the referees for their useful comments and to DRA-CMX+SPX on problem UF1, against nr values 3, 4, 5, 6, 7, 8, 9, Professor Qingfu Zhang for his help. 10, and 20. All other parameters are kept the same as in Section 4.2. It is suggested in [30] that nr should be set smaller than T in order to maintain good diversity in the population. The experi- References ments were carried out with this suggestion in mind, as can be seen in Fig. 9(right). The figure also clearly indicates that the IGD- [1] J. Branke, K. Deb, K. Miettinen, R. Slowinski (Eds.), Multiobjective Optimization: Interactive and Evolutionary Approaches, Springer, 2008. metric values of MOEA/D-DRA-CMX+SPX are decreasing from 8 to [2] Z. Cai, Y. Wang, A multiobjective optimization-based evolutionary algorithm 20 on UF1. for constrained optimization, IEEE Transactions on Evolutionary Computation Note that, since nr = 0.01N, the IGD-metric statistics reported in 10 (December (6)) (2006) 658–675. [3] C.K. Chow, S.Y. Yuen, A multiobjective evolutionary algorithm that diversifies Table 4 are obtained with nr = 2, 3, 4, 5 for each of the 2-objective population by its density, IEEE Transactions on Evolutionary Computation (99) test problems and those reported in Table 5 are with nr = 3, 5, 6, 8 (2011) 1–24. for the 3-objective ones. [4] C.A.C. Coello, A comprehensive survey of evolutionary-based multiobjec- tive optimization techniques, Knowledge and Information Systems 1 (1999) 269–308. [5] C.A. Coello Coello, G.B. Lamont, D.A. Veldhuizen, Evolutionary Algorithms for 7. Conclusion and further work Solving Multi-Objective Problems, Kluwer Academic Publishers, New York, March 2002. [6] K. Deb, Multi-Objective Optimization Using Evolutionary Algorithms, 1st ed., Different crossover operators suit different problems. Using John Wiley and Sons, Chichester, UK, 2001. multiple crossover operators should be an effective strategy for [7] K. Deb, in: S. Ross, R. Weber (Eds.), Multi-Objective Optimization Using Evolu- avoiding the use of potentially unsuitable crossovers. In this paper, tionary Algorithms, 2nd ed., John Wiley and Sons Ltd., 2002. [8] K. Deb, A. Anand, D. Joshi, A computationally efficient evolutionary algo- we have studied the combined use of two crossover operators in rithm for real-parameter optimization, Evolutionary Computation 10 (2002) MOEA/D-DRA [30], and developed a hybrid version of it, referred to 371–395. as MOEA/D-DRA-CMX+SPX using a probabilistic approach to decide [9] K. Deb, A. Pratap, S. Agarwal, T. Meyarivan, A fast and elitist multiobjective genetic algorithm: NSGA-II, IEEE Transactions on Evolutionary Computation 6 which crossover is to be used in a given generation. The latter aspect (2) (2002) 182–197. is what we refer to throughout as resource allocation. The resource [10] F.Y. Edgeworth, Mathematical Psychics, P. Keagan, London, England, 1881. allocation is in terms of the frequency with which a given crossover [11] M. Ehrgott, Multicriteria Optimization, 2nd ed., Springer, Berlin, 2005. [12] L.J. Eshelman, J.D. Schaffer, Real-coded genetic algorithms and interval- is used. This is implemented through a reward system which is just schemata, in: Foundations of Genetic Algorithms, 1993, pp. 187–202. the probability of use of each crossover. There are many possible [13] J. Horn, N. Nafpliotis, D.E. Goldberg, A niched Pareto genetic algorithm for ways to apportion probabilities to crossovers in addition to the one multiobjective optimization, in: Proceedings of the First IEEE Conference on Evolutionary Computation, CEC’94, 1994. used here. We have already carried out some experiments with [14] J. Knowles, D. Corne, Properties of an adaptive archiving algorithm for stor- an alternative approach. The results will be the subject of a future ing nondominated vectors, IEEE Transactions on Evolutionary Computation 7 paper. (April (2)) (2003) 100–116. [15] K. Deb, R. Agrawal, Simulated binary crossover for continuous search space, Results show that, based on IGD-metric values, MOEA/D-DRA- Complex System 9 (1995) 115–148. CMX+SPX is more effective on most of the CEC’09 test instances [16] H. Li, Combination of Evolutionary Algorithms with Decomposition Techniques [31], than either versions of MOEA/D-DRA using CMX or SPX as the for Multiobjective Optimization, Ph.D., Department of Computer Science, Uni- crossover operator. These test problems are more difficult than the versity of Essex, Wevehoe Park, Colchester, Essex, UK, 2007. [17] H. Li, Q. Zhang, Multiobjective optimization problems with complicated pareto ones described in [17], for instance. However, it is still sensible to sets: MOEA/D and NSGA-II, IEEE Transactions on Evolutionary Computation 13 carry out tests on these problems and others with complicated PS (April (2)) (2009) 284–302. 2780 W. Khan Mashwani, A. Salhi / Applied Soft Computing 12 (2012) 2765–2780 [18] K.M. Miettinien, Nonlinear Multiobjective Optimization, Academic Publishers [27] S. Tsutsui, M. Yamamura, T. Higuchi, Multi-parent recombination with simplex Kluwer, Norwell, MA, 1999 (ser. Kluwer’s International Series). crossover in real coded genetic algorithms, in: Proceeding of GECCO-99, 1999, [19] M.M. Raghuwanshi, P.M. Singru, U. Kale, O.G. Kakde, Simulated binary pp. 657–674. crossover with logonormal distribution, Complexity International 12 (2005) [28] D.A.V. Veldhuizen, G.B. Lamont, Evolutionary Computation and Convergence 1–10. to a Pareto Front, Morgan Kaufmann, 1998, Stanford University, California, pp. [20] N. Srinivas, K. Deb, A multiobjective optimization using nondominated sort- 221–228. ing in genetic algorithms, Journal of Evolutionary Computation 2 (3) (1994) [29] Q. Zhang, H. Li, MOEA/D: a multiobjective evolutionary algorithm based on 221–248. decomposition, IEEE Transactions on Evolutionary Computation 11 (December [21] I. Ono, S. Kobayashi, A real-coded genetic algorithm for function optimiza- (6)) (2007) 712–731. tion using the unimodal normal distribution crossover, in: 7th International [30] Q. Zhang, W. Liu, H. Li, The performance of a new version of MOEA/D on CEC’09 Conference on Genetic Algorithms, 1997, pp. 246–253. unconstrained MOP test instances, in: IEEE Congress on Evolutionary Compu- [22] I. Ono, H. Kita, S. Kobayashi, A Real-Coded Genetic Algorithm Using the Uni- tation (IEEE CEC 2009), Trondheim, Norway, May 18–21, 2009, pp. 203–208. modal Normal Distribution Crossover, Springer-Verlag New York, Inc., New [31] Q. Zhang, A. Zhou, S. Zhaoy, P.N. Suganthany, W. Liu, S. Tiwariz, Multiobjective York, NY, USA, 2003, pp. 213–237. optimization test instances for the CEC 2009 special session and competition, [23] R. Peltokangas, A. Sorsa, Real Coded Genetic Algorithm and Nonlinear Parame- Technical Report CES-487, 2009. ter Identification, Control Engineering Laboratory University of Oulu, Report A [32] E. Zitzler, M. Laumanns, L. Thiele, SPEA2: improving the strength Pareto evolu- No. 34, April 2008. tionary algorithm, in: Computer Engineering and Networks Laboratory (TIK), [24] M. Schoenauer, K. Deb, G. Rudolph, X. Yao, E. Lutton, J.J.M. Guervós, H.-P. Schwe- ETH Zurich, Zurich, Switzerland, TIK Report 103, 2001. fel (Eds.), Proceedings of the 6th International Conference on Parallel Problem [33] E. Zitzler, M. Laumanns, L. Thiele, C.M. Fonseca, V. Grunert da Fonseca, Why Solving from Nature – PPSN VI, vol. 1917, Paris, France, September 18–20, quality assessment of multiobjective optimizers is difficult, in: Genetic and Springer, 2000, Ser. Lecture Notes in Computer Science. Evolutionary Computation Conference (GECCO 2002), July, Morgan Kaufmann [25] S. Tsutsui, Multi-parent recombination in genetic algorithms with search space Publishers, New York, NY, USA, 2002, pp. 666–674. boundary extension by mirroring, in: PPSN V: Proceedings of the 5th Interna- [34] E. Zitzler, L. Thiele, An evolutionary approach for multiobjective optimiza- tional Conference on Parallel Problem Solving from Nature, Springer-Verlag, tion: the strength Pareto approach, in: Computer Engineering and Networks London, UK, 1998, pp. 428–437. Laboratory (TIK), ETH Zurich, TIK Report 43, May 1998. [26] S. Tsutsui, A. Ghosh, A study on the effect of multi-parent recombina- [35] E. Zitzler, L. Thiele, M. Laumanns, C.M. Fonseca, V.G. da Fonseca, Performance tion in realcoded genetic algorithms, in: Proc. of the IEEE ICEC, 1998, assessment of multiobjective optimizers: an analysis and review, IEEE Trans- pp. 828–833. actions on Evolutionary Computation 7 (2003) 117–132.