孫東陽(yáng)， 陳國(guó)平， 張保強(qiáng)

(1.重慶大學(xué)航空航天學(xué)院 重慶，400044) (2.南京航空航天大學(xué)機(jī)械結(jié)構(gòu)力學(xué)及控制國(guó)家重點(diǎn)實(shí)驗(yàn)室 南京，210016) (3.廈門大學(xué)航空航天學(xué)院 廈門，361005)

?

隨機(jī)和認(rèn)知不確定性量化的置信區(qū)域法*

孫東陽(yáng)1， 陳國(guó)平2， 張保強(qiáng)3

(1.重慶大學(xué)航空航天學(xué)院 重慶，400044) (2.南京航空航天大學(xué)機(jī)械結(jié)構(gòu)力學(xué)及控制國(guó)家重點(diǎn)實(shí)驗(yàn)室 南京，210016) (3.廈門大學(xué)航空航天學(xué)院 廈門，361005)

針對(duì)非線性機(jī)械系統(tǒng)中混合不確定性量化的問(wèn)題，提出了隨機(jī)和認(rèn)知不確定性量化的置信區(qū)域法。首先，分別用概率論方法和區(qū)間方法來(lái)處理混合不確定性中的隨機(jī)不確定性和認(rèn)知不確定性，得到混合不確定性的置信區(qū)域；然后，為了在時(shí)間域內(nèi)對(duì)不確定性進(jìn)行傳播，對(duì)傳統(tǒng)雙層循環(huán)蒙特卡羅抽樣方法進(jìn)行了改進(jìn)；最后，以非線性質(zhì)量-彈簧-阻尼系統(tǒng)為例討論了基于混合不確定性分析方法的有效性。結(jié)果表明，同時(shí)考慮隨機(jī)不確定性和認(rèn)知不確定性，有利于提高系統(tǒng)設(shè)計(jì)的可靠性，為非線性機(jī)械系統(tǒng)的設(shè)計(jì)與精度分析提供了理論依據(jù)。

置信區(qū)域; 隨機(jī)不確定性; 認(rèn)知不確定性; 蒙特卡羅

引 言

在非線性機(jī)械系統(tǒng)的分析中，不確定性的處理作為一個(gè)研究熱點(diǎn)，已經(jīng)受到越來(lái)越廣泛的關(guān)注[1-5]。

不確定性按照其數(shù)據(jù)信息的掌握情況通常分為隨機(jī)不確定性和認(rèn)知不確定性[6-7]。隨機(jī)不確定性[8-9]源于系統(tǒng)行為的固有隨機(jī)性，是一種不可避免的不確定性。認(rèn)知不確定性[10]是由于缺乏數(shù)據(jù)引起的，是可以隨著認(rèn)識(shí)的深入而減少的不確定性。在處理隨機(jī)不確定性時(shí)，通常選用傳統(tǒng)概率論方法[11]。然而，在處理認(rèn)知不確定性時(shí)，數(shù)學(xué)結(jié)構(gòu)的選擇一直是一個(gè)挑戰(zhàn)性問(wèn)題[12]。模糊集理論[13]、區(qū)間分析[14-15]、證據(jù)理論[16]和可能性理論[17]已經(jīng)廣泛地應(yīng)用于研究認(rèn)知不確定性。當(dāng)系統(tǒng)中同時(shí)存在隨機(jī)和認(rèn)知不確定性時(shí)，處理這樣的問(wèn)題將面臨兩方面的挑戰(zhàn)。首先，是選擇合適的數(shù)學(xué)結(jié)構(gòu)去描述這兩種不確定性；其次，需要將兩種不確定性進(jìn)行分離，以分析各種不確定性對(duì)結(jié)果的貢獻(xiàn)程度[18]。目前，處理這類混合不確定性的方法主要分為概率邊界方法，二階概率論方法和Demspster-Shafer證據(jù)理論[19]。這3類混合不確定性分析方法在分析系統(tǒng)的時(shí)域響應(yīng)時(shí)都只能分析響應(yīng)的極值或特定時(shí)刻的系統(tǒng)響應(yīng)[20]。筆者提出了隨機(jī)和認(rèn)知不確定性量化的置信區(qū)域法，該方法能夠在時(shí)域內(nèi)對(duì)系統(tǒng)響應(yīng)進(jìn)行分析;隨后，提出了改進(jìn)的雙層循環(huán)蒙特卡羅抽樣方法; 最后，以非線性質(zhì)量-彈簧-阻尼系統(tǒng)為研究對(duì)象，用隨機(jī)和認(rèn)知不確定性量化的置信區(qū)域法對(duì)模型參數(shù)中存在的不確定性進(jìn)行了分析。

