基于統(tǒng)計(jì)信息的多體系統(tǒng)區(qū)間不確定性分析

陳光宋，錢(qián)林方，王明明，林 通

(南京理工大學(xué) 機(jī)械工程學(xué)院，南京 210094)

多體系統(tǒng)動(dòng)力學(xué)被廣泛應(yīng)用于機(jī)構(gòu)、機(jī)器人、車輛及機(jī)床等機(jī)械系統(tǒng)，這些系統(tǒng)通過(guò)連接副和力元將多個(gè)零部件連接成有機(jī)的整體。盡管關(guān)于多體系統(tǒng)動(dòng)力學(xué)的建模方法已經(jīng)相對(duì)成熟，但是大多數(shù)的研究都假定系統(tǒng)內(nèi)所有參數(shù)均為確定性的。然而，現(xiàn)實(shí)中由于載荷、幾何尺寸、材料屬性、制造和裝配過(guò)程以及使用中的磨損等因素使得參數(shù)具有不確定性，這些參數(shù)的不確定性可能導(dǎo)致系統(tǒng)性能的顯著變化。因此，研究含不確定性參數(shù)的多體系統(tǒng)動(dòng)力學(xué)問(wèn)題對(duì)確保機(jī)械系統(tǒng)的安全，提高系統(tǒng)的可靠性和穩(wěn)健性等有重要意義。目前，處理不確定性問(wèn)題的方法主要包括概率方法[1-4]、模糊集方法[5-7]和區(qū)間方法[8-10]。概率方法和模糊方法所需的概率分布信息和模糊隸屬度函數(shù)在實(shí)際過(guò)程中不易獲取，需要大量的試驗(yàn)，甚至無(wú)法準(zhǔn)確獲取。相反，區(qū)間方法所需的信息只是不確定性參數(shù)的上下界限，這兩個(gè)參數(shù)獲得的難度小很多，使用起來(lái)也更靈活和方便。然而，區(qū)間方法的問(wèn)題是由“包裹效應(yīng)”引起的區(qū)間擴(kuò)大效應(yīng)，如果無(wú)法有效控制區(qū)間的擴(kuò)大，則這種方法將無(wú)法適用。Moore[11]對(duì)區(qū)間算法中用到的自然擴(kuò)張函數(shù)、中心擴(kuò)張函數(shù)和Taylor擴(kuò)張函數(shù)進(jìn)行了比較，其中Taylor擴(kuò)張函數(shù)得到的區(qū)間最窄。

動(dòng)力學(xué)問(wèn)題的數(shù)學(xué)模型通常為強(qiáng)非線性的微分方程，解算過(guò)程需要采用較小的步長(zhǎng)，將區(qū)間方法用于微分方程求解時(shí)，其“包裹效應(yīng)”將更加顯明，隨著求解時(shí)間的增加，區(qū)間方法很可能導(dǎo)致求解發(fā)散，因此如何控制動(dòng)力學(xué)問(wèn)題的區(qū)間擴(kuò)大問(wèn)題，是區(qū)間方法在動(dòng)力學(xué)問(wèn)題中應(yīng)用的關(guān)鍵。目前，國(guó)內(nèi)外學(xué)者提出了許多方案來(lái)控制區(qū)間的過(guò)估計(jì)問(wèn)題，例如區(qū)間泰勒級(jí)數(shù)方法[12]，泰勒模型方法[13]以及近年來(lái)由吳景錸提出的切比雪夫擴(kuò)張函數(shù)方法[14-16]等，這些方法能夠在一定程度上控制區(qū)間方法處理動(dòng)態(tài)不確定性問(wèn)題的過(guò)估計(jì)問(wèn)題。其中泰勒展開(kāi)的方法是侵入式的方法，需要修改動(dòng)力學(xué)方程，而且多數(shù)僅能處理區(qū)間較小的情況，其他非侵入的方法在處理高維問(wèn)題時(shí)會(huì)遇到計(jì)算量大等問(wèn)題[17]。

