CN106021943A - 一种结合gpu软硬件架构特点设计的直流故障筛选方法 - Google Patents

Abstract

本发明公开了一种结合GPU软硬件架构特点设计的直流故障筛选方法，包括：根据GPU的软硬件参数对GPU的内核函数的任务分配进行优化设计：使用一个线程块计算一个开断；CPU读取电网数据，将所需的数据传送给GPU；采用两个CUDA流，使GPU在执行内核函数的同时，实现支路开断故障筛选与发电机开断故障筛选的异步执行；CUDA流1：GPU中的内核函数1和内核函数2筛选支路开断故障集S1并回传CPU中；CUDA流2：GPU中的内核函数3和内核函数4筛选发电机开断故障集S2并回传CPU中。本发明使用一个线程块来计算一个开断，所用的线程总数量增加，单个线程的计算量减少，充分利用了GPU的硬件资源和计算能力。

Description

一种结合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)。

导纳矩阵：以系统元件的等值导纳为基础所建立的、描述电力网络各节点电压和注入电流之间关系的矩阵。

CUDA流：CUDA流表示一个GPU操作队列，并且该队列中的操作将以指定的顺序执行，每个流视为GPU的一个任务，而且这些任务之间可以并行执行。

本发明一种结合GPU软硬件架构特点设计的直流故障筛选方法，所述方法包括：

(1)根据GPU的软硬件参数对GPU的内核函数的任务分配进行优化设计：使用一个线程块计算一个开断；

(2)CPU读取电网数据，将所需的电网基础数据传送给GPU；

(3)采用两个CUDA流，使GPU在执行内核函数的同时，实现支路开断故障筛选与发电机开断故障筛选的异步执行。

(4)CUDA流1：GPU中的内核函数1计算支路开断时的节点电压相角θ₁；GPU中的内核函数2计算支路开断时各支路的有功功率并与有功功率限额比较，筛选出支路开断故障集S1并回传至CPU中；

(5)CUDA流2：GPU中的内核函数3计算发电机开断时的节点电压相角θ₂；GPU中的内核函数4计算发电机开断时各支路的有功功率并与有功功率限额比较，筛选出发电机开断故障集S2并回传至CPU中。

其中，所述步骤(1)中参与优化设定的GPU的软硬件参数包括：流处理器(SM)数量N_SM；每个SM可调度管理的最大线程数量S_thread，每个SM可调度管理的最大线程块数量S_block。

进一步，所述步骤(2)中所述电网基本数据包括电网节点数N、支路数N_branches、每条支路首尾所连节点编号df[N_branches]和dt[N_branches]、支路电抗矩阵X、支路有功功率限额P_b[N_branches]、发电机数量N_gen、发电机所连节点编号dg[N_gen]、发电机所发有功功率Pg[N_gen]、节点电压相角θ₀、基态节点注入有功功率P₀、节点电抗矩阵X₀。

优选的，所述步骤(4)中内核函数1为Kernel_1<N_blocks，N_threads>，Kernel_1的线程块大小N_threads＝S_thread/S_block，线程块数量N_blocks＝N_branches，总线程数量为：N_branches×N_threads；Kernel_1的第b号线程计算第b条支路断开时的电压相角向量θ₁＝θ₀+Δθ，设第b条支路的首尾节点编号分别为i、j，其中Δθ＝c(θ_0i-θ_0j)X_L，式中：X_L＝[X_1i-X_1j,…,X_ni-X_nj]^T，Δx_ij为断开支路的电抗值，X_ij为基态的电抗矩阵X₀的第(i,j)号元素，θ_0i为节点电压相角θ₀的第i号元素；

Kernel_1<N_blocks，N_threads>的具体计算方法为：

(1.1)CUDA自动为每个线程块和线程分配线程索引blockID，threadID分别作为线程块和线程的编号；

(1.2)分别将线程块和线程的编号赋值给变量b，t，通过变量b，t来指定b号线程块中的t号线程；

(1.3)b号线程块中的N_threads个线程计算第b条支路开断时N个节点电压的相角即：t号线程计算t+m×N_threads号节点电压相角，其中m为N_threads除以N的商；

(1.4)将第b条支路首尾节点编号分别赋值给变量i，j，即：i←df[b]，j←dt[b]；

(1.5)用公式c＝(X_ii+X_jj-2×X_ij-Δx_ij)计算c，如果c为0，置S1[b]为2，若不为0，c＝1/c；

(1.6)判断变量t的值是否小于电网节点数N，小于继续执行以下步骤，否则该线程结束运算；

(1.7)计算第b条支路开断时的各节点电压相角向量：

${\left({\mathrm{\&theta;}}_1^b\right)}_t={\left({\mathrm{\&theta;}}_0\right)}_t-c\&times;\&lsqb;{\left({\mathrm{\&theta;}}_0\right)}_i-{\left({\mathrm{\&theta;}}_0\right)}_j\&rsqb;\&times;\&lsqb;X_{ti}-X_{tj}\&rsqb;$

其中，表示第b条支路开断时第t节点的电压相角，(θ₀)_t表示基态时t节点的电压相角；

(1.8)t←t+N_threads，返回步骤1.6；

(1.9)N_blocks个线程块中的全部线程计算结束得到节点电压相角θ₁和故障集S1，供内核函数2使用。

