An Efficient Surrogate-Based Optimization Method for BWBUG Based on Multifidelity Model and Geometric Constraint Gradients

Performing shape optimization of blended-wing-body underwater glider (BWBUG) can signiﬁcantly improve its gliding performance. However, high-ﬁdelity CFD analysis and geometric constraint calculation in traditional surrogate-based optimization methods are expensive. An eﬃcient surrogate-based optimization method based on the multiﬁdelity model and geometric constraint gradient information is proposed. By establishing a shape parameterized model, deriving analytical expression of geometric constraint gradient, constructing multiﬁdelity surrogate model, the calculation times of high-ﬁdelity CFD model and geometric constraints are reduced during the shape optimization process of BWBUG, which greatly improve the optimization eﬃciency. Finally, the eﬀectiveness and eﬃciency of the proposed method are veriﬁed by performing the shape optimization of a BWBUG and comparing with traditional surrogate-based optimization methods.


Introduction
With the development of human society, the resources on land are exhausted. Because the ocean has unique advantages in resources, environment, and space, countries in the world are gradually paying attention to the development and utilization of marine resources. As a new type of underwater vehicle, Autonomous Underwater Glider (AUG) [1][2][3][4] has the advantages of long endurance, low manufacturing cost, and high reuse. It can sail in the ocean for a long time with high efficiency. At present, it has been widely used in marine environmental investigation and monitoring [5], such as sea temperature, salinity, ocean current speed, and marine environmental pollution. It can also be used to observe submarine volcanic eruptions, detect glaciers, and serve as data transmission nodes. e basic principle of the underwater glider is to convert the hydrodynamic lift generated by the wing into the driving force through adjusting its net buoyancy. erefore, the hydrodynamic shape has an important impact on the overall performance and navigation distance of the underwater glider. e blended-wing-body underwater glider (BWBUG) [6][7][8] has a fuselage with a flat airfoil profile, and the hydrofoil and the fuselage are smoothly fused together. Compared with the traditional underwater glider composed of revolving body, hydrofoil, and control surface [9,10], its outstanding feature is that it has a larger lift-drag ratio [11,12], which leads to a longer distance in the horizontal direction when it descends the same depth in the water. In recent years, it has received more and more attention and research.
To further improve the lift-drag ratio of the underwater glider, the optimization design method is widely applied in the shape design of the underwater glider. It is the top-level structure of the shape optimization design, which controls the operation of the whole optimization design and determines the development direction of the optimization. Generally, the optimization method can be divided into gradient-based optimization method [13], gradient-free optimization method [14], and surrogate-based optimization (SBO) method [15]. e gradient-based optimization method mainly determines the search direction of the next iteration by the gradient information of objective function and constraint conditions. It is efficient in finding local minima for highdimensional nonlinear problems defined by continuous smooth functions. e typical gradient-based optimization methods include the interior point method, trust region method, and sequence quadratic programming (SQP). e tests show that the computation complexity of the gradientbased optimization method is linearly proportional to the number of variables [16,17]. However, these optimization algorithms require the user to provide the gradient information and to guarantee the function smoothness of constraints and objectives. For the optimization problems of BWBUG, gradients might not be available when the lift-drag ratio is calculated by the CFD method which is a black box for most users. Although gradients can always be approximated with finite differences, these approximations suffer from potentially significant inaccuracies (a truncation error of O(h) or O(h 2 ) when second-order). erefore, the gradient-based optimization method is rarely found to be applied in the shape optimization design of BWBUG at this stage.
Different from the gradient-based optimization methods, the gradient-free optimization method does not need gradient information. erefore, it is useful when gradients are not available, such as when dealing with blackbox functions [18]. Another major advantage of gradientfree methods is that they tend to be robust to numerical noise and discontinuity, making them easier to use than gradient-based methods. ere is a wide variety of gradientfree methods which can perform a local or global search, use mathematical or heuristic criteria, and be deterministic or stochastic [19]. e gradient-free optimization methods mainly include the Nelder-Mead algorithm which is a local search algorithm based on heuristics, the generalized pattern search (GPS) which is based on mathematical criteria, and the evolutionary algorithms (such as the genetic algorithms and particle swarm optimization algorithms) which are based on the evolution of a population of designs. For the shape optimization design of BWBUG, Fu et al. [20] used NSGA2 (Nondominated Sorting Genetic Algorithm) to optimize the shape of the autonomous underwater glider. Tang et al. [21] used BESO (Bidirectional Evolutionary Structural Optimization) to optimize the shape of the underwater glider. However, the overall cost of gradient-free optimization is sensitive to the cost of the function evaluations because they require many iterations for convergence, and the number of iterations scales poorly with the number of design variables. It is found that the computational complexity of the gradient-free optimization methods increases exponentially with the number of variables [16,17]. erefore, the gradient-free optimization method is very time-consuming when it is applied to the shape optimization of BWBUG with a large number of variables.
To reduce the computational cost in the optimization while maintaining the advantages of the gradient-free optimization method, the surrogate-based optimization method [15,22] has been widely developed in recent years. It mainly uses the surrogate model to replace the time-consuming numerical simulation analysis and adds new sample points according to certain criteria in the optimization process to update the surrogate model cyclically. ere are many methods to build the surrogate model such as polynomials, radial basis functions, and kriging. For the shape optimization design of BWBUG, Wang et al. [23] established the hydrodynamic surrogate model of flying wing underwater glider using Gaussian kernel function and optimized it by particle swarm optimization algorithm. Sun and Li et al. [24,25] used the Kriging model to realize the hydrodynamic shape optimization of the underwater glider. Zhang and Li et al. [26,27] applied the surrogate-based optimization method into the shape design of BWBUG and obtained an optimal shape design. ese results show that the SBO method can effectively reduce the call times of computationally expensive functions in the optimization process, but one potential issue with the surrogate model is the curse of dimensionality [28]. With the increase of design variables, the CFD evaluations are significantly increasing to build the surrogate model. Considering the high-fidelity CFD evaluation is time-consuming for BWBUG, the SBO method still needs to be improved to reduce the computational cost when it is applied to the shape optimization design of BWBUG.
To further improve the optimization efficiency of the surrogate model and explore the high lift-drag ratio characteristics of BWBUG, an efficient surrogate-based optimization method based on the multifidelity model and geometric constraint gradient information is proposed to reduce the calculation times of high-fidelity model and geometric constraints in the shape optimization process of BWBUG. Section 2 describes the optimization model of BWBUG, and Section 3 introduces the problems existing in the traditional surrogate-based optimization method for the shape optimization design of BWBUG. Section 4 describes the proposed efficient surrogate-based method in detail. In Section 4, firstly, the parametric modeling of BWBUG is carried out, and the gradient analytical expression of geometric constraints is given. Secondly, the high-fidelity and low-fidelity CFD models of BWBUG are established, and the nested Kriging surrogate model is constructed. Finally, based on the criteria for infilling sample points, the shape optimization framework of BWBUG is built. Section 5 performs the shape optimization design of a BWBUG and compares the optimization results to demonstrate the effectiveness and efficiency of the proposed method. Section 6 summarizes the whole paper.

