互逆规划的拓宽和深化及其在结构拓扑优化的应用1)

2020-12-23隋允康彭细荣

力学学报 2020年6期

关键词：不合理约束定理

铁军隋允康彭细荣

∗(天津财经大学理工学院，天津 300222)

†(北京工业大学工程数值模拟中心，北京 100022)

∗∗(湖南城市学院土木工程学院，湖南益阳 413000)

引言

自1988 年Bendsoe 与Kikuchi[1]提出均匀化方法以来，连续体结构拓扑优化得到了快速发展，同时涌现出许多方法，包括变密度方法[2-3]、独立连续映射(independent,continuous and mapping，ICM)方法[4-6]、ESO(evolutionary structural optimization)方法[7]、水平集方法[8-9]、拓扑导数法[10-11]、相场法[12]、MMC(moving morphable components)方法[13-14]、CBS(closed B-splines)方法[15]、材料场序列展开法[16]等.对各类方法的综述可参阅Rozvany[17]及Sigmund[18]等的文献.

各种连续体结构拓扑优化方法基本上集中在建模途径和求解效率上，而对于优化模型是否合理以及模型之间的相互关系则缺乏注意.研究优化模型的合理性始于隋允康等将MCVC(minimum compliances with a volume constraint)模型同MVDC(minimum volume with displacement constraints)模型进行比较[19-21]，说明MCVC 模型约束条件值的选取和多工况时目标函数的不合理性.

需要指出的是(1)“合理与否”是各个层次的“结构优化模型”中的共同问题，只是在结构拓扑优化中显得更加突出.(2)“优化模型的合理与否”不同于“结论的正确与否”.结论的正确与否证明即可，而优化模型的合理与否需要花费较多的篇幅阐述.

为了节省读者查阅文献[19-21]时间，这里给出“结构优化设计模型合理与否”的诠释.

在结构优化设计的模型中，常常存在如下特点的问题：(1)本来可以建立单目标函数最优的问题却不适当处理成多目标最优的问题；(2)需要先验、无根据地预先给出一些约束值，从而导致最优点不恰当地依赖于约束值而变成为难以从最优解集中确定出有限个解；(3)多层次的优化问题不能按同一个标准建立模型，以至于有些不同层次选取的目标函数与约束函数是不同的，甚至出现关注没有工程意义的力学性能的情况；(4)本来应该追求满足工程规范或力学性能的经济指标最小，却不必要的追求工程规范或力学性能的最佳而很难预先设定经济指标约束限.

如果上述模型经过对于目标函数与约束条件的对调，其中的任何一个问题都不再出现，那么，就可以把原来的问题称为一个“不合理的模型”，相应地，调整后的问题成为“合理的模型”.

在连续体结构拓扑优化的背景下，作者团队进一步研究模型合理化与否，同运筹学中数学规划论的对偶规划[22-25]比对，从而提出了互逆规划的概念[26].

互逆规划的两方具有目标函数与约束条件相互交换位置的对应关系，它的意义在于：揭示了数学规划领域里，除了存在对偶理论，还存在与之表面相似却不同的互逆规划理论，也就是说，对于数学规划补充了新的理论内容；利用该理论审视结构拓扑优化的模型，揭示了不合理的模型是因为相关的研究者没有采用互逆规划的s 方(单目标函数)建模，而是采用了m 方(多目标函数)建模.研究表明，互逆规划在结构拓扑优化研究中具有明确的力学背景和广泛的应用价值.

基于文献[26]提出的互逆规划，本文发展了相关定义及理论，内容如下：

(1) 把s-m 型互逆规划拓宽出s-s 型互逆规划和m-m 型互逆规划，从而得到3 种类型的互逆规划.

(2) 对于3 种互逆规划分别推导出各自两方面的KKT 条件关系，并且揭示了两方关系的本质，提出了3 个定理.

(3) 3 种互逆规划的求解策略和具体解法.

(4)互逆规划对结构拓扑优化的应用举例.

1 将互逆规划由1 种推广为3 种

