CN105576648B - 一种基于gpu-cpu异构计算平台的静态安全分析双层并行方法 - Google Patents

Abstract

本发明涉及一种基于GPU‑CPU异构计算平台的静态安全分析双层并行方法，包括：将电力系统数据初始化，根据实际需要，设置K个预想故障；提取n个预想故障n∈K，利用GPU为n个预设故障分配多组线程，针对每个预想故障下的线程组进行潮流计算，每组线程均包含多组潮流计算中的一次迭代；并在每次迭代的结束后检查各个算例的收敛情况；如果某个线程组下的潮流计算收敛，则提取一个新的预想故障，执行潮流计算，直到所有故障计算完毕。该算法在单个GPU上实现多个预想事故的同步计算，能够将消去树同一层中的节点对应的计算任务并行处理；从而避免GPU计算资源的浪费，提高了GPU的利用效率。

Description

一种基于GPU-CPU异构计算平台的静态安全分析双层并行 方法

技术领域

本发明涉及电网仿真领域，具体涉及一种基于GPU-CPU异构计算平台的静态安全分析双层并行方法。

背景技术

静态安全分析是电力系统在线分析中的基本工具，需要检验系统出现故障后系统的运行状态。在电力系统中，为避免过负荷和电压越界引起的设备损坏，或由于过负荷设备在系统保护作用下退出运行而导致大面积连锁反应性的停电，在线地进行系统静态安全分析非常重要。静态安全分析实质上是电力系统运行的稳态分析过程，即潮流问题,不涉及元件动态特性和电力系统动态过程。静态安全在线分析对对电网元件逐个设置故障，然后进行潮流分析。由于电网元件数量巨大，其导致包含巨大的计算量。潮流计算常用PQ潮流算法和牛顿拉夫系方法。PQ法计算速度快，收敛性并不好。牛顿拉夫逊方法需要反复求解一系列的线性方程组,具有更好的鲁棒性和可预测性。

传统的在线分析通过CPU机群并行来实现，需要耗费大量的电力。相比传统CPU，GPU具备强大的浮点计算能力，却耗费更少的电力。作为加速器件，它已经越来越普遍地运用于超级计算机中。在电力系统研究领域，一些基于GPU的暂态稳定求解器以及潮流求解器已经开展。对于如何充分考虑计算过程中矩阵运算的稀疏性，充分利用GPU的计算能力，是目前研究的重点。

发明内容

为了弥补上述缺陷，本发明提供一种基于GPU-CPU异构计算平台的静态安全分析双层并行方法，该方法充分利用GPU的计算能力将多个预想事故计算中雅克比矩阵对应的消去树合并，通过这个合并后的消去树指导多个算例的粗粒度并行计算。

为了实现上述发明目的，本发明采取如下技术方案：

一种基于GPU-CPU异构计算平台的静态安全分析双层并行方法，其特征在于，所述方法包括：

(1)将电力系统数据初始化，根据实际需要，设置K个预想故障；

(2)提取n个预想故障n∈K，利用GPU为n个预设故障分配多组线程,

(3)针对每个预想故障下的线程组进行潮流计算，每组线程均包含多组潮流计算中的一次迭代；并在每次迭代的结束后检查各个算例的收敛情况；

(4)如果某个线程组下的潮流计算收敛，则提取一个新的预想故障，并根据步骤(3)执行潮流计算,直到所有预想故障计算完毕。

优选的，所述步骤(3)的潮流计算具体包括下

述步骤：

(3.1)将电力系统初始化，构建电网导纳矩阵Y；

(3.2)根据导纳矩阵Y的电网电压、功率初始值，获取雅克比矩阵的功率差额向量ΔW；

(3.3)利用三分支并行分解算法求解雅克比矩阵；

(3.4)获取电压修正量ΔV；

(3.5)确定收敛情况，若ΔV小于设定值，则结束本次计算任务；否则，返回步骤(3.2)。

进一步地，所述步骤(3.1)的初始化包括：获取电网数据，建立相应的电网导纳矩阵Y，所述初始化仅在每个算例开始计算时执行一次；

