闫　华，崔谱龙，苏永东，王　魁，张　齐

(陆军勤务学院， 重庆　401311)

多阶段任务系统(phased-mission system,PMS)是指包括多个任务阶段，且系统配置和单元参数随阶段发生变化的系统[1]。PMS是一类在军事应用领域常见的系统，如防空反导系统[2]、舰船作战系统[3]及航天测控系统等[4]。现代作战系统及其任务越来越复杂，其任务可靠性的规模也必然会越来越大，而现有的PMS任务可靠性计算方法对于大规模问题的效果都不太理想。

Markov方法是一类重要的PMS可靠性分析方法[5]，具有描述能力强，能够正确处理阶段内部和跨阶段的部件依赖性，且理论上能够求得精确解等优点。但是，当系统中部件数量较多时，Markov模型面临状态空间爆炸问题，导致模型求解困难。Markov模型的常用求解方法包括一致化方法[6]、常微分方程方法[7]和迭代计算方法[8]等。上述算法通常只适合于小规模矩阵计算，对于大规模问题计算效率较低。为提高计算效率，研究人员进一步提出了BDD和Markov相结合的层次化方法[9]以及基于任务成功路径的求解方法[10]。但是，现有的层次化方法仍然存在系统级模型的BDD节点数量随阶段数和部件数增加而迅速增加的问题；成功路径方法的主要问题是计算结果偏低，且成功路径数随阶段数增加而迅速增大。

近年来，图形处理器(graphics processing unit,GPU)的发展极其迅猛，由于具有更高的并行处理能力和存储带宽，以及成本较低等优点，GPU已经从传统的图形图像处理拓展到通用计算领域，在稀疏矩阵向量乘[11,12]、线性方程组求解[13]等科学计算领域已经有了较为广泛的应用。因此，引入并行计算思想，实现对传统可靠性求解算法的GPU并行计算，能够有效提高模型的求解效率。

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

Markov模型求解算法主要包括一致化方法(uniformization method,UM)、常微分方程方法和迭代计算方法等。由于UM算法具有简洁直观、易于实现等优点，相比Runge-Kutta和前向Euler方法，UM算法具有较高的计算效率[14]。因此，首先介绍UM算法的主要计算过程，并基于统一计算设备架构(CUDASparse)库函数，实现对该算法的并行计算。

令Q表示转移速率矩阵，Markov模型求解中的核心计算是矩阵指数运算eQt。eQt可以写为如下的无穷级数形式

(1)

若以v(t)表示系统在t时刻的状态概率向量，可得

(2)

令A=QT，则式(2)可写为

(3)

式(3)中的第k项可以由第k-1项计算得到。因此，可以利用迭代计算思想，通过计算式(3)中的前k项，得到v(t)的近似值。但是，由于转移速率矩阵Q中的对角线元素均为负数项，若直接由式(3)计算，会导致较大的舍入误差。因此，一致化方法中采用以下步骤对公式中矩阵进行处理。

A=γ(R-I)

(4)

将式(4)代入状态概率向量的基本计算公式v(t)=eAt·v(0)中，可得到

(5)

因此，连续时间Markov链可以看作是单步转移概率矩阵为R的离散时间Markov链和平均速率为γ的Poisson过程的组合。这一过程称为连续时间Markov链的一致化(uniformization)，利用该过程求解状态概率向量v(t)的方法称为一致化方法。

(6)

(7)

即当迭代次数K满足上述条件时，停止迭代，当前结果精度在可接受误差范围内。

2　一致化算法的CUDASparse并行实现

2.1　CUDA

CUDA(Compute Unified Device Architecture)是Nvidia提出的面向通用计算的并行编程开发平台，也是当前用于通用计算的主流GPU架构。CUDA包含了不同层次的调用接口，其中包括一组运行时API、一组设备驱动API以及CUDA提供的库函数API。CUDA平台函数调用的层次结构如图1所示。

CUDA驱动API能够直接控制底层硬件结构，CUDA运行时API则是对驱动API的封装。程序设计中既可以调用CUDA驱动API实现对硬件更高效的控制，也可以使用CUDA运行时API更便捷、更直观地实现CUDA平台提供的计算模式。CUDA包括一系列的数学工具库，其中CUSPARSE工具库主要针对稀疏线性代数运算进行优化[15]。

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

Markov可靠性模型求解过程中，参与运算的主要是转移速率矩阵，该矩阵是一个典型的稀疏矩阵。稀疏矩阵算法可利用压缩存储结构，只存储和处理非零元素，从而大幅降低存储空间需求及计算复杂度。但是由于在运算中引入了大量的离散间接寻址操作，导致计算访存比很低，并且计算过程中的访存轨迹与矩阵的稀疏结构相关，很难发掘其时空局部性，因此，在传统CPU处理器上稀疏矩阵算法的计算效率很低。通过引入并行计算思想，实现基于GPU的UM并行计算，能够有效提高算法效率。

首先采用按行压缩存储方案(compressed row storage,CRS)对转移速率矩阵进行压缩存储。CRS利用3个数组value、colind和rowptr，分别逐行存储矩阵中的非零元素，非零元素对应的列索引，以及每行第一个元素在数组value中的索引值。

