基于高效求矩法與最大熵原理的矩方法

姜昊銳，趙衍剛

(北京工業(yè)大學(xué)城市建設(shè)學(xué)部，北京 100124)

工程結(jié)構(gòu)在使用過程中存在著諸如材料參數(shù)、幾何參數(shù)以及外部荷載作用等不確定性，日益引發(fā)了人們對基于概率設(shè)計的重視，結(jié)構(gòu)可靠度分析在結(jié)構(gòu)設(shè)計與使用過程中起到越來越重要的作用[1]。目前常用的可靠度分析方法有蒙特卡洛模擬法[2]、一次二階矩法(first order reliability method, FORM)[3-5]與二次二階矩法(second order reliability method, SORM)[6-8]等,但蒙特卡洛模擬方法在處理小失效概率問題時需要大量的計算,效率低下；FORM與SORM因為需要計算導(dǎo)數(shù)以及迭代運算，在處理復(fù)雜工程問題時精度較低。矩方法也是可靠性分析方法之一，因其方法簡單，便于操作而受到越來越多的重視。矩方法具體包含兩個步驟：首先求解功能函數(shù)的統(tǒng)計矩，如前四階中心矩：均值，標(biāo)準(zhǔn)差，偏度，峰度；接著依據(jù)求得的統(tǒng)計矩信息，選擇合適的分布模型近似真實的概率分布，最終在失效域上進(jìn)行積分求解結(jié)構(gòu)的失效概率與可靠指標(biāo)。由上述步驟可知，矩方法不需要像蒙特卡洛模擬法進(jìn)行大量運算，也無需像FORM或SORM一樣計算導(dǎo)數(shù)以及迭代運算，但求矩的效率與精度直接影響著矩方法的精度與效率[7-9]，如何高效且準(zhǔn)確地求解功能函數(shù)的矩,目前仍舊是矩方法所研究的重點。

求解功能函數(shù)的統(tǒng)計矩，實質(zhì)上就是計算多維積分。在實際工程中，由于結(jié)構(gòu)的功能函數(shù)通常較為復(fù)雜且以隱式函數(shù)表示，直接進(jìn)行多維積分運算往往較為困難[10]。因此，許多近似計算方法相繼被提出。Zhao等[11]提出了新點估計方法，通過Rosenblatt變換[12]和一維減維方法[13]將功能函數(shù)表示成一元函數(shù)和的形式，繼而用一維點估計進(jìn)行求矩運算，很大程度上提高了計算的效率，但該方法在功能函數(shù)較為復(fù)雜的情況下計算精度較低。為了提高計算精度，Zhao等[14]提出了將Rosenblatt變換和二維減維方法[15]相結(jié)合的二維減維點估計方法，將功能函數(shù)表示成一元函數(shù)與二元函數(shù)和的形式進(jìn)行求矩運算，但隨著輸入隨機變量數(shù)目的增加，該方法往往會造成“維數(shù)災(zāi)難”問題。Xu等[16]提出了用高階無跡變換[17-19]計算二維減維點估計中的二維積分問題，該求矩方法可一定程度上提高計算精度與效率。為了進(jìn)一步提高計算的效率，F(xiàn)an等[20]提出了自適應(yīng)二維減維點估計方法，通過引入交叉項判定準(zhǔn)則，判斷二維減維后的二元函數(shù)中是否存在交叉項，若不存在交叉項，則將二元函數(shù)進(jìn)一步用一元函數(shù)表示，并用一維點估計法進(jìn)行求矩運算；若存在交叉項，則用二維點估計進(jìn)行計算，但在功能函數(shù)所含交叉項較多的情況下，該方法依舊需要大量的運算。針對此問題，Cai等[21]提出采用二維稀疏網(wǎng)格法[22]計算自適應(yīng)二維減維點估計方法中的二維積分，可一定程度上提高計算效率，但在輸入隨機變量數(shù)目較多的情況下計算量依舊較大。綜上所述，已有的求矩方法均不能在計算過程中兼顧精度與效率。

