一種自適應PC-Kriging模型的結(jié)構(gòu)可靠性分析方法

于震梁， 孫志禮， 張毅博， 王 健

(東北大學 機械工程與自動化學院， 遼寧 沈陽 110819)

隨著結(jié)構(gòu)可靠性評估和穩(wěn)健性設(shè)計概念在現(xiàn)代工程可靠性分析中的日益深入，可靠性分析方法也得到了廣泛發(fā)展和進一步研究.由于在實際工程問題中其功能函數(shù)往往是隱式的(強非線性或耗時的)，采用經(jīng)典分析技術(shù)如一階可靠度法(FORM)、二階可靠度法(SORM)[1]、蒙特卡羅模擬法(MCS)在精度和時間成本上往往難以接受.因此，采用代理模型法(如響應面法[2]、人工神經(jīng)網(wǎng)絡[3]、支持向量機[4]、Kriging[5]等)來近似功能函數(shù)計算可靠性得到了迅速發(fā)展.

Kriging代理模型[6]是一種精確的插值方法且具有隨機性，不僅能提供未采樣點的預測值，還能對預測方差進行估計.因此，該方法在結(jié)構(gòu)可靠性分析中得到了廣泛應用.為了進一步提高模型的計算精度和效率，一些學者對傳統(tǒng)的Kriging模型進行了改進研究，如在建立Kriging模型的過程中加入樣本點函數(shù)值的同時也考慮了其函數(shù)值的梯度信息值，這種模型被稱為Co-Kriging模型[7]或梯度增強Kriging(gradient-enhanced Kriging，GEK)[8]，在此基礎(chǔ)上，Han等[9]為進一步提高GEK的計算效率，提出了改進的梯度增強Kriging方法(weighted gradient-enhanced Kriging,WGEK).Schoebi等[10]提出了一種基于多項式混沌展開(polynomial-chaos-expansion，PCE)和Kriging模型相結(jié)合的可靠性方法(polynomial-chaos-based Kriging，PC-Kriging).在分析比較復雜工程結(jié)構(gòu)問題時，往往采用試驗設(shè)計的思想.因此，若干種自適應試驗設(shè)計(DoE)策略已被構(gòu)建，如Bichon等[11]通過所提的EFF函數(shù)選擇距極限狀態(tài)最近的樣本點作為新增訓練點.Echard等[12]提出U函數(shù)來衡量未測點符號預測錯誤的概率，并將U值最小的點定義為最佳樣本點.

上述這些基于自適應DoE策略的學習函數(shù)能夠有效地提高Kriging模型的準確性和效率.然而，由文獻[13]可知具有不同基函數(shù)的Kriging模型的準確性是不同的，使得Kriging模型在高階的可靠性分析中難以計算.為了解決這類問題，提出了一種基于PC-Kriging和自適應k-means相結(jié)合的結(jié)構(gòu)可靠性分析方法(APC-Kriging).首先，PC-Kriging是一種改進的Kriging算法，其回歸基函數(shù)采用稀疏多項式最優(yōu)截斷集合來近似數(shù)值模型全局行為，而用Kriging來處理模型輸出的局部變化，在保證精度的同時提高了計算效率.其次，常見的可靠性方法的采集樣本點為逐個采集，而本文的自適應k-means聚類分析將空間分成若干個區(qū)域，并從每個區(qū)域選取一個最佳樣本點，從而使多個區(qū)域同時達到提高PC-Kriging模型精度的目的，從而再次提升模型的計算效率.

1 PC-Kriging模型

在構(gòu)造結(jié)構(gòu)功能函數(shù)G(x)的近似代理模型時，Kriging將其假定為由確定性g(x)Tβ和隨機性z(x)兩部分組成，即

G(x)=g(x)Tβ+z(x).

(1)

其中：向量x為包含M個對結(jié)構(gòu)狀態(tài)(即功能函數(shù))產(chǎn)生影響的基本隨機變量的集合；g(x)為模型基函數(shù)，即在PC-Kriging模型中改進的部分；β為基函數(shù)的系數(shù)矢量；z(x)為均值為0的高斯隨機過程，其協(xié)方差為

Cov[z(xi),z(xj)]=σ2R(xi,xj;θ).

(2)

式中：σ為z(x)的標準差；R(xi,xj;θ)為z(xi)和z(xj)間的帶有參數(shù)θ的相關(guān)系數(shù)，通常為應用最廣泛的高斯相關(guān)函數(shù).

(3)

給定N個訓練樣本SDoE=[x1,x2,…,xN]及其對應的結(jié)構(gòu)狀態(tài)值Y=[y1,y2,…,yN]，G(x)預測值的無偏估計及其方差為

