隨機變量和非獨立區間變量下的可靠性序列迭代算法

潘柏松　謝少軍　Du Xiaoping　梁利華

1.浙江工業大學,杭州,3100142.Missouri University of Science and Technology,Rolla,USA

隨機變量和非獨立區間變量下的可靠性序列迭代算法

潘柏松1謝少軍1Du Xiaoping2梁利華1

1.浙江工業大學,杭州,3100142.Missouri University of Science and Technology,Rolla,USA

針對機械系統中輸入變量存在隨機變量和區間變量混合的情況，考慮區間變量的非獨立性，提出高效混合可靠性分析方法。區間變量使可靠性分析問題變為雙層優化問題。為降低雙層優化和非獨立區間變量對可靠性計算效率的影響，提出了高效序列迭代計算策略，基于橢球模型描述的非獨立區間變量，提出將非獨立區間變量轉換為獨立區間變量的方法，并利用二次泰勒近似方法，將區間分析問題轉換為易求解的二次規劃問題。算例結果表明，所提出的可靠性序列迭代算法具有較高的計算效率和精度；可靠性分析結果受區間變量獨立性假設的影響，區間變量獨立可導致較保守的可靠性分析結果。

可靠性分析;橢球模型;非獨立區間變量;序列迭代算法

0　引言

為克服傳統可靠性分析及設計在處理認知不確定性方面的不足，學者們發展了較多非概率建模理論,如可能性理論[3]、證據理論[4]、凸集模型[5],以及概率建模理論——貝葉斯理論[6]。作為凸集模型的特例，區間模型[7]在實數軸上規定了認知不確定性變量可變區間的上下限。在工程應用中,僅有的信息為不確定性變量處于某個區間內的情況十分常見,如結構幾何尺寸、運動副間隙、測量誤差、計算誤差等，基于區間模型的可靠性研究已有諸多報道。如Du等[8]提出了隨機變量和區間變量混合型可靠性設計方法;江濤等[9]基于區間模型提出了非概率可靠性指標的一維優化方法;姜潮等[10]針對分布參數存在區間變量的混合不確定問題，提出了一種時變可靠性分析方法。為提高區間模型計算效率,Du[11]基于一次二階矩法提出了一種高效的混合型可靠性分析方法;Jiang等[12]提出了序列非線性區間規劃方法。

但上述文獻均假設區間變量是相互獨立的,而在實際工程中,某些區間變量存在一定的相關性，是非獨立的，如描述結構幾何尺寸的區間變量和結構質量的區間變量一般存在相關性，較大的幾何尺寸區間變量意味著較大的結構質量區間變量,反之亦然；兩個區間變量可單獨取區間的最大值或最小值，但兩者不同時為最大值或最小值;與幾何尺寸區間變量和結構質量區間變量的相關性相反,某區間變量取值較大表明另一個區間變量較小。在一般常識意義上,本文將這三種情況的相關性分別稱為正相關性、零相關性和負相關性。考慮區間變量的非獨立性具有非常重要的研究意義,但目前關于非獨立區間變量的可靠性研究較少。Du[13]針對機構運動副,基于物理關系式推導獲得非獨立區間變量描述模型——等式與不等式約束條件,提出了一種隨機變量和非獨立區間變量的混合型可靠性設計方法。Jiang等[14]采用多維度平行六面體區間模型，考慮了區間變量為獨立或非獨立的情況,提出了一種新的非線性區間規劃方法，但該規劃方法未考慮系統中同時存在隨機變量和區間變量的混合情況。

本文針對系統輸入變量存在隨機變量和非獨立區間變量混合的情況，基于條件概率法和橢球模型，提出了混合型可靠性分析模型及高效可靠性分析算法。為解決非獨立區間變量對計算效率的影響,利用線性變換,將非獨立區間變量轉換為獨立區間變量；為解決概率分析和區間分析雙層循環計算效率低下難題，提出了序列迭代算法。

1　橢球模型

橢球模型屬于凸集模型,最早由Ben-Haim等[5,15]提出。橢球模型可方便地描述非獨立區間變量,在許多實際應用中,橢球模型較其他模型具有更多的優點[16-17]。

