A New Hypervolume-Based Evolutionary Algorithm for Many-Objective Optimization

In this article, a new hypervolume-based evolutionary multiobjective optimization algorithm (EMOA), namely, R2HCA-EMOA (R2-based hypervolume contribution approximation EMOA), is proposed for many-objective optimization. The core idea of the algorithm is to use an R2 indicator variant to approximate the hypervolume contribution. The basic framework of the proposed algorithm is the same as SMS-EMOA. In order to make the algorithm computationally efficient, a utility tensor structure is introduced for the calculation of the R2 indicator variant. Moreover, a normalization mechanism is incorporated into R2HCA-EMOA to enhance the performance of the algorithm. Through experimental studies, R2HCA-EMOA is compared with three hypervolume-based EMOAs and several other state-of-the-art EMOAs on 5-, 10-, and 15-objective DTLZ, WFG problems, and their minus versions. Our results show that R2HCA-EMOA is more efficient than the other hypervolume-based EMOAs, and is superior to all the compared state-of-the-art EMOAs.

EMOAs. The Pareto-based EMOAs (e.g., NSGA-II [3] and SPEA2 [4]) use the Pareto dominance as the main selection criterion in environmental selection. The convergence of the population is achieved by the Pareto-based selection criterion. However, the main challenge of the Pareto-based EMOAs for solving MaOPs is that the selection pressure toward the Pareto front (PF) is weakened in high-dimensional objective spaces (e.g., 10-D objective space). Almost all individuals in the population are nondominated with each other in high-dimensional objective spaces, which makes the Pareto-based selection inefficient. The decomposition-based EMOAs (e.g., MOGLS [5] and MOEA/D [6]) decompose an MOP into a set of singleobjective subproblems. Then the subproblems are optimized simultaneously in a collaborative manner. The convergence and diversity of the population are achieved by the scalarizing function criterion. However, for those decomposition-based EMOAs using a fixed weight vector grid, usually the obtained solutions cannot be uniformly distributed on the PF when the shape of the PF is irregular (e.g., the inverted triangular PF [7]). Thus, adaptively adjusting weight vectors in decomposition-based EMOAs is a hot research topic to deal with this issue [8], [9]. The indicator-based EMOAs (e.g., IBEA [10] and SMS-EMOA [11], [12]) use performance indicators (e.g., hypervolume [13], inverted generational distance (IGD) [14], R2 indicator [15]) as the main selection criterion. By properly setting the parameters of the indicators (e.g., the reference point for the hypervolume indicator and the reference point set for the IGD indicator) and optimizing the indicator values of the population, a solution set with good convergence and diversity can be obtained by the indicator-based EMOAs.
Among the indicator-based EMOAs, the hypervolume-based EMOAs (e.g., SMS-EMOA, FV-MOEA [16], and HypE [17]) adopt the hypervolume indicator within the algorithms. The general idea of the hypervolume-based EMOAs is to transform an MOP into a single-objective optimization problem, where the single-objective is to maximize the hypervolume value of the population. The hypervolume indicator is widely used for the performance evaluation of EMOAs. The advantages of the hypervolume-based EMOAs over the other EMOAs are twofold. One is that the hypervolume indicator is strictly Pareto compliant. Thus, the hypervolume-based EMOAs are sensitive to any improvement to a set with respect to Pareto dominance. The other is that the hypervolume-based EMOAs can directly optimize the hypervolume of the population. Thus, good performance comparison results are likely to be obtained by the hypervolume-based EMOAs.
The key procedure in the hypervolume-based EMOAs is the calculation of the hypervolume contribution of each individual. It has been proven that the exact calculation of the hypervolume and hypervolume contributions is NP-hard [18], [19]. Even though there exist O(N log N) algorithms to compute the hypervolume and hypervolume contributions in two-and 3-D objective spaces [20], [21], where N is the population size, the computational time of the hypervolume and hypervolume contributions increases exponentially with the increase in the number of objectives. This limited the applicability of the hypervolume-based EMOAs for solving MaOPs with many objectives (e.g., more than ten objectives). In order to overcome this drawback, the use of hypervolume approximation (instead of exact calculation) has been proposed for hypervolume-based EMOAs. HypE is a representative hypervolume-based EMOA which uses a Monte Carlo sampling method for hypervolume approximation. However, in order to achieve a good hypervolume approximation, a large number of sampling points are needed in HypE. Thus, it is still a time-consuming algorithm.
The above discussions motivate us to develop a new EMOA that inherits the advantage of the hypervolume-based EMOAs and at the same time is computationally efficient for solving MaOPs. In order to fulfil this goal, an efficient hypervolume approximation method is needed. Currently, there are mainly two methods to approximate the hypervolume. One is the Monte Carlo sampling method that is adopted in HypE. The other method is the achievement scalarizing function method proposed in [22] and [23]. The basic idea of this method is to use the average length of a set of line segments from the reference point to the hypervolume attainment surface to approximate the hypervolume. Recently, a new R2 indicator is proposed to approximate the hypervolume [24], which is an extension of the achievement scalarizing function method. Based on [24], an R2-based hypervolume contribution approximation method is proposed in [25]. The proposed method shows its superiority over the Monte Carlo sampling method [26] for the hypervolume contribution approximation in terms of both the runtime and the approximation accuracy.
In this article, a new hypervolume-based evolutionary algorithm, namely, R2HCA-EMOA, is proposed for manyobjective optimization. The proposed algorithm, which is based on the SMS-EMOA framework, employs the R2-based hypervolume contribution approximation method in [25] in environmental selection. Our algorithm has the following two advantages: 1) it inherits the advantage of the hypervolumebased EMOAs, i.e., it can directly optimize the hypervolume of the population and 2) it has a worst-case time complexity of O(N 2 | |) in one generation where N is the population size and | | is the number of the direction vectors. Since the worstcase time complexity does not increase exponentially with the number of objectives, it is suitable for solving MaOPs. The main contributions of this article are listed as follows.
1) Using the R2-based hypervolume contribution approximation method in [25], a new hypervolume-based EMOA, namely, R2HCA-EMOA, is proposed for manyobjective optimization. A normalization mechanism is incorporated into R2HCA-EMOA to enhance its performance. 2) In order to make R2HCA-EMOA computationally efficient, a utility tensor structure is introduced for the calculation of the R2 indicator variant. By efficiently updating the utility tensor, we can avoid the repeated calculation and ensure the efficiency of the algorithm. 3) R2HCA-EMOA is compared with three hypervolumebased EMOAs and several other state-of-the-art EMOAs on 5-, 10-, and 15-objective DTLZ and WFG problems and their minus versions. The results demonstrate the efficiency and effectiveness of R2HCA-EMOA for solving MaOPs. The remainder of this article is organized as follows. In Section II, we briefly review three representative hypervolumebased EMOAs. The proposed algorithm R2HCA-EMOA is elaborated in Section III. Experiments are given in Section IV. We conclude this article in Section V.