Optimization Model of BWBUG
For the BWBUG, its fuselage and wing are connected smoothly, and the cross section at each spanwise position is the airfoil section, which leads to a high lift-drag ratio. Figure 1 describes the geometric shape of BWBUG, and the blue curves show the airfoil profiles at nine spanwise locations.
e main content of the shape optimization design for the BWBUG above is to carry out parametric modeling of geometric shape and to use the efficient optimization algorithm to obtain the optimal hydrodynamic performance on the premise of meeting certain geometric constraints. e specific mathematical optimization model of BWBUG can be expressed as follows: where L/D is the lift-drag ratio of the BWBUG, and the optimization objective is to maximize the lift-drag ratio. x is the shape parameterization variable with the variable range of [x l , x u ]. g(x) denotes geometric constraints, which need to meet certain design requirements such as thickness and volume.

Problems Existing in the Traditional SBO Method
To solve the optimization model of BWBUG, the main idea of the traditional surrogate-based optimization method is to use the surrogate model to replace the time-consuming numerical simulation analysis in the optimization design. Its detailed optimization process can be described as follows: (1) e Design of Experiments (DOE) is used to generate the sample points in the design space, and the response values at these points are obtained by running accurate numerical simulations. en, the initial surrogate model is established based on these data. For a better understanding, the whole optimization process is shown in Figure 2.
Compared with the intelligent optimization method, the optimization efficiency can be greatly improved by establishing the surrogate model of lift-drag ratio and adopting the optimization process in Figure 2, but the following two problems still need to be solved. Problem 1. Figure 2 mainly uses the high-fidelity CFD model to calculate the lift-drag ratio at different sample points and builds a surrogate model of the lift-drag ratio based on these data to replace the original CFD model for optimization design, which greatly reduces the calculation consumption in the optimization process. However, if all the sample data are calculated through the high-fidelity CFD model, the calculation is still very time-consuming. Problem 2. Figure 2 often uses gradient-free optimization algorithms (such as NSGA and PSO) to solve the suboptimization problems. erefore, the calculation times of the geometric constraints have an exponential relationship with the number of optimization variables in solving the suboptimization problems. References [16][17] show that 8 optimization variables need to be calculated with thousands of geometric constraints. Even if it only takes ten seconds or a few seconds to regenerate the geometry each time, the amount of geometric constraint calculations in the whole suboptimization process is still considerable.

Efficient Surrogate-Based Optimization Method
When the surrogate model is used to optimize the shape of the BWBUG, if the high-fidelity CFD model is directly used for the hydrodynamic analysis, the calculation cost is high. In addition, the amount of geometric constraint calculation in the process of solving the suboptimization problem is large. In this section, the shape parametric model of the BWBUG is established, and the analytic expression of the geometric constraint gradient information is derived. A multifidelity CFD model composed of high-fidelity and lowfidelity CFD models is established to use the low-fidelity model to predict the results of the high-fidelity model. Finally, an efficient surrogate-based optimization method based on the multifidelity model and geometric constraint gradient information is proposed, which can greatly reduce the calculation times of the high-fidelity model and geometric constraints in the shape optimization process of BWBUG. angle, and spatial position of each airfoil section in the parametric modeling, the shape of BWBUG can be deformed in the optimization design process. In the following section, the geometric parametric modeling process is described in detail, and the related parametric variables are defined.

Shape Parameterization and Constraint Gradient
(1) Let D K ∈ R 2 , k � 0, . . . , m be the m data points of the initial airfoil, and arrange them in sequence with the airfoil trailing edge as the starting and ending points. A(ξ) is defined as the B-spline fitting curve of these m data points, which is determined by n control points P i , (i � 1, . . . , n). By changing the value of P i , the curve shape of A(ξ) can be changed. Figure 3 shows the data points of one airfoil profile, the B-spline fitting curve, and the corresponding control points. To sum up, the shape of BWBUG can be expressed by the following formula: where ξ is the direction of airfoil profile curve and η is the spanwise direction. S(η) is the 3 × 3 scale matrix and its diagonal element is S x , S y , 0 . Θ(η) represents the rotation matrix with the twist angle of θ. T(η) represents the position coordinate. Figure 4 shows an example of the left-right symmetrical BWBUG. For ease of description, the left half is taken for display. By changing the parameter values of P i , S x , S y , θ, T, the geometric deformation of the BWBUG can be realized in the optimization design process.

