Integration of Preferences in Decomposition Multiobjective Optimization

Rather than a whole Pareto-optimal front, which demands too many points (especially in a high-dimensional space), the decision maker (DM) may only be interested in a partial region, called the region of interest (ROI). In this case, solutions outside this region can be noisy to the decision-making procedure. Even worse, there is no guarantee that we can find the preferred solutions when tackling problems with complicated properties or many objectives. In this paper, we develop a systematic way to incorporate the DM’s preference information into the decomposition-based evolutionary multiobjective optimization methods. Generally speaking, our basic idea is a nonuniform mapping scheme by which the originally evenly distributed reference points on a canonical simplex can be mapped to new positions close to the aspiration-level vector supplied by the DM. By this means, we are able to steer the search process toward the ROI either directly or interactively and also handle many objectives. Meanwhile, solutions lying on the boundary can be approximated as well given the DM’s requirements. Furthermore, the extent of the ROI is intuitively understandable and controllable in a closed form. Extensive experiments on a variety of benchmark problems with 2 to 10 objectives, fully demonstrate the effectiveness of our proposed method for approximating the preferred solutions in the ROI.


Integration of Preferences in Decomposition
Multiobjective Optimization Ke Li , Member, IEEE, Renzhi Chen, Geyong Min, Member, IEEE, and Xin Yao , Fellow, IEEE Abstract-Rather than a whole Pareto-optimal front, which demands too many points (especially in a high-dimensional space), the decision maker (DM) may only be interested in a partial region, called the region of interest (ROI).In this case, solutions outside this region can be noisy to the decision-making procedure.Even worse, there is no guarantee that we can find the preferred solutions when tackling problems with complicated properties or many objectives.In this paper, we develop a systematic way to incorporate the DM's preference information into the decomposition-based evolutionary multiobjective optimization methods.Generally speaking, our basic idea is a nonuniform mapping scheme by which the originally evenly distributed reference points on a canonical simplex can be mapped to new positions close to the aspiration-level vector supplied by the DM.By this means, we are able to steer the search process toward the ROI either directly or interactively and also handle many objectives.Meanwhile, solutions lying on the boundary can be approximated as well given the DM's requirements.Furthermore, the extent of the ROI is intuitively understandable and controllable in a closed form.Extensive experiments on a variety of benchmark problems with 2 to 10 objectives, fully demonstrate the effectiveness of our proposed method for approximating the preferred solutions in the ROI.Index Terms-Decomposition-based method, evolutionary multiobjective optimization (EMO), reference points, user-preference incorporation.

