蒙特卡洛模拟.doc

是一个论文里看到的，内容如下：两元件A,B组成的串联系统,已知两元件A,B的寿命分布均为威布尔分布,则其概率密度函数和分布函数分别为:式中,m为形状参数;为尺度参数;为位置参数。由上式可得威布尔分布的直接抽样方法如下:假设元件A的参数=0小时,=5,=80小时;元件B的参数=0小时,=4,=150小时。则故障率的计算及曲线的绘制步骤为:(1)根据式(9)对元件进行抽样;(2)对系统进行50000次仿真;(3)统计寿命区间t=20小时时,每个区间的故障次数;(4)计算各个元件的故障率及系统的故障率,对于串联系统来说系统的故障率为各元件故障率的乘积;(5)绘制故障率曲线,如图2所示,其中A表示元件A的故障率,B表示元件B的故障率,S表示系统的故障率。结果是这样的：表格 1 元件和系统的故障率计算统计间隔tA的故障率aB的故障率b系统的故障率0000200.0012511e-0060.001252400.0102823e-0050.010312600.0225630.0020720.024635800.0172140.005940.0231541000.0070120.0083090.0153211200.0128450.0106960.0235411400.0176870.0096380.0273251600.0145680.0089770.0235451800.0115460.0056680.0172142000.0136980.0055150.0192132200.0156250.0058760.0215012400.0135460.0078770.0214232600.0128540.0084480.0213022800.0136250.0090070.0226323000.0140320.0088130.0228453200.0131420.0082120.0213543400.0130210.0069330.0199543600.0139850.0065550.020543800.0140320.0067530.0207854000.0137210.0073110.0210324200.0136540.0073210.0209754400.0136140.0073550.0209694600.0136320.0073350.0209674800.0136050.0073680.0209735000.0136120.0073570.0209695200.0136150.0073570.0209725400.0136140.0073490.0209635600.0136130.0073540.0209675800.0136140.0073550.020969主程序function main()n=50000; %仿真次数 %产生均匀分布随机数z1=rand(1,n); %元件A随机m1=5;eta1=80;sgma1=0;lambda1=webull(z1,m1,eta1,sgma1,n);z2=rand(1,n); %元件B随机m2=4;eta2=150;sgma2=0;figurelambda1=webull(z2,m2,eta2,sgma2,n);disp(系统的故障率)lambda=lambda1+lambda2;调用的子程序function lambda=webull(z,m,eta,sgma,n)E=eta*(-log(1-z).(1/m)+sgma;%转换成威布尔分布的随机数y=zeros(1,30);for i=1:n for t=1:30 if (E(i)20*t)&(E(i)=20*(t+1) %对于每一个随机数判断在哪个微区间内，则该区间的故障数加1 y(t)=y(t)+1; end endendy=y/n; %故障次数除以抽样次数为故障率;t=1:30;plot(t,y,.-)