基于時域響應靈敏度分析的非線性系統(tǒng)參數(shù)識別

劉 廣， 劉濟科， 呂中榮

(中山大學 航空航天學院 力學系， 廣州 510275)

非線性振動是工程中的常見振動形式[1-2]，而由于測量儀器的精度有限和部分參數(shù)測量難度，非線性系統(tǒng)中存在著大量的不確定參數(shù)[3-4]?？焖?、準確的識別工程中的非線性系統(tǒng)中不確定參數(shù)一直以來都是研究熱點[5-7]。關于此類非線性系統(tǒng)中的參數(shù)識別方法有很多，比如遺傳算法[8-9]、蜂群算法[10-11]等群智能算法，或者是神經網(wǎng)絡[12]等仿生算法。從計算角度看，這類參數(shù)識別問題其實是一個非線性尋優(yōu)問題。基于響應靈敏度分析的有限元模型修正法已廣泛應用于線性結構系統(tǒng)的局部損傷和裂紋參數(shù)等的識別[13-15]，本文將嘗試該方法推廣應用到非線性系統(tǒng)的參數(shù)識別中。對于一般的非線性振動方程，可以通過數(shù)值積分快速得到系統(tǒng)的響應[16]，本文對系統(tǒng)待識別的物理參數(shù)求導，分析得到系統(tǒng)時域響應的靈敏度，然后構造出響應靈敏度矩陣進一步用于參數(shù)識別[17]。文中以Holmes-Duffing系統(tǒng)和雙重sine Gordon系統(tǒng)為例，介紹響應靈敏度分析在非線性系統(tǒng)參數(shù)識別領域中的具體應用。算例表明，響應靈敏度法從不同的初始迭代值出發(fā)，都可以快速得到良好的結果。且進一步分析表明，該方法對工程中的測量噪聲不敏感，具備良好的應用前景。

1 非線性方程的時域靈敏度分析

考慮如下形式的動力學方程

(1)

為方便后續(xù)的靈敏度分析，對待識別的參數(shù)引入如下形式的靈敏度

(2)

則式(1)可以變?yōu)槿缦滦问?/p>

(3)

2 反問題

(4)

(5)

(6)

(7)

而位移靈敏度和加速度靈敏度可以通過式(3)求得。值得注意的是，式(3)求得的靈敏度計算量和原式(1)是一樣的，只需要計算N次。

(8)

式中：λ≥0為正則化參數(shù)，且

(9)

引入置信域限制，是為了確定更合理的正則化參數(shù)λ和迭代更新量δa。分析可知，線性化將會導致原非線性目標函數(shù)與近似目標函數(shù)出現(xiàn)偏差，根據(jù)線性化理論，要想線性化能保持對非線性函數(shù)的很好近似，迭代更新步長要足夠小，衡量一個較小的更新步長是否合適，則需要滿足如下形式的一致性條件

ρcr∈[0.25,0.75]

(10)

置信域限制保證線性化后的目標函數(shù)與原非線性目標函數(shù)足夠接近。

(11)

和

(12)

由式(10)可知，增大正則化參數(shù)可以使迭代更新步長足夠小，合理增大正則化參數(shù)，可以滿足置信域限制；此時的正則化也稱為增強的正則化。

綜上，響應靈敏度法的算法實現(xiàn)過程如圖1所示。

圖1 響應靈敏度法的算法實現(xiàn)過程Fig.1 The process of enhanced response sensitivity approach

3 算 例

3.1 Holmes-Duffing 系統(tǒng)

考慮如下形式的Holmes-Duffing方程

(13)

式中：a=[a1,a2,a3]為待識別參數(shù)。根據(jù)本文第二部分的時域靈敏度分析，可得其靈敏度方程為

(14)

待識別參數(shù)初值設定為a0=[0.5,-2,2]，圖2圖3為該系統(tǒng)正問題的響應曲線。圖2為系統(tǒng)的時程圖，圖3為對應的相圖。從圖中可以看出，在最初的時刻，由于初值的影響，系統(tǒng)并沒有做穩(wěn)定的周期振動，但是隨著時間的積累，系統(tǒng)逐漸過渡到周期振動，在相圖上的反應則是相圖逐漸過渡到一條封閉曲線上。

圖2 無噪聲時Holmes-Duffing系統(tǒng)的時程圖Fig.2 The time-displacement diagram of the Holmes-Duffing system without noise