I. INTRODUCTION
M ANY REAL-LIFE applications usually consider opti- mizing multiple conflicting objectives simultaneously.To handle such problems, called multiobjective optimization problems (MOPs), decision makers (DMs) often look for a set of Pareto-optimal solutions, none of which can be considered better than another when all objectives are of importance.Evolutionary multiobjective optimization (EMO) algorithms, which work with a population of solutions and can approximate a set of tradeoff alternatives simultaneously, have been widely accepted as a major tool for solving MOPs.Over the past two decades and beyond, many efforts have been devoted to developing EMO algorithms (e.g., elitist nondominated sorting genetic algorithm (NSGA-II) [1], indicator-based EA [2], multiobjective EA based on decomposition (MOEA/D) [3], and their variants [4]- [10]) to find a set of efficient solutions that well approximate the whole Pareto-optimal front (PF) in terms of convergence and diversity.
The ultimate goal of multiobjective optimization is to help the DM find solutions that meet, at most, his/her preferences.Supplying a DM with a large amount of tradeoff points, which approximate the entire PF, not only increases his/her workload but also provides many irrelevant or even noisy information to the decision-making procedure.Moreover, due to the curse of dimensionality, approximating a high-dimensional PF as a whole not only becomes computationally inefficient (or even infeasible) but also causes a severe cognitive obstacle for the DM to comprehend the high-dimensional data.To facilitate the decision-making procedure, it is more practical to incorporate the DM's preference information into the search process.This allows the computational efforts to concentrate on the region of interest (ROI) and, thus, has better approximation therein.In general, the preference information can be incorporated a priori, posteriori, or interactively.Note that the traditional EMO goes along the posteriori way of which the disadvantages have been discussed before.If the preference information is elicited a priori, it is directly used to guide the solutions toward the ROI.However, it is nontrivial to faithfully model the preference information before solving the MOP at hand.Eliciting the preference information in an interactive manner has been studied in the multicriterion decision-making (MCDM) field for over half a century.It enables the DM to progressively learn and understand the characteristics of the MOP at hand and adjust his/her elicited preference information.Consequently, solutions are effectively driven toward the ROI.However, since the optimization process is full of uncertainty and the DM is almost unavoidable This work is licensed under a Creative Commons Attribution 3.0 License.For more information, see http://creativecommons.org/licenses/by/3.0/ to show inconsistencies in decision-making [11], it is difficult to model the DM's behavior in an appropriate manner.
Integrating and blending the EMO and MCDM together to tailor the DM's preference information has been studied since the 1990s.Although the existing works aim at steering the search process toward the ROI, the definition of the ROI is still vague.First of all, the ROI can be any part of the PF near the DM-supplied aspiration-level vector or even subjectively determined by the DM.Second, the ROI is expected to be a partial region of the PF whereas no quantitative definition has been given to the size of this region.Although some studies (e.g., [12]- [14]) claimed to control the spread of the preferred solutions accommodating the DM's expectation of the extent of the ROI, that is, the ROI's size, the corresponding parameter setting is ad-hoc [15].In addition to the ROI, the boundary of the PF is also important for the DM to understand the underlying problem and to facilitate the further decisionmaking procedure.In particular, the boundary provides the DM general information about the PF's geometrical characteristics; and more important, it provides the information of the ideal and nadir points which facilitate the normalization of the disparately scaled objective functions.Unfortunately, how to keep solutions located in the ROI and the boundary simultaneously has rarely been studied [16].
During recent years, especially after the developments of MOEA/D and NSGA-III [17], the decomposition-based EMO methods have become increasingly popular for the posteriori multiobjective optimization.Generally speaking, by specifying a set of reference points,1 the decomposition-based EMO methods at first decompose the MOP at hand into multiple subproblems, either with scalar objective or simplified multiobjective.Then, a population-based technique is applied to solve these subproblems in a collaborative manner.Under some mild conditions, the optimal solutions of all subproblems constitute a good approximation to the PF.It is not difficult to understand that the distribution of the reference points is essential in a decomposition-based EMO method.It not only implies a priori prediction of the PF's geometrical characteristics but also determines the distribution of Pareto-optimal solutions.Das and Dennis [18] and Tan et al. [19] suggested some structured methods to generate evenly distributed reference points on a canonical simplex.To adapt to the irregular PFs, such as disconnected or mixed shapes and disparately scaled objectives, some adaptive reference point adjustment methods (e.g., [20]) have been developed to adjust the distribution of reference points on the fly.To integrate the DM's preference information into the decomposition-based EMO methods, a natural idea is to make the distribution of the reference points be biased toward the ROI.Although it sounds quite intuitive, in practice, how to obtain the appropriate reference points that accommodate the DM's preference information is far from trivial.Most recently, there have been some initiatives on adjusting the distribution of the reference points according to the DM's preference information (e.g., [21]).However, they are ad-hoc and the position and extent of the reference points around the ROI are not fully controllable.
In this paper, we present a systematic way to incorporate the DM's preference information, either a priori or interactively, into the decomposition-based EMO methods.In particular, the DM's preference information is modeled as an aspiration-level vector, which has been widely used in the EMO literature [15].Compared to three state-of-the-art preference-based EMO algorithms, the effectiveness and competitiveness of the proposed preference incorporation method for assisting three state-of-the-art decomposition-based EMO algorithms to approximate ROIs have been validated through extensive experiments on 56 test problems with 2 to 10 objectives, under both attainable and unattainable aspiration-level vector settings.Our major contributions are outlined as follows.
1) Our basic idea is a nonuniform mapping scheme (NUMS) by which the originally evenly distributed reference points on a canonical simplex can be mapped to new positions close to the DM-specified aspiration-level vector and thereby have biased distribution.
2) The mapping function is nonlinear in nature and is a function of a reference point's position with respect to the pivot point.Accordingly, the distribution of the reference points after the nonuniform mapping is biased toward the pivot point.In particular, this pivot point is representative of the ROI on the simplex and determines the ROI's position.3) Different from the existing preference-based EMO algorithms, where the extent of the approximated ROI is controlled in an ad-hoc manner, this paper provides an intuitively understandable manner to quantify this extent in a closed form.It is the ratio of the biased reference points proportional to the simplex.To a certain extent, this quantity can be used as the ratio of the ROI's size with respect to the PF. 4) Given the DM's requirements, the proposed NUMS is able to not only obtain a set of biased reference points toward the ROI, but also preserve the ones located on the boundary.This latter characteristic enables a decomposition-based EMO method to not only find the preferred solutions but also provide the global information about the PF to the DM.The rest of this paper is organized as follows.Section II overviews some state of the art related to this paper.Section III presents the technical details of our proposed NUMS.Section IV shows the empirical studies on several benchmark problems.Finally, Section V concludes this paper and provides some future directions.