A. SMS-EMOA
SMS-EMOA [11], [12] is a classical hypervolume-based EMOA. In SMS-EMOA, the population evolves in a steadystate manner. In each generation, a single offspring is generated and added to the population, then the nondominated sorting is applied to the updated population, and the individual with the smallest hypervolume contribution among the last front individuals is removed from the population. In this manner, SMS-EMOA ensures that the hypervolume of the whole population monotonically increases as the number of generations increases.
SMS-EMOA already showed its promising performance on MOPs and MaOP with 4-6 objectives [27]. However, the main drawback of SMS-EMOA is that it can hardly be applied to MaOPs with many objectives (e.g., more than ten objectives) because the exact hypervolume calculation is time-consuming in high-dimensional spaces. This drawback has been addressed in the following two ways: one is to apply a faster hypervolume calculation method, and the other way is to apply a hypervolume approximation method. These two research directions lead to the following two algorithms: 1) FV-MOEA [16] and 2) HypE [17].

B. FV-MOEA
FV-MOEA [16] is a recently proposed hypervolume-based EMOA. Differs from SMS-EMOA, FV-MOEA has the framework of a (μ + μ ) evolution strategy. In each generation, a set of offspring is generated and added to the population. Next, nondominated sorting-based environmental selection is performed to identify the last front to be included in the next generation. Then individuals in the last front are selected based on a fast hypervolume update procedure. The general idea of this procedure is to delete the individual with the smallest hypervolume contribution one by one. The hypervolume contribution of each of the remaining individuals in this front is updated once an individual is removed.
As reported in [16], FV-MOEA is much faster than SMS-EMOA on up to five-objective MaOPs. This is due to the following two reasons.
1) FV-MOEA adopts a (μ + μ ) evolution strategy while SMS-EMOA adopts a (μ + 1) evolution strategy. It has been shown that a (μ + 1) EMOA takes much longer runtime than its (μ + μ ) counterpart when they are compared under the same number of function evaluations (FEs) [28]. For this reason, FV-MOEA is more efficient than SMS-EMOA. 2) FV-MOEA has a fast hypervolume contribution update procedure. Once a solution is removed from the population, the hypervolume contribution of each individual in the remaining population is efficiently updated instead of recalculated. In this manner, FV-MOEA can achieve a fast computational performance. However, there is no such an update procedure in SMS-EMOA. FV-MOEA shows a good performance on two-and threeobjective MOPs and up to five-objective MaOPs. However, because of the NP-hardness of the exact hypervolume calculation [18], FV-MOEA becomes time-consuming in much higher dimensional spaces (e.g., 10-D space), which limits the applicability of FV-MOEA for solving MaOPs with ten or more objectives.
To overcome the limitation of the exact hypervolume calculation, while keeping the advantage of the hypervolume-based EMOAs, hypervolume approximation can be applied instead of the exact hypervolume calculation.

