A two-stage infill strategy and surrogate-ensemble assisted expensive many-objective optimization

Many optimization problems are expensive in practical applications. The surrogate-assisted optimization methods have attracted extensive attention as they can get satisfyingly optimal solutions in a limited computing resource. In this paper, we propose a two-stage infill strategy and surrogate-ensemble assisted optimization algorithm for solving expensive many-objective optimization problems. In this method, the population is optimized by a surrogate ensemble. Then a two-stage infill strategy is proposed to select individuals for real evaluations. The infill strategy considers individuals with better convergence or greater uncertainty. To calculate the uncertainty, we consider two aspects. One is the approximate variance of the current surrogate ensemble and the other one is the approximate variance of the historical surrogate ensemble. Finally, the population is revised by the recently updated surrogate ensemble. In experiments, we testify our method on two sets of many-objective benchmark problems. The results demonstrate the superiority of our proposed algorithm compared with the state-of-the-art algorithms for solving computationally expensive many-objective optimization problems.


Introduction
In industrial optimization problems, such as electric vehicle control problems [4], industrial scheduling [12], and robotics [24], there are multiple objectives to be optimized simultaneously. These multi-objective optimization problems (MOPs) can be mathematically formulated as follows:  1 subject to x ∈ R D (1) where R D represents the D-dimensional decision space, f i (x) represents the ith objective of individual x and M is the number of objectives. Objectives are often conflicting with each other, and one objective getting better may cause another to deteriorate. Pareto optimal solutions obtained by multi-objective optimization algorithms can trade off the different objectives, and they are called Pareto set (PS) in the decision space and Pareto front (PF) in the objective space, respectively. Various methods like a fast and elitist multiobjective genetic algorithm (NSGA-II) [8], a multiobjective evolutionary algorithm based on decomposition (MOEA/D) [38], and indicator-based selection in multiobjective search (IBEA) [42] have been proposed for solving MOPs with two or three objectives. These multi-objective evolutionary algorithms (MOEAs) are confronted with the lack of selection pressure when the objectives increase. To solve many-objective optimization problems (MaOPs) with more than three objectives, a lot of methods have been proposed, for example, NSGA-II/SDR [31] based on the strengthened dominance relation (SDR), AR-MOEA [30] based on an IGD non contributing solution detection (IGD-NS) indicator, reference vector guided evolutionary algorithm (RVEA) [3] which decomposes a MOP into several single-objective optimization problems. Usually, MOEAs require a lot of real evaluations of the objective function before finding a set of Pareto optimal solutions, which makes the MOEAs restricted in the expensive engineering problems, such as computing fluid kinetics (CFD) simulation [10] and engineering design optimization [11]. In these expensive problems, one simulation could take from minutes to hours [14]. Surrogate-assisted evolutionary algorithms (SAEAs), have been proposed to address the optimization of expensive problems. Commonly used surrogate models include polynomial regression model (PR) [12], radial basis function (RBF) [13], and Gaussian process (GP) [14]. For more descriptions of surrogate model, readers can refer to the review articles [13,14].
Surrogate-assisted evolutionary algorithms (SAEAs) can be divided into on-line and off-line optimization methods [33]. In the on-line surrogate-assisted optimization methods, a small number of expensive real fitness can be conducted during the optimization, and the newly generated samples can also be used to update the surrogates. Conversely, in the offline surrogate-assisted optimization algorithms, there are not new samples available [35]. In the on-line surrogate-assisted optimization algorithms, the methods to select individuals for expensive evaluation, are also known as infill sampling criterion, infill strategy or surrogate management. Infill sampling criteria, such as expected improvement (EI) [21,25], lower confidence bound (LCB) [18], and probability of improvement (PoI) [9], are widely used in the Kriging or Gaussian process (GP) assisted optimization algorithms.
On-line SAEAs are more flexible than off-line SAEAs as they have additional samples for surrogate management during the optimization process, which may have more opportunities to improve the performance of the algorithm than off-line SAEAs [14]. Therefore, in this paper, our work mainly focuses on the on-line surrogate-assisted optimization methods. A lot of on-line surrogate-assisted single-objective optimization methods have been proposed. In [16], Li selected the individuals with the best approximated fitness and maximum uncertainty for expensive evaluations and used the distance and fitness value information to calculate the uncertainty. In [29], Tian proposed a multi-objective infill strategy that considers both approximated fitness and uncertainty for solving high-dimensional expensive problems. In [22], Pan used teaching-learning-based optimization (TLBO) or differential evolution (DE) alternately to search for the best candidate solution, and the generation-based and individual-based strategies are used for surrogate management. In [37], Yu used a GP model as a coarse surrogate to learn the global landscape and an RBF model as a fine surrogate to learn the local feature of the fitness landscape. In the coarse search, approximate fitness together with the uncertainty of GP is used for environmental selection. In [19], Liu used the affinity propagation clustering technique to partition the population into several subpopulations and proposed an RBF-assisted learning strategy-based particle swarm optimizer (PSO) to update the particle in each subpopulation. For more methods of on-line single-objective SAEAs, please refer to Refs. [2,18,27,28,32,36].
Some on-line surrogate assisted-multi-objective optimization methods have been proposed. In ParEGO [15], a weight vector was selected at each iteration for optimization by the efficient global optimization (EGO). In MOEA/D-EGO [39], Gaussian process is built for each objective in the MOP and maximizing the expected improvement metrics are used to select test points for expensive evaluations. In KRVEA [5], the Kriging model is used to approximate each objective of MOPs, and approximated values or uncertainties provided by Kriging model are adaptively selected for expensive evaluations. In CSEA [23], Pan used a feedforward neural network (FNN) as a classifier to identify good solutions from the whole population. In [11], Guo used an efficient dropout neural network (EDN) to approximate the fitness of individuals and get their approximate uncertainties by randomly ignoring neurons in the neural network. In [34], Wang proposed an adaptive acquisition function in the Bayesian approach to solve expensive multi-objective optimization problems. In [26], Song used the two-archive evolutionary algorithm to optimize the population, and the differences between the individuals in the two archives are used for infill strategy. In [17], Lin used a global Kriging and several sub-models to construct the surrogate ensemble to approximate the objectives of the expensive multi-objective problems and proposed a reference vector-based infill strategy. In [10], Gu used the Kriging model to optimize the population for several generations, then the crowding degrees of the individuals in the radial space and the uncertainty information provided by the Kriging model are used for surrogate management. For ParEGO, MOEA/D-EGO, these two methods cannot obtain good performance on expensive MaOPs with more than three objectives. In KRVEA, the population is optimized by a set of uniformly distributed reference vectors, and the true evaluated individuals are selected from the population. When the shapes of reference vectors and PFs are not consistent with each other, the performance of the algorithm will be affected.
From the mentioned above, there have been some studies on multi-objective expensive optimization problems, but there are still some challenges in this field. First, different surrogate models are suitable for different types of expensive problems. Therefore, selecting an appropriate surrogate for the problems is important and difficult as there is no criterion for this work. Second, the multi-objective optimization framework also has an impact on the performance of the expensive MOPs. Third, as to surrogate management, infill sampling criteria directly affect the approximate accuracy of the surrogate and the optimal solutions finally found by the methods.
On one hand, considering the first challenge of the expensive multi-objective optimization problems as described above, a natural idea is to combine multiple base learners to form a strong learner. And in our previous work [40], ensemble constructs a strong learner by combining multiple base learners, which has proved to be superior to a single learner in terms of accuracy and robustness. Furthermore, an ensemble surrogate can provide approximate variance among different surrogates, which is important for surrogate management. On the other hand, in SAEAs, the update process of surrogate model contains some historical information. At present, there is no method to use these historical information for surrogate management. Motivated by these, we propose a method to solve the expensive many-objective optimization problems, named a two-stage infill strategy and surrogate-ensemble assisted expensive many-objective evolutionary optimization algorithm (TSEMO). The main work of this paper consists of two aspects: (1) The RBF models based on two different kernel functions are used as the base learners to construct the surrogate ensemble. And the surrogate ensemble is used to optimize the population. (2) A two-stage infill strategy is proposed to select individuals for expensive function evaluations. In the first stage, the individuals with the minimum approximate fitness or the maximum uncertainty are selected for expensive evaluations. The uncertainty is the approximate variance of the surrogate ensemble. In the second stage, the individual with the maximum uncertainty is selected for expensive evaluations and the uncertainty is the approximate variance of the historical surrogate ensemble.
The rest of this paper is organized as follows. Section "Related work" introduces the relevant work. In "A two-stage infill strategy and surrogate-ensemble assisted expensive many-objective optimization", the details of the proposed algorithm TSEMO are described. In "Numerical experiments", the numerical experiments are carried out, and the results are compared with several methods on multi-objective benchmark problems. Finally, the conclusion and future work are drawn in "Conclusion".

