CN101976835B - 一种大规模电力系统牛顿潮流的并行计算方法 - Google Patents

Abstract

一种大规模电力系统牛顿潮流的并行计算方法，属于电力系统安全分析计算技术领域。本发明利用计算机，根据并行集群，通过程序，结合雅克比矩阵更新的并行化技术，采用与准对角并行预条件子相结合的广义极小残量方法求解牛顿潮流的修正方程组，实现大规模电力系统潮流求解过程的并行化。本发明主要具有不需要任何形式的网络分割，计算简单，通用性好；算法内部不存在前推回代的依赖关系，并行效率高，计算速度快；能用常见的并行计算机、并行集群进行计算，便于推广应用等。本发明可广泛应用于电力系统并行潮流计算中，特别适用大规模复杂电力系统牛顿潮流并行计算。

Description

一种大规模电力系统牛顿潮流的并行计算方法

技术领域

本发明属于电力系统安全分析计算技术领域，具体涉及大规模电力系统在线安全分析的牛顿潮流并行计算方法。

背景技术

电力系统调度人员的主要任务是确保电力系统的安全、稳定、经济、优质运行。随着国民经济和电网互联的发展，对电力系统运行的可靠性和经济性提出了更高的要求，使调度人员面临巨大挑战：一方面，电网规模不断扩大，结构越来越复杂，运行方式灵活多变，多重故障、连锁故障造成的后果日益严重，在这种复杂的环境下，确保电力系统安全运行变得越来越困难；另一方面，随着电力市场化改革的推进，以及对节能降耗、环境保护的重视，作为指挥控制中心的电力系统调度部门在承担传统的调度任务以外，还增加了许多与电力交易、节能环保相关的工作，在运行中承担的角色和任务变得更加复杂和繁重。

《国家中长期科学和技术发展规划纲要》已经把“超大规模输配网和电网安全保障”列入了能源领域的优先主题，然而在当前电力系统网络结构和调度任务如此复杂的情况下，能量管理系统(EMS，Energy Management System)，已经不仅是辅助运行人员监视电力系统运行情况的监视系统，更是运行人员分析电力系统运行状态，做出各种决策的必备工具。而调度人员的电力系统潮流计算在能量管理系统中正处于核心地位。因此，要求电力系统潮流计算不仅是能量管理系统的必备功能，而且是一些其他功能的基础，如无功优化、静态安全分析需要反复调用潮流程序。从而，提高电力系统潮流计算的计算速度，降低潮流计算的误差，改善潮流计算的收敛性就能提高能量管理系统分析软件的计算速度、计算精度和适应性，给出正确决策，进而提高能量管理系统的效能。由于现代电力系统规模不断扩大、网络结构越来越复杂，任何模型化的计算分析问题都可能含有成千上万个方程。对于最基本的牛顿潮流计算，传统的串行计算方法已无法满足电网实时仿真和控制的需要，而并行计算为这类大规模运算提供了可行的途径。

现存电力系统牛顿潮流并行计算方法分为两类：一类是基于分块思想，即：按电力网络的地理位置或物理联接结构，将表征电力系统网络特性的系数矩阵A分解为通过协调部分互联的多个子矩阵，各子矩阵可独立求解，由此实现潮流的并行计算。如2006年2月第30卷第4期《电网技术》中的“基于支路切割方法的电力系统潮流并行协调算法”一文，公开的是一种通过切割电网支路实现网络分块的方法，该方法在计算前先从物理层面剖分电网，然后将切割支路的潮流以虚拟负荷的形式作为协调变量对子网络之间的耦合影响进行处理，再对各子网络单独进行潮流计算，最后合并。该方法的主要缺点是：由于电力网络结构千变万化，无法找到一种通用的网络分割方法来保证并行计算时的负载平衡，该类方法虽然对某些特定电网会取得比较好的并行效果，但无法推广应用到所有电力网络的潮流计算，算法的通用性较差。另一类是从线性方程组求解本身寻找算法的并行性。如2000年8月第33卷第4期《武汉水利电力大学学报》中的“求解大型稀疏线性方程组的一种并行算法及其在并行潮流计算中的应用”一文，公开的是一种适合于在消息传递并行计算机上求解大型电力网络方程组的并行方法，该方法采用稀疏矩阵因子表路径树表示LU分解法的分解过程、前推和回代过程，通过分割因子表路径树达到潮流并行计算的目的。该方法的主要缺点是：线性方程组并行求解过程中采用的均是传统的线性方程组求解算法，如LU分解法、Gauss迭代法等，由于这些算法内部存在前推回代依赖关系，无法实现粗粒度并行计算，导致该类牛顿潮流并行算法的并行效率较低。

