对复杂系统整体进行优化难度非常大。将复杂的优化问题分解为多个较简单的子问题,并根据子问题之间的耦合关系进行协调是一种有效的解决方案。Kroo等^[1]提出的协同优化方法结构简单,易于实现学科自治,是近年来受到广泛关注的一种分解和协调策略^[2-4]。在结构设计优化方面,不少学者已引入协同优化方法处理复杂结构优化问题。苏瑞意等^[5]对某全承载大客车车身骨架进行尺寸优化时,采用协同优化方法将不同的性能包括质量、刚度、强度、振动模态和翻滚等分解至多个学科中进行优化。Balling等^[6]将协同优化方法应用于轮毂框架结构截面尺寸优化中,相比于传统的整体优化,采用协同优化方法计算效率提高,且优化问题越复杂,提高越明显。

离散结构是很多工程结构的主要承载构件,如客车骨架结构、桥梁、航空航天结构等^[7-8]。近年来,很多学者对离散结构进行优化研究^[9-10]。在优化迭代过程中,需要反复进行结构有限元分析。随着结构设计规模变大,有限元分析耗费更多的时间,占用更多的资源^[11]。针对大规模结构分析问题,子结构(超级单元)方法是一种有效的方法^[12]。大模型被分解为多个规模较小、不重叠的子结构,每个子结构提供与其他子结构连接的外部节点自由度有关的缩聚矩阵,所有缩聚矩阵组装后得到系统矩阵。求解系统平衡方程即可得到结构内部边界条件。然后,将结构内部边界条件施加到相应的子结构,从而可对子结构独立进行有限元分析,恢复子结构数据(节点位移、单元应力等)。子结构方法求解静力学问题精度高,能够降低计算规模,提高计算效率。

本文提出一种基于子结构方法的分解策略,并采用基于子结构分解的协同优化框架求解强度约束下的离散结构优化问题。

1 标准协同优化算法

协同优化的核心思想是将复杂系统优化问题分解为多个子学科优化问题。双层结构是协同优化的典型结构形式,包含系统层和学科层,如图 1所示。假设原优化问题如式(1)所示：

$\begin{align} & \min f(X,Y) \\ & s.t.g\left( X,Y \right)\le 0, \\ & Y=h\left( X \right), \\ & {{X}_{\min }}\le X\le {{X}_{\max }}. \\ \end{align}$

(1)

其中： f(X,Y)为系统目标,g(X,Y)为约束函数,Y为状态变量,可分为非耦合状态变量Y_u和耦合状态变量Y_c。X为设计变量,包含共享设计变量X_sh和局部设计变量X_l。X_max和X_min分别为设计变量的上、 下限。

第i个学科的优化模型如式(2)所示：

$\begin{align} & \min {{J}_{i}}=\left\| {{X}_{\text{sh}i}}-{{Z}_{\text{sh}i}} \right\|_{2}^{2}+\left\| {{Y}_{ji}}-{{Z}_{ji}} \right\|_{2}^{2}+ \\ & \left\| {{Y}_{ik}}-{{Z}_{ik}} \right\|_{2}^{2} \\ & \text{s}\text{.t}\text{.}{{g}_{i}}\left( {{X}_{i}} \right)\le 0, \\ & {{Y}_{ik}}={{h}_{i}}\left( {{X}_{i}} \right), \\ & {{X}_{i,\min }}\le {{X}_{i}}\le {{X}_{i,\max }}. \\ \end{align}$

(2)

其中： X_i为学科i的设计变量,包括与本学科相关的共享变量X_sh、 本学科的局部设计变量X_l以及其他学科输入的耦合状态变量Y_ji。Y_ik为学科i输出的状态变量(学科i至学科k)。Z_i为系统层传递至学科i的相关变量目标值。J_i为学科i的相容性函数,g_i为学科自身的约束。

系统层的优化模型如式(3)所示：

$\begin{align} & \min f\left( Z \right) \\ & \text{s}\text{.t}\text{.}J_{i}^{*}\le \varepsilon ,i=1,2,\cdots ,{{N}_{\text{d}}}, \\ & {{Z}_{\min }}\le Z\le {{Z}_{\max }}. \\ \end{align}$

(3)

其中：J_i^*为第i个学科的相容性函数,ε为给定的允许容差。Z为系统层设计变量,包含共享设计变量X_sh和耦合状态变量Y_c。系统层的目标是优化系统目标,约束为保证各个学科优化产生的共同变量以及耦合状态变量的一致性。N_d为学科数目。

通过系统层与学科层的反复迭代,可以找到一个满足协调性的最优结果。协同优化方法求解复杂问题时,能够实现学科自治,各个学科可以并行优化,选用合适的优化算法。

