Preference-Based Evolutionary Multiobjective Optimization Through the Use of Reservation and Aspiration Points

Preference-based Evolutionary Multiobjective Optimization (EMO) algorithms approximate the region of interest (ROI) of the Pareto optimal front defined by the preferences of a decision maker (DM). Here, we propose a preference-based EMO algorithm, in which the preferences are given by means of aspiration and reservation points. The aspiration point is formed by objective values which the DM wants to achieve, while the reservation point is constituted by values for the objectives not to be worsened. Internally, the first generations are performed in order to generate an initial approximation set according to the reservation point. Next, in the remaining generations, the algorithm adapts the search for new non-dominated solutions depending on the dominance relation between the solutions obtained so far and both the reservation and aspiration points. This allows knowing if the given points are achievable or not; this type of information cannot be known before the solution process starts. On this basis, the algorithm proceeds according to three different scenarios with the aim of re-orienting the search directions towards the ROI formed by the Pareto optimal solutions with objective values within the given aspiration and reservation values. Computational results show the potential of our proposal in 2, 3 and 5-objective test problems, in comparison to other state-of-the-art algorithms.


I. INTRODUCTION
Multiobjective optimization problems (MOPs) arising in many real-world applications contain several conflicting objectives f i : S → R, with i = 1, . . . , k (k ≥ 2), that need to be optimized simultaneously over a feasible set in the decision space, S ⊂ R n , formed by decision vectors x = (x 1 , . . . , x n ) T . Generally, they can be formulated as follows: The associate editor coordinating the review of this manuscript and approving it for publication was Jagdish Chand Bansal.
All objective vectors f(x) = (f 1 (x), . . . , f k (x)) T constitute the objective set in the objective space Z = f(S) ⊂ R k . Usually, it is not possible to find a single solution where all the objectives can achieve their individual optima, given that they are in conflict. Instead, the so-called Pareto optimal solutions have to be found, in which none of the objectives can be improved without degrading, at least, one of the others.
A solution x ∈ S is said to be Pareto optimal of problem (1) if and only if there is no otherx ∈ S such that f i (x) ≤ f i (x) for all i = 1, . . . , k, and f j (x) < f j (x) for at least one index j.
The corresponding objective vector f(x) is referred to as a Pareto optimal objective vector. The set of all Pareto optimal solutions is called the Pareto optimal set (PS), denoted by E, VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ while all the Pareto optimal objective vectors form the Pareto optimal front (PF), denoted by f(E). In addition, for two vectors z,z ∈ R k , we say that z dominatesz if and only if z i ≤z i for all i = 1, . . . , k, and z j <z j for, at least, one index j. If z andz do not dominate each other, it is said that they are non-dominated. By non-dominated set we refer to a set of solutions in the decision space whose objective vectors are not dominated by the other solutions in the set. The bounds of the objective functions in the objective space are defined by their worst and their best possible values, which constitute the nadir point z nad = (z nad 1 , . . . , z nad k ) T and the ideal point z = (z 1 , . . . , z k ) T , respectively. To be more precise, the upper objective function bounds are calculated in the nadir point as z nad i = max x∈E f i (x) (i = 1, . . . , k), while the lower bounds can be obtained with the ideal point as z i = min x∈E f i (x) = min x∈S f i (x) (i = 1, . . . , k).
Pareto optimal solutions are incomparable in a mathematical sense and, to find the final solution of the MOP, a decision maker (DM) must be involved in the solution process in order to choose the Pareto optimal solution which best suits her/his preferences (the most preferred solution).
Evolutionary Multiobjective Optimization (EMO) algorithms [2], [3] are developed for solving MOPs. These algorithms generate a non-dominated set to approximate the entire PF, meaning that ideally the solutions generated should be uniformly distributed (diversity) and close enough to the PF (convergence). However, achieving convergence and diversity simultaneously is not always easy, especially when handling many-objective optimization problems (k > 3), since e.g. the percentage of non-dominated solutions increases with the number of objectives, or in problems with complicated PFs (such as e.g. non-convex, degenerated, or discontinuous). To cope with these issues, decompositionbased EMO algorithms transform the original MOP into a set of single-objective optimization sub-problems formulated using an aggregation function with one weight vector, and usually select the best individuals at each generation according to these sub-problems. Some examples of these algorithms are GWASF-GA [22], MOEA/D [14], [29], and NSGA-III [5].
The approximation found by an EMO algorithm can be useful to gain an insight into the objective functions (their ranges, their trade-offs, their degree of conflict, etc.). However, the task of identifying the most preferred Pareto optimal solution may be too demanding for a DM if (s)he has to analyze and compare all the solutions generated by an EMO algorithm, whose approximation set usually consists of a large number thereof. To facilitate this decisionmaking task, preferences can be included in the solution process, resulting in the so-called preference-based EMO algorithms [1]. Their main purpose consists of only approximating the subset of solutions which correspond to the preferences of the DM (the region of interest (ROI)) [6], [8], [25]. Within this group, note can be made of the WASF-GA algorithm [21], at which the preferences are elicited by means of a reference point q = (q 1 , . . . , q k ) T constituted by desirable values q i for the objective functions f i (i = 1, . . . , k). In practice, a reference point q can be achievable if q is equal or dominated by a Pareto optimal objective vector (that is, the desirable values can be simultaneously achieved or improved by a Pareto optimal solution), or unachievable in any other case.
Once a reference point is given as preference information, the ROI with the most interesting Pareto optimal solutions according to these preferences can be defined as suggested in [21]. To be more precise, we assume that the ROI associated with a reference point q is defined as follows: when q is achievable, the ROI is the subset of Pareto optimal objective vectors dominating q, and if q is unachievable, the Pareto optimal objective vectors dominated by q constitute the ROI.
Besides using preferences, WASF-GA is an algorithm based on decomposition. It considers a set of weight vectors to internally formulate single-objective problems using the following achievement scalarazing function [27]: where q is a reference point provided by the DM, µ = (µ 1 , . . . , µ k ) T with µ i > 0 for i = 1, . . . , k denotes the vector of weights, and ρ > 0 is a parameter to assure that the solution that minimizes (2) over the feasible set S is always a Pareto optimal solution to the original MOP (1). See [18] for further details.
Recently, IRA-EMO [23] has been proposed as an interactive EMO method based on WASF-GA. In this method, the preferences are provided in the form of preferred ranges for the objective functions. Given the multiobjective minimizing problem (1), for every i = 1, . . . , k, the DM specifies an aspiration level q a i (the lower desirable bound) and a reservation level q r i (the upper desirable bound) for the values of the objective f i . It is implicitly assumed that z i ≤ q a i < q r i ≤ z nad i (i = 1, . . . , k). With this information, two reference points can be built: an aspiration point q a = (q a 1 , . . . , q a k ) T formed by values regarded as desirable for the objective functions (i.e. values to be achieved, if possible), and a reservation point q r = (q r 1 , . . . , q r k ) T comprised of acceptable objective levels beyond which their values should be rejected. This preference scheme allows the DM to freely investigate different parts of the PF by modifying the desirable bounds for the objective functions (i.e. the aspiration and reservation points).
In IRA-EMO, Modified WASF-GA is the algorithm internally defined to generate a set of non-dominated solutions according to both points. However, this algorithm is designed to approximate the union set of the ROIs associated with q a and q r , using a small number of solutions. Nonetheless, when providing preferred objective function ranges, the DM may only be interested in the Pareto optimal solutions whose objective values are between the corresponding aspiration and the reservation levels (i.e. in the solutions lying in the intersection set of the ROIs defined by the two points).
In this paper, we focus on the use of preferences in EMO to reduce the search space to the ROI that best suits the desires of the DM. We assume that the preferences are expressed in the form of aspiration and reservation points, q a and q r . As mentioned above, when giving this information the DM wants to analyze solutions with objective function values within the desired bounds, i.e. among the aspiration and the reservation levels. Thus, the ROI to be approximated, denoted by R, includes the Pareto optimal solutions which belong, at the same time, to the ROIs associated with both points. If we denote by R a and R r the ROIs defined by the aspiration and the reservation points, respectively, the subset R with the most interesting Pareto optimal solutions for the DM can be obtained as R = R r ∩ R a . Figure 1 shows different situations for a bi-objective minimization problem, where R is the region of the PF in bold.
To converge to this ROI, we propose a decompositionbased algorithm enabling the simultaneous use of both the aspiration and reservation points, which improves the convergence procedure of the algorithms previously cited. The main purpose is to directly focus the search for new solutions onto the ROI jointly defined by the two points, avoiding producing solutions outside this region, and to do so not only with a few solutions but also with a sufficiently representative approximation set.
Our proposal is called Evolutionary algorithm based on Reservation-Aspiration Levels (ERAL). Using points q a and q r given by the DM and an initial set of weight vectors, ERAL produces an initial approximation of the ROI defined by q r in the first generations. Using the solutions obtained, the achievability of the reservation and aspiration points is investigated based on the relation of dominance between these solutions and the two points. Then, the algorithm procedure is adjusted accordingly to the scenario detected in order to converge to the subset R jointly defined by the aspiration and reservation points.
Moreover, it is well-known that the weight vectors' distribution directly affects the performance of any decompositionbased EMO algorithm, since the diversity of the population is somehow controlled by the set of weight vectors used. In problems with complicated PFs, using a well-distributed set of weight vectors does not guarantee the generation of a well-diversified approximation set of the PF. In general, the best results are obtained when the distribution of the weight vectors is consistent with the shape of the PF, but the PF's characteristics are usually unknown beforehand. To alleviate the effects of using a pre-fixed set of weight vectors, there are many proposals to update the distribution of the weights as the algorithm converges, such as e.g. [9], [11], [13], [16], [17], [20], [24]. In this line, the ERAL algorithm also incorporates the possibility of dynamically adapting the weight vectors according to the PF's shape in the last generations. Specifically, we have considered two procedures to adjust the weights: the dynamic adjustment of weight vectors proposed in [17] and the approach to update weights during an evolutionary process called AdaW [16].
The new algorithm ERAL is described below in Section II. Next, the computational experiment carried out to show the potential of our proposal is shown in Section III. Finally, the conclusions are drawn in Section IV.

