An efficient geometric constraint handling method for surrogate-based aerodynamic shape optimization

Handling a large number of geometric constraints brings a big challenge to the surrogate-based aerodynamic shape optimization (ASO) driven by computational fluid dynamics (CFD). It is not feasible to calculate the geometric constraint functions directly during the sub-optimization of a surrogate-based optimization, as the geometric constraint functions are to be evaluated thousands of times for each updating cycle and the total cost of a number of cycles can be prohibitive. This article proposes an efficient method of handling geometric constraints within the framework of a surrogate-based optimization to address this problem. The core idea is to use the Kreisselmeier-Steinhauser (KS) method to aggregate all geometric constraints into one that can be approximated by a cheap surrogate model, in order to avoid the large computational cost associated with tremendous calculation of geometric constraint values. The proposed method is verified by an analytical test case. Then, the proposed method is demonstrated and compared with the methods of building surrogate models of all geometric constraints and calculating all geometric constraints directly during each sub-optimization by drag minimizations of NACA0012 air foil and ONERA M6 wing in transonic flows. To investigate the ability of the proposed method for handling various geometric constraints, drag minimization of CRM wing in viscous transonic flow driven by CFD is performed. Results show that the proposed method can dramatically improve the optimization efficiency of ASO with the number of geometric constraints ranging from 15 to 1429 and the number of types of geometric constraints up to 3, which offers great potential for handling a larger number and more types of geometric constraints.


Introduction
Surrogate-based optimization (SBO) has become one of the most attractive methods for aerodynamic shape CONTACT Zhong-Hua Han hanzh@nwpu.edu.cn;Ke-Shi Zhang zhangkeshi@nwpu.edu.cnoptimization (ASO) driven by computationally expensive computational fluid dynamics (CFD) simulations.In a SBO, the expensive cost functions are approximated by cheap-to-evaluate surrogate models to drive the adding and evaluation of new sample point(s) towards the global optimum (Queipo et al., 2005;Forrester & Keane, 2009;Haftka et al., 2016).A typical process of a SBO includes two stages (Parr et al., 2012;Han, 2016a;Liu et al., 2019): (1) using the sample points generated by a design of experiment (DoE) method, initial surrogate models are constructed for the objective and constraint functions, respectively; (2) new samples are suggested by sub-optimization defined by certain infillsampling criteria, then their responses obtained by a solver (such as CFD solver) are used to update the surrogate models.It employs much less function evaluations to search the global optimum compared with a genetic algorithm, especially when expensive high-fidelity CFD or multidisciplinary analysis are employed (Viana et al., 2014;Han et al., 2013;Lu et al., 2021;Park et al., 2020).
SBO is free of gradient information and could be used to solve the problems with high nonlinear and multimodal design space (Haftka et al., 2016;Han, 2016b).Therefore, SBO gained much attention in the area of engineering design optimization.Despite the popularity of SBO, it is still suffering from the difficulty associated with 'curse of dimensionality', which means that the computational cost grows exponentially as the number of design variables increases.To improve the efficiency of a SBO, some advanced technologies have been proposed, such as evaluating the responses of the samples in parallel (Liu et al., 2017), variable-fidelity surrogate modeling (Courrier et al., 2016;Han et al., 2013;Han et al., 2020;Ha et al., 2014;Jo et al., 2016;Leifsson & Koziel, 2015, 2016;Zhang et al., 2018), surrogate modeling using cheap gradients (Laurenceau et al., 2012;Han et al., 2017;Mortished et al., 2016;Bouhlel & Martins, 2019), and applying dimension reduction in the optimization (Li et al., 2019;Bouhlel et al., 2018).However, handling a large number of geometric constraints presents a big challenge to the surrogate-based aerodynamic shape optimization.First, it is not feasible to calculate the geometric constraint functions directly during the sub-optimization performed by a global optimization algorithm such as genetic algorithm within the framework of a SBO, in which these functions defined by massive surface meshes and parameterization methods such as free form deformation (FFD) are evaluated thousands of times for each updating cycle, the total cost of a number of cycles can be prohibitive.Second, the method of building surrogate models for all geometric constraints during the optimization will suffer from difficulties associated with the prohibitive computational cost caused by training the large number of hyper-parameters in a one by one manner.Finally, it is improper to replace all constraints with the most violated one such as the maximum thickness constraint, since the maximum thickness t max may changes during the optimization, which leads the constraint function to be discontinuous, and may result in an increased number of iterations of the optimization and even premature convergence.At present, a large number of geometric constraints have been widely applied in gradient-based aerodynamic shape optimization of different aerodynamic configurations, such as the one-stage turbine (Xu et al., 2015), regional Jet (Bons et al., 2018), large civil aircraft wing (Lei et al., 2019;Kenway & Martins, 2016;Chen et al., 2016), strut-braced wing (Secco & Martins, 2019), 10 MW wind turbine (Madsen et al., 2019), blended-wing-body aircraft (Reist & Zingg, 2017;Lyu et al., 2017), as well as fairing (Brelje et al., 2020).Therefore, it is necessary to investigate an efficient geometric constraint handling method for SBO to solve ASO problems involving a large number of geometric constraint.
In the field of structural optimization, the KS method proposed by Kreisselmeier and Steinhauser (1979) provided a way to handle large-scale constraints, since it can make good use of the KS function to aggregate a large number of constraints into one.The KS method is developed by Raspanti et al. (2000) to handle inequality constraints in a nonlinear programming problem.Akguen et al. (2001) used it to improve the efficiency of solving the sensitivities of large-scale stress and buckling constraints in gradient-based structural optimizations.Then, the KS method was widely used in the gradientbased structural optimization (Paris et al., 2010;Lee & Martins, 2012) and aerostructural optimization (Martins, 2002;Martins et al., 2002;Hwang et al., 2014).It is found that the KS function with a small aggregation parameter ρ is overly conservative.It becomes more accurate as ρ increases, but may cause numerical difficulties and oscillations in a gradient-based optimization due to the ill-conditioned Hessian matrix.Poon andMartins (2005, 2007) addressed this problem by an adaptive KS method, more accurate and robust gradient-based optimization results were obtained compared with the KS method.The adaptive KS method was further used for solving the gradient-based aerostructural optimization (Kenway et al., 2014;Kenway & Martins, 2014) and aeroservoelastic optimization (Haghighat et al., 2012).However, it is mainly used to solve the numerical difficulties and oscillations in gradient-based optimizations, but cannot provide a higher approximate accuracy to the original constraints than the KS function with a larger value of ρ, theoretically.Different KS functions were introduced into the surrogate-based aerostructural optimization by Zhang et al. (2019).Results show that the KS method can handle large-scale structural stress constraints by aggregating all constraints to a continuous one with accurate approximation of the original constraint bound.However, the KS method has not been applied to handle strict geometric constraints of an aerodynamic shape optimization within the framework of a SBO, which motivates the study of this work.
The objective of this article is to investigate an efficient method to solve CFD-driven aerodynamic shape optimization problems involving a large number of geometric constraints within the framework of a SBO, in order to effectively improve the capability of SBO in engineering design.The core idea is to use the KS method to aggregate all geometric constraints into one that can be approximated by a cheap surrogate model, and to calculate the all geometric constraint values during the sub-optimization by using the surrogate model of the aggregated constraint, which leads SBO to avoid the large computational cost associated with building surrogate models for all geometric constraints and calculating all geometric constraint values tremendously in each iteration of optimization.This is nearly a blank area according to literature survey, which inspires the study of this article.This article continues in section 2 to give a detailed description of flow simulation method used in this article.In section 3, surrogate-based optimization for constrained aerodynamic shape optimization is described.In section 4, the proposed method for handling geometric constraints in SBO is dedicated.In section 6, first, the G9 analytical test case is used to verify the proposed method.Second, the drag minimizations of the NACA0012 airfoil and ONERA M6 wing with a small number of geometric constraints are used to compare different geometric constraint handling methods.Moreover, the ONERA M6 wing optimization example with 1000 geometric constraints is used to further study the ability of handling a large number of geometric constraints by the proposed method.Finally, to investigate the ability of the proposed method for handling various geometric constraints, drag minimization of NACA CRM wing in viscous transonic flow driven by CFD is performed.