进一步，所述步骤(4)中所述内核函数2为Kernel_2<N_blocks，N_threads>，Kernel_2的线程块大小N_threads＝S_thread/S_block，线程块数量N_blocks＝N_branches，总线程数量为：N_branches×N_threads；Kernel_2的第b号线程块计算第b条支路断开时各支路的有功功率并判断b条支路开断是否属于严重故障集；b号线程块中的计算公式为式中：P_ij ^b为第b条支路开断时首尾节点编号分别为i和j支路的有功功率，是内核函数1计算支路开断时的节点电压相角θ₁的第(i，b)和(j，b)元素；

Kernel_2<N_blocks，N_threads>使用Kernel_1的计算结果θ₁和S1，其具体计算方法为：

(2.1)CUDA自动为每个线程块和线程分配索引block ID，thread ID分别作为线程块和线程的编号；

(2.2)分别将线程块和线程的编号赋值给变量b，t，通过变量b，t来指定b号线程块中的t号线程；

(2.3)b号线程块中的N_threads个线程计算第b条支路开断时其余N_branches-1条支路的有功功率，即：t号线程计算第t+m×N_threads条支路的有功功率，其中m为N_threads除以N_branches-1的商；

(2.4)判断变量t的值是否小于支路数N_branches，并判断S1[b]不为2且t不等于b，否则该线程结束运算；

(2.5)第t条支路的首尾节点编号df[t]和dt[t]分别赋值给变量i，j，即：i←df[t]，j←dt[t]；

(2.6)计算第b条支路开断时的第t条支路的有功功率：

$P_{ij}^b=\frac{{\left({\mathrm{\&theta;}}_1^b\right)}_i-{\left({\mathrm{\&theta;}}_1^b\right)}_j}{x_{ij}}$

其中，P_ij ^b为第b条支路开断时第t条支路的有功功率，是内核函数1计算支路开断时的节点电压相角θ₁的第(i，b)和(j，b)元素，即第b条支路开断时i，j节点电压相角；并将结果与有功功率限额比较，大于限额，置S1[b]为1；

(2.7)t←t+N_threads，返回步骤2.4；

(2.8)N_blocks个线程块中的全部线程计算结束，得到S1并传回CPU处理。

优选的，所述步骤(5)中所述内核函数3为Kernel_3<N_blocks，N_threads>，

Kernel_3<N_blocks，N_threads>的线程块大小N_threads＝S_thread/S_block，线程块数量N_blocks＝N_gen，总线程数量为：N_gen×N_threads；Kernel_3的第b号线程块计算第b台发电机断开时的节点电压相角θ₂＝θ₀+Δθ，设第b台发电机所连接电网节点编号为j，其中式中，表示第b台发电机开断时t节点电压相角；

其中，Kernel_3<N_blocks，N_threads>的具体计算方法为：

(3.1)CUDA自动为每个线程块和线程分配线程索引blockID，threadID分别作为线程块和线程的编号；

(3.2)分别将线程块和线程的编号赋值给变量b，t，通过变量b，t来指定b号线程块中的t号线程；

(3.3)b号线程块中的N_threads个线程计算第b台发电机开断时N个节点电压的相角即：t号线程计算t+m×N_threads号节点电压相角，其中m为N_threads除以N的商；

(3.4)将第b台发电机所连接电网节点编号赋值j，即：j←dt[b]；

(3.5)判断变量t的值是否小于电网节点数N，小于继续执行以下步骤，否则该线程结束运算；

(3.6)计算第b台发电机开断时的各节点电压相角向量：

${\left({\mathrm{\&theta;}}_2^b\right)}_t={\left({\mathrm{\&theta;}}_0\right)}_t-{\mathrm{\&Delta;P}}_j{X}_{tj}$

其中，表示第b台发电机开断时第t节点的电压相角，(θ₀)_t表示基态时t节点的电压相角；

(3.7)t←t+N_threads，返回步骤3.6；

(3.8)N_blocks个线程块中的全部线程计算结束得到节点电压相角θ₂和故障集S2，供内核函数4使用。

再者，所述步骤(5)中所述内核函数4为Kernel_4<N_blocks，N_threads>，Kernel_4的线程块大小N_threads＝S_thread/S_block，线程块数量N_blocks＝N_gen，总线程数量为：N_gen×N_threads；Kernel_4的第b号线程块计算第b台发电机断开时各支路的有功功率并判断第b台发电机开断是否属于严重故障集。线程块中的计算公式为式中：P_ij ^b为第b台发电机开断时首尾节点编号分别为i和j支路的有功功率，是内核函数3计算发电机开断时的节点电压相角θ₂的第(i，b)和(j，b)元素，并将结果与限额比较，大于限额，置S2[b]为1；

其中，Kernel_4<N_blocks，N_threads>的具体计算方法为：

(4.1)CUDA自动为每个线程块和线程分配索引block ID，thread ID分别作为线程块和线程的编号；

(4.2)分别将线程块和线程的编号赋值给变量b，t，通过变量b，t来指定b号线程块中的t号线程；

(4.3)b号线程块中的N_threads个线程计算第b台发电机开断时各支路的有功功率，即：t号线程计算第t+m×N_threads条支路的有功功率，其中m为N_threads除以N_branches商；

(4.4)判断变量t的值是否小于支路数N_branches，否则该线程结束运算；

