Evolutionary Search with Multiple Utopian Reference Points in Decomposition-Based Multiobjective Optimization

Decomposition-based multiobjective evolutionary algorithms (MOEA/Ds) have become increasingly popular in recent years. In these MOEA/Ds, evolutionary search is guided by the used weight vectors in decomposition function to approximate the Pareto front (PF). Generally, the decomposition function will be constructed by the weight vectors and the reference point, which play an important role to balance convergence and diversity during the evolutionary search. However, in most existing MOEA/Ds, only one ideal point is used as the reference point for all the evolutionary search, which is harmful to search the entire PF when tackling the problems with difficult-to-approximate PF boundaries. To address the above problem, this paper proposes an evolutionary search method with multiple utopian reference points in MOEA/Ds. Similar to the existing MOEA/Ds, each solution is associated with one weight vector, which provides an evolutionary search direction, while the novelty of our approach is to use multiple utopian reference points, which can provide evolutionary search directions for different weight vectors. Corner solutions are used to approximate the nadir point and then multiple utopian reference points for evolutionary search can be constructed based on the ideal point and the nadir point, which are uniformly distributed on the coordinate axis or planes. The use of these utopian points can prevent solutions to gather in the same region of PF and helps to strike a good balance of exploration and exploitation in the search space. The performance of our proposed algorithm is validated on tackling 16 recently proposed test problems with difficult-to-approximate PF boundaries and empirically compared to eight state-of-the-artmultiobjective evolutionary algorithms. The experimental results demonstrate the superiority of the proposed algorithm on solving most of the test problems adopted.


Introduction
In some real-world applications [1][2][3], there are many problems involving the optimization of multiple (often conflicting) objectives, which are termed multiobjective optimization problems (MOPs).In general, a MOP can be defined mathematically as follows: min ∈Ω  () = ( 1 () , . . .,   ()) , where Ω ∈ R  is the n-dimensional decision (variable) space (n is the number of decision variables),  = ( 1 , . . .,   )  is a decision vector in Ω, and the objective vector F(x) consists of m objective functions  1 (), . . .,   ().Usually, the output of MOPs will be a set of tradeoff solutions for all the objectives, called Pareto-optimal solutions (PS).The mapping of PS on the objective space is called Paretooptimal front (PF).Traditional mathematical methods are difficult to obtain a subset of PS in a single run, while multiobjective evolutionary algorithms (MOEAs) can use a population composed of many individuals to approximate the PF, which are popular for tackling MOPs [4][5][6].Compared with traditional mathematical methods, MOEAs can well solve various kinds of MOPs without their differentiable and derivable characteristics.Due to these advantages, MOEAs have been very popular in recent years [7][8][9][10].
Most of MOEAs can be classified to three categories based on the selection mechanisms: (1) Pareto-based MOEAs [11][12][13][14]; (2) indicator-based MOEAs [15,16]; and (3) decomposition-based MOEAs [17,18].In fact, there are only two main tasks in the optimization process of MOEAs.The first one is to ensure convergence, i.e., MOEAs should have the convergence power to push its population 2 Complexity towards the PF, while the second one is to maintain diversity, i.e., MOEAs should search a set of solutions uniformly distributed along the PF.Despite the equal importance of both tasks, most state-of-the-art MOEAs follow the principle of convergence first and diversity second, such as NSGA-II [11], SPEA2 [12], MDMOEA [13], and CA-MOEA [14].Specifically, in MDMOEA [13], a new dominance relation, i.e., L-dominance, is proposed to better balance convergence and diversity when solving many-objective optimization problems.For MOPs with irregular PFs, CA-MOEA [14] uses the clustering method to detect the distribution of individuals in the objective space and then selects individuals according to cluster centers in the last nondominated front.In these MOEAs [11][12][13][14], nondominated solutions are selected in priority into the next generation for ensuring convergence, which will induce the diversity loss in sense.
One popular decomposition-based framework for MOEA (MOEA/D) proposed by Zhang et al. [17] has been widely studied in recent years.Moreover, some improved MOEA/D variants [19][20][21][22] are designed to better balance convergence and diversity.In MOEA/D-DE [19], the offspring can replace at most n r individuals in its neighborhood to prevent diversity loss, where n r is preset as a small positive integer.For example, in MOEA/D-STM [20] with a stable matching (STM) model, each individual is assigned to one subproblem based on their mutual preference information, trying to balance convergence and diversity during the evolutionary search.Moreover, MOEA/D-ATO [21] adaptively assigns the suitable combination of mating neighborhood size and reproduction operator to each subproblem at different searching stage.In IMOEA/DU [22], a crossover operator based on uniform design and selection strategy based on decomposition are used to help MOEAs to improve the search efficiency.Recently, in MOEA/HD [23], subproblems are layered into different hierarchies, aiming to better balance convergence and diversity during evolution.In MOEA/D-CPDE [24], the best individual in each generation will be regarded as a seed, which will evolve two individuals by cloud generator.In DAA [25], an efficient decomposition-based archiving approach is proposed to dealing with multiobjective optimization problems.On the other hand, recent research studies have indicated that the coevolution mechanism can significantly improve the performance of MOEAs [26,27].The basic idea in these MOEAs is to divide a complex system into several modules, which are then evolved separately [26,28].For example, MOEA/D-M2M [29] decomposes a MOP into a set of simple MOPs and then solves these simple MOPs in a collaborative manner.Similarly, E-CODBA [30] exploits coevolution, decomposition, and energy laws to come up with good solutions within an acceptable execution time, aiming to solve some real-world applications, i.e., combinatorial bilevel problems.Furthermore, some other heuristic searching algorithms, i.e., particle swarm optimization [31] and immune algorithm [32] are also developed to solving optimization problems.
However, it is found in [33][34][35][36][37][38][39][40][41][42][43][44][45] that most of MOEA/Ds adopt one single reference point for all the evolutionary search directions, which is harmful for searching the entire PF when tackling MOPs that are characterized with difficultto-approximate PF boundaries [46].To solve this kind of MOPs, in this paper, an evolutionary search method with multiple utopian reference points is proposed for MOEA/Ds.Each solution will be associated with one weight vector, which provides an evolutionary search direction.However, in our approach, multiple utopian reference points are used to provide evolutionary search directions for different weight vectors.To summarize, our approach consists of the following features.
(1) Different from the existing decomposition approaches, multiple utopian reference points are used to replace the ideal reference point in our proposed decomposition function.The use of these utopian reference points can provide multiple search directions, which will guide the solutions towards their respective utopian reference points.Therefore, our proposed approach can prevent the solutions to gather in the same region of PF and help to strike a good balance of exploration and exploitation in the search space.
(2) These utopian reference points are uniformly distributed on the coordinate axis or planes, which are used to provide different evolutionary search directions for different weight vectors.Moreover, they are adaptively adjusted according to the evaluation of the ideal point and the nadir point in each generation.By this way, the evolutionary search directions are adaptively adjusted, which guide solutions to approximate the entire PF during the evolutionary process.
The rest of this paper is organized as follows.Section 2 provides some useful preliminaries of this paper.Section 3 describes the technical details of the proposed evolutionary search method with multiple utopian reference points and the entire algorithm.The experimental results and related discussions are presented in Sections 4 and 5, respectively.Finally, Section 6 concludes this paper and provides some future directions.

