# A Kriging Model-Based Expensive Multiobjective Optimization Algorithm Using R2 Indicator of Expectation Improvement

Authors: Ding Han; Jianrong Zheng
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Most of the multiobjective optimization problems in engineering involve the evaluation of expensive objectives and constraint functions, for which an approximate model-based multiobjective optimization algorithm is usually employed, but requires a large amount of function evaluation. Aiming at effectively reducing the computation cost, a novel infilling point criterion EIR2 is proposed, whose basic idea is mapping a point in objective space into a set in expectation improvement space and utilizing the R2 indicator of the set to quantify the fitness of the point being selected as an infilling point. This criterion has an analytic form regardless of the number of objectives and demands lower calculation resources. Combining the Kriging model, optimal Latin hypercube sampling, and particle swarm optimization, an algorithm, EIR2-MOEA, is developed for solving expensive multiobjective optimization problems and applied to three sets of standard test functions of varying difficulty and comparing with two other competitive infill point criteria. Results show that EIR2 has higher resource utilization efficiency, and the resulting nondominated solution set possesses good convergence and diversity. By coupling with the average probability of feasibility, the EIR2 criterion is capable of dealing with expensive constrained multiobjective optimization problems and its efficiency is successfully validated in the optimal design of energy storage flywheel.


---

## Body

## 1. Introduction

Science and engineering practice possess a large number of multiobjective optimization problems (MOP) [1] whose objectives often conflict with each other and need to be optimized simultaneously, where the Pareto set is desired. Evolutionary Algorithm (EA) [2] has proven to be very suitable for MOP, some typical algorithms include NSGA-II [3], SPEA2 [4], MOEA/D [5], and MOPSO [6]. However, the success of these algorithms depends heavily on a large number of function evaluations and can only be applied in cases that function evaluation is cheap computationally. To handle some expensive MOPs, it is a desire to develop some efficient algorithms.To reduce the number of expensive evaluation required, approximate/surrogate model [7], which mimic the input-output relationship of expensive function but is cheap-to-evaluate, is often employed. Polynomial response surfaces [8], support vector machines [9], neural networks [10], and Kriging [11] are commonly used approximate models.Currently, there are two ways to integrate the approximate model with EA. The first one is to first construct a static global approximate model based on all sample points already evaluated, and the EA algorithm then searches for the optimal solution. However, with limited computational resources, it is hard to guarantee the approximate model’s accuracy over the entire design space. The second one utilizes an approximate model in a dynamic manner by first constructing a rough approximate model with a fraction of sample points generated by some Design of Experiments (DoE) [12] techniques, and then maximizing some infilling point criterion (IPC) [13], a new point is located by EA and is evaluated by expensive function. This process is repeated until resources are exhausted. The second manner has the advantage of high efficiency in using computation resources and is used in this paper.The infilling point criterion [14] plays a critical role in the approximate model-based optimization method. Expectation Improvement (EI), proposed in Efficient Global Optimization (EGO) [15] by Jone, may be one of the most researched method in the literature. It uses the Kriging model as an approximate model and considers not only the predicted value but also modeling uncertainty, achieving a good balance between exploration and exploitation in the optimization process. Other IPCs include statistical lower bound [16] and probability of improvement [17].However, those criterions are proposed originally only for single optimization and have to be modified for MOPs. Knowles first introduced EI into the field of multiobjective optimization and presented ParEGO [18], where the Kriging model is used and multiple objectives are converted into a single objective through a parameterized scalarizing weight vector, and then followed by the direct application of EGO, the weight vector changed randomly as iteration proceeded. Joan pointed out ParEGO is apt to favor model accuracy rather than finding nondominated solutions and demands more iterations before convergence and hence proposed Pareto Expected Improvement index(PEI) [19], which considers both EI and the probability of being Pareto optimal solution. Multi-EGO [20] builds the Kriging model and calculates separately EI for each objective and calls MOGA [21] as a solver to output a set of candidate points, from which a point is selected as an infilling point. MOEA/D-EGO [22] borrowed the ideas from ParEGO, and it generated not a single but a set of weight vectors that distributed uniformly in the objective space to form single objective optimization subproblems involving EI, optimizing jointly those subproblems produces multiple infill points and in this way parallelly infilling is achieved. Recently, some IPCs compatible with multiobjective optimization were proposed, such as MaxMin improvement [23] and expectation hypervolume improvement [24]. Svenson defined MaxMin improvement in the analytic form when the number of the objective function is two and suggested using the Monte Carlo method (MCM) [25] to calculate MaxMin value when the number of the objective is more than two.In this paper, a novel IPC is proposed for expensive multiobjective optimization. The outstanding feature of the IPC is to map one point in the objective space to a set in the expectation improvement space, where the fitness of the point is defined as a Quality indicator (QI) [26] value of the set. In this way, QI plays a role in guiding optimization towards an infilling point, rather than evaluating the nondominated solution set after optimization. There are a large number of candidate QIs in the literature, such as IGD [27, 28], ε-indicator [29], SDE indicator [30], and Hypervolume (HV) [31]. However, these QIs do not meet the requirements of easy-using and low-computation cost to some extent. For example, both IGD and ϵ-indicator demand Pareto front (PF) as the reference set, which is impossible in practical problems; HV does not require PF but a reference point, for which, if set unreasonably, will mislead the optimization search direction [32]. Besides, HV possesses high computational complexity, which increases exponentially with the number of objectives. Some QIs have a specific bias, for example, ε-indicator is more convergence-oriented, and as a result, optimization process may suffer from premature convergence, while SDE prefers diversity [33] which may make the optimization period too long [34]. R2 indicator [35, 36], having some desire properties in common with HV but is cheaper to evaluate, attracted widespread interest in recent years. In addition, it is weakly monotonic [37] and does not require PF to be provided; hence, it is adopted in this article.The rest of paper is organized as follows. In Section2, the mathematical background of Kriging, EI, and R2 indicator is reviewed. Then, the proposed EIR2 indicator and framework of the whole optimization algorithm are presented in Section 3. To evaluate the performance of the proposed IPC, a series of experiments over three benchmark test function suites against two competitors IPCs is conducted in Section 4, their results are measured by three performance metrics and then analyzed in Section 5, and followed by Section 6, where the proposed EIR2 is modified and applied to a real-world constrained expensive MOP, an optimization design of an energy storage flywheel. Lastly, Section 7 gives the conclusion and future research direction.


## 2. Background

### 2.1. Multiobjective Optimization Problem

Without loss of generality, a multiobjective optimization problem is formulated as follows:(1)minFx=f1x,⋯,fmxTs.t.gix≤0,i=1,⋯,phjx=0,j=1,⋯,qxmin≤x≤xmax,where Fx represents the objective function vector which contains m objective functions, gi and hj stand for inequality and equality constraint function with the total number p and q, respectively, and x=x1,⋯,xdT stands for design vector with lower bound xmin and upper bound xmax. A solution u is said to dominate another one v if and only if fiu≤fiv∀i=1,⋯,m and fiu<fiv∃i=1,⋯,m, written as u≺v. A solution x∗ is called Pareto optimal if there is no x satisfying x≺x∗(such x∗ is also referred to as nondominated by x). The set of all the Pareto optimal solutions is called the Pareto Set (PS), whose image in the objective space is called Pareto Front (PF).In practice, objective functions and constraint functions may come from finite element simulation; hence, they have no explicit analytic expressions and are usually computationally expensive to evaluate. Under this condition, we are interested in identifying a set of nondominated solutions, called approximated Pareto Set/Front, to represent the true Preto Set/Front, meanwhile consuming as less computing resources as possible.


