，

(南京航空航天大學(xué)機(jī)電學(xué)院，江蘇 南京210016)

基于擴(kuò)展卡爾曼濾波的磁懸浮軸承剛度阻尼辨識(shí)

鄒玥，周瑾

(南京航空航天大學(xué)機(jī)電學(xué)院，江蘇 南京210016)

提出一種磁懸浮軸承剛度阻尼的測(cè)量方法，通過(guò)對(duì)轉(zhuǎn)子施加不平衡質(zhì)量激勵(lì)，利用擴(kuò)展卡爾曼濾波算法識(shí)別出磁懸浮軸承的剛度和阻尼。仿真研究表明，此方法可有效測(cè)量磁懸浮軸承剛度阻尼。經(jīng)過(guò)試驗(yàn)辨識(shí)出使用的磁懸浮試驗(yàn)臺(tái)的支承參數(shù)，并將辨識(shí)結(jié)果帶回到仿真文件中進(jìn)行驗(yàn)證對(duì)比，結(jié)果表明，可有效辨識(shí)出磁懸浮軸承的剛度阻尼參數(shù)。

擴(kuò)展卡爾曼濾波；剛度；阻尼

0 引言

磁懸浮軸承的剛度和阻尼系數(shù)影響著控制系統(tǒng)參數(shù)優(yōu)化、臨界轉(zhuǎn)速預(yù)估和不平衡響應(yīng)計(jì)算等。因此,磁懸浮軸承的剛度和阻尼系數(shù)的辨識(shí)非常重要。針對(duì)五自由度、高速運(yùn)轉(zhuǎn)情況的下磁懸浮軸承剛度、阻尼測(cè)試與辨識(shí)還需完善。采用擴(kuò)展卡爾曼濾波方法辨識(shí)磁懸浮軸承的剛度阻尼，使用不平衡量產(chǎn)生的簡(jiǎn)諧力作為激振力，基于有限元模型的仿真和實(shí)驗(yàn)，對(duì)磁懸浮軸承的支承參數(shù)進(jìn)行辨識(shí)。

1 擴(kuò)展卡爾曼濾波算法基本理論

磁懸浮軸承-轉(zhuǎn)子系統(tǒng)試驗(yàn)臺(tái)如圖1所示。l1，l2為左右兩側(cè)磁懸浮軸承分別離轉(zhuǎn)子質(zhì)心的距離；d1,d2為左右不平衡質(zhì)量分別離質(zhì)心的距離；m1，m2分別為添加在左右兩邊圓盤(pán)上的不平衡質(zhì)量大??；r1，r2為不平衡質(zhì)量施加在2個(gè)圓盤(pán)的半徑；p1,p2為位移傳感器距轉(zhuǎn)子質(zhì)心的距離；用彈簧阻尼系統(tǒng)模擬磁懸浮軸承的支承情況。

對(duì)于磁懸浮軸承轉(zhuǎn)子系統(tǒng)，考慮運(yùn)動(dòng)方程為：

(1)

圖1 磁懸浮軸承-轉(zhuǎn)子系統(tǒng)試驗(yàn)臺(tái)

圖2 轉(zhuǎn)子不平衡質(zhì)量在轉(zhuǎn)子施加位置

z為轉(zhuǎn)子的振動(dòng)位移信號(hào)，z={zx1,zy1,zx2,zy2}T，因耦合參數(shù)對(duì)系統(tǒng)影響較小，且為了計(jì)算的簡(jiǎn)化，忽略耦合參數(shù)，則剛度阻尼矩陣為：

(2)

(3)

kxxi,cxxi(i=1,2)為徑向軸承X方向的剛度阻尼；kyyi，cyyi(i=1,2)為徑向軸承Y方向的剛度。

轉(zhuǎn)子受不平衡質(zhì)量激勵(lì)如圖2所示，與參考軸夾角分別為φ1和φ2的位置，則運(yùn)動(dòng)方程右側(cè)的F0可表示為：

(4)

ω為轉(zhuǎn)子旋轉(zhuǎn)頻率；Ω為轉(zhuǎn)子旋轉(zhuǎn)角速度，Ω=2πω；轉(zhuǎn)子軸承處的位移為zxi,zyi(i=1,2)。

