CN108804386A - 一种电力系统负荷裕度的并行化计算方法 - Google Patents

Abstract

本发明公开了一种电力系统负荷裕度的并行化计算方法，包括：采用基于负荷参数的二阶导数进行初值预估，得到临界点初值x⁰、零特征值对应的右特征向量初值g⁰、以及负荷裕度初值λ⁰；以上述参数为基础构建直接法非线性方程组对应的修正方程组，通过降阶变换将修正方程组拆分为四组同系数低维矩阵的线性方程组；基于CPU‑GPU混合架构，将Jacobi+ILU两阶段预处理和BICGSTAB迭代求解法相结合并行求解降阶之后的线性方程组，进而实现电力系统负荷裕度的快速求解。本发明保证了所求负荷裕度的正确性，计算速度较快，且能计算出一些连续潮流无法收敛系统的负荷裕度。

Description

一种电力系统负荷裕度的并行化计算方法

技术领域

本发明涉及电力系统领域，尤其涉及一种电力系统负荷裕度的并行化计算方法。

背景技术

负荷裕度作为度量电力系统电压稳定水平的性能指标，反应了系统承受负荷及故障扰动时，维持电压稳定的能力。随着电网规模不断扩大、可再生能源大规模并网、需求侧响应广泛应用，影响负荷裕度计算精度和效率的不确定因素越来越多。如何在全国联网和可再生能源大规模并网背景下，实现电力系统负荷裕度的快速、准确计算，对电力系统静态稳定在线评估具有重要实际意义。

当前，计算电力系统负荷裕度主要有两种方法：连续潮流法^[1](continuationpower flow，CPF)和直接法^[2]。CPF从系统当前运行点沿着系统负荷增长的方向，采用预测-校正方法追踪系统的静态电压稳定临界点，进而计算出系统的负荷裕度。该方法在计算过程中需不断的变换连续性参数，以处理雅可比矩阵奇异的难题，且在求解大规模系统时，计算量较大，占用内存较多，因而速度较慢；此外步长选择也是CPF的一个难点，较小的步长易提高CPF的计算耗时，而较大的步长易降低CPF的计算精度。直接法根据电压稳定临界点处潮流雅可比矩阵奇异且零特征值对应的特征向量不为0这一特点，构造一组表征电压稳定临界点性质的非线性方程组，然后采用牛顿法求解该方程组，进而计算出系统的负荷裕度。直接法原理简单，且在求取负荷裕度的同时，可根据雅可比矩阵零特征值的特征向量，有效甄别出系统的电压稳定关键节点。但直接法需求解的线性方程是原潮流方程维数的两倍，因而计算量较大。

近年来随着并行计算硬件技术的发展，图形处理器^[3](graphics processingunit,GPU)以其优异的运算性能为数值计算技术领域带来了一种新的思路。GPU作为中央处理器(central processing unit，CPU)良好的加速器，使得CPU-GPU混合架构具有超强的计算能力与并行性，已在电力系统多个领域得到了广泛应用，但尚未在负荷裕度计算中有所尝试。

由直接法计算负荷裕度的基本原理可知：采用直接法计算系统负荷裕度最为耗时的部分是求解修正方程组，若提高修正方程组的求解效率，必将大幅提升直接法的计算效率。而目前常用于求解大规模稀疏线性方程组的方法主要分为直接求解法^[4]和迭代求解法^[5]。直接求解法每求解出一个未知变量后，即认为该变量的求解已结束，无任何迭代修正过程。以基于LU分解^[6]的直接求解法为例，LU分解将一个矩阵的求解转化为对两个三角矩阵的求解，当系统规模增大时数据存储空间会以矩阵维数的平方倍增长，即算法的可伸缩性较弱，因而该方法主要应用于求解小规模稀疏线性方程组。此外，LU分解过程是逐行逐列进行的，待求变量的求解与已求变量密切相关，数据相关性强，不易实现并行计算。而迭代求解法的每次迭代都为整个未知向量产生一组解，并通过对这组解的不断修正使之逐步逼近真实解，该方法主要基于矩阵和向量的运算，数据相关性弱，适合并行计算，算法的可伸缩性也更强，所以在求解大规模稀疏线性方程组时具有明显的优势。使用迭代法求解线性方程时，可根据系数矩阵对称、正定等性质添加合适的预处理器，以改善迭代求解法收敛速度。

如何快速、准确地计算电力系统的负荷裕度，实现静态电压稳定性的快速评估仍是当前研究重点。

发明内容

本发明提供了一种电力系统负荷裕度的并行化计算方法。本发明保证了所求负荷裕度的正确性，对于大规模的系统，本发明不仅计算速度较快，还能计算出一些连续潮流无法收敛系统的负荷裕度；与其他采用单一常见预处理器的算法或者是仅基于CPU架构的算法相比，该发明同样体现出了明显的优势，能够实现电力系统负荷裕度的快速求解。详见下文描述：

一种电力系统负荷裕度的并行化计算方法，所述方法包括以下步骤：

采用基于负荷参数的二阶导数进行初值预估，得到临界点初值x⁰、零特征值对应的右特征向量初值g⁰、以及负荷裕度初值λ⁰；

以上述参数为基础构建直接法非线性方程组对应的修正方程组，通过降阶变换将修正方程组拆分为四组同系数低维矩阵的线性方程组；

基于CPU-GPU混合架构，将Jacobi+ILU两阶段预处理和BICGSTAB迭代求解法相结合并行求解降阶之后的线性方程组，进而实现电力系统负荷裕度的快速求解。

所述基于CPU-GPU混合架构具体为：

通过CPU形成预处理器矩阵L、U、J；

在第一次牛顿循环中，将系数矩阵W，线性方程等号右侧的向量b，预处理器矩阵L、U、J及其它相关参数传输至GPU中；

并行加速由GPU计算部分完成；

采用BICGSTAB迭代求解法求解α、β、γ和μ，CPU利用GPU所求结果α、β、γ和μ求得修正量Δx、Δg和Δλ，然后修正x、g和λ，并不断更新系数矩阵W，如此反复迭代直至满足收敛条件||f(x,λ)||<ε和||f_x·g||<ε，得到精确负荷裕度的过程，其实质为直接法中采用牛顿法求解负荷裕度时的循环。

所述将Jacobi+ILU两阶段预处理和BICGSTAB迭代求解法相结合并行求解降阶之后的线性方程组具体为：

1)确定初始值允许误差ε_b，计算残差初值并令循环次数i＝1，

ρ₀＝ξ＝ω₀＝1，v₀＝p₀＝0；

2)计算参数若ρ_i＝0，流程终止，输出失败信息，否则继续执行步骤3)；

3)令ψ＝(ρ_i/ρ_i-1)(ξ/ω_i-1)，p_i＝r_i-1+ψ(p_i-1－ω_i-1v_i-1)；

4)通过y＝L^-1p_i求得y，计算v_i＝JWU^-1y，

5)如果||h||≤ε_b，则退出循环，若不满足，继续执行步骤6)；

6)s_b＝r_i-1－ξv_i，令z＝L^-1s_b，t＝JWU^-1z，ω_i＝(t，s_b)/(t，t)，

7)若满足精度0.01，进入步骤8)；否则令r_i＝s_b－ω_it，i＝i+1，返回重新执行步骤3)；

8)解得后，求得α，同理依次求得解β,γ和μ。

优选地，

JWU^-1L^-1α_b＝Jb

其中，α_b为真实解，可根据α_b求取待求向量α；W为系数矩阵，α为待求向量。

本发明提供的技术方案的有益效果是：

1、本发明可实现采用直接法求解负荷裕度的并行化，提高电力系统负荷裕度的计算效率

2、与CPF计算误差相比，本发明可正确计算出系统的负荷裕度，具有极高的计算精度；

3、对系数矩阵进行Jacobi+ILU两阶段预处理，可有效降低BICGSTAB迭代求解法的累计迭代次数，改善基于GPU-CPU混合框架的电力系统负荷裕度并行求解效率；

4、相对于CPF算法，本文发明在大幅度提升电力系统负荷裕度计算效率同时，具有良好的算法稳定性，可计算出一些CPF无法收敛的测试系统的负荷裕度；

5、通过对比使用CPU架构和使用CPU-GPU混合架构所用时间可明显看出，并行计算的优势及其广泛的应用前景。

附图说明

图1是一种电力系统负荷裕度的并行化计算方法的流程图；

图2是本发明外循环流程图

图3为算例计算用时比较柱形图；

图4为另一算例计算用时比较柱形图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面对本发明实施方式作进一步地详细描述。

为了快速准确地计算出电力系统负荷裕度，本发明实施例提出一种电力系统负荷裕度的并行化计算方法。

首先，本发明实施例使用基于负荷参数的二阶导数进行初值预估，得到临界点初值x⁰、零特征值对应的右特征向量初值g⁰、负荷裕度初值λ⁰；然后构建非线性方程组对应的牛顿法修正方程组，并通过降阶变换将该方程组拆解为四组同系数低维矩阵的线性方程组；最后基于CPU-GPU混合架构，采用BICGSTAB迭代求解法协同两阶段预处理并行求解四组线性方程组，以得到精确的负荷裕度。

实施例1

本发明实施例提供了一种电力系统负荷裕度的并行化计算方法，该方法包括以下步骤：

101：采用基于负荷参数的二阶导数进行初值预估，得到临界点初值x⁰、零特征值对应的右特征向量初值g⁰、以及负荷裕度初值λ⁰；

102：以上述参数为基础构建直接法非线性方程组对应的修正方程组，通过降阶变换将修正方程组拆分为四组同系数低维矩阵的线性方程组；

103：基于CPU-GPU混合架构，将Jacobi+ILU两阶段预处理和BICGSTAB迭代求解法相结合并行求解降阶之后的线性方程组，进而实现电力系统负荷裕度的快速求解。

综上所述，本发明实施例通过上述步骤101-步骤103，基于CPU-GPU混合架构，将Jacobi+ILU两阶段预处理和BICGSTAB迭代求解法相结合并行求解降阶之后的线性方程组，提高电力系统负荷裕度的计算效率。

实施例2

本发明实施例的总体框架如图2所示，可分为CPU计算部分和GPU计算部分，并行加速由GPU计算部分完成。

其中，GPU计算部分主要包括：

采用BICGSTAB迭代求解法求解降阶后四组线性方程组，解得α、β、γ和μ的过程，本发明实施例将该过程称为内循环；

其中，CPU计算部分主要包括：

利用GPU所求结果α、β、γ和μ求得修正量Δx、Δg和Δλ，根据修正量修正参数x、g和λ，并不断更新系数矩阵W，如此反复迭代直至同时满足收敛条件||f(x,λ)||<ε和||f_x·g||<ε，得到精确负荷裕度的过程。

本发明实施例将CPU计算部分与GPU计算部分统称为外循环，其实质为直接法中求解负荷裕度时所用的牛顿法循环，因此下文称外循环为牛顿循环。

下面结合具体的计算公式、附图对实施例1中的方案进行进一步地介绍，详见下文描述：

本发明实施例以下均以线性方程组式(1)为例，介绍如何选择和采用预处理器处理系数矩阵。

Wα＝b (1)

式中，W为系数矩阵，α为待求向量，b为等号右侧的已知向量。

雅可比预处理器^[10](Jacobi preconditioner)是一种比较常用的预处理器，其通过提取式(2)所示的系数矩阵对角元素的倒数形成预处理矩阵：

式中，J为雅可比预处理器；W_(n+1)(n+1)为系数矩阵W第(n+1)行第(n+1)列的元素。

雅可比预处理后系数矩阵W_m如式(3)所示：

W_m＝JW (3)

式中，J为变换形式后的雅可比预处理器。为简化描述，本发明实施例仍然称此预处理器为雅可比预处理器。

式(1)中变量α可通过式(4)求解：

JWα＝Jb (4)

为简单描述，本发明实施例称雅可比预处理器为Jacobi预处理器。

对经过雅可比预处理的系数矩阵W_m，若将矩阵W_m分解为W_m＝LU+R，则称该分解为系数矩阵W_m的不完全LU分解预处理器(Incomplete LU decomposition preconditioner，ILU)，其中R为误差矩阵。本发明实施例选取ILU(0)预处理器^[11]对W_m进行预处理。为简化描述，下文将ILU(0)预处理器称为ILU预处理器。