Geometric Constraint Gradient Calculation.
Because the BWBUG has internal loading requirements, it is often necessary to set corresponding thickness constraints at different positions. Suppose that the thickness h between the upper and lower surfaces of the BWBUG is required to be greater than the initial thickness h 0 at position (x, 0, z), the geometric constraint g can be expressed as follows: where (ξ 1 , η 1 ) and (ξ 2 , η 2 ) are the local coordinates corresponding to the upper and lower surfaces of BWBUG at po- According to the expression formula (2) of W(ξ, η), the gradient information of the geometric constraint g relative to the parameterized variables P i , S x , S y , θ, T can be derived as follows:

Multifidelity Surrogate Model.
e lift-drag ratio of BWBUG is the primary objective of shape optimization, and its calculation by high-fidelity CFD model is the most timeconsuming. Considering that the results of the low-fidelity CFD model can reflect the change trend of lift-drag ratio, the nested Kriging model of high-fidelity and low-fidelity models are established to significantly improve the efficiency of building the Kriging model under the condition of achieving the same accuracy.

CFD Model with High and Low Fidelity.
In this section, a multifidelity CFD model composed of the highfidelity model and low-fidelity model is established. Among them, the difference of high and low-fidelity models is mainly reflected in the grid elements and the number of grids. For the grid generation of the left-right symmetrical BWBUG, a half shape is modeled to reduce the computation cost.
e computation domain is set as a box-topology shape and the size of the domain is set to , 10L] (L is the reference length) in length, width, and height direction, respectively. e trimmed mesher and prism layer mesher are used to generate the volume grids in the computational domain. e number of prism layers is set to 5, and the volume growth rate is set to 1.2. In addition, two control volumes are introduced to increase the local grid density near the BWBUG. e mesh number of the high-fidelity model is set about 350 W, and that of the low-fidelity model is set about 70 W. Figures 5 and 6 show the surface grids of the highfidelity and low-fidelity model, respectively. It can be seen that the high-fidelity model has smaller grid elements and a larger number of grids, while the low-fidelity model has larger grid elements and a smaller number of grids.
Both the high and low-fidelity models are solved using the SST k − w turbulence model, and the same boundary conditions and solving parameters are set. e boundary conditions include the following: the left surface of the computation domain is set as the velocity inlet, the right surface is set as the pressure outlet, the back surface is set as the symmetry plane, and the other three remaining surfaces are set as the slip walls. In addition, the BWBUG surface is set to the nonslip wall. e solution parameters include the following: the maximum number of iteration steps is set to 1000 steps, and the convergence residual is set to 1e − 5 . Figures 7 and 8 show the pressure distribution of high and low fidelity models at 1kn speed and 4°angle of attack. From the pressure distribution, it can be seen that the high-fidelity model exhibits smoother pressure distribution at the front edges of the BWBUG than the lowfidelity model.

Nested Kriging Surrogate Model.
e Kriging model [29,30] is an interpolation model. e interpolation result y(x 0 ) at the point x 0 is defined as the linear weighting of the response values y s of the known samples: where w(x 0 ) is the weighting coefficient.    Mathematical Problems in Engineering From formula (5), the lift-drag ratio at the estimated points can be obtained as long as the expression of the weighting coefficient w(x 0 ) is given. For the Kriging model, the unknown function is regarded as a specific realization of the Gaussian static random process, and the optimal weighting coefficient w(x 0 ) is calculated by both meeting the minimum mean square error and the unbiased condition. Based on [30], the analytical expression of w(x 0 ) is given by the following formula: where R is called the correlation matrix, which is composed of correlation function values between all known sample points. r is called the correlation vector, which is composed of correlation function values between unknown points and all known sample points. f is the vector of basis function, which is previously defined. By solving the linear equations (6) and substituting them into formula (5), the estimated value of the Kriging model can be obtained as follows: where β � (F T RF) − 1 (F T R − 1 y s ). It is assumed that the sampling points of high and lowfidelity CFD models are as follows: where the subscripts "1" and "2" represent high-fidelity and low-fidelity CFD models, respectively, and n is often much smaller than m. eir corresponding lift-drag ratio calculation results are as follows: en, the nested Kriging model of high-and low-fidelity CFD models is established by the following steps.
Step 1. Kriging model is established for the sample data of low-fidelity CFD model as follows: where y lf is the low-fidelity Kriging model.
Step 2. e difference between high-fidelity value and low-fidelity Kriging model prediction value is calculated at the sample points of the high-fidelity CFD model, as shown in the following formula: Step 3. e Kriging model for the difference value calculated by formula (11) is established by the following formula: Step 4. By summing the Kriging model established by formulas (10) and (12), the approximate results of the high-fidelity CFD model are obtained as follows: where y hf represents the approximate value of the high-fidelity CFD model.

