劉廣 劉濟科 呂中榮

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

在很多復(fù)雜科學(xué)和新興領(lǐng)域(如某些具有遺傳性和記憶性等全局相關(guān)性的領(lǐng)域[1])中，分數(shù)階導(dǎo)數(shù)相較于整數(shù)階導(dǎo)數(shù)而言有著建模準確、所需參數(shù)較少的優(yōu)勢，因此廣泛應(yīng)用于粘彈性材料[2- 3]、機械系統(tǒng)[4- 5]、控制系統(tǒng)[6]和反常擴散[7]等領(lǐng)域。在對這些系統(tǒng)進行建模的過程中，如何準確地確定系統(tǒng)的各種參數(shù)一直是一個棘手的問題，尤其是當系統(tǒng)中包含多個分數(shù)階導(dǎo)數(shù)算子的時候。確定系統(tǒng)參數(shù)的這一類問題可以被歸類為參數(shù)識別問題，所用的方法即為參數(shù)識別方法。

一般來說，參數(shù)識別方法可分為兩類，其中一類是群智能算法，如遺傳算法[8- 9]、人工蜂群體算法[10]等。這些群智能算法的基本思想是模擬自然界生物的群體行為構(gòu)建隨機優(yōu)化算法，將優(yōu)化目標轉(zhuǎn)化為個體對環(huán)境的適應(yīng)性，所以具有全局搜索能力和魯棒性強的優(yōu)點，且不要求原方程具有嚴格的連續(xù)性和可微性。但是，群智能算法也有一些缺點，比如收斂速度慢、每次計算都會得到不同的結(jié)果等。另外一類參數(shù)識別方法是基于梯度的方法，如牛頓法[11- 12]和靈敏度法[13]等。在這類方法中，頻域或時域的靈敏度法是最常用的，因為它們可以避免求解二階導(dǎo)數(shù)這一繁瑣的工作?；陬l域的靈敏度法是從大量的由傳感器得到的數(shù)據(jù)去識別系統(tǒng)的固有頻率、模態(tài)或者其他的系統(tǒng)參數(shù)，但是在實際工程中，可以測量到的模態(tài)參數(shù)總是有限的，而且往往對測量噪聲很敏感。與基于頻域數(shù)據(jù)的方法相比，基于時域的靈敏度法具有明顯的優(yōu)勢，原則上即使只有一個傳感器，只要測量時間足夠長，就能獲得足夠的數(shù)據(jù)。文中嘗試以系統(tǒng)的時域數(shù)據(jù)來識別分數(shù)階微分系統(tǒng)的參數(shù)。

1 分數(shù)階微分方程的正問題求解

考慮如下形式的單自由度線性分數(shù)階系統(tǒng)：

(1)

(2)

式中，β和γ為Newmark-β法的控制參數(shù)，Δt為時間步長。文中的分數(shù)階算子Dαx采用Caputo定義，當分數(shù)階階次0<α<1時，其有以下形式[14]：

(3)

根據(jù)Adams離散法則，式(3)可以離散為[15]

(4)

(5)

當1<α<2時，分數(shù)階算子Dαx為以下形式[15]：

(6)

同樣，根據(jù)Adams離散法則，式(6)可以被離散為

(7)

其中di,n的表達式為

(8)

當0<α<1時，將式(2)和(5)代入式(4)，式(7)變?yōu)?/p>

(9)

由此可以求得tn時刻的加速度如下：

(10)

式中，

類似地，當1<α<2時，式(7)可以離散為

(11)

則tn時刻的加速度為

(12)

式中，

值得指出的是，上文中的推導(dǎo)也同樣適用于線性分數(shù)階多自由度系統(tǒng)，為簡便起見文中從略。

2 基于增強響應(yīng)靈敏度的參數(shù)識別

繼續(xù)選用單自由度系統(tǒng)作為例子來闡述時域靈敏度的分析過程和參數(shù)識別過程?？紤]式(1)有如下初始條件的情況：

(13)

2.1 分數(shù)階方程的時域靈敏度分析

假設(shè)系統(tǒng)(13)中的μ、k和分數(shù)階階次α都是未知參數(shù)，表示為a=(μ,k,α)∈A，其中A是未知參數(shù)的可行域。用ai表示a中的第個未知參數(shù)。值得注意的是，式(13)中的響應(yīng)x是未知參數(shù)a和時間t的函數(shù)，即x=x(t,a)，則x對未知參數(shù)的響應(yīng)靈敏度為

(14)

δ是狀態(tài)變量x的系數(shù)。將式(14)代入式(13)中，靈敏度方程變?yōu)?/p>

