A Cross-Reference Line Method Based Multiobjective Evolutionary Algorithm to Enhance Population Diversity

Abstract

Multiobjective evolutionary algorithms (MOEAs) with higher population diversity have been extensively presented in literature studies and shown great potential in the approximate Pareto front (PF). Especially, in the recent development of MOEAs, the reference line method is increasingly favored due to its diversity enhancement nature and auxiliary selection mechanism based on the uniformly distributed reference line. However, the existing reference line method ignores the nadir point and consequently causes the Pareto incompatibility problem, which makes the algorithm convergence worse. To address this issue, a multiobjective evolutionary algorithm based on the adaptive cross-reference line method, called MOEA-CRL, is proposed under the framework of the indicator-based MOEAs. Based on the dominant penalty distance (DPD) indicator, the cross-reference line method can not only solve the Pareto incompatibility problem but also enhance the population diversity on the convex PF and improve the performances of MOEA-CRL for irregular PF. In addition, the MOEA-CRL adjusts the distribution of the cross-reference lines directly defined by the DPD indicator according to the contributing solutions. Therefore, the adaptation of cross-reference lines will not be affected by the population size and the uniform distribution of cross-reference lines can be maintained. The MOEA-CRL is examined and compared with other MOEAs on several benchmark problems. The experimental results show that the MOEA-CRL is superior to several advanced MOEAs, especially on the convex PF. The MOEA-CRL exhibits the flexibility in population size setting and the great versatility in various multiobjective optimization problems (MOPs) and many-objective optimization problems (MaOPs).

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Multiobjective optimization problems (MOPs), which have more than one conflicting objective to be optimized, can be defined as$$\begin{array}{}\text{(1)}& \begin{array}{cc}\text{minimize}& F\left(x\right)=\left\{f_1\left(x\right),f_2\left(x\right),\dots ,f_m\left(x\right)\right\},\\ \text{subject\hspace{0.17em}to}& x\in \Omega ,\end{array}\end{array}$$where F ⊂ R^m is the objective vector, m is the number of objectives, R^m is the objective space, x = (x₁, …, x_n)^T ∈ Rⁿ is the candidate solution, and $\mathrm{\Omega }={\displaystyle {\prod }_{i=1}^n\left[a_i,b_i\right]}\subseteq R^n$ is the n-dimensional feasible search space. $F:\mathrm{\Omega }⟶R^m$ defines m real-valued objective functions and indicates a mapping which is from the feasible search space to the objective space. Let x₁, x₂ ∈ Ω come to two solutions in the feasible search space. Next, if and only if f_i (x₁) ≤ f_i (x₂) for each i = 1, 2, …, m and f_i (x₁) ≠ f_i (x₂) for ∃ j ∈ {1, 2, …, m}, this shows that x₁ dominates x₂. If there is no any x ⊂ Ω which makes F (x) dominates F (x^∗), x^∗ will be called a global Pareto-optimal solution. The number of Pareto-optimal solutions is generally more than one in a multiobjective optimization problem, and the set of the Pareto-optimal solutions is named as the Pareto-optimal set. The Pareto-optimal set reflects the geometry of the Pareto front (PF) [1].

In recent years, many multiobjective evolutionary algorithms (MOEAs) have been proposed for solving multiobjective optimization problems in various fields [2–5]. In general, MOEAs have two main tasks: (1) ensuring the convergence pressure to drive populations to PF and (2) enhancing the diversity to spread populations evenly to PF [6]. If the two tasks can be accomplished well together, the approximation of PF will be better.

Although the existing MOEAs have been proved to be effective in many practical applications [7], they still have troubles in some complex multiobjective optimizations. The imbalance of convergence and diversity are the major issues since most MOEAs are designed on the principle of “convergence first and diversity second.” For instance, SPEA2 [8] and NSGA-II [9] reach the convergence based on nondominated relations, and the niche technology is used to decrease the crowding population for higher diversity. Meanwhile, the MOEAs based on decomposition such as MOEA/D [10] also choose the nondominated solutions firstly, instead of protecting the solutions with higher diversity in priority.

The imbalance of convergence and diversity makes it difficult to provide searches on different levels in the objective space. Most parts of the objective space are easier to be searched than the rest. If the nondominated solutions are not uniformly distributed, the candidate solutions may be remained far away from each other and the population may be in danger of losing diversity. In particular, it is difficult to generate feasible solutions in the untapped objective space by genetic factors since the remote parents cannot generate good offspring solutions efficiently in multiobjective optimization [11, 12]. In fact, some dominant feasible solutions can enhance the population diversity and use them appropriately can increase the pressure of choice in high-dimensional MOPs. In this sense, the diversity is as important as the convergence and should be emphasized in multiobjective optimization. Therefore, the studies on the method of diversity enhancement have been conducted to ensure a good approximation of PF. Some methods of diversity enhancement are used in the following three main types of EMO algorithms.

The first type of MOEAs is established on the fundamentals of Pareto dominance theory. The MOEAs selecting the population of next generation by Pareto dominance theory prefer nondominated individuals. The Pareto dominance theory itself does not promote the preservation of individuals with the diversity in the objective space, so the crowded strategy, niche theory, and other auxiliary strategies for diversity enhancement are proposed in order to gain the expansion of the objective space and enhance diversity. In NSGA-II, the diversity is improved by the crowded distance [9], and the crowded strategy has been extended to multiple strategies [13]. Corne et al. [14] proposed PESA-II based on regional selection strategies to improve diversity. Horn et al. [15] proposed NPGA based on dynamic niche strategy to enhance diversity.

The second type of MOEAs is built on the concept of decomposition. In MOEA/D [16], each candidate solution is linked with a subproblem and every subproblem is optimized by the information from its neighbors. Most of the decomposition-based MOEAs enhance the diversity by a uniform distribution of weight. There are two main forms of decomposition-based MOEAs. One is to decompose the initial MOPs into a series of single-objective optimization problems (SOPs), such as MOGLS [17], CMOGA [18], MSOPS-II [19], MOEA/D [16], and RVEA [20]. The other is to transform the initial MOPs into multiple more simplified MOPs through the way of dividing the objective space into multiple subspaces, such as MOEA/D-M2M [21], IM-MOEA [10], NSGA-III [11], and SPEA/R [22]. It is worth noting that the uniformity of weight distribution can enhance the diversity when a uniformly distributed weighting strategy is adopted, but the strategy may be failed in some special PFs. Therefore, the method of adaptive weights have been proposed in [23] to improve diversity in some special PFs. The method of adaptive weights have also been used in NSGA-III [11] to enhance the diversity.

The third category is known as the indicator-based MOEAs. The indicator-based theory guides the selection process by integrating the convergence and diversity into a single indicator. The advantage comparison method proposed by Sun et al. [24] enhances the diversity of evolutionary algorithms based on an inverted generational distance (IGD) indicator with the rank value and selection mechanism. Liu et al. [25] proposed a comparison algorithm that enhances the diversity of evolutionary algorithms based on the GD indicator. The comparison algorithm is not affected by the comparison order of individuals, so those solutions with good diversity can get more attention. Recently, AVREA [26] has adopted the adaptive reference line method and the achievement scalarizing function (ASF) as a/the secondary selection indicator, which efficiently enhance the diversity of AVREA.

In addition, many reference line methods have been widely used to enhance the diversity of MOEAs in recent years. For instance, MOEAs based on the nondominated sorting method (NSGA-III [11]) use a group of reference lines and choose those solutions being closer to the reference lines to enhance the diversity. Yuan et al. [27] offered an MOEA based on the reference line method that adopted a diversity enhancement mechanism similar to NSGA-III by measuring the distance between the origin and the projection of the candidate solution on the reference line. Furthermore, a method (RVEA) was proposed in [20] to adaptively modify the reference vector position based on the scale of the objective function to balance the diversity and the convergence. Recently, Sun et al. [28] proposed a MOEA based on reference lines. By the reference line method, the poles located on the coordinate axis in the objective space can be detected in order to construct a hyperplane. In particular, the boundary reference line is generated by linking the origin with the reference points on the axis, which improves the distribution of candidate solutions in the objective space. The boundary reference line is combined with the internal reference line, which divides the objective space into multiple subspaces. Subsequently, the selection for candidate solutions by reference lines enhances the diversity of the proposed algorithm. Due to the nature of the reference line method for the diversity enhancement and the auxiliary selection mechanism based on the uniformly distributed reference line, the reference line method is increasingly favored by MOEAs.

In general, the reference line method is constructed by ideal points and reference points. However, the existing reference line methods seldom take into account the influence of the nadir point although the objective space of each generation is often limited between the ideal point and the nadir point. Therefore, the influence of the nadir point cannot be ignored. In addition, when dealing with the problems with convex PFs, the reference line method can cause Pareto incompatibility problems, which makes the convergence of the algorithm worse. Considering the excitations and the defects of reference line, we propose the method based on adaptive cross-reference line. Compared to the existing reference line method, the main new contributions of this work can be summarized as follows:

(1)

The concept of the cross-reference line is proposed, and a new MOEA called MOEA-CRL is proposed. It inherits the advantages of the ideal point reference line in convergence, adding the nadir point reference line to enhance diversity. The ideal point reference line is combined with the nadir point reference line to divide the objective space into multiple subspaces, and the unique contributing solution preserved in each subspace to ensure the uniform distribution of the Pareto solution set. Compared with the existing reference line methods, the proposed cross-reference line method performs better in terms of diversity.

(2)

Based on the cross-reference line, a dominant penalty distance (DPD) is proposed to solve the Pareto incompatibility problem caused by the reference line method. Compared with the existing reference line indicators, the DPD indicator combines the properties of the ideal point reference line and the nadir point reference line, which not only solves the Pareto incompatibility problem but also improves the performance of the MOEA-CRL on the convex PF.

(3)

The cross-reference line adaptation method is proposed to improve the performance of MOEA-CRL for irregular PFs. The cross-reference line adaptation method not only achieves the uniform distribution by uniformly sampling points on the unit hyperplane but also adaptively adjusts the distribution of the cross-reference lines according to the contributing solutions. Compared with the existing reference line adaptive method, it adjusts the distribution of the cross-reference lines according to the contributing solutions directly defined by the DPD indicator, so the adaptation of the cross-reference line can be not affected by the population size, and the uniform distribution of the cross-reference lines can be maintained.

