CN106712029A - 小阻抗支路pq端点变雅可比矩阵的牛顿法潮流计算方法 - Google Patents

Abstract

本发明公开了一种小阻抗支路PQ端点变雅可比矩阵的牛顿法潮流计算方法，包括以下步骤：根据支路电阻和电抗的大小确定两端节点所连支路类型T；小阻抗支路的PQ端点与其他节点采用不同的方法计算雅可比矩阵；解修正方程。本发明通过在迭代过程中小阻抗支路的PQ端点采用与其他节点不同的雅可比矩阵计算方法，解决了直角坐标牛顿法潮流计算在分析含有小阻抗支路电力系统时的收敛性问题。采用常规直角坐标牛顿法潮流计算不收敛时，本发明能够可靠收敛，且比现有专利技术迭代次数少。由于本发明不仅能有效解决常规直角坐标牛顿法潮流计算分析含有小阻抗支路电力系统的收敛性问题，同时也能对正常电力系统进行潮流计算，没有不良影响。

Description

小阻抗支路PQ端点变雅可比矩阵的牛顿法潮流计算方法

技术领域

本发明涉及含小阻抗支路电力系统的潮流计算方法，特别是一种电力系统的直角坐标牛顿法潮流计算方法。

背景技术

电力系统潮流计算是研究电力系统稳态运行的一项基本计算，它根据电力系统给定的运行条件和网络结构确定整个电力系统的运行状态。潮流计算也是电力系统其他分析的基础，如安全分析、暂态稳定分析等都要用到潮流计算。由于具有收敛可靠、计算速度较快及内存需求适中的优点，牛顿法成为当前潮流计算的主流方法。牛顿法分为极坐标和直角坐标两种形式，两种形式的牛顿法潮流计算都在电力系统中得到了广泛的应用。

在直角坐标牛顿法潮流计算中，节点i的电压采用直角坐标表示为：

对正常电力网络，直角坐标牛顿法潮流计算具有良好的收敛性，但遇到含有小阻抗支路的病态网络时，直角坐标牛顿法潮流计算就可能发散。电力系统小阻抗支路可分为小阻抗线路和小阻抗变压器支路，在数学模型上线路可以看作变比为1:1的变压器，因此下面分析时仅以小阻抗变压器支路为例分析。小阻抗变压器模型见图1，变压器的非标准变比k位于节点i侧，阻抗位于标准变比侧。变压器阻抗z_ij＝r_ij+jx_ij很小，导纳为

式中，y_ij、g_ij、b_ij分别为节点i和节点j之间小阻抗支路的导纳、电导和电纳；r_ij、x_ij分别为节点i和节点j之间小阻抗支路的电阻和电抗。

由于小阻抗支路l_i-j的阻抗很小，支路的电压降也很小，因此变压器两端节点的电压应满足：

如图2所示，现有直角坐标牛顿法潮流计算方法，主要包括以下步骤：

A、输入原始数据和初始化电压

根据电力系统节点的特点，潮流计算把电力系统节点分成3类：节点有功功率和无功功率已知、节点电压幅值和电压相角未知的节点称为PQ节点；节点有功功率和电压幅值已知、节点无功功率和电压相角未知的节点称为PV节点；节点电压幅值和电压相角已知，节点有功功率和无功功率未知的节点称为平衡节点。

电压初始化采用平启动，即PV节点和平衡节点的电压实部取给定值，PQ节点的电压实部取1.0；所有电压的虚部都取0.0。这里单位采用标幺值。

B、形成节点导纳矩阵

设节点i和节点j原来的自电导与自电纳分别为G_i0、B_i0、G_j0、B_j0，在它们之间增加一条小阻抗支路后的自导纳和互导纳分别为：

式中，Y_ii、Y_jj分别为节点i和节点j的自导纳；Y_ij为节点i和节点j之间的互导纳；r_ij、x_ij分别为节点i和节点j之间小阻抗支路的电阻和电抗；k为节点i和节点j之间小阻抗支路的变比(如果是输电线支路，变比为1)；

C、计算功率及电压偏差

PQ节点的功率偏差计算公式为：

式中，P_is、Q_is分别为节点i给定的注入有功功率和无功功率，P_is为电源有功功率与负荷有功功率之差，Q_is为电源无功功率与负荷无功功率之差；a_i、b_i分别为节点i的计算注入电流相量的实部和虚部，为

式中，n为电力系统的节点数。

PV节点的有功功率及电压偏差计算公式为：

式中，V_is为节点i给定的电压幅值。

