一种基于马尔可夫过程进行安全仪表功能的功能安全验证的方法

摘要

本发明公开了一种基于马尔可夫模型进行安全仪表功能的功能安全验证的方法。采用马尔可夫模型计算单一设备采用1选1结构时的需求时失效率，继而以单一设备1选1结构的需求时失效率为基础计算多设备共失效的概率，最后以多设备共失效的概率计算冗余结构的需求时失效率。从而按照IEC61508的要求由冗余结构的平均需求时失效率验证安全仪表功能的安全完整性等级。本发明的优点：1)通过马尔可夫模型涵盖系统的多种状态，反映系统的动态行为，从而使提高失效率的计算精度；2)由1选1结构的需求时失效率计算结果间接计算冗余结构的需求时失效率，从而避免了直接对冗余结构进行马尔可夫建模和模型计算的困难；3)考虑共因失效时结合了β因子模型和二项式失效率模型，提高了对共因失效率的计算精度。

说明书

发明领域

本发明涉及一种基于马尔可夫过程的对安全仪表功能进行功能 安全验证的方法，它主要应用于化工过程，炼油过程，生物制药过程 等流程工业的功能安全验证领域。

背景技术

安全对于经济、环境和人身安全至关重要，因此安全相关系统 (Safety Related System，SRS)被广泛应用于生产过程之中。安全仪 表系统(Safety Instrumented System，SIS)是安全相关系统的一个类 别。按照IEC61511的定义，安全仪表系统是指由传感器、逻辑控制 器和执行器组成的，能够执行一项或多相安全仪表功能(Safety  Instumented Function，SIF)的仪表系统。在一定时间、一定条件下， 安全仪表系统能够成功地执行其安全功能的可能性，称为安全完整性 等级(Safety Integrity Level，SIL)，其数值代表着安全仪表系统使过 程风险降低的数量级。对于安全仪表功能的设计进行安全完整性等级 验证是十分重要的。按照IEC61508的定义，过程工业领域的安全仪 表功能多数情况下属于低要求操作模式(Low Demand)。根据 IEC61508的规定，低要求操作模式(Low Demand)下的安全完整性 等级以要求时平均失效概率(Average Probability of Failure on  Demand，PFDavg)作为指标进行划分，其对应关系如表1所示：

表1低要求操作模式下SIL与PFD_avg的对应关系

  安全完整性等级(SIL)   要求时平均失效概率(PFD_avg)   SIL4   ≥10^-5且＜10^-4  SIL3   ≥10^-4且＜10^-3  SIL2   ≥10^-3且＜10^-2  SIL1   ≥10^-2且＜10^-1

因此，对化工过程，炼油过程，生物制药过程等流程工业的功能 安全验证最为重要的就是需要时平均失效率的计算。

为了提高系统的安全完整性等级，实际应用中常利用多个设备构 成冗余结构。安全仪表系统的传感器、逻辑控制器和执行机构部分均 可以有冗余结构。冗余的设备可以采用串联降低受控系统误停车的概 率，或以并联的方式配置以降低安全仪表系统安全仪表功能失效的概 率，或者混用这两种配置方式。安全仪表系统传感器、逻辑控制器和 执行机构三部分的冗余结构通常分别以一个“N选M”形式的表达式 描述，其意义为N个设备构成该部分的冗余结构，当这N个设备中 有大于等于M个设备判断受控系统出现故障，则该部分判断受控系 统出现故障。构成冗余结构的各个设备的失效并非相互独立，冗余结 构中的设备往往会发生共因失效。因此在失效率时若将设备失效当做 独立事件考虑将会高估冗余结构的可靠性。

计算PFDavg的方法主要有故障树、可靠性框图、马尔可夫模型 法等。马尔可夫模型将系统归于若干种状态。一个状态以某种概率转 移到其他状态。利用马尔可夫模型可以分析系统的行为和可靠性。马 尔可夫模型可以分为离散时间马尔可夫模型和连续时间马尔可夫模 型。其中离散时间马尔可夫模型不但易于求解，而且更加适合具有周 期性的安全功能测试和维修的安全仪表系统。相比于其他方法，马尔 可夫模型可以包括设备的多种失效模式，通过一次建模可以得到多个 可靠性指标，可以反映系统的动态行为。但是对状态较多的具有冗余 结构的系统，马尔可夫模型的建立和求解可能十分困难。