鑒于此，現(xiàn)提出一種基于快速求矩法與最大熵原理的矩方法，將自適應(yīng)二維減維與高階無跡變換相結(jié)合求解結(jié)構(gòu)功能函數(shù)的前四階矩，并選擇最大熵原理[11]近似真實的概率分布，最終依據(jù)最大熵分布求解結(jié)構(gòu)的可靠指標(biāo)。算例部分通過兩道例題，對比所提方法與幾種已有方法的計算結(jié)果，驗證所提方法在求矩運算以及可靠指標(biāo)計算問題中在效率以及精度上的優(yōu)勢。

1 基于矩方法的可靠度分析

1.1 問題表述

結(jié)構(gòu)的功能函數(shù)[5]的表達(dá)式通常為

Z=G(X)

(1)

式(1)中：Z為功能函數(shù)；X=[X1,X2,…,Xn]為n維隨機向量，功能函數(shù)的前四階中心矩計算公式為

μG=μ1G

(2)

(3)

(4)

(5)

式中：μG、σG、α3G和α4G分別為功能函數(shù)的均值、標(biāo)準(zhǔn)差、偏度與峰度；μkG為功能函數(shù)的k階原點矩(k=1,2,3,4),具體計算表達(dá)式為

(6)

式(6)中：fx(x)為輸入隨機向量的聯(lián)合概率密度函數(shù)。由式(6)可知，求矩運算本質(zhì)上是多維積分計算問題，但實際工程中功能函數(shù)較為復(fù)雜，且隨機向量的聯(lián)合概率密度函數(shù)難以直接獲得，因此直接求解并不可行。為了高效且精確地求矩，提出基于自適應(yīng)二維減維與高階無跡變換的高效求矩法。

1.2 自適應(yīng)二維減維

結(jié)構(gòu)可靠性分析一般在標(biāo)準(zhǔn)正態(tài)空間下進(jìn)行，因此可通過Rosenblatt變換或者Nataf變換[23]，將輸入的具有相關(guān)性的非正態(tài)隨機變量轉(zhuǎn)換為獨立的標(biāo)準(zhǔn)正態(tài)隨機變量U=[U1,U2,…,Un]，此時功能函數(shù)可用標(biāo)準(zhǔn)正態(tài)隨機變量[18]表示為

Z=G(X)=G[T-1(U)]=L(U)

(7)

式(7)中：T-1為Rosenblatt逆變換或Nataf逆變換；L為G與T-1的復(fù)合函數(shù)。

基于二維減維[12]方法，功能函數(shù)可進(jìn)一步寫為

(8)

式(8)中：L2為n(n-1)/2個二元函數(shù)的和；L1為n個一元函數(shù)的和；L0為功能函數(shù)在均值點處的函數(shù)值，為一常數(shù)。L2、L1、L0的具體表達(dá)式為

(9)

(10)

L0=L(0,…,0,…,0)

(11)

式中：i,j=1,2,…,n,i

基于式(8)～式(11)可得，功能函數(shù)在二維減維后的k階原點矩的近似計算表達(dá)式為

(12)

(13)

(14)

式中：φ(u)為標(biāo)準(zhǔn)正態(tài)隨機變量的概率密度函數(shù)。由式(14)易知求解μk-Li為一維積分問題，可通過一維點估計法[15]進(jìn)行計算，即

(15)

式(15)中：ur為積分點；Pr為對應(yīng)權(quán)重；N為積分點個數(shù)。一維點估計中積分點及其對應(yīng)權(quán)重的計算式[16]為

(16)

(17)

式中：xr、ωr分別為以exp(-x2)為權(quán)函數(shù)的2N-1階高斯埃爾米特積分的積分點與權(quán)重值[24]。

式(13)中含有兩個隨機變量，為二維積分問題，文獻(xiàn)[21]采用二維點估計進(jìn)行運算，即

(18)

針對此情況，自適應(yīng)二維減維點估計方法[17]提出了交叉項判定準(zhǔn)則，通過判定準(zhǔn)則判斷函數(shù)中是否存在交叉項，將判定為不含有交叉項的二元函數(shù)進(jìn)一步用一元函數(shù)表示，對應(yīng)的矩的計算可直接利用一維積分的計算結(jié)果；判定為存在交叉項的二元函數(shù)，其矩仍通過二維點估計進(jìn)行計算。