系数矩阵W经雅可比预处理器和ILU预处理器预处理后，最终形如式(5)所示：

W_p＝JWU^-1L^-1 (5)

式中，W_p为采用两种预处理器后的矩阵形式，L、U矩阵为ILU预处理器产生的预处理矩阵。

由于基于GPU的BICGSTAB迭代求解法只支持采用一个预处理矩阵，因此本发明实施例得到预处理器J、L、U后，需将式(1)进一步变换为：

即将求解α转换为求解

本发明实施例，选取L为预处理矩阵，并通过求解式(7)得到具体过程如下：

JWU^-1L^-1α_b＝Jb (7)

最后，利用式(9)中与α的关系求得式(1)的解：

为了简化描述，本发明实施例将上述首先采用Jacobi预处理器，然后采用不完全LU分解预处理器处理系数矩阵的过程称为Jacobi+ILU两阶段预处理过程，下文统一用Jacobi+ILU两阶段预处理进行描述。

201：基于负荷参数的二阶导数进行初值预估，得到临界点初值x⁰、零特征值对应的右特征向量初值g⁰、以及负荷裕度初值λ⁰；

其中，该步骤201包括：

1)获取基础数据，包括：电力系统拓扑结构、支路参数和状态量初值预估的起点；

其中，上述获取基础数据的步骤为本领域技术人员所公知，本发明实施例对此不做赘述。

2)在电力系统电压稳定临界点处，电力系统潮流雅可比矩阵奇异，且该处零特征值对应的特征向量不为0，根据上述特点可得到一组表征电压稳定临界点性质的非线性方程组：

f(x,λ)＝0 (10)

f_x(x,λ)·g＝0 (11)

g_p-1＝0 (12)

式中，式(10)是含参数的潮流方程，f：Rⁿ×R，R为实数集，x∈Rⁿ代表各节点电压的幅值与相角；n＝npv+2×npq，npv、npq分别为电力系统中PV、PQ节点的个数；λ∈R表示电力系统的负荷裕度；式(11)是潮流雅可比矩阵的奇异方程，g∈Rⁿ为零特征值对应的右特征向量；式(12)是规范化方程，g_p表示g的第p个分量，选取初值时p为g⁰中绝对值最大分量对应的序号，在之后的过程中p为g中绝对值最大分量对应的序号。

3)选取弱节点的电压幅值变量x_p为连续性参数t，将潮流方程对连续性参数进行一阶、二阶求导，可得：

式中，f_x为潮流方程对x的一阶偏导数；f_xx为潮流方程对x的二阶偏导数；f_λ为潮流方程对λ的一阶偏导数；和分别为x对连续性参数的一阶导数和二阶导数；和分别为λ对连续性参数的一阶导数和二阶导数。

考虑到t≡x_p，故与式(13)、(14)联立，求得和

式中，和分别是x_p对自身的一阶、二阶导数。

4)采用二阶泰勒展开式近似表示各状态量初值，得到负荷裕度初值λ⁰、临界点初值x⁰和右特征向量的初值g⁰：

式中，和分别为负荷裕度、临界点状态量和零特征值对应的右特征向量初值预估的初始点，Δx_p为参数。

202：构建非线性方程组式(10)～(12)对应的牛顿法修正方程组，令牛顿循环次数k＝1；

其中，该步骤202包括：

1)假设Δx⁰、Δg⁰、Δλ⁰为各初值修正量，使其满足如下形式：

2)将式(16)分别按照泰勒级数展开，忽略包含Δx⁰、Δg⁰和Δλ⁰的二次项及以上高次项，则有：

式中，Δx、Δg和Δλ分别为x、g和λ在牛顿迭代循环过程中的修正量；表示除p个元素为1外其他均为0的单位行向量。

203：通过降阶变换将修正方程组拆分为四组同系数低维矩阵的线性方程组，减小计算量与复杂程度；

其中，该步骤203包括：

1)将式(17)展开：

f_x·Δx+f_λ·Δλ＝-f(x,λ) (18)

f_xx·g·Δx+f_x·Δg＝-f_x·g (19)