C. HypE
HypE [17] is a hypervolume estimation-based EMOA with a (μ + μ ) evolution strategy. The framework of HypE is similar to that of FV-MOEA. After a set of offspring is generated and added to the population, the nondominated sorting is performed to identify the last front to be included. Instead of using the exact hypervolume calculation for selecting individuals from the last front as in FV-MOEA, HypE uses a Monte Carlo sampling method to estimate the expected hypervolume loss (i.e., fitness value) attributed to each individual in the last front to be included. The individual with the smallest fitness value is removed one by one, and the fitness value of each of the remaining individuals is updated once an individual is removed.
HypE has been tested on up to 50-objective MaOPs and showed its effectiveness for solving MaOPs. However, HypE is still very time-consuming if a large number of sampling points are used in order to achieve a good hypervolume approximation.

A. R2-Based Hypervolume Contribution Approximation
In our previous work [25], we defined the following R2 indicator variant to approximate the hypervolume contribution of a solution s to a solution set A in an m-dimensional objective where A is a nondominated solution set and s ∈ A, is a direction vector set and each direction vector λ = (λ 1 , λ 2 , . . . , λ m ) ∈ satisfies λ 2 = 1 and λ i ≥ 0, i = 1, . . . , m, r is the reference point, α is a parameter that is chosen to be α = m as suggested in [25]. For convenience, we use R HCA 2 (s) to denote R HCA 2 (s, A, , r, α) in the remaining of this article.
The g *2tch (a|λ, s) function is defined for minimization problems as For maximization problems, it is defined as The g mtch function is defined for both minimization and maximization problems as The R2 indicator variant in (1) has a clear geometric meaning to approximate the hypervolume contribution which is illustrated in Fig. 1(a). Suppose we have three nondominated solutions a 1 , a 2 , and a 3 , a set of direction vectors , and the reference point r. If we draw different line segments with different directions starting from a 2 and ending on the boundary of the hypervolume contribution region of a 2 as shown in Fig. 1 where l i is the length of the line segment associated with the ith direction vector in .
A traditional method for the hypervolume contribution approximation that is based on the achievement scalarizing function method in [22] is illustrated in Fig. 1 [P, T] = EnvironmentalSelection(P, , T, r); 6: FEs = FEs + 1; 7: end while hypervolume of a solution set is approximated by the line segments starting from the reference point r and ending on the hypervolume attainment surface. As shown in Fig. 1(b), the hypervolume contribution of a 2 is approximated by the red lines in the hypervolume contribution region of a 2 . Comparing with the method in Fig. 1(a), the traditional method is less accurate if the same finite number (especially a small number) of direction vectors are used in both methods. This is because R HCA 2 can utilize all direction vectors whereas the traditional method can utilize only a part of them that go through the hypervolume contribution region. The approximation accuracy of these two methods are compared in [25], and the results clearly show that R HCA 2 is much more accurate than the traditional method.
The R2 indicator variant defined in (1) is the foundation for the development of the new algorithm. For more detailed explanations of this R2 indicator variant, please refer to [25]. Based on this R2 indicator variant, we will introduce R2HCA-EMOA in the following sections.

B. General Framework
The general framework of R2HCA-EMOA follows a steadystate (μ + 1) evolution strategy, which is similar to SMS-EMOA. Algorithm 1 gives the pseudocode of R2HCA-EMOA. In each generation, an offspring q is generated from the population P (line 3). Then the population P is updated by incorporating offspring q (line 4). After that, an environmental selection procedure is performed to remove one individual from the population P to keep the population size constant (line 5). Lastly, the current FEs is incremented by 1 since only one offspring is produced in each generation (line 6).
In GenerateOffspring procedure, the offspring is generated as follows. Two parents are randomly selected from population P, and the offspring is produced based on its two parents by applying simulated binary crossover (SBX) [29] and polynomial mutation [30] operators, which are commonly used in the EMO community. With respect to other procedures in R2HCA-EMOA, detailed descriptions are presented in the following sections.

