分布參數(shù)不確定性重要性分析的高效求積公式法

劉俊超,李璐祎

(西北工業(yè)大學 航空學院,陜西 西安 710072)

在工程實際中,由于生產(chǎn)水平、加工工藝等因素的制約,結(jié)構(gòu)系統(tǒng)的材料特性、結(jié)構(gòu)尺寸和制造誤差等變量中廣泛地存在多種不確定性[1]。一般按照不確定性的來源可將其分為2類:客觀的不確定性和主觀的不確定性?？陀^的不確定性與變量內(nèi)在的固有隨機性有關(guān),是不可縮減的固有屬性,而主觀的不確定性與不完善的信息及有限的數(shù)據(jù)有關(guān),可隨著認知的深入和數(shù)據(jù)的豐富而減小[2]。

靈敏度分析常常用來衡量模型輸入不確定性對模型輸出不確定性的貢獻程度[3]。靈敏度分析分為局部靈敏度分析和全局靈敏度分析。其中,全局靈敏度分析又被稱為重要性測度分析,用來衡量變量在其整個不確定性范圍內(nèi)變化時對輸出性能統(tǒng)計特征的貢獻程度[4]。目前國內(nèi)外學者們對輸入變量的重要性測度分析方法做了大量的研究,比如Saltelli[5]提出了基于非參數(shù)的方法,Sobol[6]提出了基于方差的重要性分析方法,Borgonovo[7]提出了矩獨立分析方法。由于基于方差的重要性測度分析方法具有模型獨立、重要性測度與輸入-輸出關(guān)系一一對應(yīng)等良好的特性,是目前國際上最常使用的方法之一[8]。

上述提到的重要性測度分析方法并沒有考慮輸入變量分布參數(shù)的不確定性,但實際上由于認知和測量水平的限制,輸入變量的分布參數(shù)并非是確切可知的,因此分布參數(shù)具有不確定性的重要性分析方法得到了越來越多的關(guān)注[9]。分布參數(shù)的不確定性有多種描述方式,常見的有區(qū)間模型[10]、模糊集模型[11]以及主觀概率密度函數(shù)[12]等。對于存在分布參數(shù)不確定性的情況,輸入變量和分布參數(shù)的不確定性傳遞過程可定性地描述為:分布參數(shù)→輸入變量→輸出變量→輸出性能統(tǒng)計特征。因此分布參數(shù)的不確定性會引起輸出性能統(tǒng)計特征的不確定性,這就要求在進行重要性分析時考慮分布參數(shù)在其整個變化范圍內(nèi)對模型輸出性能統(tǒng)計特征的影響。

對于分布參數(shù)具有不確定性的模型,本文借助輸入變量基于方差的重要性分析思想,建立了衡量分布參數(shù)對輸出統(tǒng)計矩影響的重要性測度指標,并針對直接求解分布參數(shù)重要性測度指標需要3層蒙特卡洛(Monte Carlo,MC)抽樣,計算成本過高的問題,提出2種基于代理抽樣概率密度函數(shù)(surrogate sampling probability density function,SSPDF)[13]的求積公式(cubature formula,CF)方法來解決此問題。所提分析方法首先將求積公式運用到分布參數(shù)不確定性基于方差的重要性測度指標求解過程中,以提高分布參數(shù)重要性測度指標中嵌套的期望和方差算子的計算效率;其次,引入代理抽樣概率密度函數(shù)方法,以降低輸入變量參數(shù)不確定性向輸出性能統(tǒng)計矩傳遞過程中的計算量;最后,將單層蒙特卡洛(quasi-Monte Carlo,QMC)方法[5]與求積公式相結(jié)合,以降低重要性測度指標求解過程中嵌套的期望和方差算子的計算量。

1 分布參數(shù)基于方差的重要性測度指標

