CN109635452B - 一种高效的多峰随机不确定性分析方法 - Google Patents

Abstract

本发明通过结合单变量降维方法（UDRM）和最大熵方法（MEM），提出了一种高效的能用于处理多峰分布随机变量的随机不确定性分析方法。该方法将MEM由四阶矩约束扩展至n阶矩约束，然后通过进行UDRM+MRM循环来确定响应分布收敛时的响应统计矩（或MEM的矩约束）的阶数，并在此基础上利用UDRM+MEM方法求出响应的概率分布和响应各点概率。本发明能够同时处理多峰分布随机变量和多峰分布响应，从而较好地解决了雅克比矩阵G接近奇异或病态的问题，增加了方程组的可解性；同时在保证结果准确性的前提下，能够大幅度降低对样本规模的需求，仅需极少的样本数量便获得了较高的计算精度。

Description

一种高效的多峰随机不确定性分析方法

技术领域

本发明属于随机不确定性分析方法领域，涉及一种高效的多峰随机不确定性分析方法。

背景技术

工程中希望产品可靠、稳定、安全，但制造误差、材料属性的变化、使用环境的变化和认识的不完备性等各种不确定性因素会对产品品质造成不利影响，因此在设计阶段对各种不确定性进行有效度量和控制对于保证产品的质量和可靠性具有十分重要的意义。不确定性通常分为随机不确定性和认知不确定性两类，而随机不确定性又称为客观不确定性，来源于物理系统或环境固有的随机性和波动性，是一种无法随着认知水平的提高而消除的不确定性，也是本发明的研究对象。

目前常用的传统随机不确定性分析方法大致有3类：第一类是基于抽样的蒙特卡洛模拟(MCS)方法，该类方法灵活易用并且都能获得良好的计算精度，但由于需要进行大规模计算，所以这类方法常被用于验证替代方法，而不直接使用其进行不确定性分析；第二类是基于MPP点的方法，这类方法包括一阶可靠性方法(FORM)、二阶可靠性方法(SORM)和一阶鞍点逼近方法(FOSA)，该类方法有较好的计算精度和计算效率，但均要求极限状态函数的非线性程度较低；第三类是矩法，这类方法包括EUDRM+Pearson、DRM+EGLD、DRM+SPA和UDRM+MEM等，相较于基于MPP点的方法，矩法由于不需要求取响应关于变量的偏导数并且无需进行变量变换而具有更高的计算效率。但是，上述传统方法(MPP和矩法)传统上只涉及单峰分布的随机变量，而在实际工程有些不确定性参数却服从多峰分布。例如，同时拥有公路和铁路两种交通的桥梁在易疲劳位置处所受的应力服从多峰分布；发电机转子在发电机启动闭合时产生的疲劳应力服从双峰分布；用于估计交通水平的轴载荷服从双峰分布。若采用传统方法，由于其自身的局限性而较难获得较好的计算精度，也很难准确估计出响应的双峰特性。

鉴于传统方法处理多峰随机不确定性问题存在上述问题而关于多峰随机不确定性的研究又很少，因此有必要提出一种高效的用于处理多峰分布随机变量的随机不确定性分析方法。

发明内容

针对上述传统方法解决多峰分布的随机变量问题的局限性，本发明通过结合单变量降维方法(UDRM)和最大熵方法(MEM)，提出了一种高效的能用于处理多峰分布随机变量的随机不确定性分析方法。该方法将MEM由四阶矩约束扩展至n阶矩约束，然后通过进行UDRM+MRM循环来确定响应分布收敛时的响应统计矩(或MEM的矩约束)的阶数，并在此基础上利用UDRM+MEM方法求出响应的概率分布和响应各点概率。

本发明提供一种高效的多峰随机不确定性分析方法，包括如下步骤：

(1)设定响应统计矩的初始阶数l＝0和初始熵值E₀，设置允许误差ε(一个较小的正数)；

(2)令l＝l+2；

(3)引入单变量降维积分方法，将高维积分转换成低维积分来计算响应Y的前l阶原点矩m_l；利用加性分解对功能函数g(X)进行近似后，响应Y的各阶原点矩公式如下：

式中μ_i表示随机变量X_i的均值，N表示随机变量的个数，通过二项式定理展开上式可重新表达如下：

作以下定义：

则式(2)可简化成：

上式

可通过下列递归公式求出：

进而问题转化成求以下一维积分：

其中，