所述步骤(3.2)获取雅克比矩阵为

简记为：

JΔV＝ΔW。

其中，J为雅克比矩阵，ΔW为功率差额向量，ΔV为电压修正量。

进一步地，所述步骤(3.2)获取雅克比矩阵的功率差额向量ΔW包括：

对于PQ节点，有功功率和无功功率为给定值，节点电压和相位为待求量；则对应的功率差额表达式为：

其中，ΔP_i为节点有功功率差额；ΔQ_i为节点无功功率差额；

对应的雅克比矩阵J和电压向量ΔV_i的表达式为：

对于PV节点，其有功功率和电压幅值为给定值，节点的无功功率和电压相位为待求量，则功率差额表达式为：

对应的雅克比矩阵和待求的电压向量为：

结合式(5)和式(9)，获得H、N、M、L、R、S的表达式：

其中，P为电力系统节点有功，Q为电力系统节点无功，e_i和f_i分别为节点电压V_i的实部与虚部；P_is和Q_is分别为节点有功与无功的初始值，V_is为节点电压的初始值；G和B分别为电网导纳矩阵Y的实部与虚部，G_ij和B_ij的下标表示节点i到j的导纳。

进一步地，所述步骤(3.3)中利用三分支并行分解算法求解雅克比矩阵包括：对雅克比矩阵J进行LU分解，并对消去树分层，获取导纳矩阵中J矩阵的依赖关系。

进一步地，所述LU分解的方法包括，将矩阵J分解为矩阵L和U的乘积，L为下三角矩阵，U为上三角矩阵；

当第i次迭代时，在矩阵L上对每一个j∈{x|i+1≤x≤n}，有

在L矩阵上，对于矩阵对角元素，有

在U矩阵上，对每一个j∈{x|i+1≤x≤n}有

最后，对U矩阵进行标幺化

U_ij＝U_ij/L_ii (13)

其中，L_ij为L矩阵上第i行第j列的元素，U_ij为U矩阵上第i行第j列的元素，J_ij代表J矩阵上第i行第j列的元素。

进一步地，所述对雅克比矩阵LU分解之前还包括：将J矩阵的下三角部分复制入L矩阵，上三角部分复制入U矩阵，并将U矩阵对角元均置为1。

进一步地，所述步骤(3.3)中，对消去树分层包括，采用消去树理论，将导纳矩阵分层；

所述潮流并行计算层，用于消去树的导纳矩阵并行计算。

进一步地，所述获得导纳矩阵中J矩阵的依赖关系包括：引入消去树结构，分析雅克比矩阵的消去顺序；

对于结构对称的n×n矩阵J，其对应的消去树上的节点j的父节点p_j满足下式：

p_j＝min{i|L_ij≠0,1≤j＜i≤n} (14)

其中，L_ij为L矩阵上第i行第j列的元素。

进一步地，所述消去树的性质包括：

性质1：若i>j，则第i行或列的元素在消去过程中受第i行或者第j列的元素直接影响的充要条件为L_ij≠0；

将消去树进行分层，使得消去树指导矩阵并行分解；对于没有子节点的节点，将其层号定义为0；对于其他节点i，其层号d_i定义为：

d_i＝max{d_j|j∈K}+1 (15)

集合K由节点i的所有子节点构成：

K＝{j|p_j＝i,1≤j＜i≤n} (16)

将节点i的消去过程定义为均按照式(1)～式(4)进行计算，并结合性质1和分层方法获得如下性质：

性质2：消去树第m层所有节点可以被消去的充分条件为第0层到第m-1层的所有节点均被消去；

性质3：消去树同一层中的所有节点消去过程之间没有依赖。

进一步地，所述根据GPU线性方程组的分层并行算法潮流任务分配方法包括：结合最小度原则和最小层原则提出基于MD-ML和改进ML-MD算法的混合算法定义轮数：

