基于馬爾可夫過程的快速SIL驗(yàn)證方法研究

付 進(jìn),張 渺,崔銘芳,杜炘潔

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

0 引言

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

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

1 Markov模型法計(jì)算平均要求時(shí)失效概率

1.1 SIL定級(jí)計(jì)算方法

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

pa=ps+pl+pfe

(1)

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

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

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

1.2 Markov模型法

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

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

安全儀表系統(tǒng)的過程周期由啟動(dòng)、運(yùn)行、失效以及修復(fù)后再運(yùn)行等一系列事件序列組成。通過Markov模型可以高效、高精度地模擬系統(tǒng)中的靜態(tài)關(guān)系和復(fù)雜的動(dòng)態(tài)變化。但模型的狀態(tài)數(shù)量會(huì)隨系統(tǒng)的復(fù)雜程度遞增。如果一個(gè)組件有5種失效模式，那么包括正常狀態(tài)在內(nèi)共有6種狀態(tài)，對(duì)應(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，兩個(gè)通道都可以正常運(yùn)行。在狀態(tài)1和狀態(tài)2，一個(gè)通道發(fā)生危險(xiǎn)失效，但另一個(gè)通道仍可執(zhí)行安全功能，系統(tǒng)仍處于安全狀態(tài)。另外，狀態(tài)1的危險(xiǎn)失效被檢測(cè)到后可以進(jìn)行在線修復(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ā)生檢測(cè)得到的危險(xiǎn)失效，狀態(tài)5則代表系統(tǒng)發(fā)生未檢測(cè)到的危險(xiǎn)失效。當(dāng)發(fā)生共因失效時(shí)，系統(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代表在一個(gè)時(shí)間段里，系統(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)過一個(gè)時(shí)間間隔，系統(tǒng)處于各狀態(tài)的概率向量S1=S0P。經(jīng)過兩個(gè)時(shí)間間隔后，系統(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í)花費(fèi)的時(shí)間也指數(shù)倍增加。

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

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

(3)

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

(4)

平均要求時(shí)失效概率為:

(5)

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

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

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

(6)

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

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

由式(6)可知，時(shí)間開銷主要來(lái)自矩陣高次冪的累計(jì)求和。而由矩陣?yán)碚摽芍绻疥嚳梢韵嗨茖?duì)角化，則存在對(duì)角矩陣D和可逆矩陣U，使P=UDU-1成立。Pi可由式(7)計(jì)算：

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

(7)

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

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

(8)

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

(9)

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

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

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

(10)

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

(11)

式(11)化簡(jiǎn)可以得到特征方程，如式(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中的一個(gè)特征值λ1=1，λ2、λ3、λ4為方程剩余項(xiàng)的解，方程中λS、λD、λDD數(shù)量級(jí)范圍為10-8～10-5，危險(xiǎn)失效恢復(fù)率μ0、安全失效恢復(fù)率μSD分別為元件平均修復(fù)時(shí)間和裝置誤動(dòng)作重啟時(shí)間的倒數(shù)，數(shù)量級(jí)范圍為10-2～10-1。將式(12)中失效參數(shù)按數(shù)量級(jí)合，并將特征方程簡(jiǎn)化為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軸相交于三個(gè)不同點(diǎn)，則λ2、λ3、λ4三個(gè)特征值互異。特征方程曲線如圖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ǎn)結(jié)果為λDD=λD。根據(jù)表2中λDD=CdλD，而危險(xiǎn)覆蓋因子Cd<1，可知λDD<λD。

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

(14)

(15)

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

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

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時(shí)，式(16)可化簡(jiǎn)為：

(17)

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

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

3 試驗(yàn)分析

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

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

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

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

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

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

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

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

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

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

表7 單元件運(yùn)算結(jié)果對(duì)比

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

3.4 結(jié)果分析

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

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

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

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

4 結(jié)論

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