CN103412997B - 一种不确定载荷下的桁架结构稳健设计方法 - Google Patents

Abstract

一种不确定载荷下的桁架结构稳健设计方法，它有九大步骤：一：建立桁架结构稳健优化模型；二：根据结构所承受的不确定载荷，选取若干基础工况；三：计算不确定载荷在基础工况坐标系中的坐标的二阶和四阶中心矩；四：初始化设计变量；五：对各个基础工况进行有限元分析得到相应的结构位移向量；六：计算结构柔度的均值与标准差；七：对结构进行灵敏度分析；八：利用渐近线方法更新设计变量；九：判断是否满足终止条件，若满足，则停止迭代；否则重复步骤五至步骤八直到满足迭代终止条件。本发明较好的解决了载荷不确定条件下以结构柔度均值与标准差的加权和最小为目标进行桁架结构设计的技术问题。

Description

一种不确定载荷下的桁架结构稳健设计方法

技术领域：

本发明涉及一种不确定载荷下的桁架结构稳健设计方法，属于桁架设计技术领域。

背景技术：

载荷不确定性是对结构系统影响最大的不确定性之一，相对于确定性优化设计方法，考虑载荷不确定性的稳健优化设计方法得到的桁架结构可以更好的抵抗外载荷的变化。不确定载荷下的桁架结构稳健设计方法的技术难点是如何准确高效的计算桁架结构在给定不确定载荷作用下的柔度均值与标准差以及目标函数对设计变量的灵敏度。

文献“MiguelCarrasco,BenjaminIvorra,AngelManuelRamos(2012)AVariance-ExpectedComplianceModelforStructuralOptimization.JOptimTheoryAppl(2012)152:136-151”公开了一种载荷不确定条件下的桁架结构稳健优化方法，该方法给出了结构柔度均值与方差的解析表达式。文献公开的方法虽然可用于结构稳健拓扑优化，但是其采用载荷分量的联合概率密度函数表示载荷的不确定性，并且假设载荷分量服从多元正态分布，此方法很难推广到其他情形，因此限制了其在实际工程中的应用；

文献“HAliciaKim,RobertAGuyer(2013)RobustTopologyOptimisationwithGeneralisedProbabilityDistributionofLoading.In:54^thAIAA/ASME/ASCE/AHS/ASCStructures,StructuralDynamics,andMaterialsConference”公开了一种载荷不确定条件下的连续体结构稳健优化方法，该方法采用载荷的大小和作用方向角的联合概率密度函数表示载荷的不确定性，并给出了结构柔度均值与方差的解析表达式。文献公开的方法虽然可用于结构稳健拓扑优化，但是其给出的柔度均值与方差的表达式冗长且繁琐，并且很难推广到三维情形和载荷相关的情形。

发明内容：

1、目的：为了避免解析方法在实际工程中的使用限制问题，本发明的目的是提供一种不确定载荷下的桁架结构稳健设计方法，该方法基于线性结构的位移叠加原理和蒙特卡罗方法计算结构的柔度均值与标准差，并且在每次迭代中基于实对称矩阵的对角化构造位移向量进行灵敏度分析。该方法每次迭代中只需要对选取的少量基础载荷工况进行结构分析，蒙特卡罗抽样彻底从结构柔度均值与标准差的计算过程中分离出来，灵敏度分析所需的位移向量亦通过叠加原理获得，无需进行额外的结构分析，因此计算量小，而且由于采用蒙特卡罗方法进行结构柔度均值与标准差的计算，因此实施简单，适用范围更广。

2、技术方案：本发明解决其技术问题所采用的技术方案是：

一种不确定载荷下的桁架结构稳健设计方法，其特征在于包括以下步骤：

步骤一：建立桁架结构稳健优化模型：

$\underset{x}{\min }:J=\mathrm{\α\μ}\left(c\right)+\mathrm{\β\σ}\left(c\right)$

s.t.:Ku(ω)=f(ω)(ω∈Θ)

$g=V-V_{\max }=\underset{e=1}{\overset{m}{\mathrm{\Σ}}}{l}_e{x}_e-V_{\max }\≤0$

x_l≤x≤x_u