(15)

2.2 反問題

(16)

且R(a)是待識別參數(shù)a的函數(shù)。我們的目標是從測量數(shù)據(jù)反過來識別出a=(a1,a2,…,aN)，這樣的反問題一般被表述為一個非線性最小二乘優(yōu)化問題，即找到參數(shù)a∈A，使得如下目標函數(shù)最?。?/p>

(17)

(18)

(19)

(20)

(21)

ρcr∈[0.25,0.75]

(22)

(23)

且

(24)

關(guān)于式(23)和(24)的詳細證明參看文獻[16]。由式(15)可知，增大正則化參數(shù)可以使迭代更新步長足夠小，合理增大正則化參數(shù)可以滿足置信域限制，此時的正則化也稱為增強的正則化。綜上，增強響應(yīng)靈敏度法的算法實現(xiàn)過程如表1所示。

表1 增強響應(yīng)靈敏度法的算法實現(xiàn)過程Table 1 Flowchart of enhanced response sensitivity approach

3 算例

接下來以一個α=0.5的分數(shù)階自由振動系統(tǒng)和α=1.6的含外激勵分數(shù)階系統(tǒng)為例，闡述文中正問題的計算方案和反問題的參數(shù)識別過程。表1中增強響應(yīng)靈敏度法的參數(shù)取值分別為tol=10-10，γ=1.414，ρcr=0.5，Nit=1 000，Ntr=20。所有算例選取的測量數(shù)據(jù)均為加速度響應(yīng)，倘若考慮測量數(shù)據(jù)中的噪聲，噪聲都按照下式進行模擬：

(25)

3.1 0<α<1的分數(shù)階自由振動系統(tǒng)

(26)

可以證明，式(26)有如下形式的解析解[18]：

(27)

在靈敏度分析中，μ和k是方程中某些項的系數(shù)，它們的靈敏度解析表達式可以根據(jù)式(14)獲得且具有如下形式：

(28)

和μ與k不同的是，分數(shù)階階次α的靈敏度是式(1)中Γ函數(shù)的偏微分形式，求解十分繁瑣，所以其靈敏度響應(yīng)用其差分響應(yīng)代替。假設(shè)α有一個攝動量Δ，則系統(tǒng)(26)的攝動系統(tǒng)為

(29)

xα為方程響應(yīng)，可通過第1節(jié)中所述計算格式獲得。因此，α的中心差分響應(yīng)為

(30)

文中所有算例的α的靈敏度都根據(jù)式(29)算得。

圖1是系統(tǒng)(25)不含噪聲的位移響應(yīng)和加速度響應(yīng)?？梢钥闯?，由第1節(jié)中所述的離散格式得到的結(jié)果和解析解吻合非常好，圖2是含有5%高斯噪聲的模擬測量結(jié)果。

圖1 系統(tǒng)(25)數(shù)值解和解析解的響應(yīng)數(shù)據(jù)對比

Fig.1 Comparison between analytical and numerical responses of system(25)

圖2 系統(tǒng)(25)含5%高斯噪聲的測量數(shù)據(jù)

Fig.2 Measured responses containing 5% Gaussian noise of system(25)

在反問題中，測量數(shù)據(jù)可以是位移數(shù)據(jù)或者是加速度數(shù)據(jù)。假設(shè)參數(shù)初值為a0=[0.8,1.2,0.5]，測量數(shù)據(jù)為不含噪聲的加速度定義為工況1，對應(yīng)的無噪聲位移響應(yīng)數(shù)據(jù)為工況2，測量加速度和測量位移數(shù)據(jù)含有5%噪聲的情況分別定義為工況3和4。圖3和4分別是工況1和2中參數(shù)的相對誤差變化曲線?？梢钥闯?，在迭代過程的前10步，參數(shù)變化很劇烈，但是大概在20步之后，每次的迭代更新量δa越來越小，識別的結(jié)果逐漸逼近真實值。在兩個工況中，都是參數(shù)α有最大的相對誤差，但是兩個工況幾乎都可以精準地識別出參數(shù)的真實值。

圖3 工況1中參數(shù)的相對誤差變化曲線Fig.3 Relative error curves of parameters of Case 1

圖4 工況2中參數(shù)的相對誤差變化曲線Fig.4 Relative error curves of parameters of Case 2

