基于多输出高斯过程的超临界翼型优化

吴宽展，刘学军，＊，吕宏强

（1.南京航空航天大学计算机科学与技术学院，江苏南京 210016；2.南京航空航天大学航空宇航学院，江苏南京 210016）

基于多输出高斯过程的超临界翼型优化

吴宽展1，刘学军1，＊，吕宏强2

（1.南京航空航天大学计算机科学与技术学院，江苏南京 210016；2.南京航空航天大学航空宇航学院，江苏南京 210016）

提出了一种基于多输出高斯过程（MOGP）回归模型和模拟退火算法相结合的翼型优化方法。采用拉丁超立方抽样方法在设计空间内构造一系列样本点，其中优化设计采用CST参数化方法对翼型的几何外形进行参数化表示。通过CFD流体计算得到其响应值来建立初始的MOGP代理模型。以阻力最小化为设计目标，考虑面积、升力等约束条件。通过单点优化和多点优化试验表明，发展的翼型优化设计方法达到了优化设计的目的，同时也说明基于MOGP模型的优化设计方法在气动优化设计中的应用是可行的。

翼型优化设计；输出高斯过程；CST参数化方法；模拟退火算法；超临界翼型

0 引 言

随着计算机计算能力的提升与计算技术的发展，CFD（计算流体动力学）软件已经成为流场分析、流场计算、流场预测的重要工具。在翼型气动优化过程中，由于需要在翼型设计空间大量调用高精度流场计算模拟，基于CFD的优化计算会受到时间、财力的限制。一个解决方案是采用代理模型［1］技术（又称响应面技术）构建近似模型来代替复杂的CFD计算。近似模型描述了目标函数（输出）与设计变量（输入）之间对应关系。通过近似模型计算目标函数值节省了大量计算时间并且可以搜索更广阔的设计空间。一些试验表明基于代理模型的优化技术不仅高效而且可靠，已逐渐在工程上得到应用［2－5］。

翼型设计领域中常用的代理模型主要有多项式模型，人工神经网络模型，Kriging模型（单输出高斯回归模型）［3－6］等。高斯过程回归模型本质上是一种基于贝叶斯框架的回归方法，有着严格的统计学理论基础。高斯模型能很好地处理小样本、非线性等复杂的问题。利用统计学方法来寻找输入变量与输出变量之间的条件分布，从而建立两者的回归方程。高斯回归模型分为单输出高斯过程（Single Output Gaussian Process，即Kriging）回归模型［7］和多输出高斯过程（Multiple Output Gaussian Process，MOGP）回归模型［8－9］。单输出模型对每个输出端口分别建模，忽略了端口之间的潜在相关性。MOGP回归模型则为各个输出端口之间建立了协方差矩阵来模拟它们之间可能存在的依赖关系，每个端口都可以利用其它输出端口来提升自己的预测性能。

在先前的工作中，我们验证了MOGP模型相比Kriging模型在多响应之间存在显著相关性时，能进一步提高代理模型在超临界翼型设计中的预测精度［10］。本文将CST参数化方法、MOGP多响应代理模型与模拟退火算法相结合设计了一种翼型优化设计方法，并将该方法应用到超临界翼型的优化设计中。通过单点和多点两个翼型优化算例证明了提出方法的有效性。

1 多输出高斯过程模型

MOGP设含有D个函数的集合，｛fd（x）｝Dd＝1，其中每一个函数都可以被表示成一个光滑核函数，｛Gd（x）｝Dd＝1，和一个隐含函数u（x）的卷积，即：

更一般的，可以认为每一个光滑核函数和隐函数都是由若干个子函数线性组合而成，

每一个子函数都拥有相同的协方差函数，且相互独立。

多输出高斯过程卷积的光滑核函数和隐函数的选择有多种，但在实际应用中，一般采取高斯形式，使光滑核函数和隐函数拥有相同的高斯形式，其可行性已由Boyle和Frean证明［8］，光滑核函数表示为：

其中，Sd，q为依赖于输出端口fd（x）和隐函数μq（x）的方差系数，P－1d为输出端口的精度矩阵。