Related work
Different multi-objective optimization algorithms have a certain degree of impact on expensive multi-objective optimization problems. To solve the expensive many-objective optimization problems, we used NSGA-II/SDR [31] as the underlying optimization framework of our method. And the NSGA-II/SDR algorithm is described below.

NSGA-II/SDR
Before giving NSGA-II/SDR, we first introduce the original dominant relation in NSGA-II, which is defined below.
Definition 1: an individual x is said to dominate another individual y (denoted as x ≺ y) if and only if ∀ i, i = 1, 2, . . . , M, f i (x) ≤ f i (y), and there is at least one j, j ∈ 1, 2, . . . , M satisfying f j (x) < f j (y).
Due to the lack of selection pressure, the traditional NSGA-II algorithm based on the dominant relation deteriorates in solving MaOPs with more than three objectives. To solve this problem, Tian [31] proposes a strengthened dominance relation (SDR), which replaces the original dominance relation in NSGA-II, and can greatly improve the selection pressure of many-objective optimization problems, thus promoting the performance of the algorithm. SDR can balance the convergence and diversity of MaOPs well by adopting the specific niche technology. The definition of strengthened dominance relation is given below.
Definition 2: an individual x strengthened dominate another individual y (denoted as x ≺ S D R y) if and only if Con(x) < Con(y) θ xy <θ where , θ xy denotes the acute angle between the two individual x, y in the objective space,θ is the parameter which is set to the (|P|/2) th minimum element of SDR is irreflexive, antisymmetric, and nontransitive, respectively. NSGA-II/SDR is based on the algorithm framework of NSGA-II. The multi-objective algorithms obtain a set of non-dominant solutions that are consistently distributed and proximate to the Pareto front by performing the environmental selection. During the environmental selection, NSGA-II/SDR first uses the strengthened dominance relation to sort the individuals into several non-dominant layers (L 1 , L 2 , . . . , L i , . . .). We suppose that the parent population size is N, and the first |i − 1| layers are put into the next generation, where|L 1 ∪ L 2 ∪ · · · ∪ L i−1 | < N and |L 1 ∪ L 2 ∪ · · · ∪ L i | > N . Crowding distance [8] which is the same as the original NSGA-II, is used to select the ith layer into the next generation one by one until the population size reaches N . More details of NSGA-II/SDR can be referred to [31].
A two-stage infill strategy and surrogate-ensemble assisted expensive many-objective optimization Radial basis function(RBF) has been widely used in SAEAs. RBF fits well for nonlinear and higher-order problems and it has a low training complexity [41]. We use two different kernel functions, cubic and Gaussian function in RBF as the basic learners to construct the surrogate ensemble. The flowchart of TSEMO is given in Fig. 1 and the pseudocode of TSEMO is in Algorithm 1. In the main loop, TSEMO first uses Latin hypercube (LHS) [20] to generate 11D − 1 individuals. These individuals are evaluated by the expensive functions and kept in archive A 1 . The nondominated individuals in A 1 are kept in A 2 . We conduct the environmental selection of NSGA-II/SDR on the individuals in A 1 , and N individuals are selected as the initial population. Then the individuals in A 1 are used to train the surrogate ensemble for each objective of the MaOPs. After that, the population is optimized for w max generations. Next, a two-stage infill strategy is used to select individuals for expensive evaluations. These new samples are added to archive A 1 and used to update archive A 2 . The surrogate ensemble is re-trained by the samples in A 1 . At last, the population is revised with the surrogate ensemble. This process is repeated until the maximum number of expensive evaluations F E max is reached. In the subsequent subsections, we will give a detailed description of the methods on the surrogate ensemble-based optimization, choosing individuals for exact function evaluations (a two-stage infill strategy), and revising the population, respectively.

