一種耦合區(qū)間隨機(jī)混合可靠性分析方法

謝少軍 潘柏松 羅路平 項涌涌

浙江工業(yè)大學(xué)機(jī)械工程學(xué)院，杭州，310014

0 引言

因試驗數(shù)據(jù)不足，某些不確定性變量采用區(qū)間模型描述，如計算誤差、試驗數(shù)據(jù)誤差、尺寸公差等，這導(dǎo)致可靠性分析既涉及傳統(tǒng)的用概率分布函數(shù)描述的隨機(jī)變量，也存在用區(qū)間模型描述的區(qū)間變量。區(qū)間變量在隨機(jī)-區(qū)間混合可靠性分析中導(dǎo)致雙層耦合循環(huán)，計算效率較低，由此，國內(nèi)外學(xué)者提出了多種高效的隨機(jī)-區(qū)間混合可靠性分析方法。如DU[1]針對隨機(jī)變量和區(qū)間變量的共存情況，利用序列迭代分析方法將耦合的可靠性問題解耦，并基于一次二階矩法提出一種高效的混合不確定性分析方法(unified uncertainty analysis method based on the first order reliability method，F(xiàn)ORM-UUA)；劉帥杰等[2]針對星載可展天線齒輪防卡問題，開展區(qū)間與概率變量下的齒輪防卡可靠性的研究；劉海波等[3]提出了一種概率-區(qū)間混合不確定性下的串并聯(lián)系統(tǒng)可靠性分析方法。但多數(shù)混合可靠性分析方法均假設(shè)區(qū)間變量是相互獨(dú)立的，區(qū)間變量可取可行域(由其形狀也稱為箱形域)內(nèi)任意值，而在工程實(shí)際中某些區(qū)間變量的可行域是箱形域的子集，即區(qū)間變量是相互耦合的。顯然，如果將耦合區(qū)間變量簡化為獨(dú)立區(qū)間變量，則會增加分析結(jié)果的不確定性，產(chǎn)生保守結(jié)果。因此，考慮區(qū)間變量耦合性，提出耦合區(qū)間-隨機(jī)混合可靠性分析方法具有重要的研究意義。

因區(qū)間變量的耦合性，在耦合區(qū)間-隨機(jī)混合可靠性分析中導(dǎo)致關(guān)于耦合區(qū)間變量優(yōu)化的區(qū)間分析更加復(fù)雜，其中，區(qū)間分析的計算效率是主要技術(shù)難點(diǎn)之一。為提高區(qū)間分析的計算效率，LUO等[4]基于單調(diào)性假設(shè)提出了一種耦合區(qū)間變量和隨機(jī)變量下的高效同步迭代算法，但當(dāng)極限狀態(tài)函數(shù)關(guān)于耦合區(qū)間變量非單調(diào)時，該迭代算法無法收斂。通過極坐標(biāo)轉(zhuǎn)換將耦合區(qū)間變量轉(zhuǎn)換為獨(dú)立區(qū)間變量，潘柏松等[5]基于投影梯度法提出了一種高效區(qū)間分析算法，但該區(qū)間分析算法需計算極限狀態(tài)函數(shù)關(guān)于區(qū)間變量的二階導(dǎo)數(shù)，當(dāng)耦合區(qū)間變量數(shù)量較多時，計算效率可能較低，并且極坐標(biāo)轉(zhuǎn)換較大地增加區(qū)間分析優(yōu)化模型的非線性，可能降低區(qū)間分析算法的精度。

為進(jìn)一步解決混合可靠性分析計算效率低的問題，本文提出了一種高效的耦合區(qū)間-隨機(jī)混合可靠性分析方法。利用序列迭代分析方法，將耦合區(qū)間-概率混合可靠性分析的雙層耦合優(yōu)化問題分解為概率分析與區(qū)間分析各自依序迭代的兩部分，并結(jié)合多項式插值法提出了一種適用于極限狀態(tài)函數(shù)關(guān)于耦合區(qū)間變量單調(diào)與非單調(diào)兩種工況的高效區(qū)間分析算法。

1 傳統(tǒng)可靠性分析方法

設(shè)系統(tǒng)響應(yīng)G關(guān)于隨機(jī)變量的函數(shù)為

G=g(X)

(1)

其中，X表示n維隨機(jī)變量向量。當(dāng)G<0時，系統(tǒng)失效，則失效概率Pf可表示為

Pf=Pr{G=g(X)<0}

(2)

其中，函數(shù)g(·)稱為極限狀態(tài)函數(shù)，Pr{·}表示概率。假設(shè)X的聯(lián)合概率分布函數(shù)為fx(X)，通過作失效域Ω={X:g(X)<0}內(nèi)的積分計算，獲得失效概率

(3)

