A Multiobjective Genetic Algorithm for the Localization of Optimal and Nearly Optimal Solutions Which Are Potentially Useful : nevMOGA

Traditionally, in a multiobjective optimization problem, the aim is to find the set of optimal solutions, the Pareto front, which provides the decision-maker with a better understanding of the problem. This results in a more knowledgeable decision. However, multimodal solutions and nearly optimal solutions are ignored, although their consideration may be useful for the decision-maker. In particular, there are some of these solutions which we consider specially interesting, namely, the ones that have distinct characteristics from those which dominate them (i.e., the solutions that are not dominated in their neighborhood). We call these solutions potentially useful solutions. In this work, a new genetic algorithm called nevMOGA is presented, which provides not only the optimal solutions but also the multimodal and nearly optimal solutions nondominated in their neighborhood. This means that nevMOGA is able to supply additional and potentially useful solutions for the decision-making stage. This is its main advantage. In order to assess its performance, nevMOGA is tested on two benchmarks and compared with two other optimization algorithms (random and exhaustive searches). Finally, as an example of application, nevMOGA is used in an engineering problem to optimally adjust the parameters of two PI controllers that operate a plant.


Introduction
In industry, there are many situations where a design problem turns into an optimization problem.Moreover, it is usual that these problems have conflicting objectives, which generates a multiobjective optimization problem (MOP) ( [1][2][3]).A MOP basically consists of three stages [4]: its definition, the optimization problem (search), and the multicriteria decision-making (MCDM).In an a priori multiobjective approach [5], the stages of optimization and decisionmaking are carried out at the same time, which results in a single solution, since the desired preferences are defined beforehand.On the contrary, in an a posteriori approach [6], the search does not supply a single solution but a set of optimal solutions (Pareto optimal solutions).This procedure is more time-consuming, but it gives the designer information about different solutions and provides them with a better understanding of the problem, which will allow them to make a well-informed decision in the decision-making stage.On the other hand, an optimization problem may have multimodal solutions ( [7][8][9]) or nearly optimal solutions ( [10,11]) that are useful for the designer.These alternative solutions, which have a similar or even the same performance as the optimal solutions, are ignored in a traditional multiobjective approach, although they offer interesting and useful information to the designer.In order not to lose this valuable information, there must be taken into account not only the optimal solutions but also the nearly optimal ones ( [12]).
Including the multimodal and nearly optimal solutions in the decision-making scenario increases the number of alternatives available to the designer.However, this poses two issues.On the one hand, the number of solutions will increase considerably and this may slow the algorithm to the point of becoming inoperative.On the other hand, having a huge amount of alternatives to choose makes the decisionmaking process more complicated.For this reason, it is important in most cases to reduce, in some way, the number of solutions.The way in which this reduction is performed will depend on the criterion or approach that the designer decides to adopt, which, in turn, depends on the application or the context.The existing algorithms usually perform this reduction by means of a discretization of the objective space.But this discretization is carried out without taking into account any consideration of the location of the solutions in the parameter space.The result of proceeding in this way is that nearly optimal solutions which have distinct characteristics could be neglected.Despite this drawback, this method is suitable for some applications.
One possible approach, however, and the one that we assume here, is the following: (1) finding a manageable number of solutions (so that the two issues mentioned above can be avoided) but, at the same time, (2) without neglecting the existing diversity in the characteristics of the nearly optimal solutions.We think that this approach is common in many applications and the existing algorithms do not serve in this case, precisely because they do not pay attention to the characteristics of the solutions.This is why, in this work, we have developed an algorithm to provide the designer with those solutions, among all the nearly optimal solutions, that fulfill that particular criterion, namely, the nondominated in their neighborhood.These solutions are the information that our designer finds more relevant.For example, a designer adopting this approach will not want to be given two nearly optimal solutions which have similar characteristics (neighboring solutions), if one of them is dominated by the other (i.e., worse for at least one of the objectives and not better for the rest), as he or she will logically choose the nondominated one.However, if these two solutions have significantly different characteristics (i.e., they are nonneighboring solutions), then both of them will be interesting for the designer.If the designer has these solutions available, they can analyze them a posteriori (for example, by including new indicators or by considering the physical sense of the solutions), in order to decide which one is the most suitable.From now on, we will call these solutions potentially useful solutions.This does not mean that the rest of nearly optimal solutions are useless, but simply that they are dispensable for a designer who has adopted the approach specified before.These potentially useful solutions are the ones that according to his or her criterion gives them the most relevant information, the information that they want to get.
An example where the specified approach could be the chosen one and, therefore, it is desirable to find the nearly optimal solutions nondominated in their neighborhood, is when one of the objectives cannot be included in the optimization process due to its high computational cost.In this case, it is possible to exclude that objective from the optimization process and to obtain the set of optimal solutions for the rest of objectives.When the optimization process has finished, the excluded objective can be evaluated (only at the optimal solutions previously obtained) and incorporated into the decision-making process.The problem here (if a classical MOP approach is followed) is that, since the excluded objective is ignored in the search, there will be interesting solutions (those with a good performance for the excluded objective and with similar performance for the rest of objectives) which will be missed.The new approach that we are proposing in this paper finds them.Another situation where it may also be interesting to examine nearly optimal solutions nondominated in their neighborhood happens when several objectives are aggregated, which is a way of defining a priori preferences.When a MOP has a lot of objectives, the decision-making stage becomes too complicated.That is why sometimes several objectives are aggregated in groups (for example, performance, cost, and robustness), in order to simplify the problem.However, this approach could miss valuable solutions for the designer.This loss of information can largely be avoided by finding the nearly optimal solutions nondominated in their neighborhood.
So, finding the nearly optimal solutions nondominated in their neighborhood solves the lack of information that typically appears in a classical multiobjective approach.Moreover, taking into account these solutions has two additional advantages: (1) it enables the designer to reassess the design objectives-if there are nearly optimal solutions which display certain characteristics that are absent in the optimal ones, this could mean a poor formulation of the objectives-and (2) it detects the existence of overparameterization (i.e., when the number of parameters is higher than needed, which could result in a man-made multimodality).
In this work, a new algorithm is presented, nevMOGA, which is aimed at finding the sets of optimal solutions and nearly optimal solutions nondominated in their neighborhood.There exist several algorithms in the literature that consider nearly optimal solutions as well ( [13,14]).However, they use these solutions simply as a means to compute an approximation of the Pareto front, whereas nevMOGA discriminates them and returns some of them (the nondominated in their neighborhood) as its outcome.nevMOGA is based on an algorithm called evMOGA, which is described in [15].This new algorithm includes additional parameters to define when a solution is considered to be a nearly optimal solution and the size of a neighborhood.

