3.1　多工况结构拓扑优化设计

3.1.1　多工况结构拓扑优化设计的扩展最优性

目前在拓扑优化领域，多数的设计对象都是围绕单工况结构展开的，也就是在结构设计域中只考虑一种载荷工况。然而在实际工程结构的设计过程中，外部工作环境条件往往比较复杂，因此需要考虑适应多种变化的工况作用，应针对不同的工况载荷目标设计不同的支撑方式以满足多种功能需求。因此，对多工况下的结构拓扑优化问题进行研究更具有工程意义和实用价值。针对多工况下结构拓扑优化问题，需要权衡多条刚度分布路径，亦可称为多刚度拓扑优化问题。需要特别注意的是，多工况拓扑优化问题与多载荷拓扑优化问题是两个不同的问题，两者最大的区别在于不同载荷的施加是否一致：多工况拓扑优化问题中各类载荷不是同时加载，而多载荷拓扑优化问题中的所有载荷则是同时施加。这便意味着多工况拓扑优化问题的本质为一种多目标优化问题。

在经典的结构拓扑优化问题中，将结构在外载荷作用下刚度最大化作为优化目标函数是比较常见的算例。除此之外，通常还会根据设计者的经验对结构减重百分比（即体积比分数）进行预估，并将其作为优化模型中的约束条件（或者以最小体积作为优化目标，将预估的结构刚度值作为约束条件），见式（3-1）和式（3-2）。也就是说，传统拓扑优化方法的设计结果实际上有部分程度还需依赖于设计者的先验知识。由此可见，其与创新设计的自动化与智能化的理念存在矛盾。

问题1　以结构刚度最大化（通常以柔度C最小化表示）为优化目标，结构体积V不超过预估体积V₀为约束条件：

问题2　以结构体积V最小化为优化目标，结构柔度C不超过预估值C₀为约束条件：

目前已经有众多学者针对单目标优化问题进行了相关研究，所采用的方法为：对结构的刚度和体积同时优化，并将该类优化问题命名为拓扑优化中的“扩展最优性（extended optimality）”。2002年，扩展最优性的概念由Rozvany等[1]首次提出，该工作主要针对的是二维结构和单工况结构。该方法的基本思想是通过优化不同结构性能指标的乘积，来实现各指标的联合优化，如式（3-3）。随后，Rozvany[2]对比了传统拓扑优化设计与具有扩展最优性拓扑优化设计结果，指出引入扩展最优性能够得到重量较轻的设计方案。相较于乘积形式的目标函数f=CV，Strömberg[3]提出了一种更加通用的目标函数f=CVη，在该目标函数中，可以通过改变参数η的值，来获取更多不同的可选优化设计结果。上述研究均采用了基于SIMP材料插值模型的拓扑优化设计方法。实际上这种优化概念同样可以应用到其他类型的拓扑优化方法，如ESO方法等：Tanskanen[4]从理论角度对拓扑优化设计中的扩展最优性展开了深入的讨论，Edwards等[5]则对比了ESO方法和SIMP方法在解决考虑扩展最优性的结构拓扑优化设计问题上的优缺点。值得注意的是，以上全部研究主要是采用了乘积形式的目标函数，且主要是用于解决单工况结构拓扑优化设计问题。

问题3　以结构柔度C与结构体积V的乘积为目标函数，同时对两项结构指标进行优化，即所谓的“扩展最优性”：

通过上述论述可以看出，扩展最优性的提出主要是为了更好地应对工程结构轻量化设计需求，这种设计概念不需要设计者在优化开始时就根据经验给出结构刚度或体积分数的预估值，而是通过求解优化问题来自动得到更加丰富且更加轻质的设计方案。可见，考虑扩展最优性的多工况结构拓扑优化问题建模与求解方法的建立就具有十分重要的理论意义和工程价值。接下来，本节将对上述问题进行重点论述。

3.1.2　多目标拓扑优化模型