发明内容

本发明的一个目的，是利用马尔可夫模型可以涵盖系统多种状 态，反映系统动态行为的特点同时，克服其对于冗余结构建模困难的 缺点，不对N选M冗余结构进行直接的马尔可夫模型建立，而是由 1选1结构的计算结果间接推出N选M冗余结构的需求时失效率， 从而避免了对复杂冗余系统建立马尔可夫模型的困难。在考虑冗余结 构的共因失效时，综合使用了β因子模型和二项式失效率模型，从而 使对多设备共因失效的计算更为准确。

根据本发明的一个方面，提供了一种基于马尔可夫模型进行安全 仪表功能的功能安全验证方法，其特征在于包括如下步骤：

确定1选1结构需求时失效率；

确定多设备共失效概率；

确定冗余结构需求时失效率，由多设备共失效概率得到冗余结构 的需求时失效率；

确定冗余结构平均需求时失效率；

确定安全完整性等级。

附图说明

图1显示了根据本发明的一个实施例的方法的流程图。

具体实施方式

本发明的目的在于利用马尔可夫模型可以涵盖系统多种状态， 反映系统动态行为的特点同时，克服其对于冗余结构建模困难的缺 点，不对N选M冗余结构进行直接的马尔可夫模型建立，而是由1 选1结构的计算结果间接推出N选M冗余结构的需求时失效率，从 而避免了对复杂冗余系统建立马尔可夫模型的困难。在考虑冗余结构 的共因失效时，利用了冲击模型，从而使对多设备共因失效的计算更 为准确。

-1选1结构需求时失效率的确定

1选1结构是最为简单的安全仪表系统结构。根据马尔可夫模型， 其离散的按小时的状态转移矩阵利用可以写作：

$P=\left(\begin{array}{ccccc}1-\Sigma & {\lambda }^{\mathrm{SD}}+{\lambda }^{\mathrm{SU}}& {\lambda }^{\mathrm{DD}}& {\lambda }^{\mathrm{DU}}{C}_{\mathrm{TI}}& {\lambda }^{\mathrm{DU}}(1-C_{\mathrm{TI}})\\ {\mu }_{\mathrm{SD}}& 1-\Sigma & 0& 0& 0\\ {\mu }_0& 0& 1-\Sigma & 0& 0\\ 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 1\end{array}\right)$

矩阵中各参数的意义为：

∑：该行其他元素之和，单位为hr^-1；

λ^SD：设备的检测到的安全失效率，单位为hr^-1；

λ^SU：设备的未检测到的安全失效率，单位为hr^-1；

λ^DD：设备的检测到的危险失效率hr^-1；

λ^DU：设备的未检测到的危险失效率，单位为hr^-1；

C_TI：周期性功能测试的覆盖率，单位为1；

μ_SD：设备的启动率，单位为hr^-1。与启动时间T_SD(单位hr)的关 系为μ_SD＝1/T_SD；

μ₀：设备的修复率，单位为hr^-1。与平均修复时间MTTR(单位hr) 的关系为μ₀＝1/MTTR。

该状态转移矩阵第i行第j列的元素为经过一个月后设备由状态i 转变为状态j的概率。矩阵对应的1选1结构的5个状态依次为：

1：正常；

2：安全失效；

3：检测到的危险失效；

4：未检测到，但周期性功能测试可以发现的危险失效；

5：未检测到，且周期性功能测试不会发现的危险失效。

一个月按744小时计(一个月31天)，按月的状态转移矩阵为：

P_m＝P⁷⁴⁴

周期性功能也对应着一个状态转移矩阵：

$W=\left(\begin{array}{ccccc}1& 0& 0& 0& 0\\ 1& 0& 0& 0& 0\\ 1& 0& 0& 0& 0\\ 1& 0& 0& 0& 0\\ 0& 0& 0& 0& 1\end{array}\right)$

其对应的状态次序与P相同。

假设系统的初始状态为正常，即初始状态的概率分布向量为：

S₀＝[1 0 0 0 0]

