一种基于GPU加速的直流故障筛选方法技术领域
本发明属于电力系统高性能计算应用领域，涉及电力系统故障筛选，尤其涉及一种
基于GPU加速的直流故障筛选方法。
背景技术
GPU是一种众核并行处理器，在处理单元的数量上要远远超过CPU。传统上的GPU
只负责图形渲染，而大部分的处理都交给了CPU。现在的GPU已经法阵为一种多核，
多线程，具有强大计算能力和极高存储器带宽，可编程的处理器。在通用计算模型下，
GPU作为CPU的协处理器工作，通过任务合理分配分解完成高性能计算。目前GPU已
经成功应用于医学影像，计算流体动力学，环境科学，石油勘测等领域。
在电力系统中，直流潮流计算，交流潮流计算，暂态稳定性分析，大型稀疏线性方
程组求解等都是比较耗时的部分。根据不同线性稀疏系统的特点，对不同的解决方法提
出了基于GPU加速策略，例如，LU分解，共轭梯度算法，雅克比矩阵迭代，在这些算
法中可以达到3到10倍的加速比。
N-1静态安全性分析用来检测电网中单一元件断开时电网的静态安全性。对于一个具有
N个元件的电力系统，严格的安全分析需要N次交流潮流计算，这是非常耗时的。在实
际应用中，只有一部分严重的故障会威胁到电网的安全性，因此，先用直流潮流来筛选
严重的故障集，再用交流潮流来进一步分析。目前广泛使用N-1静态安全性分析来检测
电网中单一元件断开时电网的静态安全性。对于一个具有N个元件的电力系统，在实际
应用中，用直流潮流来筛选一部分可能会威胁到电网安全性的严重故障，CPU处理是比
较耗时的。
发明内容
发明目的：针对现有技术的不足，本发明提供了一种可以大幅度减少直流故障的计
算时间，为在线静态安全性分析提供了可能的一种基于GPU加速的直流故障筛选方法。
技术方案：本发明提出了一种基于GPU加速的直流故障筛选方法。
直流潮流：电力学名词，指在给定电力系统网络拓扑、元件参数和发电、负荷参量
条件下，计算有功功率、电压相角在电力网中的分布。
并行计算：相对于串行运算，是一种一次可执行多个指令的算法，目的是提高计算
速度，及通过扩大问题求解规模，解决大型而复杂的计算问题。
GPU：图形处理器(英语：Graphics Processing Unit，缩写：GPU)。
导纳矩阵：以系统元件的等值导纳为基础所建立的、描述电力网络各节点电压和注
入电流之间关系的矩阵。
本发明一种基于GPU加速的直流故障筛选方法，所述方法包括：
(1)根据不同故障的特点对直流潮流算法进行优化；
(2)CPU读取电网数据，生成节点电纳矩阵B，计算节点电抗矩阵X0；利用直流
潮流算法计算基态电网的节点电压相角θ0，将直流筛选所需的电网基础数据传送给
GPU；
(3)GPU中的内核函数1计算支路开断时的节点电压相角θ1；
(4)GPU中的内核函数2计算支路开断时各支路的有功功率并与有功功率限额比
较，筛选出支路开断故障集S1并回传至CPU中；
(5)GPU中的内核函数3计算发电机开断时的节点电压相角θ2；
(6)GPU中的内核函数4计算发电机开断时各支路的有功功率并与有功功率限额
比较，筛选出发电机开断故障集S2并回传至CPU中。
其中，所述步骤(1)中所述故障包括支路开断故障和发电机开断故障，当支路
断开时，电网电抗矩阵发生变化，各节点注入功率增量ΔP为零；当发电机断开时，
电网电抗矩阵不变，即电抗矩阵增量ΔX＝0。
优选的，所述步骤(2)中所述电网基本数据包括电网节点数N、支路数Nbranches、
每条支路首尾所连节点编号df[Nbranches]和dt[Nbranches]、支路电抗X、支路有功功率限额
Pb[Nbranches]、发电机数量Ngen、发电机所连节点编号dg[Ngen]、发电机所发有功功率
Pg[Ngen]、节点电压相角θ0、基态节点注入有功功率P0、节点电抗矩阵X0。
进一步，所述步骤(3)中内核函数1为Kernel_1<Nblocks，Nthreads>，Kernel_1的线
程块大小Nthreads固定为128，线程块数量Nblocks＝(Nbranches+Nthreads–1)/Nthreads；Kernel_1的
第b个线程计算第b条支路断开时的电压相角向量θ1＝θ0+Δθ，设第b条支路的首尾节
点编号分别为i、j，其中Δθ＝c(θ0i-θ0j)ΧL，式中：XL＝[X1i-X1j,L,Xni-Xnj]T，
Δxij为断开支路的电抗值，Xij为基态的节点电抗矩阵X0的第(i,j)
号元素，θ0i为节点电压相角向量θ0的第i号元素；
Kernel_1<Nblocks,Nthreads>的具体计算方法为：
(3.1)CUDA自动为每个线程分配线程索引threadID作为线程号；
(3.2)将线程号赋值给变量tid，通过变量tid来指定线程；
(3.3)第tid号线程计算第b条支路开断时各节点电压的相角
(3.4)判断线程号变量tid的值是否小于支路数Nbranches，小于继续执行以下步骤，
否则该线程结束运算；
(3.5)将第b条支路首尾两个节点编号分别赋值给变量i，j，即：i←df[b]，j←dt[b]；
(3.6)用公式c＝(Xii+Xjj-2×Xij-Δxij)计算c，如果c为0，置S1[b]为2，若不为0，
c＝1/c；
(3.7)循环变量k从1到N-1循环，计算第b条支路开断时的电压相角向量：
其中，表示第b条支路开断时第k节点的电压相角，(θ0)k表示基态时k节点
的电压相角；
(3.8)Nblocks个线程块中的全部线程计算结束得到节点电压相角θ1和故障集S1。
再者，所述步骤(4)中所述内核函数2为Kernel_2<Nblocks，Nthreads>，Kernel_2的
线程块大小Nthreads固定为128，线程块数量Nblocks＝(Nbranches+Nthreads–1)/Nthreads；
Kernel_2<Nblocks，Nthreads>使用Kernel_1的计算结果θ1和S1，其具体计算方法为：
(4.1)CUDA自动为每个线程分配线程索引threadID作为线程的编号；
(4.2)将线程号赋值给变量tid，通过变量tid来指定线程；
(4.3)设第tid号线程计算第b条支路开断时的各支路的有功功率并与限额比较；
(4.4)判断线程号tid是否小于支路数Nbranches，小于继续执行以下步骤，否则该线
程结束运算；
(4.5)tid号线程中循环变量k从1到Nbranches循环，循环内部：a)第k条支路的
首尾节点编号df[k]和dt[k]分别赋值给变量i，j，即：i←df[k]，j←dt[k]；；b)判断S1[b]
不为2且k不等于b继续执行；否则循环结束。c)用下式计算第k条支路的有功功率：
式中：Pijb为第b条支路开断时第k条支路的有功功率，是内核
函数1计算支路开断时的节点电压相角θ1的第(i，b)和(j，b)元素，即第b条支路开断时
i，j节点电压相角；并将结果与有功功率限额比较，大于限额，置S1[b]为1；
(4.6)Nblocks个线程块中的全部线程计算结束，得到S1并传回CPU处理。
优选的，所述步骤(5)中所述内核函数3为Kernel_3<Nblocks，Nthreads>，
Kernel_3<Nblocks，Nthreads>的线程块大小Nthreads固定为128，线程块数量
Nblocks＝(Ngen+Nthreads–1)/Nthreads；
其中，Kernel_3<Nblocks,Nthreads>的具体计算方法为：
(5.1)CUDA自动为每个线程分配线程索引thread ID作为线程的编号；
(5.2)将线程号赋值给变量tid，通过变量tid来指定线程；
(5.3)第tid号线程计算第b台发电机开断时各节点电压相角
(5.4)判断线程号变量tid的值是否小于发电机数Ngen，小于继续执行以下步骤，
否则该线程结束运算；
(5.5)将第tid台发电机所连节点号赋值给变量j：j＝dg(tid)；
(5.6)循环变量k从1到N-1循环，计算第b台发电机开断时的节点电压的相角向
量：其中，表示j节点发电机开断时k节点电压相角，(θ0)k
表示稳态时k节点电压相角；
(5.7)Nblocks个线程块中的全部线程计算结束得到节点电压相角θ2，供内核函数4
使用。
进一步，所述步骤(6)中所述内核函数4为Kernel_4<Nblocks，Nthreads>，Kernel_4
的线程块大小Nthreads固定为128，线程块数量Nblocks＝(Ngen+Nthreads–1)/Nthreads；
其中，Kernel_4<Nblocks，Nthreads>的具体计算方法为：
(6.1)CUDA自动为每个线程分配线程索引thread ID作为线程的编号；
(6.2)将线程号赋值给变量tid，通过变量tid来指定线程；
(6.3)设第tid号线程计算第b台发电机开断时的各支路的有功功率并与限额比较；
(6.4)判断线程号tid是否小于发电机Ngen，小于继续执行以下步骤，否则该线程
结束运算；
(6.5)第tid号线程中循环变量k从1到Nbranches循环，循环内部：a)第k条支路
的首尾节点编号df[k]和dt[k]分别赋值给变量i，j；b)计算第k条支路的有功功率：
式中：Pijb为第b台发电机开断时第k条支路的有功功率，是内核函数4
计算发电机开断时的节点电压相角θ2的第(i，b)和(j，b)元素，并将结果与限额比较，大
于限额，置S2[b]为1；
(6.6)Nblocks个线程块中的全部线程计算结束，得到S2并传回CPU处理。
有益效果：与现有技术相比，本发明具有以下显著优点：首先该基于GPU加速的
直流故障筛选方法根据不同故障特点对直流潮流算法进行优化，避免了不必要的浮点计
算量；其次使用CPU处理电网基本数据，利用电网数据形成节点电抗矩阵，一定程度
上减少了电网数据在CPU和GPU之间的传输的次数，且效率更高；再者GPU的浮点
计算能力和内存带宽是相应CPU的5到10倍，在处理相同的浮点运算时，GPU的能耗
只有CPU的20％左右，该方法利用CPU控制程序的流程并处理基础数据和GPU处理
密集的浮点运算相结合的模式提高了直流故障筛选方法的计算效率，解决了电力系统静
态安全性分析中故障直流筛选耗时大的难题。
附图说明
图1为本发明基于GPU加速的直流故障筛选方法的算法流程图；
图2为本发明基于GPU加速的直流故障筛选方法中内核函数1的算法流程图；
图3为本发明基于GPU加速的直流故障筛选方法中电网算例基本数据；
图4为本发明基于GPU加速的直流故障筛选方法中GPU与CPU测试结果对比。
具体实施方式
本发明公开了一种基于GPU加速的直流故障筛选方法，所述方法包括：
(3)GPU中的内核函数1计算支路开断时的节点电压相角θ1；
(4)GPU中的内核函数2计算支路开断时各支路的有功功率并与有功功率限额比
较，筛选出支路开断故障集S1并回传至CPU中；
(5)GPU中的内核函数3计算发电机开断时的节点电压相角θ2；
(6)GPU中的内核函数4计算发电机开断时各支路的有功功率并与有功功率限额
比较，筛选出发电机开断故障集S2并回传至CPU中。
一、根据不同故障的特点对直流潮流算法进行优化；简化直流潮流计算公式，避免
不必要的浮点计算量；
(1)直流潮流的计算公式为：
P＝Bn×nθ或者θ＝Xn×nP (1)
其中，
n是电网节点个数；
B是节点电纳矩阵，n*n维矩阵；
P是节点有功注入功率，n维向量；
θ是节点电压角度，n维向量；
θi，θj是节点i和节点j电压相角；
X＝B-1是节点电抗矩阵，n*n维矩阵；
xij表示节点i和节点j之间的支路电抗值；
在基态时，直流潮流计算公式可以表示为：θ0＝X0P0。
当发生故障的时，公式(1)可以表示为：
θ1＝X1P1＝(X0+ΔX)(P0+ΔP)＝θ0+Δθ (3)
Δθ＝ΔXP0+X0ΔP+ΔXΔP (4)
其中，下标0和1分别表示故障前状态和故障后状态；Δ表示故障前到故障后的改
变；
如果节点i和节点j之间增加一条支路，电抗矩阵的增量ΔX可以表示为：
其中，M＝ei-ej，Δxij为增加支路的电抗值，如果是开断支
路，Δxij为负值，xij为基态的节点电抗矩阵的第(i,j)号元素，ei，ej为标准基。
(2)直流开断的算法流程如下：
1)在CPU上只计算一次电网基态下的节点电抗矩阵X0。
2)用式5计算支路开断后的电抗矩阵增量ΔX。
3)用式3和4计算节点电压相角Δθ和θ1，
4)用式2计算各支路有功功率，并检查支路是否过负荷。
当c的分母为0时，说明节点i和j之间支路断开会导致系统解列。
(3)根据支路开断故障和发电机开断故障的不同特点对直流潮流算法进行优
化。当节点i和节点j之间的支路断开时，电网电抗矩阵发生变化，各节点注入功
率增量为零，即ΔP＝0。式4可以简化为：
Δθ＝ΔXP0＝cX0MMTX0P0＝cXLMTθ0＝c(θ0i-θ0j)XL (6)
其中，XL＝X0M＝X0(ei-ej)＝[X1i-X1j,L,Xni-Xnj]T。
当连接于节点j的发电机断开时，电网电抗矩阵不变，即ΔX＝0，节点注入功率增
量ΔP是稀疏向量，只有ΔPj非零。式4可以简化为：
Δθ＝Χ0ΔP＝ΔPjXj (7)
这里，ΔP＝[0…ΔPj…0]T，Xj＝[X1j…Xnj]T，Xj为基态下的节点电抗矩阵X0的
第j列元素，Xij为基态下的节点电抗矩阵X0的第(i,j)号元素，ΔPj为连接于节点j的发
电机的有功功率。
二、算法设计：本算法将程序分为CPU和GPU两个部分。
CPU读取电网数据，生成待分析电网的节点电纳矩阵B，为GPU准备基本数据，
处理GPU的计算结果，并控制整个程序的流程。GPU主要处理密集的浮点运算，共采
用四个内核函数来完成故障筛选功能，具体包括：内核函数Kernel_1计算支路开断时的
节点电压相角θ1；内核函数Kernel_2计算支路开断时各支路的有功功率并与有功功率
限额比较，筛选出支路开断故障集S1；内核函数Kernel_3计算发电机开断时的节点电
压相角θ2；内核函数Kernel_4计算发电机开断时各支路的有功功率并与有功功率限额
比较，筛选出发电机开断故障集S2，具体流程如图1所示。
三、CPU生成基态下的节点电抗矩阵X0，传输基础数据
CPU读取电网数据，生成节点电纳矩阵B，使用KLU库函数计算其逆矩阵节点电
抗矩阵X0。利用公式：θ0＝X0P0，计算基态电网各节点电压相角θ0。将直流筛选所需
的数据传送给GPU，具体包括：开断所需的基础数据：电网节点数N，支路数Nbranches，
每条支路首尾所连节点编号df[Nbranches]和dt[Nbranches]，支路电抗X，支路有功功率限额
Pb[Nbranches]，发电机数量Ngen，发电机所连节点编号dg[Ngen]，发电机所发有功功率
Pg[Ngen]，节点电压相角θ0、基态节点注入有功功率P0、节点电抗矩阵X0。
四、计算支路开断时的节点电压相角θ1
GPU内核函数1中，电网支路数Nbranches，支路所连节点编号向量df[Nbranches]和
dt[Nbranches]，节点电压相角θ1(Nbranches×N维矩阵)，故障集S1[Nbranches]。
(1)内核函数1定义为Kernel_1<Nblocks,Nthreads>。使用一个线程thread计算一条支
路开断，共需要的线程数等于支路数Nbranches，线程块大小Nthreads固定为128，线程块数
量Nblocks为(Nbranches+Nthreads–1)/Nthreads。
(2)Kernel_1<Nblocks,Nthreads>流程：1)CUDA自动为每个线程分配线程索引thread
ID作为线程的编号；2)将线程号赋值给变量tid，通过变量tid来指定线程；3)第tid
号线程计算第b条支路开断时各节点电压的相角4)判断线程号变量tid的值是否
小于支路数Nbranches，小于继续执行以下步骤，否则该线程结束运算；5)将b条支路首
尾两个节点编号分别赋值给变量i，j，即：i←df[b]，j←dt[b]；6)用公式c＝(Xii+Xjj-
2×Xij-Δxij)计算c，如果c为0，置S1[b]为2，若不为0，c＝1/c；7)循环变量k从1到
N-1循环，用下式计算第b条支路开断时的电压相角向量：
其中，表示第b条支路开断时第k节点的电压相角，(θ0)k表示基态时k节点
的电压相角；7)Nblocks个线程块中的全部线程计算结束得到节点电压相角θ1和故障集
S1，供内核函数2使用，如图2所示。
五、计算支路开断时各支路的有功功率并筛选出支路开断故障集S1
GPU内核函数2中，电网支路数为Nbranches，支路所连节点编号向量df[Nbranches]和
dt[Nbranches]，节点电压相角矩阵θ1，节点电压相角θ0，故障集S1。
(1)内核函数2定义为Kernel_2<Nblocks,Nthreads>。使用一个线程thread计算一条支
路开断时的各支路有功功率，共需要的线程数等于支路数Nbranches。线程块大小Nthreads
固定为128，线程块数量Nblocks为(Nbranches+Nthreads–1)/Nthreads。
(2)Kernel_2<Nblocks,Nthreads>流程：使用Kernel_1的计算结果θ1和S1。1)CUDA
自动为每个线程分配线程索引thread ID作为线程的编号；2)将线程号赋值给变量tid，
通过变量tid来指定线程；3)设第tid号线程计算第b条支路开断时的各支路的有功功
率并与限额比较；4)判断线程号tid是否小于支路数Nbranches，小于继续执行以下步骤，
否则该线程结束运算；5)tid号线程中循环变量k从1到Nbranches循环，循环内部：a)
第k条支路的首尾节点编号df[k]和dt[k]分别赋值给变量i，j，即：i←df[k]，j←dt[k]；；
b)判断S1[b]不为2且k不等于b继续执行；否则循环结束。c)用下式计算第k条支
路的有功功率：
式中：为第b条支路开断时第k条支路的有功功率，是内核函数1计算
支路开断时的节点电压相角θ1的第(i，b)和(j，b)元素，即第b条支路开断时i，j节点电
压相角。并将结果与有功功率限额比较，大于限额，置S1[b]为1；6)Nblocks个线程块中
的全部线程计算结束，得到S1并传回CPU处理。
六、计算发电机开断时的节点电压的相角θ2
GPU内核函数3中，电网发电机数为Ngen，发电机所连节点编号dg[Ngen]，节点电
压相角θ2(Ngen×N维矩阵)。
(1)内核函数3定义为Kernel_3<Nblocks,Nthreads>。使用一个线程thread计算一台发
电机开断，共需要的线程数等于支路数Ngen，线程块大小Nthreads固定为128，线程块数
量Nblocks为(Ngen+Nthreads–1)/Nthreads。
(2)Kernel_3<Nblocks,Nthreads>流程：1)CUDA自动为每个线程分配线程索引thread
ID作为线程的编号；2)将线程号赋值给变量tid，通过变量tid来指定线程；3)第tid
号线程计算第b台发电机开断时各节点电压相角4)判断线程号变量tid的值是否
小于发电机数Ngen，小于继续执行以下步骤，否则该线程结束运算；5)将第tid台发电
机所连节点号赋值给变量j：j＝dg(tid)；6)循环变量k从1到N-1循环，用下式计算第b
台发电机开断时的节点电压的相角向量：
其中，表示j节点发电机开断时k节点电压相角，(θ0)k表示稳态时k节点电
压相角；7)Nblocks个线程块中的全部线程计算结束得到节点电压相角θ2，供内核函数4
使用。
七、计算发电机开断时各支路的有功功率并筛选发电机开断故障集S2
GPU内核函数4中，电网发电机数为Ngen，发电机所连节点编号dg[Ngen]，节点电
压相角θ2，故障集向量S2[Ngen]。
(1)内核函数4定义为Kernel_2<Nblocks,Nthreads>。使用一个线程thread计算一台发
电机开断时的各支路有功功率，共需要的线程数等于支路数Ngen。线程块大小Nthreads固
定为128，线程块数量Nblocks为(Ngen+Nthreads–1)/Nthreads。
(2)Kernel_4<Nblocks,Nthreads>流程：使用Kernel_3的计算结果θ2和S2。1)CUDA
自动为每个线程分配线程索引thread ID作为线程的编号；2)将线程号赋值给变量tid，
通过变量tid来指定线程；3)设第tid号线程计算第b台发电机开断时的各支路的有功
功率并与限额比较；4)判断线程号tid是否小于发电机Ngen，小于继续执行以下步骤，
否则该线程结束运算；5)第tid号线程中循环变量k从1到Nbranches循环，循环内部：a)
第k条支路的首尾节点编号df[k]和dt[k]分别赋值给变量i，j；b)用下式计算第k条支
路的有功功率：
式中：为第b台发电机开断时第k条支路的有功功率，是内核函数3计
算发电机开断时的节点电压相角θ2的第(i，b)和(j，b)元素。并将结果与限额比较，大于
限额，置S2[b]为1；6)Nblocks个线程块中的全部线程计算结束，得到S2并传回CPU处
理。
八、GPU将计算结果传回CPU进一步处理
GPU将内核函数2和4计算出来的故障集S1和S2传回CPU。若S1[b]为1，则第
b条支路断开属于故障集，若S1[b]为2，则第b条支路断开导致系统解列；若S2[b]为1，
则第b台发电机断开属于故障集。
本发明所使用的GPU计算平台配备一张Tesla K20C GPU卡和Intel Xeon
E5-2620CPU，GPU的峰值带宽可达208GB/s，单精度浮点计算量峰值可达3.52Tflops，
CPU主频为2GHz。CPU计算平台配备Intel Core i7-3520M 2.90GHz的CPU。在CPU
和GPU计算平台上分别对3个电网实例进行了测试，具体实例如图3所示。整个程序
的计算时间主要分为数据传输和核函数执行两个部分，其中数据传输超过总时间的
30％。对于规模比较小的电网，GPU内核函数节约的时间不能够弥补数据传输所耗费的
时间，因此没有加速比，如图4所示，算例1中CPU计算时间是7ms，GPU计算时间
也是7ms，但是这其中有4.76ms是数据传输所耗费的，加速比为1；算例2中CPU计
算时间为260ms，GPU计算时间为95ms，其中44.85ms用在数据传输，加速比为2.73。
当电网规模比较大时，就能达到理想的加速比，算例3中CPU计算时间为3621ms，GPU
计算时间是733ms，其中242.36ms是数据传输所耗费的，占总时间比例相对小，加速
比为4.94。虽然CPU与GPU结合增加了数据传输这部分时间，但是当电网规模比较大
的时候，加速比还是很可观的，计算速度大幅提高，可以实现电力系统的在线静态安全
分析，解决了电力系统静态安全性分析中故障直流筛选耗时大的难题。