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

劉 廣， 劉濟(jì)科， 呂中榮

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

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

1 非線性方程的時(shí)域靈敏度分析

考慮如下形式的動(dòng)力學(xué)方程

(1)

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

(2)

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

(3)

2 反問題

(4)

(5)

(6)

(7)

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

(8)

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

(9)

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

ρcr∈[0.25,0.75]

(10)

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

(11)

和

(12)

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

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

圖1 響應(yīng)靈敏度法的算法實(shí)現(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í)別參數(shù)。根據(jù)本文第二部分的時(shí)域靈敏度分析，可得其靈敏度方程為

(14)

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

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

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

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

(15)

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

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

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

(16)

圖6圖7為測(cè)量數(shù)據(jù)含10%噪聲時(shí)，由初值a0=[0.5,-2,2]得到的識(shí)別曲線。圖6為三個(gè)待識(shí)別參數(shù)隨著每次迭代更新得到的數(shù)值曲線，圖7為對(duì)應(yīng)參數(shù)的相對(duì)誤差。可以看到在初始階段，三個(gè)參數(shù)很快收斂到真值附近，約從第20步開始收斂，每次求得的迭代更新量δa很小，但是結(jié)果逐步逼近真值。其對(duì)應(yīng)的相對(duì)誤差也逐漸減小到0。

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

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

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

(17)

則非線性系統(tǒng)的參數(shù)識(shí)別問題轉(zhuǎn)化為了一個(gè)優(yōu)化問題。圖8為測(cè)量的加速度數(shù)據(jù)含不同等級(jí)噪聲時(shí)，目標(biāo)函數(shù)的進(jìn)化曲線。可以看出，人工蜂群算法差不多在100次迭代左右收斂，2%噪聲時(shí)目標(biāo)函數(shù)在第366次迭代時(shí)還有一次下降。

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

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

表1 測(cè)量噪聲等級(jí)為2%，5%和10%時(shí)，加速度響應(yīng)靈敏度法和人工蜂群算法(ABC)的參數(shù)識(shí)別結(jié)果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%

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

3.2 雙重sine-Gordon非線性方程

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

utt-uxx+f′(u)=0

(18)

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

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

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

(19)

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

(20)

(21)

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

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

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

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

(22)

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

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

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

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

4 結(jié) 論

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

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