吳蓉，劉濟(jì)科，呂中榮，汪利

（中山大學(xué)航空航天學(xué)院，廣東廣州510006）

時(shí)滯系統(tǒng)是一類存在時(shí)滯現(xiàn)象的特殊系統(tǒng)，體現(xiàn)在系統(tǒng)過去某一段時(shí)間的狀態(tài)對系統(tǒng)當(dāng)前狀態(tài)的影響存在一個(gè)時(shí)間上的滯后。時(shí)滯系統(tǒng)普遍存在于各類自然、社會和工程實(shí)際之中，比如：力學(xué)、機(jī)械學(xué)、生態(tài)學(xué)、醫(yī)學(xué)和經(jīng)濟(jì)學(xué)等，對時(shí)滯系統(tǒng)動力學(xué)的研究有利于促進(jìn)這些應(yīng)用領(lǐng)域的發(fā)展［1］。時(shí)滯系統(tǒng)的數(shù)學(xué)模型是時(shí)滯微分方程，不同于由偏微分方程描述的動力學(xué)系統(tǒng)，它是一類泛函微分方程［2］。確定一個(gè)時(shí)滯系統(tǒng)，即要確定時(shí)滯微分方程的參數(shù)。因此，參數(shù)未標(biāo)識或標(biāo)識不完全的情況下，時(shí)滯系統(tǒng)的參數(shù)識別是應(yīng)用和分析時(shí)滯系統(tǒng)的重要準(zhǔn)備工作。

時(shí)滯系統(tǒng)的參數(shù)識別問題是一類復(fù)雜的反問題，與一般微分系統(tǒng)的參數(shù)識別方法思想類似，該問題可轉(zhuǎn)化為一個(gè)非線性最小二乘優(yōu)化問題。基于靈敏度分析的模型修正法是一種常用的參數(shù)識別方法，可較好地與Tikhonov正則化結(jié)合，且便于引進(jìn)修正參數(shù)的置信域區(qū)間，以求解參數(shù)識別方程的非線性最小二乘優(yōu)化問題［3］。本文利用測量動力響應(yīng)數(shù)據(jù)作為測量數(shù)據(jù)，同時(shí)考慮時(shí)滯系統(tǒng)的特性，采用響應(yīng)靈敏度法進(jìn)行了一個(gè)簡單時(shí)滯系統(tǒng)的參數(shù)識別。算例表明，響應(yīng)靈敏度法能較好地識別時(shí)滯系統(tǒng)的參數(shù)且具有一定的抗噪能力。

1 問題描述

1.1 時(shí)滯系統(tǒng)

一般地，時(shí)滯系統(tǒng)由如下的時(shí)滯微分方程給出

其中，x=［x1， x2， …， xn］T表示n 維系統(tǒng)響應(yīng)，τ表示時(shí)滯參數(shù)，ai（i=1，2，…，m）是一般系統(tǒng)參數(shù)，t0代表初始時(shí)刻，x0（t）則是給定初始值時(shí)刻之前的解，為已知函數(shù)。當(dāng)F（·）關(guān)于x(t)，x(t - τ)是線性時(shí)，該系統(tǒng)為線性時(shí)滯系統(tǒng)。

時(shí)滯微分方程相比于一般微分方程具有更復(fù)雜的性質(zhì)。一方面，系統(tǒng)在平衡點(diǎn)附近的線性近似系統(tǒng)的特征方程由一般的有限次多項(xiàng)式代數(shù)方程變?yōu)槌椒匠蹋?］，有無窮多個(gè)特征根，解空間也成為無限維，需要數(shù)值方法進(jìn)行求解；另一方面，要保證線性時(shí)滯微分方程的穩(wěn)定性，方程的所有（無數(shù)個(gè)）特征根全部都要具有負(fù)實(shí)部，這就意味著線性時(shí)滯系統(tǒng)的穩(wěn)定性更為復(fù)雜［5］。本文將使用MATLAB 軟件中的dde23函數(shù)對目標(biāo)時(shí)滯微分方程（1）進(jìn)行數(shù)值求解。

1.2 反問題的定義

2 響應(yīng)靈敏度法識別時(shí)滯系統(tǒng)參數(shù)

2.1 時(shí)滯系統(tǒng)靈敏度分析

