基于增強(qiáng)響應(yīng)靈敏度法的多體系統(tǒng)參數(shù)識別

陳 敏， 劉 廣， 汪 利， 呂中榮

(中山大學(xué) 航空航天學(xué)院，深圳 510275)

自20世紀(jì)60年代以來，人們?nèi)找骊P(guān)注機(jī)械、車輛、機(jī)器人的動(dòng)力學(xué)特性。與此同時(shí)，電子計(jì)算機(jī)、數(shù)值計(jì)算方法的進(jìn)步為這些復(fù)雜系統(tǒng)動(dòng)力學(xué)的建模和分析提供了可能性[1]。在這樣的背景下，人們開始研究由多個(gè)部件和運(yùn)動(dòng)副組成的動(dòng)力學(xué)系統(tǒng)——多體系統(tǒng)。典型的多體系統(tǒng)包括運(yùn)載器[2-3]、航天器[4-5]、機(jī)器人[6-7]、軍事武器[8]、各種機(jī)械裝置[9]等。正是由于多體系統(tǒng)廣泛存在于工業(yè)生產(chǎn)與日常生活中[10]，對于多體系統(tǒng)的模型建立、優(yōu)化設(shè)計(jì)、測試診斷過程等，識別其系統(tǒng)參數(shù)就顯得尤為重要。多體系統(tǒng)的參數(shù)識別主要是根據(jù)系統(tǒng)的響應(yīng)識別出系統(tǒng)的幾何參數(shù)和慣性參數(shù)。從實(shí)際物理模型試驗(yàn)得到的數(shù)據(jù)可以用于識別只是已知其粗略大小的系統(tǒng)參數(shù)。

不同系統(tǒng)的組合、測量數(shù)據(jù)的要求、模型的誤差以及優(yōu)化的方法形成了不同的參數(shù)識別方法。對于多體系統(tǒng)，目前已經(jīng)有很多識別其未知參數(shù)的優(yōu)化方法。傳統(tǒng)的非線性系統(tǒng)的識別方法有Gauss最小二乘法，Young提出的工具變量法，F(xiàn)isher提出的最大似然估計(jì)法。其中，最小二乘法由于簡易以及適用性廣而被廣泛利用。該方法的主要缺陷在于它所形成的是有偏估計(jì)，為了克服這一缺點(diǎn)，Schiehlen等[11]提出了在模型噪聲是平穩(wěn)隨機(jī)過程的假設(shè)下，提出了基于協(xié)方差的參數(shù)識別方法。Serban等[12]利用最小二乘優(yōu)化算法，并且通過直接差分法進(jìn)行靈敏度分析，成功識別了簡單的擺錘模型，最后利用矩陣的條件數(shù)提出了參數(shù)可識別性的數(shù)值測試手段。Bock等[13]發(fā)展了邊值問題方法來估計(jì)微分代數(shù)方程的參數(shù)，主要是通過在有限差分或多重打靶法離散模型時(shí)引入一個(gè)可能的邊界條件，并將其作為最小二乘目標(biāo)函數(shù)的非線性約束。Uchida等[14]提出了一個(gè)新的辦法來直接識別微分代數(shù)方程形式的多體系統(tǒng)方程，主要處理手法是將Lagrange乘子表示成關(guān)于時(shí)間的分段線性函數(shù)代入到動(dòng)力學(xué)方程，通過標(biāo)準(zhǔn)的最小二乘算法就可以同時(shí)識別出系統(tǒng)的慣性參數(shù)和Lagrange乘子。Oberpeilsteiner等[15]利用伴隨靈敏度方法結(jié)合二階伴隨系統(tǒng)識別多體系統(tǒng)的摩擦參數(shù)。他們與Lauss等之后將伴隨靈敏度分析與Fourier分析相結(jié)合，提出了一個(gè)頻域識別方法用于識別振蕩多體系統(tǒng)，表現(xiàn)出很大的潛力和效率[16]。