Flow simulation method
In this article, an in-house code 'PMNS3D' is used as the CFD solver to perform 2D and 3D CFD simulations.It is a finite-volume, cell-centered multi-block solver for solving the Euler and Reynolds-averaged Navier-Stokes (RANS) equations.The Jameson-Schmidt-Turkel (JST) scheme is used for the spatial discretization.The Lower-Upper Symmetric Gauss Seidel (LU-SGS) scheme is implemented for time advancing with multi-grid technique to accelerate the iteration converging to steady states.The turbulence viscous coefficient is calculated by Spalart-Allmaras (SA) one-equation turbulence model.The in-house CFD code has been validated and used to the optimization of various airfoils (Han et al., 2019;Zhang et al., 2018) and wings (Liu et al., 2017;Han et al., 2020).

Governing equation of flow and discretization scheme
The Navier-Stokes equations in integral form for the control volume V and surface element dS can be defined as where Q = [ρ, ρu, ρv, ρw, ρE] T denotes the conservation variables, ρ a is the density.q = (u, v, w) T is the Cartesian velocity vector.E is the specific total energy.n is the surface normal of control volume boundary ∂Ω.F, F v denote the inviscid and viscous flux vector, respectively.Their detailed form can be found in reference (Han, 2021).For the Euler equations, the viscous flux vector of Eq. ( 1) is neglected.By using the cell-centered finite volume method for each grid cell i, the spatial discretization of Eq. ( 1) can be expressed as where are the residual values of the inviscid flux and the viscous flux respectively, and D i is the artificial dissipation term added to prevent non-physical oscillations.
Implicit solution algorithms have been implemented based on the LU-SGS method.The first-order backward difference is performed on the time derivative term, and the viscous flux is treated explicitly, Eq. ( 2) can be written into the following implicit form where

Rk+1
i can be written as where the superscript k represents the number of the time step.Furthermore, the JST scheme is used for the spatial discretization and the SA turbulence model is used to compute turbulence viscosity coefficient.More details can be found in reference (Han, 2021).

Wall boundary condition
The no-slip boundary condition and the adiabatic wall boundary condition are used as the boundary conditions for RANS-based CFD simulations.Since the wall is stationary, q = (u, v, w) T = 0 at the wall.
For Euler based CFD simulations, as the viscosity is neglected, the normal component of the relative velocity is q b = (u b , v b , w b ) T = 0.The no-slip boundary condition can be defined as The adiabatic wall boundary condition assumes that there is no heat exchange between the wall and the fluid, therefore, the gradient of temperature

Symmetric boundary condition
The symmetry boundary condition defines a mirror face that reflects all the flow distribution to the other side.The are no flows and scalar flux across the boundary.By using this boundary condition, the domain can essentially be halved to reduce the time to achieve a solution.