Preliminaries
. .Tchebycheff Decomposition.In decomposition approaches, the scalar subproblem is formulated as a weighted aggregation of all the objectives by decomposition function.The TCH decomposition [17] is a commonly used decomposition function in MOEA/Ds, which can be constructed by a weight vector and a reference point.It can be mathematically defined by the following formulation: where  * is used as the reference point for the weight vector  = ( 1 , . . .,   )  ,   ≥ 0 for all  ∈ {1, . . ., } and ∑  =1   = 1.Note that   is set to be a very small number, say 10 −6 , in case   = 0. Actually, in most existing MOEA/Ds [19][20][21][22][23][24][25], the ideal reference point is usually used as the initial point of all weight vectors, which is estimated by using the best value of each objective found so far during evolution.Note that all the solutions should be dominated by the ideal Figure 1: Tchebycheff decomposition with the ideal reference point for the cases: (a) reference point in the objective space.Thus, it is obvious that   () is larger than  *  for all  ∈ {1, . . ., } in the objective space.In the recent study [47], the geometric property of the subproblem's objective function in the TCH decomposition has been studied in detail.It can be summarized as follows.
In the case that F(x) is not located in L 1 , the following proposition depicts the geometric property of   (() | ,   ).
The above two propositions demonstrate the geometric property of TCH decomposition in a serious mathematical proof.More details of proofs have been systematically discussed and analyzed in [47].

. . Different Types of Reference Points.
There have been numerous studies of decomposition approaches [35][36][37][38][39][40][41][42][43][44][45] to use different types of reference points for providing evolutionary search directions.According to the position of reference point relative to the true PF in the objective space, there have three main kinds of reference points in existing MOEA/Ds.For clarity, we give the mathematical definitions of three kinds of reference points as follows.
Definition (the ideal point).The ideal point   = (  1 , . . .,    )  is assigned by the best value of each objective in the true PF, as defined by Definition (the nadir point).The nadir point   = (  1 , . . .,    )  is assigned by the worst value of each objective in the true PF, as defined by Definition (the utopian point).The utopian point   = (  1 , . . .,    )  is assigned by an infeasible objective vector that dominates the ideal point, as defined by where   > 0 for all  ∈ {1, . . ., }.As shown in Figure 2(a), the ideal point is the most commonly used reference point for search directions (as indicated by the weight vectors), which is estimated by using the best value of each objective found so far during evolution.The use of the ideal reference point can speed up convergence, but it will also cause the loss of diversity.The representative MOEA/Ds, such as MOEA/D [17] and its variants MOEA/D-DE [19], MOEA/D-STM [20], MOEA/D-ATO [21], IMOEA/DU [22], MOEA/HD [23], MOEA/D-CPDE [24], and DAA [25], all adopt the ideal point as the reference point for all the search directions.However, it can be observed from Figure 2(b) that the current population is pushed away from the nadir reference point as far as possible.Therefore, the solutions in the boundary parts of the PF can be easier to be exploited if the nadir point rather than the ideal point is used as the reference point.The existing MOEA/Ds [36][37][38][39], i.e., MOEA/D-TPN [37], MOEA/D-MR [38], and PAEA [39], adopt the nadir reference point, which show a better performance to maintain diversity.Moreover, in Figure 2(c), the utopian reference point can be used to explore the new areas in the boundary parts of the PF when tackling the MOPs with difficult-to-approximate PF boundaries, which has been validated in some studies [40][41][42][43], like MOEA/D-UP [42] and Global WASF-GA [43].Furthermore, there are some attempts [44,45] to use multiple reference points, which are used to ensure good diversity by guiding different solutions to approximate the true PF.As shown in Figure 2(d), multiple reference points rather than the ideal reference point are used to provide search directions for different solutions, which are uniformly distributed along the approximated PF.

