快速交叉验证改进的运载火箭近似建模方法

文谦,杨家伟,武泽平,杨希祥,赵海龙,王志祥

国防科技大学 空天科学学院,长沙 410073

运载火箭总体设计是涵盖多学科的复杂系统工程,在其研制过程中需综合考虑气动、载荷、结构、动力等多个因素。随着计算机技术的飞速发展,智能化设计方法逐渐成为主流,也对设计效率及精细化程度提出更高要求。同时工程仿真中常用的计算流体力学(Computational Fluid Dynamics, CFD)和有限元分析在复杂构型下计算耗时长,资源消耗多。因此亟需发展一套高精度,高效率的近似建模方法来缩短整箭的研制周期。近似建模是根据离散的训练样本,采用插值或拟合的方法对模型在全变量域的输出进行预测的技术。在复杂工程问题进行分析和优化的过程中,通常采用近似模型作为原复杂耗时模型的代替,达到减少仿真次数,提升设计效率的目标。

径向基函数(Radial Basis Function,RBF)模型由Hardy在1971年首次提出,是以待测点与样本点之间的欧氏距离为自变量的函数,通过简单基函数的线性组合来模拟复杂非线性模型的响应特性,通过求解线性方程组得到基函数的系数,实现对设计空间内任意输入参数的预测。RBF近似模型由于其原理简单、操作方便,已成为众多复杂耗时的工程设计问题中广泛使用的近似模型之一。Jin等通过14个代表不同问题的算例对包括RBF在内的4种方法进行了对比研究,结果表明RBF的精度和鲁棒性最为可靠。为快速而又准确地获取可信的近似模型,通常需要对基函数的形状参数进行合理确定,以提升RBF近似模型预测精度。常用的形状参数确定方法有按经验直接确定法和优化形状参数法等。直接确定法受样本点在空间内的分布影响,Kitayama和Yamazaki提出的非均匀形状参数确定法是常用方法之一。优化形状参数法通过定义一定的精度评估指标(通常为交叉验证误差),将各基函数的形状参数作为设计变量,寻求预测精度最高的形状参数组合,能够得到最合理的形状参数取值,但由于每个基函数都有一个形状参数,随着训练样本的增多,优化问题的规模也相应增加。

交叉验证是一种统计学和机器学习领域估计模型泛化能力的方法,利用交叉验证可以解决过拟合和欠拟合的问题,避免使用同一数据集既进行模型训练,又进行泛化误差估计,能够有效合理的评估近似模型的预测精度。交叉验证方法由Wold首次提出,其基本思想是将原始数据分组,一组作为训练集训练模型,另一组为测试集来对模型做出评价,以验证模型精度和泛化能力。在模型存在未知参数时,通过优化未知模型参数,使得交叉验证误差最小化,可实现预测模型与训练数据的匹配,从而提升预测模型性能。Van Gestel等证明了将交叉验证法用于支持向量机(Support Vector Machine, SVM)和最小二乘支持向量机分类器(Least Squares SVM classifiers, LS-SVMs)的超参数选择,两者可以获得相似的结果,为交叉验证方法在近似建模领域的应用奠定了坚实的理论基础。

目前近似建模领域常用的交叉验证方法主要有简单交叉验证(Hold-out Cross-Validation, Hold-out CV)法、K-折交叉验证方法(K-fold Cross-Validation, K-fold CV)法和留一交叉验证(Leave-One-Out Cross-Validation, LOOCV)法。Hold-out CV法是将数据集分为2个子集,一个子集为测试集,而另一个子集则为训练集。该方法对数据的处理较为简单,但对于原始数据进行随机分组,得到的结果不精确,且不同分组情况所得到的结果波动很大,验证结果不可复制。为了确保验证结果可复制,降低验证过程随机性,提高样本利用率,可选用LOOCV法。LOOCV法是将数据集中的个样本轮流作为测试集,余下的-1个样本作为训练集,按照此方法将得到个模型,计算次误差。LOOCV法的优点是几乎将所有的样本用于拟合模型,得到的结果近似无偏,不受数据集划分方法的影响。该方法虽然是对样本利用率最高的方法,但其计算成本过高,面对较大的数据集时,缺点明显。通常在样本容量较小时才会选用LOOCV法。近年来国内的张英堂、刘学艺、李佳等逐渐开展了针对留一交叉验证法求解效率提升的研究。但LOOCV法仅是K-折交叉验证法的特例,面对不同问题的响应特性和不同的样本点散布,通用K-折交叉验证方法是近似建模中不可回避的问题。目前针对任意的K-fold CV法求解效率问题鲜有研究。