當輸入變量X的分布參數(shù)具有不確定性時,輸出的統(tǒng)計特征值也將具有不確定性。衡量輸入變量的單個分布參數(shù)或多個分布參數(shù)的交互作用對模型輸出統(tǒng)計特征值(這里以均值、方差為例)方差的貢獻程度,可以為改變分布參數(shù)以減小模型輸出的不確定度提供依據(jù)[14]。

當輸入變量X=(X1,X2,…,XdX)的分布參數(shù)Θ具有不確定性時,模型輸出將同時具有主觀不確定和客觀不確定,可表示為

Y=g(X,Θ)

(1)

式中,Θ=(Θ1,Θ2,…,ΘdΘ)是輸入變量X的dΘ維相互獨立的分布參數(shù)變量。在本文中,參數(shù)的主觀不確定性用主觀概率密度函數(shù)fΘ(θ)來描述[12]。當Θ取某一特定值θ(·)時,輸入變量X的不確定性由條件概率密度函數(shù)fX(x|θ(·))來衡量。

當輸入變量X的分布參數(shù)Θ具有不確定性時,輸出的均值和方差不再只受到輸入變量X的不確定性影響,而是可以描述為分布參數(shù)Θ的函數(shù)

因此,類似于輸入變量的基于方差的重要性測度的定義,分布參數(shù)Θ分別對Μ和D的主重要性測度和總重要性測度如(4)～(5)式和(6)～(7)式所示

2 分布參數(shù)基于方差的重要性測度指標求解的高效方法

從(4)～(7)式可以看出,求解參數(shù)的重要性測度指標的關(guān)鍵仍在于求解嵌套的期望和方差算子。然而,不同于輸入變量基于方差的重要性測度指標,參數(shù)的重要性測度指標求解過程需要將參數(shù)的不確定性傳遞到輸出均值和方差,是一個分布參數(shù)Θ和輸入變量X嵌套的雙層抽樣過程。如果利用抽樣方法求解(4)～(7)式中指標,則整個過程是一個“3層嵌套”的抽樣過程,計算量太大難以為工程實際所接受[15]。求積公式利用少量的積分點和權(quán)重便能計算輸出的均值和方差,并且對更高維度的輸入仍有良好的適用性與精確度[16]。Li等[13]提出的代理抽樣概率密度函數(shù)方法可以將參數(shù)不確定性向輸出均值和方差傳遞過程中的雙層抽樣簡化為單層,在很大程度上提高了參數(shù)不確定性傳遞的效率。因此,本文將CF和SSPDF引入到參數(shù)的重要性分析中,首先提出了參數(shù)重要性測度指標求解的S-DLCF算法,然后再結(jié)合QMC方法,提出了S-SLCF算法。

2.1 求積公式

對于一個模型Z=φ(Θ)(Z可以是(2)式中的均值M或(3)式中的方差D),由于CF估計統(tǒng)計矩是在標準正態(tài)空間中進行的[17],因此首先通過Rosenblatt變換或Nataf變換將該模型轉(zhuǎn)化為獨立標準正態(tài)變量λ=(λ1,λ2,…,λdΘ)的一個函數(shù)[18]

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

(8)

式中，R-1(·)表示Rosenblatt變換或Nataf變換；ρ(λ)表示獨立標準正態(tài)變量λ的多元函數(shù)。

因此,根據(jù)(8)式,Z的一階和二階中心距,也即Z的期望和方差,可以表示為

式中,fλ(λ)表示獨立標準正態(tài)變量λ的聯(lián)合概率密度函數(shù)。

根據(jù)求積公式的理論,(9)～(10)式中的積分可利用少量合適的積分點和相應(yīng)的積分權(quán)重高效地計算得出,可以表示為[19]

式中,N是求積權(quán)重ωj與相應(yīng)的積分點λj的數(shù)目。

Xu等[18]已經(jīng)證明了CF能夠精確估計大多數(shù)問題中函數(shù)的一階矩和二階矩。表1～2中列出了目前已知的最有效的5類求積公式,從其中可以看出,CF所需的積分點個數(shù)與變量維數(shù)有關(guān)且呈二次增長關(guān)系。