文献[26]提出了s-m 型互逆规划，由s 方和m 方两个问题数学规划问题组成：(1)s 方——取单个目标函数和多个约束条件；(2)m 方——取s 方的多个约束条件作为该方的目标函数，取s 方的目标函数作为该方的约束条件.

为了完善互逆规划的理论，本文把单目标和单约束构成的一对互逆规划定义为s-s 型互逆规划，把多目标和多约束构成的一对互逆规划定义为m-m 型互逆规划.

从表达方便出发，将各种互逆规划的两方按表示的顺序分别称为正向问题(forward problem,FP)、反向问题(reverse problem,RP).若交换表示的顺序，也就将一对互逆规划的FP 问题和RP 问题换位，两种表达的本质是一样的.

(1) s-s 型互逆规划

这是FP 问题和RP 问题皆为单目标函数的互逆规划.也就是说，其FP 问题和RP 问题都只有单个约束条件.

(2) m-s 型互逆规划

这是FP 问题和RP 问题分别为多目标函数和单目标函数的互逆规划.也就是说，其FP 问题和RP 问题分别具有个单约束条件和多个约束条件.

FP 问题和RP 问题孰在前、孰在后，只是顺序的不同，并不影响问题的本质，因而m-s 型互逆规划同文献[26]提出的s-m 型互逆规划实际是同一个事物.

(3) m-m 型互逆规划

这是FP 问题和RP 问题皆为多目标函数的互逆规划.也就是说，其FP 问题和RP 问题也都具有多个约束条件.

以上提出的s-s 型互逆规划和m-m 型互逆规划，是对文献[26]提出的互逆规划的推广.推广之后的互逆规划共分上述3 类.

2 互逆规划基于KKT 条件的3 个定理

一对互逆规划，由于双方的约束条件不同，二者的可行域是不同的，一般说来，不能保证两个可行域有交集，更不可能期望这对互逆规划的解相同.本文的研究得到两个结论：其一、发掘出一对互逆规划两个KKT 条件内在关连的3 个定理；其二、利用一对互逆规划任一方的最优解调整另一方的约束条件，可以求出拟相同的最优解.

下面分别叙述s-s 型、s-m 型和m-m 型互逆规划问题的上述两个结论.限于篇幅，本文只给出定理2 的证明，略掉了相对容易的定理1 与3 的证明.

2.1 s-s 型互逆规划的两方KKT 条件及其关系

定理1对于s-s 型互逆规划的双方，其一、两个KKT 条件在形式上相像，称为拟同构；其二、部分调整参数后的互逆规划双方在数值上拟同解，即取一方的最优解或KKT 点的目标函数值，作为另一方面的约束限值，则该最优解也满足调整了上述参数后的另一方的KKT 条件.

2.2 s-m 型互逆规划的两方KKT 条件及其关系

对于s-m 型互逆规划，也有同定理1 类似却不同的定理，为了方便予以叙述和论证，首先定义多目标规划FP 问题式(2a)和单目标RP 问题式(2b)的两方Lagrange 函数如下

需要说明一下：一般的多目标权系数都有归一化条件λ1+λ2+···+λL=1，λi∈(0,1)，为了推导的方便，忽略这个不影响求解的条件，因为只要将每个权系数除以权系数之和，即可满足归一化条件.

若函数Hi(x) (i=1,2,···,L)和G(x)是可微函数，且x∗和依次是式(2a)的帕累托最优解和式(2b)的最优解，则可以给出s-m 型互逆规划的如下定理.

定理2对于m-s 型或s-m 型互逆规划的双方,其一、两个KKT 条件拟同构；其二、两种情况下互逆规划双方拟同解：情况1、先求解m 方，取其帕累托最优解的多目标函数值，作为s 方的多约束条件限值，然后求解调整上述参数后的s 方，则s 方与m 方拟同解；情况2、先求解s 方，取其最优解的单目标函数值，作为m 方的单约束条件限值，且取s 方的最优的Lagrange 乘子向量值作为m 方的多目标的权系数向量值(可以不满足归一化条件)，然后求解调整上述参数后的m 方，则m 方与s 方拟同解.