(4.5)第t条支路的首尾节点编号df[t]和dt[t]分别赋值给变量i，j，即：i←df[t]，j←dt[t]；

(4.6)计算第b台发电机开断时的第t条支路的有功功率：

${P_{ij}}^b=\frac{{\left({\mathrm{\&theta;}}_2^b\right)}_i-{\left({\mathrm{\&theta;}}_2^b\right)}_j}{x_{ij}}$

其中，P_ij ^b为第b台发电机开断时第t条支路的有功功率，是内核函数3计算发电机开断时的节点电压相角θ₂的第(i，b)和(j，b)元素，即第b条发电机开断时i，j节点电压相角；并将结果与有功功率限额比较，大于限额，置S2[b]为1；

(4.7)t←t+N_threads，返回步骤4.4；

(4.8)N_blocks个线程块中的全部线程计算结束，得到S2并传回CPU处理。

优选的，计算所述内核函数1的计算公式和内核函数3的计算公式时，先将电抗矩阵X转置形成矩阵X^T，然后采用公式完成计算。GPU在计算内核函数1和内核函数3中的计算公式时，严格同步执行的线程束中的32个线程访问电抗矩阵X的某一列中的32个元素，由于一列中元素的存储地址不连续，无法合并访问。而将电抗矩阵X转置形成矩阵X^T，再进行计算时，线程束中的32个线程将访问32个连续存储矩阵行元素，实现了高效的合并访问。

有益效果：与现有技术相比，本发明具有以下显著优点：首先本发明结合GPU软硬件架构特点设计对基于GPU加速的直流故障筛选进行任务分配优化，使用一个线程块来计算一个开断，使所用的线程总数量增加，单个线程的计算量减少，充分利用了GPU的硬件资源和计算能力；其次，根据GPU线程块中线程以线程束为单位执行及线程束访问内存的模式特点来优化合并访问设计，内存操作速度提高了接近16倍；再者，对逻辑上相互独立的支路开断和发电机开断采用CUDA流技术，使其异步执行，充分利用了GPU的硬件资源，使计算时间减少了百分之四十左右。该优化设计方法显著提高了直流故障筛选方法的计算效率，解决了电力系统静态安全性分析中故障直流筛选耗时较大的问题。

附图说明

图1为本发明结合GPU软硬件架构特点设计的直流故障筛选算法中电网算例基本数据；

图2为本发明结合GPU软硬件架构特点设计的直流故障筛选算法的内核函数1算法流程图；

图3(a)为本发明结合GPU软硬件架构特点设计的直流故障筛选算法的未合并访问时的设计图；

图3(b)为本发明结合GPU软硬件架构特点设计的直流故障筛选算法的合并访问优化设计图；

图4为本发明结合GPU软硬件架构特点设计的直流故障筛选算法的采用CUDA优化设计的算法流程图；

图5为本发明结合GPU软硬件架构特点设计的直流故障筛选算法的测试结果。

具体实施方式

本发明公开了一种结合GPU软硬件架构特点设计的直流故障筛选方法，所述方法包括：

(1)根据GPU的软硬件参数对GPU的内核函数的任务分配进行优化设计：使用一个线程块计算一个开断；

(2)CPU读取电网数据，将所需的电网基础数据传送给GPU；

(3)通过使用两个CUDA流，使GPU在执行内核函数的同时，能实现支路开断与发电机开断的同步计算。

(4)CUDA流1：GPU中的内核函数1计算支路开断时的节点电压相角θ₁；GPU中的内核函数2计算支路开断时各支路的有功功率并与有功功率限额比较，筛选出支路开断故障集S1并回传至CPU中；

(5)CUDA流2：GPU中的内核函数3计算发电机开断时的节点电压相角θ₂；GPU中的内核函数4计算发电机开断时各支路的有功功率并与有功功率限额比较，筛选出发电机开断故障集S2并回传至CPU中。

一、根据不同故障的特点对直流潮流算法进行优化；简化直流潮流计算公式，避免不必要的浮点计算量；

(1)直流潮流的计算公式为：

P＝B_n×nθ或者θ＝X_n×nP (1)

$P_{ij}=\frac{{\mathrm{\&theta;}}_i-{\mathrm{\&theta;}}_j}{x_{ij}}---2)$

其中，

n是电网节点个数；

B是节点电纳矩阵，n*n维矩阵；

P是节点有功注入功率，n维向量；

θ是节点电压角度，n维向量；

θ_i，θ_j是节点i和节点j电压相角；

X＝B^-1是节点电抗矩阵，n*n维矩阵；

x_ij表示节点i和节点j之间的支路电抗矩阵值；

在基态时，直流潮流计算公式可以表示为：θ₀＝X₀P₀。

当发生故障的时，公式(1)可以表示为：

θ₁＝X₁P₁＝(X₀+ΔX)(P₀+ΔP)＝θ₀+Δθ (3)

Δθ＝ΔXP₀+X₀ΔP+ΔXΔP (4)

其中，下标0和1分别表示故障前状态和故障后状态；Δ表示故障前到故障后的改变；

如果节点i和节点j之间增加一条支路，电抗矩阵的增量ΔX可以表示为：

$\mathrm{\&Delta;}X=\frac{X_0{\mathrm{MM}}^T{X}_0}{X_{ii}+X_{jj}-2X_{ij}+{\mathrm{\&Delta;x}}_{ij}}={\mathrm{cX}}_0{\mathrm{MM}}^T{X}_0---\left(5\right)$

