一种基于重采样粒子优化的结构拓扑优化求解方法

本发明涉及一种基于重采样粒子优化的结构拓扑优化求解方法，相比于其它算法可以更好的进行结构的拓扑优化设计，属于计算机技术领域。

结构拓扑优化的基本概念是指在给定的设计空间、支撑条件、载荷条件和某些工艺设计等要求下，确定结构构件的相互连接方式，结构内有无孔洞、孔洞的数量、位置等拓扑形式，使结构能将外载荷传递到支座，同时使结构的某种性能指标达到最优，这个过程称为结构拓扑优化，是结构优化的一种。拓扑优化相对于尺寸优化和形状优化具有更高的设计自由度，可以获得更大的设计空间，具有非常广泛的发展前景。根据优化对象的不同，结构拓扑优化可以分为离散结构拓扑优化与连续体结构拓扑优化两大类。对于拓扑优化问题主要求解方法可以分为确定性方法与不确定性方法两类，但是对于现代结构拓扑优化来说，其结构优化模型往往表现出非线性、模糊性和随机性等特点，这使得传统确定性方法难以求解，所以便出现了许多不确定性方法即智能优化方法。智能优化方法是指通过计算机编程模拟自然现象，模仿动物乃至人类的社会行为和进化机制，从而实现对复杂优化问题求解的一大类方法的统称。粒子优化方法作为一种智能优化方法，其主要思想源于对鸟捕食行为的研究。粒子优化方法具有简单易行、收敛速度快、设置参数少的优点，并且在解决实际问题中展示了其优越性。但是粒子方法也存在明显的问题，便是在求解复杂的多峰问题时易陷入局部最优解。重采样粒子优化方法在粒子方法的基础上增加了重采样过程，提出了粒子聚集度与粒子活跃值两个新概念，当粒子聚集度达到一定值时进行重采样过程，更新低权重粒子速度与位置，用于平衡全局搜索能力与搜索速度，从而解决了粒子算法易陷入局部最优解的问题。考虑到重采样粒子优化方法相比于其它同类优化算法所具有的优势，进一步丰富结构拓扑优化问题的求解方法，本发明提供了一种基于重采样粒子优化方法的结构拓扑优化求解方法。

1、目的

本发明提供一种基于重采样粒子优化的结构拓扑优化的求解方法，以提高对结构拓扑优化问题的求解能力，引入一种求解结构拓扑优化问题的新方法。

2、技术方案

为了实现上述发明目的，本发明采用以下技术方案。

本发明提供一种基于重采样粒子优化的结构拓扑优化的求解方法，主要包括以下几个步骤：

步骤一：建立具体拓扑优化问题的优化模型

所述的优化模型是指根据具体的结构拓扑优化问题所建立的数学模型，其主要包括根据具体优化问题抽象得到的目标函数、要进行优化的控制量以及控制量应满足的约束条件；在本发明中其一般数学表述如下：

min z＝f(x)

s.t.hi(x)＝0 i＝1，2，3，...，m

gj(x)≤0 j＝1，2，3，…，n

其中x为拓扑优化问题中的控制量，f(x)为拓扑优化问题目标函数，hi(x)为控制量应满足的等式约束，m为等式约束个数，gj(x)为控制量应满足的不等式约束，n为不等式约束个数；对于不同的结构拓扑优化问题，其优化目标函数以及约束形式及约束个数具有不同的形式，需根据具体问题进行建立；

步骤二：建立具体拓扑优化问题的适应度函数

所述具体优化问题的适应度函数是指由步骤一中所建立的拓扑优化模型中目标函数以及所有约束条件，采用拉格朗日方法，所建立的可以反映各个设计点优劣程度的函数关系；在本发明中其具体建立方法如下：

(1)对于仅含等式约束的优化问题采用拉格朗日乘子法，其具体建立方法如下：

设拉格朗日乘子为λ，适应度函数为F(x，λ)，则F(x，λ)表达式为：

其中f(x)为目标函数，m为等式约束个数，hi(x)为第i个等式约束；