隐含过程的协方差函数可以表示为：

在上面的表达式中，Λq是关于隐函数的精度矩阵。于是我们可以得到任意两个输出端口的协方差表达式：

给定卷积形式之后，模型的似然函数表示为：

其中，y为所有输出端口的列向量；Kf，f为卷积过程的协方差矩阵，它的每个元素都是式（5）计算得到的；Σ是一个独立的高斯白噪声过程矩阵，是为每个输出端口所加的一个扰动；θ表示该模型的超参数集合，使用最大似然估计计算出θ的值，记为θ＾。

给定新的输入集合X＊后，模型输出的预测分布为：

由于式（7）的协方差矩阵Kf，f表示的是所有输出数据的协方差，所以MOGP可以为多输出建模，能够模拟多输出端口之间的相关性。

2 CST翼型参数化方法

CST方法［11－13］是近年提出的一种可以使用较少的设计变量来精确描述复杂几何外形的翼型参数化方法。该方法通过类函数定义了要描述的一般的几何外形，型函数用来描述给定类函数的特性，具有设计变量可调节、设计空间广等优点。本文采用CST参数化方法对超临界翼型进行建模。

用CST参数化方法表示翼型的数学公式为：

其中c为翼型弦长，Δzte为翼型的后缘厚度，C（x／c）为类函数，S（x／c）为型函数，分别用以下式子表示：

对于C（x／c），当N1、N2取不同的值时，可以定义不同分类的基本几何形状。本文中取N1、N2分别为0.5和1.0，方程（8）右边的第一项和第二项分别表示翼型的圆前缘和尖后缘，以保证参数化在翼型的轴向极值处具有很好的数学特性。

在方程（10）中，多项式系数bi｛｝表示了对给定翼型的设计变量，称之为Bernstein系数，用以生成更复杂更实用的翼型，可通过最小二乘法拟合翼型得到。

3 CFD分析方法

流场控制方程采用守恒形式的N-S方程。N-S方程空间离散采用有限体积方法，对流通量采用Roe通量差分格式，湍流模型采用了κ－ε模型。采用粘性网格，网格在整个流场范围内的规模为387×104个单元，图1表示了RAE2822翼型网格拓扑结构。图2表示了该网格的Y＋分布。

为了验证计算程序的准确性，图3给出了RAE2822翼型马赫数0.729，迎角2.31°状态下计算与试验压力分布对比。可以看出，CFD计算与试验结果吻合较好，能够准确捕捉到上翼面的激波。

图1 RAE2822的网格图示Fig.1 Grid for RAE2822

图2 翼型RAE2822网格的Y＋分布Fig.2 Y＋distribution of RAE2822

图3 压力分布比较Fig.3 Comparison of pressure coefficient distributions

4 模拟退火算法

模拟退火算法［14］（Simulated Annealing，SA）的思想来源于热力学中对固体退火过程的模拟，将固体加温至融化，再让其徐徐冷却凝固成晶体。加温时，固体内部粒子随者温度升高处于不稳定状态，内能增大，而温度缓慢下降时，粒子内能减小，固体内部状态趋于稳定。在实际应用中，将温度t作模拟为控制参数，将内能模拟为目标函数值f，对于控制参数t的每一个取值，算法持续进行“产生——新解——判断——接受或舍弃”的迭代过程，并逐步衰减t值，算法终止时的当前解记为所得近似最优解。其中，使用Metroplis准则作为判断当前解是否为新解替换的标准。这是基于蒙特卡罗迭代算法的一种启发式随机搜索过程。模拟退火算法适合处理传统搜索方法难以解决的高度复杂的非线性问题。本文采用模拟退火算法作为翼型设计的优化算法。

5 翼型优化系统

5.1 试验设计

试验设计的目的是选择一些有代表性的样本点从而利用少量试验点得到较高精度的代理模型。试验设计的方法有多种，但主要有正交试验设计、均匀设计、中心复合设计和拉丁超立方设计等。拉丁超立方设计［15］（Latin Hypercube Sampling，LHS）是满足空间试验设计中一种较为优秀的试验设计方法，现在已经获得了广泛的应用。因此本文选用拉丁超立方设计来选取所需要的样本点。