因為失效域的顯式表達(dá)式往往難以獲取，并且上述積分涉及多重積分計算，故失效概率一般很難由式(3)直接求得，通常用近似方法如一次二階矩法[6-9]進(jìn)行求解。通過對極限狀態(tài)函數(shù)在標(biāo)準(zhǔn)正態(tài)分布空間內(nèi)的線性近似，一次二階矩法可高效地求得失效概率的近似值。一次二階矩法主要包括兩個步驟，首先利用以下Rosenblatt轉(zhuǎn)換，將隨機(jī)變量Xi轉(zhuǎn)換為標(biāo)準(zhǔn)正態(tài)分布變量Ui(用大寫字母表示不確定性變量，小寫字母表示對應(yīng)的不確定性變量取得某一數(shù)值)：

ui=Φ-1[FXi(xi)]

(4)

式中，F(xiàn)Xi(·)為Xi的累積分布函數(shù)；Φ-1[·]為標(biāo)準(zhǔn)正態(tài)分布變量的累積分布函數(shù)的逆函數(shù)。

在標(biāo)準(zhǔn)正態(tài)分布空間內(nèi)，由以下優(yōu)化問題求解最大概率點(diǎn)u*：

(5)

其中，‖·‖表示向量的模。一旦求得u*，則失效概率

Pf=Pr{G=g(X)<0}≈Φ(-β)

(6)

式中，β=‖u*‖為可靠性指標(biāo)。

2 混合可靠性分析方法

當(dāng)隨機(jī)變量與耦合區(qū)間變量共存時，極限狀態(tài)函數(shù)可表示為

G=g(X,Y)

(7)

其中，Y表示由耦合區(qū)間變量構(gòu)成的m維向量，其多橢球模型[10-11]可表示為

(8)

假設(shè)耦合區(qū)間變量為某一定值，即Y=y，則混合可靠性分析變?yōu)閭鹘y(tǒng)的可靠性分析問題，其失效概率可表示為

Pf(Y=y)=Pr{G=g(X,y)<0}

(9)

若耦合區(qū)間變量在其可行域內(nèi)變化，則失效概率是關(guān)于耦合區(qū)間變量的函數(shù)。顯然，可以從中找到最小和最大失效概率，分別為

Pfmin=Pr{Gmax=g(X,y*)<0}

(10)

Pfmax=Pr{Gmin=g(X,y*)<0}

(11)

由式(10)與式(11)可見，計算最小與最大失效概率是一個雙層耦合循環(huán)問題，外層循環(huán)是以隨機(jī)變量為變量的概率分析，用于計算失效概率；內(nèi)層循環(huán)是以耦合區(qū)間變量為變量的區(qū)間分析，用于求解耦合區(qū)間變量優(yōu)化點(diǎn)。為提高計算效率，采用序列迭代分析方法，將雙層耦合循環(huán)問題分解為概率分析與區(qū)間分析兩部分。

在序列迭代分析流程中，概率分析采用前文所述的一次二階矩法。為降低區(qū)間分析的求解難度，首先將耦合區(qū)間變量作如下正則化變換：

(12)

其中，i=1,2,…,Ng，vi表示正則耦合區(qū)間變量，Qi為正交矩陣，其列向量為Wi的標(biāo)準(zhǔn)特征向量，Λi為對角矩陣，其對角元素為對應(yīng)的Wi特征值，滿足WiQi=QiΛi。將式(12)代入式(8)，可得V的可行域V：

(13)

由此可見，V是量綱一的量，各組耦合區(qū)間變量的可行域轉(zhuǎn)換為中心位于坐標(biāo)原點(diǎn)、半徑為1的球模型。

計算最大失效概率的序列迭代分析流程見圖1。在第k迭代步時，標(biāo)準(zhǔn)正態(tài)分布變量點(diǎn)由最大概率點(diǎn)迭代算法更新，新的正則耦合區(qū)間變量由以下優(yōu)化問題求得：

(14)

圖1 序列迭代分析流程圖Fig.1 Flowchart of the sequential analysis procedure

為進(jìn)一步提高效率，在區(qū)間分析前，利用式(14)優(yōu)化問題的Karush-Kuhn-Tucker (KKT)條件，判斷初始值是否為優(yōu)化點(diǎn)，若是則跳過區(qū)間分析，并令v(k+1)=v(k)。KKT條件為

(15)

i=1,2,…,Ng

概率分析與區(qū)間分析的序列迭代收斂條件為

(16)

其中，ε1、ε2為給定的較小正數(shù)。當(dāng)上述條件滿足時，則迭代停止，輸出最大概率點(diǎn)u*=u(k+1)，最大失效概率

Pfmax=Φ(-‖u*‖)

(17)

在概率分析中，最大概率點(diǎn)算法采用GONG等[13]提出的迭代計算公式：

u(k+1)=β(k)α(k)

(18)