II. RELATED WORKS
In the past two decades, various methods have been developed to incorporate the DM's preference information into the EMO.This section briefly overviews the existing literature according to the ways of eliciting the DM's preference information along with the mechanisms adopted to guide the population toward the ROI.The interested reader is referred to [15] for a recent comprehensive survey.
The first one employs the weight information, that is, relative importance, to model the DM's preference information.For example, Deb [22] developed a modified fitness-sharing mechanism, by using a weighted Euclidean distance, to bias the population distribution.Branke and Deb [23] developed a weighted mapping method to modify the crowding distance calculation of NSGA-II by which the search process can be guided toward the ROI.Note that the weight-based methods become ineffective when facing a large number of objectives, because it is difficult to either specify the weights or verify the quality of the biased approximation.Moreover, it is unintuitive and challenging for the DM to steer the search process toward the ROI via the weighting scheme.
The second sort elicits the preference information by inviting the DM to make pairwise comparisons among a sample of solutions from a population.As the pioneer, Phelps and Köksalan [24] proposed using a value function model to represent the DM's preference information.Note that the precise form of the value function model is unknown a priori.It is progressively learned through a periodic interaction with the DM during the optimization process.In particular, the DM is asked to express his/her preference information about some selected alternatives, for example, their rankings, at each interaction session.Inspired by [24], many variants have been developed by using various value function models, for example, quasi-concave preference function [25], polynomial value function [26], support vector machine [27], and ordinal regression [28].Gong et al. [29] proposed using a preference polyhedron to approximate the DM's value function by choosing the best and worst solutions from the current nondominated solutions.Cvetkovic and Parmee [30] suggested a method to integrate the DM's fuzzy preference information into the EMO algorithm by converting the linguistic terms into weights.These sorts of methods are interesting but complicated, especially when the number of objectives becomes large.In addition, using such an approach interactively increases the DM's cognitive load and it is hard to control the extent of the ROI.In [16], the biased distribution of solutions is achieved by setting different territory sizes in the territory-based evolutionary algorithm [31].In particular, a smaller territory leads to a higher resolution of solutions, and vice-versa.The size of the corresponding territory is progressively adjusted by the interaction with the DM.Note that this paper is one of the few that acknowledged the importance of providing information on the extent of the solution space while converging the ROI.The major drawback of this method comes from its diversity management, especially in a high-dimensional space.Due to the same reason, it might be difficult to control the extent of the ROI precisely.
The third category transforms the DM's preference information into some modified tradeoff relationship to compare solutions [32].Greenwood et al. [33] suggested an imprecisely specified multiattribute utility theory-based weighted sum method to obtain the ranking of objectives from some candidate solutions.
The fourth category [34] invites the DM to express his/her preference information by supplying two thresholds: 1) an absolutely satisfying objective value and 2) a marginally infeasible objective value.Afterward, each objective function of the original MOP is converted into a desirability function by using these thresholds as parameters.Then, an EMO algorithm is applied to optimize the desirability functions instead of the original objective functions.
The fifth class [35], [36] uses the outranking concept [37] to incorporate the DM's preference information.Specifically, by specifying some necessary parameters, a DM develops a fuzzy predicate that models the truth degree of the predicate "solution x is at least as good as solution y".
The last one uses aspiration-level vectors to represent the DM's desired values/levels for each objective he/she would like to achieve.As the first attempt, Fonseca and Fleming [38] suggested modeling the DM's preference as a goal to achieve, that is, the aspiration-level vector.Deb et al. [12], [39] combined the reference point, that is, aspiration-level vector, related methods with NSGA-II to guide the search process toward the ROI.In particular, solutions close to the given reference point have a high priority to survive to the next generation.In [40] and [41], the aspiration-level vector is used to help select the leader swarm in the multiobjective particle swarm optimization algorithm.Molina et al. [42] suggested a modified dominance relationship, called g-dominance, where solutions satisfying either all or none of the aspiration-levels are preferred over those satisfying some aspiration-levels.Said et al. [13] developed another modified dominance relationship, called r-dominance, where nondominated solutions, according to the Pareto dominance relationship, can be distinguished by their weighted Euclidean distances toward the DM-supplied aspiration-level vector.Recently, some decomposition-based methods also used the aspiration-level vector to incorporate the DM's preference information into the search process, for example, [43] and [44]- [46].Their basic idea is to use the aspiration-level vector as the anchor around which they try to obtain some reference points.Although by specifying aspiration-level vectors, a DM is able to guide the search toward the ROI directly or interactively even when encountering a large number of objectives, existing methods cannot approximate the solutions in the ROI and the boundary simultaneously.In addition, the control of the extent of the ROI is ad-hoc.