对于采用最小度作为优先编号判据的算法，对采用同样出线度上限作为限制条件的所有节点编号过程定义为一轮；

对于采用最小层作为优先编号判据的算法，对采用相同层数作为限制条件的所有节点编号过程定义为一轮，具体包括：

设置动态参数d限制每轮编号过程中的节点出线度的最大值；

假设节点编号的轮数为R，令

d＝R (17)

或者，按分段函数的形式设置

在编号的初始阶段，采用MD-ML编号方法；在出线度为1的节点编号完成之后，改用改进的ML-MD编号方法进行编号。

进一步地，所述步骤(3.4)中电压修正量ΔV的获取如下式所示：

Lx＝-ΔW

UΔV＝x

其中，U为LU分解中的上三角矩阵，x为中间变量。

与最接近的现有技术相比，本发明达到的有益效果是：

1、使用三分支算法，可尽可能多地提升稀疏线性方程组LU分解时的GPU线程数量，提高任务并发度；

2、改进的节点编号算法，减少编号后矩阵对应的消去树的层数。由于使用三分支并行的算法计算时间和树的高度相关性很大，减少树的高度能够直接减少计算时间。

3、使用分层后的消去树，指导前代回代的并行计算。一方面加速了前代回代的计算速度；另一方面将计算数据维持在GPU上，避免了CPU和GPU之间的数据传输。

4、使用树合并技术，进一步利用消去树指导多个算例的计算任务安排，进行静态安全分析的计算。

5、采用基于GPU线程组的粗粒度并行技术，在单个GPU上实现多个预想事故的同时分析计算，避免GPU计算资源的浪费，提高了GPU的利用效率。

附图说明

图1为本发明提供的基于GPU-CPU异构计算平台的静态安全分析双层并行结构示意图；

图2为本发明提供的基于行列消去的LU分解单次迭代的示意图；

图3为本发明提供的分层消去树示意图；

图4为本发明提供的并行三分支LU分解示意图；

图5为本发明提供的一次迭代计算方法流程图。

具体实施方式

下面结合附图对本发明的具体实施方式做进一步的详细说明：

如图1所示，一种基于GPU‐CPU异构计算平台的静态安全分析双层并行方法，

GPU一类专门用于处理图形数据，加速图像输出的电子电路。近些年，将GPU作为一种面向通用计算的流处理器的做法也越来越普遍。GPU可以用于多种并行计算任务，比如分子动力学计算。他们非常适合数据输入输出量非常大的计算。大量的数据使得GPU可以充分地利用GPU的向量计算单元或者单指令多数据的结构。基于GPU的计算在大规模的计算中发挥了越来越大的作用，世界上最强的十台超级计算机中，有三台都利用了GPU的优势。GPU通过硬件实现了大量可以快速切换的线程，可供并行计算使用。

所述方法包括：

(1)将电力系统数据初始化，根据实际需要，设置K个预想故障；

(2)提取n个预想故障n∈K，利用GPU为n个预设故障分配多组线程；图1中，线程组1在第m次迭代时完成了考虑故障1的潮流计算，之后CPU为其装载了考虑故障5的潮流算例。在第m+1次迭代中，线程组1就开始进行故障5的计算。

(3)针对每个预想故障下的线程组进行潮流计算，每组线程均包含多组潮流计算中的一次迭代；如图5所示，并在每次迭代的结束后检查各个算例的收敛情况；

步骤(3)的潮流计算具体包括下述步骤：

(3.1)将电力系统初始化，构建导纳矩阵Y；

步骤(3.1)的初始化包括：获取电网数据，建立相应的电网导纳矩阵Y，所述初始化仅在每个算例开始计算时执行一次；

雅克比矩阵为：

简记为：

JΔV＝ΔW。

其中，J为雅克比矩阵，ΔW为功率差额向量，ΔV为电压修正量。

(3.2)根据导纳矩阵Y的电网电压、功率初始值，获取雅克比矩阵的功率差额向量ΔW；

对于PQ节点，有功功率和无功功率为给定值，节点电压和相位为待求量；则对应的功率差额表达式为：

