R语言实现故障树定量与定性分析——以GJB-Z 768A-1998 故障树分析指南图5.37为例

近期参与的软件开发项目，笔者直接负责FTA部分，考虑使用现有的工具包来辅助完成功能的开发，于是选用了近些年比较火热的解释性语言R来实现。

参考《GJB-Z 768A-1998 故障树分析指南》，图5.37如上图所示
下面使用R的开源包FaultTree完成故障树生成与故障树分析

1.故障树生成

首先需要加载R的开源包“FaultTree”，可以前往https://r-forge.r-project.org/R/?group_id=2125 自行下载安装。或者直接在R中输入

install.packages("FaultTree", repos="http://R-Forge.R-project.org")

可能会提示你需要加载一些相关的依赖包，例如Rcpp之类的，按照提示操作即可
接着就可以生成故障树了

library(FaultTree)
library(FaultTree.SCRAM)
library(FaultTree.widget)
#建立故障树elec <- ftree.make(type="or", name="电网失效") #E1 id=1
elec <- addLogic(elec,at=1,type = "and",tag = "E2",name = "B站输入无效")   #E2 id=2
elec <- addLogic(elec,at=1,type = "and",tag = "E3",name = "C站输入无效")  #E3 id=3
elec <- addLogic(elec, at=1,type = "or", tag="E4", name="站B和站C仅由同一单线供电") #E4 id=4
elec <- addProbability(elec, at=2, prob=.01, tag="X1", name="输电线1故障") #X1 id=5
elec <- addProbability(elec,at=2,prob = .01,tag="X2",name="输电线2故障") #X2 id=6
elec <- addLogic(elec,at=2,type = "or",tag="E5",name = "来自站C的输电线路无电") #E5 id=7
elec <- addProbability(elec,at=3,prob = .01,tag = "X3",name = "输电线3故障") #X3 id=8
elec <- addLogic(elec,at=3,type = "or",tag = "E6",name = "来自站B的输电线路无电") #E6 id=9
elec <- addLogic(elec,at=4,type = "and",tag = "E7",name = "输电线23同时故障") #E7 id=10
elec <- addLogic(elec,at=4,type = "and",tag = "E8",name = "输电线13同时故障") #E8 id=11
elec <- addLogic(elec,at=4,type = "and",tag = "E9",name = "输电线12同时故障") #E9 id=12
elec <- addDuplicate(elec,at=7,dup_id = 8) #重复事件X3 id=13
elec <- addLogic(elec,at=7,type = "and",tag = "E10",name = "输电线45同时故障") #E10 id=14
elec <- addProbability(elec,at=14,prob = .01,tag = "X4",name = "输电线4故障") #X4 id=15
elec <- addProbability(elec,at=14,prob = .01,tag = "X5",name = "输电线5故障") #X5 id=16
elec <- addDuplicate(elec,at=9,dup_id = 14) # 重复事件E10  id=17
elec <- addDuplicate(elec,at=10,dup_id = 6) # 重复事件X2 id=18
elec <- addDuplicate(elec,at=10,dup_id = 8) # 重复事件X3 id=19
elec <- addDuplicate(elec,at=11,dup_id = 5) # 重复事件X1 id=20
elec <- addDuplicate(elec,at=11,dup_id = 8) # 重复事件X3 id=21
elec <- addDuplicate(elec,at=12,dup_id = 5) # 重复事件X1 id=22
elec <- addDuplicate(elec,at=12,dup_id = 6) # 重复事件X2 id=23
elec <- addDuplicate(elec,at=9,dup_id = 12) # 重复时间E9 id=24#测试树结构
test.ftree(elec)

2.故障树图形化输出

FaultTree库提供了html来完成故障树的图形化输出

#绘制故障树
ftree2html(elec,dir = "elec",write_file = TRUE) #手动将html改为GBK编码

到这一步，故障树就以HTML文件形式输出到了当前R的工作空间下。
【注意】
笔者在测试过程发现，输出的HTML文件调用了网上的Jquery和d3js来完成图形化加载
其中，d3js控制字符，Jquery提供动态页面效果。
由于FaultTree库默认设置的Jquery引用链接失效，会直接导致输出的HTML无法正常加载。
只需要将HTML文件中头部的java script引用地址修改为国内提供的Jquery引用地址即可。笔者使用的是百度提供的Jquery【虽然版本很旧，但也够用了】
对应文件中的

