基于Kriging模型的雙變量降維方法

韓忠皓，張德權(quán)，楊美德，趙海文

（1.河北工業(yè)大學(xué)機(jī)械工程學(xué)院，天津 300401；2.湖南大學(xué)機(jī)械與運載工程學(xué)院，湖南 長沙 410082）

0 引言

功能函數(shù)的統(tǒng)計矩是機(jī)械系統(tǒng)進(jìn)行可靠性分析的前提和基礎(chǔ)，其計算效率決定著機(jī)械系統(tǒng)可靠性分析的效率。隨著機(jī)械系統(tǒng)愈加復(fù)雜，高效地獲得功能函數(shù)統(tǒng)計矩變得愈加重要，已成為近年來的研究熱點。

當(dāng)前，功能函數(shù)統(tǒng)計矩的計算方法主要包括三點估計方法[1-2]、稀疏網(wǎng)格法[3]和降維法[4-6]等。對于三點估計方法，Seo等[1]提出對于n維不確定變量問題，需要調(diào)用3n次功能函數(shù)。相關(guān)研究表明[7]：當(dāng)涉及到高維問題時，三點估計方法的功能函數(shù)調(diào)用次數(shù)會隨著不確定變量個數(shù)的增加呈指數(shù)增長，其計算效率非常低。針對三點估計方法的缺點，Smolyak[3]提出一種稀疏網(wǎng)格法，該方法具有較高的計算精度，有效地解決了高維和變量強(qiáng)相關(guān)的不確定分析問題，但仍需調(diào)用大量的功能函數(shù)。Rahman等[5]提出了降維方法，其主要包括單變量降維法[5-6]和雙變量降維法[4]，單變量降維法將一個多維問題轉(zhuǎn)換為多個一維問題進(jìn)行分析，具有較高的計算效率，但該方法忽略了二維及以上積分的殘差，在遭遇強(qiáng)非線性問題時獲得的高階統(tǒng)計矩誤差較大。相比而言，雙變量降維法在計算精度方面比單變量降維法更具有優(yōu)勢，但由于需要調(diào)用更多的功能函數(shù)導(dǎo)致其計算效率較低。

此外，考慮到代理模型[8-9]可以在保證近似精度的基礎(chǔ)上大大提高計算效率，一些研究者將代理模型應(yīng)用于功能函數(shù)的統(tǒng)計矩計算中應(yīng)對復(fù)雜問題，這類方法可以顯著改善計算效率，但其精度取決于樣本點的選取和模型參數(shù)的選取。常見的代理模型有Kriging模型[9-10]、徑向基模型[8，11]、支持向量機(jī)模型[12]、響應(yīng)面模型[13]、人工神經(jīng)網(wǎng)絡(luò)模型[14]、混沌多項式展開模型[15]等，其中Kriging 模型因其可以預(yù)測局部方差的優(yōu)點而備受關(guān)注。Won等[16]將Kriging模型與單變量降維方法相結(jié)合，提出了基于Kriging的單變量降維方法。盡管該方法改善了單變量降維方法的計算效率，但并沒有提高單變量降維方法的計算精度。范文亮等[17]將Kriging模型與雙變量降維方法相結(jié)合，提出了基于雙變量降維和Kriging近似的統(tǒng)計矩評估方法，該方法采用“米”字型節(jié)點構(gòu)建功能函數(shù)的Kriging模型，與原雙變量降維方法相比，顯著提高了統(tǒng)計矩計算效率，但該方法僅僅是簡單的使用Kriging模型代替真實功能函數(shù)。對于一般工程問題來說，該方法構(gòu)建的Kriging模型精度是足夠的，但對于復(fù)雜工程問題來說，無法判斷“米”字型節(jié)點構(gòu)建的Kriging模型是否準(zhǔn)確，可能導(dǎo)致獲得不精確甚至錯誤的統(tǒng)計矩結(jié)果。

鑒于此，本文將U學(xué)習(xí)函數(shù)和Kriging模型引入到雙變量降維方法中，提出了一種新的基于Kriging模型的雙變量降維方法。不同于僅僅采用“米”字型節(jié)點構(gòu)建功能函數(shù)的Kriging模型，提出的方法在使用Kriging模型代替真實的功能函數(shù)的基礎(chǔ)上，進(jìn)一步采用U學(xué)習(xí)函數(shù)逐步地挑選對響應(yīng)函數(shù)值影響最大的高斯積分點，減少使用對計算精度影響較小的積分點。此外，通過采用停止準(zhǔn)則，提出的方法能夠在保證精度的前提下，減少功能函數(shù)調(diào)用次數(shù)，有效克服傳統(tǒng)雙變量降維法計算效率低的問題。

1 雙變量降維方法

