大規模多階段任務系統馬爾可夫可靠性模型的存儲和計算

閆華，高黎，王魁，漆磊

（后勤工程學院后勤信息與軍事物流工程系，重慶401311）

大規模多階段任務系統馬爾可夫可靠性模型的存儲和計算

閆華，高黎，王魁，漆磊

（后勤工程學院后勤信息與軍事物流工程系，重慶401311）

由于馬爾可夫模型在進行多階段任務系統的可靠性分析時，系統狀態隨部件增加呈指數增長，從而導致大規模條件下模型求解所需的存儲量和計算量十分巨大。而根據馬爾可夫模型中轉移速率矩陣Q的取值規律和稀疏特性，給出了矩陣Q中元素qij基于狀態二進制表示的計算公式，并提出了一種Q矩陣壓縮存儲（QMCS）方法。在模型壓縮存儲的基礎上，進一步提出了基于Krylov子空間的可靠性求解算法。通過算例對比了不同壓縮存儲方案和不同求解算法的存儲量、計算時間和可靠性結果，分析表明基于QMCS和Krylov子空間的模型求解方法具有較高的存儲和計算效率，特別是在矩陣規模較大的情況下，該方法的計算耗時優于其他方法，且結果精度也能滿足可靠性計算需求。

系統評估與可行性分析；可靠性評估；多階段任務系統；壓縮存儲；Krylov子空間

0　引言

基于馬爾可夫模型進行多階段任務系統（PMS）的可靠性分析時，主要涉及到t時刻系統處于各狀態的概率向量［1］，狀態概率向量的求解需要對轉移速率矩陣（或稱無窮小生成子）Q進行運算。由于馬爾可夫模型中的狀態空間呈指數增長，當系統中單元數目較多時，Q矩陣的維數將非常大，導致存儲量和計算量都十分巨大［2-3］。針對大規模馬爾可夫可靠性模型的求解，可從模型預處理和模型求解算法兩方面進行研究，提高模型的存儲和運算效率。

馬爾可夫可靠性模型的轉移速率矩陣Q中包含大量零元素，是稀疏矩陣。因此，可利用稀疏矩陣壓縮存儲方法對矩陣Q進行預處理。稀疏矩陣壓縮儲存大致可分為兩類:通用存儲方案和特殊矩陣存儲方案［4］。通用存儲方案對矩陣中非零元素分布不作任何假設，如行壓縮存儲（CRS）［5］；特殊矩陣存儲方案針對某些具有特殊結構的矩陣，如帶狀矩陣存儲方法［6］。根據存儲的基本單元，又可分為基于元素的存儲方案和基于塊的存儲方案。CRS即為基于元素的存儲方案；基于塊的存儲方案包括固定塊存儲（FBS）和按行壓縮分塊存儲（BCRS）［7］等。上述方法各有優劣，CRS通用性較強，FBS適合于矩陣中非零塊長度相同的情形，BCRS適合于矩陣中存在很多非零塊且塊大小各不相同。同時，文獻［8-9］基于CRS、FBS和BCRS等通用壓縮存儲方案，對不同壓縮方案下的PMS任務可靠性分析方法進行了研究，但并沒有提出一種針對馬爾可夫模型特點的高效壓縮方法；文獻［10］提出了一種基于相似狀態的轉移概率矩陣壓縮方法，但該方法的前提是已知模型的轉移概率矩陣，因此并不適用于基于馬爾可夫過程的任務可靠性計算。

馬爾可夫模型的常用計算方法包括一致化方法［11］、常微分方程方法［12］和迭代計算方法［13］等。上述算法通常只適合于小規模矩陣計算，對于大規模問題計算效率較低。Lu等［14］提出了一種基于任務成功路徑的馬爾可夫可靠性模型求解方法，該方法能夠有效降低模型的計算復雜度，但主要問題是計算精度較低，且任務成功路徑數隨階段數增加而迅速增大。Krylov子空間方法是一種空間投影技術，通過將大規模問題投影至小規模子空間，得到問題的近似解。因此，在模型壓縮存儲的基礎上，可利用Krylov子空間技術推導轉移速率矩陣的近似求解算法，提高馬爾可夫模型的求解效率。

