Matlab故障树的最小割集的求解

Matlab故障树的最小割集的计算

---

一、项目描述

来源于上课时教员要求我们编一个可以自动计算一颗故障树最小割集的程序，程序的输入为矩阵模式，输出为该故障树的最小割集。

二、求解思路

故障树求解最小割集有很多种方法，主要简单介绍一下上行法和下行法。上行法是从底层事件和中间事件出发，不断的由下面的公式迭代，最终得到顶事件的表达。
下行法顾名思义，从顶事件出发，将中间事件不断向下迭代，最终得到一个全部由底事件表达关于顶事件的等式。

故障树由四类有效元素构成：顶事件、中间事件、底层事件、逻辑门。其中，顶事件、中间事件、底层事件也可以归属到一起，为了方便求解，我们分开计算。

1.代码如下：

clc
clear
prompt='How many bottom events are there?\n ';
num_variables=input(prompt);
syms variables intermediates T intermedia
for i=1:num_variables
    syms (['x',num2str(i)]); 
    variables(i)=['x',num2str(i)];
end
prompt2='How many intermediate events are there?\n';
num_mediate=input(prompt2);
intermedia(1)=T;
for j=1:num_mediate
    syms (['M',num2str(j)]); 
    variables(num_variables+j)=['M',num2str(j)];
    intermedia(j+1)=['M',num2str(j)];
end
prompt3='Please input the Adjacency Matrix of bottom and intermediate events\n';
prompt4='邻接矩阵为（（中间事件个数+顶事件）*（底事件个数+中间事件个数）），本案例中即为7*14的矩阵';
disp(prompt4);
ad_mat=input(prompt3);
prompt5='Please input the relationships of all events in order\n';
disp('输入的格式为0 1 向量的格式，1 代表 +，0 代表 - 输入数量与邻接矩阵行数一样，即为 1*7 个');
disp('输入的门的顺序为顶事件、M1、M2···案例中为[1 1 0 1 1 1 1]');
relationship=input(prompt5);
for v_num=1:num_mediate+1
    [row,col]=find(ad_mat(v_num,:)==1);
    if relationship(v_num)==1
        eq_variable_nums = (1:length(row));
        intermediates(v_num) =  sum(variables(col(eq_variable_nums)));
    else
        eq_variable_nums = (1:length(row));
        intermediates(v_num) = prod(variables(col(eq_variable_nums)));
    end
end
for mm=1:v_num
    eval([char(intermedia(v_num+1-mm)),['=',char(intermediates(v_num+1-mm))]]);
end
T=expand(T);
str_eq=char(T);
S = regexp(str_eq, '+', 'split');
S=strtrim(S);

%去除底事件的次方项数(^2,^3,……，^n)
for i = 1:length(S)
    cifang=strfind(S(i),'^');
    if isempty(cifang{1})==0
        str=cell2mat(S(i));
        start=cifang{1};
        while str(start)~='x' && str(start)~='*' && start<(length(str)+1)
            start=start+1;
        end
        S(i)=erase(S(i),str(cifang{1}:start-1));
    end
end
%去除被包含的事件
min_=[];
for i =1:length(S)-1
    for j=i+1:length(S)
        temp_str=strsplit(S{j},'*');
        temp_str_2=strsplit(S{i},'*');
        # 判断temp_str_2 是否是temp_str的子集，是则返回1，否则返回0
        judge=all(ismember(temp_str_2,temp_str));
        if judge==1
            min_(end+1)=j;
        end
    end
end
%清除S中不能作为最小割集子元素的集合并输出
S(min_)=[];
fprintf('该故障树的最小割集为：\n');
for i = 1:length(S)
    disp(S(i))
end

2.读入数据

输入的中间事件和顶事件之间的关系矩阵为
      [1 1 0 0 0 0 0 0 1 0 0 0 0 0;
       0 0 0 0 0 0 0 0 0 1 1 0 0 0;
       0 0 0 0 0 0 0 0 0 0 0 1 1 0;
       0 0 1 0 0 0 0 0 0 0 0 0 0 1;
       0 0 0 1 1 0 0 0 0 0 0 0 0 0;
       0 0 0 0 0 1 1 0 0 0 0 0 0 0;
       0 0 0 0 0 1 0 1 0 0 0 0 0 0];
   
输入的门的顺序为顶事件、M1、M2···案例中为[1 1 0 1 1 1 1]

---

结果

最后，讲讲其中的一些处理要点。
最开始我犯了几个错误：
一是没有考虑到表达式中会存在平方项，所以要把最终顶事件的表达中的平方或者更高阶项去掉以后，再去求最小割集。
二是matlab相对而言更加人性化，它会把所有的顶事件基于一定的规则进行排序，这点有利于后期对一些非最小割集进行判断。
三是最小割集判断中会存在这样的情况：比如一个割集为{‘x2’,‘x6’},另一个为{‘x2’,‘x3’,‘x4’,‘x6’}最开始的代码没有考虑这样的情况，只是简单的用了某个字符串是否被包含在另一个字符串中，面对这种情况就无法判断出来。所有更好的办法是把每个割集的底事件都单独列出，判断一个cell数组是否为另一个的子集。这里由于每个割集都不会为空，所以使用ismember不用去考虑空集是任何集合的子集的情况。