IGD Indicatorbased Evolutionary Algorithm for Manyobjective Optimization Problems
Abstract
Inverted Generational Distance (IGD) has been widely considered as a reliable performance indicator to concurrently quantify the convergence and diversity of multi and manyobjective evolutionary algorithms. In this paper, an IGD indicatorbased evolutionary algorithm for solving manyobjective optimization problems (MaOPs) has been proposed. Specifically, the IGD indicator is employed in each generation to select the solutions with favorable convergence and diversity. In addition, a computationally efficient dominance comparison method is designed to assign the rank values of solutions along with three newly proposed proximity distance assignments. Based on these two designs, the solutions are selected from a global view by linear assignment mechanism to concern the convergence and diversity simultaneously. In order to facilitate the accuracy of the sampled reference points for the calculation of IGD indicator, we also propose an efficient decompositionbased nadir point estimation method for constructing the Utopian Pareto front which is regarded as the best approximate Pareto front for realworld MaOPs at the early stage of the evolution. To evaluate the performance, a series of experiments is performed on the proposed algorithm against a group of selected stateoftheart manyobjective optimization algorithms over optimization problems with , , and objective. Experimental results measured by the chosen performance metrics indicate that the proposed algorithm is very competitive in addressing MaOPs.
Nadir point, inverted generational distance, linear assignment problem, manyobjective evolutionary optimization algorithm.
1 Introduction
\IEEEPARstartManyobjective Optimization Problems (MaOPs) refer to the optimization tasks involving (i.e., ) conflicting objectives to be optimized concurrently [1]. Generally, an MaOP is with the mathematic form represented by (1)
(1) 
where is the feasible search space for the decision variables , and is the corresponding objective vector including objectives which maps the dimensional decision space to the dimensional objective space . Without loss of generality, is assumed to be minimized since the maximization problems can be transformed into the minimization problems due to the duality principle. Because of the conflicting nature in the objective functions, there is no single perfect solution for f(x), but a set of tradeoff solutions which form the Pareto Set (PS) in the decision space and the corresponding Pareto Front (PF) in the objective space.
Optimization algorithms for addressing an MaOP aim at searching for a set of uniformly distributed solutions which are closely approximating the PF. Because the MaOPs widely exist in diverse realworld applications, such as policy management in land exploitation with objective [2] and calibration of automotive engine with objective [3], to name a few, various algorithms for solving MaOPs have been developed. Among these algorithms, the evolutionary paradigms are considerably preferable due to their populationbased metaheuristic characteristics obtaining a set of quality solutions in a single run.
During the past decades, various MultiObjective Evolutionary Algorithms (MOEAs), such as elitist Nondominated Sorting Genetic Algorithm (NSGAII) [4], advanced version of Strength Pareto Evolutionary Algorithm (SPEA2) [5], among others, have been proposed to effectively solve MultiObjective Optimization Problems (MOPs). Unfortunately, these MOEAs do not scale well with the increasing number of objectives, mainly due to the loss of selection pressure. To be specific, the number of nondominated solutions in MaOPs accounts for a large proportion of the current population because of the dominance resistance phenomenon [6] caused by the curse of dimensionality [7], so that the traditional elitism mechanism based on Paretodomination cannot effectively differentiate which solutions should survive into the next generation. As a result, the densitybased diversity promotion mechanism is considered the sole mechanism for mating and environmental selections [8]. However, the solutions with good diversity in MaOPs are generally not only distant from each other but also away from the PF. Consequently, the evolution with the solutions generated by the activated diversity promotion is stagnant or even far away from the PF [9]. To this end, various ManyObjective Evolutionary Algorithms (MaOEAs) specifically designed for addressing MaOPs have been proposed in recent years.
Generally, these MaOEAs can be divided into four different categories. The first category covers the algorithms employing reference prior to enhancing the diversity promotion which in turn improve the convergence. For example, the MaOEA using referencepointbased nondominated sorting approach (NSGAIII) [10] employs a set of reference vectors to assist the algorithm to select solutions which are close to these reference vectors. Yuan et al. [11] proposed the reference linebased algorithm which not only adopted the diversity improvement mechanism like that in NSGAIII but also introduced convergence enhancement scheme by measuring the distance between the origin to the solution projections on the corresponding reference line. In addition, a reference linebased estimation of distribution algorithm was introduced in [12] for explicitly promoting the diversity of an MaOEA. Furthermore, an approach (RVEA) was presented in [13] to adaptively revise the reference vector positions based on the scales of the objective functions to balance the diversity and convergence.
The second category refers to the decompositionbased algorithms which decompose an MaOP into several singleobjective optimization problems, such as the MOEA based on Decomposition (MOEA/D) [14] which was initially proposed for solving MOPs but scaled well for MaOPs. Specifically, MOEA/D transformed the original MOP/MaOP with objectives into a group of singleobjective optimization problems, and each subproblem was solved in its neighboring region which constrained by their corresponding reference vectors. Recently, diverse variants [15, 16, 17, 18, 19, 20] of MOEA/D were proposed for improving the performance much further.
The third category is known as the convergence enhancementbased approaches. More specifically, the traditional Pareto dominance comparison methods widely utilized in MOEAs are not effective in discriminating populations with good proximity in MaOPs. A natural way is to modify this comparison principle to promote the selection mechanism. For example, the dominance method [21] employed a relaxed factor to compare the dominance relation between solutions; Pierro et al. [22] proposed the preference order ranking approach to replace the traditional nondominated sorting. Furthermore, the fuzzy dominance methods [23, 24] studied the fuzzification of the Paretodominance relation to design the ranking scheme to select promising solutions; the optimality paradigm was proposed in [25] to pick up solutions whose objectives were with the same importance by considering their objective value improvements. In addition, Yang et al. [26] proposed the gridbased approach to select the solutions that have the higher priority of dominance, and control the proportion of Paretooptimal solutions by adjusting the grid size. Meanwhile, Antonio et al. [27] alternated the achievement function and the indicator method to improve the performance of MOEA in solving MaOPs. In [28], a modification of density estimation, termed as shiftbased density estimation, was proposed to make the dominance comparison better suited for solving MaOPs. Furthermore, the favorable convergence scheme was proposed in [29] to improve the selection pressure in mating and environmental selections. Recently, a knee pointbased algorithm (KnEA) [30] was presented as a secondary selection scheme to enhance the selection pressure. In summary, these algorithms introduced new comparison methods, designed effective selection mechanisms, or relaxed the original comparison approach to improve the selection pressure in addressing MaOPs.
The fourth category is known as the indicatorbased methods. For instance, several MOEAs based on the hypervolume (HV) were proposed in [31, 32, 33], however their major disadvantages were the costly overhead in calculating the HV values especially in solving MaOPs. To this end, Bader and Zitzler proposed the HypE method with the Monte Carlo simulation [34] to estimate the HV value. Consequently, the computational cost was largely lowered compared to its predecessors whose HV values were calculated exactly. In [35], an indicatorbased algorithm (EMOA) was proposed for solving biobjective optimization problems, and then extended further for triobjective problems [36]. Furthermore, Villalobos and Coello [37] integrated the indicator with the differential evolution [38] to solve MaOPs with up to objectives. Recently, an Inverse Generational Distance Plus (IGD) [39] indicatorbased evolutionary algorithm (IGDEMOA) was proposed in [40] for addressing MaOPs with no more than objectives. Basically, the IGD indicator is viewed as a variant of the Inverse Generational Distance (IGD) indicator.
Although the MaOEAs mentioned above have experimentally demonstrated their promising performance, major issues are easily to be identified in solving realworld applications. For example, it is difficult to choose the converting strategy of the MaOEAs from the second category, which motivates multiple variants [15, 41, 17, 42, 19, 20] to be developed further. In addition, the MaOEAs from the first and third categories only highlighted one of the characters in their designs (i.e., only the diversity promotion is explicitly concerned in the MaOEAs from the first category, and the convergence from the third category). However, both the diversity and convergence are concurrently desired by the MaOEAs. In this regard, some performance indicators, which are capable of simultaneously measuring the diversity and convergence, such as the HV and IGD indicators, are preferred to be employed for designing MaOEAs. However, the major issue of HV is its high computational complexity. Although the Monte Carlo simulation has been employed to mitigate this adverse impact, the calculation is still impracticable when the number of objectives is more than 10 [34], while the calculation of IGD is scalable without these deficiencies.
In this paper, an IGD indicatorbased ManyObjective Evolutionary Algorithm (MaOEA/IGD) has been proposed for effectively addressing MaOPs, and the contributions are outlined as follows:

A Decompositionbased Nadir Point Estimation method (DNPE) has been presented to estimate the nadir points to facilitate the calculation of IGD indicator. In DNPE, the estimation focuses only on the extreme point areas and transforms the computation of an objective optimization problem into singleobjective optimization problems. Therefore, less computational cost is required compared to its peer competitors (experiments are demonstrated in Subsection 4.5).

A comparison scheme for the nondominated sorting has been designed for improving the convergence of the proposed algorithm. In this scheme, the dominance relations of solutions are not obtained by the comparisons among all the solutions but the solutions to the reference points. Therefore, the computational complexity of the presented comparison scheme is significantly lessened compared to the traditional comparison means because the number of reference points is generally much less than that of the whole population.

Three types of proximity distance assignment mechanisms are proposed for the solutions according to their PF rank values, which make the solutions with good convergence in the same PF to have higher chances to be selected. Furthermore, these assignment mechanisms collectively assure the proposed IGD indicator to be Pareto compliance.

