基于粒子群優化算法的磁流變阻尼器多項式動力學建模方法

董致臻 馮志敏 伍廣彬 劉小鋒

1.寧波大學海運學院，寧波，315211 2.哈爾濱工業大學機器人技術與系統國家重點實驗室，哈爾濱，150001

0 引言

磁流變液作為一種智能材料在無外加磁場時黏度很低，保持Newton流體特性；而在磁場作用下，磁性顆粒在載液中迅速形成鏈狀結構，使磁流變液具有一定的屈服強度，呈現Bingham流體特性。磁流變阻尼器（magnetorheological damper，MRD）利用這一特性在工程領域得到廣泛應用［1?3］。

建立精確的MRD動力學模型是實現復雜結構控制的關鍵。針對MRD動力學建模問題，國內外學者已做了大量研究，主要根據參數是否具有物理意義分為參數模型和非參數模型［4］。在參數模型中SPENCER等［5］提出Bonc?Wen模型，能較好地反映MRD滯回特性；KWOK等［6］基于粒子群優化（particle swarm optimization，PSO）算法提出的MRD雙曲正切模型比Bonc?Wen模型所含參數少，易于推廣應用。在非參數模型中，CHANG等［7］提出神經網絡模型，建模方便且模型精度高；GAVIN等［8］最早在電流變阻尼器中應用Cheby?shev多項式模型，模型參數無需算法辨識且容易建模，模型精度取決于插值節點對應的數據誤差；CHOI等［9］運用最小二乘法辨識傳統多項式模型，并將其用于MRD動力學建模中，該方法存在高階Runge振蕩問題。近年來，為提高MRD多項式動力學模型精度和消除Runge振蕩，周鐵明等［10］提出將模型分成多個低階多項式分別擬合再拼接的分段建模思想；孔祥東等［11］采用將低階多項式模型與Bingham力學模型相結合的建模方法，提高了模型的整體擬合精度。但上述模型仍存在非線性參數多、計算冗長或模型分段處精度不足和Runge振蕩現象，因此本文提出一種基于PSO算法［12］并應用Lagrange插值理論［13］的MRD動力學建模方法。

1 MRD性能試驗

1.1 MRD試驗平臺

為研究MRD動力學特性，試驗選用SG?MRD60型兩級線圈套筒式磁流變阻尼器［14］。阻尼力檢測在中機試驗裝備股份有限公司生產的SDS?100型動靜試驗機上進行，由液壓動力源為伺服系統提供動力，以德國DOLI公司的EDC580控制器作為控制系統，并用Agilent E3632A電流源控制輸出電流。試驗平臺和工作原理見圖1。

圖1 試驗平臺和工作原理Fig.1 Test platform and working principle

1.2 試驗方法

為探究不同外部激勵與控制電流下的MRD動力學特性，通過EDC控制器在動靜試驗機輸出不同頻率、振幅的正弦激勵信號，對應電流源輸出給阻尼器的不同電流，分別進行阻尼器動力學性能試驗。由于橋梁拉索振動頻率主要在0.25～2.0 Hz之間，振幅由風雨激勵和MRD安裝位置決定［15］，故選取正弦激勵（頻率 1Hz、振幅15 mm），輸入電流為0、0.25 A、…、1.5 A的試驗工況。試驗機采集載荷（阻尼力）信號，作為模型辨識依據。

2 模型參數識別

2.1 傳統多項式模型

建立MRD傳統多項式動力學模型［9］需根據活塞運動方向，把阻尼器動力學特性曲線分為一個正向運動曲線與一個反向運動曲線疊加。其中，輸出阻尼力

式中，ai為多項式系數；v為活塞瞬時速度；n為多項式階數；I為電流；bi、ci為電流相關系數，可根據F、v、I試驗值辨識得到。

根據活塞運動方向，將v≥0定義為正向運動，v＜0定義為反向運動。

2.2 Chebyshev多項式模型

建立Chebyshev多項式模型［8］需要先確定多項式階數，再計算Chebyshev插值節點，最后將插值點試驗數據代入Lagrange插值公式。建立n階插值多項式模型需要n+1個插值節點，區間［a，b］上的Chebyshev插值節點［13］:

其中，n為多項式階數，區間[a,b]為對應工況下MRD的活塞速度變化范圍，可由下式計算得到：

