2.1 Multi-Objective Optimization

A continuous MOP can be stated in mathematical terms as:

where F is defined as the vector of the objective functions:

and where each objective fi:ℝn→ℝ, i=1,…,k, is for convenience considered to be continuously differentiable. We stress, however, that this smoothness assumption is only needed for the theoretical considerations but not for the realization of the algorithms. The inequality constraints gi, i=1,…,p, and the equality constraints hj, j=1,…,m, form the feasible set Q, i.e.:

When defining the optimal solution we cannot proceed as for the classical scalar optimization case since F is vector-valued. For MOPs, optimality is defined based on the concept of dominance [²⁹]:

Definition 1 (a) Letx,y∈Qandfx=f(x),fy=f(y)∈ℝkbe its respective images. We say that the vectorxdominatesy(denoted asx≺y) ifffix≤fiy, i=1,…,k,, and there exist an indexjsuch thatfjx<fjy.

Typically, i.e., under certain mild assumptions on the MOP, both the Pareto set and the Pareto front form (k−1)-dimensional sets [¹⁴].

2.2 Solving a MOP

So far, many different classes of methods for the numerical solution of MOPs exist such as mathematical programming techniques [³⁷, ⁶, ²⁶], continuation methods [³¹, ¹⁴, ³⁴, ³², ³⁰, ²³, ²⁴, ²²], or cell-mapping and subdivision techniques [⁹, ¹⁶, ¹², ²⁷, ¹¹].

In this study we are particularly interested in multi-objective evolutionary algorithms [⁵, ⁷] and memetic algorithms [³⁷, ²¹, ⁴⁰, ¹⁵], as they have demonstrated to be efficient in many cases.

In the following we shortly present two different MOEAs that we use within this work as the base algorithms for our GSA-based local search: the decomposition based algorithm MOEA/D [⁴²] and a memetic variant of the well-known algorithm NGSA-II that utilizes gradient information [¹⁹].

2.3 Gradient Subspace Approximation

The Gradient Subspace Approximation (GSA) [³³], is a method proposed to approximate the gradient of a function at a given candidate solution. GSA utilizes the information that is present for each generation of an evolutionary algorithm to compute such an approximation. Here, we present a brief description for the unconstrained case. For the treatment of constrained problems (inequality and equality constraints), the reader is referred to [³³].

Consider the following single objective optimization problem:

where f:ℝn→ℝ is assumed to be differentiable. We say that v∈ℝn is a descend direction of f at x if:

where ∇f(x) represents the gradient of f at x. It is known (e.g., [²⁸]) that the vector with the maximal decay at x is given by the negative of the gradient. This (normalized) vector can be expressed also as the solution of the following auxiliary SOP:

To exploit the neighboring information of the candidate solution x0, consider that r search directions vr are given along with the r directional derivatives of such search directions 〈vi,∇f(x0)〉,i=1,…,r. If we incorporate the neighborhood information such that the most greedy direction exist in v∈span(v1,…,vr), Problem (9) can be rewritten as follows:

where

The following result shows that the solution of (10) is unique under certain assumptions and that it can be expressed analytically.

Proposition 1 Let the search directionsv1,…,vr∈ℝn, r≤n, be linearly independent and:

Then, there exist a single solution ofEquation (10)given by:

and thus,

is the most greedy search direction inspan{vi,…,vr}.

Equation (14), requires that gradient information is given. To avoid such restriction we can approximate such information using the neighbors of the candidate solution. Given the candidate solution x0, a set of neighbors xi, i=1,…,r, and their respective function values f(xi), i=1,…,r, the gradient information can be approximated such that:

In the above equation, the search direction has apparently been chosen as:

For more details and for the extension to equality and inequality constrained problems we refer to [³³].

3.0.1 Applicability of GSA within MOEAs

In population based scalar optimization one typically observes that many individuals are located around the best found individual (denoted here by x0), apart from very early stages of the search. Thus, for the usage of GSA one can expect that sufficient samples are at hand near x0 to approximate ∇f(x0). The situation, however, changes when considering multi-objective optimization problems as for such problems there does typically not exist one single solution, but an entire manifold of solutions. Hence, there is no single ”best” solution any more where other solutions group around. Instead, the solutions are spread along the solution set.

Thus, the natural question arises if there exist enough samples within multi-objective evolutionary algorithms (or any other set or population based algorithms for the treatment of MOPs) such that GSA can be applied successfully. To answer this question empirically, we consider in the following the neighborhood sizes from populations obtained via MOEA/D for four different MOPs with different characteristics. In all cases, we use the 2-norm in decision to define the neighborhood. More precisely, we define the neighborhood for an individual p0 of a given population P as:

however, we stress that in principle any other norm can be chosen.

In our computations, we have used the values δ=0.25, 0.5, and 1. As MOPs we have chosen (see Table 8 for the formulations of all problems used in this work):

Figures 1-4 present the results of our experiments for population size N=100 averaged over 30 independent runs. Shown are the number of neighbors (i.e., the magnitudes of Nδ(p0)), averaged over all runs and all elements of each population generated by MOEA/D. As it can be seen, the magnitudes of the neighborhoods quickly go over 5 which has been detected as an effective value for r in order to obtain a descent direction (and which has also been used for the computations we present in the sequel). These observations are in consensus with other experiments we have done.

Based on this insight we choose our neighbor-hood N(p0) for a given candidate solution simply via selecting the r nearest elements (using the 2-norm) out of the population P. See Algorithm 3 for a pseudocode.

Hereby we assume that the population is given by μ elements, i.e.:

We stress that one potential problem when using MOEA/D is that populations can have duplicate individuals (as a point can be the best found solution for several scalarization problems at the same time). Hence, we have to avoid such solutions as else the distance between the candidate solution and the sample is zero, and GSA is not applicable. More sophisticated neighborhood selection strategies that e.g. take into account the condition number of V or the one MOEA/D is using are subject of ongoing research.

3.1 Approximating the Jacobian

Assume we are given a MOP of the form (1), and that we are given a candidate solution x0 together with the r search directions v1,…,vr∈ℝn that form the subspace.

For every objective fi, the best approximation g˜i(x0) of the gradient ∇fi(x0) at a given point x0 within S can according to the above discussion be computed via:

The Jacobian matrix J˜(x0) of F at x0 is defined by:

Thus, using (20)), a best approximation J˜(x0) of J(x0)

For the gradient free realization, consider the product of the first two matrices on the right hand side of (22) it is:

That is, every entry:

of the matrix ℱ is a directional derivative and can be approximated via the finite difference method described above, and this approximation comes for free in terms of additional function evaluations as the values fi(xj) are already known. Thus, using (22), where ℱ is computed as described, yields a Jacobian approximation without explicitly computing the objectives’ gradients. The algorithm can hence be seen as gradient free.

4.1 GSA within Laras Method

First we investigate one line search method, namely the iteration:

where vL,i is the GSA-approximation of the descent direction (25) at the current iterate xi and where ti is the step size control as described above. As example we consider the unconstrained convex bi-objective problem:

where a1=(1,…,1)T∈ℝn and a2=(−1,…,−1)T∈ℝn. For the value of n we have chosen 10.

First we investigate the influence of the number r of sampling points. Figure 6 presents the result of the following experiment: we have taken x0∈ℝ10 at random and have generated a set of r neighorhood samples from Nδ(x0) for δ=0.15 for r=3, 5, and 7.

The figure shows the points y0:=F(x0), the iterate y1:=F(x0+vL) (i.e., the image of the iterate x1:=x0+vL for the step size t=1) and further the approximations of y1 when using GSA and r samples. As it can be seen, already for r=5 a significant move toward the Pareto front can be obtained, and for r=7 the performance comes quite close to the usage of the exact gradient.

Since the consideration of one step has no significance, we consider an entire trajectory of solutions starting with x0 using r=5. In Figure 7 a sequence of 15 iterations is shown. As it can be seen, the iterations come close to the Pareto front.

As already mentioned above, more theoretical investigations are required to fully understand GSA based local search which we, however have to leave to future research. In the following we focus on the effect of the GSA based local search engines within memetic multi-objective evolutionary algorithms.

4.2 MOEA/D/GSA

For this experiment we use the original version of the MOEA/D algorithm as it was presented in [⁴²]. A comparison using two different state-of-the-art indicators is performed. The Δ2 indicator [³⁵] and the hypervolume indicator [⁴⁴] are selected for the performance assessment. Both indicators measure spread and convergence along the Pareto front to a certain extent. We propose to compare the algorithms (the standalone and the memetic version), after they have spent a certain number of function calls. We perform our comparison when the algorithms reached 30,000 function evaluations for the functions with two objectives and 50,000 for k=3. We stress that MOEA/D is not explicitly aiming at one of any existing performance indicators. Instead, it measures the success via improvements in the scalarization functions.