### 2.2. Kriging Model and Expectation Improvement

Kriging model is an unbiased estimation model with smallest estimated variance, excellent high-dimensional, and nonlinear fitting capability. It assumes the relationship between inputx and predicted output y^ as follows:(2)y^x=μ+zx,where μ is a constant and zx∼N0,σ2 is the Gaussian random variable, which represent global trend and local deviation of the output separately.Having a set of known observationsS=xi,yi,i=1,⋯,N, our job is to search for the best values of μ and σ when predicting at a new point x. In the design space, the relationship between two different points xi and xj is characterized by correlation functions, such as the Gaussian function:(3)Rzxi,zxj;θ=∏k=1dexp−θkxki−xkj,where hyperparameter θ=θ1,⋯,θdT controls the smoothness of prediction function and is usually determined through maximizing a likelihood function of R by some intelligent optimization algorithms such as particle swarm optimization. According to optimization theory, analytical form of μ and σ are obtained:(4)μ^=1TR−1y1TR−11,σ^2=y−1μ^TR−1y−1μ^N,where y=y1,⋯,ydT is known as the output vector and R is a matrix of covariance between any two points in S, that is, Ri,j=Rzxi,zxj. Finally, at the new point x, we get the predicted value and variance:(5)y^x=μ^+r⊤R−1y−1μ^,s^2x=σ^21−r⊤Rr+1−1⊤R−1r21⊤R−11,where r is a correlation vector between point x and all points in S. Unlike other approximate models such as response surface [38] and neural network [39], which can only give out the predicted value, Kriging is capable of outputting predicted variance as a byproduct, which can be regarded as a predicting uncertainty at x.With Kriging, the value at a new pointx is treated as a Gaussian random variable, that is, yx∼Nμ^,s^2. An improvement quantity is designed as Iyx,ymin=maxyx−ymin,0, where ymin is the smallest y value in S. It is noted that only those yx that are less than ymin produce improvement. The expectation of improvement (EI) [40] has a closed form, that is,(6)EIyx,ymin=ymin−y^xΦymin−y^xs^x+s^ϕymin−y^xs^x,ifs^>0,0,ifs^=0,where Φ and ϕ are the probability distribution and probability density functions, respectively.


### 2.3. R2 Indicator

The R2 indicator [41] is a unitary indicator, which was proposed to evaluate the performance of the nondominated solution set Q, given as(7)R2Q,U=−1U∑u∈Umaxq∈Quq,where U is a discrete and finite set of utility functions u, for which there are many options to choose, such as weight sums and penalty-based boundary intersection [42].


## 2.1. Multiobjective Optimization Problem

Without loss of generality, a multiobjective optimization problem is formulated as follows:(1)minFx=f1x,⋯,fmxTs.t.gix≤0,i=1,⋯,phjx=0,j=1,⋯,qxmin≤x≤xmax,where Fx represents the objective function vector which contains m objective functions, gi and hj stand for inequality and equality constraint function with the total number p and q, respectively, and x=x1,⋯,xdT stands for design vector with lower bound xmin and upper bound xmax. A solution u is said to dominate another one v if and only if fiu≤fiv∀i=1,⋯,m and fiu<fiv∃i=1,⋯,m, written as u≺v. A solution x∗ is called Pareto optimal if there is no x satisfying x≺x∗(such x∗ is also referred to as nondominated by x). The set of all the Pareto optimal solutions is called the Pareto Set (PS), whose image in the objective space is called Pareto Front (PF).In practice, objective functions and constraint functions may come from finite element simulation; hence, they have no explicit analytic expressions and are usually computationally expensive to evaluate. Under this condition, we are interested in identifying a set of nondominated solutions, called approximated Pareto Set/Front, to represent the true Preto Set/Front, meanwhile consuming as less computing resources as possible.


## 2.2. Kriging Model and Expectation Improvement

Kriging model is an unbiased estimation model with smallest estimated variance, excellent high-dimensional, and nonlinear fitting capability. It assumes the relationship between inputx and predicted output y^ as follows:(2)y^x=μ+zx,where μ is a constant and zx∼N0,σ2 is the Gaussian random variable, which represent global trend and local deviation of the output separately.Having a set of known observationsS=xi,yi,i=1,⋯,N, our job is to search for the best values of μ and σ when predicting at a new point x. In the design space, the relationship between two different points xi and xj is characterized by correlation functions, such as the Gaussian function:(3)Rzxi,zxj;θ=∏k=1dexp−θkxki−xkj,where hyperparameter θ=θ1,⋯,θdT controls the smoothness of prediction function and is usually determined through maximizing a likelihood function of R by some intelligent optimization algorithms such as particle swarm optimization. According to optimization theory, analytical form of μ and σ are obtained:(4)μ^=1TR−1y1TR−11,σ^2=y−1μ^TR−1y−1μ^N,where y=y1,⋯,ydT is known as the output vector and R is a matrix of covariance between any two points in S, that is, Ri,j=Rzxi,zxj. Finally, at the new point x, we get the predicted value and variance:(5)y^x=μ^+r⊤R−1y−1μ^,s^2x=σ^21−r⊤Rr+1−1⊤R−1r21⊤R−11,where r is a correlation vector between point x and all points in S. Unlike other approximate models such as response surface [38] and neural network [39], which can only give out the predicted value, Kriging is capable of outputting predicted variance as a byproduct, which can be regarded as a predicting uncertainty at x.With Kriging, the value at a new pointx is treated as a Gaussian random variable, that is, yx∼Nμ^,s^2. An improvement quantity is designed as Iyx,ymin=maxyx−ymin,0, where ymin is the smallest y value in S. It is noted that only those yx that are less than ymin produce improvement. The expectation of improvement (EI) [40] has a closed form, that is,(6)EIyx,ymin=ymin−y^xΦymin−y^xs^x+s^ϕymin−y^xs^x,ifs^>0,0,ifs^=0,where Φ and ϕ are the probability distribution and probability density functions, respectively.


## 2.3. R2 Indicator

The R2 indicator [41] is a unitary indicator, which was proposed to evaluate the performance of the nondominated solution set Q, given as(7)R2Q,U=−1U∑u∈Umaxq∈Quq,where U is a discrete and finite set of utility functions u, for which there are many options to choose, such as weight sums and penalty-based boundary intersection [42].


## 3. Proposed Method

### 3.1. Expectation Improvement R2 Indicator

