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

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

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

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

1 分?jǐn)?shù)階微分方程的正問題求解

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

(1)

(2)

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

(3)

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

(4)

(5)

當(dāng)1<α<2時(shí)，分?jǐn)?shù)階算子Dαx為以下形式[15]：

(6)

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

(7)

其中di,n的表達(dá)式為

(8)

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

(9)

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

(10)

式中，

類似地，當(dāng)1<α<2時(shí)，式(7)可以離散為

(11)

則tn時(shí)刻的加速度為

(12)

式中，

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

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

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

(13)

2.1 分?jǐn)?shù)階方程的時(shí)域靈敏度分析

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

(14)

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

(15)

2.2 反問題

(16)

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

(17)

(18)

(19)

(20)

(21)

ρcr∈[0.25,0.75]

(22)

(23)

且

(24)

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

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

3 算例

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

(25)

3.1 0<α<1的分?jǐn)?shù)階自由振動(dòng)系統(tǒng)

(26)

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

(27)

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

(28)

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

(29)

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

(30)

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

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

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

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

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

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

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

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

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

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

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

3.2 1<α<2的含外激勵(lì)分?jǐn)?shù)階系統(tǒng)

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

(31)

和第一個(gè)算例不同的是，假設(shè)外激勵(lì)q(t)=sint，同樣假設(shè)系統(tǒng)中的待識(shí)別參數(shù)a=[μ,k,α]，其中參數(shù)取值為μ=0.6、k=1.5，所以a=[0.6,1.5,1.6]。Newmark-β法中的參數(shù)取值和信號(hào)采樣頻率、采樣時(shí)間等都和3.1節(jié)算例相同。根據(jù)第1節(jié)中所介紹的數(shù)值方法，同樣可以求得系統(tǒng)(31)的響應(yīng)。圖5所示為系統(tǒng)(31)的位移響應(yīng)和加速度響應(yīng)。圖6為添加了5%高斯噪聲的模擬測(cè)量數(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)代替。為測(cè)試文中方法對(duì)1<α<2分?jǐn)?shù)階系統(tǒng)的參數(shù)識(shí)別效果，本算例同樣定義4種工況——測(cè)量數(shù)據(jù)為不含噪聲的加速度定義為工況1，對(duì)應(yīng)的無(wú)噪聲位移響應(yīng)數(shù)據(jù)為工況2，測(cè)量加速度和測(cè)量位移數(shù)據(jù)含有5%噪聲的情況分別定義為工況3和4。參數(shù)初值設(shè)定為a0=[2,2,2]。與算例1類似，本算例中的4種工況都可以準(zhǔn)確快速地識(shí)別出包括分?jǐn)?shù)階階次在內(nèi)的所有參數(shù)，詳細(xì)識(shí)別結(jié)果如表3所示。在不考慮測(cè)量噪聲的情況下，識(shí)別結(jié)果幾乎能夠和系統(tǒng)真值相同；在考慮噪聲的情況下，也能保持很高的識(shí)別精度，同時(shí)僅需幾十步迭代就可以收斂。

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

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

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

4 結(jié)語(yǔ)

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

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

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