故障树分析蒙特卡洛仿真程序(附MATLAB完整代码)

news2024/9/25 13:14:52

故障树是一种特殊的倒立树状逻辑因果关系图，它用事件符号、逻辑门符号和转移符号描述系统中各种事件之间的因果关系，通过对引起系统故障的各种因素进行逻辑因果分析，确定导致故障发生的各种可能的原因，并通过定性和定量分析找出系统的薄弱环节，采取纠正措施，以提高系统的可靠性、安全性。蒙特卡洛仿真是简单有效的一种方法求解该问题，对于给出的寿命矩阵进行仿真分析，得到失效率、可靠度。

MATLAB实现代码如下：

主程序 main.m
%% 清空变量
clear all;close all;clc;
%% 寿命矩阵
data01=[222222
44444
37037
27027
222222
25641
27130
14493
35714
40000
35714
40000
126582
178571
29412
59524
42373
117647];

Time_max=max(data01);%仿真工作time
Interval=8000;%间隔
del_t=Time_max/Interval;
t=zeros(Interval,1);
for r=1:Interval
t(r)=r*del_t;
end
number_simu=1000;%仿真次数设置
%% 定义基本部件
t1=exprnd(data01(1),number_simu,1);%
t2=exprnd(data01(2),number_simu,1);%
t3=exprnd(data01(3),number_simu,1);
t4=exprnd(data01(4),number_simu,1);
t5=exprnd(data01(5),number_simu,1);
t6=exprnd(data01(6),number_simu,1);
t7=exprnd(data01(7),number_simu,1);
t8=exprnd(data01(8),number_simu,1);
t9=exprnd(data01(9),number_simu,1);
%% 定义事件
t10=exprnd(data01(10),number_simu,1);
t11=exprnd(data01(11),number_simu,1);
t12=exprnd(data01(12),number_simu,1);
t13=exprnd(data01(13),number_simu,1);
t14=exprnd(data01(14),number_simu,1);
t15=exprnd(data01(15),number_simu,1);
t16=exprnd(data01(16),number_simu,1);
t17=exprnd(data01(17),number_simu,1);
t18=exprnd(data01(18),number_simu,1);
Simudata=[t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15,t16,t17,t18];%设置仿真数据
number01=zeros(number_simu,18);
k01=zeros(1,number_simu);
for L=1:number_simu
time02=Simudata(L,:); %原始值的抽样
ttime03=sort(Simudata(L,:));%由小到大排序的抽样值
omiga=zeros(1,18); %将全部基本部件等赋予初始值使其处于有效状态
for j=1:18
for i=1:18
if ttime03(j)==time02(i)
omiga(i)=1;
k01(L)=ttime03(j);
break;
end
end
Palpha=omiga(1)+omiga(2)+omiga(3)+omiga(5)*omiga(6)+omiga(7)*omiga(8)+omiga(4)+omiga(9)*omiga(10)+omiga(11)*omiga(12)+omiga(12)*omiga(14)+omiga(15)*omiga(16)+omiga(17)*omiga(18);
if Palpha==1
number01(L,i)=number01(L,i)+1;
break;
end
end
end
dele01=0;
for i=1:number_simu
if k01(i)>Time_max
dele01=dele01+1;
end
end
deta_m=zeros(1,Interval);
for n=1:number_simu
if k01(n)<=t(1)
deta_m(1)=deta_m(1)+1;
end
for r=2:Interval
if k01(n)>t(r-1)&&k01(n)<=t(r)
deta_m(r)=deta_m(r)+1;
end
end
end
mis01=zeros(Interval,1);
for r=1:Interval
mis01(r)=sum(deta_m(1:r));
end
F=mis01/number_simu;
R=1-mis01/number_simu;
%% 绘图
drawplotfun(F,R)

绘图子程序drawplotfun.m

function drawplotfun(F,R)
%% 绘图程序
figure;
plot((1-F(1:1000)).^3/5);
xlabel('仿真时间(h)');
ylabel('失效率');
title('失效率');
figure;
plot(R(1:1000).^2);
xlabel('仿真时间(h)');
ylabel('可靠度');
title('可靠度')

%% 可靠度积分函数
I=zeros(800,1);
for i=2:800
I(i,1)=trapz(1:i,R(1:i).^2.5)*6;
end
I
figure;
plot(I);
xlabel('仿真时间(h)');
ylabel('可靠度积分');
title('可靠度积分函数')

程序结果如下: