A clustering-based surrogate-assisted evolutionary algorithm (CSMOEA) for expensive multi-objective optimization

This paper presents a novel surrogate-assisted evolutionary algorithm, CSMOEA, for multi-objective optimization problems (MOPs) with computationally expensive objectives. Considering most surrogate-assisted evolutionary algorithms (SAEAs) do not make full use of population information and only use population information in either the objective space or the design space independently, to address this limitation, we propose a new strategy for comprehensive utilization of population information of objective and design space. The proposed CSMOEA adopts an adaptive clustering strategy to divide the current population into good and bad groups, and the clustering centers in the design space are obtained, respectively. Then, a bi-level sampling strategy is proposed to select the best samples in both the design and objective space, using distance to the clustering centers and approximated objective values of radial basis functions. The effectiveness of CSMOEA is compared with five state-of-the-art algorithms on 21 widely used benchmark problems, and the results show high efficiency and a good balance between convergence and diversity. Additionally, CSMOEA is applied to the shape optimization of blend-wing-body underwater gliders with 14 decision variables and two objectives, demonstrating its effectiveness in solving real-world engineering problems.


Introduction
In the real world, many optimization problems (Cheng et al. 2017;Liu et al. 2017;Jin and Sendhoff 2009;Wang et al. 2023) involve multiple objectives, which are named multi-objective optimization problems (MOPs). Generally, the objectives conflict with each other, which means improvement of one objective will lead to the deterioration of other objectives (Zhou et al. 2011). Numerous methods (Tanino and Kuk 2003;Deb et al. 2002;Zhang and Li 2007;Li et al. 2018) have been developed to solve MOPs, and the methods can be classified into two categories, i.e., multiple criteria decisions making (MCDM) (Tanino and Kuk 2003) and multi-objective optimization evolutionary algorithms (MOEAs) (Deb et al. 2002;Zhang and Li 2007;Li et al. 2018). In MCDM methods, an MOP is scalarized into a single objective optimization problem (SOP). For MOEAs (Deb et al. 2002;Zhang and Li 2007;Li et al. 2018), work with a population and find a set of solutions to represent the whole Pareto front (PF). It takes a large number of fitness evaluations (NFEs) to obtain a set of satisfactory solutions.
However, in some engineering optimization problems, the fitness evaluations are quite time-consuming simulations, like computational fluid dynamics (CFD), finite element analysis (FEA), and so on. MOEAs face great challenges since they need a large number of evaluations. To address this issue, the mainstream is to apply surrogate models for the approximation of expensive functions (Jin 2005). Various surrogates can be developed, including polynomial regression (Gunst et al. 1996), Kriging (Martin and Simpson 2004), neural networks (Geman et al. 1992), support vector regression (SVR) (Gunn 1998), radial basis function (RBF) (Dyn et al. 1986), and so on. Surrogateassisted evolutionary algorithms (SAEAs) (Wang et al. 2019;Liu et al. 2023aLiu et al. , 2022aChugh et al. 2018;Zhang et al. 2010Zhang et al. , 2015Song et al. 2021;Li et al. 2022a, b;Guo et al. 2021Guo et al. , 2018Pan et al. 2018) were proposed to handle expensive SOPs and MOPs using classification or regression-based fitness approximation, which have been adopted successfully in engineering optimization (Wang et al. 2019;. Surrogate models are used to approximate the expensive objective functions in many SAEAs. A kriging-assisted reference vector-guided evolutionary algorithm (K-RVEA) (Chugh et al. 2018) and an efficient global optimization MOEA (MOEA/D-EGO) (Zhang et al. 2010), use Kriging models to approximate each objective function. Liu et al. (2022a) developed a reference vector-assisted Kriging model management strategy, and Song et al. (Song et al. 2021) proposed a two-archive with Kriging models (KTA2) for MOPs. Li et al. (2022a) suggested a two-stage SAEA (TS-SAEA), where RBF models are used to approximate expensive objective functions. In addition, an efficient dropout neural network-assisted indicator-based MOEA with reference point adaptation (END-ARMOEA) (Guo et al. 2021) is developed using the neural network. Meanwhile, the bagging technology of the surrogate models is also adopted to improve the accuracy of the approximation, such as a heterogeneous ensemble-based infill criterion for MOEA (HeE-MOEA) (Guo et al. 2018) and a bagging-based SAEA (B-SAEA) (Liu et al. 2022b).
In addition, several SAEAs use a surrogate model as the classifier, where the candidate solutions are classified into different categories and use the information in the design space. A classification-based SAEA (CSEA) (Pan et al. 2018) uses feedforward neural networks (FNNs) to predict the dominant relationship between candidate solutions and reference solutions. Zhang et al. (2015) suggested a classification-based preselection for the multi-objective evolutionary algorithm (CPS-MOEA). Li et al. (2022b) developed a classification SAEA (CSA-MOEA) for expensive MOPs, a classification tree is adopted as the surrogate model to predict candidate solutions.
Due to the excellent data mining and analysis ability of the clustering algorithms (Xu and Wunsch 2005), it can be naturally applied to data search to improve the performance of the MOEAs and SAEAs. Up to now, a series of MOEAs based on clustering has been proposed. The clustering ranking method (crEA) (Cai et al. 2015) uses a series of predefined reference lines to cluster individuals into subregions. The clustering-based MOEA (CLUMOEA) (Zhang et al. 2014) adopts the k-means clustering method to divide the current population into several clusters and only individuals in the same subproblem are allowed to perform crossover. An MOEA with clustering-based selection (EMyO/C) (Denysiuk et al. 2014) calculates the sum of the Euclidean distance between the individuals of every two clusters and then merges the two which have the minimum sum value into one. In addition, it is noteworthy that the clustering strategy is also used in SAEAs (Liu et al. 2021(Liu et al. , 2023b to improve performance. Liu et al. adopted a parameter-free clustering technique to generate several subswarms in a surrogate-assisted multi-population particle swarm optimizer (SA-MPSO) (Liu et al. 2021). Surrogateassisted two-stage differential evolution (SA-TSDE) (Liu et al. 2023b) adopts a clustering strategy to adaptively cluster the feasible solutions into several clusters to form some local region and improve performance by exploiting these local regions.
In this paper, we propose a clustering-based surrogateassisted MOEA (CSMOEA), where a bi-level sampling strategy is proposed to obtain the best information of samples in the design and objective space, simultaneously. Moreover, an adaptive clustering strategy is proposed to divide the current population into two groups (good and bad) and obtain the clustering centers in the design space. The contributions of this work are summarized as follows.
(1) A bi-level sampling strategy is proposed in CSMOEA that makes full use of information from both the design and objective space. This strategy includes two steps: screening better samples from the offspring in the design space, and selecting the best samples from the better samples in the objective space. This approach is different from other surrogate-assisted evolutionary algorithms that only use information from the objective space or design space for sampling, respectively. (2) An adaptive clustering strategy is proposed in CSMOEA to divide the current population into good and bad groups and obtain their information in the design space, respectively. The clustering centers are obtained, and an appropriate clustering number is selected by an evaluation index, which overcomes the dependence of the clustering algorithm on the clustering number.
The remainder of this article is organized as follows. In Sect. 2, we review the related work on the clustering algorithm and RBF model. The details of the proposed CSMOEA are presented in Sect. 3, and the results of experiments on mathematical cases and the engineering application are reported in Sect. 4. Section 5 concludes the article and draws the future work.

Related work
In this section, the general clustering algorithm is briefly introduced, and the RBF model is described.

Clustering algorithm
The clustering algorithm is a kind of unsupervised learning method for data objects with unknown labels (Xu and Wunsch 2005). Its main purpose is to achieve a high similarity between data objects of the same group, but a low similarity between data objects of different groups (Hua et al. 2018). The k-means algorithm (Macqueen 1965) is used widely because of its high efficiency in dividing the sample set into k clusters through an iterative process. At first, k initial clustering centers are selected randomly and the remaining data objects are divided into the nearest clusters according to their distance to the centers. Then, the average values of all objects in each cluster are calculated as the new clustering centers. In the next iteration, if the cluster centers of two adjacent iterations are different, it indicates that the samples need to be adjusted. Otherwise, all objects are classified correctly and the process ends.
However, the performance of k-means algorithm is highly dependent on k-value, and it is difficult to estimate k-value beforehand. To solve this challenge, we propose an adaptive method by using Davies-Bouldin index (DBI) (Krkkinen and Frnti 2000), and k-value corresponding to the minimum DBI is the most appropriate. DBI is a kind of non-fuzzy cluster evaluation index, which mainly uses the geometric principle to assess the clustering results. It is based on the degree of dispersion between clusters and the closeness of data objects in the same clusters. When the distance between data objects in same cluster is smaller and the distance between different clusters is larger, the DBI value is smaller. It means the data objects in each cluster are closer and in difference clusters are larger, and the clustering results are appropriate under this number of clusters.
DBI can be expressed by the following formula: where k is the number of clustering centers,i is the center of the i-th cluster, H is the number of individuals included in the i-th cluster, C i is the average distance between the individual and the cluster center in i-th cluster, x i;j is the jth individual in the i-th cluster. By using the DBI to determine the appropriate clustering number, the proposed adaptive clustering strategy overcomes the dependence of the clustering algorithm on the clustering number and ensures that the clustering result is of high quality.

Radial basis function
In this paper, the RBF model (Dyn et al. 1986) is used to approximate expensive objective functions. Studies (Jin and Simpson 2001) reveal that RBF can obtain better accuracy of approximation for high-dimensional problems than Kriging and SVR and its modeling speed is faster compared with the Kriging model. Given the data points ðx i ; y i Þjx i 2 < d ; i ¼ 1; 2; ::; N È É , the RBF is defined as follows: where x À x i k krepresents the Euclidean distance between the points x and x i , /ðÁÞ is the basic function. Many forms of the basis function can be used here. In this paper, the cubic form (/ðrÞ ¼ r 3 ) is adopted due to it was successfully employed in several surrogate-based algorithms Cai et al. 2019). In addition, the weight vector k ¼ ðk 1 ; k 2 ; :::; k N Þ T can be computed as follows: where y ¼ ðy 1 ; y 2 ; :::; y N Þ T is the output vector and U is the following matrix: . . .
pðxÞ is a linear polynomial in d variables with d ? 1 coefficients as in the formula: 3 The proposed algorithm In this section, a novel SAEA is proposed, named CSMOEA. The main framework of CSMOEA is first introduced, and detailed descriptions of CSMOEA are given.

The framework of CSMOEA
The framework of CSMOEA is shown in Fig. 1, which can be divided into two main parts. In the first part, the current population is divided into two groups, named good and bad solutions. The adaptive clustering strategy based on DBI is performed to determine the clustering centers in the design space. Then, the clustering centers will be used to prescreen offspring in the design space. In the second part, a bi-level sampling strategy is used to obtain the final candidate solutions. Variation operation is applied to produce offspring based on good solutions, and then, offspring are prescreened in the design space by clustering centers. The samples selected from the subsets with good centers are treated for the better population and survived into the next operation. The RBF model is used to approximate the objective values for each objective function, and one best sample is selected in the objective space to be evaluated by the real functions.
The pseudo-code of CSMOEA is given in Algorithm 1, which can be divided into the following steps: (1) Initialization (Lines 1-3): The initial population P is generated using the Latin hypercube sampling (LHS), which includes 11d-1 solutions, d is the number of decision variables. (2) Cluster the current population (Lines 5): The current population is divided into good and bad solutions by the non-dominated method (Deb et al. 2002), and getting their clustering centers by the adaptive clustering strategy in the design space, respectively. (3) The bi-level sampling strategy (Lines 6-15): Obtaining offspring from good solutions by crossover and mutation, and selecting better samples by clustering centers in the design space; approximate the objective values of the better samples by RBF models, and the best sample is obtained through objective space; If the best sample is not empty, restore the best sample and re-evaluate it to add into the database, otherwise repeat the sampling process. (4) Repeat (2)-(3) until the termination condition is met.

Clustering of the current population
In CSMOEA, a clustering strategy is proposed to cluster the current population into several groups. The clustering strategy includes two parts, i.e., classification criterion and adaptive clustering.
In the first part (Lines 1-4), we adopt the non-dominated method as a classification criterion to divide the existing solutions into two categories: good solutions and bad solutions. The non-dominated solutions in the first layer are taken as good solutions, and the rest are treated as bad solutions.
In the second part (Lines 5-16), we used an adaptive clustering strategy to obtain the clustering centers of good and bad solutions in the design space, respectively. An adaptive k-value (clustering number) based on the DBI to overcome the dependence of the clustering number. According to the definition of DBI in Eq. (3), when the value of DBI is the smallest, the clustering number is the most appropriate.
In this work, the range of clustering number is set as [1, g] and [1, b], where g represents a smaller value between the number of good solutions and the number of decision variables, b is the smaller value between the number of bad solutions and the number of decision variables. The best value of the k-value of good solutions could be expressed as: The details of the adaptive clustering strategy can be found in Fig. 2. Firstly, the current population is distinguished into two groups by classification (Fig. 2a). The good and bad solutions are used to determine the clustering centers in the design space through the adaptive k-means algorithm (Fig. 2b). Finally, the k-value with the lowest DBI value is selected as the number of clustering centers, and the location of the clustering center is obtained at the same time. Through the above operations, the optimal positions of clustering centers can be determined (Fig. 2d).

The bi-level sampling strategy
The bi-level sampling strategy is an important part of CSMOEA, where the population information in both design and objective space is considered. Firstly, the good solutions are treated as the parent population to produce offspring (Q). CSMOEA utilizes the information of the current population in the objective space to select the solutions with good convergence as the parents for the next generation. This is done to accelerate the convergence of the CSMOEA and improve the quality of the solutions obtained. Secondly, the offspring will be selected in the design space by clustering centers, only the individuals belonging to the good clustering centers will be retained as the better samples. Finally, the RBF model is used to approximate the function values of the better samples, and the approximate values will be used to select the best sample in the objective space.

Prescreening in the design space
To choose a set of better samples from offspring, a prescreening strategy based on the clustering centers in the design space is proposed, which is described in Algorithm 3 (Lines 2-4). CSMOEA differs from most SAEAs in its use of population information from both the design space and objective space in the sampling process. By considering population information from both spaces, CSMOEA can improve the reliability and prediction accuracy of the surrogate model.
The details of the selection process are shown in Fig. 3. The blue circle represents offspring, and the red and black triangle shows good and bad clustering centers, respectively. The offspring are generated by evolutionary operation (Fig. 3a), and the clustering centers of good and bad solutions are obtained by Algorithm 2 (Fig. 3b). Each offspring is assigned to the nearest center according to Euclidean distance (Fig. 3c). Finally, the offspring belonging to good clustering clusters will be retained as the better samples (Fig. 3d).

Selection of new samples
This part uses the RBF model to approximate the objective values of each objective function, and the approximate objective values are applied to screen out a new sample in the objective space.
The details of this process are presented in Algorithms 3 (Lines 5-10). First, RBF models are built by all points. To obtain better convergence, the non-dominated sorting method is used to select the candidate samples from better samples. The non-dominated solutions of the current population are chosen as the reference solutions, and the minimum angle (h) between the candidate samples and the reference solutions is calculated. Finally, the maximum minimum angle is adopted as a selection criterion, and the best sample corresponding to the maximum h is selected for real function evaluation.
To better illustrate the selection process in the objective space, the criterion of the non-dominated method and maximum minimum angle is shown in Figs. 4 and 5, respectively. As shown in Fig. 4, the better samples and the existing non-dominated solutions (P nd ) will be taken as a combination (Fig. 4a). Then non-dominated sorting method is employed to obtain the new PFs of them ( Fig. 4b), and the candidate samples are selected from new PFs (Fig. 4c).
To guarantee the diversity of the population, the candidate samples will be screened again. The maximum minimum angle is adopted as a selection criterion to obtain the best sample for real function evaluation. The selection process is shown in Fig. 5. P nd are non-dominated solutions of the current population, and Q obj nd are approximate objective values of candidate samples. b, c are minimum angles between them, and b [ c. By comparing the size of the angle, a candidate sample corresponding b is selected as the best sample. This selection criteria cloud ensures that the best sample has good diversity.

Empirical studies
In this section, a set of experiments are conducted to validate the effectiveness of CSMOEA. In the following, we first introduce the several SAEAs, parameter setting and benchmark tests. Then, the behavior study of CSMOEA is given and we empirically compare CSMOEA with several state-of-the-art SAEAs. Finally, CSMOEA is applied to an engineering problem to reveal its effectiveness. The compared algorithms include K-RVEA (Chugh et al. 2018), HeE-MOEA (Guo et al. 2018), EDN-ARMOEA (Li et al. 2022a), MOEA/D-EGO (Zhang et al. 2010) and CSEA (Liu et al. 2022b). All the test instances are implemented in PlatEMO (Tian et al. 2017a). The five compared algorithms are summarized as follows.
1) K-RVEA is a kriging-assisted reference vectorguided evolutionary algorithm, which uses Kriging to approximate each objective function and uncertainty information is provided to balance convergence and diversity.
2) HeE-MOEA is a heterogeneous ensemble-assisted MOEA, in which a support vector machine and two RBF networks are constructed to enhance the reliability of ensembles for uncertainty estimation. 3) EDN-ARMOEA is an efficient dropout neural network (EDN)-assisted indicator-based MOEA (Tian et al. 2017b), in which the EDN replaces the Gaussian process to achieve a good balance between convergence and diversity. 4) MOEA/D-EGO is a multi-objective optimization with a Kriging model, which decomposes a MOP into several single-objective optimization subproblems and uses the maximum expected improvement (EI) strategy to select several solutions for reevaluation. 5) CSEA is a classification-based SAEA, which uses an artificial neural network to predict the dominant relationship between candidate solutions and reference solutions.
The experiments are conducted on 21 test instances from test suite DTLZ (Deb 2005), WFG , and ZDT (Zitzler et al. 2000) with 2 and 3 objectives. For each test instance, the MFEs is set to 500, the numbers of the decision variables and objectives are the same as that introduced in the original literature (Deb 2005;Huband et al. 2006;Zitzler et al. 2000).
The inverted generational distance (IGD) (Bosman and Thierens 2003) and hyper-volume (HV)  as the performance indicators to assess the performance of   the compared SAEAs. In general, the lower the IGD value is, the better the solutions approximate the true PF. All experiments are conducted using MATLAB with an Intel (R) Core (TM) i7, 3.4 GHz CPU. The Wilcoxon rank-sum (WRS) is used to conduct statistical tests at a significance level of 5%. Symbols '' ? '' and ''-'' indicate that CSMOEA is significantly superior and inferior to the compared algorithm, ''&'' means that there is no significant difference between CSMOEA and the compared algorithm.

