基于子系統的非線性時變系統辨識方法*

陳騰飛， 何 歡， 何 成， 陳國平

(1. 南京航空航天大學機械結構力學及控制國家重點實驗室 南京，210016)

(2.南京航空航天大學無人機研究院 南京，210016)

引 言

在時變系統辨識領域，目前對線性系統的研究比較深入?，F有各種時域、頻域方法，大多通過線性系統成熟的理論將線性時不變結構系統的模態分析理論推廣到時變系統[1-3]。但是，幾乎所有的實際工程結構都呈現出非線性特性，在系統辨識中同時考慮系統的非線性和時變特性難度較大[4-5]。近年來，在系統控制、信號處理等領域中，國內外學者已經開始研究，并做出一些有意義的嘗試[6-8]。

根據模型不同，非線性時變系統辨識方法可以分為兩大類：第1類是神經網絡模型方法，第2類是非線性參數模型方法[9]。Ahmed-Ali等[10]基于徑向基函數神經網絡發展提出一種自適應辨識方法。He等[11]基于“短時時不變”假設，將非線性自回歸滑動平均(nonlinear autoregressive moving average, 簡稱NARMA)模型應用于非線性時變系統辨識問題。一些用于線性時變系統的信號處理方法也被應用在非線性問題中。Li等[12]結合B樣條基函數和正交前向回歸(orthogonal forward regression, 簡稱OFR)算法，有效地追蹤線性系統中的慢變及快變參數，并將此方法結合時變NARMA模型，用于非線性時變系統問題，得到了較好的辨識效果[13]。

非線性時變系統辨識包括非線性的定位、非線性類型的確定以及參數估計，其中非線性定位方法已經發展成熟。目前，大多研究默認非線性類型已知，并沒有涉及這方面的討論。因此，迫切需要新的方法確定系統中時變非線性的所有信息，包括具體位置、類型和系數。正交三角分解QR是回歸量正交化常用手段，可以用于消除不同變量之間的耦合關系。本研究將QR分解應用于連續時間MDOF模型中，準確辨識出系統的所有線性及非線性時變特征。

在動力學分析中，實際工程結構通常被離散為多自由度結構系統。根據實際情況，可以采取不同的離散方法，例如集中質量系統，有限單元方法等。文中假設結構可以精確離散化為多自由度系統，并且所有自由度上的運動被精確測量。文中簡單描述非線性時變MDOF系統的數學模型；根據自由度將結構系統分解為不同子系統；結合QR分解和ERR計算，對系統非線性時變參數進行辨識；最后，通過數值算例來驗證辨識方法的精度。

1 非線性時變MDOF系統數學模型

通常，n自由度線性時變系統的運動方程可以寫為

(1)

(2)

由式(2)，系統第i個自由度上的運動方程為

(3)

在MDOF系統中，每個自由度的運動不僅由自由度上的載荷分量決定，還受到其他自由度運動連接的影響，而系統辨識的意義在于確定不同自由度之間這些未知連接關系。用一組基函數對系統中非線性恢復力展開

(4)

在一般的動力學系統中，多項式基函數可以擬合多種類型的非線性恢復力，例如三次非線性剛性力、立方根型非線性剛性力等。將式(4)代入式(3)可得

(5)

(6)

其中：βi,k(t)為時變系數。

(7)

由“短時時不變”假設可知，βtn,i,k即為式(6)中時刻tn下的參數值，即βi,k(tn)=βtn,i,k。若系統位移、速度、加速度以及外載荷均可準確測得，則系統參數辨識問題轉化為標準回歸問題。對不同時刻tn，若時刻點數大于未知參數的數量，則所有線性和非線性參數可通過求解式(7)中方程確定。在此之前，需要解決的一個重要問題是如何確定式(6)中基函數具體形式和總數。對此，文中會提出一個基于QR分解的子系統算法，有效解決此問題。

2 非線性時變MDOF系統的子系統描述

時不變集中質量MDOF系統的質量矩陣是對角矩陣，相似地，時變集中質量MDOF系統的質量矩陣可以表示成以下形式

(8)

(9)

因此，式(5)可以改寫為

(10)

盡管式(10)在形式上與式(5)類似，但是式中所有項均有具體的物理意義，包括連接的位置及形式。例如，k'i,j(t)ξi,j表示mi和mj之間存在線性彈簧連接。因此，文中將采用式(10)中的模型。

同式(6)中原理，式(10)中所有線性和非線性項可以寫成統一形式，式(11)可看作一個包括系統線性和非線性連接黑箱模型

(11)

其中：K為線性和非線性項總數；θi,k(t)為模型中時變系數；φi,k為和mi相關的相對位移和速度的任意形式的函數。

在本研究中，將采用多項式基函數展開，因此式(11)可以改寫為

(12)

引入“短時時不變”假設，系統可以表示為

(tm-Δt/2≤t≤tm+Δt/2)

(13)

其中：質量mtm,i和多項式系數γtm,j,k,p,q在此時間區間內均為常數，且