式中，A為激振振幅；f為激振頻率；φ0為初始活塞相位；t為當前時刻。

n階Lagrange插值多項式由n+1個n階插值函數構成，插值函數

構成的n階Lagrange插值多項式為

式中，yi為插值點對應的函數值。

2.3 模型辨識對比及分析

由于篇幅所限，僅選用6階、12階多項式進行分析。根據試驗工況數據，分別對上述2種模型進行參數辨識并作對比分析。其中,多項式模型用MATLAB最小二乘法工具箱作參數辨識，得到6階多項式模型參數辨識結果（表1），繪制的速度?阻尼力曲線見圖2、圖3。

表1 6階多項式模型參數識別結果Tab.1 The results of parameter identifications of the 6th order polynomial model

由于Chebyshev插值節點需根據活塞速度變化范圍確定，將工況點試驗參數代入式（4），得到活塞速度變化范圍為±94.25 mm/s。再將式（3）求得的Chebyshev插值點對應的試驗數據代入Lagrange插值多項式，并經MATLAB Simplify函數對多項式化簡。最終6階Chebyshev多項式模型插值節點與對應多項式系數見表2、表3。

表2 6階Chebyshev多項式模型參數結果Tab.2 The parameter results of the 6th order Chebyshev polynomial model

圖2 6階多項式模型及試驗曲線Fig.2 The 6th onder polynomial models and test curves

由圖2、圖3可知，高階傳統多項式模型與試驗測量的輸出阻尼力差值變化很大，并且曲線在兩個端點處發生劇烈Runge振蕩，試驗中發現Runge現象在模型提高到8階時出現，而高階Chebyshev多項式模型避免了端點處Runge振蕩。12階Che?byshev多項式模型擬合效果明顯優于6階Cheby?shev多項式模型，說明Chebyshev多項式模型精度隨多項式階數增加而提高，建模時需兼顧運算速度選擇較高的多項式階數。此外，隨MRD輸入電流的減小，速度?阻尼力曲線平順度變差，尤其在無電流工況，Chebyshev多項式模型曲線過渡段圖像嚴重失真，出現局部凸變現象。

圖3 12階多項式模型及試驗曲線Fig.3 The 12th order polynomial models and test curves

3 動力學優化建模

為優化動力學建模方法，提出利用PSO算法優化插值節點，并結合插值節點對應試驗數據，將其代入Lagrange插值方程，經MATLAB Simplify函數化簡得到所需的多項式模型。

表3 6階Chebyshev多項式模型插值節點Tab.3 The interpolation nodes of the 6th order Chebyshev polynomial model

3.1 PSO算法的插值節點求解

3.1.1 PSO算法原理

PSO算法［12］中每個粒子可以看作是解空間中的一個點，假定粒子的群體規模為m，粒子的維數為d，全部粒子的速度及位置：

式中，i為粒子索引，i=1，2，…，m；k為當前迭代次數；v(k)id為第i個粒子在第k次迭代飛行速度矢量的第d維分量的速度；x(k)id為是第i個粒子在第k次迭代飛行速度矢量的第d維分量的位置；p(k)id為第i個粒子的個體最優解的第d維分量；p(k)gd為群體最優解的第d維分量；c1、c2為學習因子；r1、r2為（0，1）的隨機數。

上述變量中，所有粒子的位置范圍都包含在所選試驗工況之內，每個粒子x的d維位置變量是多項式模型中d個插值節點在MRD動力學試驗數據中的索引值，構建n階多項式所需的粒子維度d=n+1。PSO算法優化結束后，最優插值節點的索引值保存在全局歷史最優位置pg中。

適應度函數為PSO算法的優化目標函數，在插值節點求解時，適應度函數反映了多項式計算模型輸出值與實際試驗值之間的偏差。優化通常以最小化適應度函數值為目標，表達式為

式中，f（x）為適應度函數；Fj為第j個阻尼力試驗值；Ln(xi)為xi粒子插值確定的多項式模型；vi為第j個活塞速度試驗值，dat為建模數據量。

3.1.2 算法優化流程

對PSO算法進行改進，首先加入排序算法，將所有粒子（插值向量）內各插值節點的索引值按由小到大的順序排列，消除不同索引順序的最優解更新后對其余粒子優化方向的干擾，通過保持所有粒子位置的有序性來縮短算法收斂時間；然后將全局適應度計算放在個體適應度計算之后進行，以減少運算時間。