Background
In this section, the concepts of multiobjective optimization problem and Pareto set are formally introduced, and the problem which we are dealing with is graphically explained (i.e., the problem of finding optimal and nearly optimal solutions nondominated in their neighborhood).
A multiobjective optimization problem (A maximization problem can be converted into a minimization one.For each of the objectives that have to be maximized, the transformation: max f i x = −min −f i x can be applied.)can be stated as follows: Complexity x i and x i are the lower and upper bounds of each component of x.
Definition 1 (dominance [16]).A decision vector x 1 is dominated by any other decision vector for all i ∈ 1, … , m and f j x 2 < f j x 1 for at least one j, j ∈ 1, … , m .This is denoted as x 2 ⪯x 1 Definition 2 (Pareto set).The Pareto set (denoted by P Q ) is the set of solutions in Q which are not dominated by another solution in Q: Definition 3 (Pareto front).Given a set of Pareto optimal solutions P Q , the Pareto front f P Q is defined as Figure 1 shows an example to clarify what we mean by optimal solutions, nearly optimal solutions, and nearly optimal solutions nondominated in their neighborhood in a monoobjective (Figure 1(a)) and a multiobjective (Figure 1(b)) problem.Let us consider the monoobjective case.In theory, all the nearly optimal solutions (grey lines) could be interesting for the decision-maker.However, considering all these solutions makes the decision-making stage too difficult.Moreover, many of these solutions (neighborhood 1 ) will have very similar characteristics to the optimal one (x 1 ), since they are close to it in the parameter space.This, in addition to the fact that they have a worse performance than the optimal one, leads us to consider them, under the criterion assumed, less relevant and in some cases dispensable.The rest of them (neighborhood 2 ) will have similar characteristics to x 2 , which is the best in nei ghborhood 2 and, therefore (always under the specified criterion), they could be dispensable.Consequently, the solutions which provide more relevant information (potentially useful solutions) are x 1 and x 2 , i.e., the optimal solution and the nearly optimal solution nondominated in its neighborhood.Likewise, in the multiobjective case, all the nearly optimal solutions (grey area) could be considered but, again, taking into account all of them would significantly complicate the decisionmaking stage while most of these solutions do not add much relevant information to the process.So, the sets SET n1 (optimal solutions) and SET n2 (nearly optimal solutions nondominated in their neighborhood) provide the most valuable information to the designer without unnecessarily complicating the decision-making process.
An example of the usefulness of finding the nearly optimal solutions nondominated in their neighborhood is when one of the objectives cannot be included in the optimization process due to its high computational cost.The usual way to cope with this difficulty consists of solving the MOP for the rest of objectives and then computing the excluded one only for the optimal solutions.The problem with that (if a classical MOP approach is followed) is that there will be nearly optimal solutions with an outstanding performance with respect to the excluded objective (and therefore worth considering) that will be ignored.On the other hand, if one is not willing to miss all those interesting solutions and decides to take into account all the nearly optimal solutions, then the time-consuming objective will be evaluated at a vast number of solutions, many of which are not even worth considering, as they will have similar characteristics to the optimal solutions (because they are neighboring solutions).In summary, approaching this kind of problems with a classical strategy leads to either neglecting potentially useful solutions or having to evaluate a time-expensive objective at many dispensable solutions, whose inclusion in the decisionmaking stage will, furthermore, unnecessarily complicate the Optimal and nearly optimal solutions f(x) Neighborhood 3

Optimal and nearly optimal solutions
Optimal and nearly optimal solutions Figure 1: A monoobjective example (a) and a multiobjective example (b).In the monoobjective problem, the optimal solution is x 1 (blue triangle) and the nearly optimal solution nondominated in its neighborhood is x 2 (green triangle).In the multiobjective problem, the optimal solutions are the ones in set SET n1 (blue front) and the nearly optimal solutions nondominated in their neighborhood are the ones in SET n2 (green front).
3 Complexity decision.With our new approach, these problems can be overcome.Let us look at it with a concrete example.
Let us suppose a MOP with three objectives f x = f 1 x f 2 x f 3 x to minimize.Table 1 shows five possible solutions to the problem (x 1 to x 5 ).In a classical MOP, all these solutions belong to the Pareto set (except either x 1 or x 2 , as they are multimodal solutions and, therefore, one of them would be discarded).Now, assume that the objective f 3 is computationally impracticable, i.e., its incorporation into the optimization process is impossible.For this reason, the search stage has to be solved for the first two objectives f x = f 1 x f 2 x and then f 3 will be evaluated only at the solutions previously obtained.First, note that in a classical MOP, only the Pareto front is obtained, so either x 1 or x 2 (multimodal solutions) and x 3 would be discarded, although x 1 , x 2 , and x 3 seem to be good solutions, since they are nearly optimal solutions in different neighborhoods and thus worth considering.
With our new approach, x 1 , x 2 , and x 3 would also be found, because they are nondominated in their neighborhood; in other words, they have significantly different parameter values.x 4 would be discarded since it is dominated by x 1 and they are neighboring solutions, i.e., they are expected to have similar characteristic.x 4 has a better performance in f 3 than x 1 ; however, this improvement is insignificant (if the neighborhoods are well defined) as they are very similar solutions (neighboring solutions).Through this new procedure, from all the neighboring solutions, only those with the best performance over the objectives considered in the optimization process (in this example, f 1 and f 2 ) would be selected.The solution x 3 displays a better value for f 3 than x 1 and x 2 (which are in different neighborhoods) but worse values for f 1 and f 2 (this is why it is a nearly optimal solution) and, therefore, in a classical MOP, it would be discarded, although its analysis may be interesting for the decisionmaker.Logically, x 5 would be ruled out since it is much worse for f 1 and f 2 , although it has a better value for f 3 .With the proposed approach, f 3 is only evaluated at x 1 , x 2 , and x 3 but not at x 4 (x 1 is preferred instead) and x 5 (which is discarded, since it is not a nearly optimal solution) and, consequently, the computational burden is alleviated.Therefore, thanks to this new approach; it is possible to reduce the computational cost while obtaining the potentially useful solutions for the designer.
Another example that shows the usefulness of obtaining also the nearly optimal solutions nondominated in their neighborhood can be found when several objectives are aggregated, which is a common procedure when there are a great number of them.Let us suppose that a MOP originally has four objectives f x = f 11 x f 12 x f 21 x f 22 x to minimize.Let be x 1 to x 4 (Table 2) the four solutions to this MOP.In a classical MOP, all these solutions would be selected to form the Pareto set.
If the designer wanted to reduce the number of objectives, a possibility would be to add f 11 and f 12 to obtain f 1 (f 1 = f 11 + f 12 ) and f 21 and f 22 to obtain f 2 (f 2 = f 21 + f 22 ).But now, only x 1 or x 2 (multimodal solutions) would be the optimal solution (in a classical MOP).Let us suppose that the algorithm chooses x 1 and discards x 2 .In this case, x 2 and x 3 are discarded, although they are interesting solutions for the designer.However, if the algorithm keeps also the nearly optimal solutions nondominated in their neighborhood, x 2 and x 3 would not be excluded with the aggregated objectives.Finally, x 4 presents a poor performance for both f 1 and f 2 (clearly dominated) and, therefore, it does not seem to be a good solution.
In conclusion, when the objectives are aggregated, a classical MOP finds only one solution from the three that are interesting and the other two are missed.Therefore, finding the nearly optimal solutions nondominated in their neighborhood allows the designer to aggregate objectives (and, in this way, simplifying the decision-making stage) without losing potentially useful solution.
In summary, as seen in the examples, the multimodal solutions and the nearly optimal solutions which are significantly different from those which dominate them are solutions that can be potentially useful for the designer during the decision-making stage.This additional information is not provided to the decision-maker when these kinds of MOPs are approached in a classical manner.