(2)对于包含不等式约束的优化问题采用拉格朗日对偶法，其具体建立方法如下：

设拉格朗日乘子为λ，则目标函数f(x)对应的拉格朗日对偶函数h(λ)建立公式如下：

λj≥0

其中f(x)为目标函数，n为不等式约束个数，gj(x)为第j个不等式约束；根据上式所建立的适应度函数如下：

F(x)＝max h(λ)；

步骤三：根据适应度函数利用重采样粒子优化进行求解

重采样粒子优化方法具体求解过程如下：

(1)初始化粒子体，确定粒子种中的粒子个数N、生成每个粒子的初始位置坐标x和初始速度矢量v、以及每个粒子的历史最优位置坐标即pbest和体的最优位置坐标即gbest，根据步骤二中建立的适应度函数评价每个粒子的适应度；

所述初始位置坐标x的生成方法为：

xid＝xmin(d)+rand1·(xmax(d)-xmin(d))

其中，xid是第i个粒子第d维的坐标值，xmin(d)和xmax(d)分别是粒子第d维坐标值的下限和上限，rand1是一组0～1之间的随机数；

所述初始速度矢量v的生成方法为：

vid＝xmin(d)+rand2id·(xmax(d)-xmin(d))-xid

其中，vid是第i个粒子第d维的速度值，xmin(d)和xmax(d)以及xid的含义同上，rand2是一组0～1之间的随机数；

所述每个粒子的pbest的初始化方法为：记粒子的初始位置坐标为pbest的初始值；同时根据步骤二中适应度函数求出每个粒子在pbest处的适应度函数值，称之为历史最优值；

所述gbest的初始化方法为：比较上述每个粒子的历史最优值，其中历史最优值最小的粒子的位置坐标为gbest的初始值，其对应的历史最优值为当前的全局最优值；

(2)判断粒子是否满足重采样条件，若满足则计算每个粒子权重，更新低权重粒子速度与位置；

所述重采样条件的的判别标准为粒子聚集度(PAD)是否达到给定值，主要原理为首先在粒子中给定一个随机中心PAD则为粒子中所有粒子与随机中心距离的方差的倒数，其计算公式如下：

其中：N为粒子中粒子数，DISi为第i个粒子与随机给定的中心之间的距离，xij为第i个粒子的第j维位置分量，为随机中心的第j维位置分量；

所述粒子权重值计算方法如下：

其中：qi为对第i个粒子赋予的权值，F(xi)为适应度函数，gbest为当前体最优位置，σ为以F(xi)-gbest为样本计算所得的方差，Qi为第i个粒子归一化后的权值；

所述低权重值粒子速度与位置更新方法为计算每个粒子的权重值，当某个粒子的权值小于给定的阀值qt时，便以新粒子替代原低权重值粒子，新粒子位置坐标确定方法为：

新粒子速度确定方法为：

其中：T为最大迭代次数，t为当前迭代次数，为新引入的粒子速度，vi(t)为原粒子速度，xmin和xmax分别是粒子坐标值的下限和上限，rand3与rand4均是一组0～1之间的随机数，是上述新的粒子位置；

(3)根据粒子适应度函数更新pbest与gbest，更新粒子速度与位置；

所述粒子i的第d维速度更新公式为

所述粒子i的第d维位置更新公式为

其中：为第k次迭代粒子i速度的第d维分量，为第k次迭代粒子i位置的第d维分量，c1、c2为加速度常数，由初始给定，r1、r2为两个取值在[0,1]范围内的随机数；

(4)当迭代次数达到初始给定值时终止迭代，输出gbest，迭代未达给定值时返回(2)；

步骤四：计算目标函数最优值

将步骤三中所得gbest代入到优化模型的目标函数，所得结果即为目标函数的最优值。

通过以上步骤即可利用重采样粒子优化求解得到拓扑优化问题的最优解，即迭代结束时输出的gbest，同时也可计算目标函数最优值；与其它结构拓扑优化方法相比本发明所述方法可以得到更好的结果且求解效率更高，在多次求解的结果稳定性上本发明所述方法也具有优势。