A. Overview
Reference points, as the basic components in the decomposition-based EMO algorithms, are usually generated in a structured manner, for example, the Das and Dennis's method. 2 [18] Fig. 1(a) shows an example of 91 evenly distributed reference points in a 3-D space.In this case, the DM has no preference on any particular region of the PF.These reference points are used to guide a decomposition-based EMO algorithm search for the whole PF.On the other hand, if the DM has elicited some preference information, for example, an 2 In [18], N = H+m−1 m−1 reference points, with uniform spacing δ = (1/H), are sampled from a canonical simplex m , where H > 0 is the number of divisions considered along each objective coordinate, and m is the number of objectives.aspiration-level vector, it is preferable that reference points can have a biased distribution toward the ROI accordingly.Bearing this consideration in mind, this section presents an NUMS by which we are able to change the originally evenly distributed reference points to be biased toward the ROI.Fig. 1(b) and (c) shows two examples of the biased reference points distribution after the nonuniform mapping.In the following paragraphs, we will describe the mathematical model of the NUMS in detail before showing its algorithmic implementations.

B. NUMS in 1-D Space
Let us start with a 1-D case.Considering the illustrative example shown in Fig. 2, the reference points generated in a structured manner are evenly distributed along the line starting from b 1 and ending at b 2 .Let us assume that the position of an evenly distributed reference point w obeys a uniform distribution whose probability density function (PDF) is defined as follows: where 0 ≤ ζ ≤ , = |b 2 − w p | is the distance between w p and b 2 .Here, w p is defined as the pivot point, which is the intersecting point between the reference line, connecting the DM-supplied aspiration-level vector z r and the origin, and the simplex m , to represent the ROI.When considering the DM's preference information, instead of a uniform distribution, it is preferable that the reference points have a biased distribution toward w p , that is, the closer to w p , the more reference points there are.The purpose of the NUMS is to shift w, originally generated by a structured manner, onto a new position w close to w p .Let us assume that the position of w obeys a nonuniform distribution whose PDF is defined as follows: where ξ = (δ/ ), δ = |b 2 −w | is the distance between w and b 2 .δ determines the position of w .η is a control parameter which will be further discussed in Section III-D.Note that 0 ≤ ξ ≤ 1 and δ gives the exact position of w along the line starting from b 1 and ending at b 2 .By equating the area under the probability curve of ψ e (ξ ) with that of ψ u (ζ ), we have (5)

C. NUMS in m-Dimensional Space
Now, we generalize the 1-D nonuniform mapping model into an m-dimensional case.Without loss of generality, let us consider a 2-D example shown in Fig. 3 for illustration.Similar to the 1-D case, the purpose of the NUMS in an m-dimensional case is to shift an evenly distributed reference point w onto w along the direction w p − w.For the ease of latter computation, we consider the opposite direction.That is to say, the NUMS shifts w p onto w along the direction w − w p .Accordingly, w is calculated as where • represents the 2 -norm and ρ is calculated as where = b−w p and δ is calculated based on (5) in which = w−w p .Note that w and w p are known a priori, while b is one of the intersecting points between the line connecting w p