Materials and Methods
In this section, we present in detail a novel algorithm, nev-MOGA, which is the main contribution of this work, as well as the metric and benchmarks that have been used to demonstrate its correct functioning.Let us begin with the new set of interest and its discretization, which is performed by nevMOGA.
3.1.New Set of Interest.The nevMOGA algorithm is aimed at finding not only the optimal solutions but also the nearly optimal solutions nondominated in their neighborhood.In the following, both sets of solutions are defined.
Table 1: Analysis of a MOP example with three objectives (f 1 x to f 3 x ) and five possible solutions (x 1 to x 5 ). Solutions 0.9 0.199 1.119 0.9 0.9 1.8

Complexity
Definition 4 (−ϵ-dominance [17]).Define ϵ = ϵ 1 , … , ϵ m as the maximum acceptable performance degradation.A decision vector x 1 is −ϵ-dominated by another decision vector Definition 5 (−ϵ-Pareto set [18]).The −ϵ-Pareto set (denoted by P Q,ϵ ) is the set of solutions in Q which are not −ϵ-dominated by another solution in Q: Definition 6 (neighborhood).Define n = n 1 , … , n k as the maximum distance between neighboring solutions.Two decision vectors x 1 and x 2 are neighboring solutions Definition 7 (n-dominance).A decision vector x 1 is n-dominated by another decision vector x 2 if they are neighboring solutions (Definition 6) and x 2 ⪯x 1 .This is denoted by Definition 8 (n-Pareto set).The n-Pareto set (denoted by P Q,n ) is the set of solutions of P Q,ϵ which are not n-dominated by another solution in P Q,ϵ : The set P Q,n (optimal solutions and nearly optimal solutions nondominated in their neighborhood) is the new set of interest, i.e., what nevMOGA has to find.In order to clarify its nature, we will use a graphical example (see Figure 2).In this example, the set of interest is formed by the union of P Q and P Q ′ .P Q ′ is a set of solutions nondominated in their neighborhood which belong to the set P Q,ϵ .P Q ″ is a set of solutions nondominated in their neighborhood which do not belong to the set P Q,ϵ .x 4 is not a nearly optimal solution, and, therefore, it will be discarded since it does not have the desirable performance (P Q,ϵ ).x 2 is n dominated (i.e., dominated by a neighboring solution) by x 1 and, as a result, this alternative provides less relevant information and, for this reason, is discarded.x 3 is dominated by nonneighboring solutions but not by neighboring solutions, so it belongs to P Q,n , and the same can be said of any solution in P Q ′ .x 1 belongs to the optimal set P Q .In summary, nevMOGA searches for the set P Q,n , in our example, P Q ∪ P Q ′ .
The set P Q,n clearly depends on the values given to the new parameters ϵ (Definition 4) and n (Definition 6).Given that in most problems, the objectives have a physical sense, it is possible to decide a priori how much performance we are willing to lose and to define ϵ consequently.Similarly, when the decision variables x have a physical sense, the parameter n can also be set intuitively.By setting n, the designer establishes to what extent the two solutions are considered "similar"; in other words, setting n defines the neighborhood.
If the decision variables do not have a physical sense and it is not evident how to choose n, it is still possible to follow a simple procedure (see Figure 3) which will help us to do it if ϵ has been previously set.First, a reference solution x R is chosen.Second, in the objective space, a rectangle is defined with the center at f x R and twice ϵ in width and in height.Then, starting from x R (in the search space), the value of each variable decision is increased and decreased independently, until the solution leaves the rectangle which was previously defined.Finally, each element of the vector neighborhood n is set to the minimum distance between ϵ and the first solutions (one for each direction) which left the mentioned The set of optimal solutions is P Q , and the set of nearly optimal solutions nondominated in their neighborhood is P Q ′ .x 2 is discarded since it is n dominated (i.e., dominated by a neighboring solution) by x 1 (neighborhood x 2 ).The solutions in P Q ′ ′ are discarded since their performances fall out of the area of acceptable performance (grey area).
5 Complexity rectangle (n 1 and n 2 in Figure 3).Through this procedure, we manage to quantify what excursion is needed in the search space in order to get significant changes in the objective space and, in this way, we are able to define n when no other criterion is available.

Discretization of the New Set of Interest.
Another issue is that P Q,n may contain infinite solutions, so it is necessary to discretize it and to obtain a finite set of solutions (P * Q,n ) which are adequately distributed.In order to carry out this task, the concepts of box, n_box, and their associated dominances have to be defined beforehand.