其中J为目标函数，μ(c)和σ(c)分别为结构柔度的均值与标准差，α和β为两个非负实数且满足α+β=1；K为结构整体刚度矩阵，f(ω)和u(ω)分别为结构的载荷与位移向量，ω用于表示载荷的不确定性；m为桁架基础结构中杆件的数目，l_e和x_e分别为杆件e的长度和截面积，V_max为给定的材料体积用量；x为由x_e(e=1,…，m)组成的向量，x_l和x_u分别为由x_e(e=1,…，m)的下限和上限组成的向量；

步骤二：根据结构所承受的不确定载荷，选取若干基础工况；

步骤三：计算不确定载荷在基础工况坐标系中的坐标的二阶和四阶中心矩；

步骤四：初始化设计变量；

步骤五：对各个基础工况进行有限元分析得到相应的结构位移向量；

步骤六：计算结构柔度的均值与标准差；

步骤七：对结构进行灵敏度分析；

步骤八：利用渐近线方法更新设计变量；

步骤九：判断是否满足终止条件，若满足，则停止迭代；否则重复步骤五至步骤八直到满足迭代终止条件。

其中，在步骤一中所述的“建立桁架结构稳健优化模型”的依据和方法为：在预定的载荷条件下，桁架结构的柔度越小，其承受外载荷的能力越强。在载荷不确定条件下，结构的柔度会随着载荷的变化而变化，设计桁架结构时不仅要减小结构柔度的均值μ(c)，还要减小结构柔度随载荷变化而引起的波动，即减小结构柔度的标准差σ(c)。文献“MiguelCarrasco,BenjaminIvorra,AngelManuelRamos(2012)AVariance-ExpectedComplianceModelforStructuralOptimization.JOptimTheoryAppl(2012)152:136-151”和“HAliciaKim,RobertAGuyer(2013)RobustTopologyOptimisationwithGeneralisedProbabilityDistributionofLoading.In:54^thAIAA/ASME/ASCE/AHS/ASCStructures,StructuralDynamics,andMaterialsConference”中均以桁架结构柔度均值μ(c)与方差σ²(c)的加权和最小为目标建立优化模型，其缺点是μ(c)与σ²(c)具有不同的量纲。本发明以桁架结构柔度均值μ(c)与标准差σ(c)的加权和J=αμ(c)+βσ(c)最小为目标建立优化模型，数学表达简单且物理意义更加明确。

其中，在步骤二中所述的“选取若干基础工况”的方法为：对第i个不确定载荷，选择两个位移边界条件均与实际工况相同的基础工况f_2i-1和f_2i；其中，在第一种工况中f_2i-1，仅在该不确定载荷作用点沿x方向施加一单位大小的集中力，在第二种工况f_2i中，仅在该不确定载荷作用点沿y方向施加一单位大小的集中力，由此选取出2n个基础工况f₁,…,f_2n

其中，在步骤三中所述的“计算不确定载荷在基础工况坐标系中的坐标的二阶和四阶中心矩”的方法为：对于任意的载荷f(ω)，若其中第i个不确定载荷的大小为h_i(ω)，作用方向与x方向的夹角为θ_i，则其中ξ_2i-1(ω)=h_i(ω)cos(θ_i(ω))，ξ_2i(ω)=h_i(ω)sin(θ_i(ω))。不确定载荷在基础工况坐标系中的坐标的二阶中心矩分别为ξ_ij=E(ξ_i(ω)ξ_j(ω))(i,j=1,…,2n)，四阶中心矩分别为ξ_ijkl=E(ξ_i(ω)ξ_j(ω)ξ_k(ω)ξ_l(ω))(i,j,k,l=1,…,2n)。利用蒙特卡罗方法计算ξ_ij和ξ_ijkl，给定充分大的样本数目N并根据不确定载荷的概率分布抽取N个样本ω_p(p=1,…,N)，则ξ_ij与ξ_ijkl近似值分别由下面两式给出：

${\mathrm{\ξ}}_{\mathrm{ij}}=\frac{1}{N}\underset{p=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ξ}}_i\left({\mathrm{\ω}}_p\right){\mathrm{\ξ}}_j\left({\mathrm{\ω}}_p\right)$