是输入随机变量X_j的概率密度函数，h＝j-k；

该一维积分通过基于标准矩的求积准则(NMBQR)进行计算：

为构建在X_j维坐标方向上含n个插值节点x_j,i,i＝1,2,...,n的NMBQR，定义如下函数：

其中，μ_j和σ_j分布代表X_j的均值和标准差,u_j,i表示标准插值节点，n表示积分节点的数量，且P(u_j)满足下式：

令：

则(8)式可用如下线性方程组表达：

式中μ′_j,i表示X_j的标准中心矩：

求解线性方程组(10)可得到r′_j,i，则节点u_j,i,i＝1,2,...,n即为如下多项式方程的根：

得到u_j,i后，X_j的插值节点和权重可以通过下列公式得到：

x_j,i＝σ_ju_j,i+μ_j (13)

其中

因此，式(6)的一维积分可以简化成如下代数运算：

(4)计算如下Hankel矩阵的行列式Δ_k：

如果对每一个k的取值，都满足Δ_k＞0，则继续下一步，否则直接转入第(9)步；

(5)将响应Y的原点矩m_i(i＝1,2,…,l)通过如下等式转化成非标准中心矩m′_i(i＝1,2,…,l)：

式中

是二项式系数，m₀＝∫ρ(y)dy＝1为响应Y的0阶原点矩，c为非标准中心矩m′_i给定的比例因子，为常数；

(6)假设待估响应Y的概率密度函数为ρ(y)，则其熵值为：

S(ρ)＝-∫ρ(y)logρ(y)dy (18)

使用最大熵方法(MEM)求取ρ(y)可描述为以下优化问题：

式中μ_Y是响应的均值，即μ_Y＝m₁，

使用拉格朗日乘子法求解上式，得到拉格朗日函数如下：

解上式可得ρ(y)的解析式如下

拉格朗日乘子λ_k(k＝0,1,…l)可通过求解下列非线性方程组得到：

使用标准牛顿法求解上述方程组，首先在试验值λ⁰处对G_i(λ)进行一阶泰勒展开，得到以下线性方程：

定义向量δ和v：

δ＝λ-λ⁰ (23)

v＝[m′₀-G₀(λ⁰),…,m′_n-G_n(λ⁰)]^t (24)

定义雅克比矩阵G为：

则式(22)可写成：

Gδ＝v (26)

当初值λ⁰给定时就可以求出δ，然后由λ＝λ⁰+δ不断更新λ，直至λ收敛，注意到雅可比矩阵G为对称Hankel矩阵，可以由以下2l+1个值G_i(λ)(i＝0,1,…,2l+1)来确定:

在求出拉格朗日乘子λ_k(k＝0,1,…l)并得到响应Y的概率密度函数(PDF)ρ(y)之后，可以通过在响应Y的取值范围内对ρ(y)进行积分来得到响应Y每个点的概率；

(7)使用式(18)计算待估响应Y的概率密度函数ρ(y)的熵值E_l；

(8)判断||E_l-E_l-2||/||E_l||≤ε是否成立，若成立直接进入下一步，否则重复第(2)到第(8)步；

(9)输出收敛时的响应矩的阶数l和响应Y的概率密度函数ρ(y)。

本发明的有益效果是：

1.本发明能够同时处理多峰分布随机变量和多峰分布响应。相较于传统最大熵方法(MEM)，本发明将MEM的矩约束扩展至n阶，以使其能更好地用于处理多峰分布问题，并通过使用中心矩约束代替原点矩约束和选择合适的比例因子c显著减小了G_i(λ)对不同λ_j的敏感性差异，从而较好地解决了雅克比矩阵G接近奇异或病态的问题，增加了方程组的可解性。

2.本发明对n维随机不确定问题只进行了(nN+1)次函数评价，这对于实际工程问题是可以接受的，而相较于目前广泛使用的蒙特卡洛模拟法(MCS)，本发明在保证结果准确性的前提下，能够大幅度降低对样本规模的需求，仅需极少的样本数量便获得了较高的计算精度，因此有着较高的计算效率。

3.本发明通过采用单变量降维方法对各阶响应矩进行高精度计算。本发明中提出了UDRM+MEM循环，其收敛机制可以确定最合适的矩阶来估计响应的概率密度函数，所求得的响应分布能非常好地贴合MCS求得的响应分布，因此有着较高的计算精度。

附图说明

图1是本发明的流程框图；

