分布参数不确定性重要性分析的高效求积公式法

刘俊超,李璐祎

(西北工业大学 航空学院,陕西 西安 710072)

在工程实际中,由于生产水平、加工工艺等因素的制约,结构系统的材料特性、结构尺寸和制造误差等变量中广泛地存在多种不确定性[1]。一般按照不确定性的来源可将其分为2类:客观的不确定性和主观的不确定性。客观的不确定性与变量内在的固有随机性有关,是不可缩减的固有属性,而主观的不确定性与不完善的信息及有限的数据有关,可随着认知的深入和数据的丰富而减小[2]。

灵敏度分析常常用来衡量模型输入不确定性对模型输出不确定性的贡献程度[3]。灵敏度分析分为局部灵敏度分析和全局灵敏度分析。其中,全局灵敏度分析又被称为重要性测度分析,用来衡量变量在其整个不确定性范围内变化时对输出性能统计特征的贡献程度[4]。目前国内外学者们对输入变量的重要性测度分析方法做了大量的研究,比如Saltelli[5]提出了基于非参数的方法,Sobol[6]提出了基于方差的重要性分析方法,Borgonovo[7]提出了矩独立分析方法。由于基于方差的重要性测度分析方法具有模型独立、重要性测度与输入-输出关系一一对应等良好的特性,是目前国际上最常使用的方法之一[8]。

上述提到的重要性测度分析方法并没有考虑输入变量分布参数的不确定性,但实际上由于认知和测量水平的限制,输入变量的分布参数并非是确切可知的,因此分布参数具有不确定性的重要性分析方法得到了越来越多的关注[9]。分布参数的不确定性有多种描述方式,常见的有区间模型[10]、模糊集模型[11]以及主观概率密度函数[12]等。对于存在分布参数不确定性的情况,输入变量和分布参数的不确定性传递过程可定性地描述为:分布参数→输入变量→输出变量→输出性能统计特征。因此分布参数的不确定性会引起输出性能统计特征的不确定性,这就要求在进行重要性分析时考虑分布参数在其整个变化范围内对模型输出性能统计特征的影响。

对于分布参数具有不确定性的模型,本文借助输入变量基于方差的重要性分析思想,建立了衡量分布参数对输出统计矩影响的重要性测度指标,并针对直接求解分布参数重要性测度指标需要3层蒙特卡洛(Monte Carlo,MC)抽样,计算成本过高的问题,提出2种基于代理抽样概率密度函数(surrogate sampling probability density function,SSPDF)[13]的求积公式(cubature formula,CF)方法来解决此问题。所提分析方法首先将求积公式运用到分布参数不确定性基于方差的重要性测度指标求解过程中,以提高分布参数重要性测度指标中嵌套的期望和方差算子的计算效率;其次,引入代理抽样概率密度函数方法,以降低输入变量参数不确定性向输出性能统计矩传递过程中的计算量;最后,将单层蒙特卡洛(quasi-Monte Carlo,QMC)方法[5]与求积公式相结合,以降低重要性测度指标求解过程中嵌套的期望和方差算子的计算量。

1 分布参数基于方差的重要性测度指标

当输入变量X的分布参数具有不确定性时,输出的统计特征值也将具有不确定性。衡量输入变量的单个分布参数或多个分布参数的交互作用对模型输出统计特征值(这里以均值、方差为例)方差的贡献程度,可以为改变分布参数以减小模型输出的不确定度提供依据[14]。

当输入变量X=(X1,X2,…,XdX)的分布参数Θ具有不确定性时,模型输出将同时具有主观不确定和客观不确定,可表示为

Y=g(X,Θ)

(1)

式中,Θ=(Θ1,Θ2,…,ΘdΘ)是输入变量X的dΘ维相互独立的分布参数变量。在本文中,参数的主观不确定性用主观概率密度函数fΘ(θ)来描述[12]。当Θ取某一特定值θ(·)时,输入变量X的不确定性由条件概率密度函数fX(x|θ(·))来衡量。

当输入变量X的分布参数Θ具有不确定性时,输出的均值和方差不再只受到输入变量X的不确定性影响,而是可以描述为分布参数Θ的函数

因此,类似于输入变量的基于方差的重要性测度的定义,分布参数Θ分别对Μ和D的主重要性测度和总重要性测度如(4)～(5)式和(6)～(7)式所示

2 分布参数基于方差的重要性测度指标求解的高效方法

从(4)～(7)式可以看出,求解参数的重要性测度指标的关键仍在于求解嵌套的期望和方差算子。然而,不同于输入变量基于方差的重要性测度指标,参数的重要性测度指标求解过程需要将参数的不确定性传递到输出均值和方差,是一个分布参数Θ和输入变量X嵌套的双层抽样过程。如果利用抽样方法求解(4)～(7)式中指标,则整个过程是一个“3层嵌套”的抽样过程,计算量太大难以为工程实际所接受[15]。求积公式利用少量的积分点和权重便能计算输出的均值和方差,并且对更高维度的输入仍有良好的适用性与精确度[16]。Li等[13]提出的代理抽样概率密度函数方法可以将参数不确定性向输出均值和方差传递过程中的双层抽样简化为单层,在很大程度上提高了参数不确定性传递的效率。因此,本文将CF和SSPDF引入到参数的重要性分析中,首先提出了参数重要性测度指标求解的S-DLCF算法,然后再结合QMC方法,提出了S-SLCF算法。

2.1 求积公式

对于一个模型Z=φ(Θ)(Z可以是(2)式中的均值M或(3)式中的方差D),由于CF估计统计矩是在标准正态空间中进行的[17],因此首先通过Rosenblatt变换或Nataf变换将该模型转化为独立标准正态变量λ=(λ1,λ2,…,λdΘ)的一个函数[18]

Z=φ(Θ)=φ(R-1(λ))=ρ(λ)

(8)

式中，R-1(·)表示Rosenblatt变换或Nataf变换；ρ(λ)表示独立标准正态变量λ的多元函数。

因此,根据(8)式,Z的一阶和二阶中心距,也即Z的期望和方差,可以表示为

式中,fλ(λ)表示独立标准正态变量λ的联合概率密度函数。

根据求积公式的理论,(9)～(10)式中的积分可利用少量合适的积分点和相应的积分权重高效地计算得出,可以表示为[19]

式中,N是求积权重ωj与相应的积分点λj的数目。

Xu等[18]已经证明了CF能够精确估计大多数问题中函数的一阶矩和二阶矩。表1～2中列出了目前已知的最有效的5类求积公式,从其中可以看出,CF所需的积分点个数与变量维数有关且呈二次增长关系。

表1 公式Ⅰ～Ⅳ的积分点个数[18]

表2 公式Ⅴ中的积分点个数[18]

对于中低等维度变量,积分点的个数只有几十个或几百个,体现出了高效性[19]。公式Ⅰ、Ⅱ和Ⅴ受限于变量维度,公式Ⅳ的计算误差随变量维数的增加而增加,公式Ⅲ的计算误差相比较而言更加稳定[16]。

因此,本文选择公式Ⅲ计算嵌套的期望和方差,其表示形式为

(13)

2.2 代理抽样概率密度函数方法

2.1节中的求积公式方法可以高效地求解(4)～(7)式中嵌套的期望和方差算子。然而,在求解过程中首先需要将参数Θ的不确定性传递到输出的均值和方差,该过程需要双层嵌套的输入变量X和参数Θ抽样。为解决计算成本昂贵的问题，Li等[13]提出了代理抽样概率密度函数方法。

考虑(8)式中的输出均值模型,通过引入SSPDFhX(x|θ*),可将输出均值M表示为

(14)

根据方差的计算式V(Y)=E(Y2)-E2(Y),(9)式中的输出方差模型引入SSPDFhX(x|θ*)后,可将输出方差D表示为

(15)

由(14)～(15)式可知,当在输出均值和方差函数的计算过程中引入SSPDFhX(x|θ*)时,输入变量X是由含有确定参数θ*的SSPDFhX(x|θ*)产生,不依赖于真实的分布参数。因此,在计算M和D时,分布参数Θ在外层变化,内层所产生的输入变量X的样本可重复使用。

采用(14)～(15)式计算M和D时,SSPDFhX(x|θ*)的选取是关键。文献[13]指出,当输入变量X随其具有不确定性的分布参数Θ改变时,hX(x|θ*)应该覆盖输入变量X的整个变化范围。一种简单而直接的方法是,根据分布参数Θ的变化范围确定相应的输入变量X的极限分布,然后根据输入变量X的极限分布确定SSPDFhX(x|θ*)。具体的选取方法可参考文献[13]。

2.3 基于代理抽样概率密度函数的双层求积公式(S-DLCF)算法

