基于MATLAB的蒙特卡洛方法对可靠度的计算.docx
《基于MATLAB的蒙特卡洛方法对可靠度的计算.docx》由会员分享,可在线阅读,更多相关《基于MATLAB的蒙特卡洛方法对可靠度的计算.docx(10页珍藏版)》请在冰豆网上搜索。
基于MATLAB的蒙特卡洛方法对可靠度的计算
基于MATLAB的蒙特卡洛方法对可靠度的计算
LT
摘要
对于简单的概率计算,我们可以用离散或者连续的概率分布模型进行求解;但是对于复杂的模型的近似解的求解,蒙特卡洛方法是一种非常方便的方法。
蒙特卡洛方法将最复杂的计算部分交给了电机计算机来完成,极大的方便了我们的求解过程。
本文主要是用MATLAB编写蒙特卡洛的模拟程序,然后分别验证两个正态分布的模型和两个非正态分布的模型。
非正态分布的模型中的随机变量序列都是独立同分布的,这样我们可以方便的用列维-林德伯格中心极限定理进行处理。
【关键字】:
复杂模型、蒙特卡洛、MATLAB、正太分布、独立同分布的非正态模型、列维-林德伯格中心极限定理
绪论
计算机技术的发展,促进了蒙特卡洛方法的推广、普及以及完善等。
蒙特卡洛方法诞生之初是不被重视的,因为当时的计算机技术没有达到与之匹配的程度。
蒙特卡洛模拟也称为随机模拟方法,或随机抽样技术。
它是一种以概率论和数理统计为基础,通过对随机变量的统计实验、随机模拟来求解问题近似解的数值方法。
它的主要思想是:
为了求解数学、物理、化学及工程问题,建立一个概率模型或随机过程,使它的参数等于问解;然后通过对模型或过程的观察或抽样来计算所求参数的统计特征(如均值、概率等),作为待解问题的数值解,最后给出所求解的近似值,而解的精度可用估计值的方差来表示。
蒙卡洛模拟的步骤是:
首先建立简单而又便于实现的概率分布模型,使分布模型的某些特征(如模型的概率分布或数学期望)恰好是所求问题的解;然后根据概率分布模型的特点和计算的需要改进模型,以便减少方差,降低费用,提高计算效率;再对分布模型进行随机模拟,其中包括建立产生伪随机数的方法和建立对所遇到的分布产生随机变量样本的随机抽样方法;最后建立各种统计量的估计,获得所求解的统计估计值及其方差。
蒙特卡洛模拟方法可分为直接蒙特卡洛模拟、间接蒙特卡洛模拟和蒙特卡洛积分。
(1)直接蒙特卡洛模拟采用随机数来模拟本身具有复杂随机过程的效应。
该方法是按照实际问题所遵循的概率统计规律,用计算机进行直接的抽样,然后计算其统计参数。
直接蒙卡洛模拟法能充分体现蒙特卡洛方法的特殊性和优越性,因而在物理中得到了广泛的应用,该方法也就是通常所说的“计算机实验”。
(2)间接蒙特卡洛模拟是人为地构造出一个合适的概率模型,依照该模型进行大量的统计实验,使它的某些统计参数恰好是待求问题的解。
Buffon投针实验就是运用间接蒙特卡洛模拟来求解π。
(3)蒙特卡洛积分是利用随机数系列计算积分的方法,积分维数越高,效率越高。
定积分的计算是蒙特卡洛方法被引入计算数学的开端,这里以定积分的计算说明其处理确定性问题的方法。
如计算定积分:
此时,求定积分亦即求边长为1的正方形中一个曲边梯形的面积问题,如图2所示。
可以随机地向正方形内投点,然后统计落在曲线下的点数
,当总的投点
充分大时,
就近似等于积分值s。
一、编写MonteCarlo模拟程序
1.模型的建立
本章节根据抛掷骰子编制MonteCarlo模拟程序,验证各点出现的概率均为1/6。
2.模拟流程图绘制
初始化
i=i+1
K=?
K=1
K=2
K=3
K=4
K=5
K=6
K1+1
K2+1
K3+1
K4+1
K5+1
K6+1
i<1000000
P1
P2
P3
P4
P5
P6
Y
N
图1.1流程图
3.MonteCarlo程序编写
MonteCarlo模拟程序(Matlab)
clear
N=1000000;
K_1=0;
K_2=0;
K_3=0;
K_4=0;
K_5=0;
K_6=0;
K=randi(6,N,1);
fori=1:
N
ifK(i,1)==1
K_1=K_1+1;
end
ifK(i,1)==2
K_2=K_2+1;
end
ifK(i,1)==3
K_3=K_3+1;
end
ifK(i,1)==4
K_4=K_4+1;
end
ifK(i,1)==5
K_5=K_5+1;
end
ifK(i,1)==6
K_6=K_6+1;
end
end
P_1=K_1/N
P_2=K_2/N
P_3=K_3/N
P_4=K_4/N
P_5=K_5/N
P_6=K_6/N
hist(K,6)
4.模拟结果及结论
MonteCarlo模拟得到,P_1=16.639%;P_2=16.605%;P_3=16.712%;P_4=16.710%;P_5=16.625%;P_6=16.710%。
各项约为总数的1/6,符合理论情况。
通过模拟可以得到分布直方图(图1.2)。
图1.2分布直方图
二、关于两个服从正态分布的可靠性验证
机械结构的可靠性设计中的应力-强度干涉理论的理论计算和采用蒙的卡罗方法对其进行验证。
MATLAB自带有产生正态分布的随机数,所以我们用MATLAB对N=100000实验次数进行验证。
计算次数为3次。
理论计算:
首先根据可靠度R=0.999,可得可靠度系数Z=
=3.191,然后我们确定应力XL(正态分布)的参数,均值
=200,方差
=5776,;然后再确定强度XS(正太分布)的参数,均值
=500,方差
=3062.71。
流程图的绘制
初始化
Z=S-L
i=i+1
Z<0
n=n+1
iR=n/N
Y
N
Y
N
图2.1流程图
Matlab模拟:
由理论计算的正态分布的参数进行matlab的模拟,得出的可靠度如下图:
程序如下:
N=100000;
P=[0,0,0];
R=[0.0,0.0,0.0];
forj=1:
3
S=normrnd(500,55.34,N,1);%N(500,55.34)
L=normrnd(200,76,N,1);%N(200,76)
fori=1:
N
z=S(i,1)-L(i,1);
ifz>0
P(j)=P(j)+1;
end
end
R(j)=P(j)/N
end
三、非正态分布的验证
对于非正态分布的强度-应力随机变量的可靠度计算,我们再MATLAB上用蒙的卡罗方法来验证。
验证时我们取样本值n=100000,分别验证强度服从期望为10(及λ=
)指数分布(x<0时,概率密度为0)和应力服从期望为5(及λ=
)指数分布(x<0时,概率密度为0)。
所得的可靠度如下图:
程序如下:
n=100000;
p=[0,0,0];
R=[0.0,0.0,0.0];
forj=1:
3
r1=exprnd(10,n,1);%E(0.1)
r=exprnd(5,n,1);%E(0.2)
fori=1:
n
z=r1(i,1)-r(i,1);
ifz>0
p(j)=p(j)+1;
end
end
R(j)=p(j)/n;
end
四、总结
根据强度-应力干涉模型求解系统的可靠度,对于强度和应力都服从正态分布的干涉模型,查表计算法和蒙的卡罗方法都是正确有效的。
但是当强度或者应力不再服从正态分布时,用蒙的卡罗方法可以干涉模型近似结果。
参考文献
[1]朱陆陆.蒙特卡洛方法及应用[D].华中师范大学,2014.
[2]刘岚岚.可靠性基础[D].中国质检出版社,2014.
[3]张宇.概率论与数理统计9讲[J].北京理工大学出版社.