Saúl Zapotecas-Martínez

Bilel Derbel
email: bilel.derbel@univ-lille1.fr

Arnaud Liefooghe
email: arnaud.liefooghe@univ-lille1.fr

Dimo Brockhoff
email: dimo.brockhoff@inria.fr

Hernán E Aguirre

Kiyoshi Tanaka
email: ktanaka@shinshu-u.ac.jp

Injecting CMA-ES into MOEA/D

Keywords: I.2.8 [Computing Methodologies]: Artificial Intelligence-Problem Solving, Control Methods, and Search Multi-objective Optimization, Decomposition-based MOEAs, Covariance Matrix Adaption Evolution Strategy

MOEA/D is an aggregation-based evolutionary algorithm which has been proved extremely efficient and effective for solving multiobjective optimization problems. It is based on the idea of decomposing the original multi-objective problem into several singleobjective subproblems by means of well-defined scalarizing functions. Those single-objective subproblems are solved in a cooperative manner by defining a neighborhood relation between them. This makes MOEA/D particularly interesting when attempting to plug and to leverage single-objective optimizers in a multi-objective setting. In this context, we investigate the benefits that MOEA/D can achieve when coupled with CMA-ES, which is believed to be a powerful single-objective optimizer. We rely on the ability of CMA-ES to deal with injected solutions in order to update different covariance matrices with respect to each subproblem defined in MOEA/D. We show that by cooperatively evolving neighboring CMA-ES components, we are able to obtain competitive results for different multi-objective benchmark functions.

INTRODUCTION

A multi-objective optimization problem (MOP) refers to the situation where several conflicting objectives are to be optimized simultaneously. In such a setting, solving a MOP consists in finding a whole set of solutions providing both good and diverse compromises with respect to the corresponding objectives. Computing such a set is actually a challenging task for which one can find different methodological frameworks and different algorithmic concepts. Evolutionary algorithms (EA) are particularly well-suited in order to compute an accurate approximation of the so-called Pareto set, that is the set of the best achievable objective trade-offs. Generally speaking, multi-objective EAs can be classified in different classes, ranging from dominance-based algorithms [START_REF] Deb | A Fast and Elitist Multiobjective Genetic Algorithm: NSGA-II[END_REF][START_REF] Kukkonen | GDE3: the third evolution step of generalized differential evolution[END_REF], indicator-based algorithms [START_REF] Zitzler | Indicator-based selection in multiobjective search[END_REF][START_REF] Zapotecas Martínez | Using a Family of Curves to Approximate the Pareto Front of a Multi-Objective Optimization Problem[END_REF], and aggregation-based algorithms [START_REF] Ishibuchi | Multi-Objective Genetic Local Search Algorithm[END_REF][START_REF] Zhang | MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition[END_REF]. In this paper, we are interested in the MOEA/D (multi-objective optimization based on decomposition) framework for which we can witness a growing interest from the community due to its simplicity and to its effectiveness when applied to a broad range of multi-objective optimization problems.

One can find several variants and implementations of MOEA/D, which are all based on seemingly the same concept. Instead of attempting to directly solve the MOP as a whole, all MOEA/D variants decompose the original MOP into several single-objective subproblems, and solve them cooperatively. More specifically, a set of weighting coefficient vectors is used to define different scalarized subproblems for which one solution is maintained and evolved over time. In this respect, any single-objective optimizer can be potentially plugged into MOEA/D as search engine for the singleobjective subproblems. This has been done for example with polynomial mutation and SBX crossover [START_REF] Zhang | MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition[END_REF] as well as differential evolution (DE) [START_REF] Li | Multiobjective Optimization Problems With Complicated Pareto Sets, MOEA/D and NSGA-II[END_REF], where the latter is believed to be among the best performing implementations [START_REF] Li | Multiobjective Optimization Problems With Complicated Pareto Sets, MOEA/D and NSGA-II[END_REF]. However, 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 scalarized problems.

Instead of solving every subproblem independently, MOEA/D defines a neighborhood relation among them, and allows solutions to be exchanged between neighboring problems in two ways. On the one hand, the currently best-known solution of a neighboring problem can be employed to change ones own search distribution, e.g., by using them as parents in a crossover operator. On the other hand, newly generated solutions are also allowed to be actively transferred to a neighboring subproblem in order to replace its current best solution. On an abstract level, we can say that in MOEA/D, the single-objective optimizers applied at each subproblem are not only acting selfishly to improve the solution of their own problem by "stealing information" from others, but they also behave in an altruistic way by helping to improve the solutions of their neighbors.

In this paper, we investigate new opportunities offered by the flexibility of MOEA/D in incorporating novel single-objective optimizers. More specifically, we explore the strengths of the wellestablished CMA-ES (Covariance Matrix Adaption Evolution Strategy) algorithm considered as a state-of-the-art optimizer for singleobjective blackbox continuous problems [START_REF] Hansen | Completely Derandomized Self-Adaptation in Evolution Strategies[END_REF]. Our interest in combining the CMA-ES with MOEA/D stems from two sources. On the one hand, CMA-ES has been shown experimentally to be among the best-performing single-objective blackbox algorithms on the well-established BBOB testbed, with typically superior performance to DE variants and other numerical optimizers-especially when the problems are difficult and the budgets are not too small [START_REF] Hansen | Comparing results of 31 algorithms from the black-box optimization benchmarking BBOB[END_REF]. On the other hand, a recent paper by Hansen [START_REF] Hansen | Injecting external solutions into CMA-ES[END_REF] shows that external solutions can be easily injected into the algorithm to gain from good candidate solutions that are not sampled directly from the algorithm itself. Both aspects together make the CMA-ES with injection a highly interesting candidate to be used within MOEA/D and its cooperative optimization in which solutions are exchanged between single-objective subproblems.