图2是52杆穹形桁架结构示意图；

图3是52杆穹形桁架FEM模型示意图；

图4是使用本发明方法与MCS方法求得的响应PDF结果对比；

图5是响应δ₁的PDF和熵值随l变化的演化图。

具体实施方式

本发明基于单变量降维分解法(UDRM)和最大熵原理(MEM)，提出了一种高效的能用于处理多峰分布随机变量的随机不确定性分析方法。本发明与传统的可靠性计算方法区别在于，将高维积分转换到低维空间积分，有效降低了高维高斯积分的计算难度，并能有效计算含有多峰分布随机变量的不确定性的响应统计矩问题，兼具高精度和高效率的双重优点。

下面结合附图及具体实例、采用与蒙特卡洛模拟(MCS)对比的方法对本发明作进一步详细说明：

如图2、图3所示，52杆穹形桁架，半径R＝240in，1～8杆和29～36杆的横截面积均为A₁＝2in²，9～16杆的横截面积均为A₂＝1.2in²，其余杆件的面积均为A₃＝0.6in²。桁架一共受到6个外力作用，其中在节点1处受到竖直向下的力p₁作用，在节点2和节点4受到沿穹顶内法线方向的力p₂作用，在节点3和节点5受到沿穹顶内法线方向的力p₃作用，在节点6和节点10受到沿穹顶内法线方向的力p₄作用，在节点8和节点12受到沿穹顶内法线方向的力p₅作用，在节点7、9、11和13受到沿穹顶内法线方向的力p₆作用。

假设节点1的竖直位移δ₁为兴趣值，则通过有限元分析并建立代理模型可以将节点1的竖直位移表示为：

δ₁＝g(E,P₁,P₂,P₃,P₄,P₅,P₆)

其中，弹性模量E为正态分布随机变量，所受的6个外力为多峰分布随机变量。所有变量的分布参数如表1所示。

表1随机变量P₁，P₂，P₃，P₄，P₅，P₆和E的分布参数表

如图1所示：

(1)设定响应δ₁统计矩的初始阶数l＝0和初始熵值E₀＝0.560626233634882，设置允许误差ε＝0.0006；

(2)令l＝l+2；

(3)引入单变量降维积分方法，将高维积分转换成低维积分来计算响应δ₁的前l阶原点矩m_l，利用加性分解对功能函数δ₁＝g(E,P₁,P₂,P₃,P₄,P₅,P₆)进行近似后，响应δ₁的各阶原点矩公式如下：

式中μ_j表示随机变量X_j的均值，N表示随机变量的个数，本例为N＝7，通过二项式定理展开上式可重新表达如下：

作以下定义：

则响应δ₁的各阶原点矩公式可简化成：

上式

可通过下列递归公式求出：

进而采用基于标准矩的求积准则(NMBQR)求解以下一维积分：

其中，

是输入随机变量X_j的概率密度函数，x_j,i,i＝1,2,...,n是X_j维坐标方向上的n个插值节点，w_j，i为每个插值点对应的权值，

定义如下函数：

其中，μ_j和σ_j分布代表X_j的均值和标准差,u_j,i表示标准插值节点，n表示积分节点的数量，本例取n＝12，则P(u_j)满足下式：

令：

则上式可写成如下线性方程组：

式中μ′_j,i表示X_j的标准中心矩：

求解该线性方程组可得到r′_j,i，则节点u_j,i,i＝1,2,...,12即为如下多项式方程的根：

得到u_j,i后，X_j的插值节点和权重可以通过下列公式计算：

x_j,i＝σ_ju_j,i+μ_j

其中

最终得到响应δ₁的统计矩结果及其相对MCS法计算结果的误差如下表所示：

表2响应的各阶统计矩

(4)计算如下Hankel矩阵的行列式Δ_k：

如果对每一个k的取值，都满足Δ_k＞0，则继续下一步，否则直接转入第(9)步，

(5)将响应δ₁的原点矩m_i(i＝1,2,…,l)通过如下等式转化成非标准中心矩m′_i(i＝1,2,…,l)：

式中

是二项式系数，m₀＝∫ρ(y)dy＝1为响应δ₁的0阶原点矩，c为非标准中心矩m′_i给定的比例因子，为常数；

(6)假设待估响应δ₁的概率密度函数为ρ(y)，则其熵值为：

S(ρ)＝-∫ρ(y)logρ(y)dy

