A Study of Multiobjective Metaheuristics When Solving Parameter Scalable Problems

A Study of Multiobjective Metaheuristics When Solving Parameter Scalable Problems Juan J. Durillo, Student Member, IEEE, Antonio J. Nebro, Carlos A. Coello Coello, Senior Member, IEEE, Jos´ e García-Nieto, Francisco Luna, and Enrique Alba Abstract—To evaluate the search capabilities of a multiobjective algorithm, the usual approach is to choose a benchmark of known problems, to perform a fixed number of function evaluations, and to apply a set of quality indicators. However, while real problems could have hundreds or even thousands of decision variables, current benchmarks are normally adopted with relatively few decision variables (normally from 10 to 30). Furthermore, performing a constant number of evaluations does not provide information about the effort required by an algorithm to get a satisfactory set of solutions; this information would also be of interest in real scenarios, where evaluating the functions defining the problem can be computationally expensive. In this paper, we study the effect of parameter scalability in a number of state-of-the-art multiobjective metaheuristics. We adopt a benchmark of parameter-wise scalable problems (the Zitzler–Deb–Thiele test suite) and analyze the behavior of eight multiobjective metaheuristics on these test problems when using a number of decision variables that range from 8 up to 2048. By using the hypervolume indicator as a stopping condition, we also analyze the computational effort required by each algorithm in order to reach the Pareto front. We conclude that the two analyzed algorithms based on particle swarm optimization and differential evolution yield the best overall results. Index Terms—Comparative study, efficiency, metaheuristics, multi-objective optimization, scalability. I. Introduction MANY REAL-WORLD optimization problems require the optimization of more than one objective function at the same time. These problems are called multiobjective This work was supported in part by the Consejería de Innovaci´on, Ciencia y Empresa, Junta de Andalucía under Contract P07-TIC-03044, the DIRI-COM Project (http://diricom.lcc.uma.es), and the Spanish Ministry of Sci-ence and Innovation under Contract TIN2008-06491C04-01, the M* Project (http://mstar.lcc.uma.es). The work of J. J. Durillo and F. Luna is supported by Grant AP-2006-03349, Grant ES-2006-13075, respectively, from the Spanish Government. The work of C. A. Coello Coello is supported by CONACyT Project 103570 J. J. Durillo, A. J. Nebro, J. García-Nieto, F. Luna, and E. Alba are with the Departamento de Lenguajes y Ciencias de la Computaci´on, University of M ´alaga, M´alaga 29071, Spain (e-mail: [email protected]; [email protected]; [email protected]; [email protected]; [email protected]). C. A. Coello Coello is with the Department of Computer Science, Center of Research and Advanced Studies, Mexico DF 07360, Mexico, and also with the Unité Mixte Internationale-Laboratoire Franco-Mexicaine d’Informatique et Automatique 3175 CNRS, Center of Research and Advanced Studies, Mexico DF 07360, Mexico (e-mail: [email protected]). optimization problems (MOPs). In contrast to single-objective optimization problems, the solution to MOPs is not a single solution, but a set of non-dominated solutions called the Pareto optimal set. A solution that belongs to this set is said to be a Pareto optimum and, when the solutions of this set are plotted in objective space, they are collectively known as the Pareto front. Obtaining the Pareto front is the main goal in multiobjective optimization. The fact that real-world MOPs tend to be nonlinear and with objective functions that are very expensive to evaluate has led to the use of metaheuristics [1], [3], [13] to deal with them. Metaheuristics are a family of techniques comprising evolutionary algorithms (EAs), particle swarm optimization (PSO), ant colony optimization, tabu search, differential evolution, scatter search, and many others. The most popular algorithms for multiobjective optimization based on metaheuristics in current use (NSGA-II [8] and SPEA2 [35]) adopt EAs as their search engine [5], [7]. The performance of these algorithms has been typically assessed using benchmark problems, such as the Zitzler–Deb– Thiele (ZDT) test problems [34], the Deb–Thiele–Laumanns– Zitzler (DTLZ) test problems [9], and the Walking-Fish-Group (WFG) test problems [15], [16]. These three problem families are scalable in the number of decision variables, and the last two are also scalable in the number of objectives. The methodology commonly adopted in the specialized literature is to compare several algorithms using a fixed (pre-defined) number of objective function evaluations and then to evaluate the values of different quality indicators (e.g., generational distance [32] or hypervolume [36], among others). The main motivation of this paper is that many real-world problems have hundreds or even thousands of decision variables, which contrasts with the current practice of validating multiobjective metaheuristics using the aforementioned benchmarks, but with a low number of decision variables (normally, no more than 30). Thus, the studies currently available do not consider the capability of current multiobjective evolutionary algorithms to properly scale when dealing with a very large number of decision variables. Scalability is a key issue in optimization algorithms, but it has been rarely addressed in the multiobjective domain. An example is the study presented in [14], in which the so-called intelligent multiobjective evolutionary algorithm (IMOEA) is compared to several elitist and non-elitist multiobjective evolutionary algorithms (MOEAs) in several of the ZDT test problems, adopting 63 decision variables. This comparative study was undertaken to validate the hypothesis of the authors of IMOEA, who argued the capability of such approach to deal with large parameter MOPs. Another example can be found in [33], where a study using only the ZDT1 problem with up to 100 variables is included. This contrasts with the interest in studying function scalability, which is currently an active research topic, leading to the area known as many-objective optimization [2], [10], [17], [26], [27], [29]. Another interesting issue that has been scarcely covered in the specialized literature is the analysis of the behavior of a multiobjective metaheuristic reaching the Pareto front of a problem. Typically, a fixed number of evaluations (and, in consequence, of iterations) is defined by the user, and the performance of the different algorithms studied is compared. However, this sort of comparison only measures the front aspect, and it does not provide any indication regarding the computational effort that a given algorithm requires to reach the true Pareto front of a problem, i.e., the efficiency of the algorithm. We believe that this is an important issue, because if we take into account that evaluating the objective functions of a MOP can be very time-consuming, it becomes of interest to know how expensive for a certain algorithm it is to generate the Pareto front. We carried a first analysis of these ideas in [11], where six state-of-the-art multiobjective metaheuristics were compared when solving the ZDT benchmark, considering their formulation ranging from 8 up to 2048 variables. The algorithms included in the study reported in [11] were three genetic algorithms (GAs) (NSGA-II [8], SPEA2 [35], and PESA-II [6]), one evolution strategy (PAES [18]), one PSO (OMOPSO [28]), and one cellular GA (MOCell [24]). In that paper, the number of evaluations required to provide a satisfactory solution was also analyzed. Given that the Pareto fronts of the ZDT problems are known, an algorithm was considered successful when the hypervolume of its current population (or archive, depending on the algorithm) was higher than the 95% of the hypervolume of the Pareto front. The current paper conducts further research along this line. Compared to our previous paper in [11], the main contributions of this paper are the following. 1) Two additional modern multiobjective metaheuristics have been included: generalized differential evolution (GDE3) [20] (a differential evolution algorithm) and AbYSS [23] (a scatter search algorithm), leading to a total of eight multiobjective metaheuristics, representative of the state-of-the-art. 2) We analyze the search capabilities of the algorithms when solving the scalable ZDT problems using a stronger stopping condition, by which the algorithms stop either when they find a solution set having a hypervolume higher than the 98% of the hypervolume of the true Pareto front or when they have performed 10 000 000 function evaluations (500 000 in [11]). 3) A more complete statistical analysis is performed, including pair-wise comparisons among the techniques in order to determine the significance of the results. 4) We study the behavior of the most promising techniques in order to propose mechanisms that enhance their search capabilities. The remainder of this paper is organized as follows. Section II includes basic background on multiobjective optimization. Sections III and IV describe, respectively, the problems and the metaheuristics we have used. Section V is devoted to the presentation and analysis of the experiments carried out. In Section VI, we include a discussion about the obtained results. Finally, Section VII summarizes the paper and discusses possible lines of future work. II. Multiobjective Optimization In this section, we include some background on multiobjective optimization. More specifically, we define the concepts of MOP, Pareto optimality, Pareto dominance, Pareto optimal set, and Pareto front. In these definitions we are assuming, without loss of generality, the minimization of all the objectives. A general multiobjective optimization problem (MOP) can be formally defined as follows. Definition 1 (MOP): Find a vector x∗=x∗ 1,x ∗ 2,...,x ∗ n which satisfies the minequality constraints gi(x)≥0,i = 1,2,...,m, the pequality constraints hi(x)=0,i = 1,2,...,p, and minimizes the vector function  f(x)= f1(x),f 2(x),...,f k(x)T, where x=[x1,x 2,...,x n]Tis the vector of decision variables. The set of all values satisfying the constraints defines the feasible region and any point x∈is a feasible solution. As mentioned before, we seek for the Pareto optima. Its formal definition is provided next. Definition 2 (Pareto Optimality): A point x∗∈is Pareto optimal if for every x∈and I={1,2,...,k }either ∀i∈I(fi(x)=fi(x∗))or there is at least one i∈Isuch that fi(x)>f i(x∗). This definition states that x∗is Pareto optimal if no feasible vector xexists which would improve some criterion without causing a simultaneous worsening in at least one other criterion. Other important definitions associated with Pareto optimality are the following. Definition 3 (Pareto Dominance): A vector u=(u1,...,u k)is said to dominate v=(v1,...,v k)(denoted by uv) if and only if uis partially smaller than v, i.e., ∀i∈{1,...,k },u i≤vi∧∃i∈{1,...,k }:ui<v i. Definition 4 (Pareto Optimal Set): For a given MOP  f(x), the Pareto optimal set is defined as P∗={x∈|¬∃x∈ ,  f(x) f(x)}. Definition 5 (Pareto Front): For a given MOP  f(x) and its Pareto optimal set P∗, the Pareto front is defined as PF∗= { f(x),x∈P∗}. Obtaining the Pareto front of a MOP is the main goal of multiobjective optimization. However, given that a Pareto front can contain a large number of points, a good solution must contain a limited number of them, which should be as close as possible to the true Pareto front, as well as being uniformly spread; otherwise, they would not be very useful to the decision maker. Closeness to the Pareto front ensures that we are dealing with optimal solutions, and a uniform spread TABLE I ZDT Test Functions Problem Objective Functions Variable Bounds Comments ZDT1 f1(x)=x1 f2(x)=g(x)[1 −x1/g(x)] g(x) = 1+9n i=2 xi/(n−1) 0≤xi≤1Convex ZDT2 f1(x)=x1 f2(x)=g(x)1−(x1/g(x))2 g(x) = 1+9n i=2 xi/(n−1) 0≤xi≤1Non-convex ZDT3 f1(x)=x1 f2(x)=g(x)1−x1 g(x)−x1 g(x)sin (10πx1) g(x) = 1+9n i=2 xi/(n−1) 0≤xi≤1Convex Disconnected ZDT4 f1(x)= x1 f2(x)= g(x)[1 −(x1/g(x))2] g(x) = 1 + 10(n−1)+ n i=2[x2 i−10 cos (4πxi)] 0≤x1≤1 −5≤xi≤5 i=2, ..., n Non-convex Multifrontal ZDT6 f1(x)=1−e−4x1sin6(6πx1) f2(x)=g(x)[1 −(f1(x)/g(x))2] g(x) = 1 + 9[(n i=2 xi)/(n−1)]0.25 0≤xi≤1 Non-convex Many-to-one Non-uniformly spaced of the solutions means that we have made a good exploration of the search space and no regions are left unexplored. III. Scalable Parameter-Wise Multiobjective Optimization Problems To carry out our paper, it would be helpful to use problems which are scalable in terms of the number of decision variables while keeping an invariable Pareto front. The ZDT test function family [34] fulfills this requirement. It offers, furthermore, a group of problems with different properties: convex, non-convex, disconnected, multifrontal, many-to-one problems (see Table I). These problems have been widely used in many studies in the field since they were first formulated. Table I shows the formulation of the ZDT test problem family. We omitted problem ZDT5 because it is binary encoded. The Pareto front of each problem is plotted in Fig. 1. Since we are interested in studying the behavior of the algorithms when solving scalable parameter-wise problems, we have evaluated each ZDT problem with 8, 16, 32, 64, 128, 256, 512, 1024, and 2048 variables. This way, we can study not only what techniques behave more efficiently when solving problems having many variables, but also if their search capabilities remain constant or not when the number of decision variables increases. IV. Multiobjective Optimization Algorithms In this section, we briefly describe the eight metaheuristics that we have considered in this paper. We have used the implementation of these algorithms provided by jMetal [12], a Java-based framework aimed at multi-objective optimization with metaheuristics.1 1jMetal is freely available for download at the following URL: http://jmetal.sourceforge.net. The NSGA-II algorithm was proposed by Deb et al. [8]. It is a genetic algorithm based on obtaining a new population from the original one by applying the typical genetic operators (selection, crossover, and mutation); then, the individuals in the two populations are sorted according to their rank, and the best solutions are chosen to create a new population. In case of having to select some individuals with the same rank, a density estimation based on measuring the crowding distance to the surrounding individuals belonging to the same rank is used to get the most promising solutions. SPEA2 was proposed by Zitler et al. in [35]. In this algorithm, each individual has a fitness value that is the sum of its strength raw fitness plus a density estimation. The algorithm applies the selection, crossover, and mutation operators to fill an archive of individuals; then, the non-dominated individuals of both the original population and the archive are copied into a new population. If the number of non-dominated individuals is greater than the population size, a truncation operator based on calculating the distances to the kth nearest neighbor is used. This way, the individuals having the minimum distance to any other individual are chosen. PESA-II [6] uses an internal population from which parents are selected to create new solutions, and an external population in which the non-dominated solutions found are stored. This external population uses the same hyper-grid division of phenotype (i.e., objective function) space adopted by PAES [19] to maintain diversity in which region-based selection is adopted. In region-based selection, the unit of selection is a hyperbox rather than an individual. The procedure lies in selecting (using any of the traditional selection techniques) a hyperbox and then randomly choosing an individual within that hyperbox. PAES is a metaheuristic proposed by Knowles and Corne [19]. The algorithm is based on a simple (1 + 1) evolution strategy. To find diverse solutions in the Pareto optimal Fig. 1. Pareto fronts of the ZDT test functions. set, PAES uses an external archive of nondominated solutions, which is also used to decide about the new candidate solutions. An adaptive grid is used as a density estimator in the archive. We have used a real coded version of PAES, applying a polynomial mutation operator. OMOPSO is a particle swarm optimization algorithm for solving MOPs [28]. Its main features include the use of the crowding distance from NSGA-II to filter out leader solutions and the use of mutation operators to accelerate the convergence of the swarm. The original OMOPSO algorithm makes use of the concept of -dominance to limit the number of solutions produced by the algorithm. We consider here a variant discarding the use of -dominance, and considering the leader population as the result yielded by the algorithm. GDE3 [20] is an improved version of the GDE algorithm [21]. It starts with a population of random solutions, which becomes the current population. At each generation, an offspring population is created using the differential evolution operators; then, the current population for the next generation is updated using the solutions of both, the offspring and the current population. Before proceeding to the next generation, the size of the population is reduced using non-dominated sorting and a pruning technique aimed at diversity preservation, in a similar way as NSGA-II, although the pruning used in GDE3 modifies the crowding distance of NSGA-II in order to solve some of its drawbacks when dealing with problems having more than two objectives. MOCell [24] is a cellular genetic algorithm (cGA). Like many multiobjective metaheuristics, it includes an external archive to store the non-dominated solutions found so far. This archive is bounded and uses the crowding distance of NSGA-II to keep diversity in the Pareto front. We have used here an asynchronous version of MOCell, called aMOCell4 in [25], in which the cells are explored sequentially (asynchronously). The selection is based on taking an individual from the neighborhood of the current solution (called cell in cGAs) and another one randomly chosen from the archive. After applying the genetic crossover and mutation operators, the new offspring is compared with the current one, replacing it if better; if both solutions are non-dominated, the worst individual in the neighborhood is replaced by the current one. In these two cases, the new individual is inserted into the archive. AbYSS is an adaptation of the scatter search metaheuristic to the multiobjective domain [23]. It uses an external archive similar to the one employed by MOCell. The algorithm incorporates operators from the evolutionary algorithms domain, including polynomial mutation and simulated binary crossover in the improvement and solution combination methods, respectively. V. Experimentation In this section, we describe the parameter settings used in the experiments, as well as the methodology we have followed in the tests, and the results we have obtained. A. Parameterization We have chosen a set of parameter values such that we allow a fair comparison among all the algorithms compared. All the GAs (NSGA-II, SPEA2, PESA-II, and MOCell) as well as GDE3 use an internal population size equal to 100; the size of the archive is also 100 in PAES, OMOPSO, GDE3, MOCell, and AbYSS. OMOPSO has been configured with 100 particles. For AbYSS, the population and the reference set have a size of 20 solutions. TABLE II Parameterization (L= Individual Length) Parameterization Used in NSGA-II [8] Population size 100 Individuals Selection of parents Binary tournament + binary tournament Recombination Simulated binary, pc=0.9 Mutation Polynomial, pm=1.0/L Parameterization Used in SPEA2 [35] Population size 100 individuals Selection of parents Binary tournament + binary tournament Recombination Simulated binary, pc=0.9 Mutation Polynomial, pm=1.0/L Parameterization Used in PESA-II [6] Population size 100 individuals Selection of parents Region based selection + region based selection Recombination Simulated binary, pc=0.9 Mutation Polynomial, pm=1.0/L Archive size 100 individuals Parameterization Used in PAES [18] Mutation Polynomial, pm=1.0/L Archive size 100 individuals Parameterization Used in OMOPSO [28] Particles 100 particles Mutation Uniform + non-uniform Leaders size 100 individuals Parameterization Used in GDE3 [20] Population size 100 individuals Recombination Differential evolution, CR =0.1, F=0.5 Parameterization Used in MOCell [25] Population size 100 individuals (10 ×10) Neighborhood 1-hop neighbors (8 surrounding solutions) Selection of parents Binary tournament + binary tournament Recombination Simulated binary, pc=0.9 Mutation Polynomial, pm=1.0/L Archive size 100 individuals Parameterization Used in AbYSS [23] Population size 20 individuals Reference set size 10 + 10 Recombination Simulated binary, pc=1.0 Mutation (local search) Polynomial, pm=1.0/L Archive size 100 individuals In the GAs, we have used simulated binary crossover (SBX) and polynomial mutation [7]. The distribution indices for both operators are ηc= 20 and ηm= 20, respectively. The crossover probability is pc=0.9 and the mutation probability is pm=1/L, where Lis the number of decision variables. In PAES we have also adopted a polynomial mutation operator, with the same distribution index as indicated before. AbYSS uses polynomial mutation in the improvement method and SBX in the solution combination method. GDE3 uses 0.1 and 0.5 for the two parameters CR and F, respectively [20]. OMOPSO applies a combination of uniform and non-uniform mutation to the particle swarm [28]. A detailed description of the parameter values adopted for our experiments is provided in Table II. B. Methodology We are interested in two main goals: analyzing the behavior of the algorithms when solving the scalable ZDT benchmark and their speed (efficiency) in reaching the Pareto front. Given that the Pareto fronts of the ZDT problems are known beforehand, a strategy could be to run the algorithms until they are able to produce them. However, it is possible that some of them never produce the true Pareto front, or simply take too Fig. 2. Pareto fronts with different HV values obtained for problem ZDT1. long to do it. Thus, we adopt instead a stopping condition for all the algorithms compared, based on the high quality of the Pareto front produced. For that purpose, the hypervolume [36] quality indicator is adopted. The hypervolume computes the volume (in objective function space) covered by members of a non-dominated set of solutions Qfor problems in which all objectives are to be minimized. Mathematically, for each solution i∈Q,a hypercube viis constructed with a reference point Wand the solution ias the diagonal corners of the hypercube. The reference point can simply be found by constructing a vector of worst objective function values. Thereafter, the union of all hypercubes is found and its hypervolume (HV ) is calculated HV =volume |Q|  i=1 vi.(1) Higher values of the HV performance measure imply more desirable solutions. A property of this quality indicator is that it measures both convergence to the Pareto front and diversity of the obtained fronts. Once the quality indicator we are going to use has been described, we need to establish a stopping condition to be used in the execution of the algorithms. The idea is that the metaheuristics stop when they reach a certain percentage of the HV of the Pareto front, which ensures that the obtained front represents an accurate approximation to it. To decide about that percentage, we show different approximations of the Pareto front for the problem ZDT1 with different percentages of HV in Fig. 2. We can observe that a front with a hypervolume of 98.26% represents a reasonable approximation to the true Pareto fronts in terms of convergence and diversity of solutions. This same value has been corroborated using the other test problems from the ZDT suite. So, we have taken 98% of the HV of the Pareto front as a criterion to consider that a MOP has been successfully solved. Furthermore, those algorithms requiring fewer function evaluations to achieve this termination condition can be considered to be more efficient or faster. In those situations in which an algorithm is unable to obtain a front fulfilling this condition after the maximum number of function evaluations, we consider that it has failed in solving the problem; this way, we can obtain a hit rate for the algorithms, i.e., their percentage of successful executions. We set the maximum number of evaluations to ten million. Fig. 3. Statistical analysis performed in this paper. In our experiments, we check the stopping condition every 100 evaluations (that is, each iteration in the population-based metaheuristics), where we measure the hypervolume of the non-dominated solutions found so far. Therefore, in NSGA-II, SPEA2, and GDE3 we have considered the non-dominated solutions at each generation: in PESA-II, PAES, AbYSS, and MOCell, the external population and, in MOPSO, the leaders archive. We have executed 100 independent runs for each algorithm and each problem instance. Since we are dealing with stochastic algorithms, we need to perform a statistical analysis of the obtained results to compare them with a certain level of confidence. Next, we describe the statistical test that we have carried out for assuring this [30]. First, a Kolmogorov– Smirnov test is performed in order to check whether the values of the results follow a normal (Gaussian) distribution or not. If so, the Levene test checks for the homogeneity of the variances. If samples have equal variance (positive Levene test), an ANOVA test is done; otherwise we perform a Welch test. For non-Gaussian distributions, the non-parametric Kruskal–Wallis test is used to compare the medians of the algorithms. We always consider in this paper a confidence level of 95% (i.e., significance level of 5% or p-value under 0.05) in the statistical tests, which means that the differences are unlikely to have occurred by chance with a probability of 95%. Successful tests are marked with “+” symbols in the last row in the tables containing the results (see Tables III–VII); conversely, “−” means that no statistical confidence was found (p-value >0.05). For the sake of homogeneity in the presentation of the results, all the tables include the median, ˜x, and the interquartile range, IQR, as measures of location (or central tendency) and statistical dispersion, respectively. We have performed a post-hoc testing phase using the multcompare function provided by MATLAB, which allows for a multiple comparison of samples. This way, we can make pairwise comparison between algorithms to know about the significance of their results. C. Analysis of Results Tables III–VII show the median and the interquartile range of the number of evaluations needed by the different optimizers for ZDT1, ZDT2, ZDT3, ZDT4, ZDT6, respectively. This indicates that all the 100 independent runs have been successful, which means a hit rate of 1.0. When an optimizer Fig. 4. Number of evaluations when solving ZDT1. is not able to reach an acceptable front upon performing 10 000 000 function evaluations in all the 100 independent runs, its cell in the tables includes the “−” symbol, and it is not taken into account in the statistical tests. In other words, the “−” symbol means that, in order to solve successfully the problem in all the independent runs, the optimizer may need more than 10 000 000 of function evaluations. In these cases, the hit rate would be less than 1.0. To ease the analysis of the results in these tables, the cells containing the lowest number of function evaluations have a grey colored background. There are two grey levels: the darker grey indicates the best (lowest) value, while lighter grey is used to point out the second best value. We can observe that the results in these tables are significant, as can be seen in the last row of each of them, where each cell contains a “+” symbol, except for ZDT4 with 1024 and 2048 variables (where there are no results to compare). Next, we analyze the results obtained for each of the problems. To make the results clearer, we include a figure summarizing the values, using a logarithmic scale, in addition to the corresponding table. The discussion is organized in the following order: first, we analyze the success of the algorithms when solving the different instances of the problem; second, we analyze the speed of the techniques to obtain the Pareto front when solving the problems; and, finally, we make a pairwise comparison of the algorithms, considering those problems in which there are no statistical differences between each pair of techniques (if we included the problems with statistical significance, this would lead to larger tables). 1) ZDT1: This problem presents a uniform density of solutions in the search space and a convex Pareto front. Table III and Fig. 4 show the number of evaluations needed to obtain a Pareto front with 98% of the HV of the Pareto front for this problem. Table VIII presents the hit rate indicator, where we have used the symbol √to indicate a value of 1.0 (a 100% success rate). Table IX contains the problems for which no statistical differences have been found between each pair of algorithms. We begin by analyzing which algorithms are able to solve the problems in all the independent runs carried out. The results in Table III show that all the algorithms have success in all the instances, except for PAES in the instance with 2048 TABLE III Evaluations for ZDT1 Algorithm 8 16 32 64 128 256 512 1024 2048 variables variables variables variables variables variables variables variables variables NSGA-II 4.40e+3 4.0e+2 8.10e+3 5.0e+2 1.52e+4 9.0e+2 2.88e+4 1.4e+3 5.72e+4 3.0e+3 1.21e+5 5.2e+3 2.71e+5 9.3e+3 6.31e+5 1.6e+4 1.45e+6 3.3e+4 SPEA2 5.00e+3 3.0e+2 9.30e+3 8.0e+2 1.69e+4 1.1e+3 3.14e+4 1.4e+3 6.03e+4 2.7e+3 1.25e+5 5.0e+3 2.69e+5 9.4e+3 6.01e+5 1.6e+4 1.34e+6 2.8e+4 PESA-II 4.15e+3 1.2e+3 8.40e+3 9.5e+2 1.73e+4 1.6e+3 3.74e+4 2.8e+3 8.36e+4 5.0e+3 1.96e+5 9.3e+3 4.78e+5 1.8e+4 1.18e+6 4.0e+4 2.93e+6 8.9e+4 PAES 3.40e+3 2.5e+3 7.25e+3 4.4e+3 1.34e+4 8.0e+3 2.57e+4 1.6e+4 4.70e+4 3.5e+4 9.07e+4 3.9e+4 1.73e+5 1.3e+5 3.68e+5 4.7e+5 – OMOPSO 1.40e+3 4.0e+2 3.40e+3 9.0e+2 7.40e+3 1.8e+3 1.38e+4 2.8e+3 2.80e+4 4.4e+3 6.31e+4 7.4e+3 1.58e+5 1.8e+4 4.55e+5 3.4e+4 1.41e+6 1.0e+5 GDE3 2.80e+3 2.0e+2 5.30e+3 3.0e+2 1.00e+4 4.0e+2 1.81e+4 4.5e+2 3.30e+4 9.0e+2 6.10e+4 1.2e+3 1.16e+5 1.8e+3 2.40e+5 3.0e+3 5.64e+5 7.0e+3 MOCell 1.80e+3 3.0e+2 3.80e+3 6.0e+2 9.20e+3 8.0e+2 2.18e+4 1.8e+3 4.92e+4 4.0e+3 1.13e+5 6.4e+3 2.56e+5 8.6e+3 5.62e+5 1.6e+4 1.20e+6 2.7e+4 AbYSS 3.40e+3 7.0e+2 6.80e+3 9.0e+2 1.49e+4 1.8e+3 3.30e+4 2.7e+3 7.29e+4 6.1e+3 1.58e+5 7.8e+3 3.50e+5 1.2e+4 7.06e+5 2.1e+4 1.39e+6 2.6e+4 + + + + + + + + – TABLE IV Evaluations for ZDT2 Algorithm 8 16 32 64 128 256 512 1024 2048 variables variables variables variables variables variables variables variables variables NSGA-II 7.50e+3 6.0e+2 1.37e+4 8.0e+2 2.60e+4 1.6e+3 4.98e+4 3.2e+3 1.01e+5 4.8e+3 2.08e+5 1.0e+4 4.54e+5 1.8e+4 1.01e+6 2.6e+4 2.28e+6 4.5e+4 SPEA2 7.80e+3 1.0e+3 1.46e+4 1.3e+3 2.60e+4 1.7e+3 4.86e+4 3.2e+3 9.50e+4 4.4e+3 1.90e+5 6.0e+3 3.95e+5 9.0e+3 8.51e+5 1.7e+4 1.89e+6 2.6e+4 PESA-II 1.50e+4 7.8e+3 3.51e+4 2.5e+4 9.38e+4 6.0e+4 3.98e+5 4.2e+5 −− − − − PAES 3.14e+4 4.5e+4 6.94e+4 7.1e+4 1.29e+5 1.2e+5 2.93e+5 4.2e+5 − − − − − OMOPSO 1.75e+3 4.0e+2 3.90e+3 1.6e+3 9.70e+3 3.8e+3 1.68e+4 5.0e+3 3.06e+4 5.3e+3 5.44e+4 8.4e+3 1.11e+5 1.1e+4 2.83e+5 2.7e+4 7.71e+5 7.9e+4 GDE3 3.20e+3 2.0e+2 6.10e+3 4.5e+2 1.18e+4 6.5e+2 2.26e+4 1.2e+3 4.33e+4 1.8e+3 8.18e+4 1.3e+3 1.58e+5 2.1e+3 3.27e+5 4.2e+3 7.66e+5 9.8e+3 MOCell 2.90e+3 1.0e+3 4.90e+3 1.4e+3 8.25e+3 2.2e+3 1.74e+4 1.1e+4 4.42e+4 2.5e+4 1.26e+5 4.3e+4 2.91e+5 1.1e+4 6.60e+5 1.1e+4 1.46e+6 2.0e+4 AbYSS 4.50e+3 7.5e+2 9.10e+3 1.7e+3 1.86e+4 2.5e+3 3.96e+4 3.5e+3 8.26e+4 7.2e+3 1.75e+5 8.6e+3 3.68e+5 9.6e+3 7.74e+5 1.9e+4 1.60e+6 2.4e+4 + + + + + + + + + TABLE V Evaluations for ZDT3 Algorithm 8 16 32 64 128 256 512 1024 2048 variables variables variables variables variables variables variables variables variables NSGA-II 4.20e+3 5.0e+2 7.40e+3 6.0e+2 1.36e+4 9.5e+2 2.53e+4 1.4e+3 5.06e+4 3.1e+3 1.08e+5 4.8e+3 2.38e+5 9.0e+3 5.37e+5 1.6e+4 1.20e+6 2.5e+4 SPEA2 4.90e+3 6.0e+2 9.10e+3 7.0e+2 1.62e+4 1.0e+3 3.00e+4 2.0e+3 5.89e+4 4.0e+3 1.21e+5 5.8e+3 2.56e+5 9.4e+3 5.58e+5 1.6e+4 1.22e+6 2.5e+4 PESA-II 3.75e+3 9.0e+2 7.60e+3 1.0e+3 1.59e+4 2.0e+3 3.43e+4 3.8e+3 7.76e+4 7.8e+3 1.77e+5 8.0e+3 4.14e+5 1.5e+4 9.83e+5 3.3e+4 2.28e+6 4.5e+4 PAES 6.90e+3 8.8e+3 1.21e+4 1.2e+4 2.56e+4 3.2e+4 5.68e+4 6.1e+4 9.92e+4 1.2e+5 2.03e+5 2.0e+5 3.65e+5 6.3e+5 8.17e+5 8.4e+5 – OMOPSO 2.60e+3 1.0e+3 5.40e+3 2.4e+3 1.01e+4 3.3e+3 2.20e+4 4.0e+3 5.06e+4 8.1e+3 1.26e+5 2.3e+4 3.53e+5 3.9e+4 1.03e+6 8.3e+4 3.17e+6 1.9e+5 GDE3 2.90e+3 2.5e+2 5.60e+3 3.0e+2 1.08e+4 4.0e+2 1.96e+4 8.0e+2 3.46e+4 9.5e+2 6.29e+4 1.3e+3 1.20e+5 1.8e+3 2.50e+5 2.8e+3 6.07e+5 7.2e+3 MOCell 1.90e+3 6.0e+2 4.20e+3 8.0e+2 9.90e+3 1.4e+3 2.30e+4 2.2e+3 5.24e+4 3.4e+3 1.16e+5 9.3e+3 2.57e+5 1.7e+4 5.44e+5 1.6e+4 1.15e+6 3.6e+4 AbYSS 3.35e+3 1.6e+3 6.75e+3 1.8e+3 1.40e+4 2.1e+3 2.87e+4 3.8e+3 6.02e+4 6.1e+3 1.23e+5 1.6e+4 2.57e+5 2.0e+4 5.47e+5 4.4e+4 1.12e+6 4.5e+4 + + + + + + + + + TABLE VI Evaluations for ZDT4 Algorithm 8 16 32 64 128 256 512 1024 2048 variables variables variables variables variables variables variables variables variables NSGA-II 1.62e+4 3.4e+3 4.38e+4 1.2e+4 1.37e+5 3.3e+4 4.25e+5 9.1e+4 1.20e+6 1.8e+5 3.29e+6 3.5e+5 8.85e+6 6.7e+5 − − SPEA2 2.11e+4 4.2e+3 4.54e+4 8.7e+3 1.34e+5 3.1e+4 3.69e+5 6.4e+4 1.07e+6 1.3e+5 2.98e+6 2.4e+5 8.24e+6 5.9e+5 −− PESA-II 1.84e+4 5.2e+3 5.00e+4 1.4e+4 1.51e+5 3.2e+4 4.12e+5 7.3e+4 1.06e+6 1.3e+5 2.72e+6 2.5e+5 6.69e+6 3.8e+5 − − PAES 3.08e+4 1.1e+4 8.53e+4 4.1e+4 2.17e+5 6.4e+4 5.41e+5 1.9e+5 1.28e+6 3.0e+5 3.18e+6 7.2e+5 − − − OMOPSO − − − − − − − − − GDE3 1.18e+4 8.0e+2 − − − − − − − − MOCell 8.80e+3 2.2e+3 2.04e+4 5.4e+3 5.86e+4 1.0e+4 1.65e+5 2.7e+4 4.67e+5 5.8e+4 1.23e+6 1.0e+5 3.25e+6 1.6e+5 8.14e+6 2.8e+5 − AbYSS 1.46e+4 5.4e+3 5.46e+4 2.1e+4 1.61e+5 4.3e+4 4.82e+5 1.0e+5 1.39e+6 1.5e+5 3.81e+6 4.2e+5 − − − + + +−+ + + − − TABLE VII Evaluations for ZDT6 Algorithm 8 16 32 64 128 256 512 1024 2048 variables variables variables variables variables variables variables variables variables NSGA-II 2.31e+4 1.3e+3 4.60e+4 1.4e+3 8.84e+4 2.8e+3 1.68e+5 3.5e+3 3.24e+5 6.4e+3 6.47e+5 1.1e+4 1.33e+6 1.6e+4 2.77e+6 2.7e+4 5.80e+6 4.6e+4 SPEA2 2.64e+4 1.2e+3 5.26e+4 2.0e+3 9.94e+4 2.2e+3 1.87e+5 4.2e+3 3.53e+5 7.0e+3 6.84e+5 1.1e+4 1.37e+6 1.4e+4 2.81e+6 2.3e+4 5.82e+6 3.8e+4 PESA-II 2.16e+4 1.8e+3 4.78e+4 2.4e+3 9.85e+4 4.7e+3 2.01e+5 7.2e+3 4.10e+5 1.1e+4 8.57e+5 2.1e+4 1.83e+6 2.9e+4 3.89e+6 4.4e+4 8.32e+6 8.9e+4 PAES 6.80e+3 6.8e+3 1.62e+4 1.6e+4 3.21e+4 2.8e+4 − − − − − − OMOPSO 2.90e+3 1.6e+3 4.20e+3 1.9e+3 7.70e+3 2.7e+3 1.38e+4 3.5e+3 2.88e+4 5.2e+3 6.10e+4 1.0e+4 1.26e+5 1.6e+4 2.76e+5 3.0e+4 6.12e+5 6.6e+4 GDE3 3.70e+3 5.0e+2 6.60e+3 5.0e+2 1.32e+4 1.1e+3 3.14e+4 2.6e+3 1.53e+5 6.7e+3 3.21e+5 3.8e+3 6.36e+5 4.0e+3 1.35e+6 7.5e+3 3.27e+6 1.8e+4 MOCell 1.07e+4 1.0e+3 2.58e+4 1.4e+3 5.65e+4 2.6e+3 1.19e+5 3.0e+3 2.47e+5 4.7e+3 5.18e+5 7.8e+3 1.10e+6 9.4e+3 2.35e+6 1.9e+4 5.00e+6 3.4e+4 AbYSS 1.23e+4 1.1e+3 2.57e+4 1.8e+3 5.26e+4 3.0e+3 1.08e+5 3.8e+3 2.26e+5 7.2e+3 4.67e+5 9.2e+3 9.60e+5 1.4e+4 1.96e+6 2.1e+4 4.00e+6 4.1e+4 + + + + + + + + + TABLE VIII Hit Rate for ZDT1 Algorithm/Variables 8 16 32 64 128 256 512 1024 2048 NSGA-II √ √ √ √ √ √ √ √ √ SPEA2 √√√√ √ √ √ √ √ PESA-II √√√√√ √ √ √ √ PAES √√ √ √ √ √ √ √0.95 OMOPSO √ √ √ √ √√ √ √ √ GDE3 √√ √ √ √ √ √ √ √ MOCell √ √ √ √ √ √ √ √ √ AbYSS √ √ √ √ √ √ √ √ √ Fig. 5. Number of evaluations when solving ZDT2. variables. This is corroborated considering the hit rate (see Table VIII), where all the algorithms have a hit rate of 1.0 except for PAES in the 2048 instance, which has 0.95. This means that in five out of the 100 executions carried out for this instance, PAES reached the maximum number of evaluations before obtaining a Pareto front with the desired HV value. We pay attention now to the speed, i.e., the number of function evaluations needed by the metaheuristics to find a Pareto front according to our success condition. In Fig. 4, we plot the results using a logarithmic scale. We have connected by a line the symbols of the two algorithms yielding the best values. Thus, we can observe that OMOPSO (dotted line) is the fastest algorithm up to 128 variables, while GDE3 scales better from 256 to 2048 variables. The lines clearly depict that there is a tendency change in these two algorithms at 256 variables, indicating that GDE3 tends to be faster than the other techniques as the number of decision variables Fig. 6. Number of evaluations when solving ZDT3. increases, while OMOPSO exhibits the opposite behavior. This suggests that GDE3 could be the most appropriate algorithm to solve ZDT1 with more than 2048 variables. Accordingly, we determine that GDE3 is the algorithm that scales the best in problem ZDT1. Considering the rest of techniques, MOCell is the second fastest algorithm up to 32 variables and in the case of 2048 variables. NSGA-II, SPEA2, and AbYSS tend to need a similar number of evaluations when the instances are larger. Finally, we make an analysis considering the outcome of the statistical tests included in Table IX. If we focus on OMOPSO and GDE3, the tests are non-successful only in the instances of 128 and 256 variables, which does not affect the previous analysis. It is interesting to note that PAES, the simplest of the evaluated techniques, presents no differences in many instances compared to NSGA-II and MOCell. We can also observe that NSGA-II TABLE IX Non-Successful Statistical Tests of the Number of Evaluations for ZDT1 NSGA-II 128, 256, 512, 1024 8, 16 16, 32, 64, 128 2048 – 256 32 SPEA2 16, 32 2048 − − −64 PESA-II 8 − − − 128 PAES 512, 1024 −128, 256, 512, 1024 8, 16, 32 OMOPSO 128, 256 8, 16, 32 2048 GDE3 32 − MOCell − SPEA2 PESA-II PAES OMOPSO GDE3 MOCell AbYSS TABLE X Hit Rate for ZDT2 Algorithm/Variables 8 16 32 64 128 256 512 1024 2048 NSGA-II √ √ √ √ √ √ √ √ √ SPEA2 √√ √ √ √ √ √ √ √ PESA-II √√ √ √ 0.43 0.0 0.0 0.0 0.0 PAES √ √ √ √ 0.99 0.98 0.99 0.96 0.79 OMOPSO √√ √ √ √ √ √ √ √ GDE3 √√ √ √ √ √ √ √ √ MOCell √ √ √ √ √ √ √ √ √ AbYSS √ √ √ √ √ √ √ √ √ and SPEA2 do not present statistical confidence in four cases. 2) ZDT2: This problem also presents a uniform density of solutions in the search space, but it has a non-convex Pareto front. The numbers of evaluations needed to solve the ZDT2 problem are included in Table IV and Fig. 5. Table X shows the hit rate of each algorithm. The problems in which each pair of algorithms are statistically independent appear in Table XI. We start by commenting that some optimizers have difficulties when solving this problem, as can be seen in Table IV and Fig. 5. In particular, neither PAES nor PESA-II are able to reach a hit rate of 1.0 in instances with more than 64 variables. These results indicate that problems having a nonconvex Pareto front can cause difficulties for some algorithms. Consequently, the hit rate indicator shows that six algorithms have a 100% success rate in all the instances. Among the algorithms with a hit rate smaller than 1.0, PESA-II is the algorithm having the worst value, 0.49 in the problem with 128 variables, and 0.0 in the next ones. PAES does not achieve a 100% of success after 128 variables but, in contrast to PESA-II, the values are near 100% in the instances ranging from 128 to 1024 variables. Let us examine now the speed of the algorithms. A look at Fig. 5 reveals that OMOPSO and GDE3 are the fastest algorithms. The lines connecting the values of these solvers indicate that OMOPSO requires a lower number of function evaluations than GDE3 to get the desired Pareto fronts in all but in the largest instance. In fact, when observing the lines we can see that the one of GDE3 suggests again that it could scale better than the other techniques in order to solve instances of more than 2048 variables. MOCell, AbYSS, SPEA2, and NSGA-II, in this order, are the following metaheuristics in terms of speed, although they tend to get close to each Fig. 7. Number of evaluations when solving ZDT4. other when the number of decision variables of the problem increases. The pairwise tests in Table XI reveal some interesting facts. On the one hand, the results of GDE3 and OMOPSO in the two largest instances are non-significant; on the other hand, the tests show the differences between MOCell and OMOPSO up to 64 variables and between MOCell and GDE3 up to 512 variables are also non-significant. 3) ZDT3: This problem presents a uniform density of solutions in the search space and its Pareto front is composed of several discontinuous regions; therefore, the main complexity is to find all these discontinuous regions. The evaluations required for solving the different instances are shown in Table V and Fig. 6. Table XII presents the hit rate of the algorithms. The problems in which there are not statistical differences between each pair of algorithms are shown in Table XIII. Proceeding as before, we start by analyzing which algorithms are able to solve successfully all the instances of the problem. Therefore, in Table V we see that all the algorithms Fig. 14. Results of GDE3 with polynomial mutation on the ZDT problem family. TABLE XXII GDE3 with Polynomial Mutation: Evaluations Problem 8 16 32 64 128 256 512 1024 2048 variables variables variables variables variables variables variables variables variables ZDT1 3.30e+3 3.0e+2 6.50e+3 4.0e+2 1.18e+4 5.0e+2 2.12e+4 8.0e+2 3.80e+4 9.0e+2 6.98e+4 1.4e+3 1.33e+5 2.4e+3 2.74e+5 3.5e+3 6.45e+5 8.4e+3 ZDT2 4.10e+3 3.0e+2 8.00e+3 5.0e+2 1.54e+4 7.5e+2 2.88e+4 9.0e+2 5.41e+4 1.6e+3 1.02e+5 1.8e+3 1.98e+5 3.4e+3 4.14e+5 4.4e+3 9.91e+5 1.4e+4 ZDT3 2.60e+3 1.0e+3 5.40e+3 2.4e+3 1.01e+4 3.3e+3 2.20e+4 4.0e+3 5.06e+4 8.1e+3 1.26e+5 2.3e+4 3.53e+5 3.9e+4 1.03e+6 8.3e+4 3.17e+6 1.9e+5 ZDT4 2.11e+4 1.4e+3 5.98e+4 4.8e+3 1.88e+5 1.4e+4 7.78e+5 6.9e+4 3.41e+6 4.6e+5 ––– – ZDT6 4.50e+3 4.0e+2 9.10e+3 7.5e+2 1.90e+4 2.0e+3 5.02e+4 5.2e+3 2.24e+5 8.0e+3 4.67e+5 5.0e+3 9.39e+5 7.4e+3 2.02e+6 1.4e+4 4.97e+6 3.3e+4 changes in the behavior of the technique. This, however, does not mean that there is no room for improvement; the use of different mutation and crossover operators may change the search capabilities of a GA. E. Comparison with Previous Paper In this section, we compare the results of this paper with those obtained in some of our previous paper [11]. In that paper, the main conclusions were that, considering scalability, PAES was the most competitive algorithm followed by OMOPSO, while the latter appeared as the fastest algorithm. In this paper, however, PAES appears in the last positions in the scalabily ranking. This is explained by the fact that we have not considered here the impact of not having a 100% hit rate in the tables containing the computed evaluations, which undoubtedly has penalized PAES in practically all the problems considered. If we consider OMOPSO, its results in [11] remain basically the same. The inclusion in this paper of GDE3 and AbYSS has affected OMOPSO only in the scalabity ranking, and OMOPSO is still among the fastest algorithms assessed (excluding the ZDT4 test problem). Finally, we would like to mention the effort that has been required to carry out this paper. If we consider that we have studied eight algorithms to solve nine instances of five problems, the number of experiments is 360. As we have executed one hundred independent runs per experiment, this means a total of 36 000 runs. Furthermore, the experiments with SMPSO and GDE3 with mutation have required 9000 additional runs. Taking into account also the pilot tests with the other algorithms to study the influence of the distribution indices in the mutation and crossover operators, the total number of runs have been in the order of 50 000. The fact that we have used a stronger termination condition (98% instead of 95% of the hypervolume of the Pareto front) and that the maximum number of evaluations have been raised from 500 000 in [11] to 10 000 000 in this paper, has had as a consequence the fact that the computational time of the algorithms unable to solve the problems with 1024 and 2048 variables was higher than one hour. To execute this large number of experiments, we have used the computers of the laboratories of the Department of Computer Science of the University of M´ alaga, in Spain. Most of them are equipped with modern dual core processors so, taking into account that there are more than 180 computers, that means that up to 360 cores have been available. To run all the programs, we have used Condor [31], a middleware that acts as a distributed scheduler, which has proven to be an ideal tool to cope with the large amount of tasks we have dealt with. VII. Conclusion and Future Work We have evaluated eight state-of-the-art metaheuristics over a set of parameter scalable problems in order to study the behavior of the algorithms concerning their capabilities to solve problems having a large number of decision variables. The benchmark has been composed of five problems from the ZDT family, using instance sizes ranging from 8 to 2048 variables. We have also studied the speed of the techniques when solving the problems. The stopping condition has been to reach a front with a hypervolume higher than the 98% of the hypervolume of the true Pareto front, or to compute 10 000 000 function evaluations. Our study has revealed that differential evolution and particle swarm optimization are the most promising approaches to deal with the scalable problems used in this paper. GDE3 and OMOPSO do not only scale well, but they are among the fastest algorithms. Furthermore, we have shown that their search capabilities can be improved to solve ZDT4, the problem which has appeared as the most difficult one to solve. Two modern optimizers, MOCell and AbYSS, have shown a high degree of regularity in the tests. With the exception of MOCell in ZDT4 (where it is the overall best technique), they are not in the first position in the scalabily and the speed rankings, but also they are always around the third and fifth positions. Both metaheuristics are in the group of algorithms having solved a higher number of instances. In this group we find NSGA-II and SPEA2, which are, in general, very close in the rankings, but they usually appear after MOCell. PESA-II has difficulties in ZDT2 and it normally appears among the algorithms requiring higher numbers of function evaluations to reach a front with the target HV value. Finally, PAES, the simplest of the optimizers in the paper, is the algorithm scaling the worst, due to the low hit rates it obtains in many instances. We have presented a paper about the behavior of eight multiobjective metaheuristics concerning their scalability and speed when solving a scalable benchmark. The next step is an extension including more scalable problems (e.g., DTLZ and WFG) to assess whether or not the features of these problems confirm the results obtained in this paper. Our analysis of OMOPSO and GDE3 has also shown that an open research line is to study variations and different parameter settings of the existing multiobjective metaheuristics in order to improve their scalability, and that is another path for future research that we aim to explore. References [1] C. Blum and A. Roli, “Metaheuristics in combinatorial optimization: Overview and conceptual comparison,” ACM Comput. Surveys, vol. 35, no. 3, pp. 268–308, 2003. [2] D. Brockhoff and E. Zitzler, “Are all objectives necessary? On dimensionality reduction in evolutionary multiobjective optimization,” in Proc. Ninth Int. Conf. Parallel Problem Solving Nature (PPSN), Lecture Notes in Computer Science 4193. Reykjavik, Iceland, Sep. 2006, pp. 533–542. [3] E. K. Burke and G. Kendall, Search Methodologies: Introductory Tutorials in Optimization and Decision Support Techniques. New York: Springer, 2005. [4] M. Clerc and J. Kennedy, “The particle swarm: Explosion, stability, and convergence in a multidimensional complex space,” IEEE Trans. Evol. Comput., vol. 6, no. 1, pp. 58–73, 2002. [5] C. A. Coello Coello, G. B. Lamont, and D. A. Van Veldhuizen, Evolutionary Algorithms for Solving Multiobjective Problems, 2nd ed. New York: Springer, Sep. 2007. [6] D. W. Corne, N. R. Jerram, J. D. Knowles, and M. J. Oates, “PESA-II: Region-based selection in evolutionary multiobjective optimization,” in Proc. Genetic Evol. Comput. Conf. (GECCO), San Mateo, CA: Morgan Kaufmann, 2001, pp. 283–290. [7] K. Deb, Multiobjective Optimization Using Evolutionary Algorithms. Chichester, U.K.: Wiley, 2001. [8] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan, “A fast and elitist multiobjective genetic algorithm: NSGA-II,” IEEE Trans. Evol. Comput., vol. 6, no. 2, pp. 182–197, Apr. 2002. [9] K. Deb, L. Thiele, M. Laumanns, and E. Zitzler, “Scalable test problems for evolutionary multiobjective optimization,” in Evolutionary Multiobjective Optimization: Theoretical Advances and Applications,A. Abraham, L. Jain, and R. Goldberg, Eds. New York: Springer, 2005, pp. 105–145. [10] F. di Pierro, S. Khu, and D. A. Savi´ c, “An investigation on preference order ranking scheme for multiobjective evolutionary optimization,” IEEE Trans. Evol. Comput., vol. 11, no. 1, pp. 17–45, Feb. 2007. [11] J. J. Durillo, A. J. Nebro, C. A. Coello Coello, F. Luna, and E. Alba, “A comparative study of the effect of parameter scalability in multiobjective metaheuristics,” in Proc. Congr. Evol. Comput. (CEC), Jun. 2008, pp. 1893–1900. [12] J. J. Durillo, A. J. Nebro, F. Luna, B. Dorronsoro, and E. Alba, “Jmetal: A Java framework for developing multiobjective optimization metaheuristics,” Dept. Lenguajes Ciencias Computaci´ on, Univ. M´ alaga, ETSI Inform´ atica, M´ alaga, Spain, Tech. Rep. ITI-2006-10, 2006. [13] F. W. Glover and G. A. Kochenberger, Handbook of Metaheuristics. Norvell, MA: Kluwer, 2003. [14] S. Ho, L. Shu, and J. Chen, “Intelligent evolutionary algorithms for large parameter optimization problems,” IEEE Trans. Evol. Comput., vol. 8, no. 6, pp. 522–541, Dec. 2004. [15] S. Huband, L. Barone, R. L. While, and P. Hingston, “A scalable multiobjective test problem toolkit,” in Proc. Third Int. Conf. Evol. Multicriterion Optimization (EMO), Lecture Notes in Computer Science 3410. Berlin, Germany: Springer, 2005, pp. 280–295. [16] S. Huband, P. Hingston, L. Barone, and L. While, “A review of multiobjective test problems and a scalable test problem toolkit,” IEEE Trans. Evol. Comput., vol. 10, no. 5, pp. 477–506, Oct. 2006. [17] V. Khare, X. Yao, and K. Deb, “Performance scaling of multiobjective evolutionary algorithms,” in Proc. Second Int. Conf. Evol. MultiCriterion Optimization (EMO), Lecture Notes in Computer Science 2632. Faro, Portugal, Apr. 2003, pp. 376–390. [18] J. Knowles and D. Corne, “The pareto archived evolution strategy: A new baseline algorithm for multiobjective optimization,” in Proc. Congr. Evol. Comput., 1999, pp. 98–105. [19] J. D. Knowles and D. W. Corne, “Approximating the nondominated front using the pareto archived evolution strategy,” Evol. Comput., vol. 8, no. 2, pp. 149–172, 2000. [20] S. Kukkonen and J. Lampinen, “GDE3: The third evolution step of generalized differential evolution,” in Proc. IEEE Congr. Evol. Comput. (CEC), 2005, pp. 443–450. [21] J. Lampinen, “DE’s selection rule for multiobjective optimization,” Dept. Inf. Technol., Lappeenranta Univ. Technol., Lappeenranta, Finland, Tech. Rep., 2001. [22] H. Li and Q. Zhang, “Multiobjective optimization problems with complicated pareto set, MOEA/D and NSGA-II,” IEEE Trans. Evol. Comput., vol. 2, no. 12, pp. 284–302, Apr. 2009. [23] A. J. Nebro, F. Luna, E. Alba, B. Dorronsoro, J. J. Durillo, and A. Beham, “AbYSS: Adapting scatter search to multiobjective optimization,” IEEE Trans. Evol. Comput., vol. 12, no. 4, pp. 439–457, Aug. 2008. [24] A. J. Nebro, J. J. Durillo, F. Luna, B. Dorronsoro, and E. Alba, “A cellular genetic algorithm for multiobjective optimization,” in Proc. Workshop Nature Inspired Cooperative Strategies Optimization (NICSO), Granada, Spain, 2006, pp. 25–36. [25] A. J. Nebro, J. J. Durillo, F. Luna, B. Dorronsoro, and E. Alba, “Design issues in a multiobjective cellular genetic algorithm,” in Proc. Fourth Int. Conf. Evol. Multi-Criterion Optimization (EMO), Lecture Notes in Computer Science 4403. 2007, pp. 126–140. [26] K. Praditwong and X. Yao, “How well do multiobjective evolutionary algorithms scale to large problems,” in Proc. IEEE Congr. Evol. Comput. (CEC), Sep. 2007, pp. 3959–3966. [27] R. C. Purshouse and P. J. Fleming, “On the evolutionary optimization of many conflicting objectives,” IEEE Trans. Evol. Algorithms, vol. 11, no. 6, pp. 770–784, Dec. 2007. [28] M. Reyes and C. A. Coello Coello, “Improving PSO-based multiobjective optimization using crowding, mutation and -dominance,” in Proc. Third Int. Conf. Evol. Multicriterion Optimization (EMO), Lecture Notes in Computer Science 3410. Berlin, Germany: Springer, 2005, pp. 509– 519. [29] D. K. Saxena and K. Deb, “Non-linear dimensionality reduction procedures for certain large-dimensional multiobjective optimization problems: Employing correntropy and a novel maximum variance unfolding,” in Proc. Fourth Int. Conf. Evol. Multicriterion Optimization (EMO), Lecture Notes in Computer Science 4403. Matsushima, Japan, Mar. 2007, pp. 772–787. [30] D. J. Sheskin, Handbook of Parametric and Nonparametric Statistical Procedures, 4th ed. New York: Chapman & Hall/CRC Press, 2007. [31] D. Thain, T. Tannenbaum, and M. Livny, “Distributed computing in practice: The condor experience,” Concurrency Comput. Practice Experience, vol. 17, nos. 2–4, pp. 323–356, Feb. 2005. [32] D. A. Van Veldhuizen and G. B. Lamont, “Multiobjective evolutionary algorithm research: A history and analysis,” Dept. Elec. Comput. Eng., Graduate School Eng., Air Force Inst. Technol., Wright-Patterson Afb, OH, Tech. Rep. TR-98-03, 1998. [33] Q. Zhang and H. Li, “MOEA/D: A multiobjective evolutionary algorithm based on decomposition,” IEEE Trans. Evol. Comput., vol. 8, no. 11, pp. 712–731, 2008. [34] E. Zitzler, K. Deb, and L. Thiele, “Comparison of multiobjective evolutionary algorithms: Empirical results,” Evol. Comput., vol. 8, no. 2, pp. 173–195, 2000. [35] E. Zitzler, M. Laumanns, and L. Thiele, “SPEA2: Improving the strength pareto evolutionary algorithm,” in Proc. Evol. Methods Design Optimization Control Applicat. Ind. Problems (EUROGEN), Athens, Greece, 2002, pp. 95–100. [36] E. Zitzler and L. Thiele, “Multiobjective evolutionary algorithms: A comparative case study and the strength pareto approach,” IEEE Trans. Evol. Comput., vol. 3, no. 4, pp. 257–271, Nov. 1999. Juan J. Durillo (S’08) received the Engineering degree in computer science from the University of Málaga, Málaga, Spain, in 2006. He is currently pursuing the Ph.D. degree from the Departamento de Lenguajes y Ciencias de la Computaci´ on, University of M´ alaga. His research interests include multiobjective optimization and the design of new models for parallel multiobjective metaheuristics. Antonio J. Nebro received the M.S. and Ph.D. degrees in computer science from the University of M´ alaga, M´ alaga, Spain, in 1992 and 1999, respectively. He is currently an Associate Professor of Computer Science with the Departamento de Lenguajes y Ciencias de la Computaci´ on, University of M´ alaga. Since 1999, he has published more than ten scientific papers in international indexed journals, more than 20 papers in international conferences, and coauthored 12 book chapters. His current research interests include the design and implementation of parallel evolutionary algorithms, multiobjective optimization with metaheuristics, grid computing applied to metaheuristic techniques, and the application of these techniques to real-world problems. Carlos A. Coello Coello (SM’04) received the B.S. degree in civil engineering from the Universidad Autónoma de Chiapas, Chiapas, Mexico, in 1991, and the M.S. and Ph.D. degrees in computer science from Tulane University, New Orleans, LA, in 1993 and 1996, respectively. He is currently a Professor with the Department of Computer Science, Center of Research and Advanced Studies, Mexico City, Mexico. He has authored and coauthored over 200 technical papers and several book chapters, which report over 3600 citations (half of them in the Institute for Scientific Information Web of Science). He coauthored the book Evolutionary Algorithms for Solving Multiobjective Problems (2nd ed., New York: Springer, 2007). Dr. Coello is currently an Associate Editor of the IEEE Transactions on Evolutionary Computation and serves in the Editorial Board of 11 other international journals. He also chairs the IEEE Computationally Intelligent Society Task Force on Multiobjective Evolutionary Algorithms. He received the 2007 National Research Award from the Mexican Academy of Sciences in the area of Exact Sciences. He is a member of the Association for Computing Machinery, Sigma Xi, and the Mexican Academy of Sciences. José García-Nieto received the Engineering degree in computer science and the Master degree in artificial intelligence from the University of Málaga, Málaga, Spain, in 2006 and 2008, respectively. He is currently working toward the Ph.D. degree in optimization algorithms at the Departamento de Lenguajes y Ciencias de la Computaci´ on, University of M´ alaga. His research interests include swarm intelligence and general metaheuristics, and their application to complex real-world problems in the domains of telecommunications, bioinformatics, real-parameter, and multiobjective optimization. Francisco Luna received the Ph.D. degree in computer science from the University of M´ alaga, M´ alaga, Spain, in 2008. He is currently an Assistant Researcher of computer science with the Departamento de Lenguajes y Ciencias de la Computaci´ on, University of M´ alaga. He has published more than ten scientific papers in international indexed journals and more than 20 papers in international conferences. His current research interests include the design and implementation of parallel and multiobjective metaheuristics, and their application to solve complex problems in the domain of telecommunications and combinatorial optimization. Enrique Alba received the Ph.D. degree in computer science on parallel and distributed genetic algorithms, in 1999. He is currently a Full Professor of Computer Science with Departamento de Lenguajes y Ciencias de la Computaci´ on, University of M´ alaga, M´ alaga, Spain. Since 1999, he has published more than 30 scientific papers in international indexed journals and more than 150 papers in international conferences. His research interests include the design and application of evolutionary algorithms and other metaheuristic procedures to real problems, including telecommunications, combinatorial optimization, and bioinformatics. The main focus of all his work is on parallelism. Dr. Alba has won several awards for his research quality.