静动载约束桁架结构拓扑优化的微分演化算法

胡长远，傅长荣，季光耀

（1.丽水学院工学院，浙江 丽水 323000；2.浙西南地质灾害研究所，浙江 丽水 323000）

结构拓扑优化是指以连续体开孔的数量及位置，或以离散体结构的节点布局、节点间的杆件连接关系为研究对象，通过改变结构的受力体系寻求最优的布局方案。拓扑寻优虽然可以最大程度优化结构，但求解空间和可行域在优化过程中都在不断更迭，这给结构重分析带来困难。较之尺寸优化和形状优化，拓扑优化具有更大的设计空间、经济效益和复杂程度，故逐步成为优化设计领域中的焦点和难点，受到更广泛的关注［1-3］。

目前拓扑优化方法多用于静荷载优化，而工程实际中结构常受到动荷载的作用［4］。固有频率是结构的重要特征，共振状态时固有频率的影响变得异常重要。为了保证结构的安全，须把固有频率控制在一定范围内［5］。拓扑优化设计静载约束有杆件应力、节点位移和欧拉屈曲等，动载约束有固有频率，静动荷载约束的叠加会使计算过程更为复杂、计算效率更低，而且难以收敛，无法得到较为理想的拓扑构形［6］。桁架结构拓扑优化目前主要采用基结构法，该方法以杆件截面面积为设计变量，通过逐步删除截面面积为足够小值的杆件来实现拓扑优化，求解简单方便，但存在奇异最优解问题［7］。如果用很小的断面积求解来代替零断面求解将受到各种约束的桁架拓扑优化问题的干扰，不仅会丢失可能存在的奇异最优解，而且会增加不必要的分析。以往的解决办法是二阶段尺寸和拓扑优化方法。该方法对尺寸和拓扑变量分别进行优化，求得最优拓扑后再进行尺寸优化，求解相对简单，拓扑构形容易实现。然而，由于尺寸和拓扑变量在数学上不是独立的，二阶段优化方法将真实搜索空间转化为近似域，可能会剔除包含全局最优的初始搜索域，大大增加了过早收敛的风险［8］。因此发展一些更加通用有效，能同时考虑尺寸和拓扑变量的全局优化算法来解决拓扑优化问题成为一种迫切的需要。

近年来，启发式算法逐渐被引入到工程优化领域，如遗传算法（GA）、微粒群算法（PSO）、和声搜索（HS）、教与学算法（TLBO）等［9］。它们不需要梯度信息，对目标函数和约束函数性态的要求较为宽松，且能兼顾静动态约束，因而较传统算法而言有更广泛的应用范围。计算中可以处理设计变量离散、可行域狭小的问题，能有效地得到高度复杂非线性问题的全局最优解，一些启发式优化方法已被应用于拓扑优化。Tejani 等［10］利用改进的共生生物搜索算法解决在多工况下的桁架拓扑优化问题。Dehghani 等［11］在帝国竞争算法中采用新的同化策略求解考虑尺寸、形状、拓扑变量的桁架优化问题。Dang等［12］将改进差分进化算法用于多种变量组合的桁架结构拓扑优化。Wang 等［13］用基于黄金分割系数的人工蜂鸟算法求解桁架结构的拓扑优化问题。

自Storn 和Price 于1996 年提出微分演化算法（differential evolution，DE）之后，该方法在解决复杂的优化问题上得到了很多关注［14］。微分演化算法是一种混合算法，它结合了遗传算法的更大种群概念和进化算法的自适应变异以及采用了贪婪选择策略。这些特征使微分演化算法相比进化算法和遗传算法具有更好的鲁棒性和更快的收敛速度。该算法已被广泛应用于工程优化、可靠度分析、光谱分析及无线电系统感应等领域［15-18］。

在结构工程中，微分演化算法应用于桁架拓扑优化方面的研究分析还很少。笔者将DE 算法应用于考虑尺寸和拓扑变量耦合的具有多工况静动荷载响应约束的桁架结构拓扑优化，通过对典型算例的求解并与相关文献的结果比较来评价笔者所提方法的有效性。

1 微分演化算法（DE）

一个包括t个参数的优化问题可以用一个t维的向量来描述，该向量可以表示为：xi=（xi1，xi2，…，xit）T∈S，i=1，2，3，…，NP。其中S∈Rn为优化问题的搜索空间；DE算法利用NP作为向量每一代的个体数。类似于遗传算法，DE算法通过变异、交叉和选择过程实现种群的更新进化。具体过程如下［15-16］。

1.1 变异过程

变异是为了保证种群的多样性，同时用合适的参数变化来指导已有的目标向量在合适的时间内达到一个更好的结果，从而保证了搜索的鲁棒性。

1.2 交叉过程

与GA 算法相似，DE 算法中的个体经过变异后也进行交叉操作。对于群体中第G+1代经过变异过程后的向量个体按照式（2）进行交叉，将产生新的个体：

