謝少軍 潘柏松 羅路平 項(xiàng)涌涌

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

0 引言

因試驗(yàn)數(shù)據(jù)不足，某些不確定性變量采用區(qū)間模型描述，如計(jì)算誤差、試驗(yàn)數(shù)據(jù)誤差、尺寸公差等，這導(dǎo)致可靠性分析既涉及傳統(tǒng)的用概率分布函數(shù)描述的隨機(jī)變量，也存在用區(qū)間模型描述的區(qū)間變量。區(qū)間變量在隨機(jī)-區(qū)間混合可靠性分析中導(dǎo)致雙層耦合循環(huán)，計(jì)算效率較低，由此，國內(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ū)間變量，則會(huì)增加分析結(jié)果的不確定性，產(chǎn)生保守結(jié)果。因此，考慮區(qū)間變量耦合性，提出耦合區(qū)間-隨機(jī)混合可靠性分析方法具有重要的研究意義。

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

為進(jìn)一步解決混合可靠性分析計(jì)算效率低的問題，本文提出了一種高效的耦合區(qū)間-隨機(jī)混合可靠性分析方法。利用序列迭代分析方法，將耦合區(qū)間-概率混合可靠性分析的雙層耦合優(yōu)化問題分解為概率分析與區(qū)間分析各自依序迭代的兩部分，并結(jié)合多項(xiàng)式插值法提出了一種適用于極限狀態(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時(shí)，系統(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)的積分計(jì)算，獲得失效概率

(3)

因?yàn)槭в虻娘@式表達(dá)式往往難以獲取，并且上述積分涉及多重積分計(jì)算，故失效概率一般很難由式(3)直接求得，通常用近似方法如一次二階矩法[6-9]進(jìn)行求解。通過對極限狀態(tài)函數(shù)在標(biāo)準(zhǔn)正態(tài)分布空間內(nèi)的線性近似，一次二階矩法可高效地求得失效概率的近似值。一次二階矩法主要包括兩個(gè)步驟，首先利用以下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ū)間變量共存時(shí)，極限狀態(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)可見，計(jì)算最小與最大失效概率是一個(gè)雙層耦合循環(huán)問題，外層循環(huán)是以隨機(jī)變量為變量的概率分析，用于計(jì)算失效概率；內(nèi)層循環(huán)是以耦合區(qū)間變量為變量的區(qū)間分析，用于求解耦合區(qū)間變量優(yōu)化點(diǎn)。為提高計(jì)算效率，采用序列迭代分析方法，將雙層耦合循環(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的球模型。

計(jì)算最大失效概率的序列迭代分析流程見圖1。在第k迭代步時(shí)，標(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)上述條件滿足時(shí)，則迭代停止，輸出最大概率點(diǎn)u*=u(k+1)，最大失效概率

Pfmax=Φ(-‖u*‖)

(17)

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

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

(18)

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

(19)

搜尋步長

(20)

本文提出了一種高效的區(qū)間分析算法，為節(jié)省篇幅，僅給出計(jì)算最大失效概率的計(jì)算公式，計(jì)算最小失效概率僅需將式(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ū)間分析時(shí)，迭代步初始化為零，正則耦合區(qū)間變量初始點(diǎn)設(shè)為v(l=0)=v(k)。

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

(21)

若滿足優(yōu)化點(diǎn)位于邊界點(diǎn)的假設(shè)條件，則基于式(21)的迭代式可高效地求得優(yōu)化點(diǎn)。但若優(yōu)化點(diǎn)位于可行域內(nèi)部，則式(21)不能收斂。為此，利用二階多項(xiàng)式插值法，搜索在連接兩個(gè)邊界點(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ù)迭代計(jì)算。

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

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

為驗(yàn)證提出的區(qū)間分析算法的計(jì)算效率，在算例中對本文方法與FORM-UUA方法進(jìn)行比較分析。雖然FORM-UUA方法被提出時(shí)用于處理獨(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)修改，用于驗(yàn)證本文方法的計(jì)算精度與效率。如圖4所示，滑塊受到水平力p作用，當(dāng)連桿最大應(yīng)力超出材料許用應(yīng)力S時(shí)，連桿滑塊機(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方法與本文方法計(jì)算得到的最大失效概率見表2。使用蒙特卡洛法時(shí)，將耦合區(qū)間變量的區(qū)間劃分為50等份，在每種滿足可行域的耦合區(qū)間變量組合下，對隨機(jī)變量作106次抽樣，計(jì)算失效概率，挑選出最大值為最大失效概率。以蒙特卡洛法的計(jì)算結(jié)果為參照，F(xiàn)ORM-UUA方法與本文方法在相同的迭代步下均得到了較準(zhǔn)確的結(jié)果。然而，由調(diào)用極限狀態(tài)函數(shù)的次數(shù)Nc可見，本文方法比FORM-UUA方法效率高。

表2 最大失效概率

為進(jìn)一步分析提出的區(qū)間分析算法的有效性，在迭代步k=0時(shí)，區(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)，計(jì)算效率較高。

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

4.2 懸臂圓筒

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

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

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

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

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與本文方法計(jì)算得到的結(jié)果見表5。由表5可見，F(xiàn)ORM-UUA方法與本文方法的精度均較高，但本文方法比FORM-UUA方法的效率有較大幅度的提高。

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

表5 最大失效概率

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

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

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

5 結(jié)論

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