其中，α(k)表示第k步的單位搜尋方向，即

(19)

搜尋步長

(20)

本文提出了一種高效的區(qū)間分析算法，為節(jié)省篇幅，僅給出計算最大失效概率的計算公式，計算最小失效概率僅需將式(14)的目標(biāo)函數(shù)g(u(k+1),v)替換為-g(u(k+1),v)。

3 區(qū)間分析算法

區(qū)間分析算法的任務(wù)是高效求解式(14)的優(yōu)化問題，求得正則耦合區(qū)間變量優(yōu)化點(diǎn)v*。令l表示區(qū)間分析算法的迭代步，每次實(shí)施區(qū)間分析時，迭代步初始化為零，正則耦合區(qū)間變量初始點(diǎn)設(shè)為v(l=0)=v(k)。

首先，假設(shè)優(yōu)化點(diǎn)位于邊界點(diǎn)，則可推導(dǎo)獲得以下迭代計算式：

(21)

若滿足優(yōu)化點(diǎn)位于邊界點(diǎn)的假設(shè)條件，則基于式(21)的迭代式可高效地求得優(yōu)化點(diǎn)。但若優(yōu)化點(diǎn)位于可行域內(nèi)部，則式(21)不能收斂。為此，利用二階多項式插值法，搜索在連接兩個邊界點(diǎn)線段內(nèi)部的優(yōu)化點(diǎn)。

判斷迭代前后極限狀態(tài)函數(shù)響應(yīng)值是否減小，若g(u(k+1),v(l+1))>g(u(k+1),v(l))，則表明新的迭代點(diǎn)沒有降低目標(biāo)函數(shù)響應(yīng)值，反而比舊的迭代點(diǎn)效果更差，則舍棄新點(diǎn)，而通過沿舊點(diǎn)的梯度反方向搜索新的邊界點(diǎn)：

(22)

Gvi=gvi(u(k+1),v(l))i=1,2,…,Ng

其中，系數(shù)ai為非負(fù)數(shù)，表示為

(23)

(24)

圖2 迭代點(diǎn)相互重合示意圖Fig.2 Two sequentially iterative points coincide

基于式(22)或式(24)獲得新迭代點(diǎn)后，比較目標(biāo)函數(shù)響應(yīng)值。若g(u(k+1),v(l+1))≤g(u(k+1),v(l))，則判斷是否滿足KKT條件，若滿足，則區(qū)間分析迭代停止，輸出v(k+1)=v(l+1)；否則，令l←l+1，并基于式(22)或式(24)繼續(xù)迭代計算。

若g(u(k+1),v(l+1))>g(u(k+1),v(l))，則建立二階多項式插值函數(shù)，搜索連接v(l+1)與v(l)線段內(nèi)部的極小值。基于線段端點(diǎn)處信息，建立以下二階多項式插值函數(shù)：

gq(t)=a+bt+ct2

(25)

其中，t表示線段內(nèi)的位置參數(shù)，滿足0≤t≤1。令t*為線段內(nèi)局部優(yōu)化點(diǎn)，則滿足

(26)

gq(0)=g(u(k+1),v(l))

g′q(0)=(v(l+1)-v(l))Tgv(u(k+1),v(l))

gq(1)=g(u(k+1),v(l+1))

則對應(yīng)t*的迭代點(diǎn)為

v(l+1)=(v(l+1)-v(l))t*+v(l)

(27)

檢查v(l+1)是否滿足KKT條件，若滿足則迭代停止，輸出v(k+1)=v(l+1)；否則，令l←l+1，并基于式(22)或式(24)繼續(xù)迭代。區(qū)間分析算法的流程圖見圖3。

圖3 區(qū)間分析算法流程圖Fig.3 Flowchart of the proposed interval analysis algorithm

4 算例

本文采用兩個算例來驗證提出的混合可靠性分析方法的精度與計算效率。分析方法中，采用有限差分法計算極限狀態(tài)函數(shù)的梯度，并采用調(diào)用極限狀態(tài)函數(shù)的次數(shù)來評定分析方法的計算效率。

為驗證提出的區(qū)間分析算法的計算效率，在算例中對本文方法與FORM-UUA方法進(jìn)行比較分析。雖然FORM-UUA方法被提出時用于處理獨(dú)立區(qū)間變量，但該方法在區(qū)間分析中采用了常規(guī)非線性優(yōu)化算法序列二次規(guī)劃法，用于求解關(guān)于區(qū)間變量的極值優(yōu)化問題，故該方法同樣適用于耦合區(qū)間變量。

4.1 連桿滑塊機(jī)構(gòu)

對文獻(xiàn)[1]中的連桿滑塊機(jī)構(gòu)算例作適當(dāng)修改，用于驗證本文方法的計算精度與效率。如圖4所示，滑塊受到水平力p作用，當(dāng)連桿最大應(yīng)力超出材料許用應(yīng)力S時，連桿滑塊機(jī)構(gòu)失效，則極限狀態(tài)函數(shù)可表示為

