基于二參數尋優設計點的混合結構可靠性分析算法

邱濤， 張建國， 邱繼偉， 魏娟， 游令非

(1.北京航空航天大學 可靠性與系統工程學院, 北京 100191;2.北京航空航天大學 可靠性與環境工程技術國防科技重點實驗室, 北京 100191;3.中國兵器工業標準化研究所, 北京 100089)

0 引言

在機械產品可靠性分析過程中，受加工尺寸、材料性能等條件的限制，往往需要考慮大量的不確定性。對于部分統計信息不足的不確定變量，難以得到對應的精確分布，可以采用區間變量對其進行表述。因此，研究概率-區間共存的可靠性問題具有一定的工程價值[1]。

針對概率-區間混合的可靠性問題，Kang等[2]采用單層優化模型將區間變量嵌入尋優設計點(簡稱MPP)的優化問題中。王騫等[3]通過構建響應面的方法求解最大失效概率。Xie等[4]提出一種單循環算法，將區間優化的Karush-Kuhn-Tucker(KKT)最優條件作為MPP搜索的約束條件。Du[5]提出了基于一次二階矩的統一不確定分析算法。姜潮等[6]考慮極限狀態帶的兩種不同構形，采用了基于響應面的混合可靠性算法。Elishakoff等[7]研究了不確定性的概率模型和凸模型混合問題。Xie等[8]用高維模型表示功能函數，將混合可靠性問題轉化為單步可靠性計算模型。宋向華等[9]利用非概率理論解決了混合可靠性問題。劉瞻等[10]基于Kriging模型和重要抽樣法進行了混合模型的可靠性分析。Du[11]提出了一種序列迭代算法，通過區間分析和概率分析的反復迭代求解，使隨機變量和區間變量同時收斂至最優解。Zhang等[12]利用積分方法求解了概率-區間混合的失效概率。

目前，針對極限狀態函數非線性程度較高的可靠性問題，計算結果的收斂穩定性與效率是該問題的關鍵。通過解耦雙層優化模型，對外層進行概率分析，對內層進行區間分析實現高效序列迭代。雖然一些概率分析算法(如Hasofer-Lind和Rackwitz-Fiessler(HLRF)算法等[13-16])能夠直接進行可靠度計算，但存在結果不收斂的問題。在概率分析時，通過引入兩個調整參數η和λ，分別調整搜索步長和搜尋方向，改進了傳統的尋優MPP算法；在區間分析時，將功能函數進行2次泰勒近似，并轉化為易于求解的2次規劃問題，為解決概率-區間混合的可靠性問題提供了一種新的算法。最后，通過兩個案例驗證了該算法的精度、效率及穩定性。

1 概率- 區間混合可靠性雙層優化模型

機械結構中既存在隨機變量又存在區間變量，則功能函數可以表示為

Z=g(X,Y)，

(1)

式中：X=[X1,X2,…,Xn]T為n維隨機向量；Y=[Y1,Y2,…,Ym]T為m維區間向量。當構件發生故障，即g(X,Y)≤0時，失效概率Pf定義為

Pf=Pr{g(X,Y)≤0}.

(2)

如圖1所示，X1、X2為隨機變量，功能函數因含區間變量Y，g(X,Y)=0在隨機空間中不再是一個曲面，而是由兩個邊界面構成的極限狀態帶。可以通過(3)式和(4)式來定義失效概率的范圍：

(3)

(4)

圖1 極限狀態區域Fig.1 Limit state area

通過將隨機變量轉換為獨立的標準正態隨機變量，上述問題可以轉換為求解如下兩個優化問題，最大和最小可靠度指標分別為

(5)

(6)

式中：U為隨機變量轉換到標準正態空間內的值；G(U,Y)為隨機變量轉換到標準正態空間下功能函數的值；βU為可靠度指標最大值，βL為可靠度指標最小值，對應的最大和最小失效概率分別為

Pfmax=Φ(-βL)，

(7)

Pfmin=Φ(-βU).

(8)

失效概率的范圍為Pfmin≤Pf≤Pfmax，但在實際工程問題中，人們通常關心最大失效概率。因此本文所提算法將以求解最大失效概率Pfmax為例進行混合可靠性分析。

2 混合可靠性分析

