周蘇婷， 呂震宙， 凌春燕，王燕萍

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

全局靈敏度分析可以從輸入變量的整個分布范圍來衡量輸入變量的不確定性對工程設(shè)計中所感興趣的輸出性能統(tǒng)計特征的貢獻(xiàn)程度[1-2]，通常在可靠性分析中，更關(guān)注結(jié)構(gòu)是否失效以及輸入變量對于結(jié)構(gòu)失效概率的影響程度。為了衡量輸入變量對于結(jié)構(gòu)失效概率的影響并對輸入變量的重要性進(jìn)行排序，Cui等[3]提出了可靠性全局靈敏度指標(biāo)。可靠性全局靈敏度指標(biāo)定義為某一輸入變量固定時，結(jié)構(gòu)的無條件失效概率與條件失效概率的差異在輸入變量的整個分布范圍內(nèi)的平均值，該指標(biāo)反映了輸入變量固定時對失效概率的平均影響。與可靠性局部靈敏度指標(biāo)[4-5]不同，可靠性全局靈敏度指標(biāo)可以從輸入變量的整個不確定性范圍來衡量其對于失效概率的平均影響。

本文提出了利用貝葉斯公式轉(zhuǎn)換以及元重要抽樣算法[6-7]結(jié)合自適應(yīng)Kriging代理的可靠性全局靈敏度，并組織了求解可靠性全局靈敏度指標(biāo)的高效算法。基于貝葉斯公式可以將可靠性全局靈敏度指標(biāo)的原始定義等價表示為輸入變量的無條件概率密度函數(shù)(Probability Density Fanction, PDF)與條件概率密度函數(shù)之間的差異與失效概率的乘積形式，從而避免條件失效概率的求解[8-9]?；诖说葍r形式，只需要一組樣本，便可以得到各輸入變量的可靠性全局靈敏度指標(biāo)。利用貝葉斯算法進(jìn)行可靠性全局靈敏度分析的研究已有很多，Wang等[8]利用貝葉斯公式，將可靠性全局靈敏度指標(biāo)的定義式進(jìn)行了等價轉(zhuǎn)換，并采用蒙特卡羅法來計算可靠性全局靈敏度指標(biāo)。該方法適用范圍雖廣，但在計算工程實(shí)際中的小失效概率問題時計算量過于龐大。Wang等[9]則結(jié)合重要抽樣方法來計算可靠性全局靈敏度指標(biāo)的貝葉斯等價形式。該方法只需一組樣本便可以計算得到所有輸入變量的可靠性全局靈敏度指標(biāo)，且計算量獨(dú)立于輸入變量的維數(shù)，明顯地提高了計算的效率。但是重要抽樣方法[10-11]需要構(gòu)造重要抽樣密度函數(shù)，一般來說需要將重要抽樣密度函數(shù)的中心放在設(shè)計點(diǎn)處，這就使得該方法依賴于其他方法來尋找設(shè)計點(diǎn)，并且重要抽樣方法對于多設(shè)計點(diǎn)的情況和多失效域的情況并不適用。Yun等[12]基于可靠性全局靈敏度指標(biāo)的貝葉斯公式轉(zhuǎn)換后的等價形式，利用子集模擬方法[13-16]結(jié)合重要抽樣的方法進(jìn)行可靠性全局靈敏度的計算。該方法在利用子集模擬結(jié)合重要抽樣的方法計算得到無條件失效概率的同時，通過重復(fù)利用子集模擬重要抽樣方法抽取的樣本，并利用Metropolis-Hastings(M-H)[17]準(zhǔn)則來實(shí)現(xiàn)失效樣本的轉(zhuǎn)換，進(jìn)而估計各輸入變量的失效條件概率密度函數(shù)。該方法同樣只需要一組樣本便可計算得到可靠性全局靈敏度指標(biāo)，且對于小失效概率問題是有效可行的。但該方法抽樣過程中所結(jié)合的重要抽樣方法仍然依賴于設(shè)計點(diǎn)的選取，因而該方法依舊很難適用于多設(shè)計點(diǎn)和多失效域的情況。

