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

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

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

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

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

1 PC-Kriging模型

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

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

(1)

其中：向量x為包含M個(gè)對(duì)結(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個(gè)訓(xùn)練樣本SDoE=[x1,x2,…,xN]及其對(duì)應(yīng)的結(jié)構(gòu)狀態(tài)值Y=[y1,y2,…,yN]，G(x)預(yù)測(cè)值的無偏估計(jì)及其方差為

(4)

(5)

式中：

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

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

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

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

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

(6)

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

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

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

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

(7)

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

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)造稀疏多項(xiàng)式Kriging模型的基函數(shù)時(shí)，采用LAR[14]理論定義了函數(shù)的基函數(shù)集個(gè)數(shù)，并用AIC[15]準(zhǔn)則來確定哪一個(gè)是最佳的.最后，保留中間的候選項(xiàng)以確保Kriging模型的準(zhǔn)確性，同時(shí)大大減少函數(shù)調(diào)用的數(shù)量(即保留對(duì)模型貢獻(xiàn)更多的候選項(xiàng)).選擇Kriging模型基函數(shù)的主要步驟：

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

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

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

步驟6 重復(fù)步驟5，直到滿足h=H.

步驟7 計(jì)算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},對(duì)于SDoE和Y，PC-Kriging模型的最優(yōu)基函數(shù)為

(8)

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

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

從以上幾組中醫(yī)最為常見的英譯解析中可以看出，翻譯的方式方法不是單一的，而是音譯與意譯，直譯與輔譯相結(jié)合的。中醫(yī)英譯是中醫(yī)走向世界的一道重要關(guān)口，是一項(xiàng)嚴(yán)謹(jǐn)、艱苦并富有挑戰(zhàn)性的工程，要求譯者既精通英語、又要熟悉傳統(tǒng)的中醫(yī)學(xué)理論及西醫(yī)理論，力求譯文客觀、嚴(yán)密、準(zhǔn)確、簡(jiǎ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中如果任意兩個(gè)樣本點(diǎn)之間的距離小于D0時(shí)可視為是不能接受的，此時(shí)需要將集合St-1中個(gè)別點(diǎn)位置進(jìn)行調(diào)整.

(12)

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

2.2 收斂條件

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

(13)

(14)

式中，用NU

3 算 例

3.1 算例1

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

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

3.2 算例2

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

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

表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個(gè)拉丁超立方隨機(jī)樣本點(diǎn)，再應(yīng)用所提算法評(píng)估懸臂式圓柱筒的可靠性.所提算法與其他現(xiàn)有算法所得失效概率估計(jì)值隨迭代次數(shù)的變化趨勢(shì)如圖3所示，其中橫坐標(biāo)為迭代次數(shù)Nit，縱坐標(biāo)為失效概率估計(jì)值的對(duì)數(shù)形式.此外，表3列舉了本文算法與其他方法所得結(jié)果.通過比較，不難發(fā)現(xiàn)本文算法相較于其他方法在效率和精度方面更具優(yōu)勢(shì).

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

4 結(jié) 論

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

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

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

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