本文給出了轉移速率矩陣Q在系統狀態二進制表示方法下的計算公式，總結了矩陣Q中元素的取值規律；提出了Q矩陣壓縮存儲（QMCS）方法；并在QMCS基礎上，提出了壓縮存儲和Krylov子空間相結合的PMS可靠性分析方法。通過算例對QMCS的壓縮存儲效率和Krylov子空間方法的計算效率進行了對比分析，結果證明本文所提方法能夠有效提高模型的存儲和計算效率。

1　轉移速率矩陣的壓縮存儲方法

1.1基于狀態二進制表示的Q矩陣計算

以Si表示系統的第i個狀態，則Si可以由一個二進制字符串c1c2…cn表示。其中，ck為系統中第k個單元的狀態，正常為1，失效為0.為便于描述馬爾可夫模型中轉移速率矩陣Q的取值規律，定義狀態相交權重的概念。

令Si和Sj為任意的兩個狀態（i≠j），狀態字符串的長度等于系統中單元數，記為n.令cik、cjk分別表示狀態Si和Sj中第k位（0＜k≤n）元素，若存在cik=cjk，則稱狀態Si和Sj相交。定義相交元素的個數之和為狀態相交權重，記為wij.例如，假設S3= 10011，S5=00111，則兩狀態中具有相同元素的位數分別為第2、4和5位，因此，狀態相交權重w35=3.

根據系統狀態的二進制表示和狀態相交權重，推導了轉移速率矩陣Q的計算方法。矩陣Q中元素qij表示當系統處于狀態Si，轉移至狀態Sj的速率，其一般表示方法為qij=viPij［15］.其中，vi為系統在狀態Si處的轉移速率，Pij為系統由狀態Si轉移至Sj的概率。

通常假設在極短的時間Δt內，系統只能發生一次故障或修復，因此，對于狀態Si和Sj，只有當wij=n-1時系統才有可能發生轉移。假設從Si到Sj，發生故障或修復的為第k個單元，對應其在初始狀態Si中的單元狀態為cik.并假設各單元的失效與修復時間均服從指數分布，分別以λk和μk表示單元的失效率與修復率，則有

（1）式表示在狀態Si處，可能發生的轉移有n種，因此，所有單元的的轉移速率之和就是系統在狀態Si處的轉移速率vi；（2）式表示第k個單元發生轉移的速率與總的轉移速率之比，即為系統由Si轉移至Sj的概率。

綜合（1）式、（2）式，可得到如下qij的計算公式:

式中:wij=n，表示qij為矩陣Q中的對角線元素，根據Q中對角線元素等于該行所有非零元素之和的相反數，可得qii=-∑jqij（j≠i）.

根據轉移速率矩陣Q的計算公式，總結qij取值規律如下:1）吸收態對應行中的元素全為0，吸收態即為表示系統任務失敗的狀態；2）若兩狀態相交權重wij=n或wij=n-1，則對應qij為非零元素，否則對應qij為零元素。

1.2Q矩陣壓縮存儲方案

根據轉移速率矩陣中元素qij的取值規律及其計算公式，提出一種QMCS方法。Q矩陣每行與每列均對應一個狀態，將狀態按照其二進制字符串對應的十進制值從小到大依此排列，狀態對應的十進制值即代表了矩陣元素的行索引和列索引。以4個狀態的系統為例，假設狀態00為吸收態。矩陣Q及其對應的行狀態和列狀態如下:

由狀態的二進制表示可以推出其行索引，因此，采用3個數組存儲非零元素:1）數組RStates存儲系統中的所有正常狀態的二進制字符串，并按照狀態的十進制值從小到大依此存儲，狀態的十進制值代表了該狀態所對應的行索引，數組SysStates的數據類型為字符型；2）數組UnitMTBF存儲系統中n個單元的失效率，數據類型為浮點型；3）數組UnitMTTR存儲系統n個單元的修復率，數據類型為浮點型。當qij非零時，可通過比較行狀態和列狀態，根據（3）式計算得到qij，所需數據從數組UnitMTBF和UnitMTTR中獲取。以（4）式中的矩陣Q為例，QMCS下的數據結構如表1所示。

表1　Q矩陣壓縮存儲方案下的數據結構示例Tab.1 Data structure of Q matrix in compressed storage scheme

以行狀態01為例，該狀態對應行索引為1，根據qij取值規律，該行中非零元素對應的列狀態應分別為00和11，對應矩陣元素為q10和q13，再加對角線上的非零元素，則該行中的所有非零元素為q10、q11和q13，且根據（3）式可得q10=λ2，q13=μ1，q11= -（μ1+λ2）.

為了便于對各方法的存儲量進行分析，定義以下記號:因矩陣Q為方陣，記N為Q的維數；n為系統中的單元數，由此可知N=2n；nnz為Q中非零元素的數目，r為系統中正常狀態的數量。稀疏矩陣存儲時，非零元素的值采用浮點型數據存儲，其他的輔助數組采用整型存儲。并假設浮點型數據需要8個字節，整型數據需要4個字節，字符型數據需要8個字節。若令M1表示QMCS方案下所需的存儲量，則有

由此可見，在QMCS方案下，所需的存儲量僅與正常狀態的數量和單元數量有關，具體的qij值可以在使用過程中，由（3）式進行動態計算。與其他存儲方案相比，由于不需要直接存儲非零元素值，因此，該方案能夠有效地節省存儲空間。

2　壓縮存儲模型的K rylov子空間求解方法

PMS可靠性模型的求解可以歸結為對模型中狀態概率向量d（t）的計算。對于馬爾可夫模型，根據Chapman-Kolmogorov后向方程得到d P（t）/d t= QP（t），可推得P（t）=eQt，其中P（t）為t時刻的轉移概率矩陣。若令d（0）表示系統初始時刻的狀態概率向量，根據d（0）P（t）=d（t）可得

式中:當矩陣Q較大時，其計算將十分困難。

Krylov子空間是指由形如p（A）d的向量張成的子空間［16］，A為矩陣，p（A）為由矩陣A構成的多項式，則子空間Km（A，d）可表示為

對于（6）式，若令A=QT，u（t）=dT（t），z=

u（t）的任意m-1階多項式展開都是Krylov子空間Km（A，z）中的元素，因此，基于Krylov子空間投影，可得到如下的近似計算公式［17］:

式中:p=‖z‖2；Bm=［b1，b2，…，bm］為子空間Km（A，p）中的一組標準正交基，由Arnoldi過程構造［16］；Hm為m×m階的上Hessenberg矩陣，也由Arnoldi過程得到；e1=（0，…，0，1）T.

由（9）式可以看出，其中基本運算為矩陣與向量的乘積運算。因此，根據上節QMCS方案，給出在該方案下的矩陣向量乘積運算的算法如下所示。

算法 QMCS下矩陣與向量乘積運算

Step1:For k←0 to r-1 Do

Step2: Si=SysStates［k］；

Step3: 計算行狀態Si對應的十進制索引值RIndex，并令i=RIndex；

Step4: For l←0 to n-1 Do

Step5: 根據wij=n-1的原則構造列狀態Sj；

Step6: 計算列狀態Si對應的十進制索引值CIndex，并令j=CIndex；

Step7: If cl=0

Step8: qij=μl=UnitMTTR［l］；

Step9: Else

Step10: qij=λl=UnitMTBF［l］；

Step11: End If

Step12: y［i］+=qijx［j］；

Step13: qii+=-qij；

Step14: End For