對于實(shí)際工程中普遍存在大量復(fù)雜的功能函數(shù)，自適應(yīng)Kriging(Adaptive Kriging, AK)代理模型法是處理這類問題的高效算法[18]。可靠性分析中涉及到的自適應(yīng)Kriging算法有自適應(yīng)Kriging代理模型結(jié)合Monte Carlo模擬(Monte Carlo Simulation, MCS)法(AK-MCS)[19]、自適應(yīng)Kriging代理模型結(jié)合重要抽樣(Important Sampling, IS)法(AK-IS)[20]以及元重要抽樣算法(Meta-IS)[6-7]等。AK-MCS算法將自適應(yīng)Kriging代理模型與Monte Carlo模擬法相結(jié)合，其首先由輸入變量的聯(lián)合概率密度函數(shù)產(chǎn)生MCS備選樣本池，然后利用學(xué)習(xí)函數(shù)在備選樣本池中逐步挑選對失效面擬合貢獻(xiàn)較大的點(diǎn)來更新Kriging模型，最終能夠確保Kriging模型能在一定的置信水平下識別樣本池內(nèi)樣本的功能函數(shù)值的正負(fù)號，從而較準(zhǔn)確地求得失效概率。但對于工程上常見的小失效概率問題，AK-MCS算法的樣本池容量龐大，進(jìn)而使得代理過程十分耗時。將重要抽樣法與自適應(yīng)Kriging過程相結(jié)合形成的AK-IS算法可以大大降低備選樣本池的規(guī)模，從而提高自適應(yīng)學(xué)習(xí)的效率，減少自適應(yīng)學(xué)習(xí)過程所消耗的時間。但是對于多設(shè)計點(diǎn)及多失效域問題，基于一次二階矩法構(gòu)造重要抽樣密度函數(shù)的AK-IS算法將失效，因而對于多設(shè)計點(diǎn)及多失效域的問題需要另尋解決途徑。元模型重要抽樣算法(Meta-IS)通過元模型構(gòu)造重要抽樣密度函數(shù)，進(jìn)而來抽取重要抽樣樣本點(diǎn)。該方法在降低AK-MCS樣本池規(guī)模的同時避免了設(shè)計點(diǎn)的求解，適用于多設(shè)計點(diǎn)和多失效域問題的可靠性分析。然而，Meta-IS算法在抽取重要抽樣樣本后仍需要求解重要抽樣樣本點(diǎn)處真實(shí)的功能函數(shù)值來估計失效概率，所以該算法仍有較大的計算量。如果在元重要抽樣的基礎(chǔ)上嵌入自適應(yīng)Kriging模型，形成Meta-IS-AK算法，則有可能極大程度地提高失效概率求解的效率，本文將采用這種思路求解失效概率，并在此基礎(chǔ)上組織可靠性全局靈敏度的算法。

