Zhonghua HAN, Chenzhou XU, Liang ZHANG, Yu ZHANG, Keshi ZHANG,Wenping SONG
School of Aeronautics, Northwestern Polytechnical University, Xi’an 710072, China
KEYWORDS
Abstract A variable-fidelity method can remarkably improve the efficiency of a design optimization based on a high-fidelity and expensive numerical simulation, with assistance of lower-fidelity and cheaper simulation(s). However, most existing works only incorporate‘‘two” levels of fidelity,and thus efficiency improvement is very limited.In order to reduce the number of high-fidelity simulations as many as possible,there is a strong need to extend it to three or more fidelities.This article proposes a novel variable-fidelity optimization approach with application to aerodynamic design. Its key ingredient is the theory and algorithm of a Multi-level Hierarchical Kriging(MHK), which is referred to as a surrogate model that can incorporate simulation data with arbitrary levels of fidelity.The high-fidelity model is defined as a CFD simulation using a fine grid and the lower-fidelity models are defined as the same CFD model but with coarser grids, which are determined through a grid convergence study. First, sampling shapes are selected for each level of fidelity via technique of Design of Experiments (DoE). Then, CFD simulations are conducted and the output data of varying fidelity is used to build initial MHK models for objective (e.g.CD) and constraint (e.g. CL, Cm) functions. Next, new samples are selected through infillsampling criteria and the surrogate models are repetitively updated until a global optimum is found.The proposed method is validated by analytical test cases and applied to aerodynamic shape optimization of a NACA0012 airfoil and an ONERA M6 wing in transonic flows. The results confirm that the proposed method can significantly improve the optimization efficiency and apparently outperforms the existing single-fidelity or two-level-fidelity method.
1. Introduction
During the past two decades,aerodynamic shape optimization based on high-fidelity CFD simulations is playing an increasingly important role in the area of aircraft design. However,it is still suffering from the difficulty associated with large computational cost,which can be prohibitive when a large number of computationally expensive CFD simulations are conducted.Therefore, it is of great significance to develop more efficient aerodynamic shape optimization methods that could reach an optimal design with high-fidelity and expensive CFD simulations as few as possible.
The existing aerodynamic shape optimization methods can be classified into three categories: (A) gradient-based method.It is very efficient when the gradients are computed by an adjoint approach proposed by Jameson.1The drawback is that the solution optimality can be sensitive to the initial guesses and it can become trapped into a local minimum.2(B)gradient-free heuristic method. Among the available methods, metaheuristic optimization algorithms such as Genetic Algorithms (GA),Simulated Annealing (SA), or Particle Swarm Algorithm(PSA)feature good capability of global optimization.However,when this type of algorithms is applied to an aerodynamic shape optimization, it usually requires thousands of CFD simulations3or even more,4and the overall computational cost can easily exceed the available computational budget.This situation could become even worse when dealing with complex aircraft configurations parameterized with numerous design variables.Therefore, their applications were usually limited to either 2D aerodynamic configurations3or 3D configurations using lowfidelity and fast CFD simulation methods.4(C) Surrogate-Based Optimization (SBO).5–13SBO represents a type of optimization algorithms that makes use of cheap-to-evaluate surrogate models to approximate expensive objective and constraint functions, guiding the addition and evaluation of new sample points towards global optimum.14SBO enables us to find the global optimum within much less expensive evaluations, and thus it can be viewed as a generalization of the Efficient Global Optimization (EGO) proposed by Jones et al.14Despite the growing popularity and continuous progress15–18of SBO, it is still suffering from the problem associated with the ‘‘curse of dimensionality”, which means that with increase of number of design variables,the computational cost of optimization grows rapidly and soon becomes prohibitive.To tackle this problem,a number of researchers have devoted themselves to the development of more efficient global optimization and two kinds of methods have been concerned. One is the gradient-enhanced surrogate model19,20with cheap gradients computed by an adjoint method,1and the other is the Variable-Fidelity Model(VFM) (or multi-fidelity model)21,22that uses auxiliary lowerfidelity, cheaper function to assist the prediction of a highfidelity,expensive function.Here we are mainly concerned with the VFM method,and the issue related to the use of cheap gradients is beyond the scope of this article.
VFM can be used to dramatically improve the efficiency of a design optimization driven by a high-fidelity and expensive simulation, with assistance of low-fidelity and cheaper simulation(s). The basic assumption of a VFM is that the lowerfidelity models could correctly capture the underlying variation trend of a high-fidelity objective (or constraint) function of interest, although they incorporate less physics and are thus less accurate. With this assumption, only a few high-fidelity samples are required to correct the ‘‘trend” underlying many low-fidelity samples and the resulting surrogate model is generally more accurate than the one using the high-fidelity samples alone. Fortunately, this kind of lower-fidelity models, such as generated by applying the same physical model on coarser computational grids, is always available for many engineering design problems, which inspires the engineering design community to put continuous effort studying on more efficient VFM and optimization methods.
Currently,there exist three types of techniques for building a VFM. The first type is a correction-based method that corrects low-fidelity model to approximate high-fidelity function23,24using a bridge or scaling function. The correction can be multiplicative,25,26additive27,28or hybrid.20,29The second type is space mapping.30,31The design space of a lowfidelity function is distorted to make its optimal point match that of a high-fidelity function.32The third type is variablefidelity kriging, such as co*kriging21,33,34or Hierarchical Kriging(HK).22co*kriging was originally proposed in geostatistical community35and then extended to deterministic computer experiments by Kennedy and O’Hagan.33A co*kriging introduces the assistance from a cheap low-fidelity function by constructing the cross covariance between low- and high-fidelity functions. In contrast, the HK model proposed by Han,22the first author of this article, directly takes the surrogate model built from a cheap low-fidelity function as the model trend of the kriging for an expensive high-fidelity function,and the difficulty associated with constructing the cross covariance of a co*kriging is completely avoided. It turns out that a HK model is as simple and robust as a correction-based method,and it is as accurate as a co*kriging method.Moreover,it provides more reasonable Mean-Squared-Error (MSE) estimation than any of the existing kriging and co*kriging methods.
Recently, the HK model22and its variants36,37have received much attention in the area of engineering design and optimization, such as uncertainty analysis in CFD,38multi-objective optimization,39aerodynamic shape optimization40–42and structural optimization,43as well as other areas such as film hole array optimization,44aero-servo-elastic analysis and extreme loads prediction of wind turbines.45It should be highlighted that in the well-known report by NASA,‘‘CFD Vision 2030 Study: A Path to Revolutionary Computational Aerosciences”,46HK is cited as a representative surrogate modeling method that can be used to reduce the risk of aerospace system design.
As a summary of literature review, it is observed that despite the growing popularity of HK model and other VFM methods,most existing research works only incorporate‘‘two”levels of fidelity,which restricts the efficiency improvement.In order to reduce the number of high-fidelity simulations as many as possible, there is a strong need for the community to extend it to multi-level fidelities(three or more levels),which motivates the research of this article.
The main objective of this article is to propose a novel variable-fidelity optimization method towards the application to aerodynamic design. The key point is to develop the theory and algorithm of a Multi-level Hierarchical Kriging (MHK) model, which is as an extension of the existing HK to a model that can incorporate simulation data with arbitrary levels of fidelity. The MHK is to be combined with an Expected Improvement (EI) method to obtain a more efficient variable-fidelity optimization method. The proposed method will be validated by analytical test cases and applied to benchmark aerodynamic shape optimizations of transonic airfoil and wing, in which the high-fidelity model is defined as a CFD model using a fine grid and lower-fidelity models are defined as the same CFD model but using coarser grids.
This article continues in Section 2 for the formulation and algorithm of the proposed MHK model. In Section 3, the variable-fidelity optimization method combining MHK model with an EI method will be described and analytical test cases are to be used to validate the proposed method. In Section 4,the proposed method is demonstrated by benchmark aerodynamic shape optimizations of a NACA0012 airfoil and an ONERA M6 wing in transonic flows, with number of design variables in the range from 18 to 80. At last, general conclusions will be drawn and the future work beyond the scope of this article will be discussed.
2. Formulation of multi-level hierarchical kriging model
Kriging model is a statistical interpolation method suggested by Krige47in 1951 and mathematically formulated by Matheron48in 1963. Kriging model gained popularity in design and analysis of deterministic computer experiments after the pioneer research work by Sacks et al.49in 1989. In 2012, Sacks’s kriging model was extended to incorporate multi-fidelity data by Han and Goertz,22and a VFM called‘‘HK model” was proposed. A HK model is referred to as a two-level surrogate model in which the prediction of the output of a high-fidelity simulation code (high-fidelity function)is assisted by the output of a lower-fidelity simulation code(low-fidelity function).
In this article, we aim to propose a more efficient and robust VFM technique, which is an extension of our previous work of two-level HK22to a model that can incorporate data with arbitrary (three or more) levels of fidelities. We call it MHK in this article. In what follows, the formulation of MHK will to be presented.
2.1. General assumption and descriptions
For anm-dimensional problem (ormdesign variables), suppose there are ‘‘L” simulation models with varying degree of fidelity and computational expense(e.g.potential theory,Euler equations, and Naiver-Stokes equations, or the same physical model on fine and coarse grids). Here we are concerned with the prediction of an expensive-to-evaluate high-fidelity function (e.g.CL，CD， orCm)
wherey1:Rm→R, with assistance of cheaper-to-evaluate lower-fidelity functions
whereyk:Rm→R.Please note that the index‘‘k”denotes thekth level of fidelity, with ‘‘k=1” and ‘‘k=L” representing the highest and lowest fidelity level, respectively. In other words,with increase of level index‘‘k”,the computational cost is reduced dramatically, although the model fidelity decreases.
Assume that the high-and lower-fidelity functions are sampled atn1，n2， ....，nLsampling sites respectively. In order to reduce high-fidelity evaluations as many as possible,we usually assumen1≪n2≪....≪nL. Taken thekth level of fidelity as an example, the function is observed at sites
with corresponding responses
wherenkdenotes the number of sampling sites for thekth level of fidelity. The pair （Sk，yS，k） denotes the sampled data of thekth level of fidelity in the vector space.
With the aforementioned descriptions and assumptions,our objective here is to derive a novel surrogate model for predicting the output of a high-fidelity simulation code at any untriedx∈Rm(that is to estimatey1（x）)based on the sampled datasets （S1，yS，1）（S2，yS，2），（SL，yS，L）， in an attempt to achieve the desired accuracy with the least possible number of highfidelity computations.
2.2. Mathematical formulation of MHK model
Different from the conventional kriging by Sacks et al.,49we assume that the stationary random process corresponding to thekth level of fidelity is of the form
Solving the constrained minimization problem by introducing a Lagrange multiplier μk, the weight coefficientswkin Eq.(7) can be found by solving the following linear equations:
where
Inverting the partitioned matrix, the predictor for thekth level of fidelity can be writing as
where the notation ‘‘Krig” is defined as
2.3. Strategy of Fitting an MHK model
Fig.1 Schematics of building a MHK model in a recursive manner.
The logic of such a recursive modeling process is sketched in Fig.1.
2.4. Correlation model and hyperparameter tuning
The construction of correlation matrixRand correlation vectorrfor kriging models of each level of fidelity requires the calculation of correlation functionR. Taken thekth level of fidelity as an example,here we focus on a family of correlation functions that are of the form
where ‘‘scf” denotes the spatial correlation function that only depends on the Euclidean distance between two sites (xandx′) and the hyperparameters θk.
Thus far,several types of correlation functions can be used,such as ‘‘Gaussian function”, and ‘‘cubic spline function” etc.Since the Gaussian function can lead to a correlation matrix with larger condition number,50the cubic spline function51is employed in this article:
After the correlation model is defined, we will focus on the method of tuning the hyperparameters forkth level of fidelity hereafter.
Assuming that the sampled data are distributed according to a Gaussian process,the responses at sampling sites are considered to be correlated random functions with the corresponding joint likelihood:
The optimal estimates of the scaling factor and the process variance
are obtained analytically yet depend on the unknown hyperparameters θk= [θ1，k，θ2，k，...，θm，k]. Substituting them into associated Eq. (22) and taking logarithm, we are left with maximizing the concentrated logarithmic likelihood function(neglecting constant terms):
As there is no closed form solution for θk,it has to be found by numerical optimization, which is of the form
In this article, we use an improved version of Hooke &Jeeves pattern search method (by using multi-starts search and a trust-region method)to solve the preceding optimization problem.15,19In addition, we normalize the design variablesxto the range[0,1.0]and limit the searching of optimal θkto the range [0.001, 1.0] for the cubic spline function, according to our intensive numerical experiments, to make the model more robust.
The tuning of hyperparameters for other levels of fidelity is essentially similar to above procedure. By using the strategy given by Section 2.3, an MHK model is fitted.
3. Variable-fidelity optimization based on MHK surrogate models and EI method
We are concerned with solving the following constrained optimization problem based on a high-fidelity expensive numerical analysis with assistance of low-fidelity but cheaper numerical analyses
wherey1（x） andg1，i（x） denote the objective and constraint functions, respectively, which are evaluated by a high-fidelity and expensive numerical analysis code;yk（x） andgk，i（x） are the lower-fidelity objective and constraint functions, respectively, which are evaluated by cheaper numerical analysis codes;NCis the number of constraint functions;Lis total number of fidelities;xupandxloware the upper and lower bounds of the design variablesx, respectively. Note that here we are mainly concerned with a single-objective optimization,and it can be readily extended to a multi-objective optimization. But this is beyond the scope of this article.
3.1. Optimization framework and algorithm
The optimization problem is solved based on MHK surrogate models, which are built through a few expensive high-fidelity samples and many cheaper low-fidelity samples. The MHK models are repetitively updated by the adaptive addition of expensive high-fidelity sample points suggested by maximizing the EI function (Jones et al.14in 1998), until the global optimum is reached.The framework of this MHK-based optimization is sketched in Fig.2. The basic steps are as follows:
(1) ‘‘L”sets of sample points are generated by using a DoE method and they are evaluated by the low- and highfidelity numerical analyses, respectively.
(2) Starting from the lowest level of fidelity,we sequentially build the kriging model for the objective function for each level of fidelity, until the initial MHK model for highest level of fidelity is built. The MHK model(s) for constraint function(s) is built in a similar manner. See Section 2 for the details.
(3) New sample point(s) is selected by maximizing the constrained EI function through a combination of GA,Hooke & Jeeves, and BFGS optimization methods,and then evaluated by the expensive, high-fidelity numerical analyses. When the BFGS method is used,the gradients of an EI function are computed by a central difference based on the MHK prediction. Please note that only high-fidelity sample points are chosen to update the surrogate models.
Fig.2 Flowchart of variable-fidelity optimization based on proposed MHK surrogate models and EI method.
(4) The newly selected sample point and its functional responses are augmented to the sampled dataset, and then the MHK surrogate models are updated.
(5) The steps (2)–(4) are repeated until termination condition is satisfied.15,52
For the above procedure (also see Algorithm 1), the infillsampling criterion is critical since it has significant impact on the convergence efficiency and optimal results. In this article,the EI criterion proposed by Jones et al.14is employed.Assume that the prediction of an MHK model at any untried sitexobeys a normal distribution ^Y（x）N^y1（x），s2（x）[ ],with the mean being the surrogate prediction ^y1（x） and the standard deviations（x） being its Root Mean-Squared Error (RMSE).Then the statistical improvement at any untried location w.r.t. the best high-fidelity objective function observed so faryminis defined as:
According to the derivation of Jones et al.,14the EI function can be written as
where Φ and φ are the cumulative distribution function and probability density function of standard normal distribution,respectively.
And the constrained EI function can be given by:
If there areNCconstraints in total, we should builtNCMHK models for every constraint function, and the resulting constrained EI function is
Then the following unconstrained sub-optimization problem is formulated as
A hybrid method of combing GA, Hooke & Jeeves pattern searchandBFGSgradient-basedmethodisusedtosolvethe above sub-optimization problem to suggest new sample point,15,19which is to be evaluated by high-fidelity numerical analysis code again.
Algorithm 1(ProcedureofaMHK-basedvariable-fidelity optimization).
Please note that the MHK model and the aforementioned EI method have been integrated to an in-house optimization code called ‘‘SurroOpt”.15,17–19,52,53In the next subsections, SurroOpt will be employed to run the test cases of a two-dimensional sixhump camelback unconstrained function and a seven-dimensional G9 constrained function, in order to verify and validate the proposedmethod.Thenwewilldemonstratetheproposedmethodwith twoexamplesofbenchmarkaerodynamic shape optimization.Inall of these examples, three different optimization methods will be compared: kriging-based optimization (single-fidelity method);HK-based optimization (two-level variable-fidelity method);MHK-based optimization (three- or more-level variable-fidelity method).To conduct a fair comparison,all the settings excluding the choice of surrogate models will be kept the same.
3.2.Unconstrained global optimization of a six-hump camelback function
A six-hump camelback function is taken as the first test case.The mathematical model of optimization is:
First, a single-fidelity optimization using the in-house‘‘SurroOpt” code with kriging model is conducted to verify the correctness of our code. Three initial samples are chosen by the Latin Hypercube Sampling (LHS) method to build an initial kriging model for the objective function. Then the EI infill-sampling criterion is used to repetitively select new samples to update the kriging model, until the global optimum is found. The convergence history is sketched in Fig.4, whose vertical axis is the relative error between the current objective function and the theoretical global optimum.The triangle symbol denotes all the sample points, and the solid line is plotted in a manner that shows the best observed functional value thus far versus the number of high-fidelity evaluations.The relative error is dropped below 10-7(the black dashed line) after 87 evaluations,which indicates that the optimum is very accurate and the code is very efficient.
Next,the kriging-,HK-and MHK-based methods are compared to verify the effectiveness of the proposed method. The convergence histories are sketched in Fig.5.Note that the convergence histories are plotted in a manner to only show the best observed high-fidelity functional value versus the number of evaluations,for the convenience of comparison.We assume that the optimal value is found when the relative error drops blow 10-7. As one can see that, compared with the krigingbased method(black solid line),the VFM methods remarkably speed up the optimization (reduce the number of high-fidelity evaluations), and the MHK-based method significantly outperforms the two-level HK method(blue dashed line).The reason can be explained by the fact that many cheap lower-fidelity samples have been used to explore the landscape of the objective function throughout the design space, and only very few high-fidelity samples are required to exploit the global optimum. It should be highlighted that a 3-level MHK is sufficiently efficient, and the introduction of more fidelities (4-and 5-level MHK) does not make remarkable improvement for current test case.
Fig.3 Contour of a six-hump camelback test function.
Fig.4 Convergence history of optimizing a six-hump camelback function using a kriging-based single-fidelity optimization.
3.3. Constrained global optimization test case of a G9 function
We further verify our method by using a test case of constrained seven-dimensional optimization named G9. This is a benchmark testing case intensively used for testing global optimization algorithms by the optimization community.The optimization problem is formulated as
with assistance of
This problem has four nonlinear constraints, with two of them (g1，1（x） andg4，1（x）) being active constraints. The global minimum observed by the community so far is located atx*=（2.330499，1.951372，-0.4775414，4.365726，-0.6244870-，1.038161，1.594227）, wheref（x*）=680.6300573. The lowerfidelity functions for both the objective and constraint functions are defined by adding a small offset to the design variables.
Similarly, a kriging-based single-fidelity optimization is conducted first to validate the optimization code for a constrained global optimization. The convergence history of optimization is sketched in Fig.6.It is shown that the optimum is found after 347 high-fidelity evaluations. Since this is actually a very difficult case, we consider that the optimum is found when the relative error is reduced to 0.2% (the black dashed line in Fig.6).
Fig.5 Comparison of convergence histories of kriging-, HKand MHK-based optimizations for six-hump camelback function.
Then,different methods are used to solve this optimization problem,and the convergence histories are compared in Fig.7.It is shown that the proposed MHK-based method can greatly reduce the number of high-fidelity evaluations and thus the overall optimization efficiency can be dramatically improved,if the functional evaluation is very expensive. Again,it is observed that three levels of fidelity would be the best-practice setting for this case.
Table 1 summarizes the results of these analytical test cases,including the number of initial samples, the number of highfidelity evaluations, as well as the optimum. One can see that by employing a variable-fidelity approach, the number of high-fidelity evaluations (also see Figs. 5 and 7) are significantly reduced. This enables us to greatly save computational cost in a practical problem, where the evaluation of highfidelity function is very expensive. The prediction accuracy of the different surrogate models at the optimum is also presented,which shows that all the surrogate models are accurate enough.Finally,it is found that a 3-level MHK seems to be the best practice.
Fig.6 Convergence history of optimizing a G9 function by kriging -based single-fidelity method.
Fig.7 Comparison of convergence histories of kriging-, HKand MHK-based optimizations for G9 function.
4. Examples of aerodynamic shape optimization
Two representative examples will be used to demonstrate the application of the proposed method to efficient benchmark aerodynamic shape optimizations. One is the drag minimization of a 2D airfoil, and the other is the aerodynamic shape optimization of a 3D wing.
The choice of CFD models of different fidelity is of great significance for a variable-fidelity aerodynamic shape optimization. Generally the low-fidelity CFD model can be determined through the following three ways: (A) simplified physical model, like Euler equations w.r.t. NS equations; (B)coarse discretization, like coarser grid w.r.t. fine grid; (C)relaxed convergence criteria, like insufficiently convergent results w.r.t. fully convergent results. No matter which way we choose, the low fidelity CFD model should be cheap to evaluate and features the similar variation trend with the high-fidelity model.
Here,we take RAE2822 airfoil test case as an example.The flow condition isMa=0.734,Re=6.5×106. An in-house flow solver PMNS2D54is employed to perform flow simulation. The variation of force coefficients computed with different CFD models is depicted in Fig.8.It is shown that the CFD results obtained by applying the same governing equations on grids of varying resolution are essentially very close to each other, while the CFD results obtained applying different governing equations are dramatically different, and particularly the variation trend of moment coefficients is even inconsistent.Apart from this, the pressure distributions are compared in Fig.9, which shows that the results from the same governing equations on different grids are in reasonably good agreement with each other.However,the results from different governing equations on the grid of same size significantly differ from each other.
According to the above investigation, it is obvious that using the same governing equations on grids of varying resolution is better suited for building a variable-fidelity model,compared with those using different governing equations.Therefore, in this article the lower-fidelity models are defined by the same CFD model on coarser grids.
4.1.Drag minimization of a NACA0012 airfoil in transonic flow
4.1.1. Problem statement
The first test case is a benchmark case defined by the AIAA Aerodynamic Design Optimization Discussion Group(ADODG).55,56It is the drag minimization of a NACA0012 airfoil subject to a full thickness constraint. According to the problem definition of ADODG, the airfoil is optimized at a freestream Mach number 0.85 and zero angle of attack in inviscid flow. The optimization problem is given by:
wherey≥ybaselinemeans that the local thickness should be always larger than that of the baseline, along the airfoil from the leading edge to the trailing edge.
The NACA0012 airfoil is slightly modified and features zero thickness trailing edge:
A sixteen-order CST method57(proposed by Kulfan) is used to parameterize the airfoil. Note that we force the variables for the upper surface are identical to their counterpartof the lower surface, which ensures that the airfoil is always symmetric and the lift is always zero at aero angle of attack.As a result, the number of design variables is 17 in total.The design space is defined by expanding the initial CST coefficients by 1.5 times and narrowing it by half (see in Fig.10):
Table 1 Summary of optimization results for two analytical test cases.
Fig.8 Comparison of variation of force coefficients w.r.t. angle of attack by using different CFD models for RAE2822 airfoil.
wherexbasedenotes the baseline shape.
Fig.9 Comparison of computed pressure coefficient Cp distributions by different CFD models for RAE2822 airfoil.
4.1.2. Determination of high- and low-fidelity CFD models
The high- and low-fidelity analyses are defined as CFD simulations solving the same governing equations(Euler equations)but using computational meshes of varying degree of resolutions. A grid convergence study is carried out to determine the L1, L2 and L3 grids for high-, medium- and low-fidelity CFD simulations, respectively. The result is shown in Fig.11. The black line with triangle symbols shows the drag coefficient versus the number of grid cells and the blue line with square symbols is for the computational cost. From Fig.11, we decided to use the grid L1 with 131,072 cells for the high-fidelity CFD model, and to use the L2 and L3 grids,with 32,768 cells and 8192 cells, for medium- and low-fidelity CFD models, respectively (see Fig.12). Please note all the CFD simulations in this case,either using a fine or coarse grid,are conducted by using an in-house code called ‘‘PMNS2D”,that a single CFD simulation on a modern personal computer using these grids takes around 56 min, 3.5 min and 8.5 s,respectively.
Fig.10 Geometry of baseline airfoil and definition of design space.
Fig.11 Grid convergence study for NACA0012 airfoil.
4.1.3. Results of optimization
For the kriging-based single-fidelity optimization,only 5 highfidelity samples evaluated by using the L1 grid are used to build initial surrogate models. For the HK-based variablefidelity method, additional 100 medium-fidelity samples computed by the L2 grid are used, besides the 5 high-fidelity samples; for MHK-based method, additional 2000 low-fidelity samples computed by L3 grid are employed, besides the 5 high-fidelity and 100 medium-fidelity samples. The convergence histories are sketched in Fig.13. It is obvious that both the variable-fidelity methods, either HK- or MHK-based method,perform much better than the single-fidelity optimization, and the MHK-based method (three-level) is remarkably more efficient than the HK-based method (two-level). More detailed results are listed in Table 2.One can see that although building a multi-fidelity surrogate model is more expensive,much less updating cycles are required to find the optimum.For a kriging-based method, the drag coefficient is reduced by 43.6% when 172 high-fidelity samples are added. In contrast, only 33 high-fidelity samples are needed for the HKbased method to achieve the same level of drag reduction.Furthermore, only 9 high-fidelity samples are required for the MHK-based method to achieve the same level of drag reduction.Please note that in Table 2,the results of low-fidelity optimization plus high-fidelity evaluation is also presented.40 lowfidelity samples evaluated by using the L3 grid are used to build initial surrogate model and afterwards,the optimal shape obtained by the coarsest grid is re-evaluated with high-fidelity analysis(L1 grid).The result shows that the drag coefficient is 2.7 counts larger than that of our MHK method, which confirms the effectiveness of our MHK method. The comparison of geometric shapes and pressure coefficient distributions of the baseline and the optimal airfoils are sketched in Figs. 14 and 15, respectively. One can see that both the leading- and trailing edges of optimal airfoils are much thicker than that of the baseline airfoil, and the shock waves are largely weakened. Besides, it is observed that the optimal geometries and pressure distributions obtained from different optimization methods are essentially close to each other. In Fig.15, we notice that the final design result still has a strong shock wave.Although we have obtained the best result in current design space, the shape of leading and trailing edges has reached the boundary of the design and better design may be located outside the design space. In Ref. 55 we presented a multiround optimization method, which can adaptively change the design space during the optimization. However, it is still reasonable to use a fixed design space here to fascinate a fair comparison.
Fig.12 Sketch of O-type L1,L2 and L3 grids for high-,mediumand low-fidelity CFD simulations (NACA0012 airfoil test case).
Fig.13 Comparison of convergence histories of optimizing NACA0012 airfoil using kriging-,HK-and MHK-based methods.
Table 2 Optimization results for NACA0012 case.
Fig.14 Comparison of baseline NACA0012 airfoil and optimal airfoils.
Fig.15 Comparison of surface pressure distributions of baseline NACA0012 airfoil and optimal airfoils.
4.2. Drag minimization of ONERA M6 wing in transonic flow
4.2.1. Problem statement
Fig.16 Parameterization of ONERA M6 wing with 5 control sections (80 design variables in total).
Fig.17 Determination of design space for each control section of ONERA M6 wing.
Fig.18 Grid convergence study for ONERA M6 wing.
Fig.19 Sketch of C–H type grids for high-, medium- and low-fidelity CFD simulations (ONERA M6 wing test case).
Fig.20 Comparison of computed pressure distributions with reference data at four span-wise locations.
The objective of this test case is to minimize the total drag of an ONERA M6 wing in inviscid transonic flow,subject to lift,moment and thickness constraints. The baseline shape is shown in Fig.16.We also utilize the CST method56to parameterize 5 control sections(see Fig.16)and the design space for each section is specified by expanding and shrinking it by 25%(see in Fig.17):
Note that each control section is parameterized by an eightorder CST method (18 design variables), which results in 80 design variables in total. This is very demanding for a surrogate-based optimization, since it is a high-dimensional aerodynamic optimization problem with large design space.The drag minimization is formulated as
Fig.21 Comparison of convergence histories of optimizing ONERA M6 wing using kriging-,HK-and MHK-based methods.
The drag is evaluated by a CFD simulation using the fine grid (L1) at a Mach number 0.8395 and an angle of attack 3.06°. To reduce the number of expensive CFD simulations using L1 grid, two coarse grids (L2, L3) are used to assist the optimization. Although more levels of grids can be used,three levels would be sufficient according to our experience.
4.2.2. Determination of high- and low-fidelity CFD models
Similar to the previous case, the high- and low-fidelity models are defined as CFD simulations solving the same governing equations (Euler equations) but with fine and coarser grids,which are also determined by a grid convergence study, as shown in Fig.18.Roughly 2.5 million cells are needed according to a convergence criterion. Then the L1 grid with 2.5 million cells is employed for the high-fidelity CFD model,and the L2 grid with 479,232 cells and L3 grid with 59,904 cells are selected for medium- and low-fidelity CFD simulations,respectively(see Fig.19).The computed pressure distributions of different fidelities at four span-wise locations of M6 wing are compared with reference data in Fig.20,which is obtained by solving Euler equations on a finer gird with 3.65 million cells. It is shown that with increase of grid size, the pressure distributions gradually get closer to the reference data, and all results share the similar variation trend. Please note that all the CFD simulations in this case, either using a fine or coarse grid, are conducted by using an in-house code called‘‘PMNS3D”, and that a single CFD simulation on a modern personal computer for using these grids takes about 192 min,18 min and 1.5 min, respectively.
4.2.3. Optimization results
In this case, LHS method is used to select 10, 100, and 1000 initial samples for high-,medium-and low-fidelity CFD simulations using the L1,L2,and L3 grids,respectively.Initial kriging model is built based on the high-fidelity samples alone,initial HK model is built through on high- and mediumfidelity samples and initial MHK model is built through all the samples (high-, medium- and low-fidelities). We attempt to use this manner to conduct a fair comparison of different methods.Fig.21 shows the comparison of the convergence histories of the drag coefficient. It is found that the MHK-based method significantly outperforms the kriging-based and HKbased methods.Fig.22 sketches the comparison of geometrical shapes and pressure coefficient distributions of the baseline and the optimal wings.As one can see that,the shock is weakened at different span-wise locations, for all the optimal shapes.
Fig.22 Comparison of pressure distributions and geometric shapes of baseline ONREA M6 and optimal wings.
Table 3 Optimization results for ONERA M6 wing test case.
Table 4 Geometric constraints for ONERA M6 wing test case.
A detailed comparison of different optimization methods are given in Table 3. Compared with the kriging-based method, the HK-based method is nearly 3 times faster and the MHK-based method is nearly 5 times faster (76.7% cost saving).Furthermore,the best final result,when using the proposed method, is obtained within only 31 calls of high-fidelity CFD simulation. Without the assistance of lower-fidelity model, this is definitely impossible to achieve, for such an 80-dimensional optimization problem. The prediction accuracy of surrogate models at the optimum is also checked and listed in Table 3, which proves that the surrogate models are accurate enough from an engineering designer’s point of view.Again, we conduct low-fidelity optimization using the coarse grid L3 alone and evaluate the final optimal shape by L1 grid.The drag coefficient obtained by this method is worse than that of using MHK method. The thickness constraints are checked in Table 4, which shows that all the constraints are strictly satisfied.
5. Conclusions and outlook
An efficient global optimization method based on a variablefidelity surrogate model called Multi-level Hierarchical Kriging(MHK) was proposed in this article. A self-contained derivation for the formulation of MHK was presented, and a recursive strategy of building an MHK model was proposed. The MHK surrogate model was combined with an expected improvement function and a more efficient method of variable-fidelity global optimization was obtained. The proposed method was firstly validated by two analytical test cases and then demonstrated by benchmark aerodynamic shape optimizations of a NACA0012 airfoil and an ONERA M6 wing in transonic flows, with the number of design variables in the range from 18 to 80. Some conclusions are drawn as follows:
(1) The proposed MHK model is essentially an extension of the HK model by Han et al.22to a model that can incorporate data with arbitrary levels of fidelities.
(2) The analytical examples of unconstrained and constrained global optimizations prove that the proposed method is correct and able to find the global optimum.
(3) An MHK-based optimization significantly outperforms the kriging-based single-fidelity optimization as well as the HK-based two-level optimization, with respect to number of high-fidelity evaluations,for all the test cases presented in this article.
(4) An MHK model with‘‘three”levels of fidelities seems to be the best practice. From current test cases, introduction of more levels of fidelity may not dramatically further improve the efficiency.
From an academic point of view, we were focusing on the development and validation of a novel approach for more efficient variable-fidelity optimization. Although the wing is optimized by using an Euler solver, it is readily extended to use Navier-Stokes flow solver and it will also be applied to more complex configurations. To achieve a least possible computation cost, the proposed method will be further improved by using more sophisticated infill-sampling method, such as variable-fidelity EI method40and parallel infill-sampling method. Apart from this, the cheap gradients obtained by using adjoint method could be used to improve the efficiency for a high-dimensional problem.These will be our future work beyond the scope of this article.
Acknowledgments
This research was sponsored by the National Natural Science Foundation of China (Nos.11772261 and 11972305), Aeronautical Science Foundation of China (No. 2016ZA53011) and Foundation of National Key Laboratory (No. JCKYS2019607005).