(14)

在式(13)中，如何確定多項式項成為新的難題，因為不是所有的多項式項都是模型中所需要的，這主要取決于系統中非線性特性，包括位置以及非線性的類型，這些信息正是所缺乏的。因此，引入一個新的方法以確定式(13)中多項式項，并將參數估計簡化為每個時間區間內的標準回歸問題。為方便起見，式(13)可以寫為更加簡潔的形式

(15)

3 MDOF子系統非線性時變參數辨識方法

在式(15)的辨識問題中，通過一種基于QR分解的方法，從所有的多項式項中確定模型所需的重要項，并對所選項的系數進行估計。通常，重要項的確定分為兩個步驟[11]。首先，需要確定多項式的最高階，即式(13)中的L。顯然，階次越高，多項式基數擬合非線性模型的能力越強。但是出于計算成本的考慮，常使用較低階次的多項式基函數，具體可通過比較不同階次下的擬合精度來確定。與之相比，第2個步驟是在確定最高階次的情況下，確定模型所需的重要項，是本節討論的重點。應用QR分解確定重要項的過程：①回歸量正交化，消除變量之間的相關性；②通過ERR確定模型所需的重要項；③估計參數。

式(15)中短時區間[tm-Δt/2,tm+Δt/2]內采樣點總數記為N，則方程可以改寫為如下矩陣方程的形式

F=ΦΘ+Η

(16)

其中

為了將回歸量正交化，對回歸矩陣Φ進行QR分解

Φ=QR

(17)

其中:R為K×K階的上三角矩陣且對角元素為正;Q為N×K階的上三角矩陣。

R矩陣可以作如下分解

R=DA

(18)

其中:D為一個K×K階對角矩陣，其元素為R矩陣對角線元素;A為K×K階單位上三角矩陣，即對角元素均為1。

因此，式(17)可以改寫為

Φ=WA

(19)

將式(17)代入式(16)可得

F=Wg+H

(20)

其中：g=AΘ。

式(20)可以改寫為

(21)

若對F自身作內積運算〈F,F〉，將式(21)代入可得

(22)

式(22)等號兩邊同時除以〈F,F〉可得

(23)

ERRj定義為

(24)

其中：gj=〈wj,F〉/〈wj,wj〉。

從所有候選多項式中確定模型所需的重要項，各項的ERR可以提供一個標準。在ERR計算的每一步中，具有最大EERj的候選項會被作為重要項加入模型。若給定公差常數ρ，如果以下條件滿足，此重要項選擇過程將在第Ks步停止

(25)

這里的ERR計算具有較高的效率，可以用來構造出形如式(15)的數學模型，避免所有的非重要項，有利于計算。

至此，文中提出的系統辨識過程可以總結如下。

1) 分別對MDOF系統不同的離散點施加外載荷并測量系統響應。

2) 確定多項式基函數最高階，得到含有K個候選項的模型，并確定ρ數值。

3) 計算各候選項ERR，具有最大ERR的候選項作為重要項加入式(15)中的模型。

4) 在ERR計算的第k(k≥2)步，在剩余候選項中選擇具有最大ERR的項作為模型的第k項。若條件式(25)滿足，繼續步驟5；否則，令k=k+1，重復此步驟。

5) 完成ERR計算，得到由Ks項構成的模型，系統參數估計如下

(26)

其中：As為一個Ks×Ks單位上三角矩陣;gs向量由gk(k=1,2,...,Ks)構成。

4 算例分析

4.1 3DOF集中質量系統辨識

本節通過一個非線性時變MDOF集中質量系統辨識的算例，說明提出的方法的精度和效率。若輸入、輸出數據全部測量，可以通過上一節描述的方法對形如式(5)的模型進行辨識。

考慮如圖1所示的3自由度系統。系統質量分別為m1=-0.01t+2.5 kg，m1=1 kg和m3=1.5 kg。系統中線性剛度彈簧為k1=k2=k4=1 N/m和k3=-0.03t+4 N/m，線性阻尼為c1=c2=c3=c4=0.1 N·s/m。另外，在m2和m3之間有一個三次非線性剛度knl=3+0.000 5t2N/m3。在這些系統參數中，m1,k3和knl為含時變特性。依次單獨對3個集中質量進行激勵，其中外激勵為用方差σ2=25的零均值平穩高斯白噪聲表示的集中載荷。用四階Runge-Kutta方法分別對系統響應進行求解，采樣頻率為100 Hz。在“短時時不變”假設下，每個時間區間內有1 000個時刻點參與計算，即Δt=10 s且5 s≤tm≤95 s?？紤]到測量噪聲的影響，在系統響應數據中添加3%高斯干擾。當噪聲較大情況下，通過測量加速度積分求速度及位移，會導致誤差。因此，文中假設各自由度響應數據(包括位移、速度及加速度)均可直接準確測量。

