基于CUSPARSE的多階段任務系統(tǒng)可靠性并行計算方法

閆　華，崔譜龍，蘇永東，王　魁，張　齊

(陸軍勤務學院， 重慶　401311)

多階段任務系統(tǒng)(phased-mission system,PMS)是指包括多個任務階段，且系統(tǒng)配置和單元參數(shù)隨階段發(fā)生變化的系統(tǒng)[1]。PMS是一類在軍事應用領域常見的系統(tǒng)，如防空反導系統(tǒng)[2]、艦船作戰(zhàn)系統(tǒng)[3]及航天測控系統(tǒng)等[4]。現(xiàn)代作戰(zhàn)系統(tǒng)及其任務越來越復雜，其任務可靠性的規(guī)模也必然會越來越大，而現(xiàn)有的PMS任務可靠性計算方法對于大規(guī)模問題的效果都不太理想。

Markov方法是一類重要的PMS可靠性分析方法[5]，具有描述能力強，能夠正確處理階段內部和跨階段的部件依賴性，且理論上能夠求得精確解等優(yōu)點。但是，當系統(tǒng)中部件數(shù)量較多時，Markov模型面臨狀態(tài)空間爆炸問題，導致模型求解困難。Markov模型的常用求解方法包括一致化方法[6]、常微分方程方法[7]和迭代計算方法[8]等。上述算法通常只適合于小規(guī)模矩陣計算，對于大規(guī)模問題計算效率較低。為提高計算效率，研究人員進一步提出了BDD和Markov相結合的層次化方法[9]以及基于任務成功路徑的求解方法[10]。但是，現(xiàn)有的層次化方法仍然存在系統(tǒng)級模型的BDD節(jié)點數(shù)量隨階段數(shù)和部件數(shù)增加而迅速增加的問題；成功路徑方法的主要問題是計算結果偏低，且成功路徑數(shù)隨階段數(shù)增加而迅速增大。

近年來，圖形處理器(graphics processing unit,GPU)的發(fā)展極其迅猛，由于具有更高的并行處理能力和存儲帶寬，以及成本較低等優(yōu)點，GPU已經從傳統(tǒng)的圖形圖像處理拓展到通用計算領域，在稀疏矩陣向量乘[11,12]、線性方程組求解[13]等科學計算領域已經有了較為廣泛的應用。因此，引入并行計算思想，實現(xiàn)對傳統(tǒng)可靠性求解算法的GPU并行計算，能夠有效提高模型的求解效率。

1　PMS任務可靠性的一致化算法

Markov模型求解算法主要包括一致化方法(uniformization method,UM)、常微分方程方法和迭代計算方法等。由于UM算法具有簡潔直觀、易于實現(xiàn)等優(yōu)點，相比Runge-Kutta和前向Euler方法，UM算法具有較高的計算效率[14]。因此，首先介紹UM算法的主要計算過程，并基于統(tǒng)一計算設備架構(CUDASparse)庫函數(shù)，實現(xiàn)對該算法的并行計算。

令Q表示轉移速率矩陣，Markov模型求解中的核心計算是矩陣指數(shù)運算eQt。eQt可以寫為如下的無窮級數(shù)形式

(1)

若以v(t)表示系統(tǒng)在t時刻的狀態(tài)概率向量，可得

(2)

令A=QT，則式(2)可寫為

(3)

式(3)中的第k項可以由第k-1項計算得到。因此，可以利用迭代計算思想，通過計算式(3)中的前k項，得到v(t)的近似值。但是，由于轉移速率矩陣Q中的對角線元素均為負數(shù)項，若直接由式(3)計算，會導致較大的舍入誤差。因此，一致化方法中采用以下步驟對公式中矩陣進行處理。

A=γ(R-I)

(4)

將式(4)代入狀態(tài)概率向量的基本計算公式v(t)=eAt·v(0)中，可得到

(5)

因此，連續(xù)時間Markov鏈可以看作是單步轉移概率矩陣為R的離散時間Markov鏈和平均速率為γ的Poisson過程的組合。這一過程稱為連續(xù)時間Markov鏈的一致化(uniformization)，利用該過程求解狀態(tài)概率向量v(t)的方法稱為一致化方法。

(6)

(7)

即當?shù)螖?shù)K滿足上述條件時，停止迭代，當前結果精度在可接受誤差范圍內。

2　一致化算法的CUDASparse并行實現(xiàn)

2.1　CUDA

CUDA(Compute Unified Device Architecture)是Nvidia提出的面向通用計算的并行編程開發(fā)平臺，也是當前用于通用計算的主流GPU架構。CUDA包含了不同層次的調用接口，其中包括一組運行時API、一組設備驅動API以及CUDA提供的庫函數(shù)API。CUDA平臺函數(shù)調用的層次結構如圖1所示。