2)令

Δs＝Δx-Δx_p·g (21)

对式(18)～(20)作换元变换，则：

f_x·Δs+f_λ·Δλ＝-f(x,λ)+Δx_p·(-f_x·g) (22)

3)联立式(12)与式(21)，可知Δs_p＝0，然后结合式(22)：

由上式可知：

式中，α,β∈Rⁿ⁺¹可通过求解如下两个方程获得：

4)联立式(19)与式(21)，可形成如下方程组：

f_x·Δg+f_λ·Δλ＝-f_x·g-(f_xx·g·Δs-f_λ·Δλ)-f_xx·g·Δx_p·g (27)

根据分块矩阵乘法，有：

式中，A＝[f_xx·g(f_λ)]∈R^n×(n+1)。

5)将式(24)代入式(27)，整理得：

f_x·Δg+f_λ·Δλ＝-f_x·g-A·α-Δx_p·A·(β+g₁) (29)

式中g₁＝[g；0]∈Rⁿ⁺¹。

6)联立式(20)和式(29)，可形成如下方程组：

由式(30)可得：

式中，γ、μ∈Rⁿ⁺¹通过求解如下两个方程获得：

204：形成预处理器，将数据从CPU传输到GPU中；

其中，该步骤204包括：

1)形成预处理器矩阵L、U、J。

2)在第一次牛顿循环中，将系数矩阵W，线性方程等号右侧的向量，预处理器矩阵L、U、J及其它相关参数传输至GPU中。

由于本发明实施例始终采用第一次牛顿循环中系数矩阵W生成的雅可比预处理器和ILU预处理器，故在之后计算过程中，仅需更新系数矩阵W和等号右侧的向量。

205：采用基于GPU的Jacobi+ILU两阶段预处理和BICGSTAB迭代求解法相结合并行求解式(25)、(26)、(32)和(33)四组线性方程组；

本发明实施例以求解式(25)为例，详细介绍BICGSTAB的求解过程：

其中，该步骤205包括：

1)将式(25)变换为式(6)；

2)确定初始值允许误差ε_b，计算残差初值并令循环次数i＝1，

ρ₀＝ξ＝ω₀＝1，v₀＝p₀＝0；

3)计算参数若ρ_i＝0，算法终止，输出失败信息，否则继续执行步骤4)；

4)令ψ＝(ρ_i/ρ_i-1)(ξ/ω_i-1)，p_i＝r_i-1+ψ(p_i-1－ω_i-1v_i-1)；

5)通过y＝L^-1p_i求得y，计算v_i＝JWU^-1y，

6)如果||h||≤ε_b，则退出循环，若不满足，继续执行步骤7)；

7)s_b＝r_i-1－ξv_i，令z＝L^-1s_b，t＝JWU^-1z，ω_i＝(t，s_b)/(t，t)，

8)若满足精度0.01，进入步骤9)；否则令r_i＝s_b－ω_it，i＝i+1，返回重新执行步骤3)；

9)解得后，利用式(9)求得α。同理依次求得式(26)、(32)和(33)的解β,γ和μ。

以上参数均为BICGSTAB求解式(25)、(26)、(32)和(33)过程中的参数。

具体实现时，本发明实施例将BICGSTAB迭代求解法和非精确牛顿法⁰相结合求解式(25)，以改善系统负荷裕度的求解效率。BICGSTAB迭代求解法通过比较当前相对残差和预设阈值判断何时停止迭代，故可通过修改阈值控制其计算精度。而牛顿循环中对修正方程组式(17)求解是为了逐渐逼近非线性方程组式(10)～(12)的解，故每次BICGSTAB循环的求解精度无须过高，能实现收敛即可。因此，本发明实施例采用非精确牛顿法的思想松弛内循环求解精度，降低BICGSTAB迭代求解法的迭代次数，进一步提高算法计算效率，具体精度设定为0.01，即最终所得误差小于0.01。

206：GPU部分运算完成后，将式(25)、(26)、(32)和(33)的解α、β、γ和μ与右特征向量g从GPU传回CPU中。

207：在CPU中，利用α、β、γ、μ与修正变量的关系，通过式(34)求得Δx^k、Δλ^k和Δg^k：

