復合材料板簧非線性遲滯建模和參數識別

史文庫 劉鶴龍 劉巧斌 陳志勇 閔海濤

摘要：為了描述和預測復合材料板簧的非線性遲滯特性，通過動力學分析，將滯回力分解為彈性力、阻尼力和殘余非線性力三部分組成，分別建立復合材料板簧的線性和非線性滯回力模型。進行復合材料板簧變頻和變幅激勵下的臺架試驗，以模型計算和試驗數據的均方根誤差為適應度函數，引入與迭代次數相關的非線性遞減權重系數，對螢火蟲算法進行改進，并與標準遺傳算法（GA）、標準模擬退火算法（SA）、標準螢火蟲算法（FA）進行比較，驗證了所提出的改進螢火蟲算法（MFA）具有相對高的精度和效率。將參數識別的結果分別代入線性模型和非線性模型，通過與試驗數據的對比，說明了非線性模型具有更好的適用性，為復合材料板簧的整車性能匹配和動力學分析奠定了基礎。

關鍵詞：遲滯非線性;復合材料板簧;改進螢火蟲算法;臺架試驗;幅變特性

中圖分類號：0322; U463.33+4

文獻標志碼：A

文章編號：1004-4523 （2019） 06-1050-10

DOI：10. 16 385/j. cnki. issn. 1004-4523. 2019. 06. 014

引言

節能與環保的要求日益苛刻，使得汽車輕量化技術受到了廣泛的關注，而復合材料在汽車輕量化應用上表現出了巨大的潛力，對復合材料在汽車上的應用進行深入的研究具有重要的工程價值。復合材料板簧是復合材料在汽車輕量化領域的一個應用，與傳統的鋼板彈簧相比，復合材料板簧在降低質量節省材料的同時，也降低了汽車的非簧載質量。非簧載質量的減小有助于減小車輪的動載荷，減小了車輪由于動載荷的作用而抓地能力下降的風險。因此，輕質的復合材料板簧的使用不僅可以提升汽車的燃油經濟性，同時也增加了汽車的行駛安全性。

國內外學者對復合材料板簧進行了大量的研究。主要研究集中于結構參數的優化[12]、剛度和強度的有限元分析[3-5]、制造工藝[6-7]、剛度計算理論[8-11]、疲勞損傷分析與優化[12-17]等方面。現有的研究對于復合材料板簧的遲滯特性研究的較少，而其遲滯模型在整車性能仿真方面不可或缺。因此，建立簡潔、高效的遲滯特性模型對復合材料板簧的遲滯特性進行仿真分析具有重要的意義。

傳統的鋼板彈簧遲滯非線性模型有：SAE三連桿模型[18]、有限元模型[19]、Fancher模型[20]和Bou-W en模型[21]等。復合材料是一種具有黏彈性非線性特征的材料，其遲滯特性與傳統的鋼板彈簧一樣具有幅變和頻變等特征，而上述模型對于復合材料板簧遲滯特性的建模是否適用，是否存在更簡潔實用的模型可以用于復合材料板簧的遲滯特性建模將是本文研究的一個重點。另一方面，在非線性遲滯模型的參數識別方法上，已有的研究集中于進化算法[22]、多目標遺傳算法[23]、粒子群算法[24]、人工神經網絡算法[25]、非線性狀態空間方法[26]、馬爾科夫鏈蒙特卡洛方法[27]等方法，由此可見，智能算法在非線性模型的參數識別方面展示出了極大的優勢，因此探索新型智能啟發式算法在非線性遲滯模型參數識別方面的應用是本文的另一個重點。

結合國內外研究背景和課題組前期的相關研究成果，本文在分析目前遲滯非線性建模及模型參數識別的基礎上，提出一種簡潔的非線性模型用于模擬復合材料板簧的遲滯非線性，并根據試驗測得的遲滯曲線，采用改進的螢火蟲算法對模型的參數進行識別。

合材料板簧臺架試驗

本文所研究的復合材料板簧為E玻璃纖維一聚氨酯基復合材料制作而成，采用單片拋物線結構，復合材料板簧的樣件如圖1所示。對復合材料板簧樣件進行臺架試驗，分別測試不同力激勵幅值和不同激勵頻率下的力一位移滯回曲線。試驗臺架如圖2所示。試驗時，預載力為11600 N。激勵頻率為0. 5-3 Hz，每隔0.5 Hz測取一組數據，共6組不同位移激勵幅值。最小力加載幅值為700 N，最大加載幅值為3500 N，每隔700 N測量一組數據，共5組數據。因此，本次試驗共采集到不同力加載幅值和不同激勵頻率下的數據共3 0組。