In this paper, we therefore propose a novel variant of MOEA/D where CMA-ES is used as the core single-objective evolution engine and where the injection idea allows to incorporate information from neighboring scalarizing problems into the search distributions. We show that by appropriately injecting solutions coming from the CMA-ES sampling process into neighboring subproblems, we can derive novel effective variation mechanisms which are compatible with the decomposition-based concepts used with the MOEA/D framework. To assess the validity of our approach, we conduct an experimental study where, in particular, we analyze the performance of the proposed MOEA/D-CMA algorithm compared to MOEA/D-DE using the CEC 2009 box-constrained benchmark functions. Besides being able to obtain competitive results, the investigations conducted in this paper highlight promising alternatives for designing novel efficient multi-objective optimization algorithms falling in the class of aggregation-based EAs.

Note that already several multi-objective versions of the CMA-ES algorithm exist [START_REF] Igel | Covariance matrix adaptation for multi-objective optimization[END_REF][START_REF] Igel | Steady-state selection and efficient covariance matrix update in the multi-objective cma-es[END_REF] which, however, do not resemble the framework of MOEA/D but instead aim at maximizing the hypervolume of a solution set in the framework of indicator-based algorithms. These algorithms are not in the focus of this paper although a numerical comparison would be interesting for the future.

The rest of this paper is organized as follows. In Section 2, we review some basic concepts related to multi-objective optimization as well as to the MOEA/D framework. In Section 3, we recall the main algorithmic components of the single objective CMA-ES as well as the idea of injection that will serve as a basis towards its successful incorporation into MOEA/D. In Section 4, we describe our proposed MOEA/D-CMA approach and discuss its design components in detail. In Section 5, we describe our experimental study and discuss our main findings. Section 6 finally concludes the paper and discusses related open research directions.

BASIC CONCEPTS

Multi-objective Optimization

Assuming minimization, a continuous multi-objective optimization problem (MOP), can be stated as:

minimize x∈Ω F(x) = ( f 1 (x), . . . , f k (x)) (1)
where Ω ⊂ R n defines the decision space and F is defined as the vector of the objective functions where each f j : Ω → R ( j = 1, . . . , k) is the function to be minimized. In this paper we consider the box-constrained case, i.e., Ω = n i=1 [a j , b j ]. Therefore, each variable vector x = (x 1 , . . . , x n ) ∈ Ω is such that a i ≤ x i ≤ b i for all i ∈ {1, . . . , n}. In order to describe the concept of optimality in which we are interested, the following definitions are introduced [START_REF] Miettinen | Nonlinear Multiobjective Optimization[END_REF]. Definition 1. Let x, y ∈ Ω, we say that x dominates y (denoted by x ≺ y) with respect to the problem defined in equation [START_REF] Auger | A restart CMA evolution strategy with increasing population size[END_REF] if and only if: 1) f i (x) ≤ f i (y) for all i ∈ {1, . . . , k} and 2) f j (x) < f j (y) for at least one j ∈ {1, . . . , k}. Definition 2. Let x ∈ Ω, we say that x is a Pareto optimal solution, if there is no other solution y ∈ Ω such that y ≺ x . Definition 3. The Pareto optimal set PS is defined by: PS = {x ∈ Ω|x is a Pareto optimal solution} and its image PF = {F(x)|x ∈ PS }) is called Pareto front PF.

MOEA/D

The Multi-Objective Evolutionary Algorithm Based on Decomposition (MOEA/D) [START_REF] Zhang | MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition[END_REF], transforms a MOP into several scalarizing subproblems. Therefore, an approximation of the Pareto set is obtained by solving the N scalarizing subproblems in which a MOP is decomposed.

Considering W = {w 1 , . . . , w N } as the well-distributed set of weighting coefficient vectors, MOEA/D seeks the best solution to each subproblem defined by each weight vector using the Penalty Boundary Intersection (PBI) approach [START_REF] Zhang | MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition[END_REF], which is in the form:

minimize: g pbi (x|w i , z) = d 1 + θd 2 (2)
such that

d 1 = |(F(x) -z) w| ||w|| and d 2 = (F(x) -z) -d 1 w ||w||
where x ∈ Ω ⊂ R n and z j = min{ f j (x)|x ∈ Ω} for i = 1, . . . , N and j = 1, . . . , k. Since z = (z 1 , . . . , z k ) is unknown, MOEA/D states each component z j by the minimum value for each objective f j found during the search, for j = 1, . . . , k.

In MOEA/D, a neighborhood of a weight vector w i is defined as a set of its closest weight vectors in W. Therefore, the neighborhood of the weight vector w i contains all the indexes of the T closest weight vectors to w i .

Throughout the evolutionary process, MOEA/D finds the best solution to each subproblem maintaining a population of N solutions P = {x 1 , . . . , x N } where x i ∈ Ω is the current solution to the i th subproblem. Algorithm 1 presents the general framework of MOEA/D, however, for a more detailed description of this algorithm the interested reader is referred to [START_REF] Zhang | MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition[END_REF].

CMA-ES