Based on the proposed dominance comparison scheme and the proximity distance assignments, the selection mechanism which is employed for the mating selection and the environmental selection is proposed to concurrently facilitate the convergence and the diversity.
The reminder of this paper is organized as follows. First, related works are reviewed, and the motivation of the proposed DNPE is presented in Section 2. Then the details of the proposed algorithm are documented in Section 3. To evaluate the performance of the proposed algorithm in addressing MaOPs, a series of experiments over scalable benchmark test suits are performed against stateoftheart MaOEAs, and their results are measured by commonly chosen performance metrics and then analyzed in Section 4. In addition, the performance of the proposed MaOEA/IGD is also demonstrated by solving a realworld application, and the performance of the proposed DNPE in nadir point estimation is investigated against its peer competitors. Finally, the proposed algorithm is concluded and the future works are illustrated in Section 5.
2 Related Works and Motivation
Literatures related to the nadir point estimation and IGD indicatorbased EAs are thoroughly reviewed in this section. Specifically, the Worst Crowded NSGAII (WCNSGAII) [43] and the Pareto Corner Search Evolutionary Algorithm (PCSEA) [44] would be reviewed and criticized in detail, because the insightful observations of the deficiencies of these two approaches naturally lead to the motivation of the proposed DNPE design. In addition, the IGDEMOA is reviewed as well to highlight the utilization of the IGD indicator and the reference points sampling for the calculation of IGD in the proposed MaOEA/IGD. Please note that all the discussions in this section are with the context of the problem formulation in (1).
2.1 Nadir Point Estimation Methods
According to literatures [45, 46], the approaches for estimating the nadir points can be divided into three categories including the surfacetonadir, edgetonadir, and extremepointtonadir schemes. In the surfacetonadir scheme, the nadir points are constructed from the current Paretooptimal solutions, and updated as the corresponding algorithms evolve towards the PF. MOEAs in [4, 47, 48, 49, 18] and MaOEAs in [10, 50, 13] belong to this category. However, these MOEAs are shown to perform poorly in MaOPs due to the curse of dimensionality [51]. In addition, the MaOEAs related methods are not suitable for the proposed algorithm because the MaOPs have been solved prior to the nadir point estimation, while the nadir points in this paper are targeted for addressing MaOPs.
The edgetonadir scheme covers the Marcin and Andrzej’s approach [52], Extremized Crowded NSGAII (ECNSGAII) [43], and the recently proposed Emphasized Critical Region (ECR) approach [45]. Specifically, Marcin and Andrzej’s approach decomposed an objective problem into subproblems to estimate the nadir point from the edges, in which the major issues were the poor quality in nadir point found and the impractical computation complexity beyond three objectives [45]. ECNSGAII modified the crowding distance of NSGAII by assigning large rank values to the solutions which had the minimum or maximum objective values. The ECR emphasized the solutions lying at the edges of the PF (i.e., the critical regions) with the adopted MOEAs. Although ECNSGAII and ECR have been reported to be capable of estimating the nadir points in MaOPs, they required a significantly large number of functional evaluations [45].
The extremepointtonadir approaches refer to employ a direct means to estimate the extreme points based on which the nadir points are derived, such as the Worst Crowded NSGAII (WCNSGAII) [43] in which the worst crowded solutions (extreme points) were preferred by ranking their crowding distances with large values. In WCNSGAII, it was hopeful that the extreme points were obtained when the evolution terminated. However, emphasizing the extreme points easily led to the WCNSGAII losing the diversity which inadvertently affected the convergence in turn. In addition, Singh et al. proposed the Pareto Corner Search Evolutionary Algorithm (PCSEA) [44] to look for the nadir points with the cornersort ranking method for the MaOPs whose objective values were required to be with the identical scales.
In addition, there are also various methods not falling into the above categories. For example, Benayoun et al. estimated the nadir points with the paytable [53] in which the th row denoted the objective values of the solution which had the minima on its th objective. In addition, other related works were suggested in [54, 55, 56] for the problems assuming a linear relationship between the objectives and variables. On the contrary, most of the realworld applications are nonlinear in nature.
Because the nadir point estimation is a critical part of the proposed algorithm for solving MaOPs, an approach with a high computational complexity is certainly not preferable. Furthermore, the nadir points are employed for constructing the Utopian PF, while the reference point of IGD would come from the PF. As a consequence, nadir points with a high accuracy are not necessary a guarantee. Considering the balance between the computational complexity and the estimation accuracy, we have proposed in this paper a Decompositionbased Nadir Point Estimation method (DNPE) by transforming an objective MaOP into singleobjective optimization problems to search for the respective extreme points. Thereafter the nadir point is derived.
Specifically, because the proposed nadir point estimation method (i.e., the DNPE) is based on the extremepointtonadir scheme, the WCNSGAII and the PCSEA which are with the similar scheme are detailed further. For convenience of reviewing the related nadir point estimation methods, multiple fundamental concepts of the MaOPs are given first. Then the WCNSGAII and PCSEA are discussed.
Definition 1
Generally, there are extreme points denoted as , , in an objective optimization problem, , and where and .
Definition 2
The nadir point is defined as , where .
Definition 3
The worst point is defined as , where and .
Definition 4
The ideal point is defined as , where and .
Furthermore, the ideal point, extreme point, worst point, nadir point, and the PF are plotted with a biobjective optimization problem in Fig. 1 for intuitively understanding their significance. With these fundamental definitions, a couple of nadir point estimation algorithms, WCNSGAII and PCSEA, which are in the extremepointtopoint scheme are discussed as follows.
WCNSGAII was designed based on NSGAII by modifying its crowding distance assignment. According to the definition of nadir point in Definition 2, WCNSGAII naturally emphasized the solutions with maximal objectives frontwise. Specifically, solutions on a particular nondominated front were sorted with an increasing order based on their fitness, and rank values equal to their positions in the ordered list were assigned. Then the solutions with larger rank values were preferred in each generation during the evolution. By this emphasis mechanism, it was hopeful that nadir point was obtained when the evolution of WCNSGAII was terminated. However, one major deficiency is that overemphasis on these solutions with maximal fitness leads to the lack of diversity, which in turn affects the convergence of the generated extreme points, i.e., the generated extreme points in WCNSGAII are not necessarily Paretooptimal.
PCSEA employed the cornersorting to focus on the extreme points during the evolution. Specifically, there were ascended lists during the executions of PCSEA. The first lists were about the objectives of the solutions, while the other lists were about the excluded square norm with each objective. Furthermore, the th objective of the problem to be optimized was with the excluded square norm . From these lists, solutions with smaller rank values which were equal to their positions in these lists were selected until there was no available slot. Experimental results have shown that PCSEA performs well in MaOPs due to the utilization of cornersorting other than the nondominated sorting which easily leads to the loss of selection pressure. However, the cornersorting can be viewed as to minimize the square norm of all objectives, which deteriorates the performance of PCSEA in solving the problems with different objective value scales and nonconcave PF shapes. For example in Fig. 2 which illustrates an example of biobjective optimization problem, the arc denotes the PF, points and are with different values. It is clearly observed that if the minimization of norm regarding and are emphasized, only the extreme point would be obtained while the other one (point ) would be missed. This deficiency can also be seen in the problems with nonconcave PF shapes.
Briefly, major concerns in these two nadir point estimation algorithms are summarized as 1) overemphasizing extreme points leads to the loss of diversity which in turn deteriorates the convergence of the found nadir points, and 2) simultaneously minimizing the objectives does not scale to problems with different objective value scales and nonconcave PF shapes. To this end, a natural approach is recommended by 1) decomposing the problem to be solved into several singleobjective optimization problems in which the diversity is not required, and 2) assigning different weights to the objectives. In the proposed DNPE, the respective extreme points are estimated by decomposing the objective MaOP into singleobjective problems associated with different weights. Specifically, the th extreme point estimation is with the form formulated by (2)
(2) 
where is a factor with the value greater than to highlight the priority of solving its associated term. In order to better justify our motivation, an example with biobjective is plotted in Fig. 3 in which are the extreme points, are the worse points, and the shaded region denotes the feasible objective space. To obtain the extreme point on the objective, it is required to minimize according to (2). Because is greater than , the term is optimized with a higher priority. Consequently, solutions locating in line are obtained, based on which is minimized much further, then the extreme point is obtained.
2.2 IgdEmoa
Prior to the introduction of IGDEMOA, it is necessary to compare the differences between the IGD and IGD indicators. For this purpose, we first list their respective mathematical formulations. Then the superiority of IGD indicator is highlighted. Finally, the IGDEMOA is discussed much further.
Basically, the IGD indicator is with the form formulated by (3)
(3) 
where denotes a set of reference points in the calculation of IGD, denotes the nondominated solutions generated by the algorithm, denotes the nearest distance from to solutions in , and the distance from to the solution in is calculated by . It has been pointed out in [39] that IGD cannot differentiate the quality of generated solutions when they are nondominated to the solutions in , and the IGD indicator is proposed by changing the calculation of to .
IGDEMOA employed the IGD indicator as its selection mechanism. In addition, the in IGDEMOA is sampled from the approximate PF. Specifically, it supposed that the PF is obtained by solving where is from the nondominated solutions of the current population. However, this approximate approach for generating the PF performs badly in MaOPs where multiple local Paretooptimal exist, which leads IGDEMOA only to solve MaOPs with no more than objectives [40].
In summary, we first introduce details of the proposed DNPE which is motivated by the insightful observations in deficiencies of WCNSGAII and PCSEA. With the help of the estimated nadir point, the Utopian PF is constructed and the reference points are sampled for the calculation of the proposed MaOEA/IGD. Compared to the approximation of the PF in IGDEMOA, our proposed MaOEA/IGD is capable of solving problems with many more objectives. Considering the superiority of IGD, its design principle is employed in the proposed MaOEA/IGD when the generated solutions are nondominated to the sampled reference points.
3 Proposed algorithm
In this section, the proposed Inverted Generational Distance indicatorbased evolutionary algorithm for addressing manyobjective optimization problems (in short for MaOEA/IGD) is presented. To be specific, the framework of the proposed algorithm is outlined first. Then the details of each step in the framework are documented. Next, the computational complexity of the proposed algorithm is analyzed. Finally, the mechanisms of promoting the diversity and the convergence in the proposed algorithm are discussed. Noted here that, the proposed algorithm is described within the context formulated by (1).
3.1 Framework of the Proposed Algorithm
Because the proposed algorithm is based on the IGD indicator, a set of uniformly distributed points which is generated from the PF is required. However, exact points are difficult to be obtained due to the unknown analytical form of the PF for a given realworld application. In the proposed algorithm, a set of reference points, denoted by , which are evenly distributed in the Utopian PF is generated first (Subsection 3.2). Then the population with the predefined size is randomly initialized in the feasible space, and their fitness are evaluated. Next, the population begins to evolve in pursuing the optima until the stopping conditions are satisfied. When the evolution terminates, a number of promising solutions which are hopeful to uniformly distributed in the PF with good proximity is obtained. In order to remedy the shortage caused by using the as the Paretooptimal solutions with promising diversity for IGD indicator, the rank of each solution as well as its proximity distances to all the points in are assigned (Subsection 3.3) first during each generation of the evolution. Then new offspring are generated from their parents selected based on the comparisons over their corresponding rank values and proximity distances (Subsection 3.4). Next, the fitness, ranks, and the proximity distances of the generated offspring to the solutions in are calculated. Finally, a limit number of individuals is selected from the current population to survive into the next generation with the operation of environmental selection (Subsection 3.5). In summary, the details of the framework are listed in Algorithm 3.1.
Uniformly generate reference points for IGD indicator;
Randomly initialize the population;
Fitness evaluation on ;
;
\Whilestopping criteria are not satisfied
Assign the ranks and the proximity distances for individuals in ;
Generate offspring from ;
Fitness evaluation on ;
Assign the rank and the proximity distance for each individual in ;
Environmental selection from ;
;
Return .
3.2 Uniformly Generating Reference Points
In order to obtain the , the extreme points of the problem to be optimized are calculated first. Then the ideal point and the nadir point are extracted. Next, a set of solutions is uniformly sampled from the constrained dimensional hyperplane. Finally, these solutions are transformed into the Utopian PF for the usage of the IGD indicator based on the ideal point and the nadir point. Furthermore, these steps are listed in Algorithm 3.2, and all the details of obtaining the are illustrated as follows.
\KwInOptimization problem ; the size of .
\KwOut.
Estimate the extreme points of with Algorithm 3.2;
Extract the ideal point;
Extract the nadir point;
Uniformly generate points from the constrained hyperplane;
\For to
\For to
=
Return .
To estimate the extreme points, the motivation mentioned in Subsection 2.1 is implemented, and the details are presented in Algorithm 3.2. Specifically, the extreme points are estimated individually based on the objectives of the optimization problem. Furthermore, to estimate the th extreme point , the square norm of is calculated first. Then the absolute value of is calculated. Mathematically, these two steps are formulated by and , respectively, where is the norm operator, and is the absolute value operator. Finally, the extreme point is obtained by line 3.2 in Algorithm 3.2, where is a factor with the value greater than to highlight the weight of the corresponding term in the optimization. When all the extreme points have been estimated, the nadir point and the ideal point are extracted from the extreme points based on Definitions 2 and 4, respectively. This is followed by generating a set of uniformly distributed reference points denoted by from the dimensional constrained hyperplane which is contoured by lines with the unit intercepts in the positive part of the quadrant. Noted here that the Das and Dennis’s method [57], which is widely used by some stateoftheart MaOEAs, such as MOEA/D [14] and NSGAIII [10], is employed for the generation of (line 3.2 of Algorithm 3.2). Ultimately, all the points in are transformed into the Utopian PF, which are detailed in lines 3.23.2 of Algorithm 3.2.
\KwInOptimization problem .
\KwOutExtreme points .
;
\For to
;
\For to
\If
;
;
;
;
Return .
3.3 Assigning Ranks and Proximity Distances
When has been generated, the population is randomly initialized in the feasible search space first, and then the fitness of individuals are evaluated. Next, the rank value and the proximity distances of each individual are assigned. It is noted here that, the rank values are used to distinguish the proximity of the solutions to the Utopian PF from the view of reference points, and the proximity distances are utilized to indicate which individuals are with better convergence and diversity in the subpopulation in which the solutions are with the same rank values. More details are discussed in Subsection 3.7.
Particularly, three rank values, denoted by , , and , exist in the proposed algorithm for all the individuals. Specifically, the way to rank individual is based on the definitions given as follows.
Definition 5
Individual is ranked as , if it dominates at least one solution in .
Definition 6
Individual is ranked as , if it is nondominated to all the solutions in .
Definition 7
Individual is ranked as , if it is dominated by all the solutions in , or dominated by a part of solutions in but nondominated to the remaining solutions.
With Definitions 5, 6, and 7, it is concluded that Paretooptimal solutions are all with rank values , , and , if the PF is convex, a hyperplane, and concave, respectively
Based on the conclusion mentioned above, the proximity distances of each individual with different ranks in the population are calculated. For convenience, it is assumed that there are solutions in , and individuals in the current population. Consequently, there will be proximity distances. Let denote the proximity distance of ( refers to the th individual) to the th point in , where and . Due to one individual with multiple proximity distances, the corresponding minimal proximity distance would be employed when two individuals are compared upon their proximity distances. Because the proximity distance assignment is used to differentiate the convergence of individuals with the same rank values by comparing their associate proximity distances when a prior knowledge of the PF is unknown in advance, the rank of is confirmed first in order to calculate the . Particularly, the proximity distance assignment is designed as follows. If is equal to , is set to the Euclidean distance between and ; if is equal to , is set to the negative value of the Euclidean distance between and ; if is equal to , is calculated by (4).
(4) 
For an intuitive understanding, an example is illustrated in Fig. 4 to present the motivation of the proximity distance assignment for individuals who are ranked as . In Fig. 4, the black solid circles refer to the reference points, the black triangles marked by , , and refer to the individuals , , and which are with the rank (these three individuals are nondominated to the reference points). In this situation, it is clearly shown that the individual is with the smallest minimal proximity distance if the Euclidean distance metric is employed (the minimal proximity distances of individuals , , and are , , and , respectively). Consequently, both the distance measurements of the individuals with ranks and in this situation cannot be utilized to select the desirable individual which has the most promising convergence to the PF. However, if the proximity distance quantified by (4) is employed, it is clearly that the individual is with the smallest minimal proximity distance (the minimal proximity distances of individuals , , and are , , and , respectively), which satisfies the motivation of the proximity distance assignment that a smaller value implies a better convergence.
Noted here that, the proximity distance assignment for the individuals with rank is also employed in the IGD indicator [39]. In summary, smaller proximity distance reveals the better proximity when the exact PF of the problem to be optimized is unknown. Furthermore, the algorithm of the proximity distance assignment is presented in Algorithm 4.
\KwInCurrent population with size ; reference points with size .
\KwOutProximity distances matrix .
\For to
;
Calculate the rank of ;
\For to
\uIf
\uElseIf
\Else
Return .
3.4 Generating Offspring
The process of generating offspring in the proposed algorithm is similar to that in genetic algorithms, in addition to the selection of individuals for filling up the gene pool from which the parent solutions are selected to generate offspring. In this subsection, the processes of generating offspring are elaborated in Steps 15. Then, the details for filling up the gene pool are presented in Algorithm 3.4.