Shape Optimization Framework for BWBUG.
To perform the shape optimization of BWBUG efficiently, the nested Kriging model is established according to Section 4.2 to replace the time-consuming CFD model to calculate the lift-drag ratio. e suboptimization problems for adding new sample points are solved efficiently by introducing the geometric constraint gradient information calculated in Section 4.1.
e sample points are infilled dynamically, and the nested Kriging model is updated iteratively until the optimization converges. e detailed optimization framework based on the above idea is described as follows: (1) e experimental design methods, such as full factorial design, fractional factorial, orthogonal array, central composite design, or Latin hypercube design [31,32], are used to sample the parameterized design space, and the initial sample points of high-fidelity and low-fidelity models are obtained, respectively. (2) At the corresponding sample points, the high-fidelity and low-fidelity CFD models are solved to calculate the lift-drag ratio of BWBUG, and the nested Kriging model is established according to the method in Section 4.2.2. (3) e criterion for infilling sample point is adopted and the corresponding suboptimization problem is solved to obtain the optimal solution as the new sample point. e common criteria include Minimum of Surrogate Prediction (MSP), Expected Improvement (EI), or Probability of Improvement (PI). [33,34]. Among them, MSP is the simplest, most direct, and the earliest criterion to be used. Its principle is to infill the sample point which can make the objective function of the surrogate model reach the minimum value. For the BWBUG, MSP is selected, and it mainly infills sample points by solving the following suboptimization problems with geometric constraints: where y hf (x) is the predicted value of the nested Kriging model.

Mathematical Problems in Engineering
For the above suboptimization problem, considering that the gradient analytic expression of geometric constraint g(x) is given in Section 4.1.2, gradient optimization algorithms (such as the SQP method) are used to solve the suboptimization problem efficiently. (4) At the new sample point generated by MSP, accurate CFD analysis is carried out to determine whether it converges to the optimal solution. If it converges, the optimization process ends. Otherwise, the results are added to the existing data set and the whole optimization process is repeated.
For a better understanding of the process above, an overview of shape optimization algorithm for BWBUG is shown in Algorithm 1, and the shape optimization framework for BWBUG is described in Figure 9.

Results and Discussion
In this paper, a left-right symmetrical BWBUG is taken as an example to perform the shape optimization design for verifying the effectiveness and efficiency of the proposed surrogate-based optimization method.

Optimization Problem Description.
e initial shape of the left-right symmetrical BWBUG is composed of 11 airfoil sections, and all the airfoil sections are NACA0012 airfoil. In order to reduce the number of parametric variables and the amount of calculation in the optimization process, the left half shape composed of six airfoil sections is selected for optimal design.
Considering that the internal structure and loading equipment of the BWBUG have been determined with the initial shape design, the chord length of each airfoil section and the placement position of each airfoil section are set as constant to ensure that the internal structure and equipment layout remain unchanged in the optimization design. e scale factors in y direction S y and the twist angle θ of six airfoil sections are selected as optimization design variables to improve the performance of BWBUG.
In addition, the thicknesses at 1/4 chord position of 6 airfoil sections are set to be no less than the initial value as geometric constraints to meet the loading requirements of internal equipment. And the interval constraints of optimization variables are also introduced not only to ensure sufficient optimization exploration space but also to avoid a sharp increase in the amount of calculation caused by an excessively large interval.
To sum up, the following optimization problems are obtained: where S yi 0 is the initial value of the variable S yi and h i 0 is the initial thickness value of BWBUG at 1/4 chord position of the i th airfoil section.

