基于混合螢火蟲算法的橋梁顫振分析方法

陶仕博， 湯愛平， 胡慶杰, 劉克同

(1.哈爾濱工業大學 結構工程災變與控制教育部重點實驗室，哈爾濱 150090； 2.哈爾濱工業大學 土木工程學院，哈爾濱 150090; 3.西安科技大學 建筑與土木工程學院, 西安 710054)

基于混合螢火蟲算法的橋梁顫振分析方法

陶仕博1,2， 湯愛平1,2， 胡慶杰1,2, 劉克同3

(1.哈爾濱工業大學 結構工程災變與控制教育部重點實驗室，哈爾濱 150090； 2.哈爾濱工業大學 土木工程學院，哈爾濱 150090; 3.西安科技大學 建筑與土木工程學院, 西安 710054)

在進行顫振臨界狀態分析時需要求解高次非線性方程組，為了克服傳統解法的缺陷，采用混合螢火蟲算法對方程組進行求解。使用雙參數優化模型，將橋梁顫振臨界狀態的求解問題轉化為優化問題。為彌補螢火蟲算法的不足，在螢火蟲算法基礎上，將量子遺傳算法的量子計算、交叉和變異操作與螢火蟲算法相結合，提出一種混合螢火蟲算法。最后，通過若干試驗對比分析，證實了該優化模型的可靠性及求解方法的有效性。

顫振；優化模型；量子遺傳算法;混合螢火蟲算法；最優解

在橋梁風致振動中，顫振是一種典型的氣動彈性發散現象，它會給橋梁帶來災難性的后果。因此準確判斷顫振臨界狀態對于橋梁抗風設計有著非常重要的作用[1-2]。傳統的半逆解法是最常用的顫振臨界狀態分析方法，它將問題轉化為變系數齊次方程組系數行列式的特征值問題。由于組成方程各有無窮多組解，傳統的半逆解法(如P-K法或PK-F法等)需對兩個方程的根循環對比以求得共同解[3]，不僅計算復雜，計算量大，還需要人工干預來使計算結果合理。過多的人為設置不僅影響計算效率和精度，甚至影響算法的穩定，致使分析無法進行。因此，有必要研究更為高效的、自動程度更高的方法。

文獻[4]基于矩陣奇異值理論通過計算顫振矩陣最小奇異值或條件數的倒數搜索顫振臨界點。文獻[5]采用追趕法來搜索理想平板和橋梁節段模型的顫振臨界點。文獻[6]基于準定常氣動理論對顫振行列式進行簡化，給出了顫振臨界風速的表達式。這些方法雖然避免直接解強非線性方程組，但需人工設置初始搜索值，并且對初始搜索值依賴程度大。如果初始值設置不合適，不僅增加計算量，還會使算法陷入局部最小而無法搜索到顫振臨界點。

基于上述現狀，為了避免直接求解高次非線性方程組，將顫振臨界狀態的求解轉化為優化問題，建立了搜索顫振臨界點的雙參數優化模型。為了回避傳統梯度尋優算法需求偏導數、對初始搜索值敏感及傳統螢火蟲算法(Firefly Algorithm, FA)收斂慢、易陷入局部最優的缺點，借鑒量子遺傳算法(Quantum Genetic Algorithm, QGA)[7]的思想對螢火蟲算法進行改進，采用混合螢火蟲算法(Hybrid Firefly Algorithm，HFA)來搜索該優化模型最優解。最后，通過風洞試驗實例分析，證實了該優化模型的可靠性及求解方法的有效性。

1 顫振優化模型

在進行搜索前，需要把顫振運動方程改寫成雙參數優化模型[8-9]。橋梁豎彎、側彎和扭轉方向的氣動自激力L、D和M可以通過18個無量綱的顫振(氣動)導數的線性函數來表示

(1a)

(1b)

(1c)

將氣動自激力表達式(1)代入運動方程，引入無量綱時間s=tU/B，并令K=Bω/U，并設方程的解為

(2)

式中：m和I分別為主梁單位長度的質量和質量慣矩;ζh、ζp和ζα則分別為主梁豎彎、側彎和扭轉運動的阻尼比;ωh、ωp和ωα分別為豎彎、側彎和扭轉圓頻率;h、p和α為對應自由度的廣義位移。