发明内容

本发明的目的是针对现有电力系统牛顿潮流并行计算方法的不足，提供一种大规模电力系统牛顿潮流的并行计算方法，具有不需要任何形式的网络分割，通用性较好；算法内部不存在前推回代依赖关系，加之算法中的矩阵、向量运算本身存在潜在的并行性，能实现粗粒度并行计算，且并行计算的效率高等特点。

本发明的机理是：对任一n节点电力网络，设节点电压向量为

其中

牛顿法潮流计算中，以极坐标形式表示的修正方程组如下式所示：

$\left[\begin{array}{cc}\frac{\∂P}{\∂\mathrm{\θ}}& \frac{\∂P}{\∂U}U\\ \frac{\∂Q}{\∂\mathrm{\θ}}& \frac{\∂Q}{\∂U}U\end{array}\right]\left[\begin{array}{c}\mathrm{\Δ\θ}\\ \frac{\mathrm{\ΔU}}{U}\end{array}\right]=\left[\begin{array}{c}\mathrm{\ΔP}\\ \mathrm{\ΔQ}\end{array}\right]---\left(1\right)$

式中：

为雅克比矩阵，

为各未知节点电压相角和幅值的修正量，

为功率不平衡量。

潮流计算的目的就是计算式(1)中未知的PQ节点的电压U和相角θ、PV节点的无功Q和相角θ及平衡节点的有功P和无功Q。

设所有未知量组成的向量为x∈Rⁿ(即n维列向量)，潮流计算对应于非线性方程组的求解：

f(x)＝0                                    (2)

牛顿潮流计算是将式(2)对应的非线性潮流方程组的求解转化为式(1)对应的修正方程组的迭代求解。本发明在牛顿法潮流计算中，从修正方程组系数矩阵J(即雅克比矩阵)的形成过程和修正方程组的求解过程入手，分别给出一种基于预条件广义极小残量(GMRES，Generalized Minimal Residual)方法的并行线性方程组求解方法和一种雅克比矩阵更新的并行化方法，通过这两种方法来实现牛顿潮流的粗粒度并行求解。

实现本发明目的技术方案是：一种大规模电力系统牛顿潮流的并行计算方法，利用计算机，根据并行集群，通过程序，结合雅克比矩阵更新的并行化技术，采用与准对角并行预条件子相结合广义极小残量(GMRES，Generalized Minimal Residual)方法求解牛顿潮流的修正方程组，实现大规模电力系统潮流求解过程的并行化。其具体方法步骤如下：

(1)初始化

输入电力系统网络结构参数和电气参数，如节点导纳矩阵；PQ节点的注入有功功率、无功功率；PV节点的注入有功功率、负荷无功功率、节点电压等，并对未知变量赋初值：x₀＝x_old。

(2)并行雅克比矩阵更新

第(1)步完成后，首先矢量化雅克比矩阵的更新过程，然后并行求解。具体方法如下：

设a＝[a₁，a₂，...，a_n]^T、b＝[b₁，b₂，...，b_n]^T和c＝[c₁，c₂，...，c_n]^T均为n×1阶向量。借鉴MTALB向量乘法的表示方法，定义向量点乘运算(.*)如下：

c＝a.*b＝[a₁b₁，a₂b₂，...，a_nb_n]^T                (3)

设n为电网节点总数，m为PV节点数，为节点电压列向量，其中 和f＝[f₁，f₂，...，f_n]^T分别为电压列向量的实部和虚部(依次按照PQ节点、PV节点和平衡节点顺序排列)，Y＝G+jB为n×n阶节点导纳矩阵。

节点注入复功率向量

为：

$\stackrel{\·}{S}=\stackrel{\·}{U}.*\stackrel{\‾}{\left(Y\stackrel{\·}{U}\right)}$