平衡节点不参与迭代计算，不需要计算功率偏差或电压偏差。

求各节点功率或电压偏差中绝对值最大的值，称为最大不平衡量，如果最大不平衡量的绝对值小于给定的收敛精度，转步骤F，否则执行步骤D。

D、形成雅可比矩阵J

雅可比矩阵J的元素(i≠j时)计算公式如下：

雅可比矩阵J的元素(i＝j时)计算公式如下：

PQ节点按式(15)-(18)计算雅可比矩阵元素；PV节点按式(15)、(16)、(19)、(20)计算雅可比矩阵元素。

平衡节点不计算雅可比矩阵元素。

E、解修正方程及修正电压实部e、虚部f

潮流计算的基本方程(6)和(8)是非线性方程组，通常采用逐次线性化方法迭代求解。线性化得到的方程称为修正方程，用来求电压实部和虚部的修正量。

修正方程为：

式中，J为雅可比矩阵；ΔP和ΔQ分别为有功功率和无功功率偏差列向量；ΔV²为电压幅值偏差列向量；Δe和Δf分别为电压相量的实部和虚部修正量列向量；为有功功率偏差函数列向量对电压相量实部列向量转置的偏导矩阵，上标T为转置符号。

电压修正公式为：

式中，上标t表示第t次迭代。

F、输出节点及支路数据。

对正常电力网络，牛顿法潮流计算具有良好的收敛性，但遇到含有小阻抗支路的病态网络时，牛顿法潮流计算就可能发散。而电力系统中小阻抗支路普遍存在，收敛性是电力系统潮流计算这类非线性问题的最重要指标，计算不收敛就无法得到方程的解。因此改善直角坐标牛顿法潮流计算针对含有小阻抗支路电力系统的收敛性具有非常重要的意义。

中国专利ZL201410299531.5披露了一种通过修改常规直角坐标牛顿法潮流计算雅可比矩阵的方法，改善了潮流计算的收敛性。该方法计算雅可比元素时采用由给定值P_is和Q_is计算的a_i和b_i值，有效解决了含有电阻为0的小阻抗支路电力系统潮流计算的发散问题。但当小阻抗支路的电阻不为0时，该方法迭代次数增加，收敛性变差，甚至不收敛。

中国专利ZL201410315785.1提出了一种雅可比矩阵改变的直角坐标牛顿法潮流计算方法，该方法首次迭代和后续各次迭代采用不同的雅可比矩阵计算方法，首次迭代计算雅可比元素时采用由给定值P_is和Q_is计算的a_i和b_i值，后续各次迭代计算雅可比元素时仍采用传统方法，有效解决了含电阻不为0的小阻抗支路电力系统潮流计算的发散问题，但该方法在电力系统包含多条电阻不为0的小阻抗支路时，迭代次数增加，收敛性变差。

中国专利ZL201611094297.8提出了一种随迭代和节点类型改变雅可比矩阵的牛顿法潮流计算方法，该方法首次迭代时所有PQ节点与后续各次迭代采用不同的雅可比矩阵计算方法，首次迭代时所有PQ节点计算雅可比元素时采用由给定值P_is和Q_is计算的a_i和b_i值，首次迭代时所有PV节点及后续各次迭代时所有节点计算雅可比元素时仍采用传统方法，有效解决了电力系统包含多条电阻不为0的小阻抗支路的潮流计算的发散问题，但迭代次数仍然较多，需要进一步改进。

发明内容

为解决现有技术存在的上述问题，本发明要提出一种小阻抗支路PQ端点变雅可比矩阵的牛顿法潮流计算方法，该方法可以提高其分析含有电阻不为0的小阻抗支路电力系统的收敛速度。

为了实现上述目的，本发明提出了一种直角坐标牛顿法潮流计算方法来改善潮流计算收敛性。本发明的小阻抗支路的PQ端点采用由给定值P_is和Q_is计算的a_i和b_i值计算雅可比矩阵元素，其他节点则采用传统方法计算雅可比矩阵元素。本发明的技术方案如下：小阻抗支路PQ端点变雅可比矩阵的牛顿法潮流计算方法，包括以下步骤：

A、输入原始数据和初始化电压；

B、根据支路电阻和电抗的大小确定两端节点所连支路类型T

形成节点所连支路类型数组的具体步骤如下：

B1、读入支路数据，设置小电阻阈值r_min和小电抗阈值x_min；

B2、节点所连支路类型数组T清零；

B3、令m＝1；

B4、取支路m的首末节点号i和j、电阻r、电抗x；