D. Effects and Setting of η
Fig. 4 shows six function curves with various η settings.From this figure, we can infer that η controls the gradient of the PDF curve.ψ e (ξ ) is a decreasing function of ξ when η > 0; while it is an increasing function of ξ when η < 0. From Fig. 4, we also find that the function curve is more skewed with a larger η.According to the properties of power function, it is not difficult to understand that for a given and in (5), a larger η will result in a larger δ.In summary, η has the following two effects on the NUMS.
1) To push w toward w p , we need to set η > 0; otherwise, w will be shifted away from w p .2) With a large η, which results in a large δ, w has a large probability to be closer to w p after the nonuniform mapping; on the flip side, w will be closer to b.Based on the above discussions, we realize that η is able to control the extent of the biased reference points after the nonuniform mapping.However, due to the nonlinear property of the PDF in (2), it is far from trivial to choose the appropriate η beforehand that results in the expected extent of the ROI.Instead of tweaking η, by trial-and-error, with respect to the nonlinear mapping function, here we introduce an intuitively understandable way to control the extent of the ROI.Specifically, rather than a concrete extent of the ROI, it is more plausible for the DM to specify a relative quantity in practice.Here, we use τ (0 < τ ≤ 1) as the ratio of the surface area of the biased reference points proportional to the simplex m , as this quantity.As discussed in [3], under certain smoothness assumption, each reference point is supposed to correspond to a Pareto-optimal solution.Therefore, τ can also be regarded as relative ratio of the ROI's size with respect to the PF.Given τ , collected as additional preference information elicited by the DM, Theorem 1 gives a closed form for setting the corresponding η value.
Theorem 1: Given the relative extent τ (0 < τ ≤ 1) of reference points after the nonuniform mapping, compared to the simplex m , the η value in (2) is calculated as where α = (m/H) and β = 1 − τ .The proof of Theorem 1 can be found in Appendix A of the supplementary material. 3Fig. 5 shows three examples of biased reference points after the nonuniform mapping with different τ settings.Based on Theorem 1, we have the following corollary which provides the upper and lower bounds for setting τ .
Corollary 1: To make the extent of the biased reference points shrink, we need to set 0 < τ < 1 − (m/H).
The proof of Corollary 1 can be found in Appendix B of the supplementary material.In principle, comparing to the whole PF, the relative extent of the ROI can be any number between 0 and 1.However, Corollary 1 provides a restriction on τ in order to make the evenly distributed reference points shrink to the ROI; otherwise they will expand toward the boundary.It is worth noting that Theorem 1 and Corollary 1 are derived under the condition H > m.Otherwise, all reference points generated by the Das and Dennis's method should lie on the boundary of the simplex m .How to shift the reference points lying on the boundary will be described in the next section.

E. Boundary Preservation
Note that the NUMS described so far shifts the reference points, except those lying on the boundary of the simplex m , onto the ROI.The biased reference points try to guide a decomposition-based EMO algorithm not only search for the preferred solutions but also approximate those lying on the PF's boundary.In particular, the boundary solutions provide the DM more comprehension of the PF, e.g., the PF's general shape, the ideal and nadir points which can be useful for further decision-making.Nevertheless, if the DM is not interested in the boundary any longer, we can make a simple modification on the NUMS to shift the reference points lying on the boundary toward the ROI as well.Specifically, a reference point w b is considered lying on the boundary of m if and only if the following condition is met: where = 10 −6 is a small quantity and is determined according to (9).To shift w b onto the ROI, its new position after the NUMS is calculated as where ρ = τ × w b − w p .Note that the η value derived in Theorem 1 is under the consideration that the DM is willing to keep the boundary points.If the reference points lying on the boundary are shifted onto the ROI by the NUMS as well, the η value should be calculated according to Corollary 2. Corollary 2: If all reference points are shifted onto the ROI, the η value in (2) is calculated as where α = (m/H) and The proof of Corollary 2 can be found in Appendix C of the supplementary material.Accordingly, we should have a different upper and lower bounds for η as follows.
Corollary 3: If all reference points are shifted onto the ROI, we can set 0 < τ < 1.
The proof of Corollary 3 can be found in Appendix D of the supplementary material.Fig. 1(c) gives an example that all reference points have been shifted onto the ROI.

Algorithm 1: NUMS
Input: • DM supplied aspiration level vector z • Biased reference points on a canonical simplex m by Das and Dennis's method; 2 Find the pivot point w p of z r on m ; 3 if flag = 1 then // keep the boundary

