基于多保真度代理模型的装备系统参数估计方法

刘 佳, 杨克巍, 姜 江, 夏博远

(国防科技大学系统工程学院, 湖南 长沙 410073)

0 引 言

随着科技水平的不断提高和信息技术的快速发展,为传统的武器装备系统设计和研制提出了以下3个新要求。

(1) 系统的复杂性指数级增高。复杂性是复杂系统的典型特征,随着系统和项目复杂性的增加,对系统工程的需求也随之增加。由于技术的快速发展、知识的快速更新和环境的快速变化,复杂性以指数形式增加了系统组件之间冲突的可能性,因此也增加了设计的不可靠性。在这种情况下,复杂性不仅包含工程系统,还包含逻辑上的人工数据组织。同时,由于尺度的增加以及设计中涉及的数据、变量或字段数量的增加,系统可能变得更加复杂。

(2) 需求的动态性变化。在传统的系统工程研究过程中,期望在系统分析阶段对所有需求进行尽可能完善的分析,定义和表征系统和子系统及其之间的相互作用是系统工程的目标之一。智能化技术的引入,改变了原有系统工程的固有流程,使得系统建设各个阶段的需求呈现出动态性特征。一方面,系统的复杂性增加,使得系统分析阶段难以对所有的需求进行充分分析,指数增长的复杂性使得系统方案需要不断进行调整,引发需求的动态变化;另一方面,与传统的瀑布或“V”型研发模式相比,新的螺旋式模式更强调在不断开发过程中系统性能的改进,从而降低不断增加的分析评估所消耗的资源与时间成本,因而产生了在改进过程中系统需求的持续性变化。此外,由于环境和技术的快速变化,使得系统的目标和任务具有相对不确定性,由此产生的需求变化也具有持续、动态等特点。

(3) 响应的快速性要求。智能化技术的发展和应用推广,要求大幅缩短经典系统工程的许多进程周期,尤其是针对动态性变化的多种需求,从系统整体性的角度而言,也要求系统表现出任务的快速响应特征。首先,要利用智能化技术使得系统开发中的原型系统设计和构建的进程加快,从而可以快速响应动态变化的系统需求。其次,要求缩短测试周期,使得系统方案能够得到快速的验证和评估。第三,由于信息技术的深度融合带来的系统综合程度高度集成,对经典系统工程基于文档的开发模式产生了颠覆,比如现在提出的将武器装备“试验鉴定评估”技术“左移”的策略,都是对传统系统工程设计开发模式提出改进要求的表现。

现有的与武器装备系统设计有关的文献研究大都是围绕模型驱动的武器装备体系层面展开的。文献[1]通过对作战能力和作战效能进行区分从而提出了一种武器装备体系作战能力评估的框架。文献[2]基于体系贡献率对武器装备体系的评估方法展开了详尽综述,强调评估装备对体系贡献的力量。文献[3-4]分别从能力需求的视角对武器装备体系结构的建模和评估提供了一种新思想。文献[5]详细介绍了武器装备效能评估指标,并对其关系进行了分析,从而构建装备的指标体系。文献[6]通过构建任务-能力之间的对应关系,对面向使命任务的武器装备体系能力间的关联关系进行了分析。文献[7-9]分别就武器装备试验领域中的作战试验方案设计、试验鉴定指标建立和评估体系构建展开了研究。从上述文献可以看出,目前相关研究大致可以分为以下几类:武器装备体系的结构建模、方案设计、指标体系构建、能力评估等。

但是为了应对上述3个智能化要求,基于模型的系统工程(mode based systems engineering, MBSE)方法必须融入数据的作用，不能仅聚焦于对整个武器装备体系层面的分析和研究,要下沉到单装备系统的总体设计上来。一方面,对系统底层设计参数进行估计和把握,能够更加准确地得到该装备对能力需求的满足度;另一方面,系统的总体设计不再停留在概念设计阶段,将论证往下更加推进一步,落到系统初步设计上,可以为系统后面的详细设计提供更加有力的论证依据。

