采用增广乘子法和模拟退火法的结构可靠性分析

2019-07-11高翔王林军杜义贤

西安交通大学学报 2019年7期

高翔,王林军,杜义贤

(1.三峡大学水电机械设备设计与维护湖北省重点实验室,443002,湖北宜昌;2.三峡大学机械与动力学院,443002,湖北宜昌)

在工程实际中,由于受材料属性、制造精度、安装误差等各种不确定因素的影响,结构的实际参数会偏离设计参数。为了在设计阶段有效度量及控制不确定因素对结构的影响,需要通过可靠性分析来预先评估结构的失效风险。

关于结构的可靠性分析,国内外很多学者都进行了积极地探索和研究。Dempster和Shafer提出了证据理论,并且广泛应用于非概率结构可靠性分析及可靠性优化设计[1-2]。Zadeh等对基于元模型的多目标多学科优化结构设计进行了研究,提出了求解多学科优化设计(multi-disciplinary design optimization,MDO)问题的EMOPSO法,该方法类似于MOPSO法,引入了SQP法和元模型,从而做出了位于Pareto解处的模糊逻辑决策[3]。Wang等对基于非概率集合理论的多学科可靠性优化设计问题进行了研究,提出了一种基于MDO策略的单循环算法,并且指出MDO策略主要包括MDF法、CSSO法、CO法、BLISS法[4]。Hao等对求解非概率可靠性优化设计问题的混沌控制方法进行了研究,提出了一种初期使用HL-RF法、振荡时使用ECC法进行迭代的方法,并使用混沌动力学理论加速可靠指标的收敛速度[5]。为了提高ECC法的收敛性,在内层嵌套了Wolfe-Powell准则,用来检查以及更新控制因子[5]。Hao等对桁架结构的可靠性优化设计进行了研究,提出了一种高效的IGA法,并且通过ESLA法和SSORASORM法提高了IGA法的计算效率[6]。Moon等对实验数据不足情况下的基于信度的可靠性优化设计进行了研究,通过贝叶斯模型的核密度估计法计算概率密度函数,并由基于马尔可夫链蒙特卡罗取样器的MH算法计算后验分布,最终利用Diracδ测度计算失效概率;由于Hellinger相似性的灵敏度不能用有限差分法计算,所以引入了复变量法;此外,指出DKG法是最准确的代理模型方法之一[7]。杜秀云等对基于Bregman距离函数的可靠性分析进行了研究,引入Bregman距离函数来计算可靠指标β,根据不动点同伦映射的思想构造同伦方程组进行计算,对极限状态函数为非线性方程的可靠指标计算取得了良好的效果,尤其是极限状态函数为指数函数的情况[8]。吉猛等对基于同伦分析的结构可靠性功能度量法进行了研究,使用KKT条件构造同伦方程组,并采用β-cone搜索方法追踪到最可能点(most probable point,MPP),且使用优于可靠指标法的功能度量法来分析可靠性[9]。李彬等研究了基于改进自适应混沌控制的逆可靠性分析方法,运用自适应混沌控制方法进行了逆可靠性分析[10]。黄晓旭等对基于主动学习Kriging模型和子集模拟的可靠性分析进行了研究,提出了AK-SS法,对具有隐式功能函数的小失效概率计算取得了良好的效果[11]。孟增等对基于修正混沌控制的一次二阶矩(first order second moment,FOSM)法进行了研究,运用混沌理论对迭代震荡的情况进行修正,解决了极限状态函数非线性程度较高时的迭代振荡问题[12]。

尽管FOSM法在结构可靠性分析方面的理论已经很成熟,然而在实际工程的结构可靠性分析中却很少使用,原因如下:①如果功能函数比较复杂,将会导致雅可比矩阵的计算变得烦琐,从而整个过程效率低下;②如果雅可比矩阵的条件数过大,将会导致结果误差较大;③通过极限状态方程求解雅可比矩阵的过程中涉及到单位换算,极易得到错误的雅可比矩阵,这将导致可靠指标及失效概率计算错误。因此,如何避免雅可比矩阵的求解成了一个亟待解决的问题。

