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

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

0 引言

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

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

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

1 概率- 區(qū)間混合可靠性雙層優(yōu)化模型

機械結(jié)構(gòu)中既存在隨機變量又存在區(qū)間變量，則功能函數(shù)可以表示為

Z=g(X,Y)，

(1)

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

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

(2)

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

(3)

(4)

圖1 極限狀態(tài)區(qū)域Fig.1 Limit state area

通過將隨機變量轉(zhuǎn)換為獨立的標準正態(tài)隨機變量，上述問題可以轉(zhuǎn)換為求解如下兩個優(yōu)化問題，最大和最小可靠度指標分別為

(5)

(6)

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

Pfmax=Φ(-βL)，

(7)

Pfmin=Φ(-βU).

(8)

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

2 混合可靠性分析

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

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

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

(9)

2.1 概率分析

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

Hasofer等[13]和Rackwitz等[14]提出了HLRF一次二階矩算法，其尋優(yōu)設(shè)計點為

Uk+1=Uk+dk，

(10)

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

(11)

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

Uk+1=Uk+λkdk.

(12)

根據(jù)Armijo準則：

(13)

(14)

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

(15)

2.1.2 基于二參數(shù)尋優(yōu)MPP的概率分析

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

(16)

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

(17)

圖4 在平面G=0上尋優(yōu)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)平面上尋優(yōu)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)式，得可靠度系數(shù)βk+1為

(23)

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

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

2.2 區(qū)間分析

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

(24)

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

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

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

(25)

(26)

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

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

(27)

式中：Hs為海森矩陣。

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

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

(28)

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

(29)

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

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

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

Pfmax=Φ(-βL).

(30)

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

1)將隨機變量X轉(zhuǎn)化為標準正態(tài)隨機變量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，轉(zhuǎn)步驟6，否則轉(zhuǎn)步驟5；

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

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

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

3 案例分析

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

3.1 懸臂梁失效概率計算

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

圖7 懸臂梁Fig.7 Cantilever beam

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

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

圖9為本文算法與不同混合可靠性分析算法的收斂結(jié)果對比。由圖9可知，當功能函數(shù)非線性程度較低時，相比于基于HLRF/iHLRF混合可靠性分析算法，本文所提算法同樣具有較高的收斂效率。表2為η、λ取不同值時的分析結(jié)果。由表2可知，隨著η的減小、λ的增大，本文所提算法精度基本不變，但收斂效率會逐漸提高。當η在0～1之間取任意值時，λ可取大于0的任意值，且隨著λ的增大，收斂效率會提高；當λ取大于0任意值并固定時，η可取0～1之間的任意值，且隨著η的減小，收斂效率會提高；當η=0、λ→∞時，本文算法與基于HLRF的混合可靠性分析算法具有相同的計算效率，迭代次數(shù)均為4，函數(shù)調(diào)用次數(shù)為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)計算失效概率。首先將每個區(qū)間等分為10份，在區(qū)間變量的組合值下對隨機變量抽取1×108次，則功能函數(shù)調(diào)用1×1010次。表3為懸臂梁最壞情況下的失效概率。由表3可知，本文算法相對誤差為3.9%，與基于HLRF/iHLRF[19]混合可靠性分析算法的函數(shù)調(diào)用次數(shù)為同一數(shù)量級，收斂效率基本一致。綜上所述，本文算法通過對兩個參數(shù)的調(diào)整，既具有較高的計算精度，又具有較高的收斂效率。

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

3.2 簡支梁失效概率計算

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

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

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

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

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

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

圖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 結(jié)論

本文針對機械結(jié)構(gòu)中既有隨機變量又含區(qū)間變量的混合可靠性問題，提出了基于二參數(shù)尋優(yōu)MPP的混合可靠性分析算法。所得結(jié)論如下：

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

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