為提高計算系統(tǒng)響應(yīng)統(tǒng)計矩的精度，Xu 等[18]在單變量降維方法的基礎(chǔ)上進(jìn)一步提出了雙變量降維方法。與單變量降維方法類似，雙變量降維方法的實施包括3個步驟。

1）將一個多維響應(yīng)函數(shù)等效為多個二維響應(yīng)函數(shù)與多個一維響應(yīng)函數(shù)的疊加，其數(shù)學(xué)形式如下：

式中：n表示隨機(jī)變量的維數(shù)；g(·)表示系統(tǒng)響應(yīng)的功能函數(shù)；和表示二維響應(yīng)函數(shù)的第i1和i2個隨機(jī)變量；xi表示一維響應(yīng)函數(shù)的第i個隨機(jī)變量；μi(i=1,…,n)表示隨機(jī)變量的均值。

2） 根據(jù)所有二維響應(yīng)函數(shù)與一維響應(yīng)函數(shù)的統(tǒng)計矩近似計算多維響應(yīng)函數(shù)的統(tǒng)計矩。令，根據(jù)矩的定義，其l階近似響應(yīng)函數(shù)原點矩可表示為

3）采用基于矩的求積法則[18]計算式（2）中的積分。

2 Kriging 模型和U 學(xué)習(xí)函數(shù)

2.1 Kriging 模型

Kriging模型由線性回歸項和非參數(shù)項2部分組成[19]，具體表達(dá)式為

式中：θ為相關(guān)參數(shù)；Rθ(Xi,Xj)為兩點的相關(guān)函數(shù)，其表達(dá)式為

式中：xi,q表示樣本點Xi的第q個響應(yīng)；θq為第q個相關(guān)參數(shù)。這些參數(shù)可以通過最大似然估計法[19]獲得：

式中，R為對稱相關(guān)矩陣，Rij=Rθ(Xi,Xj)，i,j=1,2,…,n。采用廣義最小二乘法[19]計算回歸系數(shù)β和σ2的估計值

式中：F是包含m1×m2個元素的矩陣，。

給定一個預(yù)測點X0，在該點預(yù)測的函數(shù)值和方差為

式中，u=FTR-1r0-f(X0)。在本研究中，采用DACE工具箱[20]建立Kriging模型并計算相應(yīng)的預(yù)測值。

2.2 U 學(xué)習(xí)函數(shù)

在保證Kriging 模型具有足夠精度的前提下，盡可能提高計算效率，Echard 等[21]提出了一種U 學(xué)習(xí)函數(shù)，并采用該學(xué)習(xí)函數(shù)尋找潛在最優(yōu)樣本點更新Kriging模型，逐步改善Kriging模型的精度。U學(xué)習(xí)函數(shù)[21]的數(shù)學(xué)表達(dá)式為

3 提出的方法

3.1 基于Gauss-Hermite 的統(tǒng)計矩評估

基于矩的求積法則在解決線性方程組時可能遭遇數(shù)值不穩(wěn)定現(xiàn)象[22-23]，受Huang 等[4]研究工作的啟發(fā)，本研究采用Gauss-Hermite數(shù)值積分求解式（2）以解決上述問題。

引入一個函數(shù)T將隨機(jī)變量轉(zhuǎn)變?yōu)榉恼龖B(tài)分布U～N(0,1 2 )的正態(tài)變量。該T函數(shù)的數(shù)學(xué)形式如下：

式中：Ti表示隨機(jī)變量xi的轉(zhuǎn)變函數(shù)；Φ-1[·]為標(biāo)準(zhǔn)正態(tài)隨機(jī)變量累積分布函數(shù)的逆函數(shù)；是隨機(jī)變量xi的累積分布函數(shù)。

采用高斯埃爾米特積分[4]，式（14）可以表示為

式中：g(·)表示系統(tǒng)響應(yīng)的功能函數(shù)；為第i1個變量關(guān)于第I1個積分點的高斯權(quán)重系數(shù)；為第i2個變量關(guān)于第I2個積分點的高斯權(quán)重系數(shù)；為第i1個變量關(guān)于第I1個積分點的高斯積分點；為第i2個變量關(guān)于第I2個積分點的高斯積分點；r，r1和r2均表示高斯積分點的個數(shù)。表1列出5階高斯埃爾米特積分節(jié)點和相應(yīng)的權(quán)重。當(dāng)求解式（15）時，為了保證足夠的精度，一般選擇5階高斯埃爾米特積分節(jié)點。

表1 高斯埃爾米特積分節(jié)點和權(quán)重[4]Tab.1 The nodes and weights of Gauss-Hermite integral

根據(jù)前4 階統(tǒng)計矩與原點矩的關(guān)系[24]，統(tǒng)計矩可表示為

式中：m1，m2，m3和m4是前4 階原點矩；D1，D2，D3和D4是前4階統(tǒng)計矩。