The main difficulty in extending EI infilling criterion to MOPs lies in the fact that we have to specifyymin before calculating the EI value, which is impractical as no such solution exists as being optimal with respect to every objective, but a set of nondominated solution found during previous optimization process. It is now highly desired to find an infilling point having large EI values with respect to all points in the current nondominated set. This motivates us to treat each nondominated solution as current optima and then calculate the EI value in each objective, resulting in a vector function:(8)Hx,p=EI1x,p,⋯,EImx,pT,where p=p1,⋯,pmT is a point in current nondominated set P and EIjx,p=EIyjx,pj. By this function, we obtain a vector consisting of all EI values for which pj is deemed as ymin in objective j. Obviously, we need to build m Kriging models, each corresponding to one objective function. After continually applying function H to all nondominated solutions, we get a set composed of P points, defined as(9)Qx=qq=Hx,p,p∈P,where q=q1,⋯,qmT. In essence, we mapped point x in design space to a set Qx in another space, called expectation improvement space, where the R2 indicator of Qx is defined as fitness of x being selected as the infilling point. The mapping from objective space to the EI space is illustrated in Figure 1.Figure 1
(a) A point with its predicted value and standard deviation in 2D objective space and (b) its corresponding mapping set in 2D EI space.
(a)(b)Combining EI and R2 indicator, a novel indicator for MOP, Expectation Improvement R2 Indicator (EIR2), is defined as follows:(10)EIR2x,P,U=−1U∑u∈Umaxq∈Qxuq.As EIR2 is a scalar indicator, some single optimization algorithm can be employed to determine the optimal solution as the infilling point. Besides, EIR2 actually defines an indicator template in which other utility functions can be specified. In this paper, Tchebyche function is employed as utility function because it can tackle MOPs with convex or concave Pareto front, defined as(11)uλq=−maxj∈1,⋯,mλjzj∗−qj,where λ=λ1,⋯,λmT is weight vector satisfying λj≥0 and ∑j=1mλj=1. These weight vectors are uniformly distributed in the objective space and are collected in the weight set Λ. z∗=z1∗,⋯,zm∗T is the ideal point, that is, zj∗=minqj,∀q∈Q. As all EI values are nonnegative, a natural selection is z∗=0. To summarize, the EIR2 indicator has the form(12)EIR2x,P,Λ=1Λ∑λ∈Λminp∈Pmaxj=1,⋯,mλj×EIjx,p.The point that maximizes the EIR2 indicator is recognized as an infilling point for expensive evaluation.


### 3.2. Weight Vectors Generating

The utility function is actually an optimization subproblem parameterized by the weight vector, by which the optimization process can be controlled. If the geometry characteristics of PF are known, such as being convex/concave, a set of weight vectors can be distributed consistent with the PF in the objective space, by which a more accurate approximate PF with uniform distribution is more likely to be obtained. However, there is often no information about the PF in priori, which is common for practical MOP; hence, a better choice in this situation is to treat all objectives equally important and generate a set of weight vectors evenly distributed in the objective space.To help generate weight vectorλ, an auxiliary set is introduced L=0/H,1/H,⋯,H/H, where H is an user-specified positive integer. We can chose randomly m−1 element from L as the first m−1 elements of λ and assign the value of λm to ensure ∑i=1mλi=1. Exhausting every possible combination, there are totally CH+m−1m−1 weight vectors generated. Figure 2 illustrates all 21 weight vectors generated using parameter m = 3 and H = 5.Figure 2
Evenly distributed weight vectors in objective space withm=3 and H=5.


### 3.3. Comparison with Other IPCs

The novelty of EIR2 is that fitness assessment of point is carried out not in the objective space but in the EI space, where the R2 indicator is calculated for the mapping set of the point. Since the EI indicator has a closed function form, so does EIR2 regardless of the number of objective functions. In addition, the computational complexity of EIR2 isOm∗P∗Λ, which is considerably less than that of hypervolume [26].ParEGO also uses the Chebyshev function; however, its purpose is to integrate multiple objectives into one single objective before applying the EGO algorithm. Therefore, in essence, the infilling point criterion used in ParEGO is EI.With regard to the MaxMin indicator, despite analytic form being existing, its calculation has to resort to dividing the objective space into many blocks on which tedious multivariable integration is carried out. Even worse, the number of blocks increases exponentially with the growth of the number of objective. So, in practice, we only apply this analytic form whenm=2 while using the Monte Carlo method when m>2.Figure3 shows contour plots of three IPCs, namely, ParEGO, MaxMin, and EIR2 applied in test function TEST2 defined in Table 1, where lighter color denotes higher value. In these plots, 20 known sample points are represented by red points, and blue points represent the true Pareto set, which are composed of three disconnected point subsets. It is easy to find that ParEGO can identify only one Pareto subset out of three, while MaxMin and EIR2 successfully locate the three Pareto subsets. In the region, where no Pareto solution exists, the function value of EIR2 is almost constant, while the function landscape of MaxMin shows a ladder-like trend.Figure 3
The landscape of infilling point criteria in 2D decision space for TEST2 by (a) ParEGO with weight vector [0.8, 0.2], (b) MaxMin, and (c) EIR2 with 20 uniformly distributed weight vectors. (d) Both the true Pareto set and 20 sample points generated by the optimal Latin hypercube sampling method.
(a)(b)(c)(d)Table 1
The definition of test functions.
NameDimensionBoundaryObjective functionComments on PFTEST12x1∈0.4,1.6x2∈2,5f1=x1−22+x2−12f2=x12+x2−62Convex, simpleTEST22x1∈0,1x2∈−2,2f1=x1f2=1+x22−x1−0.1sin5πx1Convex, disconnectedFON3xi∈−4,4f1=1−exp−∑i=1nxi−1/32f2=1−exp−∑i=1nxi+1/32NonconvexZDT15xi∈0,1f1=x1f2=z1−x1/zz=1+9/n−1∑i=2nxiConvexZDT25xi∈0,1f1=x1f2=z1−x1/z2z=1+9/n−1∑i=2nxiNonconvexZDT35xi∈0,1f1=x1f2=z1−x1/z−x1/zsin10πx1z=1+9/n−1∑i=2nxiConvex, disconnectedDTLZ25xi∈0,1f1=1+gcosπx1/2cosπx2/2f2=1+gcosπx1/2sinπx2/2f3=1+gsinπx1/2g=∑i=3,4,5xi−0.52NonconvexDTLZ55xi∈0,1f1=1+gcosπθ1/2cosπθ2/2f2=1+gcosπθ1/2sinπθ2/2f3=1+gsinπθ1/2θi=π/41+g1+2gxi,i=1,2g=∑i=3,4.5xi−0.52ConvexDTLZ75xi∈0,1f1=x1f2=x2f3=1+ghg=1+3∑i=3,4,5xih=3−∑i=12fi/1+g1+sin3πfiDisconnectedDeutz [43] recently proposed an IPC, also called EIR2, which is however very different from ours, proposed in this article. The main difference is that the R2 indicator is calculated in different spaces. Deutz’s R2 indicator is defined in objective space and the expectation of R2 indicators is carried out through multivariate integration. Therefore, the analytic form can only be obtained when m=2 but does not exist when m≥3, where the Monte Carlo method has to be adopted. In contrast, EIR2 in this paper calculates the R2 indicator in the EI space and hence having the analytic form for any m. Most importantly, our EIR2 involves no multivariable integration and hence is cheap to evaluate.


### 3.4. Overview of the EIR2-EMOA