本文提出一種區(qū)間方法用于解決含不確定性參數(shù)的動(dòng)力學(xué)問(wèn)題，利用極大極小距離的拉丁超立方采樣(Latin Hypercube Sampling,LHS)方法獲得樣本點(diǎn)，根據(jù)樣本信息統(tǒng)計(jì)獲得響應(yīng)的前四階統(tǒng)計(jì)矩，并通過(guò)Bootstrap方法提高矩的計(jì)算精度，進(jìn)而利用極大熵方法獲得響應(yīng)的概率密度函數(shù)和累積概率函數(shù)，將累積概率密度函數(shù)在樣本極大極小點(diǎn)處進(jìn)行泰勒展開(kāi)，并令分布函數(shù)在估計(jì)區(qū)間的上界和下界上的值為1和0，反求出輸出參數(shù)區(qū)間的上下界。最后，通過(guò)兩個(gè)解析函數(shù)驗(yàn)證本文方法的收斂性，通過(guò)兩個(gè)標(biāo)準(zhǔn)的機(jī)械系統(tǒng)動(dòng)力學(xué)問(wèn)題驗(yàn)證本文方法的有效性，通過(guò)一個(gè)工程算例驗(yàn)證本文方法的適用性。

1 機(jī)械系統(tǒng)不確定多體動(dòng)力學(xué)控制方程

多體系統(tǒng)動(dòng)力學(xué)控制方程通常為指標(biāo)-3的微分代數(shù)方程組(Differential Algebraic Equations,DAEs)[18]，其中微分方程部分描述狀態(tài)量對(duì)時(shí)間的導(dǎo)數(shù)，代數(shù)方程部分則是系統(tǒng)的完整約束條件。多體動(dòng)力學(xué)系統(tǒng)的控制方程可寫(xiě)成如下形式[19]

(1)

(2)

其中，

(3)

為了防止常用的積分算法在數(shù)值求解過(guò)程中的違約問(wèn)題，Baumgarte[20]基于控制反饋策略提出了一種約束穩(wěn)定方案

(4)

其中，

(5)

上述多體動(dòng)力學(xué)控制方程為確定性的，考慮系統(tǒng)參數(shù)的不確定性，則控制方程變?yōu)椴淮_定性控制方程

(6)

式中：x為不確定性參數(shù)。

為表述方便，考慮單個(gè)不確定性變量x，采用區(qū)間方法刻畫(huà)不確定性參數(shù)，則參數(shù)x的區(qū)間可定義為

(7)

定義區(qū)間[x]的區(qū)間中點(diǎn)xc、區(qū)間半徑xr和區(qū)間寬度xw為

(8)

(9)

xw=2xr

(10)

中心化的區(qū)間可表示為

[Δx]=[-xr,xr]

(11)

則真實(shí)的區(qū)間[x]可重寫(xiě)為如下形式

[x]=xc+[Δx]

(12)

2 基于統(tǒng)計(jì)信息的區(qū)間分析方法

2.1 極大極小距離的LHS

LHS具有非常靈活的特征，然而設(shè)計(jì)點(diǎn)的分布是任意的，不能很好的充滿設(shè)計(jì)空間，導(dǎo)致參數(shù)設(shè)計(jì)矩陣不能很好地代表設(shè)計(jì)空間。本文綜合計(jì)算效率和適用的靈活性，基于極大極小距離準(zhǔn)則，采用近似優(yōu)化的方案，獲得較好的樣本點(diǎn)。

定義多維空間中兩點(diǎn)a和b的歐氏距離d為

(13)

式中：p為維度；d越小表示多維空間中點(diǎn)a和b越靠近。

為使設(shè)計(jì)點(diǎn)在空間中均勻分布，可通過(guò)使多維空間中距離最小的兩個(gè)點(diǎn)a和b的距離最大來(lái)實(shí)現(xiàn)。實(shí)現(xiàn)的方法可分為優(yōu)化方法和準(zhǔn)優(yōu)化方法。優(yōu)化方法在維度p和樣本點(diǎn)數(shù)N較大的情況下，計(jì)算量將非常大，實(shí)際上可采用近似均勻的策略，通過(guò)準(zhǔn)優(yōu)化的方法獲得較均勻的樣本，實(shí)現(xiàn)步驟如下[21]：