Select solutions from the current population to fill up the gene pool, until it is full.

Select two parent solutions from the gene pool and remove them from the gene pool.

Employ the Simulated binary crossover (SBX) operator to generate offspring with the selected parent solutions.

Employ the polynomial mutation operator to mutate the generated offspring.
[htp]
\KwInCurrent population ; Gene pool size .
\KwOutGene pool .
;
\Whilethe size of is less than
Randomly select two individuals from ;
Obtain the rank of ;
Obtain the rank of ;
Obtain the proximity distances of ;
Obtain the proximity distances of ;
\uIf
;
\uElseIf
;
\Else
\uIf
;
\uElseIf
;
\Else
Randomly select one individual from ;
;
Return .
The binary tournament selection [58] approach is employed in Algorithm 3.4 to select individuals from the current population to fill up the gene pool. In other words, the binary tournament selection [58] approach is employed to select individuals from the current population. Specifically, two individuals which are denoted by and are randomly selected from the current population first (line 3.4). Then, their ranks and the proximity distances are obtained (lines 3.43.4). Next, the individual with smaller rank value is selected to be copied to the gene pool (lines 3.43.4). If and have the same rank values, the individual who has the smaller minimal proximity distance is selected (lines 3.43.4). Otherwise, an individual from and is randomly selected being as one potential parent solution to be put into the gene pool (lines 3.43.4). When the gene pool is full, two parent solutions are randomly selected from the gene pool for generating offspring, and then these selected parent solutions are removed from the gene pool until the gene pool is empty. Noted here that, the SBX [59] and the polynomial mutation [60] operators are employed for the corresponding crossover and mutation operations in the proposed algorithm. It has been reported that two solutions selected in a large search space is not necessary to generate promising offspring [7, 61]. Generally, two ways can be employed to solve this problem. One is the mating restriction method to limit the offspring to be generated by the neighbor solutions [62]. The other one is to use SBX with a large distribution index [10]. In the proposed algorithm, the latter one is utilized due to its simplicity.
3.5 Environmental Selection
When the offspring have been generated, the size of the current population is greater than that of the available slots. As a consequence, the environmental selection takes effects to select a set of representatives to survive to the next generation. In summary, the individuals are selected from the current population according to their assigned rank values and proximity distances. For convenience, it is assumed that there are available slots, the selected individuals are to be stored in , and the individuals with ranks , , as well as are grouped into , , and nondominated fronts, respectively. To be specific, the counter is increased by one until where is a countable operator. If is equal to , the individuals in are copied into and the environmental selection is terminated. Otherwise, the individuals in are copied into first, then individuals are selected from . In summary, the details of the environmental selection are presented in Algorithm 3.5. Furthermore, line 3.5 is confirmed by finding individuals who have the minimal total proximity distances to the reference points (line 3.5), which involves a linear assignment problem (LAP). In the proposed algorithm, the Hungarian method [63] is employed to solve this LAP.
\KwIn, , and ; Available slots size .
\KwOut.
;
;
\While
;
;
\uIf
;
\Else
Uniformly select reference points from ;
Select individuals from ;
;
Return .
3.6 Computational Complexity
In this subsection, the computational complexity of the proposed algorithm is analyzed. For convenience, it is assumed that the problem to be optimized is with objectives, decision variables, desired solutions for decisionmakers, and the computational complexity is analyzed in the context of Algorithm 3.1. To estimate each extreme point, the genetic algorithm is employed, and the SBX as well as polynomial mutation are used as the genetic operators. Furthermore, it is assumed that the population size for estimating extreme points is set to be , and the generation is set to be . Consequently, the total computation cost of uniformly generating reference points for IGD indicator (line 3.1) is . Furthermore, lines 3.1 and 3.1 require and computations, respectively. Because the number of the reference points is equal to that of the desired solutions, the computational complexity of assigning ranks and proximity distances in line 3.1 are and , respectively. Furthermore, generating offspring (line 3.1) needs computations because the size of gene pool is set to be . Since only the fitness, ranks and the proximity distances of the generated offspring need to be calculated, as a consequence, lines 3.1 and 3.1 consume and + , respectively. In the environmental selection, the best case scenario in computational complexity is , while the worst is given that individuals are linearly assignment to the reference points. Furthermore, it is considered common that is greater than , and in MaOPs. Therefore, lines 3.13.1 overall need computations with the generation . In summary, the computational complexity of the proposed algorithm is where is the number of the generation and is the number of solutions.
3.7 Discussions
Loss of selection pressure is a major issue for traditional MOEAs in effectively solving MaOPs because of the traditional domination comparisons between individuals giving a large proportion of nondominated solutions. In the proposed algorithm, the dominance relation of all the individuals are compared to the reference points which are employed for the calculation of IGD indicator. However, the exact reference points which are uniformly distributed in the PF are difficult to obtain. For this purpose, a set of points which are evenly distributed in the Utopian PF are sampled. Furthermore, in order to address this inefficiency given by these approximated reference points, three proximity distances are designed according to their dominance relation to the approximated reference points. This is in hope that less value of the proximity distance means that the corresponding individual is with a better proximity. Specifically, if the solutions with rank are still with the distance calculation of that with ranks or , the convergence will be lost in the proposed algorithm [39].
When the number of solutions to be selected is larger than the available slots, the representatives are chosen from a global view in the proposed algorithm. For convenience of understanding, it is first assumed that representatives need to be selected from solutions where . Then the selection of representatives is simultaneously considered by the calculation of IGD indicator, as oppose to choosing one by one. Simultaneously selecting representatives involves a linear assignment problem. By this linear assignment, each selected reference point can have one distinct individual, which improves the diversity and the convergence simultaneously and this conclusion can also be found in literatures [64, 40, 37]. If the individuals are selected by finding the individual who has the least distance to the reference points, the diversity is not necessarily guaranteed.
4 Experiments
To evaluate the performance of the proposed algorithm in solving MaOPs, a series of experiments is performed. Particularly, NSGAIII [10], MOEA/D [14], HypE [34], RVEA [13], and KnEA [30] are selected as the stateoftheart peer competitors. Although the IGDEMOA can be viewed as the peer algorithm based on IGD indicator, it is merely capable of solving MaOPs with no more than objectives. As a consequence, IGDEMOA is excluded from the list of peer competitors in our experiments.
The remaining of this section is organized as follows. At first, the selected benchmark problems used in this experiment are introduced. Then, the chosen performance metric is given to measure the quality of the approximate Paretooptimal solutions generated by the competing algorithms. Next, the parameter settings employed in all the compared algorithms are listed, and experimental results measured by the considered performance metric are presented and analyzed. Finally, the performance of the proposed algorithm in solving a realworld MaOP is shown (in Section III of the Supplemental Materials), and the performance on the proposed DNPE in estimating nadir point is empirically investigated.
4.1 Benchmark Test Problems
The widely used scalable test problems DTLZ1DTLZ7 from the DTLZ benchmark test suite [65] and WFG1WFG9 from the WFG benchmark test suite [66] are employed in our experiments. Specifically, each objective function in one given objective test problem of DTLZ has decision variables, and is set to be for DTLZ1, for DTLZ2DTLZ6, and for DTLZ7 problems. Moreover, each objective function of a given problem in WFG test suite has decision variables, and is set to be and is set to be based on the suggestion from [66].
4.2 Performance Metric
The widely used Hypervolume (HV) [67] which simultaneously measures the convergence and diversity of the MaOEAs is selected as the performance metric in these experiments. Specifically, the reference points for the calculation of HV are set to be for DTLZ1, for DTLZ2DTLZ6, and for DTLZ7 as well as WFG1WFG9 test problems. Please note that the solutions are discarded for the calculation of HV when they are dominated by the predefined reference points. Because the computational cost increases significantly as the number of objectives grows, Monte Carlo simulation [34]
4.3 Parameter Settings
In this subsection, the baseline parameter settings which are adopted by all the compared MaOEAs are declared first. Then the special parameter settings required by each MaOEA are provided.
Number of Objectives
Test problems with , , and objectives are considered in the experiments because the proposed algorithm aims specifically at effectively solving MaOPs.
Number of Function Evaluations
and Stop Criterion All compared algorithms are individually executed independent times. The maximum number of function evaluations for each compared MaOEA in one independent run is set to be , , and for , , and objective, respectively, which is employed as the termination criterion. Noted that, the parameter settings here are set based on the convention that the maximum generations for the MaOEAs with more than objectives are generally in the order of (the generation number set here is approximately to ). Because of the proposed algorithm includes the phases of nadir point estimation and the optimization for MaOPs, the function evaluations specified here will be shared by these two phases for a fair comparison.
Statistical Approach
Because of the heuristic characteristic of the peer evolutionary algorithms, all the results, which are measured by the performance metric over independent runs for each competing algorithm, are statistically evaluated. In this experiment, the MannWhitneyWilcoxon ranksum test [69] with a significance level is employed for this purpose.
# of divisions  # of population  