3.2 基于Kriging 模型的積分節(jié)點選取策略與收斂準(zhǔn)則

3.2.1 積分節(jié)點的選取策略

當(dāng)采用5 階高斯埃爾米特積分節(jié)點時，對于式（15）任意一個二維響應(yīng)函數(shù)，需要調(diào)用52次功能函數(shù)；對于任意一個一維響應(yīng)函數(shù)，需要調(diào)用5次功能函數(shù)。為了進(jìn)一步減少功能函數(shù)調(diào)用次數(shù)，本研究采用Kriging模型代替真實功能函數(shù)，并引入U學(xué)習(xí)函數(shù)逐漸增加新的積分點更新Kriging模型。

1）構(gòu)建一維響應(yīng)函數(shù)的Kriging模型

選擇坐標(biāo)軸上距離原點最近的2個積分點作為初始樣本點，采用U學(xué)習(xí)函數(shù)逐漸增加新的積分點更新Kriging模型，直到滿足收斂準(zhǔn)則。為了便于理解，圖1給出一維響應(yīng)函數(shù)Kriging 模型樣本點的選取示意圖。如圖1 所示，黑色實心圓（）表示坐標(biāo)原點，紅色空心圓（）表示建立一維響應(yīng)函數(shù)Kriging 模型的初始樣本點，黃色實心圓（）表示建立單變量響應(yīng)函數(shù)新增的樣本點，綠色實心圓（）表示建立Kriging模型未使用的樣本點。

與國外的機(jī)構(gòu)庫建設(shè)的高速發(fā)展相比，我國目前還處于起步階段。吳建中[3]2004年初發(fā)表文章探討了機(jī)構(gòu)庫對圖書館整體管理模式的沖擊，將知識庫的概念引入我國。2005年7月，北京大學(xué)圖書館率領(lǐng)國內(nèi)50多所高等院校圖書館聯(lián)合發(fā)表《圖書館合作與信息資源共享武漢宣言》，在宣言中明確指出我國高校圖書館應(yīng)“建設(shè)特色館藏，開展特色服務(wù)，建立一批特色學(xué)術(shù)機(jī)構(gòu)庫（Institutional Depository）”[4]。從那之后，機(jī)構(gòu)知識庫的建設(shè)在國內(nèi)，特別是我國高校圖書館逐步開啟[5]。

圖1 一維響應(yīng)函數(shù)Kriging 模型樣本點的選取示意圖Fig.1 Schematic diagram of sample points selection in Kriging model for one dimensional response functions

2）構(gòu)建二維響應(yīng)函數(shù)的Kriging模型

選擇坐標(biāo)軸上的9個積分點作為初始樣本點（紅色空心圓），采用U 學(xué)習(xí)函數(shù)逐漸增加新的積分點更新Kriging 模型，直到滿足收斂準(zhǔn)則。為了便于理解，圖2給出二維響應(yīng)函數(shù)Kriging模型樣本點的選取示意圖。如圖2所示，紅色空心圓（）表示建立二維響應(yīng)函數(shù)Kriging 模型的初始樣本點，黃色實心圓（）表示建立二維響應(yīng)函數(shù)新增的樣本點，綠色實心圓（）表示建立Kriging模型未使用的樣本點。

圖2 二維響應(yīng)函數(shù)Kriging 模型樣本點的選取示意圖Fig.2 Schematic diagram of sample points selection in Kriging model for two dimensional response functions

3.2.2 停止準(zhǔn)則

為了提高計算效率，當(dāng)更新的一維響應(yīng)函數(shù)和二維響應(yīng)函數(shù)Kriging模型分別達(dá)到足夠精度時，其Kriging模型則停止更新。為此，本文引入以下收斂準(zhǔn)則[21]：

式中：δ表示功能函數(shù)值的誤差；ε為允許的誤差值；表示Z(·)的預(yù)測值；Range(Z(u))=max(Z(u))-min(Z(u))，其中u為已存在的樣本點，u′為新加入的樣本點。本文中，ε取值為0.000 1。

3.2.3 計算統(tǒng)計矩

通過構(gòu)建一維響應(yīng)函數(shù)和二維響應(yīng)函數(shù)的Kriging模型，式（15）可進(jìn)一步變換為

3.3 實施步驟和流程圖

提出方法的實施步驟如表2所示，流程圖如圖3所示。

圖3 提出方法的流程圖Fig.3 Flowchart of the proposed method

表2 提出方法的實施步驟Tab.2 The implementation steps of the proposed method

4 算例驗證

4.1 算例1

采用文獻(xiàn)[5]中三維非線性功能函數(shù)驗證提出方法的有效性，其表達(dá)式為

式中，隨機(jī)變量x1，x2和x3是相互獨立且均服從均值為0.918，標(biāo)準(zhǔn)差為0.210 的韋布爾分布，其邊緣概率密度函數(shù)為。