圖3所示為0.5 Hz激勵頻率下不同力加載幅值的滯回曲線，由圖可知隨著力激勵幅值的增加，滯回曲線產生明顯的變化，加載/卸載過程的力所形成的封閉曲線包圍的面積加大。

圖4所示是在700 N力加載幅值的條件下，不同激勵頻率下的復合材料板簧滯回曲線，由圖可知，在試驗頻率范圍內，其滯回曲線變化不明顯。因此，為了簡化研究，在后文的研究中，忽略頻率對復合材料板簧滯回特性的影響，僅考慮力幅值對滯回特性的影響。

2 復合材料板簧遲滯特性建模

2.1 線性模型

復合材料板簧在受到外界載荷作用時，發生形變，產生回復力，回復力中的一部分與外界激勵的位移相關，將這一部分回復力定義為彈性力，另一部分與外界激勵的速度相關，將這一部分回復力定義為阻尼力，可以獲得如下式所示的復合材料板簧線性遲滯模型。

F =k（A，f）z+c（A，f）x+F。

（1）式中 x為位移激勵，x為位移激勵的導數，即速度，A為激勵力幅值，廠為激勵頻率，k（A，f）為彈性剛度系數，c（A，f）為阻尼系數，F0為預載荷，在線性模型中，假定F0為常數。

復合材料板簧用于支撐汽車的簧載質量，而汽車簧載質量的偏頻在1-2 Hz以內，因此，在進行復合材料板簧臺架試驗時，激勵頻率范圍取0. 5-3Hz。由復合材料板簧臺架試驗可知，在所研究的激勵頻率范圍內，頻率對復合材料的滯回特性的影響很小，為了研究的方便，忽略復合材料板簧滯回特性的頻變特性，將式（1）簡化為如下式所示

F =k（A）x-+ c（A）x+F0（2）

對于每一個力振幅激勵下的滯回曲線，線性滯回模型存在2個未知參數，即k和e。

2.2 非線性模型

根據矢量的可疊加性，對式（2）的線性滯回模型進一步完善，引入殘余非線性項Fr，建立復合材料板簧的非線性滯回特性模型，如下式

F =k（A）x-+ c（A）x-+Fr（x，x，A）

（3）式中 Fr（x，x，A）為殘余非線性力，是位移、速度和振幅的非線性函數，對于大黏彈性系統，可用Bouc-W en模型等方法對殘余非線性力進行擬合，而Bouc-Wen模型的參數較多，且參數識別的初始值不容易確定。故本文在經過對線性模型擬合后的殘余項進行分析發現，對于本文所研究的復合材料板簧，其滯回力的非線性殘余項可表示為預載力的倍數δ（A），即

Fr（x，x，A）=δ（A）F0（4）

將式（4）代人式（3），可得

F =k（A）x+c（A）x十δ（A）Fo0（5）

對于每一個力振幅激勵下的滯回曲線，非線性滯回模型存在3個未知參數，即k，c和δ。

2.3 模型參數影響分析

為分析滯回特性模型的關鍵參數對滯回環的影響，分別對k，c和δ這3個參數進行影響分析。分析中，3個參數的初始值分別為k0，C0和δ0，參數影響分析采用控制變量法，固定2個參數，每次只改變1個參數值。圖5所示是剛度參數k對滯回特性的影響，由圖可知，隨著剛度的加大，滯后曲線整體斜率增大;圖6所示是阻尼參數c對滯回特性的影響，由圖可知，隨著阻尼的增大，滯回環所占面積增大，表示復合材料板簧在加載和卸載過程中的能量損耗增加。

圖7所示是殘余力系數δ對滯回環的影響，由圖可知，隨著殘余力系數的增大，滯回環整體延y軸下移，應該指出的是非線性殘余力系數δ中包含了非線性彈性力、非線性阻尼力、預載力和回復力測量誤差等幾個部分，而δ（A）僅是在數值上將非線性殘余力等效為預載力的倍數，因此圖7的分析并不能完全體現非線性模型中δ（A）對曲線的影響，由于上述非線性彈性力和非線性阻尼力的存在，δ（A）必然對滯回環的斜率和包圍面積有影響。

3 基于改進的螢火蟲算法的參數識別