注：本定理中拟同构和拟同解的概念与定理1 中的叙述相同.

证明：

记G(x∗)=G∗，同样式(7)也可以扩展为

下面证明定理2 的后半部分，具体可以分成两种情况：

情况1：从式(2a)的解到式(2b)的解.

与式(8)同样记式(2a)的最优目标值为Hi(x∗)=(i=1,2,···,L)，用此来调整式(2b)的约束值，即取=(i=1,2,···,L).为了进一步的证明，把式(8)与式(9)分别写为

结论是：由式(2a)的解调整式(2b)的约束值之后，式(2a)的Pareto 最优解是式(2b)的KKT 点，或称为两方拟同解.

情况2：从式(2b)的解到式(2a)的解.

在式(12)的求解中，其中第2，3 式不需要求解，因它们只表示Hi(x∗)=(i=1,2,···,L).欲从式(12)的其余式子中求解，只需将(12)和(13)进行比较，可看出x∗与ˆx∗是对应的形式参数，由此解得=1，=x∗.

结论是：由式(2b)的最优解对式(2a)的两个参数予以调整，其一是用(2b)的最优化目标代替(2a)的约束值，其二是用(2b)的最优的Lagrange 乘子代替(2a)的目标函数权系数的若干倍，此时的式(2b)单目标最优解是式(2a)的带权重系数的多目标规划的KKT 点，称为拟同解.

两种情况可以概括为定理2 的后半部分，根据任一方的最优解调整第二方的约束限并取为约束式，以及在第二方为多目标时，用第一方的Lagrange 乘子代替多目标的权系数的若干倍，此时两个互逆规划拟同解.

2.3 m-m 型互逆规划的两方KKT 条件及其关系

定理3对于m-m 型互逆规划的双方，表达Pareto 最优解须满足的KKT 条件，其一、双方的KKT 条件是拟同构；其二、部分调整参数后的互逆规划与另一方是拟同解，即取任一方最优点时多目标函数值作为另一方的多约束限值，且取第一方最优的Lagrange 乘子向量值作为第二方的多目标函数权系数向量值(可以不满足归一化条件)，则按上述调整参数后的第二方与另一方拟同解.

注：本定理中拟同构和拟同解的概念与定理1 中的叙述相同.

3 互逆规划的应用和求解

3.1 互逆规划的应用

(1)鉴别结构优化模型的合理与否

文献[26]应用互逆规划揭示了结构优化中存在有不合理模型，例如：结构拓扑优化结构体积约束下多工况柔顺度最小问题，属于不合理的多目标模型，依照本文发展的互逆规划理论予以剖析，是m-s 型问题的FP(正向问题)方.

(2)使不合理的结构优化模型合理化

尽管上述m-s 型问题的FP 方属于不合理的多目标模型，但是基于互逆规划，取其RP(反向问题)方可以得到单目标模型，是合理的模型.于是，本来在FP 问题必须面对的多目标函数和体积约束限难以选取的两个致命困难，在RP 问题中就不复存在了.

单工况下，若目标函数取柔顺度最小，鉴于其结构体积约束值很难预先确定，也属于不合理问题.按照本文拓宽出的s-s 类模型互逆规划去理解和处理，取其对应的另一方，就可以得到合理的模型了.

(3)利用互逆规划的一方求解另一方

本文完善的互逆规划理论不仅对于建立模型的合理化有裨益，而且采用3 个定理，对于已经建立起的模型，可以提供有效的解法.这就是后续即将叙述的互逆规划的求解策略和具体解法.

(4)借鉴互逆规划理论求解多目标规划

多目标规划问题的求解，在数学规划中发展得比较缓慢，尤其寻求Pareto 最优解集时，前提是要选定恰当的权系数，这是较难操作的.本文提出了预选权系数的新途径，即按涉及m-m 型问题的定理3 予以处理，详细做法将在3.3 小节叙述.