The rest of this paper is organized as follows. In Section 2, the PBI reference line method is analyzed in detail, and the Pareto incompatibility problem is raised. In Section 3, new reference line methods and evaluation indicators are explored to avoid the Pareto problem on the premise of fully considering the complementarity between the nadir point and the ideal point. In Section 4, the details of the proposed MOEA-CRL are mainly described. The empirical results of MOEA-CRL compared with existing MOEAs are given in Section 5. In the end, the conclusions and the future work are set out in Section 6.

2. Related Work

2.1. Reference Method Based on MOEAs

In the MOEAs based on decomposition, the existing reference line method is more advanced compared with the reference point method [29]. Especially, for the convex PF, the reference line method can effectively enhance the diversity of candidate solutions close to the coordinate axes. As the two examples given, Figure 1 shows the objective space of the two-objective optimization problems with the concave PF and the convex PF. Both examples contain five candidate solutions located on the convex PF and the concave PF, and three of them closest to the three reference points are considered as contributing solutions. If a candidate solution is closest to a reference point and the reference point is also closest to a candidate solution, the candidate solution is defined as a contributing solution of the reference point. It should be noted that, as shown in Figure 1(b), two pole-solutions close to the coordinate axis cannot be regarded as the contributing solutions.

[figures omitted; refer to PDF]

In order to improve the evaluation limitations of the reference point method, the reference point-based method has been developed into various reference line methods [29]. The distance between a candidate solution and a reference point is replaced by the distance between a candidate solution and a reference line. Therefore, a candidate solutions close to the coordinate axis can be evaluated fairly by the reference line method. As shown in Figure 2, the reference line is generated based on each reference point and the ideal point Z^∗ separately. The distance between the candidate solution p and the reference line can be expressed as$$\begin{array}{}\text{(2)}& \text{dis}=F\left(p\right)\sin \left(\stackrel{⟶}{Z^{\ast }r},F\left(p\right)\right),\end{array}$$where $F\left(p\right)$ denotes the Euclidean distance between the candidate solution and the ideal point and $\stackrel{⟶}{Z^{\ast }r}$ denotes the vector from the ideal point to the reference point.

[figure omitted; refer to PDF]

2.2. PBI Reference Line Method Based on Ideal Points

The aggregation function is used as a fitness value function that weighs the merits of an individual. The aggregation function is usually a function of the individual x in the objective space under the condition of a given weight. The optimization of each subproblem is regarded as the optimization of the aggregation function. The PBI aggregation function is a variant of the method based on the intersection of boundaries, which aims to find the intersection point between the Pareto front and a set of lines [30]. Studies have shown that PBI aggregation functions with appropriate penalty parameter values can generate more uniform candidate solution sets, but the performance of PBI highly depends on the setting of penalty parameters that control the balance between convergence and diversity [31].

The aggregate optimization equation of the PBI function $g_{\text{pbi}}$ is [32]$$\begin{array}{}\text{(3)}& \text{Minimize}{g}_{\text{pbi}}\left(\left.x\right|\lambda \right)=d1+\theta \text{d}2,\end{array}$$where$$\begin{array}{c}\text{(4)}& d_1=\frac{{\left(f\left(x\right)-z\right)}^T\lambda }{\lambda },d_2=f\left(x\right)-\left(z-d_1\frac{\lambda }{\lambda }\right),\end{array}$$where f (x) is a candidate solution, z is the ideal point, and λ is the vector of the reference line.

PBI also uses the obtained ideal point z as the criterion to decompose the objective space. θ is the parameter of PBI, and its range is θ ≥ 0. Figure 3 shows the d₁ and d₂ of the solution x of a weight vector λ = (0.5, 0.5)^T in the two-dimensional objective space. In the PBI method, a candidate solution with a small d₁ is first considered as a better candidate solution close to the Pareto front. In addition, the distance d₂ from the weight vector λ is considered. Finally, $g_{\text{pbi}}$ is calculated by adding the value of d₂ multiplied by θ to d₁. In summary, a candidate solution with a small d₁ and d₂ is considered as a better candidate solution. The balance between d₁ and d₂ in $g_{\text{pbi}}$ is controlled by the parameter θ. Therefore, the PBI method evolves a candidate solution toward z by minimizing $g_{\text{pbi}}$.

[figure omitted; refer to PDF]

2.3. Pareto Incompatibility

The Pareto incompatibility problem means that the individual's reference line evaluation results may face conflicts with the results of the nondominated relationship during the iteration process. The reference line method of the PBI aggregation function can effectively improve the diversity of candidate solutions near the coordinate axis in convex PF by increasing the value of the parameter θ so that the influence of d₂ is much larger than the influence of d₁. This method not only maintains the diversity of candidate solutions for convex PF close to the coordinate axis, but also quickly obtains the candidate solution p with the smallest $g_{\text{pbi}}$. However, the reference line method using only ideal points causes Pareto incompatibility. Figure 4 shows an example to illustrate the problem of Pareto incompatibility.

[figure omitted; refer to PDF]

It can be clearly seen from Figure 4 that this simple example is a two-objective minimization problem, where the reference point set is {(2, 4), (5, 2)}, and the candidate solution set A = {(3, 5), (8, 4)} and the candidate solution set B = {(3, 6), (10, 4)}. The candidate solution set B = {(3, 6), (10, 4)} is the smallest distance d₂ from the reference line, so the candidate solution set B is the best. However, according to Pareto dominance theory, (3, 5) in solution set A dominates (3, 6) in candidate solution set B, and (8, 4) in solution set A dominates (10, 4), so candidate solution set A is better than candidate solution set B. Therefore, the PBI aggregate function may cause the judgment of the merits of the solution set to be contrary to Pareto dominance theory.

2.4. Complementarity of the Nadir Point and the Ideal Point

The setting of the reference point in the aggregation function plays a key role in the performance of MOEA/D. In fact, different types of reference points may have different effects on the exploration behavior of MOEA/D. Most MOEA/D improvements use ideal points as reference points. As mentioned in [33], when diversity is easy to maintain, the method using only ideal points will be effective, and only using ideal points is more helpful to promote candidate solutions to approximate PF. In [34], MOEA/D uses ideal point and a set of reference points evenly distributed along the convex PF to ensure good population diversity. In [32, 35], MOEA/D attempts to introduce the nadir point as a reference point. In [32], the reverse PBI function is proposed, and the nadir point is used to solve the reverse PBI to maximize the value of the aggregate function, which improves the search performance of MOEA/D. In [35], if the candidate solutions obtained in the boundary area after several generations are less than the PF intermediate area, the reference point will change from the ideal point to the nadir point. Recently, Wang et al. [36] studied the effect of the difference between the ideal point and the nadir point on the performance of the algorithm and showed that they can complement each other.

As Wang et al. [36] pointed out, the use of the ideal point z^∗ and the nadir point z^nad in the Chebyshev function has an important influence on the distribution of the optimal solution on the PF. In particular, in the case where the ideal point z^∗ is used as a reference point, the optimal solutions of the subproblems of convex PF and concave PF are shown in Figures 5(a) and 5(b), respectively. It can be clearly seen that the optimal solution density of the central part of the convex PF is much larger than that of the concave PF, but it is opposite near the PF boundary. Compared with the ideal point z^∗, if the nadir point z^nad is used as a reference point, the distribution directions of the optimal solutions on these PFs are reversed, as shown in Figures 5(c) and 5(d), respectively. Since the final population distribution obtained by using the ideal point z^∗ and the nadir point z^nad is complementary, using them as reference points at the same time may improve their performance, making them approximate convex PF and concave PF. In addition, if the nadir point z^nad is not used as a reference point, you may face greater diversity risks when it is not easy to maintain diversity.

[figures omitted; refer to PDF]

3. The Proposed Cross-Reference Line Method

3.1. Cross-Reference Line Method

The cross-reference line is formed by matching the ideal point reference line and the nadir point reference line one by one. As shown in Figure 6, based on the nadir point and each reference point, the nadir point reference line is constructed, and based on the ideal point and each reference point, the ideal point reference line is constructed. The ideal point reference line and the nadir point reference line corresponding to each reference point intersect at this reference point and form an angle area in the feasible area. The boundary of the angle area is composed of an ideal point reference line and a nadir point reference line and the Pareto front edge being sandwiched. The candidate solution approaches the Pareto front under the pressure of convergence. With each generation of calculation, the candidate solution tends to move into the included angle area, which is called the attraction area.

[figure omitted; refer to PDF]

It is worth noting that if a reference line is defined as the line between the ideal point and the nadir point, the distance between the candidate solution on the reference line and the ideal point reference line and the nadir point reference line is zero. Therefore, it will have an absolute advantage and break the fairness of candidate evaluation. In order to solve this problem, the connection between the ideal point and the nadir point is defined as the penalty line of the cross-reference line, and a certain additional penalty value is added to the candidate solutions that fall on the connection line. Therefore, the candidate solution on the penalty line can only be considered as a contributing solution to other cross-reference lines near the penalty line.

3.2. DPD Evaluation Indicator Based on Cross-Reference Line Method

The DPD evaluation indicator of the cross-reference line method is based on the ideal point reference line distance d_∗ and the nadir point reference line distance d_nad. Among them, the equation of the ideal point reference line distance d_∗ and the nadir point reference line distance d_nad is as follows:$$\begin{array}{c}\text{(5)}& d_{\ast }=F\left(p\right)\sin \left(\stackrel{⟶}{Z^{\ast }r}|,|F|\left(p\right)\right),d_{\text{nad}}=F\left(p\right)\sin \left(\stackrel{⟶}{Z^{\text{nad}}r}|,|F|\left(p\right)\right),\end{array}$$where $F\left(p\right)$ represents the Euclidean distance from the ideal point or reference point to the candidate solution p, $\stackrel{⟶}{Z^{\ast }r}$ represents the vector from the ideal point to the reference point, $\stackrel{⟶}{Z^{\text{nad}}r}$ represents the nadir reference point vector, and $F\left(p\right)$ represents the vector from the ideal point or the nadir point to the candidate solution p.

