Research of multi-point infill criteria based on multi-objective optimization front and its application on aerodynamic shape optimization

Abstract

In this article, efficient multi-point infill criteria based on multi-objective optimization front and its applications on analytic functions and airfoil aerodynamic shape optimization are researched. In the expensive optimization, an excellent infill algorithm should meet the demand of local exploitation and global exploration at the same time, which stand for the optimization accuracy and efficiency, respectively. In the proposed infill criteria, the predicted value and a statistic variable are chosen as objectives to implement optimization, and the obtained front stands for the points which strike a balance between local optimization accuracy and global optimization efficiency. The test results of four analytic functions, which has 2 to 20 variables, show that compared with expected improvement algorithm, the proposed multi-point infill algorithm can save 85% iteration times, so the computation efficiency could be increased remarkably. The airfoil optimization design is adopted to test the validity and efficiency of the proposed infill criteria. In the airfoil optimization design, the computational fluid dynamic method is employed to obtain the aerodynamic forces of airfoils, and the Kriging surrogate model is introduced to reduce the computational cost. The global optimization is implemented by adopting multi-objective evolutionary algorithm based on decomposition as optimization algorithm, and the drag coefficient and section area of the airfoil are taken as optimization objectives. The proposed infill method and expected improvement method are employed to implement the optimization design process. With roughly equal times of aerodynamic computation, both of the algorithms can obtain the reliable Pareto front; however, the number of infill iteration of proposed infill method is only 17% of that of expected improvement method.

Keywords

Infill criteria Kriging surrogate model airfoil optimization design expensive optimization multi-objective evolutionary algorithm based on decomposition

Introduction

The surrogate models introduced during expensive optimization design can efficiently deal with the conflict objectives, such as computational efficiency and accuracy, so they are widely used in expensive optimization design problems.^1–6 To construct the surrogate model, the first step is to choose a number of typical samples from design space, and then their responses to objective are computed by proper analysis tool with certain precision. After that, some mathematical approaches such as interpolation and fitting are employed to construct the surrogate model based on these responses of typical samples. Theoretically, with sufficient well-distributed samples, the surrogate model would be accurate enough and can be utilized to replace real analysis model to predict the response of arbitrary sample. However, for practical applications, without the knowledge of the objective function’s properties, it is hard to describe its global characteristics with few sample points. If large number of samples are used to construct the surrogate model, that would result in great computational cost. To overcome this conflict when utilizing the surrogate model to perform the global optimization design, the general approach is as follows. First, the initial surrogate model is established by a few samples, and the initial number of samples is suggested to be $10 n$ ⁷ or $(n + 1) (n + 2) / 2$ ,⁸ where n is the dimension number. Then, the surrogate model is improved by new samples via certain multi-point infill strategy. This iterative procedure carries on until the surrogate model reaches the given precision.

The EGO algorithm proposed by Jones et al.⁹ shows that infill criteria should balance the global exploration and local exploitation when determining the positions of the samples. The global exploration focuses on improving the precision of positions with relatively large uncertainty, while the local exploitation focuses on the positions near the optimal design. Mockus et al.¹⁰ proposed the famous “expected improvement” (EI) infill criteria, which try to improve the precision of the surrogate model by searching for the position with the largest EI value on the design space. Since it well balances the global exploration and local exploitation, this approach is broadly utilized.^11–13 However, the EI method converges slowly and is inefficient, which limited its application on expensive optimization design problems.¹⁴ To improve the computation efficiency, researchers propose some advanced infill criteria based on the EI criteria, such as augmented EI criteria¹⁵ and bootstrapped EI criteria.¹⁶ These new EI methods can improve the computation efficiency in certain extent; however, they make some problems such as computation stability and robustness.

Along with the development of computation hardware and parallel computation methods, the multi-point infill criteria has been well developed. Z Feng et al.¹⁷ proposed a novel multi-point infill sequential optimization method. In this method, the samples are obtained from the approximated Pareto front of two objectives which is the two terms of the EI function. By following the approach given in Feng et al.,¹⁷ efficient multi-point infill criteria based on multi-objective optimization front are proposed in this article, and its application on aerodynamic shape optimization is also researched to verify the algorithm.

The remainder of this article is arranged as follows. In the next section, the details of the Kriging surrogate model and multi-point infill criteria are illustrated. Next, four typical analytic functions are employed to test the proposed infill criteria. Then, the multi-point infill criteria proposed in this article are applied to two-dimensional airfoil aerodynamic shape optimization. This article ends with the main conclusions and recommendations for future work.