${\mathrm{\ξ}}_{\mathrm{ijkl}}=\frac{1}{N}\underset{p=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ξ}}_i\left({\mathrm{\ω}}_p\right){\mathrm{\ξ}}_j\left({\mathrm{\ω}}_p\right){\mathrm{\ξ}}_k\left({\mathrm{\ω}}_p\right){\mathrm{\ξ}}_l\left({\mathrm{\ω}}_p\right)$

其中，在步骤四中所述的“初始化设计变量”的方法为：将所有的杆件取相同的截面积 $V_{\max }/\underset{e=1}{\overset{m}{\mathrm{\Σ}}}{l}_e.$

其中，在步骤五中所述的“有限元分析”可以采用任何商业有限元分析软件或编写的程序实现，如Nastran，Ansys等。有限元分析的结果是各个基础载荷工况f₁,…,f_2n下结构的位移向量u₁,…,u_2n。

其中，在步骤六中所述的“计算结构柔度均值与标准差”的方法为：结构柔度的均值为其中表示载荷f_i在位移u_j上所作之功。结构柔度的方差为 ${\mathrm{\σ}}^2\left(c\right)=\underset{i,j,k,l=1}{\overset{2n}{\mathrm{\Σ}}}({\mathrm{\ξ}}_{\mathrm{ijkl}}-{\mathrm{\ξ}}_{\mathrm{ij}}{\mathrm{\ξ}}_{\mathrm{kl}})c_{\mathrm{ij}}{c}_{\mathrm{kl}}.$

其中，在步骤七中所述的“灵敏度分析”的方法为：根据结构柔度均值与方差的表达式，目标函数对于x_e的偏导数为：

$\frac{\∂J}{{\∂x}_e}=\underset{i,j=1}{\overset{2n}{\mathrm{\Σ}}}(\mathrm{\α}{\mathrm{\ξ}}_{\mathrm{ij}}+\frac{\mathrm{\β}}{\mathrm{\σ}\left(c\right)}\underset{k,l=1}{\overset{2n}{\mathrm{\Σ}}}({\mathrm{\ξ}}_{\mathrm{ijkl}}-{\mathrm{\ξ}}_{\mathrm{ij}}{\mathrm{\ξ}}_{\mathrm{kl}})c_{\mathrm{kl}})\frac{{\∂c}_{\mathrm{ij}}}{{\∂x}_e}$

记并以w_ij为元素构造矩阵W=(w_ij)。对W进行对角化分解得到Q^TWQ=diag{λ₁,…,λ_2n}，其中λ_i(i=1,…,2n)为W的特征值，Q为由相应的单位特征向量组成的正交阵。设q_ij为正交阵Q位于i行j列的元素，构造新的载荷工况和相应的位移向量则目标函数对设计变量x_e的偏导数由下式给出：

