A New Dominance Relation-Based Evolutionary Algorithm for Many-Objective Optimization

Many-objective optimization has posed a great challenge to the classical Pareto dominance-based multiobjective evolutionary algorithms (MOEAs). In this paper, an evolutionary algorithm based on a new dominance relation is proposed for many-objective optimization. The proposed evolutionary algorithm aims to enhance the convergence of the recently suggested nondominated sorting genetic algorithm III by exploiting the fitness evaluation scheme in the MOEA based on decomposition, but still inherit the strength of the former in diversity maintenance. In the proposed algorithm, the nondominated sorting scheme based on the introduced new dominance relation is employed to rank solutions in the environmental selection phase, ensuring both convergence and diversity. The proposed algorithm is evaluated on a number of well-known benchmark problems having 3-15 objectives and compared against eight state-of-the-art algorithms. The extensive experimental results show that the proposed algorithm can work well on almost all the test functions considered in this paper, and it is compared favorably with the other many-objective optimizers. Additionally, a parametric study is provided to investigate the influence of a key parameter in the proposed algorithm.

of objectives indeed appear widely in many real-world applications, e.g., control system design [3], [4], industrial scheduling [5], [6], and software engineering [7], [8]. Hence, the practitioners are in need of an effective optimizer to solve these problems at hand. On the other hand, the popular Pareto dominance-based multiobjective evolutionary algorithms (MOEAs), such as nondominated sorting genetic algorithm II (NSGA-II) [9], strength Pareto evolutionary algorithm 2 (SPEA2) [10], and Pareto envelope-based selection algorithm II (PESA-II) [11], have encountered great difficulties in many-objective optimization, although they have shown excellent performance on problems with two or three objectives. The primary reason is that almost all the solutions in the population become nondominated with the number of objectives increasing, which would lead to the severe loss of Pareto-based selection pressure toward the Pareto front (PF). This difficulty has been pointed out both analytically and experimentally in the early studies [12]- [14] on evolutionary many-objective optimization.
To overcome the drawback of Pareto dominance-based MOEAs, some efforts have been made in this paper. In summary, the developing techniques can be roughly classified into the following three types.
1) Adoption of New Preference Relations: Since the Paretodominance relation scales poorly in many-objective optimization, it is natural to use other preference relations, including modified Pareto dominance and different ranking schemes, so as to produce fine selection pressure toward PF. Up to now, many alternative preference relations, such as favor relation [15], dominance [16], [17], fuzzy Pareto dominance [18], [19], preference order ranking [20], and so on [21]- [25], have been proposed.

2) Adoption of New Diversity Promotion Mechanisms:
In many-objective optimization, the Pareto dominance could not provide sufficient selection pressure to make progress in a given population, so the diversity selection mechanism begins to play a key role in such cases. This phenomenon is so called active diversity promotion [26]. Some experimental observations [1], [2], [26], [27] have indicated that it has the potential detrimental effect on the convergence of MOEAs because most of the existing diversity criteria, e.g., crowding distance [9], tend to favor dominance resistant solutions [12] (i.e., the solutions with high performance in at least one of the objectives, but with especially poor performance in the rest of the objectives). Consequently, the final obtained solutions This work is licensed under a Creative Commons Attribution 3.0 License. For more information, see http://creativecommons.org/licenses/by/3.0/ may present a good diversity over the objective space but with the poor proximity to the global PF. To weaken or avoid such effect, new mechanisms promoting the diversity are thus needed. In contrast to the first type, there has not been much research in this respect yet. Adra and Fleming [28] presented two mechanisms for managing diversity and investigated their impact on the overall convergence in many-objective optimization. Deb and Jain [29] proposed an improved NSGA-II procedure, i.e., NSGA-III, which replaces the crowding distance operator in NSGA-II with a clustering operator aided by a set of well-distributed reference points.
Li et al. [30] developed a general modification of the diversity criterion for Pareto dominance-based MOEAs, i.e., the shift-based density estimation (SDE) strategy, covering both distribution and convergence information of individuals. 3) Adoption of two separate archives for convergence and diversity: This kind of approach separates the nondominated solutions into two archives, which promote convergence and diversity during the evolutionary process, respectively. The Two-Archive algorithm [31] is the first MOEA based on this idea. Recently, Wang et al. [32] suggested an improved one for many-objective optimization. It is worth noting that, unlike Pareto dominance-based MOEAs, decomposition-and indicator-based MOEAs, have been found to be very promising in many-objective optimization, although they also have their own drawbacks [25]. The former decomposes a problem with multiple objectives into a set of single-objective subproblems through aggregation functions, and then solves these subproblems simultaneously by evolving a population of solutions. MOEA based on decomposition (MOEA/D) [33], [34] is the most typical implementation of this class. So far, lots of studies based on MOEA/D have been conducted from various aspects, e.g., combing it with the swarm intelligence [35], [36], hybridizing it with local search [37], [38], and incorporating selfadaptation mechanisms [39]- [41].
The indicator-based approach employs a single performance indicator to optimize a desired property of evolutionary population. Among the current indicators available, the hypervolume [42] is probably the most popular one that is ever used in multiobjective search. This is mainly due to its good theoretical properties. It has been indicated that maximizing the hypervolume indicator is equivalent to finding the PF [43]. Nowadays, some hypervolume-based MOEAs have been well-established, such as indicator-based evolutionary algorithm [44], S metric selection evolutionary algorithm [45], and multiobjective covariance matrix adaptation evolution strategy [46]. Nevertheless, the computation cost of the hypevolume grows exponentially with the number of objectives [47], generally inhibiting the use of these algorithms for problems having more than five objectives [48]. To relieve this issue, researchers have been trying to calculate exactly the hypervolume in more efficient ways [49]- [51]. But these ways are still not efficient enough to satisfy the requirement of hypervolume-based MOEAs when solving many-objective problems. Bader and Zitzler [52] suggested the hypervolume estimation (HypE) algorithm for multiobjective optimization, where the exact hypervolume values are approximated by using Monte Carlo simulation. Since the exact calculation of hypervolume is avoided in HypE, it renders the hypervolume-based search possible in highdimensional objective space to some extent. Very recently, several other performance indicators, such as R2 [53]- [55] and the additive approximation [56], [57], have also shown potentials in guiding the many-objective search.
The reference point-based MOEAs are another class of approaches that deserve to be highlighted in many-objective optimization. Slightly different from decomposition-based MOEAs, they perform the predefined multiple targeted search by means of multiple predefined reference points instead of multiple search directions, which can effectively alleviate several difficulties in handling many objectives [29]. Figueira et al. [58] proposed a multiple reference point approach that can be divided into two consecutive phases. The first phase is the preparation phase and it is devoted to estimate the bounds of the PF, generate multiple reference points, and design a version of the solver for each reference point. The second phase is the running phase and it launches a solver for every reference point in every processor. Moen et al. [59] presented a taxi-cab surface evolutionary algorithm, where the Manhattan distance is adopted as the basis for generating the attraction points and as the single metric for selecting solutions for the next generation. Some other studies with this regard can be referred in [29] and [60]- [62]. In particular, the recently proposed NSGA-III also falls into this class.
Despite a number of recent achievements in this field, the research on evolutionary many-objective optimization is far from being fully explored. For the algorithm design, the existing state-of-the-art MOEAs specially for many-objective problems are still not powerful enough [23], [63], and the need for more effective algorithms is pressing. Moreover, although various evolutionary many-objective approaches have been proposed, there are few comparative studies of different methods available to date. In [63], eight evolutionary algorithms were compared for many-objective optimization, but the comparison was restricted to problems with only two kinds of objective space dimensions, and the latest algorithms, such as NSGA-III [29], were not investigated in their work.
Our main contribution is twofold. First of all, a simple but effective θ dominance-based evolutionary algorithm (θ -DEA), is proposed for many-objective optimization. This algorithm is motivated by the strength and weakness of two recently suggested many-objective optimizers (NSGA-III and MOEA/D). NSGA-III emphasizes population members that are Pareto nondominated but are close to the reference line of each reference point. Nevertheless, when the number of objectives is high, the Pareto-dominance relied on by NSGA-III lacks enough selection pressure to pull the population toward PF, therefore, NSGA-III indeed stresses diversity more than convergence in such cases. MOEA/D implicitly maintains the diversity via the diverse weight vectors, and it could generally approach the PF very well by means of the aggregation function-based selection operator, even in high-dimensional objective space. However, in MOEA/D, whether a new solution replaces an old solution or not is completely determined by their aggregation function values. In many-objective optimization, such replacement may lead to the severe loss of diversity in MOEA/D. The major reason is that, in highdimensional objective space, it is highly possible that a solution achieves a good aggregation function value but is far away from the corresponding weight vector. Thus, MOEA/D is at high risk of missing some search regions if aggregation function values are emphasized too much. The issue on the loss of diversity in MOEA/D has been experimentally observed in several recent studies on many-objective optimization [23], [29], [62], [63].
We aim to improve the convergence of NSGA-III in manyobjective optimization by exploiting the fitness evaluation scheme in MOEA/D, but still inherit the strength of NSGA-III in preserving the diversity. To this end, a new dominance relation, referred to as θ dominance, is introduced in this paper. In θ dominance, solutions are allocated into different clusters represented by well-distributed reference points. Only the solutions within the same cluster have the competitive relationship, where a fitness function similar to penalty-based boundary intersection (PBI) function [33] is carefully defined. When conducting the environmental selection in θ -DEA, the nondominated sorting scheme [64] based on θ dominance not only prefers solutions with better fitness values in each cluster, but also ensures that the selected solutions distribute as evenly as possible between these clusters.
The second contribution of this paper lies in the experimental aspect. We provide an extensive comparison between the proposed θ -DEA with eight state-of-the-art algorithms on 80 instances of 16 test problems taken from two well-known test suites. The results indicate that θ -DEA is a very promising algorithm for many-objective optimization.
The remainder of this paper is organized as follows. Section II introduces the background knowledge of this paper. The proposed θ -DEA is described in detail in Section III. Section IV presents the test problems, quality indicators, and algorithm settings used for performance comparison. Section V provides the extensive experimental results and discussion. Finally, the conclusion is drawn in Section VI.

