采用量子粒子群优化算法的风荷载作用下伞形膜结构形态优化设计

孙芳锦， 苏基浩， 李智达， 张大明

(1. 桂林理工大学 广西岩土力学与工程重点实验室， 广西 桂林 541004；2. 桂林理工大学 广西建筑新能源与节能重点实验室， 广西 桂林 541004；3. 桂林理工大学 土木与建筑工程学院， 广西 桂林 541004； 4. 桂林理工大学 信息科学与工程学院， 广西 桂林 541006)

膜结构是一种柔性结构体系，由于膜材不能抗压、抗弯，所以膜结构分析比刚性结构多了一个初始形态确定的过程[1].结构的“形”和“态”相互影响且一一对应，两者共同决定了膜结构的力学性能.因此，形态优化分析对于膜结构设计具有重要的指导意义.膜结构优化理论最初是以最小曲面[2]为优化目标提出的.2002年，钱基宏等[3]提出了膜结构“最优形态”的概念.随后，卫东等[4]把结构整体刚度最大，将受荷后各单元应力差值最小和结构整体投资值最小作为膜结构的优化目标，并解决了膜结构多目标优化的问题.伞冰冰等[5]以应变能、最大支座反力、应力波动系数为膜结构优化目标，展开了分析研究.

目前，膜结构优化的方法主要有传统优化算法和遗传算法[6].传统优化算法不利于解决具有多种约束或目标函数难以直接用数学式表达的优化问题，而遗传算法具有收敛速度慢、局部搜索能力差和控制变量较多等缺点.粒子群优化算法[7-8]和量子粒子群优化算法是近年来发展迅速的智能仿生优化算法，其不仅具有概念简单、执行简易、收敛迅速快等优点，且对约束条件和目标函数的要求较低，因此，被广泛应用于工程优化领域中.梁美荣[9]将粒子群优化算法应用到膜结构的找形分析中，并得到相关结论.

本文针对风荷载作用下的伞形膜结构展开形态优化分析研究，分别采用粒子群优化算法和量子粒子群优化算法[10]自编程序并对风荷载作用下的伞形膜结构进行优化分析，得到伞形膜结构最优形态的相关参数，为张拉膜结构设计提供合理的设计依据.

1 多目标优化方法

1.1 优化原理

优化的一般思路首先是提出优化目标，并以某种数量的形式表达，它是结构优化运算的一个评价标准；其次是确定约束条件和优化变量.薄膜结构的优化目标有支座反力、结构刚度和薄膜应力分布均匀程度等，由于膜结构具有多个目标函数，可以通过分配系数法将其转化为单目标函数.

1.2 目标函数的建立

1.2.1 刚度最大[11]将荷载作用下结构产生的应变能E作为刚度的表达方式，应变能大小反映了结构刚度的强弱程度，结构产生的应变能越小，刚度越大.目标函数f1表达式为

minf1=E

.

(1)

1.2.2 受荷后应力分布最均匀[11]由于膜面的面积大且厚度薄，为了避免膜面发生应力集中从而引起撕裂破坏，采用各有限单元的最大主应力均方根与其平均值的比值来衡量应力分布均匀程度，并称为应力波动系数D.应力波动系数越小，代表应力分布越均匀.目标函数f2表达式为

(2)

1.2.3 受荷后最大支座反力最小[11]为了减少膜结构给地基和下部结构的负担，选取各支座中最大支座反力Fmax为优化目标.目标函数f3表达式为

(3)

式(3)中:Fi，x，Fi，y，Fi，z分别为第i个支座x，y，z向的支座反力.

1.3 优化变量的选取

决定膜结构初始形态的参数有索预拉力、膜初始预应力、结构控制点坐标、结构跨度和高度等.由于膜预张力的可调整变化范围较小，可通过对不同膜预张力的结构进行优化分析来研究膜预张力的影响，伞形膜结构伞高度为H，平面对角线跨度为L，取矢跨比δ(δ=H/L)和脊索拉力T为优化变量.

