Covariance matrix adaptive strategy for a multi-objective evolutionary algorithm based on reference point

Abstract

In this article, an effective method, called an adaptive covariance strategy based on reference points (RPCMA-ES) is proposed for multi-objective optimization. In the proposed algorithm, search space is divided into independent sub-regions by calculating the angle between the objective vector and the reference vector. The reference vectors can be used not only to decompose the original multi-objective optimization problem into a number of single-objective subproblems, but also to elucidate user preferences to target a preferred subset of the whole Pareto front (PF). In this respect, any single objective optimizers can be easily used in this algorithm framework. Inspired by the multi-objective estimation of distribution algorithms, covariance matrix adaptation evolution strategy (CMA-ES) is involved in RPCMA-ES. A state-of-the-art optimizer for single-objective continuous functions is the CMA-ES, which has proven to be able to strike a good balance between the exploration and the exploitation of search space. Furthermore, in order to avoid falling into local optimality and make the new mean closer to the optimal solution, chaos operator is added based on CMA-ES. By comparing it with four state-of-the-art multi-objective optimization algorithms, the simulation results show that the proposed algorithm is competitive and effective in terms of convergence and distribution.

Keywords

Multi-objective optimization problem Reference point Covariance matrix adaptation evolutionary strategy Chaos operator

1 Introduction

Multi-objective optimization problems (MOPs) [37], which involve multiple objectives that need to be optimized at the same time. No matter in scientific research or engineering application [1], since multiple objectives are often coupled and conflict with each other. The optimization of one objective performance may result in the performance degradation of another or more objectives. The optimization problem can be described as follows:

$min F (x) = {(f_{1} (x), . . ., f_{m} (x))}^{T}$ (1) $s . t . x \in Ω$ where Ω is the search space and x is the decision variable vector. F: Ω → R^m, where m is the number of objective functions, and R^m is the objective space. Due to the conflict between the objectives, there is no single optimal solution for MOPs, a balanced way should be chosen to find a set of compromise solutions.

With the effort of more than twenty years from researchers and practitioners in different fields, multi-objective evolutionary algorithms (MOEAs) have become a major methodology for dealing with MOPs [2, 3]. Obviously in the process of multi-objective optimization, there are two key problems to be solved: (1) the approximate solutions obtained should be as close as possible to the true optimal front; (2) the approximate solutions obtained should be evenly distributed throughout the front (i.e. balance convergence and diversity). Over the past two decades, a large number of MOEAs have been proposed. The principal methods adopted at present can be divided into three categories [4], ranging from dominance-based framework [5 –9], indicator-based framework [10 –13], and decomposition-based framework [14 –16].

Although, the dominance-based framework and the indicator-based framework are quite popular, they have certain limitations. In the dominance-based framework, where Pareto ordering is adopted to determine the optimal front of the population. One of the main drawbacks is that most of the candidate solutions generated in a limited-scale group are non-dominated, making it difficult to distinguish candidate solutions based on the dominant selection criteria while the objective function increases [17]. Therefor, it is not suitable for many-objective optimization. Although the indicator-based algorithm is not affected by the advantages of MOEAs based on MOPs, the computational cost for the calculation of the performance becomes prohibitively expensive when the number of objectives is large [18].

In recent years, an MOEA based on decomposition (MOEA/D) proposed in [14] has attracted much attention in the area of MOEAs, which decomposes MOPs into a set of simple optimization subproblems and solves them in a collaborative manner. However, decomposition-based algorithms face their own challenges. One major challenge is that the decomposition-based method needs to specify a set of weight vectors, which significantly affect the performance of the algorithm in terms of diversity. Several variants and implementations of MOEA/D employ simplex-lattice design method [19] to generate evenly distributed weight vectors. However, in this method, the population size dramatically grows as the number of objectives increase and the setting of population size is not flexible. In MOED/D-DRA [21], a weight vector generation method based on uniform random sampling paradigm [20] is adopted. The advantage of the uniform random sampling method over simplex-lattice design is that the setting of population size is flexible. Recently, Tan et al. [22] proposed a new version of MOEA/D (UMOEA/D). UMOEA/D adopts the uniform design method to set the aggregation coefficient vectors of the subproblems. Compared with MOEA/D, distribution of the coefficient vectors is more uniform over the design space, and the population size neither increases nonlinearly with the number of objectives nor considers a formulaic setting. Deb and Jain [23] proposed a reference-point-based evolutionary algorithm (NSGA-III) with two-layer weight vector generation method. The two-layer weight vector generation method overcomes the limitation of the simplex-lattice design method in generating relatively small number of evenly spread weight vectors. In another aspect, the dynamic reference point is adopted by adjusting a parameter in [25], which can strike a good balance between exploitation and exploration. In [27], the effect of increasing population diversity by calculating the angle between the objective vector and the reference vector is achieved. To deal with the irregular PF frontier, a multi-objective evolutionary algorithm (MOEA-NAA) based on neighborhood adaptive adjustment and enhanced reference vector is proposed in [52]. The effectiveness of these algorithms in solving problem with a highly irregular PF has been proved. Due to the effectiveness and simplicity of the decompose frame based on reference point, it has been widely applied in MOPs. Therefore, a decomposition framework based on reference points will be applied in this paper.

More specifically, the decomposition mechanism decomposes MOP into many separate subproblems, and solve them cooperatively. In this respect, any single-objective operators can be potentially plugged as search engine for the single-objective subproblems. For example, polynomial mutation, SBX crossover [28, 29] and differential evolution (DE) [30], where the latter is believed to be among the better performing implementations. However, the crossover rate and mutation rate are difficult to determine in the actual calculation. The global optimal solution in the convergence domain can not be guaranteed when the crossover rate and mutation rate are small, while large rate will increase the number of iterations [14]. Sometimes the rate needs to be adjusted manually, which will reduce the objectivity of the algorithm. The underlying objective of mutative strategy parameter control is roughly to favor previously selected mutation steps in the future. If this goal is strictly pursued, a completely asymmetric adaptive scheme will be produced, which can adapt to any normal variation distribution. Moreover, this inclusion of single-objective optimizers might not always be straightforward due to two additional specificities of the MOEA/D framework related to the cooperativity among the search process of the scalar problems.

To meet the previously stated demands, CMA-ES [31] was proposed, which is one of the most successful single objective evolutionary algorithms as a core single-objective evolution engine to generate new solutions by sampling a multivariate Gaussian distribution model. In [32], CMA-ES considered as a state-of-the-art optimizer for single-objective blackbox continuous problems, which can be greatly improved through accumulation-using evolutionary paths instead of a single search step. In [33], an extension of CMA-ES for multi-objective optimization (MO-CMA-ES) is proposed and studied under the framework of NSGA-II. The algorithm uses Pareto dominance and its performance index selection. Each individual solution is associated with a different Gaussian model. Due to the invariant nature of CMA-ES, MO-CMA-ES proved to be a very favorable variable when processing MOP interactively. In MOEA/DD-CMA [34], a new algorithm based on CMA-ES and MOEA/DD [35] was proposed. It uses the selection mechanism of MOEA/DD that employs dominance and decomposition simultaneously to increase the diversity of the solutions found, although with a slower convergence rate. In [36], the covariance strategy to combine with differential evolution is used to solve the deviation problem and get good results. Since CMA-ES considers the adaptive strategy for controlling step size in sampling new solutions, it is potential to exploit the optimal solutions in a small area.

Moreover, genetic algorithms have the disadvantage of slow convergence [49]. The main reason for the premature convergence of CMA-ES is that all individuals in the group have too much selection pressure during the early evolution process, which leads to a rapid convergence in the pattern. In order to prevent local premature convergence, a chaotic operator is introduced, which uses the specific inherent randomness and ergodicity of chaotic variables to jump out of the local best. In [50], a novel genetic algorithm containing chaos operator was proposed based on the analysis of population diversity and premature convergence within the framework of Markov chain. This algorithm increases the population size dynamically so as to restore the population diversity and prevent premature convergence effectively. J. Li [51] used the properties of ergodicity, stochastic property, and regularity of chaos to lead particles’ exploration to prevent the prematurity of particles.

Based on the existing studies, covariance matrix adaptive strategy in reference point-based multi-objective evolutionary algorithm is proposed in this paper, and the main improvements of RPCMA-ES are listed as follows:

(1) A decomposition framework based on adaptively adjusted reference vectors is applied. By calculating the distance between the target vector and the relevant reference vector, which can effectively select uniformly distributed particles and search the entire boundary.

(2) CMA-ES is used as a way to generate a new population, which guides individuals to search in a better direction. The application of CMA-ES strategy improves the convergence performance and accelerates the convergence speed of the algorithm.

(3) Chaos operator is added to the original CMA-ES algorithm to prevent individuals from falling into local optimum. The new distribution mean will be generated by the chaotic equation, which can guide the search of population with better probability and find the global optimization.

The remainder of this paper is organized as follows. Section 2 describes the covariance matrix adaptive mechanism, the generation of reference points and the details of RPCMA-ES. Section 3 presents the simulation results compared with other algorithms. Finally, the conclusions and some possible paths for future work are provided in Section 4.

2 Algorithmic framework

More recently, a set of direction vectors are used to decompose the entire PF into several segments, each of which is a multi-objective subproblem. This decomposition strategy has aroused some interest. For example, in NSGA-III [23], a set of reference points or reference lines is used for niche preservation to manage diversity in each subspace for many-objective optimization, which effectively enhances convergence by giving priority to solutions closer to the reference points. In RVEA [27], a reference vector-guided EA for multi-objective optimization was proposed. The reference vectors can be used not only to decompose the original multi-objective optimization problem into a number of single-objective subproblems, but also to elucidate user preferences to target a preferred subset of the whole PF. Therefore, a framework based on reference point decomposition will be used in this paper. Fig. 1 illustrates the overall flowchart of the framework based on reference point decomposition.

Fig. 1

The overall flowchart of the framework based on reference point decomposition.

2.1 Generate reference points for area division