II. PRELIMINARIES
In this section, some basic definitions in multiobjective optimization are first given. Then, we will briefly introduce the original MOEA/D and NSGA-III, which are the basis of our proposed algorithm.

A. Basic Definitions
The multiobjective optimization problem (MOP) can be mathematically defined as consists a set of m objective functions, and is a mapping from n-dimensional decision space to m-dimensional objective space .
Definition 2: A decision vector x * ∈ is Pareto optimal iff there is no x ∈ such that x ≺ x * .
Definition 3: The Pareto set (PS) is defined as Definition 4: The PF is defined as The goal of MOEAs is to move the nondominated objective vectors toward PF (convergence), and also generate a good distribution of these vectors over the PF (diversity).

B. MOEA/D
The key idea of MOEA/D is to decompose an MOP into a number of single-objective optimization subproblems through aggregation functions. It aims to optimize these subproblems in parallel instead of trying to directly approximate the true PF. This mechanism works since the optimal solution to each subproblem is indeed a Pareto optimal solution to the given MOP. And the collection of these optimal solutions can be viewed as an approximation of the true PF. Generally, three aggregation functions, weighted sum, Chebyshev, and PBI function can well serve the purpose in MOEA/D. Just take the Chebyshev as an example, let λ 1 , λ 2 , . . . , λ N be a set of evenly spread weight vectors, then an MOP can be decomposed into N single-objective subproblems represented as where j = 1, 2, . . . , N and λ j = (λ j,1 , λ j,2 , . . . , λ j,m ) T . In MOEA/D, for each vector λ j , a set B( j) = { j 1 , j 2 , . . . , j T } is computed in the initialization phase, where {λ j 1 , λ j 2 , . . . , λ j T } is a set of T closest weight vectors to λ j according to the Euclidean distance and is also called the neighborhood of λ j . The neighborhood of the jth subproblem contains all the subproblems with weight vectors from the neighborhood of λ j . At each generation of MOEA/D, a population of N solutions x 1 , x 2 , . . . , x N are maintained, where x j is the current solution to the jth subproblem. One of the feature of MOEA/D is that the mating restriction is adopted in the reproduction phase. When producing the jth offspring, two indexes k and l are randomly selected from B( j), and a new solution y is generated from x k and x l by using genetic operators. The local replacement is another feature of MOEA/D. That is once y is obtained, it will be compared with each neighboring solution x u , u ∈ B( j), and x u is to be replaced by y, if and only if g te (y|λ u , z * ) < g te (x u |λ u , z * ). After all the N offsprings are produced in this manner one after another, a new population is formed, and the above procedure is repeated until the stopping criterion is met. Further details of MOEA/D can be referred in [33] and [34].