采用PSO算法進行插值節點優化，優化過程見圖4。主要步驟如下：

圖4 插值節點優化流程圖Fig.4 Interpolation node optimization flow chart

（1）初始化種群，即隨機生成粒子位置、速度信息,然后將各粒子位置信息代入式（9）計算，得到初始個體最優值。

（2）將式（8）更新后的各粒子位置信息代入式（6），得到Lagrange插值多項式，并計算式（9）適應度值。

（3）若粒子當前適應度值優于個體最優值，則將該粒子排序并更新個體最優值，否則進入下一輪循環；若個體最優的粒子同時達到全局最優，則更新全局最優值，否則直接進入下一輪循環。

（4）重復步驟（2）～（3），達到最大迭代次數后程序結束。

將優化后插值節點對應試驗數據代入La?grange插值多項式，得到Lagrange插值多項式模型。經MATLAB Simplify函數化簡，得到優化后的多項式動力學模型。

3.2 求解結果與比較分析

為體現PSO算法在高階多項式建模方法中的作用，選取12階PSO算法優化多項式模型與Che?byshev多項式模型，分別在正弦激勵（頻率1 Hz、振幅15 mm），輸入電流0、0.25 A、…、1.5 A工況下進行建模研究。PSO算法參數設置見表4。

表4 PSO算法參數Tab.4 PSO algorithm parameters

根據傳統多項式模型，按活塞運動方向將試驗數據分為正、反2組，并分別結合Lagrange插值多項式，用PSO算法優化插值節點。優化結果代入Lagrange插值公式，形成Lagrange插值多項式，經MATLAB化簡，所得多項式動力學模型速度-阻尼力曲線見圖5，適應度收斂軌跡見圖6。利用MATLAB最小二乘法工具箱結合輸入電流工況，根據式（2）對多項式模型系數進行辨識，得到模型參數（表5）。

圖5 2種多項式模型及試驗曲線Fig.5 Two polynomial models and test curves

圖6 適應度收斂軌跡Fig.6 Fitness convergence trajectory

表5 12階PSO優化多項式模型參數Tab.5 Parameters of the 12th order PSO optimal polynomial model

為進一步評價該建模方法的有效性，根據模型平均累計相對誤差公式：

式中，FP(vi)為動力學模型在速度點vi的阻尼力值；F(vi)為阻尼力在vi速度點的實測值。

對上述試驗結果進行比較。計算結果見表6。

表6 模型相對誤差對比Tab.6 Comparison of model relative error %

由圖5可知，經PSO算法優化插值節點建立的12階多項式模型不僅消除了Runge振蕩，而且擬合精度優于Chebyshev多項式模型。由于Che?byshev插值建模方法中插值節點無法根據實測阻尼力調整，因此，在輸入電流0～1.5 A工況下，PSO算法建模方法的相對累計誤差均小于Che?byshev多項式模型（表6），平均相對誤差下降47.0%。圖5反映出PSO算法具備良好的收斂性，能有效滿足實際工程應用需求。

4 結論

在MRD動力學特性試驗基礎上，對傳統多項式模型與Chebyshev模型進行了模型辨識及對比分析，提出基于PSO算法的多項式動力學優化建模方法，得出如下結論：

（1）傳統多項式建模方法簡便，但會出現6階精度不足、12階Runge振蕩失真問題。Chebyshev多項式模型雖然消除了高階Runge振蕩，并且模型精度隨多項式階數增加而提高，但在低電流工況下，圖像曲線過渡段模型仍然呈現凸變失真現象。

（2）由于Chebyshev多項式模型建模精度受插值點試驗數據精度的影響，并且插值節點的選取只與MRD活塞速度變化范圍有關，無法根據實測阻尼力調整，因此PSO算法建模方法能夠有效提高MRD多項式動力學模型精度。試驗表明，在正弦激勵頻率1 Hz、振幅15 mm、電流0～1.5 A工況下，相比Chebyshev多項式模型，PSO算法多項式模型平均相對累計誤差下降47.0%，能較真實地反映MRD的動力學特性。PSO算法收斂性好，具有良好的工程應用前景。