对于多工况结构拓扑优化问题来说，至关重要的一点就是通过构建合适的多目标优化模型，使其能够搜寻整个Pareto最优前端。理论上，Pareto最优的充分／必要条件是判断能否通过某种特定的多目标模型求解出整个Pareto最优解集的关键[6]。一方面，当多目标优化模型满足Pareto最优的必要条件时，通过该模型能够解出全部的Pareto最优解，但其中部分解可能不是Pareto最优解；另一方面，当多目标优化模型满足Pareto最优解的充分条件时，通过该模型解出的全部解均为Pareto最优解，但其无法求出整个Pareto解集。也就是说，当选择一种合适的多目标优化模型，使其满足Pareto最优解的充要条件，通过该模型就一定能求解出整个Pareto最优前端上的解。在多工况拓扑优化设计领域，因为优化模型的设计变量多，学者们通常采用直接加权法来处理多个载荷工况下的子目标函数，但无论如何选择子目标权重，都无法在优化问题Pareto前端非凸的情况下找到其全部Pareto最优解集[7]。这也就意味着，针对某些多工况拓扑优化问题，如果设计者改变子目标重要程度偏好，就可能会导致无法获得相应的Pareto最优解的问题。

指数加权准则（exponential weighting criterion，EWC）是一种在理论上被证明的、能够满足Pareto最优充要条件的多目标模型。[6]，[8]，[9]通过该模型，无论多目标问题的Pareto最优前端是否非凸，都可凭借连续地改变子目标权重来获取全部的Pareto最优解集。基于EWC的多目标优化模型可以表述为

其中，w_k为权重系数；m代表子目标的个数；子目标f_k(x)通过指数加权的方式转化为统一的单目标F(x)；指数常数q为实数，在这里取q=2，以保证模型能够找到多工况拓扑优化问题的整个Pareto最优前端。

需要注意的是，在采用模型（3-4）中所定义的目标函数时，q和f_k（x）通常作为指数，如果取值较大，极易造成计算数值过大甚至数值溢出（numerical overflow）。此外，当各子目标的数量级显著不同时，最优设计一定会更容易受到量级较大的目标影响，从而出现“载荷病态”的情况发生。因此，需要对基于EWC的多目标优化设计模型进行进一步的改进。本书提出一种基于归一化指数加权准则（normalized exponential weighting criterion，NEWC）的多目标拓扑优化设计模型。利用基于水平集的结构边界隐式描述方式，引入径向基函数插值中的扩展系数α作为设计变量，可以得到基于参数化水平集的优化模型：

其中，C_k(u，Φ)为子工况k作用下的结构应变能：

a_k(u，v，Φ)为子工况k的能量双线性形式：

l_k(v，Φ)为子工况k的载荷线性形式：

其中，D_pqrs为材料的弹性张量；m为所考虑的工况数量；N为有限元单元数量；V_max为体积约束值；u₀为Dirichlet边界上的位移；H(Φ)为Heaviside函数；δ(Φ)为Heaviside函数的偏导数，即Dirac函数。α_max和α_min分别为设计变量的上下限，主要用于增加优化迭代的稳定性。对于任意子工况k，u_k和v_k分别表示实位移场和运动学允许的位移空间U内的虚位移场，ε为应变场，p_k为结构体积力，τ_k为结构边界上的牵引力， 分别代表单工况k作用下的目标函数最小值与最大值，其值可通过分别求解单工况作用下的拓扑优化设计问题来预估得到。与目标函数式（3-4）相比，多目标拓扑优化设计模型式（3-5）将各个子目标值进行了归一化处理，以便有效消除由于各子工况载荷作用下的结构应变能数量级差距过大而引起的载荷病态情况。此外，归一化处理后的目标函数值被限定在一个较小的区间内，能够有效克服优化计算时的数值溢出问题。

在结构拓扑优化设计中考虑扩展最优性，其最直接的目的是充分发掘结构材料的承载潜力，以获得材料用量更少、更加轻质的结构构型。如上文所述，目前的所有研究中，学者们主要采用f=CV或f=CVη这一类乘积形式的单目标函数，而该目标函数产生可行解的多样性有限，并没能足够有效地扩大解空间。根据工程经验可知，结构的柔度和体积是一对相矛盾的指标，即提升其中一个的数值，必将会导致另外一个的数值降低。故考虑扩展最优性的拓扑优化问题本质应为一个多目标优化问题，不同柔度和体积的组合代表了解空间内所对应的Pareto最优解[10]。这里是将扩展最优性引入多工况结构拓扑优化设计问题中，并构建一种更具普遍性的多目标优化模型：

