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

韓忠皓，張德權，楊美德，趙海文

（1.河北工業大學機械工程學院，天津 300401；2.湖南大學機械與運載工程學院，湖南 長沙 410082）

0 引言

功能函數的統計矩是機械系統進行可靠性分析的前提和基礎，其計算效率決定著機械系統可靠性分析的效率。隨著機械系統愈加復雜，高效地獲得功能函數統計矩變得愈加重要，已成為近年來的研究熱點。

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

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

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

1 雙變量降維方法

為提高計算系統響應統計矩的精度，Xu 等[18]在單變量降維方法的基礎上進一步提出了雙變量降維方法。與單變量降維方法類似，雙變量降維方法的實施包括3個步驟。

1）將一個多維響應函數等效為多個二維響應函數與多個一維響應函數的疊加，其數學形式如下：

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

2） 根據所有二維響應函數與一維響應函數的統計矩近似計算多維響應函數的統計矩。令，根據矩的定義，其l階近似響應函數原點矩可表示為

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

2 Kriging 模型和U 學習函數

2.1 Kriging 模型

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

式中：θ為相關參數；Rθ(Xi,Xj)為兩點的相關函數，其表達式為

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

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

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

給定一個預測點X0，在該點預測的函數值和方差為

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

2.2 U 學習函數

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

3 提出的方法

3.1 基于Gauss-Hermite 的統計矩評估

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

引入一個函數T將隨機變量轉變為服從正態分布U～N(0,1 2 )的正態變量。該T函數的數學形式如下：

式中：Ti表示隨機變量xi的轉變函數；Φ-1[·]為標準正態隨機變量累積分布函數的逆函數；是隨機變量xi的累積分布函數。

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

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

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

根據前4 階統計矩與原點矩的關系[24]，統計矩可表示為

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

3.2 基于Kriging 模型的積分節點選取策略與收斂準則

3.2.1 積分節點的選取策略

當采用5 階高斯埃爾米特積分節點時，對于式（15）任意一個二維響應函數，需要調用52次功能函數；對于任意一個一維響應函數，需要調用5次功能函數。為了進一步減少功能函數調用次數，本研究采用Kriging模型代替真實功能函數，并引入U學習函數逐漸增加新的積分點更新Kriging模型。

1）構建一維響應函數的Kriging模型

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

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

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

2）構建二維響應函數的Kriging模型

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

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

3.2.2 停止準則

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

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

3.2.3 計算統計矩

通過構建一維響應函數和二維響應函數的Kriging模型，式（15）可進一步變換為

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

采用文獻[5]中三維非線性功能函數驗證提出方法的有效性，其表達式為

式中，隨機變量x1，x2和x3是相互獨立且均服從均值為0.918，標準差為0.210 的韋布爾分布，其邊緣概率密度函數為。

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

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

4.2 算例2

考慮含5 個隨機變量的非線性功能函數[23]，其表達式為

式中：隨機變量相互獨立且服從正態分布；均值分別為μ1=1.2，μ2=2.4，μ3=50，μ4=25和μ5=10；標準差分別為σ1=0.36，σ2=0.072，σ3=3，σ4=7.5和σ5=5。

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

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

4.3 算例3

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

圖4 SCARA 機器人示意圖和機構運動簡圖Fig.4 Schematic diagram and mechanism motion diagram of the SCARA robot

表5 SCARA 機器人D-H 參數[25]Tab.5 D-H parameters of the SCARA robot

表6 SCARA 機器人不確定參數Tab.6 Uncertain parameters of the SCARA robot

將雙變量降維方法、提出方法與蒙特卡洛模擬方法計算獲得的前4 階統計矩進行對比，表7、表8 和表9分別列出工業機器人末端位置X，Y和Z方向坐標統計矩的計算結果。

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

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

5 結論

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