設Y=(Y1,Y2,…,YNY)T為區間變量矢量,區間變量個數為NY。由于在復雜機械系統中,區間變量的維度一般較高,而且區間變量之間的非獨立關系往往不同(如某些區間變量服從正相關關系,而某些區間變量是相互獨立的),因此,需根據不同的獨立性特點,將區間變量歸入不同的組。設經分組后,區間變量矢量可表示為Y=(Y1,Y2,…,YNg)T,其中Ng為組的數量,Yi為第i組區間變量矢量。基于分組后的區間變量,橢球模型可描述為

(1)

因可行域S由Ng個橢球組成,故該模型也稱為多橢球模型,單個橢球模型的變量不超過3個。多橢球模型可方便地描述具有不同獨立特性的區間變量。如當某區間變量是獨立的，則橢球模型可退化為區間模型；當兩個區間變量存在相關性，則橢球模型可退化為橢圓模型。圖1給出了3個區間變量構成的不同幾何形狀的可行域S：在圖1a中，3個變量是相互獨立的；在圖1b中,Y3是獨立的變量,Y1和Y2存在相關性,是非獨立的;在圖1c中,3個變量存在相關性,是非獨立的。

(a)三變量相關獨立(b)一變量獨立，兩變量相關(c)三變量相關圖1　橢球模型

因各個區間變量的單位不同,區間大小不同,不便于分析計算,故將區間變量轉換為量綱一變量:

(2)

則分組后的區間變量矢量Y=(Y1,Y2,…,YNg)T轉換為V=(V1,V2,…,VNg)T。

多橢球模型可相應地表示為

(3)

(4)

(5)

2　可靠性分析模型

設系統極限狀態函數為

G=g(X,Y)

(6)

其中,X=(X1,X2,…,XNX)T為隨機變量矢量,隨機變量的個數為NX;Y=(Y1,Y2,…,YNY)T為區間變量矢量,區間變量的個數為NY。

將區間變量的變換關系式代入式(6),則極限狀態函數可寫為G=g(X,E)。設G<0時系統失效,則系統失效概率Pf可表示為Pf=Pr{g(X,E)<0}。因未知區間變量V的概率分布,故不能獲得準確的失效概率Pf。利用條件概率公式,可得失效概率Pf的最小值Pf min和最大值Pf max的計算公式:

Pf min=Pr{gmax(X,E)<0|E∈S|

(7)

Pf max=Pr{gmin(X,E)<0|E∈S}

(8)

其中,gmax(X,E)和gmin(X,E)分別表示在可行域S內極限狀態函數的全局最大值和最小值。

由式(7)和式(8)可見,當系統極限狀態函數的不確定性輸入變量存在隨機變量和區間變量時，系統失效概率的最小值和最大值分別為最大極限狀態函數和最小極限狀態函數的失效概率。相對于傳統可靠性分析問題，該可靠性分析問題涉及雙層分析循環：內循環為區間分析，在可行域S內搜尋極限狀態函數的極限值；外循環為概率分析，求解最大極限狀態函數或最小極限狀態函數的失效概率。雙層分析循環增加了可靠性分析問題的復雜性，降低了可靠性分析的計算效率，尤其當極限狀態函數由計算機數值仿真模型(如有限元模型、流體動力學模型等)隱式表述時。為提高隨機變量和非獨立區間變量混合情況下的可靠性分析計算效率，本文提出了基于一次二階矩法的高效可靠性分析方法。

3　可靠性分析方法

一次二階矩法是一種針對不確定性輸入變量均為隨機變量的近似可靠性分析方法,它包括兩個關鍵步驟:將隨機變量轉換為獨立的標準正態隨機變量;搜尋最大概率點(MPP),在最大概率點對極限狀態函數作線性近似,最后求得可靠性指標。

因一次二階矩法的計算效率和準確度令人滿意,故在實際工程中已獲大量應用，在此選用一次二階矩法進行可靠性分析。在一次二階矩法中，最大概率點u*的數學優化模型為

(9)

其中,U=(U1,U2,…,UNX)T為獨立的隨機變量矢量,服從標準正態分布,由隨機矢量X經Rosenblatt變換獲得。

e*為經變換后的區間變量最優點,由內層區間分析求得:

(10)

或

(11)

一旦尋得最大概率點,則系統失效概率的最小值和最大值分別為

Pf min=Pr{(gmax(X,E)<0|E∈S}=

Φ(-((u*)Tu*)1/2)

(12)

Pf max=Pr{gmin(X,E)<0|E∈S}=

Φ(-((u*)Tu*)1/2)

(13)

其中,Φ(·)為標準正態分布函數。

圖2給出了失效概率的計算流程,內層區間循環嵌于外層概率分析循環。由圖2可見，概率分析和區間分析的算法效率共同決定了可靠性分析的整體計算效率，為此，本文提出了高效的序列迭代算法，該算法由隨機變量迭代和區間變量迭代兩部分組成。在隨機變量迭代中，采用了高效的iHLRF迭代算法；在區間變量迭代中，將非獨立區間變量轉換為獨立變量，對極限狀態函數作二次近似，將區間分析問題轉換為二次規劃問題，最終利用梯度投影法，求得極限狀態函數的極限值。

圖2　雙層循環計算流程

3.1iHLRF迭代算法[18]

作為HLRF算法的改進算法，iHLRF算法引入了評價函數，用于調整每一步的迭代步長，解決了HLRF在處理高非線性響應函數時收斂困難的問題。因iHLRF算法收斂快速，并具有較好的穩健性,故在工程中被廣泛應用。

設當前迭代步驟為k+1，則最大概率點的迭代公式為

uk+1=uk+αkdk

(14)

其中,αk為迭代步長，dk為迭代方向，其計算式為

(15)

迭代步長αk通過求評價函數最小值獲得，評價函數為

(16)

其中,c為常數,滿足c>‖u‖/‖gu(u,e)‖。

在實際應用中，為減小計算量,在確定αk時無需準確搜尋評價函數最小值,而只需滿足評價函數足夠小條件，迭代步長αk由以下計算式確定：

(17)

在此取b=0.5,c=2(‖u‖/‖gu(u,e)‖)。

3.2區間迭代算法

為提高效率,利用Karush-Kuhn-Tucker最優化條件(KKT條件)事先判斷區間變量迭代初始點是否為優化點。若滿足KKT條件，則跳過區間迭代;若不滿足,則實施區間迭代。因在可靠性分析中工程師往往關心最壞的情況,即最大失效概率，故以下僅具體描述了求解Pf max的區間迭代算法。

基于式(5)給出多橢球模型的參數化表達式,將非獨立區間變量E轉換為相互獨立的區間變量P,轉換關系式表示為P=h(E),具體表達分以下三種情況：

(1)當i個橢球模型中有三個區間變量Ei1、Ei2和Ei3，則令Ei1=Pi1sinPi2cosPi3,Ei2=Pi1sinPi2sinPi3，Ei3=Pi1cosPi2,其中Pi1∈[0,1],Pi2∈[0,π]，Pi3∈[0,2π]；

(2)當i個橢球模型中有兩個區間變量Ei1、Ei2,即橢球模型退化為橢圓模型,則令Ei1=Pi1cosPi2,Ei2=Pi1sinPi2,其中Pi1∈[0,1],Pi2∈[0,2π]；

(3)當橢球模型中只有單個區間變量,即橢球模型退化為區間模型,表明該區間變量是獨立的,則Ei1=Pi1,其中Pi1∈[-1,1]。

將經上述變換后的區間變量代入極限狀態函數g(U,E),g(U,E)可表述為g(U,P),則式(11)可重寫為

(18)

設l為區間迭代循環次數,在迭代初始時,令l=0,初始點pl=0=h(ek)。在區間迭代循環中,隨機變量uk+1保持不變，故在以下給出的區間迭代算法中,省略隨機變量。設當前區間分析循環次數為l,在當前迭代點pl,極限狀態函數作二次泰勒近似,式(18)轉換為二次規劃問題：

(19)

其中,Hl為海森矩陣。因計算Hl需求極限狀態函數的二階偏導數,計算效率較低,故本文采用阻尼BFGS公式近似Hl。

(20)

其中,κl為迭代步長,為保證全局收斂性,利用回代法求解κl,即重復κl←κlρ,直至新迭代點滿足Armijo條件,即

(21)

在此取ρ=0.8,η=1×10-4。

若新迭代點pl+1滿足KKT條件,則區間迭代停止,令ek+1=h-1(pl+1),其中h-1(·)表示區間變量變換關系式的逆變換;否則,令l←l+1,繼續區間迭代。

本文提出的基于一次二階矩法的可靠性序列迭代算法的步驟可整理如下：

(1)輸入初始點u0和e0,初始化迭代次數k=0,初始點均取為零向量；

(2)計算uk+1=uk+αkdk;

(3)判斷ek是否滿足KKT條件，若滿足,則令ek+1←ek,轉步驟(7),若不滿足,則進入下一步,實施區間分析；

(4)初始化區間分析迭代次數l=0,利用轉換關系式,計算pl=h(ek);

(6)判斷pl+1是否滿足KKT條件,若滿足,則令ek+1=h-1(pl+1),進入下一步,若不滿足,則令l←l+1,返回步驟(5);

(7)判斷是否收斂。若|g(uk+1,ek+1)|≤ε1和‖uk+1-uk‖≤ε2(ε1和ε2為非常小的正常數),則迭代停止,令u*=uk+1,否則,令k←k+1,轉步驟(2)。

圖3給出了該算法的流程示意圖。

圖3　求解Pf max的序列迭代流程

4　算例

在MATLAB下,編寫了本文提出的序列迭代算法的可執行程序。為驗證本文提出的序列迭代算法的有效性和計算效率，在本節中給出了兩個混合型可靠性分析算例。盡管兩個算例的極限狀態函數均以顯式表達式給出，但都編寫成了可執行程序,故對于調用函數，極限狀態函數是隱式的。采用前向有限差分法計算極限狀態函數關于隨機變量和區間變量的梯度。

4.1懸臂梁算例

某懸臂梁末端受外部載荷，水平方向分量為Px,垂直方向分量為Py,如圖4所示。考慮兩種失效模式,當梁承受的最大應力超出材料屈服強度σs，則認為強度失效；當梁末端位移大于末端許用位移D0,則認為剛度失效。極限狀態函數分別為

式中,E為材料彈性模量。

圖4　懸臂梁示意圖

已知臂長L=2000 mm，矩形梁截面的寬度B和高度H均為隨機變量,服從正態分布：B～N(55,2)mm,H～N(110,5)mm。材料屈服強度σs=295 MPa,末端許用位移D0=65 mm,材料彈性模量E=210 GPa。載荷分量Px和Py為非獨立區間變量,滿足負相關關系，令Y=(Y1,Y2)T=(Px,Py)T(單位為N),其橢球模型為

Y∈S=

圖5給出了Px和Py滿足負相關和獨立關系時的不同可行域。

(a)負相關關系

(b)獨立關系圖5　Px和Py的可行域

表1給出了兩種失效工況的最大失效概率,并采用了蒙特卡羅法驗證分析結果。為比較計算效率,表1給出了各分析方法調用極限狀態函數的次數Nc。同時,表1給出Px和Py假設為獨立區間變量時,Px∈[4200,4300]N和Py∈[2200,2300]N兩種失效工況的最大失效概率。

由表1可見,本文提出的方法能較高效地求得兩種失效模式的最大失效概率。在蒙特卡羅法中，將每個區間變量的可行區間等分為50份,在區間變量的組合值下取隨機變量的抽樣數為1×106次,則極限狀態函數的調用次數為Nc=2.5×108。基于表1給出的蒙特卡羅法結果和相對誤差百分比可見，提出的方法具較高的精度。由區間變量滿足非獨立和獨立關系時的分析結果可見，區間變量的獨立性對可靠性分析結果的影響較大，假設區間變量獨立會導致較保守的分析結果。

表1　兩種工況的最大失效概率

4.2懸臂圓筒算例

某懸臂圓筒受外部載荷如圖6所示:集中力F1、F2、P和扭矩T。當最大等效von-Mises應力σmax超出材料屈服極限σs,認為懸臂圓筒強度失效,極限狀態函數可寫為

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

圖6　懸臂圓筒

最大等效von-Mises應力位于懸臂圓筒根部截面上端點,其計算式為

其中,σx為該點處的正應力,表達式為

其中,c=d/2,M為該截面處彎矩,A為截面面積,I為截面慣性矩,計算表達式分別為

M=F1L1cosθ1+F2L2cosθ2

τzx為該點的切應力,表達式為

J=2I

表2給出了各隨機參數的分布函數及其參數。角度θ1和θ2為獨立區間變量(單位：(°))，長度L1和L2為非獨立區間變量(單位：mm),滿足零相關性。

表2　隨機變量分布參數

令Y=(Y1,Y2,Y3,Y4)T=(θ1,θ2,L1,L2)T,根據區間變量的獨立性特征，將區間變量分為三組,即Ng=3,則Y=(Y1,Y2，Y3)T=(Y1,Y2,(Y3,Y4)T)T,橢球模型為

圖7給出了L1和L2滿足零相關和獨立關系時的不同可行域。

(a)零相關關系

(b)獨立關系圖7　L1和L2的可行域

表3給出了基于本文提出的方法計算獲得的最大失效概率。由表3可見，本文提出的方法能較高效地求得懸臂圓筒的最大失效概率。為驗證分析結果的正確性，在蒙特卡羅法中，將每個區間變量的可行區間等分為10份,在區間變量的組合值下取隨機變量的抽樣數為1×106，則極限狀態函數的調用次數為Nc=1.0×1010。基于表3給出的蒙特卡羅法結果和相對誤差百分比可見,本文提出的方法具有較高的精度。同樣,由區間變量滿足非獨立和獨立關系時的分析結果可見,區間變量的獨立性對可靠性分析結果的影響較大，假設區間變量獨立會導致較保守的分析結果。

表3　最大失效概率

5　結論

針對機械系統中不確定性輸入變量同時存在隨機變量和區間變量的情況，考慮非獨立性區間變量，基于混合型可靠性分析模型，利用一次二階矩法，提出了一種可靠性序列迭代算法。算例結果表明：該迭代算法的計算效率較高，計算精度較好；不考慮區間變量的非獨立性可產生較保守的可靠性分析結果，可能導致過于保守的可靠性設計結果。

[1]Elishakoff I.Uncertainties in Mechanical Structures:AMF Reudenthal’s Criticisms and Modern Convex Models[J].Journal of Applied Mechanics,1999,63(1):683-692.[2]Nikolaidis E,Ghiocel D M.Singhal S.Engineering Design Reliability Handbook[M].New York:CRC Press,2004.[3]Dubois D,Prade H.A Synthetic View of Belief Revision with Uncertain Inputs in the Framework of Possibility Theory[J].International Journal of Approximate Reasoning,1997,17(2):295-324.

[4]Shafer G.A Mathematical Theory of Evidence[M].Princeton: Princeton University Press,1976.

[5]Yakov B,Elishakoff I.Convex Models of Uncertainty in Applied Mechanics[M].Amsterdam:Elsevier,1990.

[6]Wang P,Youn B D,Xi Z,et al.Bayesian Reliability Analysis with Evolving,Insufficient and Subjective Data Sets[J].Journal of Mechanical Design,2009,131(11):111008.

[7]Moore R E.Methods and Applications of Interval Analysis[M].Philadelphia:Siam,1979.

[8]Du X,Sudjianto A,Huang B.Reliability-based Design with the Mixture of Random and Interval Variables[J].Journal of Mechanical Design,2005,127:1068.

[9]江濤,陳建軍,姜培剛,等.區間模型非概率可靠性指標的一維優化算法[J].工程力學,2007,24(7):23-27.

Jiang Tao,Chen Jianjun,Jiang Peigang,et al.A One-dimensional Optimization Algorithm for Non-probabilstic Reliability Index[J].Engineering Mechanics,2007,24(7):23-27.

[10]姜潮,黃新萍, 韓旭,等.含區間不確定性的結構時變可靠度分析方法[J].機械工程學報,2013,49(10):186-193.

Jiang Chao,Huang Xinping,Han Xu,et al.Time-dependent Structural Reliability Analysis Method with Interval Uncertainty[J].Journal of Mechanical Engineering,2013,49(10):186-193.

[11]Du X.Unified Uncertainty Analysis by the First Order Reliability Method[J].Journal of Mechanical Design,2008,130(9):091401.

[12]Jiang C,Han X,Liu G P.A Sequential Nonlinear Interval Number Programming Method for Uncertain Structures[J].Computer Methods in Applied Mechanics and Engineering,2008,197(49/50):4250-4265.

[13]Du X.Reliability-based Design Optimization with Dependent Interval Variables[J].International Journal for Numerical Methods in Engineering,2012,91(2):218-228.

[14]Jiang C,Zhang Z,Zhang Q,et al.A New Nonlinear Interval Programming Method for Uncertain Problems with Dependent Interval Variables[J].European Journal of Operational Research,2014,238(1):245-253.

[15]Ben-Haim Y.Convex Models of Uncertainty:Applications and Implications[J].Erkenntnis,1994,41(2):139-156.

[16]Kreinovich V,Beck J,Nguyen H T.Ellipsoids and Ellipsoid-shaped Fuzzy Sets as Natural Multi-variate Generalization of Intervals and Fuzzy Numbers: How to Elicit Them from Users, and Howto Use Them in Data Processing[J].Information Sciences,2007,177(2):388-407.

[17]Kreinovich V,Hajagos J G,Tucker W T,et al.Propagating Uncertainty through a Quadratic Response Surface Model[R].Albuquerque,New Mexico:Sandia National Laboratories,2008.

[18]Zhang Y,Kiureghian A.Two Improved Algorithms for Reliability Analysis[C]//Proceedings of the Sixth IFIP WG7.5 Working Conference on Reliability and Optimization of Structural Systems.Springer US,1995:297-304.

[19]Nocedal J,Wright S.Numerical Optimization[M].New York:Springer,2006.

(編輯盧湘帆)

A Sequential Iteration Algorithm for Reliability Analysis with Random and Dependent Interval Variables

Pan Baisong1Xie Shaojun1Du Xiaoping2Liang Lihua1

1.Zhejiang University of Technology,Hangzhou,310014 2.Missouri University of Science and Technology,Rolla,USA

Focusing on the mixed conditions of both random variables and dependent interval variables,an efficient hybrid reliability analysis method was proposed.Interval variables made the reliability analysis problem become a double-loop optimization problem.In order to reduce the impacts on the reliability analysis efficiency by the double-loop optimization and dependent interval variables, an efficient sequentially iterative strategy was proposed,and a formula to transform the dependent interval variables modelled by ellipsoid model was developed into independent interval variables,and the limit-state function was approximated during the inner loop in the second order form to make the optimization problem become a more solvable quadratic programming problem.The results show that the proposed sequential iteration algorithm has good efficiency and accuracy, and that the reliability analysis results are affected by the assumption of dependence between interval variables; results under the assumption of independence between interval variables can be more conservative than that of dependence.

reliability analysis;ellipsoid model;dependent interval variable;sequentially iterative algorithm

2015-01-14

國家自然科學基金資助項目(51475425,51075365)

TH122DOI：10.3969/j.issn.1004-132X.2015.12.002

潘柏松,男,1968年生。浙江工業大學機械工程學院教授、博士研究生導師。主要研究方向為可靠性設計、可靠性工程、現代設計方法。出版專著1部,發表論文40余篇。謝少軍,男,1986年生。浙江工業大學機械工程學院博士研究生。Du Xiaoping，男，1963年生。美國密蘇里科技大學機械和航空航天工程系教授、博士研究生導師。梁利華，男，1973年生。浙江工業大學機械工程學院教授、博士研究生導師。