Based on the cross-reference line method, the dominant penalty distance (DPD) indicator is defined as the maximum value of the ideal point reference line distance (d_∗) and the nadir point reference line weighted distance (d_nad). The weighting factor μ is used to verify the effectiveness, ensure effective enhancement of diversity, and improve the performance of MOEA-CRL. The equation of the cross-reference line DPD indicator of the nondominated candidate solution p is as follows:$$\begin{array}{}\text{(6)}& {\text{DPD}}_p=\max \left(d_{\ast },\mu d_{\text{nad}}\right),\end{array}$$where d_∗ is the distance from a candidate solution to the ideal point reference line, μ is the weight coefficient, and d_nad is the distance from a candidate solution to the nadir point reference line.

Taking the weighting factor μ = 1 as an example, the DPD indicator is dealing with different types of PF problems, as shown in Figure 7. According to (6), the angle area between the boundary line of μ = 1 and the vertical line of the nadir point reference line is dominated by the distance d_∗ of the nadir point reference line. The angle area between the boundary line of μ = 1 and the vertical line of the ideal point reference line is dominated by the ideal point reference line distance d_nad.

[figure omitted; refer to PDF]

The basic idea of the DPD indicator proposed in this paper is to combine the nadir point and the ideal point and use the cross-reference line as the evaluation reference. This method can not only effectively improve the diversity of candidate solutions near the coordinate axis in convex PF but also ensure the convergence under Pareto’s dominance theory. As shown in Figure 8, the evaluation method of the DPD indicator based on cross-reference lines is shown. According to the definition of the contributing solution, the nondominated candidate solutions with the smallest DPDp value in Figure 8 are the contributing solutions. This idea can be regarded as the combination of the Pareto theory of advantages and the distance evaluation between the cross-reference line and the candidate solution.

[figure omitted; refer to PDF]

3.3. Convergence and Diversity of the Cross-Reference Method

The cross-reference line method enables MOEAs to ensure good convergence and diversity when dealing with various types of PF problems. As shown in Figure 9, taking weight coefficient μ = 1 as examples, multiple sets of examples showing the two-objective minimization problem of concave PF, convex PF, and linear PF are shown. If there is a candidate solution corresponding to a reference point with the smallest DPD value and the candidate solution corresponding to the reference point has the smallest DPD value, then this candidate solution is referred to as the reference point contributing solution. In this paper, if a candidate solution p has the smallest DPD_p to a cross-reference line and the cross-reference line also has the smallest DPD_p to the candidate solution p, the candidate solution is defined as a contributing solution of the cross-reference line. According to DPD evaluation indicator equation (6) and the uniformly distributed reference point set, each of reference points has the unique contributing solution from all candidate solutions. Therefore, the contributing solution of a reference point will tend to be close to the boundary of μ = 1 as the iterative search proceeds. In short, the contributing solution tends to be close to both the nadir point reference line and the ideal point reference line.

[figure omitted; refer to PDF]

As shown in Figure 10, taking weight coefficient μ = 1 as examples, multiple sets of examples showing the two-objective minimization problem of concave PF, convex PF, and linear PF are shown. Each reference point in the reference point set matches the nadir point and the ideal point, so that the feasible search area is divided into multiple subspaces. According to DPD evaluation indicator equation (6) and the uniformly distributed reference point set, each of subspaces has unique contributing solution, thus ensuring the diversity and uniform distribution of the cross-reference line method.

[figure omitted; refer to PDF]

The cross-reference line method is an improvement to the reference line method, inheriting the advantages of the ideal point reference line in terms of convergence and adding the nadir point reference line to enhance diversity. The cross-reference line is the combination of the ideal point reference line and the nadir point reference line to divide the objective space into multiple subspaces, and unique candidate solution with the best convergence is kept in each subspace to ensure uniform distribution of the Pareto solution set.

3.4. Convergence and Diversity of Cross-Reference Method

As shown in Figure 11, the weighting factor μ = 1 is used as an example, and other parameters are the same as those in Figure 4. The distance d_∗ of the candidate solution (3, 6) to the ideal point reference line is smaller than the distance d_nad to the nadir point reference line, so the DPD indicator of (3, 6) is dominated by d_nad. The distance d_∗ from the candidate solution (3, 5) to the ideal point reference line is greater than the distance d_nad from the nadir point reference line, so the DPD indicator of (3, 5) is dominated by d_∗. It can be calculated that the DPD value of (3, 6) is greater than (3, 5), so (3, 5) in the candidate solution set A is the best. (8, 4) and (10, 4) can be compared in the same way. In summary, it can be seen that candidate solution set A is superior to candidate solution set B. According to the dominance theory, (3, 5) in solution set A dominates (3, 6) in candidate solution set B, and (8, 4) in solution set A dominates (10, 4) in candidate solution set B. Therefore, candidate solution set A is superior to candidate solution set B. Therefore, the judgment of the DPD evaluation indicator on the merits of the candidate solution set is the same as Pareto’s dominance theory.

[figure omitted; refer to PDF]

According to DPD evaluation indicator equation (6) and the uniformly distributed reference point set, each of subspaces has unique contributing solution and the contributing solutions of each subspace are jointly influenced by the ideal point reference line and the nadir point reference line. Therefore, Pareto incompatibility can be effectively avoided. The cross-reference line method not only solves the Pareto incompatibility but also can quickly obtain the nondominated candidate solution p by calculating the DPD_p indicator.

4. The Proposed Algorithm

In this section, we first describe the overall framework of the proposed MOEA-CRL in detail. Then, the case study on the implementation of the adaptive cross-reference line method is demonstrated in detail. Finally, the environment selection based on the DPD indicator is illustrated in detail, and the differences on environment selection between MOEA-CRL and other MOEAs are analyzed.

4.1. The General Framework of MOEA-CRL

In this section, the general framework of the MOEA-CRL will be elaborated through employing the cross-reference line method and the DPD indicator on the existing fundamentals of MOEAs. The cross-reference lines are formed by the intersection of the ideal point reference lines and the nadir point reference lines. As shown in Algorithm 1, the computing flow is divided into two steps. The first step is to initialize, and the second is to optimize.

Algorithm 1: General framework of MOEA-CRL.

 

Input: N (population size), N_R (number of reference point and archive size), M (number of objective)

 

Output: P (final population)

(1)

P ⟵ Ramdom Initialize (N, M);

(2)

R ⟵ Uniform Reference Point (NR, M);

(3)

A ⟵ P;

(4)

[Z^∗, Z^nad] ⟵ A;

(5)

R′ ⟵ R;

(6)

while termination condition not fulfilled

(7)

 P′ ⟵ Mating Selection (P, R′, Z^∗, Z^nad);

(8)

 O ⟵ Variation (P′, N);

(9)

 [A, R′, Z^∗, Z^nad] ⟵ Update RefLines (A ∪ O, R, Z^∗, Z^nad);

(10)

 P ⟵ Environment Selection (P∪O, R′, N, M, Z^∗, Z^nad);

(11)

end

(12)

return P;

The initialization provides preparation for MOEA-CRL. Firstly, a random initialized population is generated according to the initial parameters, which include the number of objectives M, the variable number D, the population size N, the maximum evolution generations, and so on. Secondly, the initial population is sorted according to the efficient nondominated sort (ENS) [37], and the nondominated solution set is copied to the initial archive A. Thirdly, the hyperplane is built based on initial archive A and the points on the hyperplane are uniformly sampled. Therefore, the uniform distribution of the initial reference point are insured.

The optimization is the core of the MOEA-CRL. The mating pool is established according to the tournament selection strategy. The DPD is used as an evaluation indicator to calculate the fitness value fitness_p of each candidate solution p to select individuals in the mating pool. The equation of distribution fitness_p can be expressed as$$\begin{array}{}\text{(7)}& fitnes{s}_p={\text{DPD}}_{\max }-\text{DPD}\left(\frac{\left\{p\right\}}{P,R^{\prime }}\right),\end{array}$$where DPD_max represents the maximum value of DPD in the population, R′ represents the updated reference point set, p represents a candidate solution, and P represents the total population. The establishment of the mating pool and the individual selection process based on the cross-reference line method and DPD indicator are operated in detail in Algorithm 2. It can be divided into two steps. In the first step, each individual is normalized and the fitness_p of each individual is calculated separately. In the second step, the individual selection is employed by the tournament selection strategy, which can randomly select two candidate solutions for comparison and then retain the individual with the larger fitness_p. In general, the mating pool with the number of N/2 is obtained, and a new population O with the number of N is obtained after the mutation operation.

Algorithm 2:Mating selection (P, R′, Z^∗, and Z^nad).

 

Input: P (population), R′ (set of new reference points), Z^∗ (ideal point), Z^nad (nadir point), M (number of objective)

 

Output: P′ (parents for variation)

(1)

for i = 1 to M do

(2)

 $f_i\left(p\right)\leftarrow f_i\left(p\right){\min }_{q\in P}{f}_i\left(q\right)$, $\forall p\in P$;

(3)

end

(4)

the fitness of each candidate solution is calculated;

(5)

$P^{\prime }\leftarrow \varphi$;

(6)

for i = 1 to |P| do

(7)

 Two candidate solutions p and q are selected randomly from P;

(8)

 if fitness_p > fitness_qthen

(9)

  $P^{\prime }\leftarrow P^{\prime }\cup \left\{p\right\}$;

(10)

 else

(11)

  $P^{\prime }\leftarrow P^{\prime }\cup \left\{q\right\}$;

(12)

 end if

(13)

end for

(14)

return P′;

4.2. The Cross-Reference Line Adaptation

The adaptive cross-reference line method is a key step in the optimization. As shown in Algorithm 3, the method contains six operations: (1) deleting the duplicate candidate solutions (the absolute error of each objective value of the two candidate solutions is less than the given precision ε, which is defined as the repeated candidate solution. In the paper, the precision is ε = 1e − 6.) and dominant solutions in the archive A, (2) updating the ideal points and nadir points, (3) normalizing the archive A and the reference point R and calculating the DPD indicator for each candidate solution, (4) calculating the contributing solution set and effective reference points, (5) updating the archive A, and (6) updating the reference point R′ of the cross-reference lines to update the cross-reference lines.

Algorithm 3:Update RefLines (R, A, Z^∗, and Z^nad).

 