其中，ΔP_i为节点有功功率差额；ΔQ_i为节点无功功率差额；

对应的雅克比矩阵J和电压向量ΔV_i的表达式为：

对于PV节点，其有功功率和电压幅值为给定值，节点的无功功率和电压相位为待求量，则功率差额表达式为：

对应的雅克比矩阵和待求的电压向量为：

结合式(5)和式(9)，获得H、N、M、L、R、S的表达式：

其中，P为电力系统节点有功，Q为电力系统节点无功，e_i和f_i分别为节点电压V_i的实部与虚部；P_is和Q_is分别为节点有功与无功的初始值，V_is为节点电压的初始值；G和B分别为导纳矩阵Y的实部与虚部，G_ij和B_ij的下标表示节点i到j的导纳。

(3.3)利用三分支并行分解算法求解雅克比矩阵；

步骤(3.3)中利用三分支并行分解算法求解雅克比矩阵包括：对雅克比矩阵J进行LU分解，并对消去树分层，获取线性方程组中J矩阵的依赖关系。

并行计算是一种可以同时进行多项任务的计算。其基本理念在于大的问题常常可以划分为小的问题，而这些小的问题通常可以同时解决。并行计算有多种形式：指令级并行，数据和任务并行等等。并行计算在高性能计算领域已经有多年的应用历史。在民用领域，由于单核处理器性能的限制，并行计算的重要性也越来越得到重视。

所述LU分解的方法包括，将矩阵J分解为矩阵L和U的乘积，L为下三角矩阵，U为上三角矩阵；

当第i次迭代时，在矩阵L上对每一个j∈{x|i+1≤x≤n}，有

在L矩阵上，对于矩阵对角元素，有

在U矩阵上，对每一个j∈{x|i+1≤x≤n}有

最后，对U矩阵进行标幺化

U_ij＝U_ij/L_ii (13)

其中，L_ij为L矩阵上第i行第j列的元素，U_ij为U矩阵上第i行第j列的元素，J_ij代表J矩阵上第i行第j列的元素。

雅克比矩阵LU分解之前需要进行的工作还包括：将J矩阵的下三角部分复制入L矩阵，上三角部分复制入U矩阵，并将U矩阵对角元均置为1。

如图2所示为分解过程中单次迭代的示意图。其中，点状图案覆盖区域对应的是L矩阵中的一列，是公式(1)生成的，竖线覆盖的方块对应的是对角元，是公式(2)生成的，横杠图案覆盖区域对应的是U矩阵中的一行，是公式(3)生成的。其在计算L矩阵的第i列和U矩阵的第i行时，并不修改A_mn(m>i,n>i)的值。

如图3所示，右侧所示为一个结构对称的11×11矩阵A，黑色圆圈代表的是矩阵A中的非零元，而十字叉是A在消去过程中的新增注入元。图3左侧所示为右侧矩阵对应的分层消去树。

根据消去树以及基本的分解方法，可以得到适于GPU的并行三分支LU分解方法

如图4所示为并行三分支LU分解示意图。图中十字叉代表本层节点对应的L矩阵元素，灰色正方形代表本层节点对应的对角元，加号代表本层节点对应的U矩阵元素。空心圆圈代表着已消去的元素。以计算第1层节点的计算过程为例。该层节点包含两个节点：节点3和节点6，计算时需要给这两个节点分配6个线程，分别计算L矩阵的第3列和第6列，第3个和第6个对角元，U矩阵的第3行和第6行。

对消去树分层包括，采用消去树理论，将导纳矩阵分层；

所述潮流并行计算层，用于消去树的导纳矩阵并行计算。

获得导纳矩阵中J矩阵的依赖关系包括：引入消去树结构，分析雅克比矩阵的消去顺序；

对于结构对称的n×n矩阵J，其对应的消去树上的节点j的父节点p_j满足下式：

p_j＝min{i|L_ij≠0,1≤j＜i≤n} (14)

其中，L_ij为L矩阵上第i行第j列的元素。