$\frac{\∂J}{{\∂x}_e}=\underset{i=1}{\overset{2n}{\mathrm{\Σ}}}{\mathrm{\λ}}_i\frac{\∂\left({F_i}^T{U}_i\right)}{{\∂x}_e}=\underset{i=1}{\overset{2n}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{U_{\mathrm{ie}}}^T\frac{\∂K}{{\∂x}_e}{U}_{\mathrm{ie}}$

其中U_ie为杆件e在载荷工况F_i下的单元位移向量。

约束对于设计变量的灵敏度为：

$\frac{\∂g}{{\∂x}_e}=l_e$

其中，在步骤八中所述的“利用渐近线方法更新设计变量”的方法具体实现过程如下：首先根据步骤七中灵敏度分析获得的灵敏度信息构造目标函数的近似函数：

$f\left(x\right)\≈f\left(x^{\left(k\right)}\right)+\underset{e=1}{\overset{m}{\mathrm{\Σ}}}{p}_e^{\left(k\right)}(\frac{1}{U_e^{\left(k\right)}-x_e}-\frac{1}{U_e^{\left(k\right)}-x_e^{\left(k\right)}})+\underset{e=1}{\overset{m}{\mathrm{\Σ}}}{q}_e^{\left(k\right)}(\frac{1}{x_e-L_e^{\left(k\right)}}-\frac{1}{x_e^{\left(k\right)}-L_e^{\left(k\right)}})$

其中x^(k)为当前迭代点，和为移动渐近线，其迭代格式为

$L_e^{\left(k\right)}=x_e^{\left(k\right)}-0.1(x_e^{\max }-x_e^{\min })$

$U_e^{\left(k\right)}=x_e^{\left(k\right)}-0.1(x_e^{\max }-x_e^{\min })$

为设计变量x_e的当前值，和分别为x_e的上下限； $p_e^{\left(k\right)}=\max ({(U_e^{\left(k\right)}-x_e^{\left(k\right)})}^2\∂f/{\∂x}_e,0),$ 然后利用该近似函数代替原优化问题中的目标函数而形成优化子问题；最后利用对偶方法求解该优化子问题而得到更新后的设计变量。

其中，在步骤九中所述的“终止条件”包括迭代次数和各个杆件连续两次迭代截面积变化量的最大值两方面。若迭代次数达到了最大次数或者所有杆件连续两次迭代截面积变化量均不超过给定值，则停止迭代；否则继续迭代过程。

3、优点及功效：本发明的有益效果是：基于线性结构的位移叠加原理，每次迭代只需要对少量的基础载荷工况进行结构分析；由于采用蒙特卡罗方法计算结构柔度的均值与标准差，因此易于实施并且适用范围更广，不再仅限于多元正态分布不确定性；通过实对称矩阵的对角化构造位移向量进行灵敏度分析，计算量正比于不确定载荷的数目，比参考文献中的方法计算量小，因此可应用于不确定载荷数量很大的情形。可以简单的推广到三维桁架结构的稳健设计，因此对实际工程中桁架结构设计具有更强的适用性。

附图说明

图1.不确定载荷下的桁架结构稳健设计流程图；

图2.桁架基础结构及边界与载荷条件；

图3.基础载荷工况f₁；

图4.基础载荷工况f₂；

图5.基础载荷工况f₃；

图6.基础载荷工况f₄；

图7.基础载荷工况f₅；

图8.基础载荷工况f₆；

图9.不考虑载荷不确定性的优化结果；

图10.考虑载荷不确定性（α=0.5,β=0.5）的优化结果；

图中符号说明如下：

P₁,P₂,P₃为3个大小和方向均不确定的集中载荷；

F_1x和F_1y分别为在不确定载荷P₁作用位置沿x方向和沿y方向施加的单位大小的集中力；

F_2x和F_2y分别为在不确定载荷P₂作用位置沿x方向和沿y方向施加的单位大小的集中力；

F_3x和F_3y分别为在不确定载荷P₃作用位置沿x方向和沿y方向施加的单位大小的集中力。

具体实施方式

下面结合附图和实施例对本发明作详细说明

以图2所示桁架的稳健设计为例说明本发明。桁架的基础结构由15个节点和74根杆件组成，其中结构下方的两个顶点固定，并在下端承受3个不确定载荷P₁,P₂,P₃。载荷的不确定性由载荷大小h_i和作用方向角θ_i表示，并且h_i与θ_i均服从正态分布，h_i的均值分别为5.0,30.0和5.0，标准差分别为1.0,6.0和1.0，θ_i的均值均为-π/2，标准差均为π/6。杆件材料的杨氏模量为100，泊松比为0.3，材料总体积为9.0，采用一致的单位系统。对该桁架结构进行稳健设计优化，使其在给定不确定载荷条件下的结构柔度均值与标准差的加权和最小。

见图1，本发明一种不确定载荷下的桁架结构稳健设计方法，该方法具体步骤如下：

步骤一：建立桁架结构稳健优化模型：

$\underset{x}{\min }:J=\mathrm{\α\μ}\left(c\right)+\mathrm{\β\σ}\left(c\right)$

s.t.:Ku(ω)=f(ω)(ω∈Θ)

$g=V-V_{\max }=\underset{e=1}{\overset{74}{\mathrm{\Σ}}}{l}_e{x}_e-V_{\max }\≤0$

x_l≤x≤x_u

其中J为目标函数，μ(c)和σ(c)分别为结构柔度的均值与标准差，α=β=0.5；K为结构整体刚度矩阵，f(ω)和u(ω)分别为结构的载荷与位移向量，ω用于表示载荷的不确定性；m=74为桁架基础结构中杆件的数目，l_e和x_e分别为杆件e的长度和截面积，V_max=9.0为给定的材料体积用量；x为由x_e(e=1,…,74)组成的向量，x_e(e=1,…,74)的上下限分别为0.01和0.0001，因此x_l和x_u均为由74个元素构成的向量，其中x_l的所有元素均为0.0001，而x_u的所有元素均为0.01。

步骤二：选取6个基础工况f₁,…,f₆。如图3—图8所示，对每个不确定载荷P_i(i=1,2,3)，选择两个位移边界条件均与实际工况相同的基础工况f_2i-1与f_2i，其中f_2i-1为仅在该不确定载荷作用点沿x方向施加一单位大小的集中力，f_2i为仅在该不确定载荷作用点沿y方向施加一单位大小的集中力。

步骤三：计算不确定载荷在基础工况坐标系中的坐标的二阶和四阶中心矩。三个不确定载荷P₁,P₂,P₃的大小h₁,h₂,h₃服从的概率密度函数分别为：

$p\left(h_1\right)=\frac{1}{\sqrt{2\mathrm{\π}}}\exp (-{(h_1-5)}^2/2)$

$p\left(h_2\right)=\frac{1}{36\sqrt{2\mathrm{\π}}}\exp (-{(h_2-30)}^2/72)$

$p\left(h_3\right)=\frac{1}{\sqrt{2\mathrm{\π}}}\exp (-{(h_3-5)}^2/2)$

三个不确定载荷P₁,P₂,P₃作用方向与x方向的夹角θ₁,θ₂,θ₃服从的概率密度函数分别：

$p\left({\mathrm{\θ}}_i\right)=\frac{6}{\mathrm{\π}\sqrt{2\mathrm{\π}}}\exp (-36{({\mathrm{\θ}}_i+\frac{\mathrm{\π}}{2})}^2/{\mathrm{\π}}^2)(i=\mathrm{1,2,3})$

在本例中蒙特卡罗法的抽样次数取为N=1000000，利用matlab从h₁,h₂,h₃,θ₁,θ₂,θ₃的概率分布中抽取N个样本ω_p(p=1,…,N)，ξ_ij与ξ_ijkl的近似值分别由下面两式给出：

${\mathrm{\ξ}}_{\mathrm{ij}}=\frac{1}{N}\underset{p=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ξ}}_i\left({\mathrm{\ω}}_p\right){\mathrm{\ξ}}_j\left({\mathrm{\ω}}_p\right)$