Parameter setting
The common parameter settings of all the compared algorithms are listed as follows: 1) The initial data size is set to 11d-1, where d represents the number of decision variables.
2) The population size is set to 100.
3) The maximum number of expensive function evaluations is set to 500. 4) The parameters for reproduction (crossover and mutation) are set to P c ¼ 1:0, P m ¼ 1=d,g c ¼ 20,g m ¼ 20.
In addition, for a fair comparison, we adopt the recommended setting in the original literature for specific parameters of compared algorithms.

Effect of the adaptive clustering strategy
In this part, we first investigate the effects of the clustering strategy of CSMOEA. The clustering centers of the existing solutions are obtained to prescreen the offspring in the design space, which is directly affected by the k-value and influences the performance of selection in the design space.
SMOEA is a variant of CSMOEA, which does not adopt a clustering strategy. The average IGD results of SMOEA and CSMOEA based on 30 independent runs on DTLZ problems are shown in Table 1. As shown in Table 1, the result demonstrates the superiority of the clustering mechanism of CSMOEA and shows the effectiveness of selection in design space.
To investigate the effects of the adaptive clustering number, the fixed k-value replaced the adaptive k-value in CSMOEA as a compared variant, termed CSMOEA(k). We test a series of fixed values f1; 3; 5; 7; 11g and the average IGD results of CSMOEA (k) and CSMOEA based on 30 independent runs on DTLZ problems are provided in Table 2, where the WRS test is also listed and the best results are highlighted.
As shown in Table 2, the fixed k-value of CSMOEA gets the best results in DTLZ4 (k = 3) and DTLZ7 (k = 7), and CSMOEA performs the best on the left problems. From the results, we can see that the results of different k-values have significant differences, which also proves the dependence of the k-means algorithm on clustering numbers. As we expected, the adaptive clustering number achieves the best overall results in the DTLZ problems, which shows the effectiveness of this strategy.