消去树的性质包括：

性质1：若i>j，则第i行或列的元素在消去过程中受第i行或者第j列的元素直接影响的充要条件为L_ij≠0；

将消去树进行分层，使得消去树指导矩阵并行分解；对于没有子节点的节点，将其层号定义为0；对于其他节点i，其层号d_i定义为：

d_i＝max{d_j|j∈K}+1 (15)

集合K由节点i的所有子节点构成：

K＝{j|p_j＝i,1≤j＜i≤n} (16)

将节点i的消去过程定义为均按照式(1)～式(4)进行计算，并结合性质1和分层方法获得如下性质：

性质2：消去树第m层所有节点可以被消去的充分条件为第0层到第m-1层的所有节点均被消去；

性质3：消去树同一层中的所有节点消去过程之间没有依赖。

根据GPU线性方程组分层并行算法的潮流任务分配方法包括：结合最小度原则和最小层原则提出基于MD-ML和改进ML-MD算法的混合算法定义轮数：

对于采用最小度作为优先编号判据的算法，对采用同样出线度上限作为限制条件的所有节点编号过程定义为一轮；

对于采用最小层作为优先编号判据的算法，对采用相同层数作为限制条件的所有节点编号过程定义为一轮，具体包括：

设置动态参数d限制每轮编号过程中的节点出线度的最大值；

假设节点编号的轮数为R，令

d＝R (17)

以d＝R形式为例，该混合编号算法可用以下流程描述：

或者，按分段函数的形式设置：

在编号的初始阶段，采用MD-ML编号方法；在出线度为1的节点编号完成之后，即第1轮编号完成之后，改用改进的ML-MD编号方法进行编号。

(3.4)获取电压修正量ΔV；

步骤(3.4)中电压修正量ΔV的获取如下式所示：

Lx＝-ΔW

UΔV＝x

其中，U为LU分解中的上三角矩阵，x为中间变量。

(3.5)确定收敛情况，若ΔV小于设定值，则结束本次计算任务；否则，返回步骤(3.2)。

(4)如果某个线程组下的潮流计算收敛，则提取一个新的预想故障，并根据步骤(3)执行潮流计算。

最后应当说明的是：以上实施例仅用以说明本申请的技术方案而非对其保护范围的限制，尽管参照上述实施例对本申请进行了详细的说明，所属领域的普通技术人员应当理解：本领域技术人员阅读本申请后依然可对申请的具体实施方式进行种种变更、修改或者等同替换，这些变更、修改或者等同替换，其均在其申请待批的权利要求范围之内。

Claims (11)

1.一种基于GPU-CPU异构计算平台的静态安全分析双层并行方法，其特征在于，所述方法包括：

(1)将电力系统数据初始化，根据实际需要，设置K个预想故障；

(2)提取n个预想故障n∈K，利用GPU为n个预设故障分配多组线程；

(3)针对每个预想故障下的线程组进行潮流计算，每组线程均包含多组潮流计算中的一次迭代；并在每次迭代的结束后检查各个算例的收敛情况；

(4)如果某个线程组下的潮流计算收敛，则提取一个新的预想故障，并根据步骤(3)执行潮流计算。

2.如权利要求1所述的方法，其特征在于，所述步骤(3)的潮流计算具体包括下述步骤：

(3.1)将电力系统初始化，构建电网导纳矩阵Y；

(3.2)根据导纳矩阵Y的电网电压、功率初始值，获取雅克比矩阵的功率差额向量ΔW；

(3.3)利用三分支并行分解算法求解雅克比矩阵；

(3.4)获取电压修正量ΔV；

(3.5)确定收敛情况，若ΔV小于设定值，则结束本次计算任务；否则，返回步骤(3.2)。

3.如权利要求2所述的方法，其特征在于，所述步骤(3.1)的初始化包括：获取电网数据，建立相应的电网导纳矩阵Y，所述初始化仅在每个算例开始计算时执行一次。

4.如权利要求2所述的方法，其特征在于，所述步骤(3.2)获取雅克比矩阵为：