2 离散结构优化问题描述

本文研究的对象为离散结构(如桁架、梁架等)在应力约束下的结构质量最小化问题,优化模型如式(4)所示：

$\begin{align} & \min W \\ & \text{s}\text{.t}\text{.}\left| {{\sigma }_{i}} \right|\le {{\sigma }_{u}},i=1,2,\cdots ,{{N}_{\text{L}}}, \\ & {{x}_{\min }}\le {{x}_{j}}\le {{x}_{\max }},j=1,2,\cdots ,N. \\ \end{align}$

(4)

其中： W为结构整体质量,σ_i为杆件i的应力,σ_u为允许的应力上限,N_L为杆件数目。x_j为设计变量,x_max、 x_min分别为设计变量的上、 下限,N为设计变量数目。

静力学平衡方程为

$KU=F.$

(5)

其中： K为整体刚度矩阵,U为节点位移,F为作用于结构上的外力向量。

随着结构复杂度增加,设计变量数目急剧增加,优化算法寻找最优解更加困难。因此,为了更加有效地求解复杂离散结构优化问题,本文基于有限元分析中的子结构方法,提出一种子结构分解策略,并结合该策略的特点,在标准协同优化算法的基础上,提出一种改进的协同优化框架。该框架继承了标准协同优化算法的学科自治和分布式并行的优点,但系统层和学科层的构建方法不同,主要包括以下特点： 系统层不进行优化,仅进行结构分析,以避免大量的整体结构分析; 耦合变量在学科层作为定值,以避免引入过多状态变量而导致学科层优化困难; 学科层优化目标为结构的物理性能,而不是相容性函数。

3 基于子结构分解的协同优化框架 3.1 基于子结构方法的分解策略

如图 2所示,沿着图中虚线将结构分解为2个子结构Ω₁和Ω₂,虚线构成结构内部边界。每个子结构与其他子结构连接的节点称为外部节点,其余为内部节点。在此例中,位于结构内部边界的单元属于子结构Ω₁。

用I表示子结构的内部节点集合,B表示子结构的外部节点集合,则整体刚度矩阵和第i个子结构的刚度矩阵分别为：

$K=\left[ \begin{matrix} K_{II}^{\left( 1 \right)} & K_{IB}^{\left( 1 \right)} & O \\ K_{BI}^{\left( 1 \right)} & K_{BB}^{\left( 1 \right)}+K_{BB}^{\left( 2 \right)} & K_{IB}^{\left( 2 \right)} \\ O & K_{BI}^{\left( 2 \right)} & K_{II}^{\left( 2 \right)} \\ \end{matrix} \right],$

(6)

${{K}^{\left( i \right)}}=\left[ \begin{matrix} K_{II}^{\left( i \right)} & K_{IB}^{\left( i \right)} \\ K_{BI}^{\left( i \right)} & K_{BB}^{\left( i \right)} \\ \end{matrix} \right].$

(7)

其中： K_II⁽ⁱ⁾∈R^N₁×N₁,K_IB⁽ⁱ⁾∈R^N₁×N₂,K_BI⁽ⁱ⁾∈R^N₂×N₁,以及K_BB⁽ⁱ⁾∈R^N₂×N₂,N₁为集合I中节点的自由度,N₂为集合B中节点的自由度。

子结构i的静力学平衡方程为

$\left[ {\matrix{ {K_{II}^{\left( i \right)}} & {K_{IB}^{\left( i \right)}} \cr {K_{BI}^{\left( i \right)}} & {K_{BB}^{\left( i \right)}} \cr } } \right]\left[ \matrix{ U_I^{\left( i \right)} \hfill \cr U_B^{\left( i \right)} \hfill \cr} \right] = \left[ \matrix{ F_I^{\left( i \right)} \hfill \cr F_B^{\left( i \right)} \hfill \cr} \right].$

(8)

其中： U_I⁽ⁱ⁾和F_I⁽ⁱ⁾分别为集合I中节点的位移和受到的外力; U_B⁽ⁱ⁾和F_B⁽ⁱ⁾分别为集合B中节点的位移和受到的外力,需满足与其他子结构的协调关系。使用子结构方法对整体结构进行有限元分析时,首先对内部节点位移U_I⁽ⁱ⁾进行缩聚,可以得到子结构i的缩聚方程为

${{\tilde{K}}^{\left( i \right)}}U_{B}^{\left( i \right)}={{\tilde{F}}^{\left( i \right)}}.$

(9)

其中： ${{\tilde{K}}^{\left( i \right)}}$和${{\tilde{F}}^{\left( i \right)}}$分别为缩聚刚度矩阵和缩聚载荷矩阵：