C. Reference Point Specification in R2HCA-EMOA
In fact, R2HCA-EMOA is a hypervolume-based EMOA where its objective is to maximize the hypervolume of the final population. As suggested in [31], the specification of the reference point for the hypervolume indicator should be carefully treated. The specification of the reference point has a large effect on the optimal distribution of the final population on the PF in hypervolume-based EMOAs (e.g., SMS-EMOA and FV-MOEA). If the reference point is too close to the PF, then the solutions cannot widely distribute on the PF. However, if the reference point is too far away from the PF, then the solutions tend to distribute on the boundaries of the PF when the shape of the PF is inverted triangular. In order to obtain widely and evenly distributed solutions on the PF, the reference point r = (r, r, . . . , r) can be specified in the normalized objective space with the ideal point (0, 0, . . . , 0) and the nadir point (1, 1, . . . , 1) as follows as suggested in [31]: where H is an integer satisfying C H+m−1 m−1 ≤ N < C H+m m−1 , and C n m is the total number of combinations for choosing m elements from a set of n elements, i.e., C n . In each generation of R2HCA-EMOA, the current population P will be normalized to the objective space (0, 0, . . . , 0) to (1, 1, . . . , 1), and the reference point r is specified according to (5).

D. Environmental Selection
Algorithm 2 gives the pseudocode of the environmental selection procedure. First, the population P is normalized to P (line 1), then the utility tensor T is updated based on P (line 2). After that, the nondominated sorting is performed to divide population P into different front levels (line 3). Finally, one individual with the smallest R HCA 2 value from the last front F l is removed from P (lines 4-13). If there is only one individual in F l , then this individual can be directly removed without calculating its R HCA 2 value (lines 4 and 5). If there is more than one individual in F l , then the R HCA 2 values of all individuals in F l need to be calculated first (lines 6-10).
In Algorithm 2, a utility tensor T is used to assist the calculation of the R HCA 2 values. In the next section, we will explain tensor T and its operations in detail.

1) Definition of Tensor T:
In R2HCA-EMOA, the most time-consuming part is the calculation of the R HCA 2 values. In order to maximize the computational efficiency of this part, a utility tensor T is introduced to assist the calculation of R HCA 2 .
Suppose the current population with the offspring included is P = {s 1 , s 2 , . . . , s N+1 }. T is a 3-order tensor with the following elements: Once the tensor T is obtained, its information can be used to calculate the R HCA 2 values. Suppose P is a nondominated solution set, then for a solution s k ∈ P, its R HCA 2 value can be calculated as follows: Notice that P should be a nondominated solution set. If P contains dominated solutions, then we cannot directly use T to calculate R HCA 2 values by (7). We will discuss how to deal with this situation later in this section.
Next, we give a toy example to illustrate the tensor T defined above, and how it can be used to calculate R HCA 2 .
Example 1: Suppose we have a solution set A = {a 1 , a 2 , a 3 }, a direction vector set = {λ 1 , λ 2 , λ 3 }, and the reference point r. Then the tensor T is depicted in Fig. 2. Suppose A is a nondominated solution set, then the R HCA 2 values of the solutions in A are calculated as follows: 2) Update of Tensor T: In R2HCA-EMOA, only one offspring is generated in each generation. Thus, we do not have to recalculate all elements of the tensor T in each generation. We update only a small number of elements in T for computational efficiency. Once an individual is removed from the population, a new offspring is added. This generation update mechanism can be viewed as the replacement of the removed individual with the new offspring. Thus, all elements in T related to the removed individual are updated with the new offspring (line 2 in Algorithm 2).
From the above equation we can see that totally (N+1)| |+ N| | elements of tensor T need to be updated. We illustrate this update procedure by the following toy example.
Example 2: Let us consider the case in Example 1. Suppose a 2 is removed and q is added, then the tensor T is updated as depicted in Fig. 3 where the updated elements are shaded.
3) Extraction of Tensor T: In R2HCA-EMOA, one individual in the last front F l is removed. If F l is the same as the whole population (i.e., l = 1, which means all individuals in the population are nondominated with each other), then the tensor T can be directly used to calculate the R HCA 2 values of all individuals (line 7 in Algorithm 2). However, F l usually contains only a part of the population (i.e., l > 1). In this case, a subtensor T is extracted from T (line 9 in Algorithm 2) in order to calculate the R HCA 2 values of the individuals in F l . The extraction is straightforward, i.e., elements in T that are related only to individuals in F l (i.e., those that are not related to any individual in the other fronts) are extracted to form the subtensor T .
From the above equation we can see that the size of T is M × | | × M, which means that totally M 2 | | elements need to be extracted from T.
Once the subtensor T is extracted from T, we can use it to calculate the R HCA 2 values of the individuals in F l as follows: We will use the following toy example to illustrate this procedure.
Example 3: Let us consider the case in Example 1 again. Suppose F l = {a 1 , a 3 }, i.e., a 2 dominates a 1 and a 3 while a 1 and a 3 are nondominated with each other. Thus, all elements related to a 1 and a 3 and not related to a 2 are extracted from T to form T as shown in Fig. 4, where the extracted elements in T are shaded. Then the R HCA 2 values of the solutions in F l are calculated as follows:

F. Normalization Mechanism
In R2HCA-EMOA, the normalization mechanism is incorporated before the update of the tensor T (line 1 in Algorithm 2). The normalization works as follows. First, we estimate the ideal point z * and the worst point z wst from the current population P as z * i = min s∈P s i and z wst i = max s∈P s i , i = 1, . . . , m. Then each solution s ∈ P is normalized as . . , m where s ∈ P is the normalized solution. After the normalization, the population P is within the objective space (0, 0, . . . , 0) to (1, 1, . . . , 1). The normalization mechanism is used to tackle the MOPs and MaOPs with a totally different scale in each objective. For such a problem, the hypervolume contribution approximation by R HCA 2 can be inaccurate if no normalization is performed. The normalization mechanism remedies this potential difficulty. However, the normalization will bring an issue for the update of the tensor T. In the last section, we explained that only the elements in the tensor T related to the removed individual are updated with those to the newly added offspring in order to improve the computational efficiency. If the normalization is performed, the locations of all solutions in the objective space can be changed. If this happens, we need to update all elements instead of updating only a small number of elements in the tensor T. Actually, this can happen frequently in the early generations of the optimization process. Thus, the tensor T needs to be recalculated frequently which leads to low computational efficiency.
In this article, in order to incorporate the normalization mechanism in R2HCA-EMOA while keeping the algorithm computationally efficient, we simply ignore the normalization in the update phase of the tensor T. That is, we only replace the elements in T related to the removed individual with those to the new offspring (even if other solutions are changed in the objective space after the normalization). In this manner, the computational efficiency of the proposed R2HCA-EMOA is not degraded. This simple approach is based on the following considerations. In early generations, accurate approximation of hypervolume contributions is not important since the search is mainly driven by nondominated sorting and/or the differences of the hypervolume contributions among the individuals are large. In late generations, the update of the estimated ideal and worst points does not happen frequently. In Section IV-D1, we will show the effectiveness of the normalization mechanism together with the simple update strategy in R2HCA-EMOA.

G. Computational Complexity
The computational complexity of one generation of R2HCA-EMOA is analyzed as follows. The main computational cost is from the environmental selection procedure (line 5 in Algorithm 1). The normalization of population P (line 1 in Algorithm 2) requires O(mN) computations. In the tensor update procedure, totally (N + 1)| | + N| | elements of tensor T need to be updated. Therefore, update of tensor T If the tensor T is not used in R2HCA-EMOA, that is, the R HCA 2 values of the individuals in F l are recalculated in each generation, it requires O(mN 2 | |) computations in the worst case. So the overall worst-case time complexity of R2HCA-EMOA will be O(mN 2 | |). By introducing the tensor T and its operations, we can improve the time complexity of R2HCA-EMOA by a factor of m in each generation.

IV. EXPERIMENTAL STUDIES
A. Experiment Settings 1) Platforms: The software platform we use in our experiments is PlatEMO [33], which is a MATLAB-based open source platform for EMO. The R2HCA-EMOA is implemented in PlatEMO. All the other EMOAs used in our experiments are based on their implementations in PlatEMO. 2 The hardware platform is a PC equipped with Intel Core i7-8700K CPU@3.70 GHz, 16-GB RAM.
3) Parameter Settings: The population size N is set to 100 for all test problems. For DTLZ1, DTLZ3, WFG1, and their minus versions, the maximum function evaluations FEs max is set to 100 000, while FEs max is set to 30 000 for the other test problems. These values of FEs max are used as the termination condition. SBX and polynomial mutation are used where the distribution index is specified as 20 in both operators. The crossover and mutation rates are set to 1.0 and 1/n, respectively, where n is the number of decision variables. Each algorithm is applied to each test problem 20 times (i.e., 20 independent runs).

