非獨立區間變量和隨機變量下的單步可靠性計算方法

潘柏松　謝少軍　蔣立正

浙江工業大學,杭州,310014

?

非獨立區間變量和隨機變量下的單步可靠性計算方法

潘柏松謝少軍蔣立正

浙江工業大學,杭州,310014

隨機變量和非獨立區間變量往往共存，兩種變量共存不僅導致出現雙層優化問題，而且會降低可靠性的計算效率。為解決雙層優化問題和提高可靠性計算效率，基于橢球模型描述的非獨立區間變量，利用高維模型表示方法(HDMR)解耦隨機變量和非獨立區間變量，轉換雙層優化問題為簡單的單步求解問題，基于提出的采樣方法，利用二次多項式近似HDMR展式，將隱式的單步求解問題轉化為顯式問題，提出了一種混合型單步可靠性計算方法。算例結果表明，所提出的單步可靠性計算方法具有較高的計算效率和精度；該方法僅需少量的極限狀態函數調用次數，即可獲得較高精度的計算結果。

橢球模型;非獨立區間變量;高維模型表示方法;快速可靠性方法

0　引言

在機械系統可靠性設計過程中，知識、試驗條件、時間及經費等因素的限制，使得某些不確定性變量的統計信息不足，這導致不確定性變量的分布類型及分布函數不能被精確給定。這類因信息不足引起的不確定性被稱為認知不確定性。不同于隨機不確定性，認知不確定性可隨信息量或知識的增加而減小甚至消失[1]。因此，在可靠性工程中，認知不確定性和隨機不確定性往往共存。研究表明[2]，認知不確定性對可靠性分析及設計結果的精度存在較大影響。

為定量描述認知不確定性變量，克服認知不確定性引起的可靠性分析及設計結果精度失真，目前已發展了多種不同類型的建模理論，分為概率建模理論和非概率建模理論，其中，概率建模理論主要為貝葉斯理論[3],非概率建模理論包括可能性理論[4]、證據理論[5]和凸集模型[6]。作為凸集模型特例的區間模型[7]是指在實數軸上規定認知不確定性變量可變區間的上下限。在本文中，基于區間模型描述的認知不確定性變量稱為區間變量。在工程應用中，區間變量十分常見。因此，區間變量和隨機變量共存條件下的混合型可靠性研究具有重要的學術價值及工程應用價值。

目前已有較多的處理獨立區間變量的可靠性分析方法[7-11]，但在實際工程中，某些區間變量存在一定的相關性，是非獨立的。例如描述結構幾何尺寸的區間變量和結構質量的區間變量一般存在相關性，較大的幾何尺寸區間變量意味著較大結構質量區間變量，反之亦然。為此，Du[12]針對機構運動副，基于物理關系式推導獲得非獨立區間變量描述模型——等式與不等式約束條件，提出了一種隨機變量和非獨立區間變量的混合型可靠性設計方法。Jiang等[13]采用多維度平行六面體區間模型，考慮了區間變量為獨立或非獨立的情況，提出了一種新的非線性區間規劃方法，但該規劃方法未考慮系統中同時存在隨機變量和區間變量的混合情況。Jiang等[14]引入樣本相關系數，考慮非獨立概率-區間混合情況，利用矩陣變換將非獨立變量轉換為獨立變量，提出了一種雙層迭代算法。姜潮等[15]通過引入相關角的概念定量描述了任意兩個變量之間的相關性，將不同變量之間的相關性在一個統一的框架下度量，并構建了一高效求解方法。但上述的可靠性分析方法仍為雙層優化問題，影響計算效率。故在非獨立區間變量下，提出一種單層可靠性計算方法，提高可靠性分析的計算效率仍是當前可靠性分析方法研究的一大挑戰。

為此，本文針對系統輸入變量存在隨機變量和非獨立區間變量的混合情況，基于條件概率法和橢球模型，建立了混合型可靠性分析模型，提出了一種高效的單步可靠性模型及單步可靠性計算算法。利用高維模型表示方法[16](HDMR)解耦隨機變量和非獨立區間變量，將混合型可靠性分析模型轉換為單步求解的可靠性分析模型；基于提出的采樣方案，利用二次多項式近似HDMR展式，降低極限狀態函數的調用次數，提高計算效率。

1　橢球模型

為描述非獨立區間變量，Ben-Haim等[6，17]提出了橢球模型。設Y=(Y1,Y2,…,YNY)T為區間變量的矢量，其中NY為區間變量的數量。在復雜機械系統中，區間變量的維度一般較高，非獨立關系也不盡相同，如某些區間變量服從正相關關系，某些滿足負相關關系，而某些是相互獨立的。橢球模型根據不同的非獨立關系，將區間變量歸入不同的組。