${{\tilde{K}}^{\left( i \right)}}=K_{BB}^{\left( i \right)}-K_{BI}^{\left( i \right)}{{\left[ K_{II}^{\left( i \right)} \right]}^{-1}}K_{IB}^{\left( i \right)},$

(10)

${{\tilde{F}}^{\left( i \right)}}=F_{B}^{\left( i \right)}-K_{BI}^{\left( i \right)}{{\left[ K_{II}^{\left( i \right)} \right]}^{-1}}F_{I}^{\left( i \right)}.$

(11)

然后,将缩聚矩阵组装为系统矩阵并求解,可获得结构内部边界上节点的位移,将与该子结构相关的边界条件U_B⁽ⁱ⁾代入式(8),即可对子结构i进行单独求解,恢复子结构i的内部数据。

进行离散结构协同优化时,每个子结构构成一个学科,其相关属性(如杆件截面积)在迭代过程中不断发生变化产生新的子结构,每个新的子结构性能分析相当于子结构方法中子结构的数据恢复过程,结构内部边界条件在每个学科优化开始前由系统层输入。学科层优化过程中,不涉及系统求解,因此不需要给出每个新结构如式(10)和(11)所示的缩聚矩阵。本文施加的结构内部边界条件为位移或者力,选择结构内部边界条件类型的准则是： 当子结构存在刚体位移时,选择位移作为边界条件,否则选择力。因此,U_B、 F_B是学科之间的耦合状态变量。

结构优化问题采用的子结构分解策略具有以下特点： 1) 耦合状态变量对系统目标函数,即结构质量,没有影响; 2) 系统目标和约束函数均可根据结构分解方式完全分解成只与学科设计变量相关的学科目标和约束函数。

3.2 系统层和学科层构造

本文提出的协同优化分解协调策略如图 3所示,其中设计变量X是与结构的节点或杆件相关联的属性。X⁽¹⁾和X⁽²⁾分别为两个子结构对应的设计变量。系统层起协调作用,仅进行整体结构有限元分析,不进行优化。为使各学科之间解耦,引入状态变量U_B^(i)*和F_B^(i)*,由系统层传递至各学科。状态变量在各学科中作为定值,从而避免受到耦合带宽(即结构内部边界上节点的自由度数)限制以及额外增加设计变量。初始的状态变量值从一个给定的初始结构有限元分析结果中提取。之后,每一次系统循环时,各学科将其优化获得的设计变量最优解输入系统层,经整体结构分析可获得该最优解对应的结构内部边界上节点真实的位移和力,更新状态变量值。在整个迭代过程中,每一次经过系统层协调之后,确保了相邻子结构的结构内部边界条件一致,即隐含了对状态变量的相容性约束。

原优化问题的目标和约束随着结构的分解,相应地分配到各个学科。学科以最小化对应子结构的质量W_i为目标,约束为子结构i中所有杆件的应力均小于许用应力。

3.3 算法流程

整个优化过程是一个反复分解-协调的迭代过程。用算法形式描述如下：

1) 设k=0,给定初值X₀;

2) 系统层进行整体结构分析,求得结构内部边界条件U_B^*、 F_B^*;

3) 将U_B^*/F_B^*作为定值传递至学科i,从X_k中提取X_k⁽ⁱ⁾作为学科i的初值;

4) 优化学科i,获得第k次系统迭代时学科i的最优解X_k^(i)opt;

5) 是否达到最大迭代次数？是则转7),否则转6);

6) k=k+1,X_k=[X_k-1^(1)opt X_k-1^(2)opt],将X_k作为定值传入系统层,转2);

7) 结束。

4 算 例

为了检测本文提出的离散结构协同优化算法的性能,选用经典的结构建立相应的优化模型^[13]。本文选用2维的20-杆桁架结构^[14]和3维的72-杆空间桁架结构^[15]作为测试算例。

4.1 20-杆桁架结构

算例1为图 4所示的20-杆平面桁架结构,其opt材料参数如下： Young's模量为30 000 ksi (即206.85 GPa),密度为0.1 lb/in³(即2 769.73 kg/m³),屈服极限为±25 ksi (即172.38 MPa)。节点1和5承受大小为 100 lb的力(即444.8 N),方向如图 4所示。设计变量是各杆件的截面积,可在 [0.1 15] in²(即[0.645 96.774] cm²)范围内取值,以最小化整体结构质量作为目标,所有杆件的应力不大于屈服极限。因此,该优化问题包含20个设计变量和20个约束条件。