功能测试周期为TI(单位为月)，则经过t个月(t为正整数)，对 t与TI功能的比值向下取整可得测试周期总数为：

余数为：

r＝t-n_TI·TI

于是t个月后系统的状态概率分布为：

$S\left(t\right)=\left(\begin{array}{c}{S}_0{\left(P_m^{\mathrm{TI}}W\right)}^{n_{\mathrm{TI}}-1}{P}_m^{\mathrm{TI}},r=0\\ S_0{\left(P_m^{\mathrm{TI}}W\right)}^{n_{\mathrm{TI}}}{P}_m^r,r\ne 0\end{array}\right)$

S(t)的后3个元素(即检测到的危险失效；未检测到，但周期性功 能测试可以发现的危险失效；未检测到，且周期性功能测试不会发现 的危险失效)之和即t个月后系统的需求时失效率。取

V＝[0 0 1 1 1]^T

则

λ(t)＝S(t)V

-多设备共失效的确定

多设备共失效包括独立失效和共因失效两种情况。对于共因失效， 本发明采用的共因失效模型以二项式失效率模型(Binomial Failure  Rate Model)为基础，同时常用的β因子模型的估计简化了二项式失 效率模型参数的获取。

二项式失效率模型需要考虑了冲击发生的概率μ和冲击下单一设 备失效的概率p。本发明将p作为通过经验估计获得的已知参数。

忽略独立失效和冲击同时发生的情况，那么对于N选M冗余结构， 设各个设备的共因失效因子分别估计为β_i，冲击时设备失效的概率估 计为p_i，由1选1结构的马尔可夫模型计算得到设备的需求时失效率 分别为λ_i。其中i=1，2，L，N

经过推导可得N个设备中k个设备失效的概率为

${\lambda }_{k/N}=\underset{\{i_1,i_2,L,i_k\}\subseteq \{\mathrm{1,2},L,N\}}{\Sigma }\underset{j\in \{i_1,i_2,L,i_k\}}{\Pi }({\lambda }_j-{\beta }_j{\lambda }_j/p_j)\underset{l\notin \{i_1,i_2,L,i_k\}}{\underset{l=\mathrm{1,2},L,N}{\Pi }}(1-{\lambda }_l)$

$+\mu \underset{\{i_1,i_2,L,i_k\}\subseteq \{\mathrm{1,2},L,N\}}{\Sigma }\underset{j\in \{i_1,i_2,L,i_k\}}{\Pi }{p}_j\underset{l\notin \{i_1,i_2,L,i_k\}}{\underset{l=\mathrm{1,2},L,N}{\Pi }}(1-p_l)$

其中

$\mu =\underset{i=\mathrm{1,2},L,N}{\max }\{{\beta }_i{\lambda }_i/p_i^2\}$

若N个设备相同，上式可以简化为(各设备的参数相同，故统一 用不带下标的符号表示)

${\lambda }_{k/N}=C_N^k{(\lambda -\mathrm{\beta \lambda }/p)}^k{(1-\lambda )}^{N-k}+C_N^k\mathrm{\beta \lambda }{p}^{k-2}{(1-p)}^{N-k}$

-冗余结构需求时失效率的确定

N选M冗余结构当且仅当M个或M个以上设备处于正常或安全 失效状态则系统处于正常或安全失效状态。这相当于当且仅当N-M 个或N-M个以上设备处于危险失效状态时，系统危险失效。N选M 结构的需要时失效率为

${\lambda }_{\mathrm{MooN}}=\underset{k=N-M+1}{\overset{N}{\Sigma }}{\lambda }_{k/N}$

-平均要求时失效率的计算

由已经得到的各月N选M系统的需求时失效率λ_MooN(t)，可以得到 系统运行时间LT(单位为年)内的平均需求时失效率PFD_avg

${\mathrm{PFD}}_{\mathrm{avg}}=\frac{1}{12\mathrm{LT}}\underset{t=1}{\overset{12\mathrm{LT}}{\Sigma }}{\lambda }_{\mathrm{MooN}}\left(t\right)$

-安全完整性等级SIL的确定

将PFD_avg的负对数向下取整可得系统的安全完整性等级，即

实施例

把根据本发明的方法应用到一个具体的安全仪表的安全完整性验 证上。

该设备的相关数据如表2所示

