一種新的矩獨(dú)立重要性測(cè)度分析方法及高效算法

鞏祥瑞， 呂震宙， 孫天宇， 張雷雷， 封雷

(1. 西安現(xiàn)代控制技術(shù)研究所, 西安 710065； 2. 西北工業(yè)大學(xué) 航空學(xué)院, 西安 710072)

靈敏度分析方法主要研究結(jié)構(gòu)系統(tǒng)的輸入變量對(duì)輸出響應(yīng)的影響程度[1]，已廣泛應(yīng)用于可靠性工程[2]、危險(xiǎn)分析[3-4]、環(huán)境科學(xué)[5]、計(jì)算物理化學(xué)[6-7]等領(lǐng)域。靈敏度分析一般分為局部靈敏度分析方法和全局靈敏度分析方法。局部靈敏度分析方法研究模型輸入變量的分布參數(shù)對(duì)結(jié)構(gòu)系統(tǒng)輸出響應(yīng)的影響程度[8]，只能反映輸入變量在某一名義值處對(duì)輸出響應(yīng)的影響，具有很大的局限性[9]。

全局靈敏度分析，也稱為重要性測(cè)度分析，研究一個(gè)或多個(gè)輸入變量對(duì)模型輸出響應(yīng)的影響程度，能夠反映輸入變量對(duì)結(jié)構(gòu)系統(tǒng)輸出響應(yīng)的平均影響程度，因此，該分析方法在實(shí)際工程中得到了更廣泛的應(yīng)用。重要性測(cè)度分析方法的核心思想是對(duì)影響結(jié)構(gòu)系統(tǒng)輸出響應(yīng)的輸入變量進(jìn)行重要性排序，從而得到對(duì)結(jié)構(gòu)系統(tǒng)輸出響應(yīng)影響較大的輸入變量，為結(jié)構(gòu)系統(tǒng)的設(shè)計(jì)改進(jìn)提供理論依據(jù)。近年來，研究人員已經(jīng)提出了很多重要性測(cè)度分析方法。Saltelli 和 Helton提出了無參的方法，但缺乏模型的獨(dú)立性[9-10]。Sobol[11]和 Iman[12]等提出了基于方差的重要性測(cè)度分析方法和相應(yīng)的求解算法，該方法隱含地假設(shè)了方差可以充分描述結(jié)構(gòu)系統(tǒng)輸出響應(yīng)的不確定性。然而，相關(guān)學(xué)者指出隨機(jī)變量的任何一階矩只能反映其部分信息，必然會(huì)丟失一些信息[13]。為了分析輸入變量在其整個(gè)分布區(qū)域上對(duì)結(jié)構(gòu)系統(tǒng)輸出響應(yīng)的影響程度，Borgonovo[14]、Liu[15]和 Li[16]等又分別提出了矩獨(dú)立重要性測(cè)度分析(Moment-independent Importance Measures analysis，M-IM)方法。

目前，Borgonovo[14]所提出的矩獨(dú)立重要性測(cè)度分析方法應(yīng)用最為廣泛。但在可靠性分析中，結(jié)構(gòu)系統(tǒng)的失效概率通常是工程技術(shù)人員所關(guān)注的焦點(diǎn)。Li等[16]又提出了基于失效概率的矩獨(dú)立重要性測(cè)度分析方法，該方法能夠計(jì)算輸入變量在其分布區(qū)域的任一固定點(diǎn)時(shí)，結(jié)構(gòu)系統(tǒng)的無條件失效概率和條件失效概率的平均偏差。然而，在實(shí)際工程中將輸入變量固定于某一點(diǎn)，這是很難實(shí)現(xiàn)的，只能使輸入變量在其分布區(qū)域的某一區(qū)間上縮減變化。

為了解決這一問題，Wei等[17]提出了一種新的輸入變量縮減區(qū)間的獲取方法，該方法的基本思想是：首先產(chǎn)生2個(gè)隨機(jī)數(shù)q1和q2，q1和q2服從均勻分布，q1～U(0,1)，q2～U(q1,1)，然后根據(jù)輸入變量累積分布函數(shù)的逆函數(shù)F-1(·)，得到輸入變量的縮減區(qū)間[F-1(q1),F-1(q2)]。該方法的優(yōu)點(diǎn)是很容易獲取輸入變量的縮減區(qū)間，然而它并沒有考慮輸入變量縮減區(qū)間獲取的任意性和等可能性。