$=(e+\mathrm{jf}).*\left[(G-\mathrm{jB})(e-\mathrm{jf})\right]---\left(4\right)$

$=[e.*(\mathrm{Ge}-\mathrm{Bf})+f.*(\mathrm{Gf}+\mathrm{Be})]+$

$j[-e.*(\mathrm{Gf}+\mathrm{Be})+f.*(\mathrm{Ge}-\mathrm{Bf})]$

由式(4)易知，n个节点注入有功和无功向量分别为：

P＝e.*(Ge-Bf)+f.*(Gf+Be)                         (5)

Q＝-e.*(Gf+Be)+f.*(Ge-Bf)                  (6)

令S_R＝diag(P)、S_I＝diag(Q)、E＝diag(e)和F＝diag(f)，其中diag(.)表示把向量转化为对角矩阵。设S1_R、S1_I、H′、N′、M′和L′为n×n阶中间矩阵，则有：

S1_R＝E(GE-BF)+F(GF+BE)                    (7)

S1_I＝-E(GF+BE)+F(GE-BF)                   (8)

H′＝S_I-S1_I                              (9)

N′＝S_R+S1_R                              (10)

M′＝S_R-S1_R                              (11)

L′＝S_I+S1_I                              (12)

当潮流方程用极坐标形式表示时，其雅可比矩阵为：

$J=\left[\begin{array}{cc}\frac{\∂P}{\∂\mathrm{\θ}}& \frac{\∂P}{\∂U}U\\ \frac{\∂Q}{\∂\mathrm{\θ}}& \frac{\∂Q}{\∂U}U\end{array}\right]=\left[\begin{array}{cc}H& N\\ W& L\end{array}\right]---\left(13\right)$

取n×n阶矩阵H′的前n-1行和前n-1列元素即得到矩阵H；取n×n阶矩阵N′的前n-1行和前n-m-1列元素即得到矩阵N；取n×n阶矩阵W′的前n-m-1行和前n-1列元素即得到矩阵W；取n×n阶矩阵L′的前n-m-1行和前n-m-1列元素即得到矩阵L，从而得到整个雅可比矩阵J。

得到上述雅克比矩阵的矢量化计算方法后，通过并行矩阵、向量运算实现雅克比矩阵的并行化求解。

(3)并行求解功率不平衡量

第(2)步完成后，进行并行求解有功功率不平衡量ΔP和无功功率不平衡量ΔQ，即将电力网络各节点发电机注入功率减去负荷功率。设注入功率的矩阵形式为：P_注和Q_注，用它们减去式(5)和式(6)中的有功、无功向量得到功率不平衡量为ΔP＝P_注-P和ΔQ＝Q_注-Q。

(4)生成并行预条件子矩阵

第(3)步完成后，需要生成广义极小残量(GMRES，Generalized Minimal Residual)方法求解线性方程组时并行预条件子，首先将雅克比矩阵表示为如下形式：

其中J_ij为n_i×n_j阶矩阵，且n₁+n₂+…+n_m＝n，则称准对角矩阵M＝diag(A₁₁，A₂₂，…，A_mm)为矩阵A的块雅克比预条件子，结合准对角矩阵M的结构特点，即子矩阵A₁₁，A₂₂，…，A_mm分别占据系数矩阵A独立的行和列，在并行计算过程中根据处理器数确定块雅克比预条件子的分块数，得并行块雅克比预条件子M：

(5)并行求解修正方程组

第(4)步完成后，采用与准对角并行预条件子相结合广义极小残量(GMRES，GeneralizedMinimal Residual)方法求解牛顿潮流修正方程组：

Jy＝b                        (16)

式中，J为雅克比矩阵，

为各未知节点电压相角和幅值的修正量，

为功率不平衡量。

式(16)对应的修正方程组的详细求解步骤如下(各步均通过并行矩阵、向量运算实现)：

①初始化y₀，计算r₀＝M^-1(b-Jy₀)(M为并行预条件子)，β＝||r₀||₂和v₁＝r₀/β；

②定义V_m：＝[v₁，…，v_m]，

③Arnoldi过程：

外循环开始 j＝1，...，m执行：