式中，k＝0,1,2…,i＝1,2…,n。

然后将修正变量代入式(35)中修正x、g和λ：

208：将修正后的x^k、g^k和λ^k代入式(10)和(11)，若满足精度0.01输出负荷裕度λ^k，程序结束；否则，令k＝k+1，重新执行步骤203，若已达到预设的最大迭代次数100(该取值根据实际应用中的需要进行设定，本发明实施例对此不做限制)，输出“结果不收敛”，程序结束。

综上所述，本发明实施例通过上述步骤201-步骤208以直接法为基础，提出一种电力系统负荷裕度的并行化计算方法，以提升直接法求解系统负荷裕度的计算速率。该方法首先基于负荷参数的二阶导数进行初值预估，得到临界点初值x⁰、零特征值对应的右特征向量初值g⁰、负荷裕度初值λ⁰；然后建构直接法非线性方程组对应的修正方程组，通过降阶变换将该方程组拆分为四组同系数低维矩阵的线性方程组；最后基于CPU-GPU混合架构，将Jacobi+ILU两阶段预处理和BICGSTAB迭代求解法相结合并行求解式(25)、(26)、(32)和(33)，进而实现电力系统负荷裕度的快速求解。

实施例3

下面结合具体的实例、图2、图3、图4以及表1、表2对实施例1和2中的方案进行可行性验证，详见下文描述：

本实例是以文献[13]和[14]中的case1354pegase、case2383wp、case2746wop、case2869pegase、case3012wp、case5738、case7092、case9241pegase、case11624、case13173和case13802测试系统为例，验证本方法的可行性和有效性。所有算例功率增长方向均采用所有发电机和负荷同比例增长，算法中牛顿循环的收敛阈值ε为0.01，最大外循环迭代次数N为100；BICGSTAB迭代求解法的精度ε_b为0.01。

为验证所提算法的正确性，本方法将计算结果与CPF的计算结果进行对比，结果如表1所示。由表1可知：除某些测试系统因CPF无法收敛外，其余场景下本方法计算出的系统负荷裕度与CPF计算结果进行对比，相对误差均很小，验证了本方法具有较高的计算精度。同时，需指出的是：本方法不仅具有较高的计算精度，对于一些CPF无法求解的测试系统的负荷裕度，如case3012wp、case13173与case13802测试系统，仍可计算其对应的负荷裕度，表明本方法具有较高的稳定性。

表1本方法与CPF的λ结果对比

为进一步验证本方法中Jacobi+ILU两阶段预处理的有效性，本节分别对比了采用Jacobi预处理器、ILU预处理器和Jacobi+ILU两阶段预处理对BICGSTAB迭代求解法累计内部迭代次数的影响，结果如表2所示。

表2中结果为求解α、β、γ和μ时采用不同预处理方法的累计内部迭代次数，预处理器均由第一次外循环的系数矩阵产生。表中J+ILU为采用两阶段预处理的场景，表中算法均采用第一次牛顿循环过程中系数矩阵产生的预处理器，该方法可在满足精度要求的基础上，有效减少CPU-GPU之间的数据传输量，表2中‘-’表示该算法不收敛。

表2不同预处理器对累计内部迭代次数的影响

由表3可知：当系数矩阵维度较大时，由于经Jacobi预处理后的系数矩阵特征值分布仍较为分散，使某些单独使用该预处理器的测试系统仍难以收敛，因此，在后述算例的验证分析中将不单独采用Jacobi预处理器。

进一步，对于相同BICGSTAB迭代求解法，分别采用ILU预处理器和Jacobi+ILU两阶段预处理对系数矩阵进行处理，对比结果发现：由于经Jacobi+ILU两阶段预处理过程后的系数矩阵特征值分布更为集中，采用Jacobi+ILU两阶段预处理的BICGSTAB迭代求解法累计内部迭代次数均小于采用ILU预处理器的累计内部迭代次数。此外，采用Jacobi+ILU两阶段预处理还可计算出case13173和case13802的负荷裕度，表明两阶段预处理相对ILU预处理器具有较高的稳定性。