简记为：

JΔV＝ΔW；

其中，J为雅克比矩阵，ΔW为功率差额向量，ΔV为电压修正量；

功率差额向量ΔW包括：

对于PQ节点，有功功率和无功功率为给定值，节点电压和相位为待求量；则对应的功率差额表达式为：

其中，ΔP_i为节点有功功率差额；ΔQ_i为节点无功功率差额；

对应的雅克比矩阵J和电压向量ΔV_i的表达式为：

对于PV节点，其有功功率和电压幅值为给定值，节点的无功功率和电压相位为待求量，则功率差额表达式为：

对应的雅克比矩阵和待求的电压向量为：

结合式(5)和式(9)，获得H、N、M、L、R、S的表达式：

其中，P为电力系统节点有功，Q为电力系统节点无功，e_i和f_i分别为节点电压V_i的实部与虚部；P_is和Q_is分别为节点有功与无功的初始值，V_is为节点电压的初始值；G和B分别为电网导纳矩阵Y的实部与虚部，G_ij和B_ij表示节点i到j的导纳。

5.如权利要求4所述的方法，其特征在于，所述步骤(3.3)中利用三分支并行分解算法求解雅克比矩阵包括：对雅克比矩阵J进行LU分解，并对消去树分层，获取电网导纳矩阵中J矩阵的依赖关系。

6.如权利要求5所述的方法，其特征在于，所述LU分解的方法包括，将矩阵J分解为矩阵L和U的乘积，L为下三角矩阵，U为上三角矩阵；

当第i次迭代时，在矩阵L上对每一个j∈{x|i+1≤x≤n}，有

在L矩阵上，对于矩阵对角元素，有

在U矩阵上，对每一个j∈{x|i+1≤x≤n}有

最后，对U矩阵进行标幺化

U_ij＝U_ij/L_ii (13)

其中，L_ji为L矩阵上第j行第i列的元素，U_ij为U矩阵上第i行第j列的元素，J_ji代表J矩阵上第j行第i列的元素。

7.如权利要求5所述的方法，其特征在于，所述对雅克比矩阵LU分解之前还包括：将J矩阵的下三角部分复制入L矩阵，上三角部分复制入U矩阵，并将U矩阵对角元均置为1。

8.如权利要求5所述的方法，其特征在于，所述步骤(3.3)中，对消去树分层包括，采用消去树理论，将导纳矩阵分层；

潮流并行计算层，用于消去树的导纳矩阵并行计算。

9.如权利要求5所述的方法，其特征在于，所述获得导纳矩阵中J矩阵的依赖关系包括：引入消去树结构，分析雅克比矩阵的消去顺序；

对于结构对称的n×n矩阵J，其对应的消去树上的节点j的父节点p_j满足下式：

p_j＝min{i|L_ij≠0,1≤j＜i≤n} (14)

其中，L_ij为L矩阵上第i行第j列的元素。

10.如权利要求6所述的方法，其特征在于，所述消去树的性质包括：

性质1：若i>j，则第i行或列的元素在消去过程中受第i行或者第j列的元素直接影响的充要条件为L_ij≠0；

将消去树进行分层，使得消去树指导矩阵并行分解；对于没有子节点的节点，将其层号定义为0；对于其他节点i，其层号d_i定义为：

d_i＝max{d_j|j∈K}+1 (15)

集合K由节点i的所有子节点构成：

K＝{j|p_j＝i,1≤j＜i≤n} (16)

将节点i的消去过程定义为均按照式(10)～式(13)进行计算，并结合性质1和分层方法获得如下性质：

性质2：消去树第m层所有节点可以被消去的充分条件为第0层到第m-1层的所有节点均被消去；

性质3：消去树同一层中的所有节点消去过程之间没有依赖。

11.如权利要求6所述的方法，其特征在于，所述步骤(3.4)中电压修正量ΔV的获取如下式所示：

Lx＝-ΔW

UΔV＝x

其中，x为中间变量。

Images