1 非線性質(zhì)量-彈簧-阻尼系統(tǒng)動(dòng)力學(xué)建模

(1)

其中：m為質(zhì)量;c為阻尼;k和kf分別為線性剛度和非線性剛度。

圖1 非線性質(zhì)量-彈簧-阻尼系統(tǒng)Fig.1 Nonlinear mass-spring-damper system

(2)

限定系統(tǒng)的非線性及阻尼比較弱

0≤ε?1,ζω0=εμ>0,μ=O(1)

該系統(tǒng)不具有精確解，可采用多尺度法研究解的一次近似。設(shè)

x(t)=x0(T0,T1)+εx1(T0,T1)

(3)

將式(3)代入式(2)，比較ε的同次冪得線性偏微分方程

通過(guò)對(duì)該方程組進(jìn)行求解可以得到Duffing系統(tǒng)的一次近似解。研究發(fā)現(xiàn)，該Duffing系統(tǒng)在簡(jiǎn)諧激勵(lì)下有可能出現(xiàn)主共振和次共振[22]。由于不研究該系統(tǒng)的非線性特性，因此選取的激勵(lì)頻率很低，以避開(kāi)共振頻率。在實(shí)際工程結(jié)構(gòu)中，由于加工精度問(wèn)題、試驗(yàn)條件限制和復(fù)雜環(huán)境影響，機(jī)械系統(tǒng)的結(jié)構(gòu)參數(shù)和材料參數(shù)都存在不確定性，使傳統(tǒng)基于確定性模型計(jì)算得到的響應(yīng)結(jié)果已經(jīng)不能滿足工程實(shí)際的需要，因此對(duì)機(jī)械系統(tǒng)進(jìn)行不確定性分析具有重要的意義。其中，有一部分不確定性參數(shù)可以通過(guò)試驗(yàn)找出其概率分布，這種參數(shù)不確定性為隨機(jī)不確定性。另一部分不確定性參數(shù)，難以用確定的隨機(jī)分布來(lái)描述，這就造成眾多認(rèn)知不確定性存在于系統(tǒng)中。當(dāng)系統(tǒng)同時(shí)存在隨機(jī)不確定性和認(rèn)知不確定性時(shí)，需要將這兩種不確定性進(jìn)行分離，以分析各種不確定性對(duì)總響應(yīng)的貢獻(xiàn)程度。為了便于分析，首先認(rèn)為參數(shù)不具有時(shí)變性；然后針對(duì)該系統(tǒng)存在的混合不確定性，提出了隨機(jī)和認(rèn)知不確定性量化的置信區(qū)域法。

2 置信區(qū)域法

在t時(shí)刻，包含隨機(jī)和認(rèn)知不確定性質(zhì)量塊位移的函數(shù)可以表示為

x(t)=f(a|e)=f(a1,…,anA|e1,…,enE)

(5)

其中：x(t)為t時(shí)刻質(zhì)量塊位移；a=[a1,…,anA]為隨機(jī)不確定性變量；e=[e1,…,enE]為認(rèn)知不確定性變量。

為了書寫簡(jiǎn)便，x(t)簡(jiǎn)寫為x。Ξi為e中元素ei的可能的值，所有e的可能值的集合可以表示為

Ξ=Ξ1×Ξ2×…×ΞnE

(6)

(7)

其中

(8)

選取不同的e值，將得到結(jié)果x不同的分布。因此，當(dāng)e從集合Ξ中選取不同的值時(shí)，將得到關(guān)于結(jié)果x的一組認(rèn)知不確定性分布。

下面主要討論如下集合形式的不確定性結(jié)構(gòu)

(9)

(10)

(11)