求α=[α₁，α₂，…，α_N]T

其中，Vmin和Vmax分别代表体积分数约束的最小值与最大值，在本书中设置为Vmin=0.05和Vmax=1。需要指出的是，由于Vmin取值为0代表结构域内没有任何材料，没有任何实际工程意义，故Vmin通常取值为趋近于0的一个较小值。权重向量 代表各子工况作用下的结构应变能及结构整体体积分数的权重，根据多目标优化模型的要求，所有权重系数之和设置为1。

如上所述，构建了一套考虑扩展最优性的多工况拓扑优化模型。在该多目标拓扑优化设计模型中，所有不同类型的子目标函数均实施了无量纲化操作，通过NEWC方法转化为单一优化目标，有利于多设计变量、大规模的工程优化设计问题求解。通过连续地改变权重参数，可以得到多套不同的轻量化设计方案，甚至多目标拓扑优化问题的整个Pareto前端[6]。

3.1.3　目标函数权重决策方法

在本节中，将介绍基于模糊多属性群决策（fuzzy mult-iattribute group descision making，FMAGDM）的决策方法计算多目标拓扑优化模型式（3-9）的权重系数。为了得到更可靠的方案，减少错误决策的风险，人们在进行决策活动时通常会兼顾和集中多个专家的意见，这就是所谓的群决策。由于客观事物的不确定性、不完全性以及人类思维的模糊性，导致大多数多属性群决策（mult-iattribute group decisionmaking，MAGDM）问题具有不确定性和随意性，故而形成了FMAGDM理论。随着Bellman和Zadeh[11]在1970年首次利用模糊数学知识给出模糊决策的基本模型，处理具有不确定性的参数、概念和事件，基于模糊集理论（fuzzy set theory）的FMAGDM获得了国内外学者的广泛关注[12]。在FMAGDM理论框架中，专家采用语言术语（linguis ticterms）来描述待决策对象的属性，然后将语言术语翻译成相应的模糊数并进行特定的模糊数计算、集结等操作，最终获得专家的一致性评价意见。对多目标拓扑优化问题而言，由于人的偏好以及设计方案与权重系数间的关系均具有模糊性，设计人员在确定子目标权重系数时，很难直接给出具体的数值。所以要求设计人员对各子工况下的结构柔度和结构体积指标给出诸如“比较重要”“不重要”等之类的评价语言，是一种相对切合实际的做法。FMAGDM恰好符合这种评价方式，也就是通过模糊的语言术语来描述对各个子目标重要程度的偏好。在对优化模型式（3-9）的子目标权重系数进行评估的过程中，子目标被视为待决策的“属性”，设计人员被视为“决策者”，并根据以往积累的工程实践先验知识对子目标给出语义形式的偏好。

采用图3-1所示的梯形模糊数（trapezoidal fuzzy number）[13]，[14]来处理权重系数计算过程中的模糊信息。假设有n个设计者D=(D₁，D₂，…，D_n)对m个优化子目标O=(O₁，O₂，…，O_m)进行评估。其中偏好矩阵E包含m×n个梯形模糊数 ={e_ij1，e_ij2，e_ij3，e_ij4}，且 代表设计者D_j对子目标O_i的偏好。特别地，多数FMAGDM方法并未考虑到决策者自身的知识水平对决策结果的影响，在所提出的方法中将另外采用语言术语描述决策者的等级。反映设计者重要程度及能力水平的等级矩阵记为ω={ω₁，ω₂，…，ω_m}，其中设计者D_j的等级可解释为梯形模糊数ω₁=(ω_j1，ω_j2，ω_j3，ω_j4)。这样能保证较重要的设计者偏好对最终决策结果起较大作用，反之亦然。

尽管FMAGDM理论的方法对处理含有模糊信息的语言术语有效，但这种基于专家知识的决策方法仍然带有一定的主观随意性。为了使多目标权重系数能够更加准确地反映实际工程需求，将根据客观历史数据或实验结果对子目标的取值加以限定，从而构建一种自适应权重调整机制，以确保避免因主观偏好的误差导致优化的结果不准确。综上，所提出的基于FMAGDM理论的权重计算方法基本步骤如下：