Algorithmic Details
After describing the mathematical foundations of the NUMS, this section describes its algorithmic implementation whose pseudo-code is presented in Algorithm 1.First of all, N = H+m−1 m−1 reference points w 1 , . . ., w N are initialized via the Das and Dennis's method (line 1 of Algorithm 1).Afterwards, we find the pivot point (line 2 of Algorithm 1).Then, if the DM is interested in the boundary, we use Theorem 1 to compute the exponent η of the PDF in (2); otherwise, we use Corollary 2 to do so (lines 3-7 of Algorithm 1).During the main loop, for each reference point, we use ( 9) to Fig. 6.
Estimated ideal point z * is away from the DM-supplied aspiration-level vector z r .determine the position of the corresponding boundary point for the nonuniform mapping (line 9 of Algorithm 1).If the currently investigating reference point lying on the boundary of m and the DM is not interested in the boundary, we use (12) to determine the step-length for shifting this boundary reference point onto the ROI (line 11 of Algorithm 1); otherwise we use ( 5) and (7) to serve this purpose (lines 13 and 14 of Algorithm 1).At the end of this loop, we use (6) to calculate the new position of the biased reference point (line 15 of Algorithm 1).

G. Incorporation of the NUMS into Decomposition-Based EMO Algorithm
In principle, the NUMS can be readily incorporated into any decomposition-based EMO algorithm, e.g., MOEA/D and NSGA-III, in a plug-in manner.In particular, we only need to replace the reference points with the ones generated by the NUMS.However, for MOEA/D and its variants, the commonly used subproblem formulation, e.g., the Tchebycheff function, will only result in the population that are dominated by the ideal point [47], i.e., z * = (z * 1 , . . ., z * m ) T , where z * i = min x∈PS f i (x) for all i ∈ {1, . . ., m}, which is unknown a priori.Although the ideal point can be estimated by the currently evolving population, it is highly likely that the estimated ideal point is away from the DM-supplied aspiration-level vector, as shown in Fig. 6.Note that it is fine if the DM wants to approximate the boundary and the ROI together, since the solutions on the boundary can give the appropriate ideal point.Otherwise, the algorithm will be struggling to obtain acceptable solutions if the DM is only interested in the ROI.
To overcome this aforementioned drawback, we use the following subproblem formation in MOEA/D and its variants: where is the decision space and ρ is a sufficiently small positive number, which we set as 10 −6 as suggested in [39].As discussed in [48], the optimum of ( 14) must be a Pareto-optimal solution, and ρ is able to avoid the generation of weakly Pareto-optimal solutions.By using this subproblem formulation, we can expect that the search directions are heading toward the DM-supplied aspiration-level vector.

IV. EXPERIMENTAL STUDIES
To validate the effectiveness of the NUMS, 4 we incorporate it into three state-of-the-art decomposition-based EMO algorithms (i.e., MOEA/D [3], NSGA-III [17], and our recently proposed MOEA/D variant based on stable matching model, named MOEA/D-STM [49]) and three other state-of-theart preference-based EMO algorithms (i.e., g-NSGA-II [42], R-NSGA-II [12], and r-NSGA-II [13]).Here we use the simulated binary crossover (SBX) [50] and the polynomial mutation [51] as the reproduction operators.For the SBX, the crossover probability is set as p c = 1.0 and its distribution index is set as η c = 10; for the polynomial mutation, the mutation probability is set as p m = (1/n) and its distribution index is set as η m = 20.In the experiments, we choose the popular DTLZ1 to DTLZ4, and WFG41 to WFG48 test problems [52] to form the benchmark suite.Note that WFG41 to WFG48 problems are designed to have various complex PF shapes, e.g., sharp convex/concave, mixed shape, and disconnected PF segments.For DTLZ problems, m ∈ {3, 5, 8, 10}; while for WFG problems, m ∈ {2, 3, 5, 8, 10}.The settings of aspiration-level vectors used in our experiments are given in Appendix E of the supplementary material.The population size is set as N = 100 when m = 2; N = 92 when m = 3; N = 210 when m = 5; N = 360 when m = 8; and N = 660 when m = 10, respectively.As for the NUMS, the number of divisions is set as H = 13 when m = 3 and H = 6 when m = 5.When the number of objectives is larger than 5, we use a 3-layer method suggested in our recent work [53] to generate the initially evenly distributed reference points.In particular, we set H = 3 for each layer.In our experiments, the stopping criterion of a preference-based EMO algorithm is the number of function evaluations (FEs), where the detailed settings are given in Appendix E of the supplementary material.As for R-NSGA-II, the additional parameter , used in its -clearing procedure, is set according to [12], i.e., = 0.001 when m = 2 and = 0.01 otherwise.
It is nontrivial to quantitatively compare the performance of different preference-based EMO algorithms for approximating the ROI.Here we use our recently developed preferencebased R-HV [54] to serve the purpose.The basic idea of R-HV computation is to preprocess the obtained preferred solution set S, according to the DM-supplied z r , before using the hypervolume (HV) [55] for performance assessment.Interested readers can find the technical details from [54].Similar to HV, the larger is the R-HV value, the better is the quality of S for approximating the ROI.
In the experiments, each algorithm is performed 31 independent runs.In the data tables, we show the median and the interquartile range (IQR) of metric values for different problem instances with various aspiration-level vector settings.In particular, the best median metric values are highlighted in bold face with a gray background.To have a statistically sound conclusion, we carry out the statistical analysis as suggested in [56] to validate the statistical significance of the results.More detailed description of this statistical analysis framework is provided in Appendix F of the supplementary material.
To have a visual comparison, we also show the scatter plots and the parallel coordinate plots (PCP), in the supplementary material, of the final solutions obtained by different algorithms having the best R-HV value.

