楊智春， 丁允停， 王　樂

(西北工業(yè)大學(xué)結(jié)構(gòu)動(dòng)力學(xué)與控制研究所， 陜西 西安 710072)

?

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

楊智春， 丁允停， 王樂

(西北工業(yè)大學(xué)結(jié)構(gòu)動(dòng)力學(xué)與控制研究所， 陜西 西安 710072)

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

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

引　言

在動(dòng)力學(xué)響應(yīng)分析過程中，系統(tǒng)辨識(shí)起著十分重要的作用，結(jié)構(gòu)動(dòng)力學(xué)系統(tǒng)的物理參數(shù)辨識(shí)一直是結(jié)構(gòu)動(dòng)力學(xué)領(lǐng)域的研究熱點(diǎn)。準(zhǔn)確辨識(shí)結(jié)構(gòu)的質(zhì)量矩陣、阻尼矩陣和剛度矩陣，是準(zhǔn)確地預(yù)計(jì)結(jié)構(gòu)動(dòng)力學(xué)響應(yīng)的前提。

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

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

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

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

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

黏性阻尼結(jié)構(gòu)的運(yùn)動(dòng)方程可以寫為

(1)

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

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

(2)

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

(3)

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

(4)

Padé多項(xiàng)式[9]的一般形式為

(5)

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

不失一般性，為簡明起見，令l=2，即只考慮2個(gè)修正項(xiàng)Ai，即A3和A4，

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

簡記為

(8)

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

(9)

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

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

(10)

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

(11)

(12)

以及

(13)

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

(14)

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

對(duì)非黏性阻尼結(jié)構(gòu)，其阻尼項(xiàng)一般用核函數(shù)的卷積分表示[10]，系統(tǒng)的運(yùn)動(dòng)方程可寫為

(15)

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

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

(16)

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

(17)

(18)

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

(19)

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

(20)

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

(21)

(22)

(23)

(24)

2　數(shù)值仿真算例

2.1算例1

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

圖1　帶黏性阻尼的三自由度質(zhì)量-彈簧系統(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

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

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

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

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

2.2算例2

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

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

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

2.2.1用黏性阻尼模型進(jìn)行辨識(shí)

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

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

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

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

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

2.2.2用非黏性阻尼模型進(jìn)行辨識(shí)

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

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

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

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

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

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

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

(26a)

(26b)

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

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

3　結(jié)　論

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

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

[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]王學(xué)雷, 邵惠鶴. 一種基于Padé近似的頻域辨識(shí)與頻域模型降階新方法[J]. 控制理論與應(yīng)用, 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é近似的時(shí)滯電力系統(tǒng)特征值計(jì)算方法[J]. 電力系統(tǒng)自動(dò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ī)翼的氣動(dòng)彈性響應(yīng)分析[J]. 機(jī)械科學(xué)與技術(shù), 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

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

V214.1

A

1004-4523(2016)01-0024-07

10.16385/j.cnki.issn.1004-4523.2016.01.004

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