步骤1　设计人员对子目标采用语言术语进行评价，并将给出的评价矩阵E和专家等级矩阵ω译为对应的梯形模糊数。从表述最低等级的语言术语“AP”或“AL”到描述最高等级的语言术语“AG”或“AH”，每一项语言形式的偏好都有与其唯一对应的梯形模糊数[15]，如表3-1所示。

步骤2　采用梯形模糊数混合集结算子（trapezoidal fuzzy number hybrid aggregation operator，THA）评估偏好 ，从而得到综合考虑了每位设计人员意见及其等级的子目标偏好 ：

其中， 代表设计者D_j对子目标O_i的加权偏好意见； 为 从大到小排序后的第j个模糊数；式（3-10）中的权重系数 由下述公式计算[16]：

其中，定量模糊语义算子Q由一种分段函数形式给出：

式中的参数a和b由表3-2中列出的模糊语言定量法则确定。

步骤3　计算各子目标综合偏好 的大小。将 所包含的梯形模糊数转化为参数化形式 =（σ，x₀，y₀，β），其大小可依据式（3-14）计算：

为方便计算，式（3-14）中加权函数取为g(x)=x。 的大小反映了其隶属度的高低。

步骤4　计算所有子目标的权重系数：

步骤5　利用式（3-16）给出的自适应权重调整机制修正权重系数：

其中，OV代表子目标函数值；UB和LB分别为子目标取值的上限和下限；l代表需要修正权重系数的子目标数量。较大的修正因子η会出现迭代不稳定的问题，较小的η则会延长迭代收敛的时间，故在本章中取η=2。

所提出的权重计算方法具有非常鲜明的特点：①采用语言术语描述是一种更加合理的方式，基于模糊集理论对语言术语进行处理提高了所提出方法的实用性；②基于群决策的方式有效降低了权重计算错误的风险，特殊的集结算子能够过滤掉设计人员的极端偏好；③对于采用先决策后优化方式的多目标模型在应用中更加便捷、简单，且使用范围广。

3.1.4　灵敏度分析与优化算法

与传统拓扑优化模型对比可以发现，多目标拓扑优化设计模型式（3-9）中没有任何结构响应约束，这表明无需在优化求解过程中更新拉格朗日乘子。根据形状导数概念分别对目标函数中的结构应变能子目标F_C和体积子目标F_V进行灵敏度分析，最后采用本书第2章所介绍的OC法对优化模型进行求解。

直接目标函数F_C对时间变量t的偏导可得

其中，偏导数∂C_k/∂t前的系数简记为W_k。

结构应变能的形状导数可推导为

其中：

将式（2-22）定义的参数化水平集速度场代入式（3-18）：

一方面，将式（3-20）代入式（3-17）可改写目标函数的形状导数：

另一方面，采用链式求导法则推导目标函数F_C关于时间变量t的偏导：

对比式（3-21）和式（3-22），易知子目标F_C关于设计变量α的敏度为

类似地，可推导体积子目标F_V关于设计变量α的敏度：

3.1.5　数值实例

1．二维悬臂结构

平面悬臂梁的设计空间如图3-2所示，长宽比为2∶1，左端固支，结构域受2个载荷的作用，即F₁=F₂=500kN。材料的弹性模量为201GPa，泊松比为0.3。有限元单元与水平集网格离散均为120×60个。该优化问题的目标为同时优化各子工况下的结构柔度及体积分数。

如前文所述，在加权多目标优化问题中，子目标权重系数的大小对优化设计结果有重要影响。为揭示权重系数大小与设计结果的内在联系，本例对不同权重系数下的多目标拓扑优化问题进行了求解和对比分析。首先进行11组试验（试验一），将结构体积目标的权重和各子工况下的柔度目标权重和作为试验对象，采用试验设计（design of experiments，DOE）中的均分法，以0.1为试验间隔求解不同权重系数下的多目标优化问题。为方便讨论，本例假设每个柔度子目标的权重系数相等。理论上，体积目标的权重与柔度目标权重和的取值范围为[0，1]。然而当体积目标的权重等于0时，结构会删除所有材料导致柔度无限大；当柔度目标权重和等于0时，结构不会删除任何材料导致优化的结果为实体。这两种情况在优化中均无实际应用价值，故将权重的变化范围限定为[0.01，0.99]。

