蒙特卡洛模拟.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 元件和系统的故障率计算
统计间隔△t
A的故障率λa
B的故障率λb
系统的故障率λ
0
0
0
0
20
0.001251
1e-006
0.001252
40
0.010282
3e-005
0.010312
60
0.022563
0.002072
0.024635
80
0.017214
0.00594
0.023154
100
0.007012
0.008309
0.015321
120
0.012845
0.010696
0.023541
140
0.017687
0.009638
0.027325
160
0.014568
0.008977
0.023545
180
0.011546
0.005668
0.017214
200
0.013698
0.005515
0.019213
220
0.015625
0.005876
0.021501
240
0.013546
0.007877
0.021423
260
0.012854
0.008448
0.021302
280
0.013625
0.009007
0.022632
300
0.014032
0.008813
0.022845
320
0.013142
0.008212
0.021354
340
0.013021
0.006933
0.019954
360
0.013985
0.006555
0.02054
380
0.014032
0.006753
0.020785
400
0.013721
0.007311
0.021032
420
0.013654
0.007321
0.020975
440
0.013614
0.007355
0.020969
460
0.013632
0.007335
0.020967
480
0.013605
0.007368
0.020973
500
0.013612
0.007357
0.020969
520
0.013615
0.007357
0.020972
540
0.013614
0.007349
0.020963
560
0.013613
0.007354
0.020967
580
0.013614
0.007355
0.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;
figure
lambda1=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
end
end
y=y/n; %故障次数除以抽样次数为故障率;
t=1:30;
plot(t,y,'.-')
grid on
xlabel('时间t（间隔为20小时）')
ylabel('故障率')
title('故障率曲线')
legend('元件A','元件B','系统')
grid off
lambda=y;