计算ω：＝M^-1Jv_j

内循环开始 i＝1，...，j执行：

h_ij：＝(ω，v_i)

ω：＝ω-h_ijv_i

内循环结束

计算h_(j+1)j＝||ω||₂和v_j+1＝ω/h_(j+1)j

外循环结束

④计算z_m∈R^m，使得y_m＝y₀+V_mz_m；

⑤收敛判断：如果满足收敛条件，结束；否则，令x₀＝x_m，转①重新进行循环迭代。

(6)牛顿迭代收敛判断和潮流计算结果输出

第(5)步完成后，令x_new＝x_old+y，并检查范数是否满足||x_new-x_old||＜ε：当条件满足，输出潮流计算结果，计算结束；否则返回第(2)步重新求解雅克比矩阵，再进行计算再判断，直到||x_new-x_old||＜ε条件满足为止。

本发明采用上述技术方案后，主要有以下效果：

①由于本发明基于线性方程组的并行求解和雅克比矩阵更新的并行化技术来实现整个牛顿潮流的并行求解，在整个潮流求解过程中不需要任何形式的网络分割，计算简单、通用性较好，便于推广应用，特别适合各种大规模复杂电力系统并行潮流计算；

②本发明基于预条件广义极小残量方法并行求解牛顿潮流迭代时的修正方程组，彻底摆脱传统线性方程组求解方法，如LU分解法、Gauss迭代法等，算法内部不存在前推回代的依赖关系，能实现粗粒度并行计算，并行效率较高，计算速度快；

③由于本发明涉及的并行计算过程均基于消息传递并行模式，编写出的计算机程序移植性较好，适用于常见的并行计算机、并行集群等，便于推广应用。

本发明可广泛应用于电力系统并行潮流计算中，特别适用于大规模复杂电力系统牛顿潮流并行计算。

附图说明

图1为本发明的程序流程框图。

具体实施方式

下面结合具体实施方式，进一步说明本发明。

实施例

将本发明应用到IEEE30节点、IEEE118节点标准系统及由它们合成的5个大规模电力系统。合成系统将IEEE标准系统按N×N网格结构通过高压联络线进行互联，从而得到更大规模的测试系统。合成系统的节点导纳矩阵与同等规模的实际输电网络具有接近的条件数和特征值分布。测试采用的硬件系统是由曙光TC1700服务器和Cisco千兆以太网交换机搭建的具有若干处理节点的分布式计算平台，节点处理器为Xeon 2.0GHZ/512KB，具有1G内存存储空间，操作系统为Linux v9.0，节点间的并行通信接口为：MPICH-1.2.5。

如图1所示，一种大规模电力系统牛顿潮流的并行计算方法，其具体步骤如下：

(1)初始化

输入电力系统网络结构参数和电气参数，包括节点导纳矩阵；PQ节点的注入有功功率、无功功率；PV节点的注入有功功率、负荷无功功率、节点电压；对未知变量赋初值：x₀＝x_old，其中未知电压幅值初值设为1.0，未知电压相角初值设为0.0；设置牛顿潮流迭代收敛精度ε＝0.0001。

(2)并行雅克比矩阵更新

第(1)步完成后，首先矢量化雅克比矩阵的更新过程，然后并行求解。具体方法如下：

设a＝[a₁，a₂，...，a_n]^T、b＝[b₁，b₂，...，b_n]^T和c＝[c₁，c₂，...，c_n]^T均为n×1阶向量。借鉴MTALB向量乘法的表示方法，定义向量点乘运算(.*)如下：

c＝a.*b＝[a₁b₁，a₂b₂，...，a_nb_n]^T            (1)

设n为电网节点总数，m为PV节点数，为节点电压列向量，其中

e＝[e₁，e₂，...，e_n]^T和f＝[f₁，f₂，...，f_n]^T分别为电压列向量的实部和虚部(依次按照PQ节点、PV节点和平衡节点顺序排列)，Y＝G+jB为n×n阶节点导纳矩阵。

节点注入复功率向量

为：

$\stackrel{\·}{S}=\stackrel{\·}{U}.*\stackrel{\‾}{\left(Y\stackrel{\·}{U}\right)}$