B5、判断是否满足r≤r_min且x≤x_min的条件，如果不满足，转步骤B7；

B6、令T_i＝1，T_j＝1；

B7、令m＝m+1；

B8、判断m是否大于支路数l，如果m不大于l转步骤B4；否则转步骤C；

C、形成节点导纳矩阵；

D、设置迭代计数t＝0；

E、计算功率及电压偏差，求最大不平衡量ΔW_max；

F、判断最大不平衡量绝对值|ΔW_max|是否小于收敛精度ε；如果小于收敛精度ε，执行步骤J；否则，执行步骤G；

G、形成雅可比矩阵；

雅可比矩阵元素计算方法根据支路阻抗值和节点类型采用不同方法：

潮流计算按平启动设置电压初值，正常电力系统按此电压初值计算的节点功率P_i和Q_i与实际值相差不大；当电力系统为包含小阻抗支路的病态系统时，小阻抗支路两端电压初值可能在其等值阻抗上产生较大的电压差，会计算出很大的支路功率，进而在两端节点得到与实际值相差很大的节点注入功率P_i和Q_i以及与实际值相差很大的注入电流相量的实部和虚部a_i、b_i，以此a_i、b_i计算雅可比矩阵元素，将导致潮流计算发散。

潮流计算收敛时，式(6)中ΔP_i、ΔQ_i都趋近于0，节点注入功率计算值等于节点注入功率给定值，用节点注入功率给定值计算注入电流相量的实部和虚部,记为a_is和b_is,计算公式为：

对于小阻抗支路PQ端点来说，a_is和b_is会比直接采用式(7)计算的a_i和b_i更接近实际值，更有利于改善潮流计算的收敛性。由于PV节点的电源无功功率不是给定的，输入时输入任意值，用此任意值计算的a_is和b_is，在潮流计算的后几次迭代中则与实际值差值较大，不如a_i和b_i准确，不利于提高潮流计算的收敛速度；同时后文的证明表明小阻抗支路PV端点采用传统方法计算雅可比矩阵元素，不会导致潮流计算发散。因此小阻抗支路PV端点仍然采用传统方法计算雅可比矩阵元素，即雅可比矩阵计算公式中的节点i注入电流相量的实部a_i和虚部b_i按式(7)计算。正常支路的端点也采用传统方法计算雅可比矩阵元素。

形成雅可比矩阵元素的具体步骤如下：

G1、按式(9)-(14)计算i≠j时的雅可比矩阵元素；

G2、令i＝1；

G3、判断是否同时满足(T_i＝1)和(节点i是PQ节点)的条件，如果不满足转步骤G4；如果满足，则按修改公式计算i＝j时的雅可比矩阵元素，然后转步骤G5；

计算i＝j时的雅可比矩阵元素的修改公式为：

G4、按式(15)-(20)计算i＝j时的雅可比矩阵元素；

PQ节点按式(15)-(18)计算雅可比矩阵元素；PV节点按式(15)、(16)、(19)、(20)计算雅可比矩阵元素；平衡节点不计算雅可比矩阵元素；

G5、令i＝i+1；

G6、判断i是否大于节点数n，如果i不大于n转步骤G3；否则转步骤H；

H、解修正方程及修正电压实部e、虚部f；

I、令t＝t+1，返回步骤E进行下一次迭代；

J、输出节点及支路数据。

本发明对直角坐标牛顿法潮流计算的雅可比矩阵的部分元素(i＝j时)的计算时，仅当小阻抗支路的端点是PQ节点时采用修改公式，端点是PV节点时仍然采用传统方法，其他非小阻抗支路的端点也采用传统方法。针对小阻抗支路的一个端点是PV节点，另一个是PQ节点的收敛性证明如下：

设节点i是PQ节点，节点j是PV节点。

与小阻抗支路有关的修正方程为：

式中，A_i、A_j、B_i为与Δe_k、Δf_k相关的项(k＝1,…,n且k≠i,j)。

功率及电压偏差公式为：

式中，P_i0、P_j0、Q_i0为除小阻抗支路l_i-j外节点的计算功率。

节点i按式(9)-(14)及式(24)-(27)计算雅可比矩阵元素；节点j按式(9)-(16)、(19)、(20)计算雅可比矩阵元素。

式(7)的展开式为：

式中，a_i0和b_i0分别为与节点i相连的其他非小阻抗支路计算出的注入电流相量的实部和虚部。

把式(36)代入式(15)、(16)后，节点j对应的雅可比矩阵元素(i＝j时)为