其中，M＝e_i-e_j，Δx_ij为增加支路的电抗值，如果是开断支路，Δx_ij为负值，x_ij为基态的节点电抗矩阵的第(i,j)号元素，e_i，e_j为标准基。

(2)直流开断的算法流程如下：

1)在CPU上只计算一次电网基态下的节点电抗矩阵X₀。

2)用式5计算支路开断后的电抗矩阵增量ΔX。

3)用式3和4计算节点电压相角Δθ和θ₁，

4)用式2计算各支路有功功率，并检查支路是否过负荷。

当c的分母为0时，说明节点i和j之间支路断开会导致系统解列。

(3)根据支路开断故障和发电机开断故障的不同特点对直流潮流算法进行优化。当节点i和节点j之间的支路断开时，电网电抗矩阵发生变化，各节点注入功率增量为零，即ΔP＝0。式4可以简化为：

Δθ＝ΔXP₀＝cX₀MM^TX₀P₀＝cX_LM^Tθ₀＝c(θ_0i-θ_0j)X_L (6)

其中，X_L＝X₀M＝X₀(e_i-e_j)＝[X_1i-X_1j,…,X_ni-X_nj]^T。

当连接于节点j的发电机断开时，电网电抗矩阵不变，即ΔX＝0，节点注入功率增量ΔP是稀疏向量，只有ΔP_j非零。式4可以简化为：

Δθ＝X₀ΔP＝ΔP_jX_j (7)

这里，ΔP＝[0…ΔP_j…0]^T，X_j＝[X_1j…X_nj]^T，X_j为基态下的节点电抗矩阵X₀的第j列元素，X_ij为基态下的节点电抗矩阵X₀的第(i,j)号元素，ΔP_j为连接于节点j的发电机的有功功率。

二、算法设计：本算法将程序分为CPU和GPU两个部分。

CPU读取电网数据，为GPU准备基本数据，处理GPU的计算结果，并控制整个程序的流程。GPU主要处理密集的浮点运算，共采用四个内核函数来完成故障筛选功能。本发明程序中，支路开断和发电机开断在逻辑上相互独立，将内核函数1和内核函数2放入CUDA流1，将内核函数3和内核函数4放入CUDA流2中，两个CUDA流异步执行。具体步骤为：1)创建两个流，stream1和stream2；2)stream1中包括四步：CPU到GPU传输数据，执行内核函数1，执行内核函数2，GPU到CPU传输数据；3)stream1中包括：CPU到GPU传输数据，执行内核函数3，执行内核函数4，GPU到CPU传输数据。程序流程见图4所示。其中，内核函数Kernel_1计算支路开断时的节点电压相角θ₁；内核函数Kernel_2计算支路开断时各支路的有功功率并与有功功率限额比较，筛选出支路开断故障集S1；内核函数Kernel_3计算发电机开断时的节点电压相角θ₂；内核函数Kernel_4计算发电机开断时各支路的有功功率并与有功功率限额比较，筛选出发电机开断故障集S2，具体流程如图4所示。

三、CPU生成基态下的节点电抗矩阵X₀，传输基础数据

CPU读取电网数据，生成节点电纳矩阵B，使用KLU库函数计算其逆矩阵节点电抗矩阵X₀。利用公式：θ₀＝X₀P₀，计算基态电网各节点电压相角θ₀。将直流筛选所需的数据传送给GPU，具体包括：开断所需的基础数据：电网节点数N，支路数N_branches，每条支路首尾所连节点编号df[N_branches]和dt[N_branches]，支路电抗矩阵X，支路有功功率限额P_b[N_branches]，发电机数量N_gen，发电机所连节点编号dg[N_gen]，发电机所发有功功率Pg[N_gen]，节点电压相角θ₀、基态节点注入有功功率P₀、节点电抗矩阵X₀。

四、根据所使用GPU的软硬件参数对4个GPU内核函数的任务分配进行优化设计

所述GPU的软硬件参数包括：流处理器(SM)数量N_SM；每个SM可调度管理的最大线程数量S_thread，每个SM可调度管理的最大线程块数量S_block。本发明使用使用Tesla K20C GPU进行优化设计，主要软硬件参数包括：流处理器(SM)数量N_SM：13；每个SM可调度管理管理的最大线程数量S_thread：2048；每个SM可调度管理的最大线程块数量S_block：16。

五、计算支路开断时的节点电压相角θ₁

GPU可以利用众多线程之间的自动转换来掩藏内存操作的时延。采用一个线程执行一个开断时，以图1中的算例3为例，内核函数1使用了12027个线程，仅仅是K20C约一半的计算能力。因此，本发明采用一个线程块计算一个开断，内核函数1和2所使用的线程块数N_blocks设计成等于支路数N_branches；内核函数3和4所使用的线程块数N_blocks设计成等于发电机数N_gen；线程块大小N_threads＝2048/16＝128。

内核函数1为Kernel_1<N_blocks，N_threads>，Kernel_1的线程块大小N_threads＝S_thread/S_block，线程块数量N_blocks＝N_branches，总线程数量为：N_branches×N_threads；Kernel_1的第b号线程计算第b条支路断开时的电压相角向量θ₁＝θ₀+Δθ，设第b条支路的首尾节点编号分别为i、j，其中Δθ＝c(θ_0i-θ_0j)X_L，式中：X_L＝[X_1i-X_1j,…,X_ni-X_nj]^T，Δx_ij为断开支路的电抗值，X_ij为基态的电抗矩阵X₀的第(i,j)号元素，θ_0i为节点电压相角θ₀的第i号元素；