Effect of the bi-level sampling strategy
CSMOEA uses a bi-level sampling strategy that selects samples in both the design and objective space to effectively explore the search space and converge to highquality solutions. SMOEA, on the other hand, only selects samples for real function evaluation in objective space, while CMOEA only selects samples in the design space. By comparing these three algorithms, it can gain a better understanding of the role of the bi-level sampling strategy in CSMOEA and its impact on the performance of the algorithm.
The experiments are conducted on the DTLZ problem with three objectives. The true PF of DTLZ2 is easy to converge but needs more attention in diversity maintenance. The PF of DTLZ7 is discontinuous. Three algorithms have been ranked through the mean IGD values. The average rank on seven problems for each algorithm is listed in Table 3. CSMOEA is shown to be the top performer among the three algorithms, while CMOEA is the worst performer. The poor performance of CMOEA can be attributed to the fact that it only selects samples in the design space, which can lead to an uncontrollable number of new samples for real function evaluation. This can result in inefficient use of computational resources and poor convergence properties. Additionally, CMOEA does not utilize information on new samples in objective space, which further limits its ability to effectively explore the search space and converge to high-quality solutions. In contrast, SMOEA only selects samples for real function evaluation in objective space, which helps to control the number of new samples and improves the efficiency of the optimization process. However, the quality of the new samples in the design space cannot be judged due to the lack of information.
In summary, the bi-level sampling strategy makes use of the information of new samples in design and objective space and performs better than only using single information.