Hence, along with the two state-of-the-art indicators, we propose a specialized indicator that measuares the averaged scalar value of the whole population (and thus, in a sense the overall success of the MOEA/D variants). In particular, we propose to use the indicator W that is defined as follows:

where xi∈P is the best found individual for the i-th subproblem. We introduce this new indicator, as we have observed that if W decreases (and thus, in average also the best found solutions for all sub-problems of MOEA/D), this does not necessarily yield better hypervolume nor Δp values. In other words, the sub-problems are not matched to these (or any other), performance indicators.

To confirm that the GSA method can really improve a memetic algorithm it becomes necessary to make a comparison on different test functions. In particular, for this memetic strategy we selected two of the state-of-the-art benchmark suites. The first selected set of functions is the modified version of the ZDT functions proposed in [³⁹]. For the second part of the experiments we use the DTLZ benchmark for k=3 [⁸].

Table 1 presents the parameters used for each one of these benchmarks.

We performed an experiment in order to compare the performance of the memetic algorithm based on the GSA as a local searcher. For each comparison we observed the results over 30 independent runs. The experiments were performed over the ZDT test functions and the DTLZ test functions. The effectiveness of the MOEA/D/GSA was measured comparing the memetic algorithm along with two other methods: the base variant of MOEA/D and a memetic strategy based on the Nelder-Mead method. In these experiments we use a modified version of the Nelder-Mead method as follows: we modified the algorithm to incorporate neighboring information on it (as the GSA does). In particular, we incorporate r individuals from the population into the construction of the simplex. In case that r<(n+1), we computed the remaining n−r+1 using the mechanisms proposed in the original Nelder-Mead algorithm. The resulting algorithm is termed here MOEA/D/NM.

Tables 2 to 6 present the results obtained by the three algorithms on the three different performance indicators. The nadir points for the hypervolume indicator are set as follows: (5,5)T for the ZDT test functions, (11,11)T for DTLZ 1-4 and (5,5,5)T for DTLZ 5-7. The tables present the results at a certain stage of the algorithm. That is, we measured the indicators after a certain number of function calls. For the ZDT test function we measure the results after 15,000 function evaluations. Meanwhile, for the DTLZ functions we obtained the results after 30,000 function evaluations. Bold numbers indicate that the algorithm is outperforming the other ones significantly. As it can be seen, MOEA/D/GSA wins in most cases. For instance, for the W indicator, MOEA/D/GSA significantly wins in 8 out of 12 cases, and is only outperformed on ZDT3 by MOEA/D/NM.

Figure 8 presents the computed Pareto fronts for the best individual on each of the ZDT problems.

Here, it is possible to observe that in most of the cases the GSA helps the algorithm to achieve a ”better” approximation along the entire Pareto front.

4.3 IG-NSGA-II/GSA

The next set of experiments is prepared in order to demonstrate that the GSA can be used to handle constrained problems. Since we are replacing the local search method of the IG-NSGA-II we adopt all the parameters from such algorithm. First, we consider some of the function presented in [⁵]. The definitions of such functions can be found in Table 9. The maximization problems presented on the definitions were transformed into minimization problems. The constraints of the problems were also transformed into the form a(x)≤0. We stress that all constraints are relatively easy in the sense that the feasible set is not of a highly complex structure. Table 4 presents the values for the parameters for the memetic algorithm.

The results obtained in Table 5 show that the GSA can improve the results in most of the functions according to the Δ2 indicator. Table 5 also presents the function evaluation used as stopping criteria. Moreover, the table also presents the nadir point used to compute the hypervolume indicator. It is important to mention that some indicator values obtained by the GSA outperformed significantly the standalone algorithm. Taking these results into consideration, now we are in position to confirm the efficiency of the GSA when it is used within a memetic strategy.

As a last experiment we use the constrained CF functions proposed in [⁴³] those constraints can considered to be complex. We perform a similar comparison as in the previous method. We measure the Δ2 and the hypervolume indicators at certain stage of the evolution (i.e. at 30,000 and 50,000 function evaluations). For the experiments we set the nadir point for hypervolume as (5,5)T∈ℝ2 and (5,5,5)T∈ℝ3. Table 7 presents the averaged results measured by the proposed indicators. The results are obtained taking in consideration only the feasible solutions obtained by the algorithms. By such reason the function CF10 is not in the statistical results since neither of the algorithms obtained feasible solutions.

Figure 9 presents some of the Pareto fronts obtained by the algorithms on the CF functions.