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

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

0 引言

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

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

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

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

1 雙變量降維方法

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

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

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

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

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

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

2.1 Kriging 模型

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

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

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

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

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

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

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

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

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

3 提出的方法

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

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

引入一個(gè)函數(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個(gè)變量關(guān)于第I1個(gè)積分點(diǎn)的高斯權(quán)重系數(shù)；為第i2個(gè)變量關(guān)于第I2個(gè)積分點(diǎn)的高斯權(quán)重系數(shù)；為第i1個(gè)變量關(guān)于第I1個(gè)積分點(diǎn)的高斯積分點(diǎn)；為第i2個(gè)變量關(guān)于第I2個(gè)積分點(diǎn)的高斯積分點(diǎn)；r，r1和r2均表示高斯積分點(diǎn)的個(gè)數(shù)。表1列出5階高斯埃爾米特積分節(jié)點(diǎn)和相應(yīng)的權(quán)重。當(dāng)求解式（15）時(shí)，為了保證足夠的精度，一般選擇5階高斯埃爾米特積分節(jié)點(diǎn)。

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

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

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

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

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

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

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

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

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

圖1 一維響應(yīng)函數(shù)Kriging 模型樣本點(diǎn)的選取示意圖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個(gè)積分點(diǎn)作為初始樣本點(diǎn)（紅色空心圓），采用U 學(xué)習(xí)函數(shù)逐漸增加新的積分點(diǎn)更新Kriging 模型，直到滿足收斂準(zhǔn)則。為了便于理解，圖2給出二維響應(yīng)函數(shù)Kriging模型樣本點(diǎn)的選取示意圖。如圖2所示，紅色空心圓（）表示建立二維響應(yīng)函數(shù)Kriging 模型的初始樣本點(diǎn)，黃色實(shí)心圓（）表示建立二維響應(yīng)函數(shù)新增的樣本點(diǎn)，綠色實(shí)心圓（）表示建立Kriging模型未使用的樣本點(diǎn)。

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

3.2.2 停止準(zhǔn)則

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

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

3.2.3 計(jì)算統(tǒng)計(jì)矩

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

3.3 實(shí)施步驟和流程圖

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

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

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

4 算例驗(yàn)證

4.1 算例1

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

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

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

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

4.2 算例2

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

式中：隨機(jī)變量相互獨(dú)立且服從正態(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。

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

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

4.3 算例3

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

圖4 SCARA 機(jī)器人示意圖和機(jī)構(gòu)運(yùn)動(dòng)簡(jiǎn)圖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

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

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

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

5 結(jié)論

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