判斷二元函數(shù)中是否存在交叉項的具體步驟如下。

(1)首先選取參考點Ur=(uir，uir),并帶入功能函數(shù)計算，用Lij(Ur)表示。

(2)另選取3個參考點，分別記為Ui-1=(ui-1，ujr),Uj-1=(uir，uj-1)以及Uij-1=(ui-1，uj-1)，代入式(19)計算Δ1，即

(19)

若Δ1>ε，則可判定在二元函數(shù)Lij(Ui，Uj)中Ui，Uj存在交叉項，否則執(zhí)行下一步。

(3)再選取3個計算點，分別記為Ui-2=(ui-2，ujr),Uj-2=(uir，uj-2)以及Uij-2=(ui-2，uj-2)，代入式(20)計算Δ2，即

(20)

若Δ2>ε，則可判定在二維函數(shù)Lij(Ui，Uj)中Ui，Uj存在交叉項，否則Lij(Ui，Uj)中Ui，Uj不存在交叉項。

在上述步驟中ε為一特定值，作為判別交叉項存在的標(biāo)準(zhǔn)，取ε=10-6。值得注意的是，交叉項判定選取的計算點通常是進(jìn)行一維點估計計算的積分點，因此不需額外增加計算量。

如果二維函數(shù)Lij(Ui，Uj)經(jīng)上述判別準(zhǔn)則判定為不存在交叉項時，則二元函數(shù)可用一元函數(shù)表示為

Lij(Ui,Uj)=Li(Ui)+Lj(Uj)-L0

(21)

繼而不含交叉項的二元函數(shù)矩的求解，可直接利用一元函數(shù)Li(Ui)、Lj(Uj)矩的計算結(jié)果進(jìn)行計算，即

μ1-Lij=μ1-Li+μ1-Lj-L0

(22)

(23)

(24)

(25)

式中：M2ij、M3ij、M4ij計算表達(dá)式分別為

(26)

(27)

(28)

由式(22)～式(28)可知，此時二維積分計算可直接利用一維積分計算結(jié)果，無需額外增加計算量，因此一定程度上提高了求矩的效率。

若判定二元函數(shù)中存在交叉項，則采用二維點估計進(jìn)行計算，為保證精度通常采用5點或7點估計，此時每個二元函數(shù)的求矩過程分別需計算25次或49次功能函數(shù)值。在交叉項較多的情況下，進(jìn)行二維點估計求矩的二元函數(shù)數(shù)目較多，因此仍舊需要大量的運算。針對此情況，通過采用高階無跡變換方法替換原有的二維點估計，減少了二維積分中積分點的個數(shù)，相應(yīng)減少了功能函數(shù)的計算次數(shù)，一定程度上提高了求矩的計算效率。

1.3 基于高階無跡變換的二維積分算法

高階無跡變換可以用于解決n維高斯權(quán)重積分問題[21]。高階無跡變換的基本思想是[18]：首先，選取一些特定的積分點(Sigma點)及其對應(yīng)的權(quán)重來匹配輸入隨機變量的概率信息，通常為輸入隨機變量的各階統(tǒng)計矩；其次，將選取的積分點通過非線性變換，即帶入結(jié)構(gòu)的功能函數(shù)計算對應(yīng)的輸出狀態(tài)變量的樣本點；最后，輸出狀態(tài)變量的統(tǒng)計矩可以利用輸入變量的積分點對應(yīng)權(quán)重和輸出狀態(tài)變量的樣本點進(jìn)行估計，其實質(zhì)上為一種改進(jìn)的高斯積分。

由前述可知，高階無跡變換積分點總個數(shù)為N=2n2+1，在自適應(yīng)二維減維后的二維積分問題中，n=2，則每個二元函數(shù)矩的計算僅需9次運算，其效率明顯高于二維點估計中的5點或7點估計。