利用可靠性全局靈敏指標(biāo)的貝葉斯轉(zhuǎn)換形式，只要能求得無條件失效概率和失效域條件下輸入變量的概率密度函數(shù)即可直接求得可靠性全局靈敏度，而求解無條件失效概率的數(shù)字模擬法可以同時完成失效概率的計算以及失效域條件下輸入變量概率密度函數(shù)的求解，為此本文將Meta-IS與自適應(yīng)Kriging模型結(jié)合起來，以便完成可靠性全局靈敏度的高效求解。本文所組織的Meta-IS-AK算法主要分3步來執(zhí)行。第1步是由元重要抽樣的迭代策略得到逐步逼近最優(yōu)重要抽樣函數(shù)的樣本點(diǎn)，在迭代收斂后便可以得到第1步的Kriging模型以及相應(yīng)的重要抽樣樣本點(diǎn)。第2步是基于重要抽樣樣本構(gòu)建自適應(yīng)Kriging模型，該步的Kriging模型是在第1步得到的Kriging模型的基礎(chǔ)上逐步更新得到的，其目的是為了構(gòu)建能夠?qū)Φ?步中重要抽樣樣本點(diǎn)失效與否做出準(zhǔn)確預(yù)測的Kriging模型，從而高效地求解出無條件失效概率。然而Meta-IS-AK算法抽取的失效域內(nèi)的樣本點(diǎn)的密度函數(shù)為最優(yōu)重要抽樣密度函數(shù)，因而無法直接用于估計輸入變量失效域內(nèi)的條件概率密度函數(shù)，為此第3步則對失效域內(nèi)服從于重要抽樣密度函數(shù)的樣本點(diǎn)進(jìn)行轉(zhuǎn)換，該轉(zhuǎn)換將利用Metropolis-Hastings準(zhǔn)則由重要抽樣密度函數(shù)的失效樣本點(diǎn)來得到原概率密度函數(shù)在失效域內(nèi)的樣本點(diǎn)，進(jìn)而求得各個輸入變量在失效域中的條件概率密度函數(shù)。

本文主要由以下幾部分組成。第1節(jié)簡要地回顧了已有的利用貝葉斯公式轉(zhuǎn)換的可靠性全局靈敏度的形式，第2節(jié)詳細(xì)介紹了本文所提的計算可靠性全局靈敏度的高效算法。第3節(jié)分別給出了數(shù)值與工程算例，驗證了所提算法的高效性，第4節(jié)則得出結(jié)論。

1 基于貝葉斯公式的可靠性全局靈敏度指標(biāo)的轉(zhuǎn)換形式

(1)

利用貝葉斯公式，可將條件失效概率Pf|Xi進(jìn)行如下轉(zhuǎn)換[8-9,12,23-24]：

(2)

將式(2)代入到式(1)中，即可得到基于貝葉斯公式的可靠性全局靈敏度指標(biāo)表達(dá)式：

(3)

2 求解可靠性全局靈敏度指標(biāo)的高效算法

2.1 失效概率Pf的估計

(4)

步驟1以gK1(x)作為gK2(x)的初始模型，即令gK2(x)=gK1(x)。

步驟3計算g(xu)，將(xu,g(xu))添加到訓(xùn)練集T中，形成更新的訓(xùn)練集T={T∪(xu,g(xu))}。

步驟4由T更新gK2(x)，并返回步驟2。

(5)

(6)

(7)

2.2 失效域條件下輸入變量的概率密度函數(shù)估計

步驟2計算下列比值r,并由r確定下一個馬爾可夫鏈樣本

i=1,2,…,NF-1

(8)

i=1,2,…,NF-1

(9)

式中：u為[0,1]區(qū)間上服從均勻分布的隨機(jī)數(shù)。

2.3 Meta-IS-AK算法求解可靠性全局靈敏度指標(biāo)流程圖

由2.1和2.2節(jié)的內(nèi)容可給出Meta-IS-AK算法求解可靠性全局靈敏度的流程如圖1所示。

圖1 Meta-IS-AK算法求解的流程圖Fig.1 Flowchart of estimating Meta-IS-AK algorithm

3 算例分析

為了驗證所提算法在求解可靠性全局靈敏度的效率和精度，本節(jié)給出了3個算例，分別采用準(zhǔn)Monte Carlo算法(QMC)、AK-MCS算法、Meta-IS算法及Meta-IS-AK算法求解可靠性全局靈敏度指標(biāo)。其中QMC算法結(jié)果將作為參照解。

3.1 不連續(xù)多失效域算例

考慮一個多失效域的算例，其功能函數(shù)為

(10)

(11)

圖2 算例3.1失效邊界Fig.2 Failure boundary of Example 3.1

表1 算例3.1可靠性全局靈敏度指標(biāo)的計算結(jié)果

Table 1Results of reliability global sensitivity index for Example 3.1

