劉永強,楊紹普,廖英英,張耕寧
(石家莊鐵道大學,石家莊 050043)
磁流變阻尼器是一種智能化的半主動控制實現元件,通過控制線圈中電流指令的大小可以改變磁流變液所處環境的磁場強度,磁場強度的大小直接影響著磁流變液的粘度系數,即影響著液體流經閥門的速度,從而實現阻尼器輸出阻尼力的控制。磁流變阻尼器以反應速度快,消耗能量少和控制范圍廣等優點著稱[1,2],但由于其特殊的力學特性——滯回特性[3],使得它的數學模型十分復雜。其中,Bouc-Wen模型就是較為典型的一個,該模型包含有8個未知參數,再考慮電流大小、活塞速度和激勵等因素,參數識別起來比較困難。
國內外相關學者對此進行了大量的研究。Ikhouane等人[4,5]曾采用解析的方法對 Bouc-Wen 模型的參數識別過程做過詳細的闡述,但由于解析法需要進行大量的假設和定義,過程較為復雜,因此通用性不佳。Spencer和 Dyke 等人[6,7]將 Bouc-Wen 模型中的參數均作為常值,利用有約束的非線性優化算法進行識別,雖然簡單易行,但識別精度不高,且由于不含電流項所得到的模型無法應用于半主動控制。歐進萍、Shen、Liu 以及 Jansen 等人[8~11]曾將部分參數視為電流的函數,在此假設的基礎上利用普通優化算法對參數進行識別,但值得指出的是普通優化方法在多變量識別方面局限性較大,識別精度較低,且通用性差。
本文擬通過對Bouc-Wen模型的結構和力學特性進行分析,采用擅長解決多變量優化問題和全局優化問題的遺傳算法對模型的參數進行辨識,并確定了α,c0,k0三個參數與電流指令間的函數關系和其余5個參數的值。
磁流變阻尼器的結構如圖1所示,利用MTS力學性能測試試驗臺測量其輸出阻尼力,激勵采用正弦信號x=A sin(2πft)。共進行了5組試驗,分別為:A=5 mm、f=0.5 Hz,A=10 mm、f=0.5 Hz,A=20 mm、f=0.5 Hz,A=10 mm、f=1 Hz,A=10 mm、f=1.5 Hz。每組激勵下分別取0 A~3 A范圍內7個電流指令進行測試,共得到35組數據。其中,A=10 mm、f=0.5 Hz激勵下的位移-阻尼力和速度-阻尼力曲線如圖2和圖3所示。由圖中發現,當電流為3 A時,磁流變阻尼器輸出的阻尼力達到了飽和。

圖1 磁流變阻尼器結構Fig.1 The structure form of the MR damper

圖2 位移-阻尼力試驗曲線Fig.2 Displacement-damping force test curve

圖3 速度-阻尼力試驗曲線Fig.3 Velocity-damping force test curve
Bouc-Wen模型最早由 Wen于1976年提出[12],它由滯回系統和彈簧、阻尼器并聯而成,如圖4所示。Bouc-Wen模型能夠很好地模擬滯回特性,且通用性強,易于數值處理。其力學模型描述為:

其中,滯回變量 z由下式決定


圖4 Bouc-Wen模型結構Fig.4 Bouc-Wen model structure
由于存在微分項,采用傳統方法對 Bouc-Wen模型進行數值仿真存在困難。本文在SIMULINK環境下對Bouc-Wen模型進行設計,并采用4階Runge-Kutta法對其進行數值仿真分析,如圖5所示。
Bouc-Wen模型的創始人之一 Wen[12]認為 Bouc-Wen模型中參數γ,β只影響滯回環的形狀,而不影響阻尼力的大小,參數n只影響從彈性區過渡到塑性區曲線的平順性。

圖5 Bouc-Wen模型在Simulink中的實現Fig.5 The realization of Bouc-Wen model in simulink
從試驗曲線圖1和圖2中可以發現,電流I的大小對阻尼力幅值的影響較大。那么,到底Bouc-Wen模型中哪些參數與電流有關呢?Shen和 Jansen等[9,11]認為參數α,c0與電流之間存在著如下函數關系