Surrogate-ensemble based optimization
In the surrogate-ensemble based optimization, all the individuals are evaluated by the RBF ensemble instead of the exact objective functions. Algorithm 2 gives the pseudocode of the surrogate-ensemble based optimization. In the surrogateensemble based optimization, we use simulated binary cross over [6] and polynomial mutation [7] to generate offspring population Q, and population Q is evaluated by RBF ensemble. The ith objective of individual x is calculated below: denote ith objective of the individual x evaluated by RBF models with cubic kernel function and Gaussian kernel function respectively. Then the objective functions of P ∪ Q are normalized by

Algorithm 1 The pseudocode of TSEMO
Input: Maximum number of expensive function evaluations (F E max ); population size N ; Output: Solutions in A 2 1: Use LHS to generate the 11D − 1 initial individuals and initialize the number of function evaluations F E = 0, set A 1 = ∅, A 2 = ∅, A 3 = ∅ 2: Evaluate the initial individuals with the expensive objective functions, add them to A 1 , add all nondominated individuals of the samples into the archive A 2 ; and update F E = 11D − 1; 3: Select N individuals in A 1 as the initial population. 4: Use the samples in A 1 to train surrogate ensemble 5: while F E < F E max do 6: Add the current population P into the archive A 3 7: Implement the surrogate-based optimization on the population P and get the population P (Refer to Algorithm 2) 8: Add the population P into the archive A 3 9: Select u individuals from the population P by the two-stage infill strategy and re-evaluate them with the expensive functions (Refer to Algorithm 3) 10: Add u individuals obtained from step 9 to A 1 and update A 2 11: Update the surrogate ensemble with archive A 1 12: Implement the population revising strategy on the archive A 3 (Refer to Algorithm 4) 13: Reset the archive A 3 = ∅ 14: end while where F max and F min are the vectors of the maximum and minimum objective functions of population P ∪ Q respectively. Then the environmental selection of NSGA-II/SDR is used on the combined population L = P ∪ Q to select the next generation. If the maximum number of generation w max is reached, the surrogate-ensemble based optimization is ended and the population P is as the output.