工況3和4的識別結(jié)果相較工況1和2來說，精度要低一些,這說明測量噪聲的確會造成識別精度的下降。但是無論是工況3還是工況4，其識別的最大相對誤差也都在1%以內(nèi)，所以文中提出的增強響應(yīng)靈敏度法有著良好的抗噪性。另外還可以看出，文中方法在進行參數(shù)識別時，迭代步數(shù)都只有幾十步，相較于群智能方法而言，計算量大大減少。本算例的詳細識別結(jié)果如表2所示。

表2 增強響應(yīng)靈敏度法對分數(shù)階自由振動系統(tǒng)的參數(shù)識別結(jié)果Table 2 Parameter identification results for fractional-order free vibration system obtained by enhanced response sensitivity approach

3.2 1<α<2的含外激勵分數(shù)階系統(tǒng)

不失一般性，接下來考慮一個如下形式的α=1.6的含外激勵受迫振動系統(tǒng)：

(31)

和第一個算例不同的是，假設(shè)外激勵q(t)=sint，同樣假設(shè)系統(tǒng)中的待識別參數(shù)a=[μ,k,α]，其中參數(shù)取值為μ=0.6、k=1.5，所以a=[0.6,1.5,1.6]。Newmark-β法中的參數(shù)取值和信號采樣頻率、采樣時間等都和3.1節(jié)算例相同。根據(jù)第1節(jié)中所介紹的數(shù)值方法，同樣可以求得系統(tǒng)(31)的響應(yīng)。圖5所示為系統(tǒng)(31)的位移響應(yīng)和加速度響應(yīng)。圖6為添加了5%高斯噪聲的模擬測量數(shù)據(jù)。

圖5 系統(tǒng)(31)的位移響應(yīng)和加速度響應(yīng)Fig.5 Displacement and acceleration responses of system(31)

系統(tǒng)(31)的參數(shù)靈敏度方程同樣可以按照3.1節(jié)算例的方式推導(dǎo)得到，這里不再贅述；α的靈敏度響應(yīng)根據(jù)式(30)以差分響應(yīng)代替。為測試文中方法對1<α<2分數(shù)階系統(tǒng)的參數(shù)識別效果，本算例同樣定義4種工況——測量數(shù)據(jù)為不含噪聲的加速度定義為工況1，對應(yīng)的無噪聲位移響應(yīng)數(shù)據(jù)為工況2，測量加速度和測量位移數(shù)據(jù)含有5%噪聲的情況分別定義為工況3和4。參數(shù)初值設(shè)定為a0=[2,2,2]。與算例1類似，本算例中的4種工況都可以準確快速地識別出包括分數(shù)階階次在內(nèi)的所有參數(shù)，詳細識別結(jié)果如表3所示。在不考慮測量噪聲的情況下，識別結(jié)果幾乎能夠和系統(tǒng)真值相同；在考慮噪聲的情況下，也能保持很高的識別精度，同時僅需幾十步迭代就可以收斂。

圖6 系統(tǒng)(31)含5%高斯噪聲的測量數(shù)據(jù)

Fig.6 Measured responses containing 5% Gaussian noise of system(31)

表3 增強響應(yīng)靈敏度法對分數(shù)階強迫振動系統(tǒng)的參數(shù)識別結(jié)果Table 3 Parameter identification results for fractional-order forced system obtained by enhanced response sensitivity approach

4 結(jié)語

文中基于Adams離散法則和Newmark-β法，提出了一種針對分數(shù)階系統(tǒng)正問題的求解方案。無論是分數(shù)階階次α取值為0<α<1還是1<α<2，都給出了相應(yīng)的離散格式，并用增強響應(yīng)靈敏度法對分數(shù)階微分系統(tǒng)的參數(shù)進行了識別。為了驗證所提出的正問題數(shù)值方法的有效性，闡述反問題的識別過程，文中還以一個α=0.5的分數(shù)階自由振動系統(tǒng)和α=1.6的含外激勵分數(shù)階系統(tǒng)為例進行了說明，得出以下結(jié)論：

(1)無論是位移數(shù)據(jù)還是加速度數(shù)據(jù)，都可以用于增強響應(yīng)靈敏度法的參數(shù)識別，且加速度數(shù)據(jù)的識別精度和識別效率都較位移數(shù)據(jù)要好；

(2)文中提出的識別方法有很強的抗噪性，即便是測量數(shù)據(jù)中含有5%的噪聲，仍舊可以得出很好的識別結(jié)果。因此，文中提出的增強響應(yīng)靈敏度法可以高效準確地識別出分數(shù)階微分系統(tǒng)中的各種參數(shù)。將響應(yīng)靈敏度分析方法推廣到更多的非線性領(lǐng)域，將是下一步的主要工作。