針對第1節(6)式中的概率-區間混合可靠性模型，提出一種基于二參數尋優MPP的概率-區間混合可靠性分析算法。如圖2所示，該算法分為概率分析和區間優化兩部分，二者的計算效率共同決定著混合可靠性的計算成本。在隨機變量迭代中，采用基于二參數尋優MPP的概率分析迭代算法；在區間變量迭代中，首先將區間優化問題轉換為2次規劃問題，再結合梯度投影法進行區間優化分析，求解區間最優點。

圖2 序列迭代計算流程Fig.2 Calculation process of sequence iteration

為提高雙層優化模型的計算效率，將雙層優化模型解耦為概率分析和區間分析的序列迭代可靠性模型。一旦滿足收斂條件，則得到最大失效概率點U*，最大失效概率為

(9)

2.1 概率分析

2.1.1 一次二階矩算法和改進的一次二階矩算法

Hasofer等[13]和Rackwitz等[14]提出了HLRF一次二階矩算法，其尋優設計點為

Uk+1=Uk+dk，

(10)

式中：Uk為第k次迭代的MPP；Uk+1為第k+1次迭代的MPP；dk為搜索方向，

(11)

當經過有限次迭代結果收斂后，則收斂點即為MPP. 雖然上述HLRF算法迭代效率較高，但當功能函數非線性程度較高時，此算法不能保證收斂。因此，基于不確定一維搜索準則的改進的一次二階矩(iHLRF)算法[17]被提出來，其迭代搜索設計點為

Uk+1=Uk+λkdk.

(12)

根據Armijo準則：

(13)

(14)

ck為懲罰參數，且ck>0，mk(U)梯度可表示為

(15)

2.1.2 基于二參數尋優MPP的概率分析

設當前迭代步驟為第k+1次，在隨機變量迭代過程中，區間變量保持不變，故在以下尋優MPP算法中省略區間變量。

(16)

圖3 功能函數及其在Uk點處的切平面Fig.3 Performance function and its tangent plane at point Uk

(17)

圖4 在平面G=0上尋優MPPFig.4 Optimizing MPP on plane G=0

(18)

圖5 選定η確定平面G=ηG(Uk)Fig.5 Selecting η to determine the plane G=ηG(Uk)

圖6 在G=ηG(Uk)平面上尋優MPPFig.6 Optimizing MPP on the plane G=ηG(Uk)

(19)

在平面G=ηG(Uk)上，原點連接到輔助點Hk+1的方向即為Uk+1的搜索方向，Uk+1的方向余弦矢量αk+1為

(20)

Uk+1=βk+1αk+1，

(21)

式中：βk+1為第k+1次迭代的可靠度指標。

將(21)式代入(18)式，得

(22)

求解(22)式，得可靠度系數βk+1為

(23)

根據得到的可靠度系數βk+1，可由(21)式得第k+1次迭代的MPPUk+1.

通過引入兩個調整參數分別控制搜索步長和搜索方向，通過參數η調整搜索步長，使功能函數為高度非線性時，尋優MPP收斂成為可能；通過參數λ調整搜索方向，來提高尋優MPP的效率。在選擇λ時，可以將其設定為較大的值，一般取λ最大值為50，在具體的迭代中，可以設定λ=λ/c以自適應更新迭代步長，c為步長的調整系數，是一個大于1的常數，根據數值分析，c最好取在1.1～1.3之間。如果λ的初始值較小，則不需要在整個迭代過程中改變λ.

2.2 區間分析

在第k+1次迭代中，利用概率分析求得Uk+1. 根據(6)式表示的雙層優化模型，固定隨機變量Uk+1，區間分析可表示為(24)式的優化問題：

(24)

式中：Yk+1為第k+1次混合迭代的區間最優點；Yk為第k次混合迭代的區間最優點；YU和YL分別代表區間的上、下界矢量。

為提高計算效率，在進行區間優化之前，首先檢查優化問題的KKT條件，KKT優化條件見(25)式和(26)式。若滿足KKT條件，則Yk+1=Yk；否則進行區間優化分析[18]。

對于(24)式中的優化問題，KKT條件為

(25)

(26)

式中：j=1,2,…,m表示m維區間變量的第j維。

在區間變量循環迭代中，固定隨機變量Uk+1，故在以下區間分析中省略隨機變量。設當前迭代循環次數為s，當前迭代點為Ys，功能函數進行2次泰勒展開，(24)式轉化為二次規劃問題：