4) Performance Metrics:
The hypervolume indicator is employed to evaluate the performance of the compared EMOAs. For 5-and 10-objective problems, the WFG algorithm [36] is employed to calculate the exact hypervolume of the obtained solution set in each run of each algorithm on each test problem, whereas for 15-objective problems, the Monte Carlo method is employed to estimate the hypervolume. We calculate the hypervolume as follows. First, using the true ideal point z * and the true nadir point z nad obtained from the true PF, the objective values of the obtained solutions are normalized so that two points z * and z nad are (0, . . . , 0) and (1, . . . , 1), respectively. Then the reference point r is set to (1.1, . . . , 1.1) to calculate the hypervolume. Furthermore, the results are analyzed by the Wilcoxon rank sum test with a significance level of 0.05 to determine whether one algorithm shows a statistically significant difference with the other, where "+," "−," and "≈" indicate that the compared EMOA is "significantly better than," "significantly worse than," and "statistically similar to" R2HCA-EMOA, respectively. 5) Compared Algorithms: First, three hypervolume-based EMOAs (i.e., SMS-EMOA, FV-MOEA, and HypE) are compared with R2HCA-EMOA. Both SMS-EMOA and FV-MOEA employ the WFG algorithm for the hypervolume calculation inside their implementations. In HypE, the number of sampling points is set to 10 000, which is the same as in [17]. The reference point is specified as 1 + 1/H in the normalized objective space in the three hypervolume-based EMOAs, which follows the same manner as in R2HCA-EMOA.
We also compare R2HCA-EMOA with several other stateof-the-art EMOAs. They are briefly introduced as follows.
1) AR-MOEA [37] is an indicator-based EMOA, which is based on an enhanced inverted generational distance indicator (IGD-NS). The reference points are adaptively adjusted based on the IGD-NS contributions of the candidate solutions in an external archive, thus taking both uniform distribution and PF approximation into consideration. 2) SPEA2+SDE [38] is a combination of SPEA2 and the shift-based density estimation (SDE) strategy. SDE modifies the density estimation strategy in traditional Pareto-based EMOAs and is able to cover both the distribution and convergence information of individuals. Thus, it can make Pareto-based algorithms suitable for MaOPs. 3) GFM-MOEA [39] is an EMOA based on PF modeling.
The shape of the PF is estimated by training a generalized simplex model, then the algorithm is driven by this approximated PF during the optimization process. The penalty parameter θ and the frequency f r in GFM-MOEA are set to 0.2 and 0.1, respectively, which are the same as in [39]. 4) Bi-criterion evolution (BCE)-IBEA [40] is an EMOA with a BCE framework. BCE-IBEA integrates IBEA [10] into the BCE framework, where the Pareto criterion evolution and the non-Pareto criterion evolution (i.e., IBEA) work collaboratively to facilitate each other's evolution. IBEA is based on I + indicator [41] and the scaling factor κ in IBEA is set to 0.05. The parameters involved in R2HCA-EMOA are set as follows. The number of the direction vectors | | is set to 100. The direction vectors used in R HCA 2 are generated by sampling points uniformly on the unit hypersphere [42]. 3 For all the EMOAs involved in our experiments, SMS-EMOA and R2HCA-EMOA follow the (μ + 1) evolution strategy, while the other EMOAs follow the (μ + μ ) evolution strategy where μ = μ (i.e., N offspring are generated from N parents in each generation).
First, we compare the runtime among the four algorithms to evaluate their computational efficiency. As SMS-EMOA and FV-MOEA become very time-consuming for solving 10-and 15-objective problems, we only test the four algorithms on DTLZ2 with 5, 10, and 15 objectives by a single run, and set FEs max to 3000 (i.e., 1/10 of the original FEs max setting for DTLZ2). Fig. 5 shows the runtime comparison of the four algorithms in a single run on DTLZ2 with 5, 10, and 15 objectives. From Fig. 5 we can observe that, for the five-objective DTLZ2, FV-MOEA is the fastest algorithm. The runtime of R2HCA-EMOA is comparable to that of FV-MOEA. SMS-EMOA and HypE consume much longer time than the other two algorithms. However, for the 10-and 15-objective DTLZ2, the runtime of SMS-EMOA and FV-MOEA increase exponentially while the runtime of HypE and R2HCA-EMOA increase linearly. It should be noted that SMS-EMOA and FV-MOEA take several hours for the 15-objective DTLZ2 even when we decrease the termination condition from 30 000 FEs to 3000. R2HCA-EMOA is the fastest algorithm for the 10-and 15objective DTLZ2. Because 10 000 sampling points are used in HypE while only 100 direction vectors are used in R2HCA-EMOA, HypE is much slower than R2HCA-EMOA in all the cases. From Fig. 5 we can also observe that FV-MOEA is faster than SMS-EMOA for the five-objective DTLZ2, while it is slower than SMS-EMOA for the 10-and 15-objective DTLZ2. The main reason is that with the offspring included in each generation, the population size of SMS-EMOA and FV-MOEA is N+1 and 2N, respectively. The almost 2 times larger population size of FV-MOEA will lead to much longer computational time of the hypervolume in 10-and 15-objective space. The observation suggests that FV-MOEA is more efficient than SMS-EMOA for five-objective problems which is consistent with the results in [16]. However, FV-MOEA is less efficient than SMS-EMOA for 10-and 15-objective problems due to the exponentially increasing time in the hypervolume calculation.
On account of the very long runtime of SMS-EMOA and FV-MOEA on 10-and 15-objective problems, we will only compare the four algorithms on 5-objective problems. The results are shown in Table I. In Table I, we only show the Wilcoxon rank sum test results. Detailed numerical results are provided in Appendix A in the supplementary material.
The results show that SMS-EMOA and FV-MOEA outperform HypE and R2HCA-EMOA, which is consistent with our intuition. Since SMS-EMOA and FV-MOEA are based on exact hypervolume calculation, it is likely that they can achieve better hypervolume results than HypE and R2HCA-EMOA with hypervolume approximation. However, SMS-EMOA and FV-MOEA become very time-consuming for solving 10-and 15-objective problems, which limited their applicability for solving MaOPs with more than ten objectives.
We can also observe that R2HCA-EMOA clearly outperforms HypE (+/−/≈ is 0/9/4 for DTLZ & WFG and 0/13/0 for MinusDTLZ & MinusWFG as shown in Table I). HypE employs the Monte Carlo sampling method while R2HCA-EMOA employs the R2-based hypervolume contribution approximation method. Notice that 10 000 sampling points are used in HypE while only 100 direction vectors are used in R2HCA-EMOA. Comparing with HypE, R2HCA-EMOA consumes much less runtime and is able to achieve better hypervolume results, which shows that R2HCA-EMOA is an From the figures we can observe that in general, the mean hypervolume difference decreases with the increase of FEs for all the four algorithms (except for HypE on WFG3 where the hypervolume difference increases with the increase of FEs). For most of the problems, R2HCA-EMOA has a competitive performance compared with SMS-EMOA and FV-MOEA, and clearly outperforms HypE. On some problems (e.g., DTLZ1, DTLZ3, MinusDTLZ1, and MinusWFG3), R2HCA-EMOA performs even better than FV-MOEA. Meanwhile, R2HCA-EMOA has a relatively small standard deviation on most of the problems (except for MinusDTLZ1), which shows that R2HCA-EMOA has a stable performance in general. Fig. 6 provides us a deeper insight into the behavior of R2HCA-EMOA, which clearly shows the ability of R2HCA-EMOA for achieving a good hypervolume performance.
The results are shown in Table II Table II). This is because only R2HCA-EMOA is a hypervolume-based EMOA whereas the other algorithms are not. We need to mention that in terms of the hypervolume metric, AR-MOEA outperforms MOEA/DD [43], NSGA-III [44], RVEA [45], and MOMBI-II [46] as reported in [37], and GFM-MOEA outperforms MOEA/DD, RVEA, MOEA/D-PaS [47], and VaEA [48] as reported in [39]. However, these two algorithms are both defeated by R2HCA-EMOA. Our results clearly demonstrate the effectiveness of R2HCA-EMOA for solving MaOPs in terms of the hypervolume metric.
We also show the hypervolume variations over FEs obtained by the five EMOAs on all the 10-objective problems in Fig. 7 (The figures related to 5-and 15-objective problems are provided in Appendix B in the supplementary material). Same as in Fig. 6, the hypervolume difference is shown in Fig. 7. For each figure, the y-axis is shown in a linear scale if the lines However, as will be seen in Section IV-E, R2HCA-EMOA outperforms the other algorithms on all minus problems based on a different reference point specification for hypervolume calculation. We can also observe that R2HCA-EMOA has a stable performance on most of the problems. The results in Fig. 7 clearly show the superiority of R2HCA-EMOA over the other algorithms for achieving high average hypervolume values.
D. Further Investigations on R2HCA-EMOA 1) Effect of the Normalization Mechanism: R2HCA-EMOA has introduced the normalization mechanism to enhance the performance of the algorithm. In order to investigate the effectiveness of this mechanism, we compare R2HCA-EMOA with and without the normalization mechanism.
The results are shown in Table III for five-objective DTLZ and WFG problems. From Table III we can see that R2HCA-EMOA with normalization clearly outperforms that without normalization (+/−/≈ is 1/5/7 as shown in Table III), which means that the effect of the normalization mechanism is significant on the performance of the algorithm.
We can also observe that R2HCA-EMOA with normalization clearly achieves better results than that without normalization on almost all problems (except for WFG3). Moreover, R2HCA-EMOA with normalization is significantly better than that without normalization on many WFG problems which are MaOPs with a different scale in each objective. As discussed in Section III-F, the normalization mechanism is introduced to tackle this type of problems. We also use a simple update strategy to update the tensor T in order to keep the algorithm computationally efficient. The results in Table III demonstrate the effectiveness of the normalization mechanism together with the simple update strategy in R2HCA-EMOA.   Table IV. From Table IV we can see that the results are consistent with our intuition. R2HCA-EMOA with 200 direction vectors outperforms that with 100 direction vectors, and R2HCA-EMOA with 100 direction vectors outperforms that with 50 direction vectors. However, a larger number of the direction vectors means longer runtime. There is a tradeoff between the runtime and the performance of R2HCA-EMOA. The principle to set the number of the direction vectors should depend on the computing power, available time, performance requirement, etc.
3) Influence of the Population Size: In our experiments, we set the population size N to 100 for all test problems. In order to investigate the influence of the population size on the performance of R2HCA-EMOA, we also perform our computational experiments with N = 50 and N = 200. R2HCA-EMOA is compared with AR-MOEA, SPEA2+SDE, GFM-MOEA, and BCE-IBEA on the 5-objective DTLZ and WFG problems. The other settings are the same as in Section IV-A. The results are shown in Table V. Detailed numerical results are provided in Appendix C in the supplementary material.
From Table V we can see that R2HCA-EMOA outperforms all the compared algorithms in terms of the hypervolume metric in both cases: N = 50 and N = 200. This observation further demonstrate the effectiveness of R2HCA-EMOA for solving MaOPs under different population sizes.