8  3,3  240 
15  2,2  240 
20  2,1  230 
Population Size
In principle, the population size can be arbitrarily assigned. However, the population size of the proposed algorithm, NSGAIII, MOEA/D, and RVEA depends on the number of the associated reference points or reference vectors. For a fair comparison, the sizes of population in HypE and KnEA are set to be the same as that in others. In the experiment, the reference points and reference vectors are sampled with the twolayer method [10] and the configurations are listed in Table 1.
MaOEA/IGD  NSGAIII  MOEA/D  HypE  RVEA  KnEA  
Linear 
DTLZ1  8  0.9998(2.93E4)  0.9964(6.12E4)(+)  0.9996(3.52E5)(+)  0.7213(4.31E1)(+)  0.9992(2.15E4)(+)  0.9921(4.21E5)(+) 
15  0.9990(3.18E3)  0.9984(7.23E4)(+)  0.9987(3.20E4)(+)  0.6922(5.45E1)(+)  0.9992(3.99E3)()  0.9994(2.91E4)()  
20  0.9990(2.32E3)  0.9983(3.82E4)(+)  0.9977(7.24E4)(+)  0.7672(3.88E1)(+)  0.9989(3.68E4)(=)  0.9890(3.80E4)(+)  
DTLZ7  8  0.6999(8.97E3)  0.6959(2.72E5)(=)  0.5439(6.37E5)(+)  0.2122(4.82E2)(+)  0.6894(4.24E4)(=)  0.5466(8.12E2)(+)  
15  0.3592(1.20E2)  0.2769(2.33E5)(+)  0.2119(2.42E2)(+)  0.1999(7.68E5)(+)  0.4070(9.50E4)()  0.2804(2.11E2)(+)  
20  0.5261(7.46E4)  0.2348(3.29E4)(+)  0.4325(1.12E4)(+)  0.0986(5.31E4)(+)  0.5206(7.58E2)(+)  0.5166(8.24E5)(+)  
Concave 
DTLZ2  8  0.7174(3.96E3)  0.8132(2.78E3)()  0.5221(3.83E3)(+)  0.1121(3.34E2)(+)  0.6821(3.68E3)(+)  0.7320(3.66E3)(+) 
15  0.9268(2.62E3)  0.8832(9.11E3)(+)  0.3329(1.73E2)(+)  0.0892(4.12E2)(+)  0.9020(3.48E3)(+)  0.8599(6.82E3)(+)  
20  0.8905(6.80E3)  0.9660(3.23E3)()  0.3298(2.10E2)(+)  0.0633(5.32E2)(+)  0.9443(3.19E3)()  0.9307(4.38E3)()  
DTLZ3  8  0.4664(9.25E2)  0.0055(3.80E4)(+)  0.5169(5.68E3)()  0.0085(0.76E5)(+)  0.4572(0.54E3)(=)  0.3537(5.31E4)(+)  
15  0.6984(6.68E2)  0.0091(0.78E5)(+)  0.3030(4.43E3)(+)  0.0133(1.07E5)(+)  0.7183(9.62E2)()  0.5961(0.05E2)(+)  
20  0.7476(7.52E2)  0.0002(6.48E4)(+)  0.2162(4.51E4)(+)  0.0065(5.47E4)(+)  0.6491(2.96E2)(+)  0.7317(7.45E4)(+)  
DTLZ4  8  0.8338(3.31E3)  0.8187(6.22E4)(+)  0.5322(5.87E2)(+)  0.2537(2.08E4)(+)  0.8159(3.01E4)(+)  0.8302(4.71E3)(=)  
15  0.9548(1.66E3)  0.9537(4.24E4)(=)  0.3150(5.08E3)(+)  0.1957(0.86E4)(+)  0.9267(2.62E5)(+)  0.9188(8.01E3)(+)  
20  0.9824(1.33E3)  0.9947(1.37E3)()  0.2755(7.21E5)(+)  0.2101(1.07E4)(+)  0.9854(6.54E2)()  0.9797(4.94E3)(+)  
DTLZ5  8  0.4190(0.64E3)  0.3908(7.67E3)(+)  0.3174(7.15E5)(+)  0.0451(2.24E5)(+)  0.3474(0.88E3)(+)  0.6401(2.65E4)()  
15  0.2677(9.71E3)  0.2178(5.34E5)(+)  0.1821(8.85E2)(+)  0.0418(8.99E5)(+)  0.1379(7.84E4)(+)  0.1257(3.46E3)(+)  
20  0.2101(5.57E3)  0.3390(0.44E2)()  0.1790(0.99E4)(+)  0.0423(4.90E2)(+)  0.3606(5.52E2)()  0.4139(5.49E2)()  
DTLZ6  8  0.7202(8.94E4)  0.2866(0.22E5)(+)  0.3037(7.95E3)(+)  0.0548(7.36E4)(+)  0.2467(9.72E2)(+)  0.4547(4.17E3)(+)  
15  0.7756(8.36E5)  0.4385(7.14E3)(+)  0.6748(4.56E2)(+)  0.1957(0.86E4)(+)  0.4918(9.97E2)(+)  0.7314(1.56E2)(+)  
20  0.8639(0.20E3)  0.9081(8.32E2)()  0.5170(5.20E4)(+)  0.1080(6.17E2)(+)  0.4438(2.39E4)(+)  0.3002(5.79E5)(+)  
+/=/  14/2/5  20/0/1  21/0/0  12/3/6  16/1/4 
Genetic Operators
The SBX [59] and polynomial mutation [70] are employed as the genetic operators. Moreover, the probabilities of the crossover and mutation are set to be and , respectively. The distribution indexes of mutation and crossover are set to be , in addition to NSGAIII whose mutation distribution index is specifically set to be based on the recommendation from its developers [10].
In solving the proposed nadir point estimation method for constructing the Utopian PF, evolutionary algorithm is employed. To be specific, the SBX and polynomial mutation, both of whose distribution index are set to be , and probabilities for crossover and mutation are set to be and , respectively, are utilized as the genetic operators. In addition, the population sizes are set to be the same to those in Table 1, and the numbers of generations for all are set to be . Besides, the balance parameter in (2) is specified as .
4.4 Experimental Results and Analysis
In this subsection, the results, which are generated by competing algorithms over considered test problems with specific objective numbers and then measured by the selected performance metric, are presented and analyzed to highlight the superiority of the proposed algorithm in addressing MaOPs. Specifically, the mean values as well as the standard deviations of HV results over DTLZ1DTLZ7 and WFG1WFG9 test problems are listed in Tables 2 and 3, respectively. Furthermore, the numbers with bold face imply the best mean values over the corresponding test problem with a given objective number (the second and third columns in Tables 2 and 3) against all compared algorithms. Moreover, the symbols “+,” “,” and “=” indicate whether the null hypothesis of the results, which are generated by the proposed algorithm and corresponding compared peer competitor, is accepted or rejected with the significance level by the considered ranksum test. In addition, the last rows in Tables 2 and 3 present the summarizations indicating how many times the proposed algorithm performs better than, worse than or equal to the chosen peer competitor, respectively. In order to conveniently investigate the experimental results of the welldesigned proximity distance assignments in the proposed MaOEA/IGD and the conclusion in Section 3.3, test problems are group into “Convex,” “Linear,” and “Concave” based on the respective test problem features, and displayed in the first columns of Tables 2 and 3. Noted that, although the PFs of DTLZ7 and WFG1 are mixed, they are classified into the “Linear” category due to their PF shapes being more similar to linear.
From the results measured by HV on DTLZ1DTLZ7 test problems (Table 2), it is clearly shown that MaOEA/IGD achieves the best performance among its peer competitors upon  and objective DTLZ1 and DTLZ7, while performs slightly worse upon objective DTLZ1 by KnEA and DTLZ7 by RVEA. Furthermore, MaOEA/IGD also wins the best scores on  and objective DTLZ4 and DTLZ6, but is defeated by NSGAIII upon these two problems with objective. Although NSGAIII and KnEA show better performance upon  and objective DTLZ2 and DTLZ5, MaOEA/IGD is the winner upon objective DTLZ2 and DTLZ5. In addition, MaOEA/IGD achieves the best score upon objective DTLZ3.
MaOEA/IGD  NSGAIII  MOEA/D  HypE  RVEA  KnEA  
Convex 
WFG2  8  0.9839(2.00E2)  0.9587(9.10E3)(+)  0.9474(9.09E2)(+)  0.9514(5.92E4)(+)  0.9380(3.33E3)(+)  0.9686(8.53E2)(+) 
15  0.9362(7.76E3)  0.9672(3.16E2)()  0.9402(7.00E4)()  0.6216(6.25E3)(+)  0.9475(9.71E4)()  0.9360(8.37E3)(=)  
20  0.9782(2.52E2)  0.9624(0.11E2)(+)  0.9460(5.73E2)(+)  0.8354(7.90E4)(+)  0.9657(4.04E2)(+)  0.8106(5.11E2)(+)  
Linear 
WFG1  8  0.9578(1.30E1)  0.9255(9.48E2)(+)  0.9454(0.61E4)(+)  0.6528(5.85E4)(+)  0.8383(2.85E4)(+)  0.5847(8.28E2)(+) 
15  0.9354(5.02E3)  0.9536(3.29E3)()  0.9405(6.50E2)(=)  0.6395(9.75E3)(+)  0.9538(1.84E3)()  0.9373(4.98E2)(+)  
20  0.9806(2.64E2)  0.9405(8.01E2)(+)  0.9593(1.43E2)(+)  0.6340(4.78E2)(+)  0.9020(5.43E2)(+)  0.9374(8.84E2)(+)  
WFG3  8  0.8615(1.08E1)  0.8423(7.55E2)(+)  0.8521(7.42E4)(=)  0.5660(8.31E3)(+)  0.8457(2.34E2)(+)  0.8618(1.57E2)(=)  
15  0.3346(4.58E3)  0.5091(4.10E2)()  0.3393(1.32E2)(=)  0.2852(5.41E3)(+)  0.5188(2.43E2)()  0.5076(8.26E4)()  
20  0.5750(3.13E2)  0.5313(3.89E2)(+)  0.4863(4.29E2)(+)  0.2918(9.56E3)(+)  0.3425(5.73E2)(+)  0.4579(8.50E2)(+)  
Concave 
WFG4  8  0.7800(2.64E2)  0.7877(7.02E2)(=)  0.7507(3.75E4)(+)  0.7229(9.74E3)(+)  0.7648(7.29E2)(+)  0.7715(1.74E4)(+) 
15  0.8314(5.76E3)  0.6315(0.01E2)(+)  0.8414(0.03E2)()  0.4801(0.87E2)(+)  0.8154(2.61E3)(+)  0.8690(0.23E2)()  
20  0.8108(3.86E2)  0.7885(4.68E3)(+)  0.7707(8.61E2)(+)  0.7309(4.67E2)(+)  0.8094(8.35E3)(=)  0.4794(7.43E2)(+)  
WFG5  8  0.8653(1.10E1)  0.7993(4.57E4)(+)  0.4429(6.68E2(+)  0.4915(6.99E3)(+)  0.8638(6.51E3)(+)  0.8741(3.09E2)()  
15  0.8335(6.31E3)  0.6408(1.69E2)(+)  0.3397(0.01E2)(+)  0.4351(4.18E2)(+)  0.7553(4.88E2)(+)  0.8276(1.60E3)(+)  
20  0.8905(6.28E3)  0.8022(9.87E4)(+)  0.4964(0.84E2)(+)  0.3125(2.50E2)(+)  0.7946(9.13E4)(+)  0.7258(6.64E3)(+)  
WFG6  8  0.9785(3.12E2)  0.9918(8.77E2)()  0.9488(8.06E4(+)  0.9261(4.61E2)(+)  0.9976(6.97E4)()  0.9876(3.68E2)()  
15  0.9357(8.22E3)  0.8756(5.13E3)(+)  0.8327(2.41E4)(+)  0.8406(2.60E3)(+)  0.8538(0.21E2)(+)  0.9223(8.21E2)(+)  
20  0.8854(1.82E2)  0.8189(2.62E2)(+)  0.8489(5.80E2)(+)  0.7633(8.78E4)(+)  0.7968(5.83E2)(+)  0.7502(5.00E3)(+)  
WFG7  8  0.8858(5.43E3)  0.8118(7.25E2)(+)  0.7430(8.58E4)(+)  0.7416(3.48E2)(+)  0.8192(2.51E2)(+)  0.7635(5.82E2)(+)  
15  0.8352(6.48E3)  0.8780(3.82E2)()  0.7343(7.92E4)(+)  0.4030(8.39E3)(+)  0.6366(1.79E4)(+)  0.5463(1.70E3)(+)  
20  0.7919(3.89E3)  0.8482(1.35E4)()  0.4844(9.14E3)(+)  0.6418(6.41E3)(+)  0.8706(5.71E1)()  0.8116(9.03E3)()  
WFG8  8  0.6839(1.78E2)  0.6869(1.39E2)()  0.4405(3.49E3)(+)  0.3155(1.51E3)(+)  0.5908(5.04E3)(+)  0.6850(5.72E3)()  
15  0.7340(7.08E3)  0.5470(5.14E3)(+)  0.3412(8.14E4)(+)  0.2065(0.97E2)(+)  0.6455(5.90E4)(+)  0.6246(1.24E2)(+)  
20  0.7831(2.08E2)  0.6855(5.75E3)(+)  0.4928(9.16E2)(+)  0.1117(4.95E4)(+)  0.7844(8.87E3)(=)  0.6027(4.21E2)(+)  
WFG9  8  0.7694(8.69E0)  0.7328(8.73E3)(+)  0.4488(0.55E3)(+)  0.3030(5.00E2)(+)  0.7444(3.41E2)(+)  0.7528(4.91E2)(+)  
15  0.8329(8.16E3)  0.6105(0.13E2)(+)  0.3359(7.18E4)(+)  0.2176(3.91E2)(+)  0.7294(0.34E3)(+)  0.6595(4.06E2)(+)  
20  0.7829(2.79E2)  0.7299(5.76E4)(+)  0.4881(8.07E4)(+)  0.1923(6.55E1)(+)  0.7128(8.78E3)(+)  0.7824(5.36E2)(=)  
+/=/  19/1/7  22/3/2  27/0/0  20/2/5  18/3/6 
The HV results from WFG1WFG9 test problems generated by competing algorithms are listed in Table 3. For objective WFG test problems, MaOEA/IGD shows a better performance on WFG1, WFG2, WFG7, and WFG9 than its peer competitors, and performs a little worse than that of KnEA on WFG5, RVEA on WFG6, and NSGAIII on WFG8 test problems. Although MaOEA/IGD does not show the best scores on WFG3 and WFG4, it obtains similar statistical results compared to the respective winners (i.e., KnEA and NSGAIII). For objective test problems, MaOEA/IGD shows a better performance on WFG5, WFG6, WFG8, and WFG9 than competing algorithms, while worse than RVEA on WFG2 and WFG3, NSGAIII on WFG2, and KnEA on WFG4. Although NSGAIII performs better than MaOEA/IGD on WFG7, MaOEA/IGD performs better than all other peer competitors. In addition, MaOEA/IGD wins over NSGAIII, MOEA/D, HypE, RVEA, and KnEA on objective WFG1, WFG2, WFG3, WFG4, WFG5, WFG6, and WFG9, but underperforms on WFG7 and WFG8 in which RVEA performs better.
Briefly, MaOEA/IGD wins 9 times out of the 12 comparisons upon the test problems whose PF shapes are linear (i.e., DTLZ1, DTLZ7, WFG1, and WFG3), which can be interpreted that the sampled reference points from the Utopian PF for the proposed algorithm are the Paretooptimal solutions due to the linear feature of the PF, and the proximity distance assignment for the solutions with rank value has taken effects. Furthermore, MaOEA/IGD shows competitive performance on WFG2 test problem whose feature of the PF is convex. Because the sampled reference points on the Utopian PF are all nondominated by the Paretooptimal solutions, the proximity distances for solutions with rank in MaOEA/IGD take effects in this situation. In addition, it is no strange that MaOEA/IGD obtains better results on most of other test problems whose PF features are concave because the reference points utilized to maintain the diversity and convergence of the proposed algorithm dominate the solutions uniformly generated from the PF. In summary, the proposed algorithm shows considerable competitiveness against considered competing algorithms in addressing selected MaOPs with the results measured by the HV performance metric.
Theoretically, the major shortcoming of HV indicator against IGD is its much higher computational complexity. However, noted that from Tables 2 and 3, the proposed algorithm, which is designed based on the IGD indicator, outperforms HypE, which is motivated by the HV indicator, upon all test problems with the selected numbers of objectives, although the numbers of function evaluations regarding HypE is set to be a much large number. The deficiencies of HypE in this regard are explained as follows. First, it has been reported in [71, 72, 15, 11] that the HV result is largely affected by the nadir points of the problem to be optimized. In HypE, the nadir points are determined as the evolution continues. In this way, the obtained nadir point would be inaccurate during the early evolution process (the reasons have been discussed in reviewing the nadir point estimation approaches in Section 2), which leads to the worse performance of HypE. Secondly, the HV results of HypE in solving MaOPs are estimated by Monte Carlo simulation, while the number of reference points in Monte Carlo simulation is critical to the successful performance [34]. In practice, that number is unknown and unavailable of such may lead to a poor performance.
4.5 Investigation on Nadir Point Estimation
In this subsection, we will investigate the performance of the proposed DNPE on estimating the nadir point. To be specific, two peer competitors including WCNSGAII and PCSEA which have been discussed in Section 2 are utilized to perform comparisons on selected test problems. In these comparisons, the numbers of function evaluations regarding each compared algorithm are counted until 1) the metric formulated by (5)
(5) 
where denotes the th element of the estimated nadir point derived from the extreme points generated by the compared algorithm or 2) the maximum function evaluation numbers is met. The experimental results for DTLZ1, DTLZ2, and WFG2 with , , , and objective are plotted in Fig. 5. Please note that the reason of choosing these three test problems is that they cover the various shapes of PF (i.e., DTLZ1, DTLZ2, and WFG2 are with linear, concave, and convex PF, respectively) and characteristics of objective value scales (i.e., DTLZ1 and DTLZ2 are with the same objective value scales while WFG2 is not). Specifically, the ideal points of DTLZ1, DTLZ2, and WFG2 are , and the nadir points are , , and , respectively. In addition, the population size is specified as , the probabilities of SBX and polynomial mutation are set to be and , and both of distribution index are set to be . Because the proposed DNPE is based on the decomposition to estimate the nadir point, and maximum function evaluation number with are set to be the stopping criteria for estimating each extreme point.
The results performed by compared nadir point estimation methods on , , , and objective DTLZ1, DTLZ2, and WFG2 are illustrated in Figs. a, b, and c, respectively. It is clearly shown in Fig. a that these compared algorithms find the satisfactory nadir points of the DTLZ1 which is with the linear PF within the predefined maximum function evaluation numbers, and the proposed DNPE takes the least numbers of function evaluations over the four considered objective numbers. Moreover, WCNSGAII cannot find the nadir point over DTLZ2 with concave PF and WFG2 with convex PF with , , and objective, and PCSEA cannot find the nadir point over WFG2 with different objective value scales, while the proposed DNPE performs well on both test problems with all considered objective numbers. In addition, the proposed DNPE is scalable to the objective number in the estimating nadir points of the MaOPs, which can be seen from Figs. a and b. In summary, the proposed DNPE shows quality performance in estimating nadir point of MaOPs with different PF features and objective scales.
5 Conclusion and Future Works
In this paper, an IGD indicatorbased evolutionary algorithm is proposed for solving manyobjective optimization problems. In order to obtain a set of uniformly distributed reference points for the calculation of the IGD indicator, a decompositionbased nadir point estimation method is designed to construct the Utopian PF in which the reference points can be easily sampled. For solving the deficiency of the Utopian PF being as the PF in the phase of sampling the reference points, one rank assignment mechanism is proposed to compare the dominance relation of the solutions to the reference points, based on which three types of proximity distance assignments are designed to distinct the quality of the solutions with the same front rank values. In addition, the linear assignment principle is utilized as the selection mechanism to choose representatives for concurrently facilitating the convergence and diversity of the proposed algorithm. In summary, based on the proposed nadir estimation method, the proposed dominance comparison approach, rank value and proximity distance assignment, and selection mechanism collectively improve the evolution of the proposed algorithm towards the PF with promising diversity. In order to qualify the performance of the proposed algorithm, a series of welldesigned experiments is performed over two widely used benchmark test suites with , , and objective, their results measured by the selected performance metric indicate that the proposed algorithm is with considerable competitiveness in solving manyobjective optimization problems. In addition, we utilize the proposed algorithm to solve one realworld manyobjective optimization problem, in which the satisfactory results demonstrate the superiority of the proposed algorithm. Moreover, experiments are performed by the proposed decompositionbased nadir point estimation method against a couple of competitors over three representative test problems (DTLZ1, DTLZ2, and WFG2) with challenging features in PF shapes and objective value scales, the experimental results reveal the satisfactory results obtained by the proposed nadir point estimation method. In near future, we will place our efforts mainly on two essential aspects 1) constructing more accurate PF with limited information priori to obtaining the Paretooptimal solutions to improve the development of indicatorbased algorithms which require the uniformly distributed reference points, and 2) extending the proposed algorithm to solve constrained manyobjective optimization problems.
[]Yanan Sun (S’15M’18) received a Ph.D. degree in engineering from the Sichuan University, Chengdu, China, in 2017. From 2015.082017.02, he is a jointly Ph.D. student financed by the China Scholarship Council in the School of Electrical and Computer Engineering, Oklahoma State University (OSU), USA. He is currently a Postdoc Research Fellow in the School of Engineering and Computer Science, Victoria University of Wellington, Wellington, New Zealand. His research topics are manyobjective optimization and deep learning.
[]Gary G. Yen (S’87M’88SM’97F’09) received a Ph.D. degree in electrical and computer engineering from the University of Notre Dame in 1992. Currently he is a Regents Professor in the School of Electrical and Computer Engineering, Oklahoma State University (OSU). Before joined OSU in 1997, he was with the Structure Control Division, U.S. Air Force Research Laboratory in Albuquerque. His research interest includes intelligent control, computational intelligence, conditional health monitoring, signal processing and their industrial/defense applications.
Dr. Yen was an associate editor of the IEEE Control Systems Magazine, IEEE Transactions on Control Systems Technology, Automatica, Mechantronics, IEEE Transactions on Systems, Man and Cybernetics, Parts A and B and IEEE Transactions on Neural Networks. He is currently serving as an associate editor for the IEEE Transactions on Evolutionary Computation and the IEEE Transactions on Cybernetics. He served as the General Chair for the 2003 IEEE International Symposium on Intelligent Control held in Houston, TX and 2006 IEEE World Congress on Computational Intelligence held in Vancouver, Canada. Dr. Yen served as Vice President for the Technical Activities in 20052006 and then President in 20102011 of the IEEE Computational intelligence Society. He was the founding editorinchief of the IEEE Computational Intelligence Magazine, 20062009. In 2011, he received Andrew P Sage Best Transactions Paper award from IEEE Systems, Man and Cybernetics Society and in 2014, he received Meritorious Service award from IEEE Computational Intelligence Society.
[]Zhang Yi (F’16) received a Ph.D. degree in mathematics from the Institute of Mathematics, The Chinese Academy of Science, Beijing, China, in 1994. Currently, he is a Professor at the Machine Intelligence Laboratory, College of Computer Science, Sichuan University, Chengdu, China. He is the coauthor of three books: Convergence Analysis of Recurrent Neural Networks (Kluwer Academic Publishers, 2004), Neural Networks: Computational Models and Applications (Springer, 2007), and Subspace Learning of Neural Networks (CRC Press, 2010). He was an Associate Editor of IEEE Transactions on Neural Networks and Learning Systems (2009 2012), and He is an Associate Editor of IEEE Transactions on Cybernetics (2014 ). His current research interests include Neural Networks and Big Data. He is a fellow of IEEE.
Footnotes
 This is considered in the context of a minimization problem with continuous PF, and the extreme points are excluded from the Paretooptimal solutions.
 The source code is available at: http://www.tik.ee.ethz.ch/sop/download/supplementary/hype/.
 The source code is available at: http://www.wfg.csse.uwa.edu.au/hypervolume/.