取q=0,當(dāng)x=xd時(shí)，Prbq(x)=2.5%；取q=1,當(dāng)x=xu時(shí)，Prbq(x)=97.5%。定義[xdxu]為t時(shí)刻質(zhì)量塊位移的95%置信區(qū)間，且當(dāng)q∈[0 1]時(shí)，有以下特性：a.對(duì)于Prbq(xdq)=2.5%，有xdq≥xd；b.對(duì)于Prbq(xuq)=97.5%，有xu≥xuq。在時(shí)域分析中，每個(gè)時(shí)刻都能得到一個(gè)95%置信區(qū)間[xd(t)xu(t)]，所有時(shí)刻置信區(qū)間的上(下)邊界即構(gòu)成質(zhì)量塊位移響應(yīng)的95%置信區(qū)域上(下)邊界。

3 改進(jìn)的雙層循環(huán)蒙特卡羅法

目前，常采用雙層循環(huán)蒙特卡羅法抽樣技術(shù)[21]來(lái)處理同時(shí)存在隨機(jī)和認(rèn)知不確定性的問(wèn)題。然而，傳統(tǒng)的雙層循環(huán)蒙特卡羅抽樣方法只能處理輸出結(jié)果為標(biāo)量的系統(tǒng)(如系統(tǒng)固有頻率，使用壽命等)，因此，筆者對(duì)雙層循環(huán)蒙特卡羅抽樣方法進(jìn)行了一些修改，使其能夠處理輸出為時(shí)間響應(yīng)的系統(tǒng)。

(12)

則xd?min{xd1,…，xdnSE}

張連長(zhǎng):“別攔他!誰(shuí)也別攔他!我看他想怎么樣!路上我是你們帶隊(duì),到了連隊(duì)我是你們連長(zhǎng)!想跟連長(zhǎng)打架,反教了!”

(13)

xu?max{xu1,…，xunSE}

(14)

圖2為改進(jìn)的雙層循環(huán)蒙特卡羅法抽樣技術(shù)的主要步驟，具體為：

圖2 改進(jìn)的雙層循環(huán)蒙特卡羅法抽樣技術(shù)Fig.2 Improved double-loop Monte Carlo sampling

1) 選擇處理認(rèn)知不確定性的抽樣數(shù)量M。當(dāng)采用拉丁超立方抽樣(latin hypercube sampling，簡(jiǎn)稱LHS)進(jìn)行抽樣時(shí)，最小抽樣數(shù)量M可以根據(jù)[21]M=m3+2進(jìn)行估算，其中m為認(rèn)知不確定性的數(shù)量。當(dāng)m不是很大時(shí)，為了提高計(jì)算精度應(yīng)該適當(dāng)增大M。

2) 選擇處理隨機(jī)不確定性的抽樣數(shù)量N。為了準(zhǔn)確描述系統(tǒng)響應(yīng)的分布，N一般要求較大。

3) 從每一個(gè)認(rèn)知不確定性(區(qū)間參數(shù))中選擇一個(gè)樣本。因?yàn)椴捎肔HS進(jìn)行抽樣，所以選擇的抽樣數(shù)量M相對(duì)比較少。

4) 在給定認(rèn)知不確定性樣本點(diǎn)的條件下，從隨機(jī)不確定性分布中選擇一組樣本。

5) 采用完整的抽樣序列計(jì)算系統(tǒng)響應(yīng)。

6) 判斷隨機(jī)不確定性的N個(gè)樣本是否已經(jīng)完成。如果否，回到步驟4)；如果是，繼續(xù)步驟7)。

7) 基于隨機(jī)不確定性計(jì)算得到的N個(gè)系統(tǒng)響應(yīng)，求各時(shí)刻系統(tǒng)響應(yīng)的95%(置信度2.5%～97.5%)置信區(qū)間。所有時(shí)刻置信區(qū)間上(下)邊界即構(gòu)成選定認(rèn)知不確定性樣本點(diǎn)情況下隨機(jī)不確定性樣本的95%置信區(qū)域上(下)邊界。

8) 判斷認(rèn)知不確定性的M個(gè)樣本是否已經(jīng)完成。如果是，繼續(xù)步驟9)，否則回到步驟3)。