武器装备系统的参数估计是复杂系统总体设计中的关键一环,同时参数估计的效率和精确度对系统的总体设计产生重要的影响。综述近三年国内外有关系统设计参数估计的文献,一部分是对输入输出的数值噪声进行处理从而实现估计结果精度的提高,如文献[10-11]中使用的卡尔曼滤波方法。其中绝大部分文献都是借助现有计算能力的迅速提升而开展的纯数据驱动的系统参数估计方法,通过对大量输入数据的预处理、使用智能优化算法构建参数估计模型、对算法进行改进从而提高模型的估计精度,如文献[12]中使用的是多目标遗传算法,文献[13]中的研究是基于粒子群算法,甚至文献[14]中引入了强化学习的方法。这些方法都是依赖于大量数据集的支撑,而本文的研究重点是考虑到在武器装备系统总体设计初期,缺少充足数据的支持,尤其是对于精确度高的数据,实验成本高且获取难度大。针对上述问题,本文开展了对基于多保真度代理模型的单装备系统参数估计的方法研究,旨在解决面对该种数据缺乏的情况,仍能保证系统参数估计的精确度和快速响应的要求。

1 多保真度代理模型

1.1 代理模型

代理模型[15]是一种用来获取一组独立变量与系统响应之间某种近似关系的统计技术,是试验设计、数理统计和最优化技术的一种综合应用。其基本思想是用一个简单的函数关系近似替代真实曲面函数,从而提高计算效率。设计变量和响应输出的关系可以表示为

(1)

代理模型技术的实质是以拟合精度或预测能力为约束,利用近似技术对离散数据进行回归或插值的数学模型,通过有限的已知点响应构造近似函数表达式对未知区域进行预测。通常用来近似替代在分析和优化设计过程中比较复杂和费时的数值分析的数学模型,图1为代理模型的构建过程。

图1 代理模型构建过程

目前常用的代理模型近似技术包括:响应面模型(response surface model, RSM)[16]、径向基函数(radial basis function， RBF)[17]、Kriging插值法[18]、人工神经网络(artificial neural networks, ANN)[19]、支持向量机(support vector machine, SVM)[20]等。

1.2 试验设计方法

试验设计[21]作为代理模型技术的重要一环,为代理模型技术提供了科学、经济的试验方案,使样本点能够按照不同的要求分布在参数设计空间中,更为有效地反映系统输入参数与输出响应之间的复杂函数关系。构建代理模型的第一步是选择样本点[22],即需要用试验设计的方法来决定初始化空间内样本点的位置。试验设计方法的主要目的在于对整个设计变量空间进行样本点的高效选取,使有限的样本点尽可能反映出设计变量空间变化特性[23]。进一步说,就是通过选取最少的样本点,使获取的关于未知设计空间的信息最大化。目前的试验设计方法大致分为两类:经典试验设计方法和现代试验设计方法[24]。经典试验设计方法包括全因子设计、中心组合设计、D-最优设计、Box-Behnken设计等。这些方法主要用于实验室安排仪器实验中,以减少实验随机误差。现代试验设计方法以正交试验设计、蒙特卡罗抽样、拉丁超立方设计、均匀设计为主要代表,是为确定性的计算机实验发展需要而产生的,主要采用“空间填充”的思想。目前,常用的现代试验设计方法有:蒙特卡罗抽样、正交试验设计、均匀设计和拉丁超立方抽样[25](Latin hypercube sampling,LHS)。

