一种GPU加速的电力潮流雅可比矩阵的LU分解方法与流程

本发明属于电力系统高性能计算应用领域，尤其涉及一种GPU加速的电力潮流雅可比矩阵的LU分解方法。

背景技术：

潮流计算是电力系统中应用最广泛、最基本和最重要的一种电气运算。在电力系统运行方式和规划方案的研究中，都需要进行潮流计算以比较运行方式或规划供电方案的可行性、可靠性和经济性。同时，为了实时监控电力系统的运行状态，也需要进行大量而快速的潮流计算。因此，在系统规划设计和安排系统的运行方式时，采用离线潮流计算；在电力系统运行状态的实时监控中，则采用在线潮流计算。

而实际生产过程中，无论离线潮流和在线潮流计算都对潮流的计算速度有这比较高的要求。在涉及规划设计和安排运行方式的离线潮流中，因设备落地方案等情况复杂，需要仿真运行的种类多，潮流计算量大，单个潮流计算时间影响整体仿真时长；而在电力系统运行中进行的在线潮流计算对计算时间敏感度高，需要实时给出潮流计算结果，如在预想事故、设备退出运行对静态安全的影响的潮流计算中，系统需要计算大量预想事故下潮流分布，并实时地做出预想的运行方式调整方案。

传统的牛顿拉夫逊法潮流计算中，修正方程组求解占潮流计算时间的70％，修正方程组求解的计算速度影响程序的整体性能。而随着CPU计算速度提升的放缓，现阶段的单个潮流计算计算时间已经达到一个瓶颈。目前对潮流计算的加速方法主要集中在使用集群和多核服务器对多潮流进行粗粒度加速，实际成产中对单个潮流内部运算加速的研究较少。

GPU是一种众核并行处理器，在处理单元的数量上要远远超过CPU。传统上的GPU只负责图形渲染，而大部分的处理都交给了CPU。现在的GPU已经法阵为一种多核，多线程，具有强大计算能力和极高存储器带宽，可编程的处理器。在通用计算模型下，GPU作为CPU的协处理器工作，通过任务合理分配分解完成高性能计算。

稀疏线性方程组求解计算具有并行性。对方程组系数矩阵进行LU符号分解后，得到下三角变换矩阵L和上三角阵U阵的稀疏结构，根据U阵的稀疏结构，对矩阵各列进行并行化分层。其中每层中的列的计算相互独立，没有依赖关系，天然可以被并行的计算处理，适合GPU加速。因此通过CPU和GPU之间合理的调度可以快速完成方程组系数矩阵进行LU分解，并求解稀疏线性方程组，国内外学者已经开始对GPU进行稀疏线性方程组加速求解的方法进行了研究，但是没有深入的优化线程设计，单纯从计算量的分配上研究计算线程设计，对线程计算方式，数据索引方式没有进行深入研究，无法使程序充分发挥GPU的优势。

因此，亟待解决上述问题。

技术实现要素：

发明目的：针对现有技术的不足，本发明提供一种能大幅减少电力潮流雅可比矩阵LU分解计算时间并能提升潮流计算速度的一种GPU加速的电力潮流雅可比矩阵的LU分解方法。

潮流计算：电力学名词，指在给定电力系统网络拓扑、元件参数和发电、负荷参量条件下，计算有功功率、无功功率及电压在电力网中的分布。

GPU：图形处理器(英语：GraphicsProcessingUnit，缩写：GPU)。

本发明公开了一种GPU加速的电力潮流雅可比矩阵的LU分解方法，所述方法包括：

(1)CPU中对雅可比矩阵J进行LU符号分解，得到下三角变换矩阵L和上三角矩阵U矩阵的稀疏结构，符号分解之后的J的稀疏结构等于L+U；根据U阵的稀疏结构，对矩阵J各列进行并行化分层，并将计算所需数据传输给GPU；

(2)GPU中按层次递增的顺序计算分层LU分解内核函数SparseLU。

其中，所述步骤(1)中，并行化分层将矩阵J的n列归并到MaxLevel层中，属于同一层中的列并行进行LU分解；每层包含的列的数量为Levelnum(k)，k表示层号；存储第k层中所有列号至映射表Map_k。

优选的，所述步骤(2)中，分层LU分解内核函数定义为SparseLU<N_blocks，N_threads>，其线程块大小N_thread固定为128，当对k层进行计算时，线程块数量N_blocks＝Levelnum(k)，总线程数量为：N_blocks×N_threads；按照层次递增的顺序，调用内核函数SparseLU<Levelnum(k)，N_threads>来分解属于第k层的所有列。

进一步，所述内核函数SparseLU<Levelnum(k)，N_threads>的计算流程为：

(2.1)CUDA自动为每个线程分配线程块索引blockID和线程块中的线程索引threadID；

(2.2)将blockID和threadID赋值给变量bid和t，之后通过bid和t来索引bid号线程块中的t号线程；

(2.3)第bid号线程块负责LU分解雅可比矩阵J的第j＝Map_k(bid)列；

(2.4)第bid号线程块中，变量i从1递增到j-1，如果U(i，j)≠0，采用公式J(i+1：n，j)＝J(i+1：n，j)-J(i，j)×L(i+1：n，i)更新雅可比矩阵J的第j列，具体步骤如下：

1)判断线程编号t是否小于n–i，否则线程结束执行；