$=(e+\mathrm{jf}).*\left[(G-\mathrm{jB})(e-\mathrm{jf})\right]---\left(2\right)$

$=[e.*(\mathrm{Ge}-\mathrm{Bf})+f.*(\mathrm{Gf}+\mathrm{Be})]+$

$j[-e.*(\mathrm{Gf}+\mathrm{Be})+f.*(\mathrm{Ge}-\mathrm{Bf})]$

由式(2)易知，n个节点注入有功和无功向量分别为：

P＝e.*(Ge-Bf)+f.*(Gf+Be)                (3)

Q＝-e.*(Gf+Be)+f.*(Ge-Bf)               (4)

令S_R＝diag(P)、S_I＝diag(Q)、E＝diag(e)和F＝diag(f)，其中diag(.)表示把向量转化为对角矩阵。设S1_R、S1_I、H′、N′、M′和L′为n×n阶中间矩阵，则有：

S1_R＝E(GE-BF)+F(GF+BE)                 (5)

S1_I＝-E(GF+BE)+F(GE-BF)                (6)

H′＝S_I-S1_I                           (7)

N′＝S_R+S1_R                           (8)

M′＝S_R-S1_R                           (9)

L′＝S_I+S1_I                    (10)

当潮流方程用极坐标形式表示时，其雅可比矩阵为：

$J=\left[\begin{array}{cc}\frac{\∂P}{\∂\mathrm{\θ}}& \frac{\∂P}{\∂U}U\\ \frac{\∂Q}{\∂\mathrm{\θ}}& \frac{\∂Q}{\∂U}U\end{array}\right]=\left[\begin{array}{cc}H& N\\ W& L\end{array}\right]---\left(11\right)$

取n×n阶矩阵H′的前n-1行和前n-1列元素即得到矩阵H；取n×n阶矩阵N′的前n-1行和前n-m-1列元素即得到矩阵N；取n×n阶矩阵W′的前n-m-1行和前n-1列元素即得到矩阵W；取n×n阶矩阵L′的前n-m-1行和前n-m-1列元素即得到矩阵L，从而得到整个雅可比矩阵J。

得到上述雅克比矩阵的矢量化计算方法后，以通过并行矩阵、向量运算实现雅克比矩阵的并行化求解。

(3)并行求解功率不平衡量

第(2)步完成后，进行并行求解有功功率不平衡量ΔP和无功功率不平衡量ΔQ，即将电力网络各节点发电机注入功率减去负荷功率。设注入功率的矩阵形式为：P_注和Q_注，用它们减去式(3)和式(4)中的有功、无功向量得到功率不平衡量为ΔP＝P_注-P和ΔQ＝Q_注-Q。

(4)生成并行预条件子矩阵

第(3)步完成后，需要生成广义极小残量(GMRES，Generalized Minimal Residual)方法求解线性方程组时并行预条件子，首先将雅克比矩阵表示为如下形式：

其中J_ij为n_i×n_j阶矩阵，且n₁+n₂+…+n_m＝n，则称准对角矩阵M＝diag(A₁₁，A₂₂，…，A_mm)为矩阵A的块雅克比预条件子，结合准对角矩阵M的结构特点，即子矩阵A₁₁，A₂₂，…，A_mm分别占据系数矩阵A独立的行和列，在并行计算过程中根据处理器数确定块雅克比预条件子的分块数，得并行块雅克比预条件子M：

(5)并行求解修正方程组

第(4)步完成后，采用与准对角并行预条件子相结合广义极小残量(GMRES，GeneralizedMinimal Residual)方法求解牛顿潮流修正方程组：

Jy＝b                        (14)

式中，J为雅克比矩阵，

为各未知节点电压相角和幅值的修正量，

为功率不平衡量。

式(14)对应的修正方程组的详细求解步骤如下(各步均通过并行矩阵、向量运算实现)：

①初始化y₀，计算r₀＝M^-1(b-Jy₀)(M为并行预条件子)，β＝||r₀||₂和v₁＝r₀/β；

②定义V_m：＝[v₁，…，v_m]，

③Arnoldi过程：

外循环开始 j＝1，...，m执行：

计算ω：＝M^-1Jv_j

内循环开始 i＝1，...，j执行：