使用最大熵方法(MEM)求取ρ(y)可描述为以下优化问题：

式中μ_Y是响应δ₁的均值，即μ_Y＝m₁，

使用拉格朗日乘子法求解上式，得到拉格朗日函数如下：

解上式可得ρ(y)的解析式如下

拉格朗日乘子λ_k(k＝0,1,…l)可通过求解下列非线性方程组得到：

使用标准牛顿法求解上述方程组，首先在试验值λ⁰处对G_i(λ)进行一阶泰勒展开，得到以下线性方程：

定义向量δ和v：

δ＝λ-λ⁰

v＝[m′₀-G₀(λ⁰),…,m′_n-G_n(λ⁰)]^t

定义雅克比矩阵G为：

则一阶泰勒展开后的线性方程组可写成：

Gδ＝v

当初值λ⁰给定时就可以求出δ，然后由λ＝λ⁰+δ不断更新λ，直至λ收敛，注意到雅可比矩阵G为对称Hankel矩阵，可以由以下2l+1个值G_i(λ)(i＝0,1,…,2l+1)来确定:

在求出拉格朗日乘子λ_k(k＝0,1,…l)并得到响应δ₁的概率密度函数ρ(y)之后，可以通过在响应δ₁的取值范围内对ρ(y)进行积分来得到响应δ₁每个点的概率；

(7)计算待估响应δ₁的概率密度函数ρ(y)的熵值E_l＝-∫ρ(y)logρ(y)dy；

(8)判断||E_l-E_l-2||/||E_l||≤ε是否成立，若成立直接进入下一步，否则重复第二步到第八步；

(9)输出收敛时的响应δ₁的统计矩阶数l和概率密度函数ρ(y)。

最终输出l＝12，响应δ₁的PDF如图4中实线所示，并与MCS方法所得到的响应结果进行了直观对比，可以明显发现，本发明方法得到的PDF结果能很好地与MCS结果吻合，并呈现出与MCS结果一样的双峰特性。此外，表3列出了通过计算一系列响应函数

而得到的响应δ₁的CDF结果和其相对MCS方法所得结果的误差。可以发现，本发明计算得到的相对误差在所有情况下都很小。例如，该方法最大相对误差出现在δ₁＝-1.5in时，仅为3.5088％。由此可见，本发明方法对于该算例有着很高的计算精度。

表3本发明与MCS方法所求响应CDF的结果对比

本发明之所以有如此高的精度，最重要的一点在于本发明采用的收敛机制可以通过不断增加统计矩的阶数l来逐步提高计算精度，直到满足设计精度为止。观察图5(a)可知，当响应统计矩的阶数l由2阶增加至12阶时，由本发明求得的响应δ₁的概率结果越来越接近于MCS得到的概率结果，并在l＝12时，两者达到高度吻合。而观察图5(b)容易发现，当响应统计矩的阶数l不断增加时，由本发明求得的响应δ₁的概率分布所对应的熵值不断下降并逐渐收敛，这表明响应δ₁的概率分布逐渐收敛于真实分布。

其次，使用UDRM来计算响应的统计矩也是保证本发明的计算精度的一个重要原因。表2向我们展示了分别由UDRM和MCS计算得到的响应δ₁的前12阶矩m_i＝1，2，...，12的结果以及两者的相对误差。分析可知，最大相对误差出现在m₁₂处，误差值为9.2600％，而最小相对误差仅为9.7620e-02％，出现在m₁处。值得注意的是，不同阶统计矩之间的相对误差出现了较大的差异，并容易在求解响应δ₁的PDF的过程中导致奇异或病态矩阵的产生从而影响计算精度，但本发明通过采用中心矩约束代替原点矩约束的方式可以极大程度地减小该类数值计算的难度，并保证了不确定分析的精度。

本发明也有着很高的计算效率。表4给出了本发明和MCS方法调用极限状态函数的次数对比。可以清楚地发现，由于本发明方法调用极限状态函数只发生在使用UDRM求响应统计矩的过程中，考虑到在进行UDRM+MEM循环时UDRM采用了12个积分点，所以本发明方法调用函数的次数为12×7+1＝85次，与MCS方法调用了1×10⁶次极限状态函数相比，本发明方法的计算量非常小，所以计算效率也非常高。

表4本发明与MCS方法的计算效率对比

本发明提出的处理多峰分布随机变量的随机不确定性分析方法，兼具高精度和高效率，有着较好的运用前景，值得推广。