9) 將M對(duì)置信區(qū)域上、下邊界線畫在同一張圖上，以顯示具有一定置信度的所有系統(tǒng)響應(yīng)邊界的全體。每一對(duì)邊界線顯示對(duì)應(yīng)的每個(gè)認(rèn)知不確定性樣本點(diǎn)的系統(tǒng)響應(yīng)的一個(gè)置信區(qū)域。

10) 對(duì)所有置信區(qū)域，分別找出上邊界的最大值和下邊界的最小值，該最大上邊界和最小下邊界圍成的區(qū)域即為同時(shí)考慮隨機(jī)和認(rèn)知不確定性時(shí)的95%置信區(qū)域。

11) 將系統(tǒng)響應(yīng)的最大上邊界和最小下邊界畫在同一張圖上，以顯示混合不確定性的影響。

4 數(shù)值仿真

筆者僅考慮系統(tǒng)參數(shù)的不確定性，假設(shè)m,c,k和kf相互獨(dú)立。對(duì)不確定性參數(shù)m,c,k和kf根據(jù)實(shí)驗(yàn)數(shù)據(jù)掌握情況分3種不確定性工況，分別研究質(zhì)量塊位移響應(yīng)的不確定性分布。在工況1中參數(shù)m,c,k和kf屬于隨機(jī)不確定性，假設(shè)其服從正態(tài)分布，給定均值和方差。在工況2中,參數(shù)m,c,k和kf屬于認(rèn)知不確定性，給定分布區(qū)間。在工況3中,參數(shù)k和kf屬于隨機(jī)不確定性，假設(shè)其服從正態(tài)分布，給定均值和方差，而m和c屬于認(rèn)知不確定性，給定分布區(qū)間。這3種工況的具體取值如表1所示。由于僅考慮系統(tǒng)參數(shù)的不確定性，取Y=10，ω=1。初始時(shí)刻，質(zhì)量塊位移和速度都為0。

表1 不同工況下參數(shù)的不確定性取值

在工況1中，由于參數(shù)m，c，k和kf都屬于隨機(jī)不確定性，采用傳統(tǒng)的蒙特卡羅抽樣技術(shù)對(duì)其抽樣1 000次。質(zhì)量塊位移的95%置信區(qū)域的邊界如圖3所示。在工況2中，由于參數(shù)m，c，k和kf都屬于認(rèn)知不確定性，給定參數(shù)的區(qū)間參數(shù)，所以得到的響應(yīng)也是區(qū)間數(shù)。這里采用LHS抽樣方法進(jìn)行抽樣1 000次，再將樣本點(diǎn)分別代入動(dòng)力學(xué)方程求響應(yīng)。最后，求出響應(yīng)的最大值和最小值，如圖4所示。在工況3中,由于參數(shù)m和c屬于認(rèn)知不確定性，而k和kf屬于隨機(jī)不確定性，用改進(jìn)的雙層循環(huán)蒙特卡羅法抽樣技術(shù)對(duì)該系統(tǒng)進(jìn)行不確定性分析。外層循環(huán)對(duì)應(yīng)于質(zhì)量和阻尼的認(rèn)知不確定性進(jìn)行抽樣，內(nèi)層循環(huán)對(duì)應(yīng)于線性剛度和非線性剛度的隨機(jī)不確定性進(jìn)行抽樣。質(zhì)量塊位移的95%置信區(qū)域邊界以及混合邊界如圖5所示。

圖3 質(zhì)量塊位移的置信區(qū)域邊界以及標(biāo)準(zhǔn)模型響應(yīng)Fig.3 Statistics of mass block displacement: boundaries of confidence region and response of the nominal model

圖4 質(zhì)量塊位移的響應(yīng)區(qū)域邊界以及標(biāo)準(zhǔn)模型響應(yīng)Fig.4 Statistics of mass block displacement: boundaries of the response and response of the nominal model

圖5 質(zhì)量塊位移的置信區(qū)域邊界以及混合邊界Fig.5 Statistics of mass block displacement: boundaries of confidence region for aleatory uncertainty and mixed-uncertainty

