基于Copula函數(shù)及確定優(yōu)化的結(jié)構(gòu)可靠度評估

趙維濤，李久安，祁武超

(沈陽航空航天大學(xué) 航空宇航學(xué)院， 沈陽 110136)

在工程設(shè)計中，隨著對結(jié)構(gòu)可靠性要求的不斷提高，人們將結(jié)構(gòu)設(shè)計變量間的相關(guān)性和不確定性逐步納入到考慮范疇[1-2]。而且近年來的研究結(jié)果表明，證據(jù)變量間的相關(guān)性可能對可靠性計算結(jié)果產(chǎn)生較大影響，常用的獨立性假設(shè)會對可靠性計算結(jié)果產(chǎn)生較大的誤差。因此，由Sklar定理[3]發(fā)展而來的Copula函數(shù)逐步被用于處理相關(guān)性問題中，以減小由獨立性假設(shè)帶來的誤差。近年來圍繞Copula函數(shù)的研究，著重點為最優(yōu)Copula函數(shù)的選擇方法[4]和已知Copula函數(shù)形式條件下的參數(shù)估計方法[5-6]。目前最優(yōu)Copula函數(shù)的選擇常用解析法和AIC(Akaike’s information criterion)準則法，解析法需要以二元經(jīng)驗分布函數(shù)作為二維Copula函數(shù)的估計量，AIC準則法需要計算Copula函數(shù)的密度函數(shù)，這兩種方法在使用過程中都會受到不同程度的限制。李霞[7]對Archimedean Copula函數(shù)模型選擇方法的改進，充分考慮了函數(shù)的對稱性及有效估計量，可更好地刻畫了參數(shù)之間的相依性。任仙玲[8]等在基于核密度估計的條件下，提出了選擇最優(yōu)Copula函數(shù)的核密度選擇原理，該方法避免了求解各Copula函數(shù)的密度函數(shù)。David Huard[9]等引入了一種基于Bayesian思想的最優(yōu)Copula函數(shù)的選擇方法，這種方法具有良好的理論基礎(chǔ)，并且獨立于參數(shù)的選擇，可以由較少的數(shù)據(jù)得出理想的選擇結(jié)果，在使用過程中具有一定的優(yōu)勢，本文在最優(yōu)Copula函數(shù)選擇環(huán)節(jié)沿用Bayesian思想。

確定性優(yōu)化算法[10-11]廣泛應(yīng)用于工程結(jié)構(gòu)設(shè)計中，優(yōu)化結(jié)果可靠，收斂效率高，MATLAB編寫程序簡明。而暴力組合(枚舉法)在求解極限狀態(tài)函數(shù)的極值時，可能會由于節(jié)點間距選擇不當，造成計算精度差或計算效率低等問題。鑒于確定性優(yōu)化算法中的迭代求解方式可以避免以上問題出現(xiàn)，本文將該算法嵌入到基于Copula函數(shù)的結(jié)構(gòu)可靠度評估過程中，代替原有的極值求解方法，可以避免由于節(jié)點間距選擇不當而產(chǎn)生的誤差，并且提高計算效率。本文通過算例證明所提方法在使用過程中具有一定優(yōu)勢。

1 Sklar定理及Copula函數(shù)

Sklar定理[3]是在1959年由Sklar首次提出，初步展示了一維邊緣分布和聯(lián)合分布函數(shù)的關(guān)系以及相關(guān)的變化特征，闡明Copula函數(shù)在多維分布及其一維邊緣分布的關(guān)系中所扮演的重要角色。Copula函數(shù)在統(tǒng)計學(xué)中的重要性主要體現(xiàn)在Sklar定理。以二維分布為例，設(shè)H是一個聯(lián)合分布函數(shù)，其邊緣分布函數(shù)分別是FX(x)和RY(y)，則一定有一個Copula函數(shù)C，使得對實數(shù)域中所有的x和y有式(1)成立，若FX(x)和RY(y)是連續(xù)的，則C是唯一存在的，否則C在實數(shù)域中不是唯一確定的。