Multi-point infill criteria based on multi-objective optimization front

Kriging surrogate model

The surrogate models are mathematically constructed approximation models, which have lower computation cost but with high accuracy when compared with numerical analyzing or physical experimentation results. Surrogate models are mainly divided into response surface methods (RSM), radial basis functions (RBF), artificial neural networks (ANN), and Kriging models. Compared with the RSM and RBF, the ANN is more suitable to deal with strong nonlinearity and multi-dimensional problems. However, prediction accuracy and robustness of the ANN are inferior to the Kriging method.¹⁸

The Kriging method is utilized to replace the aerodynamic analysis tool to construct the surrogate mode, which is briefly given as follows. In the Kriging model, the relationship of the attributive and independent variables is given by¹⁹

y = f {(x)}^{T} β + z (x)

(1)

where $f (x)$ is the linear combined function of functions $f_{i} (x)$ $(i = 1, 2, \dots, p)$ , called as regression model. $β$ is the regression parameter. Different combinations of $f_{i} (x)$ yield different regression models. According to the highest order of $f_{i} (x)$ , the regression model can be divided into constant model, first-order model, and second-order model, refer to Ma²⁰ for the details. $z (x)$ is the Gauss random process error function, with the mean value of 0, standard derivation $σ^{2}$ , and covariance $σ^{2} R (θ, x^{i}, x^{j})$ . For two sample points $x^{i}$ and $x^{j}$ that related by $θ$ , expression of $R (θ, x^{i}, x^{j})$ yields

R (θ, x^{i}, x^{j}) = Π_{l = 1}^{n} R_{l} (θ_{l}, x_{l}^{i} - x_{l}^{j}) = Π_{l = 1}^{n} R_{l} (θ_{l}, d_{l})

(2)

where $R_{l} (θ_{l}, d_{l})$ is the kernel function of the correlation function. $d_{l}$ denotes the distance between any two sample points, and $θ_{l}$ denotes the relative variable of the sample point along the lth direction. The kernel functions mainly the Gaussian function, exponential function, linear function, spherical function, and spline function.¹⁹

Actually, for the specified regression model and the relative model type, the correlation matrix R and the predicted values $\hat{β}$ and ${\hat{σ}}^{2}$ are determined by $θ$ .²¹ By utilizing maximum-likelihood function, solution of $θ$ can be transformed into an optimization problem, which is solved by optimization algorithms.⁹

Multi-point infill criteria algorithm

Note that the predicted optimal value $\hat{f} (x)$ based on the current Kriging model conflicts with the standard derivation $s (x)$ . Meanwhile, $\hat{f} (x)$ and $s (x)$ also reflect the local exploitation and global exploration abilities, respectively. Thus, the minimizing f and maximizing s are chosen as optimization objectives. The optimization problem is given by