C. NSGA-III
The basic framework of NSGA-III remains similar to the established NSGA-II [9] with significant changes in its selection mechanism. The main procedure of NSGA-III can be briefly described below.
NSGA-III starts with the definition of a set of reference points. Then an initial population with N members is randomly generated, where N is the population size. The next steps are iterated until the termination criterion is satisfied. At the tth generation, the current parent population P t is used to produce an offspring population Q t by using random selection, simulated binary crossover (SBX) operator and polynomial mutation [65]. The size of P t and Q t are both N. Thereafter, the two populations P t and Q t are merged together to form a new population R t = P t ∪ Q t (of size 2N). To choose the best N members from R t for the next generation, the nondominated sorting based on Pareto dominance is first used, which classifies R t into different nondomination levels (F 1 , F 2 , and so on). Then, a new population S t is constructed by filling members of different nondomination levels one at a time, starting from F 1 , until the size of S t is equal to N or for the first time becomes greater than N. Let us suppose that the last level included is the lth level. Hence, the solutions from the level l + 1 onward are simply rejected. Members in S t \ F l are already chosen for P t+1 , and the remaining population slots are chosen from F l such that a desired diversity is maintained in the population. In the original NSGA-II, the solutions in F l with the largest crowding distance values are selected. However, the crowding distance measure does not perform well on many-objective problems [66]. Thus, NSGA-III uses a new selection mechanism that conducts a more systematic analysis of members in S t with respect to the supplied reference points.
To achieve this, objective values and supplied reference points are first normalized so that they have an identical range. After normalization, the ideal point of the set S t is the zero vector. Thereafter, the perpendicular distance between a member in S t and each of the reference lines ( joining the ideal point with a reference point) is calculated. Each member in S t is then associated with a reference point having the minimum perpendicular distance. Next, the niche count ρ j for the jth reference point, defined as the number of members in S t \F l that are associated with the jth reference point, is computed for further processing. Now, a niche-preservation operation is executed to select members from F l , and it works as follows. First, the reference point set J min = { j : argmin j ρ j } having the minimum ρ j value is identified. In case of |J min | > 1, onē j ∈ J min is randomly chosen. If the level F l does not have any member associated with thejth reference point, the reference point is excluded from further consideration for the current generation, meanwhile, J min is recomputed andj is reselected.
Otherwise, the value of ρ¯j is further considered. If ρ¯j = 0, we choose the one having the shortest perpendicular distance to thejth reference line among members associated with thejth reference point in F l , and add it to P t+1 . The count of ρ¯j is then increased by one. In the event of ρ¯j ≥ 1, a randomly chosen member from level F l that is associated with thejth reference point is added to P t+1 , and the count of ρ¯j also needs to be increased by one. The above niche operation is repeated until the remaining population slots of P t+1 are filled. For more details of NSGA-III, please refer to [29].

A. Overview
The framework of the proposed θ -DEA is described in Algorithm 1. First, a set of N reference points are generated, which can be denoted as . . , m, and m k=1 λ j,k = 1. Next, the initial population P 0 with N members is randomly produced. The ideal point z * is initialized in step 3. Since it is often very time consuming to compute exact z * i , it is indeed estimated by the minimum value found so far for objective f i , and is updated during the search. The nadir point z nad is initialized in step 4, where z nad i is assigned to the largest value of f i found in P 0 , and it is updated in the normalization procedure. Steps 6-23 are iterated until the termination criterion is satisfied. In step 7, the offspring population Q t is produced by using the recombination operator. Then Q t is combined with the current population P t , and form a new population R t .
In fact, for problems having high number of objectives, S t is almost always equal to F 1 since there is a large fraction of Pareto nondominated solutions in the population. In step 11, the normalization procedure is executed to S t assisted by z * and z nad . After normalization, the clustering operator is used to split the members in S t into a set of N clusters C = {C 1 , C 2 , . . . , C N }, where the cluster C j is represented by the reference point λ j . Then, the nondominated sorting based on θ dominance (not Pareto dominance) is employed to classify S t into different θ -nondomination levels (F 1 , F 2 , and so on). θ dominance, which is the key concept in θ -DEA, would be introduced later. Once θ -nondominated sorting has been finished, the remaining steps fill the population slots in P t+1 using one level at a time, starting from F 1 . Different from both NSGA-II and NSGA-III, we just randomly select solutions in the last accepted level F l in θ -DEA, because θ dominance has stressed both convergence and diversity. Certainly, some strategies to enhance the diversity could also be used as well in step 20. In the following sections, the important procedures of θ -DEA are to be described in detail.

22:
t ← t + 1 23: end while mechanism is also used in MOEA/D, NSGA-III, and some earlier algorithms [68]. The number of reference points produced in this way depends on the dimension of objective space m and another positive integer H. Let us consider The number of solutions to (5) can be calculated as Suppose that (x j,1 , x j,2 , . . . , x j,m ) T is the jth solution. Then, the reference point λ j is obtained, where Geometrically, λ 1 , λ 2 , . . . , λ N are all located at the hyperplane m i=1 f i = 1, and H is the divisions considered along each objective axis.
Note that, if H < m, no intermediate reference point is created by this approach. However, when m is relatively larger, e.g., m = 8, H ≥ m would lead to a huge number of reference points, and hence a huge population size. To address this issue, in the proposed θ -DEA, we use two-layered reference points with small values of H as suggested in [29]. Suppose the divisions of boundary and inner layers is H 1 and H 2 , respectively, then the population size Fig. 1 illustrates the distribution of two-layered reference points using a three-objective problem with H 1 = 2 and H 2 = 1.

C. Recombination Operator
The recombination operator may be ineffective in manyobjective optimization. This is mainly because, in highdimensional objective space, there is a greater probability that solutions widely distant from each other will be selected to recombine and generate poorer performance offspring solutions, known as lethals [26], [28]. There are normally two ways to address this issue. One is to use the mating restriction scheme [69], in which two neighboring solutions are involved in the recombination operator, such as in MOEA/D. The other is to use a special recombination scheme (e.g., SBX operator with a large distribution index) [29], where near-parent solutions are emphasized, such as in NSGA-III.
θ -DEA employs the latter since it has a similar algorithm structure to NSGA-III. When performing the recombination, two parent solutions are randomly selected from the current population P t , then the child solution is created by using the SBX operator with a large distribution index and polynomial mutation.

D. Adaptive Normalization
The normalization procedure is incorporated into θ -DEA for solving problems having the PF whose objective values may be disparately scaled. In normalization, the objective f i (x), i = 1, 2, . . . , m, can be replaced bỹ As mentioned before, z * i can be estimated by the best value found so far for objective f i . However, the estimation of z nad i is a much more difficult task since it requires information about the whole PF [70], [71]. The procedure of estimating z nad i in θ -DEA is similar to that in NSGA-III but different in the identification of extreme points.
First, in the population to be normalized, i.e., S t , the extreme point e j in the objective axis f j is identified by finding the solution x ∈ S t that minimizes the following achievement scalarizing function: In (10), Thereafter, the value of z nad i is updated as a i , where i = 1, 2, . . . , m, and the population S t can be normalized using (9). In Fig. 2, we demonstrate the hyperplane construction and the intercepts formation in 3-D objective space.
Note that, if the rank of matrix E is less than m, the m extreme points will fail to constitute a m-dimensional hyperplane. And even with the hyperplane built, it is also likely to obtain no intercepts in certain directions or some intercepts a i do not satisfy a i > z * i . In all the above cases, z nad i is assigned to the largest value of f i in the nondominated solutions of S t , for each i ∈ {1, 2, . . . , m}.

E. Clustering Operator
In θ -DEA, the clustering operator is applied to population S t at each generation. The clustering works in the normalized for j ← 2 to N do 6: if d j,2 (x) < min then 7: min ← d j,2 (x) 8: n ← j 9: end if 10: end for 11: . . ,f m (x)) T is the normalized objective vector for the solution x, L is a line passing through the origin with the direction λ j , and u is the projection off(x) on L. Let d j,1 (x) be the distance between the origin and u, and d j,2 (x) be the perpendicular distance betweenf(x) and L. They can be computed, respectively, as In Fig. 3, the distances d j,1 (x) and d j,2 (x) are illustrated in the 2-D objective space. For the clustering operator, only d j,2 will be considered, d j,1 will be involved later in the definition of θ dominance. We assign a solution x to the cluster C j with the minimum d j,2 (x) value. The details of the clustering process are shown in Algorithm 2.