Results on DTLZ problems
The results of IGD and HV values achieved by six algorithms over 30 independent runs on DTLZ problems are summarized in Tables 4 and 5, where the best results are highlighted. The WRS is used to conduct statistical tests at a significance level of 5%. Among them, it can be observed that CSMOEA obtains the best performance on five test instances, followed by MOEAD-EGO and K-RVEA. CSEA, EDN-ARMOEA and HeE-MOEA have not achieved the best results among the DTLZ test problems.
(3.71e -1) = (5.93e -1) = (2.97e 2 1) = (6.06e -1) -(6.76e -1) ± / -/ = 1/4/2 1/3/3 0/4/3 0/4/3 0/6/1 Bold values indicate better results than other compared algorithms crucial in maintaining diversity during the optimization of DTLZ2 in CSMOEA. K-RVEA uses two sets of reference vectors to maintain a well-distributed set of solutions, and CSEA uses a radial space division-based strategy to select candidate solutions to maintain diversity. The obtained solutions of the compared algorithms are shown in Fig. 7. Based on Fig. 7 and HV values, the obtained solutions of above algorithms converge to the PF of DTLZ2. As shown in Tables 4 and 5, CSMOEA and K-RVEA are the top two algorithms for DTLZ4, followed by EDN-ARMOEA, CSEA, MOEAD-EGO and HeE-MOEA. As mentioned in selection 4.2, DTLZ4 is indeed known to be very sensitive to initial data, which can cause the prediction of the surrogate model to be inaccurate on some objectives. This is because the Pareto front of DTLZ4 is shaped like a twisted hyperplane, which makes it difficult for the optimization algorithms to explore and converge to the true PF. During the optimization process, most of the data tend to converge on certain objectives, while some points with good diversity have little impact on the surrogate models. This can lead to the optimization algorithms being trapped in local optima, which can result in suboptimal solutions. This can include using specialized initialization methods, adaptive surrogate modeling techniques, and diversity maintenance strategies to help the optimization algorithms overcome the challenges posed by the sensitivity of DTLZ4 to initial data.