由圖3～圖5可知，工況1中參數(shù)為隨機(jī)不確定性，因?yàn)榧僭O(shè)參數(shù)為正態(tài)分布，所以能夠得到每個(gè)時(shí)刻質(zhì)量塊位移的概率分布，從而得到每個(gè)時(shí)刻質(zhì)量塊位移的95%置信區(qū)間和質(zhì)量塊位移的95%置信區(qū)域。工況2中參數(shù)為認(rèn)知不確定性，因?yàn)榧僭O(shè)參數(shù)為區(qū)間數(shù)，所以同樣采用區(qū)間數(shù)即可描述質(zhì)量塊位移的不確定性。工況3中，由于參數(shù)同時(shí)存在隨機(jī)和認(rèn)知不確定性，無(wú)法用傳統(tǒng)的方法進(jìn)行不確定性分析，采用改進(jìn)的雙層循環(huán)蒙特卡羅法抽樣技術(shù)對(duì)其進(jìn)行不確定性分析，得到每個(gè)時(shí)刻質(zhì)量塊位移的概率邊界，再由概率邊界得到95%置信區(qū)間,進(jìn)而得到95%置信區(qū)域?？梢?jiàn)，考慮參數(shù)存在隨機(jī)不確定性或混合不確定性時(shí)，得到質(zhì)量塊位移的響應(yīng)為置信區(qū)域，該響應(yīng)結(jié)果存在一定的置信度?？紤]參數(shù)存在認(rèn)知不確定性時(shí)，得到質(zhì)量塊位移的響應(yīng)區(qū)域。通過(guò)圖5可以看出,僅考慮隨機(jī)不確定性所得到的系統(tǒng)響應(yīng)的95%置信區(qū)域包含于混合不確定性的95%置信區(qū)域內(nèi)。

5 結(jié)束語(yǔ)

為了對(duì)同時(shí)存在隨機(jī)和認(rèn)知不確定性的非線性機(jī)械系統(tǒng)進(jìn)行分析，筆者提出了隨機(jī)和認(rèn)知不確定性量化的置信區(qū)域法。該方法不僅能夠在整個(gè)時(shí)間域內(nèi)單獨(dú)分析隨機(jī)不確定性和認(rèn)知不確定性對(duì)響應(yīng)的影響，而且還能夠得到混合不確定性對(duì)響應(yīng)的影響。由數(shù)值仿真可知，同時(shí)考慮參數(shù)的隨機(jī)和認(rèn)知不確定性得到的置信區(qū)域邊界比只考慮隨機(jī)不確定性時(shí)的邊界要寬。因此，在對(duì)非線性機(jī)械系統(tǒng)進(jìn)行分析時(shí)考慮這兩種不確定性的影響，將使系統(tǒng)工作時(shí)更加安全。

[1] Schu?ller G I. On the treatment of uncertainties in structural mechanics and analysis[J]. Computers & Structures, 2007,85(5):235-343.

[2] Li L, Sandu C. On the impact of cargo weight, vehicle parameters, and terrain characteristics on the prediction of traction for off-road vehicles[J]. Journal of Terramechanics, 2007,44(3):221-238.

[3] Sandu A, Sandu C, Ahmadian M. Modeling multibody systems with uncertainties. Part I: theoretical and computational aspects[J]. Multibody System Dynamics, 2006,15(4):373-395.

[4] Sandu C, Sandu A, Ahmadian M. Modeling multibody systems with uncertainties. Part II, numerical applications[J]. Multibody System Dynamics, 2006,15(3):241-262.

[5] An D, Choi J, Schmitz T L, et al. In situ monitoring and prediction of progressive joint wear using Bayesian statistics[J]. Wear, 2011,270(11):828-838.

[6] Hofer E, Kloos M, Krzykacz-Hausmann B, et al. An approximate epistemic uncertainty analysis approach in the presence of epistemic and aleatory uncertainties[J]. Reliability Engineering & System Safety, 2002,77(3):229-238.

[7] Huang Hongzhong, Zhang Xudong. Design optimization with discrete and continuous variables of aleatory and epistemic uncertainties[J]. ASME Journal of Mechanical Design, 2009,131(3):310061-310068.