H(x,y)=C(FX(x),RY(y))

(1)

將Sklar定理推廣到n維，邊緣分布為F1(x1),F2(x2),…,Fn(xn)的n維聯(lián)合分布函數(shù)可以表示為：

F(x1,x2,…,xn)=C(F1(x1),F2(x2),…,Fn(xn))

(2)

式(2)中，C為與F對應(yīng)的Copula函數(shù)。

如果邊緣分布函數(shù)均為連續(xù)函數(shù)，則存在唯一的Copula函數(shù)C使得下式成立：

(3)

式(3)中，u1=F1(x1),u2=F2(x2),…,un=Fn(xn)，且滿足ui∈[0,1],i=1,2,…,n。

對式(2)兩側(cè)求導(dǎo)，可得X的聯(lián)合概率密度函數(shù)為：

(4)

式(4)中，c為Copula函數(shù)C的概率密度函數(shù)。其求解關(guān)系式如下：

(5)

Copula函數(shù)C(u1,u2,…,un)的上下界可由Fréchet-Hoeffding邊界[12]給出，具體為：

W(u1,u2,…,un)≤C(u1,u2,…,un)≤M(u1,u2,…,un)

(6)

式(6)中：M(u1,u2,…,un)=min(u1,u2,…,un)，W(u1,u2,…,un)=max(u1+u2+…+un-n+1,0)。

2 改進的結(jié)構(gòu)可靠性分析方法

2.1 證據(jù)理論與可靠度區(qū)間

D-S證據(jù)理論即Dempster-Shafer理論[13]，可以對各類不確定信息進行合理地描述和處理，旨在減小由多重不確定因素耦合作用造成的影響，是處理不確定信息的有力工具。將證據(jù)理論與可靠度評估結(jié)合，能夠很好地解決涉及多源不確定因素的可靠度評估問題，對可靠度區(qū)間上下限進行準確量化，而且具有較高的求解效率。

關(guān)于證據(jù)理論，最重要的基礎(chǔ)理論之一為基本可信度分配函數(shù)(Basic Probability Assignment，BPA)。令2Θ表示識別框架Θ的冪集，由Θ包含的所有命題可能組成。證據(jù)理論用基本可信度分配表示對命題A的信任程度，BPA是一個滿足如下公理的映射：

0≤m(A)≤1, ?A∈2Θ

m(φ)=0

∑A∈Θm(A)=1

(7)

式(7)中：φ表示空集，對m(A)>0的子集A稱為焦元，命題A的基本可信度分配值m(A)表示證據(jù)對命題A的支持程度。

根據(jù)證據(jù)理論，可在缺乏信息的情況下，利用區(qū)間[Bel(A),Pl(A)]來描述證據(jù)對命題的支持程度，其中可信度與似真度可以由下式得出

(8)

式(8)中：B為同一識別框架下的已知證據(jù)；Bel(A)稱為可信度，是完全支持A的證據(jù)BPA之和；Pl(A)為似真度，是所有完全或部分支持A的證據(jù)BPA之和。

2.2 基于Copula函數(shù)與證據(jù)理論的可靠度求解

2.2.1Kendall相關(guān)系數(shù)τ與Copula函數(shù)相關(guān)系數(shù)

Kendall秩相關(guān)系數(shù)[7]是在Copula函數(shù)研究中最常用到的相關(guān)性測度，主要是通過分析變量之間的變化是否協(xié)調(diào)來判斷兩變量之間的相依關(guān)系。具體以二維隨機變量為例，令(x1,y1)，(x2,y2)為隨機變量(X,Y)的兩組觀測值，若(x1-x2)(y1-y2)>0，則稱(x1,y1)與(x2,y2)是一致的，即正相關(guān)的。若(x1-x2)(y1-y2)<0，則稱(x1,y1)與(x2,y2)是不一致的，即負相關(guān)的?；谝陨险撓嚓P(guān)的討論，Kendall秩相關(guān)系數(shù)可由下式得出：