1) 膜结构形状主要由支承点的相对位置及预应力分布决定.对于伞形结构，支承点相对位置主要是指矢跨比，具体限值根据建筑功能要求给出[12]，即

δmin<δ<δmax.

(4)

式(4)中：δmin，δmax分别为最小、最大矢跨比.

2) 设定可知索截面积为As，维持曲面形状的最小拉力值为Tmin，脊索材料强度为σmax，安全系数为ε(ε<1)，则最大拉力值Tmax表达式为

Tmax=ε×σmax×As.

(5)

优化变量的取值范围通常根据工程的常用范围而定，并要求满足建筑功能[12]，T的表达式为

Tmin

(6)

2 基于量子粒子群优化算法的优化设计方法

2.1 粒子群优化算法优化策略

粒子群优化(PSO)算法主要是通过经典速度更新公式和位置更新公式实现迭代运算，即

(7)

(8)

式(7)，(8)中：V为粒子速度；r1，r2为(0，1)之间的随机数；X为粒子位置；k为算法当前迭代次数；w为粒子保持的运动惯性；c1，c2分别为局部学习因子和全局学习因子.

w一般采用从0.9到0.4线性递减的策略[13]，在此情况下，PSO算法具有较好的优化效果，计算结果已通过测试验证.w的表达式为

(9)

式(9)中：wmax，wmin分别为最大、最小惯性权重，wmax=0.9，wmin=0.4；kmax为算法的最大迭代次数.

为了确保PSO算法前期全局搜索能力较强，c1采用不断减小的策略，具体表达式为

(10)

式(10)中：cs为最大学习因子，cs=2.5；ce为最小学习因子，ce=0.5.

后期为了确保粒子具有较强的局部搜索能力，c2采用不断增加的策略[14]，具体表达式为

(11)

粒子群优化算法具体有以下5个步骤.

1) 初始化.初始化粒子种群数目，随机设定初始位置和初始速度.

3) 更新信息.通过比较每个粒子当前位置和历史最优位置更新速度和位置信息.

4) 更新gbest.通过比较上个计算迭代周期的gbest和当前所有的pbest，实现gbest更新.

5) 停止条件.如果计算结果满足设置条件，则运算结束；如不满足设置条件，则返回步骤2)继续运算，直到满足条件，最后输出搜索结果[15].

2.2 量子粒子群优化算法优化策略

量子粒子群优化(QPSO)算法是在粒子群优化算法的基础上，取消了粒子位置变化方向属性并引入平均的粒子历史最好位置mbest[16].QPSO算法将不再考虑粒子的历史运动记录，其算法步骤更新为

1) 计算平均的粒子历史最好位置mbest，即

(12)

式(12)中:S为粒子群的数目；plocal，i为当前迭代中第i个个体最优粒子.

2) 粒子位置Pi更新为

Pi=φ×plocal，i+(1-φ)pglobal.

(13)

式(13)中:Pi为第i个粒子位置的更新；pglobal为到目前为止全局最优粒子；φ为随机参数.

为了增加粒子运动的随机性，避免陷入局部收敛，将随机参数φ变为2个随机参数，粒子位置更新公式被修改为

(14)

式(14)中:φ1，φ2为(0，1)之间的随机数.

最后，得到粒子位置更新公式为

xi=Pi′±α|mbest-xi|ln(1/u).

(15)

式(15)中：xi为第i个粒子的位置；u为(0，1)中的均匀分布数值，取+和-的概率各为0.5，当u>0.5时，Pi后运算取+，反之取-；α可以根据迭代次数不断更新，从而使粒子不断趋于最优解位置.α是QPSO算法中的唯一参数，其取值一般小于1，α的表达式为

(16)

量子粒子群优化算法有以下4个步骤.