将(2)～(3)式中的输出期望和方差模型统一表示为2.1节提到的Z=φ(Θ),参数Θ的不确定性可以通过2.2节的代理抽样概率密度方法传递到输出Z。从(4)～(7)式可以看出,求解重要性测度指标的关键是计算分布参数对Z方差的贡献VΘi[EΘ～i(Z|Θi)]和VΘ～i[EΘi(Z|Θ～i)],而这些方差贡献均是期望和方差算子的嵌套。由于该算法的双层嵌套过程较长,将其分为两部分进行描述。

首先,求解SiZ,也即通过S-DLCF求解分布参数Θi对Z的主重要性测度,其基本步骤如下所示[19]

1) 估计无条件期望E(Z)和方差V(Z)

②根据2.2节方法确定SSPDFhX(x|θ*),并由hX(x|θ*)抽取输入变量X的NX个样本xk(k=1,…,NX),进而求得相应的模型输出Y。

2) 估计EΘ～i(Z|Θi)和VΘi[EΘ～i(Z|Θi)]

3) 将求得的V(Z)和VΘi[EΘ～i(Z|Θi)]代入(4)或(5)式中便可得到SiZ。

1) 和求解主重要性测度时的步骤1)相同。

2) 估计VΘi(Z|Θ～i)和EΘ～i[VΘi(Z|Θ～i)]

Θ～i)]。

2.4 基于代理抽样概率密度函数的单层求积公式(S-SLCF)算法

该算法将QMC方法[5]的思想引入到求解嵌套的期望与方差中。对于模型Z=φ(Θ),以方差贡献VΘi[EΘ～i(Z|Θi)]和VΘ～i[EΘi(Z|Θ～i)]为例,准蒙特卡洛法求解嵌套期望和方差可表示为

(16)

根据(16)式,通过将积分维度从dΘ维扩展到2dΘ-1维,可将(10)式或(11)式中的双层嵌套积分简化为单层积分,只需要分别估计期望E[G(Θ)]和无条件期望E(Z)便可得到VΘi[EΘ～i(Z|Θi)],进而可得主重要性测度,这大大简化了求解过程。

同样地,求解VΘ～i[EΘi(Z|Θ～i)]可被简化为

(17)

利用S-SLCF方法求解分布参数Θi对Z的主重要性测度和总重要性测度的基本步骤为：

3) 根据2.2节确定相应的SSPDFhX(x|θ*),用hX(x|θ*)抽取输入变量X的NX个样本xk(k=1,…,NX),进而求得相应的模型输出Y。

8) 将E[G(Θ)]和E(Z)代入到式(16)中得到VΘi[EΘ～i(Z|Θi)],将该结果及步骤5)中求得的V(Z)代入到(4)式或(5)式中得到SiZ。

2.5 2种算法的比较

3 算例分析

在本节的算例计算过程中,为了验证所提方法的效率和精度,以MC方法所得结果作为参考解[20],以QMC方法所得结果作为对比解。在求解分布参数的主重要性测度指标和总重要性测度指标时,MC方法的模型调用次数为dΘ×NΘ×NΘ×NX,QMC方法的模型调用次数为(dΘ+2)×NΘ×NX,其中MC方法中的样本量统一为NX=NΘ=2 500,QMC方法中的样本量统一为NX=NΘ=30 000,以保证所求的重要性测度指标收敛。

3.1 可加性非线性模型

根据SSPDF的选取原则,该算例中选取SSPDF为正态分布N(4,22)来抽取输入变量X的样本。

2种所提方法计算得到的重要性测度指标如表3所示,eQMC,eS-DLCF,eS-SLCF分别表示各方法相对于MC方法得到的参考解的相对误差。其中，4种方法的计算量分别为3×(2 500)3,5×(30 000)2,37×1 000,5 000。

表3 算例3.1的重要性测度指标的计算结果

3.2 铆接结构

在航空工业中,真实的铆接过程非常复杂,本文以无头铆钉为例,将铆接过程简化为如图1所示的2个阶段。第Ⅰ阶段中,铆钉从状态A(冲击前的初始状态,无变形)到状态B(中间状态,铆钉和孔隙之间无缝隙)。第Ⅱ阶段中,铆钉从状态B到状态C(铆钉受冲击后的最终状态,头部被加工成型)。

图1 简化的无头铆钉铆接过程

为了建立铆接过程中的铆钉挤压应力与几何尺寸之间的数学关系,可以假设以下几个理想条件:

1) 铆接过程中铆钉孔直径不变

2) 铆钉体积的改变忽略不计

