基于子系統(tǒng)的非線性時(shí)變系統(tǒng)辨識(shí)方法*

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

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

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

引 言

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

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

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

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

1 非線性時(shí)變MDOF系統(tǒng)數(shù)學(xué)模型

通常，n自由度線性時(shí)變系統(tǒng)的運(yùn)動(dòng)方程可以寫為

(1)

(2)

由式(2)，系統(tǒng)第i個(gè)自由度上的運(yùn)動(dòng)方程為

(3)

在MDOF系統(tǒng)中，每個(gè)自由度的運(yùn)動(dòng)不僅由自由度上的載荷分量決定，還受到其他自由度運(yùn)動(dòng)連接的影響，而系統(tǒng)辨識(shí)的意義在于確定不同自由度之間這些未知連接關(guān)系。用一組基函數(shù)對系統(tǒng)中非線性恢復(fù)力展開

(4)

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

(5)

(6)

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

(7)

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

2 非線性時(shí)變MDOF系統(tǒng)的子系統(tǒng)描述

時(shí)不變集中質(zhì)量MDOF系統(tǒng)的質(zhì)量矩陣是對角矩陣，相似地，時(shí)變集中質(zhì)量MDOF系統(tǒng)的質(zhì)量矩陣可以表示成以下形式

(8)

(9)

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

(10)

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

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

(11)

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

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

(12)

引入“短時(shí)時(shí)不變”假設(shè)，系統(tǒng)可以表示為

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

(13)

其中：質(zhì)量mtm,i和多項(xiàng)式系數(shù)γtm,j,k,p,q在此時(shí)間區(qū)間內(nèi)均為常數(shù)，且

(14)

在式(13)中，如何確定多項(xiàng)式項(xiàng)成為新的難題，因?yàn)椴皇撬械亩囗?xiàng)式項(xiàng)都是模型中所需要的，這主要取決于系統(tǒng)中非線性特性，包括位置以及非線性的類型，這些信息正是所缺乏的。因此，引入一個(gè)新的方法以確定式(13)中多項(xiàng)式項(xiàng)，并將參數(shù)估計(jì)簡化為每個(gè)時(shí)間區(qū)間內(nèi)的標(biāo)準(zhǔn)回歸問題。為方便起見，式(13)可以寫為更加簡潔的形式

(15)

3 MDOF子系統(tǒng)非線性時(shí)變參數(shù)辨識(shí)方法

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

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

F=ΦΘ+Η

(16)

其中

為了將回歸量正交化，對回歸矩陣Φ進(jìn)行QR分解

Φ=QR

(17)

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

R矩陣可以作如下分解

R=DA

(18)

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

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

Φ=WA

(19)

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

F=Wg+H

(20)

其中：g=AΘ。

式(20)可以改寫為

(21)

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

(22)

式(22)等號(hào)兩邊同時(shí)除以〈F,F〉可得

(23)

ERRj定義為

(24)

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

從所有候選多項(xiàng)式中確定模型所需的重要項(xiàng)，各項(xiàng)的ERR可以提供一個(gè)標(biāo)準(zhǔn)。在ERR計(jì)算的每一步中，具有最大EERj的候選項(xiàng)會(huì)被作為重要項(xiàng)加入模型。若給定公差常數(shù)ρ，如果以下條件滿足，此重要項(xiàng)選擇過程將在第Ks步停止

(25)

這里的ERR計(jì)算具有較高的效率，可以用來構(gòu)造出形如式(15)的數(shù)學(xué)模型，避免所有的非重要項(xiàng)，有利于計(jì)算。

至此，文中提出的系統(tǒng)辨識(shí)過程可以總結(jié)如下。

1) 分別對MDOF系統(tǒng)不同的離散點(diǎn)施加外載荷并測量系統(tǒng)響應(yīng)。

2) 確定多項(xiàng)式基函數(shù)最高階，得到含有K個(gè)候選項(xiàng)的模型，并確定ρ數(shù)值。

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

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

5) 完成ERR計(jì)算，得到由Ks項(xiàng)構(gòu)成的模型，系統(tǒng)參數(shù)估計(jì)如下

(26)

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

4 算例分析

4.1 3DOF集中質(zhì)量系統(tǒng)辨識(shí)

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

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

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

圖1 3自由度集中質(zhì)量系統(tǒng)Fig.1 The 3DOF lumped mass system

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

表1 各子系統(tǒng)中ERR結(jié)果Tab.1 ERRs of Terms in the Subsystems

圖2 子系統(tǒng)1辨識(shí)結(jié)果Fig.2 Identification result of subsystem 1

表2 子系統(tǒng)1中m1時(shí)間曲線參數(shù)Tab.2 Parameters in the Time Expression of m1 in subsystem 1

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

圖3 子系統(tǒng)2辨識(shí)結(jié)果Fig.3 Identification result of subsystem 2

表3 子系統(tǒng)2中辨識(shí)時(shí)間曲線參數(shù)Tab.3 Identified parameters in the time expression in subsystem 2

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

圖4 子系統(tǒng)3辨識(shí)結(jié)果Fig.4 Identification result of subsystem 3

圖5 機(jī)械臂結(jié)構(gòu)Fig.5 Mechanical arm structure

表4 子系統(tǒng)3中辨識(shí)時(shí)間曲線參數(shù)Tab.4 Identified parameters in the Time expression of k3 in subsystem 3

4.2 機(jī)械臂結(jié)構(gòu)辨識(shí)

圖5中機(jī)械臂結(jié)構(gòu)運(yùn)動(dòng)方程可以寫為

(27)

(28)

時(shí)變阻尼矩陣為

C(t)=

(29)

這里包括系數(shù)為0.005的比例阻尼。

非線性恢復(fù)力矩向量為

(30)

外載荷向量為

(31)

假設(shè)兩個(gè)移動(dòng)質(zhì)量塊位移的表達(dá)式為

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

(32)

兩桿的初始角度為

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

(33)

分別用方差σ2=25的零均值平穩(wěn)高斯白噪聲表示的力矩載荷對兩個(gè)子系統(tǒng)進(jìn)行激勵(lì)。用四階Runge-Kutta方法分別對系統(tǒng)10 s內(nèi)響應(yīng)進(jìn)行求解，采樣頻率為100 Hz。

表5 桿1子系統(tǒng)中kn l辨識(shí)結(jié)果Tab.5 Identification Results of kn l

圖6 立方根型非線性扭轉(zhuǎn)剛度knl辨識(shí)結(jié)果Fig.6 Identified result of the cube-root stiffness knl

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

5 結(jié)束語

基于QR分解和ERR計(jì)算，發(fā)展出一個(gè)新的非線性時(shí)變系統(tǒng)辨識(shí)方法，將MDOF系統(tǒng)分為不同的子系統(tǒng)，在不同子系統(tǒng)中，對非線性及時(shí)變連接進(jìn)行定位、估計(jì)。該方法主要優(yōu)勢在于，無需關(guān)于系統(tǒng)的先驗(yàn)知識(shí)，直接在連續(xù)時(shí)間模型中準(zhǔn)確辨識(shí)非線性時(shí)變參數(shù)。此方法需要測量模型中各自由度的響應(yīng)數(shù)據(jù)，對于大型工程結(jié)構(gòu)，首先要對模型進(jìn)行降階處理。