首先，需要在較多候選多項式中確定模型的重要項。在本算例中，多項式基函數的最高階L=9，圖1所示的3自由度系統可以由形如式(13)中的子系統模型描述，在每個子系統中，由最高階為9的多項式基函數構成最初的子系統方程。通過ERR計算，在5 s≤tm≤95 s內平均ERR大于ρc=0.005%的多項式將作為此模型的重要項，不同子系統的ERR結果如表1所示。

圖1 3自由度集中質量系統Fig.1 The 3DOF lumped mass system

表1子系統1中所有的項均為線性項，非線性項均在ERR計算中被忽略，這說明作用在m1上的恢復力均為線性恢復力，這與圖1所示一致。子系統1中的辨識結果如圖2所示，包括質量m1和m1相關的所有線性連接(紅色實線表示實際參數時間曲線，黑色點劃線表示辨識出的曲線)。圖2(a)表明質量m1辨識結果準確，包括其時變特性。為了定量分析辨識結果的精確度，將辨識出的曲線最小二乘擬合為時間的線性函數曲線，擬合曲線參數與實際質量m1比較，計算誤差。擬合的曲線在圖2(a)中用藍色虛線表示，其參數及誤差如表2所示。擬合直線與真實曲線接近，斜率和初值的誤差均小于5%。由圖2(b)～(e)可知，對子系統1中的線性剛度和線性阻尼識別精度相當高，4條辨識時間曲線與實際常參數十分接近。

表1 各子系統中ERR結果Tab.1 ERRs of Terms in the Subsystems

圖2 子系統1辨識結果Fig.2 Identification result of subsystem 1

表2 子系統1中m1時間曲線參數Tab.2 Parameters in the Time Expression of m1 in subsystem 1

在表1子系統2中，第2,5項代表m2和m1之間的線性剛度和線性阻尼，第3,6項代表m2和m3之間的線性剛度和線性阻尼，第4項代表m2和m3之間的非線性剛度連接，這與圖1所示的結構特性相吻合。子系統2的辨識結果如圖3所示。在這個子系統中，存在線性時變剛度k3和非線性時變剛度knl。為定量分析上述兩個時變參數辨識結果的精確度，將k3的辨識結果最小二乘擬合為時間的線性函數，將knl的辨識結果擬合為時間的二次函數，擬合項的系數及其誤差如表3所示。擬合曲線與真實曲線很接近，系數誤差均低于5%。

圖3 子系統2辨識結果Fig.3 Identification result of subsystem 2

表3 子系統2中辨識時間曲線參數Tab.3 Identified parameters in the time expression in subsystem 2

由表1子系統3可知，在子系統3中時變參數k3和knl也被確定，分別由第3和第6項代表，子系統3的辨識結果如圖4所示。同樣地，對時變參數k3和knl辨識結果定量分析，結果如表4所示，由圖5(c),(f)可知，用藍色虛線表示的擬合曲線幾乎與參數實際曲線重合，表明辨識精度很高。

圖4 子系統3辨識結果Fig.4 Identification result of subsystem 3

圖5 機械臂結構Fig.5 Mechanical arm structure

表4 子系統3中辨識時間曲線參數Tab.4 Identified parameters in the Time expression of k3 in subsystem 3

4.2 機械臂結構辨識

圖5中機械臂結構運動方程可以寫為

(27)

(28)

時變阻尼矩陣為

C(t)=

(29)

這里包括系數為0.005的比例阻尼。

非線性恢復力矩向量為

(30)

外載荷向量為

(31)

假設兩個移動質量塊位移的表達式為

ri=0.5+0.3sin(πt) (i=1,2)

(32)

兩桿的初始角度為

(φ1,φ2)=(π/4,0)

(33)

分別用方差σ2=25的零均值平穩高斯白噪聲表示的力矩載荷對兩個子系統進行激勵。用四階Runge-Kutta方法分別對系統10 s內響應進行求解，采樣頻率為100 Hz。

表5 桿1子系統中kn l辨識結果Tab.5 Identification Results of kn l

圖6 立方根型非線性扭轉剛度knl辨識結果Fig.6 Identified result of the cube-root stiffness knl

上述兩個數值算例結果表明，筆者提出的方法可以準確地識別出時變MDOF系統中非線性特性，確定非線性位置及類型，并估計出系統參數的時間曲線。不同自由度之間的連接可以在不同子系統中分別辨識。值得注意的是，有些系統的非線性特性取決于外載荷的大小和頻率，為了能在合理范圍內對實際結構系統進行近似，外載荷最好具有寬幅、寬頻，以激發系統的非線性特性。

5 結束語

基于QR分解和ERR計算，發展出一個新的非線性時變系統辨識方法，將MDOF系統分為不同的子系統，在不同子系統中，對非線性及時變連接進行定位、估計。該方法主要優勢在于，無需關于系統的先驗知識，直接在連續時間模型中準確辨識非線性時變參數。此方法需要測量模型中各自由度的響應數據，對于大型工程結構，首先要對模型進行降階處理。