{\begin{matrix} min F = {(\hat{f} (x), - s (x))}^{T} \\ x \in Ω \end{matrix}

(3)

The multi-objective evolutionary algorithms based on decomposition (MOEA/D) are utilized to perform the multi-objective optimization. According to the obtained Pareto front, the correlated samples are deleted.²² Details of the multi-point infill sampling procedure are given as follows (Figure 1):

Use MOEA/D to optimize the problem in equation (3). Samples on the Pareto front are taken as sampling candidates.

Two end points of the Pareto front are directly taken as samples, since they indicating the current optimal objective function and maximum uncertainty, respectively.

Delete the candidates that correlated with the current sample set.

Delete the correlated samples in the candidate sample set.

K denotes the number of samples added during each iteration. If the candidates’ number is smaller than K, then all candidates are taken as samples. Otherwise, increases the correlative estimation critical distance $ε_{d}$ and return to steps 3 and 4 until the candidates’ number is equal to or smaller than K.

Figure 1.

Algorithm procedure of multi-point infill criteria based on multi-objective optimization front.

In the previous procedure, distance between any two samples is defined as

d = \sqrt{{\sum_{i = 1}^{n} (x_{i}^{(1)} - x_{i}^{(2)})}^{2}}

(4)

where $x_{i}^{(1)}$ and $x_{i}^{(2)}$ are relatively the ith terms of random two samples $x^{(1)}$ and $x^{(2)}$ , respectively.

Then, any two samples are termed as correlated when d in equation (4) is smaller than the specified critical distance $ε_{d}$ .

Algorithm test

Test functions

Four typical analytic functions are chosen as the test functions, of which three functions have two design variables, and one function has 20 variables to test the ability of the proposed method dealing with complicated-dimensions problem:

Six-hump camel-back (SC) function

\begin{array}{l} f (x) = 4 x_{1}^{2} - 2.1 x_{1}^{4} + \frac{1}{3} x_{1}^{6} + x_{1} x_{2} - 4 x_{2}^{2} + 4 x_{2}^{4} + 2, \\ - 5 \leq x_{i} \leq 5 \end{array}

(5)

Branin (BR) function

\begin{matrix} f (x) = {(x_{2} - \frac{5.1}{4 π^{2}} x_{1}^{2} + \frac{5}{π} x_{1} - 6)}^{2} \\ + 10 (1 - \frac{1}{8 π}) \cos (x_{1}) + 10 \\ - 5 \leq x_{1} \leq 10, 0 \leq x_{2} \leq 15 \end{matrix}

(6)

Goldstein–Price (GP) function

\begin{array}{l} f (x) = (1 + {(x_{1} + x_{2} + 1)}^{2} \\ (19 - 14 x_{1} + 3 {x_{1}}^{2} - 14 x_{2} + 6 x_{1} x_{2} + 3 {x_{2}}^{2})) \\ \times (30 + {(2 x_{1} - 3 x_{2})}^{2} (18 - 32 x_{1} + 12 {x_{1}}^{2} \\ + 48 x_{2} - 36 x_{1} x_{2} + 27 {x_{2}}^{2})) - 2 \leq x_{i} \leq 2 \end{array}

(7)

Styblinski–Tang (ST) function

f (x) = \frac{1}{2} \sum_{i = 1}^{20} (x_{i}^{4} - 16 x_{i}^{2} + 5 x_{i}) (- 5 \leq x \leq 5)

(8)

The global optimal solutions and corresponding function values of the test functions are given in Table 1.

Table 1.

The optimal solutions and corresponding function values of the test functions.

Test function	SC function	BR function	GP function	ST function
Global optimal solutions	$x^{*} = {\begin{matrix} (0.0898, - 0.71261) \\ (- 0.0898, 0.71261) \end{matrix}$	$x^{*} = {\begin{matrix} (- π, 12.275) \\ (π, 2.275) \\ (3 π, 2.475) \end{matrix}$	$x^{*} = (0, - 1)$	$x^{*} = (- 2.903534, \dots, - 2.903534)$
Optimal function values	$f (x^{*}) = 0.9684$	$f (x^{*}) = 0.397887$	$f (x^{*}) = 3$	$f (x^{*}) = - 783.32$

SC: six-hump camel-back; BR: Branin; GP: Goldstein–Price; ST: Styblinski–Tang.

Test method and criterion

The computation condition during the iterations is illustrated as follows. The number of initial samples for SC and BR function is 25. Because of the nonlinear character of the GP and ST function, the number of initial samples for these two functions is relatively large, namely, 100 and 700, respectively. For each function, the initial samples obey the uniform distribution on the design space. For SC, BR, and GP functions, the number of samples added during each iteration K is 5, and for ST function, K is 10. The convergence condition is given by

ε = | \frac{f_{min} - f (x^{*})}{f (x^{*})} | \leq ε_{c}

(9)

where $ε$ is the relative error between predicted and real optimal function values according to the definition expression. $ε_{c}$ is the critical value of relative error. For SC and BR function, $ε_{c}$ is 3%, and for GP and ST function, $ε_{c}$ is 25%.

Three factors are utilized to evaluate the infill criteria. First, the relative error between predicted and real optimal function values $ε$ , defined by equation (9), which denotes the iteration precision of predicted optimal value. Second, the number of objective function computation $N_{compute}$ , which do not include the computation times of initial samples. For expensive optimization, the computation cost of objective function is huge, so $N_{compute}$ should be as small as possible. The factor $N_{compute}$ denotes the entire computation quantity. Third, the number of iteration $N_{iterate}$ . During every iteration, the computation of several samples is implemented by parallel computation, so small $N_{iterate}$ means little computation time and cost.

Test results

The classical EI infill criterion is employed to improve the accuracy of the Kriging model to compare to the proposed multi-point infill criterion. Except for the different iteration algorithms, the two infill sampling processes have the same computation condition, including the number of initial samples, the initial distribution of samples, convergence condition, and so on. The iteration results of three test functions using EI method and proposed multi-point infill criterion are depicted in Table 2. For SC and BR function, the results are similar. The relative errors $ε$ based on EI algorithm are smaller than those based on multi-point infill algorithm, and the number of objective function computation $N_{compute}$ based on EI algorithm is less. However, the number of iteration $N_{iterate}$ based on EI algorithm is much more than those based on multi-point infill algorithm. In fact, $N_{iterate}$ based on EI algorithm is 2 or 3.6 times as much as that based on multi-point infill algorithm, which means that the computation time based on EI algorithm is two to four times as much as that based on multi-point infill algorithm. For GP function, $ε$ based on EI algorithm is slightly smaller than that based on multi-point infill algorithm, and $N_{compute}$ based on EI algorithm is slightly more. Similarly, $N_{iterate}$ based on EI algorithm is 4.6 times as much as that based on multi-point infill algorithm. For ST function, $ε$ and $N_{compute}$ based on EI algorithm is slightly smaller than that based on multi-point infill algorithm, but the $N_{iterate}$ based on EI algorithm is 8.5 times as much as that based on multi-point infill algorithm.

Table 2.

The iteration results of test functions using two infill algorithms.

Test function	Results based on EI algorithm			Results based on multi-point infill algorithm
Test function	$ε$	$N_{compute}$	$N_{iterate}$	$ε$	$N_{compute}$	$N_{iterate}$
SC function	0.58%	33	33	2.5%	39	9
BR function	0.63%	10	10	2.6%	25	6
GP function	21%	55	55	25%	51	12
ST function	23.3%	502	502	24%	538	59

EI: expected improvement; SC: six-hump camel-back; BR: Branin; GP: Goldstein–Price; ST: Styblinski–Tang.

In summary, compared with EI algorithm, the proposed multi-point infill algorithm might need slightly more computation quantity to achieve the same level of predicted precision, but it can save much more iteration times (actually the number of iteration decreases about 85% on average), so the computation efficiency could be increased remarkably.

Application research on airfoil aerodynamic shape optimization design

According to the conclusion of the section above, the multi-point infill criteria can increase the computation efficiency remarkably, which is vital to expensive optimization. However, is the conclusion still effective for real complex physical problem? So in this section, an application research of the multi-point infill criteria on two-dimensional airfoil aerodynamic shape optimization problem is implemented to test its performance.

Aerodynamic shape optimization design problem

Geometry representation of the airfoil

Airfoil optimization design requires an efficient geometry representation method. The class function/shape function transformation (CST) method²³ is adopted to generate the airfoil curves. The four-order Bernstein polynomial is utilized to approximate the upper and lower airfoil outline curves. The exponential of the class function and the height of the curve’s aft end are both fixed. To approximate the real airfoil, the first four Bernstein polynomial coefficients of the upper and lower curves are set as opposite and their fifth polynomial coefficients are set as equal. The five Bernstein coefficients and their ranges are given in Table 3. Two typical airfoil curves are depicted in Figure 2. As it shows, it is very convenient for the CST method to generate abundant airfoil curves.

Table 3.

Bernstein coefficients and their ranges of the airfoil.

Bernstein coefficients	A₀	A ₁	A ₂	A ₃	A ₄
Lower bound	0.1	0.1	0.1	0.1	0.1
Upper bound	0.6	0.5	0.5	0.6	0.6

Figure 2.

Typical airfoil shape based on CST method.

Definition of optimization design problem

Computation conditions of the optimization design problem are set as Mach number $Ma = 0.707$ , steady pressure $p = 125.2 kPa$ , steady temperature $T = 300 K$ , and angle of attack $α = 1^{\circ}$ . In transonic conditions, the drag coefficient $C_{d}$ and section area $A_{s}$ are critical indexes for the airfoil aerodynamic and structure design. Meanwhile, these two indexes conflict to each other. Therefore, the MOEA/D is utilized to design the airfoil, of which the minimizing $C_{d}$ and maximizing $A_{s}$ are taken as objectives. The optimization problem can be expressed as

\begin{matrix} Ma = 0.707, α = 1^{\circ}, p = 125.2 kPa, T = 300 K \\ min {\begin{matrix} C_{d} \\ - A_{s} \end{matrix} \end{matrix}

(10)

Aerodynamic model and its validation

Aerodynamic model

The aerodynamic problem of supercritical airfoil can be considered as a two-dimensional viscous compressible turbulent flow. The commercial software FLUENT is used to evaluate the aerodynamic performance of each sample design.

Aerodynamic model will be employed many times to analyze the aerodynamic force during the optimization design process, so the aerodynamic model should be as simple as possible to make the price of every single aerodynamic analysis small. But the computation price is just one thing, the accuracy of aerodynamic analysis is also a crucial factor. Unfortunately, the computation price and accuracy of aerodynamic analysis always conflict. The main factors in aerodynamic model influencing computation price and accuracy include computing grids, difference scheme, and turbulence model. Before formal computation, some research about influence of these factors is implemented. The study subject is the typical supercritical airfoil NASA SC(2)-0518.

Three elements of computing grids, namely, the number of circumference nodes (n_c), the height of the first grid close to the wall boundary (h₁), and the increase rate of grids along normal direction (k) are chosen to research the convergence of computing grids, and three n_c values, three h₁ values, and seven k values are considered. Twelve grids with different elements are used, and results of Mesh0 whose grid number is the largest are considered as the reference value. Table 4 shows the research results. According to the table, the drag coefficient C_d is not sensitive to n_c, so for the sake of low computation price, n_c is set to be 215. The height of the first grid h₁ and the increase rate k can influence C_d obviously, and for the sake of computation accuracy, h₁ is set to be 1e−5, and k is set to be 1.3.

Table 4.

The results of convergence research of computing grids.

Grid number	n_c	h ₁	k	Total number of grids (thousand)	C_d	Remarks
Mesh0	1040	2e−6	1.15	108	0.02099	Influence element: n_c
Mesh1	420	2e−6	1.15	43.5	0.02129
Mesh2	215	2e−6	1.15	21.8	0.02189
Mesh3	215	1e−5	1.15	19.3	0.02188	Influence element: h₁
Mesh4	215	5e−5	1.15	16.9	0.02387
Mesh5	215	2.5e−4	1.15	14.4	0.02504
Mesh6	215	1e−5	1.2	15.3	0.02194	Influence element: k
Mesh7	215	1e−5	1.25	12.7	0.02198
Mesh8	215	1e−5	1.3	10.1	0.02203
Mesh9	215	1e−5	1.35	9.5	0.02465
Mesh10	215	1e−5	1.4	8.9	0.02789
Mesh11	215	1e−5	1.45	8.27	0.0331

The first-order and second-order upwind difference scheme, and Spalart–Allmaras model and Euler model are adopted to research the influence of difference scheme and turbulence model on aerodynamic character, respectively. The wall pressure distribution of airfoil with the condition of different difference scheme and turbulence model is given in Figures 3 and 4. According to the figures, both difference scheme and turbulence model can influence the wall pressure, especially the upper surface pressure, significantly. Since the computation accuracy is the first place to consider, the second-order upwind difference scheme and Spalart–Allmaras turbulence model are chosen for aerodynamic analysis.

Figure 3.

The influence of difference scheme on wall pressure.

Figure 4.

The influence of turbulence model on wall pressure.

In summary, the second-order upwind difference scheme is adopted for the convection term, and the turbulence model used in this article is Spalart–Allmaras. There are two types of boundary condition, pressure far field and wall. The former corresponds to the outer boundary of computational domain, and the latter corresponds to the surfaces of airfoil. To obtain high-fidelity aerodynamic results, structural C-type grids which cluster round surfaces of the airfoil, are constructed. The numbers of circumferential and radial nodes are 215 and 47, respectively, and the increase rate of grids along normal direction is 1.3. The total number of grid cells of each computational case is about 10,000. The typical computational grids are shown in Figure 5.

Figure 5.

Computational grids.

Validation of aerodynamic model

The flow problem of typical supercritical airfoil NASA SC(2)-0518²⁴ is employed to verify the accuracy of the results obtained by the computational fluid dynamic (CFD) simulation. Comparison of the results is depicted in Figure 6, where the hollow squares denote the ones obtained from the proposed method and hollow deltas denote the ones obtained from Zhou and Zhou.²⁴ In the figure, five angles of attack are computed at the given flow condition. As it shows, at the condition of small angle of attack, the CFD results obtained are accurate.

Figure 6.

Validation of drag characteristic of NASA SC(2)-0518 airfoil.

Optimization algorithm and design process

Optimization algorithm

The MOEA/D is utilized to perform the global optimization. Zhang and colleagues^25,26 pointed out that MOEA/D was simple to code and with a high computational efficiency, besides the obtained optimal Pareto front points distributed uniformity. MOEA/D utilizes traditional decomposing strategy, such as Tchebycheff method, to divide the multi-objective optimization problem into multiple single-objective optimization problems, which are simultaneously optimized by the evolutionary algorithms, please refer to Zhang and Li²⁵ for the details.

Optimization design process

The Kriging surrogate model and MOEA/D are utilized to optimize the two-dimensional transonic airfoil, and the optimization design process is illustrated in Figure 7. First, the Latin hypercube method is used to select 61 initial samples in the design space. CFD and computer-aided design (CAD) are then utilized to compute their objective characteristic (drag coefficient and section area). Second, the Kriging surrogate model is constructed based on the current samples. If the constructed Kriging model satisfies the required design precision, then it can be used in the MOEA/D to design the two-dimensional transonic airfoil. Otherwise, the multi-point infill criteria based on multi-objective optimization front depicted in Figure 1 are utilized to refine the Kriging model iteratively until the precision of surrogate model meets the requirement.

Figure 7.

Optimization design process.

During the optimization design procedure, the mean error is used to determine the precision of the surrogate model, which is defined as

ME = \frac{1}{N_{test}} \sum_{i = 1}^{N_{test}} | \frac{y_{i} (x) - {\hat{y}}_{i} (x)}{y_{i} (x)} |

(11)

where $y_{i} (x)$ and ${\hat{y}}_{i} (x)$ denote the real and predicted values of the samples, respectively. $N_{test}$ denotes number of the samples. Therefore, if ME is less than 9%, the corresponding surrogate model can satisfy the required precision.

Optimization design results

Improvement of the accuracy of surrogate model based on multi-point infill criteria

The modification process of the accuracy of surrogate model is a key technique for aerodynamic shape optimization design problem. Numerical results show that the multi-point infill criteria proposed in this article can improve the accuracy of surrogate model efficiently. During the iteration, the number of samples added during each iteration is set to be 6, and after 12 infill iterations, the constructed surrogate model can satisfy the required precision.

The accuracy evolution process of drag coefficient surrogate model is illustrated in Figure 8. If the model has enough accuracy, the error between predicted and true drag coefficients should be minor, so the data points on the figure should be gathered close to the 45° oblique line. Otherwise, the data points should be far away from the oblique line. It can be seen from the figure, when the initial amount of samples is 61, the surrogate model can approximately reflect the character of drag coefficient; however, a few samples with major error exist. After three times of iterations, and the amount of samples comes to 78, the distribution of data points is still dispersed, but the trend of the data points gathers to the oblique line in some certain extent. After another six times of iterations, and the amount of samples comes to 107, the data points approach to the oblique line closer. When the amount of samples increases to 121, almost all of the data points have already concentrated around the 45° oblique line, so the surrogate model should be accurate enough to replace CFD model to analyze the aerodynamic character of airfoils.

Figure 8.

Accuracy evolution process of surrogate model.

The above adding points process has revealed that the surrogate model becomes accurate and reliable piece by piece as the increasing amount and reasonable distribution of samples on the design space. This trend can be obviously illustrated from Figure 9. When the initial amount of samples is 61, the mean error is 22.1%, which is relatively large. As the amount increases, the mean error decreases gradually, and when the amount increases to 121, the error reduces to 8.95%.

Figure 9.

Evolution process of mean error of surrogate model.

Optimization design results of airfoil

The optimized two-dimensional transonic airfoils are illustrated in Figure 10, where the solid dots denote the samples and the solid line denotes the Pareto front based on proposed method. In the optimization process, the population size used in MOEA/D is 1000. As illustrated in Figure 10, the three optimized typical airfoils on the Pareto front are as follows: (1) CaseA, which has the minimum drag coefficient and the minimum section area, identified by solid circle; (2) CaseB, which has the compatible section area and drag coefficient, identified by solid pentagram, and (3) CaseC, which has the maximum section area and the maximum drag coefficient, identified by solid square. The samples with extreme objective characteristic are identified by lines with arrows, where Case1, Case2, and Case3, respectively, denote the airfoils with the minimum section area, the minimum drag coefficient, and the maximum drag coefficient and section area simultaneously. The distribution of samples in the objective space indicates that the drag coefficient and section area of the two-dimensional transonic airfoils are strongly conflict with each other, which makes necessity of the multi-objective optimization. Many design points with consideration to both of two objective characteristics are obtained, and the distribution of Pareto front in the objective space indicates that the two objective characteristics of these design points can be improved at the same time. The designers can chose the proper points from the Pareto front as the optimal design results depending on the specific circumstances.

Figure 10.

Results of airfoil aerodynamic shape optimization problem based on proposed method and EI method.

The EI method is used as infill criterion to replace the multi-point infill criteria to implement the optimization design process (as shown in Figure 7), and the corresponding Pareto front is shown in Figure 10, which is identified by red dotted line. According to the figure, the two Pareto fronts based on proposed method and EI method are similar. The discrepancy distributes nearby the CaseB, which can be obviously seen in the magnifying figure (as shown in red rectangle). As to the convergence process, taking drag coefficient surrogate model for consideration as well, the number of objective computation is 121 and 132 for proposed method and EI method, respectively, while the number of infill iteration is 12 and 71, respectively.

To validate the computation accuracy of optimization design based on Kriging model, three typical shapes on the Pareto front are chosen as the test case, namely, CaseA, CaseB, and CaseC. The computation results of two objectives with Kriging surrogate model (KG) and physics model (CFD for drag coefficient computation, and CAD for section area computation) are illustrated in Table 5. It is obvious that under the same number of samples and distribution, the errors of section area between Kriging model and CAD model are much smaller than that of drag coefficient between two types of models. The reason is that the relation between section area and design variables is much more closer, and in fact, section area can be given by exact expression of the design variables, while the nonlinear characteristics between drag coefficient and design variables are more intense. As to drag coefficient, the predicted error based on Kriging and CFD model for CaseC, CaseB, and CaseA, are namely 0.24%, −4.72%, and −7.23%. In the table, the plus sign means that the predicted value is larger than the CFD analytical one, and the minus sign means that the predicted value is smaller than the CFD analytical one. The predicted error of CaseA is the most. However, the errors of the three airfoil are both less than 8%. This level of error can satisfy the requirement of engineering design, and it shows that the optimization design results based on Kriging model and multi-point infill criteria are reliable. In addition, the exact objective values obtained by CFD and CAD tools of the three typical airfoils are plotted in Figure 10, which is identified by red hollow circle, pentagram, and square.

Table 5.

The comparison of objective predicted values based on Kriging surrogate model and physics model.

Typical shapes on the Pareto front	Drag coefficient			Section area (m²)
Typical shapes on the Pareto front	KG	CFD	Error	KG	CAD	Error
CaseA	0.01067	0.0115	−7.23%	0.102	0.101	+0.66%
CaseB	0.0318	0.0333	−4.72%	0.191	0.189	+1.06%
CaseC	0.278	0.277	+0.24%	0.264	0.263	+0.21%

KG: Kriging surrogate model; CFD: computational fluid dynamic.

Conclusion

In the present study, efficient multi-point infill criteria based on multi-objective optimization front were proposed, and its applications on analytic functions and airfoil aerodynamic shape optimization were also researched. In the proposed infill criteria, to search for several appropriate added points in single iteration, the predicted optimal value based on the current Kriging model and the standard derivation, which reflected the local exploitation and global exploration abilities, respectively, were chosen as the optimization objectives. The multi-objective optimization front obtained by MOEA/D indicated the points which striked a balance between local optimization accuracy and global optimization efficiency. So with the proposed infill criteria, the optimization process based on surrogate model could be accelerated and the precision of the surrogate model could be improved.

Four analytic functions and the transonic airfoil aerodynamic shape optimization problem were chosen to test the validity and efficiency of the proposed infill criteria. The results of three simple analytic functions with two variables showed that compared with EI algorithm, the proposed multi-point infill algorithm could save much more iteration times, so the computation efficiency could be increased remarkably. The test of analytic function with 20 variables obtained the same conclusion, and it proved the strong compatibility of the proposed algorithm dealing with the complex problem with numerous variables. According to the results of transonic airfoil aerodynamic shape optimization, the surrogate model became accurate by the proposed infill criteria, and the reliable Pareto front was obtained based on surrogate model and MOEA/D. Compared with EI algorithm, the number of infill iteration of proposed infill method is only 17% of that of EI method.

There were two problems need to be researched in the future. One, the proposed infill criteria should be applied to many more complicated nonlinear problems to test its validity and efficiency. Two, many algorithm details such as the critical distance for judging the correlation of arbitrary two samples should be improved to enhance its flexibility and robustness.

Footnotes

Academic Editor: Mohana Muthuvalu

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The first author (Y.M.) acknowledges the funding support from the Rocket Force University of Engineering scientific research project (project number: 2015QNJJ034).

References

Yang

Feng

. Hypersonic lifting body aerodynamic shape optimization based on the multiobjective evolutionary algorithm based on decomposition. J Aerospace Eng 2015; 229: 1246–1266.

Chowdhury

Mehmani

Messac

. Adaptive model refinement in surrogate-based multiobjective optimization. In: Proceedings of the 57th AIAA/ASCE/AHS/ASC structures, structural dynamics and materials conference, San Diego, CA, 4–8 January 2016. AIAA 2016–0417.

Zhang

Han

Z-H

Shi

L-X

. Multi-round surrogate-based optimization for benchmark aerodynamic design problems. In: Proceedings of the 54th AIAA aerospace sciences meeting, San Diego, CA, 4–8 January 2016. AIAA 2016–1545.

Tesfahunegn

Koziel

Leifsson

. Surrogate-based airfoil design with multi-level optimization and adjoint sensitivity. In: Proceedings of the 53rd AIAA aerospace sciences meeting, Kissimmee, FL, 5–9 January 2015. AIAA 2015–0759.

Shi

L-X

Han

Z-H

Shahbaz

. Surrogate-based robust airfoil design under aleatory operating-conditions and geometric uncertainties. In: Proceedings of the 54th AIAA aerospace sciences meeting, San Diego, CA, 4–8January 2016. AIAA 2016–0810.

Bobrowski

Ferrery

Valero

. CAD-based aerodynamic shape optimization using geometry surrogate model and adjoint methods. In: Proceedings of the 17th AIAA/ISSMO multidisciplinary analysis and optimization conference, Washington, DC, 13–17 June 2016. AIAA 2016–3831.

Loeppky

Sacks

Welch

WJ.

Choosing the sample size of a computer experiment. Technometrics 2009; 51: 366–376.

Simpson

Peplinski

Koch

. Metamodels for computer-based engineering design: survey and recommendations. Eng Comput 2001; 7: 129–150.

Jones

Schonlau

Welch

WJ.

Efficient global optimization of expensive black-box functions. J Global Optim 1998; 13: 455–492.

10.

Mockus

Tiesis

Zilinskas

The application of Bayesian methods for seeking the extremum. Towards Global Optim 1978; 2: 117–129.

11.

Hashimoto

Kenae

. Aerodynamic optimization of high-wing configuration for near future aircraft. In: Proceedings of the 10th AIAA multidisciplinary design optimization conference, National Harbor, MD, 13–17 January 2014. AIAA 2014–0291.

12.

Song

Han

Application of improved Kriging-model-based optimization method in airfoil aerodynamic design. J Northwest Polytech Univ 2010; 28: 503–510 (in Chinese).

13.

Dan

Aerodynamic optimization method based on Kriging adaptive surrogate model. Aeronaut Comput Tech 2013; 43: 13–17 (in Chinese).

14.

Jones

DR.

A taxonomy of global optimization methods based on response surface. J Global Optim 2001; 21: 345–383.

15.

Huang

Allen

Notz

. Global optimization of stochastic black-box systems via sequential Kriging meta-models. J Global Optim 2006; 34: 441–466.

16.

Kleijnen

JPC

Beers

Nieuwenhuyse

IV.

Expected improvement in efficient global optimization through bootstrapped Kriging. J Global Optim 2012; 54: 59–73.

17.

Feng

Zhang

. A multiobjective optimization based framework to balance the global exploration and local exploitation in expensive optimization. J Global Optim 2015; 61: 677–694.

18.

Ricardo

Andre

RDC

Curran

Comparison of surrogate models in a multidisciplinary optimization framework for wing design. AIAA J 2010; 48: 995–1006.

19.

Xie

Research on robust optimization of sheet metal forming based on Kriging and grey relational analysis. Shanghai, China: Shanghai Jiao Tong University, 2007 (in Chinese).

20.

Research on the aerodynamic shape optimization of vehicle based on surrogate model and MOEA/D. Changsha, China: National University of Defense Technology, 2015 (in Chinese).

21.

Chen

Guan

The optimization of parameters of Kriging correlation model based on particle swarm optimization. Microelectron Comput 2009; 26: 178–181 (in Chinese).

22.

Feng

Research on the MOEA/D and its application in the flight vehicle dynamics. Changsha, China: National University of Defense Technology, 2014 (in Chinese).

23.

Kulfan

Bussoletti

. “Fundamental” parametric geometry representations for aircraft component shapes. In: Proceedings of the 11th AIAA/ISSMO multidisciplinary analysis and optimization conference, Portsmouth, 6–8 September 2006. AIAA 2006–6948.

24.

Zhou

CFD-based studies on airfoil selection of oblique wing. Chin Quart Mech 2013; 34: 361–366 (in Chinese).

25.

Zhang

MOEA/D: a multiobjective evolutionary algorithm based on decomposition. IEEE T Evolut Comput 2007; 11: 712–731.

26.

Zhang

Feng

Liu

. Fuel-time multiobjective optimal control of flexible structures based on MOEA/D. J Natl Univ Def Technol 2009; 31: 73–76 (in Chinese).