(27)

式中：Hs為海森矩陣。

由于(27)式為區間約束的二次規劃問題，當約束條件為區間約束時，將迭代方向投影到可行域上進行計算就相當方便，因此可以采用梯度投影法[19]高效地求得功能函數的極限值。首先將迭代方向投影于可行域，再在投影后的最速下降方向上搜索功能函數的第1個局部最優點，最后利用Schur直接法求解最優點。設最優點為Ys(*)，則區間變量的迭代過程可表示為

Ys+1=Ys+γsYs(*)，

(28)

式中：γs為迭代步長，為保證全局收斂，重復γs←γsρ，直至新迭代點Ys+1滿足Arimijo條件，即

(29)

θ=1×10-4，ρ=0.8.

若新迭代點Ys+1滿足KKT條件，則區間迭代停止，Yk+1=Ys+1；否則繼續區間迭代。

按照以上算法進行迭代，如果滿足‖Uk+1-Uk‖≤ε1且|G(Uk+1,Yk+1)|≤ε2(ε1和ε2是足夠小的非負實數)，則βL=βk+1，計算最壞情況下的失效概率Pfmax：

Pfmax=Φ(-βL).

(30)

綜上所述，本文提出的基于二參數尋優MPP的混合可靠性分析算法流程如下：

1)將隨機變量X轉化為標準正態隨機變量U，得G(U,Y)=0；

2)初始迭代k=0，輸入初始值U0和Y0；

3)固定Yk，選擇η和λ的值(0≤η<1,λ>0)，在G=ηG(Uk)平面上確定Uk+1的搜索方向αk+1，并計算βk+1和Uk+1；

4)檢查KKT條件，如果滿足則Yk+1=Yk，轉步驟6，否則轉步驟5；

5)利用梯度投影法與Arimijo條件，計算Yk+1=Ys+1；

6)檢查收斂，如果滿足‖Uk+1-Uk‖≤ε1且|G(Uk+1,Yk+1)|≤ε2(ε1和ε2是足夠小的非負實數)，則βL=βk+1，轉步驟7，否則k=k+1并轉步驟3；

7)計算最壞情況下的失效概率Pfmax=Φ(-βL)。

3 案例分析

為驗證所提算法的精度、效率及穩定性，本節通過兩個案例進行分析。其中，收斂精度均取ε1=ε2=1×10-6，取λ>2，設定λ=λ/c，c=1.2以自適應更新迭代步長；N1和N2分別為混合可靠性分析循環序列迭代的迭代次數和函數調用次數。

3.1 懸臂梁失效概率計算

如圖7和圖8所示，在懸臂梁末端施加水平力Fx和垂直力Fy，懸臂梁長度為l、橫截面寬度為w、高度為h，l(mm)、w(mm)和h(mm)服從正態分布。懸臂梁的最大應力不超過屈服強度極限值S=380 MPa，Fx和Fy為區間變量，其分布參數如表1所示。

圖7 懸臂梁Fig.7 Cantilever beam

圖8 懸臂梁橫截面Fig.8 Cross section of cantilever beam

表1 懸臂梁的不確定性特征Tab.1 Uncertainty characteristics of cantilever beam

圖9為本文算法與不同混合可靠性分析算法的收斂結果對比。由圖9可知，當功能函數非線性程度較低時，相比于基于HLRF/iHLRF混合可靠性分析算法，本文所提算法同樣具有較高的收斂效率。表2為η、λ取不同值時的分析結果。由表2可知，隨著η的減小、λ的增大，本文所提算法精度基本不變，但收斂效率會逐漸提高。當η在0～1之間取任意值時，λ可取大于0的任意值，且隨著λ的增大，收斂效率會提高；當λ取大于0任意值并固定時，η可取0～1之間的任意值，且隨著η的減小，收斂效率會提高；當η=0、λ→∞時，本文算法與基于HLRF的混合可靠性分析算法具有相同的計算效率，迭代次數均為4，函數調用次數為53.

圖9 本文算法與基于HLRF/iHLRF混合可靠性分析算法的收斂過程對比(η=0.1，λ=0.9)Fig.9 Comparison of the convergence process of the algorithm proposed in this paper and the mixed reliability analysis based on HLRF/iHLRF(η=0.1，λ=0.9)

