楊永剛，丁有闖

（中國民航大學(xué)航空工程學(xué)院，天津300300）

目前，磁流變阻尼器（MRD,magneto-rheological damper）以其性能良好、耗能少、結(jié)構(gòu)簡單、響應(yīng)速度快、阻尼力大且可調(diào)等優(yōu)點(diǎn)，廣泛應(yīng)用于航空、汽車、建筑等領(lǐng)域[1]。但由于磁流變阻尼器動力學(xué)模型的滯回特性，使得其數(shù)學(xué)模型非常復(fù)雜。為精確有效地實(shí)現(xiàn)阻尼器輸出阻尼力的控制，必須保證模型的精確性和有效性。

常見磁流變阻尼器模型有Bingham 模型、修正的Dahl 模型、Bouc-Wen[2]模型、雙曲正切模型[3]等。其中Bouc-wen 模型是最典型的模型之一，但該模型的待識別參數(shù)多達(dá)8 個，因此，模型參數(shù)的精確識別具有重大意義。為此，國內(nèi)外學(xué)者對磁流變阻尼器模型的參數(shù)辨識做了大量研究。Peng 等[4]雖然利用Matlab 工具箱識別出了Bouc-Wen 模型參數(shù)，但并未對識別出的參數(shù)結(jié)果進(jìn)行檢驗(yàn)和證明。Talatahari 等[5]有效引入了一種基于帝國競爭算法（EICA,enriched imperialist competitive algorithm）的優(yōu)化方法來處理Bouc-Wen 模型非線性優(yōu)化問題，但不適用于其他模型。Yang 等[6]對磁流變阻尼器進(jìn)行了性能試驗(yàn)及力學(xué)模型分析。王修勇等[7]采用粒子群算法對磁流變阻尼器模型進(jìn)行了參數(shù)識別，但迭代次數(shù)過大，尋優(yōu)過程繁瑣。胡海剛等[8]在試驗(yàn)數(shù)據(jù)基礎(chǔ)上，采用1stopt 軟件的通用全局優(yōu)化算法（UGOA,universal global optimization algorithm）識別出Bouc-Wen 模型參數(shù)，但識別精度不高。文獻(xiàn)[9-12]首先利用遺傳算法對參數(shù)化模型進(jìn)行了參數(shù)識別，但該算法對計(jì)算機(jī)內(nèi)存和CPU 運(yùn)算速度要求較高；然后基于模型參數(shù)識別法、遺傳算法及模式搜索法對Bouc-Wen 模型進(jìn)行參數(shù)識別，但該過程過于復(fù)雜。

鑒于以上辨識方法具有過程復(fù)雜、識別精度不高等缺點(diǎn)，提出一種基于Matlab/Simulink 辨識工具箱的最小二乘法對Bouc-Wen 模型進(jìn)行參數(shù)識別。首先利用磁流變減擺器力學(xué)性能試驗(yàn)數(shù)據(jù)，在Simulink 中搭建Bouc-Wen 模型，其次利用Parameter Estimation 工具箱，選擇非線性最小二乘法迭代算法對模型參數(shù)進(jìn)行了辨識，并結(jié)合參數(shù)的物理意義對識別結(jié)果進(jìn)行分析，確定隨電流變化參數(shù)，進(jìn)而確定所有參數(shù)辨識結(jié)果。最后采用其他幅值頻率的試驗(yàn)數(shù)據(jù)對Bouc-Wen 模型進(jìn)行驗(yàn)證，結(jié)果表明，試驗(yàn)結(jié)果和仿真結(jié)果擬合度較高，該參數(shù)模型能較好地描述阻尼力的滯回特性。

1 磁流變阻尼器力學(xué)性能試驗(yàn)

磁流變阻尼器力學(xué)性能試驗(yàn)在拉伸試驗(yàn)臺上進(jìn)行，磁流變減擺器[13]結(jié)構(gòu)原理圖如圖1所示。實(shí)驗(yàn)平臺利用激振器產(chǎn)生的不同頻率幅值信號及直流電源為夾持的磁流變減擺器提供的電流，產(chǎn)生阻尼力和活塞桿位移等參數(shù)的數(shù)據(jù)變化。其中，激勵信號為正弦信號x=Msin（2πft），幅值M=10 mm，頻率f=0.4 Hz，電流強(qiáng)度分別取為0、0.25、0.50、0.75、1.00 A。各種工況下功量圖與試驗(yàn)數(shù)據(jù)處理得到的位移-阻尼力和速度-阻尼力曲線圖如圖2～圖3所示。

圖1 結(jié)構(gòu)原理圖Fig.1 Structural schematic diagram

圖2 位移-阻尼力曲線Fig.2 Curve of displacement-damping force

圖3 速度-阻尼力曲線Fig.3 Curve of velocity-damping force

2 Bouc-Wen 模型