高階無跡變換三類積分點及其對應(yīng)權(quán)重的具體表達(dá)式[13,15]如下。

第一類：

θ0=0

(29)

(30)

第二類：

(31)

(32)

(33)

式中：ej1=[0,…,1,…,0]，即只有第j1項為1。

第三類：

(34)

(35)

(36)

(37)

(38)

在計算第二類和第三類積分點時，需要引入自由參數(shù)β，其影響著高階無跡變換的計算結(jié)果的精度，文獻(xiàn)[16,18]中對β的取值進(jìn)行過研究，在此取β=7.2，經(jīng)不同例題試算結(jié)果表明，所選取的β值計算結(jié)果精度較高。

當(dāng)自由參數(shù)確定后，三類積分點及權(quán)重由式(29)～式(38)進(jìn)行計算，此時含有交叉項的二元函數(shù)的矩可直接由式(39)進(jìn)行計算，即

(39)

式(39)中：ωi、θi分別為高階無跡變換積分點與對應(yīng)權(quán)重，i=1,2,…,n。

最后，將計算所得的自適應(yīng)二維減維后的一元函數(shù)、二元函數(shù)的矩代入式(12)中，功能函數(shù)的前四階原點矩μkG即可確定(k=1,2,3,4),繼而可通過式(2)～式(5)計算功能函數(shù)G(X)的前四階中心矩。

獲得功能函數(shù)的前四階中心矩后，需選擇合適的分布模型對功能函數(shù)的真實概率分布進(jìn)行近似，并以此計算結(jié)構(gòu)可靠指標(biāo)。常用的分布模型有Pearson分布族、Burr分布族等，盡管該類分布系統(tǒng)非常靈活，但在不同分布類型的分界處往往難以實現(xiàn)分布近似。因此，采用最大熵原理近似結(jié)構(gòu)功能函數(shù)的概率分布，避免了上述問題，同時該分布模型對真實概率分布的擬合效果較好，且分布參數(shù)計算方便。

1.4 基于最大熵原理的可靠指標(biāo)求解

若功能函數(shù)Z服從概率密度函數(shù)為fZ(z)的連續(xù)分布，其熵定義[25-26]為

(40)

根據(jù)最大熵原理可知，最大熵分布為所有可能的分布形式中，在給定約束條件下,使熵最大的概率分布，其具有最小偏見，因此被認(rèn)為是構(gòu)造“最優(yōu)”概率分布的一種途徑。

為使式(40)的熵取最大值，此處以計算的功能函數(shù)前四階原點矩為約束條件。在計算最大熵分布時，為求解方便，提高收斂速度，通常先將功能函數(shù)進(jìn)行標(biāo)準(zhǔn)化處理為

(41)

此時對功能函數(shù)Z的分布近似轉(zhuǎn)變?yōu)閷?biāo)準(zhǔn)化變量Y的概率分布的近似。

最大熵分布的求解可通過引入拉格朗日算法[26]進(jìn)行計算，表達(dá)式為

(42)

令式(42)關(guān)于拉格朗日乘子的偏導(dǎo)數(shù)為0，即

(43)

(44)

利用所得到的最大熵分布，在失效域上積分可計算得到結(jié)構(gòu)的失效概率為

(45)

對應(yīng)的結(jié)構(gòu)可靠度指標(biāo)計算公式[3]為

β=Φ-1(1-pf)

(46)

式(46)中：Φ-1為標(biāo)準(zhǔn)高斯隨機變量的累計分布函數(shù)的逆函數(shù)。

綜上，以自適應(yīng)二維減維與高階無跡變換求得的功能函數(shù)前四階矩為約束條件，結(jié)合最大熵原理，最終便可計算得到結(jié)構(gòu)的可靠指標(biāo)。

1.4 算法步驟

所提出的一種新的矩方法，基于自適應(yīng)二維減維與高階無跡變換高效求解功能函數(shù)的前四階中心矩，并采用最大熵原理求解結(jié)構(gòu)的可靠指標(biāo)。具體計算步驟如下。