总之，文献[26]揭示了存在互逆规划的两方，而本文则扩充了互逆规划的类型，并且揭示了互逆规划双方的内在关连.前者有助于识别出不合理的状态，并且将其转化为合理的问题，后者则依靠互逆规划的双方关连性补充了求解的方法.

3.2 互逆规划的求解策略和具体做法

互逆规划的求解策略有如下3 点.

(1)优化模型合理与否

首先从工程实际或物理问题的自身背景角度，审视该问题是否合理？如果是合理的，就进入求解阶段；如果属于不合理的问题，建议从互逆规划出发，考虑其另一方，那可能就是一个合理的模型.

(2)合理的优化模型

对于合理的规划问题，如果是单目标函数，就按以往的有效算法；如果是多目标函数，3.3 小节将予以阐述.

(3)不合理的优化模型

对于必须求解的不合理规划问题，首先按互逆规划理论建立其对应的合理问题，然后按合理的一方求解，从而确定不合理规划一方难以确定的参数，最后就可以求解不合理规划问题了，3.4 小节将予以阐述.

求解互逆规划的具体做法按如下3 种情况进行.

(1)对于s-s 型互逆规划

存在2 种情况：其一是，待求解的一方是合理的，直接求解即可；待求解的一方是不合理的，则先求解其互逆对应的合理一方，然后借鉴定理1调整不合理一方的参数后，再予以求解.

(2)对于m-s 型或s-m 型互逆规划

存在4 种情况：其一是，待求解的s 方是合理的，直接求解这个单目标、多约束的问题即可；其二是，待求解的s 方是不合理的、而m 方是合理的规划，按3.3 小节介绍的凝聚函数方法处理目标函数，先求解m 方，用其解调整s 方的参数，最后予以求解；其三是，待求解的m 方是合理的规划，按凝聚函数方法处理目标函数，然后求解；其四是，待求解的m 方是不合理的、而s 方是合理的规划，先求解s 方，用其解调整m 方的参数，按帕累托模型求解m 方.

(3)对于m-m 型互逆规划

存在2 种情况：其一是，待求解的一方是合理的，直接求解即可，其中多目标按凝聚函数方法处理；其二是，待求解的一方是不合理的，则先按凝聚函数处理其互逆对应的合理一方的目标函数，接着求解，然后借鉴定理3，用合理方的结果调整不合理方的参数后，最后求解其帕累托模型.

3.3 多目标问题的两种解法

定理2 与定理3 都涉及到多目标规划的求解问题，本文建议如下两种解法：

(1)基于minimax 策略的凝聚函数解法

对于一个多目标优化问题，可以引用隋允康提出的统一的凝聚函数方法以及铁军、隋允康提出的抛物线凝聚函数方法[27-28]，通过适当简化以便于数值计算的方法将多目标化为一个凝聚的单目标函数.需要指出，统一的凝聚函数方法统一了模优化(norm optimization)和K-S 函数，即

(a)统一的凝聚函数

(b)模优化(norm optimization)

(d)抛物线凝聚函数

利用凝聚函数，多目标优化问题可以转化为单目标优化问题.凝聚函数Hagg(Hi(x),p)代替式(3a)的多个目标函数的原因在于：凝聚函数有一个珍贵的性质，对于适当选择的拉伸因子P，随着P的增大，凝聚函数趋于=maxHi(x) (i=1,2,···,L).由此式得

式(18)左边表达了L个多目标函数，右边表达了1 个单目标函数.式(18)的本质是：它表述了多目标优化的minimax 解法，而凝聚函数则是实施该解法的途径之一.

将式(18)和凝聚函数代入式(3a)中，得到了用单目标函数求解多目标优化式(3a)的如下问题

特别需要指出的是，在采用式(14)∼式(17)把多目标处理为一个凝聚函数时，为了满足被凝聚的每个函数都大于0 的条件，可以对于每一可能会小于0 的函数加上适合自己的充分大正数.