將雙變量降維方法、提出的方法與蒙特卡洛模擬方法計算得到的前4 階統(tǒng)計矩進(jìn)行對比如表3 所示，可以看出傳統(tǒng)雙變量降維方法和提出方法獲得的前4階矩結(jié)果基本一致，且都接近蒙特卡洛模擬計算結(jié)果，表明計算精度能夠得到保證。在計算效率方面，提出方法調(diào)用功能函數(shù)次數(shù)比傳統(tǒng)雙變量降維方法少16次，從而提高了統(tǒng)計矩的計算效率。

表3 算例1 不同方法計算結(jié)果Tab.3 Computational results of different methods in example 1

4.2 算例2

考慮含5 個隨機(jī)變量的非線性功能函數(shù)[23]，其表達(dá)式為

式中：隨機(jī)變量相互獨立且服從正態(tài)分布；均值分別為μ1=1.2，μ2=2.4，μ3=50，μ4=25和μ5=10；標(biāo)準(zhǔn)差分別為σ1=0.36，σ2=0.072，σ3=3，σ4=7.5和σ5=5。

雙變量降維方法、提出的方法與蒙特卡洛模擬方法的計算精度和效率對比如表4 所示。從表4 可以看出，傳統(tǒng)雙變量降維方法和提出方法獲得的前4 階矩結(jié)果基本一致，且都接近蒙特卡洛模擬計算結(jié)果，再次證明提出方法具有較高的計算精度。在計算效率方面，提出方法功能函數(shù)調(diào)用次數(shù)約為傳統(tǒng)雙變量降維方法的一半，計算效率得到了顯著提升。與算例1相比，算例2隨機(jī)變量維數(shù)增多，提出方法在節(jié)省調(diào)用功能函數(shù)次數(shù)方面具有顯著的效果。

表4 算例2 不同方法計算結(jié)果Tab.4 Computational results of different methods in example 2

4.3 算例3

SCARA機(jī)器人[25]有4個關(guān)節(jié)，其中3個旋轉(zhuǎn)關(guān)節(jié)的軸線相互平行，在空間內(nèi)進(jìn)行定位和定向，1個移動關(guān)節(jié)用于實現(xiàn)末端件在垂直于水平面的運動。SCARA 機(jī)器人示意圖和機(jī)構(gòu)運動簡圖如圖4a）和圖4b）所示。表5列出了SCARA機(jī)器人的D-H參數(shù)[25]，其中L1=200，θ1=0和d3=250為確定值，其他不確定變量的信息如表6所示。

圖4 SCARA 機(jī)器人示意圖和機(jī)構(gòu)運動簡圖Fig.4 Schematic diagram and mechanism motion diagram of the SCARA robot

表5 SCARA 機(jī)器人D-H 參數(shù)[25]Tab.5 D-H parameters of the SCARA robot

表6 SCARA 機(jī)器人不確定參數(shù)Tab.6 Uncertain parameters of the SCARA robot

將雙變量降維方法、提出方法與蒙特卡洛模擬方法計算獲得的前4 階統(tǒng)計矩進(jìn)行對比，表7、表8 和表9分別列出工業(yè)機(jī)器人末端位置X，Y和Z方向坐標(biāo)統(tǒng)計矩的計算結(jié)果。

表7 不同方法計算SCARA 機(jī)器人X 方向統(tǒng)計矩Tab.7 The computational statistical moments by different methods in X direction for the SCARA robot

根據(jù)計算結(jié)果可以看出，提出的方法無論是在均值，標(biāo)準(zhǔn)差，偏度還是峰度，計算精度都與原雙變量降維方法相近，且都接近蒙特卡洛模擬方法的計算結(jié)果。此外，與原雙變量降維方法相比，提出方法在3個方向的功能函數(shù)調(diào)用次數(shù)都明顯減少，Z方向功能函數(shù)調(diào)用次數(shù)減少了46次，X和Y方向功能函數(shù)調(diào)用次數(shù)約為原雙變量降維方法的一半，證明提出方法在保證計算精度的前提下，實現(xiàn)統(tǒng)計矩的高效計算。

5 結(jié)論

本文將Kriging模型和U學(xué)習(xí)函數(shù)引入到雙變量降維方法中，提出一種改進(jìn)的雙變量降維方法，與傳統(tǒng)雙變量降維方法和蒙特卡洛模擬法進(jìn)行精度和效率的比較。兩個數(shù)值算例和SCARA機(jī)器人算例的計算結(jié)果表明：提出的方法能夠在保證計算精度的前提下，提高傳統(tǒng)雙變量降維方法的計算效率，減少功能函數(shù)的調(diào)用次數(shù)。此外，提出方法可應(yīng)用于其他復(fù)雜機(jī)械裝備的不確定性分析和可靠性優(yōu)化設(shè)計。