(9)

式(9)中：τ是由觀測值計算得出的秩相關(guān)系數(shù)；sign(·)為符號函數(shù)。

(10)

2.2.2Bayesian方法求最優(yōu)Copula函數(shù)

(11)

式(11)中：Cl為第l個備選Copula函數(shù)；P(Cl)是在得出觀測數(shù)據(jù)之前，Copula函數(shù)為Cl的先驗概率；P(Q|Cl)為在備選Copula函數(shù)Cl條件下，數(shù)據(jù)Q的似然。

在備選Copula函數(shù)給定的情況下，P(Q)為常數(shù)，所以式(11)等價于：

P(Cl|Q)∝P(Q|Cl)P(Cl)

(12)

在實際工程設(shè)計工作中，由于對結(jié)構(gòu)參數(shù)的未知性，缺乏具有說服力的先驗認識，可以將P(Cl)取為相同的概率，因此最優(yōu)Copula函數(shù)完全由P(Q|Cl)確定。當參數(shù)θl為離散值時，根據(jù)似然函數(shù)的定義，P(Q|Cl,θl)具體表達式為：

(13)

綜上關(guān)于Bayesian思想的敘述，最優(yōu)Copula函數(shù)的選擇標準為：

(14)

式(14)中，arg max(·)表示求解具有最大函數(shù)值的參量。

2.3 可靠度評估方法的改進

2.3.1嵌入確定優(yōu)化的可靠度評估

基于上述理論基礎(chǔ)，本文創(chuàng)新性地將確定性優(yōu)化算法嵌入到可靠度評估過程中，將原有的可靠度評估方法進行改進。首先考慮，一個含有n維證據(jù)矢量X的不確定性結(jié)構(gòu)，其可靠域G可定義為[15]

G{g∶g(x1,x2,…,xn)≥0}

(15)

式(15)中，g為功能函數(shù)，定義一個n維識別框架如下：

D=X1×X2×…×Xn=

{dk=[x1i,x2j,…,xnl];

x1i∈X1,x2j∈X2,…,xnl∈Xn}

(16)

式(16)中：[x1i,x2j,…,xnl]構(gòu)成第k個焦元dk，i,j,…,l為證據(jù)變量所分區(qū)間數(shù)量，且滿足k=i×j×…×l。

證據(jù)變量聯(lián)合BPA和可靠域G，結(jié)構(gòu)安全的可信度Bel和似真度Pl可以通過式(8)得出。

(17)

理論上，真實的結(jié)構(gòu)可靠度P=P{g(x1,…,xn)≥0}應(yīng)該屬于區(qū)間[Bel(G),Pl(G)]。為計算上述兩個可靠性測度Bel和Pl則需要判斷dk?G(焦元dk完全位于可靠域內(nèi))或dk∩G≠φ(焦元d完全位于或者部分位于可靠域內(nèi))，為此需要計算極限狀態(tài)方程在每個焦元dk上的極值。本文將利用確定性優(yōu)化算法求解極值，該方法具有可行性，并且相比較簡單極值求解，具有一定優(yōu)勢確定性優(yōu)化模型如下：

(18)

(19)

對于焦元dk，如果gmin和gmax均為正，則dk?G，焦元dk的BPA同時計入Bel(G)和Pl(G)中；如果gmin和gmax均為負，則dk∩G=φ焦元dk的BPA既不計入Bel(G)，也不計入Pl(G)中；如果gmin為負，gmax為正，則dk∩G≠φ，焦元dk的BPA計入Pl(G)但不計入Bel(G)。

2.3.2本文流程

基于以上所提的理論基礎(chǔ)，總結(jié)得出本文方法的總體流程如下：

1) 設(shè)計變量X樣本值，邊緣BPA。

2) 將樣本值轉(zhuǎn)換為[0,1]上的均勻變量，并計算Kendall秩相關(guān)系數(shù)。

