基于馬爾可夫過程的快速SIL驗證方法研究

付 進,張 渺,崔銘芳,杜炘潔

(中國石油西南油氣田分公司安全環(huán)保與技術(shù)監(jiān)督研究院,四川 成都,610041)

0 引言

近年來，社會對工業(yè)過程安全日益重視，安全儀表系統(tǒng)的應(yīng)用也越來越普遍。為保證新建或改造的安全儀表系統(tǒng)能夠滿足現(xiàn)有的生產(chǎn)工藝需求，將生產(chǎn)風險控制在合理、可接受的范圍內(nèi)，需要對安全儀表系統(tǒng)開展安全完整性等級(safety integrity level，SIL)評價[1]。這也是完整性管理體系中的重要環(huán)節(jié)[2]。目前，常用的方法主要包括可靠性框圖[3]、故障樹[4]、馬爾可夫(Markov)模型法[5-6]等。其中，Markov模型法對安全儀表系統(tǒng)的各方面因素考慮全面。但該方法運算量大，時間成本相對較高。尤其在安全性能要求高的高冗余系統(tǒng)(例如大型高含硫天然氣凈化廠的安全儀表系統(tǒng))中，Markov模型法的缺陷將被進一步放大，因此并未得到廣泛應(yīng)用。

郭海濤等[7]通過對狀態(tài)轉(zhuǎn)移矩陣預(yù)先迭代相乘來降低計算量；Brissaud等[8]通過部分測試和全面測試的試驗分析，提出一種簡化的計算方法；王海清等[9]提出一種多階段馬爾科夫模型簡化算法。這些方法在一定程度上減少了計算量，但其實質(zhì)都是通過數(shù)學近似來簡化運算，從而引入了計算誤差。因此，本文以Markov過程理論為基礎(chǔ)，提出了快速馬爾可夫方法(rapid Morkov,R-Markov)。該方法在保證計算精度的前提下顯著降低時間成本，能夠應(yīng)用于大規(guī)模的安全儀表系統(tǒng)的SIL驗證。

1 Markov模型法計算平均要求時失效概率

1.1 SIL定級計算方法

安全儀表系統(tǒng)的SIL定級本質(zhì)是計算系統(tǒng)中安全儀表功能(safety instrumentation function，SIF)的平均要求時失效概率(probability of failure on demand,PFD)。該概率值pa為回路中傳感器、邏輯單元、最終元件的平均要求時失效概率ps、pl、pfe之和[10]。

pa=ps+pl+pfe

(1)

SIL與PFD的關(guān)系如表1所示。

表1 SIL與PFD的關(guān)系

根據(jù)每個安全儀表功能的平均要求時失效概率大小，驗證其是否滿足所需 SIL要求。如果計算所得的概率小于或等于相應(yīng)SIL等級的失效概率，說明按現(xiàn)有的配置狀態(tài)滿足SIL的要求；反之，則說明按現(xiàn)有的配置狀態(tài)不滿足所需SIL的要求，要對現(xiàn)有配置中的某一個元件或某幾個元件的冗余配置進行改進，例如逐級增加冗余配置。對于改進后的安全儀表功能作再次計算，直到最終概率滿足要求為止[11]。

1.2 Markov模型法

Markov模型法基于Markov過程理論，將系統(tǒng)劃分為若干獨立狀態(tài)，狀態(tài)間會以特定概率相互轉(zhuǎn)移。系統(tǒng)將來所處的狀態(tài)和系統(tǒng)的歷史狀態(tài)無關(guān)，只和當前狀態(tài)有關(guān)。Markov模型的這個特性恰好能夠解決電子電氣系統(tǒng)失效的指數(shù)概率密度問題。因此，Markov模型很適合用來進行安全儀表系統(tǒng)失效分析，包括多個影響可靠性的因素，例如結(jié)構(gòu)冗余、共因失效、自診斷、在線或離線測試維修、非理想的維修測試、周期性功能測試、單一或多個維修隊伍等。本文涉及的可靠性參數(shù)及計算方法如表2所示。表2中，M、S、Cd、CS、β分別代表平均修復(fù)時間、裝置誤動作重啟時間、危險覆蓋因子、安全覆蓋因子以及共因失效因子。

表2 可靠性參數(shù)及計算方法