E. Additional Results
In our experiments, the hypervolume indicator is employed to evaluate the performance of R2HCA-EMOA. Since R2HCA-EMOA is a hypervolume-based EMOA, the hypervolume indicator is able to directly evaluate the ability of R2HCA-EMOA for maximizing the hypervolume value of the population. In addition to the hypervolume indicator, the IGD indicator [14] is another widely used performance indicator for the performance evaluation. We provide the IGD results of all the compared EMOAs in this article in Appendix D in the supplementary material. Interested readers can refer to the detailed discussions there.
In addition, the hypervolume indicator is evaluated based on the reference point specification r = 1.1 as explained in Section IV-A4. This reference point specification is commonly used in the EMO community for the hypervolume performance evaluation of EMOAs. However, different reference point specifications may affect the hypervolume evaluation results [31]. Thus, in order to fairly compare the hypervolume performance of different EMOAs, we provide the hypervolume results of all the compared EMOAs based on the reference point specification r = 1 + 1/H in Appendix E in the supplementary material. Detailed discussions can be found there.
Finally, in order to visually explain the behavior of R2HCA-EMOA, we show the final solution distributions on some 3-objective test problems obtained by R2HCA-EMOA and the other three hypervolume-based EMOAs. The results are provided in Appendix F in the supplementary material.