算法X1/10-3X2/10-3NcallQMC14.5[0.30]14.5[0.46]1×105AK-MCS14.2[0.23]14.1[0.35]440Meta-IS14.2[0.41]14.0[0.65]42+5 000Meta-IS-AK14.2[0.43]13.9[0.65]42+254

注：數(shù)據(jù)右上標(biāo)為可靠性全局靈敏度指標(biāo)估計值的標(biāo)準(zhǔn)差,Ncall為功能函數(shù)的調(diào)用次數(shù)。

3.2 鋼架結(jié)構(gòu)算例

對于一鋼架結(jié)構(gòu)，有4個失效模式：

(12)

因為是串聯(lián)系統(tǒng)，所以系統(tǒng)的功能函數(shù)g為g=min{g1,g2,g3,g4}。在各個失效模式的功能函數(shù)中，Mi(i=1,2,3)與S是獨(dú)立的，均值和標(biāo)準(zhǔn)差分別為:μMi=2(i=1,2,3),μS=1,σMi=2(i=1,2,3),σS=0.25。各輸入變量所對應(yīng)的可靠性全局靈敏度指標(biāo)的計算結(jié)果如表2所示。

表2 算例3.2可靠性全局靈敏度指標(biāo)的計算結(jié)果

由表2中結(jié)果可以看出，Meta-IS-AK算法得到的計算結(jié)果與QMC算法計算所得結(jié)果基本一致，且模型的調(diào)用次數(shù)低于QMC算法與AK-MCS算法以及Meta-IS算法。這說明Meta-IS-AK算法能夠很好地適用于多個失效模式的情況。對于此算例，輸入變量的重要性排序為S>M3>M1>M2。此外，輸入變量S對失效概率的影響較其他輸入變量的影響大很多。因此在進(jìn)行可靠性設(shè)計時通過調(diào)整輸入變量S的不確定性可以更有效地滿足可靠性要求。此外，M2對失效概率影響的靈敏度指標(biāo)較其他的輸入變量小很多，因此可以忽略M2的不確定性對于失效概率的影響，可以通過將M2固定在均值處來簡化失效概率的求解模型。

3.3 簡化的翼盒結(jié)構(gòu)算例(多設(shè)計點(diǎn)算例)

圖3(a)是簡化翼盒模型的示意圖。這種翼盒結(jié)構(gòu)由64個桿和42個板組成。64個桿根據(jù)它們的方向分為3組，x、y和z方向上桿的長度分別是2L、L和 3L。所有桿的截面積均為A，所有板的厚度均為TH，E和P分別是所有板和桿的彈性模量和外部載荷。泊松比為0.3。假設(shè)輸入變量是獨(dú)立的正態(tài)變量，其分布參數(shù)如表3所示。

圖3 翼盒結(jié)構(gòu)Fig.3 Wing box structure

表3 算例3.3輸入變量的分布參數(shù)

Table 3Distribution parameters of input variables of Example 3.3

變量分布均值變異系數(shù)A/m2正態(tài)1×10-40.1L/m正態(tài)0.20.1E/Pa正態(tài)7.1×10100.12P/N正態(tài)1 5000.1TH/m正態(tài)2.5×10-30.15

由于翼盒結(jié)構(gòu)承擔(dān)了來自機(jī)翼的大部分垂直方向的載荷，其自由端的位移可能相對較大。大的位移可能會引起桿的變形并導(dǎo)致桿或板的破壞，這在飛行過程中是不允許的。因此計算翼盒結(jié)構(gòu)的位移是很有必要的，且位移的計算可以通過有限元分析完成。圖3(b)顯示了在ANSYS 14.0中構(gòu)建的翼盒結(jié)構(gòu)的有限元模型，其中輸入變量固定在它們的均值處，此外，翼盒的變形如圖3(c)所示。