代入雅可比矩阵元素和功率偏差及电压偏差后，与小阻抗支路有关的修正方程为：

式(41)-(43)中，小阻抗支路r_ij、x_ij很小，因而g_ij、b_ij很大。与g_ij、b_ij相比，式中A_i、A_j、a_is、b_is、a_j0、b_j0、B_i、G_i0、G_j0、B_i0、P_is、P_js、Q_is、P_i0、P_j0、Q_j0各项较小，忽略这些较小量，得：

下面分析首次迭代的情况。式(45)-(48)考虑到首次迭代时，电压为电压初值，即e_i＝1.0，f_i＝0.0，e_j＝V_js，f_j＝0.0，得：

-(g_ij/k²)Δe_i+(g_ij/k)Δe_j+(b_ij/k²)Δf_i-(b_ij/k)Δf_j≈-g_ij/k²+g_ijV_js/k (49)

(b_ij/k²)Δe_i-(b_ij/k)Δe_j+(g_ij/k²)Δf_i-(g_ij/k)Δf_j≈b_ij/k²-b_ijV_js/k (51)

-2e_jΔe_j＝0 (52)

由式(52)，得Δe_j＝0，即

式(49)乘以b_ij与式(51)乘以g_ij相加，得

式(53)中由于得

Δf_i≈kΔf_j (54)

由于初值则电压虚部修正后满足公式(2)。

把式(52)和式(54)代入式(49)，得

Δe_i≈1-kV_jS (55)

式(55)中，考虑电压实部初值首次迭代后电压实部为

式(56)满足公式(2)。

式(41)考虑到首次迭代时，电压为电压初值，并考虑Δe_j＝0，得

把式(54)和式(55)代入式(57)，得

(-a_is-G_i0-g_ij/k²)(1-kV_js)+(B_i0-b_is)Δf_i+A_i≈P_is-(G_i0+g_ij/k²)+g_ijV_js/k-P_i0 (58)

式(58)整理，得

(a_is+G_i0)(kV_js)-a_is+(B_i0-b_is)Δf_i+A_i＝P_is+g_ijV_js/k-P_i0 (59)

式(59)不存在因小阻抗出现较大的Δf_i问题，同时从式(52)和式(55)可见，电压实部的修正量也很小。且满足小阻抗支路两端电压关系式(2)。由于小阻抗的影响已经不存在了，因此首次迭代时小阻抗不会对收敛有影响。

第2次迭代时，式(41)-(44)中电压为上次迭代得到电压，满足 式(41)-(44)代入上述电压关系，得

式(60)-(62)整理，得

式(63)-(65)中，忽略相对g_ij、b_ij较小的量，得：

(-g_ije_j-b_ijf_j)Δe_i/k+(g_ije_j+b_ijf_j)Δe_j+(b_ije_j-g_ijf_j)Δf_i/k+(g_ijf_j-b_ije_j)Δf_j≈0 (66)

(-g_ije_j-b_ijf_j)Δe_j+(g_ije_j+b_ijf_j)Δe_i/k+(b_ije_j-g_ijf_j)Δf_j+(g_ijf_j-b_ije_j)Δf_i/k≈0 (67)

(b_ije_j-g_ijf_j)Δe_i/k+(g_ijf_j-b_ije_j)Δe_j+(g_ije_j+b_ijf_j)Δf_i/k-(g_ije_j+b_ijf_j)Δf_j≈0 (68)

式(66)乘以b_ij与式(68)乘以g_ij相加，得

式(70)中由于得

-f_jΔe_i/k+f_jΔe_j+e_jΔf_i/k-e_jΔf_j≈0 (71)

式(68)乘以b_ij，再与式(66)乘以g_ij相减，得

式(72)中由于得

e_jΔe_i/k-e_jΔe_j+f_jΔf_i/k-f_jΔf_j≈0 (73)

式(71)乘以e_j与式(73)乘以f_j相加，得

式(74)中由于得

Δf_i≈kΔf_j (75)

由于首次迭代后有则修正后满足公式(2)。

式(75)代入式(73)，得

Δe_i≈kΔe_j (76)

由于首次迭代后有则修正后满足公式(2)。

式(63)加式(64)，得

这样式(63)-(65)和式(69)经过变换得到式(75)-(77)和式(69)，而式(75)-(77)和式(69)已经不存在小阻抗了，且满足小阻抗支路两端电压关系式(2)。由于小阻抗的影响已经不存在了，因此第2次迭代时小阻抗不会对收敛有影响。

同理可证第2次后各次迭代时小阻抗不会对收敛有影响。