Step15: y［i］+=qiix［i］；

Step16:End For

其中，Si為Q矩陣中第i行對應的狀態，Sj為第j列對應的狀態。對于正常狀態Si，只有Si與Sj的相交權重wij=n或wij=n-1時，對應qij為非零元素。wij=n時，即有Si=Sj，對應qij是矩陣Q中對角線上的元素；wij=n-1即與行狀態Si相比只有一位發生變化的狀態。

QMCS實現了對模型中轉移速率矩陣的高效壓縮存儲，同時，上述算法實現了QMCS下的矩陣向量計算，為利用（9）式進行投影后的模型近似求解提供了運算基礎。

3　算例分析

CRS是最為通用的方法［18］，對矩陣中非零元素的分布不作任何假設；Haque［7］通過實驗得出結論:在各種稀疏矩陣存儲方案中，FBS的性能最優。FBS通常可記為FBS l，其中l為規定的非零塊長度，l的長度取決于具體的計算平臺，常用的有FBS2和FBS3.因此，算例分析部分采用CRS方案、FBS2方案和FBS3方案與QMCS方案進行對比。

令M2、M3分別表示CRS方案和FBS l方案下的存儲量，βl表示矩陣Q中長度為l的非零塊的個數，則有［7］

采用文獻［9］中算例對QMCS壓縮存儲與該存儲方案下的Krylov子空間求解方法進行分析。記該多階段任務為T，包括兩個任務階段P1和P2，且各階段任務持續時間均為10min.假設系統中各單元的平均修復時間和平均故障間隔時間均為30min和30 000min.任務T各階段的系統故障樹分別如圖1和圖2所示，其中A1、A2，…，A13為系統部件。

圖1　任務T中階段P1的系統故障樹Fig.1 System fault tree of Phase P1in Task T

根據文獻［9］中的實驗結果，相比Runge-Kutta方法和前向Euler方法，一致化方法（UM）具有較高的計算效率。因此，為了驗證Krylov子空間算法的有效性，采用UM算法與Krylov子空間算法進行對比。

任務T中兩個階段參與任務的單元數分別為7和13，階段矩陣Q1和Q2的規模分別為128×128和8 192×8 192，非零元素分別為208和37 128，兩個階段系統中的正常狀態數分別為26和2 652.矩陣Q1和Q2中長度為2的非零塊數目為26和2 652；長度為3的非零塊數目為17和1 478.可得矩陣Q1和Q2在CRS方案、FBS2方案、FBS3方案和QMCS方案下的存儲量如表2所示。

圖2　任務T中階段P2的系統故障樹Fig.2 System fault tree of Phase P2in Task T

表2　Q1和Q2在不同存儲方式下所需的存儲量比較Tab.2 Storage space of Q1and Q2

Krylov子空間算法和UM算法在CRS方案、FBS2方案、FBS3方案和QMCS方案存儲方式下各階段的可靠性計算時間如表3所示，兩種算法的可靠性計算結果如表4所示。由于采用矩陣壓縮存儲對模型進行預處理，因此，表3中的可靠性計算時間分別為:t1表示模型預處理時間，即對矩陣的壓縮處理；t2為求解算法運行時間，不包括模型預處理時間；t3表示可靠性計算的整個耗時，包括預處理時間和算法運行時間，即t3=t1+t2.表3中CRS-Krylov方案表示模型預處理采用CRS壓縮方案，任務可靠性求解采用Krylov子空間算法；CRS-UM方案表示模型壓縮采用CRS方法，并利用UM算法求解任務可靠性。表3中其他符號具有類似的含義。

表3　不同存儲方式下利用Krylov子空間算法和UM算法的可靠性計算時間比較（s）Tab.3 Computation times of Krylov subspace and UM algorithms in different storage modes（s）

表4　基于Krylov子空間算法和UM算法的任務可靠性計算結果Tab.4 Calculated task reliabilities of Krylov subspace and UM algorithms