Input: R (set of initial reference points), A (Initial archive), Z^∗ (ideal point), Z^nad (nadir point)

 

Output: R′ (set of new reference points), A′ (new archive), Z^∗′ (new ideal point), Z^nad′ (new nadir point)

(1)

Duplicate candidate solutions are deleted in A;

(2)

Dominated candidate solutions are deleted in A;

(3)

for i = 1 to M do

(4)

 if min_p∈Af_i (p) < Z^∗

(5)

  Z^∗′ ← min_p∈Af_i (p);

(6)

 else

(7)

  Z^∗′ ← Z^∗;

(8)

 end if

(9)

 if max _p∈Af_i (p) > Z^nadthen

(10)

  Z^nad′ ← max_p∈Af_i (p);

(11)

 else

(12)

  Z^nad′ ← Z^nad;

(13)

 end if

(14)

 f1_i (p) ← f_i (p) − Z^∗, $\forall p\in A\cup P$;

(15)

 f2_i (p) ← f_i (p) − Z^nad, $\forall p\in A\cup P$;

(16)

 $R_i\leftarrow R_i^{\ast }\left(Z^{nad\prime }-Z^{\ast \prime }\right)$;

(17)

 DPD (f (p), R_i) ←max (Caldis (f_i (p), R_i), μ^∗Caldis (f_i (p), R_i));

(18)

 end for

(19)

A^con ← {p ∈ A| ∃r ∈ R: DPD (f (p), r) = min_p∈A DPD (f (p), r)};

(20)

A′ ← A^con;

(21)

while |A′| < min (|R|, |A|) do

(22)

 p ← argmax_p1∈A\A′ min_p2∈A′ arccos (f (p1), f (p2))

(23)

 A′ ← A′ ∪ {p};

(24)

end

(25)

R^valid ← {r ∈ R| ∃r ∈ A^con: DPD (f (p), r) = min_r′∈R DPD (f (p), r′)};

(26)

R′ ← R^valid;

(27)

while |R′| < min (|R|, |A|) do

(28)

 p ← argmax_{p∈A′\R′} min_r∈R′ arccos (f (p), r);

(29)

 p′ ← projection (p, hyperplane);

(30)

 R′ ← R′ ∪ {f (p′)};

(31)

end

(32)

return A′, R′, Z^∗ and Z^nad;

In the second step of Algorithm 3, the update of the ideal point and the nadir point depends on the archive A for each generation, which provides the support for the normalization of different objective functions and calculation of the DPD indicator. In the third step of Algorithm 3, the archive A and the reference points are normalized to the same interval ${\displaystyle {\prod }_{i=1}^M\left[0,z_i^{\text{nad}}-z_i^{\ast }\right]}$, so the influence of the difference objective functions is eliminated, which is convenient to compare. In addition, the calculation of the reference line distance according to the DPD indicator is shown in Algorithm 4, and the specific equation is shown (6).

Algorithm 4:Caldis (P and R).

 

Input: P (population), R (set of reference points)

 

Output: Distance (perpendicular distance from the candidate solution to the reference line)

(1

)/^∗O denotes the coordinate origin of P/

(2

)for i = 1 to |N| do

(3

) for j = 1 to |R| do

(4

)  Distance (i, j) ←$F\left(p_i\right)\sin \left(\stackrel{⟶}{O{r}_j}|,|F|\left(p_i\right)\right)$;

(5

) end for

(6

)end for

(7

)return Distance;

In the fourth step of Algorithm 3, the contributing solutions will be calculated, and finally the contributing solution set is obtained. The contributing solution must satisfy that the solution p has the smallest DPD_p for a cross-reference line. Through the calculation of the DPD evaluation indicator, all contributing solutions are copied from A^con to the new archive A′, and the remaining space of A′ is filled up by candidate solutions from A\A′ one by one until A′ reaches its maximal size of min (|R|, |A|), where at each time, the candidate solution p having the maximum value of min_{p1∈A\A;p2∈A′} arccos (f (p1), f (p2)) in A\A′ is copied to A′, with arccos (f (p1), f (p2)) indicating the acute angle between p1 and p2 in objective space. In this way, the archive always contains a number of nondominated solutions with good distribution. It is worth noting that these nondominated solutions with good distribution, which are noncontributing and selected into A′, will serve as supplement for subsequent reference point updates.

The fifth and sixth steps of Algorithm 3 are the key of the cross-reference line adaptation method. Firstly, the valid reference point set R^valid is obtained. The valid reference points must satisfy two conditions at the same time: (1) the solution p has the smallest DPD_p for the cross-reference line and (2) the cross-reference line with the solution p has the smallest DPD_p. The calculation of the DPD indicator will be affected by the distance between the candidate solution p to the ideal point reference line and the nadir point reference line at the same time and is calculated according to (6). Subsequently, the valid reference point set R^valid is copied into the reference point set R′. Finally, the remainder of R′ is complemented by the candidate solutions in the new archive A′ until |R′| = min (|R|, |A′|) is satisfied. The complementary strategy is to calculate the maximum value of min_r∈R′ arccos (f (p), r), and the remaining part of R′ is the candidate solutions in A′ corresponding to these maximum values.

Figure 12 shows the adaptive update process of the reference point set R and the archive A. Firstly, four contributing solutions are obtained by calculating the DPD indicator of the candidate solutions to each reference point as indicated in Figure 12(a). Secondly, the four contributing solutions and the other two noncontributing solutions are copied into the new archive A′ as illustrated in Figure 12(b), where the two noncontributing solutions are nondominated solutions with good distribution and elected into A′ according to the maximum value of min_{p1∈A\A; p2∈A′} arccos (f (p1), f (p2)) in A\A′. Thirdly, the cross-reference lines which the contributing solution have the smallest DPD for are copied into the valid reference point set R^valid as shown in Figure 12(c). Finally, the projection points of the two candidate solutions in the new file A′ on the hyperplane and four valid reference points are copied into R′ as shown in Figure 12(d).

[figures omitted; refer to PDF]

The adaptation of reference points not only ensure their own uniformity but also reflect the geometric property of the PF. Therefore, after updating the reference points, the cross-reference lines can also adaptively update to improve the performances of MOEA-CRL for irregular PF.

4.3. Environmental Selection Based on DPD Indicator

The environment selection based on the DPD indicator presents as shown in Algorithm 5. Being similar to most MOEAs, MOEA-CRL uses an elite strategy to make environmental choices for each generation. It is worth noting that after normalization and ENS, the smallest k-th generation whose individual number reaches N needs to be selected to enhance diversity, and the DPD indicator is employed for the selection.

Algorithm 5: Environmental selection (P, R′, N, M, Z^∗, and Z^nad).

 Input: P (population), R′ (set of new reference points), N (size of

population), M (number of objective), Z^∗ (ideal point), Z^nad (nadir point)

 

Output: Q (population for next generation)

(1)

for i = 1 to M do

(2)

 f_i (p) ⟵ f_i (p) − min_q∈Af_i (q), $\forall p\in P$;

(3)

end

(4)

Front ← Non-Dominated Sort (P);

(5)

k ← Satisfied_min$\left|{\displaystyle {\cup }_{i=1}^{k}Fron{t}_i}\right|\ge N$;

(6)

$Q\leftarrow {\displaystyle {\cup }_{i=1}^{k-1}Fron{t}_i}$;

(7)

while |Front_k| > N − |Q|

(8)

 p ← argmax_p ∈ Front_k DPD (Front_k\{p}, R′);

(9)

 Front_k ← Front_k\{p};

(10)

end

(11)

Q ← Q ∪ Front_k;

(12)

return Q;

Although the selection of most decomposition-based evolutionary algorithms is guided by a set of reference points, the reference lines in MOEA-CRL have different purposes. In the MOEA-CRL, the cross-reference lines are adopted to calculate the DPD indicator to evaluate candidate solutions, but each candidate solution is associated with unique reference point in the decomposition-based MOEAs. Therefore, the population size of MOEA-CRL can be unequal to the number of reference points and is not necessarily the same as the method proposed by Das and Dennis [37].

In addition, MOEA-CRL adjusts the distribution of the cross-reference lines according to the contributing solutions directly defined by the DPD indicator, so the adaptation of the cross-reference line can be not affected by population size, and the uniform distribution of cross-reference lines can be maintained. Regardless of the size of the population, MOEA-CRL is always able to obtain uniformly distributed candidate solutions, providing the flexibility for population size settings. This conclusion is further evidenced by the empirical results in Section 5.4.

5. Experimental Results and Analysis

In this section, the sensitivity analysis of the DPD weight coefficients is firstly conducted. It not only proves the validity of the cross-reference line method but also offers the best weight coefficient μ of the MOEA-CRL. Subsequently, the parameters are set in detail. In Section 5.3, the proposed MOEA-CRL is compared with four existing popular MOEAs, including MOEA/D [16], NSGA-III [11], RVEA [20], and KnEA [38]. Finally, the sensitivity analysis of the population size of the MOEA-CRL was performed.

In the experiment, 19 test problems from three widely used test suites, including DTLZ1-DTLZ7 [39], WFG1-WFG9 [40], MaF3, MaF11, and MaF15 [41], were used to verify the algorithm in this paper. DTLZ1-DTLZ7 and WFG1-WFG9 are the problems of the quantity of scalable objectives, which are used to test the performances of the MOEAs on various MOPs and MaOPs. MaF3, MaF11, and MaF15 possess highly irregular “convex” PFs which can be used to test the performances of algorithms on the highly irregular “convex” PF.

5.1. Sensitivity Analysis of Weight Coefficients of MOEA-CRL

In the MOEA-CRL, the maximum of the ideal point reference line and the nadir point reference line is selected as the DPD indicator, so the choice of the DPD weight coefficient μ significantly affects the performance of the MOEA-CRL. The ideal point reference line determines the convergence performance of the MOEA-CRL, and the supplement of the nadir point reference line not only solves the Pareto incompatibility but also enhances diversity. The weighting coefficient μ changes the fairness of the ideal point reference line distance and the nadir point reference line penalty distance, which will affect the convergence and population diversity of MOEA-CRL. Therefore, the weight coefficient μ is an important factor that determines the performance of MOEA-CRL. A suitable μ is set to meet requirements of the convergence and diversity, which helps to enhance the flexibility of the MOEA-CRL for different MOPs and MaOPs.