往下的具体求解可以采用研究者青睐的各种行之有效的约束规划算法，本文倾向于应用DSQP(对偶序列二次规划)算法[29]或DP-EM 方法[30-31].

(2)依据互逆规划定理求解帕累托模型

一般地，多目标规划问题并不存在唯一的最优解，故以往解决多目标规划问题通常多依赖于寻求解帕累托模型的过程.该解法的前提是，设取L个目标函数的权系数，其本质也是将多个目标函数化为单个目标函数.众所周知，权系数绝非是理性的预判，而是感性偏好的预设，最多也只能说是以往经验的积累.

顺便阐述帕累托解法与minimax 解法的关系.由=maxHi(x)(i=1,2,···,L)得

进而得

式(21)与式(18)比较，表明minimax 解与帕累托解的关系，类似于minimax 解与一般多目标解的关系.二者的区别在于：前者落脚于右式的单目标的选取方式，后者落脚于左式的多目标的权系数的选取上.

利用帕累托解法求解一般多目标问题解的前提是恰当地选取权系数，为了提出解决这一关键问题的新方法，前面在3.2 小节介绍具体做法时已经提及，这里再详细阐明其中的两种途径.

途径1——利用定理2.

3.2 小节的具体做法(2)对应于m-s 型或s-m 型互逆规划，其中情况四是：待求解的m 方是不合理的模型，而它对应的s 方是合理的模型.

此时先求解s 方问题(2b)，根据定理2，用其解来调整m 方的参数：取问题(2b)的Lagrange 乘子的最优值作为式(2a)的多目标权系数，取问题(2b)的最优目标值作为式(2a)的约束限.此时作为m 方的式(2a)就可以按帕累托模型进行求解.

途径2——利用定理3.

3.2 小节的具体做法(3)对应于m-m 型互逆规划，其中情况二是：待求解的一方是不合理的模型，而它对应的另一方是合理的模型，不妨假定它们依次为问题(3a)和问题(3b).

此时先着手合理方对应的问题(3b)，先按凝聚函数处理式(3b)目标函数为单目标函数，接着求解.接着根据定理3，用式(3b)的解来调整式(3a)的参数：取问题(3b)的Lagrange 乘子的最优值作为问题(3a)的多目标权系数，问题(3b)的最优的多目标值取作问题(3a)的多约束限.此时的式(3a)可按帕累托模型进行求解.

上述两个途径的共性是：利用各自先求解的互逆规划第一方的最优Lagrange 乘子向量取作第二方的多目标权系数向量，这是本文对于帕累托模型的贡献.两个途径的个性只是在于第一方模型不同而造成的第一方求解差异.

另外，处理后的每一个单目标问题均可以利用DSQP 算法求解.

3.4 利用互逆规划理论求出不合理问题的“合理结果”

文献[26]从互逆规划的角度，揭示了结构拓扑优化里存在有不合理模型的情况，例如给定体积约束求柔顺度最小化的MCVC 问题，本文3.2 小节的求解策略和具体做法也提到不合理模型的情况，并且在具体做法里从3 种类型互逆规划逐一进行了讨论.

本文作者一如既往的观点是：尽量建立和求解合理的规划；不过，也有不得不遭遇计算不合理模型的尴尬，例如：为了呼吁学者们抛弃不合理模型，采用合理的模型，就需要耐心地进行数值比较，既要计算出合理的结果，也要计算出不合理的结果.

为此，本文进行如下的工作：首先，在多工况条件下，取式(2a)代表给定体积约束的多个柔顺度皆最小化的MCVC 问题，取式(2b)代表多个柔顺度约束的体积最小化MVCC 问题.利用ICM 方法分别对于这个互逆规划的双方建立起优化模型.

设单元刚度矩阵和单元体积的过滤函数分别为fk(xi)=和fv(xi)=xi(i=1,2,···,N).

MCVC 模型