參數識別是遲滯特性建模的核心之一，其精度和效率直接關系到建模的成敗。現有的非線性最小二乘法等傳統算法用于參數識別常有初始值不易確定、收斂速度慢和目標函數確定困難等弊端。新興的智能算法在非線性參數識別領域展現了巨大的潛力，各種算法不斷涌現，算法的適用性在解決實際工程問題中不斷得到驗證，各算法也針對不同的問題得到了不斷的改進和完善。本文嘗試將螢火蟲算法這種智能算法用于所提出的非線性遲滯模型的參數識別，并有針對性的對算法進行了改進，在提升了算法收斂精度的同時，保證了算法的運行效率。為驗證改進的螢火蟲算法的相對優越性，將改進的螢火蟲算法與遺傳算法、模擬退火算法這兩種相對成熟的智能算法比較;另外，也與改進前的螢火蟲算法做了比較。

3.1 參數初始值和取值范圍的確定

采用人工智能算法進行參數識別時，參數的初始值和取值范圍關系到參數識別的效率和精度，合適的初始值和取值范圍可大幅提升參數識別結果的可信度，減小參數識別的時間成本。對于非線性和線性模型中涉及到的3個參數k，c和δ，k的初始值設為滯回特性曲線一次線性回歸結果的斜率，如圖8所示，即k0=130 N/mm。c的初始值根據下式計算式中 △E為滯后環所包圍的面積，E為加載過程所包圍曲線的面積，這兩個面積表征的是能量，f為位移激勵的頻率。

為分析方便，取最小激勵振幅下的滯回環，即700 N力激勵下的位移一力滯回曲線表征的阻尼值為初始值，通過計算得到c0=5 N·s/mm。通過分析各組不同振幅下的試驗數據，發現非線性殘余力與預載力基本相當，故將非線性殘余項系數δ的初始值設為δ0=1。

各參數的取值范圍為

3.2 改進螢火蟲算法

螢火蟲算法（Firefly Algorithm，FA）是2 008年由劍橋大學學者Yang Xinshe所提出的一種新型啟發式群智能優化算法[28-30]，其基本思想借鑒了螢火蟲通過發光而實現信息溝通的功能。通過發光的隨機粒子之間距離的迭代計算從而實現函數尋優。

從數學的角度對螢火蟲算法相關概念的定義如下：

螢火蟲的亮度定義如下式

I（rij）=I0e-rr2j

（8）式中 IO為螢火蟲的初始亮度，γ為光強吸收系數，rij為第i個和第j個螢火蟲之間的笛卡爾距離，如下式所示式中 xi和xi分別為第i和j個螢火蟲的坐標，d為坐標的維度。

螢火蟲之間的吸引度定義如下式所示

β=β0e-yr2ij

（1 0）式中 β0為螢火蟲的初始吸引度。

亮度較小的螢火蟲朝較亮螢火蟲的位置移動的位置更新公式如下式所示

xi（t+1）=xi（t） +β（Xi （t） -xj （t））+

a（rand-0.5）

（11）式中 t為迭代次數，a為步長因子，rand為O到1之間的隨機數。

式（11）的位置更新迭代過程中，容易出現早熟而陷入局部最優，且在迭代后期螢火蟲聚集在最優點附近，螢火蟲之間的吸引力過大，可能導致解的震蕩。為解決以上問題，學者們提出了很多種改進方法[31-39]，這些方法實現的目的就是在保證算法的全局搜索能力的同時，提升算法的局部尋優能力，希望在迭代初期，算法具有較強的全局搜索能力，而在迭代后期，算法具有較強的局部搜索能力。現有文獻報道中對螢火蟲算法的改進基本都是線性遞減權重的，很少有非線性遞減權重的報道，本文根據所優化問題的特性，經過大量的測試，發現使用1/2次根值遞減的權重可實現更好的優化效果。

改進后的位置更新公式如下式

xi（t+1）=φ（t）xi（t）+β（xi（t）-xi（t））+

a（rand-0.5）

（12）式中 φ（t）為螢火蟲當前迭代位置的權重，與迭代次數成非線性遞減關系，如下式所示

由式（13）獲得的權重隨迭代過程的變化曲線如圖9所示。由圖可知，在迭代的初始階段，權重φ（t）較大，這保證了初始迭代時算法的全局尋優能力，而在迭代的終了階段，權重φ（t）較小，這保證了算法在后期的局部搜索能力。

3.3 適應度函數的確定

將適應度函數定義為模型計算的回復力Fsi與試驗測得的回復力F矗之間相對誤差平方和的均值，即如下式所示的均方誤差（RootMean Square Er-ror，RMSE）式中 X為待識別的參數變量，n為參與識別的樣本點總數，i為樣本點序號。參數識別的目的是在可行域內獲得一組使得目標函數值最大（RMSE最小）的變量。

