楊智春， 丁允停， 王　樂

(西北工業(yè)大學結構動力學與控制研究所， 陜西 西安 710072)

?

用Padé多項式擬合法辨識動力學系統(tǒng)的物理參數(shù)*

楊智春， 丁允停， 王樂

(西北工業(yè)大學結構動力學與控制研究所， 陜西 西安 710072)

提出了一種動力學系統(tǒng)的物理參數(shù)辨識方法。應用Padé多項式對動力學系統(tǒng)的動剛度曲線進行擬合，通過最小二乘法確定Padé多項式中的系數(shù)矩陣，利用遺傳算法對Padé擬合式中的參數(shù)進行優(yōu)化，從而得到系統(tǒng)的質量矩陣、阻尼矩陣和剛度矩陣。數(shù)值算例表明該方法具有較高的辨識精度且適用于黏性阻尼系統(tǒng)和非黏性阻尼系統(tǒng)。

參數(shù)識別； 系統(tǒng)辨識； 結構動力學系統(tǒng)； Padé擬合； 最小二乘法

引　言

在動力學響應分析過程中，系統(tǒng)辨識起著十分重要的作用，結構動力學系統(tǒng)的物理參數(shù)辨識一直是結構動力學領域的研究熱點。準確辨識結構的質量矩陣、阻尼矩陣和剛度矩陣，是準確地預計結構動力學響應的前提。

Phan[1]等利用系統(tǒng)的輸入輸出信號，通過狀態(tài)空間模型辨識系統(tǒng)的質量矩陣、阻尼矩陣和剛度矩陣。Chen和Tsuei[2]同時考慮了黏性阻尼和結構阻尼來對系統(tǒng)的物理參數(shù)進行了辨識。Lee和Kim[3]對Chen和Tsuei的方法進行了改進，將原來方法拓展到多輸入多輸出系統(tǒng)，并在實驗驗證中發(fā)現(xiàn)，Tsuei等人的方法若從動剛度的角度出發(fā)，辨識過程將得到很大簡化，且辨識結果受測量誤差和噪聲的影響較小。但是，正如Lee和Kim[3]在文中所說，利用結構動剛度進行動力學參數(shù)辨識的研究還很少。

廣泛應用于系統(tǒng)降階及參數(shù)擬合的Padé多項式是一種曲線擬合方法。Chazot[4]等將Padé多項式用于黏彈性結構降階，其計算效率與直接計算方法相比，得到很大提高。王學雷[5]提出了一種基于Padé近似的頻域辨識方法，研究了基于積分最小二乘指標的SISO時滯系統(tǒng)頻域辨識問題。葉華[6]等利用Padé多項式來逼近時滯環(huán)節(jié)，提出了一種時滯電力系統(tǒng)特征值的計算方法。Fournodavlos和Nestoridis[7]從數(shù)學角度也研究了Padé在參數(shù)擬合方面的應用。作者[8]在之前的研究中，曾研究過利用Padé多項式對頻域廣義氣動力擬合，得到時域氣動力表達式，進而研究帶遲滯非線性環(huán)節(jié)二元機翼的氣動彈性響應問題。

本文從線性結構動力學系統(tǒng)的動剛度出發(fā)，采用Padé多項式擬合，對動力學系統(tǒng)的物理參數(shù)進行辨識。首先分別從黏性阻尼和非黏性阻尼兩種動力學系統(tǒng)介紹了系統(tǒng)參數(shù)辨識方法，并通過數(shù)值仿真算例對兩種動力學系統(tǒng)的物理參數(shù)進行辨識，驗證了該方法具有較高的辨識精度。

1　系統(tǒng)物理參數(shù)辨識方法

1.1黏性阻尼系統(tǒng)

黏性阻尼結構的運動方程可以寫為

(1)

式中M，C，K分別為質量、阻尼和剛度矩陣，x(t)為位移向量，f(t)為外激勵向量。

假設外激勵為簡諧的，即f(t)=F(ω)ejωt，則響應也是簡諧的，即x(t)=X(ω)ejωt，代入式(1)，得到

(2)

式中X(ω)，F(xiàn)(ω)分別為x(t)，f(t)的Fourier變換。則根據(jù)系統(tǒng)的傳遞函數(shù)的定義，得