V. CONCLUSION
In this article, we proposed R2HCA-EMOA, a new hypervolume-based EMOA for many-objective optimization. By using the R2-based hypervolume contribution approximation method in environmental selection, R2HCA-EMOA is able to inherit the advantage of the hypervolume-based EMOAs. At the same time, it is efficient for solving MaOPs. We demonstrated the efficiency and the effectiveness of the proposed algorithm through comparative experiments. First, we showed that R2HCA-EMOA is slightly outperformed by SMS-EMOA and FV-MOEA, both of which need much longer computation time than R2HCA-EMOA. Comparing with HypE, R2HCA-EMOA is superior in terms of both the runtime and the hypervolume metric. Finally, we compared R2HCA-EMOA with several state-of-the-art EMOAs. Our results showed that R2HCA-EMOA outperforms all the compared algorithms in terms of the hypervolume metric.
As future research, we have the following three research directions: 1) it has been shown that a dynamic reference point specification mechanism can enhance the performance of the hypervolume-based EMOAs [49]. In the future, we will investigate the effect of the dynamic reference point mechanism on the performance of the proposed algorithm; 2) it is interesting to test the proposed algorithm on real-world problems; and 3) it is also interesting to extend the proposed algorithm for solving large-scale MOPs [50] and multimodal MOPs [51].