4 參數識別結果與討論

4.1 改進螢火蟲算法參數識別的優點

為驗證所提出的算法進行參數識別的優越性，分別采用改進螢火蟲算法（MFA）、標準遺傳算法（GA）、標準模擬退火算法（SA）、和標準螢火蟲算法（FA）進行0.7 kN激勵力幅值T的非線性滯回力模型的3個關鍵參數識別。4種算法的種群數量統一設置為30，最大迭代步數為80。螢火蟲算法和標準螢火蟲算法的初始亮度I0=1，初始吸引度β0=1，步長因子a=0.5，光強吸收系數r設為0.9。標準遺傳算法的交叉概率為0.7，變異概率為0.1，采取輪盤賭法選擇機制。模擬退火算法的參數設置按照MATLAB最優化工具箱的默認設置。

分別記錄以上4種算法的運行時間和最終適應度值，如表1所示。由表可知，模擬退火算法在運行時間上略有優勢，但是其優化效果較差。標準遺傳算法和標準螢火蟲算法在優化效果上相近，標準螢火蟲算法的運行時間稍長，而經過改進后的螢火蟲算法，無論在求解精度上，還是在運行效率上，都得到了一定的提升，其優化精度優于標準遺傳算法，均方根誤差控制在1%以內，且運行效率略高于標準遺傳算法。

圖1 0所示是這4種算法隨迭代次數的收斂曲線，由圖可以直觀地看出各種算法的精度和效率。由圖可知，改進以后的螢火蟲算法的精度最高，改進后的螢火蟲算法在初始迭代階段，其適應度函數呈現了階梯狀的上升，而標準粒子群算法出現了早熟，因而陷入局部最優。

4.2 參數識別最終結果

在驗證了改進螢火蟲算法的優越性后，采用改進的螢火蟲算法分別對各力振幅激勵下的線性滯回曲線和非線性滯回曲線進行參數識別。

表2和3分別為采用改進的螢火蟲算法對不同力載荷下的線性模型和非線性模型的參數進行識別的結果。由表2可知，隨著力激勵幅值的增加，線性模型的誤差逐漸增大，而線性模型的等效剛度和阻尼略有減小。對比表2和3的結果可知，采用非線性模型后，滯回力的平均預測誤差的最大值從7. 8133%下降到1.3866%，而參數識別的時間基本不變。說明引入非線性殘余力項可以大大提高模型的精度。圖11所示是不同預載力下的位移一滯回力和時間一滯回力曲線。由圖可知，非線性模型與實測數據吻合的很好，而線性模型誤差較大，且激勵力幅值越大，線性模型的誤差也越大，說明大幅值的激勵力增大了復合材料板簧的非線性。

在驗證了本文所提出的非線性滯回模型的準確性的基礎上，進一步對非線性模型中3個參數的幅變特性進行分析，發現使用二次多項式可以實現這3個參數的良好擬合，下式是擬合結果

圖12 -14所示，分別是剛度k，阻尼c和非線性殘余力系數δ的二次多項式擬合結果和基于改進的螢火蟲算法參數識別結果的對比。由圖可知，隨著力振幅的增大，復合材料板簧的滯回模型中的剛度、阻尼和非線性殘余項均有所下降，說明隨著力振幅的增大，板簧的非線性增加，耗散能量增加。

5 結 論

本文在分析復合材料板簧滯回力物理意義的基礎上，將滯回力分解為線性彈性力、線性阻尼力和非線性殘余力三部分，建立起了復合材料板簧滯回力的非線性模型，并采用改進螢火蟲算法對非線性模型中的模型參數進行識別，將識別結果與線性模型和實測結果對比，驗證了所提出的非線性模型具有很好的適用性，且參數識別方法具有相對較高的效率和精度。主要結論如下：

（1）復合材料板簧的滯回力存在非線性因素，其滯回特性的幅變特性明顯，而頻變特性不明顯，采用線性模型對滯回力進行擬合的誤差較大;

（2）采用本文所提出的非線性滯回力模型可以較好地描述和預測復合材料板簧的滯回特性，且效率和建模精度都相對較高;

（3）通過位置坐標權重非線性改進后的螢火蟲算法參數識別方法，與經典螢火蟲算法、經典粒子群算法和經典遺傳算法相比，具有更好的精度和相對較高的效率，可以為相關非線性模型的參數識別提供參考;