Combining the devised EIR2 infilling point criterion with the Kriging model, optimal Latin hypercube sampling method [44], and particle swarm optimization [45], an optimization algorithm for multiobjective problems with expensive functions is proposed, referred to as EIR2-MOEA. The detailed procedure is shown in Algorithm 1.Algorithm 1:EIR2-MOEA.
Input:Nmax: Maximum number of evaluations by expensive functionα: Ratio of the number of initial samples to TmaxH: Integer used to generate weight vectorOutput:P: Nondominated solution set(1) Initialization: Obtain design variable limitsxmax and xmin, number of objective m, etc. according to MOP to be solved. Set S=∅ and P=∅.(2) GenerateαNmax sample points in design space xmin,xmax using optimal Latin hypercube sampling.(3)fork = 1 to αNmaxdo(4) Calculate the expensive objective function valuesyk=Fxk for sample point xk(5)S=S∪xk,yk, i.e., add sample point to S(6)end for(7) Generatec=CH+m−1m−1 weight vectors in objective space and save them in set Λ=λ1,⋯,λc.(8) Find non-dominated solutions inS and put them into P(9)fork=1to1−αNmaxdo(10)fork=1tomdo(11) Using PSO to find the best hyperparameterθ(12) Build Kriging model of them-th objective function based on S(13)end for(14) According to EIR2 indicator, apply PSO to find the best infilling pointxinf(15) Calculate expensive objective functionyinf=Fxinf(16) UpdateS=S∪xinf,yinf(17) Find the non-dominated solutions inS and put them into P(18)end for return Output P as approximated Pareto setIn the initial stage of EIR2-MOEA, a set of sample points of small size need to be selected for establishing the initial Kriging model. As a design of experiment method, the optimal Latin Hypercube Sampling (LHS) method is utilized as it has flexibility in setting the number of sampling points and can ensure these points filling in the entire design space.When maximizing the EIR2 value and looking for optimal hyperparameters of the Kriging model, an optimization method is required. However, the landscapes of these two functions are nonlinear, for which traditional gradient-based optimization algorithms are prone to falling into local minima if the initial point is not selected wisely and has to be rerun multiple times. Particle swarm optimization algorithm, as an intelligent optimization algorithm, does not need information about gradient and has fewer parameter settings, high efficiency of search, and lower complexity, hence is employed to find the optimal parameters.


## 3.1. Expectation Improvement R2 Indicator

The main difficulty in extending EI infilling criterion to MOPs lies in the fact that we have to specifyymin before calculating the EI value, which is impractical as no such solution exists as being optimal with respect to every objective, but a set of nondominated solution found during previous optimization process. It is now highly desired to find an infilling point having large EI values with respect to all points in the current nondominated set. This motivates us to treat each nondominated solution as current optima and then calculate the EI value in each objective, resulting in a vector function:(8)Hx,p=EI1x,p,⋯,EImx,pT,where p=p1,⋯,pmT is a point in current nondominated set P and EIjx,p=EIyjx,pj. By this function, we obtain a vector consisting of all EI values for which pj is deemed as ymin in objective j. Obviously, we need to build m Kriging models, each corresponding to one objective function. After continually applying function H to all nondominated solutions, we get a set composed of P points, defined as(9)Qx=qq=Hx,p,p∈P,where q=q1,⋯,qmT. In essence, we mapped point x in design space to a set Qx in another space, called expectation improvement space, where the R2 indicator of Qx is defined as fitness of x being selected as the infilling point. The mapping from objective space to the EI space is illustrated in Figure 1.Figure 1
(a) A point with its predicted value and standard deviation in 2D objective space and (b) its corresponding mapping set in 2D EI space.
(a)(b)Combining EI and R2 indicator, a novel indicator for MOP, Expectation Improvement R2 Indicator (EIR2), is defined as follows:(10)EIR2x,P,U=−1U∑u∈Umaxq∈Qxuq.As EIR2 is a scalar indicator, some single optimization algorithm can be employed to determine the optimal solution as the infilling point. Besides, EIR2 actually defines an indicator template in which other utility functions can be specified. In this paper, Tchebyche function is employed as utility function because it can tackle MOPs with convex or concave Pareto front, defined as(11)uλq=−maxj∈1,⋯,mλjzj∗−qj,where λ=λ1,⋯,λmT is weight vector satisfying λj≥0 and ∑j=1mλj=1. These weight vectors are uniformly distributed in the objective space and are collected in the weight set Λ. z∗=z1∗,⋯,zm∗T is the ideal point, that is, zj∗=minqj,∀q∈Q. As all EI values are nonnegative, a natural selection is z∗=0. To summarize, the EIR2 indicator has the form(12)EIR2x,P,Λ=1Λ∑λ∈Λminp∈Pmaxj=1,⋯,mλj×EIjx,p.The point that maximizes the EIR2 indicator is recognized as an infilling point for expensive evaluation.


## 3.2. Weight Vectors Generating

The utility function is actually an optimization subproblem parameterized by the weight vector, by which the optimization process can be controlled. If the geometry characteristics of PF are known, such as being convex/concave, a set of weight vectors can be distributed consistent with the PF in the objective space, by which a more accurate approximate PF with uniform distribution is more likely to be obtained. However, there is often no information about the PF in priori, which is common for practical MOP; hence, a better choice in this situation is to treat all objectives equally important and generate a set of weight vectors evenly distributed in the objective space.To help generate weight vectorλ, an auxiliary set is introduced L=0/H,1/H,⋯,H/H, where H is an user-specified positive integer. We can chose randomly m−1 element from L as the first m−1 elements of λ and assign the value of λm to ensure ∑i=1mλi=1. Exhausting every possible combination, there are totally CH+m−1m−1 weight vectors generated. Figure 2 illustrates all 21 weight vectors generated using parameter m = 3 and H = 5.Figure 2
Evenly distributed weight vectors in objective space withm=3 and H=5.


## 3.3. Comparison with Other IPCs

The novelty of EIR2 is that fitness assessment of point is carried out not in the objective space but in the EI space, where the R2 indicator is calculated for the mapping set of the point. Since the EI indicator has a closed function form, so does EIR2 regardless of the number of objective functions. In addition, the computational complexity of EIR2 isOm∗P∗Λ, which is considerably less than that of hypervolume [26].ParEGO also uses the Chebyshev function; however, its purpose is to integrate multiple objectives into one single objective before applying the EGO algorithm. Therefore, in essence, the infilling point criterion used in ParEGO is EI.With regard to the MaxMin indicator, despite analytic form being existing, its calculation has to resort to dividing the objective space into many blocks on which tedious multivariable integration is carried out. Even worse, the number of blocks increases exponentially with the growth of the number of objective. So, in practice, we only apply this analytic form whenm=2 while using the Monte Carlo method when m>2.Figure3 shows contour plots of three IPCs, namely, ParEGO, MaxMin, and EIR2 applied in test function TEST2 defined in Table 1, where lighter color denotes higher value. In these plots, 20 known sample points are represented by red points, and blue points represent the true Pareto set, which are composed of three disconnected point subsets. It is easy to find that ParEGO can identify only one Pareto subset out of three, while MaxMin and EIR2 successfully locate the three Pareto subsets. In the region, where no Pareto solution exists, the function value of EIR2 is almost constant, while the function landscape of MaxMin shows a ladder-like trend.Figure 3
The landscape of infilling point criteria in 2D decision space for TEST2 by (a) ParEGO with weight vector [0.8, 0.2], (b) MaxMin, and (c) EIR2 with 20 uniformly distributed weight vectors. (d) Both the true Pareto set and 20 sample points generated by the optimal Latin hypercube sampling method.
(a)(b)(c)(d)Table 1
The definition of test functions.
NameDimensionBoundaryObjective functionComments on PFTEST12x1∈0.4,1.6x2∈2,5f1=x1−22+x2−12f2=x12+x2−62Convex, simpleTEST22x1∈0,1x2∈−2,2f1=x1f2=1+x22−x1−0.1sin5πx1Convex, disconnectedFON3xi∈−4,4f1=1−exp−∑i=1nxi−1/32f2=1−exp−∑i=1nxi+1/32NonconvexZDT15xi∈0,1f1=x1f2=z1−x1/zz=1+9/n−1∑i=2nxiConvexZDT25xi∈0,1f1=x1f2=z1−x1/z2z=1+9/n−1∑i=2nxiNonconvexZDT35xi∈0,1f1=x1f2=z1−x1/z−x1/zsin10πx1z=1+9/n−1∑i=2nxiConvex, disconnectedDTLZ25xi∈0,1f1=1+gcosπx1/2cosπx2/2f2=1+gcosπx1/2sinπx2/2f3=1+gsinπx1/2g=∑i=3,4,5xi−0.52NonconvexDTLZ55xi∈0,1f1=1+gcosπθ1/2cosπθ2/2f2=1+gcosπθ1/2sinπθ2/2f3=1+gsinπθ1/2θi=π/41+g1+2gxi,i=1,2g=∑i=3,4.5xi−0.52ConvexDTLZ75xi∈0,1f1=x1f2=x2f3=1+ghg=1+3∑i=3,4,5xih=3−∑i=12fi/1+g1+sin3πfiDisconnectedDeutz [43] recently proposed an IPC, also called EIR2, which is however very different from ours, proposed in this article. The main difference is that the R2 indicator is calculated in different spaces. Deutz’s R2 indicator is defined in objective space and the expectation of R2 indicators is carried out through multivariate integration. Therefore, the analytic form can only be obtained when m=2 but does not exist when m≥3, where the Monte Carlo method has to be adopted. In contrast, EIR2 in this paper calculates the R2 indicator in the EI space and hence having the analytic form for any m. Most importantly, our EIR2 involves no multivariable integration and hence is cheap to evaluate.