(4)

(5)

式中：

u(x)=GTR-1r(x)-g(x)，

R=(R(xi,xj;θ))N×N，

G=[g(x1),…,g(xN)]T.

其中，參數(shù)向量θ需采用極大似然估計獲取.

PC-Kriging采用多項式混沌展開替代傳統(tǒng)Kriging模型的回歸基函數(shù)來增強預測模型的全局近似精度，并利用Kriging模型來捕捉預測模型局部特征的能力.采用最小角回歸(LAR)構(gòu)建回歸基函數(shù)的最優(yōu)多項式數(shù)量集，同時用Akaike信息準則(AIC)來確定最優(yōu)的截斷集合.

(6)

式中，fXi(x)為x的第i個邊緣概率密度函數(shù).

一個完備的Hilbert空間L2(R,f)的標準正交基(f表示x的聯(lián)合概率密度函數(shù))是

其中，α=[α1,α2,…,αM]是一個自然數(shù)的M維向量.這里，滿足總項數(shù)|α|=α1+α2+…+αM不超過給定閾值T0的項將被保留作為Kriging模型基函數(shù)的候選項.基函數(shù)的候選項ΑM,T0可表示為

AM,T0={α∈NM,|α|≤T0}，AM,T0中項數(shù)的數(shù)量可由P表示:

(7)

然后，考慮基函數(shù)的所有候選項的“完備”設(shè)計矩陣：

GM,T0=[ψα0,ψα2,…,ψαP-1]，

ψαi=[ψαi(x1),ψαi(x2),…,ψαi(xN)]T.

式中：αi∈A(i=0,1,…,P-1)；xnSDoE(n=1,2,…,N).

根據(jù)式(7)，如果將ΑM,T0中的所有候選函數(shù)用作Kriging模型的基函數(shù)，則函數(shù)調(diào)用的次數(shù)將隨著T0的增加而急劇增加.為了避免這類問題，在構(gòu)造稀疏多項式Kriging模型的基函數(shù)時，采用LAR[14]理論定義了函數(shù)的基函數(shù)集個數(shù)，并用AIC[15]準則來確定哪一個是最佳的.最后，保留中間的候選項以確保Kriging模型的準確性，同時大大減少函數(shù)調(diào)用的數(shù)量(即保留對模型貢獻更多的候選項).選擇Kriging模型基函數(shù)的主要步驟：

步驟1 設(shè)置相關(guān)參數(shù)值，即保留多項式T0的最大階次和基函數(shù)pmax的最多項數(shù)，本文設(shè)T0=3，pmax=0.5card(SDoE).H=min{P,pmax}.

步驟2 初始化所有候選項中的系數(shù)aαi=0 (i=0,1,…,P-1).根據(jù)LAR理論，初始化后的剩余項等于結(jié)構(gòu)響應Y.

步驟5h=h+1.共同移動a朝向當前殘差為Gh-1的聯(lián)合最小二乘系數(shù)，直到向量ψαh與當前殘差Gk-1具有相同的相關(guān)系數(shù).

步驟6 重復步驟5，直到滿足h=H.

步驟7 計算Gh(h=1,…,H)的AIC值：

AICh=Nln(SSEh)+2h.

其中,

SSEh=[Gh(GhTGh)-1GhY-Y]T×
[Gh(GhTGh)-1GhY-Y].

步驟8 找到最小的AICh(h=1,…,H),p=argminh{AICh;h=1,…,H},對于SDoE和Y，PC-Kriging模型的最優(yōu)基函數(shù)為

(8)

2 自適應PC-Kriging方法

構(gòu)造完近似代理模型后，結(jié)構(gòu)的失效概率可通過式(9)近似計算：

從以上幾組中醫(yī)最為常見的英譯解析中可以看出，翻譯的方式方法不是單一的，而是音譯與意譯，直譯與輔譯相結(jié)合的。中醫(yī)英譯是中醫(yī)走向世界的一道重要關(guān)口，是一項嚴謹、艱苦并富有挑戰(zhàn)性的工程，要求譯者既精通英語、又要熟悉傳統(tǒng)的中醫(yī)學理論及西醫(yī)理論，力求譯文客觀、嚴密、準確、簡潔，既符合英語的思維習慣又表達中醫(yī)的豐富思想，在中西方文化鴻溝中建起一座橋梁。

(9)

2.1 自適應PC-Kriging方法

所提出的k-means聚類分析自適應策略選取樣本點方法步驟：

Ω0={(x0,i,y0,i),i=1,2,…,M0}，