方程組(2)有非零解的條件是系數行列式為0，即

(3)

式中,Aij為式(3)未知量對應的系數。其具體表達式可見文獻[10]。

將式(3)展開，得到一個關于X的六次復系數多項式。臨界狀態下X是實數，則該多項式的實部和虛部均應為0，由此整理得出臨界狀態顫振方程組的表達式為

(4)

顫振臨界狀態的求解實質上是在可行域尋找一組ω和U，使方程組(4)的實部和虛部同時為0，則該組ω和U即為待求臨界圓頻率ωcr和臨界風速Ucr。對于實際橋梁來說，在顫振臨界狀態，圓頻率ωcr處在max{ωα,ωh,ωp}和min{ωα,ωh,ωp}之間，臨界風速Ucr可由公式[11]

(5)

粗略估計后給出一個搜索范圍。

式中:ε=ωα/ωh為橋梁扭彎頻率比；r為回轉半徑,m，即橋梁慣性矩除以截面積所得商的平方根值；μ為結構相對空氣的質量比；b為半寬,m。

方程組(4)的求解可以轉化為這樣一個優化問題：

(6)

式(6)即為搜索臨界圓頻率和臨界風速的顫振優化模型。

2 顫振優化模型的求解

2.1 傳統螢火蟲算法

FA是模擬自然界中螢火蟲的行為而衍生出的啟發式全局優化算法，它利用螢火蟲發光特性在搜索空間中尋找伙伴，并向位置較優的螢火蟲移動，從而達到進化的目的。在該算法中螢火蟲彼此吸引的原因取決于兩個要素：自身發光亮度和吸引度。其中，螢火蟲發出熒光的亮度取決于自身所在位置的目標值(適應度)，亮度越高表示所處的位置越好，即目標值越佳。

在描述具體的 FA之前，需進行如下假設[12-15]:

(1)螢火蟲不分性別，發光較強的螢火蟲可以無差別吸引其他發光較弱的螢火蟲。

(2)每個螢火蟲的吸引度β與其發光強度I成正比。對于任意兩只螢火蟲，發光較弱的螢火蟲會朝發光較強的螢火蟲移動，且β與I隨著螢火蟲之間的距離r的增大而減小。最亮的螢火蟲(即I最大的螢火蟲)是隨機飛行的。如果螢火蟲發光亮度相同，則螢火蟲隨機移動。

(3)螢火蟲的發光強度I與目標函數值有關。

在滿上述三個假設的前提下，螢火蟲算法的基本步驟如下。

步驟1 設置算法參數。種群規模N，最大吸引度β0，吸收系數γ，隨機步長α，最大迭代次數。在解空間中隨機初始化螢火蟲的位置X，令t=0。

步驟2 每只螢火蟲的發光強度Ii(i=1，2，…，N)，Ii作為適應度 -f(Xi) (Xi表示問題的一組可能解)，即Ii=-f(Xi)，1≤i≤N。

步驟4 更新螢火蟲的位置。螢火蟲i被發光強度更亮j吸引而發生位置移動。位置移動方式由式(7)決定

(7)

步驟5 判斷算法是否滿足終止條件，若滿足，則算法結束，輸出最優解；否則，令t=t+1，返回步驟2。根據參考文獻[16]本文設置的螢火蟲算法參數，α=0.2，β0=1，γ=0.1。

2.2 螢火蟲算法的改進

FA已被證明在求解精度和穩定性方面都超過了許多其他進化算法，如粒子群優化算法，遺傳算法等。但由于自身缺乏變異機制，一旦受到局部極值的束縛將很難擺脫，特別是在進化初期，算法對初始解分布的依賴性較強。而在進化后期，由于收斂速度慢、求解精度低和易陷入局部最優等缺點。

量子遺傳算法是量子計算理論與遺傳算法結合的產物，此方法具有良好的魯棒性和廣泛的適應性。具體表現在，繼承了遺傳算法不受問題性質和優化準則形式等因素限制的優勢，有很強的變異機制，克服了遺傳算法收斂速度慢、易陷入局部極值和對初始值比較依賴等缺點。而且目標函數在概率引導下能夠自適應的進行全局搜索。因此本文作者采用量子遺傳算法的思想對傳統的螢火蟲算法進行改進。具體改進方式如下：