安全儀表系統(tǒng)的過程周期由啟動、運行、失效以及修復(fù)后再運行等一系列事件序列組成。通過Markov模型可以高效、高精度地模擬系統(tǒng)中的靜態(tài)關(guān)系和復(fù)雜的動態(tài)變化。但模型的狀態(tài)數(shù)量會隨系統(tǒng)的復(fù)雜程度遞增。如果一個組件有5種失效模式，那么包括正常狀態(tài)在內(nèi)共有6種狀態(tài)，對應(yīng)的是6維狀態(tài)轉(zhuǎn)移矩陣。1oo2系統(tǒng)狀態(tài)轉(zhuǎn)移模型如圖1所示。

圖1 1oo2系統(tǒng)狀態(tài)轉(zhuǎn)移模型

在狀態(tài)0，兩個通道都可以正常運行。在狀態(tài)1和狀態(tài)2，一個通道發(fā)生危險失效，但另一個通道仍可執(zhí)行安全功能，系統(tǒng)仍處于安全狀態(tài)。另外，狀態(tài)1的危險失效被檢測到后可以進行在線修復(fù)，以修復(fù)率μ0回到狀態(tài)0。狀態(tài)3、狀態(tài)4、狀態(tài)5都是系統(tǒng)失效狀態(tài)：狀態(tài)3表示系統(tǒng)發(fā)生安全失效，狀態(tài)4表示系統(tǒng)發(fā)生檢測得到的危險失效，狀態(tài)5則代表系統(tǒng)發(fā)生未檢測到的危險失效。當發(fā)生共因失效時，系統(tǒng)從狀態(tài)0直接跳轉(zhuǎn)到狀態(tài)4或狀態(tài)5。狀態(tài)轉(zhuǎn)移矩陣P如式(2)所示。其中，∑表示矩陣元素所在行除該元素外其他元素之和[12]。

(2)

轉(zhuǎn)移矩陣P的元素Pij代表在一個時間段里，系統(tǒng)由狀態(tài)i轉(zhuǎn)移到狀態(tài)j的概率。定義初始狀態(tài)向量S0。該向量在1oo2結(jié)構(gòu)中為6維，代表系統(tǒng)共有6種狀態(tài)。若系統(tǒng)初始狀態(tài)正常，則S0=[1 0 0 0 0 0]。經(jīng)過一個時間間隔，系統(tǒng)處于各狀態(tài)的概率向量S1=S0P。經(jīng)過兩個時間間隔后，系統(tǒng)處于各狀態(tài)的概率向量S2=S1P。以此類推，Sn=Sn-1P或Sn=S0Pn。顯然，系統(tǒng)越復(fù)雜，存在的狀態(tài)數(shù)越多，狀態(tài)向量S的維數(shù)越多，轉(zhuǎn)移矩陣P越復(fù)雜，求解P的n次冪時花費的時間也指數(shù)倍增加。

假設(shè)第6個狀態(tài)為未檢測到的危險失效狀態(tài)、第5個狀態(tài)為檢測到的危險失效狀態(tài)、第4個狀態(tài)為安全失效狀態(tài)，則定義危險失效向量Vd=[0 0 0 0 1 1]T，安全失效向量Vs=[0 0 0 1 0 0]T。假設(shè)測試周期T為一年，約為8 760 h。經(jīng)過i小時后，系統(tǒng)狀態(tài)向量Si、要求時失效概率di為：

Si=S0Pi,i=1,2,...,8 760

(3)

di=S0PiVd,i=1,2,...,8 760

(4)

平均要求時失效概率為:

(5)

通過式(5)，分別計算出安全儀表功能中傳感器、邏輯單元、最終元件的平均要求時失效概率。由式(1)可求得所在SIF平均要求時失效概率，從而得到SIL等級。

2 Markov模型法優(yōu)化分析

2.1 Markov模型法時間復(fù)雜度分析

(6)

由式(6)可知，隨著T和n取值的增大，計算時間將呈指數(shù)上升。而T的取值通常大于8 760。這也是目前Markov模型法在SIL驗證中應(yīng)用較少的原因。對于大型安全儀表系統(tǒng)，出于時間成本考慮，一般會選擇其他方法。

2.2 Markov轉(zhuǎn)移矩陣冪運算優(yōu)化

由式(6)可知，時間開銷主要來自矩陣高次冪的累計求和。而由矩陣理論可知，如果方陣可以相似對角化，則存在對角矩陣D和可逆矩陣U，使P=UDU-1成立。Pi可由式(7)計算：

