基于通用生成函数的结构可靠性优化策略

周金宇 朱达伟 王志凌

1.金陵科技学院机电工程学院，南京，211169 2.江苏理工学院机械工程学院，常州，213001

0 引言

不确定性是客观世界的重要特征，普遍存在于工业产品的设计、制造和使用过程。工程问题中，加工尺寸、外载荷、材料等往往存在不确定性，设计时若不考虑这些不确定性，可能得到不可靠的设计结果[1]。近年来，计算技术飞跃发展，考虑不确定性的优化设计方法逐步应用于航空航天、远洋深海、高速机车、特种装备等工程领域。其中，基于概率理论和随机不确定性的可靠性设计优化(reliability-based design optimization，RBDO)是理论较成熟、应用最广泛的考虑不确定性的优化设计方法[2]。

根据优化流程的迭代结构，可将RBDO方法分为双循环法、单循环法和解耦法。双循环法是解决RBDO问题的最基本方法，具有外层设计变量优化与内层可靠性分析相互耦合的双层嵌套结构。双循环法计算成本高、效率低，因此学者相继提出了单循环法和解耦法。单循环法利用近似等效条件替代可靠性约束，避免优化循环内嵌的可靠性分析循环。为进一步提高优化求解效率，LI等[3]基于可靠设计空间，将RBDO问题完全转换为确定性优化问题；ZHOU等[4]提出的基于顺序逼近的两阶段RBDO方法显著减小了计算量。对于解耦法，DU等[5]引入偏移向量，将设计变量优化与可靠性分析解耦分离，提出了序列优化与可靠性评估(sequential optimization and reliability assessment, SORA)方法；TORII等[6]提出一种适用于各类可靠性分析方法的通用解耦方法；HAO等[7]提出一种增强步长调整(enhanced step length adjustment, ESLA)迭代算法和基于二阶可靠性分析的逐步增强顺序优化与可靠性评估方法(SSORA-SORM)。

工程结构RBDO通常面临的两大技术瓶颈是复杂功能函数的调用成本高和可靠度指标的求解精度低。一方面，实际工程计算中含复杂功能函数的问题时，往往需借助有限元分析等成本高昂的学科工具。近年来，研究人员借助多项式响应面(RSM)、人工神经网络(ANN)、支持向量机(SVM)、克里金(Kriging)模型等代理模型，通过试验样本的高效训练建立分析变量与复杂功能函数之间的简明映射，在一定程度上解决了复杂功能函数调用成本高的难题，因此功能函数调用次数不总是算法评估的核心指标。另一方面，针对工程结构RBDO中常见的变量非正态分布(如指数分布、均匀分布、多峰分布)和功能函数非线性场合，可靠度指标求解的常用方法(FORM、SORM)存在较大误差，甚至因误差不可控、迭代不收敛导致方法失效。Monte Carlo模拟(MCS)法是高精度求解方法，但计算成本极高，一般难以胜任工程结构RBDO的可靠性分析。因此，在可控成本下实施可靠度指标的高精度求解，已成为工程结构RBDO亟待攻克的现实难题。

上世纪80年代USHAKOV[8]提出通用生成函数(universal generating function, UGF)概念以来，UGF法被引入工程领域并取得丰硕成果。文献[9-11]在多状态系统研究领域中应用并进一步发展该方法，使之逐步成为系统概率分析的重要工具。LISNIANSKI[12]基于变量离散化思想，利用UGF计算连续状态系统可靠性指标的界后，UGF法逐步扩展到具有连续型随机变量的结构体系概率分析[13]。

笔者在结构RBDO中引入UGF，针对随机变量为非正态分布和功能函数非线性的RBDO问题，提出一种基于UGF的可靠性设计优化策略。通过系列响应面建立分析变量与概率指标的动态映射，消除传统双循环法的内层循环；利用UGF法替代传统矩法完成响应面动态更新过程中的可靠性分析。该方法能在可控的计算成本内显著提高随机变量/参数非正态、功能函数非线性场合的工程结构RBDO的求解精度。

1 结构可靠性分析的UGF法

1.1 连续型随机变量的UGF

对于连续型随机变量s，累积分布函数及概率密度函数分别为FS(s)和fS(s)，将s在其定义域(smin,smax)内近似均匀离散化为m个点，分别记作s1、s2、…、sm。离散点si(i=1,2,…,m)对应的概率为

(1)

式中，δ为离散步长，δ=(smax-smin)/m。

这样，可根据离散数据集{(si,pi)|i=1，2，…，m}定义连续型随机变量s的UGF：

(2)

式中，离散值si为随机变量的第i个状态值，si=smin+(i-0.5)δ;z为UGF模型中的默认字符，仅用于表示函数结构，无实际含义。

1.2 基于UGF的可靠性分析

对于涉及n维连续型随机向量S=(S1,S2,…,Sn)的工程结构，可靠性分析可利用式(1)、式(2)获得S第j个分量的UGF：