II. EVOLUTIONARY ALGORITHM BASED ON RESERVATION AND ASPIRATION LEVELS A. MOTIVATION
As mentioned above, according to the aspiration and reservation levels given by a DM, Modified WASF-GA [23] generates an approximation set of the region R r ∪ R a , instead of directly approximating the intersection set R = R r ∩ R a . This may not entail a significant waste of computational effort if the purpose is to generate a small number of solutions, as in IRA-EMO [23], although this also depends on the complexity of the problem considered. Nevertheless, there are more efficient ways to approximate R and avoid producing solutions outside it. This motivated us to propose the Evolutionary algorithm based on Reservation and Aspiration Levels (ERAL), as opposed to Modified WASF-GA, directly focuses the search for new non-dominated solutions on the subset R. VOLUME 9, 2021 Our proposal is a decomposition-based algorithm that internally considers the ASF given in (2) to converge to the desired region R. Due to the practical meaning of minimizing this ASF over the feasible set, 1 it is important that the weight vectors considered are defined so that they produce projection directions that allow finding a well-diversified approximation of the PF [22], or of the subset of solutions to be approximated [21]. Therefore, in ERAL, a set of weight vectors are initially generated so that they define uniformly distributed projection directions.
The first generations of ERAL are performed using the weight vectors initially generated and the reservation point. Once an approximation set has been produced, the algorithm's working procedure is reconsidered depending on the dominance relation of the non-dominated solutions obtained so far with respect to both the reservation and aspiration points. This dominance information allows us to know the achievability of the aspiration and reservation points: whether both points are of the same type (achievable or unachievable) or if the aspiration point is unachievable and the reservation point is achievable. Let us remark that this information is unknown before the algorithm starts; indeed, even the DM is not able to foresee the achievability of both points when providing them. Then, the algorithm adopts different strategies to perform the remaining generations according to the knowledge gained in this respect, with the aim of guiding the search as best as possible towards the subset R to be approximated, which changes depending on the achievability of the two points.
In addition, the ERAL algorithm also allows the weight vectors to be dynamically adapted according to the PF's characteristics. As mentioned above, in decomposition-based EMO algorithms, the diversity of the final population is somehow controlled by the weight vectors, and a pre-fixed set of weight vectors may not lead to a uniformly distributed set of solutions in the region to be approximated [11], [13]. Several papers propose internal mechanisms to adapt the weight vectors along the convergence process [11], [16], [17], [20], [24]. Following a similar idea, the vectors of weights initially generated in ERAL can be internally updated. For this purpose, the weight vector adjustment procedures suggested in [16] and [17] can be executed in the last generations of the algorithm to better adjust the search directions towards R, according to the specific features of the PF associated with the problem being solved.