本文利用样本点局部密度实现将多个形状参数确定转换为缩放系数单一参数确定,建立通用快速交叉验证方法实现缩放系数的高效求解,应用于运载火箭智能设计过程中,提高近似模型求解效率及预测精度。

1 基于样本局部密度的径向基函数形状参数计算方法

1.1 径向基函数近似模型

基本径向基函数的数学模型为

(1)

式中：为待预测样本；为采样点个数；为每个基函数的权系数；()是以未知样本到已知样本的欧式距离为自变量的基函数；表示样本点距样本点的欧氏距离,即

(2)

如式(3)所示,基函数()选用高斯函数

(3)

式中：为形状参数。的值确定后,通过训练集合及其输出代入式(1)中,求解线性方程组即可计算出当前训练样本对应的基函数权系数矩阵

=

(4)

(5)

基函数系数决定后,针对任意样本,利用式(1)可以得到相应的预测输出。任意不为0的代入上述过程,均可得到经过所有训练样本的近似模型,但不合理的会导致模型泛化能力差,所以需要合理确定。本文在分析形状参数直观意义的基础上,提出了基于局部密度的形状参数归一化表征方法。

1.2 高斯径向基函数形状参数归一化表征方法

RBF的核心过程在于确定形状参数的值,直接影响数值试验的误差结果。而Rippa指出交叉验证是估算形状参数的一种有效方法,形状参数的合理选择能够有效提高近似模型的精度,减小训练样本个数,指导下一步的采样优化。李响和王晓鹏指出形状参数越大,基函数越光滑,近似模型也越光滑,有利于提高近似模型精度。Wu等的研究表明,随着值的不断增大,模型的精度呈现先上升后下降的趋势,且由于不同真实模型性质不同,因此也难以通过选择统一的形状参数保证不同近似模型的精度。

对于RBF而言,未知样本点的预测结构由所有已知样本点的输出共同影响决定,而形状参数则控制着当前样本点在设计空间内向周围影响衰减的快慢。如图1所示,形状参数越小,样本点对周围空间的影响衰减越快,反之,则对周围的影响能够传播更远。

图1 不同形状参数下Gauss基函数曲线Fig.1 Gaussian radial basis function curves with different shape parameters

样本点在空间内并不是均匀分布的,所以不能采用单一的形状参数。应根据样本点分布情况,选择合适的形状参数。在样本密集区域,应采用较小的形状参数,使每个样本的影响相对较小；在样本稀疏区域,应采用较大的形状参数,适当增大每个样本的影响范围。采用样本局部密度表征样本疏密程度,形状参数的选择就应与样本局部密度相关,样本局部密度的表达式如式(6)和式(7) 所示

(6)

(7)

式中：与形状参数类似,表示样本点对样本局部密度响应的影响程度。

在将样本归一化到单位超立方体后,通常可取

(8)

式中：为样本数量；为设计维度。

采用式(6)～式(8)计算出样本对应的局部密度,得到如图2所示结果。图中局部密度曲线表明,样本越密集,局部密度值越大,样本越稀疏则值越小。表明该方法可以合理地对离散样本点的局部密度进行量化。

理想的形状参数应在样本密集处取值相对小、样本稀疏处较大,因此得到样本点局部密度之后,将径向基函数形状参数视为该样本点影响半径,令每个样本点影响范围(影响体积)与该点局部密度成反比,即

(9)

(10)

(11)

图2 一维样本局部密度示意图Fig.2 Diagram of one-dimensional sampling points local density

1.3 缩放系数选择方法

(12)

式中：为缩放系数。

如图3所示,当取较大值时,每个基函数的影响范围被放大,适用于目标函数变化较平滑的情况；取较小值时,基函数的影响范围缩小,适用于目标函数局部变化剧烈的情况。图中实心点为采样点,细实线代表基函数,虚线区域代表基函数影响范围,粗实线是叠加后的输出预测。

图3 不同缩放系数影响范围示意图Fig.3 Diagram of influence range of different scale factors

为合理确定,建立如式(13)所示的优化问题,并采用黄金分割法等一维优化算法对其进行求解以获得最优缩放系数。

min()

(13)

式中：(·)为缩放系数对应的交叉验证误差,用来评估所构建的近似模型精度。

交叉验证基本思想是将样本集中的样本轮流作为测试样本,以评估近似模型的预示精度。相对于Hold-out CV法和LOOCV法,K-折交叉验证是最通用的方法,其基本思想是将数据集分为个子集,每个子集轮流作为测试集,其余-1个子集则作为训练集。LOOCV法是K-fold CV法中=时的特殊情况。所有的样本都作为训练集和测试集,且每个样本都被验证了一次。同时,K-fold CV法中选取不同的值可能会得到具有差异的结果。其迭代过程如图4所示。