## 3.4. Overview of the EIR2-EMOA

Combining the devised EIR2 infilling point criterion with the Kriging model, optimal Latin hypercube sampling method [44], and particle swarm optimization [45], an optimization algorithm for multiobjective problems with expensive functions is proposed, referred to as EIR2-MOEA. The detailed procedure is shown in Algorithm 1.Algorithm 1:EIR2-MOEA.
Input:Nmax: Maximum number of evaluations by expensive functionα: Ratio of the number of initial samples to TmaxH: Integer used to generate weight vectorOutput:P: Nondominated solution set(1) Initialization: Obtain design variable limitsxmax and xmin, number of objective m, etc. according to MOP to be solved. Set S=∅ and P=∅.(2) GenerateαNmax sample points in design space xmin,xmax using optimal Latin hypercube sampling.(3)fork = 1 to αNmaxdo(4) Calculate the expensive objective function valuesyk=Fxk for sample point xk(5)S=S∪xk,yk, i.e., add sample point to S(6)end for(7) Generatec=CH+m−1m−1 weight vectors in objective space and save them in set Λ=λ1,⋯,λc.(8) Find non-dominated solutions inS and put them into P(9)fork=1to1−αNmaxdo(10)fork=1tomdo(11) Using PSO to find the best hyperparameterθ(12) Build Kriging model of them-th objective function based on S(13)end for(14) According to EIR2 indicator, apply PSO to find the best infilling pointxinf(15) Calculate expensive objective functionyinf=Fxinf(16) UpdateS=S∪xinf,yinf(17) Find the non-dominated solutions inS and put them into P(18)end for return Output P as approximated Pareto setIn the initial stage of EIR2-MOEA, a set of sample points of small size need to be selected for establishing the initial Kriging model. As a design of experiment method, the optimal Latin Hypercube Sampling (LHS) method is utilized as it has flexibility in setting the number of sampling points and can ensure these points filling in the entire design space.When maximizing the EIR2 value and looking for optimal hyperparameters of the Kriging model, an optimization method is required. However, the landscapes of these two functions are nonlinear, for which traditional gradient-based optimization algorithms are prone to falling into local minima if the initial point is not selected wisely and has to be rerun multiple times. Particle swarm optimization algorithm, as an intelligent optimization algorithm, does not need information about gradient and has fewer parameter settings, high efficiency of search, and lower complexity, hence is employed to find the optimal parameters.


## 4. Experiment Setup

In order to assess the performance of the proposed algorithm, we applied it to three sets of standard test function suites and compared it with another two algorithms, ParEGO and MaxMin.


### 4.1. Test Problems

Three benchmark function suites, namely, simple set, ZDT set, and DTLZ set, are selected to verify the performance of the IPCs in different aspects, whose definitions are shown in Table1. The simple set includes three biobjective optimization problems TEST1, TEST2, and FON with the number of variables less than three, whose corresponding PFs are convex, concave, and discontinuous, respectively. The ZDT set contains three biobjective test functions, ZDT1, ZDT2, and ZDT3, with convex, concave, and discontinuous PF, respectively. We set the number of variables to 5 to test IPCs’ performance on problems having more decision variables. DTLZ2, DTLZ5, and DTLZ7 form the DTLZ test set, which all have three objectives and five decision variables, and they are selected to assess the performance of the IPCs when the number of objectives is large.


### 4.2. Performance Metric

We hope that the approximated Pareto front found close enough to the true Pareto front and is evenly distributed in the objective space, that is, to meet the requirement on convergence and diversity. For this purpose, we choose three performance metrics to quantify the performance of the approximated Pareto front.(1)
Inverted Generational Distance (IGD) [28]:(13)IGDT,P=1T∑x∈Tminp∈Pdx,p,whereT and P are true and approximated Pareto set, respectively, T is the number of solutions in T, and dx,p represents the Euclidean distance between x and p in objective space. IGD measures the average distance between the true and approximated Pareto set. Therefore, the smaller IGD value is preferred.(2)
HyperVolume (HV) [31]:(14)HVP=Vol∪x∈Pf1x,r1×…fmx,rm,wherer=r1,⋯,rmT represents reference point which is dominated by all solutions in P and Vol stands for Lebesgue measure. Geometrically, HV is the volume enclosed by the approximated Pareto set P and the reference point r. HV can reflect both convergence and uniformity. If the approximate Pareto set performances well with respect to both aspects, then the HV value will be high.(3)
Nondominated solution ratio (NR):(15)NRP=PNmax,whereP is the number of nondominated solutions and Nmax the maximum number of evaluations by expensive function. NR is used to measure the proportion of the nondominated solutions out of the total sample points evaluated by expensive functions. In other words, NR represents the efficiency of IPC in using computing resources. The larger the NR value, the more nondominated solutions found by the IPC and the higher the utilization efficiency of computing resources.


### 4.3. Parameter Setting

We hope to evaluate the performance of each IPC with limited computing resources, that is, for each test function we set the maximum number of evaluations involving expensive functionsNmax, shown in Table 2. The size of the initial sample set is specified to be a fraction of Nmax, namely, αNmax, where α is a proportional coefficient which is 1/3 in this paper. To ensure fairness, all three IPCs start with the same initial sample point set generated by optimal LHS in each run, and ten runs are executed for each test function.Table 2
Maximum number of evaluations involving expensive functions.
TEST1TEST2FONZDT1ZDT2ZDT3DTLZ2DTLZ5DTLZ7Nmax4560909090120210210210Reference points are required for calculating the HV value. All biobjective optimization problems share the same reference point1.1,1.1 except TEST1 for which 20.5,20.5 is used; all three-objective optimization problems are assigned the same reference point 1.1,1.1,1.1 except DTLZ7, where 0.95,0.95,6.6 is utilized.As for particle swarm optimization, the population size is set to 100, the evolution generation is 100, and the scaling factor is 1.2. Other parameters of MaxMin and ParEGO are taken as the recommended values in original papers.