(3)

動剛度矩陣為系統(tǒng)傳遞函數(shù)矩陣的逆矩陣

(4)

Padé多項式[9]的一般形式為

(5)

式中Ai+2為修正項，βi共有l(wèi)項，一般事先給出。下面將推導如何擬合出Padé多項式的系數(shù)矩陣Ai，i=0,1,2,…,l+2。

不失一般性，為簡明起見，令l=2，即只考慮2個修正項Ai，即A3和A4，

將式(5)的實部和虛部分離：

簡記為

(8)

式(8)為矛盾方程組，求解過程中，頻率點的選取至少需要2個，其最小二乘解為

(9)

從而得到系數(shù)矩陣Ai，i=0,1,2,…,l+2。其中，(STS)-1ST為S的Moore-Penrose逆，當l>2時，修正項由更多項組成，但系數(shù)矩陣的求解過程與上述相同。

將式(5)表示為復數(shù)形式

(10)

比較式(4)，(10)，根據(jù)復數(shù)相等的充要條件并比較兩式的同類項，得到：

(11)

(12)

以及

(13)

由式(11)，(13)可見，βi的取值會影響參數(shù)辨識的結果，因此，在對動剛度矩陣進行擬合時，需要對βi的取值進行優(yōu)化，即βi值的確定為一個尋優(yōu)過程。本文利用遺傳算法對優(yōu)化變量βi值的選取進行優(yōu)化，優(yōu)化目標為使得重構后的動剛度矩陣與原始的動剛度矩陣在關心的頻率范圍內(nèi)其誤差的范數(shù)最小，其中，重構的動剛度矩陣通過對重構的頻響函數(shù)求逆獲得，目標函數(shù)為

(14)

1.2非黏性阻尼系統(tǒng)

對非黏性阻尼結構，其阻尼項一般用核函數(shù)的卷積分表示[10]，系統(tǒng)的運動方程可寫為

(15)

式中c(t)為核函數(shù)。

令c(t)=C0g(t)，C0為對稱的正定系數(shù)矩陣，g(t)為核函數(shù)的類型。顯然，當g(t)=δ(t)，δ(t)為狄拉克函數(shù)(Dirac delta function)時，式(15)退化為黏性阻尼系統(tǒng)。本文重點研究一種特殊的非黏性阻尼模型，其核函數(shù)為

(16)

對這種指數(shù)型阻尼模型，其核函數(shù)也稱為“松弛函數(shù)”，該阻尼模型廣泛用于表征黏彈性阻尼結構[11]。對式(16)進行傅里葉變換，并乘以系數(shù)矩陣C0，即可得到阻尼矩陣，顯然此時阻尼矩陣為一復數(shù)矩陣，即C=CR+jCI，其中：

(17)

(18)

根據(jù)結構動力學系統(tǒng)的動剛度定義，由式(4)得到黏彈性阻尼系統(tǒng)的動剛度為

(19)

可見，黏彈性阻尼同時也具有剛度效應。將式(5)表示為

(20)

比較式(19)，(20)，由實部、虛部相等的條件以及同冪次項系數(shù)相等，得到：

(21)

(22)

(23)

(24)

2　數(shù)值仿真算例

2.1算例1

如圖1所示的三自由度質量-彈簧系統(tǒng)，假設阻尼為黏性阻尼。

圖1　帶黏性阻尼的三自由度質量-彈簧系統(tǒng)Fig.1　3-DOF mass-spring system with viscous damping

系統(tǒng)的物理參數(shù)矩陣如下：

表1　系統(tǒng)的3階固有頻率

Tab.1　The three natural frequencies of the system

Modenumber123Naturalfrequency/Hz1.03452.96234.3948

識別的頻率修正項接近于零，由前文可知選用黏性阻尼模型進行辨識是合理的。由式(11)，(12)，(13)，得到：

可見，對黏性阻尼系統(tǒng)的物理參數(shù)矩陣進行了準確地辨識。

由式(11)，(12)，(13)，得到：

可見，當模態(tài)不完備時，本文方法對黏性阻尼系統(tǒng)的物理參數(shù)矩陣同樣能夠準確辨識。

