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

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

(陸軍勤務(wù)學(xué)院， 重慶　401311)

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

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

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

1　PMS任務(wù)可靠性的一致化算法

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

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

(1)

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

(2)

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

(3)

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

A=γ(R-I)

(4)

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

(5)

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

(6)

(7)

即當(dāng)?shù)螖?shù)K滿足上述條件時(shí)，停止迭代，當(dāng)前結(jié)果精度在可接受誤差范圍內(nèi)。

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

2.1　CUDA

CUDA(Compute Unified Device Architecture)是Nvidia提出的面向通用計(jì)算的并行編程開(kāi)發(fā)平臺(tái)，也是當(dāng)前用于通用計(jì)算的主流GPU架構(gòu)。CUDA包含了不同層次的調(diào)用接口，其中包括一組運(yùn)行時(shí)API、一組設(shè)備驅(qū)動(dòng)API以及CUDA提供的庫(kù)函數(shù)API。CUDA平臺(tái)函數(shù)調(diào)用的層次結(jié)構(gòu)如圖1所示。

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

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

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

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

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

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

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

算法1　Kernel1算法核心步驟

算法2　Kernel2算法核心步驟

3　算例分析

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

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

任務(wù)T中3個(gè)階段參與任務(wù)的設(shè)備數(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算法，可得任務(wù)T的可靠性計(jì)算時(shí)間如表1所示，相應(yīng)的各階段任務(wù)可靠性計(jì)算結(jié)果如表2所示。

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

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

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

同時(shí)，由表2可知，并行UM與傳統(tǒng)UM算法的計(jì)算結(jié)果完全一致。這是因?yàn)?，并行UM算法只是將部分運(yùn)算轉(zhuǎn)移至GPU中進(jìn)行，并沒(méi)有改變算法具體計(jì)算過(guò)程和步驟。同時(shí)也可看出，利用式(7)，能夠有效控制算法的計(jì)算精度。

4　結(jié)論

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

參考文獻(xiàn)：

[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]劉斌,武小悅.基于多階段貝葉斯網(wǎng)絡(luò)的反導(dǎo)系統(tǒng)任務(wù)可靠性建模[J].裝備指揮技術(shù)學(xué)院學(xué)報(bào),2012,23(1):75-78.

[3]狄鵬.考慮多階段任務(wù)的艦船任務(wù)可靠性分配方法[J].海軍工程大學(xué)學(xué)報(bào),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方法的多階段任務(wù)系統(tǒng)可靠性分析綜述[J].兵器裝備工程學(xué)報(bào),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ùn)算優(yōu)化研究[D].武漢：華中科技大學(xué),2012.

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

[13] 柳有權(quán),尹康學(xué),吳恩華.大規(guī)模稀疏線性方程組的GMRES-GPU快速求解算法[J].計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào),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.