## 4.1. Test Problems

Three benchmark function suites, namely, simple set, ZDT set, and DTLZ set, are selected to verify the performance of the IPCs in different aspects, whose definitions are shown in Table1. The simple set includes three biobjective optimization problems TEST1, TEST2, and FON with the number of variables less than three, whose corresponding PFs are convex, concave, and discontinuous, respectively. The ZDT set contains three biobjective test functions, ZDT1, ZDT2, and ZDT3, with convex, concave, and discontinuous PF, respectively. We set the number of variables to 5 to test IPCs’ performance on problems having more decision variables. DTLZ2, DTLZ5, and DTLZ7 form the DTLZ test set, which all have three objectives and five decision variables, and they are selected to assess the performance of the IPCs when the number of objectives is large.


## 4.2. Performance Metric

We hope that the approximated Pareto front found close enough to the true Pareto front and is evenly distributed in the objective space, that is, to meet the requirement on convergence and diversity. For this purpose, we choose three performance metrics to quantify the performance of the approximated Pareto front.(1)
Inverted Generational Distance (IGD) [28]:(13)IGDT,P=1T∑x∈Tminp∈Pdx,p,whereT and P are true and approximated Pareto set, respectively, T is the number of solutions in T, and dx,p represents the Euclidean distance between x and p in objective space. IGD measures the average distance between the true and approximated Pareto set. Therefore, the smaller IGD value is preferred.(2)
HyperVolume (HV) [31]:(14)HVP=Vol∪x∈Pf1x,r1×…fmx,rm,wherer=r1,⋯,rmT represents reference point which is dominated by all solutions in P and Vol stands for Lebesgue measure. Geometrically, HV is the volume enclosed by the approximated Pareto set P and the reference point r. HV can reflect both convergence and uniformity. If the approximate Pareto set performances well with respect to both aspects, then the HV value will be high.(3)
Nondominated solution ratio (NR):(15)NRP=PNmax,whereP is the number of nondominated solutions and Nmax the maximum number of evaluations by expensive function. NR is used to measure the proportion of the nondominated solutions out of the total sample points evaluated by expensive functions. In other words, NR represents the efficiency of IPC in using computing resources. The larger the NR value, the more nondominated solutions found by the IPC and the higher the utilization efficiency of computing resources.


## 4.3. Parameter Setting

We hope to evaluate the performance of each IPC with limited computing resources, that is, for each test function we set the maximum number of evaluations involving expensive functionsNmax, shown in Table 2. The size of the initial sample set is specified to be a fraction of Nmax, namely, αNmax, where α is a proportional coefficient which is 1/3 in this paper. To ensure fairness, all three IPCs start with the same initial sample point set generated by optimal LHS in each run, and ten runs are executed for each test function.Table 2
Maximum number of evaluations involving expensive functions.
TEST1TEST2FONZDT1ZDT2ZDT3DTLZ2DTLZ5DTLZ7Nmax4560909090120210210210Reference points are required for calculating the HV value. All biobjective optimization problems share the same reference point1.1,1.1 except TEST1 for which 20.5,20.5 is used; all three-objective optimization problems are assigned the same reference point 1.1,1.1,1.1 except DTLZ7, where 0.95,0.95,6.6 is utilized.As for particle swarm optimization, the population size is set to 100, the evolution generation is 100, and the scaling factor is 1.2. Other parameters of MaxMin and ParEGO are taken as the recommended values in original papers.


## 5. Result and Discussion

For each test function, ten independent runs are performed, and the HV, IGD, and NR indicators of the resulting approximated Pareto set are calculated. Finally, their mean and standard deviation of the ten runs are calculated and listed, respectively, in Tables3–5. For each test function, among the three IPCs, the winning one is shown in boldface with a gray background.Table 3
The mean and standard deviation of the HV value of approximated Pareto set found by ParEGO, MaxMin, and EIR2-EMOA on all test functions.
ParEGOMaxMinEIR2MeanStdMeanStdMeanStdTEST1272.20241.71E + 00278.86261.03E − 01279.02288.64E − 02TEST20.69008.02E − 030.72271.49E − 030.72337.44E − 04FON0.45833.10E − 020.53242.21E − 030.53643.08E − 03ZDT10.84685.53E − 030.84635.01E − 030.86431.08E − 03ZDT20.50572.14E − 020.51214.63E − 030.53121.04E − 03ZDT31.29258.52E − 031.29464.72E − 021.32301.14E − 03DTLZ20.64999.69E − 030.70785.60E − 030.73782.85E − 03DTLZ50.40985.02E − 030.41513.90E − 030.43583.45E − 04DTLZ71.54102.31E − 021.63332.72E − 021.68583.39E − 03Table 4
The mean and standard deviation of the IGD value of approximated Pareto set found by ParEGO, MaxMin, and EIR2-EMOA on all test functions.
ParEGOMaxMinEIR2MeanStdMeanStdMeanStdTEST10.43474.88E − 020.19395.19E − 030.19084.49E − 03TEST20.04151.02E − 020.00975.06E − 040.00975.49E − 04FON0.06403.47E − 020.01211.05E − 030.01091.98E − 03ZDT10.02112.42E − 030.02324.17E − 030.00968.15E − 04ZDT20.02851.51E − 020.02493.26E − 030.01261.16E − 03ZDT30.05919.69E − 030.05441.21E − 020.01431.79E − 03DTLZ20.09995.40E − 030.09926.36E − 030.07004.21E − 03DTLZ50.02893.01E − 030.03344.89E − 030.00844.72E − 04DTLZ70.13931.70E − 020.08511.25E − 020.05764.23E − 03Table 5
The mean and standard deviation of the NR value of approximated Pareto set found by ParEGO, MaxMin, and EIR2-EMOA on all test functions.
ParEGOMaxMinEIR2MeanStdMeanStdMeanStdTEST10.66895.78E − 020.87112.52E − 020.86442.86E − 02TEST20.28832.95E − 020.61832.66E − 020.60172.54E − 02FON0.18674.50E − 020.42003.22E − 020.43677.33E − 02ZDT10.26111.76E − 020.25563.19E − 020.45223.31E − 02ZDT20.24786.87E − 020.20332.10E − 020.43113.77E − 02ZDT30.12422.17E − 020.12831.77E − 020.29582.95E − 02DTLZ20.27521.81E − 020.27332.49E − 020.43812.00E − 02DTLZ50.16521.27E − 020.12621.64E − 020.35291.09E − 02DTLZ70.19242.43E − 020.36866.87E − 020.54382.53E − 02


### 5.1. Simple Tests

TEST1, TEST2, and FON are all biobjective test functions but have, respectively, convex, piecewise continuous, and concave Pareto front. The comparison of mean and std value for HV, IGD, and NR metric can be found at the top part of Tables3–5, respectively. It is shown that ParEGO achieved the worst performance on all three performance metrics, while EIR2 and MinMax achieve similar results but both significantly better than ParEGO. Specifically, for TEST1 and TEST2, EIR2 is better than MaxMin in terms of HV and IGD but is slightly worse concerning NR. For FON, EIR2 ranked first on all three performance metrics.Figure4 plots the approximated Pareto front obtained by the three IPCs on TEST1, TEST2, and FON test function, respectively, with its HV value being median in ten runs. It is evident that the number of nondominated solutions obtained by ParEGO is small and they are not evenly distributed, which is more obvious on TEST2 and FON. This phenomenon can be foreseen because ParEGO lacks a mechanism to guarantee the uniform distribution of the solutions, that is, its weight vector is randomly generated at each iteration. In contrast, the approximated Pareto fronts obtained by EIR2 and MaxMin seem to have completely covered the true Pareto fronts and have good distribution performance, regardless of the nature of the true Pareto front of the test function.Figure 4
Approximated Pareto front found by three infilling criterions ParEGO(left), MaxMin(middle), and EIR2(right) on test function TEST1 (a∼c), TEST2 (d∼f), and FON (g∼i).
(a)(b)(c)(d)(e)(f)(g)(h)(i)Overall, EIR2 and MaxMin have comparable performance and are much better than ParEGO on simple test functions with a low number of variables.