将有约束模型转换为无约束模型的常用方法为罚函数法。内点罚函数法的罚因子逐步减小,常用于不等式约束;外点罚函数法的罚因子逐步增大,常用于等式约束。由于罚函数法是一种序列无约束极小化方法,所以收敛较慢,且两种罚函数法的收敛性均依赖于罚因子的初值[13]。如果初始罚因子选取不当,就可能会导致两种罚函数法中构造的目标函数均不收敛。因此,如何避免罚函数法初始罚因子的选取成了一个亟待解决的问题。

针对以上两个技术难点,本文提出了一种采用增广乘子法与模拟退火法的结构可靠性分析方法,通过增广乘子法将等式约束优化模型转换为无约束优化模型,采用模拟退火法求解结构的可靠指标。此外,本文研究了各参数的不确定度对可靠指标的影响,并与FOSM法和蒙特卡罗模拟(Monte Carlo simulation,MCS)法进行了对比,数值算例和悬臂梁算例的结果表明,本文方法迭代次数比MCS法少,且所得结果比FOSM法更接近于MCS法。

1 可靠性理论

正态分布的概率密度函数f(x)和累积分布函数F(x)的公式[14]分别为

(1)

(2)

式中μx和σx分别为变量x的均值和标准差。

拉科维茨·菲斯莱法可以实现正态分布到标准正态分布的转换,公式[14]为

y=(x-μx)/σx

(3)

μy=x-φ-1[F(x)]σx

(4)

σy=φ{φ-1[F(x)]}/f(x)

(5)

(6)

(7)

式中:y为转换后的x;μy和σy分别为变量y的均值和标准差;φ(y)为y的概率密度函数;φ(y)为y的累积分布函数;t为积分变量。

误差函数erf(y)和互补误差函数erfc(y)的公式[15]分别为

(8)

(9)

引入误差函数后,累积分布函数φ(y)变为

(10)

可靠指标β和失效概率Pf的公式[14]分别为

β=(x-μx)/σx=μy/σy

(11)

(12)

变异系数V的公式[14]为

V=σx/μx

(13)

由大数定理中的Bernoulli定理可知,进行N次模拟之后,若功能函数g(x)<0的次数为nf,则直接抽样蒙特卡罗法的失效概率Pf[14]为

(14)

式中:xu为第u次模拟时的变量;I[x]为示性函数,公式为

(15)

(16)

为提高直接抽样蒙特卡罗法的精度,将抽样中心改为MPP。重要抽样蒙特卡罗法的失效概率Pf[14]为

(17)

式中fMPP(x)表示将概率密度函数f(x)中变量x的均值替换为MPP后的概率密度函数。

标准化正态空间中坐标原点到极限状态面的最短距离就是可靠指标β[14],此时对应的极限状态面上的点就是MPP[14]。

(18)

2 优化算法理论

2.1 等式约束的增广乘子法理论

根据可靠性理论,可建立优化模型为

subject togj(X)=0,j=1,2,…,m

(19)

式中:M为适应度;X=[X1,X2,…,Xn];μv和σv分别为Xv的均值和标准差。

式(19)为具有m个等式约束gj(X)的模型,可以由增广乘子法转化为无约束优化模型,公式[13]为

(20)

式中:Mλ为无约束优化模型的适应度;r为外罚函数法的罚因子;λ=[λ1,λ2,…,λm],为拉格朗日乘子。式(20)右端第二项为惩罚项,第三项为乘子项。

使用增广乘子法时,并不要求罚因子趋于无穷大,只需取一个比较大的值或按照一定的比例递增[14]。此方法同时应用于外点罚函数法及拉格朗日乘子法,避开了外罚函数法的初始罚因子选取。

2.2 模拟退火法理论

1983年,IBM公司的Kirkpatrick等提出了一种基于物理学正则系统的模拟退火法,该算法使用了辅助分布和多重马尔科夫链,且类似于吉布斯采样器[16-18]。吉布斯分布又称玻尔兹曼分布,表达式为

(21)