經分組后，Y可表示為Y=(Y1,Y2,…,YNg)T，其中，Ng為組的數量，Yi為第i組區間變量矢量，則橢球模型為

i=1,2,…,Ng)

(1)

多橢球模型可描述不同非獨立關系的區間變量：如當某區間變量是獨立的，則橢球模型可退化為區間模型；當兩個區間變量存在相關性，則橢球模型可退化為橢圓模型。圖1給出了3個區間變量構成的不同幾何形狀的可行域S：在圖1a中，3個變量是相互獨立的；在圖1b中，Y3是獨立的變量，Y1和Y2存在相關性，是非獨立的；在圖1c中，3個變量存在相關性，是非獨立的。

(a)Y1、Y2、Y3相互獨立　　　(b)Y3獨立，Y1、Y2相關(c)Y1、Y2、Y3相關圖1　橢球模型

由于各個區間變量的單位不同，區間大小不同，不利于數值計算的穩定性，因此將區間變量Yi轉換為量綱一變量Vi，轉換關系式為

(2)

則分組后的區間變量矢量Y=(Y1,Y2,…,YNg)T轉換為V=(V1,V2,…,VNg)T，多橢球模型相應地表示為

(3)

式(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}，其中Pr{·}表示概率。因未知區間變量E的概率分布，不能獲得準確的失效概率。利用條件概率公式，可得失效概率的最小值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內搜尋極限狀態函數的極限值；外循環為概率分析，求解最大極限狀態函數或最小極限狀態函數的失效概率。

雙層循環增加了可靠性分析問題的復雜性，會大幅度增加極限狀態函數的調用次數。對于復雜機械系統，極限狀態函數一般由計算機數值仿真模型(如有限元模型、流體動力學模型等)隱式表述，極限狀態函數調用次數的增加,會導致計算效率低下，增大可靠性分析及設計的難度。

為提高隨機變量和非獨立區間變量混合情況下的可靠性計算效率，本文采用高維模型表示方法，將雙層循環問題轉化為單步問題，提出了一種單步快速的可靠性計算方法。

2.1高維模型表示方法(HDMR)

高維模型表示方法是一種用于模型近似的處理方法，它常用于近似高維度輸入系統的響應函數。研究表明，低階的高維模型表示方法可精確描述大部分工程中的極限狀態函數[18]。

本節主要介紹如何采用HDMR描述極限狀態函數g(X,E)。設Z=(X,E)T，則g(X,E)可寫為g(Z)?；诟呔S模型表示方法，g(Z)可表示為

(9)

其中，NZ為Z向量的元素個數，NZ=NX+NY；g0為0階分量函數，為常量；gi(Zi)為1階分量函數，表示輸入變量Zi單獨作用時對輸出響應g(Z)的影響；gij(Zi,Zj)為2階分量函數，表示輸入變量Zi和Zj共同作用時對輸出響應g(Z)的影響；更高階的分量函數表示多個輸入變量共同作用時對輸出響應g(Z)的影響；最后一項g12…NZ(Z1,Z2,…,ZNZ)表示所有殘余的耦合輸入變量對輸出響應的影響。

切割法(cut-HDMR)是一種確定高維模型各個分量函數的常用方法。在切割法中，首先選定參考點c，再計算經過參考點c的線、平面、體積等切割幾何上的響應值，分別確定各個分量函數。實際使用中，參考點c一般選定為輸入變量可行空間內最感興趣的點。

利用切割法，各分量函數可表示為

g0=g(c)

(10)

gi(Zi)=g(Zi,ci)-g0

(11)

gij(Zi,Zj)=g(Zi,Zj,cij)-gi(Zi)-gj(Zj)-g0

(12)

其中，g(Zi,ci)=g(c1,c2,…,ci-1,Zi,ci+1,…,cl) ，表示除了分量Zi，其余所有輸入變量均固定在參考點c處，它是一個一元函數；類似地，g(Zi,Zj,cij)為二元函數；最后項g12…l(Z1,Z2,…,Zl)由真實響應值和基于高維模型表示方法的預測值的殘差確定。

