概率盒全局灵敏度和活跃子空间跨层降维

2021-10-20胡政文张保强邓振鸿

航空学报 2021年9期

胡政文，张保强，邓振鸿

厦门大学航空航天学院，厦门 361000

由于制造、测量、计算及模型本身中存在的各种误差，航空航天系统参数化建模和多学科优化需要考虑不确定性的影响[1-4]。在不确定性分析时，不确定性参数按照类别大致可分为随机不确定性参数和认知不确定性参数[5]。不同类型不确定性参数混合作用与单种不确定性参数作用导致的模型输出不确定性存在明显的差异，因此模型分析过程中须同时考虑随机和认知不确定性的影响[6-8]。Williamson和Downs[9]在累积分布函数的基础上引入概率盒来处理随机和认知不确定性混合量化问题。概率盒可以让随机和认知不确定性在模型中分别传播[10]，并保留2种不确定性各自的特性，因此在混合不确定性量化分析中得到了广泛应用[11-13]。

同时由于涉及多层次和多学科，航空航天系统参数的维数通常较多[14]。学科间耦合作用下系统通常呈现高度非线性，难以直接判断输入和输出之间的关系，高维参数给系统研究分析带来了挑战[15]。因此系统建模前通常须考虑减少输入参数的维数，以往研究中主要有两种思路处理维数简化问题：一种是参数选择，这种思路认为在所有的输入参数中只有少数真正与模型输出相关，其余参数对模型输出不确定性影响较小可以忽略，参数选择通常采用灵敏度分析的方法[16]；第二种是降维方法，与参数选择相比，降维方法假设所有输入参数都影响模型的输出结果，但只在参数的某些特定组合下才能最大程度体现出来，而降维方法的目标就是构造这类组合[17]。

灵敏度分析主要研究输入不确定性对输出不确定性的贡献大小。随着模型复杂性的增加，输入参数的灵敏度信息逐渐成为建模时需要考虑的因素之一[18-20]。灵敏度分析方法主要分为局部灵敏度分析和全局灵敏度分析，全局灵敏度分析方法主要有Sobol’指标、矩灵敏度指标等[21-22]。传统的灵敏度分析方法大都基于概率框架，然而实际中模型通常包含认知不确定，须同时考虑区间的灵敏度分析[21, 23]。Oberguggenberger等[24]阐述了概率和非概率框架下的灵敏度分析方法。近年来，随着概率盒在混合不确定性量化中的应用，概率盒下的全局灵敏度研究越来越多。Ferson和Troy Tucker[25]研究了在混合不确定性概率边界分析时通过缩减输入参数的不确定性进行灵敏度分析。Song等[26]提出一种扩展蒙特卡罗模拟方法，计算概率盒下各输入的基于方差的全局灵敏度指标。Bi等[27]使用Bhattacharyya距离定量描述概率盒下的不确定性参数灵敏度。不确定性缩减法分析概率盒全局灵敏度具有易于理解、工程中实现较为容易的优势，逐渐得到研究人员的关注[28]。

另一方面，降维方法在航空发动机、卫星、汽车、芯片、图像处理等许多领域中得到了广泛应用[14]。Fodor[29]回顾了常用的降维方法，包括主成分分析、因子分析、随机投影法、连续拓扑映射、回归等。主成分分析是普遍使用的线性降维方法之一[30-31]，该方法基于参数的协方差矩阵，可以有效地降低均方误差，也被称为本征正交分解、Karhunen-Loève (KL) 变换、霍特林变换和经验正交函数法等[32-33]。Constantine等[34]通过构造输出的梯度协方差矩阵来判断输入参数空间中的主要方向，提出了活跃子空间降维方法。活跃子空间法推广了传统的主成分分析思想，对处理高维复杂系统具有一定的优势[35]。Lukaczyk等[36]将活跃子空间法应用于ONERA-M6跨声速机翼的形状设计优化中，在设计空间区域内找到了最佳设计变量。Jefferson等[37]将活跃子空间法应用在综合水文模型的灵敏度分析。Hu等[38]使用活跃子空间法对卫星系统进行多学科优化设计。

以往处理高维参数的简化问题的研究中，一般使用灵敏度分析或降维方法两者之一，对于高维混合不确定性系统单种方法降维效果不明显。本文将灵敏度分析和降维方法结合在一个框架中，提出一种结合概率盒全局灵敏度和活跃子空间的跨层降维方法。案例分析中将跨层降维方法应用到NASA多学科不确定性量化挑战问题中，验证提出的跨层降维方法的有效性。

1 跨层降维方法基本理论

1.1 概率盒全局灵敏度分析