Shape Optimization and Results
Analysis. e proposed optimization method is used to solve the optimization problem (15). Firstly, the Latin hypercube method is used to select 30 and 70 sample points as the initial sample points of high-fidelity and low-fidelity CFD models, respectively. Secondly, the corresponding lift-drag ratio is calculated at each sample point, and the nested Kriging model is built according to Section 4.2.2. Finally, the gradient-based optimization algorithm is used to solve the MSP suboptimization problem and the obtained optimal points are infilled dynamically to update the nested Kriging model until the optimization converges. e convergence condition is that the difference between the two optimal sample points is less than 10 − 4 .  Figure 9: e shape optimization framework for BWBUG. e shape optimization was carried out in the workstation, whose processor is Intel ® Core ™ i7-7700K and memory size is 16 GB. e whole optimization converged after running 362 hours. In the process, 23 sample points were infilled dynamically to update the nested Kriging model. A total of 53 high-fidelity CFD analyses and 93 lowfidelity CFD analyses were performed to obtain the optimal shape for BWBUG. Figure 10 shows the convergence history of the whole optimization process. Table 1 lists all the initial values from which the optimization was started, and all the optimal values obtained after the shape optimization of BWBUG. Figure 11 compares the pressure distribution on the upper surface between the initial shape and the optimal shape of BWBUG. Figure 12 compares the pressure distribution on the lower surface of the initial shape and the optimal shape.
It can be seen from Table 1 that the thicknesses at 1/4 chord position of six airfoil sections are all greater than the initial values, which meets the requirements of geometric constraints. Compared with the initial shape of BWBUG, the lift-drag ratio has an increase of 23.9% from 13.24 to 16.41. Figures 11 and 12 give the reason why the lift-drag ratio increases. As Figure 11 shows, the low-pressure area on the upper surface of the optimal shape is larger than that of the upper surface of the initial shape, especially at the leading edge. While Figure 12 shows that the area of the high-pressure zone on the lower surface of the optimal shape is larger than that of the lower surface of the initial shape. Considering the total pressure distribution on the upper and lower surfaces, the optimal shape produces a higher lift force than the initial shape, which indirectly increases the high lift-drag ratio.
Taken together, the proposed optimization method tries to find a better design shape of BWBUG at each iteration ( Figure 10) by dynamically infilling new sample points to update the nested Kriging model of the multifidelity CFD model, and finally, it obtains an optimal shape which has a better hydrodynamic distribution meeting the design requirements than the initial shape (Table 1). It is demonstrated that the proposed optimization is effective in improving the hydrodynamic performance of BWBUG.