式中:i为迭代次数;P为概率;Z为配分函数或正则化常数;U(i)为势能;k为玻尔兹曼常数;T为温度;E(i)为系统的能量。

Creutz在研究物理学的伊辛模型时,提出了基于微正则系统的微正则退火法,其配分函数[18-19]为

(22)

式中:E0为初始能量值;ED为热系统中具有能量交换能力的“妖”(Demon)的能量,更新规则为

(23)

式(21)中的温度T有很多种计算方法。经过查阅相关资料,收集了以下3种常见的第i次迭代时的温度Ti的计算方法。

(1)文献[17]指出,若目标函数值的标准差为σf(x)、接受概率p>3σf(x),则温度Ti为

(24)

(2)文献[17]指出,为保证接受新解的概率大于设定值a0,应设定温度Ti为

(25)

式中:Δ+为目标函数值上升的平均值;m1、m2分别为先前实验中使目标函数下降、上升的解的数量。

(3)文献[20]指出,设定一个足够大的常数γ,使γ等于或者大于函数图形的深度,则温度Ti为

(26)

分析上述方法可知:方法(1)需要统计目标函数值的标准差,比较烦琐;方法(2)需要设定初始概率值a0;方法(3)需要知道函数图像的深度。由于以上3种方法的参数设置均比较困难,所以本文没有采用。

模拟退火法利用μ-1原理[21]计算迭代步长,公式为

(27)

式中:gμ-1=[gμ-1,1,gμ-1,2,…,gμ-1,n]为计算迭代步长的中间变量;μ0=10100η,η=(i/imax)q,q为退火因子,必须大于0,q越大则退火速度越快[21];变量Xrand=[Xrand,1,Xrand,2,…,Xrand,n]与变量X的元素数量相同,且所有元素均是位于[-1,1]内的随机数;sgn为符号函数。

模拟退火法利用Metropolis准则[13]来判断是否接受新解,Metropolis准则为:当适应度的变化ΔMλ<0或随机数p满足一定的条件时,接受新解。p需要满足的具体条件为

(28)

式中:ε是算法的精度;ε0是极小数常量,分母加上ε0是为了防止分母为零。在MATLAB软件中,p由RAND函数生成,ε0由eps函数生成。

为避免式(28)中的指数函数计算,文献[22]提出了Demon算法,具体如下:若ΔE≤D,则接受新解,同时更新Demon值D,即D=D-ΔE。但是,文献[22]提出的Demon算法的初始Demon值如何设置,文献[17]中并未记载,故本文未采用。

采用模拟退火法求解无约束优化问题的步骤如下。

步骤1 初始化X,并求出函数Mλ的函数值。

步骤2 开始迭代,设i为迭代次数。

步骤3 对每个Xν,计算迭代步长dXν=gμ-1,ν·(Bu,ν-Bl,ν),式中Bu,ν和Bl,ν分别为变量Xν的上下界。

步骤4 通过迭代步长求得新解X′=X+dX,如果X′不在范围内则随机赋予新值。

步骤5 首先求函数的变化值ΔMλ,然后根据Metropolis准则选择是否接受新解。

步骤6 保留适应度最小的解。

步骤7 判断是否达到循环终止条件。当i小于迭代上限imax时,返回步骤2;当i达到imax时,循环终止。

文献[23]指出,为了在后期能够跳出局部最优解,有学者提出了在后期提高温度(即回火)的回火退火法,也有学者提出反复执行退火降温和回火升温的计算方法。若采用文献[23]的回火退火法,则式(28)变为

(29)

本文将文献[13]的增广乘子法与文献[13,21]的模拟退火法和文献[23]的回火退火法结合,并将其应用于可靠指标的计算,分析结构的可靠性,本文方法的流程图如图1所示。

图1 本文算法的流程图

3 数值算例分析

设变量X=[X1,X2]服从正态分布,X1和X2均值分别为10和2.5,标准差分别为2和0.375,某个结构的功能函数g(X)为

(30)