### 5.2. ZDT Tests

The middle part of Tables3–5 presents the IGD, HV, and NR indicator comparison on ZDT test suites. Figure 5 illustrates the approximated Pareto fronts obtained by each IPC on ZDT1∼ZDT3 whose HV value being median in ten runs. It is clear that EIR2 is the best in terms of HV, IGD, and NR indicators for all ZDT test functions. Specifically, the HV value and IGD value of EIR2 are much better than that of ParEGO and MaxMin. It is worth noting that the advantages of EIR2 are quite obvious for IGD and NR indicators. For example, for ZDT3 EIR2 gets an NR value of 0.5438, while the other two algorithms only have 0.1924 and 0.3686; meanwhile, it gets an IGD value of 0.0143 while the other two 0.0544 and 0.0591.Figure 5
Approximated Pareto front found by three infilling criterions ParEGO(left), MaxMin(middle), and EIR2(right) on test function ZDT1 (a∼c), ZDT2 (d∼f), and ZDT3 (g∼i).
(a)(b)(c)(d)(e)(f)(g)(h)(i)According to Figure5, it is seen that both ParEGO and MaxMin successfully converged to the true Pareto front of ZDT1 and ZDT2; however, they failed in making them evenly distributed. This phenomenon is even worse on ZDT3 whose true Pareto front consists of five segmented curves. The majority of solutions found by ParEGO and MaxMin are clustered in one segment, while only a few solutions located in other segments. In contrast, the approximated Pareto set obtained by EIR2 are uniformly and equally distributed among all five segments.In summary, EIR2 is capable of finding more nondominated solutions than ParEGO and MaxMin using the same amount of computing resources and keeping them evenly distributed along the approximated Pareto front.


### 5.3. DTLZ Tests

The DTLZ test function suites are used aiming at testing an algorithm’s ability to three-objective optimization problems. The metric comparisons are listed at the bottom part of Tables3–5, respectively. Figure 6 shows each IPC’s results on DTLZ2, DTLZ5, and DTLZ7 with their HV value being median in ten runs. Same as the ZDT case, on the three test functions, the EIR2 is still significantly ahead of its two competitors in terms of all three performance metrics. To be more specific, take DTLZ2 as an instance, the NR value of EIR2 is 0.3529, which is nearly three times 0.1262 of MaxMin and more than twice 0.1652 of ParEGO, which shows again EIR2 is more efficient in using computation resources compared to other IPCs. The larger HV value and smaller IGD value obtained by EIR2 also indicate its capacity in achieving strong convergence ability.Figure 6
Approximated Pareto front found by three infilling criterions ParEGO(left), MaxMin(middle) and EIR2(right) on test function DTLZ1 (a∼c), DTLZ3 (d∼f), and DTLZ5 (g∼i).
(a)(b)(c)(d)(e)(f)(g)(h)(i)The same conclusion is also reached by inspecting Figure6. Specifically, ParEGO has difficulty in converging to the true Pareto front, while MaxMin, despite converged to PF, suffered from being not able to spread its approximated Pareto front uniformly distributed along the true Pareto front which is a 1/8 sphere of radius 1, and there are only a few solutions in the middle part and most of the solutions are found in three edges. In contrast, the approximated PF obtained by EIR2 is adequate to converge to and cover entirely and evenly the 1/8 spheres surface.


## 5.1. Simple Tests

TEST1, TEST2, and FON are all biobjective test functions but have, respectively, convex, piecewise continuous, and concave Pareto front. The comparison of mean and std value for HV, IGD, and NR metric can be found at the top part of Tables3–5, respectively. It is shown that ParEGO achieved the worst performance on all three performance metrics, while EIR2 and MinMax achieve similar results but both significantly better than ParEGO. Specifically, for TEST1 and TEST2, EIR2 is better than MaxMin in terms of HV and IGD but is slightly worse concerning NR. For FON, EIR2 ranked first on all three performance metrics.Figure4 plots the approximated Pareto front obtained by the three IPCs on TEST1, TEST2, and FON test function, respectively, with its HV value being median in ten runs. It is evident that the number of nondominated solutions obtained by ParEGO is small and they are not evenly distributed, which is more obvious on TEST2 and FON. This phenomenon can be foreseen because ParEGO lacks a mechanism to guarantee the uniform distribution of the solutions, that is, its weight vector is randomly generated at each iteration. In contrast, the approximated Pareto fronts obtained by EIR2 and MaxMin seem to have completely covered the true Pareto fronts and have good distribution performance, regardless of the nature of the true Pareto front of the test function.Figure 4
Approximated Pareto front found by three infilling criterions ParEGO(left), MaxMin(middle), and EIR2(right) on test function TEST1 (a∼c), TEST2 (d∼f), and FON (g∼i).
(a)(b)(c)(d)(e)(f)(g)(h)(i)Overall, EIR2 and MaxMin have comparable performance and are much better than ParEGO on simple test functions with a low number of variables.


## 5.2. ZDT Tests

The middle part of Tables3–5 presents the IGD, HV, and NR indicator comparison on ZDT test suites. Figure 5 illustrates the approximated Pareto fronts obtained by each IPC on ZDT1∼ZDT3 whose HV value being median in ten runs. It is clear that EIR2 is the best in terms of HV, IGD, and NR indicators for all ZDT test functions. Specifically, the HV value and IGD value of EIR2 are much better than that of ParEGO and MaxMin. It is worth noting that the advantages of EIR2 are quite obvious for IGD and NR indicators. For example, for ZDT3 EIR2 gets an NR value of 0.5438, while the other two algorithms only have 0.1924 and 0.3686; meanwhile, it gets an IGD value of 0.0143 while the other two 0.0544 and 0.0591.Figure 5
Approximated Pareto front found by three infilling criterions ParEGO(left), MaxMin(middle), and EIR2(right) on test function ZDT1 (a∼c), ZDT2 (d∼f), and ZDT3 (g∼i).
(a)(b)(c)(d)(e)(f)(g)(h)(i)According to Figure5, it is seen that both ParEGO and MaxMin successfully converged to the true Pareto front of ZDT1 and ZDT2; however, they failed in making them evenly distributed. This phenomenon is even worse on ZDT3 whose true Pareto front consists of five segmented curves. The majority of solutions found by ParEGO and MaxMin are clustered in one segment, while only a few solutions located in other segments. In contrast, the approximated Pareto set obtained by EIR2 are uniformly and equally distributed among all five segments.In summary, EIR2 is capable of finding more nondominated solutions than ParEGO and MaxMin using the same amount of computing resources and keeping them evenly distributed along the approximated Pareto front.


## 5.3. DTLZ Tests