式中：j=1，2，…，t；rand（j）∈［0，1］，是t个0～1 相互独立的随机数中的第j个；randn（i）是随机从集合{1，2，…，t} 中取得个体向量维度的序号；CR为交叉因子，将决定个体之间交叉的概率。

1.3 选择过程

DE 算法采用和GA 算法不同的贪婪准则：通过比较由变异和交叉产生的子代个体和父代个体，选择适应值好的变量，即如果父代个体适应值更优将继续保留在种群中，否则保留子代个体。选择过程如式（3）所示：

2 桁架结构拓扑优化

2.1 桁架结构拓扑优化的数学模型

结合拓扑优化的设计思想，构建以结构构件尺寸和拓扑变量为设计变量，以结构重量极小化为目标函数，受多工况荷载作用，考虑固有频率、杆件应力、节点位移、欧拉屈曲以及动态失稳等约束的离散型结构拓扑优化的数学模型，表达式为：

式中：X是设计向量，A1，A2，…，Am为杆件截面积，B1，B2，…，Bm为拓扑变量；Ai、Ei、Li、ρi、σi、σcri和Bi分别是第i类设计变量的截面面积、弹性模量、杆件长度、密度、杆件应力、欧拉屈曲应力和拓扑值。设Alim为杆件的临界面积，当Ai小于Alim时，Bi等于0，表示杆件删除；当Ai大于等于Alim时，Bi等于1，表示杆件保留。bj和δj分别表示节点j的质量和位移；ki为欧拉屈曲系数；fr是第r阶固有频率。上标 “comp” 表示抗压， “max” 和 “min” 分别表示上下限值。

数学模型中的约束函数还包括结构动态稳定性检查g6和有效性检查g7，前者是为了避免不稳定结构的出现，后者则是为了保证荷载作用点和支承节点等基本节点不被删除。由于结构的拓扑构形是由DE算法迭代时随机操作产生的，难免会出现一些不稳定的拓扑结构，因此有必要对新生成的结构进行稳定性检查。结构动态稳定性检查分两步进行：1）检查结构自由度是否违反格鲁布里尔准则；2）检查结构的总刚度矩阵是否为正定矩阵。

采用惩罚函数法处理上述约束，如果没有违反约束，惩罚函数F（X）为零；否则，按照以下标准执行处罚：

其中：

式（12）中：q为违反约束的个数；pi和分别表示违反的约束值及其限值；β1和β2为系数，本文均取1.5。

2.2 基于DE算法的桁架结构拓扑优化程序

为了提高DE算法的拓扑优化效率，采用以下改进策略和技巧：1）对于离散型结构拓扑优化，取截面下限Amin=-Amax，将搜索空间转为约两倍的设计变量限值，这样会丰富拓扑构形，有利于算法快速寻找到全局最优拓扑方案。2）对DE 初始随机生成以及经迭代交叉变异后生成的拓扑构形进行动态稳定性和有效性检查。若不满足要求，则更新个体，不再进入有限元计算，减少计算量。3）计算过程中，一旦杆件被删除，接续删除其连接，重新构建有限元模型，消除应力矩阵的奇异性，避免不必要的分析。

结合DE算法及前述拓扑优化数学模型，可以给出图1所示的离散型结构拓扑优化设计流程图。由算法原理及优化设计流程图可知，DE算法原理简单，能方便地应用到结构的拓扑优化设计中。

图1 基于DE算法的离散型结构拓扑优化设计流程图

3 算例分析

为验证笔者方法对桁架结构拓扑优化设计的有效性，以下对3 个典型桁架结构进行分析，通过改变杆件尺寸以及拓扑构形获得最小重量，并与其他算法的结果进行比较。对所有问题，欧拉屈曲系数（ki，i=1，2，…，m）和节点质量（bj，j=1，2，...，n）分别取为4.0和5 kg。DE算法的参数选择如下：F1=0.75，F=0.6，CR=0.85，20 杆和24 杆平面桁架算例种群规模NP 均取为100，72 杆空间桁架算例则取150。

3.1 24杆平面桁架

图2 所示为24 杆平面桁架的基结构，材料的弹性模量E=6.9×1010Pa，密度ρ=2 740 kg/m3，节点3上有双向集中质量500 kg。频率约束为f1≥30 Hz，应力约束为±172.43 MPa，节点5、6 竖向位移约束为δmax=10 mm。设计变量为连续变量，其下限为-40 cm2，上限为40 cm2，临界截面面积为1 cm2。两工况下节点荷载见表1。

表1 24杆平面杆桁架荷载工况

图2 24杆平面桁架的基结构

24杆平面桁架拓扑优化的目标函数收敛曲线见图3，优化后桁架最优拓扑布局见图4。为了验证笔者方法的有效性，与文献［8，12，19-22］相比较，对比结果见表2。

表2 24杆平面桁架拓扑优化结果比较