概率盒定义为一组包含变量所有可能累积分布函数(Cumulative Distribution Function, CDF)曲线的边界，可以有效表示混合不确定性影响下的参数特征，图1为典型概率盒示意图。概率盒不确定性分析提供了一种较为直观的方法来描述受到随机和认知不确定性综合作用的输出特性。

图1 典型概率盒示意图Fig.1 Schematic of typical probability box

不确定性缩减法可以用于概率盒框架下不确定性参数的全局灵敏度评估。输入参数的灵敏度取决于参数自身不确定性的大小以及对模型输出不确定性的影响，概率盒全局灵敏度可以表示为[25]

(1)

式中：B为原始基准输入参数；T为缩减不确定性后的输入参数；unc(·)表示输出不确定性的度量;S为全局灵敏度指标。若unc(·)为标量，则灵敏度指标S也是标量，由此可以对输入参数进行灵敏度排序。输入参数不确定性缩减后，式(1)中分子unc(T)的值会下降。与基于方差的Sobol’指标不同，由于不确定性的减少，对于所有输入参数，unc(T)均小于unc(B)，因此灵敏度指标S范围是[0,1]。对各输入参数依次应用不确定性缩减法，由式(1)可以得到各输入参数的S值，将所有的S按照大小进行排序即得到对应参数的灵敏度排序。

输入参数不确定性的缩减有多种方法，不同缩减方式得到的不确定性度量unc(·)的结果不同。通常有以下2种方式：① 使用特定的概率分布替代不确定性输入参数，特定的概率分布法只消除参数的认知不确定性，不影响其随机不确定性；② 使用固定值替代不确定性输入参数，固定值法同时消除参数的随机不确定性和认知不确定性。

不确定性度量unc(·)也有多种不同的定义方法(如方差、矩等)，不同的定义方式可以用于解决不同的灵敏度分析问题[25]。概率盒是参数所有可能累计分布函数的包络线，沿着概率盒的上界和下界进行积分可以得到概率盒的面积。本文结合概率盒的特征，采用概率盒的面积作为输出的不确定性度量unc(·)，此时不确定性缩减法的流程见图2，式(1)转化为

图2 概率盒灵敏度分析的不确定性缩减法流程图Fig.2 Flow chart of pinching method for probability box sensitivity analysis

(2)

式中：A(·)表示概率盒面积。

不确定性缩减法主要比较缩减输入参数不确定性后输出概率盒的面积占比下降值，进而判断参数的灵敏度信息。

1.2 活跃子空间法