步驟1令I(lǐng)=mN，一般m=5，N為需要生成隨機(jī)矩陣的行數(shù)，在[0,1]區(qū)間上產(chǎn)生I×p的均勻隨機(jī)數(shù)，構(gòu)成隨機(jī)矩陣Rd；

步驟2對(duì)I中的第i行元素，計(jì)算歐氏距離dij，j∈1,2,…,I,i≠j，取出其中距離最小的兩個(gè)并取平均值；

步驟3存儲(chǔ)第i行元素的最小平均距離，i=i+1，返回步驟2，直到i=I；

步驟4刪除平均距離最小的一組數(shù)據(jù)，因此I=I-1；

步驟5i=I，返回步驟2，直到I=N，形成均勻隨機(jī)數(shù)矩陣R；

步驟6根據(jù)矩陣R，結(jié)合輸入?yún)?shù)的區(qū)間范圍，產(chǎn)生最終輸入?yún)?shù)的設(shè)計(jì)矩陣D。

2.2 各階統(tǒng)計(jì)矩的Bootstrap求解方法

將設(shè)計(jì)矩陣D中的輸入?yún)?shù)代入多體系統(tǒng)動(dòng)力學(xué)方程，通過(guò)解算獲得輸出參數(shù)樣本Y，則樣本Y的前四階統(tǒng)計(jì)矩為：

樣本均值

(14)

樣本方差

(15)

樣本偏度系數(shù)

(16)

樣本峰度系數(shù)

(17)

通常在樣本數(shù)量較多時(shí)，通過(guò)式(14)～式(17)可獲得較準(zhǔn)確的解，然而在樣本數(shù)量較少時(shí)，對(duì)低階統(tǒng)計(jì)矩的計(jì)算較準(zhǔn)確，對(duì)高階統(tǒng)計(jì)矩的估算可能產(chǎn)生較大的誤差。為此，本文通過(guò)Bootstrap重采樣方法來(lái)減小各階統(tǒng)計(jì)矩估算的偏差。根據(jù)初始獲得的樣本，采用有放回的抽樣方法，獲得B組Boostrap樣本，繼而根據(jù)式(14)～式(17)計(jì)算獲得每組Bootstrap樣本的各階統(tǒng)計(jì)矩，并取均值獲得最終輸出參數(shù)的統(tǒng)計(jì)矩

(18)

(19)

(20)

(21)

2.3 參數(shù)分布的極大熵估計(jì)

最大熵原理[22]可以描述為在給定的條件下，在所有可能的分布中存在一個(gè)使信息熵取極大值的分布，該分布是最小偏見(jiàn)的，由此可以得到一種構(gòu)造最佳分布的方法。輸出參數(shù)Y的信息熵為

(22)

式中：fY(y)為Y的概率密度函數(shù)；H為信息熵；Ω為積分空間。

假定fY(y)為如下形式

(23)

式中：ai為待定系數(shù)。我們的目標(biāo)是尋找到一個(gè)分布函數(shù)(即一組系數(shù)ai)使得信息熵的值最大，并將隨機(jī)變量Y的前四階統(tǒng)計(jì)矩作為約束條件，可建立優(yōu)化模型如下

(24)

基于梯度的優(yōu)化方法能夠得到較高精度的解，初值的選擇對(duì)優(yōu)化求解的收斂性有很大影響。采用基于正態(tài)分布的初始點(diǎn)則能更快更準(zhǔn)確地搜索到最大熵對(duì)應(yīng)的系數(shù)，初始點(diǎn)可設(shè)置為

(25)

通過(guò)求解上述優(yōu)化問(wèn)題，可得系數(shù)ai(i=1,2,3,4)，進(jìn)而確定概率密度函數(shù)fY(y)，通過(guò)對(duì)fY(y)進(jìn)行積分可得隨機(jī)變量Y的分布函數(shù)FY(y)為

(26)

由概率密度函數(shù)的定義可知，F(xiàn)Y(y)∈[0,1]，并且FY(y)單調(diào)遞增，F(xiàn)Y(y)的導(dǎo)數(shù)即為fY(y)，且fY(y)始終為正，這些性質(zhì)是本文求解參數(shù)區(qū)間的基本保證。