圖3 無噪聲時Holmes-Duffing系統(tǒng)的相圖Fig.3 The phase diagram of the Holmes-Duffing system without noise

圖4圖5分別為反問題中的測量數(shù)據(jù)，采樣頻率均為100 Hz，圖4為測量數(shù)據(jù)含2%高斯白噪聲的相圖，圖4為含5%噪聲系統(tǒng)的相圖。本文中的含噪聲的測量數(shù)據(jù)通過式(15)模擬得到

(15)

圖4 測量數(shù)據(jù)含2%噪聲時Holmes-Duffing系統(tǒng)的相圖Fig.4 Measured data with noise level 2%, the phase diagram of the Holmes-Duffing system

圖5 測量數(shù)據(jù)含5%噪聲時Holmes-Duffing系統(tǒng)的相圖Fig.5 Measured data with noise level 5%, the phase diagram of the Holmes-Duffing system

本算例中采用的測量數(shù)據(jù)是加速度，對應的是加速度靈敏度矩陣。根據(jù)式(14)可得系統(tǒng)加速度靈敏度矩陣式(16)所示，代入每一步得到的a然后用四階Runge-Kutta法可得系統(tǒng)加速度靈敏度矩陣，按照本文第二部分流程，可最終得到識別結果。

(16)

圖6圖7為測量數(shù)據(jù)含10%噪聲時，由初值a0=[0.5,-2,2]得到的識別曲線。圖6為三個待識別參數(shù)隨著每次迭代更新得到的數(shù)值曲線，圖7為對應參數(shù)的相對誤差?？梢钥吹皆诔跏茧A段，三個參數(shù)很快收斂到真值附近，約從第20步開始收斂，每次求得的迭代更新量δa很小，但是結果逐步逼近真值。其對應的相對誤差也逐漸減小到0。

圖6 測量數(shù)據(jù)含10%噪聲時Holmes-Duffing系統(tǒng)參數(shù)識別結果Fig.6 Measured data with noise level 10%, the value of parameters in Holmes-Duffing system

圖7 測量數(shù)據(jù)含10%噪聲時Holmes-Duffing系統(tǒng)參數(shù)的相對誤差Fig.7 Measured data with noise level 10%, the relative error of parameters in Holmes-Duffing system

下面用人工蜂群算法(Artificial Bee Colony Algorithm, ABC)對系統(tǒng)式(13)進行識別并作對比。人工蜂群算法是模仿蜜蜂采蜜行為的一種元啟發(fā)式算法，和增強的響應靈敏度一樣，測量數(shù)據(jù)采用加速度數(shù)據(jù)。人工蜂群算法的核心是目標函數(shù)的建立，本算例中利用測量的加速度數(shù)據(jù)和計算的加速度響應誤差來建立目標函數(shù)，表達式為

(17)

則非線性系統(tǒng)的參數(shù)識別問題轉化為了一個優(yōu)化問題。圖8為測量的加速度數(shù)據(jù)含不同等級噪聲時，目標函數(shù)的進化曲線?？梢钥闯?，人工蜂群算法差不多在100次迭代左右收斂，2%噪聲時目標函數(shù)在第366次迭代時還有一次下降。

表1為不同噪聲等級下加速度響應靈敏度法和人工蜂群算法的識別結果。可以看到，測量噪聲無論是2%，5%還是10%，兩種算法都可以得到結果，但是不同噪聲等級下，本文方法識別到的參數(shù)最大相對誤差也低于1%，而人工蜂群算法的相對誤差則基本比本文方法大一個量級，且本文方法的所用時間遠小于人工蜂群算法?？梢灶A見，倘若非線性系統(tǒng)進一步復雜，本文方法還將將進一步節(jié)省計算資源。

圖8 測量數(shù)據(jù)含10%噪聲時，用蜂群算法(ABC)識別的目標函數(shù)變化曲線Fig.8 Measured data with level 10%, the curves of objective function by artificial bee colony algorithm(ABC)

表1 測量噪聲等級為2%，5%和10%時，加速度響應靈敏度法和人工蜂群算法(ABC)的參數(shù)識別結果Tab.1 The results of identification by acceleration response sensitivity approach and artificial bee colony algorithm (ABC), meastred data with noise level 2%,5% and 10%