在高維模型表示方法中，1階、2階、3階等分量函數與泰勒級數展開式中相應的分量函數有著本質區別。經證明，高維模型表示方法中的1階分量函數gi(Zi)是泰勒級數展開式中僅含有變量Zi分量函數的集合；類似地，2階分量函數gij(Zi,Zj)是泰勒級數展開式中僅含有變量Zi和Zj分量函數的集合。1階分量函數gi(Zi)可以是非線性的。因此，較截斷的泰勒級數展式，任意相應截斷的高維模型表示方法的展式具有較高的精度。

2.2單步可靠性計算模型

若極限狀態函數可描述為關于隨機變量和區間變量相互分離的兩部分函數，則雙層循環問題可變為單層問題，提高計算效率。

利用1階HDMR展式，極限狀態函數可近似為

(13)

如前所述，1階HDMR展式gi(Xi)和gi(Ei)分別是泰勒展式內僅含有變量Xi和Yi所有分量函數的集合,因此，相對于1階泰勒展式，1階HDMR展式沒有限定近似表達式的非線性。這提高了1階HDMR展式的近似精度。

基于1階HDMR展式的近似表達式，失效概率的最小值和最大值的計算模型可寫為

pf，min=Pr{Gmax=(1-NX-NY)g0+

(14)

pf,max=Pr{Gmin=(1-NX-NY)g0+

(15)

(16)

(17)

3　計算失效概率

使用單步可靠性計算模型計算失效概率上下限時，需確定參考點c和各個一元函數表達式。參考點c對可靠性計算結果的精度具有一定的影響。為提高精度，參考點c一般選定為輸入變量可行空間內最感興趣的點。區間變量的可行區間往往較小，區間變量可行區間內的中點可兼顧兩個邊界點，為此，將區間變量可行區間內的中點及區間變量固定在中點時的最大概率點(MPP)對應的隨機變量值設定為參考點c。

最大概率點u*的數學計算模型為

(18)

s.t.g(U,0)=0

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

一旦求得MPP，則參考點c為

(19)

基于求得的參考點c，使用二次多項式近似表達各個一元函數表達式。設gi(Xi,ci)和gi(Ei,cn+i)的近似式分別為

(20)

(21)

其中，ai0、ai1、ai2、bi0、bi1和bi2分別為二次多項式的待定系數。

采用最小二乘法求解各個二次多項式的待定系數。在最小二乘法中，對于隨機變量Xi，沿過參考點c的Xi軸，在[μi-3σi,μi+3σi]區間內均勻分布k(k=5，7，9)個樣本點，如圖2所示，其中μi和σi分別為隨機變量Xi的均值和標準差。各個樣本點的坐標值為μi-3σi，μi-3σi+6σi/(k-1)，…，μi，μi+6σi/(k-1)，…，μi+3σi。對于變換后的區間變量Ei，沿過參考點c的Ei軸，在[-1,1]區間內均勻分布k(k=5，7，9)個樣本點，如圖2所示。各個樣本點的坐標值為-1，-1+2/(k-1)， …，0，2/(k-1)，…，1。

(a)隨機變量的樣本點分布

(b)區間變量的樣本點分布圖2　樣本點分布示意圖

將式(20)和式(21)代入式(14)～式(17)，可得失效概率的最小值和最大值的計算式分別為

pf，min=Pr{Gmax=(1-NX-NY)g(c)+

(22)

pf，max=Pr{Gmin=(1-NX-NY)g(c)+

(23)

(24)

(25)

一旦最大概率點(MPP)和各個二次多項式系數確定，基于HDMR的混合型可靠性計算方法就無須再調用原始的狀態極限函數。這可大幅度提高計算效率，尤其當極限狀態函數以隱式的計算機仿真模型表述時，調用一次狀態極限函數的用時一般較長。為高效求得最大概率點，已有學者提出了較多的數值算法，如HLRF法[19-20]、iHLRF法[21]等。作為HLRF法的改進，iHLRF法引入了價值函數，在處理非線性度較高的極限狀態函數時仍具較好的收斂性，故被廣泛應用。為此，采用iHLRF法求解最大概率點。

4　算例

在MATLAB下，編寫了提出的單步可靠性算法可執行程序，計算了兩個混合型可靠性分析算例，其中第一個算例的非線性較低，輸入變量的維度較低，相比于第一個算例，第二個算例的非線性度較高，輸入變量的維度也較高，用于驗證提出的單步可靠性算法在處理較高非線性和高維度時的計算效率及精度。在算例中，使用調用原始極限狀態函數的次數Nc評定計算效率。盡管兩個算例的極限狀態函數均以顯式表達式給出，但都編寫成了可執行程序，故對于調用函數，極限狀態函數是隱式的。在用于求解最大概率的iHLRF算法中采用向前有限差分法計算極限狀態函數關于隨機變量的梯度。

4.1懸臂梁

某懸臂梁末端受外部載荷，其中，水平方向分量為Px，垂直方向分量為Py，如圖3所示。當梁末端位移大于末端許用位移D0時，認為剛度失效，則極限狀態函數為

其中，L為懸臂梁長度；b和h分別為矩形梁截面的寬度和高度；E為材料彈性模量。

圖3　懸臂梁

已知，末端許用位移D0=65 mm。表1給出了各個隨機變量的分布參數及區間變量的特征矩陣，其中，L、b和h均服從正態分布；彈性模量E為獨立區間變量，載荷分量Px和Py為非獨立區間變量。

表1　不確定性變量的參數

為研究不同樣本數量對計算結果精度的影響，在使用提出的可靠性計算方法計算失效概率時，令k值分別為5、7和9。表2給出了不同k值時的失效概率結果。為驗證計算結果的精度，同時使用了蒙特卡羅法計算失效概率。因蒙特卡羅法僅能計算系統中不確定性都為隨機變量的工況，故在使用蒙特卡羅法時，將各個區間變量在可行域內均布取樣30個點，對滿足多橢球模型約束條件的區間樣本點，調用107次原始極限狀態函數，計算失效概率，最后挑選出失效概率的最小值和最大值，作為失效概率的上下限。由表2可見，當k值等于9時，計算結果與蒙特卡羅法獲得的結果較接近，具有較高的計算精度；根據Nc可知，提出的基于HDMR的單步可靠性計算方法可較少地調用原始極限狀態函數，求得較高精度的失效概率上下限值。

表2　失效概率上下限

4.2懸臂圓筒

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

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

圖4　懸臂圓筒

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

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

c=d/2

M=F1L1cosθ1+F2L2cosθ2

其中，M為該截面處彎矩，A為截面面積，I為截面慣性矩，τzx為該點的切應力。

表3給出了各個隨機變量的分布參數及區間變量的特征矩陣，其中θ1和θ2為獨立區間變量，其余不確定性變量均為隨機變量。

表3　隨機變量分布參數

比較研究了k值對計算結果精度的影響，表4給出了基于提出的方法，令k值分別為5、7和9時的計算結果。為驗證計算結果的正確性，在蒙特卡洛法中，將兩個區間變量在可行域內均布取樣50個點，對滿足多橢球模型約束條件的區間樣本點，調用106次原始極限狀態函數，計算失效概率，最后挑選出失效概率的最小值和最大值，作為失效概率的上下限。由表4可見，k值對該算例的計算結果幾乎沒有影響，并在k值較小時，失效概率上下限已接近基于蒙特卡羅法獲得的值；根據Nc可得提出的基于HDMR的單步可靠性計算方法計算效率高。

表4　失效概率上下限

5　結語

針對機械系統中隨機變量和非獨立區間變量共存的常見工況，基于橢球模型，利用HDMR法，提出了單步可靠性計算模型；使用多項式近似，提出了一種快速可靠性計算算法。由算例結果表明：該算法僅利用少量的原始極限狀態函數的響應信息，或較少的調用次數，即可快速地計算獲得較高精度的失效概率上下限。

在處理極限狀態方程關于輸入變量在可行區間內高度非線性情況時，基于二次多項式函數近似的高維模型的一階分量函數可能會存在較大的誤差，影響計算精度，提出多項式函數階數自適應極限狀態方程非線性的近似方法是一種可行的改進方法。

[1]NikolaidisE,GhiocelDM,SinghalS.EngineeringDesignReliabilityHandbook[M].NewYork:CRCPress, 2004.

[2]ElishakoffI.UncertaintiesinMechanicalStructures:AMFReudenthal’sCriticismsandModernConvexModels[J].JournalofAppliedMechanics, 1999, 63(1):683-692.

[3]WangP,YounBD,XiZ,etal.BayesianReliabilityAnalysiswithEvolving,Insufficient,andSubjectiveDataSets[J].JournalofMechanicalDesign, 2009, 131(11):111008.

[4]DuboisD,PradeH.ASyntheticViewofBeliefRevisionwithUncertainInputsintheFrameworkofPossibilityTheory[J].InternationalJournalofApproximateReasoning, 1997, 17(2):295-324.

[5]ShaferG.AMathematicalTheoryofEvidence[M].Princeton:PrincetonUniversityPress, 1976.

[6]Ben-HaimY,ElishakoffI.ConvexModelsofUncertaintyinAppliedMechanics[M].Amsterdam:Elsevier, 1990.

[7]DuX,SudjiantoA,HuangB.Reliability-basedDesignwiththeMixtureofRandomandIntervalVariables[J].JournalofMechanicalDesign, 2005, 127:1068-1076.

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

JiangTao,ChenJianjun,JiangPeigang,etal.AOne-dimensionalOptimizationAlgorithmforNon-probabilisticReliabilityIndex[J].EngineeringMechanics, 2007, 24(7):23-27.

[9]DuX.UnifiedUncertaintyAnalysisbytheFirstOrderReliabilityMethod[J].JournalofMechanicalDesign, 2008, 130(9):091401.

[10]JiangC,HanX,LiuGP.ASequentialNonlinearIntervalNumberProgrammingMethodforUncertainStructures[J].ComputerMethodsinAppliedMechanicsandEngineering, 2008, 197(49/50):4250-4265.

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

[12]DuX.Reliability-basedDesignOptimizationwithDependentIntervalVariables[J].InternationalJournalforNumericalMethodsinEngineering, 2012, 91(2):218-228.

[13]JiangC,ZhangZG,ZhangQF,etal.ANewNonlinearIntervalProgrammingMethodforUncertainProblemswithDependentIntervalVariables[J].EuropeanJournalofOperationalResearch, 2014, 238(1):245-253.

[14]JiangC,ZhengJ,NiBY,etal.AProbabilisticandIntervalHybridReliabilityAnalysisMethodforStructureswithCorrelatedUncertainParameters[J].InternationalJournalofComputationalMethods, 2015, 12(4)：No.4.

[15]姜潮, 鄭靜, 韓旭,等. 一種考慮相關性的概率-區間混合不確定模型及結構可靠性分析[J]. 力學學報, 2014, 46(4):591-600.JiangChao,ZhengJing,HanXu,etal.AProbabilityandIntervalHybridStructuralReliabilityAnalysisMethodConsideringParameters’Correlation[J].EngineeringMechanics, 2014, 46(4):591-600.

[16]RabitzH,AliF,ShorterJ,etal.EfficientInput-outputModelRepresentations[J].ComputerPhysicsCommunications, 1999, 117(1/2):11-20.

[17]Ben-HaimY.ConvexModelsofUncertainty:ApplicationsandImplications[J].Erkenntnis, 1994, 41(2):139-156.

[18]LiG,RosenthalC,RabitzH.HighDimensionalModelRepresentations[J].TheJournalofPhysicalChemistryA, 2001, 105(33):7765-7777.

[19]HasoferAM,LindNC.ExactandInvariantSecond-momentCodeFormat[J].JournaloftheEngineeringMechanicsDivision, 1974, 100(1):111-121.

[20]RackwitzR,FlesslerB.StructuralReliabilityUnderCombinedRandomLoadSequences[J].Computers&Structures, 1978, 9(5):489-494.

[21]ZhangY,KiureghianA.TwoImprovedAlgorithmsforReliabilityAnalysis[C]//Proceedingsofthe6thIFIPWG7.5WorkingConferenceonReliabilityandOptimizationofStructuralSystems.US:Springer, 1995:297-304.

(編輯王艷麗)

An One-step Reliability Analysis Method with Random and Dependent Interval Variables

Pan BaisongXie ShaojunJiang Lizheng

Zhejiang University of Technology,Hangzhou,310014

Random variables and dependent interval variables often coexisted, and this made the reliability analysis problem into a double-loop optimization problem and reduces the efficiency of the reliability analysis. So an one-step reliability analysis method was proposed. Specifically, the double-loop optimization problem where the dependent interval variables were modelled by the ellipsoid model, was decoupled into a simple one-step problem by using the HDMR. Based on the proposed sampling strategy, a quadratic polynomial was applied to approximate each of the HDMR expression, making the implicit one-step problem into the explicit one. The example results show that the proposed method has good efficiency and accuracy, and may calculate the reliability results with good accuracy at a little cost of calling the origin limit-state function.

ellipsoid model; dependent interval variable; high dimensional model representation(HDMR); rapid reliability analysis method

2015-09-21

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

TH122

10.3969/j.issn.1004-132X.2016.18.003

潘柏松(通信作者)，男，1968年生。浙江工業大學機械工程學院教授、博士研究生導師。主要研究方向為可靠性設計、可靠性工程、現代設計方法。出版專著1部，發表論文40余篇。謝少軍，男，1986年生。浙江工業大學機械工程學院博士研究生。蔣立正，男，1979年生。浙江工業大學機械工程學院博士研究生。