<script src="http://d3js.org/d3.v3.min.js" charset="GBK"></script>
<script src="https://libs.baidu.com/jquery/2.1.4/jquery.min.js"></script>

注意到，笔者将charset也改成了GBK，这样可以兼容中文

页面展示效果如下：

中间节点可以实现自由折叠，上图是将右侧折叠后的效果，被折叠的节点将以蓝色框显示。

3.定性分析

FaultTree库提供了最小割集的一键计算函数，通过设置cutsets函数的参数method，可以选用不同的计算方法，默认是下行法MOCUS

#最小割集
cutsets(elec) #默认下行法MOCUS

R中计算结果如下：

[[1]]
NULL[[2]][,1] [,2]
[1,] "X2" "X3"
[2,] "X1" "X3"
[3,] "X1" "X2"[[3]]
[1] "X3" "X4" "X5"

可以发现，二阶最小割集有3个，分别是X2X3，X1X3和X1X2；三阶最小割集有1个 X3X4X5
通过最小割集，我们可以轻松写出系统的结构函数
X2X3+X1X3+X1X2+X3X4*X5
也可以自定义一个结构函数求解函数struct()

#根据最小割集写出结构函数
struct<-function(cuts){level<-length(lengths(cuts))SF <- NULLfor(i in 1:level){length <- length(cuts[[i]]) #该阶段最小割集数量for(j in 1:length){if(i==1 & !is.null(cuts[i][j])){}else if(i==level & j==length){SF <- paste0(SF,t(cuts[[i]])[j]) #结尾}else if(i==1){SF <- paste0(SF,t(cuts[[i]])[j],"+") #一阶}else if(i>1 & j%%i==0){SF <- paste0(SF,t(cuts[[i]])[j],"+") #二阶及以上割集尾}else {  # i>1 & j%%1 != 0SF <- paste0(SF,t(cuts[[i]])[j],"*")}}}return(SF)
}

调用函数前记得先source加载一下

source("D:/Reliability/R-project/FaultTreeExtend/struct.R") #路径是你存储该R文件的绝对路径，建议不要有中文

接着就可以直接调用该函数了，输入参数是刚刚求得的最小割集，返回值是结构函数表达式