Far-field boundary condition
In this paper, the local one-dimensional Riemannian invariants are used to deal with the far-field boundary.the local one-dimensional Riemannian invariants are defined as where q n denotes the normal component of the velocity at the far-field boundary, a is the sound velocity.

Framework of surrogate-based constrained optimization
The purpose of this article is to explore an efficient method that can be used to handle a large number of geometric constraints of an aerodynamic shape optimization problem in the framework of a SBO.The generic constrained single-objective optimization problem can be defined as where f (x) is the objective function, x is the vector of design variables, c j (x), j = 1, • • •, N c are the constraint functions, N c denotes the number of the constraints.x low and x up denote the lower and upper boundary of the design space, respectively.An in-house surrogate-based optimizer, SurroOpt, is applied to address the optimization problem mentioned above.It has been tested with numerous benchmark numerical examples and practical engineering problems (Han, 2016b;Liu et al., 2017;Han et al., 2020).Figure 1 sketches the framework of SurroOpt.The optimization process could be described as

Infill-sampling criteria (ISC) and constraint handling methods
As we described in previous section, new feasible samples are selected by ISC.This step is called the suboptimization in SBO, which turns out to be a constrained problem that involves the maximization or minimization of a function defined by the ISC.In SurroOpt, ISC involving EI, MSP, PI, LCB and MSE are combined to generate five new sample points.Note that each infill criterion can only be used to select one point.Since each infill criterion has its advantages and disadvantages, the simultaneous implementation of all criteria can complement each other.According to these ISC, the constraints handling methods can be divided into two categories.

MSP and LCB infill criteria with constraint handling
MSP considers adding a new sample at the current optimum.It is not guaranteed that MSP can converge to the global optimum due to the lack of exploration.The sub-optimization problem of MSP is where ' ∧ ' means the surrogate model of the corresponding item of Eq. ( 8).LCB is developed by taking the uncertainty of surrogate models into account, the following sub-optimization problem would be solved where A denotes a user-defined parameter.ŝ(x)is the estimated standard deviation of the surrogate model prediction.
For these two ISC, in order to address a constrained problem, the global optimization algorithm GA is used.The constrained problem is transformed into an unconstrained one by adding a one pass penalty function to the model prediction (Deb, 2000), where J(x) is the objective function value of the worst feasible solution in the population, ĵ is the objective function of ISC for Eq. ( 9) and Eq. ( 10), max(•) returns the largest value.ĉi (x) is the normalized constraint model.

PI, EI and MSE infill criteria with constraint handling
Using the statistical feature of kriging model, both the exploitation and exploration of surrogate model can be accounted to select new sample points by estimating the regions with a high probability of improvement.The probability of improvement function is defined as Correspondingly, the sub-optimization of PI is given by max P(I(x)) PI does not indicate the magnitude of improvement.It only identifies areas where some improvement may be expected.The magnitude of improvement is exposed in the concept of EI, Thus, the sub-optimization to be solved is max In addition, the new update samples can be added at the location with maximum MSE, to improve the global accuracy of surrogate model, The constrained optimization method using penalty function mentioned in section 3.2.1 can also be used for PI, EI and MSE method.Making use of the statistical feature of kriging model, a probabilistic approach is proposed (Schonlau, 1997).Taking the EI criterion as an example, the probabilistic approach is implemented by using a product of the expected improvement of the objective function and the PoF calculated from the constraint functions.The PoF is calculated in the same manner as the probability of improvement, which is given by where ĉi (x) is the surrogate model of the constraint function, s c,i is the root mean square error (RMSE) of ĉi (x).
Then a new sample can be obtained by solving the following optimization problem: where J(x) is P(I(x)), E(I(x)), MSE( f (x)) respectively for PI, EI and MSE infill criteria.Parr et al. (2010) found that the probabilistic method of constraint handling performs better than the penalty method on a number of problems.

Methods of handling geometric constraints
In ASO, a large number of geometric constraints are needed due to the structure or assembly requirements, e.g.thickness constraints on the wing or cross-sectional area constraints on fuselage/nacelle, which becomes a great challenge for SBO.Generally, the geometric parameter (e.g.thickness or cross-sectional area) should be no less than a certain value at required locations, ) where c i (x) denotes the remaining constraints, such as the lift and moment constraints, and N c is the number of c i (x).g i (x) denotes the geometric constraint.gp 0,i ,gp i stand for the geometric parameter of the baseline aerodynamic shape and the new shape, respectively.N t is the number of geometric constraints.

Evaluating all geometric constraints directly (EGCD)
Solving geometric constraints directly during the suboptimization of a SBO can provide true values of the geometric constraints for the optimization.Then, Eq. ( 19) becomes where g i (x) is the real geometric constraint value, which is calculated directly in the sub-optimization.The values of f (x) and c(x) are still obtained by predicting their surrogate models.

Building surrogate models for all geometric constraints (BSMGC)
This method is to build surrogate models for the objective and constraints in a one-by-one manner, Eq. ( 19) becomes where ' ∧ ' stands for the surrogate model of the corresponding item of Eq. ( 19).
The probability of feasibility of a sample point, satisfying all constraints, is the product of the PoF of each constraint, which is calculated using the surrogate models of the constraint functions.Then, the new updates can be obtained via solving the sub-optimization problem defined by the ISC.