表3的对比结果表明：采用本方法中Jacobi+ILU两阶段预处理方法对系数矩阵进行处理可大幅度降低BICGSTAB迭代求解法的累计内部迭代次数，提高线性方程求解效率，在大规模电力系统的负荷裕度求解中具有显著优势。

本方法进一步对比了CPF、基于CPU架构的Jacobi+ILU两阶段预处理和BICGSTAB迭代求解法相结合的方法、基于CPU-GPU混合架构的ILU预处理和BICGSTAB迭代求解法相结合的方法以及本方法的计算耗时，结果如图3和4所示。图中若没有显示出的计算耗时，表明该测试系统的计算结果不收敛，J+ILU为Jacobi+ILU两阶段预处理。

进一步，为了更直观比较上述算法的计算耗时，根据图3和4中的结果，本方法以CPF的计算耗时为基准，即CPF收敛时加速比为1，将其分别除以其他3种算法的计算耗时定义为加速比，结果如表4所示，表中加速比越大表示该算法的计算耗时越少。若某些测试系统CPF算法无法收敛，而其它算法可计算出该测试系统的负荷裕度，即该测试系统计算加速比无基准值时，因此表4中只给出该算法的收敛情况。

表4四种算法加速比比较

由表4可知：①CPF计算速度和稳定性均不太理想，对于某些测试系统，本方法虽可计算其对应的负荷裕度，而CPF却无法收敛，如case3012wp、case13802系统；②某些场景下，CPF虽可收敛，但其计算耗时非常大，如在case7092测试系统中，本方法相对于CPF的加速比可达42.799；③对于同一测试系统时，本方法中Jacobi+ILU两阶段预处理加速比均大于采用ILU预处理的加速比，说明Jacobi+ILU两阶段预处理可使矩阵特征值分布更为集中，有利于加快线性方程组的求解速度。

此外需要指出的是，当测试系统节点数较小时，如case1354pegase系统，基于CPU-GPU混合架构的负荷裕度求解算法并不占优，主要原因是：修正方程组维数较小时，GPU并行求解带来的加速性不足以抵消CPU和GPU之间数据交互式的耗时，因而导致采用CPU-GPU混合架构的负荷裕度求解算法计算时间较长。而当系统规模大于case2383wp后，CPU-GPU混合架构的计算优势将逐步显现。

上述分析表明：本方法基于CPU-GPU混合架构，采用Jacobi+ILU两阶段预处理和BICGSTAB迭代求解法相结合的并行求解方法在计算较大规模测试系统的负荷裕度时，可有效提高负荷裕度的计算效率，且具有良好的算法稳定性，较传统负荷裕度计算方法的优势更为明显。

本发明实施例对各器件的型号除做特殊说明的以外，其他器件的型号不做限制，只要能完成上述功能的器件均可。

本领域技术人员可以理解附图只是一个优选实施例的示意图，上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

参考文献

[1]赵晋泉，钱天能.计及发电机励磁电流约束和电枢电流约束的连续潮流[J].中国电机工程学报，2012，32(22)：118-125.

[2]A.B.Neto，D.A.Alves.Singularities analysis of the jacobian matrixmodified in the continuation power flow:performance evaluation[J].IEEE LatinAmerica Transactions,2017,15(11):2137-2143.

[3]X.Li，F.Li.GPU-based power flow analysis with Chebyshevpreconditioner and conjugate gradient method[J].Electric Power SystemsResearch，2014，116(11)：87-93.

[4]张永杰，孙秦.大型稀疏线性方程组符号LU分解法[J].计算机工程与应用，2007，43(28)：29-30.

[5]Kulkarni A Y，Pai M A，Sauer P W.Iterative solver techniques in fastdynamic calculations of power systems[J].International Journal ofElectricalPower&Energy Systems，2001，23(3)：237-244.

[6]陈颖，林锦贤，吕暾.LU分解和Laplace算法在GPU上的实现[J].计算机应用，2011，31(3)：851-855.

[7]江伟，王成山，余贻鑫，等.直接计算静态电压稳定临界点的新方法[J].中国电机工程学报，2006，26(10)：1-6.

[8]X.Li，F.Li.GPU-based two-step preconditioning for conjugategradient method in power flow[C]//IEEE Power&Energy Society General Meeting，2015：1-5.