设算法的终止条件为‖X(i)-X0‖/‖X0‖≤1×10-6,X0为最优解,X(i)为第i次迭代时的解。分别采用设计点法、JC法和简化加权分位值法3种FOSM法[14]求解此问题,所得结果如下。

(1)3种方法的迭代次数均为12次,且功能函数g均等于8.689 9×10-12。

(2)3种方法所得可靠性参数几乎相同,MPP约等于(11.185 490,1.654 912),可靠指标β≈2.330 217,失效概率Pf≈0.009 897。

(3)设计点法和JC法所得失效概率相同,而简化加权分位值法的失效概率略小。文献[14]指出,简化加权分位值法的精度低于JC法,但是此次的计算结果表明,简化加权分位值法与JC法的失效概率误差仅为1×10-17。

取退火因子q=1,M0为当前函数值,Mv为最优函数值,设|M0-Mv|<1×10-6或迭代次数小于等于12作为终止条件。本文方法迭代4次得:MPP为(11.746 19,1.683 011),g=0.029 489,β=2.348 96,Pf=0.009 413。

设计点法、JC法和简化加权分位值法3种FOSM法所得的可靠指标β=2.330 217,本文所得可靠指标β=2.348 96。两种方法的可靠指标几乎相同,证明本文方法可行。

经过1×107次模拟,直接抽样蒙特卡罗法求得参数如下:β=2.346 2,Pf=0.009 5,σPf=3.064 7×10-5。

设重要抽样蒙特卡罗法的抽样中心为设计点(11.185 5,1.654 9),经过1×107次模拟,重要抽样蒙特卡罗法求得参数如下:β=2.344 3,Pf=0.009 5,σPf=4.942 0×10-6。由此可见,与直接抽样蒙特卡罗法相比,在模拟次数相同时,重要抽样蒙特卡罗法的失效概率标准差σPf明显偏小,说明重要抽样蒙特卡罗法的改进效果较好。

将FOSM法、MCS法和本文方法的计算量和计算结果进行对比,结果见表1和表2,其中,FOSM法取JC法的数据,MCS法取重要抽样蒙特卡罗法的数据。

表1 数值算例下FOSM法、MCS法和本文方法

表2 数值算例下FOSM法、MCS法和本文方法

由表1和表2可以看出:本文方法的迭代次数少于MCS法,计算结果比FOSM法更接近MCS法。由此可得出,本文方法比MCS法效率更高,比FOSM法更精确。

在不同确定度下，X1、X2的均值μX1、μX2以及标准差σX1、σX2会有不同的变化范围，此时通过本文方法计算得到的对应可靠指标β也会有不同的变化范围和不确定度。不同不确定度下的μX1、μX2、σX1、σX2、β的变化范围以及β的不确定度如表3和表4所示。

表3 不同均值不确定度下μX1、μX2、β的变化范围及β的不确定度

表4 不同标准差不确定度下σX1、σX2、β的变化范围及β的不确定度

由表3和表4可知均值和标准差的不确定度与β的不确定度之间的关系,如图2所示。

(a)均值的不确定度与β的不确定度的关系

(b)标准差的不确定度与β的不确定度的关系图2 均值和标准差的不确定度与β的不确定度的关系

对表3、表4、图2进行分析,可得如下结论:

(1)β与均值呈正比关系,与标准差呈反比关系;

(2)β的不确定度与均值、标准差的不确定度分别呈现线性关系;

(3)均值的不确定度为8%时,β的不确定度为43.34%,标准差的不确定度为20%时,β的不确定度为40.07%,证明均值对β的影响比标准差对β的影响大一倍。

4 悬臂梁的可靠性分析

悬臂梁[24]长度为L,矩形横截面的宽和高分别为b和h。悬臂梁末端承受的水平载荷和竖直载荷分别为Ph、Pv,末端许可挠度[ω]=3 mm。悬臂梁的材料为45号钢,弹性模量E=210 GPa,屈服极限σs=350 MPa。此时,变量X=[b,h,L,Ph,Pv]。