Claims (1)

1.一种高效的多峰随机不确定性分析方法，包括如下步骤：

(1)设定响应统计矩的初始阶数l＝0和初始熵值E₀，设置允许误差ε；

(2)令l＝l+2；

(3)引入单变量降维积分方法，将高维积分转换成低维积分来计算响应Y的前l阶原点矩m_l；利用加性分解对功能函数g(X)进行近似后，响应Y的各阶原点矩公式如下：

式中μ_i表示随机变量X_i的均值，N表示随机变量的个数，通过二项式定理展开上式可重新表达如下：

作以下定义：

则式(2)可简化成：

上式

可通过下列递归公式求出：

进而问题转化成求以下一维积分：

其中，fX_j(x_j)是输入随机变量X_j的概率密度函数，h＝j-k；

该一维积分通过基于标准矩的求积准则NMBQR进行计算：

为构建在X_j维坐标方向上含n个插值节点x_j,i,i＝1,2,...,n的NMBQR，定义如下函数：

其中，μ_j和σ_j分布代表X_j的均值和标准差,u_j,i表示标准插值节点，n表示积分节点的数量，且P(u_j)满足下式：

令：

则(8)式可用如下线性方程组表达：

式中μ′_j,i表示X_j的标准中心矩：

求解线性方程组(10)可得到r′_j,i，则节点u_j,i,i＝1,2,...,n即为如下多项式方程的根：

得到u_j,i后，X_j的插值节点和权重可以通过下列公式得到：

x_j,i＝σ_ju_j,i+μ_j (13)

其中

因此，式(6)的一维积分可以简化成如下代数运算：

(4)计算如下Hankel矩阵的行列式Δ_k：

如果对每一个k的取值，都满足Δ_k＞0，则继续下一步，否则直接转入第(9)步；

(5)将响应Y的原点矩m_i，i＝1，2，...，l通过如下等式转化成非标准中心矩m'_i，i＝1，2，...，l：

式中

是二项式系数，m₀＝∫ρ(y)dy＝1为响应Y的0阶原点矩，c为非标准中心矩m′_i给定的比例因子，为常数；

(6)假设待估响应Y的概率密度函数为ρ(y)，则其熵值为：

S(ρ)＝-∫ρ(y)logρ(y)dy (18)

使用最大熵方法MEM求取ρ(y)可描述为以下优化问题：

式中μ_Y是响应的均值，即μ_Y＝m₁，

使用拉格朗日乘子法求解上式，得到拉格朗日函数如下：

解上式可得ρ(y)的解析式如下

拉格朗日乘子λ_k，k＝0，1，...，l可通过求解下列非线性方程组得到：

使用标准牛顿法求解上述方程组，首先在试验值λ⁰处对G_i(λ)进行一阶泰勒展开，得到以下线性方程：

定义向量δ和v：

δ＝λ-λ⁰ (23)

v＝[m′₀-G₀(λ⁰),…,m′_n-G_n(λ⁰)]^t (24)

定义雅克比矩阵G为：

则式(22)可写成：

Gδ＝v (26)

当初值λ⁰给定时就可以求出δ，然后由λ＝λ⁰+δ不断更新λ，直至λ收敛，注意到雅可比矩阵G为对称Hankel矩阵，可以由以下2l+1个值G_i(λ),i＝0,1,...,2l+1来确定:

在求出拉格朗日乘子λ_k，k＝0，1，...，l并得到响应Y的概率密度函数ρ(y)之后，可以通过在响应Y的取值范围内对ρ(y)进行积分来得到响应Y每个点的概率；

(7)使用式(18)计算待估响应Y的概率密度函数ρ(y)的熵值E_l；

(8)判断||E_l-E_l-2||/||E_l||≤ε是否成立，若成立直接进入下一步，否则重复第(2)到第(8)步；

(9)输出收敛时的响应矩的阶数l和响应Y的概率密度函数ρ(y)；

在步骤(1)中，定义穹形桁架的半径R、穹形桁架的杆的横截面积、穹形桁架的节点所受的外力作用，以及节点处的竖直位移，则通过有限元分析并建立代理模型可以将节点的竖直位移δ₁表示为δ₁＝g(E,P₁,P₂,P₃,P₄,P₅,P₆)，其中，弹性模量E为正态分布随机变量，所受的外力P₁，P₂，P₃，P₄，P₅，P₆作用为多峰分布随机变量。

Images