References
 V. Khare, X. Yao, and K. Deb, “Performance scaling of multiobjective evolutionary algorithms,” in International Conference on Evolutionary MultiCriterion Optimization. Springer, 2003, pp. 376–390.
 O. Chikumbo, E. Goodman, and K. Deb, “Approximating a multidimensional pareto front for a land use management problem: A modified moea with an epigenetic silencing metaphor,” in Evolutionary Computation (CEC), 2012 IEEE Congress on. IEEE, 2012, pp. 1–9.
 R. J. Lygoe, M. Cary, and P. J. Fleming, “A realworld application of a manyobjective optimisation complexity reduction process,” in Evolutionary MultiCriterion Optimization. Springer, 2013, pp. 641–655.
 K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan, “A fast and elitist multiobjective genetic algorithm: Nsgaii,” IEEE Transactions on Evolutionary Computation, vol. 6, no. 2, pp. 182–197, 2002.
 E. Zitzler, M. Laumanns, L. Thiele, E. Zitzler, E. Zitzler, L. Thiele, and L. Thiele, “Spea2: Improving the strength pareto evolutionary algorithm,” 2001.
 C. M. Fonseca and P. J. Fleming, “Multiobjective optimization and multiple constraint handling with evolutionary algorithms. i. a unified formulation,” Systems, Man and Cybernetics, Part A: Systems and Humans, IEEE Transactions on, vol. 28, no. 1, pp. 26–37, 1998.
 R. C. Purshouse and P. J. Fleming, “On the evolutionary optimization of many conflicting objectives,” IEEE Transactions on Evolutionary Computation, vol. 11, no. 6, pp. 770–784, 2007.
 B. Li, J. Li, K. Tang, and X. Yao, “Manyobjective evolutionary algorithms: A survey,” ACM Computing Surveys (CSUR), vol. 48, no. 1, p. 13, 2015.
 T. Wagner, N. Beume, and B. Naujoks, “Pareto, aggregation, and indicatorbased methods in manyobjective optimization,” in Evolutionary multicriterion optimization. Springer, 2007, pp. 742–756.
 K. Deb and H. Jain, “An evolutionary manyobjective optimization algorithm using referencepointbased nondominated sorting approach, part i: Solving problems with box constraints,” IEEE Transactions on Evolutionary Computation, vol. 18, no. 4, pp. 577–601, 2014.
 Y. Yuan, H. Xu, B. Wang, and X. Yao, “A new dominance relationbased evolutionary algorithm for manyobjective optimization,” IEEE Transactions on Evolutionary Computation, vol. 20, no. 1, pp. 16–37, 2016.
 Y. Sun, G. G. Yen, and Z. Yi, “Reference linebased estimation of distribution algorithm for manyobjective optimization,” KnowledgeBased Systems, p. doi:10.1016/j.knosys.2017.06.021, 2017.
 R. Cheng, Y. Jin, M. Olhofer, and B. Sendhoff, “A reference vector guided evolutionary algorithm for manyobjective optimization,” IEEE Transactions on Evolutionary Computation, vol. 20, no. 5, pp. 773–791, 2016.
 Q. Zhang and H. Li, “Moea/d: A multiobjective evolutionary algorithm based on decomposition,” IEEE Transactions on Evolutionary Computation, vol. 11, no. 6, pp. 712–731, 2007.
 Y. Yuan, H. Xu, B. Wang, B. Zhang, and X. Yao, “Balancing convergence and diversity in decompositionbased manyobjective optimizers,” IEEE Transactions on Evolutionary Computation, vol. 20, no. 2, pp. 180–198, 2016.
 Z. Wang, Q. Zhang, M. Gong, and A. Zhou, “A replacement strategy for balancing convergence and diversity in moea/d,” in 2014 IEEE Congress on Evolutionary Computation (CEC). IEEE, 2014, pp. 2132–2139.
 K. Li, Q. Zhang, S. Kwong, M. Li, and R. Wang, “Stable matchingbased selection in evolutionary multiobjective optimization,” IEEE Transactions on Evolutionary Computation, vol. 18, no. 6, pp. 909–923, 2014.
 K. Li, S. Kwong, Q. Zhang, and K. Deb, “Interrelationshipbased selection for decomposition multiobjective optimization,” IEEE Transactions on Cybernetics, vol. 45, no. 10, pp. 2076–2088, 2015.
 M. Asafuddoula, T. Ray, and R. Sarker, “A decompositionbased evolutionary algorithm for many objective optimization,” IEEE Transactions on Evolutionary Computation, vol. 19, no. 3, pp. 445–460, 2015.
 S. B. Gee, K. C. Tan, V. A. Shim, and N. R. Pal, “Online diversity assessment in evolutionary multiobjective optimization: A geometrical perspective,” IEEE Transactions on Evolutionary Computation, vol. 19, no. 4, pp. 542–559, 2015.
 M. Laumanns, L. Thiele, K. Deb, and E. Zitzler, “Combining convergence and diversity in evolutionary multiobjective optimization,” Evolutionary computation, vol. 10, no. 3, pp. 263–282, 2002.
 F. di Pierro, S.T. Khu, and D. A. Savic, “An investigation on preference order ranking scheme for multiobjective evolutionary optimization,” IEEE Transactions on Evolutionary Computation, vol. 11, no. 1, pp. 17–45, 2007.
 G. Wang and H. Jiang, “Fuzzydominance and its application in evolutionary many objective optimization,” in Computational Intelligence and Security Workshops, 2007. CISW 2007. International Conference on. IEEE, 2007, pp. 195–198.
 Z. He, G. G. Yen, and J. Zhang, “Fuzzybased pareto optimality for manyobjective evolutionary algorithms,” IEEE Transactions on Evolutionary Computation, vol. 18, no. 2, pp. 269–285, 2014.
 X. Zou, Y. Chen, M. Liu, and L. Kang, “A new evolutionary algorithm for solving manyobjective optimization problems,” Systems, Man, and Cybernetics, Part B: Cybernetics, IEEE Transactions on, vol. 38, no. 5, pp. 1402–1412, 2008.
 S. Yang, M. Li, X. Liu, and J. Zheng, “A gridbased evolutionary algorithm for manyobjective optimization,” IEEE Transactions on Evolutionary Computation, vol. 17, no. 5, pp. 721–736, 2013.
 A. Lâopez, C. A. C. Coello, A. Oyama, and K. Fujii, “An alternative preference relation to deal with manyobjective optimization problems,” in International Conference on Evolutionary MultiCriterion Optimization. Springer, 2013, pp. 291–306.
 M. Li, S. Yang, and X. Liu, “Shiftbased density estimation for paretobased algorithms in manyobjective optimization,” IEEE Transactions on Evolutionary Computation, vol. 18, no. 3, pp. 348–365, 2014.
 J. Cheng, G. G. Yen, and G. Zhang, “A manyobjective evolutionary algorithm with enhanced mating and environmental selections,” IEEE Transactions on Evolutionary Computation, vol. 19, no. 4, pp. 592–605, 2015.
 X. Zhang, Y. Tian, and Y. Jin, “A knee pointdriven evolutionary algorithm for manyobjective optimization,” IEEE Transactions on Evolutionary Computation, vol. 19, no. 6, pp. 761–776, 2015.
 M. Emmerich, N. Beume, and B. Naujoks, “An emo algorithm using the hypervolume measure as selection criterion,” in International Conference on Evolutionary MultiCriterion Optimization. Springer, 2005, pp. 62–76.
 C. Igel, N. Hansen, and S. Roth, “Covariance matrix adaptation for multiobjective optimization,” Evolutionary computation, vol. 15, no. 1, pp. 1–28, 2007.
 D. Brockhoff and E. Zitzler, “Improving hypervolumebased multiobjective evolutionary algorithms by using objective reduction methods,” in 2007 IEEE Congress on Evolutionary Computation. IEEE, 2007, pp. 2086–2093.
 J. Bader and E. Zitzler, “Hype: An algorithm for fast hypervolumebased manyobjective optimization,” Evolutionary computation, vol. 19, no. 1, pp. 45–76, 2011.
 K. Gerstl, G. Rudolph, O. Schütze, and H. Trautmann, “Finding evenly spaced fronts for multiobjective control via averaging hausdorffmeasure,” in Electrical Engineering Computing Science and Automatic Control (CCE), 2011 8th International Conference on. IEEE, 2011, pp. 1–6.
 H. Trautmann, G. Rudolph, C. DominguezMedina, and O. Schütze, “Finding evenly spaced pareto fronts for threeobjective optimization problems,” in EVOLVEA Bridge between Probability, Set Oriented Numerics, and Evolutionary Computation II. Springer, 2013, pp. 89–105.
 C. A. Rodríguez Villalobos and C. A. Coello Coello, “A new multiobjective evolutionary algorithm based on a performance assessment indicator,” in Proceedings of the 14th annual conference on Genetic and evolutionary computation. ACM, 2012, pp. 505–512.
 R. Storn and K. Price, Differential evolutiona simple and efficient adaptive scheme for global optimization over continuous spaces. ICSI Berkeley, 1995, vol. 3.
 H. Ishibuchi, H. Masuda, Y. Tanigaki, and Y. Nojima, “Modified distance calculation in generational distance and inverted generational distance,” in International Conference on Evolutionary MultiCriterion Optimization. Springer, 2015, pp. 110–125.
 E. M. Lopez and C. A. C. Coello, “Igd+emoa: A multiobjective evolutionary algorithm based on igd+,” in Evolutionary Computation (CEC), 2016 IEEE Congress on. IEEE, 2016, pp. 999–1006.
 Y. Sun, G. G. Yen, H. Mao, and Z. Yi, “Manifold dimension reduction based clustering for multiobjective evolutionary algorithm,” in Evolutionary Computation (CEC), 2016 IEEE Congress on. IEEE, 2016, pp. 3785–3792.
 Y. Sun, G. G. Yen, and Z. Yi, “Global viewbased selection mechanism for manyobjective evolutionary algorithms,” in Evolutionary Computation (CEC), 2017 IEEE Congress on. IEEE, 2017, pp. 427–434.
 K. Deb, S. Chaudhuri, and K. Miettinen, “Towards estimating nadir objective vector using evolutionary approaches,” in Proceedings of the 8th annual conference on Genetic and evolutionary computation. ACM, 2006, pp. 643–650.
 H. K. Singh, A. Isaacs, and T. Ray, “A pareto corner search evolutionary algorithm and dimensionality reduction in manyobjective optimization problems,” IEEE Transactions on Evolutionary Computation, vol. 15, no. 4, pp. 539–556, 2011.
 H. Wang, S. He, and X. Yao, “Nadir point estimation for manyobjective optimization problems based on emphasized critical regions,” Soft Computing, pp. 1–13, 2015.
 K. Deb and K. Miettinen, “A review of nadir point estimation procedures using evolutionary approaches: A tale of dimensionality reduction,” Indian Instit. Technol., Kanpur, India, KanGAL Rep, vol. 2, no. 008, p. 004, 2008.
 L. Ke, Q. Zhang, and R. Battiti, “Moea/daco: a multiobjective evolutionary algorithm using decomposition and antcolony,” IEEE Transactions on Cybernetics, vol. 43, no. 6, pp. 1845–1859, 2013.
 X. Ma, F. Liu, Y. Qi, X. Wang, L. Li, L. Jiao, M. Yin, and M. Gong, “A multiobjective evolutionary algorithm based on decision variable analyses for multiobjective optimization problems with largescale variables,” IEEE Transactions on Evolutionary Computation, vol. 20, no. 2, pp. 275–298, 2016.
 N. Chen, W.n. Chen, Y.J. Gong, Z.H. Zhan, J. Zhang, Y. Li, and Y.S. Tan, “An evolutionary algorithm with doublelevel archives for multiobjective optimization,” IEEE Transactions on Cybernetics, vol. 45, no. 9, pp. 1851–1863, 2015.
 Y. Yuan, H. Xu, B. Wang, B. Zhang, and X. Yao, “Balancing convergence and diversity in decompositionbased manyobjective optimizers,” IEEE Transactions on Evolutionary Computation, vol. 20, no. 2, pp. 180–198, 2015.
 K. Praditwong and X. Yao, “How well do multiobjective evolutionary algorithms scale to large problems,” in 2007 IEEE Congress on Evolutionary Computation. IEEE, 2007, pp. 3959–3966.
 M. Szczepański and A. P. Wierzbicki, “Application of multiple criteria evolutionary algorithms to vector optimisation, decision support and reference point approaches,” Journal of Telecommunications and Information Technology, pp. 16–33, 2003.
 R. Benayoun, J. De Montgolfier, J. Tergny, and O. Laritchev, “Linear programming with multiple objective functions: Step method (stem),” Mathematical programming, vol. 1, no. 1, pp. 366–375, 1971.
 M. Dessouky, M. Ghiassi, and W. Davis, “Estimates of the minimum nondominated criterion values in multiplecriteria decisionmaking,” Engineering Costs and Production Economics, vol. 10, no. 2, pp. 95–104, 1986.
 H. Isermann and R. E. Steuer, “Computational experience concerning payoff tables and minimum criterion values over the efficient set,” European journal of operational research, vol. 33, no. 1, pp. 91–97, 1988.
 P. Korhonen, S. Salo, and R. E. Steuer, “A heuristic for estimating nadir criterion values in multiple objective linear programming,” Operations Research, vol. 45, no. 5, pp. 751–757, 1997.
 I. Das and J. E. Dennis, “Normalboundary intersection: A new method for generating the pareto surface in nonlinear multicriteria optimization problems,” SIAM Journal on Optimization, vol. 8, no. 3, pp. 631–657, 1998.
 B. L. Miller and D. E. Goldberg, “Genetic algorithms, tournament selection, and the effects of noise,” Complex systems, vol. 9, no. 3, pp. 193–212, 1995.
 K. Deb and R. B. Agrawal, “Simulated binary crossover for continuous search space,” Complex Systems, vol. 9, no. 3, pp. 1–15, 1994.
 K. Deb, Multiobjective optimization using evolutionary algorithms. John Wiley & Sons, 2001, vol. 16.
 S. F. Adra and P. J. Fleming, “Diversity management in evolutionary manyobjective optimization,” IEEE Transactions on Evolutionary Computation, vol. 15, no. 2, pp. 183–195, 2011.
 K. Deb and D. E. Goldberg, “An investigation of niche and species formation in genetic function optimization,” in Proceedings of the 3rd international conference on genetic algorithms. Morgan Kaufmann Publishers Inc., 1989, pp. 42–50.
 R. Jonker and T. Volgenant, “Improving the hungarian assignment algorithm,” Operations Research Letters, vol. 5, no. 4, pp. 171–175, 1986.
 J. A. M. Berenguer and C. A. C. Coello, “Evolutionary manyobjective optimization based on kuhnmunkresâ algorithm,” in International Conference on Evolutionary MultiCriterion Optimization. Springer, 2015, pp. 3–17.
 K. Deb, L. Thiele, M. Laumanns, and E. Zitzler, Scalable test problems for evolutionary multiobjective optimization. Springer, 2005.
 S. Huband, P. Hingston, L. Barone, and L. While, “A review of multiobjective test problems and a scalable test problem toolkit,” IEEE Transactions on Evolutionary Computation, vol. 10, no. 5, pp. 477–506, 2006.
 E. Zitzler and L. Thiele, “Multiobjective evolutionary algorithms: a comparative case study and the strength pareto approach,” IEEE Transactions on evolutionary computation, vol. 3, no. 4, pp. 257–271, 1999.
 L. While, L. Bradstreet, and L. Barone, “A fast way of calculating exact hypervolumes,” IEEE Transactions on Evolutionary Computation, vol. 16, no. 1, pp. 86–95, 2012.
 R. G. Steel, D. JH Dickey et al., Principles and procedures of statistics a biometrical approach. WCB/McGrawHill, 1997, no. 519.5 S813 1997.
 K. Deb and M. Goyal, “A combined genetic adaptive search (geneas) for engineering design,” Computer Science and Informatics, vol. 26, pp. 30–45, 1996.
 A. Auger, J. Bader, D. Brockhoff, and E. Zitzler, “Theory of the hypervolume indicator: optimal distributions and the choice of the reference point,” in Proceedings of the tenth ACM SIGEVO workshop on Foundations of genetic algorithms. ACM, 2009, pp. 87–102.
 H. Ishibuchi, Y. Hitotsuyanagi, N. Tsukamoto, and Y. Nojima, “Manyobjective test problems to visually examine the behavior of multiobjective evolution in a decision space,” in Parallel Problem Solving from Nature, PPSN XI. Springer, 2010, pp. 91–100.