与现有技术相比，本发明具有以下有益效果：

1、本发明的雅可比矩阵元素计算方法根据支路阻抗值和节点类型采用不同方法。如果小阻抗的端点是PQ节点，因潮流计算收敛时，计算功率P_i和Q_i与给定值P_is和Q_is相等，故采用a_is和b_is来计算雅可比矩阵元素效果较好，如果采用电压初值计算出的电流相量的实部和虚部a_i和b_i来计算雅可比矩阵元素则会出现很大的值导致潮流计算发散；由于PV节点的电源无功功率不是给定的，原始数据输入时该值输入是任意值，因此Q_is与潮流计算收敛时的计算无功功率Q_i相差很大，采用a_is和b_is来计算雅可比矩阵元素效果可能较差，影响收敛速度，应该采用a_i和b_i计算，并且PV节点采用传统方法计算雅可比矩阵元素会提高潮流计算后几次迭代的收敛速度。本发明通过在迭代过程中小阻抗支路的PQ端点采用与其他节点不同的雅可比矩阵计算方法，解决了直角坐标牛顿法潮流计算在分析含有小阻抗支路电力系统时的收敛性问题。采用常规直角坐标牛顿法潮流计算不收敛时，本方法能够可靠收敛，且比现有专利技术迭代次数少。

2、由于本发明不仅能有效解决常规直角坐标牛顿法潮流计算分析含有小阻抗支路电力系统的收敛性问题，同时也能对正常电力系统进行潮流计算，没有不良影响。

附图说明

本发明共有附图6张。其中：

图1是电力系统小阻抗变压器模型示意图。

图2是直角坐标牛顿法潮流计算的流程图。

图3是专利方法1直角坐标牛顿法潮流计算的流程图。

图4是专利方法2直角坐标牛顿法潮流计算的流程图。

图5是本发明直角坐标牛顿法潮流计算的流程图。

图6是本发明形成节点所连支路类型数组的流程图。

具体实施方式

下面结合附图对本发明进行进一步地说明。根据图1所示的小阻抗变压器模型，采用图5-6所示的直角坐标牛顿法潮流计算的流程图，对一个实际大型电网进行了潮流计算。该实际大型电网有445个节点，含有大量的小阻抗支路。其中，x≤0.001的小阻抗支路有49条，x≤0.0001的小阻抗支路有41条，x≤0.00001的小阻抗支路有22条。其中阻抗值最小的是节点118和节点125之间的小阻抗支路l_118-125为x＝0.00000001，变比k＝0.9565，k位于节点118侧。潮流计算的收敛精度为0.00001。为了验证本发明计算含电阻不为0的小阻抗支路电力系统的收敛性，把小阻抗支路l_118-125、l_60-122及l_287-310的电阻改为r＝0.0001。

作为对比，同时采用以下3种对比方法对该实际大型电网进行了潮流计算：

常规方法：常规的直角坐标牛顿法潮流方法；

专利方法1：专利号为ZL201410315785.1的专利方法；

专利方法2：申请号为ZL201611094297.8的专利方法。

迭代次数结果见表1。

表1不同潮流方法的迭代结果

方法	常规方法	专利方法1	专利方法2	本专利方法
					迭代结果	不收敛	7次收敛	6次收敛	5次收敛

由表1可见，对于修改后的445节点实际电力系统算例，常规直角坐标牛顿法潮流方法不收敛，本专利方法和现有专利方法都能够收敛，但本专利方法的迭代次数比现有专利方法1少2次，比现有专利方法2少1次。

不同潮流计算方法各次迭代最大不平衡量见表2。单位为标幺值。

表2不同潮流方法各次迭代最大不平衡量

由表2可知，4种方法首次迭代前最大不平衡量相同且很大。首次迭代后，现有专利方法和本专利方法最大不平衡量明显减少，现有专利方法1迭代7次收敛；现有专利方法2迭代6次收敛；本专利方法后几次迭代最大不平衡量减少速度更快，迭代5次收敛；而常规方法的最大不平衡量则变大，最终发散。

为了验证小阻抗支路端点为PV节点时采用传统方法计算雅可比矩阵元素的有效性，作为对比方法，也用小阻抗支路的所有端点都采用修改方法计算雅可比矩阵元素的方法进行了潮流计算。各次迭代最大不平衡量结果见表3。

表3不同潮流方法各次迭代最大不平衡量