将该结构沿着节点3和4所在直线(图 4中虚线)分为子结构1(虚线左侧)和子结构2(虚线右侧),位于结构内部边界上的杆件属于子结构1。相应地,原结构静力学求解时需要求解包含16个线性方程的方程组; 分解后,需并行求解两个包含8个线性方程的方程组。两个子结构对应的学科分别包含10个设计变量和10个约束。子结构1输入的结构内部边界条件是节点3和4的力,子结构2输入的结构内部边界条件是节点3和4的位移。

在Isight软件中搭建协同优化模型,使用Dell Precision T7910(16核,64 GB)工作站进行计算。两个学科均采用序列二次规划法求解,选择4个不同的初值,其中2个初值为不可行解,2个初值为可行解。从不同的初值出发,均收敛到非常相近的最优解,如表 1所示,两个子结构的质量迭代过程如图 5所示,均快速收敛。1 in²=6.45 cm²。1 lb=453.59 g。

作为对比,对该算例进行整体优化。同样从表 1所示的4个初值出发,采用序列二次规划法求解,协同优化和整体优化的结果对比列于表 2。可见,整体优化收敛到不同的最优解,稳定性不如协同优化算法。另外,协同优化方法需要的整体结构分析次数明显减少,且找到更好的解。有限元分析花费的总时间见表 2。由于协同优化时各学科并行优化,因此协同优化的总时间为学科层中花费的最多时间与系统层花费的时间之和。虽然学科层子结构分析的次数较多,但除了初值3外,协同优化的总时间与整体优化相差不大; 而整体优化从初值3开始优化时,出现了明显的早熟收敛现象。

4.2 72-杆空间桁架结构

算例2 为72-杆空间桁架结构,如图 6所示,由4个相同的模块组成。材料特性为： Young's模量为 10 000 ksi (即68.95 GPa),密度为 0.1 lb/in³(即 2 769.73 kg/m³),所有杆件的拉压应力均应小于25 ksi(即172.38 MPa)。承受的两种载荷工况列于表 3 (1 kips=4.448 kN)。设计变量为杆件截面积,允许的最大截面积为4 in²(即25.81 cm²),最小截面积为0.1 in²(即0.645 cm²)。所有杆件分成16组： 1) A₁—A₄,2) A₅—A₁₂,3) A₁₃—A₁₆,4) A₁₇—A₁₈,5) A₁₉—A₂₂,6) A₂₃—A₃₀,7) A₃₁—A₃₄,8) A₃₅—A₃₆,9) A₃₇—A₄₀,10) A₄₁—A₄₈,11) A₄₉—A₅₂,12) A₅₃—A₅₄,13) A₅₅—A₅₈,14) A₅₉—A₆₆,15) A₆₇—A₇₀,16) A₇₁—A₇₂。每组中杆件的截面积要求相同。该结构优化问题共有16个设计变量和72个约束。

每层结构形成一个子结构,从下至上分别为子结构1、 子结构2、 子结构3和子结构4。位于结构内部边界的杆件属于相邻子结构中的下层子结构。原结构静力学平衡方程包含48个线性方程; 分解后为4个包含12个线性方程的方程组,且4个方程组并行求解。每个学科有4个设计变量和18个约束。子结构输入的结构内部边界条件列于表 4。

计算环境同4.1节算例1。4个学科均采用修正可行方向法求解。分别从一个不可行初值和一个可行初值出发,最终收敛到相近的最优解附近,如表 5所示,子结构的质量收敛曲线如图 7所示,均以非常快的速度收敛到最优解。

同样对该结构进行整体优化,采用修正可行方向法求解。表 6列出协同优化与整体优化的求解结果。其中协同优化的总时间计算方法同算例1,由于从初值1出发整体优化未找到最优解,因此未列出其时间; 而从初值2出发,协同优化的寻优效果优于整体优化,需要的整体结构分析次数非常少。求解效率两者相差不大。

5 结论

本文针对复杂结构优化时存在设计变量多、寻优困难等问题,提出基于子结构分解策略的协同优化算法。该方法继承了传统协同优化方法学科层并行优化、学科自治等优点,同时降低了结构分析的复杂度,且不会因为耦合带宽增加,额外引入设计变量而导致优化问题更加复杂。

根据优化对象的结构形式将其分解为多个子结构,目标函数、约束函数以及设计变量均相应地分配至各个学科。学科之间的状态耦合变量为结构内部边界上的位移和力边界条件。

结合基于子结构分解策略,每次系统迭代时,系统层进行一次结构分析来协调各个学科之间耦合状态变量的差异,显著减少了整体结构分析次数。学科层优化时,只需要对该学科对应的子结构进行有限元分析,降低了结构分析的复杂度。

最后,2维的20-杆桁架结构和3维的72-杆桁架结构算例结果表明： 相比于采用相同求解算法的整体优化,本文提出的基于子结构分解策略的协同优化算法更加稳定,能够找到更好的解。