本文设计的多保真度代理模型方法拟采用的试验设计方法是LHS方法。LHS[26]是一种在现代试验设计中比较流行的能够对大型参数设计空间中进行均匀高效采样的试验设计方法,受到基于计算机仿真分析设计领域的广泛研究与关注。LHS方法约束随机地生成均匀样本点,通过控制样本点的位置,避免抽样点在小领域内重合。LHS的基本原理是如果在设计空间内抽n个样本点,那么就把m个随机设计变量的抽样范围都分成等距离或等概率的n个区间,然后对于每个变量分别在每个区间随机取一个值,这样对于每个变量就有对应的n个水平值,然后将m组的n个变量值随机组合配对就构成了n个样本点。对于每个设计变量来说,n个样本点一定分别落在其每个小区间中,因而得到的实际抽样点就会等概率地分散在整个设计空间中。LHS法不仅具有有效的空间填充能力,还能够拟合非线性响应。与其他现代试验设计方法中另一常用的正交试验法相比,在同样的采样点数目下LHS设计可以研究更多的样本点组合[27]。LHS方法不限制问题的维数、样本数目的多少,因为LHS方法产生的样本点可以确保其代表向量空间中的所有部分并且具有相当大的随意性。同时,LHS方法还具有样本记忆功能,抽样效率高,能够避免重复抽取已经抽过的样本点,并且对于分布在抽样空间边界处的样本点也能保证使其参与抽样[28]。因此，拉丁超立方试验设计方法可以保证在抽样较少的情况下获得较高的计算精度。

1.3 多保真度问题描述

多保真度问题可以描述为:对于一个m维问题,要想获得一个难以评估的函数(计算成本大或者函数响应值难以获取)y1的预测值,可以通过借助一个易于评估的函数(计算成本低或者函数响应值容易得到)y2的响应结果,以及函数y1上的几个样本点值进行预测。首先需要对不同保真度的样本空间(高保真度样本空间S1,低保真度样本空间S2)进行采样,得到

S1=(x1,1,x1,2,…,x1,n1)T∈Rn1×m

(2)

S2=(x2,1,x2,2,…,x2,n2)T∈Rn2×m

(3)

式中，n1是高保真度数据的采样数量;n2是低保真度数据的采样数量。在实际情况中,精确度高的高保真度采样点的响应值计算复杂性高、计算耗时长导致难以获取,甚至在有些情况下受客观原因的限制无法获得大量的高保真度数据,而低保真度采样点的响应值没有精度要求的约束,其获取的途径有很多,比如通过经验统计、基本原理公式、快速建模仿真软件等[29],所以获取成本相对较低,容易获取。由此可见,高保真度数据的获取难度远远高于低保真度数据的获取难度,因此通常情况下,n2≫n1。

接着,通过独立的不同复杂性的计算途径分别获取高、低保真度的采样点响应值:

y1=[y1(x1,1),y1(x1,2),…,y1(x1,n1)]∈Rn1

(4)

y2=[y2(x2,1),y2(x2,2),…,y2(x2,n2)]∈Rn2

(5)

因此,该问题的输入就是低保真度数据集(S2,y2)和高保真度数据集(S1,y1),输出目标是得到高保真函数y2。有了该函数,给出任意参数样本点,都能快速得到该样本点对应的高保真响应值。通过上述方法得到的函数响应值不仅计算成本低、响应快、而且数据精度高。

2 基于多保真度代理模型的参数估计方法

2.1 多保真度代理模型构建流程

多保真度代理模型的构建流程如图2所示。

图2 多保真度代理模型构建流程

具体可分为以下6个步骤。

步骤 1使用试验设计方法从初始化参数空间中进行低保真度数据的样本采样。

步骤 2通过低保真度数值分析获得低保真样本点的响应值,从而得到一组低保真数据集。

步骤 3基于低保真度数据集构建低保真代理模型。

步骤 4再次使用试验设计方法从初始化参数空间中进行高保真度数据的样本采样。

步骤 5通过使用复杂的高保真数值分析获得高保真度数据样本点的响应值,从而获得一组高保真度数据集(高保真数据集的数量远低于低保真数据集)。

步骤 6通过输入一组高保真度数据集,在步骤3构建的低保真度代理模型的基础上建立多保真度代理模型。该多保真度代理模型就可以用来预测初始化参数空间中任意一点的响应值,快速且精确。

2.2 基于原始Kriging构建低保真度代理模型

Kriging[30]是一种基于统计理论的插值技术,是以已知样本信息的动态构造为基础充分考虑到变量在空间的相关特征,建立对象问题的近似函数关系来模拟某一点未知信息的代理模型方法,是一种基于随机过程的估计方差最小的无偏估计模型。Kriging模型与其他常用代理模型的区别在于,该模型不仅能给出对未知函数的预估值,还能给出预估值的误差估计。Kriging模型对非线性函数具有良好的近似能力和独特的误差估计功能。