假设原始输入参数的维数为n，表示为x=[x1,x2,…,xn]T，降维方法的目的是在满足损失数据信息最少的条件下，通过某种方法将原始的n维参数x用m维的s=[s1,s2,…,sm]T(m

(3)

对于参数x，其标准参数空间可以视为边界为1的超立方体。非标准的参数空间可以通过正则化转化为标准参数空间。假设在x的参数空间中随机抽样，样本数为N，则输入参数x可以表示为矩阵X={xi,j:1≤i≤n,1≤j≤N}。对于任意的xi,j，可正则化为

(4)

目标函数f梯度的协方差矩阵定义为

(5)

(6)

C=WΛWT

(7)

式中：W为n×n的特征向量矩阵；Λ=diag(λ1,λ2, …,λn)是非负特征值组成的对角矩阵。特征值大小反映了目标函数f在x参数空间中沿着对应特征向量微小扰动的平均变化。将特征值按照降序进行排列λ1≥λ2≥…≥λn，特征值越小表示对应的特征向量在梯度协方差矩阵C中的贡献度越低，降维时可以适当舍弃这些向量。

活跃子空间确定基向量维数有多种方法：通过设定固定阈值λ0，只保留大于λ0的特征值及对应的特征向量；定义κ为前m项特征值占所有特征值的累加比例：

(8)

(9)

(10)

式中：U为活跃子空间的基向量，包含W的前m列；Λ1为对应的特征值对角矩阵。确定活跃子空间的基向量后，原输入的参数空间可以投影到该活跃子空间：

s=UTx

(11)

此时目标函数f(x)可以近似表示为活跃子空间内降维变量s的函数：

f(x)≈g(UTx)=g(s)

(12)

g(s)的定义域为γ={s=UTx,x∈Rn}∈Rm。在降维变量s的基础上构建低维代理模型：

g(s)≈g*(s)≡R(s;g1,g2,…,gM)

(13)

式中：M为代理模型训练样本数；g*(s)≡R为活跃子空间中定义的代理模型，可以降低计算模型的复杂性和减少计算时间和成本。原函数和代理模型的误差可以表示为[34]

(14)

系数C1取决于x和ρ(x)。由式(14)可以看出，舍弃的特征值之和越小，即累加比例κ的值越大，原函数和新函数误差也越小，即活跃子空间降维法所产生的误差也就越小。若前m项特征值累加比例κ>90%，则可以选择前m项特征向量组成基向量用于降低输入参数空间的维度。

1.3 跨层降维方法

在不确定性缩减法和活跃子空间的基础上，提出结合概率盒全局灵敏度和活跃子空间的跨层降维方法。

图3是跨层降维方法的流程图，基本步骤如下：

图3 跨层降维方法流程图Fig.3 Flow chart of cross-layer method for dimension reduction

步骤1整理系统的不确定性参数，将其按照随机、认知和混合不确定性进行分类。

步骤2使用概率盒法对初始模型进行不确定性分析，得到模型输出的初始概率盒，计算概率盒的面积。

步骤3缩减待分析的输入参数的不确定性，对模型进行分析得到缩减不确定性后输出概率盒的面积。

步骤4比较缩减参数不确定性前后输出概率盒的面积改变量，计算参数的灵敏度指标S值。

步骤5选取另一输入参数，重复步骤3和步骤4，直至完成所有不确定性参数的灵敏度分析。

步骤6对各输入参数的S值按照降序排列，将S值较小的输入参数固定。

步骤7对保留不确定性的输入参数进行抽样，取样本数N，计算模型输入并得到梯度协方差矩阵。

步骤8对梯度协方差矩阵进行特征分解，按照特征值大小对特征值进行排序，确定保留的活跃子空间的维数和基向量。

步骤9对输入参数使用活跃子空间的基向量进行降维。

2 NASA多学科不确定性量化挑战问题

为了验证提出的概率盒全局灵敏度和活跃子空间的跨层降维方法，本文以NASA Langley研究中心提出的多学科不确定性量化挑战问题为案例[40]。

2.1 挑战问题描述

挑战问题的数学模型描述了地面基站无线电波操控遥控飞机的通用传输动力学特性，如图4所示[40]。要求对不确定性系统模型进行综合评估，解决不确定性量化、全局灵敏度分析、极限状态分析和稳健优化设计中的一系列关键问题。

图4 NASA多学科不确定性量化挑战问题模型示意图[40]Fig.4 Model schematic of NASA multidisciplinary uncertainty quantification challenge problem[40]

模型中的不确定性参数按照类型划分为随机、认知以及混合不确定性参数3类。Ψ表示所研究的系统模型，p表示Ψ中不确定性输入参数(p∈R21)，d表示设计参数(d∈R14)，g表示Ψ性能约束(g∈R8)，g取决于p和d。不确定性输入参数p、设计参数d和约束g之间的关系可以表示为

x=h(p)

(15)

g=f(x,d)

(16)

式中：x是参数为p的中间函数(x∈R5)；g和x假定为连续可微函数。式(15)展开为

(17)

x的每个分量都表示系统一个学科的输出，x1～x5均为标量。若系统Ψ满足对于g的所有分量均满足gi≤0,i=1,2,…,8，则可以认为该系统是可靠的，符合设计要求。对于确定设计参数d，满足约束向量g≤0的输入参数p的集合称为安全域，反之称为失效域。g的分量有一个不满足约束条件，对应的参数p就落在失效域内。图5为挑战问题系统模型Ψ中各参数变量的层次关系示意图。

图5 NASA挑战问题模型参数层次关系Fig.5 Hierarchical relationship of model parameters in NASA challenge problem

2.2 全局灵敏度分析

由于系统模型Ψ的输入参数中包含随机、认知以及混合不确定性，首先对Ψ进行概率盒不确定性分析，研究系统的不确定性传递和量化的相关问题，确定不确定性输入参数p对学科输出x的全局灵敏度，以便后续活跃子空间降维过程中固定灵敏度低的参数。

输入参数p的不确定性通过函数h传播，得到x1、x2、x3、x4和x5的概率盒，缩减p中的认知和混合不确定性参数时，对应x概率盒的面积也会改变。

由式(17)可知，子学科Ψ1的输出x1是p1、p2、p3、p4和p5的函数，表1给出了参数p1～p5的具体分布信息，其中p1为单峰Beta分布族组成的混合不确定性参数，p2是区间不确定性参数，p3服从[0,1]均匀分布，p4和p5为相关的二元正态分布族组成的混合不确定性参数。表中:Δ表示参数的区间范围;E(·)、V(·)和ρ分别代表均值、方差和相关系数。随机不确定性是系统模型的固有属性，对模型输出影响较大，为减少计算成本，只对认知和混合不确定性参数按照S值大小排序。

表1 x1对应不确定性输入参数p1～p5

对于子学科Ψ1，对初始条件下的x1进行概率盒不确定性分析。取1 000条x1的样本累积概率分布曲线，得到初始条件下x1的概率盒，如图6 所示，此时x1概率盒的面积为0.229。

图6 初始条件下的x1概率盒Fig.6 Probability box of x1 under initial conditions

使用不确定性缩减法分别缩减参数p1、p2、p4、p5的不确定性。由于p1为混合不确定性参数，为计算p1对x1的概率盒全局灵敏度，不确定性缩减法的固定值取为p1的期望值0.7。p2为认知不确定性参数，对应的缩减值取为区间中值0.5。p4～p5服从二维正态分布，相关系数ρ的绝对值在0～1范围内，在计算p4的灵敏度时，取p4为期望值0，p5仍然服从正态分布，得到p4的灵敏度，同理计算p5的概率盒全局灵敏度。综合考虑计算精度以及效率成本，结合x1概率盒边界情况，每个参数均取1 000条累积分布函数曲线，概率盒边界趋于稳定，进而得到缩减不确定性后各参数对应的x1概率盒，如图7所示。运用式(2)计算得到各参数的灵敏度指标S值，计算结果列于表2之中。由图7和表2可以看出，子学科Ψ1的参数灵敏度排序依次是p1>p5>p4>p2。由灵敏度指标S值的计算结果，p2和p4对x1概率盒的灵敏度均小于0.1，可以将p2和p4视为固定的参数。

表2 x1概率盒的面积及对应输入参数灵敏度排序

图7 缩减各参数不确定性对应的x1概率盒Fig.7 Probability boxes of x1 after pinching corresponding uncertain parameters

其余子学科x2、x3、x4和x5对应的16个不确定性输入参数p6～p21的具体分布信息如表3所示。

表3 x2～x5对应不确定性输入参数p6～p21Table 3 Description of uncertain input parametersp6-p21 corresponding to x2-x5

表4 x2概率盒的面积及对应输入参数灵敏度排序

表5 x3概率盒的面积及对应输入参数灵敏度排序

表6 x4概率盒的面积及对应输入参数灵敏度排序

图8 缩减参数p21不确定性对应的x5概率盒Fig.8 Probability boxes of x5 after pinching p21

综合概率盒灵敏度分析的计算结果，筛选出p1、p5、p6、p7、p12、p16、p17、p18、p21共9个影响程度较大的认知或混合不确定性参数以及p3、p9、p11、p19共4个随机不确定性参数。

2.3 活跃子空间降维

根据全局灵敏度分析的结果，共有13个影响较大的不确定性参数，应用概率盒全局灵敏度和活跃子空间的跨层降维方法，采取活跃子空间进一步降低不确定性参数的维度。

为描述系统模型的可靠性，定义如下系统性能度量函数：

(18)

性能度量函数ζ描述了系统模型最容易超出约束的g分量，间接反映了系统失效的最大可能性。为构建合适的活跃子空间基向量，将ζ作为目标函数。由拉丁超立方抽样选取300个输入参数样本，分别计算对应目标函数ζ的值，并根据式(6) 构造ζ的近似梯度协方差矩阵C。

对梯度协方差矩阵进行特征分解，将特征值按照从大到小排序。表7列出了C的13个特征值，图9为特征值随着维度的衰减曲线图。由表7可知，矩阵C的一维对应的特征值为0.063 84，远远大于二维的特征值1.400×10-17，且一维特征值占所有特征值和的比例κ超过99.99%。由式(8) 和式(14)可知，对于目标函数ζ，活跃子空间维数为一，选取一维对应的特征向量作为活跃子空间的基向量。此时式(11)中基向量，U=[0.256，-0.087 4，-0.010 1，0.015 9，0.025 5，-0.082 6，-0.033 6，0.026 3，-0.019 0，-0.038 2，0.088 8，0.103，0.94 7]。

图9 梯度协方差矩阵的特征值衰减曲线Fig.9 Curve of eigenvalues of gradient covariance matrix

表7 梯度协方差矩阵的特征值Table 7 Eigenvalues of gradient covariance matrix

为进一步验证活跃子空间选取的基向量的准确性，可以使用统计学中Bootstrap自助法重复抽取样本组成新的梯度协方差矩阵进行活跃子空间的基向量计算。取Bootstrap自助法重复次数为1 000，图10和图11分别为活跃子空间基向量各分量分布柱状图和区间，由图可知活跃子空间基向量的各分量分布较为集中，可以认为求得的U是准确的。