The KS method
If the surrogate model of the aggregated constraint is built, the optimization problem becomes where g KS denotes the aggregated constraint of original constraints 19).The subscript denotes the KS function, which could aggregate the constraints into one, is defined as where ρ is the aggregation parameter.g max (x) denotes the maximum value among all constraints at the current sample x, which is defined as The range of the KS function at a sample point x is From Eq. ( 25), one can observe that, (1) the upper bound value of the KS function decrease with the increase of ρ, and when ρ gets close to infinity, the value of the KS function just becomes equivalent to g max (x); (2) If KS [g(x)] ≤ 0 is satisfied, g max (x) must be less than zero, which may lead to a conservative evaluation of the original active constraints.Therefore, the feasible region defined by the KS function can be a subset of the original feasible region.
In the following single-variable example, we use the KS function to aggregate the inequality constraints.Consider the following convex inequality constraints: The expression of the aggregation constraint is Figure 2 shows that the aggregated constraint changes with different values of ρ.Two numerical corners (corner 1 and 2) are located at where the constraints intersect in the feasible region, and the details of these two numerical corners on the right side of Figure 2. It can be observed clearly that the value of the KS function approaches the maximum constraint g max (x) as ρ increases.If ρ is large enough, as ρ = 1000, the aggregated constraint is almost coinciding with the original active constraints.

Discussion of computational time of sub-optimization by using different constraint handling methods
To evaluate the contribution of the KS constraint method to the optimization efficiency, we compare the computational time of three constraints handling methods mentioned above by a 40-design-variable wing ASO case.The objective is to minimize the drag of the wing, 100 thickness constraints and another aerodynamic constraint are used to limit the wing internal volume and lift coefficient.To obtain new samples that satisfy all constraints, performing the sub-optimization is essential.In our case, GA is used for sub-optimization, with the maximum evolutional generation set as 500 and the population size set as 300.The computational cost of once suboptimization performed by GA in each iteration of SBO can be expressed as where n ME , n PS denote the maximum evolutional generation and the population size set for GA respectively.Their product n ME • n PS denotes the number of iterations of sub-optimization.t * is the computational cost of calculating all real geometric constraint values for once by generating a new aerodynamic shape.n SM denotes the number of surrogate models, and t is the computational cost of evaluating the objective or a constraint value for once by using corresponding surrogate models.
(t * + n SM • t ) denotes the computational cost of evaluating the values of the objective and all constraints for once time.
The computational time of sub-optimization for optimizations by using different constraint handling methods at the 1000th sample are evaluated according to Eq. ( 28) and list in Table 1.All above calculations are performed with Intel Xeon CPU e5-2690 at 2.6 GHz.From Table 1 one can see that, (1) for the optimization via EGCD method, t * = 0.35(s),t = 1.4 × 10 −3 (s).Although calculating 100 real geometric constraint values for once only takes 0.35s, the computational time of the suboptimization via EGCD can be achieved according to Eq. (28), i.e. t subopt = 500 × 300 × (0.35 + 2 × 0.0014) = 52,920 (s) = 14.70 (h); (2) For the optimization via BSMGC method, it takes about 26.5 minutes to train a surrogate model with 1000 samples.However, training 102 surrogate models using BSMGC costs up to 45.73 hours, and so many surrogate models result in 5.95 hours

Examples and results
In this section, the proposed method is verified by G9 test case.Then, it is demonstrated by three typical CFDdriven aerodynamic shape optimization examples.First, the drag minimizations of the NACA0012 airfoil and ONERA M6 wing with a small number of geometric constraints are used to compare different geometric constraint handling methods.Second, the ONERA M6 wing optimization with 1000 geometric constraints is used to further study the ability of handling a large number of geometric constraints by the proposed method.Finally, drag minimization of NACA CRM wing in viscous transonic flow is performed to investigate the ability of the proposed method for handling various geometric constraints.
The optimizations are performed using the advanced TH-1A computing system at National Supercomputer Center (NSCC) in Tianjin.Each computing node in this cluster has two 2.6 GHz Intel Xeon CPUs (e5-2690 V4) and 128 GB RAM in total.Please note that, direct comparisons between SBO and other optimization methods are beyond the scope of this article, since we mainly concern about investigating a method that can efficiently handle geometric constraints within the framework of a SBO.

Problem statement
This is a benchmark case for strongly constrained global optimization.The mathematical model is This case contains 7 design variables and 4 constraints, in which g 1 and g 4 are active constraints.It is very difficult to optimize since the feasible region is only 0.5% of the design space and the order of the objective function value ranges from 10 2 to 10 7 .Although the real optimum is unknown, the optimization community has observed a best solution of f (x * ) = 680.6300573.

Aggregation parameter setting
As mentioned in the section 4, the approximation accuracy of the KS function, which influences the optimal results to a large extent, significantly depends on the aggregation parameter ρ.Therefore, it is necessary to find an appropriate aggregation parameter of the KS function for this case.
We select two sample points from the sampled datasets generated by DoE randomly.One is feasible that satisfy all constraint, while the other is infeasible that does not satisfy any constraint.Since the aggregated constraint KS (x, ρ) can provide a conservative approximation to the boundary of original constraints, we define the relative error as where g max (x) is the maximum value of geometric constraint functions.Figure 3 shows the relative errors (REs) of both the feasible and infeasible samples get close to   zero with the increase of ρ.The values of g max and relative errors of the selected feasible and infeasible designs are shown in Table 2.One can see that the setting ρ = 1000 is reasonable for this case and the aggregation constraints are accurate enough.