Bouc-Wen 模型雖然在數(shù)學(xué)表達(dá)式上復(fù)雜，但其易于處理且能較好地描述磁流變阻尼器在低速區(qū)的非線性滯回特性和能耗性，因此，該模型可用來做磁流變阻尼器阻尼力模型。Bouc-Wen 模型簡圖如圖4所示。

圖4 Bouc-Wen 模型Fig.4 Bouc-Wen model

其微分方程表達(dá)式為

式中：F 為磁流變阻尼器產(chǎn)生的阻尼力；c0為磁流變液阻尼系數(shù)；k0為剛度系數(shù)；x 為位移；x0為初始位移；α為滯回變量的系數(shù)；滯回變量為

式中：n 為曲線圓滑系數(shù)；γ 為滯回模型的寬度調(diào)節(jié)系數(shù)；β 為滯回模型的高度調(diào)節(jié)系數(shù)；A 為比例系數(shù)。

在Simulink 環(huán)境下搭建的Bouc-Wen 模型如圖5所示。

圖5 Simulink 環(huán)境下搭建的Bouc-Wen 模型Fig.5 Bouc-Wen model built in Simulink environment

3 參數(shù)辨識

3.1 非線性最小二乘法

非線性最小二乘法是最小二乘法的分支，可用來估計(jì)非線性模型參數(shù)。其原理是通過求取誤差平方和的目標(biāo)函數(shù)Q 的最小值，從而得到對所估計(jì)目標(biāo)函數(shù)的最優(yōu)解。

設(shè)所辨識模型的數(shù)學(xué)表達(dá)式為

式中：y 為系統(tǒng)的輸出；x′，x″，…為輸入；θ′，θ″，…為參數(shù)。在估計(jì)參數(shù)時模型的數(shù)學(xué)表達(dá)式f 是已知的，經(jīng)過實(shí)驗(yàn)取得數(shù)據(jù)為（x1′，x1″，…，y1），（x2′，x2″，…，y2），…，（xn′，xn″，…，yn）。

非線性模型誤差平方和的目標(biāo)函數(shù)Q 為

3.2 Bouc-Wen 模型的參數(shù)辨識

Bouc-Wen 模型的參數(shù)達(dá)8 個，使用遺傳算法、粒子群優(yōu)化算法等識別方法代碼繁瑣，比較復(fù)雜。使用基于Matlab 最小二乘法的迭代算法，利用試驗(yàn)數(shù)據(jù)自動調(diào)用迭代算法，使試驗(yàn)值和仿真值無限接近，達(dá)到辨別所有參數(shù)的目的。辨識流程圖如圖6所示，具體步驟如下：

圖6 辨識流程圖Fig.6 Identification flow chart

1）利用已搭建的Bouc-Wen 模型，通過Parameter Estimation 工具箱導(dǎo)入試驗(yàn)數(shù)據(jù)，輸入數(shù)據(jù)為位移和速度，輸出數(shù)據(jù)為阻尼力；

2）根據(jù)經(jīng)驗(yàn)，設(shè)定初次辨識的參數(shù)初始值分別為：α=30 N·mm-1、β=-500、γ=800、n=2、x0=-26 mm、k0=300 N·mm-1、A=1 000、c0=80 N·s·mm-1。確定辨識范圍，最后數(shù)據(jù)上傳模型；

3）點(diǎn)擊option 選擇優(yōu)化選項(xiàng)，選擇最小二乘法迭代算法，代價函數(shù)選擇和平方差，最后點(diǎn)擊Estimation進(jìn)行參數(shù)估計(jì)。計(jì)算機(jī)會根據(jù)試驗(yàn)曲線和仿真曲線的接近度根據(jù)判定規(guī)則停止參數(shù)估計(jì)。辨識結(jié)束后的試驗(yàn)值與仿真值擬合圖如圖7所示。參數(shù)辨識結(jié)果數(shù)據(jù)如表1所示。

圖7 試驗(yàn)值和仿真值的曲線擬合Fig.7 Curve fitting between test and simulation value

由表1可知，不論是試驗(yàn)值還是理論值，n 和x0都近似為常數(shù)值，因此取其平均值n=1.01、x0=-25.64 mm作為下次參數(shù)識別的初始值。由于參數(shù)A、β、γ 是滯回模型的調(diào)節(jié)系數(shù)，在對誤差平方和最小的追求中，A、β、γ 相較于α、c0、k0更容易實(shí)現(xiàn)調(diào)節(jié)。同時在表1中可發(fā)現(xiàn)，由于采用第一次參數(shù)估計(jì)后得出的參數(shù)值作為下次參數(shù)估計(jì)的初始值，存在著迭代關(guān)系，所以表1中后3 組的A、β、γ、k0值受電流影響較小，更具參考意義，因此將結(jié)合后3 組的數(shù)據(jù)采用逐步縮小邊界法再次進(jìn)行參數(shù)識別。因此接下來要識別6 個參數(shù)，因此減小了識別難度。以首次辨識電流為0.5 A 時識別出的參數(shù)為初始值，辨識結(jié)果如表2所示。