1) 参数初始化，设置粒子尺寸、粒子种群数目、运算终止条件和最大迭代次数.

2) 通过比较粒子当前位置适应度和到目前为止最优位置的适应度，决定是否更新最优位置信息和局部最佳权重.

3) 更新迭代mbest，根据mbest和粒子位置更新公式(15)，从而更新粒子信息.

4) 若迭代次数达到上限或达到设定终止条件，则停止搜索输出优化结果；若未达到，则回到步骤2)继续搜索[16].

2.3 算法适应度函数的构造

适应度函数f是算法进行迭代计算的衡量标准，其表达式为

(17)

β1+β2+β3+β4=1.

(18)

式(17)，(18)中:βi为优化目标的权重分配系数；fi为优化目标，各个优化目标之间存在不同量级的问题，通过优化目标乘以系数ζi，则各个fiζi成为统一数量级关系，ζi表达式为

(19)

3 风荷载下伞形膜结构的形态优化分析

3.1 优化目标分析和权重确定

以伞形膜结构为算例[12]，采用刚性直线性边界结构，脊索沿对角线布置，对角线跨度为20 m，结构顶端开洞，洞口直径为1 m，膜材厚度为0.001 m，薄膜经、纬向刚度均为1 100 kN·m-1，结构剪切刚度为100 kN·m-1，泊松比为0.2，脊索截面积为3.14×10-4m2，结构受竖向均布荷载为-0.48 kN·m-2，风荷载为20 m·s-2.伞形膜结构示意图，如图1所示.

(a) 俯视图 (b) 主视图 (c) 风向角示意图 (d) 结构模型图

与雪荷载不同，风荷载不仅竖直向下作用于结构，还会垂直作用于膜表面.为了消除膜表面积的影响，设单位面积膜应变能E1、单位长度索应变能E2、应力波动系数D和结构最大支座反力Fmax为膜结构优化目标[4]，设矢跨比和脊索拉力为优化变量.E1，E2的表达式为

E1=E膜/A，

(20)

E2=E索/Lt.

(21)

式(20)，(21)中:E1为单位面积膜应变能；E膜为膜总应变能；A为膜表面积总和；E2为单位长度索应变能；E索为索总应变能；Lt为脊索总长度.

为了得到优化目标的权重分配系数，首先进行简单的参数分析，脊索拉力T的取值要在具体设计荷载作用下保证结构内部具有维持曲面形状的应力值区间，T取10，16，22，28 kN；矢跨比δ采用伞形膜结构的常用取值范围[5]，δ取1/4.20，1/2.50，1/1.70，1/1.34.各优化目标与矢跨比和脊索拉力的关系，如图2所示.

(a) 单位面积膜应变能 (b) 应力波动系数

由图2可知：随着矢跨比的增大，应力波动系数和最大支座反力逐渐增大，但单位长度索应变能和单位面积膜应变能逐渐减小；单位长度索应变能、单位面积膜应变能在脊索拉力较小时达到最小值，同时，最大支座反力随着应变能的减小而减小.因此，最大支座反力的权重分配系数可以适当减小；相对于单位长度索应变能，单位面积膜应变能的变化量级更大，变化范围更广，优化的目的是使总应变能尽可能小，所以，单位面积膜应变能的权重分配系数应比单位长度索应变能的大.

综上所述，无论矢跨比和脊索拉力趋于更大或更小，都会导致有些优化目标向不好的方向发展，从而使结构因为某些性能被忽略而提前产生破坏.

引进层次分析法[17]给每个优化目标分配权重系数，4个目标函数则可集合成一个目标函数.优化目标层次结构模型，如表1所示.表1可反映优化目标之间的重要程度关系.表1中：数值1表示横坐标和纵坐标重要性同等重要；数值3表示横坐标优化目标的重要性比纵坐标的强；数值1/3表示横坐标优化目标的重要性比纵坐标的弱.经过反复多次试算，得到优化目标应力波动系数、最大支座反力、单位面积膜应变能、单位长度索应变能对应的权重分配系数β1，β2，β3，β4分别为0.35，0.15，0.35，0.15.