2.4 基于泰勒展開(kāi)的參數(shù)區(qū)間估計(jì)

假定ymin和ymax分別為“2.1”節(jié)和“2.2”節(jié)中通過(guò)極大極小距離LHS得到的有限樣本的極小值和極大值。根據(jù)“2.3”節(jié)獲得的累積概率密度函數(shù)FY(y)，在ymin和ymax處對(duì)其進(jìn)行一階泰勒展開(kāi)，可得

(27)

(28)

(29)

(30)

(31)

圖1 計(jì)算流程圖Fig.1 Flowchart of computation process

3 算 例

首先通過(guò)兩個(gè)解析函數(shù)算例驗(yàn)證本文方法的收斂性，然后通過(guò)兩個(gè)典型的多體機(jī)械系統(tǒng)(雙擺機(jī)構(gòu)和曲柄滑塊機(jī)構(gòu))來(lái)驗(yàn)證本文提出區(qū)間分析方法的有效性，最后將本文方法應(yīng)用于一個(gè)工程算例。

3.1 解析算例

考慮解析函數(shù)表達(dá)式如下

(32)

y2=-sinx1-2(sinx2)2

(33)

式中：x1∈[-π,π]，x2∈[-π/2,π/2]，輸出參數(shù)在輸入?yún)?shù)定義域內(nèi)的曲線和曲面如圖2(a)和圖2(b)所示，其中輸出參數(shù)的上下界均可通過(guò)優(yōu)化的方法獲得。為驗(yàn)證本文方法的收斂特性，在輸入?yún)?shù)空間內(nèi)分別采樣20，40，60，80和100個(gè)樣本點(diǎn)，計(jì)算輸出參數(shù)上下界如圖3(a)和圖3(b)所示。點(diǎn)劃線為本文計(jì)算出的上下界，虛線為通過(guò)優(yōu)化方法獲得的上下界。中間過(guò)程的各階統(tǒng)計(jì)矩，如表1和表2所示。

圖2 解析函數(shù)Fig.2 The analytical function

表1 函數(shù)1參數(shù)統(tǒng)計(jì)矩Tab.1 The statistical moment of the output of funtion 1

圖3 不同采樣點(diǎn)下的區(qū)間Fig.3 The interval evaluated by different number of samples

表2 函數(shù)2參數(shù)統(tǒng)計(jì)矩Tab.2 The statistical moment of the output of funtion 2

從圖3可知，本文方法不僅可獲得較緊的包含輸出響應(yīng)的區(qū)間，而且在較少樣本的情況下也能獲得較好的解，具有很好的收斂性。

3.2 雙擺機(jī)構(gòu)

雙擺機(jī)構(gòu)示意圖如圖4所示。其動(dòng)力學(xué)方程為典型的多體動(dòng)力學(xué)模型。參照Wu等研究中的參數(shù)，兩根桿的長(zhǎng)度分別為l1=l2=1 m，質(zhì)心分別為O1和O2，質(zhì)量分別為m1=m2=1 kg，重力加速度為9.8 m/s2。

圖4 雙擺機(jī)構(gòu)Fig.4 Schematic of the double pendulum

(34)

式中：

(35)

式中：diag([])為對(duì)角矩陣，對(duì)角元素為[]；慣性矩

(36)

考慮兩桿長(zhǎng)均為不確定性參數(shù)，其變化范圍為名義值的5%，則有

(37)

式中：ξ1,ξ2∈[-1,1]，相應(yīng)的質(zhì)量和慣量矩也為不確定性參數(shù)，可表示為

(38)

(39)

通過(guò)四階龍哥庫(kù)塔法解算動(dòng)力學(xué)方程，步長(zhǎng)為0.001 s，采用本文的方法和優(yōu)化方法計(jì)算從0～10 s時(shí)間段內(nèi)兩桿的質(zhì)心位置，本文方法的樣本點(diǎn)數(shù)取為100，為更有效地對(duì)比計(jì)算結(jié)果增加了區(qū)間中點(diǎn)以及任意一組確定性參數(shù)的響應(yīng)值。計(jì)算結(jié)果如圖5～圖8所示。