Algorithm 2 surrogate-ensemble based optimization
Input: population P Output: population P 1: while w < w max do 2: Generate offspring Q using simulated binary crossover and polynomial mutation 3: Combine parent and offspring populations L = P ∪ Q, and evaluate the objective values of L by RBF ensemble 4: Normalize the population L 5: Use environmental selection on L 6: Update w = w + 1 7: end while

A two-stage infill strategy
After the surrogate-ensemble assisted optimization, u individuals are selected for exact evaluations. In this paper, we propose a two-stage infill strategy for selecting individuals, and the pseudocode of the infill strategy is given in Algorithm 3.
In the first stage, the convergence of population or the uncertainty of the surrogate ensemble is adaptively con- Fig. 1 The flow chart of TSEMO sidered according to the distribution of the population. By considering the diversity of the population, we use k-means to divide the population P into K clusters in objective space. An individual from each cluster with the minimum Euclidean distance (ED) to the origin or the maximum approximate variance ε 1 is selected for expensive function evaluations and ε 1 is calculated as follows: To select individuals from each cluster, the population P is sorted into several non-dominated fronts by the original dominant relation in NSGA-II. For each cluster, if there is only one non-dominated front, which means the individuals are non-dominated with each other, the individual with the maximum uncertainty will be selected. Otherwise, the individual with the minimum Euclidean distance (ED) to the origin will be selected. Then the distance d 1 , which is the minimum distance from each of the K individuals to the samples in A 1 in the decision space is calculated. If the distance d 1 is greater than the threshold δ, the individual selected from the first stage will be evaluated by the exact objective function. After that, the individuals calculated by the expensive function are added to A 1 . Else if the distance d 1 is less than or equal to δ, the individual will not be reevaluated. Here we set the value of δ to 10 −6 . A simple example is given in Fig. 2 to show the criterion to select individuals in each cluster. In Fig. 2, there are three clusters, c 1 , c 2 , c 3 , and one individual from each cluster will be selected. In cluster c 1 , because there are two non-dominated fronts, the individual a 1 with the small- est ED value will be selected. In cluster c 2 , there is only one non-dominated front, thus the individual a 3 with the largest ε 1 will be selected. In cluster c 3 , a 6 will be selected as in cluster c 2 .
In the second stage, the u − K individuals with the maximum uncertainty of ε 2 are selected. In this stage, the surrogate ensemble is updated by A 1 and the population P is reevaluated by the surrogate ensemble. The uncertainty ε 2 is calculated below: where f 1 i (x) denotes the ith objective value of the individual x approximated by the surrogate ensemble in the first stage and f 2 i (x) denotes the ith objective value of the individual x approximated by the surrogate ensemble in the second stage. The distance d 2 of u − K individuals is calculated in the same way as d 1 . If d 2 is greater than the threshold δ, the individuals will be evaluated by the exact objective function and they will be added to A 1 . The approximated variance from the first and the second stage is used to select individuals for exact evaluation and the aim is to exploit the historical information of the surrogate ensemble.