表1 各參數(shù)的辨識結(jié)果Tab.1 Identification result of each parameter

表2 縮小辨識邊界范圍的辨識結(jié)果Tab.2 Identification result with reduced boundary range

由表2可以看出參數(shù)A、β、γ、k0變化幅度不大，可以其平均值作為其模型的識別后參數(shù)，A = 3 507.8、β=-1 013.6、γ=1 235.8、k0=68.31 N·mm-1。α 和c0隨電流變化而有規(guī)則地變化，因此即可認(rèn)為α 和c0與電流有以下關(guān)系

借助Matlab 曲線擬合工具箱（curve fitting tool），擬合結(jié)果如圖8所示，可得出α1、α2、c1、c2、c3。

圖8 α、c0 曲線擬合Fig.8 Curve fitting between α and c0

因此，由擬合結(jié)果可知α1= 20.204 N·mm-1、α2=16.248 N·mm-1、c1=-30.651 N·s·mm-1、c2=50.675 N·s·mm-1、c3= 34.049 N·s·mm-1。至此最終辨識結(jié)果為：α1= 20.204 N·mm-1、α2= 16.248 N·mm-1、c1= -30.651 N·s·mm-1、c2=50.675 N·s·mm-1、c3=34.049 N·s·mm-1、A=3 507.8、β=-1 013.6、γ=1 235.8、k0=68.31 N·mm-1、n=1.01、x0=-25.64 mm。

4 模型驗(yàn)證

為了驗(yàn)證辨識參數(shù)結(jié)果的正確性，在辨識模型基礎(chǔ)上增加一個可發(fā)出隨時間變化正弦信號的模塊組代替輸入模塊x，輸入數(shù)據(jù)由變量改為辨識后得出的參數(shù)常量，c2、α 作為電流I 的函數(shù)輸入，將輸出模塊數(shù)據(jù)導(dǎo)入工作空間，得到了仿真模型如圖9所示。

圖9 參數(shù)驗(yàn)證仿真模型Fig.9 Parameter verification simulation model

將上述辨識結(jié)果代入仿真模型，首先驗(yàn)證同幅值頻率的任意電流下的仿真輸出數(shù)據(jù)與試驗(yàn)數(shù)據(jù)的吻合程度。位移-阻尼力的仿真試驗(yàn)對比圖和速度-阻尼力的仿真試驗(yàn)對比圖如圖10所示。

圖10 0.4 Hz、10 mm、0.5 A 條件下的仿真試驗(yàn)Fig.10 Simulation experiment diagram under 0.4 Hz,10 mm,0.5 A

由圖10可知，位移-阻尼力仿真試驗(yàn)圖以及速度-阻尼力仿真試驗(yàn)圖吻合度較好，驗(yàn)證了參數(shù)辨識結(jié)果的正確性。為了驗(yàn)證基于Matlab/Simulink 辨識工具箱——最小二乘法方法的正確性和通用性，再隨機(jī)選取不同幅值頻率的數(shù)據(jù)來進(jìn)行驗(yàn)證，如圖11所示。

圖11 0.7 Hz、7.5 mm、0.5 A 條件下的仿真試驗(yàn)Fig.11 Simulation experiment diagram under 0.7 Hz、7.5 mm、0.5 A

由圖11可知，選取0.5 A，頻率為0.7 Hz，幅值為7.5 mm 的數(shù)據(jù)進(jìn)行驗(yàn)證，也能較好地反映辨識結(jié)果的正確性以及該動力學(xué)模型的通用性。

5 結(jié)語

首先在Simulink 中搭建Bouc-Wen 模型，然后利用在拉伸試驗(yàn)臺上做的阻尼器力學(xué)性能試驗(yàn)數(shù)據(jù)，基于Parameter Estimation 工具箱使用非線性最小二乘法迭代算法對Bouc-Wen 模型進(jìn)行參數(shù)識別，最后搭建仿真驗(yàn)證模型對參數(shù)辨識結(jié)果進(jìn)行驗(yàn)證。得出結(jié)論如下：

1）辨識出的Bouc-Wen 模型參數(shù)不僅能很好地與參與辨識的數(shù)據(jù)相吻合，而且還能較好地表達(dá)阻尼器其他幅值頻率下的動力學(xué)特性；

2）辨識出的模型各參數(shù)物理意義清晰，有助于進(jìn)行下一步的減擺器控制工作。使用非線性最小二乘法，算法效率高，易于實(shí)現(xiàn)，對于其他阻尼器模型參數(shù)辨識也具有通用性。