In this section, in order to study the effect of the weight coefficient μ on the performance of MOEA-CRL, the different μ values are employed for performance comparison. In the experiment, μ was set to 1e − 6, 0.25, 0.5, 0.75, 1, 2.5, 5, and 7.5, respectively. μ was set to 1e − 6 to prove the role of the nadir point reference line and further prove the effectiveness of the cross-reference line method. In order to research the effect of the weight coefficient μ on the proposed MOEA-CRL, the test problems select three kinds of test problems “linear,” “convex,” and “concave” according to the feature of the PF, such as DTLZ1, DTLZ2, and MaF3. The DTLZ1, DTLZ2, and MaF3 were tested with 3 objectives. In addition, the other parameter settings of MOEA-CRL are the same as in Section 5.2. A box-plot of the DPD indicator obtained for the eight μ cases among the three test problems is indicated in Figure 13.

[figures omitted; refer to PDF]

It can be seen from Figure 13 that the mean DPD with the weight coefficient μ = 1e − 6 is larger than the mean DPD with μ = 0.25. The results show that the complement of the nadir point reference line even can increase the convergence pressure, which can improve the convergence performance of MOEA-CRL, in the case of a specific weight coefficient μ. It can be seen from Figure 13(a) that the mean DPD does not change significantly as the weight coefficient μ changes, which can be found that the weight coefficient μ has little influence for linear PF. As illustrated in Figure 13(b), the mean DPD changes drastically. When the weight coefficient μ > 0.25, the mean DPD increases as μ increases. It is particularly noteworthy that when μ > 1, the mean DPD increases dramatically. This indicates that when the weighting coefficient μ is so large that the convergence performance of the MOEA-CRL is significantly deteriorates for concave PF. The reason is that the population convergence pressure will gradually decrease for concave PF as the weight coefficient μ increases. It is shown in Figure 13(c) that the mean DPD with the weight coefficient μ = 1e − 6 is the biggest. This indicates that the complement of the nadir point reference line can enhance convergence for convex PF. It proves that it is valid to employ the nadir point reference line in the cross-reference line method as a strategy for evaluation. Furthermore, Figure 13(c) illustrates that the convergence performance of MOEA-CRL with the weight coefficient μ = 0.25 is optimal compared with other cases. Therefore, in the work of this paper, the weight coefficient μ is set to 0.25.

5.2. Experimental Settings

In order to compare fairly with existing advanced algorithms, this article uses general parameter settings, as follows:

(1)

Setting of the reference point. The reference point generation of MOEA/D, NSGA-III, and RVEA is on the basis of the two-layer method proposed by Das and Dennis [37]. Table 1 makes a list of the number of reference points in the test experiment for each objective quantity, in which p1 and p2 represent the number of divisions of each objective of the boundary layer and the inner layer, respectively. For fair comparison, MOEA-CRL also uses the same number of preset reference points listed in Table 2, and the population size of all MOEAs is the same as the number of reference points.

(2)

Relevant parameter settings of the competition algorithm. In MOEA/D, the size of neighborhood T is set to 1/10 of the population size, and the aggregate function used by the algorithm is the Chebyshev method. The penalty parameter α of RVEA is set to 2, and the reference point adaptive frequency fr is set to 0.1. The preset parameter T of KnEA is 0.5. There are no additional parameters for NSGA-III.

(3)

Genetic operation. The crossover operators in all experiments in this experiment are analog binary crossovers, and the mutation operators are polynomial mutations [42]. The distribution indicators of the crossover operators and the mutation operators are both set to 20, and the crossover probability and the mutation probability are set to 1.0 and 1/D, respectively, in which D represents the number of decision variables.

(4)

Performance indicators. The convergence and the diversity of the solution sets are indicated by the IGD and the hyper-volume (HV). In the HV calculation, all individuals of the population have been normalized, then the normalized HV value is calculated with a reference point (1.1, 1.1, …, 1.1). The MOEA with a larger HV value has better performance than the other. In addition, in order to reduce the computational complexity and improve the computational efficiency, the Monte Carlo estimation method is adopted for problems with the objective number is 5 and 10, and the number of sampling points required for the calculation is set to 1,000,000. In the DPD calculation, approximately 5,000 uniformly distributed points are sampled at the PF by the two-layer method proposed by Das and Dennis [37]. All tests were run 30 times independently, and the mean and standard deviation of each metric were recorded. The results of the experiment were statistically analyzed by the Wilcoxon rank sum test with a significance level of 5%, as Tables 3 and 4, where the symbol “+” indicates that the result of the other MOEA is significantly better, the “−” indicates that the result of the other MOEA is significantly worse, and “≈” indicates the similar performance of MOEA-CRL.

Table 1

Settings of the number of reference points for each number of objectives, where p1 and p2 denote the numbers of divisions on each objective for the boundary layer and the insider layer, respectively.

Number of objectives (M)	Parameter (p1, p2)	Number of reference points/population size (N)
3	13, 0	105
5	6, 0	210
10	3, 2	275

Table 2

Settings of the number of objectives, the number of decision variables, and the maximal number of generations for each test problem.

Test problem	M	D	Pareto front
Regular Pareto front
DTLZ1	3, 5, 10	M −1 + 5	Linear
DTLZ3	3, 5, 10	M −1 + 10	Concave
DTLZ2, 4	3, 5, 10	M −1 + 10	Concave
WFG4-9	3, 5, 10	M −1 + 10	Concave
Irregular Pareto front
DTLZ5-6	3, 5, 10	M −1 + 10	Mostly degenerate
DTLZ7	3, 5, 10	M −1 + 20	Disconnected
WFG1	3, 5, 10	M −1 + 10	Sharp tails
WFG2	3, 5, 10	M −1 + 10	Disconnected
WFG3	3, 5, 10	M −1 + 10	Mostly degenerate
MaF3	3, 5, 10	M −1 + 10	Convex
MaF11	3, 5, 10	M −1 + 10	Convex disconnected
MaF15	3, 5, 10	20∗M	Convex large scale

Table 3

Statistical results (mean values and standard deviations) of IGD value obtained by MOEA/D, NSGA-III, RVEA, KnEA, and MOEA-CRL on DTLZ1-DTLZ7, WFG1-WFG9, MaF3, MaF11, and MaF15 with 3 objectives.

Problem	M	D	MOEA/D	NSGA-III	RVEA	KnEA	MOEA-CRL
DTLZ1	3	7	2.8508e − 2 (4.30e − 6) −	1.8979e − 2 (4.28e − 6) ≈	1.8978e − 2 (5.63e − 6) ≈	5.0030e − 2 (2.51e − 2) −	1.8977e − 2 (2.29e − 5)
DTLZ2		12	6.9661e − 2 (5.68e − 5) −	5.0301e − 2 (4.52e − 7) −	5.0301e − 2 (4.21e − 7) −	6.6663e − 2 (2.44e − 3) −	4.6814e − 2 (5.07e − 5)
DTLZ3		12	1.0106e − 1 (1.20e − 1) −	5.0394e − 2 (1.67e − 4) −	5.0355e − 2 (7.33e − 5) −	1.0817e − 1 (2.89e − 2) −	4.7075e − 2 (1.21e − 4)
DTLZ4		12	2.3255e − 1 (3.36e − 1) −	1.3215e − 1 (1.86e − 1) ＋	5.0300e − 2 (4.47e − 7) ≈	1.5318e − 1 (2.69e − 1) ＋	2.2800e − 1 (2.42e − 1)
WFG4	3	12	2.7142e − 1 (6.38e − 4) −	2.0405e − 1 (3.67e − 5) −	2.0800e − 1 (2.21e − 3) −	2.4949e − 1 (7.34e − 3) −	1.9326e − 1 (4.61e − 4)
WFG5			2.8760e − 1 (9.10e − 4) −	2.1444e − 1 (3.75e − 5) −	2.1580e − 1 (6.17e − 4) −	2.6120e − 1 (1.07e − 2) −	2.0682e − 1 (2.34e − 4)
WFG6			2.9642e − 1 (1.27e − 2) −	2.1871e − 1 (8.11e − 3) ≈	2.2604e − 1 (1.00e − 2) −	2.8667e − 1 (1.28e − 2) −	2.1759e − 1 (1.29e − 2)
WFG7			2.7206e − 1 (6.18e − 4) −	2.0414e − 1 (5.96e − 5) −	2.0612e − 1 (9.83e − 4) −	2.4178e − 1 (9.90e − 3) −	1.9289e − 1 (3.86e − 4)
WFG8			3.1847e − 1 (6.38e − 3) −	2.6527e − 1 (2.75e − 3) −	2.8098e − 1 (5.23e − 3) −	3.3436e − 1 (7.90e − 3) −	2.6066e − 1 (2.48e − 3)
WFG9			2.7920e − 1 (3.28e − 2) −	2.0537e − 1 (7.59e − 4) −	2.0722e − 1 (1.72e − 3) −	2.2643e − 1 (6.75e − 3) −	1.9470e − 1 (8.91e − 4)
＋/−/≈			0/10/0	1/8/1	1/8/1	1/9/0
DTLZ5	3	12	1.2417e − 2 (1.51e − 6) −	1.1730e − 2 (1.24e − 3) −	5.8469e − 2 (1.08e − 3) −	1.0462e − 2 (2.65e − 3) −	5.7501e − 3 (3.73e − 4)
DTLZ6		12	1.2419e − 2 (7.06e − 7) −	1.7437e − 2 (2.69e − 3) −	5.9099e − 2 (3.01e − 3) −	4.7828e − 3 (3.60e − 4) −	4.2868e − 3 (3.45e − 5)
DTLZ7		22	2.4700e − 1 (8.06e − 2) −	7.0580e − 2 (2.35e − 3) ＋	1.0489e − 1 (3.58e − 3) ＋	8.4509e − 2 (7.48e − 2) ＋	2.0368e − 1 (1.79e − 1)
WFG1	3	12	2.1051e − 1 (1.22e − 2) −	1.3620e − 1 (2.39e − 3) ＋	1.5393e − 1 (6.17e − 3) −	1.8632e − 1 (7.92e − 3) −	1.4766e − 1 (2.95e − 3)
WFG2		12	2.1674e − 1 (1.53e − 3) −	1.5105e − 1 (1.13e − 3) ≈	1.6491e − 1 (5.59e − 3) −	1.8398e − 1 (8.28e − 3) −	1.5058e − 1 (9.94e − 4)
WFG3		12	4.1686e − 2 (3.12e − 4) ＋	8.8550e − 2 (6.34e − 3) ＋	2.1668e − 1 (8.51e − 3) −	9.6049e − 2 (8.33e − 3) −	8.9822e − 2 (8.87e − 3)
Convex
MaF3	3	12	1.5083e − 1 (9.68e − 2) −	4.3564e − 2 (3.22e − 4) −	3.8402e − 2 (4.17e − 4) −	1.4269e − 1 (6.04e − 2) −	3.4876e − 2 (8.30e − 4)
MaF11		12	2.1673e − 1 (1.60e − 3) −	1.5147e − 1 (9.71e − 4) ≈	1.6463e − 1 (4.37e − 3) −	1.8721e − 1 (9.46e − 3) −	1.5082e − 1 (1.11e − 3)
MaF15		5	4.5845e − 1 (2.10e − 1) −	3.6301e − 1 (1.82e − 1) −	9.0510e − 1 (2.23e − 1) −	5.1952e − 1 (1.64e − 1) −	2.6890e − 1 (3.57e − 2)
＋/−/≈			1/8/0	3/4/2	1/8/0	1/8/0