Kriging模型是一种插值模型,其插值结果是由已知样本函数响应值的线性加权[31],即

(6)

式中,ω(i)是加权系数。

基于式(6)的描述,只要给出加权系数ω=[ω(1)，ω(2)，…，ω(n)]T的表示式,就可以得出设计空间中任意样本点的预估值。为了求出加权系数,传统Kriging模型引入统计学假设[32]:将未知函数看成是某个高斯静态随机过程的具体实现,即对于任意位置x,对应的函数响应值y(x)被一个随机函数Y(x)替代,y(x)只是Y(x)可能的结果之一。该静态随机过程写成:

Y(x)=β0+Z(x)

(7)

式中,β0是全局趋势模型,表示Y(x)的数学期望值,是一个未知常数;Z(·)是均值为0,方差为σ2(σ2(x)≡σ2),协方差为cov[z(x),z(x′)]=σ2R(x,x′)的静态随机过程;R(x,x′)是空间相关函数，该函数只与空间中两点之间的欧式距离有关。

传统Kriging模型要寻找最优加权系数ω(i),须使得均方差(mean squared error,MSE)最小[29],即

(8)

同时根据无偏估计要求,还须满足:

(9)

经求解,得到传统Kriging模型的预估值为

(10)

式中,

β0=(FTR-1F)-1FTR-1ys

(11)

同时也能得到传统Kriging模型预估值的MSE估计为

(12)

对于一个有m个设计变量的优化问题,假定需要建立未知性能函数y:Rm→R(目标函数或约束函数)对设计变量x=[x1,x2,…,xm]T∈Rm的近似模型。

基于上述问题描述,首先需要选择如拉丁超立方采样等试验设计方法,在设计空间中进行抽样得到n个样本:

S=[x(1),x(2),…,x(n)]T∈Rn×m

(13)

这n个样本点也可以看作是n个设计方案。接着,需要对这n个设计方案进行数值分析(计算流体力学或者计算结构力学),从分析结果中得到n个函数响应值:

ys=[y(1),y(2),…,y(n)]T=[y(x(1)),y(x(2)),…,y(x(n))]T∈Rn

(14)

要为高保真函数建立代理模型,首先需要为低保真函数建立一个代理模型,用于后面的辅助预测。根据前文对传统Kriging模型方法的介绍,将对低保真函数的随机过程响应值定义为

Ylf(x)=β0,lf+Zlf(x)

(15)

式中,β0,lf是未知常数;Zlf(x)是静态随机过程。按照前文介绍的建立Kriging模型的方法,基于低保真数据集(Slf,ys,lf)建立传统Kriging模型。在拟合出Kriging模型后,由于F=[1,1,…,1]T,可参照式(10)将低保真函数在任意未知点x的预估值写为

(16)

式中,

(17)

其中，Rlf∈Rnlf×nlf是由所有已知样本点之间的函数值组成的“相关矩阵”;I∈Rnlf是单位列向量;rlf∈Rnlf由未知点与所有已知样本点之间的相关关系组成的“相关向量”。

2.3 基于改进Kriging构建多保真度代理模型

(18)

(19)

参照传统Kriging模型的求解过程,可以将改进Kriging模型的预估值写为

(20)

除此以外,改进Kriging模型还给出了预估值的MSE:

(21)

不论是在建立传统Kriging代理模型还是改进Kriging代理模型的过程中,代表已知样本点函数之间相关性和已知样本点与未知样本点之间相关性的“相关矩阵”和“相关矢量”的构造都与“相关函数”的选择和计算有关。这里考虑到满足高斯假设和计算简便的要求[34],建议采用一类三次样条函数,其形式为

(22)

式中,

(23)

假设采样数据是根据高斯过程分布的,采样点的响应值被认为是相关的随机函数,其对应的似然函数为

(24)

比例因子β0和过程方差σ2的最优估计值可解析求出,但均由未知超参数θ=[θ1,θ2,…,θm]决定。

β0(θ)=(FTR-1F)-1FTR-1ys

(25)

(26)