(1)以基于量子位的二進制編碼表示螢火蟲位置。一個量子位狀態可用基本量子比特狀態|0〉和|1〉的疊加描述

(8)

(9)

(2)將量子旋轉門的概念引入到螢火蟲算法中，使螢火中算法在保持原有局部尋優能力基礎上，進一步提高其全局搜索能力。采用量子旋轉門更新群體Q(t)，量子旋轉門的操作公式為

(10)

定義旋轉門更新策略如下

(11)

(3)增加交叉及變異操作,交叉及變異的目的是為了獲取新的信息，以保持種群的多樣性。交叉通過交換部分量子位的編碼實現。變異通過改變部分量子位的編碼實現。

改進后的螢火蟲算法的演化流程為

(1)初始化種群Q(t),t=0，產生量子種群。

(2)對初始種群Q(t)中的每個個體進行測量,以得到Q(t)的二進制解P(t)。

(3)對P(t)進行適應度評價，記算P(t)中的個體對應的發光強度。

(4)判斷是否滿足終止條件：是，終止；否，繼續。

(5)t=t+1；通過式(11)量子旋轉門更新策略改變種群位置，得到新種群Q′(t)。

(6)實施量子交叉、變異操作。

(7)計算個體對應的發光強度。

(8)將種群保存為Q(t)，轉至(4)。

2.3 顫振臨界狀態的HFA求解

根據上述演化流程，采用HFA搜索顫振臨界點，對下面關鍵步驟進行說明：

(1)初始化種群Q(t0)，令種群的螢火蟲位置為顫振優化模型的解向量[ω,U]。

(2)取式(6)中f(X)為適應度函數，計算其相對熒光亮度。

(3)量子旋轉門更新策略。旋轉角的調整策略采用式(11)來實現。

3 算例分析

本文用三個算例進行比較，實驗環境為Matlab，機器配置為AMD A8-4500M APU with Radeon(tm) HD Geaphics 1.90 GHz處理器，4 GB內存。三種算法采用相同種群規模100，交叉概率為0.2，變異概率為0.1，最大迭代次數為500。前后兩代平均適應度的差值小于10-6時，認為算法開始收斂。分別采用三種算法對優化模型運行100輪。

3.1 仿真分析(算例1)：二自由度—理想平板

均勻流作用下的薄平板，半寬度b=0.225 m，每延米質量m=11.25 kg/m，慣性矩I=0.282 8 kg·m2，豎彎頻率ωh=12.11 rad/s，扭轉頻率ωα=19.0 rad/s，豎彎阻尼比ζh=0.005，扭轉阻尼比ζα=0.008，氣動導數采用Theodorsen理論解[18]。采用顫振優化模型搜索顫振臨界頻率和臨界風速，ωcr∈[12,19]，Ucr∈[13,28]。

計算結果如表1所示。表中最大解是指連續進行100輪計算，取搜索結果中的最大值；相應的，最小解指的是100論計算中搜索結果中的最小值。而平均值指的是100輪計算所有搜索結果的平均值。由表1可知，雖然各種算法每次運行的演化過程不盡相同，但收斂后演化趨于一致。

表1 三種算法優化結果的對比(平板)

許福友等采用追趕法得到的計算結果為ω=15.186 rad/s，U=16.91 m/s，從表1中可以看出，HFA的計算結果較FA和QGA與其吻合更好。QGA存在不收斂的情況，FA和HFA所有運算全部收斂。表1中， HFA的收斂速度較FA有較明顯的提高。HFA每次運行得到的最優值基本收斂于平均值，QGA和FA仿真結果的離散性較大，尤其是QGA。在算例一中，HFA無論是在數值穩定性還是仿真精度上均優于QGA和FA。在計算所需要的時間上，HFA略大于FA及QGA。

3.2 算例2：二自由度自由振動顫振分析風洞試驗

為驗證HFA的有效性，選取我國某懸索橋主梁模型為研究對象。大橋主梁的節段模型按1∶40的縮尺比制作，其材料為有機玻璃，模型具體尺寸的如圖1所示。試驗完成地點為哈爾濱工業大學大氣邊界層風洞與浪槽聯合實驗室。試驗時沒有考慮欄桿和防撞墻等附屬物。