Pi=(UDU-1)i=UDU-1UDU-1×...×UDU-1=UDiU-1

(7)

對角矩陣D的元素對應(yīng)轉(zhuǎn)移矩陣P的特征值λ，因此Pi的計算就簡化為求P的n個特征值λ1、λ2、...、λn的i次冪。式(8)為優(yōu)化后算法計算平均要求時失效概率時的計算式。

i=1,2,...,8 760

(8)

該算法的時間復(fù)雜度為:

(9)

相比式(6)中的O(T2n3)，優(yōu)化后的時間復(fù)雜度顯著降低。隨著系統(tǒng)冗余程度的提高和測試周期的提高，即n和T越大，算法效率提升越明顯。但要使上面的結(jié)論成立，需要證明Markov轉(zhuǎn)移矩陣可以相似對角化。

2.3 證明轉(zhuǎn)移矩陣可相似對角化

矩陣可相似對角化的定義為：n階方陣存在n個線性無關(guān)的特征向量。該問題可轉(zhuǎn)化為求證方陣特征值互異。

(10)

考慮式(10)所示1oo1系統(tǒng)的轉(zhuǎn)移矩陣P，求解滿足特征方程|λE-P|=0的特征值λ，如式(11)所示。

(11)

式(11)化簡可以得到特征方程，如式(12)所示。

(λ-1)[λ3+λ2(λDD+λS+μ0+μSD-3)+λ(λDμ0+λDμSD-2λD-λDμ0+λSμ0-2λS+μ0μSD-2μ0-2μSD+3)+λDμ0μSD-λDμ0-λDμSD+λD-λDDμ0μSD+λDDμ0-λSμ0+λS-μ0μSD+μ0+μSD-1]=0

(12)

由式(12)可知：P中的一個特征值λ1=1，λ2、λ3、λ4為方程剩余項的解，方程中λS、λD、λDD數(shù)量級范圍為10-8～10-5，危險失效恢復(fù)率μ0、安全失效恢復(fù)率μSD分別為元件平均修復(fù)時間和裝置誤動作重啟時間的倒數(shù)，數(shù)量級范圍為10-2～10-1。將式(12)中失效參數(shù)按數(shù)量級合，并將特征方程簡化為f(λ)=λ3+Aλ2+Bλ+C。令f(λ)=0，A、B、C近似取值范圍分別為(-2.927,2.833)、(2.671,2.855)、(-0.838,-0,928)。網(wǎng)格分區(qū)生成特征方程曲線族，均與x軸相交于三個不同點，則λ2、λ3、λ4三個特征值互異。特征方程曲線如圖2所示。

圖2 特征方程曲線

將λ=1代入三次特征方程，得：

1+λDD+λS+μ0+μSD-3+λDμ0+λDμSD-2λD-λDμ0+λSμ0-2λS+μ0μSD-2μ0-2μSD+3+

λDμ0μSD-λDμ0-λDμSD+λD-λDDμ0μSD+λDDμ0-λSμ0+λS-μ0μSD+μ0+μSD-1=0

(13)

式(13)化簡結(jié)果為λDD=λD。根據(jù)表2中λDD=CdλD，而危險覆蓋因子Cd<1，可知λDD<λD。

對于1oo2系統(tǒng)的六階轉(zhuǎn)移矩陣P，|λE-P|=0的特征值λ如式(14)和式(15)所示。

(14)

(15)

式中：∑0為第0行元素之和；∑1為平行元素之和；∑2為第2行元素之和。

為求證特征方程不存在特征值為1的重根，將λ=1代入式(15)中的行列式，化簡得到特征多項式為：

Y=μ0μSD∑0∑1∑2-μ0μSD(2λDDNλD∑2+

2λDUNλDD∑1+∑2λDDC)

(16)

由表2可知:λDD<λD，λDDC<λD，λDDN<λD，λDUN<λD。因此，在λD<0.1時，式(16)可化簡為：

(17)

綜上可知，特征方程式(12)、式(14)不存在值為1的重根，即λ1、λ2、λ3、λ4互異。

由此證得上述Markov狀態(tài)轉(zhuǎn)移矩陣P可相似對角化。

3 試驗分析

3.1 天然氣集氣站安全儀表系統(tǒng)數(shù)據(jù)

