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

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

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

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

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

1 PC-Kriging模型

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

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

(1)

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

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

(2)

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

(3)

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

(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ù)來增強(qiáng)預(yù)測模型的全局近似精度，并利用Kriging模型來捕捉預(yù)測模型局部特征的能力.采用最小角回歸(LAR)構(gòu)建回歸基函數(shù)的最優(yōu)多項式數(shù)量集，同時用Akaike信息準(zhǔn)則(AIC)來確定最優(yōu)的截斷集合.

(6)

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

一個完備的Hilbert空間L2(R,f)的標(biāo)準(zhǔn)正交基(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]準(zhǔn)則來確定哪一個是最佳的.最后，保留中間的候選項以確保Kriging模型的準(zhǔn)確性，同時大大減少函數(shù)調(diào)用的數(shù)量(即保留對模型貢獻(xiàn)更多的候選項).選擇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īng)Y.

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

步驟6 重復(fù)步驟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 自適應(yīng)PC-Kriging方法

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

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

(9)

2.1 自適應(yīng)PC-Kriging方法

所提出的k-means聚類分析自適應(yīng)策略選取樣本點(diǎn)方法步驟：

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

(12)

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

2.2 收斂條件

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

(13)

(14)

式中，用NU

3 算 例

3.1 算例1

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

其中，x1,x2為服從標(biāo)準(zhǔn)正態(tài)N(0,1)的獨(dú)立同分布隨機(jī)變量.

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

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

3.2 算例2

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

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

表2 隨機(jī)變量的分布特征

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

g(x)=S-σmax.

其中：σmax表示在原點(diǎn)處筒上表面所受的最大等效應(yīng)力；σx為正應(yīng)力；τzx為扭應(yīng)力.

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

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

4 結(jié) 論

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

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

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

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