Table 4

Statistical results (mean values and standard deviations) of HV value obtained by MOEA/D, NSGA-III, RVEA, KnEA, and MOEA-CRL on DTLZ1-DTLZ7, WFG1-WFG9, MaF3, MaF11, and MaF15 with 5 objectives and 10 objectives.

Problem	M	MOEA/D	NSGA-III	RVEA	KnEA	MOEA-CRL
DTLZ1	5	9.0873e − 1 (8.52e − 2) −	9.7979e − 1 (1.56e − 4) ≈	9.7984e − 1 (1.55e − 4) ≈	6.2616e − 1 (1.64e − 1) −	9.7988e − 1 (1.63e − 4)
	10	9.7273e − 1 (5.32e − 3) −	9.8648e − 1 (3.26e − 2) −	9.9967e − 1 (1.88e − 5) ≈	0.0000e + 0 (0.00e + 0) −	9.9973e − 1 (1.10e − 4)
DTLZ2	5	7.1112e − 1 (5.08e − 4) −	8.1269e − 1 (4.52e − 4) −	8.1252e − 1 (4.48e − 4) −	7.9064e − 1 (3.54e − 3) −	8.1626e − 1 (9.62e − 4)
	10	6.2665e − 1 (2.70e − 2) −	9.4539e − 1 (3.51e − 2) −	9.6963e − 1 (1.72e − 4) −	9.5650e − 1 (4.53e − 3) −	9.7102e − 1 (1.70e − 3)
DTLZ3	5	4.3566e − 1 (1.27e − 1) −	7.7543e − 1 (2.88e − 3) ＋	7.7757e − 1 (1.36e − 3) ＋	3.9573e − 1 (1.31e − 1) −	7.3975e − 1 (3.36e − 3)
	10	6.4046e − 1 (1.97e − 2) −	3.3590e − 1 (4.26e − 1) −	9.6443e − 1 (9.40e − 3) ＋	0.0000e + 0 (0.00e + 0) −	9.6234e − 1 (2.68e − 1)
DTLZ4	5	2.7690e − 1 (1.65e − 1) −	7.3262e − 1 (6.50e − 2) ＋	7.7683e − 1 (1.68e − 2) ＋	7.6308e − 1 (4.31e − 3) ＋	7.2000e − 1 (5.17e − 2)
	10	6.2169e − 1 (2.11e − 2) −	9.6721e − 1 (1.24e − 2) −	9.6984e − 1 (1.99e − 4) −	9.5637e − 1 (4.37e − 3) −	9.7165e − 1 (1.51e − 3)
WFG4	5	6.2968e − 1 (8.71e − 4) −	8.0469e − 1 (8.35e − 4) −	8.0565e − 1 (1.10e − 3) ≈	7.8730e − 1 (2.17e − 3) −	8.0575e − 1 (1.74e − 3)
	10	5.8107e − 1 (5.61e − 2) −	9.4578e − 1 (3.92e − 3) ＋	9.4326e − 1 (3.57e − 3) ＋	9.5750e − 1 (1.56e − 3) ＋	9.1518e − 1 (5.08e − 3)
WFG5	5	5.9294e − 1 (1.98e − 2) −	7.6126e − 1 (3.68e − 4) ＋	7.6092e − 1 (4.04e − 4) ＋	7.4495e − 1 (2.82e − 3) ＋	7.3465e − 1 (2.49e − 3)
	10	5.3478e − 1 (1.75e − 2) −	8.9907e − 1 (6.43e − 4) ＋	8.9790e − 1 (1.17e − 3) ＋	8.9664e − 1 (8.14e − 4) ＋	8.4566e − 1 (3.56e − 3)
WFG6	5	5.5065e − 1 (3.13e − 2) −	7.4242e − 1 (1.58e − 2) ＋	7.4662e − 1 (1.44e − 2) ＋	7.2041e − 1 (1.03e − 2) ≈	7.2008e − 1 (1.20e − 2)
	10	4.7917e − 1 (9.01e − 2) −	8.6935e − 1 (1.66e − 2) ≈	8.6291e − 1 (2.01e − 2) ≈	8.6854e − 1 (1.59e − 2) ≈	8.7014e − 1 (1.72e − 2)
WFG7	5	6.2952e − 1 (1.76e − 3) −	8.0771e − 1 (6.15e − 4) ＋	8.0685e − 1 (6.00e − 4) ＋	7.9537e − 1 (2.41e − 3) ＋	7.9284e − 1 (1.74e − 3)
	10	6.1632e − 1 (3.69e − 2) −	9.4638e − 1 (2.02e − 2) ＋	9.4552e − 1 (2.93e − 3) −	9.5737e − 1 (6.09e − 3) ＋	9.4268e − 1 (2.70e − 3)
WFG8	5	3.2668e − 1 (1.10e − 2) −	6.9440e − 1 (3.44e − 3) −	6.9749e − 1 (1.55e − 3) ≈	6.6088e − 1 (4.04e − 3) −	6.9786e − 1 (2.26e − 3)
	10	5.2900e − 1 (2.26e − 2) −	8.3115e − 1 (2.65e − 2) −	7.4198e − 1 (7.45e − 2) −	8.1031e − 1 (5.88e − 2) −	8.5967e − 1 (2.09e − 2)
WFG9	5	4.1801e − 1 (7.01e − 2) −	7.6559e − 1 (4.98e − 3) ＋	7.6952e − 1 (3.01e − 3) ＋	7.6796e − 1 (3.00e − 3) ＋	7.4665e − 1 (4.89e − 3)
	10	5.3897e − 1 (5.22e − 2) －	8.6420e − 1 (4.42e − 2) ＋	8.7342e − 1 (1.19e − 2) ＋	9.0581e − 1 (3.21e − 2) ＋	8.3182e − 1 (1.92e − 2)
＋/−/≈		0/20/0	10/8/2	10/5/5	8/10/2
DTLZ5	5	1.0968e − 1 (8.30e − 3) ≈	1.0483e − 1 (1.62e − 2) ≈	9.1890e − 2 (1.22e − 3) −	7.1580e − 2 (2.93e − 2) −	1.0808e − 1 (5.72e − 3)
	10	9.7922e − 2 (3.20e − 4) ＋	2.8809e − 2 (3.32e − 2) −	9.0906e − 2 (1.34e − 4) ＋	3.5304e − 2 (3.14e − 2) −	8.2365e − 2 (2.16e − 2)
DTLZ6	5	9.4634e − 2 (6.13e − 3) ≈	5.7651e − 2 (4.09e − 2) −	9.9670e − 2 (4.73e − 3) ≈	9.1029e − 2 (2.57e − 3) −	9.2650e − 2 (2.20e − 2)
	10	9.8124e − 2 (2.62e − 4) ＋	3.0230e − 3 (1.66e − 2) −	9.2019e − 2 (9.23e − 4) ＋	0.0000e + 0 (0.00e + 0) −	9.1491e − 2 (3.74e − 2)
DTLZ7	5	1.7168e − 1 (5.33e − 2) −	2.3157e − 1 (9.31e − 3) ＋	2.1343e − 1 (8.52e − 3) −	2.4843e − 1 (1.07e − 2) ＋	2.2732e − 1 (3.03e − 3)
	10	4.5806e − 3 (7.31e − 3) −	1.7116e − 1 (7.07e − 3) ＋	1.3371e − 1 (2.42e − 2) ＋	9.1243e − 2 (2.89e − 2) ＋	7.4082e − 2 (2.95e − 2)
WFG1	5	9.4334e − 1 (7.14e − 2) −	9.7370e − 1 (2.11e − 2) −	9.8276e − 1 (2.89e − 2) ≈	9.9246e − 1 (1.41e − 3) −	9.9835e − 1 (1.60e − 2)
	10	5.0234e − 1 (1.69e − 1) −	9.4008e − 1 (4.39e − 2) ＋	9.9015e − 1 (2.42e − 2) ＋	9.9756e − 1 (8.35e − 4) ＋	7.4040e − 1 (5.17e − 2)
WFG2	5	9.6566e − 1 (3.42e − 2) −	9.9555e − 1 (8.56e − 4) −	9.9404e − 1 (1.15e − 3) −	9.9320e − 1 (7.41e − 4) −	9.9683e − 1 (7.33e − 4)
	10	9.9648e − 1 (2.10e − 3) ＋	9.9700e − 1 (1.46e − 3) ＋	9.8471e − 1 (4.33e − 3) −	9.9306e − 1 (1.47e − 3) ＋	9.9191e − 1 (2.27e − 3)