其中，x为设计变量向量，Cj，Kj，Fj和Uj分别为第j号载荷工况的结构柔顺度、结构刚度矩阵、总载荷向量和总位移向量，ki j和ui j分别为第j号载荷工况的单元位移列向量和单元刚度矩阵，vi为单元i的体积，V为结构总体积，为指定的体积上限约束.

MVCC模型

其中，Cj是结构柔顺度约束值.

式(22)相比式(23)，是不合理的.原因不仅在于多目标上，更在于所给出作为结构体积上限为指定的值，不具有理性的根据.因此，本文作者多次呼吁以MVCC 问题替代MCVC 问题，而且文献[26]给出了结构柔顺度约束值的物理意义和计算公式.

从学术角度考虑，应当与同行作者的MCVC 问题结果进行比较，导致本文作者不得不求解式(22)，从而探究不合理一方解法，进而利用互逆规划的定理2，计算出不合理问题的“合理结果”.亦即给出如下的建议：

首先求出式(23)的解V∗，然后取式(22)的=V∗，且取式(23)最优的Lagrange 乘子向量作为式(22)的多目标权系数向量，加权求和使之成为单目标问题，当然也可以采用凝聚函数将式(22)化为单目标问题.也就是说，经过这一番由式(23)的解调整式(22)的参数，最后才有可能求解计算出不合理问题的“合理结果”.

上面探讨了对于m-s 型(或s-m 型)互逆规划，借助于定理2，从求解合理的s 方的途径得到不合理的m 方的“合理解”的过程.类似地，对于m-m 型互逆规划，借助于定理3，也可以从求解合理的m 方的途径，调整不合理的m 方的参数，最后求出“合理解”.具体做法相似，此处不再赘述了.

4 数值算例

这里以两个数学问题和两个结构拓扑优化问题，作为本文应用的数值算例.

4.1 算例1:简单的s-s 型数学规划问题

一对s-s 型互逆规划如式(24)所示，其中式(20b)中〈RP 问题〉中的约束值，有待根据定理1，由式(24a)即〈FP 问题〉中的最优目标值确定.

得到KKT 条件

得到退化为鞍点条件形式的KKT 条件

图1 G 与H 函数关系Fig.1 Relationship between G and H function

为了从直观上把握互逆规划，从3 个方面进行如下讨论.

(1) FP 问题的几何意义

先考虑更一般的情况

单约束问题若不退化为无约束问题，约束必定取等式.按KKT 条件叙述，即Lagrange 乘子必不为0，而为正数.因此，式(29)是式(24a)的更一般情况.

目标函数是由大圆向圆心下降的等值圆族，这个结论可以由式(29)的目标函数负梯度的指向得出，也可以从点(3,0),(2.5,1.5)和(3,2)目标函数值−9,−12.5 和−13 看到，它们分别是棕色圆、紫色圆和圆心的目标值.

可行域是绿色直线斜下方区域.

最优点是在紫色圆与绿色直线的相切点，最优目标值等于−12.5.想象在设计变量--目标函数的里，目标函数构成了一个倒置的旋转面，旋转面的锥点高度为−13，在高度为−12.5 处的旋转面的水平截线圆上一点，有一直线过该点，而这个直线的投影与过(4,4)点，该直线就是可行域的边线.

顺便指出：式(29)的半无穷大可行域退化式(24a)等式表示的直线.

(2) RP 问题的几何意义

同上，先考虑比式(24b)更一般的情况

图2 FP 问题的几何直观Fig.2 Geometry of the FP program

目标函数等值线是由与绿色直线平行的直线族组成，垂直于平行直线向左下方下降；可行域是包括紫色圆边线的圆区域.

最优点是在紫色圆与绿色直线的相切点，最优目标值等于4.

在设计变量--目标函数的里，目标函数等值线构成了过(0，0，0),(4，0，4)和(0，4，4)三个点的无限大的平面.

图3 RP 问题的几何直观Fig.3 Geometry of the RP program

顺便指出：式(29)的圆面积可行域退化式(24b)等式表示的圆周线可行域.

(3)参数变动导致结果改变的几何背景

当H∗=−12.5 时，求解式(24b)，可以得到如下结论:

该结论可如下推证.取R为圆心坐标为(3,2)的圆的半径，令

解得

于是，得到下面3 种情况：

4.2 算例2:简单的m-s 型数学规划问题

为求解FP，构造Lagrange 函数

得到KKT 条件

图4 FP 问题的目标函数及可行域Fig.4 Objective function and feasible region of the FP program

因式(28)中仅一个约束，也可以按等式约束处理.为了求解FP，构造Lagrange 函数

得到KKT 条件

为了方便比较结果，图5 中画了式(31b)中的两个目标函数G1及G2最优点对应的直线线，它们对应的目标函数最优值分别为=−3 及=−3.还画了[w1,w2]=[0.1,0.9]情况下最优点(1.25,0.66)对应的直线，其对应的目标函数最优值为G∗=−3.036.它们的最优点与式(28)的最优点在图5 上以2 个黑点表示.

图5 RP 问题的目标函数及可行域Fig.5 Objective function and feasible region of the RP program

4.3 算例3:单工况MBB 梁问题拓扑优化

材料弹性模量E=2.1×105MPa，泊松比为0.3.基结构宽为300 mm、高50 mm，取一半结构进行分析，左、右下角分别为铰支和滑动铰支约束，计算简图如图6 所示，采用150×50 网格.一个集中载荷200 N 载荷工况，作用于跨中顶点，为避免应力集中，F=100 N 分散作用在相邻的4 个结点上，右下角竖向约束也分散在2 个结点上.过滤半径为3 mm，收敛精度取0.001.为了进行结果的比较，涉及到MCVC 模型，惩罚因子取p=3,分别采用20%与60%的体积比约束值，进行两种情况的计算；涉及到MVCC 模型，采用许用应力[σ]=125 MPa，对应的结构应变能约束值为=10.012 N·mm[26].此例计算按如下3 步进行.

图6 MBB 梁问题取一半结构进行分析及优化的模型Fig.6 Half of structure for analysis and optimization of the MBB beam

其一，MCVC 模型计算.

分别采用20%和60%的体积比约束值，最优拓扑图如图7 所示.体积比约束20%时，最大应力148.71 MPa，最优柔顺度13.242 N·mm；体积比约束60%时，最大应力77.60 MPa，最优柔顺度3.889 N·mm.

图7 MCVC 模型得到的最优拓扑Fig.7 Optimized topology obtained by the MCVC model

可以看到MCVC 模型对不同的体积比约束值，不仅各自对应的最大Mises 应力不同，而且更重要的是二者最优拓扑构型也不一样.拓扑构型的不同与最大应力的差别皆源于体积比约束值各异.究竟如何做，才能避免难以选取适当体积比约束值的困扰？为此，进行第二步计算.

其二，MVCC 模型计算.

以体积极小为目标、全局应变能约束计算所得最优拓扑如图8(a)所示.对应的Mises 应力分布如图8(b)所示.最优体积比为23.88%.最大应力为125.07 MPa，满足应力约束条件[26].为了便于第三步的数值比较，算出最优点处结构柔顺度为=10.012 N·mm.

图8 MVCC 模型得到的结果Fig.8 The results obtained by the MVCC model

其三，MCVC 模型再计算.

依据定理1，将MVCC 模型得到的最优体积比23.88%，作为修正的MCVC _模型的体积约束，再计算得到最优点结构柔顺度=10.121 N·mm，最优拓扑与Mises 应力分布如图9 所示，最大应力为124.75 MPa.

图9 体积比23.88%时MCVC 模型得到的结果Fig.9 The results obtained by the MCVC model with 23.88%volume ratio constraint

计算结果的拓扑构型、最大应力和结构柔顺度，都接近MVCC 模型得到的结果.以上三步计算，从MVCC 模型—MVCC 模型—MCVC 模型，恰似一个否定之否定的过程.这个过程与依据“以往MCVC 模型”求解的不同，在于“新MCVC 模型”插入了上述的第二步.可见，依据MVCC 模型可以确定MCVC 模型一个恰当的体积比约束值.也就是说，经过定理1 的改造，不合理的模型求出了合理的结果.可见，这是本文对于MCVC 模型求解的一个贡献.