Definition 9 (box [19]). Let δ
for δ i > 0, and where f max Definition 10 (box dominance [19]).Given two decision vectors x 1 and x 2 whose boxes are box x 1 and box x 2 , respectively, Definition 11 (n-box dominance).Given the neighboring decision vectors x 1 and x 2 whose boxes are box x 1 and box x 2 , respectively, The algorithm only keeps non-n-box-dominated solutions.A box cannot contain two neighboring solutions.When a neighboring solution is found in the same box, the nearest solution to the ideal corner (the lower left one) will be chosen (see Figure 4, d 1 and d 2 ) and the rest will be discarded.Figure 4 shows how the algorithm performs the discretization.x 1 is n-box dominated (i.e., box dominated by a neighboring solution) by x 2 and, therefore, it is discarded.x 3 is box dominated but not n-box dominated (since x 2 and x 3 are nonneighboring solutions), so it is an alternative, as it is a nearly optimal solution.x 4 , x 5 , and x 6 are in the same box.Two neighboring solutions cannot be in the same box, so either x 4 or x 5 has to be eliminated.In this case, x 4 is eliminated because it has a greater distance to the ideal corner than x 5 (d 2 < d 1 ).x 6 is a good alternative because it is a nearly optimal solution, it is not n-box dominated and it does not belong to either the neighborhood of x 4 or x 5 .
The algorithm searches a discretization of the set of interest.This new discrete set is P * Q,n .The maximum number of solutions (in the worst case ).Similarly, P * Q and P * Q,ϵ are the discrete sets of P Q and P Q,ϵ , respectively.
x j − x j n i 7 3.3.Description of nevMOGA.The algorithm nevMOGA is based on the algorithm evMOGA [15].The main difference between them is that nevMOGA computes an additional population, namely, the set of nearly optimal solutions Pareto set x R (reference solution) Nearly optimal olutions Solutions leaving the rectangle  (a) Figure 3: Determination of the vector neighborhood n when no other criterion is available.Each element of the vector n (n 1 and n 2 ) is set to the minimum distance between x R and the first solutions (one for each direction) which leave the rectangle.The algorithm manages four populations (see Figure 5): (1) P t : by means of the population P t , the search space is explored, with the aims of obtaining the solutions of P Q,n and having diversity in the set of solutions.The number of individuals in this population is constant and equal to Nind P (2) Front t is the archive where P * Q is stored, i.e., a discrete approximation of the Pareto front.The size of this population varies but is always less than or equal to a given maximum size which depends on the number of boxes previously defined by the user (3) Subfront t is the archive where P * Q,n /P * Q is stored, i.e., the nearly optimal solutions nondominated in their neighborhood.Its size is variable but bounded, depending on the number of boxes (4) G t is an auxiliary population where the new individuals generated by the algorithm in each iteration are stored.The number of individuals of this population is Nind G , which must be multiple of 4 Now that the archives front t and subfront t have been presented, it is possible, with the help of the prior definitions, to establish the conditions that a solution must fulfill in order to enter them.This is the aim of the following two additional definitions: Definition 12 (inclusion of x in front t ).Given a solution x and the archive front t , x will be included in front t if and only if where dist z is the distance from z to the ideal corner (lower left) of the box (see Figure 3, d 1 and d 2 ).Additionally, if x is included in front t , then all the solutions x * that fulfill the following condition will be eliminated from front t : Neighborhood x 1 x 5 x 1 x 6 x 3 x 4 x 2 (b) 7 Complexity Then, it will be determined whether any of these solutions x * which have been eliminated from front t is included in subfront t or not, by applying Definition 13.Furthermore, all the solutions x′ that fulfill the following condition will be eliminated from subfront t : Definition 13 (inclusion of x in subfront t ).Given a solution x such that x ∉ f ront t , x will be included in subfront t if and only if Additionally, if x is included in subfront t , then all the solutions x * that fulfill the following condition will be eliminated from subfront t : In order to maintain the diversity of solutions in P * Q,n , it is necessary that the population P t be diverse too.For this, P t is permanently ordered by using the niche count, which is an indicator of how densely populated a solution is in P t (Definition 14).In this way, during the processes of selection and substitution, the solutions in the less populated areas have a higher probability of being chosen to generate new individuals, whereas those in the more populated areas have a higher probability of being substituted.Definition 14 (niche count [20]).
where sh d ij is the sharing function, which gives a measure of how similar two elements in a population are and is defined by The algorithm pseudocode is the following: Lines 3, 6, 8, 10, 13, 14, and 15 of the previous pseudocode constitute the changes which have been added to evMOGA, the algorithm on which ours is based.Next, the most important steps in nevMOGA are explained: x G1 and x G2 are obtained by crossing over x P1 and x F 6: else 7: x G1 and x G2 are obtained by mutating x P1 and x F 8: end if 9: Choose x SF ∈ Sub-Front (t) ⊳ Random selection 10: Choose x P2 ∈ P (t) ⊳ According to an exponential distribution (μ = Nind P /10) 11: u≔ random(1) 12: if u>P c /m then 13: x G3 and x G4 are obtained by crossing over x P2 and x SF 14: else 15: x G3 and x G4 are obtained by mutating x P2 and x SF 16: end if Algorithm 4: Create population G t 9 Complexity probability of being chosen) and another from P t (x P ) (where each one has a different probability of being chosen).In effect, the solutions of P t are ordered according to their niche count (from least to greatest) and its individual is chosen according to an exponential distribution, which makes it more likely for a solution with a lower niche count to be chosen.This procedure favors the uniform distribution of the solutions (2) A random number u ∈ 0 1 is generated.If u > P c /m (probability of crossover/mutation), then a crossover will be performed (step 3); otherwise, a mutation will be performed (step 4) (3) x P and x F are crossed by means of an extended intermediate crossover operator, which produces two new individuals x h1 and x h2 : The parameter α i t is a random value uniformly distributed which belongs to interval −d t 1 + d t , and d t is a parameter which is adjusted by using an exponential decreasing function, as in simulated annealing [19]: Figure 6 shows d t from t = 0 to t = t max − 1, where d 0 = d ini and d t max − 1 = d fin .
(4) x P and x F are mutated by using a random mutation with Gaussian distribution: where the variances β1 i t and β2 i t are expressed as a percentage of x i max − x i min .These variances are calculated by a function which is similar to the one previously used for the parameter d t : (5) x h1 and x h2 are included in G t These five steps are repeated again from step one but now with a solution from subfront t instead of one from front t .The whole process is executed again and again until G t is full.In order to set the parameters which appear in the equations of the steps 3 and 4, the default values suggested by [19] for the original algorithm (evMOGA) are taken.
(iv) Line 14 (Algorithm 5).In the updating of P t , it is determined whether a new individual from G t (x G ) substitutes any other existing individual in P t or not.The search for the individual to be substituted (x P ) starts from an individual chosen by an exponential distribution which is the inverse of the one used for the choice of x P in the creation of G t (line 10, Algorithm 4) and, therefore, now the more populated solutions will be more likely to be chosen.This is aimed at achieving a more uniform distribution of the solutions.During this updating, there may occur three distinct cases: (1) The new individual x G is ϵ dominated by some member of front t .In this case, an individual x P dominated by x G will be searched to be substituted (2) The new individual x G is not ϵ dominated by any member of front t , and there exists an individual x P dominated by x G which is ϵ dominated by some member of front t .In this case, it will be substituted (3) The new individual x G is not ϵ dominated by any member of front t , and there does not exist an individual x P dominated by x G which is not ϵ dominated by any member of front t .In this case, the initial solution x P is chosen.If x P and x G are neighboring solutions, a solution n dominated by x G will be searched to be substituted.If they are not, a solution n dominated by x P will be searched to be substituted.In this last case, if that solution is not found, x G will be randomly substituted by a solution from the neighborhood of x P 3.4.Performance Metrics.Next, we present one metric which will serve to assess the performance of the different algorithms to be compared.This metric measures both the convergence and the diversity between two sets (the outcome 10 Complexity set and the target set).The outcome set of each algorithm is the set of solutions that it returns, whereas the target set results from discretizing the search space with a fine grain.We use this metric to measure the performance of the algorithms in the objective space (convergence toward the front) and in the parameter space (diversity of solutions).