${\mathrm{\ξ}}_{\mathrm{ijkl}}=\frac{1}{N}\underset{p=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\ξ}}_i\left({\mathrm{\ω}}_p\right){\mathrm{\ξ}}_j\left({\mathrm{\ω}}_p\right){\mathrm{\ξ}}_k\left({\mathrm{\ω}}_p\right){\mathrm{\ξ}}_l\left({\mathrm{\ω}}_p\right)$

步骤四：初始化设计变量。取每个杆件的截面积x_e=0.000025(e=1,…,74)；

步骤五：利用有限元软件Nastran对基础工况f₁,…,f₆进行有限元分析得到相应的结构位移向量u₁,…,u₆；

步骤六：结构柔度的均值为其中结构柔度的方差为 ${\mathrm{\σ}}^2\left(c\right)=\underset{i,j,k,l=1}{\overset{6}{\mathrm{\Σ}}}({\mathrm{\ξ}}_{\mathrm{ijkl}}-{\mathrm{\ξ}}_{\mathrm{ij}}{\mathrm{\ξ}}_{\mathrm{kl}})c_{\mathrm{ij}}{c}_{\mathrm{kl}}.$

步骤七：进行结构灵敏度分析。目标函数对于x_e的偏导数为：

$\frac{\∂J}{{\∂x}_e}=\underset{i,j=1}{\overset{6}{\mathrm{\Σ}}}(\mathrm{\α}{\mathrm{\ξ}}_{\mathrm{ij}}+\frac{\mathrm{\β}}{\mathrm{\σ}\left(c\right)}\underset{k,l=1}{\overset{6}{\mathrm{\Σ}}}({\mathrm{\ξ}}_{\mathrm{ijkl}}-{\mathrm{\ξ}}_{\mathrm{ij}}{\mathrm{\ξ}}_{\mathrm{kl}})c_{\mathrm{kl}})\frac{{\∂c}_{\mathrm{ij}}}{{\∂x}_e}$