The DTLZ test function suites are used aiming at testing an algorithm’s ability to three-objective optimization problems. The metric comparisons are listed at the bottom part of Tables3–5, respectively. Figure 6 shows each IPC’s results on DTLZ2, DTLZ5, and DTLZ7 with their HV value being median in ten runs. Same as the ZDT case, on the three test functions, the EIR2 is still significantly ahead of its two competitors in terms of all three performance metrics. To be more specific, take DTLZ2 as an instance, the NR value of EIR2 is 0.3529, which is nearly three times 0.1262 of MaxMin and more than twice 0.1652 of ParEGO, which shows again EIR2 is more efficient in using computation resources compared to other IPCs. The larger HV value and smaller IGD value obtained by EIR2 also indicate its capacity in achieving strong convergence ability.Figure 6
Approximated Pareto front found by three infilling criterions ParEGO(left), MaxMin(middle) and EIR2(right) on test function DTLZ1 (a∼c), DTLZ3 (d∼f), and DTLZ5 (g∼i).
(a)(b)(c)(d)(e)(f)(g)(h)(i)The same conclusion is also reached by inspecting Figure6. Specifically, ParEGO has difficulty in converging to the true Pareto front, while MaxMin, despite converged to PF, suffered from being not able to spread its approximated Pareto front uniformly distributed along the true Pareto front which is a 1/8 sphere of radius 1, and there are only a few solutions in the middle part and most of the solutions are found in three edges. In contrast, the approximated PF obtained by EIR2 is adequate to converge to and cover entirely and evenly the 1/8 spheres surface.


## 6. Engineering Application

In addition to boundary constraints, most MOPs also contain equality or inequality constraints. We assume that constraint functions are expensive; hence, the Kriging model need to be established for each constraint. On the purpose of measuring the degree to which the solution satisfies the constraints, probability of feasibility (PoF) [46, 47] is introduced as follows:(16)PoFix=Pgix≤0=Φ−g^ixs^ix,where g^ix and s^ix is the predicted value and standard deviation of constraint function gix. Based on PoF, the average probability of feasibility (APoF) is defined as follows:(17)APoF=1v∑i=1vPoFix,where v is the number of all inequality constraints. The APoF value of an infeasible solution is low, but it will increase as more constraints get satisfied. Clearly, the value of APoF varies in interval 0,1. The closer the value of APoF is to 1, the higher the probability that the solution is feasible, and hence more preferred.In combination withAPoF, EIR2 is extended to constrained MOPs(18)CEIR2x=EIR2x,P,Λ×APoFx=1Λ∑λ∈Λminp∈Pmaxj=1,⋯mλj×EIjx,p×1v∑i=1vΦ−g^ixs^ix.In this section, the proposed CEIR2 indicator is applied to the optimization design of energy storage flywheel [48] to test its effectiveness.The energy storage flywheel [49] uses a flywheel in high-speed rotating to store kinetic energy and convert it into electrical energy when needed. Given a constant speed, the maximum energy storage is proportional to the moment of inertia of flywheel. Larger moment of inertia can be obtained by adjusting the geometric parameters of the flywheel, but at the same time on the risk of increasing the total mass. In this paper, our job is to maximize the flywheel rotational inertia as well as minimize total mass.Figure7 is a cross-sectional profile of an energy storage flywheel. There are six decision variables: r1∼r4 and t1andt2, representing the radial sizes and thickness of the hub and spoke of the flywheel, respectively. Three inequality constraints involve restrictions on radial stress, radial deformation, and hoop stress. The overall optimization problem is expressed as follows:(19)minFx=−Ix,Mxs.t.g1x=σrmaxx−σr≤0g2x=σθmaxx−σθ≤0g3x=Δrmaxx−Δr≤0,where x=r1,r2,r3,r4,t1,t2, whose upper and lower boundary are listed in Table 6. Ix and Mx represent the rotational inertia and total mass of the flywheel, respectively. σrmax, σθmax, and Δrmax, respectively, represent the maximum radial stress, the maximum hoop stress, and the maximum radial deformation of the outer ring when the flywheel rotates, accompanied by their corresponding allowable values, σr, σθ, and Δr.Figure 7
Diagram of energy storage flywheel.Table 6
Range of decision variables.
Range (mm)x1x2x3x4x5x6xmax114.22128.19159.94203.1236.758.35xmin104.22118.19149.94193.1226.754.35Although formulas of the objective functionsIx and Mx can be obtained directly from the geometric parameters, the calculation of the constraint functions has to resort to finite element analysis, considering the complexity of flywheel structure and external force conditions; hence, this problem is a constrained expensive multiobjective optimization problem.The material used for the flywheel is Ti-6Al-4V titanium alloy, with corresponding allowable stressσr=σθ=800MPa. We set Δr=10mm at the working rotation speed ω=20,000rpm. ANSYS is employed for calculating values of objective and constraint functions.As for parameter setting, the maximum number of function evaluations is 360,α is set to 1/4, and other parameters remain unchanged as for the test function case.The resulting feasible nondominated solutions are obtained by applying CEIR2 to this engineering problem, as shown in Figure8. For comparison, we use the optimal LHS method to get 900 sample points at once and then evaluate values of the expensive objective and constraint functions, out of which the feasible nondominated solutions are also shown in Figure 8.Figure 8
Nondominated solutions obtained by CEIR2 and optimal LHS method for the design of the energy storage flywheel.CEIR2 found a total of 44 feasible nondominated solutions, which have a large distribution range in the objective space, namely,Ix∼16.1,24.7kg and Mx∼0.68,0.89kg⋅m2, and uniformly distributed along the approximated PF. In comparison, the optimal LHS spend 900 expensive function evaluations but found only 11 feasible nondominated solutions, which are inferior to that of CEIR2 both in terms of convergence and diversity, indicating CEIR2’s advantage in using computing resources.


## 7. Conclusion and Future Work

In this paper, a novel infilling point criterion EIR2 is proposed, which is combined with the Kriging approximate model, optimal Latin hypercube sampling, and PSO algorithm for dealing with expensive multiobjective optimization problems. The basic idea of EIR2 is to map a point in the objective space to a set in the EI space, and the R2 indicator of the set is defined as the fitness of this point. Due to its analytic form, the computation required for this criterion is very low. The proposed algorithm is applied to three suites of standard test functions. The results show that the algorithm can get more nondominated solutions with higher convergence and dispersion under the equal computing resources. Combined with the average probability of feasibility indicator, EIR2 can also be applied to constrained expensive MOPs. In the engineering optimization of energy storage flywheel, the nondominated solutions having a wider range and diversity are obtained exhausting relatively low computation cost.Future work will have to extend EIR2 to a parallel version that is capable of finding multiple infilling points [50–54] in one iteration with the help of parallel computers to further accelerate the speed of the optimization process. Besides, regarding the parameterization of the utility function, we use a fixed set of uniformly distributed weight vectors. This strategy runs smoothly if the geometry of PF is relatively regular, such as a simplex-like shape, but may struggle [55] if the true PF turns out to be irregular, such as having a sharp peak of a low tail or discontinuous, because multiple weight vectors cloud corresponds to one same solution. As a result, the diversity of the nondominated solution set cannot be guaranteed. A potential promising strategy [56–58] is to dynamically adjust the weight vectors during the optimization process so that the solutions are directed towards each part of the PF. Future work will be incorporating this strategy to our IPC to further improve its ability to maintain diversity.



---

*Source: 9474580-2020-06-27.xml*