Complexity
. .e Motivations of Our Method .In decomposition approaches, each subproblem is built by using a weight vector and a reference point, which provides an evolutionary search direction for its associated solution.As it is impossible to know the exact ideal point beforehand, an approximate ideal point is often used as the reference point for all the search directions in MOEA/Ds [17,[19][20][21][22][23][24][25], which is obtained by using the best value of each objective found so far during evolution.For example, in MOEA/D [17], a set of uniform weight vectors decomposes a MOP into a number of subproblems, which are assigned with the same computational resource in evolutionary search.However, as pointed out in MOEA/D-DRA [48] and MOEA/D-IRA [49], these subproblems are not equally important and the computational resources should be allocated according to the difficulties of subproblems.Especially, when tackling the problems with difficult-to-approximate PF boundaries, solutions in the middle part of the PF are easier to be found, while that in the PF boundaries are difficult to be searched.Thus, the solutions with good convergence in the middle part of the PF will replace the ones close to the PF boundaries with a high probability, which will cause overcrowdedness at the middle part of the PF.As observed in Figure 3(a), the newly generated offspring close to the ideal reference point have a better aggregated function value for their corresponding subproblems, which will replace these parent solutions with poor convergence into the next generation.Once solutions are all trapped into the middle part of the PF, the boundary solutions will be almost impossible to be found due to the restriction of search directions in decomposition approaches.Similarly, as shown in Figure 3(b), although the nadir point is used as the reference point for search directions, the loss of diversity in the PF boundaries still exists when tackling the MOPs with difficult-to-approximate PF boundaries, while its middle part of the PF can be easily found.
Therefore, for the problems with difficult-to-approximate PF boundaries, this paper proposes an evolutionary search method with multiple utopian reference points (MUP) in MOEA/Ds, which can avoid overcrowdedness at the middle part of the PF and search the entire PF.Similar to the existing MOEA/Ds, each solution is associated with one weight vector, which can provide an evolutionary search direction.However, the main difference of our approach with the existing MOEA/Ds is that multiple utopian reference points are used to provide evolutionary search directions for different weight vectors.As these utopian reference points are uniformly distributed along all the coordinate axis or planes that are constructed by the ideal reference point and the nadir reference point, a much wider breadth-diversity [46], which indicates how widely the solutions are distributed over the PF, can be obtained in the evolutionary process.Our approach true PF approximate solutions in T 2 (T 2 ≫T 1 ) generation approximate solutions in T 1 generation in T 2 (T 2 ≫T 1 ) generation true PF approximate solutions approximate solutions in T 1 generation can help to explore the boundary areas on the PF due to the use of multiple utopian reference points.At the same time, all the solutions are guided by their respective utopian reference points to approximate the true PF, which prevents them to gather at the same region of the PF and helps to strike a good balance of exploration and exploitation in the search space.
. .Related Work of MOEA/Ds.As recently revealed in [35], the reference point providing evolutionary search directions plays a crucial role in balancing convergence and diversity of MOEA/D.Actually, different types of reference points, such as the ideal point, the nadir point, and the utopian point, could have different impacts on the search behaviors of MOEA/Ds.Recently, the nadir point is used in [36][37][38][39] as the reference point for strengthening diversity, while the utopian point is adopted in [40][41][42] as the reference point for all the evolutionary search directions in MOEA/Ds.Both of a utopian point and a generalized nadir point are used simultaneously as the reference points in [43], which indicates that MOEA/Ds can achieve a much wider distribution to approximate the PF if the utopian point rather than the ideal point is treated as the reference point.In the above-mentioned MOEA/Ds, there are at most two reference points adopted to provide search directions.However, in R-MOEA/D [44], a new scalarizing function using a series of new reference points derived from a given reference point in the preference model is proposed, which ensures a much wider spread of approximate solutions.However, this approach only works well on MOPs with two objectives, which cannot be extended for MOPs with more objectives.Similarly, a new scalarization method is proposed in [45] to get a more uniform distribution on the PF.To summarize, the appropriate selection of reference points for evolutionary search directions can improve the performance of MOEA/Ds, which will better balance exploration and exploitation during the evolutionary process.

The Proposed Method
In this section, the generation of multiple utopian points (MUP) is first introduced in Section 3.1.Then, the decomposition approach of using MUP is described in Section 3.2, while the incorporation of MUP method into MOEA/D is clarified in Section 3.3.[40][41][42][43], a utopian point is usually set as Definition 5.It is clear that the utopian point dominates the ideal point according to this definition.As illustrated in Figure 4, a utopian point is only located in the narrow dominated areas between the ideal point and the origin point, as named by the narrow utopian area in Figure 4.

. . Generation of Multiple Utopian Points. In existing MOE-A/Ds
In this paper, one more generalized formula of the utopian point is given, as mathematically defined by It can be observed from ( 9) that the utopian area can be extended to a much wider objective space using this general definition.Thus, in other words, the area except for the dominated area of PS can be served as the utopian area.
In our proposed decomposition method with MUP, the weight vectors and the reference points are still the critical components.Given a set of uniformly distributed weight vectors  = { 1 , . . .,   } that is generated by using Das and Dennis's [50] systematic approach and an ideal reference point that is estimated by using the best value of each objective found so far during evolution, a MOP is decomposed into a set of scalar optimization subproblems by using a decomposition function.Unlike the existing MOEA/Ds [19][20][21][22][23][24][25] with only an ideal point, a set of utopian reference points is derived from the ideal point in this paper.Each of these utopian reference points { shifting the ideal reference point along the coordinate axis, as described in ( 10) where The ideal point and the nadir point decide the range of these newly generated utopian points as illustrated in (11).In this paper, we use the corner solutions [50] to compute the nadir point.In fact, these M corner solutions lying on the objective axes have a good ability in representing the boundary of each objective on the PF, which can be obtained by simultaneously minimizing −1 objectives.As suggested in [51], the Euclidean distance to the ideal point is considered in the paper, as follows: Then, the corner solution   can be obtained as follows: Thus, the nadir point can be computed by In our proposed method, each utopian point is associated with a weight vector as illustrated in (10).Although the ideal reference point is inaccurate, the case that solutions are overcrowded around the ideal reference point will not occur as multiple utopian reference points are used to provide evolutionary search directions for different weight vectors.As shown in Figure 5, for MOPs with two objectives, these newly generated utopian reference points are uniformly distributed in the coordinate axis, which are bounded between the ideal point and the nadir point.However, for MOPs with more than two objectives, these utopian reference points will be distributed in the coordinate planes.According to (11), these utopian reference points are adaptively adjusted according to the change of the ideal point and the nadir point during the evolutionary process.
. .Decomposition with MUP.In our proposed decomposition approach, a set of uniformly distributed utopian reference points replaces the ideal reference point to guide the search behaviors of individuals.Due to the replacement of reference point, the classical TCH decomposition function in (2) is not suitable for balancing convergence and diversity in our approach.Therefore, to strike a good balance of exploration and exploitation in the search space, a novel decomposition function is proposed by applying some mild modifications on the classical TCH decomposition function, as mathematically defined by min where the utopian point is used as the starting point for each search direction generated by (10).Unlike the classical TCH decomposition function in (2), the utopian point rather than the ideal point is used to provide evolutionary search direction for each associated individual, and the absolute operator is removed to allow the objective value of the solution to be smaller than the utopian reference point in some dimensions, i.e.,   () <    , where ∃ ∈ {1, . . ., }.Without using the absolute operator in (15), it is obvious that there are two situations in our proposed decomposition function.Actually, there is no difference with the classical TCH decomposition function when the solution is dominated by the utopian reference point.However, if the solution is not dominated by the utopian reference point (i.e., their Figure 6: Tchebycheff decomposition with the utopian reference point. relationship is nondominated), the objective value of the solution that is smaller than the utopian reference point in one dimension will not play any effect on the aggregated function value of (15), whereas the objective values larger than the utopian reference point in other dimensions will decide the aggregated function value.Of course, according to the definition of utopian point, solutions are impossible to dominate the utopian reference point.Similar to the study in [47], the geometric property of the subproblems' objective functions is investigated when using (15) for decomposition.
Proposition 6.Let   = (  1 , . . .,    ) be a utopian reference point of ( ) and  = ( 1 , . . .,   )  be a positive weight vector.If a given objective vector () = ( 1 (), . . .,   ()) is in the line: In this proposition, there are two different cases according to the dominance relation between the solution and the utopian reference point in the objective space.As shown in Figure 6(b),   ( 1 ) ≥    for all  ∈ {1, . . ., }.In this case, there is no difference between our method and the classical TCH approach on the geometric property.Thus, Proposition 7 can be tenable according to Proposition 2. However, there is another case that the solution is nondominated by the utopian reference point.As shown in Figure 6(b),   ( 2 ) <    for at least one dimension ∃ ∈ {1, . . ., }.The proof is given as follows.
Proof.According to the conditions of (), we get →            () −            1 Remark .In geometry,   (( 2 ) | ,   ) is equal to the l 1 -norm of () −   as shown in Figure 6(b).() is the intersection point between the line L 2 and the contour line of   (( 2 ) | ,   ) in the objective space, as follows: since →   ( () | ,   ) (2) For example, given ( 2 ) = (1, 0.5),   = (0, 1), and  = (1/3, 2/3), we can calculate . .Incorporation of MUP Method into MOEA/D.In this section, we introduce the framework of our evolutionary search method with MUP in MOEA/Ds, i.e., MOEA/D-MUP, whose pseudocode is given in Algorithm 1.At the beginning, a population with N solutions  = { 1 , . . .,   } is initialized in line (1), and a set of uniformly distributed weight vectors  = { 1 , . . .,   } is generated by using the Das and Dennis [50] systematic approach in line (2).Then, these N solutions are randomly allocated to associate each of N subproblems.Moreover, for each weight vector   ( ∈ {1, . . ., }), its T closest weight vectors are assigned as its neighbors () in lines (3)-( 5), where 1 ≤  ≤ .Then, the ideal point will be initialized by finding the best value of each objective and the nadir point will be initialized using the corner solutions in lines (6)- (7), which have been introduced in Section 3.1.Afterward, like other MOEA/Ds [19][20][21][22][23][24][25], parent solutions are randomly selected from the neighborhood with a high probability  in line (11).In this case, the more similar child solutions will be produced to have a good exploitation on the local region.Otherwise, parent solutions are randomly selected from the whole population with a low probability 1 −  in line (13).In this case, a very wide range of child solutions could be generated by evolutionary operators due to the dissimilarity among these parent solutions.After that, two parent solutions are randomly selected in line (15), which are used to generate a new offspring by evolutionary operators in lines ( 16)- (17), i.e., differential evolution (DE) and polynomial-based mutation [19].Next, the newly generated offspring will be used to update the ideal point in line (18), and then the neighborhood of i-th subproblem is also updated by the offspring in line (19).As shown in Algorithm 2, each subproblem is associated with a unique utopian reference point in line (2), which will be used to compute the aggregated function values of solutions.Then, as shown in lines (3)-( 5) of Algorithm 2, the compared parent solution will be replaced by the newly generated offspring if the offspring has the smaller aggregated function value.After updating all the neighborhoods of subproblems, the nadir point will be reevaluated by recalculating the corner solutions of the updated population in line (21).If the stopping criterion is not satisfied in line (8), the above evolutionary process will be repeated in lines (9)- (21); otherwise, the final population will be outputted in line (23).

Experimental Settings
. .Test Instance.Sixteen test problems [46] (i.e., F1-F16) featured by difficult-to-approximate PF boundaries are considered as the benchmark problems in our experimental studies.F1-F8 are two-objective optimization problems, while F9-F16 are three-objective optimization problems.As suggested in

Input:
B(i): the neighboring of the i-th subproblem y: a new generated solution z ide : the ideal point z nad : the nadir point (1) for each index  ∈ () do (2) Compute the utopian point   associated with j-th subproblem using Eq. ( 10) =  (5) end (6) end Algorithm 2: UpdateNeighboringSolutions.[46], the number of decision variables is set as 31 in twoobjective optimization problems (F1-F8).However, in threeobjective optimization problems (F9-F16), the number of decision variables is set as 32.The true PFs of F1-F7 are the same, which can be represented by a straight line in the objective space.The true PF of F8 is convex.Similar to F1-F8, the geometry of the true PFs of F9-F15 is a hyperplane except that the PF of F16 is convex.The design details and control parameters for F1-F16 are introduced in [46].
. .e Compared Algorithms.Three classical MOEAs (NSGA-II [11], SMS-EMOA [15], and MOEAD/D-DRA [48]) and a modified variant of MOEA/D-DRA (MOEA/D-DRA-UT [46]) are chosen for comparison to assess the performance of the proposed algorithm on F1-F16.NSGA-II is a famous dominance-based MOEA, while SMS-EMOA is a well-known indicator-based MOEA using the hypervolume metric for population diversity maintenance.Moreover, MOEA/D-DRA is a typical decomposition-based MOEA, in which all subproblems are dynamically allocated with Table 1: Parameter settings of all the compared algorithms.
(2) Reproduction operators: To be fair, we use DE operator and Polynomial-based mutation for offspring generation in all the compared algorithms except that MOEA/D-M2M employs another recombination operator [53] as suggested in [29].As recommended in [19], we set  = 1.0 and  = 0.5 in DE operation,  = 20 and   = 1/ in Polynomial-based mutation operator, where  is the number of decision variables.
(3) The stopping condition and the number of runs: The maximal number of function evaluation is set to 300,000 for all test problems.All the algorithms are run 30 times independently on each test problem.
. .Performance Metrics.To compare the performance of all algorithms on test instances, the inverted generational distance (IGD) indicator and the hypervolume indicator (I H ) are used in our experiments.A set of solutions uniformly sampled along the true PF is needed for computing the former indicator, i.e., IGD.However, for the later, i.e., I H , a reference point is necessary.In our experiment, all the nondominated solutions are selected from the final population, which are used to compare the performance of all algorithms in terms of IGD values and I H values.
(1) IGD [54]: Let  * be a set of solutions which are uniformly sampled along the PF in the objective space.The approximate set (denoted as A in ( 23)) composes of a set of nondominated solutions, which are selected from the final population obtained by the corresponding algorithm.Then, the IGD value can be computed as follows: where and m is the number of objectives.As suggested in [46], 500 representative points are uniformly sampled from the PF to form  * for F1-F8, and 1,000 points are chosen as  * for F9-F16.
(2) I H [55]: Let  = { 1 , . . .,   }  be a reference point in the objective space.The approximate set (denoted as A in ( 24)) composes of a set of nondominated solutions, which are obtained by the corresponding algorithm in the objective space.Then, the I H value is computed as follows: (24) subject to  ≺ , ∀ ∈ , where V(⋅) is the Lebesgue measure.As suggested in [29], the reference point  = (1.0, . . ., 1.0)  is considered in our Complexity 11 experiments.Please note that all the nondominated solutions in A should dominate the reference point in the objective space.

Experimental Studies
In this section, we present and discuss the comparisons of MOEA/D-MUP with other state-of-the-art algorithms.Three classical algorithms (NSGA-II, SMS-MOEA, and MOEA/D-DRA), and a competitive algorithm with the utopian reference point (MOEA/D-DRA-UT) are chosen for comparison.Experimental results obtained by the above four algorithms in terms of the two common indicators, i.e., IGD and I H , are shown in Tables 2 and 3 4 and  5, respectively.Please note that the final nondominated solutions with the median IGD values obtained by different algorithms in 30 runs are plotted in Figures 8 and  9.
. .Quantitative Comparisons.According to Wilcoxon's rank sum test with a significant level of 0.05, the experimental results obtained by the nine MOEAs in terms of IGD metric and I H metric are shown in Tables 2-5, respectively.In general, it can be observed that MOEA/D-MUP performs better than the other eight algorithms on most test problems, i.e., F1-F16.More specifically, it can be seen from Tables 2 As shown in Tables 2 and 3, we can conclude that MOEA/D-MUP is significantly better than NSGA-II and SMS-EMOA in terms of mean IGD and I H values on all test problems adopted, i.e., F1-F16.In fact, the fast nondominated sorting approach always tends to choose these solutions in the middle part of the PF, which is caused by the characteristics of the problems that the solutions located in the middle part of the PF always can dominate other boundary solutions during the evolutionary process.Thus, as shown in Tables 2 and 3, NSGA-II always obtains very poor performance on all test problems, i.e., F1-F16.Similarly, the update strategy of neighborhood used in most MOEA/Ds is executed based on the aggregated function value to replace parent solutions, which also follows the strategy of convergence first and diversity second.Thus, overcrowdedness in the middle part of the PF still occurs in most MOEA/Ds.For example, although a dynamic computational resource allocation strategy is adopted in MOEA/D, MOEA/D-DRA still fails to cover the boundaries of these PFs when solving F1-F16.As shown in Tables 2 and  3, in terms of mean IGD and I H values, MOEA/D-MUP achieves better performance than MOEA/D-DRA on all test problems except that they have no significant difference on F7 regarding the mean IGD value.For MOEA/D-DRA-UT, more boundary areas can be found in the evolutionary search process, while MOEA/D-DRA-UT is still beaten by our proposed MOEA/D-DRA-UT on most test problems.These results validate that solutions in the middle part of the PF are likely maintained in the environmental selection mechanisms of these traditional MOEAs when tackling the problems with difficult-to-approximate PF boundaries.It can be observed from Tables 2 and 3, MOEA/D-MUP significantly outperforms MOEA/D-DRA-UT on most test problems.For example, on F1-F6, MOEA/D-MUP is significantly better than MOEA/D-DRA-UT in terms of mean IGD and I H values.However, on F7, no algorithm can approximate the entire PF due to the extreme difficulty to search the PF boundaries.Quantitatively, MOEA/D-MUP is slightly beaten by MOEA/D-DRA-UT in terms of mean IGD and I H values on F7.In addition, the performance of MOEA/D-DRA is similar to MOEA/D-MUP in terms of IGD on F7.This is mainly because the used dynamic resource allocation strategy can allocate more computational resources to search the undiscovered boundaries of the true PF on F7.Quantitatively, for F8, MOEA/D-MUP is slightly beaten by MOEA/D-DRA-UT in terms of IGD metric, while they have similar results in terms of I H metric.As discussed above, the use of the utopian reference point instead of the ideal reference point is helpful to exploit the boundaries of the PF.However, in MOEA/D-MUP, multiple utopian reference points are used to provide different evolutionary search directions, which fails to cover the entire PF of F8.This is mainly because the inaccuracy of the nadir point leads to the wrong judgment about the upper bounds of the entire PF of F8.Furthermore, For F9-F16, MOEA/D-MUP achieves a better performance than MOEA/D-DRA-UT in terms of mean IGD metric, except that they are similar on F15.In terms of mean I H metric, MOEA/D-MUP performs much better than MOEA/D-DRA-UT on F10-F12 and F16, and they obtain statistically similar results for the rest three-objective test problems.In general, these results confirm the superiority of MOEA/D-MUP for tackling these problems with difficult-to-approximate PF boundaries.
As shown in Tables 4 and 5 4 and 5, MOEA/ D-MUP is only beaten by MOEA/D-M2M on F7, and they obtain statistically similar results on F13 and F14 in terms of I H . On other problems, MOEA/D-MUP is significantly better than MOEA/D-M2M.In terms of IGD, MOEA/D-STM is similar to MOEA/D-MUP on F7, while MOEA/D-STM is significantly worse than MOEA/D-MUP on the rest problems.When compared to MOEA/D-AWA and MOEA/D-PaS, MOEA/D-MUP performs much better on all test problems.Furthermore, when compared to MOEA/D-STM, MOEA/D-M2M, MOEA/D-AWA, and MOEA/D-PaS, MOEA/D-M2M attains the best mean IGD and I H values on all two-objective test problems.From these comparison results, it can be concluded that our proposed MOEA/D-MUP obtains the best performance, and the decomposition strategy used in MOEA/D-M2M achieves some success for some problems with difficult-toapproximate PF boundaries, as it can effectively maintain the boundary solutions into the next population.To conclude, our proposed MOEA/D-MUP can attain significantly better convergence and breadth-diversity than other compared MOEAs on most test problems adopted.For the problems with difficult-to-approximate PF boundaries, diversity should be more considered during the evolution.
To visually show the convergence speeds of all the compared algorithms during the evolution, the curves of IGD value over the number of evaluations are depicted in Figure 7, when solving two biobjective optimization problems (F1 and F6) and three triple-objective optimization problems (F10, F13, and F16).As shown in Figures 7(a . .Qualitative Comparisons.Due to page limitations, we select two competitive algorithms, i.e., MOEA/D-DRA-UT and MOEA/D-M2M, to qualitatively demonstrate the superiority of our proposed algorithm on solving most of the test problems adopted.For example, some representative twoobjective problems, i.e., F1, F3, F5, F6, and F8 are shown in Figure 8 and some three-objective problems are shown Figure 9, respectively.As shown in Figure 8, it is obvious that our proposed algorithm MOEA/D-MUP performs best on these test problems.For F8 with the convex PF, MOEA/D-MUP is only beaten by MOEA/D-DRA-UT, while MOEA/D-MUP is obviously better than MOEA/D-M2M.As shown in Figure 8, MOEA/D-MUP obtains a better approximation of the PF on F8 than MOEA/D-M2M.For example, boundary solutions obtained by MOEA/D-DRA-UT have not completely converged to the true PF, while MOEA/D-MUP performs much better.
As shown in Figure 9, the final nondominated sets obtained by the three algorithms on F9, F12 and F16 are plotted.It is obvious that MOEA/D-MUP performs significantly better than MOEA/D-DRA-UT and MOEA/D-M2M for these three-objective problems, i.e., F9, F12 and F16.Moreover, it can be seen from Figure 9 that more boundary solutions are obtained by MOEA/D-MUP, while the other two algorithms cannot approximate the PF boundaries sufficiently.
. .Parameter Sensitivity Analysis.In our algorithm, the maximal number of solutions replaced by a newly generated solution could be set as large as T (the neighborhood size).However, in most MOEA/Ds [19], the maximal number of solutions replaced by a newly generated solution is bounded by nr, which is set to be much smaller than T. To study the impact of nr in MOEA/D-MUP, different values (i.e., 2, 10, 20 and 30) are used for performance comparison.Due to page limitations, only the curves of IGD values obtained by four nr values on some test problems, i.e., F6 and F16, are provided in Figures 10(a) and 10(b).As shown in Figure 10, it is found that the performance of MOEA/D-MUP with a large nr value is superior to that with a small nr value.Thus, in MOEA/D-MUP, the maximal number of solutions replaced by a newly generated solution is set to be as large as the neighborhood size, which can better balance convergence and diversity during the evolutionary process.

Conclusion
In this paper, we have proposed an evolutionary search method with multiple utopian reference points for MOEA/D, called MOEA/D-MUP, which uses multiple utopian reference points to provide evolutionary search directions for different weight vectors.In order to effectively exploit the undiscovered boundaries of PF and enhance the breadthdiversity, the utopian reference point of each search direction is got by transforming the ideal point at a specified distance.Moreover, the nadir point is used to restrict the maximal translational distance.In other words, the positions of multiple utopian reference points in the objective space are bounded by the ideal point and the nadir point.Comparing our proposed algorithm with eight state-of-the-art MOEAs in terms of both IGD and I H metrics, we have witnessed the effectiveness and competitiveness of MOEA/D-MUP on F1-F16 used in this paper.As indicated by the experimental results, most existing MOEA/Ds cannot approximate the entire PFs on these test problems because the boundaries of PF cannot be effectively exploited.For the problems with difficult-to-approximate PF boundaries, there is a little effort that has been made to (1) effectively exploit the boundary areas and (2) balance exploration and exploitation in the search space.Therefore, there is much room for further development.
In the future work, we would like to investigate more effective selection mechanisms to better balance convergence and diversity during the evolutionary search process when tackling the problems with difficult-to-approximate PF boundaries.Specifically, we will further study a more  accurate generation method for utopian reference points and design recombination operators with stronger search ability.We also will put more efforts to develop mechanisms that can adaptively approximate the complex PFs and to extend our algorithm for solving many-objective optimization problems.Moreover, it is valuable to apply our proposed algorithm to solve some real-world applications [56,57].was supported by the National Engineering Laboratory for Big Data System Computing Technology.

Figure 2 :
Figure 2: Four different cases of reference points in decomposition functions: (a) the ideal point; (b) the nadir point; (c) the utopian point; (d) multiple reference points.

Figure 3 :
Figure 3: The missing of boundary PF parts in the cases: (a) the ideal reference point: (b) the nadir reference point.

2 Figure 5 :
Figure 5: An illustration of correlation between a utopian point and its corresponding search direction for two objectives.

)Proposition 7 .
which is shown in Figure (a).en,   ( () ,   ) =            () −            1 .(17)In the above equation, since   () ≥    for all  ∈ {1, . . ., }, the geometric property of the subproblems' objective functions is the same with that of the classical TCH decomposition function.According to Proposition 1 Complexity introduced in Section 2, Proposition 6 can be tenable by replacing the ideal reference point with the utopian reference point.Let   = (  1 , . . .,    ) be a utopian reference point of ( ) and  = ( 1 , . . .,   )  be a positive weight vector.Given an objective vector () = ( 1 (), . . .,   ()), if () = ( 1 (), . . .,   ()) is generated to satisfy two constraints: ( ) () and () are located in the same contour line, i.e.,   (() | ,   ) =   (() | ,   ); ( ) () lies in the line  2 as shown in Figure (b); we have   ( () | ,   ) =            () −            1 , respectively.In addition, four representative MOEA/D variants such as MOEA/D-STM, MOEA/D-M2M, MOEA/D-AWA and MOEA/D-PaS are also chosen to further validate the superiority of our proposed algorithm.The values of IGD and I H obtained by the four MOEA/D variants are shown in Tables -5 that our proposed MOEA/D-MUP achieves the best results in terms of the mean IGD values in all compared algorithms on F1-F6 and F9-F16.In addition, compared with the other algorithms except MOEA/D-DRA-UT and MOEA/D-M2M, MOEA/D-MUP obtains the best values for both mean IGD metric and I H metric on F7.For F8, MOEA/D-MUP only performs slightly worse than MOEA/D-DRA-UT in terms of mean IGD metric, while they have no significant difference in terms of mean I H metric. Next, we specifically analyze and discuss the superiority and drawback of MOEA/D-MUP with four classical MOEAs and four representative MOEA/Ds, respectively.The comparison results of MOEA/D-MUP with four classical MOEAs in terms of mean IGD and I H values are presented in Tables2 and 3, respectively.The best metric values are highlighted in bold face and the metric values which are significantly similar to the metric value of MOEA/D-MUP are marked with italic font.
, the experimental results of MOEA/D-MUP with four representative MOEA/Ds (MOEA/D-STM, MOEA/D-M2M, MOEA/D-AWA, and MOEA/D-PaS) also demonstrate the superiority of MOEA/ D-MUP on these test problems.It is clear from the two tables that MOEA/D-MUP significantly outperforms MOEA/D-STM, MOEA/D-AWA, and MOEA/D-PaS on F1-F16 and

Table 2 : 1 \
Mean and standard deviation of IGD values obtained by four classical algorithms and MOEA/D-MUP.According to Wilcoxon's rank sum test with a significant level  = 0.05, +, − and = indicate that the corresponding MOEA is significantly better than, worse than or similar to MOEA/D-MUP.

Table 3 :
Mean and standard deviation of I H values obtained by four classical algorithms and MOEA/D-MUP.on most test problems than MOEA/D-M2M.As shown in Tables ) and 7(b), although MOEA/D-DRA-UT shows a faster convergence speed than MOEA/D-MUP in the early stage, MOEA/D-MUP can finally obtain the best IGD values on F1 and F6.For F10, F13, and F16, respectively, in Figures 7(c), 7(d), and 7(e), MOEA/D-MUP also shows a much faster convergence speed during the evolution and achieves the best IGD values.