3) 運用Bayesian方法求出最優(yōu)Copula函數(shù)。

4) 利用優(yōu)化算法式(18)求解出每個焦元上極限狀態(tài)函數(shù)極值[gmin,gmax]。

5) 獲得結(jié)構(gòu)可靠度區(qū)間[Bel(G),Pl(G)]。

本文所提出的改進方法流程圖如圖1所示。

3 算例

本文所提出的改進可靠度評估改進方法涉及到確定性優(yōu)化算法，因此在極限狀態(tài)函數(shù)極值計算時調(diào)用MATLAB優(yōu)化工具箱中的fmincon[16-17]函數(shù)求解。

3.1 算例1

矩形截面梁[18]載荷位置與截面尺寸如圖2所示，本算例根據(jù)對該梁固定端最大應(yīng)力處進行強度校核，得出極限狀態(tài)方程如下：

g(x1,x2)=S-3 000x1-1 500x2

(20)

式(20)中：x1和x2為證據(jù)變量，也是確定性優(yōu)化算法中的優(yōu)化變量；S為懸臂梁極限強度。

在本算例中，證據(jù)變量離散點數(shù)據(jù)來源于文獻[18]，并對材料參數(shù)進行調(diào)整后，將自編程序計算結(jié)果作為對比參考。將證據(jù)變量x1,x2轉(zhuǎn)換為在[0,1]上的均勻變量。然后由式(9)計算得出Kendall秩相關(guān)系數(shù)τ=0.683 3。根據(jù)式(10)中給出的Copula函數(shù)相關(guān)性系數(shù)θ與Kendall秩相關(guān)系數(shù)τ的微積分關(guān)系，求得各Copula函數(shù)的相關(guān)性系數(shù)θ見表1。對于該算例，沿用Bayesian方法選取最優(yōu)Copula函數(shù)，各Copula函數(shù)的權(quán)重比由表1給出，符合證據(jù)變量的最優(yōu)Copula函數(shù)為Gumbel Copula函數(shù)。證據(jù)變量x1,x2的區(qū)間分布與邊緣BPA由表2給出。當充分考慮證據(jù)理論相關(guān)性模型與最優(yōu)Copula函數(shù)的內(nèi)在關(guān)系后，本文將確定性優(yōu)化算法的表達式(18)引入可靠性評估過程中，得出極限狀態(tài)函數(shù)在各焦元中的極值。在此基礎(chǔ)上，懸臂梁可靠度評估結(jié)果如表3所示，從結(jié)果中可以明顯看出，本文算法在保證計算精度的條件下，很大程度地提高了計算效率。

表1 最優(yōu)Copula函數(shù)選擇

表2 邊緣BPA

表3 可靠度評估結(jié)果

3.2 算例2

某齒輪減速箱的齒輪組及箱體如圖3所示，本算例中例包含三個證據(jù)變量，分析后得出的極限狀態(tài)函數(shù)如下：

(21)

式(21)中：φ1為輸出軸直徑；φ2為輸入軸直徑；x1為輸出軸固定軸承的間距。

本算例因為無法進行試驗操作得出數(shù)據(jù)，所以通過MATLAB隨機生成符合Gaussian Copula函數(shù)分布特征的離散數(shù)據(jù)點，作為可靠度評估的數(shù)據(jù)來源。通過后續(xù)的計算分析，也間接驗證了本文最優(yōu)Copula函數(shù)選擇方法的有效性與正確性。

在本算例中，將進一步計算包含三維證據(jù)變量x1，φ1和φ2的齒輪減速機構(gòu)的可靠度區(qū)間。首先根據(jù)證據(jù)變量離散數(shù)據(jù)點的一致性概率求出Kendall秩相關(guān)系數(shù)為τ=0.508 0，然后由Copula相關(guān)性系數(shù)θ與Kendall秩相關(guān)系數(shù)τ的微積分關(guān)系，求得各Copula函數(shù)的相關(guān)性系數(shù)見表4，同時表4也給出了由Bayesian方法求出的各Copula函數(shù)的權(quán)重比，本算例最優(yōu)Copula函數(shù)為Gaussian Copula函數(shù)。證據(jù)變量x1，φ1和φ2的區(qū)間分布及邊緣BPA，由表5給出。在此基礎(chǔ)上，齒輪減速機構(gòu)可靠度評估過程及結(jié)果如表6所示。