F. θ Dominance
The proposed θ dominance is defined on population S t with the supply of a set of reference points . And each solution in S t is associated with a cluster among a set of clusters C by the clustering operator.
where θ is a predefined penalty parameter. The form of F j (x) is the same as that of the PBI function [33]. But here, the distances d j,1 and d j,2 are both computed in the normalized objective space. Generally, d j,2 (x) = 0 ensures that f(x) is always in L, resulting in perfect diversity, and smaller d j,1 (x) value under the condition d j,2 (x) = 0 means better convergence. With the definition of F j and C, several concepts related to θ dominance can be defined as follows.
Definition 7: Given two solutions x, y ∈ S t , x is said to θ -dominate y, denoted by x ≺ θ y, iff x ∈ C j , y ∈ C j , and Definition 9: All solutions that are θ -optimal in S t form the θ -optimal set (θ -OS), and the corresponding mappings of θ -OS in the objective space form the θ -optimal front.
Based on the definition of θ dominance, we have the following three properties, which respectively illustrate that the relation ≺ θ is irreflexive, asymmetric, and transitive.
. So, the supposition is invalid, and the proposition is true.
Property 3: If three solutions x, y, z ∈ S t satisfy x ≺ θ y and y ≺ θ z, then x ≺ θ z.
Proof: By x ≺ θ y, we have that ∃j ∈ {1, 2, . . . , N}, x ∈ C j , y ∈ C j , and F j (x) < F j (y). Then according to y ≺ θ z, we have that z ∈ C j and F j (y) < F j (z). Overall, x ∈ C j , z ∈ C j , and F j (x) < F j (z). Thus, x ≺ θ z.
Due to the above three properties, the θ dominance defines a strict partial order on S t . Hence, the fast nondominated sorting approach [9] can be immediately adopted in θ -dominance sense, and the population S t would be partitioned into different θ -nondomination levels.
Note that, there is no competitive relationship between clusters in θ dominance, and thus it can indeed use different θ values in different clusters. To explore this characteristic, we do a bit more work on the proposed θ -DEA. In the normalized objective space, if λ j is the axis direction, then we assign a large θ value (θ = 10 6 is used) in the cluster C j , otherwise assign a normal θ value. This is just to cooperate with the normalization procedure presented in Section III-D, and if normalization is disabled, it may be unnecessary. The large θ values in the clusters represented by axis directions would make θ -DEA more likely to capture the nadir point in highdimensional objective space, and thus conduct a more stable normalization.

G. Computational Complexity of θ -DEA
The computational complexity of θ -DEA in one generation is dominated by the clustering operator that is described in Algorithm 2 under the general condition. In Algorithm 2,

IV. EXPERIMENTAL DESIGN
This section is devoted to the experimental design for investigating the performance of the proposed θ -DEA. First, the test problems and the quality indicators used in our experiments are given. Then, we briefly introduce eight state-of-the-art algorithms that are employed for comparison. Finally, the experimental settings adopted in this paper are provided.

A. Test Problems
As a basis for the comparisons, two well-known test suites for many-objective optimization, Deb-Thiele-Laumanns-Zitzler (DTLZ) [72] and Walking Fish Group (WFG) [73], are involved in the experiments. To compute the quality indicators reliably, we only consider DTLZ1-4 and DTLZ7 problems for DTLZ test suite, since the nature of DTLZ5 and DTLZ6s PFs is unclear beyond three objectives [73]. Moreover, we also use two scaled test problems, i.e., scaled DTLZ1 and DTLZ2 problems [29], which are the modifications of DTLZ1 and DTLZ2 problems, respectively. To illustrate, if the scaling factor is 10 i , the objectives f 1 , f 2 , and f 3 for three-objective scaled DTLZ1 problem are multiplied by 10 0 , 10 1 , and 10 2 , respectively. In our experiments, we call scaled DTLZ1 and DTLZ2 problems SDTLZ1 and SDTLZ2 for short, respectively.
All these problems can be scaled to any number of objectives and decision variables. We consider the number of objectives m ∈ {3, 5, 8, 10, 15}. For DTLZ1-4, DTLZ7, SDTLZ1, and SDTLZ2 problems, the total number of decision variables is given by n = m + k − 1, unless otherwise specified, k is set to 5 for DTLZ1 and SDTLZ1, 10 for DTLZ2-4 and SDTLZ2, and 20 for DLTZ7 as recommended in [29] and [72]. As for all WFG problems, unless otherwise stated, the number of decision variables is set to 24 and the position-related parameter is set to m − 1 according to [54] and [73]. The scaling factors for SDTLZ1 and SDTZL2 problems with different number of objectives are shown in Table I.
These test problems have a variety of characteristics, such as having linear, mixed (convex/concave), multimodal, disconnected, degenerate, and disparately scaled PFs, which  Table II, we summarize the main features of all the adopted test problems.

B. Quality Indicators
The quality indicators are needed to evaluate the performance of the concerned algorithms. In the EMO literature, the inverted generational distance (IGD) [74] is one of the most widely used indicators, which could provide a combined information about convergence and diversity of a solution set. The calculation of IGD requires a set of uniformly distributed points along the known PF. Surely, the true PFs of most benchmark problems are known, and it is relatively easy to uniformly sample the points on the 2-or 3-D PF. However, for a high-dimensional PF, how to make sure the sampled points are uniformly distributed and how many points are enough for representing the true PF, are themselves difficult questions. In fact, numerous studies [23], [30], [63] relating to many-objective optimization, where IGD is used as the quality indicator, failed to indicate how they sampled those points along the PF.
Recently, Deb and Jain [29] suggested a way to compute IGD for MOEAs in which the reference points or reference directions are supplied, e.g., MOEA/D and NSGA-III. This way works as follows. For each reference direction λ j , j = 1, 2, . . . , N, we can exactly locate its targeted point v j on the known PF in the normalized objective space. All N targeted points constitute the set V = {v 1 , v 2 , . . . , v N }. For any algorithm, let A be the set of final nondominated points obtained in the objective space. Then IGD is computed as where d(v i , f) is the Euclidean distance between the points v i and f. Note that, for the scaled problems, the objective values in the set A should be first normalized using the ideal and nadir points of the exact PF before computing IGD. The set A with smaller IGD values is better.
It makes sense to compute IGD using the above method for reference point/direction-based MOEAs. This is mainly because that the many-objective optimization task for these algorithms can be seen as finding the Pareto-optimal points close to the supplied reference points to some extent. Since the proposed θ -DEA is also based on reference points, we will evaluate and compare its performance using IGD defined by (14) in the experiments.
However, such IGD is not applicable to MOEAs without using reference points/directions, e.g., HypE [52] and SDE [30]. For these algorithms, the many-objective optimization task is to search for sparsely distributed Pareto-optimal points over the entire PF [75]. In this scenario, another popular indicator, i.e., hypervolume [42], is adopted to evaluate the performance. The hypervolume is strict Pareto-compliant [74], whose nice theoretical qualities make it a rather fair indicator. Let A be the set of final nondominated points obtained in the objective space by an algorithm, and r = (r 1 , r 2 , . . . , r m ) T be a reference point in the objective space which is dominated by any point in the set A. Then the hypervolume indicator value of A with regard to r is the volume of the region dominated by A and bounded by r, and can be described as HV can measure both convergence and diversity of a solution set in a sense. Given a reference point r, larger HV value means better quality. In the calculation of HV, the choice of reference point is a crucial issue. It has been found that choosing r i that is slightly larger than z nad i is suitable since the balance between convergence and diversity of the solution set is well emphasized [76], [77]. In our experiments, we set r to 1.1z nad , where z nad can be analytically obtained for all the adopted test problems. Following the practice in [2] and [22], the points that do not dominate reference point are discarded for the HV calculation. Considering the problems that have PFs with differently scaled objective values, we first normalize the objective values of points in A and the reference point r using z nad and z * (z * is 0 for all the adopted test problems) prior to computing HV by (15). Thus, the presented HV value for a m-objective instance in the experiments is between 0 and 1.1 m −V m , where V m is the hypervolume of the region enclosed by the exact normalized PF and the coordinate axes. In addition, for problems with no more than ten objectives, we calculate HV exactly using the recently proposed WFG algorithm [50]. As for problems having 15 objectives, we approximate the HV by the Monte Carlo simulation proposed in [52], and 10 000 000 sampling points are used to ensure the accuracy.