5.2 翼型优化设计算法框架

基于MOGP多响应代理模型的翼型优化设计方法的步骤为：首先确立使用数学语言描述的优化目标，通过CST参数化方法确定合理的翼型设计空间（设计变量的个数及其取值范围）；然后采用拉丁超立方抽样在设计空间内选取合适的初始样本点，并使用CFD软件计算出初始翼型所对应的响应值（升力系数、阻力系数等）；根据初始样本点信息构造MOGP模型；根据初始样本点构造MOGP代理模型。对代理模型进行优化分析，翼型的几何约束采用罚函数法［16］处理，即对违反约束的个体适应度施加一个惩罚项，使用模拟退火算法寻找最佳的设计变量取值，对搜索到的局部最优值，将它加入到初始样本点中，重新构造模型，如此迭代计算，直到收敛到全局最优点。图4是基于MOGP模型与模拟退火算法的优化流程图。

图4 优化流程图Fig.4 Flow chart of optimization

6 算例分析

6.1 单点翼型优化设计

算例1。以超临界翼型RAE2822为初始翼型，以面积、升力等作为约束条件，进行单目标减阻优化设计，设计状态为迎角2.79°，马赫数0.734，雷诺数6.5×106。优化问题可以描述为：

Minimize CD

Subject to CL＝0.7983

Area≥0.99Area0

在优化设计过程中，取12维Bernstein系数作为设计变量，在给定的设计变量约束范围内，采用拉丁超立方抽样选取121个初始试验点，建立初始MOGP模型，通过模拟退火算法得到单次搜索最优解并加入原样本集中，经过21次循环后得到最终最优解。表1给出了设计变量范围与最终优化结果，表2给出了基准翼型和优化后翼型的气动力参数比较。由表可知，在满足约束条件的前提下，优化后的翼型阻力系数比基准翼型减小了31%。图5和图6分别给出了原始外形与最终优化得到的翼型的形状以及压力分布的比较。由图可知，优化翼型有着较小的前缘半径，并且翼型的激波有着明显的减小。

表1 设计变量约束范围及优化结果Table 1 Scope of design variable constraint and optimization result

表2 初始翼型和优化翼型气动力参数比较Table 2 Comparison of airfoil character between initial and optimized airfoil

图5 优化前后翼型形状Fig.5 Comparison between initial and optimized airfoil configurations

图6 优化前后翼型表面压力分布比较Fig.6 Comparison of airfoil pressure distribution between initial and optimized airfoil

6.2 多点翼型优化设计

飞行器需要在一定的马赫范围内均有良好的性能。通常只有一个设计状态的单点优化很难满足要求，而多目标优化可以有效地将多个设计要求统一起来处理。这里我们对RAE2822翼型在马赫数0.70和0.75进行两点优化，并考察优化翼型在马赫数范围0.68至0.76内，阻力系数的发散情况。

对于多目标的优化问题，采用“统一目标函数法”［17－19］处理多目标函数。统一目标函数法，即将各个分目标函数统一到一个总的目标函数中去。在两点翼型优化中，假设第一设计点和第二设计点的目标函数分别为F1（x）、F2（x），使用线性插值法将各个分目标函数规格化为：

其中，αi、βi为目标函数的最大值与最小值。为使各个分目标能均匀一致的趋向各自的最优目标，根据各个目标函数的重要性分别给予一个加权因子ω，于是，统一的目标函数可表示为：

本文中，目标函数F1（x）、F2（x）分别代表待优化翼型在低马赫数与高马赫数下的气动性能参数，其权重是相同的，我们令加权因子ω1、ω2均为0.5。

算例2。设计状态（1）：迎角2.79°，马赫数0.70，雷诺数6.5×106；设计状态（2）：迎角2.79°，马赫数0.75，雷诺数6.5×106。

Minimize CD1，CD2