而 Liu 等人[10]則將參數 α,c0,k0視為電流的多項式函數。因此,本文假設參數α,c0,k0與電流指令間存在著函數關系,但具體函數表達式待定。
Bouc-Wen 模型中待辨識的參數有 α,c0,k0,γ,β,A,n,x0共8個之多,辨識起來比較復雜。本文擬采用擅長多變量優化和全局搜索技術的遺傳算法工具箱進行參數的識別工作。
1962 年 Holland[13,14]首次提出一種嶄新的全局優化算法——遺傳算法,它借用了生物遺傳學的觀點,通過自然選擇、遺傳、變異等作用機制,體現了自然界中“物競天擇、適者生存”進化過程,且計算速度快,效果比較理想,從而被迅速推廣到優化、搜索、機器學習等方面,亦在復雜數學模型的參數識別方面得到廣泛應用。MATLAB中遺傳算法的命令為:


其中,fitnessfiunc為適應度函數;N var s為待求未知變量的個數;options為GA的結構參數。
GA的線性約束表達為:

GA的非線性約束用Nonlcon函數來表達,它包括:

定義Bouc-Wen模型的適應度函數為:

其中,m為試驗數據點個數;Fsim和Fexp分別為通過仿真和試驗得到的阻尼力大小;和分別為試驗數據中最大和最小阻尼力值。
遺傳算法是一種隨機搜索、逐漸尋優的算法,其結果具有隨機性。對任意一個MR阻尼器來說,由于Bouc-Wen模型不存在真實解,也就無法判斷所識別的參數是否全局最優解,因此判斷識別結果是否準確的唯一標準就是仿真結果與試驗數據之間的吻合程度,即適應度函數值。另外,經過大量的仿真試驗發現在保持其余參數固定不變的條件下,各參數的單調變化(單調遞增或遞減)都會引起滯回環形狀或阻尼力幅值的單調變化。Bouc-Wen模型參數的這種特殊的單調特性使得只要目標函數足夠小,得到的結果就應該在全局最優解的附近。
基于以上認識,本文擬采用逐漸縮小參數取值范圍的方法來提高識別精度。利用幅值為A=10 mm,頻率為f=0.5 Hz的正弦激勵下電流為0 A~3 A的試驗數據參與模型參數辨識,對辨識結果進行數值仿真,對比其他激勵狀況下的試驗數據與仿真結果間的吻合情況,以驗證辨識結果的正確性。
Bouc-Wen 模型共包括 c0,k0,x0,α,γ,β,A,n 等 8個未知參數需要識別。從圖2和圖3可以看到,磁流變阻尼器的滯回曲線不存在偏移現象,因此可知x0=0。對于固定的磁流變液,參數n為固定值,本文所用磁流變液為 n=2。因此,Bouc-Wen 還有 c0,k0,α,γ,β,A等6個參數待識別。
本文采用初始群體個數為N=50,采用實數編碼方式,選擇算法采用隨機均勻分布選擇;交叉采用分散交叉的方法,交叉概率取0.8(變異概率為0.2);變異采用高斯函數方法,生成服從高斯分布、均值為0的隨機數加到父代上,尺度取0.5,壓縮取0.7。
由文獻[6,12] 可知Bouc-Wen模型中各參數均取正值,但由于無具體參數取值范圍的界定,進行識別時應取較寬的取值范圍,如:分別在電流指令為0 A~3 A條件下,對Bouc-Wen模型各參數進行GA辨識運算,迭代次數為500。

由于在參數取值范圍較寬的條件下得到的參數適應度函數值較大,I=0 A時適應度值只能達到f1=5.41e-2。并且由于遺傳算法在局部尋優能力較弱,因此需要縮小參數取值范圍來提高識別精度。
從初步識別得到的參數值中選擇各自的最小值作為低限值LB,最大值作為高限值UB,再次進行GA運算。不同的電流指令作用下,各參數的取值范圍會稍有不同。以 I=0 A 時為例,c0,k0,α,γ,β,A 的線性約束取為:

此時得到的適應度值減小為 f2=1.38e-0.02。以此類推,逐漸縮小參數的取值范圍,參數識別的精度也逐漸提高。經過10組運算后,適應度值減小到f10=1.36e-0.03。
運算得到的各參數值及誤差如表1所示。繪制各參數隨電流的變化曲線如圖6所示。

表1 GA運算得到的參數值及誤差值Tab.1 The parameter and error values of GA

圖6 各參數隨電流的變化趨勢Fig.6 The change trends of parameters under different
從圖6可以發現γ,β,A為隨機變化,無明顯規律可循,α,c0,k0隨電流的變化趨勢則可以近似表達為3階多項式形式:

因此,需要辨識的參數由原來的6個增加為12個:

根據式(9)修改Bouc-Wen仿真模型,并根據式(6)中的適應度函數,定義包含所有電流的總適應度函數為:

式中,p為試驗中電流的個數,本文取p=7;fitnrssfunci為單個電流作用下的適應度值。
采用遺傳算法對所有12個參數進行參數識別,采用逐步縮小參數取值范圍的方法提高識別精度,式(9)中衍生出的9個參數可在[-1e8,1e8] 范圍內取值,其余參數的取值范圍為[0,1e8] ,最終得到總適應度值為fitnrssfuncall=1.107e-3時的參數識別結果為:

因此,Bouc-Wen模型最終可以表示為:

其中:

模型驗證分為擬合結果驗證和預測驗證兩部分。擬合結果驗證是為了驗證對參數 α,c0,k0與電流間函數關系進行的3階多項式假設的正確性;預測驗證則是為了驗證采用遺傳算法辨識得到的Bouc-Wen模型能否真正表達該磁流變阻尼器的力學特性。
首先驗證利用逐漸縮小參數取值范圍的方法得到的Bouc-Wen模型的識別結果是否為全局最優解。
遺傳算法的全局搜索能力較強,但局部搜索能力較弱,得到的結果無法保證為全局最優解。而模式搜索(Pattern Search,PS)則在求解全局最優解方面表現優異,但它對初始值的設定依賴性較強。因此本文利用遺傳算法的結果作為模式搜索的初始值,較寬的取值范圍內搜索全局最優解,以驗證遺傳算法結果的正確性。
按照式(7)設定為參數的取值范圍,將式(11)作為初始值,采用模式搜索法得到的近似的全局最優解為:

對比式(11)和式(13)中的結果,可以發現二者的結果基本相同,由此可以證明利用逐步縮小取值范圍的方法提高Bouc-Wen模型參數識別精度的方法是可行的,由此得到的解在全局最優解的附近。
對式(12)所示的Bouc-Wen模型進行數值仿真分析,采用 Runge-Kutta法進行定步長求解,步長取為1 e-2 s,仿真時間為2 s,以獲得與試驗數據點個數相同的仿真值。當正弦激勵幅值為A=10 mm,頻率為f=0.5 Hz時,分別繪制仿真與試驗數據的位移-阻尼力和速度-阻尼力曲線圖,如圖7所示。

圖7 A=10 mm,f=0.5 Hz時仿真值與試驗數據對比Fig.7 Comparation between simulation results and experiment data when A=10 mm,f=0.5 Hz
從圖中可以看出,仿真結果與試驗數據基本吻合,由此可見本文對函數關系的假設是合理的,所得到的曲線擬合和辨識結果均滿足要求。
由于Bouc-Wen模型的參數是通過振幅為A=10 mm,頻率為f=0.5 Hz的正弦激勵下的試驗數據獲得的,因此需要用其他幅值和頻率下的試驗數據進行驗證才能證明該模型的正確性。

圖 8 A=5 mm,f=0.5 Hz,I=0.5 A時仿真值與試驗數據對比Fig.8 Comparation between simulation results and experiment data when A=5 mm,f=0.5 Hz,I=0.5 A
預測結果驗證,即將數值仿真結果分別與 A=5 mm,f=0.5 Hz,A=20 mm,f=0.5 Hz,A=10 mm,f=1 Hz,A=10 mm,f=1.5 Hz正弦激勵中不同電流指令作用下的試驗數據進行對比,如圖8~圖11所示。