2)J(t+i+1，j)＝J(t+i+1，j)-J(i，j)×L(t+i+1，i)；

3)t＝t+128，返回1)；

(2.5)计算L阵的第j列变换向量：

采用公式L(j：n，j)＝J(j：n，j)/J(j，j)计算下三角阵L阵的第j列，具体步骤如下：

1)判断线程编号t是否小于n–j+1，否则线程结束执行；

2)L(j+t，j)＝J(j+t，j)/J(j，j)；

3)t＝t+128，返回1)。

(2.6)更新J阵第j列：J(j+1：n，j)＝0。

有益效果：与现有技术比，本发明的有益效果为：首先本发明采用CPU对电力潮流的雅可比矩阵J进行LU符号分解，根据U阵的稀疏格式，可以减少不必要的浮点计算；其次根据U阵的稀疏结构将J阵的各列分到可以并行计算的不同层次，并将分层结果传给GPU；再者GPU中按层次递增的顺序启动分层LU分解内核函数SparseLU；最后本发明利用CPU控制程序的流程并处理基础数据和GPU处理密集的浮点运算相结合的模式提高了电力潮流雅可比矩阵LU分解的效率，解决了电力系统静态安全性分析中潮流计算耗时大的问题。

附图说明：

图1为本发明的实例电网的数据表；

图2为本发明的实例计算时间；

图3为本发明的流程示意图。

具体实施方式：

如图3所示，本发明一种GPU加速的电力潮流雅可比矩阵的LU分解方法，所述方法包括：

(1)CPU中对雅可比矩阵J进行LU符号分解，得到下三角变换矩阵L和上三角矩阵U矩阵的稀疏结构，符号分解之后的J的稀疏结构等于L+U；根据U阵的稀疏结构，对矩阵J各列进行并行化分层。

(2)GPU中按层次递增的顺序启动分层LU分解内核函数SparseLU。

其中LU符号分解原理参见：Algorithm 907:KLU,A Direct Sparse Solver for Circuit Simulation Problems,Timothy A.Davis,EkanathanPalamadai Natarajan,ACM Transactions on Mathematical Software,Vol 37,Issue 6,2010,pp 36:1-36:17。LU分层原理参见：陈德扬，李亚楼，江涵，徐得超.基于道路树分层的大电网潮流并行算法及其GPU优化实现[J].电力系统自动化，2014，38(22)：63-69。

一、CPU中对电力潮流雅可比矩阵J进行LU符号分解方法

首先，在CPU中对雅可比矩阵J进行LU符号分解，得到下三角变换矩阵L和上三角阵U阵的稀疏结构，符号分解之后的J的稀疏结构等于L+U；然后，并行化分层将矩阵J的n列归并到MaxLevel层中，属于同一层中的列并行进行LU分解；每层包含的列的数量为Levelnum(k)，k表示层号；映射表Map_k存储第k层中所有列的具体列号。最后，CPU将GPU计算所需数据传输给GPU，数据包括：雅可比矩阵J，其维度n，上三角阵U阵，层数MaxLevel，每层包含的列数Levelnum以及映射表Map。

二、GPU中按层次递增的顺序启动分层LU分解内核函数SparseLU

分层LU分解内核函数定义为SparseLU<N_blocks，N_threads>，其线程块大小N_thread固定为128，当对k层进行计算时，线程块数量N_blocks＝Levelnum(k)，总线程数量为：N_blocks×N_threads；按照层次递增的顺序，调用内核函数SparseLU<Levelnum(k)，N_threads>来分解属于第k层的所有列。

SparseLU<Levelnum(k)，N_threads>的计算流程为：

(1)CUDA自动为每个线程分配线程块索引blockID和线程块中的线程索引threadID；

(2)将blockID和threadID赋值给变量bid和t，之后通过bid和t来索引bid号线程块中的t号线程；

(3)第bid号线程块负责LU分解雅可比矩阵J的第j＝Map_k(bid)列；

(4)第bid号线程块中，变量i从1递增到j-1，如果U(i，j)≠0，采用公式J(i+1：n，j)＝J(i+1：n，j)-J(i，j)×L(i+1：n，i)更新雅可比矩阵J的第j列，具体步骤如下：

1)判断线程编号t是否小于n–i，否则线程结束执行；

2)J(t+i+1，j)＝J(t+i+1，j)-J(i，j)×L(t+i+1，i)；

3)t＝t+128，返回1)；

(5)计算L阵的第j列变换向量：

采用公式L(j：n，j)＝J(j：n，j)/J(j，j)计算下三角阵L阵的第j列，具体步骤如下：

1)判断线程编号t是否小于n–j+1，否则线程结束执行；

2)L(j+t，j)＝J(j+t，j)/J(j，j)；

3)t＝t+128，返回1)。

(6)更新J阵第j列：J(j+1：n，j)＝0。

本发明所使用的GPU计算平台配备一张TeslaK20CGPU卡和IntelXeonE5-2620CPU，GPU的峰值带宽可达208GB/s，单精度浮点计算量峰值可达3.52Tflops，CPU主频为2GHz。CPU计算平台配备IntelCorei7-3520M2.90GHz的CPU。GPU计算平台上对图1中的五个实例电网的雅可比矩阵进行了测试，图2为对应五个雅可比矩阵LU分解的测试时间。