Comparison with Traditional SBO Method.
In this section, the traditional SBO method and the proposed optimization method are both applied to optimize the optimization problem (15), and their optimized results are compared with each other to verify the efficiency of the proposed optimization method.
In order to make the optimized results statistically significant, five initial points are selected evenly within the variable bounds to start the optimization, as shown in Table 2. As for the parameter settings of the proposed optimization method, they are the same as those in Section 5.2. As for the parameter settings of the traditional SBO method, the Latin hypercube method is also adopted to select 70 sample points to establish the initial surrogate model, but the high-fidelity CFD model is used to calculate the lift-drag ratio at all sample points. e remaining parameter settings are the same as those in Section 5.2.
In order to shorten the time to complete the whole optimization above, the ten optimizations above were performed on a massively parallel supercomputer and each optimization was carried out concurrently with the same computational resources. Table 3 lists the detailed comparison data between the traditional SBO method and the proposed optimization method.
It can be seen from Table 3 that the average values of L/D obtained by the traditional SBO method and the proposed optimization method are 16.39 and 16.40, respectively, which are roughly the same as a relative error of only 0.06%. And for each case of five initial points, the maximum relative error of L/D is 0.18%. All these results imply that the proposed optimization method can get almost the same optimal objective as the traditional SBO method. Considering that the traditional SBO method has been widely proved to be an effective optimization method [35], therefore the proposed optimization method is also proved to be effective indirectly.
In addition, Table 3 shows that the average time cost by the traditional SBO method and the proposed optimization method is 261.0 hours and 185.7 hours, respectively. Compared with the traditional SBO method, the average time cost by the proposed optimization method is reduced by 28.9%, which demonstrates that the proposed optimization method has higher optimization efficiency than the traditional SBO method. After the in-depth analysis, the calls to the low-fidelity CFD model and high-fidelity CFD model (Table 3) by the two methods reveal the reasons why the proposed optimization method is more efficient. For the traditional SBO method, just the high-fidelity CFD model is called for the calculation of L/D and the average calls are 98. While for the proposed optimization method, the average calls to the low-fidelity CFD model are 91 and the average calls to the high-fidelity CFD model are 51. e proposed optimization method calls 91 more low-fidelity CFD models and 47 less high-fidelity CFD models than the traditional SBO method. What is more, the computational time of the high-fidelity CFD model is almost 4-6 times that of the low-fidelity CFD model for the optimization problem (15). erefore, the proposed optimization method mainly improves the optimization efficiency by using the low-fidelity CFD model to reduce the computational time of high-fidelity model, which agrees well with the hypotheses before.
However, it should be noted in Table 3 that the optimal values of L/D obtained from different initial points are different with a small deviation for both the proposed optimization method and the traditional SBO method. It implies that both methods cannot strictly converge to the global optimality. Although it is difficult and uneconomical for surrogate-based optimization methods to find a strict global optimal solution to complex engineering problems, it should be fully considered in our future work.

Inputs:
n, m: number of initial samples for high-and low-fidelity CFD models [x l , x u ]: lower and upper bounds τ: Minimum Expected Improvement Outputs: x * : best point identified y * : corresponding function value Calculate by multifidelity CFD models While k < k max and (y hf − y * )/y * < τ do y lf � kriging(S 2 , Y 2 ) Build Kriging for high-fidelity model Calculate the difference values d � kriging(S 1 , d) Build Kriging for difference values y hf � y lf + d Build the nested Kriging model Solve the MSP problem y hf � High_CFD(x * ), y lf � Low_CFD(x * ) Evaluate at predicted optimum y * � min(y hf , Y 1 ) Update best point if necessary Add new point to training data

Conclusions
When the traditional SBO method is used for the shape optimization design of BWBUG, there are some problems such as time consumption of high-fidelity CFD analysis and large amount of calculation of geometric constraints in the suboptimization process. In this paper, a shape parametric model of the BWBUG is established, the analytic expression of geometric constraint gradient is derived, and the multifidelity CFD surrogate model is constructed. Finally, an efficient surrogate-based optimization method based on the multifidelity model and geometric constraint gradient information is proposed, which reduces the calculation times of the high-fidelity model and geometric constraints in the shape optimization process, and greatly improves the optimization efficiency.
e shape of one left-right symmetrical BWBUG was optimized to verify the effectiveness and efficiency of the proposed optimization method. e results show that the lift-drag ratio of the optimal shape obtained by the proposed method is 23.9% higher than that of the initial shape, and the optimization efficiency of the proposed method is 28.9% higher than that of the traditional SBO method, which proves that the proposed method has a good application prospect in the shape optimization design of the BWBUG.
In addition, the main theories of the proposed optimization method, such as the parametric modeling method, the gradient calculation method of geometric constraints, and the multifidelity surrogate model method, are not just suitable for the shape of BWBUG, but also for that of the aircraft, ships, cars, and so on. erefore, the proposed optimization method has a broad applications prospect in Table 2: Initial points from which the optimization is started.

Case
S y1 /S y10 S y2 /S y20 S y3 /S y30 S y4 /S y40 S y5 /S y50 S y6 /S y60 θ 1 θ 2 θ 3 θ 4 θ 5    Data Availability e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that they have no conflicts of interest.