[9]张芡.大规模稀疏线性系统的稀疏近似逆预处理技术[D].清华大学,2013.

[10]Benzi M.Preconditioning Techniques for Large Linear Systems:ASurvey[J].Journal of Computational Physics,2002,182(2):418-477.

[11]M.Benzi.Preconditioning Techniques for Large Linear Systems：ASurvey[J].Journal of Computational Physics，2002，182(2)：418-477.

[12]Dembo R S，Eisenstat S C，Steihaug T.Inexact Newton methods[J].SiamJournal on Numerical Analysis，1982，19(2)：400-408.

[13]R.D.Zimmerman，C.E.Murillo-Sánchez，R.J.Thomas，"MATPOWER：Steady-State Operations，Planning and Analysis Tools for Power Systems Research andEducation，"Power Systems,IEEE Transactions on，2011，26(1)：12-19.

[14]X.Li，F.Li，H.Yuan，et al.GPU-based Fast Decoupled Power Flow withPreconditioned Iterative Solver and InexactNewtonMethod[J].IEEE Transactionson Power Systems，2016，PP(99)：1-1

Claims (4)

1.一种电力系统负荷裕度的并行化计算方法，其特征在于，所述方法包括以下步骤：

采用基于负荷参数的二阶导数进行初值预估，得到临界点初值x⁰、零特征值对应的右特征向量初值g⁰、以及负荷裕度初值λ⁰；

以上述参数为基础构建直接法非线性方程组对应的修正方程组，通过降阶变换将修正方程组拆分为四组同系数低维矩阵的线性方程组；

基于CPU-GPU混合架构，将Jacobi+ILU两阶段预处理和BICGSTAB迭代求解法相结合并行求解降阶之后的线性方程组，进而实现电力系统负荷裕度的快速求解。

2.根据权利要求1所述的一种电力系统负荷裕度的并行化计算方法，其特征在于，所述基于CPU-GPU混合架构具体为：

通过CPU形成预处理器矩阵L、U、J；

在第一次牛顿循环中，将系数矩阵W，线性方程等号右侧的向量b，预处理器矩阵L、U、J及其它相关参数传输至GPU中；

并行加速由GPU计算部分完成；

采用BICGSTAB迭代求解法求解α、β、γ和μ，CPU利用GPU所求结果α、β、γ和μ求得修正量Δx、Δg和Δλ，然后修正x、g和λ，并不断更新系数矩阵W，如此反复迭代直至满足收敛条件||f(x,λ)||<ε和||f_x·g||<ε，得到精确负荷裕度的过程，其实质为直接法中采用牛顿法求解负荷裕度时的循环。

3.根据权利要求2所述的一种电力系统负荷裕度的并行化计算方法，其特征在于，所述将Jacobi+ILU两阶段预处理和BICGSTAB迭代求解法相结合并行求解降阶之后的线性方程组具体为：

1)确定初始值允许误差ε_b，计算残差初值并令循环次数i＝1，ρ₀＝ξ＝ω₀＝1，v₀＝p₀＝0；

2)计算参数若ρ_i＝0，流程终止，输出失败信息，否则继续执行步骤3)；

3)令ψ＝(ρ_i/ρ_i-1)(ξ/ω_i-1)，p_i＝r_i-1+ψ(p_i-1－ω_i-1v_i-1)；

4)通过y＝L^-1p_i求得y，计算v_i＝JWU^-1y，

5)如果||h||≤ε_b，则退出循环，若不满足，继续执行步骤6)；

6)s_b＝r_i-1－ξv_i，令z＝L^-1s_b，t＝JWU^-1z，ω_i＝(t，s_b)/(t，t)，

7)若满足精度0.01，进入步骤8)；否则令r_i＝s_b－ω_it，i＝i+1，返回重新执行步骤3)；

8)解得后，求得α，同理依次求得解β,γ和μ。

4.根据权利要求3所述的一种电力系统负荷裕度的并行化计算方法，其特征在于，

JWU^-1L^-1α_b＝Jb

其中，α_b为真实解，可根据α_b求取待求向量α；W为系数矩阵，α为待求向量。