(3)

设结构的功能函数为G(S)，G(S)>0时，结构可靠，否则失效。对结构进行可靠性分析时，需对各随机变量UGF进行复合运算，获得描述总体性能分布的结构UGF：

(4)

其中，⊗G为复合算子，由复合算子根据结构物理内涵定义不同变量UGF之间的运算规则。

结构UGF包含总体性能分布信息，可用来计算各类可靠性指标[9]。设结构的性能分布用UGF表示为

(5)

则依据该UGF的概率项进行条件求和，得到结构可靠度

(6)

其中，M为结构UGF的离散状态组合总数，可考虑对性能值接近的状态组合进行同类项合并，且通常M≤mn；ψ(·)为条件求和算子；I(·)为示性函数，当其自变量大于0时取1，否则取0。

由式(6)可知，基于UGF的可靠度计算原理是，对各随机变量离散状态组合进行条件枚举。借助UGF复合算子内禀的普适性、递推性、可分性、互换性等优异特性[9]，该方法适用于任意分布的随机变量，以及任意形态的功能函数。

1.3 连续型随机变量的非均匀离散化

UGF法已在多状态系统分析中得到成功应用，因存在组合爆炸的壁垒，因此在连续状态结构系统分析中的应用受限。对于连续状态的结构系统，通过低密度离散化将连续型变量描述为UGF往往达不到预期精度，高密度离散化将在后继复合运算中消耗较高计算成本，甚至导致组合爆炸，而借助同类项合并的增效技术在较多场合下无效。因此，本文针对结构RBDO问题，对随机变量进行非均匀离散化，在减少离散状态组合数的情况下保证可靠度指标的求解精度。

非均匀离散化的基本思路是当给定的离散状态数m因计算成本原因而限定为较小值时，选择结构RBDO当前迭代步的最大概然失效点(most probable failure point, MPP)为敏感点SMPP。以该点为中心，依照等比级数对各随机变量进行敏感点密集、边缘点稀疏的离散化，使SMPP邻域内的离散点密集化，以保证极限状态敏感区的概率分析精度。随机变量Si在敏感区可实现的最小离散步长为

(7)

式中，simax、simin分别为Si的最小值和最大值；SMPP(i)为SMPP的第i个分量；q为等比级数的公比，通常取1

这里的SMPP为极限状态超曲面上距离随机空间的均值点μS最近的点，通过求解下面的数学模型可得敏感点SMPP：

(8)

2 基于UGF的RBDO方法

RBDO求解的代表性数学模型为

(9)

2.1 指标函数的响应面模型

分析RBDO数学模型不难发现，不确定性功能函数Gq(d,X,P)与设计点D密切相关，每个迭代步的可靠性分析均基于当前设计点Dk。因此，Dk和该点对应的当前可靠度指标βk之间存在如下映射关系：

Dk→βk

(10)

确定指标函数的响应面模型时，可选用不含交叉项的一次和二次响应面，这两种响应面模型结构简单、运算量小且精度满足求解要求，各自表达式分别为

(11)

(12)

式中，nd为设计变量的个数；Dk,i为设计点Dk第i个分量；a、bi、ci为待定系数。

初始响应面的构建依赖于初始设计点D0邻域的N对输入输出数据：

B0=[D0,1D0,2…D0,N]
β0=(β0,1,β0,2, …,β0,N)

式中，B0为初始迭代步的响应面输入数据矩阵；D0,i为D0邻域的第i个试验点向量；β0为初始迭代步的响应面输出数据向量；β0,i为对应输入数据的第i个可靠度指标值。

对于一次响应面模型和二次响应面模型，N=2nd+1。

首先，确定初始设计点D0与初始迭代步的响应面输入数据矩阵B0。具体做法是：通过等距抽样在设计空间中生成若干点样本，对这些样本进行约束函数评估，选定位于可行域且目标函数值较小的点为初始设计点D0。以此为中心，通过下列规则确定B0中的试验点向量：

(13)

其中，α为调节系数，可取0.05～1。θ(t)(t=1,2,…,n)为扩展向量，其分量若与确定性设计变量相对应，则取值为初始设计点到该分量定义域最近边界的距离；其分量若与不确定性设计变量相对应，则取值为相应随机变量的标准差。

(14)

式中，βT为许用可靠度指标。

2.2 指标函数响应面的动态更新

图1 基于UGF的RBDO算法流程

3 算例分析

为比较不同算法的求解精度，计算最优点处可靠度指标的相对误差：

(15)

式中，βMCS为在不同算法最优点处采用Monte Carlo模拟求得的可靠度指标。

3.1 算例1

本例为通用数值算例[14]，涉及2个确定性设计变量d1、d2和2个随机参数X1、X2，且功能函数为非线性。RBDO数学模型为

(16)