测试结果表明，当体积目标的权重系数大于一定值时，例如wV>0.3时，优化迭代过程会随着材料的急剧减少而不断振荡，从而直接导致优化难以收敛（>400步迭代），优化结果中某些结构刚度路径甚至会因为优化迭代过程振荡发生断裂。因此仅给出如表3-3和图3-3所示的6组测试结果。为方便讨论，在此定义一种材料性能指标来反映最优结构拓扑中每一单位数量的材料减少所导致的结构刚度丧失，即刚度降低率（decreasing ratio，DR）：

其中，ΔC和ΔVF分别代表相邻两次测试中结构柔度的差值和体积分数的差值；T代表不同试验的编号。DR的理论意义与经济学中的“边际效益”类似，在此可称之为“边际刚度”(marginal stiffness)：较大的DR意味着每减小一单位的结构体积会导致结构刚度的剧烈下降，此时以牺牲结构刚度为代价来降低结构重量是“不经济”的，反之亦然。

从表3-3可以看到，当结构的最优体积分数在0.2以下时，DR值急剧增加，表明结构的最优刚度会随着体积分数的减小而剧烈下降。因此在二维情况下，推荐将考虑扩展最优性的多工况拓扑优化问题的体积分数下限设为0.2，以保证优化结果材料使用的“经济性”。需要特别提出的是，三维空间的结构拓扑优化问题不同于二维问题。测试结果表明，由于三维结构具有更强的承载能力，体积分数下限可设置为0.1。从图3-3可以看出，不同权重系数的组合会得到不同拓扑形式的优化结果，但本方法始终能保证结构边界的清晰和光滑。

此外，数值算例表明当体积目标的权重系数在[0.01，0.99]间取值时，最优设计的拓扑构型对权重系数十分敏感。表3-4和图3-4分别给出了试验二的数值结果和结构拓扑形式。优化结果显示，体积目标的权重系数的微小改变就会导致优化设计的体积分数发生较大差异，尤其是当权重系数wV<0.05时最优设计的体积分数变化更加明显。通过控制体积分数上限可一定程度避免权重系数的微小改变带来体积分数的剧烈波动。而且为实现结构的轻量化设计，工程上通常将优化设计方案的体积分数设置在0.5附近。在本实例中体积分数0.5将被设置为体积目标值的上限，以保证优化结果的减重效果。

将各工况下结构柔度的加权和作为一种特殊的结构柔度目标，结构的体积分数作为体积目标，可将考虑扩展最优性的拓扑优化问题简化视为一种两目标优化问题。经简化后的优化问题的部分Pareto最优前端如图3-5所示。

2．二维桥梁结构

桥梁的拓扑优化属于大型连续稀疏结构的设计问题，在实际工程中该问题会涉及十分复杂的工况载荷。桥梁结构的设计空间如图3-6所示，桥面以下左右端固支，桥面上综合作用5种载荷工况，且F₁=F₂=F₃=F₄=F₅=50kN。材料的弹性模量为180GPa，泊松比为0.3。整个结构设计域采用240×60个四边形单元进行离散，其中桥面（受载区域）为非设计区域，由240×2个被动单元（passive element）构成，与被动单元关联的CSRBF扩展系数不参与优化迭代。该优化问题的目标为同时最小化各子工况下的结构柔度及体积分数。利用结构对称性，仅7200个单元参与计算。

首先采用基于FMAGDM的权重计算方法确定多目标拓扑优化问题的权重。5个设计者对6个目标的评价，以及5位设计者的水平及经验的综合等级如表3-5所示。为说明方法的稳定性，特别令设计者1给出与其他设计者完全相反的极端评价意见。通过计算，可以快速得到规范化的子目标权重系数w=[0.2356，0.1514，0.0982，0.0952，0.0911，0.3285]。计算结果表明，个别极端的评价意见对最终的权重系数确定没有明显影响。