在系統(tǒng)與控制理論中，常常用狀態(tài)變量的形式來(lái)描述一個(gè)系統(tǒng)的運(yùn)動(dòng)。因此，可以用它來(lái)描述轉(zhuǎn)子系統(tǒng)的運(yùn)動(dòng)。擴(kuò)展卡爾曼濾波算法是非線性狀態(tài)濾波的遞推算法，可以將系統(tǒng)的狀態(tài)和未知參數(shù)同時(shí)辨識(shí)出來(lái)[1-3]，需將軸承參數(shù)作為狀態(tài)變量的一部分。在擴(kuò)展卡爾曼濾波算法中就是把轉(zhuǎn)子的支承參數(shù)也看作狀態(tài)向量，對(duì)于時(shí)不變的系統(tǒng)，該向量就是一個(gè)常數(shù)向量。將其與系統(tǒng)的真實(shí)狀態(tài)向量s組合，構(gòu)成增廣狀態(tài)變量，然后利用擴(kuò)展卡爾曼濾波器對(duì)這個(gè)增廣狀態(tài)進(jìn)行最優(yōu)或者次優(yōu)估計(jì)[4]。那么，可以設(shè)全局狀態(tài)變量為：

(5)

r={kxx1,kyy1,kxx2,kyy2,cxx1,cyy1,cxx2,cyy2}T。

由于剛度阻尼參數(shù)未知，可將它們看成系統(tǒng)的另外8個(gè)狀態(tài)，這樣全部16個(gè)狀態(tài)如下：

si為第i個(gè)自由度相對(duì)于基礎(chǔ)的位移；

si+4為第i個(gè)自由度相對(duì)于基礎(chǔ)的速度；

si+8為第i個(gè)自由度剛度系數(shù)

si+12為第i個(gè)自由度阻尼系數(shù)

i=1,2,3,4，系統(tǒng)的一階微分方程組為：

=f(s,F,t)+w

(6)

f(s,F,t)中不為0的元素可以通過(guò)運(yùn)動(dòng)方程(1)轉(zhuǎn)化得到。 由位移傳感器測(cè)得并轉(zhuǎn)換至磁懸浮軸承處的位移向量，即系統(tǒng)的輸出量為Yk，Yk向量中也包含了測(cè)量誤差，那么Yk與位移的實(shí)際值szk之間存在一個(gè)線性關(guān)系，即

Yk=Hszk+vk

(7)

Yk中的變量與各軸承處的位移有以下關(guān)系，則有觀測(cè)方程為：

(8)

狀態(tài)方程和觀測(cè)方程都用離散時(shí)間的形式為：

s(k+1)=f[s(k),F0,t]·τ+s(k)+w

(9)

Y(k)=h[s(k),F0,t]+v(k)

(10)

s(k)為16維狀態(tài)向量；Y(k)為4維軸承處位移；w(k),v(k)為零均值的正態(tài)白噪聲序列;w(k)表示系統(tǒng)的描述誤差;v(k)是測(cè)量中的噪聲[5]。

預(yù)測(cè)階段為[6]：

u(tk),tk]·τ

(11)

(12)

增益陣的計(jì)算為：

(13)

濾波階段為：

(14)

P(tk+1|tk+1)=(I-Kk+1Hk+1)·P(tk+1|tk)

(15)

τ為采樣時(shí)間間隔；Φ為轉(zhuǎn)移矩陣。

(16)

Hk+1=H

(17)

(18)

(19)

估計(jì)準(zhǔn)則是使?fàn)顟B(tài)向量的誤差方差最小[7]。根據(jù)式(11)～式(19)，只要給出系統(tǒng)參數(shù)及狀態(tài)初始值的估計(jì)及初始誤差的協(xié)方差陣，就可以對(duì)量測(cè)數(shù)據(jù)進(jìn)行遞推處理和計(jì)算。實(shí)際表明，對(duì)初始值的要求是非常低的，只要給出大約的量級(jí)就可以。這對(duì)實(shí)際應(yīng)用是十分有利的，采用全局迭代的離線遞推辨識(shí)算法，可以在有限數(shù)據(jù)量測(cè)的條件下能得到較高的參數(shù)估計(jì)精度，將預(yù)測(cè)階段與濾波階段反復(fù)進(jìn)行后，狀態(tài)變量中的支承參數(shù)會(huì)逐漸逼近真實(shí)值。