對比表2～表4中存儲量、計算時間和任務可靠性等計算結果，可以得到如下結論:

1）CRS、FBS2、FBS3和QMCS 4種存儲方案中，QMCS方案所需存儲量最?。ㄒ姳?），有效提高了模型中轉移速率矩陣的存儲效率。

2）在QMCS方案下，并未直接存儲轉移速率矩陣Q中的非零元素qij，qij是在模型求解過程中動態生成。因此，若只考慮算法運行時間（即表3中t2），由于QMCS方案下需要額外的qij計算過程，不論Krylov算法或是UM算法，QMCS方案下的算法運行時間t2都大于同階段的其他存儲方案（見表3）；但是，CRS方案和FBS l方案下均包含繁瑣的矩陣壓縮處理過程，該過程耗時（表3中t1）相對算法運行時間，通常不可忽略。因此，若對比不同存儲方案下可靠性計算的整個耗時（表3中t3），同一任務階段中QMCS方案的計算耗時均要優于其他存儲方案，且矩陣規模越大時其優勢越明顯（見表3）。綜上所述，相比CRS方案、FBS2方案和FBS3方案，QMCS方案在存儲和計算耗時方面均更高效。

3）以任務階段P2為例，采用Krylov子空間算法求解時，其運行時間（表3中t2）的最大值和最小值分別為0.671 s和0.024 s；而UM算法運行時間的最大值和最小值分別為2.028 s和0.078 s（見表3）。從中可以看出，相比UM算法，Krylov子空間算法的計算速度更快。對比階段P1的相關數據，也可以得到同樣的結論。

4）對比Krylov子空間和UM兩種算法的可靠性計算結果，任務階段P1和P2的可靠性結果誤差分別為7.00×10-11和8.20×10-11.由此可知，Krylov子空間算法也具有較高的計算精度。

4　結論

在馬爾可夫模型中，當系統單元數量較多時，其狀態空間呈指數增長，如單元數量為20時，轉移速率矩陣Q的維數就達到了百萬級。本文針對大規模馬爾可夫模型的存儲和計算問題，提出了基于矩陣壓縮存儲技術的可靠性求解方法。實例分析表明，相比CRS、FBS2和FBS3 3種壓縮存儲方案，QMCS方案所需的存儲空間最??；同時，對比分析了不同壓縮存儲方案下Krylov子空間算法和UM算法的運行耗時和計算精度，結果表明QMCS方案下的Krylov子空間方法具有較高的計算效率，特別是在矩陣規模較大的情況下，其計算耗時優于UM算法。綜上所述，QMCS方案和Krylov子空間技術相結合的求解方法較好地解決了大規模馬爾可夫模型的存儲和計算問題。

（References）

［1］ Alam M，Al-Saggaf U M.Quantitative reliability evaluation of repairable phased-mission systems using Markov approach［J］.IEEE Transactions on Reliability，1986，35（5）:498-503.

［2］ Peng R，Zhai Q，Xing L D.Reliability of demand-based phasedmission systems subject to fault level coverage［J］.Reliability Engineering and System Safety，2014，121（1）:18-25.

［3］ Wu X Y，Wu X Y.Extended object-oriented Petri net model for mission reliability simulation of repairable PMS with common cause failures［J］.Reliability Engineering and System Safety，2015，136（4）:109-119.

［4］ Shahnaz R，Usman A，Chughtai IR.Review of storage techniques for sparse matrices［C］∥2005 Pakistan Section Multitopic Conference.Karachi，Pakistan:IEEE，2005.

［5］ Zhang J，Wan J，Li F.Efficient sparse matrix-vector multiplication using cache oblivious extension quadtree storage format［J］. Future Generation Computer Systems，2016，54（1）:490-500.

［6］ Golub G H，Loan C F V.Matrix computations［M］.4th ed.Baltimore，US:Johns Hopkins University Press，2013.