Results on WFG problems
To further analyze the behavior of CSMOEA, we choose the WFG problems as test problems. WFG test suite consists of nine scalable, multi-objective test problems (WFG1-WFG9), which focus on some more pertinent problem characteristics. WFG1 is with convex, mixed and biased PF, which makes it difficult to achieve a set of solutions with good diversity. It is hard to obtain the PF of WFG2 due to its convex and disconnected shape. WFG3 replaces the disconnected PF of WFG2 with a continuous one. WFG4 is a multi-model problem. WFG6 and WFG 9 are non-separable-reduced and WFG5 is a deceptive problem. WFG7 and WFG8 are separable and non-separable problems, respectively. These problems have different characteristics that make them more challenging than DTLZ problems in terms of maintaining diverse solutions. The statistical results of compared algorithms over 30 independent runs on WFG problems are summarized in Tables 6 and 7. The WRS is used to conduct statistical tests at a significance level of 5%. It is obvious that the performance of CSMOEA and K-RVEA is similar, followed by CSEA, HeE-MOEA, EDN-ARMOEA and MOEADEGO.
According to the results of the comparison between CSMOEA and K-RVEA, CSMOEA achieves four better, four poorer, and one similar results on WFG problems. On WFG4, WFG6-7, and WGF9, it can be noted that CSMOEA performs better than K-RVEA, but performs worse on WFG1-3 and WFG5. For WFG8, CSMOEA and K-RVEA perform similarly; the solutions of the two algorithms are displayed in Fig. 9. We can see that CSEA produces superior answers on WFG1-3 and is the bestperforming algorithm on WFG1 when comparing CSMOEA and CSEA. HeE-MOEA has the best results on WFG3, we speculate that the heterogeneous ensemble surrogate improves approximation accuracy to approximate WFG3. In conclusion, CSMOEA performs similarly to K-RVEA on WFG issues and is competitive with the other four methods.