X0={x0,1,x0,2,…,x0,M0}.

(10)

(11)

其中，i=1,2,…,k.

步驟4 調(diào)整集合St-1中各點位置.定義距離D0，如式(12)所示.假設(shè)在集合St-1中如果任意兩個樣本點之間的距離小于D0時可視為是不能接受的，此時需要將集合St-1中個別點位置進行調(diào)整.

(12)

式中e為給定常數(shù)．

2.2 收斂條件

收斂條件采用文獻[16]中的學習停止條件，其基本思想為隨著迭代過程進行，符號預測錯誤的樣本點數(shù)占總失效樣本點數(shù)的比重很小時，失效概率估計值滿足精度要求，學習過程停止，其表達式為

(13)

(14)

式中，用NU

3 算 例

3.1 算例1

選取文獻[12]中具有多個設(shè)計點的狀態(tài)函數(shù)，其表達式為

其中，x1,x2為服從標準正態(tài)N(0,1)的獨立同分布隨機變量.

適當抽取構(gòu)建初始PC-Kriging模型時所需樣本點數(shù)，設(shè)N0=6.為取得分布較為均勻的初始樣本點，采用拉丁超立方抽樣(Latin hypercube sampling, LHS)在[-5, 5]區(qū)域內(nèi)抽取初始樣本點.通過Matlab中的工具箱根據(jù)所提算法建立PC-Kriging預測模型，并通過主動學習，更新樣本空間DoE，并不斷提高模型精度，重復該過程，直到滿足迭代停止條件，結(jié)果見圖1、表1.其中，Nit表示迭代次數(shù)，Ncall為調(diào)用樣本點數(shù)，ε為失效概率估計值與標準值的相對誤差，失效概率標準值由MCS方法計算得到.

通過圖1和表1的結(jié)果對比，可以發(fā)現(xiàn)在Ncall項所提方法需要最少的樣本點就可以達到足夠高精度，而通過對Nit項，因為應用到聚類分析，使得在滿足精度的要求下迭代次數(shù)得到大幅度的降低，從而提高計算效率.

3.2 算例2

選用文獻[17]懸臂式圓柱筒結(jié)構(gòu)(圖2)，它是一個擁有9個隨機變量的高維非線性的工程結(jié)構(gòu)實例，該結(jié)構(gòu)的9個輸入隨機變量分別為t，d，L1，L2，F(xiàn)1，F(xiàn)2，P，T，S.其分布特征如表2所示.

表1 二維算例結(jié)果對比

表2 隨機變量的分布特征

圖2中懸臂式圓筒結(jié)構(gòu)受到外力F1，F(xiàn)2，P和扭矩T的作用，其功能函數(shù)表示為屈服強度S和最大應力σmax的差：

g(x)=S-σmax.

其中：σmax表示在原點處筒上表面所受的最大等效應力；σx為正應力；τzx為扭應力.

首先，采用Nataf變換將上述隨機變量映射到標準正態(tài)空間，并在[-5，5]9立方體內(nèi)抽取N0=11個拉丁超立方隨機樣本點，再應用所提算法評估懸臂式圓柱筒的可靠性.所提算法與其他現(xiàn)有算法所得失效概率估計值隨迭代次數(shù)的變化趨勢如圖3所示，其中橫坐標為迭代次數(shù)Nit，縱坐標為失效概率估計值的對數(shù)形式.此外，表3列舉了本文算法與其他方法所得結(jié)果.通過比較，不難發(fā)現(xiàn)本文算法相較于其他方法在效率和精度方面更具優(yōu)勢.

表3 九維算例結(jié)果對比

4 結(jié) 論

1) 本文提出了一種新的改進Kriging基函數(shù)與k-means聚類分析相結(jié)合的可靠性分析方法.使得PC-Kriging方法在基函數(shù)回歸方面具有更大的靈活性.結(jié)合k-means聚類分析方法的自適應策略在每次迭代時添加多個新增樣本點來提高計算效率.

2) 通過與兩個算例中MCS抽樣計算的結(jié)果對比可知，所提出的算法與MCS模擬較為接近，從而驗證了本文所提算法計算精度的正確性.

3)通過對比的可靠性分析結(jié)果可以發(fā)現(xiàn)，所提方法相對于AK-MCS+U和AK-MCS+EFF將迭代次數(shù)降低同時提高了估計精度，說明所提方法具有更快的效率以及更高的精度.

4) 所提方法亦可適用于解決功能函數(shù)為隱式多維非線性問題，這為解決實際工程中的可靠性評估問題提供了重要的參考價值.