基于混合萤火虫算法的桥梁颤振分析方法

陶仕博， 汤爱平， 胡庆杰, 刘克同

(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