C. Other Algorithms in Comparison
To verify the proposed θ -DEA, the following eight state-of-the-art algorithms are considered as the peer algorithms.
1) Grid based evolutionary algorithm (GrEA) [23]: It exploits the potential of a grid to strengthen the selection pressure toward the PF while maintaining an extensive and uniform distribution of solutions. To this end, two concepts (i.e., grid dominance and grid difference), three grid-based criteria, (i.e., grid ranking, grid crowding distance, and grid coordinate point distance), and a fitness adjustment strategy are incorporated into GrEA. 2) Preference ordering genetic algorithm (POGA) [20]: It uses the preference order-based approach as an optimality criterion in the ranking stage of MOEAs. By exploiting the definition of efficiency of order in the subsets of objectives, a ranking procedure is employed within the framework of NSGA-II, which exerts the higher selection pressure over objective spaces of different dimensionality compared with the traditional Pareto dominance-based ranking scheme. 3) NSGA-III [29]: It is a reference point-based evolutionary algorithm for many-objective optimization, whose framework is still similar to the original NSGA-II. But unlike NSGA-II, the maintenance of diversity among population members is aided by providing and adaptively updating a number of well-distributed reference points. Overall, NSGA-III emphasizes population members which are nondominated yet close to a set of supplied reference points. 4) SDE [30]: It is a general modification of the density estimation strategies that could make the Pareto dominance-based MOEAs suitable for many-objective optimization. Its basic idea is that, given the preference of density estimators for solutions in sparse regions, the solutions with poor convergence are put into crowded regions by SDE, so that they will be assigned a high density value and then be eliminated easily during the evolutionary process. SDE is simple in implementation and can be applied to any specific density estimator with negligible computational cost and no additional parameters. In this paper, the version that integrates SDE into SPEA2 (SPEA2+SDE) is used, since it shows the best overall performance among all the three considered versions (NSGA-II+SDE, SPEA2+SDE, and PESA-II+SDE) in [30]. 5) MOEA/D [33]: It is a representative algorithm that belongs to the decomposition-based approach. In [34], a new version of MOEA/D (MOEA/D-DE) based on differential evolution (DE) [78] was proposed to deal with problems with complicated PSs. In this paper, the original MOEA/D with the PBI function is selected, since it has been reported in [29] that MOEA/D-DE shows poor performance on many-objective problems and PBI is more suitable for solving problems having a high-dimensional objective space. 6) Decomposition-based multiobjective particle swarm optimizer (dMOPSO) [35]: It is an MOEA that extends the particle swarm optimization [79] technique to the decomposition-based multiobjective approach. It updates the position of each particle using a set of solutions considered as the global best according to the decomposition-based approach. And it is mainly characterized by the use of a memory reinitialization process aiming to provide diversity to the swarm. Similar to MOEA/D, the PBI function is also chosen in dMOPSO. HypE makes a tradeoff between the accuracy of the estimates and the available computing resources, making the hypervolume-based search more easily applied to many-objective problems.

8) Many-objective metaheuristic based on the R2 indicator (MOMBI) [54]:
It is a many-objective metaheuristic based on the R2 indicator. To use R2 in the selection mechanism, a nondominated scheme based on the adopted utility functions is designed for that purpose. Its main idea is that the solutions which achieve the best values on any of the chosen utility functions are given the first rank. Such solutions are then removed and a second rank will be identified in the same manner. The process will continue until all the solutions have been ranked. These eight algorithms have covered most categories of techniques mentioned in Section I for many-objective optimization. Table III summarizes the list of selected algorithms for comparison. All the concerned algorithms, including the proposed θ -DEA are implemented in the jMetal framework [80], and run on an Intel 2.83 GHz Xeon processor with 15.9 Gb of RAM. For all the algorithms except SDE, the obtained final population is used for computing quality indicators, whereas for SDE, the final archive is used.

D. Experimental Settings
The experimental settings include general settings and parameter settings. The general settings are listed as follows.
1) Number of Runs: Each algorithm is run 20 times independently for each test instance. 2) Termination Criterion: The termination criterion of an algorithm for each run is specified in the form of the maximum number of generations (MaxGen). Since the used test problems are of varying computational complexity, we use different MaxGen for different problems.  [81] at a 5% significance level is carried out on the assessment results obtained by two competing algorithms. As for the parameter settings, several common settings for algorithms are first given as follows.
1) Population Size: The setting of population size N for NSGA-III, MOEA/D, and θ -DEA cannot be arbitrarily specified, where N is controlled by a parameter H [see (6)]. Moreover, as presented in Section III-B, we use two-layered reference points for problems having 8, 10, and 15 objectives to create intermediate reference points. For the other algorithms, the population size can be set to any positive integer, but for ensuring a fair comparison, the same population size is adopted. Table IV lists the population sizes used in this paper for the problem with different number of objectives. Noting that, POGA, HypE, MOMBI, and θ -DEA have similar framework to that of NSGA-III, we slightly adjust the population size of these algorithms to the multiple of four just as in the original NSGA-III study [29], i.e., 92, 212, 276, and 136 for 3-, 5-, 10-, and 15-objective problems. 2) Penalty Parameter θ : Since PBI function is involved in MOEA/D, dMOPSO, and θ -DEA, the penalty parameter θ needs to be set for them. In this paper, θ is just set to 5 for both MOEA/D and dMOPSO as suggested in [33], whereas some studies [75], [82] have indicated that the good specification of θ for MOEA/D may depend on the problem to be solved and its number of objectives. For the proposed θ -DEA, θ is also set to 5, but we will investigate the influence of θ on the performance of the proposed θ -DEA in Section V-C.