本文提出了一種新的基于失效概率的矩獨(dú)立重要性測(cè)度分析方法，并給出了一種更加合理的輸入變量縮減區(qū)間的獲取方法。同時(shí)，引入自適應(yīng)超球重要抽樣(Adaptive Radial-Based Importance Sampling，ARBIS)方法[18-19]，構(gòu)建所提新指標(biāo)的高效求解算法。應(yīng)用數(shù)值算例和工程算例，驗(yàn)證本文所提新指標(biāo)的合理性和求解算法的高效性。

1 新的矩獨(dú)立重要性測(cè)度指標(biāo)

1.1 新指標(biāo)的定義

EQ(E(IF)-E(IF|Xi∈Ui(Q)))2

(1)

式中：Pf為非條件失效概率；Q為分位數(shù)矩陣，用來生成輸入變量Xi的縮減區(qū)間；當(dāng)Xi的分布區(qū)間縮減到其所有可能的子區(qū)間Ui(Q)時(shí)，就能夠得到條件失效概率Pf|Xi∈Ui(Q)；EQ(·)為對(duì)Q求期望；IF為指示函數(shù)，當(dāng)X位于失效域(即g(X)≤0)時(shí)，IF(x)=1，否則IF(x)=0。

1.2 輸入變量縮減區(qū)間獲取方法

為了等概率地獲取輸入變量所有可能的縮減區(qū)間，首先定義一個(gè)分位數(shù)矩陣Hk(3×2)：

(2)

式中：hk(1)和hk(2)為[0,1]區(qū)間上的2個(gè)隨機(jī)數(shù)，滿足0≤hk(1)

Hk(k=1,2,…,N)組成一個(gè)(3N×2)的分位數(shù)矩陣Q：

(3)

式中：Hk(k=1,2,…,N) 是被等概率獲取的，例如P{Hk}=1/N。

為了更好地理解式(3)，可表示為

(4)

式中：q3k-2(1)=0，q3k-2(2)=q3k-1(1)=hk(1)，q3k-1(2)=q3k(1)=hk(2)，q3k(2)=1 (k=1,2,…,N)。

因此，根據(jù)分位數(shù)矩陣Q就可以定義縮減區(qū)間矩陣U為

(5)

對(duì)于輸入變量Xi，其縮減區(qū)間可表示為

(6)

上述是一種新的區(qū)間劃分技術(shù)來獲取輸入變量在其分布區(qū)域上所有可能的縮減區(qū)間。本文所提出的這種新的區(qū)間獲取方法具有2個(gè)優(yōu)點(diǎn)：①輸入變量所有可能的縮減區(qū)間能夠被等可能地獲取，從而就能夠更加合理地計(jì)算當(dāng)輸入變量在其分布區(qū)域的縮減區(qū)間上變化時(shí)對(duì)結(jié)構(gòu)系統(tǒng)失效概率的平均影響程度；②全期望公式在這種區(qū)間劃分技術(shù)中成立，基于全期望公式，第2節(jié)中就能夠?qū)⒈疚奶岢龅男碌幕谑Ц怕实木鬲?dú)立重要性測(cè)度指標(biāo)轉(zhuǎn)化為更加便于計(jì)算的基于方差的重要性測(cè)度指標(biāo)。

2 新指標(biāo)與基于方差的指標(biāo)的關(guān)系

由于1.2節(jié)中提出的輸入變量縮減區(qū)間獲取方法中全期望公式成立，就能夠推導(dǎo)出新的矩獨(dú)立重要性測(cè)度指標(biāo)與基于方差的重要性測(cè)度指標(biāo)間的關(guān)系，同時(shí)得到一種新的基于方差的矩獨(dú)立重要性測(cè)度指標(biāo)表示方式。

輸入變量在縮減區(qū)間上的全期望公式表示為

E(IF)=EQ(E(IF|Xi∈U(Q)))

(7)

證明過程詳見附錄A。

VQ(E(IF|Xi∈Ui(Q)))

(8)

從式(8)中可以看出，新的矩獨(dú)立重要性測(cè)度指標(biāo)能夠很容易地轉(zhuǎn)換為基于方差的重要性測(cè)度指標(biāo)，從而新的矩獨(dú)立重要性測(cè)度指標(biāo)就能夠采用基于方差的重要性測(cè)度分析方法來進(jìn)行求解。

3 新指標(biāo)的高效求解算法

新的矩獨(dú)立重要性測(cè)度指標(biāo)傳統(tǒng)的計(jì)算方法是雙層重復(fù)抽樣蒙特卡羅(Double-Loop-Repeat-Set Monte Carlo，DLRS MC)[17]方法，該方法先用N個(gè)樣本點(diǎn)計(jì)算非條件失效概率，然后在每一個(gè)子區(qū)間上，再產(chǎn)生N個(gè)新的樣本點(diǎn)來計(jì)算條件失效概率，最后計(jì)算非條件失效概率和條件失效概率差異的平均值。

