可靠性和可靠性灵敏度分析的Monte-Carlo数值模拟法

东北大学机械设计及理论研究所机械系统可靠性和可靠性灵敏度分析的MonteCarlo数值模拟法授课教师：黄贤振东北大学机械设计及理论研究所由构件和运动副组成的机械系统若干元件组成的承受外部作用力并有特定功能的整体构机系统多体系统是指由多个刚性或柔性体通过一定的联接方式互相关联起来的完成预期动作的一个整体机械系统机械结构系统体多系统图1机械系统的分类机械系统是指由多个机械基本要素（如机械零、部件，机构，机器等）组成的，可以完成所需动作(或动作过程)，实现动作的传递和力的变换以及机械能的转化和利用的装置。1机械系统东北大学机械设计及理论研究所2机械可靠性应力分析应力分布参数强度分析强度分布参数载荷分布参数机械系统构件分布参数机械系统性能分布参数机械系统构件分布参数f(x)fR(xR)fS(xS)x机械可靠性模型应力是指对产品功能有影响的各种因素，比如机械应力、变形、磨损等。强度是指产品承受应力的能力，比如机械强度、刚度、抗裂度等。东北大学机械设计及理论研究所假定机械系统的强度随机变量xR，应力随机变量为xS，机械系统极限状态函数（功能函数）为(,)RSRSZgxxxx极限状态的定义为：整个机械系统的一部分超过某一特定状态就不能满足设计规定的某一功能要求，此特定状态为该功能的极限状态。0(,)00RSZgxx失效状态极限状态可靠状态2机械可靠性东北大学机械设计及理论研究所f(x)fR(xR)fS(xS)xxSdxS应力xS落在附近dxS小区间内的概率为d()dSSSSSPxxfxx在应力xS落在小区间dxS内的条件下，强度xR小于应力xS的概率为d()SxRSSSRRRPxxxxfxdx,d()d()SxRSSSSSSRRRPxxxxfxxfxdx应力xS落在小区间dxS内，同时强度xR小于应力xS的概率为根据可靠度的定义，对于应力xS所有的可能值强度xR均小于应力xS的概率，及事件(xRxS)的边缘概率就是机械系统的可靠度(对上式的求和)2机械可靠性东北大学机械设计及理论研究所(0)()()dd()()()()SxfRSSSRRRSrsRSSSSPPZfrfsrsfxfxdxdxFxfxdx(0)()()1()()RfRRSSSRSRRRRxPPZfxfxdxdxFxfxdx同理，也可以求得失效概率的另一种表达式2机械可靠性具有n设计变量的机械系统的功能函数可以表示为12()(,,...,)nZggxxxx则极限状态方程g(x1,x2,…,xn)=0将结构的基本随机变量空间分为失效域和可靠区域两部分。东北大学机械设计及理论研究所2机械可靠性12f1212()0()()()dddnXXXnngPfxfxfxxxxx式中,f(xi)(i=1,2,…,n)为随机变量xi的概率密度函数。式中，f(x1,x2,…,xn)是基本随机变量x=[x1,x2,…,xn]T的联合概率密度函数。若基本随机变量是相互独立的，则有f1212()0(,,,)dddnngPfxxxxxxXx失效概率Pf可表示为东北大学机械设计及理论研究所MonteCarlo可靠性分析方法又称随机抽样法、概率模拟法或统计试验法。该方法是通过随机模拟或者说统计试验来进行结构可靠性分析的。由于它是以概率和数理统计理论为基础的，故被无理学家以赌城MonteCarlo来命名。3MonteCarlo可靠性分析东北大学机械设计及理论研究所3MonteCarlo可靠性分析根据概率论的大数定律，若有来自同一母体且有相同分布的n个相互独立的随机样本x1,x2,…,xn,他们具有相同的有限均值μ和方差σ2，则对于任意的ε0,有11lim1niniPxn上式表明样本均值是依概率收敛于母体的均值μ的。11niixn3.1MonteCarlo模拟的理论基础东北大学机械设计及理论研究所3MonteCarlo可靠性分析另外，设随机事件A发生的概率P(A),在n次独立试验中，事件A发生的频率为m，则随机事件A发生的频率W(A)=m/n,对于任意ε0，有上式表明事件发生的频率是依概率收敛于事件发生的概率的。lim()1nmPPAn东北大学机械设计及理论研究所3MonteCarlo可靠性分析3.2MonteCarlo可靠性分析的原理与计算公式1.MonteCarlo方法求解失效概率估计值的计算公式失效概率可以改写成下式所示的失效域指示函数IF(x)的数学期望形式：f1212()01212(,,,)ddd()(,,,)dddE[()]nnngFnnFRPfxxxxxxIfxxxxxxIXxXxx式中，为失效域指示函数；Rn为n维变量空间；E[.]为数学期望算子。1,()0,FxFIxxF东北大学机械设计及理论研究所3MonteCarlo可靠性分析失效概率为失效域指示函数的数学期望，依据大数定律，失效域指示函数的数学期望可以有失效域指示函数的样本均值来近似。即随机变量的联合概率密度函数fX(x)抽取N个样本xj(j=1,2,…,N),落入失效域F内样本点的个数Nf与总样本点的个数N之比即为失效概率的估计值。f11ˆ()NfFjjNPINNx(1)东北大学机械设计及理论研究所3MonteCarlo可靠性分析2.MonteCarlo失效概率估计值的方差分析由式(1)可以看出，失效概率估计值的随机样本xj(j=1,2,…,N)的函数，因此也是一个随机变量。为了对计算出的的收敛性有一个清楚的认识，有必要对的方差进行分析。对式(1)两边求数学期望，可得失效概率估计值的期望如下所示：fˆPfˆEPf11ˆ()NFjjEPEINxff11ˆ()()NFjFjjEPEIEIPNxx样本xj与母体x独立同分布东北大学机械设计及理论研究所3MonteCarlo可靠性分析由上式可知，，即为Pf的无偏估计ffˆEPPff11ˆˆ()NFFjjEPIIPNx在数值模拟的过程中，以指示函数IF(x)的样本均值近似代替E[IF(x)],则失效概率估计值的期望可以近似表达为FIff11ˆ()()NFjFjjEPEIEIPNxx()FFEIIx东北大学机械设计及理论研究所3MonteCarlo可靠性分析失效概率估计值的方差可以通过对式(1)两边求方差得如下：f21111ˆVarVar()Var()jxNNFjFjjjPIINNxx独立由于样本方差依概率收敛于母体的方差，所以可以用IF(.)的样本方差22211()1NFjFjSININxf1ˆVarVar()FPINx(2)xj与母体x独立同分布东北大学机械设计及理论研究所3MonteCarlo可靠性分析将上式代人式(2)，可得失效概率估计值的方差估计为2f1ˆˆˆVar()1ffPPPN进而得到估计值的变异系数为fˆˆVar[]1ˆCov()ˆˆ[](1)ffffPPPEPNP222211122fff1111Var()()()()11ˆˆ()1ˆ()11NNNFFjFFjFkjjkNFjjNIINIIINNNNNPPNIPNNNxxxxx指示函数IF(x)方差的无偏估计可以进一步表达为？东北大学机械设计及理论研究所3MonteCarlo可靠性分析f11ˆ()NFjjPINxf11ˆVarVar()NFjjPINxf211ˆVarVar()NFjjPINx2221221122fff11()111()()1ˆˆ1()ˆ()11NFjFjNNFjFkjkNFjjSININNIINNNNNPPIPNNNxxxxIF(xj)方差的无偏估计f1ˆVarVar()FjPINx2f1ˆˆˆVar()1ffPPPNfˆˆVar[]1ˆCov()ˆˆ[](1)ffffPPPEPNPxj相互独立xj与母体x独立同分布开始设置抽样样本数N，初始值m=0,j=0结束否IF(xj)=0IF(xj)=1是否是j=j+1由概率密度函数fX(x)参数随机样本点xj将样本xj代人功能函数，计算功能函数值g(xj)g(xj)≤0?m=m+IF(xj)j=N?2ffffˆˆˆˆ,Var1PPmPPNN依据随机变量的分布形式和参数，由随机样本的产生方法产生N组随机样本向量的样本xj=(xj1,xj2,…,xjn)(j=1,2,…,N)。将随机向量样本xj代人极限状态方程，并根据状态指示函数IF(xj)进行累加。求得失效概率估计值fˆP估计失效概率估计值的方差及变异系数。fˆP东北大学机械设计及理论研究所3MonteCarlo可靠性分析常见分布随机数生成函数的调用格式东北大学机械设计及理论研究所3MonteCarlo可靠性分析,0.577...6XXXXXabb,0.577...6XXXXXabbⅠ型极小值分布由均值和方差求Ⅰ型极小值分布的分布参数mu--meanparametersigma--standarddeviationparametermu--locationparametersigma--scaleparameter东北大学机械设计及理论研究所,0.577...6XXXXXabb5MonteCarlo可靠性分析东北大学机械设计及理论研究所5MonteCarlo可靠性分析采用向量运算代替循环提高Matlab计算效率东北大学机械设计及理论研究所4可靠性灵敏度可靠性灵敏度定义为基本随机变量分布参数的变化引起失效概率变化的比率，在数学上可靠性灵敏度是由失效概率Pf对基本随机变量分布参数θx的偏导数予以表达，即。可靠性灵敏度反应了基本变量分布参数对失效概率的影响程度，无量纲正则化的可靠性灵敏度可以给出基本变量分布参数对可靠度的重要排序。东北大学机械设计及理论研究所5MonteCarlo可靠性灵敏度分析可靠性灵敏度为失效概率Pf对基本随机变量xi的分布参数θ(i=1,2,…,n;k=1,2,…,mi)，其中mi为第i个变量xi的分布参数的总数)的偏导数，将失效概率的积分公式对分布参数求导数，便可以得到可靠性灵敏度如下所示：()ihx()ifkxPf()(),,()dhhFXiXiPfXxxf()0()dyPfXxxx5.1MonteCarlo可靠性灵敏度估计的计算公式将可靠性灵敏度定义式做如下变换，可使可靠性灵敏度变成数学期望的形式，之后就可以采用MonteCarlo数值模拟来估计可靠性灵敏度。东北大学机械设计及理论研究所5MonteCarlo可靠性灵敏度分析()()()()()()()1()()()()()1()()()()iiiniifkkkxxxFFkkRxxPffdfdffIfIfdEffXXXXXXXXXxxxxxxxxxxxxxxf1212()01212(,,,)ddd(,,,)dddnnngFnnRPfxxxxxxfxxxxxxXxXI()()1ˆ()1()()iijNfFjkkjxjxPIfNfXXxxxxx采用样本均值代替总体均值xj是按联合概率密度函数fX(x)