對于式（2）的極小值問題，一般采用迭代優(yōu)化的方法進(jìn)行求解。而迭代法的關(guān)鍵在于如何在已有參數(shù)結(jié)果-a的基礎(chǔ)上快速找到合理的更新值δa使得g(-a + δa)盡可能地小。將R?在-a 處進(jìn)行一階泰勒展開并取其線性部分進(jìn)行近似，則非線性目標(biāo)函數(shù)（2） 近似轉(zhuǎn)化為如下線性最小二乘函數(shù)。

而，系統(tǒng)響應(yīng)x（t）關(guān)于時(shí)滯參數(shù)τ 的靈敏度則應(yīng)滿足

2.2 響應(yīng)靈敏度法

式（3）中所得的線性最小二乘目標(biāo)函數(shù)可直接求得最小值，為改善反問題求解的適定性，可采用正則化方法，Tikhonov 正則化［7］是一種常見的正則化方法。

確定合理的置信域半徑η，使得相應(yīng)的更新量δa滿足一致性條件（12），這即為置信域約束。至此，響應(yīng)靈敏度法可列出具體的算法，實(shí)現(xiàn)步驟如表1。

3 數(shù)值算例

考慮一個(gè)機(jī)械元件加工時(shí)常見的顫振方程［4］，是一個(gè)單自由度系統(tǒng)的線性時(shí)滯微分方程：

表1 響應(yīng)靈敏度法的具體步驟Table 1 Procedure of response sensitivity approach

真 實(shí) 參 數(shù) 值a =(-17.6π2，- 0.4π2，1.6 π2，0.5)的系統(tǒng)，噪聲水平分別為0、2%、5%和10%時(shí)所測得的加速度響應(yīng)如圖1所示。

圖1 不同級別測量誤差下測得的加速度響應(yīng)圖Fig.1 Acceleration response of example 1 under different noise level

固定參數(shù)a1=(-17.6π2，- 0.4π，1.6π2，0.5）的系統(tǒng)，分別按式（7）靈敏度方程和差分法計(jì)算時(shí)滯參數(shù)T 的靈敏度，結(jié)果見圖2。兩種方法的計(jì)算結(jié)果十分吻合，證明了時(shí)滯參數(shù)靈敏度方程的正確性。圖2 中，時(shí)滯參數(shù)T 的靈敏度在［0，0.5］的區(qū)間上恒為零，說明要能夠識別時(shí)滯參數(shù)，測量數(shù)據(jù)的采樣時(shí)間長度必須大于時(shí)滯參數(shù)。

圖2 加速度響應(yīng)對時(shí)滯參數(shù)T的靈敏度（T=0.5），F(xiàn)ig. 2 Diagram of sensitivity to delay parameter T

固定測量數(shù)據(jù)誤差水平為2%，考慮4 種初值情況，由響應(yīng)靈敏度法計(jì)算的參數(shù)識別結(jié)果見表2。不同的初始參數(shù)下，參數(shù)識別結(jié)果基本一致，最大相對誤差為1.23%，最大迭代步為102。也就是說，對于該系統(tǒng)響應(yīng)靈敏度法精度高且收斂快。

固 定 初 始 參 數(shù)a1=(-16π2，- 0.3π，π2，0.4)，考慮4 種噪聲水平0、1%、5%和10%，由響應(yīng)靈敏度法計(jì)算的參數(shù)識別結(jié)果見表3。在10%的最大噪聲水平下，最大相對誤差為1.89%，最大迭代步為36。

表2 不同初始值的參數(shù)識別結(jié)果（2%噪聲）Table 2 Identification result of different initial values under noise level of 2%

4 結(jié) 論

基于靈敏度分析，結(jié)合Tikhonov正則化和置信域約束，將時(shí)滯系統(tǒng)的參數(shù)識別反問題轉(zhuǎn)化為一個(gè)標(biāo)準(zhǔn)的非線性最小二乘優(yōu)化問題，提出一種基于響應(yīng)靈敏度法的時(shí)滯系統(tǒng)的參數(shù)識別方法，并對一個(gè)單自由度的線性時(shí)滯系統(tǒng)模型進(jìn)行參數(shù)識別數(shù)值分析，結(jié)果表明：采樣時(shí)長必須大于時(shí)滯參數(shù)才能成功識別時(shí)滯參數(shù)；響應(yīng)靈敏度法能在一定的初值范圍內(nèi)得到準(zhǔn)確的參數(shù)識別結(jié)果；即使在10%的噪聲水平下，響應(yīng)靈敏度法仍能正確地識別時(shí)滯系統(tǒng)的參數(shù)，具有較好的抗噪性；針對該算例，響應(yīng)靈敏度法精度高、收斂快、有較好的魯棒性。