Averaged Distance Hausdorff.
A single indicator is used to measure the convergence and diversity between two sets (averaged distance Hausdorff Δ p , [21]).The sets Y and X have n y and n x solutions, respectively.The Δ p is calculated by where This indicator is the average of the Hausdorff distance, and we use it here in order to measure the performance between an outcome set H and a target set P * Q,n .This performance is measured both in the objective space and in the parameter space.In this way, it is possible to measure convergence toward the front and diversity of solutions.We use p = 2.It is desirable that this indicator has the smallest possible value.Ideally, when H = P * Q,n , then Δ p = 0.
3.5.Benchmarks.Our algorithm nevMOGA has been tested on two distinct benchmarks.In this section, we describe them.

Benchmark 1.
The first benchmark considered in this work corresponds to an academic example [18] and is stated as follows:

21
where  The second benchmark is an adaptation of the modified Rastrigin [22].This benchmark has been modified in order to turn it into a MOP with one optimal set and nearly optimal solutions lying in different neighborhoods.This is its formulation:

Results and Discussion
In this section, we present the results of nevMOGA on the two benchmarks introduced and on a control design problem.These results will be analyzed on the basis of the metric which was introduced in Section 3.4.In addition, we compare nevMOGA with two other algorithms, namely, random and exhaustive searches.In order to make this comparison possible, all three algorithms evaluate the objective functions the same number of times.Although there do exist algorithms which search the set of nearly optimal solutions ( [17,18,23]), none of them has the same set of interest as nevMOGA (P * Q,n ).For this reason, nevMOGA cannot directly be compared with any of them.In order to obtain statistical results, each algorithm is tested 50 times on each benchmark.For the exhaustive search case, the evaluated points are uniformly distributed over the search space, forming a hypergrid.In order for it to obtain different results, the hypergrid is slightly lifted each time.
4.1.Benchmark 1.In this benchmark, since the global Pareto set and the local Pareto sets are known (Section 3.5.1),we can set the parameters ϵ and n more easily.Here, we set ϵ = 0 15 0 15 so that the eight sets described in (25) are regarded as nearly optimal solutions.We get n by applying the procedure described in Section 3.1 with a reference solution x R = 0 0 and the parameter ϵ previously set, resulting in n = 0 13 0 38 .Once these parameters have been determined, the set to be found (target set P * Q,n , see Figure 7) is Figure 8 and Table 3 show the statistical results of the three algorithms (each one is tested 50 times on this benchmark) that have been compared (nevMOGA, random search, and exhaustive search) using the metric Δ p .Each algorithm evaluates the objective function 5000 times in each test run.The outcome sets of the three algorithms are compared with the target set P * Q,n .In Figure 8(a), there are the results against Δ p f H , f P * Q,n , whereas in Figure 8(b) against Δ p H, P * Q,n , using boxplots and density estimation.The values yielded by nevMOGA for the metric in the objective space and in the parameter space are good (near 0) and better than those provided by random and exhaustive searches.For the random search, the approximation quality to the front varies enormously.For the exhaustive search, the results display less variability, although it performs worse than nevMOGA in both spaces.Moreover, a multiple comparison test ( [24,25]) leads to the conclusion that the results of the three algorithms are significantly different with respect to the metric in both spaces.For these reasons, it can be concluded that nevMOGA significantly outperforms random and exhaustive searches on benchmark 1.
In order to graphically compare the three algorithms, a representative solution (from the 50 test runs) of each one has been plotted in Figure 9.This representative solution has been chosen as the one whose values of Δ p H, P * Q,n and are closest to the average values of that metric in both spaces.The outcome set obtained by nevMOGA has 70 solutions, whereas those obtained by random and exhaustive searches have 78 and 48 solutions, respectively.The great number of solutions found by nevMOGA in addition to the results shown in Figure 9 proofs that nevMOGA is capable of characterizing the target set P * Q,n better than its competitors (random and exhaustive searches) on benchmark 1.This result, plus the fact that nevMOGA achieves lower values   3.   13 Complexity of the metric Δ p , leads us to conclude that nevMOGA accomplishes a better general performance on this benchmark.This benchmark has also been solved by [18], where the set P Q,ϵ is characterized.The algorithm tested in that work is called P Q,ϵ -MOEA.A comparison of the results of nev-MOGA with those of P Q, -MOEA shows that nevMOGA obtains a significantly smaller number of solutions (approximately four times smaller) without missing any  Q,n in the parameter space (right) and in the objective space (left).This representative solution has been chosen as the one whose values of Δ p are closest to the average values of that metric in the objective space and in the parameter space.14 Complexity neighborhood containing nearly optimal solutions.This was our goal, namely, finding only those solutions that provide the most relevant information, and this with two aims: (1) not to slow down the algorithm and (2) to simplify the decision-making stage.