3) Parameters for Crossover and Mutation:
The SBX and polynomial mutation [65] are used in all the considered algorithms except dMOPSO. For GrEA, POGA, SDE, MOEA/D, HypE, MOMBI, the parameter values for crossover and mutation are presented in Table V. As for NSGA-III and θ -DEA, the settings are only a bit different according to [29], where η c is set to 30. Besides the parameters mentioned above, GrEA, SDE, MOEA/D, dMOPSO, and HypE have their specific parameters. These parameters are set mainly according to the suggestions given by their developers, which are shown below.
1) Parameter Setting in GrEA: The grid division (div) needs to be set. Since the population size and the  [23], we adjust div according to the guidelines provided in [23] for each problem instance as shown in Table VI, aiming to well balance the convergence and diversity.

2) Parameter Setting in SDE:
The archive size is just set as the same as the population size.

3) Parameter Setting in MOEA/D: The neighborhood size
T is set to 20. 4) Parameter Setting in dMOPSO: The age threshold T a is set to 2.

5) Parameter Setting in HypE:
The bound of the reference point is set to 200, and the number of sampling points M is set to 10 000.

V. EXPERIMENTAL RESULTS
In this section, the performance of θ -DEA is to be validated according to the experimental design described in Section IV. Our experiments can be divided into three parts. The first one is to compare θ -DEA with the other two MOEAs with reference points/directions, i.e., NSGA-III and MOEA/D. The aim is to demonstrate the superiority of θ -DEA in achieving the desired convergence and diversity as a reference point-based algorithm. The second one is to compare θ -DEA with various types of many-objective techniques. The aim is to show the great ability of θ -DEA in searching for the sparsely distributed nondominated points over the entire PF as a general many-objective optimizer. The third one is to investigate the influence of parameter θ on the performance of the proposed algorithm. Moreover, we will also make some further comments on our experimental results in Section V-D.

A. Comparison With NSGA-III and MOEA/D
In this section, IGD is used to evaluate the algorithms. Since all the experimental settings including the way to compute IGD are consistent with those in the original NSGA-III study [29], we compare the IGD results of θ -DEA with those of NSGA-III and MOEA/D-PBI taken from [29].
First, the normalized test problems, DTLZ1-4, which have an identical range of values for each objective over the PF, are employed for comparison. The DTLZ1 problem has a simple, linear PF ( m i=1 f i = 0.5) but with 11 5 − 1 local optima in the search space. The difficulty in this problem is to converge to the hyperplane. The PFs of DTLZ 2-4 problems have the same geometrical shape ( m i=1 f 2 i = 1), but they are designed to challenge different capacities of an algorithm. The DTLZ2 problem is a relatively easy problem with a spherical PF. The DTLZ3 problem introduces a huge number of local PFs paralleling to the global PF based on DTLZ2, which poses a stiff challenge for algorithms to converge to the global PF. The DTLZ4 problem challenges the ability of an algorithm to maintain the diversity in the objective space by introducing a variable density of solutions along the PF. Table VII shows the results of θ -DEA and NSGA-III on the four problems, where the best, median, and worst IGD values are reported. From Table VII, θ -DEA performs consistently better than NSGA-III on all the instances except 3, 15-objective DTLZ1 and 3-objective DTLZ2. For 15-objective DTLZ1 and DTLZ3, θ -DEA is unable to get close to the PF in some runs, as evident from the large worst IGD value. We suspect that the reason is that θ -DEA sometimes fails to capture the nadir point in high-dimensional objective space and does a wrong normalization. But for 15-objective DTLZ4, θ -DEA performs very well in all runs, whereas NSGA-III sometimes struggles to maintain good convergence and diversity on DTLZ4 problems having more than five objectives. It is interesting to note that, for three-objective DTLZ4, neither of θ -DEA and NSGA-III could achieve good performance all the time.
MOEA/D-PBI does not incorporate the normalization technique. To compare against it more reasonably, we remove the normalization procedure from θ -DEA, and refer this version as θ -DEA * . Table VIII shows the comparison results between θ -DEA and MOEA/D-PBI. It can be seen that θ -DEA * outperforms MOEA/D-PBI on DTLZ1, DTLZ3, and DTLZ4 problems, whereas MOEA/D-PBI wins on DTLZ2 problem. The advantage of θ -DEA * is particularly obvious on DTLZ4 problem. We also notice that the results of θ -DEA * are generally better than those of θ -DEA shown in Table VII, which indicates that the normalization is not quite necessary for the normalized problems. In addition, unlike θ -DEA, for 15-objective DTLZ1 and DTLZ3, θ -DEA * works well in all runs. This illustrates that the normalization is the bottleneck of the performance of θ -DEA on these two instances to some extent. Next, two WFG problems, WFG6 and WFG7, are used to test the performance of θ -DEA, NSGA-III, and MOEA/D-PBI. The two problems have the same PF shape, which is part of  NSGA-III performs best on WFG6 problems and threeobjective WFG7, whereas the proposed θ -DEA performs best on WFG7 with more than three objectives. MOEA/D-PBI is outperformed by θ -DEA on both WFG6 and WFG7 problems.
To further investigate the performance of θ -DEA on problems with disparately scaled objective values, SDTLZ1 and SDTLZ2 problems are considered. The comparison results are shown in Table X. For SDTLZ1 problem, the situation is similar to DTLZ1 problem. That is, NSGA-III performs better on 3-and 15-objective instances, whereas θ -DEA performs better on the remaining instances. As for SDTLZ2 problem, θ -DEA clearly outperforms NSGA-III.
Based on the above comparisons, it can be concluded that the proposed θ -DEA can generally maintain a good balance between convergence and diversity assisted by structured reference points. Indeed, through the test on problems with varying features, θ -DEA outperforms NSGA-III and MOEA/D-PBI on most of them in terms of IGD values.