表1 公式Ⅰ～Ⅳ的積分點個數(shù)[18]

表2 公式Ⅴ中的積分點個數(shù)[18]

對于中低等維度變量,積分點的個數(shù)只有幾十個或幾百個,體現(xiàn)出了高效性[19]。公式Ⅰ、Ⅱ和Ⅴ受限于變量維度,公式Ⅳ的計算誤差隨變量維數(shù)的增加而增加,公式Ⅲ的計算誤差相比較而言更加穩(wěn)定[16]。

因此,本文選擇公式Ⅲ計算嵌套的期望和方差,其表示形式為

(13)

2.2 代理抽樣概率密度函數(shù)方法

2.1節(jié)中的求積公式方法可以高效地求解(4)～(7)式中嵌套的期望和方差算子。然而,在求解過程中首先需要將參數(shù)Θ的不確定性傳遞到輸出的均值和方差,該過程需要雙層嵌套的輸入變量X和參數(shù)Θ抽樣。為解決計算成本昂貴的問題，Li等[13]提出了代理抽樣概率密度函數(shù)方法。

考慮(8)式中的輸出均值模型,通過引入SSPDFhX(x|θ*),可將輸出均值M表示為

(14)

根據(jù)方差的計算式V(Y)=E(Y2)-E2(Y),(9)式中的輸出方差模型引入SSPDFhX(x|θ*)后,可將輸出方差D表示為

(15)

由(14)～(15)式可知,當在輸出均值和方差函數(shù)的計算過程中引入SSPDFhX(x|θ*)時,輸入變量X是由含有確定參數(shù)θ*的SSPDFhX(x|θ*)產(chǎn)生,不依賴于真實的分布參數(shù)。因此,在計算M和D時,分布參數(shù)Θ在外層變化,內(nèi)層所產(chǎn)生的輸入變量X的樣本可重復使用。

采用(14)～(15)式計算M和D時,SSPDFhX(x|θ*)的選取是關(guān)鍵。文獻[13]指出,當輸入變量X隨其具有不確定性的分布參數(shù)Θ改變時,hX(x|θ*)應(yīng)該覆蓋輸入變量X的整個變化范圍。一種簡單而直接的方法是,根據(jù)分布參數(shù)Θ的變化范圍確定相應(yīng)的輸入變量X的極限分布,然后根據(jù)輸入變量X的極限分布確定SSPDFhX(x|θ*)。具體的選取方法可參考文獻[13]。

2.3 基于代理抽樣概率密度函數(shù)的雙層求積公式(S-DLCF)算法

將(2)～(3)式中的輸出期望和方差模型統(tǒng)一表示為2.1節(jié)提到的Z=φ(Θ),參數(shù)Θ的不確定性可以通過2.2節(jié)的代理抽樣概率密度方法傳遞到輸出Z。從(4)～(7)式可以看出,求解重要性測度指標的關(guān)鍵是計算分布參數(shù)對Z方差的貢獻VΘi[EΘ～i(Z|Θi)]和VΘ～i[EΘi(Z|Θ～i)],而這些方差貢獻均是期望和方差算子的嵌套。由于該算法的雙層嵌套過程較長,將其分為兩部分進行描述。

首先,求解SiZ,也即通過S-DLCF求解分布參數(shù)Θi對Z的主重要性測度,其基本步驟如下所示[19]

1) 估計無條件期望E(Z)和方差V(Z)

②根據(jù)2.2節(jié)方法確定SSPDFhX(x|θ*),并由hX(x|θ*)抽取輸入變量X的NX個樣本xk(k=1,…,NX),進而求得相應(yīng)的模型輸出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 基于代理抽樣概率密度函數(shù)的單層求積公式(S-SLCF)算法

該算法將QMC方法[5]的思想引入到求解嵌套的期望與方差中。對于模型Z=φ(Θ),以方差貢獻VΘi[EΘ～i(Z|Θi)]和VΘ～i[EΘi(Z|Θ～i)]為例,準蒙特卡洛法求解嵌套期望和方差可表示為