Benchmark 2.
The statement of this second optimization problem was presented in Section 3.5.2.First, and with the aim of setting the two parameters of nevMOGA (ϵ and n), a tentative search for solutions is carried out within the defined range.This leads us to choose ϵ.Once ϵ = 7 7 0 3 has been set, we take a reference solution x R = 0 66 0 5 and follow the procedure detailed in Section 3.1 to obtain n.This results in n = 0 15 0 15 .The selection of ϵ and n leads to a particular target set P * Q,n , which is shown in Figure 10.As can be seen in the figure, P * 4Q is dominated by the set P * 3Q but they are in distinct neighborhoods.The set P * 4Q does not have any solution dominated in its neighborhood and, therefore, it is part of the target set P * Q,n .Figure 11 is similar to Figure for benchmark 2. Each algorithm is run 50 times and can evaluate the objective function 5000 times in each test run.Again, nevMOGA achieves better results than random and exhaustive searches with regard to the metric Δ p , since it reaches values which are closer to zero.As we did in benchmark 1, we have also performed a multiple comparison test on these results.This comparison shows that the results obtained by the three search strategies are significantly different in both spaces.For these reasons, it can be concluded that nevMOGA  4. 15 Complexity significantly outperforms random and exhaustive searches on benchmark 2.
Figure 12 shows a representative solution (chosen in the same way as we did in benchmark 1) of each algorithm (nev-MOGA, random, and exhaustive searches).A glance at the figure is sufficient to show that nevMOGA accomplishes a better approximation to the target set P * Q,n .Moreover, the outcome set obtained by nevMOGA has 73 solutions, whereas those obtained by random and exhaustive searches have 62 and 48 solutions, respectively.These two last facts (a better approximation to the target set P * Q,n and a greater number of solutions) confirm the conclusion drawn before, when the results were statistically examined, namely, that nevMOGA outperforms both random and exhaustive searches.

Multiobjective Control Design Problem.
In this section, we present the application of nevMOGA to a real design problem in the field of control engineering.The problem consists of optimizing the performance of two PI controllers, which operate a multivariable plant with two inputs and two outputs.
The model of the plant is [26] Y s = y 1 s where y 1 and y 2 are the outputs of the plant, and u 1 and u 2 are the inputs (control actions).The time constants are all in seconds.
The control structure (also taken from [26]) is defined as a multiloop PI control which uses an off-diagonal pairing scheme, i.e., output y 1 is controlled by u 2 and y 2 by u 1 .That is to say, where Kc 1 and Kc 2 are the proportional gains, Ti 1 and Ti 2 are the integral time constants, and e 1 = r 1 − y 1 and e 2 = r 2 − y 2 are the errors, where r 1 and r 2 are the setpoints for y 1 and y 2 , respectively (see Figure 13).
The optimization problem is stated as follows: where dt, dt, We launch two different simulations, sequentially.In the first one, a unitary step is applied to r 1 , whereas r 2 is kept to zero.In the second one, the unitary step is now applied to r 2 and r 1 is kept to zero.In each simulation, the metric ISE (Integral Squared Error) is calculated for both outputs y 1 and y 2 .Our first objective f 1 is then defined as the sum of the values of the ISE metrics of each simulation.So, the objective f 1 measures the setpoint tracking performance of the control system.The second objective f 2 is calculated in a similar manner but using the control action signals instead of the errors.Therefore, f 2 measures the control effort.Note that each objective (f 1 and f 2 ) can be seen as a sum of 4 subobjectives (see 32).
Once the optimization problem has been defined, the two parameters of nevMOGA (ϵ and n) have to be set.It is important to choose these parameters correctly, especially n.First of all, we have to decide what performance loss (with respect to  Figure 12: A representative solution of each one of the three algorithms is compared with the target set P * Q,n in the parameter space (right) and in the objective space (left).This representative solution has been chosen as the one whose values of Δ p are closest to the average values of that metric in the objective space and in the parameter space.17 Complexity the Pareto set) we are willing to accept (i.e., the parameter ϵ), in other words, what we are still going to regard as a nearly optimal solution.Choosing a high value would lead to a great number of possible solutions, whereas a low value would provide very few.In this application example, although the design objectives do have a physical sense, it is not self-evident how to choose ϵ.In order to solve this difficulty, we have proceeded in the following way.First, we have taken the values of f 1 and f 2 corresponding to the control designed in [26] as a reference solution.Then, we have pinpointed that solution on the Pareto front (which we have taken from [27]).Finally, on the basis of that information, we have decided to choose ϵ = 10 5 .Note that there is an unavoidable subjective element concerning the way in which ϵ is chosen.Now, we have to assign a value to n, which is the most critical parameter of the algorithm.On the one hand, choosing an excessively low value would make most solutions be nonneighboring solutions, which in turn would cause the algorithm to tend to obtain the whole set P * Q,ϵ .On the other hand, choosing an excessively high value would result in many solutions being neighboring solutions, which would cause the algorithm to tend to obtain only the set P * Q .This is why it is crucial that n be chosen carefully, supporting the choice with good reasons.Next, we describe how we have set n. First, we choose a reference controller x = − 0 36,7 86,0 64,628 4 (taken from [26]).Second, starting from the solution corresponding to this reference controller, the value of each of its parameters is increased and decreased independently, until reaching significant changes in the system response (see Figure 3).In this way, we finally get n = 0 25 50 0 1 200 ; these are the values from which the system response is significantly different from the one provided by the reference controller.
As mentioned before, the objective functions include several terms which measure the performance and the control effort of each control loop (as if there were 4 subobjectives added for each objective).For this reason, there may be solutions with the same objective value which are, however, completely different.This way of defining the objective functions is useful when the number of control loops is high, since it reduces the number of objectives and simplifies the decision-making stage.Yet, this approach has a shortcoming, namely, many interesting solutions may be Nearly optimal solutions (subfront(t)) Optimal solutions (front(t)) Figure 14: Set P * Q,n found by nevMOGA.The optimal solutions (set P Q ) are shown in red and the rest (nearly optimal solutions nondominated in their neighborhood) in blue.x 1 , x 3 , and x 5 belong to P Q (Pareto optimal solutions), whereas x 2 , x 4 , and x 6 are nearly optimal solutions.18 Complexity neglected, because they may become nearly optimal solutions or multimodal solutions.Figure 14 shows the set of solutions P * Q,n found by nevMOGA (for a given discretization of the objective space), where x 1 , x 3 , and x 5 are the optimal solutions, and x 2 , x 4 , and x 6 are the nearly optimal solutions.In order to show the decision variables, we employ a visualization tool called level diagrams ( [28,29]).For the ordinate axis, we have used the 2-norm.For the graphical representation of the objectives, we have not used this tool though, since there are only two and, in this case, plotting them on a single graph, where they can directly be compared, is preferable.Let us draw attention to the decision variables x 1 and x 2 (Figures 14(b)).Note how nevMOGA has been able to identify new areas in the search space containing solutions which may be interesting for the decision-maker and that would have been neglected in a traditional optimization approach, where only the set of Pareto optimal solutions P Q is characterized.
Figure 15 shows the system responses for the solutions x 1 (optimal solution) and x 2 (nearly optimal solution), which have very similar values of f 1 and f 2 (see Figure 14(a)).These two solutions lie in distinct neighborhoods, and their corresponding responses are different too.This should not cause surprise since their corresponding controllers are significantly dissimilar from each other.
A more careful comparison of these two transient responses shows that x 1 presents a much better response than x 2 with respect to y 2 when the step is applied to r 1 , but x 2 presents a better response than x 1 with respect to y 2 and y 1 (better in the sense that it displays less overshoot) when the step is applied to r 2 .Therefore, it is not clear at all which solution should be chosen.In these cases, the decision-maker will have to include new criteria or preferences, in order to make the decision.Anyway, note that the nearly optimal solution x 2 could ultimately be preferred over the optimal one x 1 .This evidences that it would not have been wise to discard x 2 in the search stage and proofs the convenience of finding the nearly optimal solutions besides the optimal ones.This situation is very common when the objective functions are defined as an aggregation of several subobjectives, as is the case here.Additionally, it is worth noting that if this MOP had been addressed with a traditional approach and eight objectives (one for each term in 32, without adding objectives)-which is the ideal approach if one wants to find all the optimal solutions, although impractical in the decision-making stage-then x 2 would not be dominated by x 1 and both x 1 and x 2 could belong to the Pareto front.This is an additional evidence of the relevance of x 2 , which nevMOGA has been able to identify.
Next, we assess the robustness of every solution found by the algorithm, that is to say, how well each controller maintains its performance when identification errors are contemplated.This analysis is typical in control design.In order to do that, first, we get 50 sets of plant parameter (gains, time constants, and time delays) by generating random values for each parameter within a range of ±20%.Each of this sets represents a deviation from the nominal plant.Then, we evaluate the performance (f 1 and f 2 ) of each controller for each one of those deviations.Including an additional objective which measures robustness and computing it in the optimization stage, along with the two other objectives, would have been very expensive computationally.This is why this robustness analysis is shifted to the decision-making stage.As an example of the results of this robustness analysis, Figure 16 shows the boundaries of the performance degradation (determined by the values of f 1 and f 2 when they are evaluated for each one of the 50 sets of plant parameters) experienced by two controllers (x 3 and x 4 ).Note that x 4 (a nearly optimal solution) experiences less degradation than x 3 .This is a very valuable information and could make the decision-maker opt for x 4 instead of x 3 .