图3 24杆平面桁架拓扑优化的目标函数收敛曲线

图4 24杆平面桁架的最优拓扑形状

由图3可知：考虑静动荷载多种约束时，24杆平面桁架拓扑优化在迭代约150次后已经收敛，收敛速度快，不易陷入局部最优解。在表2 中，第1阶固有频率和节点6的竖向位移均已达到限值，优化后的桁架总重量为118.895 5 kg，优于文献［8，12，19-22］的结果。根据表2中的优化结果绘制其最优拓扑布局，如图4 所示，最优拓扑布局与文献［8，12，20，22］的结果相一致，杆件1～6，10～11，14，17～21，23以及节点1、4被删除，均为稳定性结构。

3.2 20杆平面桁架

本例的基结构由9 个节点和20 根杆件组成，如图5所示，其中节点1和9为支座，节点4上作用有荷载F1和F2。材料弹性模量E=6.9×1010Pa，密度ρ=2 740 kg/m3，两工况下节点荷载见表3。频率约束为f1≥60 Hz，f2≥100 Hz，应力约束为±172.43 MPa，位移约束以节点4 的竖向位移等于10 mm 和60 mm两种情形进行计算。设计变量为连续变量，其下限为-100 cm2，上限为100 cm2，临界截面面积为1 cm2。

表3 20杆平面杆桁架荷载工况

图5 20杆平面桁架的基结构

考虑不同位移约束的20杆平面桁架拓扑优化目标函数曲线见图6，优化的最终拓扑形状见图7。笔者计算结果与相关文献的比较见表4和表5。

表4 δ4y=10 mm时20杆平面桁架拓扑优化结果比较

表5 δ4y=60 mm时20杆平面桁架拓扑优化结果比较

图6 20杆平面桁架拓扑优化的目标函数收敛曲线

图7 20杆平面桁架的最优拓扑形状

由图6 可知：节点4 的竖向位移约束为10 mm和60 mm时，计算分别在迭代约40次和120次后完全收敛。表4 中，节点4 的竖向位移已达到限值。由表4和表5可知，两种情况下笔者优化后的桁架总质量分别为315.913 0 kg和150.637 8 kg，略好于文献［12］的结果，优于文献［4，19-22］的结果。对比表4和表5还可知，位移约束严苛时结构的总重量会增加。图7所示的最优拓扑形状中节点3、6、7被删除，布局上具有对称性，仅与文献［4］不同。

3.3 72杆空间桁架

如图8 所示的72 杆空间桁架基结构，由20 个节点和72 根杆件组成，依据结构对称特性，杆件被分为16 组。材料弹性模量E=6.895×1010Pa，密度ρ=2 767.99 kg/m3，节点1～4 上有三向集中质量2 270 kg。频率约束为f1≥4 Hz，f3≥6 Hz，应力约束为±172.375 MPa，节点1～4 的X和Y向位移限值为6.35 mm。设计变量为连续变量，其下限为-30 cm2，上限为30 cm2，临界截面面积为1 cm2。两工况下节点荷载见表6。

表6 72杆空间杆桁架荷载工况

图8 72杆空间桁架的基结构

静动荷载约束下72杆空间桁架拓扑优化目标函数收敛曲线见图9，与相关文献的比较列于表7，最优拓扑形式见图10。

表7 72杆空间桁架拓扑优化结果比较

图9 72杆空间桁架拓扑优化的目标函数收敛曲线

图10 72杆空间桁架的最优拓扑形状

图9 结果显示迭代约180 次时已找到最优拓扑布局方案，评价函数趋于收敛。由表7 可知：笔者DE 算法得出的桁架总质量为434.894 9 kg，优于其他所有结果，最优解受固有频率约束控制。图10最优拓扑形状保留了2层外围系杆，4层水平系杆及全部竖杆和腹杆，竖杆截面积向下逐层增大，与其他文献结果一致。

4 结论

DE 算法是一种启发式搜索的混合进化算法，和其他进化算法相比，DE算法不易陷入局部最优解，具有原理简单、易于实现以及鲁棒性强等优点。以往用二阶段优化方法求解桁架拓扑优化，虽然降低了求解难度与计算量，但割裂了解空间而难以保证得到全局最优解。笔者在DE 算法基础上，通过对每个杆件引入拓扑变量，便于实现同时考虑两种不同性质设计变量的拓扑优化问题。考虑到桁架结构在实际使用时可能会因动载作用产生噪声过大、共振等现象，需要在约束条件中引入固有频率约束。虽然考虑静动载约束耦合会造成求解难度的增加，但是更符合实际情况。此外，在优化过程中对有限元模型进行修正，重新构建刚度、质量和荷载矩阵，解决了矩阵奇异和不必要分析给计算带来的困难。数值分析结果表明：笔者方法很好地保持了所得最优解的稳定性、有效性以及较强的收敛性，可以成功地进行桁架结构的拓扑优化设计。