為驗證本文算法的精度，采取蒙特卡洛法(MCS)計算失效概率。首先將每個區間等分為10份，在區間變量的組合值下對隨機變量抽取1×108次，則功能函數調用1×1010次。表3為懸臂梁最壞情況下的失效概率。由表3可知，本文算法相對誤差為3.9%，與基于HLRF/iHLRF[19]混合可靠性分析算法的函數調用次數為同一數量級，收斂效率基本一致。綜上所述，本文算法通過對兩個參數的調整，既具有較高的計算精度，又具有較高的收斂效率。

表3 懸臂梁最壞情況下的失效概率Tab.3 Failure probabilities cantilever beam in the worst case

3.2 簡支梁失效概率計算

如圖10所示，對簡支梁施加均布載荷q和集中載荷F，梁的長度為l，楊氏模量為E. 簡支梁在C點的撓度不超過極限值δ=6 mm，q和F為區間變量，E、l、w和h服從正態分布，其分布參數特征如表4所示。本例為概率-區間混合的結構可靠性問題，求解簡支梁在C點的撓度不超過極限值的最大可能失效概率。

圖10 簡支梁Fig.10 Simple beam表4 簡支梁的不確定性特征Tab.4 Uncertainty characteristics of simple beam

參數分布均值標準差下界上界l/mm正態81612w/mm正態1107h/mm正態904E/GPa區間200220q/N區間840895F/N區間16801825

表5為不同η、λ值時的失效概率分析結果。由表5可知，當功能函數非線性程度較高時，隨著η取值減小，λ取值增大，分析結果存在不收斂的情況；當η取值趨于0時，本文算法和基于HLRF的混合可靠性分析算法類似，會存在計算收斂不穩定的問題；當η取值趨于1時，即收斂迭代步長減小，能解決計算結果收斂不穩定的問題，使計算結果收斂；當λ取值增大時，隨著沿梯度方向的步長增加，收斂速度變快，但會引起收斂不穩定的問題，此時可以增加步長調整系數c，以自適應調整迭代步長，解決計算收斂不穩定的問題。

表5 簡支梁不同η、λ取值的分析結果Tab.5 Analyzed results of different η and λ values for simple beam

為驗證本文算法的精度，將每個區間等分為10份，采用MCS在區間變量的組合值下對隨機變量抽取1×108次，則功能函數調用1×1011次。根據圖11的不同分析方法收斂結果對比圖及表6分析結果誤差與效率對比可知，與MCS相比，本文所提算法計算結果相對誤差為4.3%. 基于HLRF的混合可靠性分析算法，結果不收斂；基于iHLRF的混合可靠性分析算法，在迭代循環47次、函數調用526次后達到收斂；本文所提算法可在迭代循環50次、函數調用638次后達到收斂。以上對比結果表明，當功能函數非線性程度較高時，通過選定合適的η和λ的值，本文算法能解決收斂不穩定的問題，同時具有較高的計算精度和效率。

圖11 非線性程度較高時本文算法(η=0.6，λ=10，c=1.2)與基于HLRF/iHLRF混合可靠性分析的收斂過程對比Fig.11 Comparison of convergence processes of the proposed algorithm (η=0.6，λ=10，c=1.2) and the mixed reliability analysis based on HLRF/iHLRF when the performance function having a high degree of nonlinearity 表6 簡支梁最壞情況下的失效概率Tab.6 Failure probabilities of simple beam in the worst case

分析算法Pfmax相對誤差/%N1N2MCS0.011710111011基于HLRF混合可靠性分析算法基于iHLRF混合可靠性分析算法0.01224.345526本文算法0.01224.349638

4 結論

本文針對機械結構中既有隨機變量又含區間變量的混合可靠性問題，提出了基于二參數尋優MPP的混合可靠性分析算法。所得結論如下：

1) 針對概率-區間混合的可靠性問題，可將其轉化為雙層優化模型，實現概率分析和區間分析共同作用的高效序列迭代算法，使計算結果具有較高的計算精度和效率。

2) 在概率分析時，通過引入兩個分別控制搜索方向和搜索步長的調整參數，當功能函數非線性程度較高時，同樣能使計算結果收斂，且具有較高的計算效率。