圖 9 A=10 mm,f=1 Hz,I=1 A時仿真值與試驗數據對比Fig.9 Comparation between simulation results and experiment data when A=10 mm,f=1 Hz,I=1 A

圖 10 A=20 mm,f=0.5 Hz,I=1.5 A 時仿真值與試驗數據對比Fig.10 Comparation between simulation results and experiment data when A=20 mm,f=0.5 Hz,I=1.5 A

圖 11 A=10 mm,f=1.5 Hz,I=2 A 時仿真值與試驗數據對比Fig.11 Comparation between simulation results and experiment data when A=10 mm,f=1.5 Hz,I=2 A
從圖8~圖11中可以看出,本文通過遺傳算法識別出的Bouc-Wen模型不僅與參與辨識的試驗數據吻合較好,而且對其他幅值和頻率的正弦激勵下的MR阻尼器的動力學響應也能準確反映。
本文在基于Bouc-Wen模型部分參數與電流指令間存在函數關系的假設上,通過遺傳算法對模型的參數進行識別,采用逐步縮小參數取值范圍的方法提高識別精度,并采用模式識別方法對該方法識別出的結果進行驗證。最后,采用數值仿真的方法對函數關系的假設、曲線擬合及參數識別結果進行驗證。結果顯示,由于Bouc-Wen模型參數特殊的單調特性,通過逐步縮小參數取值范圍來提高識別精度的方法可以找到近似的全局最優解。識別出的Bouc-Wen模型不僅與參與識別的試驗數據相吻合,而且還可以準確表達其他幅值和頻率正弦激勵下的MR阻尼器動力學響應。
[1] 陳進新,王 宇.一種可應用于小阻尼力領域的新型磁流變液阻尼器[J] .振動與沖擊,2009,28(2):155 -161,166.
[2] 施 亮,何 琳.磁流變阻尼器參數辨識方法研究[J] .振動與沖擊,2009,28(1):131 -133.
[3] Yang S P,Li S H,et al.A hysteresis model for magnetorheological damper[J] ,Int J Nonlinear Sci Numer Simul,2005,6(2):139 -144.
[4] Ikhouane F,Rodellar J.Systems with hysteresis:analysis,identification and control using the Bouc-Wen model[M] .Chichester:John Wildy& Sons,Ltd,2007.
[5] Sireteanu T,Giuclea M,Mitu A M.An analytical approach for approximation of experimental hysteretic loops by Bouc-Wen model[J] .Proceedings of the Romanian Academy,Series A,2009,10(1):1 -12.
[6] Spencer B F,Dyke SJ,Sain M K,et al.Phenomenological model of a agnetorheological damper[J] .J.Engrg.Mech.,ASCE,1997(123):230-238.
[7] Dyke SJ.Acceleration feedback control strategies for active and semi-active control systems:modeling,algorithm development,and experimental verification[D] .University of Notre Dame,Indiana,1996.
[8] 歐進萍.結構振動控制——主動、半主動和智能控制[M] .北京:科學出版社,2003.
[9] Shen Y. Vehicle suspension vibration control with magnetorheological dampers[D] .University of Waterloo,Canada,2005.
[10] Liu Y Q,Matsuhisa H,Utsuno H,et al.Controllable vibration of the car-body using magnetorheological fluid damper[J] .Vehicle System Dynamics Supplement,2004(41):627-636.
[11] Jansen L M,Dyke SJ.Semiactive control strategies for MR dampers:comparative study[J] .Journal of Engineering Mechanics,2000(8):795 ~803.
[12] Wen Y K.Method of random vibration of hysteretic systems[J] .Journal of Engineering Mechanics Division,ASCE,1976,102(EM2):249 -263.
[13] Kristinsson K.System identification and control using genetic algorithms[J] .System,Man and Cybernetics,1992,22(5):1033-1046.
[14] Schneider T R.A genetic algorithm for the identification of conformationally invariant regions in protein molecules[J] .Acta Crystallographica Section D:Biological Crystallography,2002,58(2):195 -208.