接下來進一步研究不同的迭代初值a0對計算過程的影響。一般來說，不同的迭代初值會對計算的結果產生一定的影響，也會影響到計算的時間(迭代步數(shù)不同)，好的初值可以很快收斂到精確解。選取如下三組迭代初值進行計算，a0=[0.5,-2,2]，a0=[0.3,-20,40]和a0=[0.5,-8,60]。圖9為測量數(shù)據(jù)含有10%噪聲時的識別結果。圖中x坐標為a1，y坐標為a2，z坐標為a3。從圖9中可以看到，無論是哪組初值，最終都可以收斂到正確的結果；證明了響應靈敏度法有著良好的初值適應性。

3.2 雙重sine-Gordon非線性方程

考慮如下形式的雙重sine Gordon非線性方程，其在在許多物理工程中有廣泛應用[18]

utt-uxx+f′(u)=0

(18)

圖9 測量數(shù)據(jù)含10%噪聲時系統(tǒng)中待識別參數(shù)由不同初值得到的迭代曲線Fig.9 Measured data with noise level 10%, the iterative curve of parameters begin with different initial values

α,β參數(shù)為任意正數(shù)。假設

u=u(θ),θ=kx-ωt

(19)

式中：u為周期波函數(shù)；θ為自變量；k為波數(shù)；x為空間坐標；ω為頻率；t為時間。將式(19)代入式(18)，有

(20)

(21)

圖10圖11為無噪聲時和含5%測量噪聲時雙重sine Gordon非線性系統(tǒng)的振動圖像。從時程圖中可以看出，系統(tǒng)一直做周期振動，且振幅的最小值即為U0。 由于雙重sine Gordon系統(tǒng)的特殊性，該方程無論如何選取初值，系統(tǒng)都會以初值為振幅的最低點做周期振動。同理，系統(tǒng)的相圖是一個封閉的環(huán)。且加入噪聲后的振動圖像變?yōu)檎劬€圖。

圖10 無噪聲時雙重sine-Gordon系統(tǒng)的時程圖Fig.10 The time-displacement diagram of the double sine-Gordon system without noise

圖11 含5%測量噪聲時雙重sine-Gordon系統(tǒng)的時程圖Fig.11 The time-displacement diagram of the double sine-Gordon system with noise level 10%

本算例中，測量數(shù)據(jù)仍舊假定為加速度，由式(21)可得系統(tǒng)加速度靈敏度矩陣如式(22)所示，代入每一步計算得到的a，然后用四階Runge-Kutta法對式(22)進行數(shù)值積分，可得系統(tǒng)加速度靈敏度矩陣用于反問題的識別當中。

(22)

如圖12所示，為測量數(shù)據(jù)含5%噪聲時，待識別參數(shù)α和β由初值a0=[90,8]出發(fā)的識別結果，約迭代20次左右結果開始收斂。從相對誤差曲線可以看到，隨著迭代步數(shù)的增加，相對誤差逐漸趨于0。

圖12 測量數(shù)據(jù)含5%噪聲時，由初值a0=[90,8]出發(fā)系統(tǒng)參數(shù)的相對誤差變化曲線Fig.12 Measured data with noise level 5%, the relative error of parameters about system with initial values a0=[90,8]

圖13和圖14分別為為測量數(shù)據(jù)含5%噪聲時，α和β由初值a0=[85,15]和a0=[75,20]出發(fā)的識別結果。從相對誤差變化曲線可以看出，即便是和真值相差較大的初值，本文算法也可以得到較好的結果，具有較強的初值適應性。

圖13 測量數(shù)據(jù)含5%噪聲時，由初值a0=[85,15]出發(fā)系統(tǒng)參數(shù)的相對誤差變化曲線Fig.13 Measured data with noise level 5%, the relative error of parameters about system with initial values a0=[85,15]

4 結 論

本文算例表明，基于時域的響應靈敏度分析在非線性參數(shù)識別問題上有著快速、準確的優(yōu)點。且進一步分析表明，該方法在非線性系統(tǒng)的參數(shù)識別上，有著對測量噪聲不敏感，不依賴于初始迭代值的特點。和人工蜂群算法相比，其計算量大為減少，僅需幾十步即可收斂到正確結果。推廣響應靈敏度分析方法到更多的非線性領域，將是下一步的主要工作。

圖14 測量數(shù)據(jù)含5%噪聲時，由初值a0=[75,20]出發(fā)系統(tǒng)參數(shù)的相對誤差變化曲線Fig.14 Measured data with noise level 5%, the relative error of parameters about system with initial values a0=[75,20]