［7］ Haque SA.Acomputational study of sparse matrix storage schemes［D］.Lethbridge，Canada:University of Lethbridge，2008.

［8］ 黎麗榮.基于Markov模型的大型PMS任務可靠性分析方法［D］.長沙:國防科學技術大學，2011. LI Li-rong.Method of mission reliability analysis of large PMS based on Markov model［D］.Changsha:National University of Defense Technology，2011.（in Chinese）

［9］ Wu X Y，Yan H，Li L R.Numerical method for reliability analysis of phased-mission system using Markov chains［J］.Communications in Statistics—Theory and Methods，2012，41（21）: 3960-3973.

［10］ 石磊，姚瑤.馬爾可夫預測模型中轉移概率矩陣的壓縮與應用［J］.計算機應用，2007，27（11）:2746-2749. SHI Lei，YAO Yao.Compression and application of transiton probability matrix in Markov prediction model［J］.Computer Applications，2007，27（11）:2746-2749.（in Chinese）

［11］ Moorsel A PV，Sanders W H.Transient solution of Markov models by combining adaptive and standard uniformization［J］.IEEE Transactions on Reliability，1997，46（3）:430-440.

［12］ Stewart W J.Introduction of thenumerical solution of Markov chains［M］.Princeton，NJ:Princeton University Press，1994.

［13］ Rauzy A.An experimental study on iterative methods to compute transient solutions of large Markov models［J］.Reliability Engineering and System Safety，2004，86（1）:105-115.

［14］ Lu JM，Wu X Y.Reliability evaluation of generalized phased-mission systems with repairable components［J］.Reliability Engineering and System Safety，2014，121（1）:136-145.

［15］ Ross SM.Introduction to probability models［M］.Burlington，US:Academic Press，2006.

［16］ Saad Y.Iterative methods for sparse linear systems［M］.New York，US:PWS Publishing Company，1996.

［17］ 閆華，武小悅.航天測控通信系統可靠性分析的改進Krylov投影算法［J］.系統工程與電子技術，2012，34（10）:2180-2186. YAN Hua，WU Xiao-yue.Improved Krylov subspace projection algorithm for reliability analysis of TT＆amp;C and communication system［J］.Systems Engineering and Electronics，2012，34（10）:2180-2186.（in Chinese）

［18］ Simecek I，Langr D，Tvrdik P.Space-efficient sparse matrix storage formats for massively parallel systems［C］∥Proceedings of14th IEEE International Conference on High Performance Computing and Communications.Liverpool:IEEE，2012:54-60.

Storage and Computation of Markov Reliability Model for Large-scale Phased-mission System

YAN Hua，GAO Li，WANG Kui，QI Lei
（Department of Logistics Information＆amp;Logistics Engineering，Logistic Engineering University of PLA，Chongqing 401311，China）

When Markov model is used to analyzed the reliability of phased-mission system，the system state grows exponentially with the increase in the number of components，thus resulting in a huge storage space and calculated amount resolved by the model.According to the element value rules and sparsity of the transition rate matrix Q in Markovmodel，the formula of computing the elements qijis derived based on binary description of states，and a Q-matrix compressed storage scheme（QMCS）is proposed.A reliability computing algorithm using Krylov subspace method is proposed based on the model compressed storage scheme.Taking a practical phased-mission system for example，the required storage spaces，computation times and reliability results of different compressed storage schemes and different algorithms are compared.The analysis results show that the method combining QMCS and Krylov subspace method has higher efficiency in storage and computation.Especially in the case of a large matrix，the QMCS-Krylov method is superior to other methods both in computation time and accuracy.

system assessment and feasibility；reliability evaluation；phase-mission system；compressed storage；Krylov subspace

TP202+.1；N945.17

A

1000-1093（2016）09-1715-06

10.3969/j.issn.1000-1093.2016.09.023

2016-02-02

國家自然科學基金項目（71401172）

閆華（1983—），男，講師。E-mail:yanhua_8304@163.com