h_ij：＝(ω，v_i)

ω：＝ω-h_ijv_i

内循环结束

计算h_(j+1)j＝||ω||₂和v_j+1＝ω/h_(j+1)j

外循环结束

④计算z_m∈R^m，使得

y_m＝y₀+V_mz_m；

⑤收敛判断：如果满足收敛条件，结束；否则，令x₀＝x_m，转①重新进行循环迭代。

(6)牛顿迭代收敛判断和潮流计算结果输出

第(5)步完成后，令x_new＝x_old+y，并检查范数是否满足||x_new-x_old||＜ε：当条件满足，输出潮流计算结果，计算结束；否则返回第(2)步重新求解雅克比矩阵，再进行计算再判断，直到||x_new-x_old||＜ε条件满足为止。

基于本发明的并行牛顿潮流计算性能如下表所示：

为了说明本发明的并行性能，基于LU分解法的并行牛顿潮流计算性能如下表所示：

从上述结果可知，运用本发明进行并行潮流计算时，当系统规模超过3000节点时，与并行LU分解法潮流计算相比，本发明对应的并行潮流计算方法在程序执行时间、加速比和效率方面均优于前者，并且随着计算规模的增大，两者之间的差别增大。对于12000节点系统，采用2处理器、4处理器和8处理器的执行时间分别比前者减少50.7％、53.0％和34.2％，并行效果较好。

Claims (1)

1.一种大规模电力系统牛顿潮流的并行计算方法，利用计算机，根据并行集群，通过程序，结合雅克比矩阵更新的并行化技术，采用与准对角并行预条件子相结合广义极小残量方法求解牛顿潮流的修正方程组，实现大规模电力系统潮流求解过程的并行化，其具体方法步骤如下：

(1)初始化

输入电力系统网络结构参数和电气参数，包括节点导纳矩阵；PQ节点的注入有功功率、无功功率；PV节点的注入有功功率、负荷无功功率、节点电压，并对未知变量赋初值：x₀＝x_old；

(2)并行雅克比矩阵更新

第(1)步完成后，首先矢量化雅克比矩阵的更新过程，然后并行求解，具体方法如下：

设a＝[a₁，a₂，...，a_n]^T、b＝[b₁，b₂，...，b_n]^T和c＝[c₁，c₂，...，c_n]^T均为n×1阶向量，借鉴MTALB向量乘法的表示方法，定义向量点乘运算如下：

c＝a.*b＝[a₁b₁，a₂b₂，...，a_nb_n]^T                    (1)

设n为电网节点总数，m为PV节点数，

为节点电压列向量，其中e＝[e₁，e₂，...，e_n]^T和f＝[f₁，f₂，...，f_n]^T分别为电压列向量的实部和虚部，Y＝G+jB为n×n阶节点导纳矩阵，

节点注入复功率向量

为：

$\stackrel{\·}{S}=\stackrel{\·}{U}.*\stackrel{\‾}{\left(Y\stackrel{\·}{U}\right)}$

$=(e+\mathrm{jf}).*\left[(G-\mathrm{jB})(e-\mathrm{jf})\right]---\left(2\right)$

$=[e.*(\mathrm{Ge}-\mathrm{Bf})+f.*(\mathrm{Gf}+\mathrm{Be})]+$

$j[-e.*(\mathrm{Gf}+\mathrm{Be})+f.*(\mathrm{Ge}-\mathrm{Bf})]$

由式(2)易知，n个节点注入有功和无功向量分别为：

P＝e.*(Ge-Bf)+f.*(Gf+Be)                             (3)

Q＝-e.*(Gf+Be)+f.*(Ge-Bf)                            (4)

令S_R＝diag(P)、S_I＝diag(Q)、E＝diag(e)和F＝diag(f)，设S1_R、S1_I、H′、N′、M′和L′为n×n阶中间矩阵，则有：

S1_R＝E(GE-BF)+F(GF+BE)                              (5)

S1_I＝-E(GF+BE)+F(GE-BF)                             (6)

H′＝S_I-S1_I                                        (7)

N′＝S_R+S1_R                                    (8)

M′＝S_R-S1_R                                    (9)

L′＝S_I+S1_I                                    (10)