接着对基于NEWC的多目标优化模型进行求解。通过对各工况单独作用时拓扑优化问题的快速求解，可以预估子目标的最大和最小值： =25.8242和 =0.2584； =18.8834和 =0.1881； =13.5942和 =0.1356； =8.7989和 =0.0879； =5.4642和 =0.0549。经过优化迭代，本问题的最优目标函数值为3.1235，其对应的各子工况下的最优结构柔度值为C=[1.3292，0.9653，0.5900，0.3686，0.3038]，结构体积分数为0.2479，迭代步数为151，对应的结构最优拓扑形式如图3-7（f）。结构拓扑、水平集函数及CSRBF扩展系数分布的演化过程分别如图3-7、图3-8和图3-9所示。可以看到结构边界在优化初期经历剧烈的拓扑和形状变化，最终收敛到光滑的最优拓扑形状，说明该方法很好地保留了水平集方法的优点。扩展系数作为设计变量在优化过程中动态地调整，优化后的扩展系数沿结构边界均匀化分布，表明优化过程已趋于稳定，此时的结构边界Γ_d应为最优边界 。从优化后的结构拓扑来看，任一工况载荷作用处均有相应的刚度传递路径延伸出来，并最终与主传递路径交汇，保证了刚度的有效传递，说明本节所涉及的方法可以有效地克服多工况拓扑优化问题中载荷病态的缺陷。

优化问题的迭代收敛曲线如图3-10所示。各子目标的收敛曲线如图3-11所示。从图3-10中看到，整个优化迭代过程高效且稳定，目标函数值在第20步迭代时已逼近最优值，此时结构的拓扑基本形成，体积约束已接近设定的体积分数，后续131步迭代过程主要是针对结构局部拓扑的微调以及进行结构形状优化。从图3-11可以看到，在初始的20步迭代内，各子目标值均发生了一定程度的变化，部分目标值甚至出现局部振荡，主要是因为迭代初期结构的材料去除速度较快，导致拓扑变化较大，各子目标通过不断变化以便获得较为稳定的折中值。在这一迭代阶段完成后，各子目标值均开始稳定收敛，最终得到相应的Pareto最优解。说明该方法对二维拓扑优化设计问题的有效性。

考虑到拓扑优化中存在的“边际刚度”随体积分数的减少而快速递减的现象，最优设计的体积分数会直接影响材料的使用效率以及结构的性能。在考虑扩展最优性的拓扑优化问题中，结构的体积分数在优化前是未知的，当设计者对体积目标给出较为过大或过小的权重系数时，结构的体积分数也会相应地达到某种过高或过低的水平。因此引入自适应权重调整机制将二维拓扑优化问题的最优体积分数维持在区间[0.2，0.5]内，以保证材料使用的经济性。根据式（3-16）可得到关于体积目标权重系数的自适应调整算子，用于进行体积目标和柔度目标权重系数的修正。

为说明自适应权重调整机制的作用，本例给出另外5位设计者对6个子目标的偏好矩阵，其中所有设计者均对体积目标较为偏重，见表3-6。通过计算可获得规范化的权重系数w=[0.1503，0.1036，0.0819，0.0702，0.0819，0.5121]。

设计者第二种偏好下的桥梁结构最优拓扑形式如图3-12所示。结构的体积分数和体积目标权重的收敛曲线如图3-13所示。从图中可知，前15次迭代中结构的体积分数总是高于所设定的下限值（图中红色虚线），此时的结构体积分数稳定下降，体积目标权重值保持不变，材料使用的“经济性”良好。从第16次迭代开始，体积分数开始低于下限值0.2，材料的使用变为“不经济”的状态，自适应权重调整机制开始起作用，体积目标权重开始逐渐下降，直到第54次迭代结构的体积分数被拉升到下限值以上，此时体积目标权重维持不变。各子工况下的最优柔度值为2.0338、1.4036、0.6774、0.5646和0.4039，体积分数为0.2006，迭代步数为232。特别指出，虽然结构初始化的体积分数可能高于上限值0.5，但随着OC算法的迭代，体积分数值会迅速下降至上限值以下，故可以假设在这一阶段上限值是不起作用的。本实例说明所论述的优化方法能通过调整权重值来保证材料的使用效率。