其中，l1為曲柄長度；l2為連桿長度；d1、d2分別為圓筒連桿的內(nèi)外徑；S、p、l1、l2為隨機(jī)變量，其分布參數(shù)見表1。

表1 隨機(jī)變量分布參數(shù)

圖4 連桿滑塊機(jī)構(gòu)Fig.4 A crank-slider mechanism

因機(jī)構(gòu)安裝地點(diǎn)不定，偏移量、滑塊與地面的摩擦因數(shù)被設(shè)為區(qū)間變量，而圓筒連桿的內(nèi)外徑設(shè)為耦合區(qū)間變量，則多橢球模型為

Y=[Y1Y2Y3Y4]T=[eμd1d2]T

Y3=[Y3Y4]T

由蒙特卡洛法、FOMR-UUA方法與本文方法計算得到的最大失效概率見表2。使用蒙特卡洛法時，將耦合區(qū)間變量的區(qū)間劃分為50等份，在每種滿足可行域的耦合區(qū)間變量組合下，對隨機(jī)變量作106次抽樣，計算失效概率，挑選出最大值為最大失效概率。以蒙特卡洛法的計算結(jié)果為參照，F(xiàn)ORM-UUA方法與本文方法在相同的迭代步下均得到了較準(zhǔn)確的結(jié)果。然而，由調(diào)用極限狀態(tài)函數(shù)的次數(shù)Nc可見，本文方法比FORM-UUA方法效率高。

表2 最大失效概率

為進(jìn)一步分析提出的區(qū)間分析算法的有效性，在迭代步k=0時，區(qū)間分析算法的歷史迭代記錄見表3，表3中最后一列表示迭代點(diǎn)是否滿足KKT條件，數(shù)值為0則表示不滿足，為1則表示滿足。可見，正則耦合區(qū)間變量優(yōu)化點(diǎn)位于可行域邊界，滿足優(yōu)化點(diǎn)位于邊界點(diǎn)的假設(shè)，因此，本文區(qū)間分析算法僅需較少的迭代就可快速求得優(yōu)化點(diǎn)，計算效率較高。

表3 區(qū)間分析算法歷史迭代記錄

4.2 懸臂圓筒

通過修改文獻(xiàn)[1]中的懸臂圓筒算例，驗證本文可靠性分析方法。如圖5所示，懸臂圓筒受到外部力F1、F2、p與扭矩T。當(dāng)圓筒最大von-Mises應(yīng)力σmax超出材料屈服強(qiáng)度Sy時懸臂圓筒失效，則極限狀態(tài)函數(shù)為

G=g(X,Y)=Sy-σmax

圖5 懸臂圓筒Fig.5 A cantilever tube

最大von-Mises應(yīng)力σmax位于懸臂圓筒支撐附近，計算公式為

M=F1L1cosθ1+F2L2cosθ2

τzx=Td/(2J)J=2I

式中，L1為圓筒長度；L2為外力F2作用點(diǎn)距懸臂圓筒根部的距離。

隨機(jī)變量的分布函數(shù)參數(shù)見表4。因信息量不足，力的角度θ1與θ2(單位為度)設(shè)為耦合區(qū)間變量，其多橢球模型為

Y=[Y1Y2]T=[θ1θ2]T

由蒙特卡洛法、FORM-UUA與本文方法計算得到的結(jié)果見表5。由表5可見，F(xiàn)ORM-UUA方法與本文方法的精度均較高，但本文方法比FORM-UUA方法的效率有較大幅度的提高。

表4 隨機(jī)變量分布參數(shù)

表5 最大失效概率

同理，在迭代步k=0時，區(qū)間分析算法的歷史迭代記錄見表6。由表6可見，正則耦合區(qū)間變量優(yōu)化點(diǎn)位于可行域內(nèi)部，未滿足優(yōu)化點(diǎn)位于邊界點(diǎn)的假設(shè)，因此，區(qū)間分析算法在迭代步l為3、5時調(diào)用了二階多項式插值法，用于搜索位于可行域內(nèi)部的優(yōu)化點(diǎn)，并在l=5時求得優(yōu)化點(diǎn)。

表6 區(qū)間分析算法歷史迭代記錄

注：*表示由二階多項式插值函數(shù)獲得。

5 結(jié)論

算例結(jié)果表明，相比已有算法，本文耦合區(qū)間-隨機(jī)混合可靠性分析方法具有更高的計算效率，提出的區(qū)間分析算法可明顯提高混合可靠性分析的整體計算效率，并且該方法可適應(yīng)于耦合區(qū)間變量與獨(dú)立區(qū)間變量，應(yīng)用范圍較廣。