劉俊超,李璐祎

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

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

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

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

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

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

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

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

Y=g(X,Θ)

(1)

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

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

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

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

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

2.1 求積公式

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

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

(8)

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

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

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

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

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

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

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

表2 公式Ⅴ中的積分點(diǎn)個(gè)數(shù)[18]

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

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

(13)

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

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

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

(14)

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

(15)

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

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

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

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

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

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

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

2) 估計(jì)EΘ～i(Z|Θi)和VΘi[EΘ～i(Z|Θi)]

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

1) 和求解主重要性測(cè)度時(shí)的步驟1)相同。

2) 估計(jì)VΘi(Z|Θ～i)和EΘ～i[VΘi(Z|Θ～i)]

Θ～i)]。

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

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

(16)

根據(jù)(16)式,通過將積分維度從dΘ維擴(kuò)展到2dΘ-1維,可將(10)式或(11)式中的雙層嵌套積分簡(jiǎn)化為單層積分,只需要分別估計(jì)期望E[G(Θ)]和無條件期望E(Z)便可得到VΘi[EΘ～i(Z|Θi)],進(jìn)而可得主重要性測(cè)度,這大大簡(jiǎn)化了求解過程。

同樣地,求解VΘ～i[EΘi(Z|Θ～i)]可被簡(jiǎn)化為

(17)

利用S-SLCF方法求解分布參數(shù)Θi對(duì)Z的主重要性測(cè)度和總重要性測(cè)度的基本步驟為：

3) 根據(jù)2.2節(jié)確定相應(yīng)的SSPDFhX(x|θ*),用hX(x|θ*)抽取輸入變量X的NX個(gè)樣本xk(k=1,…,NX),進(jìn)而求得相應(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é)的算例計(jì)算過程中,為了驗(yàn)證所提方法的效率和精度,以MC方法所得結(jié)果作為參考解[20],以QMC方法所得結(jié)果作為對(duì)比解。在求解分布參數(shù)的主重要性測(cè)度指標(biāo)和總重要性測(cè)度指標(biāo)時(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,以保證所求的重要性測(cè)度指標(biāo)收斂。

3.1 可加性非線性模型

根據(jù)SSPDF的選取原則,該算例中選取SSPDF為正態(tài)分布N(4,22)來抽取輸入變量X的樣本。

2種所提方法計(jì)算得到的重要性測(cè)度指標(biāo)如表3所示,eQMC,eS-DLCF,eS-SLCF分別表示各方法相對(duì)于MC方法得到的參考解的相對(duì)誤差。其中，4種方法的計(jì)算量分別為3×(2 500)3,5×(30 000)2,37×1 000,5 000。

表3 算例3.1的重要性測(cè)度指標(biāo)的計(jì)算結(jié)果

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

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

圖1 簡(jiǎn)化的無頭鉚釘鉚接過程

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

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

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

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

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

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

(18)

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

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

(19)

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

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

(20)

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

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

σmax=K(εy)nSHE

(21)

式中：K為強(qiáng)度系數(shù)；nSHE為鉚釘材料的硬化因子；εy為鉚釘頭在鉚接過程中y方向上的真實(shí)應(yīng)變。真實(shí)應(yīng)變?chǔ)舮由第Ⅰ階段中產(chǎn)生的應(yīng)變?chǔ)舮1和第Ⅱ階段中產(chǎn)生的應(yīng)變?chǔ)舮2組成,因此,ε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ù)材料手冊(cè),鉚釘?shù)臄D壓強(qiáng)度為σsq=565 MPa,如果最大擠壓應(yīng)力大于擠壓強(qiáng)度,鉚釘就可能會(huì)失效,因此可以建立極限狀態(tài)方程如(24)式所示

g=σsq-σmax

(24)

在整個(gè)鉚接過程中,假設(shè)各隨機(jī)變量之間相互獨(dú)立且服從表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種方法得到的重要性測(cè)度指標(biāo)以及每種方法所需調(diào)用功能函數(shù)的次數(shù),從中可以得出,4種方法的計(jì)算量分別為5×(2 500)3,7×(30 000)2,181×10 000,18 000。從表中可以看出，相比于MC方法和QMC方法,S-DLCF方法和S-SLCF方法在滿足計(jì)算精度的前提下,2種新方法都能有效降低分布參數(shù)重要性分析過程中的計(jì)算量,并且S-SLCF方法比S-DLCF方法的計(jì)算量更少,表現(xiàn)了S-SLCF方法在重要性測(cè)度指標(biāo)求解過程中更高的計(jì)算效率。

表6 算例3.2的重要性測(cè)度指標(biāo)的計(jì)算結(jié)果

圖2 算例3.2在不同方法下計(jì)算得出的輸出均值重要性測(cè)度對(duì)比圖

圖3 算例3.2在不同方法下計(jì)算得出的輸出方差重要性測(cè)度對(duì)比圖

4 結(jié) 論

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

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

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