由表3可知，本专利方法和对比方法的首次迭代前最大不平衡量相同且很大，之后最大不平衡量明显减少，前2次迭代时最大不平衡量变化都较快。但小阻抗支路的所有端点都采用修改方法计算雅可比矩阵元素时，以后各次收敛速度较慢，迭代17次收敛；本专利方法最大不平衡量减少速度一直较快，迭代5次收敛。可见小阻抗支路的端点是PV节点时不宜采用修改方法计算雅可比矩阵元素。

PV节点的电源无功功率输入值与潮流计算收敛时的计算值及第3次迭代后的计算值见表4。单位为标幺值。

表4PV节点的电源无功功率输入值与收敛计算值及第3次迭代后的计算值

由表4可见，PV节点的电源无功功率输入值与潮流计算收敛时的计算值有较大的差值；第3次迭代后的计算值则与潮流计算收敛时的计算值就比较接近了，此后迭代再用PV节点的电源无功功率输入值求出给定值P_is和Q_is计算的a_i和b_i值计算雅可比元素就不合适了，影响潮流计算的收敛速度；同时前面的证明过程表明，小阻抗支路的PV端点使用传统方法计算雅可比矩阵元素不影响收敛性。

本发明可以采用任何一种编程语言和编程环境实现，如C语言、C++、FORTRAN、Delphi等。开发环境可以采用Visual C++、Borland C++Builder、Visual FORTRAN等。

本发明不局限于本实施例，任何在本发明披露的技术范围内的等同构思或者改变，均列为本发明的保护范围。

Claims (1)

1.小阻抗支路PQ端点变雅可比矩阵的牛顿法潮流计算方法，其特征在于：包括以下步骤：

A、输入原始数据和初始化电压

根据电力系统节点的特点，潮流计算把电力系统节点分成3类：节点有功功率和无功功率已知、节点电压幅值和电压相角未知的节点称为PQ节点；节点有功功率和电压幅值已知、节点无功功率和电压相角未知的节点称为PV节点；节点电压幅值和电压相角已知，节点有功功率和无功功率未知的节点称为平衡节点；

电压初始化采用平启动，即PV节点和平衡节点的电压实部取给定值，PQ节点的电压实部取1.0；所有电压的虚部都取0.0；这里单位采用标幺值；

B、根据支路电阻和电抗的大小确定两端节点所连支路类型T

形成节点所连支路类型数组的具体步骤如下：

B1、读入支路数据，设置小电阻阈值r_min和小电抗阈值x_min；

B2、节点所连支路类型数组T清零；

B3、令m＝1；

B4、取支路m的首末节点号i和j、电阻r、电抗x；

B5、判断是否满足r≤r_min且x≤x_min的条件，如果不满足，转步骤B7；

B6、令T_i＝1，T_j＝1；

B7、令m＝m+1；

B8、判断m是否大于支路数l，如果m不大于l转步骤B4；否则转步骤C；

C、形成节点导纳矩阵

设节点i和节点j原来的自电导与自电纳分别为G_i0、B_i0、G_j0、B_j0，在它们之间增加一条小阻抗支路后的自导纳和互导纳分别为：

$Y_{ii}=(G_{i0}+\frac{r_{ij}}{k^2(r_{ij}^2+x_{ij}^2)})+j(B_{i0}-\frac{x_{ij}}{k^2(r_{ij}^2+x_{ij}^2)})---\left(1\right)$

$Y_{jj}=(G_{j0}+\frac{r_{ij}}{(r_{ij}^2+x_{ij}^2)})+j(B_{j0}-\frac{x_{ij}}{(r_{ij}^2+x_{ij}^2)})---\left(2\right)$

$Y_{ij}=-\frac{r_{ij}}{k(r_{ij}^2+x_{ij}^2)}+j\frac{x_{ij}}{k(r_{ij}^2+x_{ij}^2)}---\left(3\right)$

式中，Y_ii、Y_jj分别为节点i和节点j的自导纳；Y_ij为节点i和节点j之间的互导纳；r_ij、x_ij分别为节点i和节点j之间小阻抗支路的电阻和电抗；k为节点i和节点j之间小阻抗支路的变比，如果是输电线支路，变比k为1；

D、设置迭代计数t＝0；

E、计算功率及电压偏差，求最大不平衡量ΔW_max；

PQ节点的功率偏差计算公式为：