本文利用Lü等[17]提出的增強(qiáng)響應(yīng)靈敏度法用于識別多體系統(tǒng)問題。增強(qiáng)相應(yīng)靈敏度方法首先通過時(shí)域測量數(shù)據(jù)與理論系統(tǒng)計(jì)算響應(yīng)數(shù)據(jù)的誤差的最小二乘估計(jì)，將參數(shù)識別問題轉(zhuǎn)化為數(shù)學(xué)上典型的反問題，然后通過引入Tikhonov正則化解決該問題可能存在的不適定性，并通過置信域流程選取得到合適的正則化參數(shù)，再利用有限差分獲得的一階參數(shù)靈敏度數(shù)據(jù)，最終通過迭代法求解反問題得到系統(tǒng)的真實(shí)參數(shù)。由此可見，增強(qiáng)相應(yīng)靈敏度方法只需利用到系統(tǒng)的動(dòng)力學(xué)響應(yīng)，在實(shí)際中只要試驗(yàn)測量時(shí)長足夠，響應(yīng)數(shù)據(jù)很容易通過傳感器直接大量獲取，這是其最明顯的優(yōu)勢；再者，增強(qiáng)響應(yīng)靈敏度法屬于梯度法，收斂速度快，并且不像Newton法那樣需要二階靈敏度數(shù)據(jù)，再利用置信域流程，該方法計(jì)算效率更高，且對測量數(shù)據(jù)噪聲不敏感；最后，該方法也已經(jīng)成功地應(yīng)用于振動(dòng)系統(tǒng)[18]、混沌系統(tǒng)[19]、滯回系統(tǒng)[20]等的識別，均表現(xiàn)出了良好的識別精度、抗噪性和穩(wěn)定性。

1 多體系統(tǒng)的動(dòng)力學(xué)建模及求解

利用絕對坐標(biāo)方法可以建立起多體系統(tǒng)動(dòng)力學(xué)方程

(1)

對于位置φ∈Rn-m和速度φ′∈Rn-m的系統(tǒng)的初始條件可以表示為如下形式

(2)

利用傳統(tǒng)增廣法的Baumgarte方法求解方程式(1)。由于式(1)是指標(biāo)3的微分代數(shù)方程，先對式(1)中的約束方程關(guān)于時(shí)間t求二階導(dǎo)，將其化為指標(biāo)1的微分代數(shù)方程，并寫成如下的等價(jià)形式

(3)

其中

由式(3)中后兩式直接增廣計(jì)算，得到

(4)

(5)

由式(3)的第一式和式(4)即構(gòu)成常規(guī)的一階常微分方程

(6)

式(6)可以利用MATLAB軟件的ode進(jìn)行求解。

(7)

從而，式(3)的第三式變?yōu)?/p>

(8)

再利用MATLAB軟件求解即可。

對于阻尼參數(shù)的選取，參考文獻(xiàn)[22]提出了參數(shù)α和β的一種自動(dòng)選取方法，能得到很好的數(shù)值穩(wěn)定性，即

(9)

式中，h為求解微分方程式(3)所用的離散差分格式的步長。

2 多體系統(tǒng)的參數(shù)識別

2.1 反問題的建立

(10)

(11)

(12)

其中各項(xiàng)如下所示

至此，即可得到關(guān)于Δb的一個(gè)線性最小二乘問題，如下所示

(13)

然而，問題式(13)通常是不適定的，因此需要引入Tikhonov正則化進(jìn)行處理，可得

(14)

(15)

2.2 時(shí)域靈敏度分析

基于梯度法解決反問題式(13)，靈敏度分析是不可或缺的。由于直接利用傳統(tǒng)增廣法求解多體系統(tǒng)的靈敏度復(fù)雜且耗時(shí)巨大，故采用有限差分法求取。本文中心差分求取響應(yīng)靈敏度，即

(16)

式中：ε為一個(gè)充分小的正數(shù)；而q(t,bi+ε)和q(t,bi-ε)可以通過求解正問題式(1)得到。本文ε取0.000 1。

2.3 引入置信域

(17)

(18)

從而需要選取合適的δ使得一致性條件式(18)滿足。

(19)

以及

(20)

從上可知，通過適當(dāng)?shù)脑黾诱齽t化參數(shù)μ，可以使迭代更新Δbμ變小，并且滿足置信域條件式(18)。

2.4 增強(qiáng)響應(yīng)靈敏度法

通過以上幾部分的分析，至此，建立起了用于識別多體系統(tǒng)參數(shù)的增強(qiáng)響應(yīng)靈敏度法的實(shí)現(xiàn)過程。具體的算法流程，如表1所示。

表1 增強(qiáng)響應(yīng)靈敏度的算法流程Tab.1 Algorithmic procedure for enhanced response sensitivity approach

3 數(shù)值算例

這個(gè)部分以平面連桿機(jī)構(gòu)和空間連桿機(jī)構(gòu)為例進(jìn)行參數(shù)識別。系統(tǒng)的動(dòng)力學(xué)方程利用MapleSim軟件的Multibody Library建模提取，通過第1章所述的方法利用MATLAB軟件求解。