Subject to CL1＝0.7526

CL2＝0.7858

算例2采用了242个初始样本点构造MOGP代理模型，表3、表4给出了优化结果。经过优化后，优化翼型在两个设计状态下，阻力系数比基准翼型分别减小了5.5%和23%。

表3 设计状态1优化结果Table 3 Optimization result of design state 1

表4 设计状态2优化结果Table 4 Optimization result of design state 2

图7给出了优化前后的阻力发散特性。从图7中可以看到，由于单点优化只是局部地优化了一个设计点（Ma＝0.73）的气动性能，所以优化后的翼型虽然在设计点处和高马赫区域取得了良好的减阻效果，却较为严重地恶化了低马赫数条件下翼型的阻力特性。而多目标优化在设计点Ma＝0.70（低马赫数区域）和设计点Ma＝0.75（高马赫数区域）进行了全局的优化，所以优化后的翼型在给定的马赫数范围0.68至0.76内都有良好的减阻效果，极大地提高了阻力发散特性。

图7 翼型单点与多点优化阻力特性比较Fig.7 Drag character of airfoil at different optimization result

7 结 论

本文设计了一种基于机器学习算法和传统优化算法相结合的优化方法，并用于超临界翼型的优化设计。其核心部分为采用训练后的MOPG模型作为气动性能评估的代理模型，优化算法采用的是模拟退火算法。根据所采用的优化算例可以看出，基于MOGP模型的优化设计方法有效地改善了初始基准翼型的气动性能。从优化设计压力分布来看，所采用的优化方法能给出合理的翼型，并给出期望的压力分布。优化结果表明，所发展翼型优化设计方法可以有效地应用于超临界翼型的多目标气动优化设计，有着较好的工程应用前景。

［1］ Pena F L，Casás V D，Gosset A，et al.A surrogate method based on the enhancement of low fidelity computational fluid dynamics approximations by artificial neural networks［J］.Computers and Fluids，2012，58：112-119.

［2］Hacioglu A.Fast evolutionary algorithm for airfoil design via neural network［J］.AIAA Journal，2007，45（9）：2196-2203.

［3］ Tao Y，Li Y H，Huang Y.Multi-point optimization of supercritical airfoil based on agent model［J］.Scientific Research，System Simulation Technology ＆Application，2012，（14）：77-81.

［4］Han Z，Grtz S，Zimmermann R.Improving variable-fidelity surrogate modeling via gradient enhanced kriging and a generalized hybrid bridge function［J］.Aerospace Science and Technology，2013，25（1）：177-189.

［5］ Jeong S，Murayama M，Yamamoto K.Efficient optimization design method using Kriging model［J］.Journal of Aircraft 2005，42（2）：413-420.

［6］ Ren Q Z，Song W P.Multi-objective aerodynamic design optimization for airfoil using Kriging model［J］.Aeronautical Computing Technique，2009，39（3）：78-82.（in Chinese）任庆祝，宋文萍.基于Kriging模型的翼型多目标气动优化设计研究［J］.航空计算技术，2009，39（3）：78-82.

［7］Rasmussen C E，Williams C K I.Gaussian processes for machine learning［M］.Cambridge，MA：MIT Press，2006.

［8］ Boyle P，Frean M.Dependent gaussian processes［C］.Advances in Neural Information Processing Systems 17.Saul L K，Weiss Y，Bottou L，edit.Cambridge，MA：MIT Press，2005：217-224.

［9］ Alvarez M A，Lawrence N D.Computationally efficient convolved multiple output Gaussian processes［J］.Journal of Machine Learning Research，2011，（12）：1425-1466.

［10］Liu X，Zhu Q，Lu H.Modeling multi-response surfaces for airfoil design with multiple output gaussian process regression［J］.Journal of Aircraft，2014，51：740-747.doi：10.2514／1.C032465.

［11］Jacob C H，Willam A C，Gosset A.A parametric approach to supercritical airfoil design optimization［R］.AIAA 2007-62.