DLRS MC方法可以獲得高精度的結(jié)果值，但其計(jì)算量很大，計(jì)算效率較低。本節(jié)提出了一種新的求解算法——ARBIS方法，旨在對(duì)新的矩獨(dú)立重要性測(cè)度指標(biāo)進(jìn)行高效求解。

3.1 自適應(yīng)超球重要抽樣方法

ARBIS方法最初是用來計(jì)算失效概率的，而計(jì)算新的基于失效概率的矩獨(dú)立重要性測(cè)度指標(biāo)的核心就是計(jì)算非條件失效概率和條件失效概率。因此，本節(jié)基于ARBIS方法構(gòu)建新指標(biāo)的求解算法。

ARBIS方法[18-19]的基本思想是自適應(yīng)地尋找到一個(gè)落在功能函數(shù)安全域內(nèi)的最大的超球。當(dāng)抽樣點(diǎn)落入超球內(nèi)部時(shí)，就將該抽樣點(diǎn)歸類為安全點(diǎn)，不必再代入結(jié)構(gòu)系統(tǒng)功能函數(shù)進(jìn)行求解。因此，該方法在計(jì)算結(jié)構(gòu)系統(tǒng)輸出響應(yīng)的失效概率時(shí)，計(jì)算量將大幅度降低。同時(shí)， ARBIS方法采用自適應(yīng)的方法來獲取滿足條件的最大超球半徑，相比于傳統(tǒng)的梯度搜索算法，效率更高，穩(wěn)健性也更好。

任何一個(gè)隨機(jī)變量通過合適的變換方法，都可以轉(zhuǎn)換到標(biāo)準(zhǔn)正態(tài)空間中[20-21]。因此，以下步驟都假定在標(biāo)準(zhǔn)正態(tài)空間中進(jìn)行。為了更好地理解ARBIS方法，給出該方法的一些關(guān)鍵步驟。

步驟1設(shè)定初始化超球半徑β0。 為了確保超球與失效域相交，初始半徑β0可通過式(9)來確定：

(9)

步驟2根據(jù)隨機(jī)變量的分布類型，通過轉(zhuǎn)化得到標(biāo)準(zhǔn)正態(tài)空間的樣本點(diǎn)uk(k=1,2,…,N)。

(10)

圖1 自適應(yīng)策略獲取最優(yōu)化半徑Fig.1 Adaptive strategy for obtaining optimal radius

步驟6重復(fù)步驟4和步驟5，直到收斂條件被滿足[22]，即可得到最優(yōu)的超球半徑βopt。

3.2 采用自適應(yīng)超球重要抽樣方法求解新指標(biāo)

當(dāng)輸入變量Xi縮減到區(qū)間Ui(Q)時(shí)，Xi就是一個(gè)截?cái)嘧兞?，則所有輸入變量的聯(lián)合概率密度函數(shù)可表示為

(11)

式中：K為截?cái)嘞禂?shù)，可表示為

(12)

對(duì)于式(12)的估計(jì)，并不需要計(jì)算結(jié)構(gòu)系統(tǒng)的功能函數(shù)值，其計(jì)算量可以忽略。因此，條件失效概率Pf|Xi∈Ui(Q)可重新表示為

(13)

從式(13)中可以看出，條件失效概率Pf|Xi∈Ui(Q)的計(jì)算已轉(zhuǎn)化為一個(gè)不含截?cái)嘧兞康亩嗄B(tài)并聯(lián)系統(tǒng)問題。經(jīng)過這種轉(zhuǎn)化，構(gòu)建ARBIS方法計(jì)算新的矩獨(dú)立重要性測(cè)度指標(biāo)的算法如下：

步驟1自適應(yīng)搜索獲取超球半徑βopt，并計(jì)算無條件失效概率Pf；同時(shí)得到半徑為βopt的超球外失效點(diǎn)組成的樣本點(diǎn)矩陣B。

步驟2當(dāng)輸入變量Xi縮減到區(qū)間Ui(Q)時(shí)，在樣本點(diǎn)矩陣B中搜索輸入變量Xi落在縮減區(qū)間Ui(Q)的樣本點(diǎn)，記樣本點(diǎn)個(gè)數(shù)為m。

步驟3計(jì)算條件失效概率：

(14)

式中：N為采用ARBIS方法的初始樣本點(diǎn)個(gè)數(shù)。

4 算例分析

4.1 數(shù)值算例

某一結(jié)構(gòu)系統(tǒng)的功能函數(shù)可表示為

(15)