In order to screen out the non-dominant solutions with good diversity and convergence effectively, reference points are introduced into objective space and the objective space is divided according to the reference vectors. Thus, a set of uniformly distributed reference vectors on a hyperplane need to be generated by the method proposed in [41]. ${\begin{matrix} ɛ_{i}^{j} & \in {\frac{0}{P}, \frac{1}{P}, . . ., \frac{P}{P}}, \sum_{j = 1}^{m} ɛ_{i}^{j} = 1 \\ ɛ_{i} & = (ɛ_{i}^{1}, ɛ_{i}^{2}, . . ., ɛ_{i}^{m}), i = 1, 2, . . ., N \end{matrix}$ (2) where m is the number of objective functions, i is the number of reference points, i = 1, 2, . . . , N, and P is a predefined parameter. The reference points are uniformly distributed on the hyperplane that can be mapped to the hypersphere by the Eq. (3):

$w_{i} = \frac{ɛ_{i}}{‖ ɛ_{i} ‖}$ (3)

In this method, the number of reference points is calculated by $H = (\begin{matrix} P + m - 1 \\ P \end{matrix})$ , and it can be seen that the number of reference points depends on the objective number m and integer P of the objective. Take Fig. 2 as an example, if P is set to 4 for a 3-objective problem, 15 reference points will be generated.

Fig. 2

15 reference points on the normalized hyperplane for 3-objective problem with P = 4.

If the population size is N, N reference vectors will be generated to divide the objective space into N sub-regions. Thus, the population is divided by calculating the angle of individuals in the population to each reference vector and the reference vectors of the minimum angle are correlated. As shown in Fig. 3, for each objective vector t_i, the angle between t_i and the reference vector w_i is calculated by the Eq. (4): $cos θ = \frac{t_{i} \cdot w_{i}}{‖ t_{i} ‖ ‖ w_{i} ‖}$ (4)

Fig. 3

A schematic diagram of the reference vector associated with the objective vector. t_i is a objective vector, w₁ and w₂ are the the reference vectors. θ₁ and θ₂ are the angles between t_i from w₁ and w₂.As can be seen from the Fig. 3, θ₂ is smaller than θ₁, so t_i is related to w₂.

The quality of the solution is evaluated by a series of reference points so as to assist in controlling the distribution of population in the objective space. This evaluation method will be described in detail in the Environmental Selection.

2.2 Covariance matrix adaptive mechanism

Covariance matrix adaptive is a kind of random numerical optimization algorithm that calculate gradient is not needed, which increases the probability of correct search direction by changing the covariance matrix. Note that the evolutionary strategy is different from the genetic algorithm, but both are important variants of the evolutionary algorithm.

The parameters involved in CMA-ES are mean value of distribution m^(g) ∈ Rⁿ, which also called the center of a population or expectation, step size σ^(g) > 0, covariance matrix C^(g) ∈ R^n×n and evolutionary path at generation g. The basic idea is that the probability of producing excellent individuals gradually increased by adjusting these pre-set parameters (the probability of searching in a good direction increases). As a way to generate offspring, CMA-ES is different from genetic algorithm in that crossover mutation operation is not implemented, instead, gaussian distribution in the solution space of optimization problem is adopted. Therefore, the purpose of generating λ sample individuals and evaluating the fitness of new samples in the population can be better achieved. The evolution direction of the next generation population is adjusted by using the evolution strategy parameters so as to generate the next generation of individuals.

More specifically, the detailed operations of CMA-ES are presented in Algorithm 1. As described in Algorithm 1, each subproblem generates subpopulation through the CMA-ES operation. The CMA-ES adopted in this paper follows the literature [42]. The step1 is the initialization process. Firstly, the initial parameters are set, which m^(g) is a sample randomly selected from the search space, the covariance matrix C^(g) is an identity matrix I_n×n, the step size σ^(g) is set to 0.35, and the evolutionary path vector $P_{c}^{(g)}$ , $P_{σ}^{(g)}$ are set to 0 ∈ Rⁿ. Although CMA-ES identifies the parameters involved, it is adaptive in subsequent iterations. The main aim of step2 is elitism strategy. In each subsequent iteration, λ ≥ 1 samples are taken for each individual, which are sampled in the multivariate normal variation distribution following Eq. (5):

Algorithm 1 The process of CMA-ES.

Input: parent population P_t;

Output: offspring population Q;

1: Step1: Initialization

2: for i = 1 to N do

3: mⁱ← the ith individual in P_t;

4: Initialize ${P_{c}^{i}, P_{σ}^{i}, C, m^{i}, σ^{i}, g}$ ;

5: end for

6: Step2: Update and Select

7: while stopping criterion not satisfied do

8: for j = 1 to λ do

9: $X_{j} \sim N (m_{i}, σ_{i}^{2} C_{i})$ ;

10: end for

11: Generate the average candidate distribution using (11);

12: SetΛ = X ⋃ mⁱ;

13: Λ = TruncationSelection (Λ);

14: Update Parameters{P_c, P_σ, C, m, σ, g};

15: Q = Select (Λ);

16: end while

17: return Q

$x_{i}^{(g + 1)} = m^{(g)} + σ^{(g)} N ({0, C}^{(g)}), i = 1, 2, . . ., λ$ (5)

After the sampling process is completed, the central point of the new population is updated by selecting the optimal solution from the generated population as a sub-population for recombination. Each candidate solution (i.e. the number of individuals is combined with the mean value of the extracted sample and the birth distribution) is sorted according to the value of the objective function where the annotation x_i:λ for the ith best solution(i.e. f (x_1:λ) ≤ . . . ≤ f (x_u:λ) ≤ f (x_u+1:λ) ≤ . . . ≤ f (x_λ:λ)) is used, then the first μ (μ = λ/2) individuals are selected to update distribution parameters by truncation selection.

The new average distribution is calculated by the Eq. (6): $\begin{matrix} z_{k} = (x_{k} - m^{(g)}) / σ, k = 1, 2, . . ., u \\ z_{w} = \sum_{k = 1}^{μ} (w_{k} \times z_{k}) \\ m^{(g + 1)} = m^{(g)} + σ z_{w} \end{matrix}$ (6) where x_k is the k-th beat solution among {m^(g)} ∪ {x_i ∣ i = 1, 2, . . . , λ}, w_k > 0, is weighted coefficient and $\sum_{k = 1}^{μ} w_{k} = 1$ . The result of the distribution center of the next generation population will shift to the sub-population until tangent to the optimal solution is presented by the Eq. (6). In this paper, chaos operator is used to generate the mean value of candidate distribution to find the optimal mean, which will be introduced in detail in the next section. In addition, the update of covariance matrix is the next most important step of the population, which the updating principle is to increase the variance along the successful search direction, that is, to increase the probability of sampling along these directions.

The Method of combining rank-1-update and rank-μ-update is used in update operation, which is described by the Eq. (7):

$\begin{matrix} C^{(g + 1)} = & (1 - c_{1} - c_{μ}) C^{(g)} \\ + c_{1} (\frac{1}{μ_{cov}}) (P_{c}^{(g + 1)} P_{c}^{(g + 1)} δ (h_{σ}) C^{(g)}) \\ + c_{μ} \sum_{k = 1}^{μ} w_{k} y_{k} y_{k}^{T} \end{matrix}$ (7) the first term on the right-hand side is rank-1-update [38], which can be interpreted as directly viewing the evolutionary path as a successful search direction; the second is rank-μ-update [39, 40] that is actually the weighted maximum likelihood estimation using the μ solutions selected. Rank-μ-update can also be interpreted as the geometric optimization, stochastic optimization, and evolutionary strategy for the natural gradient information of C. c₁ and c_μ are the learning rates of the above two methods respectively, and μ_eff is the variance effective selection mass. The μ_cov is the weight parameter for the rank-1-update and rank-μ-update. δ (h_σ) = (1 - h_σ) c_c (2 - c_c) ≤ 1 is inessential and can be set to 0 in the usual case, h_σ is a Heaviside function.

Covariance requires fully excavating a series of relevant information between evolutionary generations to be updated. Hence, the concept of evolutionary path is introduced in [31, 38]. Evolutionary path P_c is the sum of the variation steps of continuous evolutionary generation, which is updated by the Eq. (8): $P_{c}^{(g + 1)} = (1 - c_{c}) P_{c}^{(g)} + h_{σ} \sqrt{c_{c} (2 - c_{c}) μ_{eff}} z_{w}$ (8) where c_c is the backward time horizon of the path P_c. The Heaviside function can stall the update of P_c, which is defined as h_σ = 1 if $\frac{‖ P_{σ} ‖}{\sqrt{1 - {(1 - c_{σ})}^{2 (g + 1)}}} < (1.4 + \frac{2}{n + 1}) E ‖ N (0, I) ‖$ or 0 otherwise. The fast increase of axis of C in a linear surrounding (when the step size is far too small) can be prevented by P_c.

In the updating of covariance matrix, only the variation direction is controlled by P_c, which makes the overall scaling efficiency is not high. In order to control the step size, another evolution path P_σ (i.e. a sum of successive steps) is introduced. This control method is independent of the covariance matrix, which is updated by Eq. (9): $P_{σ}^{(g + 1)} = (1 - c_{σ}) P_{σ}^{(g)} + \sqrt{c_{σ} (2 - c_{σ}) μ_{eff}} C^{(g) - \frac{1}{2}} z_{w}$ (9) where the learning rate for the cumulation for the step-size control is $c_{σ} = \frac{μ_{eff} + 2}{n + μ_{eff} + 3}, and \frac{1}{c_{σ}}$ is the backward time horizon of the evolution path. The transformation $C^{(g) - \frac{1}{2}}$ makes the expected length of $P_{σ}^{(g)}$ independent of its direction.

The evolution path P_σ is calculated to control the step size σ^(g) by the Eq. (10): $σ^{(g + 1)} = σ^{(g)} exp (\frac{c_{σ}}{d_{σ}} (\frac{‖ P_{σ}^{(g)} ‖}{E ‖ N (0, I) ‖}))$ (10) where $E ‖ N (0, I) ‖ = \sqrt{n} (1 - \frac{1}{4 n} + \frac{1}{21 n^{2}})$ is the expectation Euclidean norm of a distributed random vector, and d_σ is a damping parameter.

Most of the parameters in the above formula are self-contained. When the birth parameter is set, the algorithm will iterate and loop in accordance with the above operations and gradually search for the global optimal solution until the stop condition is reached.

Termination conditions play an important role in the optimization of restart subproblems and the reduction of resource waste in CMA-ES. The part of the descriptions in the literature are applied [42], which are described as follows:

NoEffectCoord. Reset if adding 0.2-standard deviations in any single coordinate does not change m^(g) (i.e. m^(g) equals m^(g) + 0.2σ^(g)c_j,j for any j = 1, . . . , n).

NoEffectAxis. Reset if adding a 0.1-standard deviation vector in any principal axis direction of C^(g) does not change m^(g). More formally, stop if m^(g) equals $m^{(g)} + 0.1 σ^{(g)} \sqrt{d_{jj}} b_{j}$ , where j = (g modn) +1 and d_jjand b_j are the j^th eigenvalue and eigenvector of C^(g), with ||b_j||=1, respectively.

TolXUp. Stop/Reset if σ^(g) · max (diag (D)) increased by more than 10⁴, where D ∈ Rⁿ, is a diagonal matrix which the diagonal elements are square roots of eigenvalues of C^(g) and correspond to the respective columns of eigenvectors of C^(g).

ConditionCov. Stop/Reset if the condition number of the covariance matrix exceeds 10¹⁴.

Therefore, it means that the current optimization of the problem is over and the next sub-optimization process is started as long as any of the above conditions are met.

2.3 Chaos operator

Most CMA-ES algorithms iterate according to Eq. (5)-(10) and search for the optimal solution. In order to make the new average distribution mean more likely to approach the optimal solution and avoid falling into the local optimal solution, chaos operator [43] is added on the basis of CMA-ES.

The chaos operator takes the current mean as the initial point after each update of the distribution mean, and uses the Eq. (11) to generate the candidate mean: $m_{a}^{(n + 1)} = α m_{a}^{(n)} (1 - m_{a}^{(n)})$ (11) where m_a is the distribution mean in the candidate sequence. The initial m_a expectation is the mean updated by the candidate solution for the first time after each iteration. n represents the n-th candidate mean, and α is the chaos control parameter. Moreover, the current distribution mean is replaced by the distribution mean in the candidate sequence.

The new generation of distribution mean can not only ensure the population to search in a better direction, but also jump out of the local optimum effectively, which improves the convergence of the population.

2.4 Environmental selection

After the generation of offspring population, the current population size is 2N, and it is necessary for an effective selection mechanism to select excellent N individuals from 2N individuals as the population of the next iteration. The process uses the reference vectors to divide the objective space into subspaces with the same number as the initial population size. The best individuals will be chosen in these subspaces. The process is described in Algorithm 2

Algorithm 2 RP-Select Strategy.

Input: the number of iterations (g), population N, the current population (P), the set of reference vectors (V);

Output: population P_t+1;

1: for i = 1 to |P| do

2: R = Normalization (P);

3; end for

4: Partition the Population: calculate angle between x_i ∈ R and v_j ∈ V using (16),associate all individuals with reference vectors;

5: Calculate the Angle-Distance of each reference vector and the individuals associated with it which will operate independently in each of the sub-regions by using (17);

6: Select the closest individual to each reference vector through AD;

7: return P_t+1

The main framework of the environmental selection is divided into five parts. First of all, since the non-dominant solutions may be different in the range of each objective value, scaling operation is needed to convert the objective value to facilitate the next association operation, the initial population can be standardized as

$f_{i}^{'} (x) = \frac{f_{i} (x) - z_{i}^{min}}{z_{i}^{max} - z_{i}^{min}}, i = 1, 2, . . ., m$ (12) where the ideal point $Z^{min} = {(z_{1}^{min}, z_{2}^{min}, . . ., z_{m}^{min})}^{T}$ is the minimum value of each objective dimension i in M objective functions for all individuals in P and the nadir point $Z^{max} = {(z_{1}^{max}, z_{2}^{max}, . . ., z_{m}^{max})}^{T}$ can be constructed by finding the maximum value for each objective. It is convenient for the next step to use the normalization process of the reference vector and the individual so as to make the ideal point becomes the origin.

The angles between individuals x_i ∈ R and each reference vector v_j ∈ V are calculated as

$θ (x_{i}, v_{j}) ≜ arccos | \frac{f_{i}^{'} (x) \cdot v_{j}}{‖ f_{i}^{'} (x) ‖} |$ (13) then the individuals with the smallest angle (i.e. largest cosine value) of each reference vector are found to be correlated (lines 6-13 in Algorithm 2).

The objective space is divided into N independent subspaces by reference vectors. Through the numerical ordering of Eq. (4), the individuals of the population can be determined which reference vectors to associate with.

Eq. (4) calculates the cosine value between the objective vector and all reference vector. By cosine value sorting, the reference vector with the largest cosine value (i.e., the smallest angle) of the objective vector can be determined. The reference vector can be considered to be associated with the objective vector. Traverse all reference vectors associated with individuals, then the method mentioned in literature [27] will be used to calculate the angle-distance (AD) which be explained by the Eq. (14):

$D_{i, j} = (1 + P (θ_{i, j})) \cdot ‖ f_{i}^{'} (x) ‖$ (14) the penalty function P (θ_i,j) is associated with g and θ_i,j, α is a user-defined parameter that controls the rate of change. The individual with the smallest AD value is the optimal solution in this sub-region.

Finally, Z^min and Z^max are recalculated to adjust the reference vector adaptively by the Eq. (15): $v_{i}^{'} = \frac{v_{i} \cdot (Z^{max} - Z^{min})}{‖ v_{i} \cdot (Z^{max} - Z^{min}) ‖}$ (15)

The basic idea is use Z^min and Z^max to zoom in or out according to the shape of the leading edge so as to change the position of the reference vector for guidance and provide a better search direction for the next iteration.

2.5 RPCMA-ES

The pseudocode of the method is shown by Algorithm 3. The strengths of the selection mechanism according to the position of reference points and population individuals are used in RPCMA-ES. The individuals with good convergence and diversity are selected to guide the population to search in a better direction. It can be seen that the proposed algorithm includes initialization, CMA-ES evolution process, selection of reference vectors and adaptive adjustment of reference vectors. First, randomly generate an initial population P of size N. Second, a set of uniformly distributed reference vectors V (|V| = |P|) are generated to evaluate the diversity of these solutions. It is well-known that a set of uniformly distributed reference vectors can divide the objective space into several independent sub-regions. Next, covariance adaptive adjustment strategy is adopted to ensure the generation of more excellent sub-population, which is the main contribution of this paper. Finally, the optimal individual is selected from the combined population according to the relationship between the objective value vector and the reference vector.

Algorithm 3 General Framework of the RPCMA-ES.

Input: stopping criterion MaxGen, the number of population N, the number of objectives m;

Output: Pareto-optimal solution P;

1: Generate a random set of solutions: P = {x₁, x₂, . . . , x_N};

2: Generate a set of uniform reference vectors from reference points: V = Generate (N, m);

3: while gen < MaxGen do

4: Q = CMA - ES (P);

5: R = P ⋃ Q;

6: P = RP - Select (gen, R, V);

7: V = AdaptiveAdjustment (V);

8: gen = gen + 1;

9: end while

10: return P

3 Experimental results

Aim at testing the effectiveness of the proposed algorithm, four popular algorithms are selected for comparison including RVEA, NSGA-II, R2MOPSO, and MOEA/D-STM. In this section, an introduction to the test function and performance index is given. Then, the parameter settings of each algorithm are described, and the experimental results are analyzed finally.

3.1 Test problems

To test the performance of the algorithm which is proposed in this paper, WFG1-WFG10 benchmarks are used. The PF of WFG test problems has the characteristics of linear, convex, concave, mixed, multi-modal or discontinuous [44]. Besides, any number of objectives and decision variables can be scaled to in WFG test problems. They present a significant challenge for an evolutionary algorithm to find a well-converged and well-distributed solution set. Moreover, the true Pareto optimal front is known for these problems. The number of decision variables is set to D = k + l for WFG problem instances, where the position-related variable k = 18 and the distance-related variable l = 14. In this work, 2-objective and 3-objective of these problems are focused.

3.2 Performance indicators

In this paper, the performance indicators of inverted generational distance (IGD) [45, 46] and HV are widely used to evaluate the performance of the algorithm, on account of no one indicator can assess a MOEA [47, 48] comprehensively.

(1) IGD indicator [26]: IGD is the comprehensive performance evaluation index of a algorithm, reflecting the distribution and convergence of a algorithm. Let S be a set of approximate solutions in the objective space, F is a solution set uniformly distributed on the real PF. The PF contains 1,000 points which is knew from the start. The IGD value is defined as:

$IGD (S, F) = \frac{1}{| F |} \sum_{λ \in F} dist (λ, S)$ (16) where dist (λ, S) is the the Euclidean distance between the individual λ ∈ F and its nearest point in S, |F| is the cardinality of F. The good performance of the algorithm can be reflected by a small IGD value.

(2) HV indicator [24]: The hypervolume measures the volume of dimensional regions in the objective space enclosed by the non-dominant solutions set and reference points obtained by the multi-objective optimization algorithm. The mathematical expression of hypervolume is as follows: $HV (S) = V (⋃_{μ \in S} [f_{1} (μ), z_{1}] \times . . . \times [f_{m} (μ), z_{m}])$ (17) where set S is the obtained non-dominated solutions, V (·) is the Lebesgue measure which can be measured volume. Let z = (z₁, . . . , z₂) ^T be the reference point in the objective space. In this experiment, z is used (3, 5) ^T for 2-objective test problems, and (3, 5, 7) ^T for 3-objective test instances. HV is an effective unitary quality measurement index, which is strictly monotonous in Pareto dominance. The higher the value of HV is, the better the performance of the corresponding algorithm will be. In addition, the calculation of HV index is not needed to test the ideal PF of the problem, which greatly facilitates the use of HV in practical applications.

3.3 Parameter settings

The parameters of the proposed algorithm are set as follows: the population size is set to 50 for 2-objective and 105 for 3-objective, which is determined by the objective number m. For WFG1-WFG3, WFG8 and WFG10, the number of generations is set to 1000. The number of generations is set to 800 for WFG4, is set to 500 for WFG5, WFG7, and is set to 300 for WFG6, WFG9.

The parameter settings for the CMA-ES strategy are shown in Table 1. λ represents the number of samples and λ ≥ 2. μ is the number of selected search points in the population, which is set to λ/2 by truncation selection. ω_k is positive weight coefficients for recombination, Eq.(6) calculates the mean value of µ selected points. c₁ and c_μ are the learning rates of the rank-one and rank-µ update respectively. μ_eff is the variance effective selection mass. μ_cov is the parameter for weighting between rank-one and rank-µ update, which determines their relative weighting and is chosen as μ_eff is most appropriate. d_σ is damping parameter for step-size update, which scales the change magnitude of step size. c_c and c_σ are the backward time horizon of the path P_c and the learning rate for cumulating for the step-size control respectively.

Table 1
Parameter settings for CMA-ES

λ = 4 +3 ln(n) μ = (λ/2)

$ω_{k} = log ((λ + 1) / 2) - Σ_{k = 1}^{μ} log (k), k = 1, . . . μ$

$c_{1} = \frac{2}{(n + 1.3)^{2} + μ_{eff}}$

$c_{μ} \approx \min (1 - c_{1}, \frac{μ_{eff} - 2 + 1 / μ_{eff}}{(n + 2)^{2} + μ_{eff}})$

$μ_{eff} = \frac{1}{Σ_{k = 1}^{μ} ω_{k}^{2}}$ μ_cov ≈ μ_eff

$d_{σ} = 1 + 2 \max (0, \sqrt{\frac{μ_{eff} - 1}{n + 1}} - 1) + c_{σ}$

$c_{σ} = \frac{μ_{eff} + 2}{n + μ_{eff} + 3}$ $c_{c} = \frac{4 + μ_{eff} / n}{n + 4 + 2 * μ_{eff} / n}$

In order to be fair, the number of iterations and population size in the comparison algorithm are consistent with the RPCMA-ES. For the crossover and mutation operator, the distribution index is set to η_c = 20. In addition, the crossover probability p_c = 1.0 is used in the compared algorithms except R2MOPSO. For the polynomial mutation, the distribution index and the mutation probability are set to η_m = 20 and p_m = 1/n, respectively in the four algorithms. The neighborhood size T is set to 20 in decomposition framework, and the particle of agemax T_a is set to 2 in R2MOPSO. For each test function, all algorithms run 20 times independently.

3.4 Significant test

To verify the statistical difference between RPCMA-ES and other compared algorithms, an useful statistical test is recommended to be conducted. The Wilcoxon signed-rank test can be used to make simple pairwise comparisons. Therefore, this statistical test is used herein.

Wilcoxon signed-rank test is a non-parametric method used to test whether there is a significant difference in the distribution of the population from which two pairs of samples come from. First, the observations of the corresponding samples of the first group are subtracted from the observations of the second group of samples, respectively. If the difference is positive, it is recorded as ^′ + ′; if it is negative, it is recorded as ^′ - ′, and the difference data is also saved. Subsequently, the difference data are sorted in ascending order, and the rank of the difference is obtained. Next, the positive sign rank sum W⁺ and the negative sign rank sum W^- are calculated in steps.

Eventually, if the p-value does not meet the acceptable criterion, then, the null hypothesis could not be held. In this paper, the significance level is set as 0.05. When p-value is greater than 0.05, it means that the compared HV mean values or IGD mean values are similar in statistics.

As can be seen from Tables 2-5, the oblique bold of p-value in tables shows that there is no significant difference between RPCMA-ES and other algorithms. For WFG2, there is no significant difference between RPCMA-ES and NSGAII in the HV values of 2-objective, 3-objective, and IGD values of 3-objective; similarly, there is no significant difference from RVEA in HV values of 3-objective Sexual differences. For WGF9, RPCMA-ES is not significantly different from NSGAII and RVEA in the index values of the 2-objective and 3-objective, respectively. In addition, there is no significant difference between the algorithm of this paper and NSGAII in the IGD value of WFG4 with a 3-objective. In addition to the above-mentioned relationship with no significant difference, there are obvious differences in the remaining test functions. These tables clearly indicate that RPCMA-ES received more significant in most test functions than other comparison algorithms.

Table 2
HV values obtained by RPCMA-ES, RVEA, NSGA-II, R2MOPSO and MOEA/D-STM on WFG1-WFG10 with 2-objective

Problem RPCMA-ES RVEA NSGA-II R2MOPSO MOEA/D-STM

WFG1 mean 2.715E-01 3.841E-01 4.423E-01 8.317E-01 1.664E-01

std 3.100E-02 2.030E-02 1.400E-02 2.300E-03 2.520E-02

p 6.796E-08 6.796E-08 6.796E-08 2.960E-07

WFG2 mean 6.128E-01 5.995E-01 6.127E-01 4.899E-01 7.070E-02

std 2.000E-03 1.780E-02 9.500E-03 6.700E-03 1.000E-03

p 1.375E-06 9.204E-02 6.786E-08 6.748E-08

WFG3 mean 7.076E-01 6.494E-01 7.074E-01 4.497E-01 9.670E-02

std 9.600E-03 1.670E-02 1.220E-02 9.900E-03 2.700E-03

p 6.796E-08 4.922E-02 6.796E-08 6.796E-08

WFG4 mean 4.626E-01 4.893E-01 5.095E-01 3.669E-01 7.550E-02

std 2.250E-02 1.290E-02 1.790E-02 6.600E-03 2.900E-03

p 6.015E-07 3.416E-07 6.796E-08 6.786E-08

WFG5 mean 5.013E-01 4.768E-01 5.193E-01 3.655E-01 1.999E-01

std 1.130E-02 9.400E-03 1.490E-02 7.600E-03 3.800E-03

p 2.352E-06 1.600E-05 6.796E-08 6.796E-08

WFG6 mean 5.144E-01 4.619E-01 5.019E-01 4.828E-01 1.774E-01

std 6.500E-03 1.420E-02 1.600E-02 1.730E-02 2.800E-03

p 6.796E-08 5.091E-04 6.796E-08 6.796E-08

WFG7 mean 4.844E-01 4.136E-01 4.567E-01 2.872E-01 4.270E-02

std 6.600E-03 1.580E-02 1.470E-02 9.200E-03 5.500E-03

p 6.796E-08 5.166E-06 6.796E-08 6.796E-08

WFG8 mean 4.108E-01 3.115E-01 4.102E-01 2.599E-01 1.281E-01

std 1.340E-02 1.170E-02 1.530E-02 5.900E-03 6.900E-03

p 6.796E-08 4.672E-02 6.796E-08 6.796E-08

WFG9 mean 4.873E-01 4.741E-01 4.848E-01 4.396E-01 1.664E-01

std 1.410E-02 6.800E-03 1.090E-02 7.700E-03 7.700E-03

p 8.292E-05 3.942E-01 6.796E-08 6.796E-08

WFG10 mean 8.149E-01 8.672E-01 8.642E-01 7.196E-01 2.965E-01

std 1.230E-02 6.500E-03 4.700E-03 8.100E-03 9.000E-03

p 6.796E-08 6.796E-08 6.796E-08 6.796E-08

+/=/- 7/0/3 6/0/4 9/0/1 10/0/0

Problem		RPCMA-ES	RVEA	NSGA-II	R2MOPSO	MOEA/D-STM
WFG1	mean	2.715E-01	3.841E-01	4.423E-01	8.317E-01	1.664E-01
	std	3.100E-02	2.030E-02	1.400E-02	2.300E-03	2.520E-02
	p		6.796E-08	6.796E-08	6.796E-08	2.960E-07
WFG2	mean	6.128E-01	5.995E-01	6.127E-01	4.899E-01	7.070E-02
	std	2.000E-03	1.780E-02	9.500E-03	6.700E-03	1.000E-03
	p		1.375E-06	9.204E-02	6.786E-08	6.748E-08
WFG3	mean	7.076E-01	6.494E-01	7.074E-01	4.497E-01	9.670E-02
	std	9.600E-03	1.670E-02	1.220E-02	9.900E-03	2.700E-03
	p		6.796E-08	4.922E-02	6.796E-08	6.796E-08
WFG4	mean	4.626E-01	4.893E-01	5.095E-01	3.669E-01	7.550E-02
	std	2.250E-02	1.290E-02	1.790E-02	6.600E-03	2.900E-03
	p		6.015E-07	3.416E-07	6.796E-08	6.786E-08
WFG5	mean	5.013E-01	4.768E-01	5.193E-01	3.655E-01	1.999E-01
	std	1.130E-02	9.400E-03	1.490E-02	7.600E-03	3.800E-03
	p		2.352E-06	1.600E-05	6.796E-08	6.796E-08
WFG6	mean	5.144E-01	4.619E-01	5.019E-01	4.828E-01	1.774E-01
	std	6.500E-03	1.420E-02	1.600E-02	1.730E-02	2.800E-03
	p		6.796E-08	5.091E-04	6.796E-08	6.796E-08
WFG7	mean	4.844E-01	4.136E-01	4.567E-01	2.872E-01	4.270E-02
	std	6.600E-03	1.580E-02	1.470E-02	9.200E-03	5.500E-03
	p		6.796E-08	5.166E-06	6.796E-08	6.796E-08
WFG8	mean	4.108E-01	3.115E-01	4.102E-01	2.599E-01	1.281E-01
	std	1.340E-02	1.170E-02	1.530E-02	5.900E-03	6.900E-03
	p		6.796E-08	4.672E-02	6.796E-08	6.796E-08
WFG9	mean	4.873E-01	4.741E-01	4.848E-01	4.396E-01	1.664E-01
	std	1.410E-02	6.800E-03	1.090E-02	7.700E-03	7.700E-03
	p		8.292E-05	3.942E-01	6.796E-08	6.796E-08
WFG10	mean	8.149E-01	8.672E-01	8.642E-01	7.196E-01	2.965E-01
	std	1.230E-02	6.500E-03	4.700E-03	8.100E-03	9.000E-03
	p		6.796E-08	6.796E-08	6.796E-08	6.796E-08
+/=/-			7/0/3	6/0/4	9/0/1	10/0/0

Table 3

HV values obtained by RPCMA-ES, RVEA, NSGA-II, R2MOPSO and MOEA/D-STM on WFG1-WFG10 with 3-objective

Problem		RPCMA-ES	RVEA	NSGA-II	R2MOPSO	MOEA/D-STM
WFG1	mean	4.089E-01	4.982E-01	4.494E-01	4.193E-01	3.024E-01
	std	1.460E-02	1.130E-02	6.000E-03	3.800E-03	4.000E-02
	p		6.786E-08	6.796E-08	4.300E-03	6.796E-08
WFG2	mean	7.322E-01	7.303E-01	6.567E-01	6.276E-01	5.120E-02
	std	5.740E-02	6.270E-02	3.410E-02	2.510E-02	1.000E-03
	p		9.620E-02	2.977E-01	1.431E-07	6.786E-08
WFG3	mean	8.041E-01	7.894E-01	7.856E-01	5.189E-01	1.057E-01
	std	1.350E-02	1.840E-02	7.300E-03	8.900E-03	2.800E-03
	p		6.674E-06	7.948E-07	6.796E-08	6.796E-08
WFG4	mean	5.773E-01	6.828E-01	6.078E-01	4.260E-01	9.780E-02
	std	1.910E-02	2.000E-02	2.250E-02	1.180E-02	4.400E-03
	p		6.796E-08	1.105E-05	6.796E-08	6.796E-08
WFG5	mean	5.754E-01	6.393E-01	6.186E-01	3.838E-01	2.069E-01
	std	2.510E-02	2.030E-02	1.680E-02	7.500E-03	5.300E-03
	p		6.796E-08	7.898E-08	6.796E-08	6.796E-08
WFG6	mean	6.281E-01	6.141E-01	5.848E-01	5.798E-01	1.643E-01
	std	2.470E-02	1.870E-02	2.310E-02	9.800E-03	3.800E-03
	p		2.100E-03	2.563E-07	6.796E-08	6.796E-08
WFG7	mean	5.667E-01	5.287E-01	6.132E-01	3.491E-01	7.190E-02
	std	2.080E-02	1.900E-02	1.500E-02	1.240E-02	9.100E-03
	p		5.896E-05	1.235E-07	6.796E-08	6.796E-08
WFG8	mean	5.366E-01	4.430E-01	5.059E-01	2.992E-01	9.390E-02
	std	2.090E-02	1.970E-02	2.730E-02	6.600E-03	4.300E-03
	p		6.796E-08	1.997E-04	6.796E-08	6.796E-08
WFG9	mean	5.734E-01	5.671E-01	5.427E-01	4.955E-01	1.585E-01
	std	8.500E-03	1.490E-02	1.540E-02	8.900E-03	1.090E-02
	p		6.390E-02	1.429E-07	6.796E-08	6.796E-08
WFG10	mean	9.722E-01	9.387E-01	9.570E-01	7.976E-01	3.999E-01
	std	8.000E-03	6.600E-03	8.900E-03	6.500E-03	4.800E-03
	p		6.796E-08	6.796E-08	6.796E-08	6.796E-08
+/=/-			7/0/3	7/0/3	9/0/1	10/0/0

Table 4

IGD values obtained by RPCMA-ES, RVEA, NSGA-II, R2MOPSO and MOEA/D-STM on WFG1-WFG10 with 2-objective

Problem		RPCMA-ES	RVEA	NSGA-II	R2MOPSO	MOEA/D-STM
WFG1	mean	1.418E+00	1.058E+00	1.286E+00	1.285E+00	1.695E+00
	std	8.710E-02	4.690E-02	4.070E-02	2.600E-03	1.238E-01
	p		6.796E-08	6.796E-08	6.796E-08	1.201E-06
WFG2	mean	6.128E-01	6.358E-01	9.531E-01	7.266E-01	2.978E+00
	std	2.400E-03	8.240E-02	6.700E-03	3.280E-02	2.900E-03
	p		1.201E-06	1.440E-02	6.796E-08	6.796E-08
WFG3	mean	3.570E-02	1.719E-01	7.763E-01	6.635E-01	2.451E+00
	std	4.500E-03	3.900E-02	1.110E-02	3.060E-02	1.440E-02
	p		6.796E-08	3.939E-07	6.796E-08	6.796E-08
WFG4	mean	2.141E-01	1.339E-01	9.464E-01	5.010E-01	1.722E+00
	std	5.650E-02	2.730E-02	3.450E-02	2.190E-02	2.910E-02
	p		2.960E-07	3.417E-07	6.796E-08	6.796E-08
WFG5	mean	1.278E-01	1.737E-01	9.833E-01	6.025E-01	1.528E+00
	std	1.350E-02	2.030E-02	6.200E-03	4.080E-02	1.910E-02
	p		9.748E-06	9.173E-08	6.796E-08	6.796E-08
WFG6	mean	1.120E-01	2.257E-01	1.034E+00	2.010E-01	1.677E+00
	std	1.240E-02	3.790E-02	2.600E-02	1.140E-02	1.380E-02
	p		6.796E-08	1.230E-02	6.796E-08	6.796E-08
WFG7	mean	1.792E-01	4.070E-01	8.613E-01	8.356E-01	1.862E+00
	std	1.440E-02	6.800E-02	5.300E-02	5.020E-02	9.380E-02
	p		6.796E-08	1.807E-05	6.796E-08	6.796E-08
WFG8	mean	3.448E-01	6.719E-01	6.729E-01	8.703E-01	1.963E+00
	std	3.370E-02	7.310E-02	1.489E-01	5.390E-02	2.540E-02
	p		6.796E-08	1.807E-05	6.796E-08	6.796E-08
WFG9	mean	1.733E-01	2.023E-01	1.195E+00	3.006E-01	1.613E+00
	std	3.810E-02	1.350E-02	2.810E-02	2.300E-02	5.690E-02
	p		7.406E-05	1.475E-01	6.796E-08	6.796E-08
+/=/-			8/0/1	8/0/1	8/0/1	9/0/0

Table 5

IGD values obtained by RPCMA-ES, RVEA, NSGA-II, R2MOPSO and MOEA/D-STM on WFG1-WFG10 with 3-objective

Problem		RPCMA-ES	RVEA	NSGA-II	R2MOPSO	MOEA/D-STM
WFG1	mean	1.491E+00	1.259E+00	1.732E+00	1.512E+00	1.766E+00
	std	1.690E-02	3.040E-02	1.680E-02	3.000E-03	7.460E-02
	p		6.768E-08	2.563E-07	3.293E-05	6.796E-08
WFG2	mean	5.602E-01	5.030E-01	1.701E+00	8.454E-01	3.832E+00
	std	2.229E-01	1.799E-01	1.544E-01	3.620E-02	5.200E-03
	p		4.570E-01	3.942E-01	3.500E-06	6.796E-08
WFG3	mean	1.688E-01	1.994E-01	1.243E+00	7.754E-01	3.089E+00
	std	1.320E-02	3.260E-02	1.880E-02	2.230E-02	2.530E-02
	p		1.000E-03	3.416E-07	6.796E-08	6.796E-08
WFG4	mean	3.818E-01	2.172E-01	1.371E+00	7.638E-01	2.309E+00
	std	4.600E-02	4.600E-03	2.330E-02	2.400E-02	6.870E-02
	p		6.796E-08	7.353E-01	6.796E-08	6.796E-08
WFG5	mean	3.816E-01	2.522E-01	1.192E+00	1.159E+00	2.843E+00
	std	3.580E-02	1.690E-02	1.300E-02	7.040E-02	4.020E-02
	p		6.796E-08	3.705E-05	6.796E-08	6.796E-08
WFG6	mean	2.918E-01	3.007E-01	1.195E+00	4.414E-01	3.200E+00
	std	1.700E-02	2.100E-02	2.500E-02	2.420E-02	3.820E-02
	p		2.256E-02	6.796E-08	6.796E-08	6.796E-08
WFG7	mean	4.813E-01	5.362E-01	1.073E+00	1.210E+00	2.305E+00
	std	4.470E-02	8.610E-02	1.720E-02	5.700E-02	8.870E-02
	p		4.680E-02	1.294E-04	6.796E-08	6.796E-08
WFG8	mean	5.186E-01	6.845E-01	1.081E+00	1.123E+00	3.743E+00
	std	3.100E-02	7.790E-02	7.380E-02	2.770E-02	2.210E-02
	p		9.173E-08	1.657E-07	6.796E-08	6.796E-08
WFG9	mean	3.585E-01	3.691E-01	1.444E+00	6.086E-01	3.034E+00
	std	1.780E-02	2.210E-02	3.060E-02	3.290E-02	1.028E-01
	p		3.366E-02	6.796E-08	6.796E-08	6.796E-08
+/=/-			5/0/4	8/0/1	9/0/0	9/0/0

3.5 Results

Aim at evaluating the performance of the proposed algorithm, the other four compared algorithms are RVEA, NSGA-II, R2MOPSO and MOEA/D-STM.

Table 2 and Table 3 summarize the mean, standard deviation (std) and p-value (p) of the HV performance indicators for 2-objective and 3-objective problems. Similarly, Table 4 and Table 5 record the mean, std and p of the IGD performance indicators for 2-objective and 3-objective problems. In these tables, the best results of the mean or std for each test function are highlighted. The better/similar/worse (+/=/-) show that RPCMA-ES is better than, similar to, and worse than the performance of these compared algorithms.

As can be seen from the tables and figures, RPCMA-ES is significantly superior to R2MOPSO and MOEA/D-STM in the mean value of all test functions. It may be that the mechanism in the particle swarm optimization algorithm is not suitable for dealing with test functions like WFG with complex fronts. MOEA/D-STM is the worst performer among all the comparison algorithms. The lack of fairness of the stable matching strategy may lead to the reduction of convergence and diversity of the algorithm in dealing with complex problems.

Compared with NSGAII, Table 2 shows that the HV values of NSGAII are slightly better than that of RPCMA-ES in the two-objective problems of WFG1, WFG4, and WFG5. The reason is that NSGA-II is a typical algorithm for solving MOPs which the non-dominant solutions are selected according to the crowding distance. However, The selection pressure of NSGAII decreases with the increase of objective dimension. It can be seen from Tables 3 and 5 that the HV and IGD values of NSGAII are worse than that of RPCMA-ES. It can be clearly seen from Fig. 7, the uniformity and diversity of NSGAII are far inferior to that of RPCMA-ES. RPCMA-ES performed well mainly because it adopted the mechanism of reference point in the objective space, which can effectively screen out candidate solutions with good diversity and convergence according to the angle of non-dominant solution and reference point to reduce complexity.

Figs. 4-7 show the trend of performance indicators of five algorithms on WFG3 and WFG5-WFG9 test functions. It can be seen that, RVEA performs slightly better on the 3-objective WFG5 problem. Table 2 and Table 4 also show that the performance of RVEA is slightly better on the 3-objective problems of WFG1, WFG4, and WFG5. It may be due to WFG1, WFG4, and WFG5 are designed with different complexity criteria (e.g., WFG4 is a multimodal problem, WFG5 is deceptive, WFG1 is a problem with complex structure of the PF). The ability of RPCMA-ES to deal with such problems needs to be improved. However, whether it is HV value or IGD value, RPCMA-ES is the most effective algorithm in terms of the number of optimal outcome. Although both RPCMA-ES and RVEA adopted the strategy of reference point, the effect of RVEA on the optimization of 2-objective or 3-objective is not very good. The main reason for this result is that CMA-ES strategy combines the reliability of ES, the overall situation and the efficient guidance of adaptive covariance matrix. The adaptability is strong in CMA-ES to solve non-convex optimization problems and makes up for the shortcomings of more human factors in crossover and mutation operations in RVEA. The introduction of chaotic operator enables the population to jump out of local optimum.

Fig. 4

HV indicator value on WFG3,WFG5-WFG9 with 2-objective.

Fig. 5

IGD indicator value on WFG3,WFG5-WFG9 with 2-objective.

Fig. 6

HV indicator value on WFG3,WFG5-WFG9 with 3-objective.

Fig. 7

IGD indicator value on WFG3,WFG5-WFG9 with 3-objective.

Fig. 8 and Fig. 9 describe candidate solutions achieved by each compared algorithm on WFG3 and WFG5-WFG9 in the run resulting the best IGD among 20 runs of 2-objective and 3-objective test functions, respectively. It may be known from Fig. 8 that MOEA/D-STM performs poorly on convergence and diversity in the 2-objective, especially on the WFG3 and WFG8 problems. Besides, in addition to the better convergence on WFG9, R2MOPSO has poor convergence on other test problems. At the same time, RVEA does not perform as well as RPCMA-ES on WFG3 and WFG6-WFG8. Although the convergence of RVEA is good on WFG5 and WFG9, the search for boundary solutions is inferior to that of RPCMA-ES. By comparison, it is found that the approximate solutions of RPCMA-ES algorithm can cover and converge to the real front in a wider range. An interesting observation is that, although the convergence of the approximate PF obtained by RPCMA-ES looks slightly worse, the distribution of these solutions is better than other comparison algorithms, and thus the results obtained by the RPCMA-ES are very encouraging.

Fig. 8

The solutions set of each compared algorithm on the PF of WFG3,WFG5-WFG9 with 2-objective.

Fig. 9

The solutions set of each compared algorithm on the PF of WFG5, WFG6, WFG7, WFG8 with 3-objective.

As can be seen from Fig. 9, performance of NSGAII in the 3-objective is worse than that in the 2-objective. The reason for this may be the reduction of selection pressure. Whether in terms of convergence or diversity, R2MOPSO and MOEA/D-STM obtain relatively poor performance compared with other algorithms by the results demonstrated. Similarly, the ability to search the boundary in the 3-objective test functions of RVEA is still worse than RPCMA-ES. The RPCMA-ES algorithm performs very well in solving many problems, especially in the continuous objective function. Not only the convergence speed is fast, but also the optimal result is better than other algorithms, because in the experiment, a chaotic mutation operator is introduced into the mutation operation of the better individual. The results show that RPCMA-ES present a clear advantage over the other algorithms on the majority of the test instances.

4 Conclusion

In this paper, an adaptive covariance strategy based on reference points is proposed to solve multi-objective optimization problems which named RPCMA-ES. This algorithm is mainly composed of the adaptive covariance strategy of the probability distribution model based on the potential candidate solutions to generate new individuals by random sampling. The difference with traditional CMA-ES is that chaos operator is used to updating the distribution mean value to search for good individuals with higher probability, which enhances the convergence speed and reliability of the population. In the solution selection mechanism, the entire search space is divided into independent sub-regions by reference points, and the AD value is calculated to effectively screen out diverse candidate solutions. In this paper, several test functions are simulated and compared with the current popular algorithms. In fact, although RPCMA-ES has achieved good results in 2-objective and 3-objective test functions, the algorithm has many challenges. Next, we will pay more attention to the high-dimensional multi-objective evolutionary algorithm with more objective functions.

Footnotes

Acknowledgement

This work was supported by the Natural Science Foundation of Hebei Province (No. F2016203249) and the Youth Fund of Hebei Natural Science Foundation (No. E2018203162).

References

, Wei

, Sun

, Yang

and Wei

, Optimization of Metal Rolling Control Using Soft Computing Approaches: A Review, Archives of Computational Methods in Engineering (2019), 1–17.

Lamont

G.B.

, Evolutionary Algorithms for Solving Multiobjective Problems, Springer US (2017).

Zhou

, Qu

B.Y.

, Li

, et al., Multiobjective evolutionary algorithms: A survey of the state of the art, Swarm and Evolutionary Computation 1(1) (2011), 32–49.

Trivedi

, Srinivasan

, Sanyal

and Ghosh

, A Survey of Multi-objective Evolutionary Algorithms based on Decomposition, IEEE Transactions on Evolutionary Computation 21(3) (2017), 440–462.

Deb

, Pratap

, Agarwal

and Meyarivan

, A fast and elitist multiobjective genetic algorithm: NSGA-II, IEEE Transactions on Evolutionary Computation 6(2) (2002), 182–197.

Laumanns

, Thiele

, Deb

and Zitzler

, Combining convergence and diversity in evolutionary multiobjective optimization, Evolutionary Computation 10(3) (2002), 263–282.

Hadka

and Reed

, Borg: an auto-adaptive many-objective evolutionary computing framework, Evolutionary Computation 21(2) (2013), 231–259.

Wang

and Jiang

, Fuzzy-Dominance and its Application in Evolutionary Many Objective Optimization, International Conference on Computational Intelligence and Security Workshops (2007), 195–198.

Elarbi

, Bechikh

, Gupta

, Ben Said

and Ong

Y.S.

, A New Decomposition-Based NSGA-II for Many-Objective Optimization, IEEE Transactions on Systems, Man, and Cybernetics: Systems 48(7) (2018), 1191–1210.

10.

, Kwonga

, Li

, Zheng

and Shen

, Achieving balance between proximity and diversity in multi-objective evolutionary algorithm, Information Sciences 182(1) (2012), 220–242.

11.

Zitzler

and K"unzli

, Indicator-Based Selection in Multiobjective Search, Lecture Notes in Computer Science 3242 (2004), 832–842.

12.

Tian

, Cheng

, Zhang

, Cheng

and Jin

, An Indicator Based Multi-Objective Evolutionary Algorithm with Reference Point Adaptation for Better Versatility, IEEE Transactions on Evolutionary Computation 22(4) (2018), 609–622.

13.

Wei

, Li

, Fan

, Sun

and Hu

, A Hybrid Multiobjective Particle Swarm Optimization Algorithm Based on R2 Indicator, IEEE Access 6 (2018), 14710–14721.

14.

Zhang

and Hui

, MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition, IEEE Transactions on Evolutionary Computation 11(6) (2007), 712–731.

15.

Wang

, Zhang

, Zhou

, Gong

and Jiao

, Adaptive Replacement Strategies for MOEA/D, IEEE Transactions on Cybernetics 46(2) (2017), 474–486.

16.

, Yongshan

, Xiaobo

, Zhihua

and Yaoming

, A Multiobjective Optimization-Based Sparse Extreme Learning Machine Algorithm, Neurocomputing 317 (2018), 88–100.

17.

Giagkiozis

and Fleming

, Methods for multi-objective optimization: An analysis, Information Sciences 293 (2015), 338–350.

18.

While

, Hingston

, Barone

, et al., A faster algorithm for calculating hypervolume, IEEE transactions on evolutionary computation 10(1) (2006), 29–38.

19.

Das

and Dennis

J.E.

, Normal-boundary intersection: A new method for generating the pareto surface in nonlinear multicriteria optimization problems, SIAM J on Optimization 8(3) (1998), 631–657.

20.

Zhang

, Liu

and Li

, The performance of a new version of MOEA/D on CEC09 unconstrained MOP test instances, in Evolutionary Computation, 2009. CEC ’09. IEEE Congress on (2009), 203–208.

21.

Jaszkiewicz

, On the performance of multiple-objective genetic local search on the 0/1 knapsack problem – a comparative experiment, IEEE transactions on evolutionary computation 6(4) (2002), 402–412.

22.

Tan

Y.Y.

, Jiao

Y.C.

, Li

, et al., MOEA/D + uniform design: A new version of MOEA/D for optimization problems with many objectives, Computers operations research 40(6) (2013), 1648–1660.

23.

Deb

and 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 18(4) (2014), 577–601.

24.

Zitzler

and Thiele

, Multiobjective evolutionary algorithms: A comparative case study and the strength Pareto approach, IEEE Transactions on Evolutionary Computation 3(4) (1999), 257–271.

25.

Rui

, Jian

, Ishibuchi

, Wu

and Tao

, On the effect of reference point in MOEA/D for multi-objective optimization, Applied Soft Computing 58 (2017), 25–34.

26.

Bosman

P.A.N.

and Thierens

, The balance between proximity and diversity in multiobjective evolutionary algorithms, 7(2) (2003), 174–188.

27.

Cheng

, Jin

, Olhofer

and Sendhoff

, A Reference Vector Guided Evolutionary Algorithm for Many-Objective Optimization, IEEE Transactions on Evolutionary Computation 20(5) (2016), 773–791.

28.

Agrawal

R.B.

, Deb

and Agrawal

R.B.

, Simulated Binary Crossover for Continuous Search Space, Complex Systems 9(3) (1994), 115–148.

29.

Fan

, Wei

, Li

and Hu

, A novel multi-objective PSO algorithm based on completion-checking, Journal of Intelligent & Fuzzy Systems 34(1) (2018), 321–333.

30.

and Zhang

, Multiobjective optimization problems with complicated Pareto sets, MOEA/D and NSGA-II, IEEE Transactions on Evolutionary Computation 13(2) (2009), 284–302.

31.

Hansen

and Ostermeier

, Adapting arbitrary normal mutation distributions in evolution strategies: the covariance matrix adaptation, IEEE International Conference on Evolutionary Computation (1996), 312–317.

32.

Hansen

and Ostermeier

, Completely Derandomized Self-Adaptation in Evolution Strategies, Evolutionary Computation 9(2) (2001), 159–195.

33.

Igel

, Hansen

and Roth

, Covariance Matrix Adaptation for Multi-objective Optimization, Evolutionary Computation 15(1) (2007), 1–28.

34.

Castro

O.R.

, Santana

, Lozano

J.A.

, et al., Combining CMA-ES and MOEA/DD for many-objective optimization, IEEE Congress on Evolutionary Computation (CEC) (2017).

35.

, Deb

, Zhang

and Kwong

, An Evolutionary Many-Objective Optimization Algorithm Based on Dominance and Decomposition, IEEE Transactions on Evolutionary Computation 19(5) (2014), 694–716.

36.

, Zhang

and Deng

, Biased Multiobjective Optimization and Decomposition Algorithm, IEEE Transactions on Cybernetics 47(1) (2016), 52–66.

37.

Kalyanmoy

, Multi objective optimization using evolutionary algorithms, John Wiley and Sons (2001), 124–124.

38.

Hansen

and Ostermeier

, Completely derandomized self-adaptation in evolution strategies, Evolutionary Computation 9(2) (2001), 159–195.

39.

Hansen

, M"uller

S.D.

and Koumoutsakos

, Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES), Evolutionary Computation 11(1) (2014), 1–18.

40.

M"uller

S.D.

, Hansen

and Koumoutsakos

, Increasing the Serial and the Parallel Performance of the CMA-Evolution Strategy with Large Populations, International Conference on Parallel Problem Solving from Nature (2002), 422–431.

41.

Cornell

J.A.

, Experiments With Mixtures: Designs, Models, and the Analysis of Mixture Data, John Wiley and Son (2011).

42.

Auger

and Hansen

, A restart CMA evolution strategy with increasing population size, Proc IEEE Congress on Evolutionary Computation 2 (2005), 1769–1776.

43.

De Castro

L.N.

and Von Zuben

F.J.

, A Learning and optimization using the clonal selection principle, IEEE Transactions on Evolutionary Computation 6(3) (2001), 239–251.

44.

Huband

, Hingston

, Barone

and While

R.L.

, A review of multiobjective test problems and a scalable test problem toolkit, IEEE Transactions on Evolutionary Computation 10(5) (2006), 477–506.

45.

Zhou

, Jin

, Zhang

, Sendhoff

and Tsang

, Combining Model-based and Genetics-based Offspring Generation for Multi-objective Optimization Using a Convergence Criterion, IEEE Congress on Evolutionary Computation (2006), 892–899.

46.

Zhou

, Zhang

, Jin

, Tsang

and Okabe

, A model-based evolutionary algorithm for Bi-objective optimization, IEEE Congress on Evolutionary Computation (2005), 2568–2575.

47.

Zitzler

, Thiele

, Laumanns

, Fonseca

C.M.

and Fonseca

V.G.D.

, Performance Assessment of Multiobjective Optimizers: An Analysis and Review, IEEE Transactions on Evolutionary Computation 7(2) (2003), 117–132.

48.

Schutze

, Esquivel

, Lara

and Coello

C.A.C.

, Using the Averaged Hausdorff Distance as a Performance Measure in Evolutionary Multiobjective Optimization, IEEE Transactions on Evolutionary Computation 16(4) (2012), 504–522.

49.

Rudolph

, Convergence Analysis of Canonical Genetic Algorithms, IEEE Transactions on Neural Networks 5(1) (1994), 96–101.

50.

Juan

, Zixing

and Jianqin

, Premature convergence in genetic algorithm: analysis and prevention based on chaos operator, Proceedings of the 3rd World Congress on Intelligent Control and Automation (Cat. No. 00EX393) 1 (2000), 495–499.

51.

, A modified particle swarm optimization based on genetic algorithm and chaos, Proceeding of the 11th World Congress on Intelligent Control and Automation (2014), 509–512.

52.

Fan

, Wei

, Sun

, et al., An enhanced reference vectors-based multi-objective evolutionary algorithm with neighborhood-based adaptive adjustment, Neural Computing and Applications (8) (2019).