3、本发明的优点和功效

本发明使用了重采样粒子优化解决了结构拓扑优化问题，提供了一种基于重采样粒子优化的结构拓扑优化问题求解方法。其主要优点在于相比于其它结构拓扑优化方法可以得到更好的结果，而且具有更高的求解效率，除此之外在多次求解中本发明使用方法具有更高的结果稳定性。本方法简单实用，实施容易，具有推广应用价值。

图1一种基于重采样粒子优化的结构拓扑优化求解方法。

图2重采样过程流程图。

图3重采样粒子优化流程图。

图4算例一九杆桁架结构图。

图5算例一九杆桁架拓扑构型图。

图6算例一三种方法迭代过程图。

图7算例二梁结构初始设计区域及载荷图。

图8算例二结构拓扑构型图。

图9算例二两种方法迭代过程图。

图中序号、符号、代号说明如下：

pbest为粒子最优位置坐标；gbest为粒子体最优位置坐标；A1～A9为杆件横截面积代号；RPSO为重采样粒子优化；PSO为传统粒子优化；GA为遗传优化。

以下结合附图和两个结构拓扑优化算例对本方法作进一步描述，为了体现本方法在求解能力上的优势，对于相同算例除本方法外还采用了遗传优化方法或传统粒子优化方法进行比较。

本发明提出了一种基于重采样粒子优化的结构拓扑优化问题的求解方法，本方法主要流程见图1，重采样过程具体流程见图2，重采样粒子优化流程如图3所示，其算例如下：

算例一、九杆桁架结构

九杆桁架结构如附图4所示，该问题有九个设计变量，除应力约束外还需保持结构对称关系，具体参数如下：

密度：ρ＝0.283Lbs/in3；

对称关系：A1＝A2；A3＝A4；A5＝A6；A7＝A8；

许用应力：-18≤σi≤20psi，i＝1，2，9；

-12≤σi≤20psi，i＝3，4，5，6，7，8；

载荷：

(1)节点4上，Px＝100kip；节点2上，Py＝-100kip

(2)节点5上，Px＝-100kip；节点2上，Py＝-100kip

根据本发明所述求解方法具体过程如下：

步骤一：建立具体拓扑优化问题的优化模型

根据上述问题描述，所建立的拓扑优化模型如下：

s.t.gj(ti，Ai)≤0 j＝1，2，3，...，J

0≤Ai≤AU

ti＝0/1

其中W为结构质量，i为杆件编号，ti为一个0或1的常数，用于表示对应杆件是否存在，ρi为杆件密度，Ai为杆件截面积，AU为截面积上限，li为杆件长度，gj(ti，Ai)为结构应满足的约束，J为约束个数。

步骤二：建立具体拓扑优化问题的适应度函数

设拉格朗日乘子为λ，则根据步骤一中模型建立的拉格朗日对偶函数h(λ)公式如下：

λj≥0

根据上式建立的适应度函数如下：

F(x)＝max h(λ)

步骤三：根据适应度函数利用重采样粒子优化进行求解

(1)初始化粒子体，确定粒子种中的粒子个数N、生成每个粒子的初始位置坐标x和初始速度矢量v、以及每个粒子的历史最优位置坐标pbest和体的最优位置坐标gbest，根据步骤二中建立的适应度函数评价每个粒子的适应度；

所述初始位置坐标x的生成方法为：

xid＝xmin(d)+rand1·(xmax(d)-xmin(d))

其中，xid是第i个粒子第d维的坐标值，xmin(D)和xmax(D)分别是粒子第d维坐标值的下限和上限，rand1是一组0～1之间的随机数；

所述初始速度矢量v的生成方法为：

vid＝xmin(d)+rand2id·(xmax(d)-xmin(d))-xid

其中，vid是第i个粒子第d维的速度值，xmin(D)和xmax(D)以及xid的含义同上，rand2是一组0～1之间的随机数；

所述每个粒子的pbest的初始化方法为：记粒子的初始位置坐标为pbest的初始值；同时根据步骤二中适应度函数求出每个粒子在pbest处的适应度函数值，称之为历史最优值。

所述gbest的初始化方法为：比较上述每个粒子的历史最优值，其中历史最优值最小的粒子的位置坐标为gbest的初始值，其对应的历史最优值为当前的全局最优值；

(2)判断粒子是否满足重采样条件，若满足则计算每个粒子权重，更新低权重粒子速度与位置；

所述重采样条件的的判别标准为粒子聚集度(PAD)是否达到给定值，主要原理为首先在粒子中给定一个随机中心PAD则为粒子中所有粒子与随机中心距离的方差的倒数，其计算公式如下：

其中：N为粒子中粒子数，DISi为第i个粒子与随机给定的中心之间的距离，xij为第i个粒子的第j维位置分量，为随机中心的第j维位置分量；

所述粒子权重值计算方法如下：

其中：qi为对第i个粒子赋予的权值，F(xi)为适应度函数，gbest为当前体最优位置，σ为以F(xi)-gbest为样本计算所得的方差。Qi为第i个粒子归一化后的权值；

所述低权重值粒子速度与位置更新方法为计算每个粒子的权重值，当某个粒子的权值小于给定的阀值qt时，便以新粒子替代原低权重值粒子，新粒子位置坐标确定方法为：

新粒子速度确定方法为：

其中：T为最大迭代次数，t为当前迭代次数，为新引入的粒子速度，vi(t)为原粒子速度，xmin和xmax分别是粒子坐标值的下限和上限，rand3与rand4均是一组0～1之间的随机数，是上述新的粒子位置。

(3)根据粒子适应度函数更新pbest与gbest，更新粒子速度与位置；

所述粒子i的第d维速度更新公式为

所述粒子i的第d维位置更新公式为

其中：为第k次迭代粒子i速度的第d维分量，为第k次迭代粒子i位置的第d维分量，c1、c2为加速度常数，由初始给定，r1、r2为两个取值在[0,1]范围内的随机数；

(4)当迭代次数达到初始给定值时终止迭代，输出gbest，迭代未达给定值时返回(2)；

步骤四：求解结构质量并使用其它方法求解该问题

根据(3)中所得gbest代入目标函数求得结构质量，并将该适应度函数代入传统粒子优化(PSO)与遗传优化(GA)进行计算，得到的结构拓扑构型如附图5所示，三种方法迭代过程曲线如附图6所示，为提高计算的精确性对每种方法进行了十次重复计算取得的设计变量数值结果最优值如下表所示：

算例二、梁柔度优化问题

在该拓扑优化问题中，假设设计区域为矩形且被平方有限元离散化，边界形式与载荷条件如附图7所示，网格划分方式为x＝60，y＝20，主要约束为优化结果的体积为原结构体积的一半，即体积分数为f＝0.5。

根据本发明所述求解方法具体过程如下：

步骤一：建立具体拓扑优化问题的优化模型

根据上述问题描述，所建立的拓扑优化模型如下：

s.t.V(x)/V0＝f

KU＝F

0≤xmin≤x≤1

其中：U和F为整体位移和外载荷矢量，K为整体刚度矩阵，x是设计变量即相对密度，p为惩罚因子，在此例中取p＝3，ue和ke为各元素的位移矢量和刚度矩阵，N为设计区域离散化元素数目，在此例中N＝60×20，xmin为设计变量最小值，V(x)与V0为材料体积和设计区域总体积，f为给定体积分数，在此例中f＝0.5。

步骤二：建立具体拓扑优化问题的适应度函数

从优化模型中可以看出，该优化问题中外力约束即KU＝F在求解目标函数时可自动满足，控制量x的范围约束0≤xmin≤x≤1也可以在重采样粒子优化方法中通过给定粒子搜索范围来进行约束，所以在建立适应度函数时可以仅考虑体积约束，因为体积约束为等式约束所以采用拉格朗日乘子法进行建立。体积约束变形得到下式：

V(x)-fV0＝0

设拉格朗日乘子为λ，建立的拉格朗日函数F(x，λ)如下式：

原模型转化为求解函数F(x，λ)最小值问题，该拓扑优化问题的适应度函数即为F(x，λ)。

步骤三：根据适应度函数利用重采样粒子优化进行求解

(1)初始化粒子体，确定粒子种中的粒子个数N、生成每个粒子的初始位置坐标x和初始速度矢量v、以及每个粒子的历史最优位置坐标pbest和体的最优位置坐标gbest，根据步骤二中建立的适应度函数评价每个粒子的适应度；

所述初始位置坐标x的生成方法为：

xid＝xmin(d)+rand1·(xmax(d)-xmin(d))

其中，xid是第i个粒子第d维的坐标值，xmin(d)和xmax(d)分别是粒子第d维坐标值的下限和上限，rand1是一组0～1之间的随机数；

所述初始速度矢量v的生成方法为：

vid＝xmin(d)+rand2id·(xmax(d)-xmin(d))-xid

其中，vid是第i个粒子第d维的速度值，xmin(d)和xmax(d)以及xid的含义同上，rand2是一组0～1之间的随机数；

所述每个粒子的pbest的初始化方法为：记粒子的初始位置坐标为pbest的初始值；同时根据步骤二中适应度函数求出每个粒子在pbest处的适应度函数值，称之为历史最优值。

所述gbest的初始化方法为：比较上述每个粒子的历史最优值，其中历史最优值最小的粒子的位置坐标为gbest的初始值，其对应的历史最优值为当前的全局最优值；

(2)判断粒子是否满足重采样条件，若满足则计算每个粒子权重，更新低权重粒子速度与位置；

所述重采样条件的的判别标准为粒子聚集度(PAD)是否达到给定值，主要原理为首先在粒子中给定一个随机中心PAD则为粒子中所有粒子与随机中心距离的方差的倒数，其计算公式如下：

其中：N为粒子中粒子数，DISi为第i个粒子与随机给定的中心之间的距离，xij为第i个粒子的第j维位置分量，为随机中心的第j维位置分量；

所述粒子权重值计算方法如下：

其中：qi为对第i个粒子赋予的权值，F(xi)为适应度函数，gbest为当前体最优位置，σ为以F(xi)-gbest为样本计算所得的方差。Qi为第i个粒子归一化后的权值；

所述低权重值粒子速度与位置更新方法为计算每个粒子的权重值，当某个粒子的权值小于给定的阀值qt时，便以新粒子替代原低权重值粒子，新粒子位置坐标确定方法为：

新粒子速度确定方法为：

其中：T为最大迭代次数，t为当前迭代次数，为新引入的粒子速度，vi(t)为原粒子速度，xmin和xmax分别是粒子坐标值的下限和上限，rand3与rand4均是一组0～1之间的随机数，是上述新的粒子位置。

(3)根据粒子适应度函数更新pbest与gbest，更新粒子速度与位置；

所述粒子i的第d维速度更新公式为

所述粒子i的第d维位置更新公式为

其中：为第k次迭代粒子i速度的第d维分量，为第k次迭代粒子i位置的第d维分量，c1、c2为加速度常数，由初始给定，r1、r2为两个取值在[0,1]范围内的随机数；

(4)当迭代次数达到初始给定值时终止迭代，输出gbest，迭代未达给定值时返回(2)；

步骤四：求解结构质量并使用其它方法求解该问题

根据(3)中所得gbest将全局最优解以图形的形式输出，得到的拓扑构型如附图8所示，同时将该适应度函数代入传统粒子优化(PSO)与优化准则法(OC)进行计算，两种方法十次计算取最优值迭代过程如附图9所示，最终结果如下表所示：

通过以上所述方法可以有效处理结构拓扑优化问题，而且通过以上两个算例结果显示，使用本发明所述方法求解拓扑优化问题相比于遗传优化方法和传统粒子优化方法可以取得更好的结果。

上面对本专利的实施方式作了详细说明，但是本专利并不限于上述实施方式，在本领域的技术人员所具备的知识范围内，还可以在不脱离本专利宗旨的前提下做出各种变化。