Complexity
Let us look at one last proof of the convenience of considering the nearly optimal solutions.If we compare x 5 (optimal solution) and x 6 (nearly optimal solution), we see that they have a similar performance with respect to both f 1 and f 2 .However, x 6 has lower values (in absolute value) of the gains Kc 1 and Kc 2 .From a design perspective, this fact may be an advantage, since this solution will present less sensibility to any measurement noise in the outputs of the plant and, therefore, it could be preferred over the optimal solution.

Conclusions
In this paper, we have presented a new multiobjective genetic algorithm, nevMOGA, which is capable of characterizing not only the Pareto set but also the set of nearly optimal solutions nondominated in their neighborhood.We think that this algorithm will be useful for designers who want to get a manageable number of solutions without neglecting the existing diversity in the characteristics of the solutions.In comparison with other algorithms that only return the Pareto front, our algorithm provides, in addition, new alternatives to be analyzed in the decision-making stage.Given that these solutions have a very similar performance to the optimal ones and significantly different values in the parameter space, they will be regarded by a designer adopting the specified approach as solutions worth studying.This is why obtaining these solutions, which we have called potentially useful solutions, can be very valuable in many multiobjective optimization problems.
nevMOGA has been compared with two other algorithms (random search and exhaustive search) on two benchmark problems, in order to assess its performance.It has not been compared to any other because no genetic algorithm of similar characteristics (one which is capable of finding nearly optimal solutions nondominated in their neighborhood) is available in the literature.In order to perform this comparison, the metric Δ p has been employed, which measures the approximation quality of the outcome set of an algorithm to the target set P * Q,n (the set of solutions to be found) in the objective space and in the parameter space.In addition, a statistical analysis has been conducted.The results reveal the outperformance of nevMOGA in characterizing the set of optimal and nearly optimal solutions in both spaces.
nevMOGA has also been applied to a real design problem in the field of process control engineering.This design problem consists of optimally adjusting the parameters of two PI controllers (K c and T i ) which operate a plant with two inputs and two outputs.The control performance is evaluated by two widely used, conflicting metrics: f 1 , which measures the setpoint tracking performance of the control system, and f 2 , which measures the control effort.Each of these objective functions has been formulated as a sum of four subobjectives (see 32), which means that eight different objectives have been incorporated in the problem definition.Dealing with eight independent objective functions would have complicated enormously the analysis of the Pareto front and, consequently, the decision-making stage.This is why the eight objectives have been aggregated to form only two objective functions, whose Pareto front can be plotted on a twodimensional graph and, thus, easily analyzed.
Aggregating objectives has, however, a major drawback if the optimization algorithm searches only for the optimal solutions: the resulting set of Pareto optimal solutions will lack many solutions that may be interesting and useful for the designer and, in some cases, even preferable to the optimal ones.nevMOGA, as it searches for nearly optimal solutions nondominated in their neighborhood too, is capable of finding those valuable solutions which are neglegted by a classical optimization algorithm.This means that the decision-maker will have a greater number of alternatives, all of them worth considering.For instance, if we compare x 2 (a nearly optimal solution) and x 1 (an optimal solution), we see that they have similar performances (f x 2 ≃ f x 1 ) but their corresponding system responses are different, since they lie in distinct neighborhoods.The designer could choose the nearly optimal solution x 2 instead of the optimal solution x 1 , since the former achieves a better system response when the step is applied to the setpoint r 2 .A comparison between x 5 (optimal solution) and x 6 (nearly optimal solution) leads to the same conclusion: they have similar performances (f x 6 ≃ f x 5 ) but x 6 has lower values of Kc 1 and Kc 2 (the gains of the PI controllers) than x 5 and, therefore, x 6 will be less sensitive to the measurement noises in the outputs, which is a considerable advantage in some applications.This could make the designer prefer x 6 to x 5 .
For this same design problem, a robustness analysis has been performed on the optimal and nearly optimal solutions obtained in the optimization stage.The results of this analysis show how a nearly optimal solution (x 4 ) can be more robust than an optimal one (x 3 ) while providing a similar performance with respect to the objective  (optimal solution) and x 4 (nearly optimal solution) when the plant parameters vary in ±20% from their nominal values.It must be highlighted that a traditional MOP approach would not have found the solution x 4 .20 Complexity function (f x 3 ≃ f x 4 ).This is due to the fact that they are nonneighboring solutions and, therefore, likely to present distinct characteristics.If we had included robustness as a metric in the definition stage in the same way as we did with f 1 and f 2 , the computational cost would have increased drastically (around 50 times higher).These examples demonstrate the pertinence of the criterion assumed and, consequently, the great convenience of taking into account the potentially useful solutions, that is to say, the nearly optimal solutions that are not dominated in their neighborhood.This approach is specially advisable when the MOP presents multimodal solutions, either because they are intrinsic to the problem or because they have arisen artificially as a result of aggregating several objectives into one single-objective function.Likewise, having the nearly optimal solutions available represents an advantage when there exist objectives which cannot be incorporated within the optimization process due to its high computational cost.When this happens, these objectives are excluded from the optimization stage and evaluated afterwards only for the optimal and nearly optimal solutions resulting from the search.

