一种基于gpu的多波前潮流计算方法和系统的制作方法
【专利摘要】本发明公开了一种基于GPU的多波前潮流计算方法，该方法包括：1)波前矩阵链异步并发执行；2)CPU-GPU间任务的分配；3)GPU上矩阵乘法、减法算法的优化。波前矩阵链采用异步并发执行的模式处理，GPU的资源得到了充分的利用；根据矩阵的大小决定采用CPU还是GPU处理，使得单个矩阵的处理时间最少，结合波前矩阵链的异步并发执行，能尽可能的减少CPU、GPU的空闲等待时间；在GPU上对矩阵乘法、减法进行优化，并适时的利用shared memory，使程序的性能达到最优。上述3种方法的结合，可以显著地提升分解矩阵的性能。
【专利说明】一种基于GPU的多波前潮流计算方法和系统
【技术领域】
[0001] 本发明属于计算机高性能计算领域，更具体地，涉及一种基于GPU的多波前潮流 计算方法和系统，该方法能够在GPU上利用多波前算法高效求解潮流计算。
【背景技术】
[0002] 如今，各行各业都已离不开电，在一个大城市中，一次意外的断电事故将会造成巨 大的损失，因此，对电网能否安全稳定的运行进行预测显得尤为重要。电力系统潮流计算是 对电力系统进行安全分析的基础，其实质是求解一组多元非线性方程组。牛顿-拉夫逊法 因收敛速度快的特性，是求解潮流方程最常用的方法，它通过引入雅可比矩阵将非线性方 程组的求解转变为相应的线性方程组的求解。由于雅可比矩阵非常的稀疏，当电网的规模 非常大时，如何对雅可比矩阵进行快速的分解是一个具有挑战性的研宄课题。
[0003] 目前，求解上述线性方程组主要有两种方案：⑴稀疏三角分解技术：利用电力系 统中方程组的稀疏特性，尽量减少不必要的计算以提高求解的效率，但该方案不太易并行 化；（2)多波前方法：将大的稀疏的雅可比矩阵转化为一系列小的稠密矩阵（波前矩阵）， 然后处理这些稠密矩阵，该方案可以很好的并行化。由于稀疏三角分解技术的原理较简单， 如今还是稀疏三角分析技术应用的比较多，因为不易并行化的特点，也就无法与GPU很好 的结合。而多波前方法因并行性好，可以很好的发挥GPU的计算能力，但目前研宄的还比较 少。
【发明内容】
[0004] 针对现有技术的以上缺陷或改进需求，本发明提供了一种基于GPU的多波前潮流 计算方法和系统，其目的在于解决现有方法中存在的并行性程度不高及GPU资源没充分利 用的技术问题。
[0005] 为实现上述目的，按照本发明的一个方面，提供了一种基于GPU的多波前潮流计 算方法，包括以下步骤：
[0006] (I)CPU端读取潮流计算的原始数据，对其进行处理得到线性方程组的系数矩阵J 和常数项b，并构建波前矩阵链的拓扑关系图TG;
[0007] (2)在CPU和GPU上分配计算资源处理线性方程组J·ΔX=b:在CPU上对波前 矩阵作预处理；在GPU上作矩阵-矩阵乘法、减法的数值运算；
[0008] (3)根据TG为每个CPU端线程分配一个待处理的波前矩阵链，之后每个CPU端线 程开始处理其对应波前矩阵链中的波前矩阵；(4)-(9)表示第cthread-i个CPU线程处理 波前矩阵链的过程，CPU中各个线程并行执行（4)-(9)步；
[0009] (4)cthread-i号线程判断波前矩阵链中的波前矩阵是否都已被处理，如果没有则 转入步骤（6);否则转入步骤（5);
[0010] (5)cthread-i号线程判断TG中的所有波前矩阵链是否都已被处理，如果是则转 入步骤（10);否则根据TG判断目前是否有可以被处理的波前矩阵链j，如果有则将波前矩 阵链j分配给cthread-i，然后进入步骤（6)，否则再次转（5)，进行循环查询；
[0011] (6)cthread-i号线程从为其分配的波前矩阵链中取一个波前矩阵F处理，将F分 为4块小矩阵：13、4、0、(3，其维度分别为1^*1^、1^*11、111*1^、111*11，1^、111、11为任意值，？的分解过程 如下所5，
[0012]
【权利要求】
1. 一种基于GPU的多波前潮流计算方法，其特征在于，所述方法包括以下步骤： (1) CPU端读取潮流计算的原始数据，对其进行处理得到线性方程组的系数矩阵J和常 数项b，并构建波前矩阵链的拓扑关系图TG; (2) 在CPU和GPU上分配计算资源处理线性方程组J?AX=b:在CPU上对波前矩阵 作预处理；在GPU上作矩阵-矩阵乘法、减法的数值运算； (3) 根据TG为每个CPU端线程分配一个待处理的波前矩阵链，之后每个CPU端线程开 始处理其对应波前矩阵链中的波前矩阵；(4)-(9)表示第cthread-i个CPU线程处理波前 矩阵链的过程，CPU中各个线程并行执行（4)-(9)步； (4) cthread-i号线程判断波前矩阵链中的波前矩阵是否都已被处理，如果没有则转入 步骤（6);否则转入步骤（5); (5)Cthread-i号线程判断TG中的所有波前矩阵链是否都已被处理，如果是则转入步 骤（10);否则根据TG判断目前是否有可以被处理的波前矩阵链j，如果有则将波前矩阵链 j分配给cthread-i，然后进入步骤（6)，否则再次转（5)，进行循环查询； (6) cthread-i号线程从为其分配的波前矩阵链中取一个波前矩阵F处理，将F分为4 块小矩阵：13、4、0、(]，其维度分别为1^*1^、1^*11、111*1^、111*11，1^、111、11为任意值，？的分解过程如下 所示：
(7) 确定计算量为Cal=m*k*n,在cthread-i号线程中判断Cal与阈值Threshold的 大小关系；若Cal〈 =Threshold,则转入步骤（8)，否则转入步骤（9); (8) cthread-i号线程执行矩阵L2、U2的乘法运算，利用得到的结果与矩阵C做减法运 算，再将结果存储在矩阵C中，之后转入步骤（4); (9) cthread-i号线程将矩阵L2、U2、C存放在SQ的第sq-i个存储单元中，然后GPU中 的第gblock-i个线程块从第sq-i个存储单元中取数据执行，处理完成后将结果传回给CPU 端的第cthread-i个线程；之后转入步骤（4); (10) 终止GPU程序的运行； (11) 利用分解得到的矩阵L、U进行前代回代计算求得AX，然后利用AX去更新X，之 后判断AX中绝对值最大的|AX」是否小于10,如果小于则转入步骤（12);否则转入步 骤⑶； (12) 根据所求得的X去更新潮流计算中的功率变量。
2. 根据权利要求1所述的多波前潮流计算方法，其特征在于，所述步骤（1)包括以下子 步骤： (1-1)读取潮流计算的原始数据，根据电路理论中的基尔霍夫第一定律得到非线性方 程组I=YV，I、V、Y分别为电流、电压、导纳矩阵； (1-2)利用牛顿-拉夫逊法通过求导将非线性方程组转化为线性方程组，转化的公式
(1-3)对J作行列变换，之后根据多波前算法理论对J进行符号分解，得到消去树T、每 个波前矩阵的信息、波前矩阵链的结构及波前矩阵之间的依赖关系； (1-4)根据波前矩阵链的结构、波前矩阵之间的依赖关系得到波前矩阵链的拓扑关系 图TG。
3. 根据权利要求1或2所述的多波前潮流计算方法，其特征在于，所述步骤（2)具体包 括： (2-1)在CPU端分配n个线程，每个线程处理一个波前矩阵链； (2-2)在GPU端分配n个线程块，并在GPU显存的全局内存上分配大小为n*unit的存 储队列SQ，其中每段的长度为unit，共有n段，n段是连续的； (2-3) CPU中的第cthread-i个线程与GPU中的第gblock-i个线程块、存储队列中的第sq_i个存储单元--对应，其中 1〈 =cthread-i、gblock-i、sq_i〈=n; (2-4)利用GPU显存上的纹理内存处理SQ的同步操作。
4. 根据权利要求1或2所述的多波前潮流计算方法，其特征在于，所述步骤（6)具体包 括： (6-1)利用LU分解法将矩阵B分解为上三角矩阵L1，其维度为k*k，和下三角矩阵U1， 其维度为k*k; (6-2)利用（6-1)得到的1^、％，对矩阵D、A进行矩阵-矩阵乘法运算，得到矩阵L2，其 维度为m*k，以及U2，其维度为k*n; (6-3)利用（6-2)得到的L2、U2，通过矩阵-矩阵的乘法与减法，对矩阵C更新，再将结 果存储在C中，C=C-L2U2。
5. 根据权利要求1或2所述的多波前潮流计算方法，其特征在于，所述步骤（7)中的阈 值Threshold具体为： 设CPU处理矩阵所花费的时间为Tepu，GPU处理矩阵所花费的时间为，下式表示了T^、U具体的值
&=mnk，N2=mk+nk+2mn，1表示做矩阵运算操作的次数，1表示数据传输的大小；aagpu分别为在CPU、GPU上执行矩阵-矩阵乘法、减法操作的平均性能；0为CPU-GPU 间传输矩阵的平均带宽
6. 根据权利要求4所述的多波前潮流计算方法，其特征在于，所述步骤（9)包括以下子 步骤： (9-l)cthread-i号线程判断SQ的sq-i个存储单元是否为空，如不为空则再次转 (9-1)，进行循环查询；否则转（9-2); (9-2)cthread_i号线程将L2、U2、C拷贝到SQ的第sq-i个存储单元中，之后通知GPU的第gblock-i个线程块有数据需要处理； (9-3)当GPU中第gblock-i个线程块的线程处于空闲状态时，判断SQ的第sq-i个存 储单元是否为空，如果为空，则转（9-3)，进行循环查询；否则转（9-4); (9-4)判断Cal的大小，如果Cal>=Gthreshold(阈值），则将第sq-i个存储单元中 的数据转存到GPU的共享内存中，然后执行矩阵乘法、减法运算；否则直接对第sq-i个存储 单元的数据执行矩阵乘法、减法运算；处理完后，将结果传回给CPU端的cthread-i号线程。
7. -种基于GPU的多波前潮流计算系统，其特征在于，包括： 第一模块，用于在CPU端读取潮流计算的原始数据，对其进行处理得到线性方程组的 系数矩阵J和常数项b，并构建波前矩阵链的拓扑关系图TG; 第二模块，用于在CPU和GPU上分配计算资源处理线性方程组J?A X=b :在CPU上 对波前矩阵作预处理；在GPU上作矩阵-矩阵乘法、减法的数值运算； 第三模块，用于根据TG为每个CPU端线程分配一个待处理的波前矩阵链，之后每个CPU端线程开始处理其对应波前矩阵链中的波前矩阵； 第四模块，用于使cthread-i号线程判断波前矩阵链中的波前矩阵是否都已被处理， 如果没有则转入第六模块；否则转入第五模块； 第五模块，用于使cthread-i号线程判断TG中的所有波前矩阵链是否都已被处理，如 果是则转入第十模块；否则根据TG判断目前是否有可以被处理的波前矩阵链j，如果有则 将波前矩阵链j分配给cthread-i，然后进入第六模块，否则再次转第五模块，进行循环查 询； 第六模块，用于使cthread-i号线程从为其分配的波前矩阵链中取一个波前矩阵F处 理，将F分为4块小矩阵：B、A、D、C，其维度分别为k*k、k*n、m*k、m*n，k、m、n为任意值，F 的分解过程如下所示：
第七模块，用于确定计算量为Cal=m*k*n,在cthread-i号线程中判断Cal与阈值Threshold的大小关系；若Cal〈 =Threshold,则转入第八模块，否则转入第九模块； 第八模块，用于使cthread-i号线程执行矩阵L2、U2的乘法运算，利用得到的结果与矩 阵C做减法运算，再将结果存储在矩阵C中，之后转入第四模块； 第九模块，用于使cthread-i号线程将矩阵L2、U2、C存放在SQ的第sq-i个存储单元 中，然后GPU中的第gblock-i个线程块从第sq-i个存储单元中取数据执行，处理完成后将 结果传回给CPU端的第cthread-i个线程；之后转入第四模块； 第十模块，用于终止GPU程序的运行； 第十一模块，用于利用分解得到的矩阵L、U进行前代回代计算求得AX，然后利用AX去更新X，之后判断AX中绝对值最大的|AX」是否小于10,如果小于则转入第十二模 块；否则转入第三模块； 第十二模块，用于根据所求得的X去更新潮流计算中的功率变量。
【文档编号】G06F9/50GK104484234SQ201410670758
【公开日】2015年4月1日 申请日期:2014年11月21日 优先权日:2014年11月21日
【发明者】徐得超, 陈勇, 王伟, 江涵, 郑然
申请人:中国电力科学研究院, 华中科技大学