圖1 節段模型(mm)Fig.1 Section model (mm)

試驗時,模型有豎彎及扭轉兩個自由度,試驗風速為0 m/s、4 m/s、6 m/s、8 m/s、10 m/s、12 m/s、14 m/s及16 m/s。自由振動系統主要參數如表2所示。加速度信號由懸掛在節段模型四角上的4個加速度計采集。模型自由衰減振動信號經數據采集、濾波并進行適當的線性組合得到豎彎和扭轉振動的時程曲線。顫振導數識別采用MITD法[19]，識別結果如圖2所示。圖2中，f為橋梁振動頻率。優化算法搜索范圍，ωcr∈[0.6,1.3]，Ucr∈[60,120]。

圖2 氣動導數Fig.2 Flutter derivatives

采用傳統的Scanlan法[20]及QGA、FA和HFA三種優化算法計算臨界頻率及臨界風速，結果如表3所示。表3中可以看出，三種仿真算法計算結果的平均值相差不多，但QGA未收斂的情況最多，需要更多的迭代步數才能收斂。HFA收斂所需迭代次數最少，且計算全部收斂。HFA較QGA和FA在計算效率和穩定性方面表現出明顯的優勢。但在平均CPU時間上，HFA略大。

表2 模型系統的主要參數

表3 優化計算結果對比(二自由度)

3.3 算例3：三自由度強迫振動顫振分析

本試驗所選用的橋梁主梁模型為大貝爾特東橋主梁，其截面作為經典模型被廣泛應用于檢驗各種識別技術的正確性。其模型的縮尺比為1∶40。模型截面尺寸如圖3所示。剪切中心位于截面左右對稱線上，距橋面0.465倍截面高度處，試驗時沒有考慮欄桿和防撞墻等附屬物。

采用全解耦三自由度橋梁節段模型強迫振動系統，使用強迫振動法進行顫振試驗。顫振導數提取依然采用MITD法。優化算法搜索范圍，ωcr∈[0.5,1.2]，Ucr∈[35,80]。

圖3 強迫振動試驗橋梁尺寸(mm)Fig.3 The size of the model in the forced vibration tests(mm)

分別采用QGA，FA和HFA搜索最優解。取種群規模為150。搜索時，各算法參數取值同算例1。每種算法運行100輪。計算結果列于表4中。

表4 優化結果對比(三自由度)

從表4可以看出，QGA的開始收斂代數早于HFA和QGA。但由于三自由度問題復雜度增加，QGA及FA不收斂的情況很多。在平均CPU時間上，HFA需要的時間最長。但在數值穩定性方面，HFA要好于QGA和FA。

將HFA計算結果與前人研究結果的對比列于表5中。由于Poulsen風洞試驗考慮了橋面人行道欄桿和中央防撞欄等附屬物, 所以其臨界風速小于其他未考慮附屬設施的試驗結果對比。另外，HFA計算結果與表5中其他研究者的計算結果相差不大，這說明采用HFA算法搜索得到的結果是可以接受的。因此，對于二維三自由度顫振問題，HFA同樣具有較好的適用性。

表5 顫振臨界速度和臨界頻率

4 結 論

本文將量子遺傳算法與螢火蟲算法相結合，形成一種混合螢火蟲算法。并將其算法應用于搜索橋梁顫振臨界頻率及風速，主要結論如下:

(1) 利用量子遺傳算法的思想，對螢火蟲算法進行改進。采用二進制編碼，將螢火蟲算法中的種群量子化，并引入交叉變異操作，定義了旋轉角的調整策略，得到了一種混合螢火蟲算法，給出了采用該混合算法搜索顫振優化模型最優解的步驟。

(2) 采用改進的螢火蟲算法，螢火蟲算法和量子遺傳算法，對兩個二自由度顫振問題和一個三自由度顫振問題進行了分析。分析結果表明，本文的混合螢火蟲算法較螢火蟲算法和量子遺傳算法尋優結果更好，穩定性更強。算法所用的計算時間稍多于螢火蟲算法和量子遺傳算法。

(3)本文的混合螢火蟲算法計算過程中最大程度減少人為參與對計算結果的影響，實現了對顫振臨界點的全域自動搜索，不僅數值穩定性好，而且計算精度高，具有較好的實用價值。

