在工程优化设计领域，目标函数的获得往往十分耗时，发展高效的优化方法(即减少目标函数的计算次数)是缩短设计周期的关键。此外，为了提高设计效果，优化方法应尽可能地具有好的全局性，即应尽可能地找到全局最优解。
基于代理模型的优化方法被认为是一类能够实现全局优化的高效优化方法。Kriging模型具有较好的预测非线性、多峰值的能力，且能预测出未知点的误差，因而成为最受欢迎的代理模型之一，并被越来越多地应用于多学科优化设计中[1, 2, 3]。尽管如此，Kriging模型的效率和精度有时仍然不能满足工程实际的需求，这是由于在很多实际问题中，建立一个具有合理精度的近似模型仍然需要大量的高精度数值模拟，而这一问题将随着设计变量的个数和范围的增加变得更加突出。为此，近些年来，人们正寻求比传统代理模型更为高效的新型代理模型。
提高代理模型效率的一种方法是使用大量低可信度样本数据来辅助高可信度样本数据来建立一种用于预测高可信度分析程序输出结果的变可信度模型[4, 5, 6, 7]，如Cokriging模型。该方法的思想是在设计空间内同时采用高可信度分析程序(如Navier-Stokes方程程序)和低可信度程序(如Euler方程程序)进行抽样，通过大量的低可信度样本数据获得设计空间内响应值的全局趋势，通过少量高可信度样本数据来校正预测值。
另一种提高Kriging模型效率的方法是在Kriging模型中引入梯度信息，即在建立代理模型时不仅利用样本点处的函数值，同时利用样本点处的梯度值。引入梯度的Kriging模型被称为梯度增强的Kriging模型(gradient-enhanced Kriging,GEK)。Liu等[8]提出了一种间接GEK模型。Chung等[9]提出了一种直接GEK模型。Laurenceau等[10]比较了Kriging、直接GEK、间接GEK模型在气动问题中的预测精度，发现在相同样本点数下2种GEK的精度都高于Kriging模型。Han等[1]发展了一种更一般化的GEK模型，允许梯度信息在任意位置引入，且允许在某些样本点处不使用梯度信息。
目前，对于Kriging模型及基于该模型的优化方法已得到大量研究，对于GEK模型本身的研究也不少见[9, 10, 11, 12]，而鲜见对GEK模型和Kriging模型在应用于优化设计时的优化效率和优化效果的比较仔细的对比研究。为了探讨GEK模型在优化设计中的效率及效果，本文首先在Han等[11]发展的GEK模型基础上，结合多种加点准则[13]及多种传统优化算法，发展了一套基于代理模型的高效优化设计框架；其次，在统一的优化框架下，采用典型数值算例将GEK模型和Kriging模型进行了对比研究，最后采用工程中的气动设计问题将二者进行了对比。结果显示，对于大多数问题，使用GEK时收敛更快且优化结果更好。
1 Kriging模型/梯度增强的Kriging模型对于Kriging模型的详细介绍，本文不再赘述，见相关参考文献[1, 14, 15]。GEK模型是在Kriging模型基础上发展起来的，其思想、原理和建模方法与Kriging模型基本一致，只是在模型中不仅利用了样本点处的函数值，还利用了样本点处的梯度值。下面介绍本文所采用的由Han等[11]发展的一种直接GEK模型。
首先,在Kriging模型中,对于一个m维问题,假设对目标函数(或约束函数)y抽样获得n个函数值,抽样位置及相应的响应值可写成
对于GEK,将上述样本数据的响应值列向量加以拓展,把样本点处的梯度信息以偏导数的形式加入到响应值中。于是得到含梯度信息的样本及响应值矩阵:
式中,xg(j),j∈｛1,…,n′｝表示第j个提供偏导数的样本点(每个偏导数看做一个样本点),向量yS中的
表示第j个偏导数信息(对第k个独立变量的偏导数,k∈｛1,…,m｝)。
GEK模型对未知函数的预测值定义为所有抽样函数值和所有抽样偏导数值的加权,即
式中,w(i)代表第i个函数值的加权系数,λ(j)代表第j个偏导数信息的加权系数。与Kriging模型类似,引入静态随机过程假设:
式中，β0是未知的常数,Z(·)为均值为0的静态随机过程,且协方差满足关系式
式中,σ2 是随机过程Z(·)的方差,R是仅与空间距离有关的“相关函数”,
和
分别为相关函数相对于x(i) 和x(j) 处第k个分量的偏导数。
与Kriging模型推导过程类似,可得到GEK模型的预测值:
式中,R为样本点相关矩阵,r为样本点与未知点x的相关向量:
均方误差MSE由下式计算:
1.2 相关函数模型
相关矩阵R和相关向量r的建立需要计算相关函数及其一阶导数、二阶导数值。此处考虑以下形式的相关模型:
本文采用与文献[11]类似的3次样条相关函数:
1.3 模型参数的获得
由(8)式可知,GEK模型中存在超参数θ,与Kriging模型相同,θ值采用最大似然估计法求解下列优化问题获得:
2 多点加点准则及其约束处理
基于代理模型的优化方法的一个关键技术是加点准则。由于建立一个较精确的全局近似模型需要大量的样本点,特别是设计变量较多时。因此,为了减少总的样本点数,首先使用较少的样本点建立初始代理模型,再寻找一定准则下的最优点并将其作为新的样本点来更新代理模型重新优化。加点准则的使用直接决定了优化的成败及优化效率。本文发展一种多点加点准则:同时采用局部收敛性好的目标函数最小值准则和全局收敛性好的EI最大准则,并分别将其推广到含约束的形式。
2.1 目标函数值最小准则及其约束处理该准则认为代理模型足够精确,直接寻找代理模型上目标函数的最小值,即建立代理模型后,求解下列子优化问题:
式中，NG是约束个数。本文使用遗传算法(GA)和梯度优化相结合的方法求解最小化问题。
该准则没有考虑模型的误差,新加的样本点主要集中在当前真实最优解附近,因而它是一种局部加点准则。本文称之为minimizing the constrained prediction,简称CMP。
2.2 EI最大值准则及其约束处理(CEI)设样本集中的最优目标函数值为ymin,在x处目标函数相对的改进量为I=ymin-Y(x)>0,于是I的期望值可由下式计算:
式中,φ表示标准正态分布概率密度函数,Φ表示标准正态分布累积函数。
当存在约束g(x)>0时,约束函数也建立代理模型,并也假定G(x)服从正态分布,G(x)∈
。于是G(x)满足约束的概率如下:
满足约束的EI值(CEI)如下:
若存在多个约束,则将(12)式中乘以每一个约束函数满足约束的概率。
该准则以当前真实最优解为参照,同时考虑了代理模型的误差,因而是一种全局加点准则。
上述2种加点准则的结合兼顾了优化的局部收敛性和全局收敛性,因而可以更好地适用局部优化问题和全局优化问题。
3 基于代理模型的优化设计方法本文结合代理模型、试验设计方法(本文采用LHS)、多点加点准则、遗传算法、梯度优化算法等发展了1套基于GEK模型的优化设计方法。该方法的流程如下:
1) 使用试验设计在设计空间内生成一些初始样本点。
2) 调用真实分析模型,获得每个样本点的真实目标函数值和约束函数,对于GEK模型,还需计算二者的梯度值。
3) 针对目标函数及每一个约束函数,分别建立目标函数和约束函数的代理模型。
4) 分别使用CEI和CMP进行子优化,获得代理模型上的2个最优解。
5) 检查2个新的最优点是否重合或与已有的样本点出现重合,若重合,剔除新出现的重叠点。
6) 调用真实函数校验两种加点准则所得的最优点的真实目标函数值和约束函数,对于GEK模型,还需计算二者的梯度。
7) 检验优化是否达到终止标准,若是,停止优化,否则将上步中所得最优点作为新样本点加入样本点集,回到步骤3)。
4 算例测试在本文优化算例中,为了与Kriging模型进行对比，假设计算1次目标函数的梯度所需的计算量等同于1次目标函数的计算量(类似在气动优化设计中,若使用Adjoint方法[12]计算梯度,其计算量大致等于1次流场求解的计算量。由于初始样本点的生成具有随机性,且在子优化过程中使用了具有随机性的遗传算法,因此,本文基于代理模型的优化算法也会存在随机性。因此本文对所有算例重复优化30次来进行研究。
4.1 无约束算例(局部优化)该算例是求解Rosenbrock函数的最小值。优化模型如下:
其真实最优解为x*=(1.0,…,1.0),f(x*)=0。图 1给出了Rosenbrock函数的云图,图中的“+”为最优点。
 |
| 图 1 Rosenbrock函数云图 |
分别使用Kriging模型和GEK模型进行优化,在优化过程中,使用LHS方法生成4个初始样本点。使用GEK模型进行优化时,最大样本点数取100,加上100次梯度所需的计算量,总计200个计算量。在使用Kriging模型进行优化时,最大样本点数取200。图 2给出了使用Kriging模型和GEK模型运行30次得到的平均收敛曲线的比较。可以看出,在相同的样本点数情况下,使用GEK模型优化的结果明显优于Kriging模型优化的结果。然而,计算梯度需要额外的计算量,在相同计算量情况下,使用GEK模型优化的结果好于Kriging模型优化的结果。表 1给出了使用GEK模型、Kriging模型优化结果的比较。作为参照,表中还给出了直接使用GA优化得到的结果。从表 1还可以看出使用2种代理模型方法的优化结果都明显优于GA所得的结果,尽管GA调用了10 000次目标函数。(表中“计算量”一列中“f”表示真实目标函数计算，“g”表示真实目标函数的梯度计算作为参照。)
 |
| 图 2 使用Kriging和GEK优化Rosenbrock函数的目标函数平均收敛曲线对比 |
| 方法 | 最好值 | 平均值 | 标准差 | 计算量 |
| 真实值 | 0 | ／ | ／ | ／ |
| GEK | 1.429 3×10 -7 | 5.1020×10 -5 | 9.706 1×10 -5 | 200(100f+100g) |
| Kriging | 1.655 2×10 -7 | 9.2790×10 -4 | 4.324 4×10 -3 | 200 |
| GA | 8.230 8×10 -5 | 4.552 9×10 -2 | 5.842 9×10 -2 | 10 000 |
4.2 无约束算例(全局优化)
本算例是一个全局优化算例,目标是求解Rastrigin函数的最小值,优化模型如下:
函数的真实最优解为x*=(0,…,0),f(x*)=0。图 3给出了Rastrigin函数图,可以看出,该函数是一个复杂多极值函数。
 |
| 图 3 Rastrigin函数 |
使用LHS方法获得10个初始样本点。使用GEK模型进行优化时,最大样本点数取150;使用Kriging模型进行优化时,最大样本点数取300。图 4给出了使用Kriging模型和GEK模型运行30次得到的平均收敛曲线。可以看出,基于Kriging模型的优化方法和基于GEK模型的优化方法都表现很好。然而,GEK模型优化的收敛精度要高于Kriging模型优化的收敛精度。表 2给出了使用本文GEK方法、基于Kriging模型的优化方法优化结果的比较。该表也反映了图 4所显示的结果。作为参考,表中也给出了直接使用GA优化得到的结果。
 |
| 图 4 使用Kriging和GEK优化Rastrigin函数的目标函数平均收敛曲线对比 |
| 方法 | 最好值 | 平均值 | 标准差 | 计算量 |
| 真实值 | 0 | ／ | ／ | ／ |
| GEK | 1.579 2×10 -12 | 1.008 8×10 -10 | 1.697 8×10 -10 | 300(150f+150g) |
| Kriging | 4.502 5×10 -11 | 2.440 9×10 -8 | 7.944 2×10 -8 | 300 |
| GA | 0 | 3.552 7×10 -16 | 1.945 9×10 -15 | 30 000 |
本算例是一个压力容器设计问题。该优化问题的数学模型如下所示:
图 5给出了2种代理模型方法重复优化的平均收敛曲线的对比。表 3给出了3种方法分别优化30次所得最好最优值、平均最优值、最优值的标准差,以及计算量的比较。与上例类似,GEK模型几乎每次都精确地收敛到相同的值,30次重复优化最优值的平均值和标准差都小于Kriging,这说明GEK模型效果好于Kriging且表现更稳定。该表也给出了直接使用GA优化得到的结果,可以看出,使用基于GEK和Kriging的优化方法所得的优化结果都明显优于GA,而计算量却远小于GA。
 |
| 图 5 使用Kriging和GEK优化pressure vessel design的目标函数平均收敛曲线对比 |
| 方法 | 最好值 | 平均值 | 标准差 | 计算量 |
| 真实值 | 未知 | ／ | ／ | ／ |
| GEK | 5 885.333 | 5 885.336 | 0.005 | 100(50f+50g) |
| Kriging | 5 885.789 | 5 902.444 | 19.731 | 100 |
| GA | 5 910.793 | 6 671.307 | 623.203 | 40 000 |
为了比较2种代理模型应用于工程优化设计的的优化效率和效果,以翼型压力反设计进行测试。
设计状态为Ma=0.73,Re=6.5×105,α=2.0°。首先给定RAE2822翼型在该设计状态的压力分布,以该压力分布作为目标压力分布,将反设计问题转化为优化问题求解,优化问题的数学表达式如下:
式中,p为计算压力,pd为目标压力,B表示沿物面的积分。
翼型的几何参数化采用“CST”参数化方法[17],18个设计变量。流场求解采用课题组发展的求解器PMNS2D;计算网格采用C型结构化网格。目标函数的梯度使用Adjoint方法计算。因此,通过求解一次RANS方程和一次与之计算量相当的Adjoint方程即可得到目标函数及其梯度值。
分别使用GEK模型、Kriging模型进行翼型的反设计,重复运行30次。作为参考,这里还给出了使用基于Adjoint的梯度优化算法所得的结果,其中梯度优化算法采用BFGS拟牛顿法,初始翼型为NACA0012翼型。图 6给出了基于代理模型的优化方法的平均收敛曲线及Adjoint方法收敛曲线的对比。从图 6可以看出,使用Kriging时花费300个样本点已基本收敛,而GEK方法仅花费50个样本点(100个计算量)目标函数就已收敛,优化结果也好于使用Kriging模型。
 |
| 图 6 基于GEK方法和Kriging方法的优化方法的翼型反设计平均收敛曲线及Adjoint翼型反设计收敛曲线对比(一次Adjoint方程计算等同于一次流场计算) |
图 7给出了分别使用GEK、Kriging所得的某次反设计结果对应的翼型几何形状及压力分布与目标翼型几何形状及压力分布的对比。可以看出,使用2种模型反设计所得的翼型的几何形状及压力分布与目标翼型的几何形状及压力分布都吻合得很好,但是仔细比较,GEK表现比Kriging更好,这可在图中圆圈部分可以看出。
 |
| 图 7 使用GEK和Kriging模型反设计前后翼型的压力分布及几何形状与目标翼型的对比 |
在气动设计领域，由于当前已具备高效计算梯度的技术——Adjoint方法，因此，将目标函数与其梯度相结合建立代理模型来进行优化设计，可能有着很大的潜力。为了研究梯度的引入对优化设计的影响，本文首先将Kriging模型、梯度增强的Kriging模型、多点加点准则、多种传统优化算法相结合，发展了一套基于代理模型的高效、通用优化设计框架。其中，本文发展了一种多点加点准则，并将加点准则推广到含约束的形式。再次，基于相同的优化设计框架，对使用梯度增强的Kriging模型和使用Kriging模型的优化结果进行了比较。在假设梯度的计算花费与目标函数的计算花费相同的情况下，对于绝大多数问题，在总的计算量一致的前提下，使用GEK时收敛更快且优化结果更好。此外，还使用工程实际中的翼型压力反设计对Kriging和GEK的效果进行了对比，结果显示，当使用GEK模型时，优化结果和优化效率都高于使用Kriging模型。
| [1] | Hoyle N, Bressloff N W, Keane A J. Design Optimization of a Two-Dimensional Subsonic Engine Air Intake[J]. AIAA Journal, 2006, 44(11):2672-2681 |
| Click to display the text | |
| [2] | Ahmed Mym, Qin N. Surrogate-Based Multi-Objective Aerodynamic Design Optimization of Hypersonic Spiked Bodies[J]. AIAA Journal, 2012, 50(4):797-810 |
| Click to display the text | |
| [3] | Yim J, Lee B J, Kim C. Exploiting Multi-Stage Shape Optimization Strategy of Multi-Body Geometries Using Kriging-Based Model and Adjoint Method[J]. Computers & Fluids, 2012, 68:71-87 |
| Click to display the text | |
| [4] | Kennedy M C, O'Hagan A. Predicting the Output from a Complex Computer Code When Fast Approximations Are Available[J]. Biometrika, 2000, 87(1):1-13 |
| Click to display the text | |
| [5] | Forrester A I J, Soberster A, Keane A J. Multi-Fidelity Optimization via Surrogate Modeling[J]. Proceedings of the Royal Society A:Mathematical, Physical and Engineering Sciences, 2007, 463:3251-3269 |
| Click to display the text | |
| [6] | Han Z H, Zimmermann R, Goertz S. An Alternative CoKriging Model for Variable-Fidelity Surrogate Modeling[J]. AIAA Journal, 2012, 50(5):1205-1210 |
| Click to display the text | |
| [7] | Han Z H, Goertz S. A Hierarchical Kriging Model for Variable-Fidelity Surrogate Modeling[J]. AIAA Journal, 2012, 50(9):1885-1896 |
| Click to display the text | |
| [8] | Liu W, Batill S M. Gradient-Enhanced Response Surface Approximations Using Kriging Models[R]. AIAA-2002-5456 |
| Click to display the text | |
| [9] | Chung H S, Alonso J J. Using Gradients to Construct Cokriging Approximation Models for High-Dimensional Design Optimization Problems[R]. AIAA-2002-0317 |
| Click to display the text | |
| [10] | Laurenceau J, Sagaut P. Building Efficient Response Surfaces of Aerodynamic Functions with Kriging and CoKriging[J]. AIAA Journal, 2008, 46(2):498-507 |
| Click to display the text | |
| [11] | Han Z H, Goertz 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 |
| Click to display the text | |
| [12] | Jameson A. Aerodynamic Design via Control Theory[J]. Journal of Scientific Computation, 1988, 3(3):233-260 |
| Click to display the text | |
| [13] | Jones D R. A Taxonomy of Global Optimization Methods Based on Response Surfaces[J]. Journal of Global Optimization, 2001, 21:345-383 |
| Click to display the text | |
| [14] | Han Z H, Zhang K S, Song W P, et al. Optimization of Active Flow Control Over an Airfoil Using a Surrogate-Management Framework[J]. Journal of Aircraft, 2010, 47(2):603-612 |
| Click to display the text | |
| [15] | Han Z H, Zhang K S, Liu J, et al. Surrogate-Based Aerodynamic Shape Optimization with Application to Wind Turbine Airfoils[R]. AIAA-2013-1108 |
| Click to display the text | |
| [16] | 许瑞飞,宋文萍, 韩忠华. 改进Kriging模型在翼型气动优化设计中的应用研究[J]. 西北工业大学学报, 2010,28(4):503-510 Xu R F, Song W P, Han Z H. Application of Improved Kriging-Model-Based Optimization Method in Airfoil Aerodynamic Design[J]. Journal of Northwestern Polytechnical University, 2010, 28(4):503-510 (in Chinese) |
| Cited By in Cnki (19) | Click to display the text | |
| [17] | 卜月鹏,宋文萍,韩忠华,许建华. 基于CST参数化方法的翼型气动优化设计[J]. 西北工业大学学报,2013,31(5):829-836 Bu Yuepeng, Song Wenping, Han Zhonghua, et al. Aerodynamic Optimization Design of Airfoil Based on CST Parameterization Method[J]. Journal of Northwestern Polytechnical University, 2013, 31(5):829-836(in Chinses) |
| Cited By in Cnki (0) | Click to display the text |