2.2算例2

如圖2所示的二自由度質量-彈簧系統(tǒng)，假設阻尼為黏彈性阻尼。

圖2　帶黏彈性阻尼的二自由度質量-彈簧系統(tǒng)Fig.2　2-DOF mass-spring system with viscoelastic damping

系統(tǒng)的物理參數(shù)矩陣為：

2.2.1用黏性阻尼模型進行辨識

此時，系統(tǒng)的阻尼矩陣、質量矩陣和剛度矩陣由式(11)，(12)，(13)得到：

可見，質量矩陣得到準確辨識，但識別得到的系統(tǒng)剛度矩陣不是常數(shù)陣，阻尼矩陣為實數(shù)矩陣，由前文所述可知，選用黏性阻尼模型對該系統(tǒng)進行辨識是不合理的。這里，僅給出在1～100 rad/s頻率帶寬范圍內(nèi)，辨識得到的阻尼矩陣(或剛度矩陣)與原始阻尼矩陣(或剛度矩陣)中的一些元素隨頻率的變化曲線對比，如圖3和4所示。

圖3　用黏性阻尼模型識別的阻尼矩陣Fig.3　Identified damping matrix by viscous damping model

圖4　用黏性阻尼模型識別的剛度矩陣Fig.4　Identified stiffness matrix by viscous damping model

由圖3和4可見，雖然剛度矩陣中的元素K11和K22的最大相對誤差分別為1.99%和3.32%，但已表現(xiàn)出隨頻率變化的特性，而且阻尼矩陣的虛部信息明顯缺失，所以用于辨識的阻尼模型選用黏性阻尼模型是不合理的，應按非黏性阻尼模型進行辨識。

2.2.2用非黏性阻尼模型進行辨識

當辨識阻尼模型選用非黏性阻尼模型時，采用前述針對非黏性阻尼系統(tǒng)的Padé多項式擬合法，對系統(tǒng)的物理參數(shù)矩陣進行辨識，取修正項數(shù)l=2。如圖5所示為對βi的取值優(yōu)化前，取不同βi值得到的辨識結果，其中實線表示的是松弛因子μ取100時的原始阻尼矩陣中的元素隨頻率的變化曲線。顯然，需要按前一節(jié)所述對βi的取值進行優(yōu)化。

由式(21)，(24)可得：

可見，質量矩陣和剛度矩陣得到準確辨識。根據(jù)式(22)，(23)并與式(17)，(18)比較，得到阻尼矩陣為：

可見，松弛因子和系數(shù)矩陣得到了精確地辨識。在1～100 rad/s頻率帶寬范圍內(nèi)，如圖6所示為辨識得到的阻尼矩陣與原始阻尼矩陣的各個元素隨頻率的變化曲線對比(根據(jù)阻尼矩陣對稱性，C21=C12，C22=C11)。顯然，阻尼矩陣的辨識精度也相當高。

圖5　不同βi取值下識別的C22的實部與虛部Fig.5　Real and imaginary part of identified C22 with different values of βi

由式(21)，(24)可得：

可見，質量矩陣和剛度矩陣得到準確辨識。根據(jù)式(22)，(23)并與式(17)，(18)比較，得到阻尼矩陣為：

(26a)

(26b)

可見，當模態(tài)不完備時，本文方法對非黏性阻尼系統(tǒng)的物理參數(shù)矩陣同樣能夠準確辨識。

圖6　用非黏性阻尼模型識別的阻尼矩陣Fig.6　Identified damping matrix by non-viscous damping model

3　結　論

(1) 本文利用Padé多項式對系統(tǒng)的動剛度進行擬合，提出了動力學系統(tǒng)參數(shù)辨識的一種新方法，該方法同時適用于黏性阻尼系統(tǒng)和非黏性阻尼系統(tǒng)。并且，本文方法的辨識結果能夠反映一定的阻尼機理，當頻率修正項較小或接近為零時，用于辨識的阻尼模型應按黏性阻尼模型進行辨識；當頻率修正項較大時，用于辨識的阻尼模型應按非黏性阻尼模型進行辨識。