WFG3	5	9.1998e − 2 (3.47e − 4) −	1.6765e − 1 (1.32e − 2) ≈	1.5967e − 1 (1.56e − 2) ≈	7.7484e − 2 (1.82e − 2) −	1.6083e − 1 (1.27e − 2)
	10	7.5376e − 2 (7.87e − 3) ＋	2.4516e − 4 (1.34e − 3) −	0.0000e + 0 (0.00e + 0) −	0.0000e + 0 (0.00e + 0) −	4.0606e − 2 (7.87e − 3)
MaF3	5	9.9646e − 1 (1.07e − 4) −	9.9870e − 1 (2.01e − 3) −	9.9895e − 1 (4.31e − 4) −	8.8686e − 1 (1.27e − 1) −	9.9975e − 1 (7.14e − 4)
	10	9.9993e − 1 (4.54e − 5) ≈	2.9125e − 1 (4.53e − 1) −	9.8343e − 1 (5.59e − 2) −	0.0000e + 0 (0.00e + 0) −	9.9997e − 1 (1.30e − 1)
MaF11	5	9.7508e − 1 (2.07e − 2) −	9.9583e − 1 (5.38e − 4) −	9.9397e − 1 (1.11e − 3) −	9.9299e − 1 (8.84e − 4) −	9.9751e − 1 (5.92e − 4)
	10	9.3320e − 1 (6.36e − 2) −	9.9641e − 1 (2.18e − 3) −	9.8256e − 1 (4.48e − 3) −	9.9284e − 1 (8.14e − 4) −	9.9818e − 1 (2.03e − 3)
MaF15	5	2.9995e − 2 (1.11e − 2) ＋	0.0000e + 0 (0.00e + 0) −	2.4759e − 2 (1.06e − 2) ＋	0.0000e + 0 (0.00e + 0) −	1.0571e − 2 (5.05e − 3)
	10	3.0314e − 11 (1.16e − 10) −	0.0000e + 0 (0.00e + 0) −	1.9696e − 7 (2.48e − 7) −	0.0000e + 0 (0.00e + 0) −	2.6803e − 7 (2.42e − 10)
＋/−/≈		5/10/3	4/12/2	5/10/3	4/14/0

5.3. Comparisons between MOEA-CRL and Existing MOEAs

Table 3 lists the comparison of the mean values of the IGD results on the test problems with 3 objectives between MOEA-CRL and four popular MOEAs. It can be seen from the evaluation results of the mean value of the IGD in Table 3 that the MOEA-CRL proposed in this paper is superior to the other four MOEAs in dealing with test problems for 3 objectives. Among 10 test problems with regular PFs, 9 test problems except DTLZ4 obtained the best solutions by MOEA-CRL. The mean value of the IGD of MOEA-CRL was slightly larger than RVEA on DTLZ1 and NSGA-III on WFG6. The MOEA-CRL also shows the good performance on the test problems with 9 irregular PFs. Especially, on the MaF3, MaF11, and MaF15 problems with three concave PF, MOEA-CRL shows the best performance.

Figure 14 plots the nondominated solution sets for each algorithm that obtained on DTLZ1, DTLZ2, and MaF3 problems with three objectives, which are the mean values of the IGD after 30 runs. It can be further observed from Figure 14 that the MOEA-CRL obtains a uniformly distributed nondominated solution set on DTLZ1, DTLZ2, and MaF3. It can be seen that the proposed MOEA-CRL can not only perform well on the linear PF and concave PF, but also on the convex PF in MOPs. Specifically, most popular MOEAs perform well for the linear and concave PF such as DTLZ1 and DTLZ2, but the population diversity will deteriorate significantly on the test problems with convex PF and the Pareto solution sets lost the uniform distribution. As shown in Table 3 and Figure 14, it indicates that MOEA-CRL performs good convergence and diversity for MOPs, especially on the convex PF.

[figure omitted; refer to PDF]

Table 4 lists the HV values gained by the MOEA-CRL and the four popular MOEAs on the test problems with 5 objectives and 10 objectives. Overall, MOEA-CRL achieved the best performance 16 times in a total of 38 experiments, while MOEA/D, NSGA-III, RVEA, and KnEA achieved the best performance 5, 6, 6, and 5 times, respectively. The evaluation results show that the overall performance of MOEA-CRL in dealing with MaOPs is better than the other four MOEAs, but the performance of MOEA-CRL deteriorates with the increase of dimension. For the regular Pareto test problems with 10 objectives, the advantage of MOEA-CRL is not significant. For the 9 irregular Pareto test problems, MOEA-CRL is competitive. Especially, MOEA-CRL performs better on the three concave PF test problems, MaF3, MaF11, and MaF15. Since the introduction of the nadir point will bring the performance gain on the convex PF more obvious, the cross-reference line method performs better on the convex PF than other multiobjective algorithms.

Figure 15 plots the nondominated solution sets for each algorithm that obtained on DTLZ1, DTLZ2, and MaF3 with 10 objectives, which are the mean values of the HV after 30 runs. The parallel coordinate is a way of data visualization. Multiple vertical and parallel coordinate axes represent multiple dimensions, and the scale on the dimension represents the corresponding value on the objective. Each sample corresponds to a value in each dimension, and a connected polyline represents the sample. It can be further observed from Figure 15 that the MOEA-CRL obtains the uniformly distributed nondominated solution sets on DTLZ1, DTLZ2, and MaF3. The MOEA-CRL can still effectively enhance the population diversity. Specifically, MOEAs other than MOEA/D can maintain diversity and obtain the good uniformity of nondominated solution sets in test problems on DTLZ1 and DTLZ2 with the 10 objectives. However, the MOEA-CRL obtained a uniform distribution in test problem on MaF3 with the 10 objectives and the performance of other four MOEAs significantly deteriorates. Therefore, MOEA-CRL can maintain better population diversity on MaOPs as shown in Table 4 and Figure 15, especially on the convex PF. It should be noted that the convergence of MOEA-CRL deteriorates significantly as the dimension increases. The reason is that the nadir point reference lines are employed to enhance the diversity, which results in a decline of the convergence pressure.

[figure omitted; refer to PDF]

Through the sensitivity analysis of the weight coefficient of the 3-objective problem, it can be concluded that the use of the nadir point does not weaken the convergence but increases the convergence pressure. However, as the objective dimension increases, the use of the nadir point will indeed cause the deterioration of convergence and even the failure to converge.

5.4. Sensitivity Analysis of Population Size

In the experiment, the population size N is set as same as the number of reference points since the reference points are generally associated with each candidate solution in most decomposition-based MOEAs. The number of reference points depends on the method of Das and Dennis [37].

The population size setting of MOEA-CRL proposed in this paper is flexible, and the number of reference points has less influence on it. The number of candidate solutions can be less than the number of reference points and also can be larger.

MOEA-CRL with different population sizes was tested on DTLZ1, DTLZ2, and MaF3 with three objectives. Figure 16 exhibits the nondominated solution sets of MOEA-CRL with population sizes of 35, 70, 105, 140, and 175, and the number of reference points is always set to 105. It can be shown that the Pareto solution sets obtained by MOEA-CRL are always uniformly distributed regardless of the population sizes. Therefore, MOEA-CRL provides greater flexibility for population size setting.

[figure omitted; refer to PDF]

6. Conclusions and Remarks

In this paper, an evolutionary algorithm based on the adaptive cross-reference line method, called MOEA-CRL, is proposed to inherit the advantages of the ideal point reference line for better convergence and add the nadir point reference line for higher diversity. Especially, on the convex PF, MOEA-CRL solves the Pareto incompatibility problem and significantly enhances the population diversity. Furthermore, this paper proposed the DPD indicator based on the cross-reference lines. The properties of the ideal point reference line and the nadir point reference line are combined to solve the Pareto incompatibility problem as well as improve the performance of the MOEA-CRL on the convex PF. Based on the DPD evaluation strategy of the cross-reference line method, MOEA-CRL retains unique solution with the best convergence in each attraction region as a nondominated solution, which ensures that the Pareto solution set is distributed evenly. Finally, this paper proposed a cross-reference line adaptation method in order to enhance the performance of MOEA-CRL in dealing with the irregular problems.

The experimental results show the superiority of MOEA-CRL on the convex PF. It also has the competitiveness due to the adaptability of cross-reference lines while solving those MOPs and MaOPs with other types of PFs. Remarkably, the cross-reference line method is only used to calculate the DPD indicator. Therefore, the population size is irrelated to the number of the cross-reference lines, and subsequently, the population size setting is flexible. The proposed MOEA-CRL proves that the adaptive cross-reference line method is prospective for significantly improving the diversity especially in the convex PF.

In fact, the experimental results also clearly illustrate that the performance of MOEA-CRL deteriorates significantly with the increase of dimensions. That means that the cross-reference line method still poses the challenges in dealing with some research issues such as high-dimensional deterioration and more complex convex PF problems.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (nos. 51875419 and 51605345).

References

[1] C. A. C. Coello, "Evolutionary multi-objective optimization: a historical view of the field," IEEE Computational Intelligence Magazine, vol. 1 no. 1, pp. 28-36, DOI: 10.1109/mci.2006.1597059, 2006.

[2] C. A. C. Coello, M. S. Lechuga, "MOPSO: a proposal for multiple objective particle swarm optimization," Proceedings of the 2002 Congress on Evolutionary Computation. CEC’02 (Cat. No.02TH8600),DOI: 10.1109/cec.2002.1004388, .

[3] L. Zuo, S. Lei, S. Dong, C. Zhu, T. Hara, "A multi-objective optimization scheduling method based on the ant colony algorithm in cloud computing," IEEE Access, vol. 3, pp. 2687-2699, DOI: 10.1109/access.2015.2508940, 2015.

[4] H. M. Alsaket, R. Mahmoud, H. M. Elattar, M. A. Aboul-Dahab, "Exploring evolutionary multi-objective techniques in self-organizing networks," IEEE Access, vol. 5, pp. 12049-12060, DOI: 10.1109/access.2017.2706188, 2017.

[5] D. Alanis, J. Hu, P. Botsinis, Z. Babar, S. X. Ng, L. Hanzo, "Quantum-assisted joint multi-objective routing and load balancing for socially-aware networks," IEEE Access, vol. 4, pp. 9993-10028, DOI: 10.1109/ACCESS.2016.2629671, 2016.