Kernel_1<N_blocks，N_threads>的具体计算方法为：

(1.1)CUDA自动为每个线程块和线程分配线程索引blockID，threadID分别作为线程块和线程的编号；

(1.2)分别将线程块和线程的编号赋值给变量b，t，通过变量b，t来指定b号线程块中的t号线程；

(1.3)b号线程块中的N_threads个线程计算第b条支路开断时N个节点电压的相角θ₁ ^b，即：t号线程计算t+m×N_threads号节点电压相角，其中m为N_threads除以N的商；

(1.4)将第b条支路首尾节点编号分别赋值给变量i，j，即：i←df[b]，j←dt[b]；

(1.5)用公式c＝(X_ii+X_jj-2×X_ij-Δx_ij)计算c，如果c为0，置S1[b]为2，若不为0，c＝1/c；

(1.6)判断变量t的值是否小于电网节点数N，小于继续执行以下步骤，否则该线程结束运算；

(1.7)计算第b条支路开断时的各节点电压相角向量：

${\left({\mathrm{\&theta;}}_1^b\right)}_t={\left({\mathrm{\&theta;}}_0\right)}_t-c\&times;\&lsqb;{\left({\mathrm{\&theta;}}_0\right)}_i-{\left({\mathrm{\&theta;}}_0\right)}_j\&rsqb;\&times;\&lsqb;X_{ti}-X_{tj}\&rsqb;$

其中，表示第b条支路开断时第t节点的电压相角，(θ₀)_t表示基态时t节点的电压相角；

计算所述内核函数1的计算公式时，先将电抗矩阵X转置形成矩阵X^T，然后采用公式完成计算。GPU在计算内核函数1中的计算公式时，严格同步执行的线程束中的32个线程访问电抗矩阵X的某一列中的32个元素，由于一列中元素的存储地址不连续，无法合并访问，如图3(a)所示。而将电抗矩阵X转置形成矩阵X^T，再进行计算时，线程束中的32个线程将访问32个连续存储矩阵行元素，实现了高效的合并访问，如图3(b)所示。

(1.8)t←t+N_threads，返回步骤1.6；

(1.9)N_blocks个线程块中的全部线程计算结束得到节点电压相角θ₁和故障集S1，供内核函数2使用，如图2所示。

六、计算支路开断时各支路的有功功率并筛选出支路开断故障集S1

内核函数2为Kernel_2<N_blocks，N_threads>，Kernel_2的线程块大小N_threads＝S_thread/S_block，线程块数量N_blocks＝N_branches，总线程数量为：N_branches×N_threads；Kernel_2的第b号线程块计算第b条支路断开时各支路的有功功率并判断b条支路开断是否属于严重故障集；b号线程块中的计算公式为式中：P_ij ^b为第b条支路开断时首尾节点编号分别为i和j支路的有功功率，是内核函数1计算支路开断时的节点电压相角θ₁的第(i，b)和(j，b)元素；

Kernel_2<N_blocks，N_threads>使用Kernel_1的计算结果θ₁和S1，其具体计算方法为：

(2.1)CUDA自动为每个线程块和线程分配索引block ID，thread ID分别作为线程块和线程的编号；

(2.2)分别将线程块和线程的编号赋值给变量b，t，通过变量b，t来指定b号线程块中的t号线程；

(2.3)b号线程块中的N_threads个线程计算第b条支路开断时其余N_branches-1条支路的有功功率，即：t号线程计算第t+m×N_threads条支路的有功功率，其中m为N_threads除以N_branches-1的商；

(2.4)判断变量t的值是否小于支路数N_branches，并判断S1[b]不为2且t不等于b，否则该线程结束运算；

(2.5)第t条支路的首尾节点编号df[t]和dt[t]分别赋值给变量i，j，即：i←df[t]，j←dt[t]；

(2.6)计算第b条支路开断时的第t条支路的有功功率：

${P_{ij}}^b=\frac{{\left({\mathrm{\&theta;}}_1^b\right)}_i-{\left({\mathrm{\&theta;}}_1^b\right)}_j}{x_{ij}}$

其中，P_ij ^b为第b条支路开断时第t条支路的有功功率，是内核函数1计算支路开断时的节点电压相角θ₁的第(i，b)和(j，b)元素，即第b条支路开断时i，j节点电压相角；并将结果与有功功率限额比较，大于限额，置S1[b]为1；

(2.7)t←t+N_threads，返回步骤2.4；

(2.8)N_blocks个线程块中的全部线程计算结束，得到S1并传回CPU处理。

七、计算发电机开断时的节点电压的相角θ₂

内核函数3为Kernel_3<N_blocks，N_threads>，

Kernel_3<N_blocks，N_threads>的线程块大小N_threads＝S_thread/S_block，线程块数量N_blocks＝N_gen，总线程数量为：N_gen×N_threads；Kernel_3的第b号线程块计算第b台发电机断开时的节点电压相角θ₂＝θ₀+Δθ，设第b台发电机所连接电网节点编号为j，其中式中，表示第b台发电机开断时t节点电压相角；