(2) 本文以Padé多項式修正項中的參數(shù)為變量，求得辨識得到的動剛度矩陣與原始的動剛度矩陣之間的誤差矩陣，以誤差矩陣的范數(shù)為目標函數(shù)，通過遺傳算法對修正項中的參數(shù)進行優(yōu)化，從而提高Padé多項式曲線擬合的精度，辨識得到的質量矩陣、阻尼矩陣和剛度矩陣具有較高的準確度。

[1]Phan M Q, Longman R W. Extracting mass, stiffness, and damping matrices from identified state-space models[C]. AIAA Guidance, Navigation, and Control Conference and Exhibit, 2004, Providence, Rhode Island.

[2]Chen S Y, Ju M S, Tsuei Y G. Estimation of mass, stiffness and damping matrices from frequency response functions[J]. American Society of Mechanical Engineers Journal of Vibration and Acoustics, 1996, 118:78—82.

[3]Lee J H, Kim J. Development and validation of a new experimental method to identify damping matrices of a dynamic system[J]. Journal of Sound and Vibration, 2001, 246(3):505—524.

[4]Chazot J D, Nennig B, Chettah A. Harmonic response computation of viscoelastic multilayered structures using a ZPST shell element[J]. Computers and Structures, 2011, 89:2522—2530.

[5]王學雷, 邵惠鶴. 一種基于Padé近似的頻域辨識與頻域模型降階新方法[J]. 控制理論與應用, 2003, 20(1): 54—58.

Wang Xuelei, Shao Huihe. New method of frequency-domain identification and model reduction based on Pagé approximation[J]. Control theory and Applications, 2003,20(1):54—58.

[6]葉華, 霍健, 劉玉田. 基于Padé近似的時滯電力系統(tǒng)特征值計算方法[J]. 電力系統(tǒng)自動化, 2013, 37(7): 25—30.

Ye Hua, Huo Jian, Liu Yutian. A method for computing eigenvalue of time-delayed power system based on Pagé approximation[J]. Automation of Electric Power Systems, 2013,37(7):25—30.

[7]Fournodavlos G, Nestoridis V. Generic approximation of functions by their Padé approximants[J]. Journal of Mathematical Analysis and Applications,2013, 408:744—750.

[8]谷迎松，楊智春. 帶遲滯非線性環(huán)節(jié)二元機翼的氣動彈性響應分析[J]. 機械科學與技術, 2006, 25(8): 901—904.

Gu Yingsong, Yang Zhichun. Aeroelastic response analysis of an airfoil with hysteresis nonlinearity[J]. Mechanical Science and Technology, 2006,25(8):900—904.

[9]Padé H. Sur la représentation approachée d'une fonction par des fractions rationnelles[J]. Annales de l'Ecole Normale Sup, 1892, 9 (3): 3—93.

[10]Woodhouse J. Linear damping models for structural vibration[J]. Journal of Sound and Vibration, 1998, 215(3): 547—569.

[11]Biot M A. Variational principles in irreversible thermodynamics with application to viscoelasticity[J], Physical Review, 1955, 97(6):1463—1469.

Identifying physical parameters of structural dynamical system using Padé approximation

YANGZhi-chun,DINGYun-ting,WANGLe

(Institute of Structural Dynamics and Control, Northwestern Polytechnical University, Xi′an 710072, China)

A new identification method for the physical parameters of structural dynamical system is proposed. The Padé approximants is used to fit the dynamic stiffness curve of the structural dynamical system, and the coefficient matrices in the Padé polynomial are determined by the least squares method.In addition, genetic algorithms is adopted to optimize the parameters in Padé polynomial. Then the mass, damping and stiffness matrices in the physical space can be extracted from the Padé polynomial. Numerical examples illustrate that the proposed method has good accuracy and is effective for viscous or non-viscous damped systems.

parameters identification; system identification; structural dynamical system; Padé approximants; least squares method

2014-07-17；

2015-06-26

高等學校學科創(chuàng)新引智計劃資助項目(B07050)；國家自然科學基金資助項目(11402205)

V214.1

A

1004-4523(2016)01-0024-07

10.16385/j.cnki.issn.1004-4523.2016.01.004

楊智春(1964—)，男，教授，博士生導師。電話：(029)88460461；E-mail:yangzc@nwpu.edu.cn