A. Experimental Results
Due to the page limit, we only discuss some results on DTLZ problems here.More comprehensive results and the discussion on WFG problems are in Section H of the supplementary material.From the results shown in Table I, we can clearly see that the decomposition-based EMO algorithms, i.e., MOEA/D-STM, MOEA/D, and NSGA-III, assisted by the NUMS are the best candidates for approximating the ROI of various test problems.Their superiority becomes more evident with the increase of the number of objectives.In the following paragraphs, we explain the results instance by instance.
Let us start from the DTLZ1 problem which has a linear PF shape, i.e., a hyper-plane intersects with each coordinate at 0.5.Note that DTLZ1 also has many local optima in its search space, which obstruct the convergence toward the global PF.
In the 3-objective case, all algorithms, except g-NSGA-II, are able to drive solutions to converge toward the PF.As the DM expects the ROI to be 20% of the whole PF, solutions found by the decomposition-based EMO algorithms assisted by the NUMS are the best candidates with respect to the DM's expectation.Fig. 7 shows the scatter plots of solutions obtained by all algorithms with respect to z r = (0.05, 0.05, 0.2).From this figure, we can see that solutions obtained by three decomposition-based EMO algorithms assisted by the NUMS are consistent and well distributed.In contrast, although the solutions found by R-NSGA-II are in the ROI, they crowd in a narrow region.In this case, R-NSGA-II cannot provide as many tradeoff alternatives as the decomposition-based EMO algorithms assisted by the NUMS.As shown in Fig. 7, solutions found by r-NSGA-II do not converge to the ROI.With the increase of the number of objectives, g-NSGA-II and r-NSGA-II have difficulty in driving solutions toward the PF due to the multimodal property of DTLZ1.As for R-NSGA-II, solutions are even more focused in the high-dimensional space as shown in Fig. 8, an 8-objective example with z r = (0.01, 0.02, 0.07, 0.02, 0.06, 0.2, 0.1, 0.01) T .Although the spread of the preferred solutions obtained by R-NSGA-II can be controlled by its parameter, there is no rule-of-thumb for tuning it to adapt to the DM' expected extent of the ROI.
DTLZ2 is a relatively simple test problem, where the objective functions of a Pareto-optimal solution x * satisfies:    Scatter plots of solutions on 3-objective DTLZ2 where z r = (0.2, 0.5, 0.6).m i=1 f 2 i (x * ) = 1.All algorithms do not have too much difficulty in driving solutions toward the PF.As the examples shown in Figs. 9 and 10, the performance of three decomposition-based EMO algorithms assisted by the NUMS are consistent.In contrast, the spread of the solutions obtained by the other three preference-based EMO algorithms is not fully controllable.In particular, solutions found R-NSGA-II and r-NSGA-II are very focused while those found by g-NSGA-II scattered in a wide region.The PF of DTLZ3 is the same as DTLZ2.But its search space contains many local optima which can make an EMO algorithm get stuck at any local PF before converging to the global PF.Similar to the observations in DTLZ1, g-NSGA-II cannot find any converged solutions in all 3-to 10-objective cases.The performance of the decomposition-based EMO algorithms assisted by the NUMS is very robust.It is interesting to note that, as the example shown in Fig. 11, solutions found by r-NSGA-II do not converge well to the ROI in the 3-objective case.This might be caused by the failure of its adaptive parameter control given a limited number of FEs.As shown in Fig. 12, we also notice that solutions found by r-NSGA-II do not converge to the PF when the number of objectives becomes large.
DTLZ4 also has the identical PF shape as DTLZ2.However, in order to investigate an EMO algorithm's ability to maintain a good distribution of solutions, DTLZ4 introduces a parametric variable mapping to the objective functions of DTLZ2.This modification allows a biased density of points away from f m (x) = 0.It is interesting to note that the performance of all these algorithms are similar to the DTLZ2.As the examples shown in Figs. 13 and 14, g-NSGA-II cannot drive all solutions converge to the PF due to the biased density of solutions.As shown in Fig. 13, some solutions found by r-NSGA-II are still drifted away from the PF when encountering an attainable aspiration-level vector, i.e., z r = (0.7, 0.8, 0.5), in the 3-objective case.