CUDA驅動API能夠直接控制底層硬件結構，CUDA運行時API則是對驅動API的封裝。程序設計中既可以調用CUDA驅動API實現(xiàn)對硬件更高效的控制，也可以使用CUDA運行時API更便捷、更直觀地實現(xiàn)CUDA平臺提供的計算模式。CUDA包括一系列的數(shù)學工具庫，其中CUSPARSE工具庫主要針對稀疏線性代數(shù)運算進行優(yōu)化[15]。

2.2　基于CUSPARSE的一致化并行算法

Markov可靠性模型求解過程中，參與運算的主要是轉移速率矩陣，該矩陣是一個典型的稀疏矩陣。稀疏矩陣算法可利用壓縮存儲結構，只存儲和處理非零元素，從而大幅降低存儲空間需求及計算復雜度。但是由于在運算中引入了大量的離散間接尋址操作，導致計算訪存比很低，并且計算過程中的訪存軌跡與矩陣的稀疏結構相關，很難發(fā)掘其時空局部性，因此，在傳統(tǒng)CPU處理器上稀疏矩陣算法的計算效率很低。通過引入并行計算思想，實現(xiàn)基于GPU的UM并行計算，能夠有效提高算法效率。

首先采用按行壓縮存儲方案(compressed row storage,CRS)對轉移速率矩陣進行壓縮存儲。CRS利用3個數(shù)組value、colind和rowptr，分別逐行存儲矩陣中的非零元素，非零元素對應的列索引，以及每行第一個元素在數(shù)組value中的索引值。

UM算法的并行實現(xiàn)，可通過調用CUDA平臺中CUSPARSE工具庫中的函數(shù)實現(xiàn)。主要利用CUSPARSE庫中的函數(shù)cusparseDcsrmv(a,b,x,y,M)，實現(xiàn)了形如v=a*op(M)*v1+b*v2的表達式運算。函數(shù)cusparseDcsrmv的各參數(shù)含義如下：M為以CRS格式存儲的稀疏矩陣，v1和v2為向量，a和b為標量；op(M)表示是否對矩陣M進行轉換，有3種形式：M、MT和MH，M表示不轉換，MT表示對矩陣M進行轉置，MH表示將M轉換為共軛矩陣。

由公式(6)及其中y(k)的定義可知，UM運算主要涉及對矩陣R及相應向量的乘積運算，這部分運算的并行實現(xiàn)可直接調用CUSPARSE庫函數(shù)cusparseDcsrmv()完成；同時，根據(jù)R=I+A·Δt可知，矩陣R通過對轉移速率矩陣A的運算得到，因此，可編寫相應的核函數(shù)，實現(xiàn)對矩陣R的并行計算。綜上所述，基于CUSPARSE的UM算法并行實現(xiàn)，主要包括兩個核函數(shù)Kernel1和Kernel2，且Kernel1的計算結果為Kernel2提供輸入，因此，兩個核函數(shù)實際上為串行關系。Kernel1主要根據(jù)公式R=I+A·Δt，實現(xiàn)由矩陣A計算得到矩陣R；Kernel2通過調用cusparseDcsrmv()函數(shù)，根據(jù)式(6)實現(xiàn)對v(t)的并行計算，并在每次迭代過程中，根據(jù)式(7)判斷是否達到停止條件。Kernel1和Kernel2的算法步驟分別如下。

核函數(shù)Kernel1中，由于矩陣A采用CRS壓縮存儲方案，算法中對應該矩陣的存儲數(shù)組分別為value_A，colind和rowptr；由于矩陣R與矩陣A相比，只有非零元素值發(fā)生了變化(矩陣A與單位矩陣相加，不會影響非零元素的位置，因為轉移速率矩陣A中，對角線元素總是非零元素)，因此，算法中對應矩陣R的存儲數(shù)組分別為value_R，colind和rowptr。

算法1　Kernel1算法核心步驟

算法2　Kernel2算法核心步驟

3　算例分析

并行算法實現(xiàn)基于Nvidia推出的CUDA 7.5版本，實驗平臺CPU為Inter Core i5-4210 @ 1.70 GHz，GPU為Nvidia Geforce GT 730M顯卡。計算過程中，將UM算法的可接受誤差設置為ε=1.0×10-10。

記有多階段任務T，該任務可劃分為3個任務階段P1、P2和P3，持續(xù)時間分別為90 s、150 s和90 s。假設所有設備的平均修復時間(mean time to repair,MTTR)和平均故障間隔時間(mean time between failures,MTBF)均為30 min和300 min。任務T各階段的系統(tǒng)故障樹如圖2所示。

任務T中3個階段參與任務的設備數(shù)分別為10、13和16，各階段矩陣規(guī)模分別為1 024×1 024，8 192×8 192和65 536×65 536，非零元素分別為4 275，51 282和514 233。

采用并行UM算法和傳統(tǒng)UM算法，可得任務T的可靠性計算時間如表1所示，相應的各階段任務可靠性計算結果如表2所示。

表1　并行UM和傳統(tǒng)UM的可靠性計算時間　　s

