基于Copula函數及確定優化的結構可靠度評估

趙維濤，李久安，祁武超

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

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

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

1 Sklar定理及Copula函數

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

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

(1)

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

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

(2)

式(2)中，C為與F對應的Copula函數。

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

(3)

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

對式(2)兩側求導，可得X的聯合概率密度函數為：

(4)

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

(5)

Copula函數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 改進的結構可靠性分析方法

2.1 證據理論與可靠度區間

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

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

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

m(φ)=0

∑A∈Θm(A)=1

(7)

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

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

(8)

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

2.2 基于Copula函數與證據理論的可靠度求解

2.2.1Kendall相關系數τ與Copula函數相關系數

Kendall秩相關系數[7]是在Copula函數研究中最常用到的相關性測度，主要是通過分析變量之間的變化是否協調來判斷兩變量之間的相依關系。具體以二維隨機變量為例，令(x1,y1)，(x2,y2)為隨機變量(X,Y)的兩組觀測值，若(x1-x2)(y1-y2)>0，則稱(x1,y1)與(x2,y2)是一致的，即正相關的。若(x1-x2)(y1-y2)<0，則稱(x1,y1)與(x2,y2)是不一致的，即負相關的。基于以上正負相關的討論，Kendall秩相關系數可由下式得出：

(9)

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

(10)

2.2.2Bayesian方法求最優Copula函數

(11)

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

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

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

(12)

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

(13)

綜上關于Bayesian思想的敘述，最優Copula函數的選擇標準為：

(14)

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

2.3 可靠度評估方法的改進

2.3.1嵌入確定優化的可靠度評估

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

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

(15)

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

D=X1×X2×…×Xn=

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

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

(16)

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

證據變量聯合BPA和可靠域G，結構安全的可信度Bel和似真度Pl可以通過式(8)得出。

(17)

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

(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本文流程

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

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

2) 將樣本值轉換為[0,1]上的均勻變量，并計算Kendall秩相關系數。

3) 運用Bayesian方法求出最優Copula函數。

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

5) 獲得結構可靠度區間[Bel(G),Pl(G)]。

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

3 算例

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

3.1 算例1

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

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

(20)

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

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

表1 最優Copula函數選擇

表2 邊緣BPA

表3 可靠度評估結果

3.2 算例2

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

(21)

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

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

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

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

表4 最優Copula函數選擇

表5 邊緣BPA

表6 可靠度評估結果

3.3 算例3

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

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

(22)

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

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

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

表7 最優Copula函數選擇

表8 邊緣BPA

表9 可靠度評估結果

4 結論

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

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