UM算法的并行实现，可通过调用CUDA平台中CUSPARSE工具库中的函数实现。主要利用CUSPARSE库中的函数cusparseDcsrmv(a,b,x,y,M)，实现了形如v=a*op(M)*v1+b*v2的表达式运算。函数cusparseDcsrmv的各参数含义如下：M为以CRS格式存储的稀疏矩阵，v1和v2为向量，a和b为标量；op(M)表示是否对矩阵M进行转换，有3种形式：M、MT和MH，M表示不转换，MT表示对矩阵M进行转置，MH表示将M转换为共轭矩阵。

由公式(6)及其中y(k)的定义可知，UM运算主要涉及对矩阵R及相应向量的乘积运算，这部分运算的并行实现可直接调用CUSPARSE库函数cusparseDcsrmv()完成；同时，根据R=I+A·Δt可知，矩阵R通过对转移速率矩阵A的运算得到，因此，可编写相应的核函数，实现对矩阵R的并行计算。综上所述，基于CUSPARSE的UM算法并行实现，主要包括两个核函数Kernel1和Kernel2，且Kernel1的计算结果为Kernel2提供输入，因此，两个核函数实际上为串行关系。Kernel1主要根据公式R=I+A·Δt，实现由矩阵A计算得到矩阵R；Kernel2通过调用cusparseDcsrmv()函数，根据式(6)实现对v(t)的并行计算，并在每次迭代过程中，根据式(7)判断是否达到停止条件。Kernel1和Kernel2的算法步骤分别如下。

核函数Kernel1中，由于矩阵A采用CRS压缩存储方案，算法中对应该矩阵的存储数组分别为value_A，colind和rowptr；由于矩阵R与矩阵A相比，只有非零元素值发生了变化(矩阵A与单位矩阵相加，不会影响非零元素的位置，因为转移速率矩阵A中，对角线元素总是非零元素)，因此，算法中对应矩阵R的存储数组分别为value_R，colind和rowptr。

算法1　Kernel1算法核心步骤

算法2　Kernel2算法核心步骤

3　算例分析

并行算法实现基于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，持续时间分别为90 s、150 s和90 s。假设所有设备的平均修复时间(mean time to repair,MTTR)和平均故障间隔时间(mean time between failures,MTBF)均为30 min和300 min。任务T各阶段的系统故障树如图2所示。

任务T中3个阶段参与任务的设备数分别为10、13和16，各阶段矩阵规模分别为1 024×1 024，8 192×8 192和65 536×65 536，非零元素分别为4 275，51 282和514 233。

采用并行UM算法和传统UM算法，可得任务T的可靠性计算时间如表1所示，相应的各阶段任务可靠性计算结果如表2所示。

表1　并行UM和传统UM的可靠性计算时间　　s

由表1可知，在任务T的各阶段计算中， 并行UM比传统UM算法分别能达到1.00，4,32和7.24倍的加速效果。由表1也可以看出，虽然并行UM算法具有一定的加速效果，但并不明显，主要原因一是硬件性能限制；二是由于算法基于CUSPARSE函数库实现，UM并行计算效率依赖于所调用库函数；同时，由于CUSPARSE库中并没有直接实现UM算法的函数，导致Kernel2算法过程不够简洁高效。比如，由Kernel2算法步骤可以看出，Step8和Step9中两次调用了函数cusparseDcsrmv()，因为每次迭代中y(k)必须单独保存，但在Step8中，y(k)作为中间计算结果无法保存，因此，Step9中再次调用函数cusparseDcsrmv()，单独计算并保存y(k)。

表2　并行UM和传统UM的可靠性计算结果

同时，由表2可知，并行UM与传统UM算法的计算结果完全一致。这是因为，并行UM算法只是将部分运算转移至GPU中进行，并没有改变算法具体计算过程和步骤。同时也可看出，利用式(7)，能够有效控制算法的计算精度。

4　结论

针对大规模条件下PMS任务可靠性求解困难的问题，将并行计算思想引入传统UM算法，提出了基于CUSPARSE的并行UM算法。该算法利用Nvidia提出的并行计算平台CUDA，通过调用平台工具库CUSPARSE中的函数，实现了对UM算法的并行计算，以及算法的误差控制。实例分析表明，并行UM算法与传统UM算法相比，其计算结果完全一致，且并行UM计算效率具备优势，但并不明显。这一方面是由于实验平台中GPU的性能限制，另一方面因为调用CUDA函数库导致性能受限。综上所述，基于CUSPARSE的PMS任务可靠性并行算法，比传统算法具有较高的效率，但由于受到调用CUDA库函数的限制，下一步可考虑基于CUDA平台自主实现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]刘斌,武小悦.基于多阶段贝叶斯网络的反导系统任务可靠性建模[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方法的多阶段任务系统可靠性分析综述[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的稀疏矩阵运算优化研究[D].武汉：华中科技大学,2012.

[12] 王迎瑞,任江勇,田荣.基于GPU的高性能稀疏矩阵向量乘及CG求解器优化[J].计算机科学,2013,40(3):46-49.

[13] 柳有权,尹康学,吴恩华.大规模稀疏线性方程组的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.