[ 1 ] 王凱,廖海黎,李明水. 基于風洞試驗的大跨度鋼桁梁懸索橋顫振性能研究[J]. 振動與沖擊, 2015, 34(15):175-180. WANG Kai, LIAO Haili, LI Mingshui. Flutter stability study of a long-span suspension bridge with steel truss based on wind tunnel testing [J]. Journal of Vibration and Shock, 2015, 34(15): 175-180.

[ 2 ] MISHRA S S, KUMAR K, KRISHNA P. Multimode flutter of long-span cable-stayed bridge based on 18 experimental aeroelastic derivatives [J]. Journal of Wind Engineering & Industrial Aerodynamics, 2008, 96(1):83-102.

[ 3 ] 劉克同, 湯愛平,劉玥君,等.基于浸沒邊界和格子玻爾茲曼的流固耦合算法[J].華中科技大學學報(自然科學版),2015, 43(1):61-66. LIU Ketong, TANG Aiping, LIU Yuejun, et al. Fluid-structure interaction method using immersed boundary and Lattice Boltzmann method [J]. Journal of Huazhong University of Science and Technology (Natural Science Edition), 2015, 43(1):61-66.

[ 4 ] 谷迎松,楊智春,王巍.一類基于顫振行列式的顫振分析新方法[J].航空動力學報, 2009, 24(4):815-818. GU Yingsong, YANG Zhichun,WANG Wei. New type of flutter solution based on flutter determinant [J]. Journal of Aerospace Power，2009, 24(4): 815-818．

[ 5 ] 許福友, 陳艾榮,王達磊．搜索平板顫振臨界風速的追趕法[J].工程力學,2005, 22(5): 48-53. XU Fuyou, CHEN Airong, WANG Dalei. Pursuit method for searching critical flutter wind velocity of flat plates [J]. Engineering Mechanics, 2005, 22(5):48-53.

[ 6 ] LEE H E, VU T V, YOO S Y, et al. A simplified evaluation in critical frequency and wind speed to bridge deck flutter [C]//Procedia Engineering: The Twelfth East Asia-Pacific Conference on Structural Engineering and Construction,2011: 1784-1790.

[ 7 ] 楊宇,李紫珠,何知義,等. QGA-VPMCD智能診斷模型研究[J]. 振動與沖擊, 2015, 34(13): 31-35. YANG Yu, LI Zizhu, HE Zhiyi, et al. The intelligent diagnosis model research of QGA-VPMCD[J]. Journal of Vibration and Shock, 2015, 34(13): 31-35.

[ 8 ] 許福友, 陳艾榮, 梁艷. 平板顫振參數研究[J].同濟大學學報(自然科學版), 2005, 32(10):1365-1370. XU Fuyou, CHEN Airong, LIANG Yan. Research on flutter parameters of flat plate[J]. Journal of Tongji University (Natural Science), 2005, 32(10):1365-1370.

[ 9 ] 許福友.橋梁結構顫振導數識別與顫振分析[D].上海:同濟大學，2006.

[10] 許福友, 陳艾榮,張哲, 等. 確定橋梁模型顫振臨界風速實用方法[J].振動與沖擊,2008, 27 (12):97-100. XU Fuyou, CHEN Airong, ZHANG Zhe, et al. Practical technique for determining critical flutter wind speed of bridge model [J].Journal of Vibration and Shock, 2008, 27(12):97-100.

[11] PUT V D. Rigidity of structures against aerodynamic forces [C]//10th Congress of IABSE. Zurich: International Association of Bridges and Strutural Engineering，1976.

[12] 黃凱, 周永權. 帶交尾行為的混沌人工螢火蟲優化算法[J]. 計算機科學, 2012, 39(3):231-234. HUANG Kai, ZHOU Yongquan. Chaotic artificial glowworm swarm optimization algorithm with mating behavior[J]. Computer Science, 2012, 39(3): 231-234.

[13] 劉長平，葉春明. 具有混沌搜索策略的螢火蟲優化算法[J].系統管理學報,2013, 22(4):538-543. LIU Changping, YE Chunming. Firefly algorithm with chaotic search strategy[J]. Journal of Systems & Management, 2013, 22(4): 538-543.