Results of optimizations
LHS method is used to generate 30 initial sample points.
The ISC in conjunction with EI, MSP are applied to add 2 new sample points sequentially at each iteration of a SBO.In order to obtain new samples that satisfy geometric constraints, GA is used for sub-optimization.The maximum evolutional generation is set as 200 and the population size is set as 300.Since G9 case is difficult to optimize, the acceptable maximum number of sample points is set as 1200.To eliminate the randomness in the process of optimization, the optimization using each method of handling geometric constraints is repeated 5 times.
Figure 4 shows a typical convergence history of an optimization, in which the constant parameter ρ is set as 1000.The feasible and infeasible points are denoted by the blue squares and red triangles, respectively.The convergence history of the objective function is shown in Figure 5.The result of the optimization is listed in Table 3.It is observed that the value of objective function decreases as the increase of ρ.The best objective functional of 680.66319 is obtained when ρ = 1000 due to the accurate aggregated constraint, which is very close to the  reference value of 680.63006.Besides, all constraints are satisfied strictly.

Problem statement
The design point of NACA0012 airfoil is that the angle of attack AoA = 0°and the freestream Mach number Ma = 0.85 in inviscid flow.The optimization problem is defined as where C l = 0 means the lift coefficient of the airfoil should be equal to zero in the whole optimization process.y ≥ y baseline means the local thicknesses should be no less than the ones of the baseline geometry, at all the positions from leading edge to the trailing edge.The NACA0012 airfoil is modified slightly, which leads to a zero-thickness trailing edge: We use the free-form deformation method to parameterize the airfoil (see Figure 6).There are 19 control points on the upper and lower surface of the FFD control volume respectively.To ensure the airfoil shape to be symmetric, the displacements of the active control points on the upper surface are forced to be the same as their counterparts on the lower surface.The design variables are the y-coordinate of the active control points.As a result, there are 19 variables practically.In addition, as shown in Figure 6, 15 thickness constraints are distributed on NACA0012 airfoil and can be probed uniformly along the chord.The design space is shown in Figure 7. Please note that the optimization is solved in the fixed design space to avoid the failure of optimization due to the emergence of odd geometries.

Grid convergence study
The computational grid of NACA0012 airfoil with 51,200 cells (L2 grid) is generated by using our in-house O-type grid generator '2d-Grid'.To determine the resolution accuracy of this grid, a grid convergence study is performed.All grids (see Figure 8) with deficient size are calculated by our in-house CFD solver 'PMNS3D'.Figure 9 shows that the pressure distributions computed using all grids are in good agreement with the reference data (Batina 1991), which turns out that the in-house CFD solver has sufficient accuracy.The simulation results are listed in Table 4, one can see that the L2 grid is of sufficient accuracy, the difference in the drag coefficient for the L2 grid and the finest grid (L0 grid) is only about one drag count.

Aggregation parameter setting
We select two sample points from the sampled datasets generated by DoE randomly.One is feasible, and the other is infeasible (see Figure 10) in the design space.
Figure 11 shows how the REs of both the feasible and infeasible samples get close to zero with the increase of ρ.The values of g max and REs of the selected feasible and infeasible designs are listed in Table 5.One can see that the setting ρ = 1000 is reasonable for our optimization problem and the aggregation constraints are accurate enough.

Results of optimizations
Initial surrogate models are built based on 50 sample points generated by the LHS method.The ISC in conjunction with EI, MSE, MSP, PI and LCB are applied to add 5 new sample points sequentially at each iteration.
In order to obtain new samples that satisfy geometric constraints, GA is used for the sub-optimization.The maximum evolutional generation is set as 100 and the population size is set as 150.To compare different methods of handling geometric constraints, we set the acceptable maximum number of sample points to 600.Besides,   to eliminate the randomness in the process of optimization, the optimization using each method of handling geometric constraints is repeated 5 times.Figure 12 shows a typical convergence history of an optimization, in which the geometric constraints are calculated directly.The feasible and infeasible points are denoted by the blue squares and red triangles, respectively.The average convergence histories of the objective over time are shown in Figure 13.The results of optimizations are listed in Table 6.The time costs of NACA0012 Table 5.The maximum constraint g max and relative errors of the selected feasible and infeasible designs.airfoil designs using different constraint handling methods within the framework of a SBO are shown in Table 7.
From Figure 13, Tables 6 and 7, one can see that, (1) The using of BSMGC results in the most timeconsuming optimization, since the large number of surrogate models are built and used to predict new samples during the optimization.On one hand, the cost of training surrogate models for the objective and all geometric constraints becomes higher and higher as the number of sample points increases.On the other hand, predicting the values of the objective and constraints with so many surrogate models result in the high-cost sub-optimization.In addition, the optimization takes 51.55 hours, which is too expensive for a simple airfoil aerodynamic shape optimization driven by Euler-based CFD simulations.Therefore, BSMGC is only suitable for solving optimization problems with a very small number of geometric constraints.
(2) The minimum C d is obtained by EGCD, because the real values of geometric constraints are calculated directly and applied for the sub-optimization to select new samples.However, if it is expensive to solve the geometric constraint values, the application of this method will be limited due to the unacceptable cost of the sub-optimization solved by GA.
(3) As ρ increases, the approximate accuracy of the KS aggregated constraint to the original geometric constraints improves.As a result, the average C d decreases as ρ changes from 50 to 1000.When ρ = 1000, the accuracy of the KS function is high enough, both the minimum and the average C d obtained by the KS method are very close to those obtained by EGCD.(4) The KS method reduces the number of surrogate models by aggregating all geometric constraints into one, and the aggregated constraint is used to select new samples that satisfy all geometric constraints during the sub-optimization, which makes it more efficient than the other methods.The average total time cost of the optimization via the KS method is 11.57hours, which is only 22.44% and 35.11% of that via BSMGC and EGCD respectively.
The airfoil shapes and pressure distributions of the baseline airfoil and the best designs obtained from repeated optimizations are compared and sketched in Figures 14 and 15, respectively.One can see that the optimal airfoils become thicker at the leading and trailing edges compared with the baseline airfoil, the shock    waves are significantly weakened.The comparison of Mach number contours of baseline and optimal airfoil achieved via the KS method with ρ = 1000 is shown in Figure 16, it is clearly shown that the strongly shock wave on the baseline airfoil surface is broken into an attached shock wave and a detached shock wave.Since the attached shock wave is much weaker, the shock wave drag is greatly reduced.The thicknesses of the optimal airfoils are compared with the baseline airfoil, as listed in Table 8.One can see that the thickness constraints marked in red are active constraints, as well as all constraints are satisfied strictly because the thicknesses of these optimal shapes are bigger than those of the baseline shape.