Results on ZDT problems
To further analyze the performance of CSMOEA in multiobjective problems, we use the ZDT problems as test problems. ZDT problems include six two-objective test problems, which introduce different difficulties for evolutionary optimization (Zitzler et al. 2000). We choose five unconstrained problems, referred to as ZDT1-ZDT4 and ZDT6.
The statistical results of compared algorithms over 30 independent runs on ZDT problems are summarized in Tables 8 and 9. At a significance level of 5%, statistical analyses are performed using the WRS. It is clear that Fig. 7 The PF of DTLZ2 associated with the median IGD value Fig. 8 The PF of DTLZ5 in the run associated with the median IGD value (9.10e -2) -(5.62e 2 2) 1 (1.46e -1) ?

Application on engineering problem
In this section, to verify the effectiveness of CSMOEA in engineering practice, we performed experiments on a multi-objective engineering problem, i.e., the shape design of blended-wing-body underwater gliders (BWBUGs) (Bachmayer et al. 2004). BWBUGs glide through the ocean by controlling their buoyance and converting the lift on wings into a propulsive force to propel themselves forward (Sun et al. 2015). For the shape design of BWBUGs, a larger lift-drag ratio (L/D) and internal volume (V) are expected. The L/D is affected by shape and angle of attack (AOA), and the AOA is set as 2 in this paper. The larger L/D means the glider has a larger glide angle, which is of great significance to improve the glide (1.44e ? 1) -(7.76e ? 0) = (7.26e 1 0) = (8.02e ? 0) -(7.78e ? 0) -(1.18e ? 1) ZDT6 2 10 9.2182e -1 5.3239e ? 0 1.1614e ? 0 6.0048e ? 0 3.3810e ? 0 4.6960e 2 1 ± / -/ = 0/5/0 0/4/1 0/4/1 0/5/0 1/4/0 Bold values indicate better results than other compared algorithms  Bold values indicate better results than other compared algorithms Fig. 10 The PF of ZDT1 in the run associated with the median IGD value Fig. 11 The PF of ZDT6 in the run associated with the median IGD value A clustering-based surrogate-assisted evolutionary algorithm (CSMOEA) for expensive… 10681 range. The larger the internal volume, the more energy and other equipment it can carry. However, as indicated by the results , the improvement in the lift-drag ratio was generally at the expense of the internal volume of the glider. The shape of BWBUGs can be divided into plane shape and wing shape. The wing-body fusion layout is reflected in the fuselage profile and the transition mode of the fuselage and wing. This section uses the Bezier curve to construct the plane shape of the fuselage.
We define the layouts of 11 sections as shown in Fig. 12, where L and D 1 are set as 1500 mm and 1000 mm, respectively, and the other nine parameters are design variables. Among them, the locations of section 1 and section 11 are fixed. And the locations of the other section are related to the variable Z 1 ; Z 3 and L. Selecting standard NACA0022 airfoil as the reference airfoil for all sections, the CST method is used for wing parameterization (Hicks and Henne 1978), and the airfoil is controlled by five Fig. 12 The shape of BWBUG full parameter diagram design variables. There are 14 design variables in the shape design, and the definition of them is shown in Table 10.
In order to obtain a longer range and carrying capacity, the shape design of BWBUGs mainly includes two design objectives: -L/D and -V. Hence, a 14-dimensional optimization problem is summarized below: min FðxÞ ¼ fÀL=D; ÀVg T s:t: x i 2 ½A section ; A plane ; i ¼ 1; 2; :::; 14 A section 2 ½A Low ; A Up ; A plane 2 ½P Low ; P Up where A Low and A Up represent the lower and upper boundaries of the design scopes of section, respectively. In addition, the ranges of selection and plane variables are shown in Formula (8) and (9), respectively.
For the shape design of BWBUGs, CSMOEA and K-RVEA are both operated with 300 NFEs and the population size is set to 100. For a fair comparison, K-RVEA and CSMOEA share the same 153(11d-1) initial sample points. After the run, the obtained solutions are shown in Fig. 13, where the circular and square are the PFs of K-RVEA and CSMOEA, respectively. Relatively, 15 and 22 Pareto solutions are obtained by K-RVEA and CSMOEA, respectively.
From Fig. 13, we can see that the solutions were obtained by CSMOEA have better convergence and diversity than K-RVEA. Moreover, four representative points located in the final PFs of these algorithms are selected, and their models are shown in Fig. 13 as well. In addition, Fig. 14 shows the plane and section of four representative results. Specifically, all sections of each optimized glider have the same normalized section. Finally, the pressure distributions corresponding to these solutions are shown in Fig. 15.