表2实施例的安全仪表相关数据

  数据名称/单位   数值   检测到的安全失效率/小时^-1  1e-5   未检测到的安全失效率/小时^-1  5e-6   检测到的危险失效率/小时^-1  2e-5

  未检测到的危险失效率/小时^-1  1e-5   平均修复时间/小时   24   启动时间/小时   24   功能测试周期/月   12   周期性功能测试覆盖率   0.95   运行时间/年   3

冗余结构选择3选2结构。

根据数据可以得到该设备按小时的状态转移矩阵为

$P=\left(\begin{array}{ccccc}9.99955e-1& 1.50000e-5& 2.00000e-5& 9.50000e-6& 5.00000e-7\\ 4.16667e-2& 9.58333e-1& 0& 0& 0\\ 4.16667e-2& 0& 9.58333e-1& 0& 0\\ 0& 0& 0& 1.00000& 0\\ 0& 0& 0& 0& 1.00000\end{array}\right)$

求该矩阵的744次幂可得按月的状态转移矩阵为

$P_m=\left(\begin{array}{ccccc}9.91760e-1& 3.57119e-4& 4.76159e-4& 7.03611e-3& 3.70321e-4\\ 9.91998e-1& 3.57205e-1& 4.76273e-4& 6.80993e-3& 3.58417e-4\\ 9.91998e-1& 3.57205e-1& 4.76273e-4& 6.80993e-3& 3.58417e-4\\ 0& 0& 0& 1.00000& 0\\ 0& 0& 0& 0& 1.00000\end{array}\right)$

由按月的状态转移矩阵可得设备采用1选1结构时在运行时间内 各月的需求时失效率如表3所示

表3实施例的安全仪表采用1选1结构时各月的需求时失效率

  月序   需求时失效率   月序   需求时失效率   月序   需求时失效率   1   7.883e-3   13   1.212e-2   25   1.633e-2   2   1.523e-2   14   1.943e-2   26   2.362e-2   3   2.252e-2   15   2.669e-2   27   3.085e-2   4   2.976e-2   16   3.390e-2   28   3.803e-2   5   3.695e-2   17   4.106e-2   29   4.515e-2   6   4.408e-2   18   4.816e-2   30   5.222e-2   7   5.116e-2   19   5.521e-2   31   5.924e-2   8   5.819e-2   20   6.221e-2   32   6.621e-2   9   6.516e-2   21   6.915e-2   33   7.313e-2   10   7.209e-2   22   7.605e-2   34   7.999e-2   11   7.896e-2   23   8.289e-2   35   8.680e-2   12   8.578e-2   24   8.968e-2   36   9.357e-2

取共因失效因子β＝0.1，本例中冲击下单一设备失效概率取为 p＝0.75，则

λ_2/3＝0.075λ+2.253λ²-2.253λ³

λ_3/3＝0.075λ+0.651λ³

故λ_2oo3＝λ_2/3+λ_3/3＝0.15λ+2.253λ²-1.602λ³

按该等式换算2oo3结构在运行时间内各月的需求时失效率如表4 所示

表4实施例的安全仪表采用2oo3结构时各月的需求时失效率

  月序   需求时失效率   月序   需求时失效率   月序   需求时失效率   1   1.322e-3   13   2.145e-3   25   3.044e-3   2   2.802e-3   14   3.754e-3   26   4.778e-3   3   4.503e-3   15   5.579e-3   27   6.724e-3   4   6.418e-3   16   7.613e-3   28   8.874e-3   5   8.537e-3   17   9.846e-3   29   1.122e-2   6   1.085e-2   18   1.227e-2   30   1.375e-2   7   1.336e-2   19   1.488e-2   31   1.646e-2   8   1.604e-2   20   1.766e-2   32   1.934e-2   9   1.890e-2   21   2.062e-2   33   2.239e-2   10   2.192e-2   22   2.373e-2   34   2.559e-2   11   2.510e-2   23   2.700e-2   35   2.895e-2   12   2.845e-2   24   3.042e-2   36   3.245e-2

对各月的需求时失效率取平均值，得到平均需求时失效率

PFD_avg＝1.465e-2

故该设备采用2oo3冗余结构时的安全完整性等级为1。