Algorithm 3 A two-stage infill strategy
Calculate the number of non-dominance layers τ in the c i cluster 5: if τ = 1 then 6: Select the individual with the maximum ε 1 value in each cluster 7: else 8: Select the individual with the minimum Euclidean distance to origin in each cluster 9: end if 10: Calculate the minimum distance d 1 from the individual to the samples in archive A 1 11: if d 1 > δ then 12: Evaluate the individual using the exact objective functions and add it to A 1 . 13: end if 15: end for 16: Update the surrogate ensemble with archive A 1 17: Select u − K individuals in P with the maximum ε 2 value 18: for j = 1 to u − K do 19: Calculate the minimum distance d 2 from the individual to the samples in archive A 1 20: if d 2 > δ then 21: Evaluate the individual using the exact objective functions and add it to A 1 . 22:

Revising the population
In surrogate-assisted optimization, the approximated errors of objectives may lead to the wrong direction of search. And the approximated objectives of the surrogate ensemble may become increasingly accurate as the training samples in A 1 increase. Considering the above two points, we propose a strategy to revise the population. First, the population P before the surrogate-based optimization and the population P are combined into A 3 . After the two-stage infill strategy, the surrogate ensemble will be updated by the samples in A 1 . Then the objectives of the individuals in A 3 are re-evaluated by the updated surrogate ensemble. Finally, environmental selection of NSGA-II/SDR is performed on A 3 and N individuals with better convergence and diversity are retained. Figure 3 gives an example to explain the effect of revising the population. Figure 3 shows a two-dimensional decision space. In Fig. 3, region 1 represents the region occupied by individuals with worse convergence and smaller approximate error, region 2 represents the region occupied by individuals with better convergence and larger approximate error, region 3 represents the region occupied by individuals with better convergence and smaller approximate error, and region 4 represents the region occupied by individuals with worse convergence and larger approximate error. The five-pointed star represents the individuals whose exact function values are convergent. In the process of surrogate-ensemble based optimization, the optimization of the population is based on approximated value of the surrogate ensemble. Besides, due to the limited training samples at the initial stage of optimization, approximated error of the surrogate ensemble is large, the individuals in region 4 with worse convergence will be discarded. However, region 4 may contain some individuals with good convergence. By revising the population, the initial population is re-evaluated with the updated surrogate ensemble after the two-stage infill strategy. At this time, the approximate accuracy of the surrogate ensemble is improved due to the increased samples. In this way, the individuals represented by the five-pointed star in region 4 can be retained. And the effectiveness of this step will be validated in the experimental section.