Problem statement
In this example, drag minimization of ONERA M6 wing is performed in a fixed design space.This optimization problem subjects to lift and thickness constraints, the design condition is that the angle of attack AoA = 3.06°and the freestream Mach number Ma = 0.8395 in inviscid transonic flow.We apply the FFD method to parameterize the wing with a control volume which contains five control sections.Each control section contains five control points at the upper and lower surface respectively, as shown in Figure 17.Note that the control points (shown in black color) at the trailing edge are also fixed in this case, which results in 40 design variables in total.The drag minimization problem is formulated as (33)   The optimization problem is defined to subject to 20 (4 spanwise × 5 chordwise) full thickness constraints (see Figure 18).

Grid convergence study
The structured C-H grid with 1,245,816 cells (L2 grid) is generated by our in-house 3D grid generator '3D_CH_grid' for ONERA M6 wing.A grid convergence  study is performed to determine the resolution accuracy of this grid.All grids (see Figure 19) are calculated by our in-house CFD solver 'PMNS3D'.Figure 20 shows that the pressure distributions computed using all grids are in good agreement with the reference data (Batina, 1991), which turns out that the in-house 3D CFD solver has sufficient accuracy.The simulation results are listed in Table 9, one can see that the difference in the drag coefficient for the L2 grid and the finest grid (L0 grid) is only 1.2 drag counts, which turns out that the L2 grid is of sufficient accuracy.

Aggregation parameter setting
The feasible and infeasible sample points are selected from the sampled datasets generated by DoE randomly.Table 10.The maximum constraint g max and relative errors for the selected feasible and infeasible designs.10 lists detail information about g max and REs of the selected sample points.It is obvious that the aggregated constraints with ρ = 2000 is accurate enough for this optimization case.

Results of optimizations
For this case, 100 initial sample points are generated to build the initial surrogate models for the objective and constraint functions.EI, MSE, MSP, PI and LCB are used to select 5 new samples sequentially, the responses are solved in parallel to update the surrogate models.GA is used for sub-optimization to obtain new samples that satisfy such large number of geometric constraints, the maximum evolutional generation is set as 500 and the population size is set as 300.The maximum number of sample points is set as 750 and the optimization is repeated for 3 times to eliminate the randomness in the process of optimization.
A typical optimization convergence history is shown in Figure 22, in which the geometric constraints are handled by the KS method.The average convergence history of the objective over the same time are shown in Figure 23.The pressure distributions of the baseline and optimal shapes obtained from repeated optimizations are compared and sketched in Figure 24.The results of optimizations are listed in Table 11.The time costs of ONERA M6 wing designs using different constraint handling methods within the framework of a SBO are shown in Table 12.From Figures 23 and 24, Tables 11 and 12, one can  see that, compared with the other geometric constraint handling methods, the KS method makes SBO achieve the preset maximum number of sample points first after 130 iterations, and the optimal wing shape with the minimum drag coefficient is obtained by the KS method within the same run time.It is obvious that without such an efficient geometric constraint handling method, performing a wing shape optimization that subjects to so many geometric constraints would be very intractable.
The thicknesses of the optimal wings are compared with the baseline, as listed in Table 13.We can see that most thickness constraints are active constraints.Moreover, all geometric constraints are satisfied strictly.

Further study of the KS method for handling a large number of geometric constraints
In order to explore the ability for handling large-scale geometric constraints of the KS method, the optimization  of M6 wing is defined to subject to 1000 (20 panwise × 50 chordwise) full thickness constraints, which is shown in Figure 25.Two sample points are selected from the sampled datasets generated by DoE randomly.One is feasible and the other one is infeasible.Figure 26 shows how relative errors get close to zero with the increase of ρ.Table 14 lists the values of g max and relative errors of the selected sample points.One can see that the aggregated constraints with ρ = 2000 is accurate enough for this optimization case.
For this case, the maximum number of sample points is set as 1000 to get a convergent optimization and the optimization is repeated for 3 times.The other optimization parameters for SBO are set the same as those mentioned in section 6.3.4.A typical optimization convergence history is shown in Figure 27.The average convergence history is shown in Figure 28.It can be seen that all optimizations are converged by using the KS method.The pressure distribution of the optimal wing shape is compared with the baseline, as shown in Figure 29.One can see that for the optimal wing shape, the shock is weakened significantly, as well as all constraints are satisfied.
The optimization results are listed in Table 15.We can see that, compared with the baseline, large drag reduction is achieved.Please note that, since the number of the geometric constraints is much more than 20, the minimum   objective value is slightly larger than that obtained via solving the problem with 20 geometric constraints by using the KS method.We analyze the results of the 20constraint-problem and find that although all the geometric constraints are strictly satisfied, the thicknesses are smaller than the original shape at many locations that are not constrained.It is obvious that the KS method is an efficient geometric constraint handling method and can be used to solving aircraft design problems with large number of constraints.