i
and f min i are the maximum and minimum values of f x .

Figure 6 :
Figure 6: Graph of d t from t = 0 to t = t max − 1.The parameters d ini and d fin are the initial and final values of d t .

Figure 8 :
Figure 8: Boxplots for benchmark 1, comparing nevMOGA, random search, and exhaustive search.(a) Δ p f H , f P * Q,n measures the approximation of the outcome set toward the Pareto front.(b) Δ p H, P * Q,n measures the approximation of the outcome set toward the Pareto set.The most significant values are shown in Table3.

Figure 7 :
Figure 7: Target set P * Q,n , formed by P * Q ∪ P * 1Q .It consists of nine subsets, one of them contains the optimal solutions, and the rest the nearly optimal ones.
nevMOGA outcome set compared with the target set P * Random search outcome set compared with the target set P * Exhaustive search outcome set compared with the target set P * Q,n

Figure 9 :
Figure9: A representative solution of each one of the three algorithms is compared with the target set P * Q,n in the parameter space (right) and in the objective space (left).This representative solution has been chosen as the one whose values of Δ p are closest to the average values of that metric in the objective space and in the parameter space.

Figure 10 :Figure 11 :
Figure 10: The target set P * Q,n is formed by P * Q ∪ P * 1Q .The set P * 4Q is dominated by the set P * 3Q ; however, they lie in distinct neighborhoods and, for this reason, both sets belong to the target set.

1 ( 1 ( 1 (
a) nevMOGA outcome set compared with the target set P * b) Random search outcome set compared with the target set P * c) Exhaustive search outcome set compared with the target set P * Q,n

Figure 13 :
Figure 13: The control structure: a multiloop PI control which uses an off-diagonal pairing scheme.

1 Figure 15 :
Figure15: System responses for x 1 (optimal solution) and x 2 (nearly optimal solution).Although f x 2 ≃ f x 1 , both transient responses are different, since x 1 and x 2 are nonneighboring solutions.

Boundary for x 3
Boundary for x4

Figure 16 :
Figure16: Boundaries of the performance degradation for x 3 (optimal solution) and x 4 (nearly optimal solution) when the plant parameters vary in ±20% from their nominal values.It must be highlighted that a traditional MOP approach would not have found the solution x 4 .

Table 2 :
Analysis of a MOP example with four objectives (f 11 x , f 12 x , f 21 x , and f 22 x ) and four possible solutions (x 1 to x 4 ) where the MOP is reduced by aggregating f 11 and f 12 (f 1 = f 11 + f 12 ) and f 21 and f 22 (f 2 = f 21 + f 22 ).
Lines 7 and 12 (Algorithm 2).Here, it is determined whether a new individual is included in f ront t or not.For this, the individual has to fulfill Definition 12.If the new individual is finally included in front t , then it is analyzed whether its inclusion eliminates other individuals from front t (in which case it will be established whether they must be included in subfront t ) and from subfront t Algorithm 2: Inclusion of x in front t .⊳ Definition 12 1: Choose x F ∈ Front (t) ⊳ Random selection 2: Choose x P1 ∈ P (t) ⊳ According to an exponential distribution (μ = Nind P /10) 3: u≔ random(1) 4: if u> P c /mthen⊳P c /m is the probability of crossover and mutation 5: (ii) Lines 8 and 13 (Algorithm 3).Here, it is determined whether a new individual which was not included in front t has to be included in subfront t or not.For this, the individual has to fulfill Definition 14.If the new individual is included in subfront t , then it is analyzed whether its inclusion eliminates other individuals from subfront t (iii) Line 10 (Algorithm 4).The population G t is created by using the following procedure:(1) Two individuals are chosen at random, one from front t (x F ) (where any individual has the same 1: if ∃ x * ∈ Front t x * ⪯ −ϵ x then 2: return 3: else if ∃ x * ∈ Front t ⋃ Sub − Front t Box x * ⪯ n Box x then 4: return 5: else if ∃ x * ∈ Front t Box x * = n Box x & dist x * ≤ dist x then 6: return 7: end if 8: Include x in Sub-Front (t) 9: if ∃ x * ∈ Sub − Front t Box x ⪯ n Box x * then ⊳ Remove non-nearly-optimal solutions 10: Remove x * from Sub-Front (t) 11: else if ∃ x * ∈ Sub − Front t Box x = n Box x * & dist x ≤ dist x * then 12: Remove x * from Sub-Front (t) then 13: end if Algorithm 3: Inclusion of x in subfront t .⊳ Definition 13 1: if ∃ x * ∈ Front t Box x * ⪯Box x then 2: return 3: else if ∃ x * ∈ Front t Box x * = Box x & dist x * ≤ dist x then * ∈ Front t Box x = Box x * & dist x ≤ dist x *

Table 3 :
Numerical values of the statistical results presented in Figure8.

Table 4 :
Numerical values of the statistical results presented in Figure11.