2 磁懸浮軸承剛度和阻尼的仿真研究

建立試驗(yàn)所用磁懸浮軸承-轉(zhuǎn)子系統(tǒng)有限元模型，附加不平衡質(zhì)量，計(jì)算出磁懸浮軸承處的位移響應(yīng)。根據(jù)上述理論，可辨識(shí)出磁懸浮軸承的剛度和阻尼系數(shù)。有限元仿真時(shí)，先假定左側(cè)磁懸浮軸承的剛度阻尼系數(shù)分別為：

K1=5×106N/m，C1=5 000N·s/m

右側(cè)磁懸浮軸承的剛度阻尼系數(shù)分別為：

K2=6×106N/m，C2=6 000N·s/m

辨識(shí)出磁懸浮軸承的剛度和阻尼系數(shù)隨旋轉(zhuǎn)頻率變化的曲線如圖3和圖4所示。

圖3 20～200Hz剛度變化仿真值曲線

圖4 20～200Hz阻尼變化仿真值曲線

由圖可見(jiàn)參數(shù)識(shí)別精度較高，辨識(shí)結(jié)果在200Hz以內(nèi)與原值誤差較低，與初始給定值極為相符，辨識(shí)結(jié)果較為可靠。隨著轉(zhuǎn)速的升高，辨識(shí)值與真實(shí)值之間的誤差也越來(lái)越大，主要考慮為隨著轉(zhuǎn)子轉(zhuǎn)速的升高，轉(zhuǎn)子的柔性逐漸顯現(xiàn)出來(lái)，采用的動(dòng)力學(xué)方程不能完全表達(dá)其運(yùn)動(dòng)規(guī)律。但在辨識(shí)范圍內(nèi)，即20～200Hz之間，仿真計(jì)算的誤差基本處于10%以下，辨識(shí)結(jié)果較為可靠，該辨識(shí)方法用于辨識(shí)磁懸浮軸承的剛度阻尼參數(shù)具有一定的可行性。

3 磁懸浮軸承剛度和阻尼的實(shí)驗(yàn)研究

與仿真的算法相同，得到兩側(cè)磁懸浮軸承處的位移后，帶入編制好的Matlab程序進(jìn)行計(jì)算。對(duì)磁懸浮軸承的徑向剛度阻尼值進(jìn)行求解。各路方向各轉(zhuǎn)速下磁懸浮軸承的剛度與阻尼曲線如圖5所示。磁懸浮軸承在轉(zhuǎn)子頻率較低時(shí)剛度相對(duì)較小，隨著轉(zhuǎn)速的增加，剛度逐漸上升，阻尼在低速時(shí)也相對(duì)較小，隨著轉(zhuǎn)速的增加，阻尼逐漸上升。

將辨識(shí)所得的參數(shù)記錄成txt文件，用Matlab編制好的程序?qū)⑵渲械谋孀R(shí)出的磁懸浮軸承的剛度阻尼系數(shù)，代入上一節(jié)的仿真文件中，替換掉其中轉(zhuǎn)子的支承參數(shù)，計(jì)算對(duì)應(yīng)不同頻率，磁懸浮軸承處轉(zhuǎn)子的位移響應(yīng)，最后將此位移響應(yīng)與實(shí)驗(yàn)采集的位移響應(yīng)進(jìn)行對(duì)比，以驗(yàn)證辨識(shí)結(jié)果的可靠性，EKF方法辨識(shí)效果驗(yàn)證如圖6所示。

圖5 各路方向各轉(zhuǎn)速下磁懸浮軸承的剛度與阻尼系數(shù)

由驗(yàn)證結(jié)果圖來(lái)看，在200Hz以下，將辨識(shí)出的數(shù)值反帶回仿真文件中的計(jì)算出的幅值，與試驗(yàn)值基本吻合，相差較小，辨識(shí)結(jié)果可信。而200Hz轉(zhuǎn)速以上，兩者差值比較大，并在260Hz左右達(dá)到最大，說(shuō)明在200Hz以上的辨識(shí)結(jié)果誤差逐漸增大，辨識(shí)結(jié)果可靠程度降低。所以，在200Hz以下的范圍，辨識(shí)結(jié)果相對(duì)可靠。