当潮流方程用极坐标形式表示时，其雅可比矩阵为：

$J=\left[\begin{array}{cc}\frac{\∂P}{\∂\mathrm{\θ}}& \frac{\∂P}{\∂U}U\\ \frac{\∂Q}{\∂\mathrm{\θ}}& \frac{\∂Q}{\∂U}U\end{array}\right]=\left[\begin{array}{cc}H& N\\ W& L\end{array}\right]---\left(11\right)$

取n×n阶矩阵H′的前n-1行和前n-1列元素即可得到矩阵H；取n×n阶矩阵N′的前n-1行和前n-m-1列元素即可得到矩阵N；取n×n阶矩阵W′的前n-m-1行和前n-1列元素即可得到矩阵W；取n×n阶矩阵L′的前n-m-1行和前n-m-1列元素即可得到矩阵L，从而得到整个雅可比矩阵J，得到上述雅克比矩阵的矢量化计算方法后，可以通过并行矩阵、向量运算实现雅克比矩阵的并行化求解；

(3)并行求解功率不平衡量

第(2)步完成后，进行并行求解有功功率不平衡量ΔP和无功功率不平衡量ΔQ，即将电力网络各节点发电机注入功率减去负荷功率，设注入功率的矩阵形式为：P_注和Q_注，用它们减去式(3)和式(4)中的有功、无功向量可得到功率不平衡量为ΔP＝P_注-P和ΔQ＝Q_注-Q；

(4)生成并行预条件子矩阵

第(3)步完成后，需要生成广义极小残量方法求解线性方程组时并行预条件子，首先将雅克比矩阵表示为如下形式：

其中J_ij为n_i×n_j阶矩阵，且n₁+n₂+…+n_m＝n，则称准对角矩阵M＝diag(J₁₁，J₂₂，…，J_mm)为雅克比矩阵J的块雅克比预条件子，结合准对角矩阵M的结构特点，即子矩阵J₁₁，J₂₂，…，J_mm分别占据雅克比矩阵J独立的行和列，在并行计算过程中根据处理器数确定块雅克比预条件子的分块数，可得并行块雅克比预条件子M：

式(13)中，子矩阵J₁₁，J₂₂，…，J_mm分别属于处理器1，处理器2，...，处理器m；

(5)并行求解修正方程组

第(4)步完成后，采用与准对角并行预条件子相结合广义极小残量方法求解牛顿潮流修正方程组：

Jy＝b                                (14)

式中，J为雅克比矩阵， $y=\left[\frac{\stackrel{\mathrm{\Δ\θ}}{\mathrm{\ΔU}}}{U}\right]$ 为各未知节点电压相角和幅值的修正量， $b=\left[\begin{array}{c}\mathrm{\ΔP}\\ \mathrm{\ΔQ}\end{array}\right]$ 为功率不平衡量，式(14)对应的修正方程组的详细求解步骤如下：

①初始化y₀，计算r₀＝M^-1(b-Jy₀)，β＝||r₀||₂和v₁＝r₀/β；

②定义V_m：＝[v₁，…，v_m]， ${\stackrel{\‾}{H}}_m={\left\{h_{\mathrm{ij}}\right\}}_{1\≤i\≤j+\mathrm{1,1}\≤j\≤m};$

③Arnoldi过程：

外循环开始j＝1，...，m执行：

计算ω：＝M^-1Jv_j

内循环开始i＝1，...，j执行：

h_ij：＝(ω，v_i)

ω：＝ω-h_ijv_i

内循环结束

计算h_(j+1)j＝||ω||₂和v_j+1＝ω/h_(j+1)j

外循环结束

④计算z_m∈R^m，使得

y_m＝y₀+V_mz_m；

⑤收敛判断：如果满足收敛条件，结束；否则，令y₀＝y_m，转①重新进行循环迭代；

(6)牛顿迭代收敛判断和潮流计算结果输出

第(5)步完成后，令x_new＝x_old+y_m，并检查范数是否满足||x_new-x_old||＜ε：当条件满足，输出潮流计算结果，计算结束；否则返回第(2)步重新求解雅克比矩阵，再进行计算再判断，直到||x_new-x_old||＜ε条件满足为止。

Images