B. THE ALGORITHM
Next, the ERAL algorithm is described in detail. Let q r and q a denote the reservation and the aspiration points given by the DM, respectively. As said before, it may happen that q r is achievable and q a is unachievable, but it is also possible that q r and q a are both achievable, or that both of them 1 Minimizing the ASF given in (2) over S means, in practice, to project q onto the PF in the projection direction defined by the inverse components of the weight vector used (i.e. by 1 are unachievable, but this information is unknown at the beginning. As already said, a set of weight vectors is used to internally decompose the original MOP into various sub-problems. Let W = {µ 1 , . . . , µ N µ } be the set of weight vectors in (0, 1) k used, where N µ > 0 is the number of vectors considered. These weight vectors are generated so that their inverse components (which define their search or projection directions in the objective space) are evenly distributed in the weight vector space. Besides, let us denote by N the population size and by G T the total number of generations.
The first generations of ERAL are performed as in the original WASF-GA [21], using the weight vectors in W and q r as reference point. If G p denotes the number of generations run using q r as reference point, we can define it as G p = p · G T , where 0 < p < 1 is a real number. In practice, the purpose of executing these G p generations is to produce an initial approximation of the ROI only associated with q r . Thus, G p cannot be too small in order not to stop the algorithm too prematurely, nor too large so that there are still generations to improve the convergence towards the subset R (delimited by both q r and q a ) before the process finalizes.
Once these G p generations have been completed, we check the dominance relation of the non-dominated solutions obtained with respect to both the aspiration and reservation points. Let P G p denote the population of non-dominated solutions at generation G p . Internally, we build the following subsets: • Subset of solutions whose objective vectors dominate the reservation point: • Subset of solutions whose objective vectors dominate the aspiration point: • Subset of solutions whose objective vectors dominate the reservation point and, at the same time, are dominated by the aspiration point: Depending on whether there are solutions or not in D a , D r and D r,a , three different scenarios can be distinguished that tell us if q a and q r are achievable or not. Accordingly, the algorithm adapts its working procedure to perform the remaining number of generations (i.e. the last G 1−p = (1 − p) · G T = G T − G p generations as described next: • Scenario 1: If D r = ∅, that is, none of the nondominated solutions generated so far dominates the reservation point, then we can assume that both q r and q a are unachievable (D r = ∅ implies that D a = ∅). In this case, the ROI R to be approximated is defined just by the reservation point q r and the last G 1−p generations of ERAL are performed using in WASF-GA the same set of weight vectors W and q r a as reference point.
• Scenario 2: If D a = ∅, then, there exists at least one solution in P G p which dominates the aspiration point, meaning that both q r and q a are achievable (D a = ∅ implies that D r = ∅). In this situation, the subset R to be approximated is defined just by the aspiration point q a and the last G 1−p generations of ERAL are executed using in WASF-GA the weight vectors in W and q a as a reference point.
• Scenario 3: Whether D r,a = ∅ (and this implies that D r = ∅) and D a = ∅, then, among the solutions generated at generation G p , we find at least one solution dominating the reservation point that is also dominated by aspiration point q a ; at the same time, there are no solutions dominating q a . Thus, this implies that q r is achievable and q a is unachievable. In this case, the way to proceed with the remaining G 1−p generations is the following one. Let us consider that q min,i is the (best) minimum value achieved by the objective function f i among the solutions in the subset D r,a (i = 1, . . . , k). That is, q min,i = min x∈D r,a f i (x). Then, a new aspiration point is defined asq a = (q min,1 , q min,2 , . . . , q min,k ) and the last G 1−p generations of ERAL are performed as in WASF-GA using the weight vectors in W andq a as a reference point. Figure 2 shows several situations corresponding to this scenario. As can be observed, point q a is calculated so that it adjusts to the delimited region R (in bold). To some extent, the non-dominated solutions produced at each generation of ERAL are the best individuals produced so far for minimizing the ASF given in (2), using the weight vectors in W , and with respect to either q r , q a orq a , depending on the generation being executed, as previously described in the three scenarios.
Last but not least, the weight vectors in the set W can be updated, if desired, so that the convergence towards R is performed taking into account the PF's characteristics. In the last G 1−p generations (corresponding to the three scenarios), the procedures suggested in [16], [17] can be applied to re-calculate the weight vectors during the optimization process based on the non-dominated solutions already obtained, which are likely to reflect the shape of the PF since they are gradually converging towards it. The idea underlying these weight vectors' adaptations consist of replacing the weight vectors generating solutions in overcrowded areas by new ones pushing the search towards the parts with a lack of solutions. The main differences between both procedures lye in the way to identify the vectors that will be removed and to calculate the new vectors that will replace the removed ones. For more details, see [16], [17]. In ERAL, these approaches are applied to adapt all the weight vectors in W only once, when 10% of the remaining generations have been performed (i.e. at generation G p + 0.1 · G 1−p ). Note that any other procedure allowing us to update the weight vectors can also be used in ERAL.
To summarize, Algorithm 1 contains ERAL's main steps. We denote by h the current generation, P h the population of individuals at generation h, Q h the offspring population at generation h, Z h the combined population of parents and offspring at generation h, F h n the n-th front at generation h and P F the final set of non-dominated solutions generated. Algorithms 2, 3, 4 individually describe the steps followed in the three different scenarios.
Finally, we would like to say that, in case the DM is only interested in providing one of the two points, ERAL can also be applied just with the aspiration point or the reservation one. This option may be suitable if one needs to relax the cognitive burden required by using our algorithm to solve a specific MOP. When the DM only specifies an aspiration point, an approximation to the nadir point can be considered as the reservation point, while if only a reservation point is provided, an approximation to the ideal point can be used as the aspiration point.

III. COMPUTATIONAL EXPERIMENTS A. EXPERIMENTAL DESIGN
To the best of our knowledge, none of the existing reference point-based EMO algorithms (see e.g. [6], [8], [25]) allows working with aspiration and reservation points as our proposal does (in the sense that these two points define desirable lower and upper bounds for the objective function values). Note that algorithms enabling the use of several reference points, such as e.g. [6], [25], can be executed using the

Algorithm 1 Main Steps of ERAL
Require: The aspiration and reservation points q a and q r , a set of weight vectors W = {µ 1 , . . . , µ N µ }, the population size N , the number of generations G T and the percentage p of generations initially executed with q r . Ensure: A approximation set P F of the region R.
1: Set G p = p · G T , G 1−p = (1 − p) · G T and h = 0. 2: Randomly generate N solutions to form the initial population P 0 . 3: while h ≤ G p do 4: Set Z h = ∅ and n = 0.

5:
Apply the selection, crossover and mutation operators to P h to create a population Q h with N new solutions. Set Z h = P h ∪ Q h .

9:
for j = 1, . . . , N µ do 10: Update h = h + 1. 22: end while 23: Using the population P G p , obtain the subsets D r , D a and D r,a as in (3), (4) and (5), respectively. 24: if D r = ∅ then 25: Execute Algorithm 2 using q r , the set of weight vectors W , the remaining number of generations G 1−p and P G p as the population for the next generation. Set P F as the outcome of this algorithm. 26: end if 27: if D a = ∅ then 28: Execute Algorithm 3 using q a , the set of weight vectors W , the remaining number of generations G 1−p and P G p as the population for the next generation. Set P F as the outcome of this algorithm. 29: end if 30: if D r,a = ∅ and D a = ∅ then 31: Execute Algorithm 4 usingq a obtained from the set of non-dominated solutions D r,a , the set of weight vectors W , the remaining number of generations G 1−p and P G p as the population for the next generation. Set P F as the outcome of this algorithm. 32: end if aspiration and reservation points as two different reference points, but would approximate the ROIs determined by the two points independently. That is, their working procedures are not designed to converge to the region R = R r ∩R a , as our algorithm does. Therefore, the existing reference point-based For every x ∈ Z h , calculate s(q r , f(x),μ j ) with j = 1, . . . , N µ . 6: while Z h = ∅ do 7: Set n = n + 1 and initialize F h n = ∅.
EMO algorithms are not valid to perform a fair comparison against ERAL. In addition, we have not considered Modified WASF-GA in the computational tests because this approach only generates a few non-dominated solutions to approximate the region R r ∪ R a , and not to directly converge to R, as already explained.
Owing to this, we have considered the EMO algorithms NSGA-III [5] and MOEA/D-DE [12] but, given that these algorithms directly approximate the whole PF, we have introduced the preference information in q a and q r as constraints bound to the objective functions in the problems' formulation in order to perform a fair comparison. That is, when running these algorithms (hereinafter denoted by cNSGA-III and cMOEA/D-DE, respectively), the following constraints are introduced into the formulation of each test problem considered: q a i ≤ f i (x) ≤ q r i , for each i = 1, . . . , k. This assures that they generate non-dominated solutions with objective function values within the desired values; that is, cNSGA-III and cMOEA/D-DE will approximate the region R defined by these two points.

Algorithm 3 Scenario 2
Require: An aspiration point q a , a set of weight vectors W = {µ 1 , . . . , µ N µ }, the population size N , the remaining number of generations G 1−p , and an initial population P with N solutions. Ensure: A approximation set P F of the region R.

4:
Apply the selection, crossover and mutation operators to P h and create a population Q h with N new individuals. Set Z h = P h ∪ Q h .
Moreover, in the comparison, we also execute two variants of ERAL at which dynamic adaptations of the weight vectors are carried out while the algorithm is converging. These two variants are referred to as Adapt-ERAL and AdaW-ERAL, depending on whether the weight vectors' adaptation is done following [17] or [16], respectively.
The test problems involved in our study are from the families DTLZ [7], UF [30], WFG [10], and ZDT [31]. We have considered problems with two, three and five objective functions. These test problems have various properties, such as having a linear, multi-modal, concave, discontinuous, or degenerate PF, having a multi-modal, biased, or deceptive search space, and/or having strong-linkage decision variables [15]. Therefore, approximating their PFs represents a challenge for any of the algorithms included in the comparison. In Table 1, we give a detailed description of the main features of these problems, and we also indicate the number of objectives (k) considered for each problem. In addition, we use k + 4, k + 9, and k + 19 decision variables for the DTLZ1, DTLZ2-DTLZ4, and DTLZ7

Algorithm 4 Scenario 3
Require: A set of non-dominated solutions D r,a , a set of weight vectors W = {µ 1 , . . . , µ N µ }, the population size N , the remaining number of generations G 1−p , and an initial population P with N solutions. Ensure: A approximation set P F of the region R.
problems, respectively. For the WFG problems, we set the position-and distance-related parameters to k − 1 and 10, respectively. Our experiment has been implemented in Java using the jMetal framework [19]. For each problem, we have randomly generated six pairs of reservation and aspiration points that have been used to test each algorithm. Note that using such number of reservation and aspiration points ensures that the three scenarios described in Section II are considered in our computational study. Let us remark that, before running the algorithms, it is not possible to know if the generated points are achievable or unachievable. Besides, to have statistically significant results, we have executed 10 independent runs for each algorithm, each test problem and each pair of aspiration and reservation points. Thus, a total of 60 independent runs have been executed for each problem and each algorithm. For the two-, three-, and five-objective problems, respectively, the population sizes used are N = 50, 91 and 210, and the maximum number of generations considered (stopping criterion) is G T = 400, 600 and 800.
In ERAL, Adapt-ERAL, AdaW-ERAL and cNSGA-III, we use the SBX crossover operator [4] (distribution index η c = 30 and probability P c = 0.9), and the polynomial mutation operator [4] (distribution index η m = 20 and probability P m = 1/n, where n is the number of variables). In the three ERAL variants, the value p has been set to 0.6, which means that G p = 240, 360 and 480 for the two-, three-, and five-objective problems, respectively. In cMOEA/D-DE, the neighborhood size used is 0.1 · N and the probability of choosing the mate sub-problem from the neighborhood is 0.9. The crossover ratio and the scale factor are both set to 0.5 for the DE operator. For a fair comparison, all algorithms initially use the same set of weight vectors.
The algorithms' performances have been evaluated with the hypervolume (HV) metric, which is obtained using the WFG calculation method suggested in [26] for all test problems. For the HV calculation, the reference point needed has been set depending on the achievability of the aspiration and reservation points considered in each case, so that it is suited as best possible to the ROI to be approximated. Specifically, the nadir point, the aspiration point and the reservation point have been used as the reference point required, respectively, in scenarios 1, 2 and 3 described in Section II.
Note that the HV values obtained by each algorithm have been normalized to have comparable values. For each problem and for each pair of aspiration and reservation points (each region of interest), the HV value of each algorithm has been divided by the sum of the HV values of the five algorithms considered. Then, we have calculated the mean of the normalized HV values obtained in all the runs for each problem. Hereinafter, for simplicity, when we talk about the HV mean value, we refer to the mean of the HV values normalized in this way.
In addition, we apply a Wilcoxon rank-sum test [28] to check if the HV mean values achieved by ERAL, Adapt-ERAL and AdaW-ERAL are significantly different from the ones of the other algorithms. For each problem, the null hypothesis is that the distribution of their HV mean values in the 60 runs differ by a value α, assuming that the difference is significant if the p-value obtained is lower than α = 0.05. We use the wilcox.test function from the R software. 2 Tables 2, 3 and 4 show the HV mean values obtained in the 60 independent runs, for the two-, three-and five-objective problems, respectively. For each problem, we emphasize the algorithm with the best HV value in dark gray color, and the one with the second best value in light gray color. Besides, in Table 5, we indicate the numbers of problems for which the HV mean values of ERAL, Adapt-ERAL or AdaW-ERAL are significantly better than ( ), equal to ( ), or worse than ( ) the HV mean value of the other algorithms, according to the Wilcoxon test with a significance level α = 0.05.

B. RESULTS
In relation to the two-objective problems (see Table 2), our results indicate that the three variants of ERAL perform better than cMOEA/D-DE and cNSGA-III in all problems. Observe that ERAL reaches the best HV mean value in 15 out of the 26 problems (hereafter indicated as 15/26), Adapt-ERAL is the best in 1/26 problems, AdaW-ERAL achieves the best results in 10/26 cases, while cMOEA/D-DE and cNSGA-III do not win in any problem. Moreover, the Wilcoxon test (see Table 5) reveals that, in the problems with two objectives, ERAL shows a statistically better performance in 24/26 cases compared to cNSGA-III, and it produces significantly better results in 16/26 cases with respect to cMOEA/D-DE. Similar results are found by the two ERAL variants. Indeed, Adapt-ERAL is significantly better in 25/26 cases compared to both cNSGA-III and cMOEA/D-DE, and AdaW-ERAL obtains significantly better results in 24/26 and 25/26 problems when compared to cNSGA-III and cMOEA/D-DE, respectively. The ROI approximation of the ZDT4 problem is shown in Figure 3. This representation corresponds to the ROI obtained for each algorithm in ten independent runs for a fixed level of preferences (an aspiration point and reservation point).
Concerning the problems with three objectives (see Table 3), in general, ERAL, Adapt-ERAL and AdaW-ERAL also outperform the other two algorithms. Specifically, each one wins in 6/16, 1/16 and 4/16 problems, respectively.   In turn, cNSGA-III achieves the best HV mean in 3/16 problems, and cMOEA/D-DE wins in 2/16 cases. Regarding the significance of the results for the 3-objective problems (see Table 5), both ERAL and Adapt-ERAL obtain significantly better results in 9/16 problems compared to cNSGA-III, and in 8/16 problems when compared with  cMOEA/D-DE, while AdaW-ERAL is significantly better than both cNSGA-III and cMOEA/D-DE in 9/16 problems.
Finally, let us analyze the results for the five-objective problems (see Table 4). Overall, ERAL and its variants achieve the best results once more (they win in 12/14 problems). We can see that ERAL is the best one in 4/14 problems, Adapt-ERAL in 5/14 problems and AdaW-ERAL in 3/14 problems. However, cNSGA-III obtains the best results in just 2/14 cases, while cMOEA/D-DE does not win in any problem. Moreover, in relation to the Wilcoxon test (see Table 5), the HV mean values of ERAL, Adapt-ERAL and AdaW-ERAL are statistically better than those of cNSGA-III in 11/14 problems, and better than the cMOEA/D-DE's in 14/14 problems.
With this computational study, our proposal's performance has been demonstrated in a wide range of benchmark problems. Overall, ERAL and its two variants have outperformed the other two algorithms in most of the problems.
Moreover, although we have not analyzed in detail which algorithm has achieved the second best HV mean value, it can be seen at a glance that, whenever they do not obtain the best results, ERAL, Adapt-ERAL and AdaW-ERAL are ranked as the second best ones in most cases of the two-, threeand five-objective problems. This indicates that our proposal seems to achieve very promising results, especially as the number of objective functions increases.
Nevertheless, we cannot clearly state that any of the variants of the proposed algorithm is more adequate than the others for the problems considered, given that the results they reached were very similar (at least for the HV metric). In the problems with two objectives, we observe that the fact of internally adapting the weight vectors while ERAL is converging may imply a significant improvement in the approximation set obtained, since ERAL stood out from Adapt-ERAL and AdaW-ERAL in these cases. However, in the comparison using the problems with three and five objectives, the performances of the three algorithms were very similar. Anyway, since the purpose of performing an adaptation of the weight vectors is to orient the search for new non-dominated solutions according the PF's features, Adapt-ERAL and AdaW-ERAL may be more suitable when solving problems with complicated PFs.

IV. CONCLUSION
In this paper, a new preference-based EMO algorithm has been proposed that manages preference information in the form of aspiration and reservation points. The new algorithm, named ERAL, is designed to approximate the ROI of the PF formed by the Pareto optimal solutions whose objective function values are within the desired aspiration and reservation levels. To achieve this, a number of generations are first executed using the reservation point and an initial set of weight vectors. Next, according to the relation of dominance between the solutions generated so far and the two points given by the DM, the algorithm analyzes the aspiration and reservation points to see what type they are (achievable or unachievable). Finally, the other generations are performed following different procedures depending on their achievability, which determines the ROI to be approximated.
Our proposal also incorporates the possibility of adapting the weight vectors used internally, with the aim of adjusting the convergence towards the desired ROI as much as possible, taking into account the PF's characteristics. Two procedures have been incorporated in ERAL for the adaptation of the weight vectors based on [17] and [16].
For the two-, three-, and five-objective benchmark problems considered in the computational experiments, we concluded that ERAL and its two variants, Adapt-ERAL and AdaW-ERAL (which perform the adaptation of the weights following [16], [17], respectively), have generated very promising results. They have shown a better performance in comparison to other well-known EMO algorithms, according to the HV metric.
As future research, we plan to study how the algorithm performs in higher dimension benchmark problems. In addition, we want to investigate the potential of ERAL to solve real-world problems using preferences given by real DMs, especially in areas such as engineering or economics.