圖5 桿1質(zhì)心水平位置Fig.5 Horizontal displacement of the top rod center

圖6 桿1質(zhì)心垂直位置Fig.6 Vertical displacement of the top rod center

圖7 桿2質(zhì)心水平位置Fig.7 Horizontal displacement of the bottom rod center

圖8 桿2質(zhì)心垂直位置Fig.8 Vertical displacement of the bottom rod center

從圖5～圖8可知，本文方法所求區(qū)間能夠較好地包含優(yōu)化方法獲得的結(jié)果，且能夠長(zhǎng)期抑制區(qū)間過(guò)估計(jì)問(wèn)題，從推導(dǎo)過(guò)程可知，與傳統(tǒng)區(qū)間方法不同，本文所提出的方法無(wú)導(dǎo)致區(qū)間放大的區(qū)間運(yùn)算過(guò)程，可避免區(qū)間的“包裹效應(yīng)”，其區(qū)間估計(jì)的準(zhǔn)確程度僅取決于估計(jì)輸出參數(shù)統(tǒng)計(jì)矩所需樣本點(diǎn)的數(shù)量，不會(huì)隨著輸入?yún)?shù)個(gè)數(shù)的增加而大幅增加。

3.3 曲柄滑塊

曲柄滑塊機(jī)構(gòu)示意圖如圖9所示。圖9中A點(diǎn)、B點(diǎn)和C點(diǎn)分別為曲柄、連桿和滑塊的質(zhì)心。Oxy為慣性坐標(biāo)系，局部坐標(biāo)系O1x1y1，O2x2y2和O3x3y3與3個(gè)構(gòu)件質(zhì)心固接，θ1和θ2為曲柄和連桿的轉(zhuǎn)動(dòng)角。在O點(diǎn)處施加一轉(zhuǎn)矩τ，滑塊末端與彈簧阻尼機(jī)構(gòu)相連，剛度為k，阻尼為c。當(dāng)θ1和θ2為0時(shí)，該彈簧處于平衡位置。圖中l(wèi)1和l2分別表示曲柄和連桿的長(zhǎng)度，m1，m2和m3分別表示曲柄、連桿和滑塊的質(zhì)量，參照Wu等研究的參數(shù)，各參數(shù)取值如表3所示。

圖9 曲柄滑塊機(jī)構(gòu)Fig.9 Schematic of slider crank

表3 曲柄滑塊參數(shù)Tab.3 Parameters of slider crank

選擇廣義坐標(biāo)q=[x1,y1,θ1,x2,y2,θ2,x3]T，則控制方程為

(40)

其中，

(41)

慣性矩為

(42)

考慮曲柄長(zhǎng)度存在不確定性，其不確定范圍為名義值的1%，即

(43)

質(zhì)量和慣性矩的不確定性可表示為

(44)

考慮力矩τ在其名義值附近有1%的不確定性，即

(45)

初始條件如下

(46)

通過(guò)四階龍哥庫(kù)塔法解算動(dòng)力學(xué)方程，步長(zhǎng)為0.001 s，采用本文的方法和優(yōu)化方法計(jì)算從0～2 s時(shí)間段內(nèi)的響應(yīng)，本文方法的樣本點(diǎn)數(shù)取為200，并增加了區(qū)間中點(diǎn)以及任意一組確定性參數(shù)的響應(yīng)值作為比較。計(jì)算結(jié)果如圖10～圖12，結(jié)果顯示本文所得的結(jié)果能夠完全包含實(shí)際解，即使在多個(gè)周期的往復(fù)運(yùn)動(dòng)后，估計(jì)的區(qū)間依然可靠，并且區(qū)間的放大現(xiàn)象得到很好的抑制。

圖10 曲柄轉(zhuǎn)角Fig.10 Rotation angle of crank

圖11 連桿轉(zhuǎn)角Fig.11 Rotation angle of connecting rod

圖12 滑塊位置Fig.12 Displacement of piston

3.4 火炮搖架俯仰運(yùn)動(dòng)