(1)根據(jù)式(8)將結(jié)構(gòu)功能函數(shù)進(jìn)行二維減維，表示成二元函數(shù)與一元函數(shù)和的形式。

(2)一元函數(shù)的矩直接通過式(15)的一維點估計進(jìn)行運算求解。

(3)引入交叉項判定準(zhǔn)則：①通過準(zhǔn)則判定為不含有交叉項的二元函數(shù)，利用式(21)進(jìn)一步用一元函數(shù)表示，其對應(yīng)的矩可利用已求得的一元函數(shù)矩的結(jié)果直接進(jìn)行計算；②判定為含有交叉項的二元函數(shù)，其矩采用高階無跡變換方法進(jìn)行計算。

(4)將求得的一元函數(shù)、二元函數(shù)的矩的計算結(jié)果代入式(12)中，即可計算功能函數(shù)的前四階原點矩，繼而由式(2)～式(4)計算前四階中心矩。

(5)以計算得到的前四階中心矩為約束條件，基于最大熵原理求解功能函數(shù)的最大熵分布，繼而通過最大熵分布在失效域上的積分，可最終求解結(jié)構(gòu)的失效概率與可靠指標(biāo)。

2 算例分析

為驗證所提方法在計算精度與效率上的優(yōu)勢，選取兩道算例，將所提高效求矩方法計算得到的功能函數(shù)前四階中心矩結(jié)果，與蒙特卡洛模擬方法、一維減維點估計方法、二維減維點估計方法、自適應(yīng)二維減維點估計、二維減維無跡變換及自適應(yīng)二維減維稀疏網(wǎng)格積分方法進(jìn)行了對比。此外，將基于最大熵原理計算得到的可靠指標(biāo)與蒙特卡洛方法結(jié)果進(jìn)行了對比。

2.1 算例1：軸向壓力環(huán)形柱可靠度分析

某根受軸向壓力的環(huán)形柱，其受力圖與平面幾何尺寸如圖1所示，其中環(huán)形柱的彈性模量為E，長度為L，內(nèi)徑為D，壁厚為T，受到軸向荷載P的作用。其中，E、L、D、T、P為相互獨立的隨機變量，具體對應(yīng)的概率分布及分布參數(shù)如表1所示。

圖1 軸向壓力環(huán)形柱示意圖

表1 各隨機變量的分布信息

此處受軸向壓力環(huán)形柱的失效模式為屈曲破壞，其對應(yīng)的極限狀態(tài)函數(shù)表達(dá)式為

(47)

表2中分別列出了蒙特卡洛模擬方法，一維減維點估計方法、二維減維點估計方法，自適應(yīng)二維減維點估計，二維減維無跡變換及自適應(yīng)二維減維稀疏網(wǎng)格積分方法的求矩結(jié)果與計算次數(shù)。同時以蒙特卡洛模擬法結(jié)果為標(biāo)準(zhǔn)，不同方法計算結(jié)果的相對誤差也在表內(nèi)括號里列出。其中，相對誤差計算公式為

(48)

式(48)中：e為相對誤差；rMCS為蒙特卡洛模擬方法計算結(jié)果；r為各求矩方法計算結(jié)果。

由表2可知，一維減維點估計方法計算效率較高，但偏度與峰度的計算誤差較大。所提求矩方法相較于其他方法，計算效率有一定程度的提升，僅需計算59次，且提高效率的同時，計算的精度也能保證，且誤差更小。其中，本文方法計算得到的前四階中心矩中偏度的相對誤差最大，但僅有1.31%。

表2 各求矩方法計算結(jié)果

求解得到前四階中心矩后，便可按照1.4節(jié)步驟計算對應(yīng)的最大熵分布，繼而進(jìn)行可靠指標(biāo)計算。表3中列出了本文方法結(jié)合最大熵分布計算得到的可靠指標(biāo)結(jié)果，并與蒙特卡洛方法計算結(jié)果進(jìn)行對比，以驗證所求結(jié)果的精度。同樣，采用式(12)計算本文方法計算得到的可靠指標(biāo)相對于蒙特卡洛模擬方法結(jié)果的相對誤差。