 struct(cuts)

输出效果如下：

[1] "X2*X3+X1*X3+X1*X2+X3*X4*X5"

至此便基本完成了故障树的定性分析

4.定量分析

这里我们考虑计算顶事件的发生概率，以及三种重要度（概率重要度，结构重要度和关键重要度）
【说明】
1、底事件概率重要度：这里的概率重要的是指伯恩鲍姆提出的概率重要度（即伯恩重要度MIF），其计算式：

2、底事件结构重要度：这里的结构重要度是简单地将伯恩重要度式中的概率全部视为0.5求得
3、底事件关键重要度：体现的是底事件概率发生微小变化时导致顶事件发生概率变化的相对变化率（CIF），其计算式：

4.1 顶事件的发生概率

假设底事件彼此之间相互独立，且发生概率p都为0.01。但由于最小割集之间存在着交集，即不同的割集内存在相同的底事件，故最小割集之间不相互独立，不可简单计算。
若考虑精确计算，往往有两种方法，一种是“不交化算法（即Sum of Disjoint Products(SDP)）”
一种是Inclusion and Exclusion Principle算法（IEP）

4.1.1 不交化算法SDP

在《GJB-Z 768A-1998 故障树分析指南》的附录A中，使用了不交化算法完成计算，由于该系统比较简单，且底事件数量相对较少，故计算过程不太复杂。笔者摘录了其计算过程如下，供参考学习。

4.1.2 Inclusion and Exclusion Principle算法（IEP或容斥原理）

图片摘自电子科技大学刘宇教授可靠性工程教学课件《Chapter06_Reliability Analysis(FTA)》

下面使用R来手动完成计算

p<-.01
#Inclusion and Exclusion Principle [精确计算]
F1 <- 3*p^2+p^3
F2 <- 3*p^3+2*p^4+p^5
F3 <- 3*p^5+p^3
F4 <- p^5
prob2 <- F1-F2+F3-F4
prob2

计算结果与不交化算法一致：

[1] 0.0002989801

4.1.3 近似计算

在工程应用上往往没有必要精确计算，该计算结果仅供参考，不属于设计指标。另外，精确计算的计算量大，且底事件发生概率值本身的误差就可能大到使得精确计算没有意义。在《GJB-Z 768A-1998 故障树分析指南》中推荐了集中近似计算方法，一般都可以满足工程要求。

容斥原理取部分项近似计算【以首项近似为例】

当最小割集事件互斥，该近似计算就等于精确计算值
容斥原理上下限平均近似计算

其实就是IEP的前两项的变形
独立近似计算
当各个最小割集中相同的底事件比较少，且发生概率很低时，可以假设各个最小割集之间相互独立

这里只用R实现独立近似计算

#结构函数：X2*X3+X1*X3+X1*X2+X3*X4*X5
#独立近似计算
#顶事件概率函数  1-(1-q2*q3)*(1-q1*q3)*(1-q1*q2)*(1-q3*q4*q5)
fun=expression(1-(1-q2*q3)*(1-q1*q3)*(1-q1*q2)*(1-q3*q4*q5))
imp <- deriv(fun, c("q1", "q2","q3","q4","q5"), func = TRUE)
prob <- imp(.01,.01,.01,.01,.01)[1] #独立近似计算顶事件概率
prob

计算结果如下：

[1] 0.0003009697

综上对比可以发现，即使使用独立近似计算，与精确计算的相对误差仅0.67%，在工程上是可以接受的

4.2 重要度计算

由于重要度计算在《GJB-Z 768A-1998 故障树分析指南》中的计算过程已经非常详细了，这里就只使用R对其实现近似计算并排序。
【注意】这里使用的顶事件概率函数（故障树的故障函数）是独立近似计算的函数，不是精确函数式，故与国军标指南上结果有一点点偏差，但不影响排序
首先，计算概率重要度

fun=expression(1-(1-q2*q3)*(1-q1*q3)*(1-q1*q2)*(1-q3*q4*q5))
imp <- deriv(fun, c("q1", "q2","q3","q4","q5"), func = TRUE)
MIF<- imp(.01,.01,.01,.01,.01) #独立近似计算概率重要度0.01999598 0.01999598 0.02009595 9.997e-05 9.997e-05
MIF

输出结果：

[1] 0.0003009697
attr(,"gradient")q1         q2         q3        q4        q5
[1,] 0.01999598 0.01999598 0.02009595 9.997e-05 9.997e-05

概率重要度排序：X3>X1=X2>X4=X5

接着计算结构重要度（将概率全部更改为0.5）

SIF<- imp(.5,.5,.5,.5,.5) #独立近似计算概率重要度0.01999598 0.01999598 0.02009595 9.997e-05 9.997e-05
SIF

输出结果：

[1] 0.6308594
attr(,"gradient")q1        q2        q3        q4        q5
[1,] 0.4921875 0.4921875 0.5976562 0.1054688 0.1054688

结构重要度排序：X3>X1=X2>X4=X5
最后，计算关键重要度

MIF<- c(0.01999598,0.01999598,0.02009595,9.997e-05,9.997e-05);
CIF<- MIF*.01/prob #关键重要度CIF
CIF

输出结果

[1] 0.664385150 0.664385150 0.667706747 0.003321597 0.003321597

关键重要度排序：X3>X1=X2>X4=X5

综上，三种重要度排序都一致，故在制定系统故障诊断时的核对清单时，应考虑参考该顺序进行排序，系统运行期间的监测部件优先级也可以参考该顺序。
重要度排序：X3>X1=X2>X4=X5

至此，故障树分析就全部完成了。

笔者在持续跟进故障树的R实现以及C++软件应用方面的工作，需要进一步交流可以给我留言。

参考资料

GJB-Z 768A-1998 故障树分析指南
电子科技大学 Professor Yu Liu 可靠性工程教学课件《Chapter06_Reliability Analysis(FTA)》

转载请注明出处，谢谢