记并以w_ij为元素构造矩阵W=(w_ij)。对W进行对角化分解得到Q^TWQ=diag{λ₁,…,λ₆}，其中λ_i(i=1,…,6)为W的特征值，Q为由相应的单位特征向量组成的正交阵。设q_ij为正交阵Q位于i行j列的元素，构造新的载荷工况和相应的位移向量则目标函数对设计变量x_e的偏导数由下式给出：

$\frac{\∂J}{{\∂x}_e}=\underset{i=1}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\λ}}_i\frac{\∂\left({F_i}^T{U}_i\right)}{{\∂x}_e}=\underset{i=1}{\overset{6}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{U_{\mathrm{ie}}}^T\frac{\∂K}{{\∂x}_e}{U}_{\mathrm{ie}}$

其中U_ie为杆件e在载荷工况F_i下的单元位移向量。

约束对于设计变量的灵敏度为：

$\frac{\∂g}{{\∂x}_e}=l_e$

步骤八：利用渐近线方法更新设计变量；首先根据步骤七中灵敏度分析获得的灵敏度信息构造目标函数的近似函数：

$f\left(x\right)\≈f\left(x^{\left(k\right)}\right)+\underset{i=e}{\overset{74}{\mathrm{\Σ}}}{p}_e^{\left(k\right)}(\frac{1}{U_e^{\left(k\right)}-x_e}-\frac{1}{U_e^{\left(k\right)}-x_e^{\left(k\right)}})+\underset{e=1}{\overset{74}{\mathrm{\Σ}}}{q}_e^{\left(k\right)}(\frac{1}{x_e-L_e^{\left(k\right)}}-\frac{1}{x_e^{\left(k\right)}-L_e^{\left(k\right)}})$

其中x^(k)为当前迭代点，和为移动渐近线，其迭代格式为

$L_e^{\left(k\right)}=x_e^{\left(k\right)}-0.1(x_e^{\max }-x_e^{\min })$

$U_e^{\left(k\right)}=x_e^{\left(k\right)}-0.1(x_e^{\max }-x_e^{\min })$

为设计变量x_e的当前值，和分别为x_e的上下限； $p_e^{\left(k\right)}=\max ({(U_e^{\left(k\right)}-x_e^{\left(k\right)})}^2\∂f/{\∂x}_e,0),$ 然后利用该近似函数代替原优化问题中的目标函数而形成优化子问题；最后利用对偶方法求解该优化子问题而得到更新后的设计变量。

步骤九：判断是否满足终止条件，若满足，则停止迭代；否则重复步骤二至步骤五直到迭代终止。本例中迭代终止的条件是相邻两次迭代所有杆件横截面积的变化量均不超过0.000001或者迭代次数超过1000.

本例中进行了两种条件下的桁架优化，其中第一种不考虑载荷不确定性，即按照确定性优化进行桁架设计；第二种利用本发明的桁架稳健设计方法，其中α=0.5,β=0.5。两种情况下桁架设计的结果分别如图9，图10所示。两种结构中各杆件的截面积是不同的，其中在图10中考虑载荷不确定性的稳健设计得到的桁架结构最下方的四根杆件的截面积比较大。这是由于载荷的不确定性，三个载荷都可能存在沿x方向的分量，该分量的存在使得确定性优化结果变得不稳定，本发明的稳健设计方法很好的考虑了这一点，增强了下方四根具有较强的水平方向承载能力的杆件，因此所得到的结构也具有更好的抵抗外载荷变化的能力。本例中每次迭代只需要对6个基础工况进行结构分析，结构柔度的均值与方差计算以及灵敏度分析所需要的位移向量均通过对以上6个基础工况的位移向量进行线性运算得到，因此额外的计算量很小，计算效率很高。

Claims (3)

1.一种不确定载荷下的桁架结构稳健设计方法，其特征在于：它包括以下步骤：

步骤一：建立桁架结构稳健优化模型：

$\underset{x}{Min}:J=\mathrm{\α}\mathrm{\μ}\left(c\right)+\mathrm{\β}\mathrm{\σ}\left(c\right)$

s.t.:Ku(ω)＝f(ω)(ω∈Θ)

$g=V-V_{max}=\underset{e=1}{\overset{m}{\Σ}}{l}_e{x}_e-V_{max}\≤0$

x_l≤x≤x_u

其中，J为目标函数，μ(c)和σ(c)分别为结构柔度的均值与标准差，α和β为两个非负实数且满足α+β＝1；K为结构整体刚度矩阵，f(ω)和u(ω)分别为结构的载荷与位移向量，ω用于表示载荷的不确定性；m为桁架基础结构中杆件的数目，l_e和x_e分别为杆件e的长度和截面积，V_max为给定的材料体积用量；x为由x_e组成的向量，x_l和x_u分别为由x_e的下限和上限组成的向量；其中，e＝1,…,m

步骤二：根据结构所承受的不确定载荷，选取若干基础工况；

步骤三：计算不确定载荷在基础工况坐标系中的坐标的二阶和四阶中心矩；

步骤四：初始化设计变量；

步骤五：对各个基础工况进行有限元分析得到相应的结构位移向量；

步骤六：计算结构柔度的均值与标准差；

步骤七：对结构进行灵敏度分析；

步骤八：利用渐近线方法更新设计变量；

步骤九：判断是否满足终止条件，若满足，则停止迭代；否则重复步骤五至步骤八直到满足迭代终止条件；

其中，在步骤二中所述的选取若干基础工况的方法为：对第i个不确定载荷，选择两个位移边界条件均与实际工况相同的基础工况f_2i-1和f_2i；其中，在第一种工况中f_2i-1，仅在该不确定载荷作用点沿x方向施加一单位大小的集中力，在第二种工况f_2i中，仅在该不确定载荷作用点沿y方向施加一单位大小的集中力，由此选取出2n个基础工况f₁,…,f_2n；

其中，在步骤三中所述的计算不确定载荷在基础工况坐标系中的坐标的二阶和四阶中心矩的方法为：对于任意的载荷f(ω)，若其中第i个不确定载荷的大小为h_i(ω)，作用方向与x方向的夹角为θ_i，则

其中ξ_2i-1(ω)＝h_i(ω)cos(θ_i(ω))，ξ_2i(ω)＝h_i(ω)sin(θ_i(ω))；不确定载荷在基础工况坐标系中的坐标的二阶中心矩分别为ξ_ij＝E(ξ_i(ω)ξ_j(ω))，其中，i,j＝1,…,2n，四阶中心矩分别为ξ_ijkl＝E(ξ_i(ω)ξ_j(ω)ξ_k(ω)ξ_l(ω)，)其中i,j,k,l＝1,…,2n，利用蒙特卡罗方法计算ξ_ij和ξ_ijkl，给定充分大的样本数目N并根据不确定载荷的概率分布抽取N个样本ω_p，其中，p＝1,…,N，则ξ_ij与ξ_ijkl近似值分别由下面两式给出：

${\mathrm{\ξ}}_{ij}=\frac{1}{N}\underset{p=1}{\overset{N}{\Σ}}{\mathrm{\ξ}}_i\left({\mathrm{\ω}}_p\right){\mathrm{\ξ}}_j\left({\mathrm{\ω}}_p\right)$

${\mathrm{\ξ}}_{ijkl}=\frac{1}{N}\underset{p=1}{\overset{N}{\Σ}}{\mathrm{\ξ}}_i\left({\mathrm{\ω}}_p\right){\mathrm{\ξ}}_j\left({\mathrm{\ω}}_p\right){\mathrm{\ξ}}_k\left({\mathrm{\ω}}_p\right){\mathrm{\ξ}}_l\left({\mathrm{\ω}}_p\right);$

其中，在步骤四中所述的初始化设计变量的方法为：将所有的杆件取相同的截面积

其中，在步骤五中所述的有限元分析，采用商业有限元分析软件Nastran或Ansys实现，有限元分析的结果是各个基础载荷工况f₁,…,f_2n下结构的位移向量u₁,…,u_2n；

其中，在步骤六中所述的计算结构柔度均值与标准差的方法为：结构柔度的均值为其中表示载荷f_i在位移u_j上所作之功，结构柔度的方差为 ${\mathrm{\σ}}^2\left(c\right)=\underset{i,j,k,l=1}{\overset{2n}{\Σ}}({\mathrm{\ξ}}_{ijkl}-{\mathrm{\ξ}}_{ij}{\mathrm{\ξ}}_{kl})c_{ij}{c}_{kl};$

其中，在步骤七中所述的灵敏度分析的方法为：根据结构柔度均值与方差的表达式，目标函数对于x_e的偏导数为：

$\frac{\∂J}{\∂x_e}=\underset{i,j=1}{\overset{2n}{\Σ}}({\mathrm{\α\ξ}}_{ij}+\frac{\mathrm{\β}}{\mathrm{\σ}\left(c\right)}\underset{k,l=1}{\overset{2n}{\Σ}}({\mathrm{\ξ}}_{ijkl}-{\mathrm{\ξ}}_{ij}{\mathrm{\ξ}}_{kl}\left)c_{kl}\right)\frac{\∂c_{ij}}{\∂x_e}$

记并以w_ij为元素构造矩阵W＝(w_ij)；对W进行对角化分解得到Q^TWQ＝diag{λ₁,…,λ_2n}，其中λ_i为W的特征值，其中，i＝1,…,2n，Q为由相应的单位特征向量组成的正交阵；设q_ij为正交阵Q位于i行j列的元素，构造新的载荷工况和相应的位移向量则目标函数对设计变量x_e的偏导数由下式给出：

$\frac{\∂J}{\∂x_e}=\underset{i=1}{\overset{2n}{\Σ}}{\mathrm{\λ}}_i\frac{\∂\left(F_i^T{U}_i\right)}{\∂x_e}=\underset{i=1}{\overset{2n}{\Σ}}{\mathrm{\λ}}_i{U_{ie}}^T\frac{\∂K}{\∂x_e}{U}_{ie}$

其中，U_ie为杆件e在载荷工况F_i下的单元位移向量；

约束对于设计变量的灵敏度为：

$\frac{\∂g}{\∂x_e}=l_e;$

其中，在步骤八中所述的利用渐近线方法更新设计变量的方法具体实现过程如下：首先根据步骤七中灵敏度分析获得的灵敏度信息构造目标函数的近似函数：

$f\left(x\right)\≈f\left(x^{\left(k\right)}\right)+\underset{e=1}{\overset{m}{\Σ}}{p}_e^{\left(k\right)}(\frac{1}{U_e^{\left(k\right)}-x_e}-\frac{1}{U_e^{\left(k\right)}-x_e^{\left(k\right)}})+\underset{e=1}{\overset{m}{\Σ}}{q}_e^{\left(k\right)}(\frac{1}{x_e-L_e^{\left(k\right)}}-\frac{1}{x_e^{\left(k\right)}-L_e^{\left(k\right)}})$

其中x^(k)为当前迭代点，和为移动渐近线，其迭代格式为

$L_e^{\left(k\right)}=x_e^{\left(k\right)}-0.1(x_e^{\max }-x_e^{\min })$

$U_e^{\left(k\right)}=x_e^{\left(k\right)}+0.1(x_e^{\max }-x_e^{\min })$

为设计变量x_e的当前值，和分别为x_e的上下限； $p_e^{\left(k\right)}=max({\left(U_e^{\left(k\right)}-x_e^{\left(k\right)}\right)}^2\∂f/\∂x_e,0),$ 然后利用该近似函数代替原优化问题中的目标函数而形成优化子问题；最后利用对偶方法求解该优化子问题而得到更新后的设计变量。

2.根据权利要求1所述的一种不确定载荷下的桁架结构稳健设计方法，其特征在于：在步骤一中所述的建立桁架结构稳健优化模型的依据和方法为：在预定的载荷条件下，桁架结构的柔度越小，其承受外载荷的能力越强；在载荷不确定条件下，结构的柔度会随着载荷的变化而变化，设计桁架结构时不仅要减小结构柔度的均值μ(c)，还要减小结构柔度随载荷变化而引起的波动，即减小结构柔度的标准差σ(c)；这里以桁架结构柔度均值μ(c)与标准差σ(c)的加权和J＝αμ(c)+βσ(c)最小为目标建立优化模型，数学表达简单且物理意义更加明确。

3.根据权利要求1所述的一种不确定载荷下的桁架结构稳健设计方法，其特征在于：在步骤九中所述的终止条件包括迭代次数和各个杆件连续两次迭代截面积变化量的最大值两方面；若迭代次数达到了最大次数或者所有杆件连续两次迭代截面积变化量均不超过给定值，则停止迭代；否则继续迭代过程。