图4 K-fold CV迭代图Fig.4 Iteration diagram of K-fold CV

将原始系数矩阵的行数据划分为组,每组内有行数据。K-fold CV法误差Δ为近似模型预测值与真实值之差,即

Δ=pre,-

(14)

式中：pre,为针对-个训练样本构建的近似模型在个测试样本处的预测值；为测试样本真实输出。

(15)

由式(15)可见,在K-fold CV中每折都需要训练一个新近似模型,即对一个-阶矩阵进行求逆,整个过程需求解个-阶矩阵的逆。和越大、越小,计算速度越慢。线性方程组式(15) 的求解占据了交叉验证过程大部分的计算时间,消耗了大量的计算资源。

为提升K-fold CV运算效率,节约计算成本,本文提出了逆矩阵快速求解算法和交叉验证误差快速求解算法,有效解决了反复求解高阶矩阵的逆的问题。采用K-fold CV进行缩放系数的确定,由于一次交叉验证误差评估需要对训练集进行次近似建模,增加了RBF近似建模的计算耗时。针对该问题,提出快速K-折交叉验证(Fast K-fold Cross-Validation, FK-fold CV)方法以提升近似建模效率与精度。

2 快速交叉验证方法

2.1 逆矩阵快速求解算法

针对K-折交叉验证计算效率低下的问题,本文提出的FK-fold CV法通过矩阵初等变换与分块矩阵求逆来改进常规K-fold CV法,降低计算复杂度,提升直接求解交叉验证误差的效率。

如图5所示,K-fold CV法通过将样本划分为组,每次取一组样本进行测试,剩下的所有样本作为训练样本,循环次之后即得到近似模型在每一组样本点上的误差。虚线框中将第组与第组交换的步骤是FK-fold CV法中独有的。

图5 快速K-折交叉验证误差求解示意图Fig.5 Flow chart of FK-fold CV

FK-fold CV法以K-fold CV法为基础,将数据集分为个子集,求解每一子集的交叉验证误差时,通过初等变换将测试集置于训练数据集的末端,变换过程如图6所示。

图6 FK-fold CV迭代图Fig.6 Iteration diagram of FK-fold CV

通过左乘变行,右乘变列的原则对原系数矩阵进行行列交换,将应作为测试集的第组与第组换位后,得到新的系数矩阵。(,)定义为阶单位阵第组与第组交换后得到的矩阵,与关系如式(16)所示。

(,)(,)=

(16)

(17)

系数矩阵被划分为组后,设当前测试组内有个样本,将第组与第组数据交换位置,式(17)矩阵中的每个元素代表一个×的小矩阵,例,表示数据集中第组数据。对式(16) 求逆得

(18)

由(,)(,)=,可以得到式(19)。

(19)

()=

(20)

(21)

、、、满足如下数学关系时,可以计算得出

(22)

其中:是(-)阶方阵;是(-)×阶矩阵;是×(-)阶矩阵;是阶方阵。采用式(23)即可直接计算出(-),运用阶矩阵求逆替代了(-)阶矩阵直接求逆的过程,由于≪,因此转化后的矩阵求逆计算量会大幅降低。

(23)

式中：-表示矩阵去掉最后行、列所得矩阵的逆矩阵。

上述推导表明,当LOOCV法是K-fold CV法的特殊情况,即等于样本点的数量,每次取1个样本进行测试,其余作为训练样本集。当采用留一法进行推导时,式(21)中为一常数,此时近似模型的系数矩阵的逆矩阵可通过式(24) 计算

(24)

2.2 交叉验证误差快速求解算法

根据式(15)可得,FK-fold CV法误差Δ为

Δ=pre,-=,-(-)--

(25)

式中：,-为系数矩阵中测试样本对应的行和前-列元素；-为训练样本预测输出。由式(23)和式(25),可得式(26)

(26)

(27)

式中：为所有样本的真实输出,且测试样本的个输出对该项没有任何贡献；,为系数矩阵中测试样本对应的行元素。式(27)展开得

(28)

根据逆矩阵定义有

(29)

(30)

因此,式(28)转化为

(31)

根据RBF的定义可知

(32)

将式(31)、式(32)代入式(25),可得

(33)

式中：为测试样本点对应的个权系数。