其中，Kernel_3<N_blocks，N_threads>的具体计算方法为：

(3.1)CUDA自动为每个线程块和线程分配线程索引blockID，threadID分别作为线程块和线程的编号；

(3.2)分别将线程块和线程的编号赋值给变量b，t，通过变量b，t来指定b号线程块中的t号线程；

(3.3)b号线程块中的N_threads个线程计算第b台发电机开断时N个节点电压的相角即：t号线程计算t+m×N_threads号节点电压相角，其中m为N_threads除以N的商；

(3.4)将第b台发电机所连接电网节点编号赋值j，即：j←dt[b]；

(3.5)判断变量t的值是否小于电网节点数N，小于继续执行以下步骤，否则该线程结束运算；

(3.6)计算第b台发电机开断时的各节点电压相角向量：

${\left({\mathrm{\&theta;}}_2^b\right)}_t={\left({\mathrm{\&theta;}}_0\right)}_t-{\mathrm{\&Delta;P}}_j{X}_{tj}$

其中，表示第b台发电机开断时第t节点的电压相角，(θ₀)_t表示基态时t节点的电压相角；

计算内核函数3的计算公式时，先将电抗矩阵X转置形成矩阵X^T，然后采用公式完成计算。GPU在计算内核函数3中的计算公式时，严格同步执行的线程束中的32个线程访问电抗矩阵X的某一列中的32个元素，由于一列中元素的存储地址不连续，无法合并访问，如图3(a)所示。而将电抗矩阵X转置形成矩阵X^T，再进行计算时，线程束中的32个线程将访问32个连续存储矩阵行元素，实现了高效的合并访问，如图3(b)所示。

(3.7)t←t+N_threads，返回步骤3.6；

(3.8)N_blocks个线程块中的全部线程计算结束得到节点电压相角θ₂和故障集S2，供内核函数4使用。

八、计算发电机开断时各支路的有功功率并筛选发电机开断故障集S2

内核函数4为Kernel_4<N_blocks，N_threads>，Kernel_4的线程块大小N_threads＝S_thread/S_block，线程块数量N_blocks＝N_gen，总线程数量为：N_gen×N_threads；Kernel_4的第b号线程块计算第b台发电机断开时各支路的有功功率并判断第b台发电机开断是否属于严重故障集。线程块中的计算公式为式中：为第b台发电机开断时首尾节点编号分别为i和j支路的有功功率，是内核函数3计算发电机开断时的节点电压相角θ₂的第(i，b)和(j，b)元素，并将结果与限额比较，大于限额，置S2[b]为1；

其中，Kernel_4<N_blocks，N_threads>的具体计算方法为：

(4.1)CUDA自动为每个线程块和线程分配索引block ID，thread ID分别作为线程块和线程的编号；

(4.2)分别将线程块和线程的编号赋值给变量b，t，通过变量b，t来指定b号线程块中的t号线程；

(4.3)b号线程块中的N_threads个线程计算第b台发电机开断时各支路的有功功率，即：t号线程计算第t+m×N_threads条支路的有功功率，其中m为N_threads除以N_branches商；

(4.4)判断变量t的值是否小于支路数N_branches，否则该线程结束运算；

(4.5)第t条支路的首尾节点编号df[t]和dt[t]分别赋值给变量i，j，即：i←df[t]，j←dt[t]；

(4.6)计算第b台发电机开断时的第t条支路的有功功率：

${P_{ij}}^b=\frac{{\left({\mathrm{\&theta;}}_2^b\right)}_i-{\left({\mathrm{\&theta;}}_2^b\right)}_j}{x_{ij}}$

其中，P_ij ^b为第b台发电机开断时第t条支路的有功功率，是内核函数3计算发电机开断时的节点电压相角θ₂的第(i，b)和(j，b)元素，即第b条发电机开断时i，j节点电压相角；并将结果与有功功率限额比较，大于限额，置S2[b]为1；

(4.7)t←t+N_threads，返回步骤4.4；

(4.8)N_blocks个线程块中的全部线程计算结束，得到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。在GPU计算平台上分别对3个电网实例进行了测试，具体实例计算结果，如图5所示。相对电网规模比较小的算例，优化所得到的加速比并不显著，主要是因为GPU的线程计算能力没有得到充分地利用。以算例3为例，没有经优化之前加速比是4.94，经过任务分配优化之后，取得的加速比是15.1，是优化之前的3.05倍；经过内存合并访问优化之后，取得的加速比是27.6，是优化之前的5.40倍；经过CUDA流优化之后，取得的加速比是47.6，是优化之前的9.64倍，也即是，经三步优化之后，算例3的计算时间减少了9倍之多。而且从3个算例取得的加速比可见，随着电网规模的扩大，本发明的优化设计所取得的效果更加显著，计算速度大幅提高，可以实现电力系统的在线静态安全分析，很好地解决了电力系统静态安全性分析中故障直流筛选耗时较大的问题。

Claims (8)

1.一种结合GPU软硬件架构特点设计的直流故障筛选方法，其特征在于：所述方法包括：

(1)根据GPU的软硬件参数对GPU的内核函数的任务分配进行优化设计：使用一个线程块计算一个开断；

(2)CPU读取电网数据，将所需的电网基础数据传送给GPU；