将式(25)和式(26)代入到式(24)中,取对数,优化问题就转化为求解:

ln[L(θ)]=-nlnσ2(θ)-ln|R(θ)|

(27)

由于无法解析求出θ的最优值,就需要采用数值优化算法,模型参数θ的求解目标为

(28)

由此可以得出,θ是该代理模型优化的关键参数。θ最优值选取的好坏会直接影响该模型的效果和精度。考虑到遗传算法与其他的一般优化算法的不同之处,在于其他算法往往是先选定一个初始值,然后再进行迭代,最终选出最优解,但是这样得到的最优解,其最优性很有可能被限制在局部范围内,而遗传算法会预先选择一个初始值集合,对一群初始值进行迭代,最终再选出最优解,这种选择方法相对其他算法在最优解的选择范围上更大,更有利于寻找全局最优解[35]。因此，建议使用遗传算法对模型参数θ最优值进行求解。由于现有文献中对遗传算法的基本原理和算法流程都有详尽地说明,这里就不对其进行展开叙述。

3 应用研究

为了验证所提出的系统参数估计方法,应用研究选择飞机这一单航空装备系统设计为示例进行分析。首先给出模型的数据输入,这里包括低保真度和高保真度两种不同保真度数据,接着借助建立的多保真度代理模型对采样点的响应值进行求解,然后对求解结果进行说明分析,最后引入了几个模型评价指标进一步地验证该系统参数估计方法的有优越性。

飞机系统总体设计的主要参数包括翼载荷、推重比和最大升力系数[36]。由于翼载荷和推重比主要是通过对飞机自身重力、机翼面积和发动机种类的选择决定的,因此该示例着重研究最大升力系数这一主要参数。飞机的升力系数主要由飞机的机翼提供,即要求飞机的升力系数就是求飞机机翼的升力系数,而机翼的升力系数是通过对机翼的结构参数进行设计来改变的,在机翼面积、展弦比、后掠角等一般参数确定的情况下,机翼弯度的改变将对机翼气动特性产生重大的影响[37],机翼弯度是产生升力的最基本要素。飞机机翼的弯度是指机翼剖面的弯曲程度,具体是指中弧线与翼弦之间的距离[38],平直翼的机翼弯度一般为零。综合以上分析,将本示例研究的自变量定为机翼弯度,因变量是最大升力系数。

3.1 模型数据输入

示例研究实验中,首先需要对低保真样本点进行抽样,采用的试验设计方法是LHS方法,样本机翼弯度变量的抽样范围是(20 cm, 40 cm),采样点数为10。由于LHS方法是归一化的,将抽样结果转化为对应的抽样区间,采样得到的10个点是38 cm, 27 cm, 29 cm, 22 cm, 36 cm, 33 cm, 25 cm, 31 cm, 20 cm, 35 cm。然后,本文将使用建模仿真工具作为低保真数值分析程序,计算时间很快,得到一组低保真点及其响应值(机翼弯度和及其对应的最大升力系数),具体数值如表1所示。

表1 低保真度样本点及响应值

同样,采样LHS方法在区间(20 cm, 40 cm)中进行高保真样本点采样,由于获取高保真样本点响应值的成本高,因此只选取5个采样点2 cm, 25 cm, 30 cm, 33 cm, 38 cm。本文获取高保真点响应值的计算程序是求飞机气动特性的专业求解器,由于这种求解器一般都不是开源的,且计算耗时较长。得到一组高保真样本点及其响应值(机翼弯度及其对应的最大升力系数),具体数值如表2所示。

表2 高保真度样本点及响应值

3.2 模型求解