［12］Li J，Gao Z H，Huang J T，et al.Aerodynamic optimization system based on CST technique［J］.Acta Aerodynamica Sinica，2012，30（4）：443-449.（in Chinese）李静，高正红，黄江涛，等.基于CST参数化方法气动优化设计研究［J］.空气动力学学报，2012，30（4）：443-449.

［13］Kulfan B，Bussoletti J E.Fundamental parametric geometry representations for aircraft component shapes［R］.AIAA 2006-6498.

［14］Xia L，Chang Y X，Zhang L.The application of an improved simulated annealing algorithm to airfoil design［J］.Flight Dynamics，2008，26（1）：71-74.（in Chinese）夏露，常彦鑫，张龙.改进的模拟退火算法在翼型设计中的应用［J］.飞行力学，2008，26（1）：71-74.

［15］Bai J Q，Wang B，Sun Z W.The research of robust supercritical airfoil design optimization［J］.Acta Aerodynamica Sinica，2011，29（4）：459-463.（in Chinese）白俊强，王波，孙智伟.超临界翼型稳健型优化设计研究［J］.空气动力学学报，2011，29（4）：459-463.

［16］Zhang L，Chen H Q.Research on genetic algorithm for aerodynamic shape optimization based on CST［J］.Aeronautical Computing Technique，2011，41（6）：53-57.（in Chinese）张磊，陈红全.基于CST参数化的翼型优化遗传算法研究［J］.航空计算技术，2011，41（6）：53-57.

［17］Zhong X P，Ding J，Li W，et al.Robust airfoil optimization with multi-objective estimation of distribution algorithm［J］.Chinese Journal of Aeronautics，2008，21（4）：289-295.

［18］Ursache N M，Bressloff N W，Keane A J.Aircraft roll enhancement via multi-objective optimization using surrogate modeling［J］.AIAA Journal，2011，49（7）：1525-1541.

［19］Yu G，Li D.Optimization design of airfoil based on hybrid genetic algorithm and compound form method［J］.Sciece Technology and Engineering，2007，7（10）：2292-2295.

A supercritical airfoil design based on multi-output surrogate model

Wu Kuanzhan1，Liu Xuejun1，＊，Lyu Hongqiang2
（1.College of Computer Science and Technology，Nanjing University of Aeronautics and Astronautics，Nanjing 210016，China；2.College of Aerospace Engineering，Nanjing University of Aeronautics and Astronautics，Nanjing 210016，China）

An airfoil optimization method based on a MOGP（Muliple Output Gaussian Process）model and the simulated annealing algorithm is developed.The Latin hypercube method is used to construct a series of sample points in the design space，where the CST parameterization method is used to define the airfoil geometry.The initial MOGP model is established based on the response values obtained by CFD（Computational Fluid Dynamics）calculation.The proposed numerical test is to minimize the drag coefficient with consideration of some constrained conditions such as the area of airfoil，the lift coefficient and so on.Experimental results show that the drag coefficient is significantly reduced and the aera of airfoil and the lift coefficient remains，which is exactly meet the expectation of the optimization.Furthermore，the proposed method can be extended to 3Dcases straightforwardly and is applicable to practical engineering problems.

airfoil design；MOGP；CST method；simulated annealing；supercritical airfoil

V211.3

：Adoi：10.7638／kqdlxxb-2014.0059

2014-06-30；

2014-10-10

国家自然科学基金（61170152）；航空基金重点实验室类（20151452021）

吴宽展（1988-），男，硕士，研究方向：机器学习.E-mail：514400413＠qq.com

刘学军＊，教授，研究方向：计算机科学与技术.E-mail：xuejun.liu＠nuaa.edu.cn

吴宽展，刘学军，吕宏强.基于多输出高斯过程的超临界翼型优化［J］.空气动力学学报，2015，33（6）：728-732.

10.7638／kqdlxxb-2014.0059 Wu K Z，Liu X J，LV H Q.A supercritical airfoil design based on multi-output surrogate model［J］.Acta Aerodynamica Sinica，2015，33（6）：728-732.

0258-1825（2015）06-0728-06