試驗采用四川地區(qū)某天然氣集氣站的安全儀表系統(tǒng)可靠性參數(shù)。該集氣站共劃分為15個安全儀表功能，分別執(zhí)行火災(zāi)關(guān)斷、進站高高壓力切斷、進站低低壓力切斷和污水系統(tǒng)低液位關(guān)斷的安全功能。本次試驗選取其中最具代表性的5個安全儀表功能。傳感器、最終元件、邏輯單元的可靠性參數(shù)分別如表3、表4及表5所示。

表3 傳感器的可靠性參數(shù)

表4 最終元件的可靠性參數(shù)

表5 邏輯單元的可靠性參數(shù)

3.2 兩種方法結(jié)果對比

分別使用原始Markov模型法和R-Markov模型法計算5組SIF的平均要求時失效概率，分析比較兩種方法花費的時間和計算結(jié)果。每組數(shù)據(jù)分別進行10次試驗，時間取平均值。為消除干擾因素，只計算兩種方法核心算法耗時，不考慮構(gòu)建Markov轉(zhuǎn)移矩陣等公共計算開銷。兩組試驗均在相同軟件、硬件環(huán)境下進行。兩種方法計算結(jié)果如表6所示。

表6 兩種方法計算結(jié)果

3.3 系統(tǒng)冗余度和測試周期T的影響

從式(6)、式(9)所示兩種方法的計算復(fù)雜度可知，相比R-Markov模型法，Markov模型法耗時對轉(zhuǎn)移矩陣階數(shù)(即系統(tǒng)冗余度)和測試周期更為敏感。通過一組試驗，分別驗證系統(tǒng)冗余結(jié)構(gòu)和測試周期對計算時間的影響。試驗數(shù)據(jù)采用表3中SIF02的傳感器可靠性參數(shù)，僅計算傳感器的要求時失效概率。設(shè)定測試周期為8 760 h、13 140 h和17 520 h，分別在1oo1、1oo2、2oo3條件下進行運算。單元件運算結(jié)果對比如表7所示。冗余結(jié)構(gòu)和測試時間對效率提升的影響如圖3所示。

表7 單元件運算結(jié)果對比

圖3 冗余結(jié)構(gòu)和測試時間對效率提升的影響

3.4 結(jié)果分析

分析表6、表7數(shù)據(jù)和圖3可得如下結(jié)果。

①R-Markov模型法在5組SIF可靠性數(shù)據(jù)上提高約48～80倍效率。該計算結(jié)果和Markov模型法完全一致。

②結(jié)合表6中的時間比和表3～表5中各SIF回路的可靠性數(shù)據(jù)可知，SIF02、SIF03、SIF05中部分元件使用了更長的測試周期和更復(fù)雜的冗余結(jié)構(gòu)，此時優(yōu)化后算法的執(zhí)行效率有更顯著的提高。

③由表7和圖3可知，Markov模型法計算效率對系統(tǒng)冗余結(jié)構(gòu)和測試周期非常敏感，而R-Markov模型法的計算耗時始終能控制在非常低的水平。隨著冗余結(jié)構(gòu)和測試周期的提高，兩者計算時間的比值也顯著增加。當測試周期為一年時，R-Markov模型法在3種冗余結(jié)構(gòu)上計算效率分別提高至51.1、51.9和77.2倍。當測試周期為一年半時，計算效率則提高至70.7、72.4和118.6倍。當測試周期為兩年時，計算效率達到了91.5、94.6和164.3倍。隨著測試周期的進一步增加，這一差距將會更加明顯。

4 結(jié)論

本文在Markov模型法的基礎(chǔ)上，提出了一種新的驗算安全儀表系統(tǒng)SIL的方法，解決了Markov模型法計算量大、計算時間長的缺陷。通過數(shù)學原理，證明了該方法的可行性。本文采用真實天然氣集氣站的安全儀表系統(tǒng)可靠性數(shù)據(jù)進行了兩組對比試驗。試驗結(jié)果表明，在不同的冗余結(jié)構(gòu)和測試周期條件下，R-Markov模型法將SIL驗算的計算效率提高了48～160倍。分析試驗數(shù)據(jù)，R-Markov模型法即使用于大型復(fù)雜系統(tǒng)的SIL驗算任務(wù)，最終運算時間也可控制在10 s內(nèi)，為大型化工處理廠批量SIL驗算工作提供了高效的解決方案。