B. Comparison With State-of-the-Art Algorithms
In this section, we compare the proposed θ -DEA with all the eight algorithms mentioned in Section IV-C. The version of θ -DEA without normalization, i.e., θ -DEA * , will also be involved in the comparison. The HV indicator is used to evaluate the concerned algorithms.   Table XIII gives a summary of the significance test on HV results between the proposed θ -DEA (θ -DEA * ) and the other algorithms. In this table, Alg 1 versus Alg 2 , "B" ("W") means that the number of instances on which the results of Alg 1 are significantly better (worse) than those of Alg 2 , and "E" indicates the number of instances where there exists no statistical significance between the results of Alg 1 and Alg 2 .
To describe the distribution of obtained solutions in highdimensional objective space, we use 15-objective WFG7 instance as an illustration. Fig. 4 plots the final solutions of four competitive algorithms, i.e., θ -DEA, GrEA, NSGA-III, and SDE, in a single run by parallel coordinates. This particular run is associated with the result closest to the average HV value. It is clear from Fig. 1 that θ -DEA and NSGA-III are able to find a good approximation and coverage of the PF, whereas GrEA and SDE can only converge to a portion of the PF.
To quantify how well each algorithm performs overall, the performance score [52] is introduced to rank the algorithms. For a specific problem instance, suppose there are l algorithms Alg 1 , Alg 2 , . . . , Alg l involved in the comparison, let δ i,j be 1, if Alg j is significantly better than Alg i in terms of HV, and 0 otherwise. Then, for each algorithm Alg i , the performance score P(Alg i ) is determined as This value reveals how many other algorithms significantly outperform the corresponding algorithm on the considered test instance. So, the smaller the value, the better the algorithm. In Fig. 5, the average performance score is summarized for different number of objectives and different test problems. Fig. 6 shows the average performance score over all 80 problem instances for the selected ten algorithms, and the overall rank of each algorithm according to the score is also given in the corresponding bracket. Based on the above results, we can obtain some observations for each algorithm. The proposed θ -DEA works well on nearly all the considered instances. In particular, it shows the best overall performance on DTLZ4, SDTLZ1, SDTLZ2, and WFG4-9 problems. Relatively speaking, θ -DEA does not show such outstanding performance on WFG1-3 problems. For WFG2 problem, θ -DEA exhibits an interesting search behavior, it remains competitive on 3-, 5-, 8-, and 10-objective instances, but performs worst on 15-objective instance. θ -DEA * shows the advantage on the normalized test problems, i.e., DTLZ1-4 problems. But it cannot be compared to θ -DEA when handling scaled problems. Indeed, θ -DEA * is significantly outperformed by θ -DEA on 53 out of 60 scaled problem instances, verifying the effectiveness of the normalization procedure in θ -DEA.
GrEA can effectively deal with scaled problems in general, which is verified by its competitive results on DTLZ7, SDTLZ1, SDTLZ2, and most of WFG problem instances. This is mainly because that GrEA divides each dimension of the objective space into the same number of divisions. So, GrEA indeed does the objective normalization implicitly during the evolutionary process. In addition, it is worth mentioning that the performance of GrEA is sensitive to the parameter div, the overall excellent performance of GrEA in our experiments is obtained by suitably setting div for each test instance. In this regard, GrEA takes advantage of the other compared algorithms.
POGA shows overall competitive performance on WFG2 and WFG3 problems, and it even achieves the best performance on ten-objective WFG2 instance and 5-, 8-, and 10-objective WFG3 instances. But it generally cannot obtain very satisfying results on the other problems. It is interesting to note that, for DTLZ1, DTLZ3, and SDTLZ1 problems, POGA obtains zero HV values on 5-, 8-, and 10-objective instances, but nonzero values for 15-objective instance. Since the PFs of the three problems have a huge number of local PFs, it is not surprising that POGA always fails to converge to the PFs and obtains zero HV values. But the reason for nonzero HV values (especially poor) for 15-objective instances is waiting to be explored.
NSGA-III shows the closest overall performance to the proposed θ -DEA, and it can perform very well over a wide range of test problems. There is no significant difference between It is interesting to find that, for WFG1 problem, NSGA-III performs poorly on 3-, 5-, 8-, and 10-objective instances, but it takes the second place on 15-objective instance. SDE generally has the medium-high performance on most of the considered problems among the compared algorithms, leading to its fourth rank in the average performance score over all the instances as shown in Fig. 6. It is worth noting, for DTLZ7 problem, SDE performs best on three-and five-objective instances, but its performance scales poorly in instances with a higher number of objectives.
MOEA/D achieves good performance on the normalized test problems except DTLZ4, but it cannot produce satisfactory results on almost all the scaled test problems. This is why it ranks poorly in Fig. 6. Nevertheless, from this, we do not claim that MOEA/D is a poor many-objective optimizer, since a naive normalization procedure may even enhance its ability to deal with scaled problems [33]. dMOPSO is a relatively poor algorithm, and it does not perform well even on the normalized problems.
HypE is very competitive on three-objective instances, which is reflected in Fig. 5(a). Indeed, it performs best on 12 out of the 16 three-objective instances. However, it does not show advantage over the other algorithms on problems having more than three objectives except on WFG3, where it performs very well on 8-, 10-, and 15-objective instances. Note that, HypE computes exactly the hypervolume-based fitness values when m ≤ 3, otherwise it estimates the fitness values using Monte Carlo simulation. Thus, we suspect that its relatively poor performance on problems with high number of objectives is mainly due to its inaccurate fitness estimation. Increasing the number of sampling points may improve the situation, but the computational effort will soon become unacceptable. Although, HypE is a popular many-objective optimizer, our experimental results indicate that it is not compared favorably with some newly proposed many-objective algorithms.
The similar observation can also be found in several recent studies [23], [30], [63]. MOMBI performs well on most of WFG1-3 problem instances, especially on WFG1, where it attains the best results on 5-, 8-, 10-, and 15-objective instances. But it does not behave quite well on the other problems.
Since, we have allocated more computational efforts for problems having higher number of objectives, it is not possible to analyze performance scalability. In fact, it is not clear from our experimental results that the performance of each algorithm decays as the number of objectives is increased. For the proposed θ -DEA, it performs as well even on 15-objective instances, which can be seen from the results on the problems whose PFs have regular geometrical shapes, i.e., DTLZ1-4, SDTLZ1, SDTLZ2, and WFG4-9 problems. The average HV values of θ -DEA on these instances are close to 1.1 15 ≈ 4.177, indicating that it can achieve a good approximation of PF even in high-dimensional objective space.   Next, we would like to briefly investigate what will happen if we set a larger number of decision variables (n) for test problems, although this is not the focus of this paper.
As an illustration, we select a normalized problem DTLZ2 and a scaled problem WFG7. For DTLZ2, k is reset to 98, thus n = m + 97. For WFG7, n is reset to 100, and the Fig. 6. Ranking in the average performance score over all test problem instances for the selected ten algorithms. The smaller the score, the better the overall performance in terms of HV. position-related parameter is still set to m − 1. We run three competitive algorithms θ -DEA, NSGA-III, and SDE on the two problems, and the algorithm parameters and the termination criterion are kept unchanged. Table XIV shows the average HV results. It can be seen that all the three algorithms obtain smaller HV values than on the original problem instances, which indicates that larger n would pose greater difficulty to all of them. However, the comparison situation changes to some extent. For example, θ -DEA performs significantly better than SDE on the original 3-, 5-, and 10-objective DTLZ2 instances, but it is significantly outperformed by SDE on these instances with larger n. It seems that, for DTLZ2 problem, SDE is less influenced by the increasing of n than θ -DEA. In addition, it is worth mentioning that a cooperative coevolution technique [83] has been developed specially for solving MOPs with large number of decision variables.
In the end, we intend to gain some insights based on the extensive experimental results provided in this section. First, it should be pointed out that the performance of an algorithm not only depends on its ability to cope with specific problem features but also depends on its ability to handle high number of objectives (e.g., θ -DEA on WFG2 and SDE on DTLZ7). Second, for certain problems, the increasing of the number of objectives would exert different levels of impact on different algorithms (e.g., NSGA-III on WFG1). Third, θ -DEA, NSGA-III, and GrEA show strong competitiveness on most of the problem instances considered in this paper. θ -DEA may have difficulty in handling problems with high-dimensional convex disconnected PFs (e.g., 15-objective WFG2); NSGA-III may struggle on relatively low-dimensional biased problems having mixed PFs (e.g., 3-, 5-, 8-, and 10-objective WFG1); GrEA may not be good at solving problems with a huge number of local PFs (e.g., DTLZ1 and DTLZ3). The other concerned algorithms show advantage on specific problem instances, and their main characteristics can be summarized as follows.

C. Influence of Parameter θ
In this section, we investigate the effect of parameter θ on the performance of the proposed algorithm. The variation of θ would influence the normalization process presented in Section III-D, and hence influence the search behavior of θ -DEA. To observe the pure effect of θ , we hypothesize that the ideal and nadir points are known a priori and the accurate normalization can be done. Therefore, here we decide to show the influence of θ by running θ -DEA * on the normalized test problems, i.e., DTLZ1-4. Similar observation can be obtained from the other test cases.
Figs. 7 and 8 present how the performance of θ -DEA * varies with the change of θ on DTLZ1-4 problems in terms of average IGD and average HV, respectively. We vary θ between 0 and 50 with a step size of 5, and θ = 0.1 is also examined. Based on the two figures, we obtain the following observations. 1) θ = 0 almost always leads to the worst performance.
2) θ = 0.1 is also usually a bad choice, but it seems to be the most suitable setting for three-objective DTLZ3.
3) The most appropriate setting of θ depends on the problem instance to be solved. 4) The performance of θ -DEA * on most of the problem instances is robust over a wide range of θ values, which is beneficial to the practical use of the proposed algorithm. Note that, for eight-objective DTLZ4, the IGD values fluctuate with the increase in θ , but the corresponding HV values remain stable. The way we compute IGD in our experiments makes IGD a more sensitive indicator than HV. For example, it may lead to much larger IGD even if there is only one reference point that cannot be associated well with the found solutions, but the corresponding HV may only get slightly smaller in this case.
It is also interesting to understand how well θ -DEA * performs when θ tends to positive infinity. In this case, only the closeness to the corresponding reference points for each solution is emphasized in the environmental selection phase, which makes the search behavior of θ -DEA * more like NSGA-III. To investigate this extreme case, we set θ to a large value 10 6 . Table XV compares the average IGD and average HV values of θ -DEA * (θ = 5) and θ -DEA * (θ = 10 6 ). It can be seen clearly that θ = 5 generally achieves a better balance between convergence and diversity than θ = 10 6 by means of reference points.

D. Further Discussion
In this section, we further discuss three issues. The first is about the quality indicators adopted in our experiments; the second is about the computational effort of the concerned algorithms; and the third is about the comparison of our experimental findings to the existing ones in this paper.
In our experiments, we adopt two quality indicators for the purpose of comparison, i.e., IGD and HV, both of which can provide a combined information about convergence and diversity of a solution set. IGD used in this paper only applies to MOEAs based on reference points/directions, whereas HV applies to all MOEAs. We find that IGD and HV can usually result in consistent conclusions, but there exist a few exceptions. For example, in Table XV, θ -DEA * (θ = 5) significantly outperform θ -DEA * (θ = 10 6 ) on 8-, 10-, and 15-objective DTLZ1 in terms of IGD, but the HV results show that θ -DEA * (θ = 10 6 ) performs significantly better. The main reason is that the task to find Pareto-optimal points close to the supplied reference points is not completely equivalent to the task to maximize the hypervolume, which strongly depends on the distribution of the reference points. But it is indeed desirable to evaluate reference points/directions based MOEAs using IGD indicator, since IGD can better measure the "closeness" between the outcome and goal for these algorithms. In addition, IGD can also be used in the scenario that the user is interested in only a preferred part of PF [29], where only a few representative reference points are used. It is worth pointing out that if we use more uniformly spread reference points in θ -DEA, it would be beneficial for the task to search for sparsely distributed Pareto-optimal points over the entire PF, and thus may enable θ -DEA to achieve better HV results.
Among the concerned algorithms, we find that MOEA/D and dMOPSO require the least computational effort. θ -DEA, MOMBI, and NSGA-III are also efficient enough to cope with the problem instances under the parameter specifications in this paper. GrEA, POGA, and SDE are generally much more computation expensive than the above mentioned algorithms, whose computational time increases sharply with increasing number of objectives. We would like to specially mention HypE. In the original HypE study [52], the authors used a small population size, i.e., 50, and claimed that HypE is a fast algorithm. However, we find that the computation time of HypE would increase severely not only with the increase in the  TERMS OF AVERAGE IGD AND HV VALUES.  THE PERFORMANCE THAT IS SIGNIFICANTLY BETTER  THAN THE OTHER IS SHOWN IN BOLD number of objectives but also with the increase in population size. Hence, under the parameter settings in our experiments, HypE is indeed the most time-consuming algorithm. The similar observation can be referred in [75], where the computation time of HypE is reported. It should be mentioned that, not just HypE, the efficiency of the other considered algorithms is also influenced by the population size. However, HypE seems to be more affected than the others based on our experimental observations. Although there have been a number of experimental studies on the comparison of many-objective optimizers, the parameter settings, the adopted algorithm variants (e.g., MOEA/D has several variants), the termination criterion, and the problem settings (e.g., the number of objectives and the number of decision variables) usually vary from study to study, making a rigorous comparison between the experimental findings impractical. However, this paper represents one of the most comprehensive experimental comparisons in the literature so far. A number of different research issues have been investigated, as presented in Section V.

VI. CONCLUSION
In this paper, we have presented a new many-objective evolutionary algorithm, called θ -DEA, whose environmental selection mechanism is based on θ dominance. Given the complementary advantages of NSGA-III and MOEA/D, θ -DEA is expected to enhance the convergence ability of NSGA-III in high-dimensional objective space by utilizing the aggregation function-based fitness evaluation scheme in MOEA/D, thereby achieving a better compromise between convergence and diversity in many-objective optimization. To achieve the goal, a new dominance relation, θ dominance, is introduced into the proposed algorithm, managing to emphasize both convergence and diversity.
In the experimental investigation of θ -DEA, we have shown the evidence of the robustness of the algorithm with respect to the key parameter θ . We have experimentally verified that the θ -DEA, in general, performs better than NSGA-III and MOEA/D in searching for the Pareto-optimal points close to the supplied reference points as a reference point-based algorithm. It has also been found that the embedded normalization procedure enables θ -DEA to handle scaled problems more effectively.
To demonstrate the strong competitiveness, we have made an extensive experimental comparison of θ -DEA with eight state-of-the-art algorithms that belong to five different categories of technologies. A number of well-known benchmark problems are chosen to challenge different abilities of the algorithms. The comparison results reveal that the proposed θ -DEA works well on almost all the problem instances considered in this paper, and it is compared favorably with state-of-the-art many-objective optimizers. However, we have also observed that none of the algorithms is capable of beating any of the other algorithms on all the instances, which indicates that a careful choice of algorithms is sill needed at present when solving a many-objective problem at hand.
In the future, we will have a deeper insight into the search behavior of θ -DEA, so as to further improve its performance. It is also necessary to evaluate θ -DEA further by comparing it to the two most recent many-objective optimization algorithms [25], [32], so that the strength and weakness of different algorithms can be better understood. It would also be interesting to extend our θ -DEA to solve constrained many-objective problems by incorporating constraint handling techniques [84]. Moreover, we would apply θ -DEA to real-world problems in order to further verify its effectiveness.