式中：X1、X2和X3為相互獨(dú)立的輸入隨機(jī)變量，均服從標(biāo)準(zhǔn)正態(tài)分布，即Xi～N(0,1),i=1,2,3。

表1 DLRS MC方法和ARBIS方法計(jì)算的新的矩獨(dú)立重要性測(cè)度指標(biāo)值Table 1 New moment-independent importance measure indices of numerical example computed by DLRS MC and ARBIS methods

圖2 DLRS MC方法和ARBIS方法計(jì)算的新的矩獨(dú)立重要性測(cè)度指標(biāo)的收斂曲線Fig.2 Convergence curves of new moment-independent importance measure indices of numerical example computed by DLRS MC and ARBIS methods

4.2 工程算例

在汽車結(jié)構(gòu)中，車橋承載著大部分汽車重量，其通過懸臂與車架相接，將來自車輪的牽引力和制動(dòng)力，還有側(cè)向力經(jīng)過懸架傳遞給車架，起主要承載作用的就是汽車前軸[23]。由于工字型截?cái)嗔耗軌蛱岣呖箯潖?qiáng)度，因此前軸通常采用工字結(jié)構(gòu)梁。

圖3為前軸的結(jié)構(gòu)示意圖。危險(xiǎn)截面處的最大正應(yīng)力為σ=M/Wx，最大切應(yīng)力為τ=T/Wρ，M為彎矩，T為扭矩，Wx為截面系數(shù)，Wρ為極截面系數(shù)，且有

(16)

Wρ=0.8bt2+0.4[a3(h-2t)/t]

(17)

考慮前軸結(jié)構(gòu)靜強(qiáng)度失效，有以下極限狀態(tài)方程：

(18)

式中：σS為靜強(qiáng)度屈服極限，根據(jù)前軸材料特性有σS=460 MPa;g為裕度值。

前軸結(jié)構(gòu)尺寸參數(shù)和承受外載看作獨(dú)立正態(tài)隨機(jī)變量，分布參數(shù)如表2所示。

采用DLRS MC方法和ARBIS方法計(jì)算的新的矩獨(dú)立重要性測(cè)度指標(biāo)結(jié)果如表3所示，同時(shí)也給出了計(jì)算功能函數(shù)的總次數(shù)。可以看出，2種算法計(jì)算結(jié)果相近。圖4給出了DLRS MC方法和ARBIS方法計(jì)算結(jié)果標(biāo)準(zhǔn)差的收斂曲線?？梢?，估計(jì)值的標(biāo)準(zhǔn)差都比較小，因此計(jì)算結(jié)果可靠穩(wěn)健。DLRS MC方法的總次數(shù)為446 000，而ARBIS方法計(jì)算總次數(shù)僅為4 475，因此ARBIS方法的計(jì)算效率很高。

2種算法計(jì)算的新的矩獨(dú)立重要性測(cè)度指標(biāo)的排序均相同：t>T>b>a>h>M。這表明當(dāng)分別固定6個(gè)輸入隨機(jī)變量在其各自的分布區(qū)域內(nèi)縮減變化時(shí)，t對(duì)汽車前軸結(jié)構(gòu)的失效概率影響程度最大，而M對(duì)其失效概率的影響幾乎可以忽略不計(jì)。

圖3 汽車前軸結(jié)構(gòu)示意圖Fig.3 Schematic of automobile front axle structure

表2 汽車前軸結(jié)構(gòu)各輸入變量分布參數(shù)Table 2 Distribution parameters of input variables of automobile front axle structure

表3 DLRS MC方法和ARBIS方法計(jì)算汽車前軸指標(biāo)結(jié)果Table 3 Results of indices of automobile front axle computed by DLRS MC and ARBIS methods

圖4 DLRS MC方法和ARBIS方法計(jì)算汽車前軸的收斂曲線Fig.4 Convergence curves of automobile front axle computed by DLRS MC and ARBIS methods

5 結(jié) 論

本文提出了一種新的基于失效概率的矩獨(dú)立重要性測(cè)度分析方法來計(jì)算輸入隨機(jī)變量對(duì)結(jié)構(gòu)系統(tǒng)輸出響應(yīng)失效概率的平均影響程度。

1) 建立了一種新的區(qū)間劃分技術(shù)來等可能地獲得輸入隨機(jī)變量所有可能的縮減區(qū)間，并給出了相應(yīng)證明。

2) 引入了自適應(yīng)超球重要抽樣方法來進(jìn)行新指標(biāo)的求解，提高了計(jì)算效率。

3) 給出了一個(gè)數(shù)值算例和一個(gè)工程算例，說明了本文所提新指標(biāo)的意義，同時(shí)也驗(yàn)證了新算法的高效性。