本文中，該翼盒的隱式功能函數(shù)是根據(jù)翼盒的最大位移Dmax不超過臨界值構(gòu)造的(這里臨界值等于0.01)，其表達(dá)式為

g(X)=0.01-|Dmax|

(13)

式中：Dmax=D(L,A,E,P,TH)，是關(guān)于輸入變量的隱式函數(shù)，由有限元分析確定。算例3.3的可靠性全局靈敏度計算結(jié)果如圖4所示，計算結(jié)果的標(biāo)準(zhǔn)差如表4所示。

圖4 算例3.3可靠性全局靈敏度指標(biāo)的計算結(jié)果Fig.4 Results of reliability global sensitivity index for Example 3.3

表4 算例3.3可靠性全局靈敏度指標(biāo)結(jié)果

Table 4Reliability global sensitivity index forExample 3.3

變量SDi/10-5QMCAK-MCSMeta-ISMeta-IS-AKA/m20.270.200.170.18L/m0.260.110.140.11E/Pa5.363.875.045.05P/N0.350.420.330.32TH/m8.178.856.376.30Ncall5×10512033+5 00033+59

由表4中的結(jié)果可以看出，Meta-IS-AK算法得到的計算結(jié)果精度較高且模型的調(diào)用次數(shù)低于QMC算法與AK-MCS算法以及Meta-IS算法，這說明了Meta-IS-AK算法能夠很好地適用于隱式功能函數(shù)的情況。對于算例3.3，輸入變量的的不確定性對于失效概率影響的重要性排序為TH>E>P>L>A。從圖4的結(jié)果可以看出，輸入變量TH對失效概率的影響較其他輸入變量的影響大很多，因此在進(jìn)行可靠性設(shè)計時調(diào)整輸入變量TH的不確定性可以更有效地滿足可靠性要求。從圖4的結(jié)果還可以看出，L和A的不確定性對于失效概率的影響較其他的輸入變量小很多，因此忽略L和A的不確定性對于失效概率的影響可以簡化失效概率的求解模型。

4 結(jié) 論

1) 可靠性全局靈敏度指標(biāo)可以有效地衡量輸入變量對于結(jié)構(gòu)失效概率的影響。為了高效地計算可靠性全局靈敏度指標(biāo)，本文提出了可靠性全局靈敏度指標(biāo)計算的Meta-IS-AK算法。

2) Meta-IS-AK算法在計算可靠性全局靈敏度指標(biāo)時充分利用了可靠性全局靈敏度指標(biāo)貝葉斯算法的維度獨(dú)立性，也使得該指標(biāo)的求解轉(zhuǎn)化為輸出樣本的分類問題，這就使可靠性全局靈敏度指標(biāo)可以通過嵌入式代理模型結(jié)合數(shù)字模擬法來求解。相比于傳統(tǒng)的重要抽樣函數(shù)，元重要抽樣法構(gòu)造的準(zhǔn)最優(yōu)抽樣密度函數(shù)能夠很好地適用于多設(shè)計點(diǎn)和多失效域的情況。在此基礎(chǔ)上構(gòu)造重構(gòu)的功能函數(shù)的代理模型能夠準(zhǔn)確地對準(zhǔn)重要抽樣密度函數(shù)抽取的樣本進(jìn)行分類，進(jìn)一步地提高了失效概率的計算效率。M-H準(zhǔn)則無需額外調(diào)用功能函數(shù)便可以將準(zhǔn)重要抽樣密度函數(shù)的失效樣本轉(zhuǎn)化成原概率密度函數(shù)的失效樣本，大大提高了求解輸入變量的失效條件概率密度函數(shù)的效率。

3) 所提方法能夠很好地適用于多設(shè)計點(diǎn)和多失效域的情況，也適用于隱函數(shù)問題的可靠性靈敏度分析。所提方法的計算效率遠(yuǎn)遠(yuǎn)高于QMC算法，同時較AK-MCS算法以及Meta-IS算法也有提高。所給算例充分驗證了以上結(jié)論。