(16)

根據(jù)(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方法求解分布參數(shù)Θi對Z的主重要性測度和總重要性測度的基本步驟為：

3) 根據(jù)2.2節(jié)確定相應(yīng)的SSPDFhX(x|θ*),用hX(x|θ*)抽取輸入變量X的NX個樣本xk(k=1,…,NX),進而求得相應(yīng)的模型輸出Y。

8) 將E[G(Θ)]和E(Z)代入到式(16)中得到VΘi[EΘ～i(Z|Θi)],將該結(jié)果及步驟5)中求得的V(Z)代入到(4)式或(5)式中得到SiZ。

2.5 2種算法的比較

3 算例分析

在本節(jié)的算例計算過程中,為了驗證所提方法的效率和精度,以MC方法所得結(jié)果作為參考解[20],以QMC方法所得結(jié)果作為對比解。在求解分布參數(shù)的主重要性測度指標和總重要性測度指標時,MC方法的模型調(diào)用次數(shù)為dΘ×NΘ×NΘ×NX,QMC方法的模型調(diào)用次數(shù)為(dΘ+2)×NΘ×NX,其中MC方法中的樣本量統(tǒng)一為NX=NΘ=2 500,QMC方法中的樣本量統(tǒng)一為NX=NΘ=30 000,以保證所求的重要性測度指標收斂。

3.1 可加性非線性模型

根據(jù)SSPDF的選取原則,該算例中選取SSPDF為正態(tài)分布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的重要性測度指標的計算結(jié)果

3.2 鉚接結(jié)構(gòu)

在航空工業(yè)中,真實的鉚接過程非常復雜,本文以無頭鉚釘為例,將鉚接過程簡化為如圖1所示的2個階段。第Ⅰ階段中,鉚釘從狀態(tài)A(沖擊前的初始狀態(tài),無變形)到狀態(tài)B(中間狀態(tài),鉚釘和孔隙之間無縫隙)。第Ⅱ階段中,鉚釘從狀態(tài)B到狀態(tài)C(鉚釘受沖擊后的最終狀態(tài),頭部被加工成型)。

圖1 簡化的無頭鉚釘鉚接過程

為了建立鉚接過程中的鉚釘擠壓應(yīng)力與幾何尺寸之間的數(shù)學關(guān)系,可以假設(shè)以下幾個理想條件:

1) 鉚接過程中鉚釘孔直徑不變

2) 鉚釘體積的改變忽略不計

3) 鉚接過程結(jié)束后,鉚釘?shù)氖芰γ鏋閳A柱面

4) 鉚釘材料具有各向同性

鉚接前,圖1a)中的狀態(tài)A下鉚釘?shù)某跏俭w積V0可以由(18)式表示

(18)

式中,d和h分別表示狀態(tài)A下鉚釘?shù)闹睆胶烷L度。

第Ⅰ階段后,狀態(tài)B下鉚釘?shù)捏w積可由(19)式表示為

(19)

式中,D0和h1分別表示狀態(tài)B下鉚釘?shù)闹睆胶烷L度。

第Ⅱ階段后,假設(shè)鉚釘在狀態(tài)C下的上下表面的尺寸相同,則狀態(tài)C下鉚釘?shù)捏w積可由(20)式表示為

(20)

式中：t為薄壁件的整體厚度；D1和H分別表示狀態(tài)C下鉚釘頭的直徑和長度。

根據(jù)硬化強度理論,y方向上的最大擠壓應(yīng)力σmax可以由(21)式表示為

σmax=K(εy)nSHE

(21)

式中：K為強度系數(shù)；nSHE為鉚釘材料的硬化因子；εy為鉚釘頭在鉚接過程中y方向上的真實應(yīng)變。真實應(yīng)變εy由第Ⅰ階段中產(chǎn)生的應(yīng)變εy1和第Ⅱ階段中產(chǎn)生的應(yīng)變εy2組成,因此,εy可由(22)式表示為