显然,采用FK-fold CV法计算预测误差,仅需要对小矩阵进行求逆,极大地节省了计算量。当每折样本个数为1时,交叉验证误差按照以上过程推导得

(34)

本文针对LOOCV法计算误差的推导结果,Rippa所推导的快速误差计算公式是本文取时的特殊情况,表明本方法具有更强的通用性。在计算时间上,对采样点个数为的近似建模问题,普通K-折交叉验证的时间复杂度为,快速K-折交叉验证的时间复杂度为()。采用500个样本点的数值算例分别对K-fold CV法和FK-fold CV法求解效率进行测试,在不同折数下时间节省幅度如图7所示。折数越多时间节省幅度越大,尤其在与常规留一交叉验证对比时,节省了超过90%的时间,算例表明快速K-折交叉验证具有更高的求解效率。

图7 时间节省幅度Fig.7 Graph of time saving percentage

3 数值算例校验

选取4个多峰非线性函数验证本文近似建模方法的效率、精度和稳定性,算例形式如表1所示,对应的函数图像如图8～图11所示。

图8 测试函数1图像Fig.8 Graph of test function 1

图9 测试函数2图像Fig.9 Graph of test function 2

图10 测试函数3图像Fig.10 Graph of test function 3

图11 测试函数4图像Fig.11 Graph of test function 4

为验证本方法的有效性,将本文方法与Kitayama、Rippa等提出的近似建模方法以及Kriging方法在相同条件下进行性能对比。采用4种近似模型训练方法分别针对4个测试函数进行训练,选择样本预测输出与样本真实输出的均方根误差(Root Mean Square Error, RMSE)来表征近似模型的精度。RMSE指标定义如下：

(35)

表1 数值测试函数Table 1 Numerical test function

为避免样本选取带来的影响,针对随机取样的情况独立运行30次对方法性能进行测试。基于快速交叉验证法选取最优缩放系数,对于一维函数在设计域内随机选取20个样本点,对于二维函数在设计域内随机选取40个样本点,通过样本点局部密度计算不同样本点相对形状参数大小,然后采用一维搜索法针对缩放系数进行寻优,得到最优缩放系数及其对应的预测模型。之后采用相同的样本点,用其他3种方法构造近似模型并计算对应的RMSE,不同算法得到的均方根误差箱型图分别如图12～图15所示,其中测试函数3、4由于其误差较大,将其结果取对数处理后,再用于表征其趋势。结果表明,基于FK-fold CV法构建的近似模型在4个测试函数中,都展示出比其他3种方法更优越的性能。

本文近似模型构建算法每一次测试都能达到与已知方法相近甚至更优的性能。针对测试函数1、2,如图12、图13所示,本文算法展现出针对复杂多峰问题近似建模的潜力,可以看出Rippa提出的FLOOCV法在处理多峰值问题时,没有考虑样本在设计空间内的分布,仅使用单一的形状参数对所有基函数进行统一处理,其精度不够理想。对于测试函数3,在处理角落处有尖峰出现的二维多极值问题时,Kitayama方法和Rippa方法均体现出一定的局限性。当测试函数存在大量极值时,例如测试函数4,Rippa方法的性能较差。

图12 测试函数1均方根误差箱型图Fig.12 Boxplot of RMSE of test function 1

图13 测试函数2均方根误差箱型图Fig.13 Boxplot of RMSE of test function 2

图14 测试函数3均方根误差箱型图Fig.14 Boxplot of RMSE of test function 3

图15 测试函数4均方根误差箱型图Fig.15 Boxplot of RMSE of test function 4

4 快速交叉验证方法工程应用

针对工程仿真中常用的计算流体力学(Computational Fluid Dynamics, CFD)和有限元分析模型,以运载火箭减阻特性和和结构承载能力近似建模为例,验证本文方法预测精度。

4.1 超声速减阻杆阻力特性近似建模

针对火箭飞行过程中存在的激波阻力问题,常在火箭头部安装减阻杆以进行流场重构,减小气动阻力。减阻杆作为被动减阻技术的一种,具有结构简单、减阻效果好等优点,在实际工程系统中得到广泛应用,如图16、图17所示的三叉戟潜射洲际弹道导弹及Igla单兵肩扛式防空导弹均采用减阻杆降低飞行器激波阻力。

图16 三叉戟潜射洲际弹道导弹Fig.16 UGM-133A ballistic missile

图17 Igla单兵肩扛式防空导弹Fig.17 Igla air defense missile