$\left\{\begin{array}{c}{\mathrm{\ΔP}}_i=P_{is}-P_i=P_{is}-e_i{a}_i-f_i{b}_i\\ {\mathrm{\ΔQ}}_i=Q_{is}-Q_i=Q_{is}-f_i{a}_i+e_i{b}_i\end{array}\right.---\left(4\right)$

式中，P_is、Q_is分别为节点i给定的注入有功功率和无功功率，P_is为电源有功功率与负荷有功功率之差，Q_is为电源无功功率与负荷无功功率之差；e_i和f_i分别为节点i的电压相量的实部和虚部；a_i、b_i分别为节点i的计算注入电流相量的实部和虚部，为

$\left\{\begin{array}{c}{a}_i=\underset{m=1}{\overset{n}{\Σ}}(G_{im}{e}_m-B_{im}{f}_m)\\ b_i=\underset{m=1}{\overset{n}{\Σ}}(G_{im}{f}_m+B_{im}{e}_m)\end{array}\right.---\left(5\right)$

式中，n为电力系统的节点数；

PV节点的有功功率及电压偏差计算公式为：

$\left\{\begin{array}{c}{\mathrm{\ΔP}}_i=P_{is}-P_i=P_{is}-e_i{a}_i-f_i{b}_i\\ {{\mathrm{\ΔV}}_i}^2=V_{is}^2-(e_i^2+f_i^2)\end{array}\right.---\left(6\right)$

式中，V_is为节点i给定的电压幅值；

平衡节点不参与迭代计算，不需要计算功率偏差或电压偏差；

求各节点功率或电压偏差中绝对值最大的值，称为最大不平衡量ΔW_max；

F、判断最大不平衡量绝对值|ΔW_max|是否小于收敛精度ε；如果小于收敛精度ε，执行步骤J；否则，执行步骤G；

G、形成雅可比矩阵J

雅可比矩阵元素计算方法根据支路阻抗值和节点类型采用不同方法：

潮流计算按平启动设置电压初值，正常电力系统按此电压初值计算的节点功率P_i和Q_i与实际值相差不大；当电力系统为包含小阻抗支路的病态系统时，小阻抗支路两端电压初值可能在其等值阻抗上产生较大的电压差，会计算出很大的支路功率，进而在两端节点得到与实际值相差很大的节点注入功率P_i和Q_i以及与实际值相差很大的注入电流相量的实部和虚部a_i、b_i，以此a_i、b_i计算雅可比矩阵元素，将导致潮流计算发散；

潮流计算收敛时，式(4)中ΔP_i、ΔQ_i都趋近于0，节点注入功率计算值等于节点注入功率给定值，用节点注入功率给定值计算注入电流相量的实部和虚部,记为a_is和b_is,计算公式为：

$\left\{\begin{array}{c}{a}_{is}=\frac{a_i{P}_{is}+f_i{Q}_{is}}{e_i^2+f_i^2}\\ b_{is}=\frac{f_i{P}_{is}-e_i{Q}_{is}}{e_i^2+f_i^2}\end{array}\right.---\left(7\right)$

对于小阻抗支路PQ端点来说，a_is和b_is会比直接采用式(5)计算a_i和b_i更接近实际值，更有利于改善潮流计算的收敛性；由于PV节点的电源无功功率不是给定的，输入时输入任意值，用此任意值计算a_is和b_is，在潮流计算的后几次迭代中则与实际值差值较大，不利于提高潮流计算的收敛速度，小阻抗支路PV端点采用传统方法用a_i和b_i计算雅可比矩阵元素，不会导致潮流计算发散；因此小阻抗支路PV端点仍然采用传统方法计算雅可比矩阵元素；正常支路的端点也采用传统方法计算雅可比矩阵元素；

形成雅可比矩阵元素的具体步骤如下：

G1、计算i≠j时的雅可比矩阵元素；

当i≠j时，雅可比矩阵J的元素计算公式如下：

$\frac{\∂{\mathrm{\ΔP}}_i}{\∂e_j}=-G_{ij}{e}_i-B_{ij}{f}_i---\left(8\right)$

$\frac{\∂{\mathrm{\ΔP}}_i}{\∂f_j}=B_{ij}{e}_i-G_{ij}{f}_i---\left(9\right)$

$\frac{\∂{\mathrm{\ΔQ}}_i}{\∂e_j}=B_{ij}{e}_i-G_{ij}{f}_i---\left(10\right)$

$\frac{\∂{\mathrm{\ΔQ}}_i}{\∂f_j}=G_{ij}{e}_i+B_{ij}{f}_i---\left(11\right)$

$\frac{\∂{{\mathrm{\ΔV}}_i}^2}{\∂e_j}=0---\left(12\right)$

$\frac{\∂{{\mathrm{\ΔV}}_i}^2}{\∂f_j}=0---\left(13\right)$

G2、令i＝1；

G3、判断是否同时满足以下条件：

T_i＝1；

节点i是PQ节点；

如果不满足转步骤G4；如果满足，则按修改公式计算i＝j时的雅可比矩阵元素，然后转步骤G5；

计算i＝j时的雅可比矩阵元素的修改公式为：

$\frac{\∂{\mathrm{\ΔP}}_i}{\∂e_i}=-a_{is}-G_{ii}{e}_i-B_{ii}{f}_i---\left(14\right)$

$\frac{\∂{\mathrm{\ΔP}}_i}{\∂f_i}=-b_{is}+B_{ii}{e}_i-G_{ii}{f}_i---\left(15\right)$

$\frac{\∂{\mathrm{\ΔQ}}_i}{\∂e_i}=b_{ia}+B_{ii}{e}_i-G_{ii}{f}_i---\left(16\right)$

$\frac{\∂{\mathrm{\ΔQ}}_i}{\∂f_i}=-a_{is}+G_{ii}{e}_i+B_{ii}{f}_i---\left(17\right)$

G4、按传统计算i＝j时的雅可比矩阵元素；

传统方法计算i＝j时的雅可比矩阵元素的公式为：

$\frac{\∂{\mathrm{\ΔP}}_i}{\∂e_i}=-a_i-G_{ii}{e}_i-B_{ii}{f}_i---\left(18\right)$

$\frac{\∂{\mathrm{\ΔP}}_i}{\∂f_i}=-b_i+B_{ii}{e}_i-G_{ii}{f}_i---\left(19\right)$

$\frac{\∂{\mathrm{\ΔQ}}_i}{\∂e_i}=b_i+B_{ii}{e}_i-G_{ii}{f}_i---\left(20\right)$

$\frac{\∂{\mathrm{\ΔQ}}_i}{\∂f_i}=-a_i+G_{ij}{e}_i+B_{ii}{f}_i---\left(21\right)$

$\frac{\∂{{\mathrm{\ΔV}}_i}^2}{\∂e_i}=-2e_i---\left(22\right)$

$\frac{\∂{{\mathrm{\ΔV}}_i}^2}{\∂f_i}=-2f_i---\left(23\right)$

PQ节点按式(18)-(21)计算雅可比矩阵元素；PV节点按式(18)、(19)、(22)、(23)计算雅可比矩阵元素；平衡节点不计算雅可比矩阵元素；

G5、令i＝i+1；

G6、判断i是否大于节点数n，如果i不大于n转步骤G3；否则转步骤H；

H、解修正方程及修正电压实部e、虚部f

潮流计算的基本方程(4)和(6)是非线性方程组，采用逐次线性化方法迭代求解；线性化得到的方程称为修正方程，用来求电压实部和虚部的修正量；

修正方程为：

$\left[\begin{array}{c}\mathrm{\Δ}P\\ \mathrm{\Δ}Q\\ {\mathrm{\ΔV}}^2\end{array}\right]=J\left[\begin{array}{c}\mathrm{\Δ}e\\ \mathrm{\Δ}f\end{array}\right]=\left[\begin{array}{cc}\frac{\∂\mathrm{\Δ}P}{\∂e^T}& \frac{\∂\mathrm{\Δ}P}{\∂f^T}\\ \frac{\∂\mathrm{\Δ}Q}{\∂e^T}& \frac{\∂\mathrm{\Δ}Q}{\∂f^T}\\ \frac{\∂{\mathrm{\ΔV}}^2}{\∂e^T}& \frac{\∂{\mathrm{\ΔV}}^2}{\∂f^T}\end{array}\right]\left[\begin{array}{c}\mathrm{\Δ}e\\ \mathrm{\Δ}f\end{array}\right]---\left(24\right)$

式中，J为雅可比矩阵；ΔP和ΔQ分别为有功功率和无功功率偏差列向量；ΔV²为电压幅值偏差列向量；Δe和Δf分别为电压相量的实部和虚部修正量列向量；为有功功率偏差函数列向量对电压相量实部列向量转置的偏导矩阵，上标T为转置符号；

电压修正公式为：

$\left\{\begin{array}{c}{e}_i^{(t+1)}=e_i^{\left(t\right)}-{\mathrm{\Δe}}_i^{\left(t\right)}\\ f_i^{(t+1)}=f_i^{\left(t\right)}-{\mathrm{\Δf}}_i^{\left(t\right)}\end{array}\right.---\left(25\right)$

式中，上标t表示第t次迭代；

I、令t＝t+1，返回步骤E进行下一次迭代；

J、输出节点及支路数据。