（4）進一步的研究要通過對復合材料板簧的彈塑性力學機理進行深入的探索，分析滯回特性中的非線性殘余項與位移和速度等物理量的關系;另一方面，采用不同的參數識別方法對模型的參數進行識別，以進一步提高參數識別的效率和精度。

參考文獻：

[1] Rajendran I， Vijayarangan S.Optimal design of acomposite leaf spring using genetic algorithms[J].Computers & Structures， 2001， 79（11）：1121-1129.

[2] Sancaktar E，Gratton M. Design， analysis， and opti-mization of composite leaf springs for light vehicle ap-plications[J]. Composite Structures， 1999， 44（2-3）：195-204.

[3] Mithari R，Patil A， Aitavade E N.Analysis of com-posite leaf spring by using analytical&FEA[J]. In-ternational Journal of Engineering Science & Technol-ogy， 2012， 4（12）：4809-4814.

[4] Shokrieh M M， Rezaei D.Analysis and optimization ofa composite leaf spring[J]. Composite Structures，2003， 60（3）：317-325.

[5] Hou J P，Cherruault J Y， Nairne I，et al.Evolution ofthe eye-end design of a composite leaf spring for heavyaxle loads[J]. Composite Structures， 2007， 78（3）：3 51- 358.

[6] Banka H， Muluka R，Reddy V. Fabrication and experi-mental analysis of epoxy-glass fiber composite leaf spring[R]. SAE Technical Paper Series， 2017-28-1985.

[7]Rai V， Saxena G. Development of a composite leafspring for a light commercial vehicle （Tata Magic）[J].International Journal of Engineering Research& Ap-plications， 2013， 3（5） ：110-114.

[8]Jadhao K K， Dalu D R S. Experimental investigation& numerical analysis of composite leaf spring[J]. In-ternational Journal of Engineering Science& Technol-ogy， 2011， 3（6） ：445-477.

[9] Shankar G S S， Vijayarangan S. Mono composite leafspring for light vehicle-design， end joint analysis andtesting[J]. Materials Science， 2006， 12（3） ：220-225.

[10] Azzam B S N， Mokhtar M O A， Naga S A R. Theo-retical and numerical analysis of fibrous composite C-springs[R]. SAE Technical Paper Series， # 2001-01-2 710.

[11] Shi W， Qian C， Chen Z， et al. Establishment of theo-retical model of composite leaf springs by using themechanics of composite materials[J]. Journal of Rein-forced Plastics& Composites， 2017， 36（18）：1316-1326.

[12] Kumar M S， Vijayarangan S. Static analysis and fa-tigue life prediction of steel and composite leaf springfor light passenger vehicles[J]. Journal of Scientific&Industrial Research， 2007， 66（2） ：128-134.

[13] Kueh J T J， Faris T. Finite element analysis on thestatic and fatigue characteristics of composite multi-leaf spring[J]. Journal of Zhejiang University-ScienceA（Applied Physics& Engineering）， 2012， 13（3）：15 9-16 4.

[14] Manjunath H N， Manjunath K， Rangaswamy T. Stat-ie analysis and fatigue life prediction of composite leafspring for a light commercial vehicle （TATA ACE）[J]. International Journal of Engineering Research，2014， 3（7） ：422-425.

[15] Shamsaei N， Rezaei D. Comparing fatigue life reliabil-ity of a composite leaf spring with a steel leaf spring[C]. ASME， Biennial Conference on Engineering Sys-tems Design and Analysis， 2004：371-374.

[16] Kumar M S， Vijayarangan S. Analytical and experi-mental studies on fatigue life prediction of steel andcomposite multi-leaf spring for light passenger vehiclesusing life data analysis[J]. Materials Science， 2007，13（2） ：141-146.

[17] Qian C， Shi W， Chen Z， et al. Fatigue reliability designof composite leaf springs based on ply scheme optimization[J]. Composite Structures， 2017， 168： 40-46.

[18]官欣，齊海政，詹軍，等，基于SAE三連桿模型的鋼板彈簧遲滯特性建模[J].吉林大學學報（工學版） ，2010， 40（S1） ：76-0080.

Guan X， Qi H Z， Zhan J， et al. Modeling for leafspring hysteretic characteristic based on SAE three-link model[J]. Journal of Jilin University（EngineeringEdition） ， 2010， 40（S1） ：76-0080.

[19]周孔亢 ，陸建輝 ，侯永濤 ，等.基于 RecurDyn的鋼板報， 2014， 50（4） ：128-134.