[8] 張保強(qiáng)，陳國(guó)平，郭勤濤.不確定性熱彈耦合梁的固有振動(dòng)分析[J].振動(dòng)與沖擊,2012,31(19):160-164.

Zhang Baoqiang, Chen Guoping, Guo Qintao. Free vibration analysis of a thermoelastic coupled beam with material uncertainty[J]. Journal of Vibration and Shock, 2012,31(19):160-164. (in Chinese)

[9] 趙寬，陳建軍，閻彬，等.含隨機(jī)參數(shù)的多體系統(tǒng)動(dòng)力學(xué)分析[J].力學(xué)學(xué)報(bào),2012,44(4):802-806.

Zhao Kuan, Chen Jianjun, Yan Bin, et al. Dynamic analysis of multibody systems with probabilistic parameters[J]. Chinese Journal of Theoretical and Applied Mechanics, 2012,44(4):802-806. (in Chinese)

[10]曾開(kāi)春，向錦武.高超聲速飛行器飛行動(dòng)力學(xué)特性不確定分析[J].航空學(xué)報(bào),2012,33(4): 1-11.

Zeng Kaichun, Xiang Jinwu. Flight dynamic characteristics of hypersonic vehicles with uncertain parameters[J]. Acta Aeronautica et Astronautica Sinica, 2012,33(4):1-11. (in Chinese)

[11]Nilsson N J. Probabilistic logic[J]. Artificial Intelligence, 1986,28(1):71-87.

[12]Oberkampf W L, Helton J C, Joslyn C A, et al. Challenge problems: uncertainty in system response given uncertain parameters[J]. Reliability Engineering & System Safety, 2004,85(1):11-19.

[13]Ross T J. Fuzzy logic with engineering applications[M]. 2nd ed. New York: Wiley, 2004:6-20.

[14]Merlet J P. Solving the forward kinematics of a Gough-type parallel manipulator with interval analysis[J]. The International Journal of Robotics Research, 2004,23(3):221-235.

[15]杜永峰，李萬(wàn)潤(rùn)，李慧.基于測(cè)量數(shù)據(jù)不確定性的結(jié)構(gòu)參數(shù)識(shí)別[J].振動(dòng)、測(cè)試與診斷,2012,32(4):629-633.

Du Yongfeng, Li Wanrun, Li Hui. Structural parameters identification based on uncertainty of measurement data[J]. Journal of Vibration, Measurement & Diagnosis, 2012,32(4):629-633. (in Chinese)

[16]Beynon M, Curry B, Morgan P. The Dempster-Shafer theory of evidence: an alternative approach to multicriteria decision modelling[J]. Omega, 2000,28(1):37-50.

[17]Dubois D. Possibility theory and statistical reasoning[J]. Computational Statistics & Data Analysis, 2006,51(1):47-69.

[18]Helton J C, Johnson J D, Oberkampf W L, et al. Representation of analysis results involving aleatory and epistemic uncertainty[J]. International Journal of General Systems, 2010,39(6):605-646.

[19]Sentz K, Ferson S. Probabilistic bounding analysis in the quantification of margins and uncertainties[J]. Reliability Engineering & System Safety, 2011,96(9):1126-1136.

[20]Helton J C. Quantification of margins and uncertainties: conceptual and computational basis[J]. Reliability Engineering & System Safety, 2011,96(9):976-1013.

[21]Oberkampf W L, Roy C J. Verification and validation in scientific computing[M]. New York: Cambridge University Press, 2010:606-610.

[22]胡海巖.應(yīng)用非線性動(dòng)力學(xué)[M].北京:航空工業(yè)出版社,2000:68-72.

10.16450/j.cnki.issn.1004-6801.2015.05.017

*重慶大學(xué)科研啟動(dòng)基金資助項(xiàng)目(0240001104412)

2013-10-25；

2013-11-21

TB114; TH113

孫東陽(yáng)，男，1985年6月生，博士、講師。主要研究方向?yàn)槎囿w系統(tǒng)動(dòng)力學(xué)建模及控制。曾發(fā)表《Model reduction of a multibody system including a very flexible beam element》(《Journal of Mechanical Science and Technology》2014,Vol.28,No.28)等論文。 E-mail：sundongyang@cqu.edu.cn