(3)采用两个CUDA流，使GPU在执行内核函数的同时，实现支路开断故障筛选与发电机开断故障筛选的异步执行；

(4)CUDA流1：GPU中的内核函数1计算支路开断时的节点电压相角θ₁；GPU中的内核函数2计算支路开断时各支路的有功功率并与有功功率限额比较，筛选出支路开断故障集S1并回传至CPU中；

(5)CUDA流2：GPU中的内核函数3计算发电机开断时的节点电压相角θ₂；GPU中的内核函数4计算发电机开断时各支路的有功功率并与有功功率限额比较，筛选出发电机开断故障集S2并回传至CPU中。

2.根据权利要求1所述的基于GPU加速的直流故障筛选方法，其特征在于：所述步骤(1)中所述GPU的软硬件参数包括：流处理器(SM)数量N_SM；每个SM可调度管理的最大线程数量S_thread，每个SM可调度管理的最大线程块数量S_block。

3.根据权利要求1所述的基于GPU加速的直流故障筛选方法，其特征在于：所述步骤(2)中所述电网基本数据包括电网节点数N、支路数N_branches、每条支路首尾所连节点编号df[N_branches]和dt[N_branches]、支路电抗矩阵X、支路有功功率限额P_b[N_branches]、发电机数量N_gen、发电机所连节点编号dg[N_gen]、发电机所发有功功率Pg[N_gen]、节点电压相角θ₀、基态节点注入有功功率P₀、节点电抗矩阵X₀。

4.根据权利要求1所述的基于GPU加速的直流故障筛选方法，其特征在于：所述步骤(4)中内核函数1为Kernel_1<N_blocks，N_threads>，Kernel_1的线程块大小N_threads＝S_thread/S_block，线程块数量N_blocks＝N_branches，总线程数量为：N_branches×N_threads；Kernel_1的第b号线程计算第b条支路断开时的电压相角向量θ₁＝θ₀+Δθ，设第b条支路的首尾节点编号分别为i、j，其中Δθ＝c(θ_0i-θ_0j)Χ_L，式中：X_L＝[X_1i-X_1j,…,X_ni-X_nj]^T，Δx_ij为断开支路的电抗值，X_ij为基态的电抗矩阵X₀的第(i,j)号元素，θ_0i为节点电压相角θ₀的第i号元素；

Kernel_1<N_blocks，N_threads>的具体计算方法为：

(1.1)CUDA自动为每个线程块和线程分配线程索引blockID，threadID分别作为线程块和线程的编号；

(1.2)分别将线程块和线程的编号赋值给变量b，t，通过变量b，t来指定b号线程块中的t号线程；

(1.3)b号线程块中的N_threads个线程计算第b条支路开断时N个节点电压的相角即：t号线程计算t+m×N_threads号节点电压相角，其中m为N_threads除以N的商；

(1.4)将第b条支路首尾节点编号分别赋值给变量i，j，即：i←df[b]，j←dt[b]；

(1.5)用公式c＝(X_ii+X_jj-2×X_ij-Δx_ij)计算c，如果c为0，置S1[b]为2，若不为0，c＝1/c；

(1.6)判断变量t的值是否小于电网节点数N，小于继续执行以下步骤，否则该线程结束运算；

(1.7)计算第b条支路开断时的各节点电压相角向量：

${\left({\mathrm{\&theta;}}_1^b\right)}_t={\left({\mathrm{\&theta;}}_0\right)}_t-c\&times;\&lsqb;{\left({\mathrm{\&theta;}}_0\right)}_i-{\left({\mathrm{\&theta;}}_0\right)}_j\&rsqb;\&times;\&lsqb;X_{ti}-X_{tj}\&rsqb;$

其中，表示第b条支路开断时第t节点的电压相角，(θ₀)_t表示基态时t节点的电压相角；

(1.8)t←t+N_threads，返回步骤1.6；

(1.9)N_blocks个线程块中的全部线程计算结束得到节点电压相角θ₁和故障集S1，供内核函数2使用。

5.根据权利要求1所述的基于GPU加速的直流故障筛选方法，其特征在于：所述步骤(4)中所述内核函数2为Kernel_2<N_blocks，N_threads>，Kernel_2的线程块大小N_threads＝S_thread/S_block，线程块数量N_blocks＝N_branches，总线程数量为：N_branches×N_threads；Kernel_2的第b号线程块计算第b条支路断开时各支路的有功功率并判断b条支路开断是否属于严重故障集；b号线程块中的计算公式为式中：P_ij ^b为第b条支路开断时首尾节点编号分别为i和j支路的有功功率，是内核函数1计算支路开断时的节点电压相角θ₁的第(i，b)和(j，b)元素；

Kernel_2<N_blocks，N_threads>使用Kernel_1的计算结果θ₁和S1，其具体计算方法为：

(2.1)CUDA自动为每个线程块和线程分配索引block ID，thread ID分别作为线程块和线程的编号；

(2.2)分别将线程块和线程的编号赋值给变量b，t，通过变量b，t来指定b号线程块中的t号线程；

(2.3)b号线程块中的N_threads个线程计算第b条支路开断时其余N_branches-1条支路的有功功率，即：t号线程计算第t+m×N_threads条支路的有功功率，其中m为N_threads除以N_branches-1的商；

(2.4)判断变量t的值是否小于支路数N_branches，并判断S1[b]不为2且t不等于b，否则该线程结束运算；