表3 可靠指標(biāo)計算結(jié)果

從表3可知，本文方法求得的可靠指標(biāo)與蒙特卡洛較為吻合，在保留三位有效數(shù)字情況下結(jié)果與蒙特卡洛方法結(jié)果基本一致，相對誤差僅有0.7%。由此可知本文方法結(jié)合最大熵原理進(jìn)行可靠指標(biāo)計算時精度較高，符合誤差范圍。

2.2 算例2：空間剛架可靠度分析

如圖2所示，有一個由72根桿件構(gòu)成的空間桁架，其中各個桿件的橫截面面積均為34.849 cm2,F1、F2、F3、F4、F5、F6、F7、F8分別為作用在節(jié)點5、8、9、12、13、16、17、20處的水平作用力，所有作用力與構(gòu)成空間剛架的桿件的彈性模量E互為獨立隨機變量，對應(yīng)的概率分布與分布參數(shù)如表4所示。

圖2 空間72桿框架示意圖

表4 各隨機變量的分布信息

此處考慮該空間桁架的極限狀態(tài)為所有節(jié)點的水平位移絕對值最大值超過界限值，該失效模式所對應(yīng)的結(jié)構(gòu)功能函數(shù)表達(dá)式為

Z=G(X)=0.05-max|μi(x)|,i=5,6,…,20

(49)

式(49)中：0.05為沿x方向水平位移限值，m。本例題計算過程中，空間桁架的各個節(jié)點水平位移通過MATLAB軟件進(jìn)行有限元分析計算。表5中列出了各個求矩方法計算得到的功能函數(shù)前四階中心矩，同算例1，表中括號內(nèi)也列出了各個求矩方法計算結(jié)果相較于蒙特卡洛模擬方法結(jié)果的相對誤差。

在計算得到功能函數(shù)前四階中心矩后，依據(jù)最大熵原理計算相應(yīng)的結(jié)構(gòu)可靠指標(biāo)并與蒙特卡洛方法對比，具體結(jié)果如表6所示。

表6 可靠指標(biāo)計算結(jié)果

由表5可知，在求解隱式功能函數(shù)(14)的前四階中心矩時，本文方法計算功能函數(shù)次數(shù)較少，僅需計算211次功能函數(shù)值，且效率提升的同時也能保證計算的精度，其中本文方法計算得到的前四階中心矩中峰度結(jié)果的相對誤差最大，但也僅有3.9%。

表5 各求矩方法計算結(jié)果

最后通過最大熵分布擬合功能函數(shù)的分布函數(shù)并進(jìn)行結(jié)構(gòu)可靠度分析，計算得到該空間桁架的結(jié)構(gòu)可靠指標(biāo)為2.934，蒙特卡洛模擬法計算結(jié)果為3.051，相對誤差為3.8%，可知本文方法在效率提升的同時能保證一定的精度。

3 結(jié)論

提出了一種新的矩方法，基于自適應(yīng)二維減維與高階無跡變換進(jìn)行高效求矩，并結(jié)合最大熵原理計算結(jié)構(gòu)的可靠指標(biāo)。通過算例部分的兩道例題，對比了本文方法與既有方法的求矩結(jié)果以及可靠指標(biāo)計算結(jié)果，最終能得到以下的結(jié)論。

(1)求矩方面，所提高效求矩方法在運算過程中，計算結(jié)構(gòu)功能函數(shù)的次數(shù)少于其他求矩方法，因此在計算效率方面優(yōu)于既有的求矩方法。通過不同的求矩結(jié)果與蒙特卡洛模擬方法結(jié)果的對比，可知所提求矩方法在提高效率的同時，在計算精度上也有明顯優(yōu)勢，誤差更小。

(2)可靠指標(biāo)計算方面，通過結(jié)合高效求矩法與最大熵原理計算得到可靠指標(biāo)，并與蒙特卡洛方法結(jié)果對比，可知所提矩方法計算的結(jié)構(gòu)可靠指標(biāo)與蒙特卡洛模擬方法的結(jié)果較為吻合，相對誤差滿足工程要求。