由表6所示的計算過程和結(jié)果對比可以看出，當節(jié)點間距較大時，暴力組合(枚舉法)可以提高計算效率，但是會引起可靠度評估結(jié)果不穩(wěn)定。隨著節(jié)點間距逐漸減小，可靠度評估結(jié)果逐漸趨于穩(wěn)定，但是計算量較大，計算耗時長。當引進確定性優(yōu)化迭代計算過程代替極值求解方法，可以通過少量迭代計算得出較為理想的評估結(jié)果，該方法在保證計算精度的條件下，很大程度地提高了計算效率。

表4 最優(yōu)Copula函數(shù)選擇

表5 邊緣BPA

表6 可靠度評估結(jié)果

3.3 算例3

十桿桁架結(jié)構(gòu)[19-20]如圖4所示，圖4中水平方向和豎直方向桁架長度L=9.144 m，材料彈性模量E=68 948 MPa。其中1～6桿橫截面積為70 cm2，7～10桿橫截面積為65 cm2，節(jié)點4受到豎直作用力x1，節(jié)點2受到豎直作用力x2和水平作用力x3。本算例將節(jié)點2處的豎向位移d2max作為最大位移極限，隱式極限狀態(tài)方程如下：

g(x1,x2,x3)=d2max-δ2(x1,x2,x3)

(22)

式(22)中：δ2(x1,x2,x3)為節(jié)點2處的豎向位移。

在本算例中，x1,x2,x3離散數(shù)據(jù)生成原理同算例2，此處不再贅述。對十桿桁架結(jié)構(gòu)進行可靠度評估，著重研究本文算法與其他算法在計算精度及計算效率方面的比較。首先，由證據(jù)變量x1,x2,x3的離散數(shù)據(jù)點的一致性概率得出秩相關(guān)系數(shù)τ=0.505 1，各備選Copula函數(shù)的相關(guān)性系數(shù)θ見表7，同時表7也給出了各備選Copula函數(shù)的權(quán)重比，本算例的最優(yōu)Copula函數(shù)為t Copula。證據(jù)變量x1,x2,x3的區(qū)間分布及邊緣BPA由表8給出。

本文方法與其他方法的計算結(jié)果見表9，從表9中可以看出，本文方法與文獻[18]方法求得的可靠度區(qū)間一致，說明本文方法具有一定的計算精度。文獻[20]中的響應(yīng)面法選取106個樣本點進行分析，得出可靠概率Pr=0.879 1，文獻[21]靜態(tài)響應(yīng)矩法計算得出的可靠概率Pr=0.888 2?？梢钥闯?，所涉及文獻中計算得出的可靠度都介于本文方法計算得出的可靠度區(qū)間內(nèi)。通過該算例可知，本文方法在求解可靠度區(qū)間問題是可行的。

表7 最優(yōu)Copula函數(shù)選擇

表8 邊緣BPA

表9 可靠度評估結(jié)果

4 結(jié)論

1) 本文提出了一種高效的可靠度評估方法，將確定性優(yōu)化算法巧妙地嵌入到可靠度評估過程中，將原有的可靠性分析方法進行改進，避免由于暴力組合(枚舉法)過程中節(jié)點間距選擇不當帶來的誤差，可在保證計算精度條件下提高計算效率。

2) 本文方法適用于極限狀態(tài)函數(shù)為線性函數(shù)、非線性函數(shù)和隱函數(shù)的各種情況。具有一定的求解精度，而且在很大程度上提高了計算效率。