εy=εy1+εy2

(22)

已經(jīng)假設(shè)鉚釘在鉚接過程中體積不變,由(18)～(22)式可得鉚釘?shù)淖畲髷D壓應(yīng)力可以由(23)式表示為

(23)

本文選取的鉚釘材料為2017-T4,其硬化指數(shù)nSHE=0.15。為保證鉚釘上下端有一定的余量,假設(shè)鉚釘頭高度H=2.2 mm。根據(jù)材料手冊,鉚釘?shù)臄D壓強度為σsq=565 MPa,如果最大擠壓應(yīng)力大于擠壓強度,鉚釘就可能會失效,因此可以建立極限狀態(tài)方程如(24)式所示

g=σsq-σmax

(24)

在整個鉚接過程中,假設(shè)各隨機變量之間相互獨立且服從表4所示的正態(tài)分布,其均值存在主觀不確定性且服從表5所示的正態(tài)分布。

表4 無頭鉚釘?shù)妮斎胱兞康姆植紖?shù)

表5 無頭鉚釘?shù)妮斎胱兞烤档姆植紖?shù)

根據(jù)選取原則,該算例中輸入變量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種方法得到的重要性測度指標以及每種方法所需調(diào)用功能函數(shù)的次數(shù),從中可以得出,4種方法的計算量分別為5×(2 500)3,7×(30 000)2,181×10 000,18 000。從表中可以看出，相比于MC方法和QMC方法,S-DLCF方法和S-SLCF方法在滿足計算精度的前提下,2種新方法都能有效降低分布參數(shù)重要性分析過程中的計算量,并且S-SLCF方法比S-DLCF方法的計算量更少,表現(xiàn)了S-SLCF方法在重要性測度指標求解過程中更高的計算效率。

表6 算例3.2的重要性測度指標的計算結(jié)果

圖2 算例3.2在不同方法下計算得出的輸出均值重要性測度對比圖

圖3 算例3.2在不同方法下計算得出的輸出方差重要性測度對比圖

4 結(jié) 論

本文研究了分布參數(shù)具有不確定性時的重要性分析問題,首先定義了衡量分布參數(shù)不確定性對輸出均值和方差影響的重要性測度指標。其次,針對傳統(tǒng)MC方法在進行參數(shù)重要性分析時計算量大、工程難以接受的缺點,以求積公式為基礎(chǔ),結(jié)合SSPDF,提出了參數(shù)不確定性重要性測度求解的兩種新算法,大大降低了參數(shù)不確定性重要性分析的成本。

所提新算法的效率和精度取決于兩方面:一方面是CF的選擇,另一方面是SSPDF的選取。對于CF的選擇,依賴于問題的維數(shù),本文選擇的是求積公式Ⅲ,該公式在求解低等到中等維度的問題時有較好的精度及效率;對于SSPDF的選取,一個合適的SSPDF應(yīng)該覆蓋分布參數(shù)確定的輸入變量的整個變化范圍,本文采用了與原輸入變量具有相同分布的方法確定SSPDF。

通過數(shù)值算例和工程算例的分析結(jié)果可以得出,所提的新算法在求解主重要性測度和總重要性測度時都體現(xiàn)出了較高的精度,并能準確對分布參數(shù)的重要性進行排序。與QMC方法比較,2種新算法在大幅度減少計算量的同時表現(xiàn)出了更優(yōu)的計算精度;而S-DLCF與S-SLCF相比,S-DLCF所得結(jié)果的相對誤差更小,計算結(jié)果更有效;而S-SLCF通過引入準蒙特卡洛法,真正意義上將求解分布參數(shù)不確定性下的重要性測度指標時的3層MC抽樣方法簡化為單層抽樣,計算過程更加簡潔。在工程實際中,可以根據(jù)具體的問題選擇不同的算法。