[14] KRISHNANAND K N, GHOSE D. Glowworm swarm optimization for simultaneous capture of multiple local optima of multimodal function[J]. Swarm Intelligence, 2009, 3(2): 87-124.

[15] YANG X S. Nature-Inspired Metaheuristic Algorithms[M]. Luniver Prsee,2010.

[16] ZOUACHE D, NOUIOUA F, MOUSSAOUI A. Quantum-inspired firefly algorithm with particle swarm optimization for discrete optimization problems[J]. Soft Computing, 2015:1-19.

[17] RIGAS I, ECONOMOU G, FOTOPOULOS S. Efficient modeling of visual saliency based on local sparse representation and the use of hamming distance[J]. Computer Vision and Image Understanding, 2015, 134: 33-45.

[18] THEODORSON T. General theory of aerodynamic instability and the mechanism of flutter [R]. NACA Report,1935.

[19] KORLIN R, STAROSSEK U. Wind tunnel test of an active mass damper for bridge decks [J]. Journal of Wind Engineering and Industrial Aerodynamics, 2007,95(4): 267-277.

[20] 陳政清. 橋梁風工程[M]. 北京：人民交通出版社，2005.

[21] POULSEN N K, DAMSGAARD A, REINHOLD T A. Determination of flutter derivatives for the Great Belt Bridge [J]. Journal of Wind Engineering and Industrial Aerodynamics,1992, 41 (1/2/3):153-164.

[22] LARSEN A, WALTHER J H. Aeroelastic analysis of bridge girder sections based on discrete vortex simulations [J]. Journal of Wind Engineering and Industrial Aerodynamics, 1997, 67:253-265.

[23] LARSEN A, WALTHER J H. Discrete vortex simulation of flow around five generic bridge deck sections[J]. Journal of Wind Engineering and Industrial Aerodynamics, 1998, 77:591-602.

[24] 祝志文，陳政清. 數值模擬橋梁斷面氣動導數和顫振臨界風速[J].中國公路學報, 2004, 17(3):41-50. ZHU Zhiwen, CHEN Zhengqing. Numerical simulations for aerodynamic derivatives and critical flutter velocity of bridge deck [J]. China Journal of Highway and Transport，2004, 17(3): 41-50.

[25] 崔益華,陳國平. 橋梁顫振導數的耦合強迫振動仿真識別[J]. 振動工程學報, 2007,20(1):35-39. CUI Yihua, CHEN Guoping. Simulations for identification of flutter derivatives of bridge section using the coupled-forced-vibration method [J]. Journal of Vibration Engineering, 2007, 20 (1):35-39.

Flutter analysis for bridges based on a hybrid firefly algorithm

TAOShibo1, 2,TANGAiping1,2,HUQingjie1,2,LIUKetong3

(1. Key Laboratory of Structures Dynamic Behavior and Control Ministry of Education, Harbin Institute of Technology, Harbin 150090, China; 2.School of Civil Engineering, Harbin Institute of Technology, Harbin 150090, China; 3.College of Architecture and Civil Engineering, Xi’an University of Science and Technology, Xi’an 710054, China)

When performing flutter analysis, high-order strong nonlinear equations need to be solved. For overcoming the difficulties encountered by traditional methods, a hybrid firefly algorithm was used to solve the equations. The solution of the critical flutter state problem was converted to an optimization problem by using a double-parameter optimization model. Therefore the optimization model can be employed to calculate the critical velocity and the critical frequency of two or three degrees of freedom flutter. To compensate shortcoming of the firefly algorithm, a hybrid firefly algorithm was proposed and used for searching the optimal solution of the optimization model. Finally, the reliability and the validity of the optimization model as well as its solution were confirmed by numerical and experimental examples.

flutter; optimization model; quantum genetic algorithm; hybrid firefly algorithm; the optimal solution

國家高技術研究發展計劃(863計劃)項目(2008AA11Z104)；國家國際科技合作項目(2011DFA21460)

2015-09-23 修改稿收到日期：2015-12-04

陶仕博 男，博士生，1985年生

湯愛平 男，博士，教授，1968年生 E-mail:taoshibo1985@163.com

U441.3;V211.3

A

10.13465/j.cnki.jvs.2017.04.023