Numerical experiments
In this section, we first describe the performance metrics and parameter settings. Then we test the performance of TSEMO by comparing it with several state-of-the-art MaOPs algorithms on DTLZ and MaF test problems. After that, we conduct experiments on different values of K to determine the thresholds used in the two-stage infill strategy in TSEMO. Then the effectiveness of the two-stage infill strategy is studied. At last, the efficiency of revising the population is examined. Wilcoxon rank-sum test is used to compare the And the symbol '+' indicates that the compared algorithm performs statistically better than TSEMO. The '−' represents that the compared algorithm performs statistically worse than TSEMO. And the '≈' means that there is no significant difference between the performance of TSEMO and the compared algorithm.

Performance metrics and experimental settings
In this paper, inverted generational distance (IGD) [1] is used as the performance indicator to compare the performance of different algorithms. The IGD indicator can provide information about convergence and diversity of the non-dominated solutions, and it is defined as follows: where P * denotes the evenly distributed point set on the true PF, Ω is the non-dominated solutions obtained by the algorithm and dis(x, P * ) represents the minimum Euclidean distance between x and the point in P * . The number of the reference points in the true PF is set as recommended in KRVEA [5].
We also use HV as the performance indicator to compare the performance of our method with other algorithms. The HV indicator measures the volume of the hypercube dominated by the non-dominated solutions obtained by the algorithm, and it is defined below: where V ol(·) indicates the Lebesgue measure, and z r = (z r 1 , z r 2 , . . . , z r M ) is a reference point in the objective space that is dominated by Ω. z r is set as recommended in AR-MOEA [30].
The experimental settings are given below: (1) Genetic operators: in this paper, the simulated binary crossover [6] and polynomial mutation [7] are used to generate offspring population in all the algorithms. The distribution index of crossover is set to 20, and the distribution index of mutation is set to 20. The crossover probability p c is set to 1.0, and the mutation probability p m is set to 1/D.