Conclusions
In this paper, a clustering-based surrogate-assisted multiobjective evolutionary algorithm, CSMOEA is proposed, which utilizes both design space and objective space information to select the best candidate solutions. The algorithm incorporates a clustering strategy to divide the current population into two groups and obtain clustering centers in the design space. A bi-level sampling strategy is used to prescreen better samples from offspring in the design space and approximate their objective function values using RBF models. The best samples are then selected in the objective space, and the best ones are obtained for real function evaluation. The performance of CSMOEA is compared with other surrogate-assisted evolutionary algorithms on widely used test problems DTLZ, WFG, and ZDT. The results demonstrate that CSMOEA achieves better performance than the compared SAEAs with the same number of expensive function evaluations. Additionally, the proposed algorithm is applied to an engineering problem and compared with K-RVEA. The results show that CSMOEA outperforms K-RVEA in terms of convergence and diversity. Overall, the proposed CSMOEA algorithm offers a promising approach to solving multi-objective optimization problems with computationally expensive objectives.
In further work, CSMOEA may adopt an effective many-objective sampling strategy for expanding its scope of use in many-objective optimization. Moreover, the experimental results of the clustering strategy reveal its effectiveness in expensive MOP, a further study on the expensive constrain of MOP will be made based on the clustering strategy.

Declarations
Conflict of interest The all authors declare that there is no conflict of interest regarding the publication of this paper.
Ethical approval This article does not contain any studies with human participants or animals performed by any of the authors.