(2.5)第t条支路的首尾节点编号df[t]和dt[t]分别赋值给变量i，j，即：i←df[t]，j←dt[t]；

(2.6)计算第b条支路开断时的第t条支路的有功功率：

${P_{ij}}^b=\frac{{\left({\mathrm{\&theta;}}_1^b\right)}_i-{\left({\mathrm{\&theta;}}_1^b\right)}_j}{x_{ij}}$

其中，P_ij ^b为第b条支路开断时第t条支路的有功功率，是内核函数1计算支路开断时的节点电压相角θ₁的第(i，b)和(j，b)元素，即第b条支路开断时i，j节点电压相角；并将结果与有功功率限额比较，大于限额，置S1[b]为1；

(2.7)t←t+N_threads，返回步骤2.4；

(2.8)N_blocks个线程块中的全部线程计算结束，得到S1并传回CPU处理。

6.根据权利要求1所述的基于GPU加速的直流故障筛选方法，其特征在于：所述步骤(5)中所述内核函数3为Kernel_3<N_blocks，N_threads>，

Kernel_3<N_blocks，N_threads>的线程块大小N_threads＝S_thread/S_block，线程块数量N_blocks＝N_gen，总线程数量为：N_gen×N_threads；Kernel_3的第b号线程块计算第b台发电机断开时的节点电压相角θ₂＝θ₀+Δθ，设第b台发电机所连接电网节点编号为j，其中式中，表示第b台发电机开断时t节点电压相角；

其中，Kernel_3<N_blocks，N_threads>的具体计算方法为：

(3.1)CUDA自动为每个线程块和线程分配线程索引blockID，threadID分别作为线程块和线程的编号；

(3.2)分别将线程块和线程的编号赋值给变量b，t，通过变量b，t来指定b号线程块中的t号线程；

(3.3)b号线程块中的N_threads个线程计算第b台发电机开断时N个节点电压的相角即：t号线程计算t+m×N_threads号节点电压相角，其中m为N_threads除以N的商；

(3.4)将第b台发电机所连接电网节点编号赋值j，即：j←dt[b]；

(3.5)判断变量t的值是否小于电网节点数N，小于继续执行以下步骤，否则该线程结束运算；

(3.6)计算第b台发电机开断时的各节点电压相角向量：

${\left({\mathrm{\&theta;}}_2^b\right)}_t={\left({\mathrm{\&theta;}}_0\right)}_t-{\mathrm{\&Delta;P}}_j{X}_{tj}$

其中，表示第b台发电机开断时第t节点的电压相角，(θ₀)_t表示基态时t节点的电压相角；

(3.7)t←t+N_threads，返回步骤3.6；

(3.8)N_blocks个线程块中的全部线程计算结束得到节点电压相角θ₂和故障集S2，供内核函数4使用。

7.根据权利要求1所述的基于GPU加速的直流故障筛选方法，其特征在于：所述步骤(5)中所述内核函数4为Kernel_4<N_blocks，N_threads>，Kernel_4的线程块大小N_threads＝S_thread/S_block，线程块数量N_blocks＝N_gen，总线程数量为：N_gen×N_threads；Kernel_4的第b号线程块计算第b台发电机断开时各支路的有功功率并判断第b台发电机开断是否属于严重故障集。线程块中的计算公式为式中：P_ij ^b为第b台发电机开断时首尾节点编号分别为i和j支路的有功功率，是内核函数3计算发电机开断时的节点电压相角θ₂的第(i，b)和(j，b)元素，并将结果与限额比较，大于限额，置S2[b]为1；

其中，Kernel_4<N_blocks，N_threads>的具体计算方法为：

(4.1)CUDA自动为每个线程块和线程分配索引block ID，thread ID分别作为线程块和线程的编号；

(4.2)分别将线程块和线程的编号赋值给变量b，t，通过变量b，t来指定b号线程块中的t号线程；

(4.3)b号线程块中的N_threads个线程计算第b台发电机开断时各支路的有功功率，即：t号线程计算第t+m×N_threads条支路的有功功率，其中m为N_threads除以N_branches商；

(4.4)判断变量t的值是否小于支路数N_branches，否则该线程结束运算；

(4.5)第t条支路的首尾节点编号df[t]和dt[t]分别赋值给变量i，j，即：i←df[t]，j←dt[t]；

(4.6)计算第b台发电机开断时的第t条支路的有功功率：

${P_{ij}}^b=\frac{{\left({\mathrm{\&theta;}}_2^b\right)}_i-{\left({\mathrm{\&theta;}}_2^b\right)}_j}{x_{ij}}$

其中，P_ij ^b为第b台发电机开断时第t条支路的有功功率，是内核函数3计算发电机开断时的节点电压相角θ₂的第(i，b)和(j，b)元素，即第b条发电机开断时i，j节点电压相角；并将结果与有功功率限额比较，大于限额，置S2[b]为1；

(4.7)t←t+N_threads，返回步骤4.4；

(4.8)N_blocks个线程块中的全部线程计算结束，得到S2并传回CPU处理。

8.根据权利要求1所述的基于GPU加速的直流故障筛选方法，其特征在于：计算所述内核函数1的计算公式和内核函数3的计算公式时，先将电抗矩阵X转置形成矩阵X^T，然后采用公式完成计算。