圖6 EKF方法辨識(shí)效果驗(yàn)證

4 結(jié)束語(yǔ)

提出了一種磁懸浮軸承剛度和阻尼系數(shù)的測(cè)試方法。進(jìn)行了不平衡質(zhì)量試驗(yàn)仿真，根據(jù)擴(kuò)展卡爾曼濾波算法建立迭代方程，對(duì)磁懸浮轉(zhuǎn)子系統(tǒng)支撐參數(shù)進(jìn)行了辨識(shí)。對(duì)比設(shè)定值表明，識(shí)別方法對(duì)磁懸浮軸承參數(shù)辨識(shí)的有效性。并進(jìn)行不平衡響應(yīng)實(shí)驗(yàn)，對(duì)磁懸浮轉(zhuǎn)子實(shí)驗(yàn)臺(tái)進(jìn)行辨識(shí)支承參數(shù)，辨識(shí)結(jié)果帶回仿真文件的結(jié)果表明，辨識(shí)結(jié)果比較可靠。

[1] Brad A.Miller,Samuel A.Howard.Identifying bearing rotor-dynamic coefficients using an extended Kalman filter[J].Tribology Transactions,2009,52(5):671-679.

[2] 寶志雯,史文月.基于推廣卡爾曼濾波算法的結(jié)構(gòu)模型的參數(shù)識(shí)別[J].振動(dòng)與沖擊,1990(1):11-21.

[3] 童余德,周永余,陳永冰,周崗. 基于Matlab的卡爾曼濾波法參數(shù)辨識(shí)與仿真[J].船電技術(shù),2009(8):47-50.

[4] 樊素英,李忠獻(xiàn).基于廣義卡爾曼濾波的橋梁結(jié)構(gòu)物理參數(shù)識(shí)別[J].計(jì)算力學(xué)學(xué)報(bào),2007(4):472-476.

[5] 吳子燕,丁蘭,劉書(shū)奎.基于廣義卡爾曼濾波的橋梁結(jié)構(gòu)物理參數(shù)識(shí)別的子結(jié)構(gòu)法[J].西北工業(yè)大學(xué)學(xué)報(bào),2010(3):425-428.

[6] Gelb A.Applied Optimal Estimation.Cambridge：MIT Press, 1974.

[7] 仲衛(wèi)進(jìn),艾芊.擴(kuò)展卡爾曼濾波在動(dòng)態(tài)負(fù)荷參數(shù)辨識(shí)中應(yīng)用[J].電力自動(dòng)化設(shè)備,2007(2):47-50.

Identification of the Magnetic Bearings Support Parameters Using an Extend Kalman Filter

ZOUYue,ZHOUJin

(College of Mechanical and Electrical Engineering,Nanjing University of Aeronautics and Astronautics,Nanjing 210016,China)

An identify measurement of stiffness and damping to AMB system is introduced.By applying unbalance excitation to the rotor,stiffness and damping can be identified use extended Kalman filter(EKF) algorithm.The simulation result shows that it is efficient for this method to measure the stiffness and damping of AMB.Then used the method to identify the supporting parameters in experiment and brought the identified results to simulate files to compare the displacement response.The results indicated that the EKF method can efficiently identify the stiffness and damping.

extend Kalman filter;stiffness;damping；Magnetic bearing

2014-02-24

國(guó)家自然科學(xué)基金資助項(xiàng)目(51075200);江蘇省自然科學(xué)基金資助項(xiàng)目(BK2011070)

TH999

A

1001-2257(2014)11-0012-04

鄒玥(1989-),女，安徽蕪湖人，碩士研究生，研究方向?yàn)榇艖腋≥S承技術(shù)；周瑾(1972-)，女，江蘇徐州人，博士,教授，研究方向?yàn)榇艖腋〖夹g(shù)、轉(zhuǎn)子動(dòng)力學(xué)和機(jī)電系統(tǒng)控制。