結構動態可靠性分析的最大熵方法

趙恩康，李洪雙

(南京航空航天大學 航空宇航學院，江蘇 南京 210016)

0 引言

結構動態可靠性分析的任務是量化計算一個結構或結構系統在工作了一定的時間t或者時段[0,t]后的可靠度。盡管針對結構動態可靠性分析已經有了許多研究工作，但是兼顧高精度和高效率仍然是一個挑戰。傳統的蒙特卡洛法(monte carlo simulation, MCS)可以用于結構動態可靠性分析，但在實際應用中由于需要大量抽樣計算，所需的計算成本太高，阻礙了它在實際工程中的應用?，F有的方差縮減技術同樣面臨大計算量困難，例如重要抽樣法[1]。因此結構動態可靠性分析領域目前主要采取近似方法，主要包括首次超越法[2-3]和極值法[4-6]。首次超越法認為極限狀態函數超過安全閾值的上界或低于其下界，結構就會發生失效。通過假設不同“超越事件”之間是相互獨立的，結構動態可靠度就可以通過在時間域上對穿越率進行積分來獲得。大多數首次超越法的局限性在于忽略了超越時間之間的相關性及線性近似帶來的誤差。極值法的基本思想是考慮極限狀態函數在目標時段內的極值，如果極值大于給定的閾值，則結構發生失效。目前，基于極值思想已經發展出多種結構動態可靠性分析方法，例如，WANG提出了一種嵌套極值響應面法[6]，DU和HU提出了一種混合高效全局優化法[7]，CHEN和WANG則提出了一種自適應極值響應面法[8]。但是現有極值法的性能依賴于代理模型對極限狀態函數的逼近精度和范圍。若不能高精度地擬合結構極限狀態函數及其極值，相應的動態可靠性分析結果是不可信的。

本文基于極值思想，提出一種結構動態可靠性分析的最大熵方法，可以直接處理一般形式的動態極限狀態函數。首先采用拓展最優線性估計(expansion optimal linear estimation, EOLE))[9]方法，將極限狀態函數中的輸入隨機過程離散化為一系列隨機變量之和的形式。對所研究結構進行結構分析，計算得出目標時間區間內極限狀態函數的極值；進而利用極值的前四階統計矩和最大熵原理擬合極值分布?；谧畲箪胤植迹蠼饨Y構動態可靠度。文中利用工程算例驗證了所提方法的有效性。

1 結構動態可靠性分析

1.1 動態可靠性定義

結構動態可靠性分析的極限狀態函數的一般形式為

g(X,Y(t),t)

(1)

其中X=[X1,X2,…,Xn] 是n個輸入隨機變量的向量；Y(t)=[Y1(t),Y2(t),…,Ym(t)] 是m個與時間有關的輸入隨機過程。如果在時間區間[t1,t2] 內存在某一時間節點t，極限狀態函數值超出安全閾值b，結構發生失效，即有

?t∈[t1,t2],g(X,Y(t),t)>b

(2)

結構失效概率可以定義為

Pf(t1,t2)=Pr(g(X,Y(t),t)>b,?t∈[t1,t2])

(3)

相應的動態可靠度為

R(t1,t2)=1-Pf(t1,t2)

(4)

由于極限狀態函數的復雜性以及時間參數的存在，得到式(3)的解析解是非常困難的。蒙特卡洛法能夠得到精確的結果，但由于需要大量樣本來評估極限狀態函數，計算成本極高。因此本文研究基于極值思想的方法。

1.2 極值法原理

極值法關注極限狀態函數一段時間內的極值函數或者極值分布。那么基于極值思想，式(3)中的結構失效概率可以表示為

Pf(t1,t2)=Pr(gmax>b)

(5)

顯然，式(5)中不含有時間參數，這就將時變問題轉化為時不變問題。但是由于極限狀態函數的復雜性，在實際工程中解析地得到其極值往往非常困難，所以目前已有研究工作主要是使用代理模型方法和近似方法得到結構響應的極值或者極值分布。

2 動態可靠性分析的最大熵方法

2.1 所提方法

為了解決現有方法的一些缺點，基于最大熵原理，本文提出一種新的結構動態可靠性分析方法。該方法的實施步驟如下：

1) 采用EOLE方法對輸入隨機過程進行離散化；

2) 采用拉丁質心Voronoi網格抽樣生成高質量抽樣樣本；

3) 計算每條樣本軌跡的極值，得到響應極值的樣本；

4) 計算極值的前四階矩，分別是均值、標準差、偏度和峰度；

5) 由極值的前四階矩計算Lagrange乘子，然后擬合極值的最大熵分布；

6) 計算動態可靠度。

本文提出的動態可靠性分析的最大熵方法流程如圖1所示。

圖1 最大熵方法流程圖

如圖2所示，陰影部分的面積即為結構失效概率。

圖2 極值分布計算失效概率示意圖

下面針對步驟中關鍵理論和技術進行詳細討論。

2.2 隨機過程的離散

在實際工程應用中，許多結構參數都具有不確定性，例如材料屬性。大多數情況下，這種不確定性都可以通過隨機變量來表征。但是仍然有一些參數隨結構服役時間的累積呈現出時變隨機特性，例如結構上的隨機載荷，需要用隨機過程來表征。若想對隨機過程Y(t)進行抽樣計算，在此之前需要對該隨機過程Y(t)進行離散化處理。本文討論的隨機過程限定為高斯隨機過程，該過程的均值函數、標準差函數和自相關系數函數分別為μY(t)、σY(t) 和ρY(t1,t2)。非高斯隨機過程可以通過轉換方法轉化為高斯隨機過程。

EOLE方法首先將時間[0,T]離散成s個離散時間點ti,i=1,2,…,s?？紤]到計算上的簡便，時間步Δt一般取定值。在時間點處多維隨機變量的相關矩陣可以表示為

(6)

(7)

式中：Ui為相互獨立的標準正態隨機變量；ρY(t)=ρY(t,ti),i=1,2,…,s。

式(7)表明：在每一個時間點ti處，隨機過程Y(ti)可以表示為一組標準正態隨機變量Ui的線性和。基于以上結論，結構動態可靠性的極限狀態函數可以轉化為g(X,Y(t),t)=g(X,U,t) ，使得極限狀態函數中不再隱式含有時間參數t。

2.3 拉丁質心Voronoi網格抽樣

抽樣樣本點的質量會影響擬合的精度。目前常用的拉丁超立方抽樣(latin hypercube sampling, LHS)是一種全空間非重疊填充的隨機抽樣方法，能夠兼顧可行域與不可行域，但難以保證抽樣的空間均勻性；BURKARDT[10]基于Voronoi網格提出了質心Voronoi網格抽樣方法 (centroidal voronoi tessellation, CVT)。它將整個抽樣空間劃分為指定數目的Voronoi區域，Voronoi區域的質心即為抽樣樣本點，具有均勻性與穩定性的特點。

ROMERO[11]等人結合LHS方法和CVT抽樣方法，提出了一種新的抽樣方法，稱為拉丁質心Voronoi網格(‘latinized’ centroidal voronoi tessellation，LCVT)抽樣方法。對質心Voronoi網格進行特殊處理，使其具有拉丁超立方體的性質，此處理過程稱為拉丁化。拉丁質心Voronoi網格抽樣的具體實施步驟如下：

(8)

式中Uji表示在第j維空間內均勻分布的隨機數。

(9)

對質心Voronoi網格進行拉丁化，獲得的新網格稱為拉丁質心Voronoi網格。LCVT抽樣方法兼具CVT方法的均勻性與LHS方法一定的隨機性，能更好地反映出隨機變量在設計空間的分布情況。

基于極值思想的結構動態可靠性分析方法的關鍵步驟是準確地計算極限狀態函數g(X,Y(t),t)的極值。獲得高質量抽樣點后，帶入結構分析模型，計算每個抽樣點對應的軌跡，并獲得結構響應在感興趣時段內的極值。進一步計算結構響應極值的前四階統計矩，為擬合最大熵分布提供約束條件。

2.4 最大熵分布

1948年，SHANNON[12]將熱力學熵的概念引入到信息科學中，來定量描述一個隨機事件的不確定性或信息量。基于信息熵的概念，JAYNE[13]在1957年提出了最大熵原理。最大熵原理可以表述為：在所有滿足給定的約束條件的概率密度函數中, 使信息熵最大的概率密度函數就是最佳，即偏差最小的概率密度函數。

若X為連續型隨機變量，其概率密度函數為f(x)，則其信息熵HX可定義為

(10)

根據最大熵原理，應該求使信息熵HX最大的概率密度函數f(x)。通常將隨機變量的前n階原點矩μi(i=1,…,n)作為約束條件，相應的優化問題為：

(11)

式(11)是一個等式約束極值問題，可以用拉格朗日乘數法求解，經過簡單推導，可以得到最大熵分布

(12)

式中λi(i=1,…,n)為Lagrange乘子。根據概率密度函數的數學定義，求得標準化因子：

在估計最大熵分布時，隨著使用的統計矩μi階次增大，對樣本中的信息利用也更加充分，最大熵分布擬合的精度也會提高。文中采用極值數據的前四階矩作為約束條件擬合極值分布，它們分別是均值、標準差、偏度和峰度。

由于本文中使用了結構響應極值的前四階矩擬合最大熵分布，無法獲得最大熵分布的顯示解，此時需要采用諸如標準牛頓法[14]的數值算法求解Lagrange乘子。但是標準牛頓法需要求解n個帶有積分的等式，過程繁瑣，計算量大，效率較低，因此文獻中使用一種無約束最小化方法求解Lagrange乘子，基本過程如下：

定義一個以λi(i=1,…,n)為自變量的勢函數Q(λ1,…,λn)

(13)

勢函數Q(λ1,…,λn)對自變量λi的梯度為：

(14)

基于式(14)，進一步計算Q(λ1,…,λn) 的Hessian矩陣元素：

(15)

3 算例分析

水動力渦輪機可以將水的動能轉化為葉片的轉動機械能。此算例為時變河流載荷下的水動力渦輪葉片[18]。渦輪葉片的簡化橫截面如圖3所示。

圖3 渦輪葉片根部橫截面

為了表征河流流速的時變特性，河流流速采用一個高斯隨機過程v(t)表示，其均值μv(t)、標準差σv(t)和自相關函數ρv(t1,t2)分別為：

(16)

(17)

ρv(t1,t2)=cos(2π(t2-t1))

(18)

式中系數a、b和c的值如下：

在河流流速v(t)作用下，葉片根部的彎矩可以表征為

(19)

式中：ρ=1×103kg/m3為水的密度；Cm=0.3422為力矩系數。

該葉片的極限狀態函數定義為

(20)

(21)

表1 渦輪葉片中隨機變量

在[0,12]個月區間上計算渦輪葉片的失效概率，通過LCVT方法生成500條樣本曲線，并與樣本數目為106的MCS計算結果進行對比，失效概率結果如表2和圖4所示，[0,12]個月區間極值的概率密度函數曲線如圖5所示。

表2 渦輪葉片失效概率

圖4 渦輪葉片失效概率

圖5 [0,12]個月極值概率密度函數曲線

由表2和圖4可知：隨著時間的推移，渦輪葉片的失效概率逐漸增大，從第7個月開始保持不變，這是由于極限狀態函數的極值落在[6,7]個月之間。在整個時間區間內，最大熵方法的計算結果都和MCS計算結果十分接近，并且相比MCS需要106次抽樣，最大熵方法只需要500條樣本曲線，說明最大熵方法兼具較好的計算精度和計算效率。

4 結語

本文基于極值思想，提出了一種結構動態可靠性分析的最大熵方法，能夠處理形如g(X,Y(t),t)的一般形式的動態可靠性極限狀態函數。該方法首先將隨機過程離散為一組隨機變量，再對極限狀態函數中的隨機變量進行抽樣并計算樣本點處的極值，最后通過樣本極值數據利用最大熵原理擬合極值分布并計算動態可靠度。從文中算例的計算結果可以看出：該方法在精度上接近蒙特卡洛方法且兼具較高計算效率，有一定的工程應用價值。