Zhou Kongkang， LU Jianhui， Hou Yongtao， et al.Dynamics modeling and parameter identification of leafspring based on recurdyn[J]. Journal of MechanicalEngineering， 2014， 50（4） ：128-134.

[20] Liu X， Chen W， Zhao J， et al. Fitting precision analy-sis on the dynamic stiffness of leaf spring based on theFancher Model[J]. Automobile Technology， 2017，（3） ：38-41.

[21] Xie Q X， Zhou S Y， Zhang N， et. al. Simulation andexperimental study on Bouc-Wen hysteresis model ofleaf spring[J]. Automobile Engineering， 2016， 38（5） ：632-637.

[22] Charalampakis A E， Koumousis V K. Identification ofBouc-Wen hysteretic systems by a hybrid evolutionaryalgorithm[J]. Journal of Sound& Vibration， 2008，314（3） ：571-585.

[23] Ortiz G A， Alvarez D A， Bedoya-Ruiz D. Identifica-tion of Bouc-Wen type models using multi-objectiveoptimization algorithms[J]. Computers& Structures，2013 ， s114-115 （1） ：121-132.

[24] Charalampakis A E， Dimou C K. Identification ofBouc-Wen hysteretic systems using particle swarm op-timization[J]. Computers& Structures， 2010， 88（ 21） ：119 7-12 0 5.

[25] Xie S L， Zhang Y H， Chen C H， et al. Identificationof nonlinear hysteretic systems by artificial neural net-work[J]. Mechanical Systems& Signal Processing，2013， 34（1-2） ：76-87.

[26] Noel J P， Esfahani A F， Kerschen G， et al. A nonlin-ear state-space approach to hysteresis identification[J]. Mechanical Systems& Signal Processing， 2016，84 ：171-184.

[27] Ortiz G A， Alvarez D A， Bedoya-Ruiz D. Identifica-tion of Bouc-Wen type models using the TransitionalMarkov Chain Monte Carlo method[J]. Computers&Structures， 2015， 146：252-269.

[28] Yang X S. Firefly algorithms for multimodal optimiza-tion[J]. Mathematics， 2012， 5792：169-178.

[29] Yang X S， He X. Firefly algorithm： Recent advancesand applications[J]. International Journal of SwarmIntelligence， 2013， 1（1） ：36-50.

[30] Fister I， Jr I F， Yang X S， et al. A comprehensive re-view of firefly algorithms[J]. Swarm& EvolutionaryComputation， 2013 ， 13 （1） ：34-46.

[31] Hong CO ， Tilahun S L， Tang S S. A comparativestudy on standard， modified and chaotic firefly algo-rithms[J]. Pertanika Journal of Science& Technolo-gy， 2015， 23（2） ：251-269.

[32] Fister I， Yang X S， Brest J， et al. Modified firefly al-gorithm using quaternion epresentation [J]. ExperSystems with Applications， 2013， 40（18） ：7220-7230.

[33] Niknam T， Azizipanah-Abarghooee R， Roosta A. Re-serve constrained dynamic economic dispatch： A newfast self-adaptive modified firefly algorithm[J]. IEEESystems Journal， 2012 ，6 （4） ：635-646.

[34] Verma O P， Aggarwal D， Patodi T. Opposition and di-mensional based modified firefly algorithm[Jl. ExpertSystems with Applications， 2016， 44（C） ：168-176.

[35] Xiao L， Wang C， Liang T， et al. Combined modelbased on multiple seasonal patterns and modified fire-fly algorithm for electrical load forecasting[J]. Ap-plied Energy， 2016， 167：135-153.

[36] Li M， Zhang Y， Zeng B， et al. The modified fireflyalgorithm considering fireflies' visual range and its ap-plication in assembly sequences planning[J]. Interna-tional Journal of Advanced Manufacturing Technology，2016 ， 82 （ 5-8） ：1381-1403.

[37] Wang B， Li D X， Jiang J P， et al. A modified fireflyalgorithm based on light intensity difference[J]. Jour-nal of Combinatorial Optimization， 2016， 31（3） ：1045-10 60.

[38] Yelghi A， Kose C. A modified firefly algorithm forglobal minimum optimization[J]. Applied Soft Com-puting， 2018， 62：29-44.

[39] Sedaghati R， Hakimzadeh M， Fotoohabadi H， et al.Study of network reconfiguration in distribution sys-tems using an adaptive modified firefly algorithm[J].Automatika， 2016， 57（1） ：27-36.