表1中的各項(xiàng)參數(shù)的取值分別為tol=1×10-6，γ=1.414,ρcr=0.5,Nmax=1 000,Ntr=20。測量數(shù)據(jù)均由MapleSim仿真得到，若考慮噪聲，均以服從標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)噪聲進(jìn)行模擬，表示為如下形式

(21)

3.1 平面曲柄連桿機(jī)構(gòu)

考慮如圖1所示的平面曲柄連桿機(jī)構(gòu)，l1和l2分別是曲柄和連桿的長度，m1,m2和m3分別是曲柄、連桿和滑塊的質(zhì)量。取系統(tǒng)的狀態(tài)變量為q=[x1,y1,θ1,x2,y2,θ2,x3]T。除了考慮該機(jī)構(gòu)受重力作用外，為了充分激勵(lì)參數(shù)靈敏度，在曲柄假設(shè)受到如下的轉(zhuǎn)動(dòng)副摩擦模型

圖1 平面曲柄連桿機(jī)構(gòu)Fig.1 A planar slider-crank system

(22)

從而，系統(tǒng)的廣義質(zhì)量矩陣為

M=diag[m1,m1,J1,m2,m2,J2,m3]

(23)

式中，J1和J2分別為曲柄和連桿的轉(zhuǎn)動(dòng)慣量。

待識別的系統(tǒng)參數(shù)取為b=[J1,J2]T。系統(tǒng)所受的外力矢量為

Q=[0,-m1g,τ(t),0,-m2g,0,0]T

(24)

系統(tǒng)的約束方程為

(25)

相應(yīng)的Jacobian矩陣為

(26)

利用傳統(tǒng)增廣法，結(jié)合MATLAB軟件的ode113求解結(jié)果與MapleSim軟件仿真結(jié)果比較，如圖2所示；其誤差如圖3所示。可以看出，引入Baumgarte穩(wěn)定性理論的傳統(tǒng)增廣法得到的結(jié)果與仿真結(jié)果相比的誤差相當(dāng)小，誤差也沒有發(fā)散的趨勢，說明正問題的求解流程是正確。接下來采用MapleSim軟件的仿真數(shù)據(jù)加上不同水平的噪聲來模擬實(shí)際的測量數(shù)據(jù)，選取連桿的平動(dòng)位移數(shù)據(jù)x2和y2進(jìn)行識別。如圖4所示是噪聲水平為20%的情況下的測量數(shù)據(jù)。取參數(shù)初始值b0=[1,1]T，利用增強(qiáng)響應(yīng)靈敏度法識別系統(tǒng)參數(shù)，結(jié)果如表2所示。20%噪聲水平下識別過程中目標(biāo)函數(shù)的變化情況,如圖5所示。可以看出，增強(qiáng)響應(yīng)靈敏度法識別系統(tǒng)參數(shù)均在二十幾個(gè)迭代步內(nèi)收斂，且即使在20%噪聲污染的情況下識別結(jié)果的誤差在可以接受的范圍內(nèi)，該方法表現(xiàn)出了良好的收斂速度和抗噪性。

(a) 連桿水平方向的位移響應(yīng)

(b) 連桿豎直方向的位移響應(yīng)圖2 平面曲柄連桿機(jī)構(gòu)的MATLAB求解結(jié)果(ode113)與MapleSim仿真結(jié)果Fig.2 The calculation results of calculation(ode113) and the simulation results of MapleSim for planar crank-slider system

(a)

(b)圖3 平面曲柄連桿機(jī)構(gòu)的計(jì)算結(jié)果(ode113)與仿真結(jié)果的誤差Fig.3 The error between the calculation results(ode113) and the simulation results of planar crank-slider system

(a) 受噪聲污染的連桿水平方向的位移響應(yīng)

(b) 受噪聲污染的連桿豎直方向的位移響應(yīng)圖4 20%高斯噪聲污染的測量數(shù)據(jù)Fig.4 Measured response polluted by Gaussian noise at level 20%

圖5 20%高斯噪聲水平下目標(biāo)函數(shù)隨迭代步數(shù)的變化曲線Fig.5 Objective function curves with respect to the iteration steps at 20% Gaussian noise level

表2 平面曲柄連桿機(jī)構(gòu)的參數(shù)識別結(jié)果Tab.2 Parameter identification results of the planar crank-slider system

3.2 空間曲柄連桿機(jī)構(gòu)

圖6 空間曲柄連桿機(jī)構(gòu)Fig.6 Spatial slider-crank system

表3 空間曲柄連桿機(jī)構(gòu)的各項(xiàng)參數(shù)及其取值Tab.3 The parameters of the spatial crank-slider system and their values