Problem statement
In this example, drag minimization of NASA CRM wing is performed in a fixed design space.The design    points, resulting 90 control points in total.In order to ensure that the locations of the leading edge and the trailing edge of the wing does not change, the two control points at the leading edge and trailing edge of each control section share one design variable to move the same distance along the opposite direction.We also choose ydirection displacements of the other 70 control points and the angle of attack as design variables.So there are 81 wing shape variables in total.The mathematical model of the optimization is Min. C D s.t.C L = 0.5 C MZ ≥ −0.17 t i − 0.85t i,0 ≥ 0, i = 0, . . ., 1400 area j ≥ area j,0 , j = 0, . . ., 28 v ≥ v 0 (34) Three types of geometric constraints including thickness constraints, cross-sectional area constraints and the volume constraint are enforced on the wing (see Figure 31).Meanwhile, there are 1400 thickness constraints imposed at 50 chordwise and 28 spanwise locations of the wing, they are enforced to be no less than 85% of ones of the baseline wing at each location.The cross-sectional areas of these 28 wing spanwise locations are enforced to be no less than those of the baseline wing, as well as the volume is enforced to be no less than that of the baseline wing.
As a result, there are 1429 geometric constraints in total.
Another 2 aerodynamic constraints are that the lift coefficient of the wing is required to be equal to 0.5, and the pitch moment of the wing C MZ is enforced to be no less than −0.17.The mathematical model of the optimization is

Grid convergence study
The multi-block CFD structured grid with 2,442,240 cells (L1 grid) is generated by ANAYS ICEM CFD for NASA CRM wing.A grid convergence study is performed to determine the resolution accuracy of this grid.All grids (see Figure 32) are calculated by our in-house CFD solver 'PMNS3D'.Figure 33 shows that the pressure distributions computed using all grids are in good agreement with the reference data (Kenway & Martins, 2016), which turns out that the in-house CFD solver has sufficient accuracy.The simulation results are listed in Table 16, one can see that the difference in the drag coefficient for the L1 grid and the finest grid (L0 grid) within 1 drag count, which turns out that the L1 grid has sufficient accuracy.

Aggregation parameter setting
The feasible and infeasible sample points are selected from the sampled datasets generated by DoE randomly.Figure 34 shows how REs get close to zero with the increase of ρ.Table 17 lists detail information about g max and REs of the selected sample points.It is obvious that the aggregated constraints with ρ = 5000 is accurate enough for this optimization case.

Results of optimizations
For this case, 100 initial sample points are generated to build the initial surrogate models for the objective and constraint functions.EI and MSP are used to select 2 new samples sequentially, the responses are solved in parallel to update the surrogate models.GA is used for suboptimization to obtain new samples that satisfy such large number of geometric constraints, the maximum evolutional generation is set as 500 and the population size is set as 300.The maximum number of sample points is set as 700 and the optimization is repeated for 2 times to eliminate the randomness in the process of optimization.A typical optimization convergence history is shown in Figure 35, in which the geometric constraints are handled by the KS method.The average convergence history of the objective over the same time are shown in Figure 36.The pressure distributions of the baseline and optimal shapes obtained from repeated optimizations are compared and sketched in Figure 37.The results of optimizations are listed in Table 18.From Figure 37 and Table 18, one can see that, compared with the baseline, the shock is weakened significantly for the optimal wing shape, the drag coefficient drops from 205.9 counts to 197.8 counts, as well as all constraints are satisfied.It is obvious that the KS method is an efficient geometric constraint handling method and can be used to solving aircraft design problems with a large number and multi-types of geometric constraints.

Conclusions
In this article, the constraint handling methods for solving CFD-driven aerodynamic shape optimization problems with geometric constraints were investigated within the framework of a SBO, in order to improve the capability of SBO for the engineering design optimizations.First, G9 numerical test case was performed to verify whether the KS method can handle the constraints effectively.Second, a comparative study of different methods of handling geometric constraints was presented via the aerodynamic shape optimizations of NACA0012 airfoil and M6 wing in transonic flows.Moreover, the M6 wing aerodynamic shape optimization with up to 1000 geometric constraints was conducted to explore the ability of the KS method for handling a large number of geometric constraints.Finally, drag minimization of NACA CRM wing in viscous transonic flow is performed to investigate the ability of the proposed method for handling various geometric constraints.According to the current test cases, following conclusions can be drawn as (1) The cost of SBO via building surrogate models for all geometric constraints (BSMGC) is very expensive, which induces a long design cycle.As the number of constraints increases, the cost becomes higher and higher, even prohibitive.Therefore, it is only suitable for solving optimization problems with very small number of constraints.(2) Evaluating all geometric constraints directly (EGCD) is not efficient, since for each of sub-optimizations using genetic algorithm (GA) the geometric modeling of aerodynamic configuration has to be conducted thousands of times or even more, and the total cost associated with a number of updating cycles can be extremely large.Therefore, the method EGCD will be limited due to aerodynamic shape optimization of simple configuration and with small number of updating cycles.
(3) The proposed method of this article can be used to handle constraints effectively within the framework of a SBO by aggregating all inequality constraints into one.As ρ increases, the approximate accuracy of the KS constraint to the original geometric constraints improves.As a result, the KS method can dramatically improve the optimization efficiency of CFD-driven aerodynamic shape optimizations with the number of geometric constraints in the range  from 15 to 1000 and the number of types of geometric constraints up to 3, which improves the usability of SBO in engineering design optimization and offers great potential for handling a larger number and more types of geometric constraints.(4) The parameter ρ needs to be selected according to a specific example, because the accuracy of the KS function may vary from a sample point to another and it is sensitive to the number of the geometric constraints.The method to obtain the optimal parameter ρ should be investigated to make the KS function has high accuracy for any sample.
The future work beyond the scope of this article will focus on improving the method of achieving the optimal parameter ρ to make the KS function has high accuracy for any sample, as well as the optimization of the complex aerodynamic shape such as wing-body-tail configuration with a large number and many types of geometric constraints.

Figure 2 .
Figure 2. Schematics of the aggregated constraint changing with different values of ρ.

Figure 3 .
Figure 3. Schematics of the relative errors changing with the increase of ρ.

Figure 4 .
Figure 4.The initial DoE and infill-sampling process of a typical SBO.

Figure 5 .
Figure 5. Convergence histories of the optimizations.

Figure 6 .
Figure 6.Design variables and thickness constraints on NACA0012 airfoil.

Figure 7 .
Figure 7. Illustration of the design space.

Figure 8 .
Figure 8. O-type grids of NACA0012 airfoil for the grid convergence study.

Figure 9 .
Figure9.Comparison of pressure distributions computed using all grids with reference data(Poole et al., 2015).

Figure 10 .
Figure 10.The selected feasible and infeasible designs compared with baseline airfoil.

Figure 11 .
Figure 11.Schematics of the relative errors changing with the increase of ρ.

Figure 12 .
Figure 12.The initial DoE and infill-sampling process of a typical SBO.

Figure 13 .
Figure 13.Average convergence histories of aerodynamic shape optimization of NACA0012 airfoil.

Figure 15 .
Figure 15.Comparison of surface pressure distributions of NACA0012 airfoil with optimal airfoils.

Figure 16 .
Figure 16.Comparison of Mach number contours of baseline and optimal airfoil.

Figure 19 .
Figure 19.C-H grids of ONERA M6 wing for the grid convergence study.

Figure 20 .
Figure20.Comparison of pressure distributions computed using all grids with reference data(Batina, 1991) at two typical cross-sections.

Figure 21 .
Figure 21.Schematics of the relative errors changing with the increase of ρ.
Figure21shows how REs get close to zero with the increase of ρ.Table10lists detail information about g max and REs of the selected sample points.It is obvious that the aggregated constraints with ρ = 2000 is accurate enough for this optimization case.

Figure 22 .
Figure 22.The initial DoE and infill-sampling process of a typical SBO.

Figure 23 .
Figure 23.Average convergence histories of wing aerodynamic shape optimization.

Figure 24 .
Figure 24.Comparison of pressure distributions and geometric shapes of the baseline and optimal wings.

Figure 26 .
Figure 26.Schematics of the relative errors changing with the increase of ρ.

Figure 27 .
Figure 27.The initial DoE and infill-sampling process of a typical SBO.

Figure 28 .
Figure 28.Average convergence history of wing aerodynamic shape optimization.

Figure 29 .
Figure 29.Comparison of pressure distribution of the baseline and the optimal wing.

Figure 31 .
Figure 31.Geometric constraints enforced on the NASA CRM.

Figure 32 .
Figure 32.O-H grids of NACA CRM wing for the grid convergence study.

Figure 33 .
Figure 33.Comparison of pressure distributions computed using all grids with reference data (Kenway & Martins, 2016) at typical crosssections.

Figure 34 .
Figure 34.Schematics of the relative errors changing with the increase of ρ.

Figure 35 .
Figure 35.The initial DoE and infill-sampling process of a typical SBO.

Figure 36 .
Figure 36.Average convergence histories of wing aerodynamic shape optimization.

Figure 37 .
Figure 37.Comparison of pressure distributions and geometric shapes of the baseline and optimal wing.

Table 1 .
Time cost of the wing design with different constraint handling methods at the 1000th sample.In contrast, large benefit can be obtained by the proposed method.It uses KS function to aggregate all geometric constraints into one, and then the cheap surrogate model of the aggregated constraint is used to search for new shapes, which can avoid the large computational cost associated with building surrogate models for all geometric constraints and tremendous calculation of real geometric constraint values caused by sub-optimization.It is obvious that only the proposed method can be applied to the wing design instead of two other methods.

Table 2 .
The maximum constraint g max and relative errors for the selected feasible and infeasible design.

Table 3 .
Results of G9 test case with the KS method.

Table 4 .
Mesh-convergence study for the baseline NACA0012 airfoil.

Table 6 .
Results of NACA0012 airfoil design with different constraint handling methods.

Table 7 .
Time cost of NACA0012 airfoil design within the framework of a SBO (hours).
BSMGC: building surrogate models for geometric constraints.EGCD: evaluating all geometric constraints directly.KS: the KS method.

Table 8 .
Thicknesses of optimal airfoils compared with the baseline airfoil.

Table 9 .
Mesh-convergence study for the baseline ONERA M6 wing.

Table 11 .
Results of M6 wing optimization with 20 geometric constraints.

Table 12 .
Time cost of M6 wing design within about 210 h.

Table 13 .
Thicknesses of optimal wings compared with the baseline.

Table 14 .
The maximum constraint g max and relative errors of the selected feasible and infeasible designs.

Table 15 .
Optimization results of ONREA M6 wing optimization with 1000 geometric constraints.

Table 16 .
Mesh-convergence study for the baseline NASA CRM wing.

Table 17 .
The maximum constraint g max and relative errors for the selected feasible and infeasible designs.

Table 18 .
Results of NASA CRM wing optimization with 3 types of geometric constraints.