本文实验主要分为两步,首先是输入经过采样得到的10对低保真数据组,构造低保真度代理模型;再输入经过采样得到的5对高保真数据组,在前一步构造的低保真度代理模型的基础上得到最终的多保真度代理模型。通过两次数据输入,完成两部分实验,可以得到基于保真度代理模型的参数预测结果。本文的研究中是想要求得机翼弯度在任意一点对应的最大升力系数响应值,因此通过运行编写的算法程序,可以求解机翼弯度在(20 cm, 40 cm)内每一点的响应值,代码运行时间即结果计算时间是秒级,可见算法结果响应十分迅速。为了将实验结果进行对比,同时还分别采用低保真数值分析程序和高保真数值分析程序将机翼弯度在(20 cm, 40 cm)中的其他点的响应值补充完整,这里基本可以将高保真数值分析得到的响应结果真实值来进行对比。这样就获得了3组机翼弯度在(20 cm, 40 cm)区间中的全数据集:低保真度数据、高保真度数据和通过本文算法实验得到的多保真度数据,每组数据的自变量都是机翼弯度,因变量是最大升力系数。因此，以机翼弯度为横坐标,最大升力系数为纵坐标作图,得到图3中的3条曲线。其中,黑色圆点代表的是低保真度数据结果,橙色三角形代表的是高保真度数据结果,蓝色方块就是本文实验的预测结果,即多保真度数据结果。图3中5个醒目的绿色方表示第二步实验中输入的5对高保真度数据。这5个点多保真度预测曲线和高保真度数据结果曲线的交点,并且观察图3可见横坐标即机翼弯度从30 cm开始到40 cm,多保真度代理模型的预测结果与高保真度数据结果十分接近,两条曲线的拟合效果很好。同时,从图3中也能够很明显地看出由于大部分低保真数据结果都是基于数值统计或者经验公式等粗糙的数值分析程序快速求得的,导致与真实数据(即高保真度数据结果)相比,差距略大。

图3 不同保真度模型的估计值结果

3.3 结果评价与分析

显然,本文研究问题想要得到的结果估计模型是一个回归模型,而不是一个分类模型,因此本文选择均方根误差(root mean square error, RMSE)、最大差值(max error, MAX)和R-平方值这3个常用的回归模型性能评价指标[39]作为评价本文求解结果模型好坏的指标。

(29)

(30)

(31)

以高保真度数据为真实数据作对照(即RMSE=0,MAX=0,R2=1),分别计算低保真度数据模型和多保真度代理模型的上述3个指标值,得到结果如表3所示。

表3 模型性能指标结果值

由表3的结果可以看出,多保真度代理模型的4项指标的表现,无论是整体精度还是局部精度都是优于低保真度数据的。尤其是多保真度预测模型的R-平方值为0.986,几乎接近1,说明模型拟合效果很好,而低保真度数据的R-平方值只有0.086,拟合效果不尽人意。

4 结 论

针对由信息技术的快速发展给复杂系统设计和研制带来的一系列智能化要求:复杂性激增、需求动态性变化、响应的快速性;并且,系统参数估计又是武器装备系统设计研制中的重要一环,同时又考虑到不同保真度数据获取难度差异的情况,本文提出了一种基于多保真度代理模型的参数估计方法。该方法首先基于传统Kriging模型对低保真度数据建立代理模型,然后将该模型作为全局趋势模型,再基于改进的Kriging模型插入高保真度数据建立多保真度代理模型,得到同时满足高精度和高效率双重要求的参数估计方法。最后以飞机机翼系统设计中的最大升力系数参数估计问题为例,验证了所提出的方法在装备设计研制中的可行性和有效性。

值得一提的是,虽然低保真度数据计算方便且快捷,相对于高保真度数据获取难度占明显优势,甚至有时不仅仅是因为计算耗时成本高,还包括一些客观原因致使高保真度数据的获取难度大且获取数量少,因此想要获取全样本点的高保真度响应值耗费巨大成本,在有些情况下是不可能完成的。但是通过本文给出的基于多保真度代理模型参数估计方法,所需输入的数据获取难度仅仅等于多个低保真度数据获取难度加上少量高保真度获取难度(重点是可实现的且可行性高),就可以获取十分接近真实数据的精度和拟合效果,可见性价比十分高。同时也为一些特殊情况:无法获得采样空间中所有点的高保真度数据响应值但是想要获取高精度的结果响应值提供了一种解决办法,即揭示了在实际情况下用来指导武器装备总体设计的作用和意义,为装备系统设计人员在研制总体设计方案时面对动态性的需求提供了可操作、科学性且能够快速响应的方法支撑。