［1］ VULCU C，DUBIN? D，POPA N，et al.Hybrid SeismicProtection System：Buckling Restrained Brace of Nano?Micro Composite Magneto Rheologi?cal Damper［J］.Ce/papers，2017，1（2/3）：2936?2945.

［2］ 胡紅生，肖平，江明，等.基于魚群算法的永磁體?電磁閥式磁流變阻尼器半主動懸架系統［J］.中國機械工程，2017，28（21）：2526?2533.HU Hongsheng，XIAO Ping，JIANG Ming，et al.Experimentally Reducing Vibration and Noise of Gear System Using Magneto?rheological Damper［J］.China Mechanical Engineering，2017，28（21）：2526?2533.

［3］ WANG Zhihao，CHEN Zhengqing，GAO Hui，et al.Development of a Self?powered Magnetorheological Damper System for Cable Vibration Control［J］.Ap?plied Sciences，2018，8（1）：118.

［4］ JIANG Z，CHRISTENSON R E.A Fully Dynamic Magneto ?rheological Fluid Damper Model［J］.Smart Materials&Structures，2012，21（6）：65002?65013.

［5］ SPENCER B F Jr，DYKE S J，SAIN M K，et al.Phenomenological Model for a Magnetorheological Damper ［J］.Journal of Engineering Mechanics，1997，123（3）：230?238.

［6］ KWOK N M，HA Q P，NGUYEN T H，et al.A Nov?el Hysteretic Model for Magnetorheological Fluid Dampers and Parameter Identification Using Particle Swarm Optimization［J］.Sensors and Actuators A，2006，132：441?451.

［7］ CHANG C，ROSCHK P E.Neural Network Model?ing of a Magnetorheological Damper［J］.Journal of Intelligent Material Systems&Structures，1998，9（9）：351?358.

［8］ GAVIN H P，HANSON R D，FILISKO F E.Elec?trorhelogical Dampers，Part 2：Testing and Modeling［J］.Applied Mechanics，1996，63（3）：676?682.

［9］ CHOI S B，LEE S K，PARK Y P.A Hysteresis Model for the Field?dependent Damping Force of a Magnetorheological Damper［J］.Journal of Sound and Vibration，2001，245（2）：375?383.

［10］ 周鐵明，陳恩偉，陸益民，等.磁流變液阻尼器的改進多項式模型及驗證［J］.振動與沖擊，2014，33（7）：221?226.ZHOU Tieming，CHEN Enwei，LU Yimin，et al.Modified Polynomial Model and Its Verification for a MR Damper［J］.Journal of Vibration and Shock，2014，33（7）：221?226.

［11］ 孔祥東，李斌，權凌霄，等.磁流變液阻尼器Bingham?多項式力學模型研究［J］.機械工程學報，2017，53（14）：179?186.KONG Xiangdong，LI Bin，QUAN Lingxiao，et al.Study on Dynamic Bingham?polynomial Model of a MRF Damper［J］.Journal of Mechanical Engineer?ing，2017，53（14）：179?186.

［12］ KENNEDY J，EBERHART R.Particle Swarm Op?timization［C］//Proceedings of IEEE International Conference on Neural Networks.Perth，1995：1942?1948.

［13］ 李慶揚，王能超，易大義.數值分析［M］.5版.北京：清華大學出版社，2008：23?25.LI Qingyang，WANG Nengchao，YI Dayi.Numeri?cal Analysis［M］.5th ed.Beijing：Tsinghua Univer?sity Press，2008：23?25.

［14］ 李宏偉，冷志鵬，徐慧.一種新型磁流變阻尼器：中國，CN201110378234.6［P］.［2012?04?25］.LI Hongwei，LENG Zhipeng，XU Hui.A Novel Mag?netorheological Damper：China，CN201110378234.6［P］.［2012?04?25］.

［15］ 馮志敏，張興軍，張剛，等.斜拉索?阻尼器系統建模與減振控制研究［J］.農業機械學報，2013，44（增1）：282?287.FENG Zhimin，ZHANG Xingjun，ZHANG Gang，et al.System Modeling and Vibration Reduction Control for Stayed Cable?Magnetorheological Dam?per［J］.Transactions of the Chinese Society for Ag?ricultural Machinery，2013，44（s1）：282?287.*