Experimental results
In this section, we compare TSEMO with three surrogateassisted algorithms namely, KTA2, CSEA, KRVEA, and two MOEAs called NSGA-II/SDR, RVEA on DTLZ1-DTLZ7, and the comparison results are given in Tables 1 and 2 Table 1 show that TSEMO obtained the best results on all the test problems of DTLZ2, DTLZ5, and DTLZ6, and TSEMO obtained the best results on most of the test problems of DTLZ4 and DTLZ7. For two multi-objective algorithms, NSGA-II/SDR and RVEA, TSEMO obtained 25 better results than NSGA-II/SDR and RVEA among 35 test problems, which demonstrate that surrogate models are effective on most test problems. For a few test problems like DTLZ1 and DTLZ3, TSEMO needs more expensive function evaluations to find the true PF. TSEMO got 20 better results and 12 worse results than KTA2. TSEMO got 18 better results and 6 worse results than KRVEA. Compared to CSEA, TSEMO obtained 23 better results and 11 worse results. Then the statistical results of HV in Table 2 also show that TSEMO performed the best among all the algorithms. Therefore, for DTLZ test suite, TSEMO is more efficient than the three representative surrogate-assisted many-objective optimization algorithms, KTA2, KRVEA, and CSEA. Next, the parallel coordinates plot of the final  2823e-1 (1.99e-3 1403e-2 (6.38e-3) -1 (3.46e-1)           The best results are highlighted  -1 (1.84e-3) -1 (1.83e-3) -1 (6.72e-3) -1 (6.59e-3 -1 (1.23e-2) -1 (1.06e-1 -1 (4.48e-2) 10 8.3660e−1 (6.01e−2) 9470e-1 (1.59e-3 3543e-1 (2.25e-3) 0801e-1 (1.44e-3) 0051e-1 (9.49e-4)       The best results are highlighted Table 3  Statistical results for IGD values obtained by KTA2                       non-dominated solutions obtained by TSEMO and the compared algorithms for DTLZ2 with 10 objectives is presented in Fig. 4. From it, we can see that the maximum value of objectives for all the non-dominated solutions obtained by TSEMO is smaller than that of the compared algorithms and the non-dominated solutions of TSEMO have a higher den-sity than the two non-surrogate assisted MOEAs algorithms. These mean TSEMO can achieve a set of well convergent and evenly distributed solutions and has a better capability on the convergence and diversity. As to MOEAs, NSGA-II/SDR and RVEA have a lower density of solutions which means they have a worse distribution. Furthermore, to verify the effectiveness of TSEMO on many-objective optimization problems, we conduct experiments on MaF1-MaF7, which are the modified versions of the DTLZ. And the results are presented in Table 3. For MOEAs without surrogates, TSEMO obtained 15 better results than RVEA and NSGA-II/SDR among 21 test problems, which further verified the effectiveness of adding the surrogate in MOEAs. We can see that our proposed TSEMO can obtain 11, 10 better results among 21 test problems than KRVEA and CSEA, and lose to win on 6 and 7 problems than KRVEA and CSEA, respectively, which shows that our proposed method is more efficient than both of these algorithms on MaF1-MaF7 test problems.

Parameter analysis
In the two-stage infill strategy, a total of u individuals are selected for expensive objective evaluations. For fairness, we set u to the same value of 5 as in the KRVEA [5]. In the first stage of the infill strategy, K individuals are selected for real evaluations. And the remaining u−K individuals are selected in the second stage. The number of individuals selected in the first stage will affect individuals with large uncertainty found in the second stage and ultimately affect finding a good Pareto front. Therefore, we first conducted experiments with u = 5, and different numbers of individuals selected in the first stage for expensive evaluations, i.e., K = 5, 4, 3, 2, 1. When K = 5, it means that individuals are selected for real evaluations only in the first stage. Figure 5 shows the mean IGD values of 20 independently runs obtained by TSEMO with the different numbers of individuals selected for exact evaluations in the first stage. From Fig. 5, we can see that when the number of individuals selected for exact evaluation in the first stage is set to 4, the performance of TSEMO is the best.
Then we further analyze the performance variation of the algorithm, when u, K increase further, and the number of infill samples in the second stage remains unchanged. Fig- The best results are highlighted ure 6 shows the mean IGD values of 20 independently runs obtained by TSEMO with the different numbers of u individuals selected for exact evaluations. From Fig. 6, we can see that there are no significant differences among the performance of TSEMO on 3, 4, 6, 8, 10 objectives of DTLZ2 when u increases further. On the low-dimensional 3 objectives of DTLZ1, the algorithm performs best when u = 6, K = 5. However, on the high-dimensional DTLZ1 of 4, 6, 8, 10 objectives, the algorithm performs best when u = 5, K = 4. Thus, we set u = 5 and K = 4 in our method for solving the expensive many-objective optimization problems.

Effectiveness of the two-stage infill strategy
In this part, we investigate the effects of the two-stage infill strategy by comparing TSEMO with another variant TSEMO-I, which only selects individuals for exact evaluations in the first stage. Table 4 gives the statistical results of our proposed TSEMO, and TSEMO-I. From Table 4, The best results are highlighted we can see that TSEMO obtained 7/35 better results than TSEMO-I. Compared to TSEMO-I, our TSEMO obtained 0 worse and 28 comparable results than TSEMO-I. TSEMO-I without the second stage sampling strategy has worse performance mainly on DTlZ1 with 8 objectives, DTLZ2 with 3, 4, 6, 8 and 10 objectives, DTLZ3 with 8 objectives, and DTLZ5 with 10 objectives. The reason we analyze is that DTLZ2 is used to investigate the diversity and distribution of the algorithm and the second stage sampling strategy has a potential benefit to improve the diversity of algorithms. For other problems, like DTLZ4 and DTLZ7, TSEMO-I showed comparable performance with TSEMO. Therefore, the effectiveness of the second stage infill sampling strategy cannot be ignored and the two-stage infill strategy is adopted to be used in our method.

Effectiveness of the two different kernel functions used in the surrogate ensemble
In this paper, we adopt two different kernel functions of RBF models to construct a surrogate-ensemble and utilize the ensemble to approximate each objective function. To see the efficiency of using surrogate-ensemble, we compare TSEMO with another variant TSEMO-II, where the surrogate-ensemble is replaced with GP, the uncertainty is provided by GP, and the two-stage infill strategy and the strategy of revising the population are employed. The average IGD results of TSEMO-II and TSEMO based on 20 independent runs are presented in Table 5. As can be seen in Table 5, TSEMO performed 12 better, 3 comparable, and 0 worse results among a total of 15 problems than TSEMO-II.
TSEMO performs well on most test problems. The results showed that TSEMO can achieve a better balance between convergence and diversity. This may be attributed to that the two different kernels of RBF can provide more accurate fitness predictions and the uncertainty information provided by our method is useful. Thus, we use the surrogate ensemble to approximate each objective function.

Effectiveness of revising the population
In this section, we investigate the effects of revising the population by comparing TSEMO with another variant TSEMO-R, which removes the strategy of revising the population. The average IGD results of TSEMO-R and TSEMO based on 20 independent runs on DTLZ1-DTLZ7 problems with 3, 4, 6, 8, and 10 objectives are presented in Table 6. As can be seen in Table 6, TSEMO performed 5 better , 25 comparable, and 0 worse results among a total of 35 problems than TSEMO-R. TSEMO performs better on DTLZ2 with 10 objectives, DTLZ5 with 8 and 10 objectives and DTLZ6 with 3 and 4 objectives. The results showed that the strategy of revising the population is effective for the surrogate-assisted multi-objective optimization problems.

Conclusion
In this paper, we propose a two-stage infill strategy and surrogate-ensemble assisted expensive many-objective optimization (TSEMO) method. To improve the robustness of the surrogate model, we use two kernel functions of RBF to con- The best results are highlighted struct a surrogate ensemble. After optimization, a two-stage infill strategy is proposed to select individuals for expensive evaluations. The approximated variance of the surrogate ensemble and approximated fitness are considered adaptively in the first stage. And historical information of the approximated variance by the surrogate ensemble is considered in the second stage. To avoid the approximated errors from misleading the search direction, a strategy of revising the population is proposed.
We compare our method with three surrogate-assisted MOEAs and two non-surrogate-assisted MOEAs on two sets of test problems. The results show that TSEMO achieves better performance than the compared algorithms on the limited number of expensive function evaluations. In our work, the surrogate-ensemble, two-stage infill strategy, and strategy of revising the population are effective for solving expensive MOPs. However, for multimodal problems, like DTLZ3 and MaF4, TSEMO still needs more expensive function evaluations to find the true PFs. Therefore, a more effective optimization process and infill strategy should be developed for solving these complex problems.

Declarations
Competing interests The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.