设宽度b/mm～N(100,52),高度h/mm～N(200,52),Ph/kN∈[45,75],Pv/kN∈[22,28]。长度L的均值为1 000 mm,变异系数V=0.01,由式(13)可知L/mm～N(1 000,102)。根据正态分布的3σ准则可知,Ph/kN～N(60,52),Pv/kN～N(25,12)。

4.1 考虑挠度失效的可靠性分析

矩形截面的梁弯曲变形时,横截面对中性轴的惯性矩IZ为

(31)

当悬臂梁的一端承受集中力Pv时,最大挠度ωmax为

(32)

所以,当悬臂梁的一端承受Ph、Pv两个集中载荷时,产生的最大挠度为

(33)

考虑挠度失效时,功能函数g1为

(34)

由于此功能函数比较复杂,求导烦琐,雅可比矩阵较难获取,此处不再与FOSM法进行对比。

经过50次迭代,本文方法求得参数如下:MPP为(111.27,202.32,994.48,51.34,25.01),g1=-0.488 3,β=3.303 6,Pf=4.77×10-4。求解过程中β和Pf的变化如图3所示。

图3 挠度失效求解过程中β和Pf的变化

由图3可知,本文方法在前10次迭代中收敛较快,而在后40次迭代中收敛较慢。

4.2 考虑应力失效的可靠性分析

矩形截面的梁弯曲变形时,梁的抗弯截面系数WZ为

(35)

式中ymax为悬臂梁中性层到横截面两端的最大距离。

当悬臂梁一端承受Ph、Pv两个集中载荷时,矩形截面上产生的最大应力σmax为

(36)

考虑应力失效时,功能函数g2为

(37)

经过1 000次迭代,JC法求得参数如下:MPP为(84.97,195.93,100 2.76,60.0,25.0),g2=48.756 6,β=3.126 1,Pf=8.857 2×10-4。

经过500次迭代,本文方法求得参数如下:MPP为(85.4,192.9,100 4.0,70.7,25.0),g2=-0.024 0,β=3.904 1,Pf=4.728 8×10-5。求解过程中β和Pf的变化如图4所示。

图4 应力失效求解过程中β和Pf的变化

由图4可知,本文方法在前70次迭代中收敛较快,而在后430次迭代中收敛较慢。

采用式(28)的模拟退火法及式(29)的回火退火法计算可靠指标,分别运行6次后得到两组结果,分别为:3.896 6、3.920 5、3.983 8、3.911 4、3.931 4、3.949 6,以及3.914 0、3.988 3、3.916 4、3.955 7、3.952 9、3.889 2。可靠指标的统计参数见表5。

表5 可靠指标的统计参数

分析表5可知,在相同的迭代次数下,回火退火法的算术平均值、标准差和极差都较大,回火退火法的改进效果较差。这是因为回火退火法虽然提高了全局搜索能力,但是牺牲了局部搜索能力。如果将其应用于强欺骗性、多模态、多漏斗的函数,情况可能相反。

经过1×107次模拟,直接抽样蒙特卡罗法求得参数如下:β=3.828 0,Pf=6.460 0×10-5,σPf=2.541 6×10-6。

设重要抽样蒙特卡罗法的抽样中心为设计点(84.97,195.93,1 002.7,60,25),经过1×107次模拟,重要抽样蒙特卡罗法求得参数如下:β=3.835 4,Pf=6.267 7×10-5,σPf=4.120 8×10-8。与直接抽样蒙特卡罗法相比,在模拟次数相同时,重要抽样法的失效概率标准差σPf明显偏小,说明重要抽样蒙特卡罗法的改进效果较好。

FOSM法、MCS法和本文方法的计算量见表6,计算结果见表7。

分析表6和表7可知,本文方法的迭代次数比MCS法更少,所以本文方法比MCS法效率更高。FOSM法的失效概率与MCS法相比相差一个数量级,误差非常大,而本文算法的失效概率4.8×10-5更接近MCS法的失效概率6.3×10-5,说明本文方法比FOSM法更精确。由于本文方法和MCS法都不需要计算雅可比矩阵,所以FOSM法的结果误差较大极有可能是由于雅可比矩阵条件数过大导致的。