某大口徑火炮上裝部分示意圖，如圖13所示。其主要部件包括身管、搖架、上架、反后坐裝置、炮尾等，火炮的高低射角通過(guò)高平機(jī)液壓油缸的伸縮實(shí)現(xiàn)，此處為分析問(wèn)題方便，假定在火炮發(fā)射過(guò)程中上架固定。由于高平機(jī)液壓油缸受溫度、容氣量等的影響，使得火炮每次發(fā)射過(guò)程，搖架運(yùn)動(dòng)規(guī)律不一，對(duì)于長(zhǎng)身管而言，搖架俯仰運(yùn)動(dòng)情況能夠間接反映炮口擾動(dòng)情況，雖無(wú)指標(biāo)要求，但在火炮設(shè)計(jì)過(guò)程中，通常需要通過(guò)動(dòng)力學(xué)模型進(jìn)行檢驗(yàn)，在射擊精度試驗(yàn)時(shí)，也需要測(cè)試搖架的運(yùn)動(dòng)情況。因此，基于本文中的區(qū)間方法，對(duì)火炮發(fā)射過(guò)程搖架的運(yùn)動(dòng)進(jìn)行不確定性分析。

圖13 火炮上裝部分Fig.13 The upper part of howitzer

高平機(jī)油缸示意圖如圖14所示，其中C腔與蓄能器相連，A腔和B腔為工作腔，初始?jí)簭?qiáng)分別為PA和PB，初始液壓油體積分別為VA和VB，工作面積分別為

圖14 高平機(jī)液壓缸Fig.14 The hydraulic cylinder of elevating-equilibrium mechanism

(47)

(48)

液壓油體積模量定義為

(49)

式中：Eh為不確定性參數(shù)，受溫度和含氣量影響，這里取值范圍用區(qū)間表示為

(50)

(51)

另外還需考慮結(jié)構(gòu)參數(shù)D1，d，D2以及液壓油體積的誤差，即

(52)

(53)

(54)

(55)

(56)

由此得到A和B腔液壓油等效剛度分別為

(57)

(58)

則高平機(jī)的力學(xué)模型為

(59)

(60)

圖15 搖架俯仰角速度Fig.15 The rotation angle of cradle

由圖15可知，在火炮發(fā)射過(guò)程中搖架發(fā)生了周期性運(yùn)動(dòng)，而且考慮高平機(jī)參數(shù)不確定性的情況下，搖架的運(yùn)動(dòng)是動(dòng)態(tài)不確定的，通過(guò)本文估計(jì)的區(qū)間能夠較緊的包含優(yōu)化方法得到的結(jié)果，說(shuō)明即使在輸入?yún)?shù)較多的情況下，本文方法仍能通過(guò)較少的樣本點(diǎn)獲得可靠的輸出參數(shù)區(qū)間。

4 結(jié) 論

本文針對(duì)工程中多體系統(tǒng)的動(dòng)態(tài)不確定性問(wèn)題，提出了一種區(qū)間分析方法。該方法能有效避免由于區(qū)間方法固有的“包裹效應(yīng)”導(dǎo)致的過(guò)估計(jì)問(wèn)題，能夠用較少的樣本得到可靠的輸出參數(shù)區(qū)間，且樣本的數(shù)量對(duì)輸入?yún)?shù)的數(shù)量不敏感。數(shù)值算例表明，本文方法對(duì)典型的機(jī)械系統(tǒng)動(dòng)力學(xué)問(wèn)題均有效，且能應(yīng)用于復(fù)雜的工程實(shí)例。另外，值得注意的是：

(1)本文中計(jì)算輸出參數(shù)的統(tǒng)計(jì)信息僅是為了獲得輸出參數(shù)區(qū)間而進(jìn)行的數(shù)學(xué)處理手段，并非真實(shí)輸出參數(shù)的統(tǒng)計(jì)信息，因此，本質(zhì)上本文的方法還是一種利用統(tǒng)計(jì)信息的區(qū)間方法。

(2)采用本文的方法在樣本數(shù)量很少的情況下，可能估計(jì)到的區(qū)間不能完全包含參數(shù)的真實(shí)區(qū)間，可通過(guò)增加樣本數(shù)量解決，而且增加樣本的過(guò)程不受區(qū)間估計(jì)算法的限制。