4.4 算例4:多工况正方形平板问题拓扑优化

材料弹性模量E=2.1×105MPa，泊松比为0.3.基结构宽为100 mm、高50 mm，结构左边为固定约束，计算简图如图10 所示，采用100×50 网格.工况1 为一集中载荷F1=200 N 竖直向上作用于上边中点，工况2 为一集中载荷F2=150 N 竖直向上作用于上边右角，为避免应力集中，F1分散作用在相邻的4 个结点上，F2分散作用在相邻的3 个结点上.过滤半径为3 mm，收敛精度取0.001.为了进行结果的比较，涉及到MCVC 模型，惩罚因子取p=3；涉及到MVCC 模型，采用许用应力[σ]=125 MPa，对应于两个工况的结构应变能约束值为=3.841 4 N·mm 和=7.412 3 N·mm[26].此例计算按如下三步进行.

图10 正方形平板拓扑优化问题模型Fig.10 Definition of the topology optimization problem of a square plate

其一，MCVC 模型计算.

MCVC 模型因工况不同产生的多目标，其加权系数按载荷大小的比值4:3 计算，得到归一化的权系数w1=0.571 4 和w2=0.428 6，分别计算20%和60%的体积比约束值两种情况.体积比约束20%时，两工况下最大应力分别为=163.34 MPa 和=195.40 MPa，均不满足应力约束条件.最优点处结构柔顺度为=5.902 6 N·mm 和=18.688 5 N·mm.体积比约束60%时，两工况下最大应力分别为=56.39 MPa 和=80.58 MPa，均比应力约束值小很多，未能充分发挥材料的强度性能.最优点处结构柔顺度为=1.157 7 N·mm 和=3.297 4 N·mm.最优拓扑如图11 所示.

图11 MCVC 模型得到的最优拓扑Fig.11 Optimized topology obtained by the MCVC model

为了比较，再计算一种情况：保持60%的体积比约束，使两工况的加权系数相同w1=0.5 及w2=0.5，计算得到两工况下最大应力为=58.60 MPa 和=80.58 MPa，同样比应力约束值小很多.最优点处结构柔顺度为=1.206 0 N·mm和=3.264 1 N·mm.

可见，MCVC 模型对不同的体积比约束值，对应与不同的工况的权系数，不仅各自对应的最大Mises 应力不同，而且最优拓扑构型也不一样.原因在于，体积比约束值和多目标加权系数的预估，都没有找到理性的取值方法.为此，导致下一步计算.

其二，MVCC 模型计算.

取结构体积极小化为目标、应力约束对应的全局应变能上限为约束条件[26]，计算所得最优拓扑如图12 所示.对应的Mises 应力分布如图13 所示.最优解的体积比为29.44%.两工况下最优点处结构柔顺度皆正好满足各自的应变能约束条件为=3.841 4 N·mm 和=7.412 3 N·mm；对应的最大应力为=125.00 MPa 和=125.00 MPa，恰好满足应力约束条件[26].此时两工况约束对应的Lagrange 乘子为=22.2135 和=67.318 4.

其三，MCVC 模型再计算.

依据定理2，确定再计算的MCVC 模型体积约束值和多目标的权系数值：将MVCC 模型得到的最优体积比29.44%，作为MCVC 模型的体积约束；依据上述得到的Lagrange 乘子进行归一化处理，得到多目标的权系数w1=0.248 1 和w2=0.751 9.再进行寻优计算，得到两工况下最大应力为=120.65 MPa 和=125.12 MPa，基本满足应力约束条件.最优点处结构柔顺度为=3.777 6 N·mm 和=7.381 0 N·mm，与MVCC 模型的约束值接近.最优拓扑如图14 所示，与MVCC 模型得到的最优拓扑构型相同.对应的Mises 应力分布如图15 所示.