The Covariance Matrix Adaptation Evolution Strategy (CMA-ES, [START_REF] Hansen | Completely Derandomized Self-Adaptation in Evolution Strategies[END_REF]) is one of the state-of-the-art numerical blackbox optimization algorithms available-outperforming other algorithms especially for difficult functions and larger budgets [START_REF] Hansen | Comparing results of 31 algorithms from the black-box optimization benchmarking BBOB[END_REF]. To minimize a function f : R n → R, CMA-ES iteratively samples λ solutions from a multivariate normal distribution, parameterized by a mean vector m τ ∈ R n , a step size σ τ > 0, and a covariance matrix C τ ∈ R n×n at time step τ. After the evaluation and ranking of the λ candidate solutions, their relative steps (between the sampled points and the old mean, sorted according to the solution's objective function value) are used to update the parameters of the sampling distribution. This loop of sampling and sample distribution update is repeated until any of several stopping criteria is met.

More concretely, the CMA-ES version, employed in this paper and following [START_REF] Hansen | The CMA evolution strategy: a comparing review[END_REF], works as follows. In the initialization step at τ = 0, the mean m 0 is typically sampled uniformly at random in a bounded search domain [a, b] ∈ R n , the initial covariance matrix C 0 is chosen as the identity matrix, and the evolution paths p ρ and p c are set to 0 ∈ R n .

At each iteration, CMA-ES then samples λ ≥ 1 candidate solutions

x i ∈ R n (1 ≤ i ≤ λ) according to a multivariate normal distribution (x i ∼ N(m τ , σ 2 τ C τ ) ∼ m τ + σ τ • N(0, C τ )
). The candidate solutions x i are evaluated with respect to the objective function f : R n → R and ranked according to f -values where we use the notation x i:λ for the ith best solution (i.e. f (

x 1:λ ) ≤ • • • ≤ f (x μ:λ ) ≤ f (x μ+1:λ ) ≤ • • • ≤ f (x λ:λ ).
The new mean vector is computed via weighted recombination of the μ best candidate solutions as:

m τ+1 = μ i=1 ω i x i:λ = m τ + μ i=1 ω i (x i:λ -m τ ) (3)
where the positive (recombination) weights ω i > 0 with μ i=1 ω i = 1 are typically chosen linearly on the log-scale with μ ≤ λ/2.

The step size σ τ is updated using cumulative step size adaptation (CSA). The evolution path (or search path) p σ is thereby firstly updated as:

p σ ← (1 -c σ )p σ + 1 -(1 -c σ ) 2 √ μ eff C -1/2 τ m τ+1 -m τ σ τ (4)
with the help of which the step size is finally updated as:

σ τ+1 = σ τ × exp c σ d σ × ||p σ || E||N(0, I)|| -1 (5)
where c -1 σ ≈ n/3 is the backward time horizon for the evolution path p σ and larger than one,

μ eff = μ i=1 ω 2 i -1
is the variance effective selection mass and 1 ≤ μ eff ≤ μ by definition of ω i ,

C -1/2 τ = C -1
τ is the unique symmetric square root of the inverse of C τ , and d σ is a damping parameter usually close to one. For d σ = ∞ or c σ = 0 the step size remains unchanged. Note that the step size σ τ is increased if and only if p σ is larger than the expected step length of a fully random sample:

E N(0, I) = √ 2 Γ((n + 1)/2)/Γ(n/2) ≈ √ n (1 -1/(4 n) + 1/(21 n 2 )) (6)
and decreased if it is smaller. For this reason, the step size update tends to make consecutive steps C -1 τ -conjugate, in that after the adaptation has been successful,

m τ+2 -m τ+1 σ τ+1 λ C -1 τ m τ+1 -mτ στ
≈ 0 holds. Finally, the covariance matrix is updated by means of rank-one and rank-μ updates for which, again, the respective evolution path is firstly updated as:

p c ← (1-c c )p c +1 [0,α √ n] ( p σ ) 1 -(1 -c c ) 2 √ μ eff m τ+1 -m τ σ τ (7)
which is used to update the covariance matrix as:

C τ+1 = (1 -c 1 -c μ + c s )C τ + c 1 p c p c rank-one update + c μ μ i=1 ω j x i:λ -m τ σ τ x i:λ -m τ σ τ rank-μ update (8)
where denotes the transpose and c -1 c ≈ n/4 is the backward time horizon for the evolution path p c and larger than one, α ≈ 3/2 and the indicator function

1 [0,α √ n] ( p σ ) evaluates to one iff p σ ∈ [0, α √ n] or, in other words, p σ ≤ α √ n, which is usually the case. The constant c s = (1 -1 [0,α √ n] ( p σ ) 2 ) c 1 c c (2 -c c
) makes partly up for the small variance loss in case the indicator is zero, c 1 ≈ 2/n 2 is the learning rate for the rank-one update of the covariance matrix and c μ ≈ μ eff /n 2 is the learning rate for the rank-μ update of the covariance matrix and must not exceed 1c 1 .

This completes an iteration of CMA-ES which continues with the sampling of new solutions and updates of the sample distribution parameters until a stopping criterion is met (see Sec. 4 for details). For a more detailed description, the interested readers are referred to [START_REF] Hansen | The CMA evolution strategy: a comparing review[END_REF].

CMA-ES with Injection. If external solutions are available to CMA-ES, for example from a gradient step or from evaluating a meta-model of the (expensive) objective function, these solutions can be directly used in the update of the sampling distribution's parameters. The only change to make the above algorithm work when handling such "injections" of solutions-as argued in [START_REF] Hansen | Injecting external solutions into CMA-ES[END_REF]-is to restrict the distance between injected solution and previous mean and to rescale the corresponding search step accordingly before taking it into account in the updates of mean, step size, and covariance matrix.

MOEA/D-CMA

The proposed MOEA/D with Covariance Matrix Adaption Evolution Strategy (MOEA/D-CMA) is given in the template of Algorithm 2. As we can see, the proposed approach consists in an initialization step, recombination step and a specific CMA evolution process. In the following, we will provide a step-by-step descrip-

Procedure Initialize( p i c , p i σ , C i , m i , σ i , τ i , x i , σ init ) p i c ← p i σ ← 0 ; // Initial cumulative paths C i ← I ; // Initial covariance matrix m i ← x i , σ i ← σ init ;
// Initial mean and sigma τ i ← 0 ; // Initial iteration tion of the different components involved in our proposed approach.

Preliminary Considerations and Initialization. MOEA/D-CMA decomposes the problem (1) into N single-objective optimization subproblems. It is in fact based on MOEA/D in the sense that it updates a neighborhood of solutions which solves a set of neighboring subproblems. However, instead of using genetic operators (crossover and mutation), each subproblem evolves the parameters of a multivariate normal distribution and samples according to this distribution, as presented in Section 3. Therefore, analogously to MOEA/D, MOEA/D-CMA employs a well-distributed set of weight vectors W = {w1 , . . . , w N } to define the set of singleobjective subproblems to be optimized cooperatively. Each weight vector w i defines a scalarized function by using the PBI approach (equation ( 2)). However, the use of other scalarized functions is also possible, see for example those presented in [START_REF] Miettinen | Nonlinear Multiobjective Optimization[END_REF]. Each of the subproblems is then solved by a CMA-ES instance, involving a mean vector, step size, covariance matrix, and the corresponding cumulation paths as presented in Section 3. The i th individual is thereby denoted by the six-tuple p i c , p i σ , C i , m i , σ i , τ i with τ i being an iteration counter.

Let us consider P = {x 1 , . . . , x N } as the set of initial random solutions generated in the feasible search space Ω. At the beginning, the reference point z is set with infinite values and then, its jth component is updated with the minimum value of the jth objective function f j found along the optimization process. Analogously to MOEA/D, a neighborhood B i which contains the indexes of the T closest weight vectors to w i is defined. Each solution in P defines the initial mean of each CMA-ES's search distribution. The covariance matrices, step sizes and evolution paths are initialized the same for each subproblem as shown in procedure Initialize.

Recombination. In Step 1 of Algorithm 2, the ith subproblem generates a set of solutions V i by using the multivariate normal distribution which is in the form N(m i , σ 2 i C i ) 1 . It is worth noticing that the new solutions in V i could be located outside of the feasible region. In this case, we substitute each infeasible solution by its closest vector in the feasible region, which is denoted by V i rep . If the solution is feasible we set v rep = v, this procedure is referred to as Repair in Algorithm 2.

After sampling new candidate solutions, we update the neighboring solutions maintained at each subproblem. For this purpose we adopt the update mechanism proposed in MOEA/D-DE [START_REF] Li | Multiobjective Optimization Problems With Complicated Pareto Sets, MOEA/D and NSGA-II[END_REF], where a maximum number of replacements n r (lines 20-21) and a dynamic selection of neighborhood (line 18) are employed. In this way, the probability of losing diversity is reduced and the preservation of the best solutions to each neighboring subproblem is maintained during the sampling procedure.

It is important to note that in order to avoid stagnation and allow for restarts of the single optimization runs in each subproblem, we adopt a reset procedure denoted by ResetCriteria in line 12. This procedure first checks the following four criteria taken from the literature (e.g., see [START_REF] Auger | A restart CMA evolution strategy with increasing population size[END_REF]): Therefore, if the components of the ith tuple satisfies one of the above criteria, the tuple p i c , p i σ , C i , m i , σ i , τ i is reset using the same Initialize procedure as introduced previously. It is worth noticing that with each reinitialization, we reset the step size σ i to the initial value divided by a factor of two, and using the current best solution for this subproblem as the new initial mean before continuing. More precisely, the ith mean is set as m i = x i and

σ i = σ init /2.
Evolution Strategy. After performing the recombination and replacement stage, we consider to evolve and update the components of every single CMA-ES optimizer at each single-objective subproblem. This is the main aim of Step 2. Generally speaking, we adopt a variant of the original CMA-ES presented in [START_REF] Hansen | Injecting external solutions into CMA-ES[END_REF], where the injection of external solutions into the evolution of the covariance matrix is possible. As it was pointed out by Hansen [START_REF] Hansen | Injecting external solutions into CMA-ES[END_REF], the injection of solutions in the evolutionary process of CMA-ES could improve the adaptation of the matrix if the injected solutions provide sufficient information about the fitness landscape. Our main concern is then to design an accurate injection mechanism which can deal in a proper way with our multi-objective setting.

We argue that, under some assumptions, from the Karush-Kuhn-Tucker conditions, it can be deduced that the PS of a continuous MOP with k objectives forms a (k -1)-dimensional piecewise continuous manifold in decision variable space [START_REF] Miettinen | Nonlinear Multiobjective Optimization[END_REF]. It means that an optimal solution of the scalarizing problem defined by a weight vector w p is close to the one defined by another weight vector w q (p q), if w p and w q are close to each other.

We hypothesize that during the search, the sample distribution of each subproblem converges to the region in which the optimal solution to each subproblem is located. Therefore, considering continuous MOPs and having the reference to the neighboring candidate solutions to each subproblem, the solutions to be injected are chosen precisely from this neighborhood. With that, a set of promising solutions are injected while eventually optimizing each separate subproblem.

Injecting Solutions. In Step 2 of Algorithm 2, the set Q denotes the set of candidate solutions to be injected in the adaptation of the covariance matrix. This set of solutions, in fact, contains the samples, generated by the ith subproblem, and the best solution for the ith scalarized subproblem (defined by w i ) found among the samples V j (for j ∈ B i ). Since we are injecting (hypothetically) promising solutions, we consider the repaired solutions V j rep . Note q j:|Q|m i σ i where:

ResetCriteria( p i c , p i σ , C i , m i , σ i , τ i , x i , σ init /2); 13 V i = V i rep = ∅; 14 for j = 1, . . . , λ do 15 V i ← V i ∪ {v j }, where v j ←∼ N(m i , σ 2 i C i )}; 16 V i rep ← V i rep ∪ {v j rep }, where v j rep ← Repair(v j );
f fit (q 1:|Q| |w i , z) ≤ • • • ≤ f fit (q μ:|Q| |w i , z) ≤ • • • and q j ∈ Q; 31 y j ← α clip c y , ||(C i ) -1/2 y j || × y j , if q j:|Q| was injected; 32 Δm ← μ j=1 ω j y j ; 33 m i ← m i + c m σ i Δm; 34 m i rep ← Repair(m i ); 35 UpdateStepSize(C i , Δm, p i σ , σ i ); 36 UpdateCovarianceMatrix(C i , p i c , Δm, y 1 , . . . , y μ ); 37 τ i ← τ i + 1; 38 for l ∈ π do 39 if g pbi (m i rep |w l , z) < g pbi (x l |w l , z) then x l ← v j rep ;
40 return P ← {x 1 , . . . , x N };

that the samples V i are the samples given by the ith subproblem and i j.

To obtain a ranking among the to-be-injected solutions, the sorting is carried out by using the definition of an auxiliary fitness function adopted from [START_REF] Hansen | A method for handling uncertainty in evolutionary optimization with an application to feedback control of combustion[END_REF], which allows us, to deal with the boxconstrained case, and it is stated as follows.

f fit (q j |w i , z) = g pbi (q j rep |w i , z) + α||q j -q j rep || 2 (9)
It is worth noticing that line 31 is introduced as part of the updated process of the CMA-ES with injection [START_REF] Hansen | Injecting external solutions into CMA-ES[END_REF], where α clip (c, x) = 1 ∧ c

x and the notation a ∧ bc + d refers to the minimum of a and bc + d. In this way, the tuple p i c , p i σ , C i , m i , σ i , τ i is updated by means of equations provided by [START_REF] Hansen | Injecting external solutions into CMA-ES[END_REF], which are referred to as the UpdateStepSize. More precisely, the step size σ i and the corresponding evolution path p i c are updated according to the following equations.

p i σ ← (1 -c σ )p i σ + c σ (2 -c σ )μ eff (C i ) -1/2 Δm ( 10
)
σ i ← σ i × exp Δ max σ ∧ c σ d σ ||p i σ || E||N(0, I)|| -1 (11)
The covariance matrix is updated by the procedure UpdateCova-rianceMatrix by using Equation [START_REF] Hansen | A method for handling uncertainty in evolutionary optimization with an application to feedback control of combustion[END_REF]. However, because we are injecting a set of solutions, the learning path p i c is updated by means of the following equation.

p i c ← (1 -c c )p i c + h σ c c (2 -c c )μ eff Δm (12)
In fact, as it is shown in [START_REF] Hansen | Injecting external solutions into CMA-ES[END_REF], Equations ( 10), ( 11) and ( 12) differ from Equations ( 4), ( 5) and [START_REF] Hansen | Comparing results of 31 algorithms from the black-box optimization benchmarking BBOB[END_REF], by introducing the parameter Δm and Δ max σ . Note however, that using Δ max σ = +∞ and Δm =

mτ i -m τ i -1 σ i
, the original equations of CMA-ES are recovered. It is also possible to inject an arbitrary mean, which shifts the current mean m i by means of additional equations. In the study presented herein, we focused only on the injection of already evaluated solutions. However, the injection of a determined mean, is indeed, a possible path for future research.

In the remainder, our proposed approach is denoted as MOEA/D-CMA+I as exactly described in Algorithm 2; however, we shall also consider a second variant of this algorithm, denoted by MOEA/D-CMA, by taking off the injection mechanism, i.e., Q = V i and Δ max σ = +∞.

EXPERIMENTAL STUDY

In order to analyze the efficiency of the proposed approach with and without injection, i.e. MOEA/D-CMA+I and MOEA/D-CMA, we compare its performance against two competing algorithms: (i) the conventional Multi-Objective Evolutionary Algorithm based on Decomposition (MOEA/D) [START_REF] Zhang | MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition[END_REF], and (ii) an improved variant of MOEA/D which uses differential evolution operators and a dynamic neighborhood selection. The corresponding algorithm is called MOEA/D-DE [START_REF] Li | Multiobjective Optimization Problems With Complicated Pareto Sets, MOEA/D and NSGA-II[END_REF]. In this section, the benchmark problems and the performance assessment design adopted in our analysis are presented.

Experimental Setup

We consider the continuous MOPs with complicated Pareto sets proposed in [START_REF] Zhang | Multiobjective Optimization Test Instances for the CEC 2009 Special Session and Competition[END_REF], and extracted from the CEC 2009 special session and competition on the performance assessment of multi-objective optimization algorithms. This benchmark test suite has been specifically designed to resemble complicated real-life optimization problems. The MOPs therein present different properties in terms of separability, multi-modality, and shape of the Pareto front, i.e. convexity, concavity, discontinuities, gaps, etc. More particularly, we consider all box-constrained functions UF1-F10 under their original setting [START_REF] Zhang | Multiobjective Optimization Test Instances for the CEC 2009 Special Session and Competition[END_REF], with UF1-UF7 being two-objective problems and UF8-UF10 being three-objective problems. Notice that, for all of them, the number of variables is n = 30, and the Pareto fronts lie in the hyper-box [0, 1] k , where k denotes the number of objective functions for the problem under consideration.

All the competing algorithms considered in this study were compared by following the performance assessment experimental setup recommended in [START_REF] Knowles | A Tutorial on the Performance Assessment of Stochastic Multiobjective Optimizers[END_REF].

Relative Hypervolume Deviation. The first performance measure indicates the relative hypervolume achieved by the Pareto set approximation given by an algorithm. This Relative Hypervolume (RHV) deviation can be computed by:

RHV(A) = HV(R) -HV(A) HV(R) ( 13
)
where HV is the hypervolume indicator [START_REF] Zitzler | Multiobjective Optimization Using Evolutionary Algorithms -A Comparative Case Study[END_REF], A is an approximation set and R is a reference set for the instance under consideration. The reference vector is set to (2, . . . , 2) . In this quality indicator, a lower value means a better approximation set.

Inverted Generational Distance. The Inverted Generational Distance (IGD, [START_REF] Coello Coello | Solving Multiobjective Optimization Problems using an Artificial Immune System[END_REF]) indicates how far a given Pareto front approximation is from a reference set. Let R be the true Pareto front, the IGD for a set of approximated solutions A is calculated as:

IGD(A) = 1 |R| |R| i=1 d i ( 14
)
where 2 and k is the number of objective functions. A value of zero in this performance measure, indicates that all the solutions obtained by the algorithm are on the true Pareto front.

d i = min |A| j=1 k l=1 ( f l (i) -f l ( j))
Both performance measures (i.e. RHV and IGD) are computed by using the reference sets available at: http://dces.essex.ac.uk/ staff/qzhang/moeacompetition09.htm.

Parameter Setting

As we mentioned before, we consider two variants of the proposed MOEA/D-CMA paradigm. The first version exactly maps to the one presented in Algorithm 2 (i.e. MOEA/D-CMA+I). The second version, denoted by MOEA/D-CMA, corresponds to the same framework, but without performing the phase where solutions are injected to the neighboring subproblems. That is, Step 2.1 is not performed in Algorithm 2, i.e. Q = V i . This shall allow us to appreciate by how much the search process can actually benefit from injecting solutions sampled from neighboring subproblems.

Besides these two versions, we also consider the conventional MOEA/D [START_REF] Zhang | MOEA/D: A Multiobjective Evolutionary Algorithm Based on Decomposition[END_REF] as well as the more sophisticated MOEA/D-DE [START_REF] Li | Multiobjective Optimization Problems With Complicated Pareto Sets, MOEA/D and NSGA-II[END_REF] in our comparative study. For all competing algorithms, the set of weighting coefficient vectors is generated following a simplexlattice design [START_REF] Scheffé | Experiments With Mixtures[END_REF]. The settings of N and W = {w 1 , . . . , w N } are controlled by a parameter H. More precisely, let {w 1 , . . . , w N } be the set of weight vectors. Each individual weight w i j , such that i = 1, . . . , N and j = 1, . . . , k, can take a value from 0 H , 1 H , . . . , H H . Therefore, the number of vectors in W is given by N

= C k-1 H+k-1 ,
where k is the number of objective functions. In this work, we set H = 99 for two-objective problems, and H = 19 for three-objective problems, that is 100 and 210 weight vectors for MOPs having two and three objectives, respectively.

To fix the parameters of MOEA/D-CMA+I and MOEA/D-CMA, we essentially get inspired by the standard setting values suggested in [START_REF] Hansen | The CMA evolution strategy: a comparing review[END_REF], which are summarized in Table 1. However, following our empirical observations, we reduce the population size λ to be the same as μ. Nonetheless, the adjustment of the λ parameter, is indeed, one important open issue that would deserve to be investigated. The remaining parameters required by the MOEA/D framework are set as follows: T = 20, η c = η m = 20, P c = 1 and P m = 1/n; which represent respectively the neighborhood size, crossover index (for Simulated Binary Crossover (SBX)), mutation index (for Polynomial-Based Mutation (PBM)), crossover rate and

= 4 + 3 ln n μ = λ 2 ω i = ln λ+1 2 -ln i μ j=1 ln λ+1 2 -ln j for i = 1, . . . , μ μ eff = μ i=1 (ω 2 i ) -1 c y = √ n + 2n/(n + 2) c m = 1 c σ = μ eff +2 n+μ eff +3 c c = 4 n+4 d σ = 1 + c σ + 2 max 0, μ eff -1 n+1 -1 c 1 = αcov min(1,λ/6) (n+1.3) 2 +μ eff c μ = 1 -c 1 ∧ α cov μ eff -2+1/μ eff (n+2) 2 +αcovμ eff /2 α cov = 2 Δ max σ = 1
mutation rate. Finally, the parameter θ in the PBI approach was set to θ = 5. For MOEA/D-DE, we adopted the set of parameters given in [START_REF] Li | Multiobjective Optimization Problems With Complicated Pareto Sets, MOEA/D and NSGA-II[END_REF]. More precisely, the differential factor was set as F = 0.5, the crossover ratio was set as CR = 1, the maximum number of replacements was set as n r = 2, and δ = 0.9. We define the initial step size as

σ init = 1 4 × (U b -L b )
, where U b and L b are the upper-and the lower-bounds in the search space. Since we consider the problems with the same boundaries in all dimensions, the decision variables of the original CEC 2009 test functions are simply rescaled without modifying the shape of the PS or the PF. Finally, the auxiliary fitness function was computed with α = 1 × 10 -5 .

For each MOP, we performed 30 independent runs and we measure the performance of the algorithms after N ×1 000 and N ×2 000 fitness function evaluations.

Numerical Results

Our results are summarized in Table 2, where we can see the relative performance of the four competing algorithms using two different stopping conditions and for the two quality indicators described above. Notice that we can read two informations in Table 2. We put in bold the best average indicator-value obtained among all algorithms. We also perform a Mann-Whitney non-parametric statistical test between each pair of algorithms in order to determine whether a given algorithm is significantly outperformed by 2: Comparison of the competing algorithms respect to the relative hypervolume deviation (RHV) and to the inverted generational distance (IGD). The first number stands for the average indicator-value (lower is better). The number in brackets stands for the standard deviation. Bold values correspond to the best average indicator-value for the instance and the indicator under consideration. Underline values correspond to algorithms that are not statistically outperformed by any other algorithm for the instance and the indicator under consideration with respect to a Mann-Whitney non-parametric statistical test with a p-value of 0.05 by using a Bonferroni correction [START_REF] Bonferroni | Teoria statistica delle classi e calcolo delle probabilita[END_REF]. any other. If no other algorithm is statistically better than the algorithm under consideration, then the corresponding indicator-value is underlined in the table. Several interesting observations can be extracted from Table 2. First, we can see that there are no significant differences between the behavior of the algorithms when considering different stopping conditions. This hopefully informs that all algorithms are rather ranking in the same manner independently of the available function-evaluation budget. More interestingly, the results of Table 2 allows us to validate the accuracy of the introduced approach from several perspectives.

0 0 1 -1 -0 1 -1 0 1 x 3 x 1 x 2 x 3 (a) UF2 0 0 1 0 0 1 0 0 1 x 3 x 1 x 2 x 3 (b) UF3 0 0 1 0 0 1 -2 0 2 x 3 x 1 x 2 x 3 (c) UF8 0 0 1 0 0 1 -2 0 2 x 3 x 1 x 2 x 3 (d) UF9

RHV IGD

M O E A / D -C M A + I M O E A / D -C M A M O E A / D M O E A / D -D E M O E A / D -C M A + I M O E A / D -C M A M O E A / D M O E A / D -
When comparing MOEA/D-CMA+I with MOEA/D-CMA, we can see that MOEA/D-CMA is outperformed by MOEA/D-CMA+I for all instances except for UF6 when using the IGD indicator. This indicates that the way in which solutions are injected from neighboring subproblem into the single-objective CMA-ES engine, and the way in which CMA-ES is taking care of those solutions when adapting its components, is drastically improving the search when compared to the straightforward version where CMA-ES is plugged within MOEA/D to simply sample new points independently at every subproblem. Notice however, that the MOEA/D-CMA version exhibits, in the worst case, very comparable performance with respect to conventional MOEA/D, and the more advanced MOEA/D-CMA+I algorithm is outperforming the conventional version of MOEA/D in almost all instances expect for instance UF5 and UF10. This indicates that the variation operator in-duced by the MOEA/D-CMA algorithm is performing well. This is confirmed when looking at the relative performance of MOEA/D-CMA+I with respect to MOEA/D-DE which is known to perform extremely well on these benchmark functions. We can see that there is no algorithm that is performing better in all the considered instances and that the difference is most of the time not statistically significant. Notice also that the difference between the two algorithms is more pronounced for the last 3 instances UF8, UF9 and UF10, in favor of MOEA/D-DE-these instances are actually the three-objective problems considered in our experiments. For the bi-objective instances, MOEA/D-CMA+I appears to be relatively competitive compared to MOEA/D-DE.

In Fig. 1, we show the Pareto set approximations obtained by our proposed MOEA/D-CMA+I for the UF2, UF3, UF8 and UF9 problems. As we can see, the benchmark functions exhibit complicated shapes and depending on the considered instance, the algorithm has some potential in accurately approaching the true Pareto set. This is actually a common behavior for the other competing algorithms as well which is to be attributed to the particular difficulty of some of the considered instances.

We recall that the main purpose of this paper is to investigate at what extent the CMA-ES single-objective optimizer could be appropriately integrated within the MOEA/D framework and what are the benefits one can obtain. As such, we can conclude from the previous set of experiments that injecting accurately solutions from neighbors plays a crucial role to this We can also conclude the newly proposed MOEA/D-CMA+I algorithm is a promising candidate to solve hard and complicated optimization problems. In fact, the CMA-ES optimizer enjoys several attractive properties, such as invariance and robustness, which are not verified by several other evolutionary operators. Hence, the experimental study presented in this paper, and the relatively good performance that MOEA/D-CMA+I is able to obtain can be viewed as a first promising step towards the establishment of highly accurate algorithms.

CONCLUSIONS AND FUTURE WORK

In this paper, we consider the single-objective CMA-ES algorithm as a plausible optimizer to be incorporated into the MOEA/D framework. We proposed to assign to each subproblem, obtained by decomposition, a single CMA-ES engine and to profit from the solutions of the neighboring subproblems in order to adapt the CMA-ES components accordingly. We thereby rely on the ability of CMA-ES to deal with external solutions in order to derive novel variation operators for the MOEA/D framework. The experimental study conducted in this paper allows us to validate the proposed approach and to show its effectiveness. Besides, the presented algorithm opens the road to further challenging research questions. In fact, the proposed approach is to be viewed as a first step enlightening how the CMA-ES algorithm can specifically be used as a single-objective optimizer when solving multi-objective problems; and this is thanks to the basic concepts introduced by MOEA/D and the ability of CMA-ES of handling external solutions. It is worth noticing that different strategies to take into account the sampled solutions from neighbors could be investigated in the future. Furthermore, it would be insightful to study more deeply the behavior of the so-obtained strategies when considering other benchmark functions. This would enable to fully appreciate the strength of the CMA-ES algorithm on a broad range of problems with different properties. It is our hope that the contribution of this paper can serve as a starting point to derive more powerful aggregation-like EMO methods based on the well established CMA-ES algorithm.

17 zl 19 c ← 0; 20 foreach l ∈ π do 21 if 22 Step 2 . 2 . 2 .

 1719202122222 ← min(z l , f l (v j rep )) ; // for l = 1, . . . , k 18 if rand() < δ then π ← B i else π ← {1, . . . , N}; g pbi (v j rep |w l , z) < g pbi (x l |w l , z) and c < n r then x l ← vj rep and c ← c + 1; Covariance Matrices Adaptation 23 for i = 1, . . . , . . . , N do 24 if rand() < δ then π ← B i else π ← {1, . . . , N}; 25 Q ← V i ; 26 2.1. Injecting Solutions 27 foreach j ∈ π do 28 Q ← Q ∪ {v }, such that: v = arg min v∈V j rep g(v|w i , z); 29 Covariance matrix adaptation 30 y j ←

Figure 1 :

 1 Figure 1: Pareto set approximations given by MOEA/D-CMA+I in UF2, UF3, UF8 and UF9

Algorithm 1 :

 1 General Framework of MOEA/D Input: N: the number of subproblems to be decomposed; W: a well-distributed set of weight vectors {w 1 , . . . , w N }; T : the neighborhood size.Reproduction: Randomly select two indexes k, l from B i , and then generate a new solution y from x k and x l by using genetic operators.

		Output:
		P: the final approximation to PS.
	9	Mutation: Apply a mutation operator on y to produce y .
	10	Update of z: z j ← min(z j , f j (x i ));	// for j = 1, . . . , k
	11	Update of Neighboring Solutions: For each index j ∈ B i , if
		g(y |w j , z) < g(x i |w j , z), then set x j = y ;

1 z = (z 1 = +∞, . . . , z k = +∞) ; 2 Generate a random set of solutions P = {x 1 , . . . , x N } in Ω; 3 for i = 1, . . . , N do 4 B i ← {i 1 , . . . , i T }, such that: w i 1 , . . . , w i T are the T closest weight vectors to w i ; 5 z j ← min(z j , f j (x i )); // for j = 1, . . . , k 6 while stopping criterion is not satisfied do 7 for x i ∈ P do 8 12 return P = {x 1 , . . . , x N };

1 .

 1 NoEffectCoord. Reset if adding 0.2-standard deviations in any single coordinate does not change m i (i.e. m j equals m i j + 0.2σc j, j for any j = 1, . . . , n). 2. NoEffectAxis. Reset if adding a 0.1-standard deviation vector in any principal axis direction of C i does not change m i . More formally, stop if m i equals m i +0.1σ d j j b j , where j = (τ mod n) + 1 and d j j and b j are the j th eigenvalue and eigenvector of C i , with ||b j || = 1, respectively. 3. TolXUp. Stop/Reset if σ • max(diag(D)) increased by more than 10 4 . 4. ConditionCov. Stop/Reset if the condition number of the covariance matrix exceeds 10 14 .

  ← {i 1 , . . . , i T }, such that: w i 1 , . . . , w i T are the T closest weight vectors to w i ; C i , m i , σ i , τ i , x i , σ init )

		Algorithm 2: MOEA/D-CMA+I
	7 σ , 8 Evolution Strategy Initialize ( p i c , p i
	9 while stopping criterion is not satisfied do
	10	Step 1. Reproduction and Replacement
	11	for i = 1, . . . , N do
	12

Input:

N: the number of subproblems to be decomposed; W: a well-distributed set of weight vectors W = {w 1 , . . . , w N }; T : the neighborhood size. Output: P: the final approximation to PS ; 1 Initialization 2 z = (z 1 = +∞, . . . , z k = +∞) ; 3 Generate a random set of solutions P = {x 1 , . . . , x N } in Ω; 4 for i = 1, . . . , N do 5 z l = min(z l , f l (x i )) ; // for l = 1, . . . , k 6 B i

Table 1 :

 1 Parameters for MOEA/D-CMA

	λ

Note that the sampling is implemented via an Eigen decomposition of the covariance matrix into C i = BDB -1 and then sampling according to m i + σ i BDN(0, I) with B an orthonormal basis of eigenvectors and D a diagonal matrix containing the corresponding positive eigenvalues.