减阻杆的减阻效果受到包括减阻杆顶端半径及减阻杆杆长等外形尺寸参数的影响,本节将采用FLUENT对圆盘形减阻杆与球形头部结合减阻杆的减阻效果受到包括减阻杆顶端半径及减阻杆杆长等外形尺寸参数的影响,本节将采用FLUENT对圆盘形减阻杆与球形头部结合后外形的阻力系数进行数值计算,几何外形如图18 所示。设计参数为减阻杆顶端半径和减阻杆杆长,见图19。飞行器头部半径=100 mm,减阻盘厚度=10 mm。数值模拟中边界条件设置为压力远场条件及壁面无滑移绝热壁条件,来流马赫数=2,气体模型设置为完全气体模型,气流参数设置为标准海平面大气参数,湍流模型选用SST-模型。由于该几何外形为轴对称外形,采用轴对称条件进行计算可提高计算效率,因此选用轴对称条件进行数值模拟。

图18 超声速飞行器外观及减阻结构示意图Fig.18 Diagram of supersonic vehicle appearance and drag reduction structure

图19 减阻结构设计参数示意图Fig.19 Diagram of drag reduction structure design parameters

和的取值范围如表2所示,采用10×10个均匀分布的样本点作为测试集,进行数值模拟得到100 组不同顶端半径及杆长的减阻杆与球形头部组合外形的阻力系数。运用本文所提出的方法进行近似建模,实验设计方法生成20个点作为快速交叉验证的数据集。选用快速5-折交叉验证,每折样本容量为4,在单次验证过程中,16个为训练数据,4个为测试数据。

表2 减阻结构设计参数Table 2 Drag reduction structure design parameters

4种方法分别计算得出的均方根误差值如图20 所示。在样本点不均匀的2维气动工程问题中,由于Kitayama法中形状参数为一个固定值,不会根据实际情况进行调参,其所构建近似模型精度缺点明显,误差达到了FK-fold CV法的7倍左右。

图20 工程算例1均方根误差值Fig.20 RMSE of engineering example 1

4.2 重型运载火箭级间段加筋壳结构承载能力近似建模

采用复杂加筋壳结构设计问题对本文方法的有效性进行进一步验证。运载火箭及导弹舱段等航天结构中常用的薄壁加筋柱壳结构,具有较高的轴压、弯曲和扭转承载效率的特点。作为运载火箭的主要承力部段,薄壁加筋柱壳的轻量化设计可以在提高其运载能力的同时降低成本。蒙皮桁条结构是典型框桁加强薄壁加筋柱壳结构之一,轴压承载性能优越,广泛应用于运载火箭结构设计。

图21 加筋柱壳结构及中间框布局Fig.21 Diagram of cylindrical stiffened shell structure and intermediate frame

上述复杂结构设计工程问题具有19维设计变量,各参数变化范围如表3所示,输出响应为结构轴向承载极限能力。采用优化拉丁超立方得到设计空间内均匀分布的400个样本点,作为训练样本来构建近似模型。在选择缩放系数时,采用快速10-折交叉验证,每折中应有40个数据。在单次交叉验证过程中,360个样本点为训练数据,40个样本点为测试数据。同样采用优化拉丁超立方得到设计空间内均匀分布的5 000个样本点,仿真得到的极限载荷作为真实输出,计算出真实输出与近似模型预测输出的RMSE值,再对比本文方法与常用近似建模技术所得结果的误差大小。

表3 蒙皮桁条结构设计参数Table 3 Skinned purlin structure design parameters

4种方法分别计算得出的均方根误差值如图22 所示。在高维结构问题中,Kitayama法与Kriging法都有一定的局限性,其误差分别为FK-fold CV法的194%和158%。

图22 工程算例2均方根误差值Fig.22 RMSE of engineering example 2

5 结 论

1) 提出基于样本局部密度的形状参数表征方法,将多个形状参数确定转化为缩放系数确定的单变量优化问题,降低径向基函数形状参数确定的计算复杂度,解决现有RBF形状参数确定中待求变量多和计算复杂等问题。基函数中形状参数的合理确定大幅提升了RBF近似模型预测精度。采用4种方法分别对运载火箭进行气动力和结构承载能力近似建模,结果表明采用本文方法在复杂工程问题上预测精度远高于现有方法且相对稳定,通用性强且能满足工程实际需求。

3) 由于Rippa提出的FLOOCV法是本文所提出的FK-fold CV法的特例,2种方法思路大致相同但本方法普适性更强。在针对上述算例进行对比时发现,这2种方法应用于低维测试函数时差异较为明显,所以面向不同问题时,应对的取值进行讨论。对于FK-fold CV算法本身而言,综合时间与近似模型精度考虑,选取的最优值,是下一步值得探索的问题。