On the Velocity Update in Multi-Objective Particle Swarm Optimizers

On the Velocity Update in Multi-Objective Particle Swarm Optimizers Juan J. Durillo, Antonio J. Nebro, Jos´eGarcı´a-Nieto, and Enrique Alba Dept. Lenguajes y Ciencias de la Computaci´on, ETSI Inform´atica, University of M´alaga, Campus de Teatinos, 29071 M´alaga, Spain {durillo,antonio,jnieto,eat}@lcc.uma.es Summary. Since its appearance, Particle Swarm Optimization (PSO) has become a very popular technique for solving optimization problems because of both its simplicity and its fast convergence properties. In the last few years there has been a variety of proposals for extending it to handle with multiples objectives. Although many of them keep the same properties of the original algorithm, they face difficulties when tackling the optimization of some multi-modal problems, i.e., those having more than one suboptimal front of solutions. Recent studies have shown that this disadvantage could be related to the velocity of the particles: uncontrolled high velocities may have no effect in particles movements. While many of the contributions on the specialized literature have focused on the selection of the leaders of the swarm, studies about different schemes for controlling the velocity of the particles are scarce in the multi-objective domain. In this work, we study different mechanisms in order to update the velocity of each particle with the idea of enhancing the search capabilities of multi-objective PSO algorithms. Our experiments show that some modifications help to over-coming the difficulties observed in previous proposals when dealing with hard optimization problems. 3.1 Introduction Particle Swarm Optimization (PSO) has become a popular algorithm due to its relative simplicity and competitive performance when solving a wide range of optimization problems in the continuous domain [11]. A considerable number of these problems has to optimize more than one objective function at the same time which are in conflict with respect to each other, so many proposals of Multi-Objective PSO (MOPSO) techniques have been developed [16]. In the survey presented in [16], the authors considered that the main characteristics of all the existing MOPSOs are the existence of an external archive of non-dominated solutions, the selection strategy for non-dominated solutions as leaders for guiding the swarm, the neighborhood topology, and the existence or not of a mutation operator. In this chapter, our ap-proach is to study a different feature: the scheme for controlling the velocity of the particles. This issue has been studied in single-objective optimization [1, 2, 15, 18], and we are interested in investigating how it can affect the search effectivity in MOPSOs. In our previous research, we analyzed the performance of six representative MOPSO metaheuristics in [8], concluding that a variant of the OMOPSO algorithm [17] provided the best overall performance over a set of 21 benchmark problems. OMOPSO also shown to be very fast in a comparison with other state- of-the-art multi-objective metaheuristics [14]; however, we realized that OMOPSO had difficulties when facing the solution of some multi-modal problems. We deeply studied this issue in [8], finding that the velocity of the particles in MOPSOs can become too high, and hence resulting in erratic movements towards the upper and lower limits of the positions of the particles. This is an example of the so-called swarm explosion [2], and we found out that it can be prevented by using a velocity constriction mechanism [8]. Our motivation then is, taking OMOPSO as our baseline MOPSO algorithm, and to study different velocity update schemes in order to have an insight of the potential improvements in the search capacity they can lead to MOPSO metaheuristics. The contributions of this chapter can be summarized as follows: •We propose four velocity schemes to be applied to MOPSO algorithms. •The resulting techniques are evaluated on a benchmark of 21 problems (those comprising the ZDT [20], DTLZ [5], and WFG [10] problem families). •Three quality indicators are used to assess the performance of the algorithms (additive unary epsilon (I1 ε +) [13],spread ( Δ ) [4], and hypervolume (HV) [19]). •We use a statistical analysis methodology to ensure the significance of the obtained results. The remainder of this chapter is organized as follows. Section 3.2 includes basic background about PSO and MOPSO algorithms. Section 3.3 is aimed at describing OMOPSO, the baseline approach considered, and the different velocity schemes that we have applied, leading to four resulting versions of OMOPSO. Section 3.4 is devoted to the experimentation, including the benchmark problems, the quality indicators, the parameter setting, and the methodology adopted in the statistical tests. In Section 3.5, we analyze the obtained results regarding the three used quality indicators. The results are discussed in Section 3.6. Finally, Section 3.7 contains the conclusions and some possible lines for future work. 3.2 PSO Background In a PSO algorithm, each potential solution to the problem is called particle and the population of solutions is called swarm. A basic PSO updates the particle xiat the generation twith the formula: xi(t)=xi(t−1)+vi(t)(3.1) where the factor vi(t)is known as velocity and is given by vi(t)= χ [w·vi(t−1)+C1·r1·(xpi−xi)+C2·r2·(xgi−xi)] (3.2) In this formula, xpiis the best solution that xihas viewed (pbest), xgiis the best particle (gbest, also known as the leader) that the entire swarm has viewed, wis the inertia weight of the particle and controls the trade-off between global and local experience, r1and r2are two uniformly distributed random numbers in the range [0,1], and the parameters C1and C2are specific parameters which control the effect of the personal and global best particles. χ is a constriction coefficient introduced to control the particle’s velocity [2]. Algorithm 3.1. Pseudocode of a general PSO algorithm. 1: initializeSwarm() 2: locateLeader() 3: generation = 0 4: while generation <maxGenerations do 5: for each particle do 6: updatePosition() // flight (Formulas 1 and 2) 7: evaluation() 8: updatePbest() 9: end for 10: updateLeader() 11: generation ++ 12: end while Algorithm 3.1 describes the pseudocode of a general single-objective PSO. The algorithm starts by initializing the swarm (Line 1), which includes both the positions and velocities of the particles. The correspondingpbest of each particle is initialized, as well as the leader (Line 2). Then, during a maximum number of iterations, each particle flies through the search space updating its position (Line 6). Then, it is evaluated (Line 7), and its pbest is also calculated (Lines 6-8). At the end of each iteration, the leader is updated. As commented before, the leader is usually the best particle in the swarm (i.e., gbest). However, it can be a different particle depending on the social structure of the swarm (i.e., the topology of the neighborhood of each particle) [12]. To apply a PSO algorithm in multi-objective optimization the previous scheme has to be modified to cope with the fact that the solution of a problem with multiple objectives is not a single one but a set of non-dominated solutions. Therefore, issues that have to be considered now are [16]: 1. How to select the particles to be used as leaders? 2. How to retain the non-dominated solutions found during the search? 3. How to maintain diversity in the swarm in order to avoid convergence to single solutions? The pseudocode of a general MOPSO is included in Algorithm 3.2. After initialising the swarm (Line 1), the typical approach is to use an external archive to store the leaders, which are taken from the non-dominated particles in the swarm. After initializating the leaders archive (Line 2), some kind of quality measure has to be calculated (Line 3) for all the leaders to select usually one leader for each particle of the swarm. In the main loop of the algorithm, the flight of each particle is performed after a leader has been selected (Lines 7-8) and, optionally, a mutation or turbulence operator can be applied (Line 9); then, the particle is evaluated and its corresponding pbest is updated (Lines 10-11). After each iteration, the set of leaders is updated and the quality measure is calculated again (Lines 13-14). After the termination condition, the archive is returned as the result of the search. For further details about the operations contained in the MOPSO pseudocode, please refer to [16]. Algorithm 3.2. Pseudocode of a general MOPSO algorithm. 1: initializeSwarm() 2: initializeLeadersArchive() 3: determineLeadersQuality() 4: generation = 0 5: while generation <maxGenerations do 6: for each particle do 7: selectLeader() 8: updatePosition() // flight (Formulas 1 and 2) 9: mutation() 10: evaluation() 11: updatePbest() 12: end for 13: updateLeadersArchive() 14: determineLeadersQuality() 15: generation ++ 16: end while 17: returnArchive() 3.3 Velocity Schemes and Resulting Algorithms In this section we describe the MOPSO variants we have developed for our study. We start by giving details of OMOPSO, the baseline approach. 3.3.1 OMOPSO As commented in the introduction, our base MOPSO algorithm is a variant of OMOPSO (Optimized MOPSO), proposed by Reyes-Sierra and Coello Coello in [17]. This algorithm is characterized by using the crowding distance of NSGAII to filter out leader solutions and the combination of two 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, but in our experiments we always discard this feature, being the leaders archive the result of the execution of the technique. The velocity scheme in OMOPSO, following Equation 3.2, is defined as follows: •The inertia weight wis a uniformly distributed random number in the range [0,1]. •The coefficientsC1andC2are two uniformly distributed random numbers in the range [1.5,2.0]. •The constriction coefficient χ takes the value 1.0. In addition, we introduce a mechanism in such a way that the accumulated velocity of each variable j(in each particle) is further bounded by means of the following velocity constriction equation: vi,j(t)=⎧ ⎪ ⎨ ⎪ ⎩ deltajif vi,j(t)>deltaj −deltajif vi,j(t)≤−deltaj vi,j(t)otherwise (3.3) where deltaj=(upper limitj−lower limitj) 2(3.4) After applying Equation 3.1, OMOPSO checks whether the resulting positions are out of the bounds of the variables of the problem. In that case, the positions are assigned the corresponding upper or lower bound value; additionally, the direction of the velocity is reversed by multiplying it by −1.0. Finally, OMOPSO applies a combination of uniform and non-uniform mutation to the particle swarm (uniform mutation to the 30% of the swarm, non-uniform to other 30%, and no mutation to the rest of the particles). Once we have defined our base MOPSO algorithm, we present next the four velocity schemes which will lead to the same number of algorithms. Each scheme affects to each different component of Equation 3.2: the constriction coefficient χ , the inertia weigth w, the coefficients C1and C2, and the component vi(t−1)(the current velocity). 3.3.2 SMPSO Our previous research in [8] indicated that OMOPSO had difficulties when solving some multi-modal problems (e.g., ZDT4, DTLZ1, and DTLZ3). Our analysis of this issue showed that including a constriction coefficient similar to the one proposed in [2] the resulting algorithm could successfully solve these problems. We called this algorithm SMPSO (Speed contrained Multi-objective PSO). The constriction coeffient applied in SMPSO is defined as follows: χ =2 2− ϕ − ϕ 2−4 ϕ (3.5) where ϕ =C1+C2if C1+C2>4 4ifC1+C2≤4(3.6) Besides using this velocity scheme, the coeffientsC1and C2are random numbers in the range [1.5,2.5]; the range used in OMOPSO, [1.5,2.0]would lead Equation 3.6 to always return a value of 4. 3.3.3 MOPSO TVAC Ratnaweera et al. proposedin [15] the use of Time-Varying Acceleration Coefficients (TVAC) with the idea of enhancing the search in the early part of the optimization and to encourage the particles to converge toward the global optima at the end of the search process. This can be carried out by linearly changing the coefficients C1and C2through the time (number of iterations). They suggest the following definitions of C1and C2: C1=(C1f−C1i)iter MAXITR +C1i(3.7) C2=(C2f−C2i)iter MAXITR +C2i(3.8) where C1i,C1f,C2i,andC2fare constants, iter is the current iteration number and MAXITR is the maximum number of iterations of the PSO algorithm. As suggested in [15], we use the values of C1andC2changing from 2.5 to 0.5 and from 0.5 to 2.5, respectively. The inertia weight wis not considered, so it takes the value 1.0. The resulting algorithm after applying this velocity scheme to OMOPSO is named MOPSO TVAC. In [15], a mutationoperator is introduced,but we have omitted it due to the fact that OMOPSO already includes its own mutation mechanism. 3.3.4 MOHPSO The self-Organizing Hierarchical Particle Swarm Optimization optimizer (HPSO) was also proposed in [15]. The authors of this work observed that in the absence of the previous velocity term, particles rapidly rush to a local optimum and stagnate due to the lack of momentum.To cope this issue, they proposed a reinitialization scheme proportional to the maximum allowable velocity. Taking these ideas, the MOHPSO algorithm is characterized by the following equation defining the velocity: vi(t)=C1·r1·(xpi−xi)+C2·r2·(xgi−xi)(3.9) If a given velocity element vi,j(t)gets 0, then it is reinitialized according to: vi,j(t)=rand1·deltajif rand2<0.5 rand3·−deltajif rand2>=0.5(3.10) where rand1, rand2, and rand3 are separately generated uniformly distributed random numbers in [0,1]. The limit deltajresults from Equation 3.4. 3.3.5 MOPSO TVIW The last variant we have considered is based on the Time-Varying Inertia Weight (TVIW) proposed by Shi and Eberhart in [18]. They found that the performance of a PSO method could improve by linearly varying the inertia weigh w. We adopt this scheme in OMOPSO and the result is MOPSO TVIW. The inertia weight is defined as follows: w=(w1−w2)MAXITR−iter MAXITR +w2(3.11) where w1and w2are the initial and final values of the inertia weight, respectively, iter is the current iteration number, and MAXITR is the maximum number of iterations of the PSO algorithm. The study carried out in [18] yielded that the most promising results were obtained by varying wfrom 0.9 (w1) at the beginning of the search to 0.4 (w2) at the end for most of the studied problems; here, we use the same values of w1and w2. 3.4 Experimentation In this section we explain the benchmark problems used to evaluate the algorithms, the quality indicators used to assess their performance, the parameter settings used, and the statistical tests carried out. 3.4.1 Benchmark Problems Here, we describe the different sets of problems addressed in this work. These problems are well-known, and they have been used in many studies in this area. The problems families are the following: •Zitzler-Deb-Thiele (ZDT): This benchmark is composed of five bi-objective problems [20]: ZDT1 (convex), ZDT2 (nonconvex), ZDT3 (nonconvex, disconnected), ZDT4 (convex, multimodal), and ZDT6 (nonconvex, nonuniformly spaced). These problems are scalable according to the number of decision variables. •Deb-Thiele-Laumanns-Zitzler (DTLZ): The problems of this family are scalable both in the number of variables and objectives [5]. It is composed of the following seven problems: DTLZ1 (linear), DTLZ2-4 (nonconvex), DTLZ5-6 (degenerate), and DTLZ7 (disconnected). •Walking-Fish-Group (WFG): This set is composed of nine problems, WFG1 - WFG9, that have been constructed using the WFG toolkit [10]. The properties of these problems are detailed in Table 3.1. They all are scalable both in the number of variables and the number of objectives. In this work we have used the bi-objective formulation of the DTLZ and WFG problem families. A total of 21 MOPs are used to evaluate the six metaheuristics. 3.4.2 Quality Indicators To assess the search capabilities of multi-objective metaheuristics on the test problems, two different issues are normally taken into account: the distance between the Table 3.1. Properties of the MOPs created using the WFG toolkit Problem Separability Modality Bias Geometry WFG1 separable uni polynomial, flat convex, mixed WFG2 non-separable f1uni, f2multi no bias convex, disconnected WFG3 non-separable uni no bias linear, degenerate WFG4 non-separable multi no bias concave WFG5 separable deceptive no bias concave WFG6 non-separable uni no bias concave WFG7 separable uni parameter dependent concave WFG8 non-separable uni parameter dependent concave WFG9 non-separable multi, deceptive parameter dependent concave solution set generated by the proposed algorithm to the optimal Pareto front should be minimized (convergence) and the spread of solutions should be maximized in order to obtain as smooth and uniform a distribution of solutions as possible (diversity). To measure these two criteria it is necessary to know the exact location of the optimal Pareto front; the benchmark problems used in this work have known Pareto fronts. The quality indicators can be classified into three categories depending on whether they evaluate the closeness to the Pareto front, the diversity in the solutions obtained, or both [3]. We have adopted one indicator of each type. •Unary Epsilon Indicator (I1 ε +). This indicator was proposed by Zitzler et al. [21] and makes direct use of the principle of Pareto-dominance. Given an approximation set, A, of a problem, the I1 ε +indicator is a measure of the smallest distance one would need to translate every point in A so that it dominates the optimal Pareto front of the problem. More formally, given z1=(z1 1,..., z1 n)and z2=(z2 1,..., z2 n),wherenis the number of objectives: I1 ε +(A)=inf ε ∈R∀z2∈ParetoOptimal Front ∃z1∈A:z1≺ ε z2(3.12) where, z1≺ ε z2if and only if ∀1≤i≤n:z1 i< ε +z2 i. •Spread ( Δ ). The Spread indicator [4] measures the extent of spread achieved among the obtained solutions. This indicator (illustrated in Fig. 3.1) is defined as: Δ =df+dl+∑N−1 i=1di−¯ d df+dl+(N−1)¯ d,(3.13) where diis the Euclidean distance between consecutive solutions, ¯ dis the mean of these distances, and dfand dlare the Euclidean distances to the extreme (bounding) solutions of the optimal Pareto front in the objective space (see [4] for the details). Δ takes a value of zero for an ideal distribution, pointing out a Fig. 3.1. Calculating the Spread quality indicator f1 f2 Pa r eto-optimalf r ont W A B C 1 0 Fig. 3.2. The hypervolume enclosed by the non-dominated solutions perfect spread out of the solutions in the Pareto front. We apply this indicator after a normalization of the objective function values. •Hypervolume (HV). The HV indicator calculates the volume, in the objective space, covered by members of a non-dominated set of solutions Qfor problems where all objectives are to be minimized [19]. In the example depicted in Fig. 3.2, the HV is the region enclosed within the discontinuous line, where Q={A,B,C}(in the figure, the grey area represents the objective space that has been explored). Mathematically, for each solution i∈Q, a hypercube voliis constructed with a reference point Wand the solution ias the diagonal corners of the hypercube. The reference point can be found simply by constructing a vector of worst objective function values. Thereafter, a union of all hypercubes is found and its hypervolume (HV) is calculated: HV =volume⎛ ⎝ |Q|  i=1 voli⎞ ⎠.(3.14) Table 3.6. SMHPSO vs SMPSO vs MOHPSO Problem SMPSO MOHPSO I1 ε + Δ HV I1 ε + Δ HV ZDT1       ZDT2       ZDT3       ZDT4 –    ZDT6     – DTLZ1 –     DTLZ2 – – –    DTLZ3 – –     DTLZ4 – – – – –  DTLZ5 – – –    DTLZ6       DTLZ7       WFG1       WFG2       WFG3 –     WFG4       WFG5   –  – WFG6 –    WFG7 –    WFG8      – WFG9 –    corresponding algorithm is significantly better than SMHPSO in the indicator specified by the column containing that symbol, a symbol “” means that SMHPSO is perform better, and a symbol “–” means that no statistical differences were found. The results in Table 3.6 indicate that SMHPSO outperforms SMPSO in the WFG family and, in a similar way, it improves the values obtained by MOHPSO also in the DTLZ problems. Thus, we can state that SMHPSO has improved the results of SMPSO and MOHPSO in those problems in which these two algorithms encounter difficulties, but it has not been able to improve the results in those problems in which SMPSO and MOHPSO are the best algorithms. This leads usto proposea researchline, relatedto investigatehow to use theproper velocityupdateschemein orderto design a MOPSO ableof improvethe performance of those studied in this chapter. Besides analyzingother update strategies (e.g., those proposedin [1]),finding outhybridapproachesaswell as designingadaptativemechanisms to vary the update scheme during the search are promising ideas. 3.7 Conclusions and Future Work In this chapter we have studied the effect of applying different velocity schemes to OMOPSO, a multi-objective PSO algorithm which has proven to be competitive against a set of state-of-the-art multi-objective optimizers in previous works. Concretely, we have developed and evaluated four alternatives called SMPSO, MOPSO TVAC, MOHPSO, and MOPSO TVIW, each of them characterized by modifying a different component of the formula defining the velocity scheme in a PSO. As bechmark problems, we have used 21 instances corresponding to the well-known ZDT, DTLZ, and WFG test suites. The obtained results have shown that, in the context of the problems, the quality indicators, and the parameter settings considered, SMPSO, MOPSO TVAC, and MOPSO TVIW improves the results obtained by OMOPSO in the ZDT and DTLZ families, whereas MOHPSO is well suited for solving the problems composing the WFG benchmark. We have also carried out a first attempt to hybridize two of the most promising MOPSOs trying to combine their search habilities but, although the results are promising, it is a matter of further research. Other future research topics in this line are related to study how the velocity update scheme may affect the MOPSOs concerning issues such as their speed to converge faster to the true Pareto front, or their ability to solve scalable problems in the number of variables and/or objectives. Acknowledgement. Authors acknowledge funds from the Spanish Ministry of Sciences and Innovation European FEDER under contract TIN2008-06491-C04-01 (M* project, available at http://mstar.lcc.uma.es) and CICE, Junta de Andaluc´ıa under contract P07-TIC-03044 (DIRICOM project, http://diricom.lcc.uma.es). Juan J. Durillo is supported by grant AP-2006-003349 from the Spanish Ministry of Education and Science. Jos´eGarc´ıa-Nieto is supported by grant BES-2009-018767 from the Spanish Ministry of Sciences and Innovation. References 1. Bui, L., Soliman, O., Abbass, H.: A modified strategy for the constriction factor in particle swarm optimization. In: Randall, M., Abbass, H.A., Wiles, J. (eds.) ACAL 2007. LNCS (LNAI), vol. 4828, pp. 333–344. Springer, Heidelberg (2007) 2. Clerc, M., Kennedy, J.: The particle swarm - explosion, stability, and convergence in a multidimensional complex space. IEEE Transactions on Evolutionary Computation 6(1), 58–73 (2002) 3. Deb, K.: Multi-Objective Optimization Using Evolutionary Algorithms. John Wiley & Sons, Chichester (2001) 4. Deb, K., Pratap, A., Agarwal, S.,Meyarivan, T.: A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II. IEEE Transactions on Evolutionary Computation 6(2), 182–197 (2002) 5. Deb, K., Thiele, L., Laumanns, M., Zitzler, E.: Scalable Test Problems for Evolutionary Multiobjective Optimization. In: Abraham, A., Jain, L., Goldberg, R. (eds.) Evolutionary Multiobjective Optimization. Theoretical Advances and Applications, pp. 105–145. Springer, USA (2005) 6. Demˇsar, J.: Statistical Comparisons of Classifiers over Multiple Data Sets. J. Mach. Learn. Res. 7, 1–30 (2006) 7. Durillo, J., Nebro, A., Luna, F., Dorronsoro, B., Alba, E.: jMetal: a Java Framework for Developing Multi-objective Optimization Metaheuristics. Tech. Rep. ITI-2006-10, Departamento de Lenguajes y Ciencias de la Computaci´on, University of M´alaga, E.T.S.I. Inform´atica, Campus de Teatinos (2006) 8. Durillo, J., Garc´ıa-Nieto, J., Nebro, A., Coello Coello, C., Luna, F., Alba, E.: Multiobjective particle swarm optimizers: An experimental comparison. Accepted for publication in EMO 2009 (2009) 9. Hochberg, Y., Tamhane, A.C.: Multiple Comparison Procedures. Wiley, Chichester (1987) 10. Huband, S., Barone, L., While, R.L., Hingston, P.: A scalable multi-objective test problem toolkit. In: Coello Coello, C.A., Hern´andez Aguirre, A., Zitzler, E. (eds.) EMO 2005. LNCS, vol. 3410, pp. 280–295. Springer, Heidelberg (2005) 11. Kennedy, J., Eberhart, R.: Particle swarm optimization. In: Proceedings of IEEE International Conference on Neural Networks, pp. 1942–1948 (1995) 12. Kennedy, J., Eberhart, R.C.: Swarm Intelligence. Morgan Kaufmann Publishers, San Francisco (2001) 13. Knowles, J., Thiele, L., Zitzler, E.: A Tutorial on the Performance Assessment of Stochastic Multiobjective Optimizers. Tech. Rep. 214, Computer Engineering and Networks Laboratory (TIK), ETH Zurich (2006) 14. Nebro, A.J., Durillo, J.J., Coello Coello, C., Luna, F., Alba, E.: A study of convergence speed in multi-objective metaheuristics. In: Rudolph, G., Jansen, T., Lucas, S., Poloni, C., Beume, N. (eds.) PPSN 2008. LNCS, vol. 5199, pp. 763–772. Springer, Heidelberg (2008) 15. Ratnaweera, A., Halgamuge, S., Watson, H.: Self-organizing hierarchical particle swarm optimizer with time-varying acceleration coefficients. International Journal of Computational Intelligence Research 8(3), 240–255 (2004) 16. Reyes-Sierra, M., Coello Coello, C.: Multi-Objective Particle Swarm Optimizers: A Survey of the State-of-the-Art. International Journal of Computational Intelligence Research 2(3), 287–308 (2006) 17. Reyes-Sierra, M., Coello Coello, C.A.:Improving PSO-Based Multi-objective Optimization Using Crowding, Mutation and ε -Dominance. In: Coello Coello, C.A., Hern´andez Aguirre, A., Zitzler, E. (eds.) EMO 2005. LNCS, vol. 3410, pp. 505–519. Springer, Heidelberg (2005) 18. Shi, Y., Eberhart, R.: Empirical study of particle swarm optimization. In: Proceedings of the 1999 Congress on Evolutionary Computation, 1999. CEC 1999, pp. 1945–1950 (1999) 19. Zitzler, E., Thiele, L.: Multiobjective Evolutionary Algorithms: A Comparative Case Study and the Strength Pareto Approach. IEEE Transactions on Evolutionary Computation 3(4), 257–271 (1999) 20. Zitzler, E., Deb, K., Thiele, L.: Comparison of multiobjective evolutionary algorithms: Empirical results. IEEE Transactions on Evolutionary Computation 8(2), 173–195 (2000) 21. Zitzler, E., Thiele, L., Laumanns, M., Fonseca, C., da Fonseca, V.G.: Performance assessment of multiobjective optimizers: An analysis and review. IEEE Transactions on Evolutionary Computation 7, 117–132 (2003)