3) 铆接过程结束后,铆钉的受力面为圆柱面

4) 铆钉材料具有各向同性

铆接前,图1a)中的状态A下铆钉的初始体积V0可以由(18)式表示

(18)

式中,d和h分别表示状态A下铆钉的直径和长度。

第Ⅰ阶段后,状态B下铆钉的体积可由(19)式表示为

(19)

式中,D0和h1分别表示状态B下铆钉的直径和长度。

第Ⅱ阶段后,假设铆钉在状态C下的上下表面的尺寸相同,则状态C下铆钉的体积可由(20)式表示为

(20)

式中：t为薄壁件的整体厚度；D1和H分别表示状态C下铆钉头的直径和长度。

根据硬化强度理论,y方向上的最大挤压应力σmax可以由(21)式表示为

σmax=K(εy)nSHE

(21)

式中：K为强度系数；nSHE为铆钉材料的硬化因子；εy为铆钉头在铆接过程中y方向上的真实应变。真实应变εy由第Ⅰ阶段中产生的应变εy1和第Ⅱ阶段中产生的应变εy2组成,因此,εy可由(22)式表示为

εy=εy1+εy2

(22)

已经假设铆钉在铆接过程中体积不变,由(18)～(22)式可得铆钉的最大挤压应力可以由(23)式表示为

(23)

本文选取的铆钉材料为2017-T4,其硬化指数nSHE=0.15。为保证铆钉上下端有一定的余量,假设铆钉头高度H=2.2 mm。根据材料手册,铆钉的挤压强度为σsq=565 MPa,如果最大挤压应力大于挤压强度,铆钉就可能会失效,因此可以建立极限状态方程如(24)式所示

g=σsq-σmax

(24)

在整个铆接过程中,假设各随机变量之间相互独立且服从表4所示的正态分布,其均值存在主观不确定性且服从表5所示的正态分布。

表4 无头铆钉的输入变量的分布参数

表5 无头铆钉的输入变量均值的分布参数

根据选取原则,该算例中输入变量d,h,D0,t和K的SSPDF分别为N(5,0.62),N(20,2.52),N(5.1,0.1022),N(5,0.12)和N(547.2,16.4162)。

表6展示了4种方法得到的重要性测度指标以及每种方法所需调用功能函数的次数,从中可以得出,4种方法的计算量分别为5×(2 500)3,7×(30 000)2,181×10 000,18 000。从表中可以看出，相比于MC方法和QMC方法,S-DLCF方法和S-SLCF方法在满足计算精度的前提下,2种新方法都能有效降低分布参数重要性分析过程中的计算量,并且S-SLCF方法比S-DLCF方法的计算量更少,表现了S-SLCF方法在重要性测度指标求解过程中更高的计算效率。

表6 算例3.2的重要性测度指标的计算结果

图2 算例3.2在不同方法下计算得出的输出均值重要性测度对比图

图3 算例3.2在不同方法下计算得出的输出方差重要性测度对比图

4 结 论

本文研究了分布参数具有不确定性时的重要性分析问题,首先定义了衡量分布参数不确定性对输出均值和方差影响的重要性测度指标。其次,针对传统MC方法在进行参数重要性分析时计算量大、工程难以接受的缺点,以求积公式为基础,结合SSPDF,提出了参数不确定性重要性测度求解的两种新算法,大大降低了参数不确定性重要性分析的成本。

所提新算法的效率和精度取决于两方面:一方面是CF的选择,另一方面是SSPDF的选取。对于CF的选择,依赖于问题的维数,本文选择的是求积公式Ⅲ,该公式在求解低等到中等维度的问题时有较好的精度及效率;对于SSPDF的选取,一个合适的SSPDF应该覆盖分布参数确定的输入变量的整个变化范围,本文采用了与原输入变量具有相同分布的方法确定SSPDF。

通过数值算例和工程算例的分析结果可以得出,所提的新算法在求解主重要性测度和总重要性测度时都体现出了较高的精度,并能准确对分布参数的重要性进行排序。与QMC方法比较,2种新算法在大幅度减少计算量的同时表现出了更优的计算精度;而S-DLCF与S-SLCF相比,S-DLCF所得结果的相对误差更小,计算结果更有效;而S-SLCF通过引入准蒙特卡洛法,真正意义上将求解分布参数不确定性下的重要性测度指标时的3层MC抽样方法简化为单层抽样,计算过程更加简洁。在工程实际中,可以根据具体的问题选择不同的算法。