B. Summary of the Experimental Results
Based on the observations in Section IV-A, we summarize the comparisons between the decomposition-based EMO algorithms assisted by the NUMS and the other preference-based EMO algorithms as follows.
1) Solutions found by three decomposition-based EMO algorithms assisted by the NUMS are consistent.This is because the search directions of a decompositionbased EMO algorithm is determined by the reference points.By using the NUMS, the reference points are    Scatter plots of solutions on 3-objective DTLZ4 where z r = (0.7, 0.8, 0.5).
transformed from an even distribution to a biased distribution toward the DM-supplied aspiration-level vector.In contrast, the driving force of the other preferencebased EMO algorithms is to find solutions close to the DM-supplied aspiration-level vector.Since this "closeness" by itself is vague, it brings uncertainty to the search process.
2) For the decomposition-based EMO algorithms assisted by the NUMS, the extent of the approximated ROI is controlled by the DM in an intuitively understandable manner.For the other preference-based EMO algorithms, the approximated ROI can be any crowd of solutions "close" to the DM-supplied aspiration-level vector.Although there are some parameters that control this extent, there is no rule-of-thumb to tweak those parameters.
3) As shown in the proof-of-principle results shown in the supplementary material, the NUMS can help a decomposition-based EMO algorithm not only find solutions in the ROI but also those lying on the boundary of the PF.In contrast, The other preference-based EMO algorithms can only approximate the ROI.4) Different from the other preference-based EMO algorithms, the NUMS does not incur additional computations to the baseline algorithm.As introduced in Section III-G, we only need to change the distribution of the reference points to be biased toward the ROI.As shown in Table V of the supplementary material, the average CPU time costs of the NUMS assisted algorithms are almost the same as the baseline algorithms.
V. CONCLUSION This paper present a systematic way to incorporate the DM's preference information into the decomposition-based EMO methods in either a priori or interactive manner.In particular, the DM's preference information is modeled as an aspiration-level vector which represents the DM's expected value on each objective.Our basic idea is an NUMS that transforms the originally evenly distributed reference points into a biased distribution.In particular, the closer to the DM specified aspiration-level vector, the more reference points in view of their higher relevance to the DM's preference information.Different from the existing literature, the ROI's size is fully controllable and intuitively understandable according to a quantitative definition.To facilitate the interactive decision-making process, our proposed NUMS is able to preserve the ones located on the boundary as given the DM's requirements.By incorporating the NUMS into some decomposition-based EMO algorithms, i.e., MOEA/D-STM, MOEA/D, and NSGA-III, its effectiveness is validated by proof-of-principle experiments and comparative studies with other state-of-the-art preference-based EMO algorithms on a variety of benchmark problems with 2 to 10 objectives.
It is clear that the distribution of the biased reference points is determined by the transformation function defined in (2).One direct extension of this paper is to use some other distribution functions that are tailored according to the DM's requirements.As discussed in Section II, there are several other ways of eliciting the DM's preference information.The other extension of this paper is the adaptation of the NUMS to other types of preference model.To further facilitate the interactive process, it is worth considering the combination of human computer interaction techniques [57] and the preference-based EMO.Moreover, discrete and mixed variable optimization problems are ubiquitous in real-world applications, e.g., scheduling [58], [59].It is interesting to study the application of the NUMS for finding DM preferred solutions in those cases.

TABLE I COMPARISON
RESULTS OF MEDIAN R-HV VALUES AND THE IQR OBTAINED BY SIX PREFERENCE-BASED EMO ALGORITHMS ON DTLZ1 TO DTLZ4 PROBLEMS WITH UNATTAINABLE AND ATTAINABLE ASPIRATION LEVEL VECTORS