[6] H.-L. Liu, L. Chen, K. Deb, E. Goodman, "Investigating the effect of imbalance between convergence and diversity in evolutionary multi-objective algorithms," IEEE Transactions on Evolutionary Computation, vol. 21 no. 3, pp. 408-425, DOI: 10.1109/tevc.2016.2606577, 2017.

[7] Y. Tian, R. Cheng, X. Zhang, Y. Jin, "PlatEMO: a MATLAB platform for evolutionary multi-objective optimization [educational forum]," IEEE Computational Intelligence Magazine, vol. 12 no. 4, pp. 73-87, DOI: 10.1109/mci.2017.2742868, 2017.

[8] E. Zitzler, M. Laumanns, L. Thiele, "SPEA2: improving the strength pareto evolutionary algorithm," TIK-report, vol. 103, 2001.

[9] K. Deb, A. Pratap, S. Agarwal, T. Meyarivan, "A fast and elitist multiobjective genetic algorithm: NSGA-II," IEEE Transactions on Evolutionary Computation, vol. 6 no. 2, pp. 182-197, DOI: 10.1109/4235.996017, 2002.

[10] R. Cheng, Y. Jin, K. Narukawa, B. Sendhoff, "A multiobjective evolutionary algorithm using Gaussian process-based inverse modeling," IEEE Transactions on Evolutionary Computation, vol. 19 no. 6, pp. 838-856, DOI: 10.1109/tevc.2015.2395073, 2015.

[11] K. Deb, H. Jain, "An evolutionary many-objective optimization algorithm using reference-point-based nondominated sorting approach, part I: solving problems with box constraints," IEEE Transactions on Evolutionary Computation, vol. 18 no. 4, pp. 577-601, DOI: 10.1109/tevc.2013.2281535, 2014.

[12] N. Kowatari, A. Oyama, H. E. Aguirre, K. Tanaka, "A study on large population MOEA using adaptive ε -box dominance and neighborhood recombination for many-objective optimization," Proceedings of the International Conference on Learning and Intelligent Optimization, pp. 86-100, DOI: 10.1007/978-3-642-34413-8_7, .

[13] S. Kukkonen, K. Deb, "Improved pruning of non-dominated solutions based on crowding distance for bi-objective optimization problems," Proceedings of the 2006 IEEE International Conference on Evolutionary Computation, pp. 1179-1186, DOI: 10.1109/CEC.2006.1688443, .

[14] D. W. Corne, N. R. Jerram, J. D. Knowles, M. J. Oates, "PESA-II: region-based selection in evolutionary multiobjective optimization," pp. 283-290, .

[15] J. Horn, N. Nafpliotis, D. E. Goldberg, "A niched Pareto genetic algorithm for multiobjective optimization," Proceedings of the First IEEE Conference on Evolutionary Computation, IEEE World Congress on Computational Intelligence, pp. 82-87, DOI: 10.1109/icec.1994.350037, .

[16] Q. Zhang, H. Li, "MOEA/D: a multiobjective evolutionary algorithm based on decomposition," IEEE Transactions on Evolutionary Computation, vol. 11 no. 6, pp. 712-731, DOI: 10.1109/tevc.2007.892759, 2007.

[17] H. Ishibuchi, T. Murata, "A multi-objective genetic local search algorithm and its application to flowshop scheduling," IEEE Transactions on Systems, Man and Cybernetics, Part C (Applications and Reviews), vol. 28 no. 3, pp. 392-403, DOI: 10.1109/5326.704576, 1998.

[18] T. Murata, H. Ishibuchi, M. Gen, "Specification of genetic search directions in cellular multi-objective genetic algorithms," Proceedings of the International Conference on Evolutionary Multi-Criterion Optimization, pp. 82-95, DOI: 10.1007/3-540-44719-9_6, .

[19] E. J. Hughes, "MSOPS-II: a general-purpose many-objective optimiser," Proceedings of the 2007 IEEE Congress on Evolutionary Computation, pp. 3944-3951, DOI: 10.1109/cec.2007.4424985, .

[20] R. Cheng, Y. Jin, M. Olhofer, B. Sendhoff, "A reference vector guided evolutionary algorithm for many-objective optimization," IEEE Transactions on Evolutionary Computation, vol. 20 no. 5, pp. 773-791, DOI: 10.1109/tevc.2016.2519378, 2016.

[21] H.-L. Liu, F. Gu, Q. Zhang, "Decomposition of a multiobjective optimization problem into a number of simple multiobjective subproblems," IEEE Transactions on Evolutionary Computation, vol. 18 no. 3, pp. 450-455, DOI: 10.1109/tevc.2013.2281533, 2014.

[22] S. Jiang, S. Yang, "A strength Pareto evolutionary algorithm based on reference direction for multiobjective and many-objective optimization," IEEE Transactions on Evolutionary Computation, vol. 21 no. 3, pp. 329-346, DOI: 10.1109/tevc.2016.2592479, 2017.

[23] Y. Qi, X. Ma, F. Liu, L. Jiao, J. Sun, J. Wu, "MOEA/D with adaptive weight adjustment," Evolutionary Computation, vol. 22 no. 2, pp. 231-264, DOI: 10.1162/evco_a_00109, 2014.

[24] Y. Sun, G. G. Yen, Z. Yi, "IGD indicator-based evolutionary algorithm for many-objective optimization problems," IEEE Transactions on Evolutionary Computation, vol. 23, pp. 173-187, 2018.

[25] Y. Liu, J. Wei, X. Li, M. Li, "Generational distance indicator-based evolutionary algorithm with an improved niching method for many-objective optimization problems," IEEE Access, vol. 7, pp. 63881-63891, DOI: 10.1109/access.2019.2916634, 2019.

[26] J. Li, G. Chen, M. Li, H. Chen, "An adaptative reference vector based evolutionary algorithm for many-objective optimization," IEEE Access, vol. 7, pp. 80506-80518, DOI: 10.1109/access.2019.2923422, 2019.

[27] Y. Yuan, H. Xu, B. Wang, X. Yao, "A new dominance relation-based evolutionary algorithm for many-objective optimization," IEEE Transactions on Evolutionary Computation, vol. 20 no. 1, pp. 16-37, DOI: 10.1109/tevc.2015.2420112, 2016.

[28] Y. Sun, G. G. Yen, Z. Yi, "Reference line-based estimation of distribution algorithm for many-objective optimization," Knowledge-Based Systems, vol. 132, pp. 129-143, DOI: 10.1016/j.knosys.2017.06.021, 2017.

[29] Y. Tian, R. Cheng, X. Zhang, F. Cheng, Y. Jin, "An indicator-based multiobjective evolutionary algorithm with reference point adaptation for better versatility," IEEE Transactions on Evolutionary Computation, vol. 22, pp. 609-622, DOI: 10.1109/tevc.2017.2749619, 2018.

[30] Y. Yuan, H. Xu, B. Wang, B. Zhang, X. Yao, "Balancing convergence and diversity in decomposition-based many-objective optimizers," IEEE Transactions on Evolutionary Computation, vol. 20 no. 2, pp. 180-198, DOI: 10.1109/TEVC.2015.2443001, 2015.

[31] A. Mohammadi, M. N. Omidvar, X. Li, K. Deb, "Sensitivity analysis of penalty-based boundary intersection on aggregation-based EMO algorithms," Proceedings of the 2015 IEEE Congress on Evolutionary Computation (CEC), pp. 2891-2898, DOI: 10.1109/CEC.2015.7257248, .

[32] H. Sato, "Inverted PBI in MOEA/D and its impact on the search performance on multi and many-objective optimization," Proceedings of the 2014 Annual Conference on Genetic and Evolutionary Computation, pp. 645-652, DOI: 10.1145/2576768.2598297, .

[33] H. Ishibuchi, K. Doi, H. Masuda, Y. Nojima, "Relation between weight vectors and solutions in MOEA/D," Proceedings of the 2015 IEEE Symposium Series on Computational Intelligence, pp. 861-868, DOI: 10.1109/ssci.2015.127, .

[34] Q. Zhang, H. Li, D. Maringer, E. Tsang, "MOEA/D with NBI-style Tchebycheff approach for portfolio management," Proceedings of the IEEE Congress on Evolutionary Computation,DOI: 10.1109/cec.2010.5586185, .

[35] S. Jiang, S. Yang, "An improved multiobjective optimization evolutionary algorithm based on decomposition for complex Pareto fronts," IEEE Transactions on Cybernetics, vol. 46, pp. 421-437, DOI: 10.1109/tcyb.2015.2403131, 2016.

[36] Z. Wang, Q. Zhang, H. Li, H. Ishibuchi, L. Jiao, "On the use of two reference points in decomposition based multiobjective evolutionary algorithms," Swarm and Evolutionary Computation, vol. 34, pp. 89-102, DOI: 10.1016/j.swevo.2017.01.002, 2017.

[37] I. Das, J. E. Dennis, "Normal-boundary intersection: a new method for generating the pareto surface in nonlinear multicriteria optimization problems," SIAM Journal on Optimization, vol. 8 no. 3, pp. 631-657, DOI: 10.1137/s1052623496307510, 1998.

[38] X. Zhang, Y. Tian, Y. Jin, "A knee point-driven evolutionary algorithm for many-objective optimization," IEEE Transactions on Evolutionary Computation, vol. 19 no. 6, pp. 761-776, DOI: 10.1109/tevc.2014.2378512, 2015.

[39] K. Deb, L. Thiele, M. Laumanns, E. Zitzler, "Scalable test problems for evolutionary multiobjective optimization," Evolutionary Multiobjective Optimization, pp. 105-145, DOI: 10.1007/1-84628-137-7_6, 2005.

[40] S. Huband, P. Hingston, L. Barone, L. While, "A review of multiobjective test problems and a scalable test problem toolkit," IEEE Transactions on Evolutionary Computation, vol. 10 no. 5, pp. 477-506, DOI: 10.1109/tevc.2005.861417, 2006.

[41] R. Cheng, M. Li, Y. Tian, "Benchmark functions for the CEC’2018 competition on many-objective optimization," 2018. Technical report

[42] K. Deb, M. Goyal, "A combined genetic adaptive search (GeneAS) for engineering design," Computer Science and Informatics, vol. 26, pp. 30-45, 1996.