由表1可知，在任務T的各階段計算中， 并行UM比傳統(tǒng)UM算法分別能達到1.00，4,32和7.24倍的加速效果。由表1也可以看出，雖然并行UM算法具有一定的加速效果，但并不明顯，主要原因一是硬件性能限制；二是由于算法基于CUSPARSE函數(shù)庫實現(xiàn)，UM并行計算效率依賴于所調用庫函數(shù)；同時，由于CUSPARSE庫中并沒有直接實現(xiàn)UM算法的函數(shù)，導致Kernel2算法過程不夠簡潔高效。比如，由Kernel2算法步驟可以看出，Step8和Step9中兩次調用了函數(shù)cusparseDcsrmv()，因為每次迭代中y(k)必須單獨保存，但在Step8中，y(k)作為中間計算結果無法保存，因此，Step9中再次調用函數(shù)cusparseDcsrmv()，單獨計算并保存y(k)。

表2　并行UM和傳統(tǒng)UM的可靠性計算結果

同時，由表2可知，并行UM與傳統(tǒng)UM算法的計算結果完全一致。這是因為，并行UM算法只是將部分運算轉移至GPU中進行，并沒有改變算法具體計算過程和步驟。同時也可看出，利用式(7)，能夠有效控制算法的計算精度。

4　結論

針對大規(guī)模條件下PMS任務可靠性求解困難的問題，將并行計算思想引入傳統(tǒng)UM算法，提出了基于CUSPARSE的并行UM算法。該算法利用Nvidia提出的并行計算平臺CUDA，通過調用平臺工具庫CUSPARSE中的函數(shù)，實現(xiàn)了對UM算法的并行計算，以及算法的誤差控制。實例分析表明，并行UM算法與傳統(tǒng)UM算法相比，其計算結果完全一致，且并行UM計算效率具備優(yōu)勢，但并不明顯。這一方面是由于實驗平臺中GPU的性能限制，另一方面因為調用CUDA函數(shù)庫導致性能受限。綜上所述，基于CUSPARSE的PMS任務可靠性并行算法，比傳統(tǒng)算法具有較高的效率，但由于受到調用CUDA庫函數(shù)的限制，下一步可考慮基于CUDA平臺自主實現(xiàn)UM并行計算。

參考文獻：

[1]WU Xinyang,WU Xiaoyue.Extended Object-Orient Petri Net Model for Mission Reliability Simulation of Repairable PMS with Common Cause Failures[J].Reliability Engineering & System Safety,2015,136:109-119.

[2]劉斌,武小悅.基于多階段貝葉斯網絡的反導系統(tǒng)任務可靠性建模[J].裝備指揮技術學院學報,2012,23(1):75-78.

[3]狄鵬.考慮多階段任務的艦船任務可靠性分配方法[J].海軍工程大學學報,2014,26(3):108-112.

[4]LU Jimin,WU Xiaoyue.A new reliablility evaluation method for spaceflight telemetry,tracking and control system[C]//International Conference on Quality,Reliability,Risk,Maintenance,and Safety Engineering,2013.

[5]閆華.基于Markov方法的多階段任務系統(tǒng)可靠性分析綜述[J].兵器裝備工程學報,2016,37(6):92-96.

[6]AAD P A.VAN MOORSEL,SANDERS WILLIAM H.Transient Solution of Markov Models by Combining Adaptive and Standard Uniformization[J].IEEE Transactions on Reliability,1997,46(3):430-440.

[7]STEWART W J.Introduction of the Numerical Solution of Markov Chains[M].Princeton,NJ:Princeton University Press,1994.

[8]ANTOINE RAUZY.An Experimental Study on Iterative Methods to Compute Transient Solutions of Large Markov Models[J].Reliability Engineering & System Safety,2004,86(1):105-115.

[9]DAZHI WANG,TRIVEDI KISHOR S.Reliability Analysis of Phased-Mission System with Independent Component Repairs[J].IEEE Transactions on Reliability,2007,56(3):540-551.

[10] LU Ji-min.Reliability analysis of large phased-mission systems with repairable components based on success-state sampling[J].Reliability Engineering & System Safety,2015,142:123-133.

[11] 梁添.基于GPU的稀疏矩陣運算優(yōu)化研究[D].武漢：華中科技大學,2012.

[12] 王迎瑞,任江勇,田榮.基于GPU的高性能稀疏矩陣向量乘及CG求解器優(yōu)化[J].計算機科學,2013,40(3):46-49.

[13] 柳有權,尹康學,吳恩華.大規(guī)模稀疏線性方程組的GMRES-GPU快速求解算法[J].計算機輔助設計與圖形學學報,2011,23(4):553-560.

[14] WU Xiaoyue,YAN Hua,LI Liriong.Numerical Method for Reliability Analysis of Phased-Mission System using Markov Chains[J].Communications in Statistics-Theory and Methods,2012,41(21):3960-3973.

[15] Nvidia Corporation.NVIDIA C.CUSPARSE Library[Z].

2015.