孫東陽， 陳國平， 張保強

(1.重慶大學航空航天學院 重慶，400044) (2.南京航空航天大學機械結構力學及控制國家重點實驗室 南京，210016) (3.廈門大學航空航天學院 廈門，361005)

?

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

孫東陽1， 陳國平2， 張保強3

(1.重慶大學航空航天學院 重慶，400044) (2.南京航空航天大學機械結構力學及控制國家重點實驗室 南京，210016) (3.廈門大學航空航天學院 廈門，361005)

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

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

引 言

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

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

1 非線性質量-彈簧-阻尼系統(tǒng)動力學建模

(1)

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

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

(2)

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

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

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

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

(3)

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

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

2 置信區(qū)域法

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

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

(5)

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

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

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

(6)

(7)

其中

(8)

選取不同的e值，將得到結果x不同的分布。因此，當e從集合Ξ中選取不同的值時，將得到關于結果x的一組認知不確定性分布。

下面主要討論如下集合形式的不確定性結構

(9)

(10)

(11)

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

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

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

(12)

則xd?min{xd1,…，xdnSE}

張連長:“別攔他!誰也別攔他!我看他想怎么樣!路上我是你們帶隊,到了連隊我是你們連長!想跟連長打架,反教了!”

(13)

xu?max{xu1,…，xunSE}

(14)

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

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

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

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

3) 從每一個認知不確定性(區(qū)間參數(shù))中選擇一個樣本。因為采用LHS進行抽樣，所以選擇的抽樣數(shù)量M相對比較少。

4) 在給定認知不確定性樣本點的條件下，從隨機不確定性分布中選擇一組樣本。

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

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

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

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

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

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

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

4 數(shù)值仿真

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

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

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

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

圖4 質量塊位移的響應區(qū)域邊界以及標準模型響應Fig.4 Statistics of mass block displacement: boundaries of the response and response of the nominal model

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

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

5 結束語

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

[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] 張保強，陳國平，郭勤濤.不確定性熱彈耦合梁的固有振動分析[J].振動與沖擊,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] 趙寬，陳建軍，閻彬，等.含隨機參數(shù)的多體系統(tǒng)動力學分析[J].力學學報,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]曾開春，向錦武.高超聲速飛行器飛行動力學特性不確定分析[J].航空學報,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]杜永峰，李萬潤，李慧.基于測量數(shù)據(jù)不確定性的結構參數(shù)識別[J].振動、測試與診斷,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]胡海巖.應用非線性動力學[M].北京:航空工業(yè)出版社,2000:68-72.

10.16450/j.cnki.issn.1004-6801.2015.05.017

*重慶大學科研啟動基金資助項目(0240001104412)

2013-10-25；

2013-11-21

TB114; TH113

孫東陽，男，1985年6月生，博士、講師。主要研究方向為多體系統(tǒng)動力學建模及控制。曾發(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