取q=[s,θ,α,β]T為廣義坐標(biāo)，由MapleSim軟件建模提取得到系統(tǒng)的質(zhì)量矩陣為

(27)

系統(tǒng)受到的廣義外力矩陣為

(28)

系統(tǒng)的約束矩陣為

(29)

其對應(yīng)的Jacobian矩陣為

(30)

由此，結(jié)合方程式(1)便建立起了空間連桿機(jī)構(gòu)的動(dòng)力學(xué)方程。同樣地，利用傳統(tǒng)增廣法，結(jié)合MATLAB軟件的ode113求解結(jié)果與MapleSim軟件仿真結(jié)果比較，如圖7所示；其誤差如圖8所示，可以看出計(jì)算結(jié)果和仿真結(jié)果相當(dāng)吻合，說明正問題的求解方法是正確的。利用MapleSim軟件在0～10 s內(nèi)的仿真數(shù)據(jù)q作為測量數(shù)據(jù)，采樣率為200 Hz，利用式(21)加入不同水平的高斯噪聲模擬，待識別的參數(shù)取為b=[Ix1,Iy1,Iy2,Iz2]T，參數(shù)初值取為[1,4.2,1,4.2]T，同樣地，利用增強(qiáng)響應(yīng)靈敏度法識別系統(tǒng)參數(shù)，結(jié)果如表4所示。可以看出，增強(qiáng)響應(yīng)靈敏度法識別系統(tǒng)參數(shù)能夠成功識別出系統(tǒng)的慣性參數(shù)，且即使在20%噪聲污染的情況下識別結(jié)果的誤差依舊在可以接受的范圍內(nèi)，同樣地表現(xiàn)出了良好抗噪性。

(a) 滑塊的位移響應(yīng)

(b) 曲柄的位移響應(yīng)

(c) 萬向鉸α角位移響應(yīng)

(d) 萬向鉸β角位移響應(yīng)圖7 空間曲柄連桿機(jī)構(gòu)的計(jì)算(ode113)與仿真結(jié)果對比Fig.7 The comparison between the calculation results(ode113) and the simulation results of spatial crank-slider system

(a)

(b)

(c)

(d)圖8 空間曲柄連桿機(jī)構(gòu)的計(jì)算結(jié)果(ode113)與仿真結(jié)果的誤差Fig.8 The error between the calculation results(ode113) and the simulation results of spatial crank-slider system

表4 不同噪聲水平下空間連桿機(jī)構(gòu)的識別結(jié)果Tab.4 The parameter identification results of spatial crank-slider system on different noise levels

值得注意的是，由迭代步數(shù)可以看出，本算例在不同噪聲水平下的收斂步數(shù)差異較大。對于無噪聲的情況，收斂速度很慢，是因?yàn)橄到y(tǒng)的參數(shù)靈敏度并沒有被外力充分激勵(lì)，導(dǎo)致靈敏度過小，使得收斂速度大大降低。合適的系統(tǒng)激勵(lì)可以通過參考文獻(xiàn)[24-25]等的方法得到，但只適用于系統(tǒng)響應(yīng)可以視作與參數(shù)線性相關(guān)的情況，并不存在適用于一般系統(tǒng)的情況[26]。

有噪聲的情況下，收斂速度較快得益于置信域的引入。實(shí)際上，若不引入置信域，5%噪聲的情況下，一般的靈敏度算法會(huì)在47步循環(huán)后得到結(jié)果[0.145 2,3.956 9,0.359 5,3.423 7]T，此時(shí)算法遠(yuǎn)遠(yuǎn)未滿足循環(huán)終止的容限tol，與真實(shí)參數(shù)的差值將會(huì)越來越大，這是由于缺少了置信域流程的抗噪性，且系統(tǒng)的參數(shù)靈敏度過低，而噪聲污染的影響相比之下很大，從而無法正確地收斂到真實(shí)參數(shù)。

4 結(jié) 論

本文利用了增強(qiáng)靈敏度方法識別多體系統(tǒng)的參數(shù)，通過兩個(gè)算例表現(xiàn)出了良好的識別精度、收斂速度以及抗噪性。然而，實(shí)際實(shí)現(xiàn)過程其實(shí)存在著很多挑戰(zhàn)，多體系統(tǒng)的特性會(huì)帶來很多問題，諸如系統(tǒng)的最優(yōu)激勵(lì)問題，系統(tǒng)參數(shù)的獨(dú)立性和可識別性等，未來的工作將著力解決以上問題，使得該方法能更好地應(yīng)用于多體系統(tǒng)的參數(shù)識別。