其中，X1、X2均服从正态分布，其均值μ1、μ2分别为5和3，标准差σ1、σ2分别为1.5和0.9；许用可靠度指标βT=2.32；随机参数的离散状态数取30。

根据式(11)构建迭代开始时的响应面和迭代收敛时的响应面，如图2所示。不同算法的目标函数值迭代过程如图3所示，优化结果见表1，其中，N为功能函数调用次数；算法1为基于FORM的双循环法；算法2为SORA法；算法3为本文所提方法，其可靠性分析环节由UGF法完成；算法4为本文所提方法的变形，其可靠性分析环节由MCS法完成。

(a)迭代开始时响应面

图3 不同算法目标函数迭代过程

由表1可以看出，所提算法与传统算法均可收敛到最优解。算法3的精度高于传统方法(算法1和算法2)，算法4的精度更高，但对功能函数的调用过多(5×106次)。因此，与传统方法相比，算法3在计算成本可控的条件下获得了更高精度，算法可行。值得一提的是，本例中的随机变量均为正态变量，所以传统算法也能有效实施并获得一定精度。

表1 不同RBDO算法的算例1优化结果

3.2 算例2

本算例涉及设计随机变量X3和X4，其中，X3服从标准差为0.6的正态分布，X4服从标准差为0.6的均匀分布，两设计随机变量的均值μ3和μ4为不确定性设计变量。RBDO数学模型为

(17)

其中，许用可靠度指标βT=2。设计随机变量的离散状态数取30，不同算法的优化结果见表2。

由表2数据可以发现，优化求解问题含高度非正态(如均匀分布)的随机变量时，若采用矩法完成可靠性分析，则当量正态化处理极易导致显著误差[15]，因此算法1、算法2均无法收敛到最优解，方法失效。算法3和算法4均可收敛到最优点，且精度较高。兼顾求解效率，算法3最优，解决了传统算法对高度非正态随机变量无法求解的普遍难题。

表2 算例2不同RBDO算法优化结果

3.3 算例3

短梁结构如图4所示。矩形截面的高为h、宽为b，自由端受到双轴弯矩M1、M2和轴向力F作用，材料的屈服强度为Y。

图4 短梁结构示意图

确定性设计变量d=(b,h)，随机参数P=(M1,M2,F,Y)。各随机参数相互独立，其中，M1、M2、Y服从正态分布，F服从均匀分布，各自统计信息见表3。该结构RBDO的数学模型为

表3 随机参数统计信息

(18)

其中，许用可靠度指标βT=3，随机参数的离散状态数取10，不同算法的优化结果见表4。

分析表4数据可知，RBDO问题中的随机变量增多且存在非正态(如均匀分布)随机变量时，算法1可以收敛，但在最优解处的可靠度不满足概率约束的要求，求解精度极低；算法2无法收敛，方法失效；算法3以较高精度收敛到最优解。算法3为抑制多变量UGF复合运算引发的组合爆炸，将随机变量的离散状态数减为10，但由于算法采用了非均匀离散化技术，故仍可保证较高的求解精度。

表4 不同RBDO算法的算例3优化结果

3.4 算例4

本算例对单级齿轮传动减速箱的局部结构进行可靠性优化设计，要求在满足齿轮强度可靠度R≥0.99的条件下结构总质量最小。已知传动比为3.2，指定齿宽系数d1、模数d2、主动轮齿数d3、主动轴的最小直径d4、从动轴的最小直径d5为确定性设计变量。齿轮传动的输入转矩T服从均值为108.8 N·m、标准差为8.704 N·m的均匀分布；齿轮的弯曲疲劳强度极限σF,lim服从均值380 MPa、标准差30.4 MPa的对数正态分布；齿轮的接触疲劳强度极限σH,lim服从均值500 MPa、标准差44 MPa的对数正态分布。该优化问题的数学模型为

(19)

其中，确定性设计变量d=(d1,d2,d3,d4,d5)，随机参数P=(T,σF,lim,σH,lim)，许用可靠度指标βT=2.33，随机参数的离散状态数取20，不同算法的优化结果见表5。

表5 不同RBDO算法的算例4优化结果

由表5所示的优化结果可知，由于随机变量为均匀分布和对数正态分布，且同时存在2个概率约束，因此算法1、算法2均无法得到收敛解而失效，算法3可较好地收敛到最优解。与算法4相比，算法3的求解精度满足工程需求并具有更高的求解效率。

4 结论

(1)针对通常的结构RBDO，本文方法比传统方法具有更高的求解精度和可控的计算效率。

(2)RBDO问题涉及非正态随机变量、高度非线性功能函数时，传统方法存在求解精度低或无法收敛的劣势；本文方法的可靠性分析由UGF法完成，不受随机变量非正态、功能函数高度非线性的影响，具有较高的鲁棒性。

(3)RBDO问题存在高维随机空间时，所提方法可借助非均匀离散化技术，在保持随机变量离散状态数不变的条件下，提高了优化求解的精度。