表1 优化目标层次结构模型

3.2 优化过程分析

在伞形优化模型确定的基础上，分别以PSO和QPSO算法编制程序并对其进行优化分析.采用表1的优化目标权重分配系数，优化过程如图3，4所示，优化结果如表2所示.

表2 优化结果

(a) 单位面积膜应变能 (b) 应力波动系数

每个粒子对应一个适应度和一组目标函数值，每组目标函数值是4个优化目标加权相加后的结果.优化目标由于设置了权重，没法保证每个优化目标值在迭代过程中都能比上一代更小，只要优化目标值的主要变化趋势随着适应度的减小而趋于减小，即是合理的.由图3，4可知:QPSO算法的适应度、单位面积膜应变能和应力波动系数的变化趋势保持一致；两种算法的适应度一直朝着更小方向变化，表示搜索到的目标函数越来越小，目标函数越小，则结构越优化，最后适应度趋于平稳，即两种算法已找到最优解位置.

由图4可知：QPSO算法在迭代更新到63代时得到最小适应度5.325 8，而PSO算法在迭代更新到72代时得最小适应度5.327 1.与PSO算法相比，QPSO算法能以更少的迭代次数和更小的适应度找到最优解，证明了QPSO算法能以更快的速度找到最优解且优化效果更好.

图4 QPSO算法和PSO算法的适应度

3.3 结果分析

为了探究伞形膜结构最优形态相关参数之间的规律，采用量子粒子群优化算法对伞形膜结构展开大量的优化分析.风荷载取值为风荷载规范中5～11级区间内的风速；跨度和矢跨比取值范围为伞形膜结构常用参数[5]；膜预张力和脊索拉力大小由设计荷载确定，基本原则是在设计荷载作用下保证结构内部具有维持曲面形状的预应力值，同时随着跨度变化，结构需要的膜预应力和索拉力会发生变化[11-12].16，20，25 m跨度结构的优化结果，如表3～5所示.表3～5中：St为膜预应力；v为风速.

表3 16 m跨度结构的优化结果

表4 20 m跨度结构的优化结果

表5 25 m跨度结构的优化结果

由表3～5可知：随着膜预应力的增大，脊索拉力最优值也会随之增大，矢跨比最优值出现微弱减少；随着风速的增大，矢跨比和脊索拉力最优值几乎不受影响，即最优结构参数与风荷载大小无关；随着结构跨度的增大，脊索拉力最优值增大，矢跨比最优值变化细微且稳定在一个定值附近.

为探究最优结构形态的跨度、脊索拉力和膜预应力之间的关系，定义比值系数γ=T/(St×0.5×L)[5].最优形态的比值系数γ和矢跨比δ，如图5，6所示.

图5 最优形态的比值系数 图6 最优形态的矢跨比

通过以上分析可知：在风荷载作用下，最优伞形膜结构的脊索拉力、膜预应力和跨度之间的比值系数γ在0.77附近波动，最优结构矢跨比出现在0.33～0.41，大部分最优值分布在0.35附近.将程序计算结果与系数测算结果进行对比验证，结果基本一致.

4 结论

根据膜结构的特点建立伞形膜结构多目标优化模型，并采用粒子群优化算法和量子粒子群优化算法对伞形膜结构展开优化设计研究，得到以下2点结论.

1) 相对于粒子群优化算法，量子粒子群优化算法的收敛速度快且具有较强的全局搜索能力，在伞形膜结构形态优化分析中能够得到更精准的全局最优解.

2) 在风荷载作用下，最优伞形膜结构的脊索拉力、膜预应力和跨度之间的比值系数γ在0.77左右，最优结构矢跨比为0.33～0.41，最优值主要集中在0.35附近，结构最优形态与风荷载大小无关.