第3章 P—Q分解法潮流计算基本原理
3.1 牛顿—拉夫逊法的基本原理
设有单变量非线性方程
f(x)=0 （3-1）
求解此方程时，先给出解的近似值x(0)，它与真解的误差为Δx(0)，则将满足方程(3-1)，即
f(x(0)+Δx(0))=0 （3-2）
将式(3-2)左边的函数在x(0)附近展成泰勒级数，于是便得
（3-3）
式中f’(x(0)),…，f(n)(x(0))分别为函数f(x)在x(0)处的一阶导数，…，n阶导数。
如果差值Δx(0)很小，式(4-3)右端的二次及以上阶次的各项均可略去。于是，式(3-3)可简化为
（3-4）
这是关于修正量Δx(0)的线性方程式，亦称为修正方程式。解此方程式可得修正量 （3-5）
用所求得的Δx(0)去修正近似解，便得
由于式(3-4)是略去了高次项的简化式，因此所解出的修正量Δx(0)也只是近似值。修正后的近似解x(1)与真解仍然有误差。但是，这样的迭代计算可以反复进行下去，迭代计算的通式是
（3-6）
迭代过程的收敛判据为
（3-7）
或 （3-8）
式中的ε1，ε2是预先给定的小正数。
图（3.1）牛顿法的几何解释
这种解法的几何意义可以从图(3.1)中得到证明。由此可见，牛顿—拉夫逊法实质上就是切线法，是一种逐步线性化的方法。设有n个联立的非线性代数方程：
（3-9）
应用牛顿法求解多变量非线性方程组(3-9)时，假定已给出各变量的初值x1(0)，x2(0)，…，xn(0)，令Δx1(0)，Δx2(0)，…，Δxn(0)分别为各变量的修正量，使其满足方程，即
（3-10）
将上式中的n个多元函数在初始值附近分别展成泰勒级数，并略去含有Δx1(0)，Δx2(0)，…，Δxn(0)的二次及以上阶次的各项，使得
方程式可以写成矩阵形式
（3-12）
方程式(3-12)是对于修正量Δx1(0)，Δx2(0)，…，Δxn(0)的线性方程组，称为牛顿法的修正方程式。利用高斯消去法或三角分解法可以解出修正量Δx1(0)，Δx2(0)，…，Δxn(0)。然后对初始近似解进行修正
xi(1)=xi(0)+Δxi(0) (i=1,2,…,n) （3-13）
如此反复迭代，在进行第次迭代时，从求解修正方程式
（3-14）
得到修正量Δx1(k)，Δx2(k)，…，Δxn(k)，并对各变量进行修正
xi(k+1)=xi(k)+Δxi(k) (i=1,2,…,n) （3-15）
式(4-14)和(4-15)也可以缩写为
F(x(k))=-J(k)Δx(k) （3-16）
和
x(k+1)=x(k)+Δx(k) （3-17）
式中，x和Δx分别是有n个变量和修正量组成的n维列向量；F(x)是由n个多元函数组成的n维列向量；J是n×n阶方阵，称为雅可比矩阵，它的第Jij个元素
迭代过程一直进行到满足收敛判据
（3-18）
或 （3-19）
为止，ε1和ε2为预先给定的小正数。
将牛顿—拉夫逊法进行潮流计算，要求将潮流方程写成形如（3-9）的形式。由于节点电压可以采用不同的坐标系来表示，牛顿—拉夫逊法潮流计算也将相应的采用不同的计算公式。
3.2 节点电压用极坐标牛顿—拉夫逊法潮流计算
采用极坐标时，节点电压表示为
节点的功率方程(4-20)将写成
（3-20）
式中，δij=δi-δj是两节点电压的相角差。
方程式(4-33)把节点功率表示为节点电压的幅值和相角的函数。在有n个节点的系统中，假定第1～m号节点为节点，第m+1～n-1号节点为节点，第n号节点为平衡节点。Vn和δn是给定的，PV节点的电压幅值Vm+1～Vn-1也是给定的。因此，只剩下n-1个节点的电压相角δ1，δ2，…，δn-1和m个节点的电压幅值V1，V2，…，Vm是未知量。
实际上，对于每一个PQ节点或每一个PV节点都可以列写一个有功功率不平衡量方程式
(i=1,2,…,n) （3-21）
而对于每一个PQ节点还可以在列写一个无功功率不平衡量方程式
(i=1,2,…,m) （3-22）
式（3-21）和（3-22）一共包含了n-1+m个方程式，正好同未知数的数目相同，而比直角坐标形式的方程少了n-1-m个。
对于方程式（3-21）和（3-22）可以写出修正方程式如下
（3-23）
式中
（3-24）
H是(n-1)×(n-1)阶方阵，其元素为 ；N是(n-1)×m阶矩阵，其元素为 ；K是m×(n-1)阶矩阵，其元素为 ；L是m×m阶方阵，其元素为 。
在这里把节点不平衡功率对节点电压幅值的偏导数都乘以该节点电压，相应地把节点电压的修正量都除以该节点的电压幅值，这样，雅可比矩阵的表达式就具有比较整齐的形式。
对式（3-21）和式（3-22）求偏导数，可以得到雅可比矩阵元素的表达式如下：
（当i≠j时） （3-25）
（当i=j时） （3-26）
其计算步骤和程序框图在这里不给出。
3.3 P—Q分解法的原理
采用极坐标形式表示节点电压，能够根据电力系统实际运行状态的物理特点，对牛顿潮流计算的数学模型进行合理的简化。
在交流高压电网中，输电线路的电抗要比电阻大得多，系统中母线有功功率的变化则主要受母线电压幅值变化的影响。在修正方程式的系数矩阵中，偏导数 和 数值是相当小的。作为简化的第一步，可以将方程式（3-24）中的子块N和K略去不计，即认为它们的元素为零。这样，n-1+m阶的方程式（3-24）便分解为一个n-1阶和一个m阶的方程
△P = - H△δ （3-27）
△Q= - LVD-1△V （3-28）
这一简化大大地节省了机器内存和解题时间。方程式（3-27）和（3-28）表明，节点的有功功率不平衡量只用于修正电压的相位，节点的无功功率不平衡量只用于修正电压的幅值。这两组方程轮流迭代，这就是所谓的有功—无功功率分解法。
但是矩阵H和L的元素都是节点电压幅值和相角差的函数，其数值在迭代过程中是不断变化的。因此，最关键的一步简化就在于，把系数矩阵H和L简化为常数矩阵。它的根据是什么呢？在一般情况下，线路两端电压的相角差是不大的（不超过１０o～２０o），因此可以认为
cosδij≈1, Gijsinδij<<Bij
此外，与系统各节点无功功率相适应的导纳ＢLDi必远小于该节点自导纳的虚部，即
ＢＬＤi＝ 〈〈Ｂii和Ｑi〈〈 Ｂii
考虑到以上的关系，矩阵Ｈ和Ｌ的元素的表达式便简化成
Ｈij＝ＶiＶjＢij（i，j=1，2， …，n-1） （3-29）
Ｌij＝ＶiＶjＢij（i，j=1，2， …，m） （3-30）
而系数矩阵Ｈ和Ｌ则可以分别写成
V1B11V1 V1B12V2 V1B1，n-1V1
H = V2B21V1 V2B22V2 V2B2，n-1V2
Vn-1Bn-1,1V1 Vn-1Bn-1,2V2 Vn-1Bn-1，n-1Vn-1
V1 B11 B12 B1，n-1
= V2 B21 B22 B2，n-1
Vn-1 Bn-1,1 Bn-1,2 Bn-1，n-1
V1
V2 =VD1B/VD1 （3-31）
Vn-1
V1B11V1 V1B12V2 V1B1，n-1V1
L = V2B21V1 V2B22V2 V2B2，n-1V2
VmBm,1V1 VmBm,2V2 VmBm，mVm
V1 B11 B12 B1，m
= V2 B21 B22 B2，m
Vm Bm,1 Bm,2 Bm，m
V1
V2 =VD2B//VD2 （3-32）
Vm
将（3-32）和（3-32）分别代入（3-27）和（3-28），便得到
△P= -VD1B/VD1△δ
△Q= -VD2B//△V
用 和 分别左乘以上两式便得
△P= - B/VD1△δ （3-33）
△Q= - B//△V （3-34）
这就是简化了的修正方程式，它们也可展开写成
B11 B12 B1，n-1 V1△δ1
= B21 B22 B2，n-1 V2△δ2
（3-35）
Bn-1,1 Bn-1,2 Bn-1，n-1 Vn-1△δn-1
B11 B12 B1，m V1△δ1
= B21 B22 B2，m V2△δ2
（3-36）
Bm,1 Bm,2 Bm，m Vm△δm
在这两个修正方程式中，系数矩阵都由节点导纳矩阵的虚部构成，只是阶次不同，矩阵B/为n-1阶，不含平衡节点对应的行和列，矩阵B//为m阶，不含平衡节点和PV节点对应的行和列。由于修正方程式的系数矩阵为常数矩阵，只要作一次三角分解，即可反复使用，结合使用稀疏技巧，还可以进一步的节省机器内存和计算时间。
利用公式（3-21）和（3-22）计算节点功率的不平衡量，用修正方程（3-35）和（3-36）解出修正量，并按下述条件
max{ | △P i(k)| }〈εp，max{ | △Q i(k)| }〈εQ
检验收敛，这就是分解法的主要计算内容（其程序框图将在第五章给出）。
需要说明，分解法所作的种种简化只涉及到解题过程，而收敛条件的校验仍然是以精度的模型为依据，所以计算结果的精度是不受影响的。单要注意，在各种简化条件中，关键是输电线路的r/x比值的大小。110kV及以上电压等级的架空线r/x比值较小，一般都符合PQ分解法的简化条件。在35kV及以下电压等级的电力网中，线路的r/x比值比较大，在迭代计算中可能出现不收敛的情况。
顺便指出，P—Q分解法在实际应用中还有一些改进。最常用的是，在形成P—δ迭代用的矩阵B/时，将一些对有功功率和电压相位影响较小的因素略去不计，即在计算B/的对角线元素时，忽略输电线路和变压器Π型等值电路中的对地电纳支路。实验表明，这样的处理能加快P—δ迭代的收敛过程。
3.4 P—Q分解法的特点
P—Q分解法与牛顿法潮流程序的主要差别表现在它们的修正方程式上。P—Q分解法通过对电力系统具体特点的分析，对牛顿法修正方程式的雅可比矩阵进行了有效的简化和改进，最后得到（3-35）和（3-36）所示的修正方程式，这两组方程式和牛顿法修正方程式（3-23）相比，有以下三个特点：
（1） 式（3-35）、（3-36）用二个n阶线形方程组代替了一个2n阶线形方程组；
（2） 式（3-35）、（3-36）中系数矩阵的所有元素在迭代过程中维持常数；
（3） 式（3-35）、（3-36）中系数矩阵是对称矩阵。
第一个特点在提高计算速度和减少内存方面的作用是很明显的，在这里不再述及。
第二个特点使我们得到以下好处。首先，因为修正方程式的系数矩阵就是导纳矩阵的虚部，因此在迭代过程中不必像牛顿法那样进行形成雅可比矩阵的计算，这样不仅减少了运算量，而且大大简化了程序。其次，由于系数矩阵在迭代过程中维持不变，因此在求解修正方程式时，不必每次都对系数矩阵进行消去运算，只需要在进入迭代过程以前，将系数矩阵用三角分解形成因子表，然后反复利用因子表对不同的常数项△P/V或△Q/V进行消去运算和回代运算，就可以迅速求得修正量，从而显著提高了迭代速度。
第三个特点可以使我们减少形成因子表时的运算量，而且由于对称矩阵三角分解后，其上三角矩阵和下三角矩阵有非常简单的关系，所以在计算机中可以只存储上三角矩阵或下三角矩阵，从而也进一步节约了内存。
由于P—Q分解法大大提高了潮流计算的速度，所以不仅可以用于离线计算，而且也可以用于电力系统在线静态安全监视。
P—Q分解法所采用的一系列简化假定只影响修正方程式的结构，也就是说只影响了迭代过程，但不影响起最终结果。因为P—Q分解法和牛顿法都采用同样的数学模型，最后计算功率误差和判断收敛条件都严格按照精确公式进行的，所以P—Q分解法和牛顿法一样都可以达到很高的精度。
第4章 P—Q分解中用到的技术和方法
4.1 网络节点编号的优化
目前电力系统计算程序中，在解电力网络节点方程I=YV时，大多采用直接解法。为了对网络方程反复求解，往往首先采用顺序消去法对导纳矩阵进行三角分解。然后将分解以后的三角矩阵及对角矩阵存储在计算机内，就可以用对不同的右端常数项进行前代和回代运算，从而得到网络方程的解。
如前所述，导纳矩阵是零元素很多的稀疏矩阵，分解后得到的三角矩阵一般也是稀疏矩阵。为了节约内存及避免在计算机内对非零元素的运算，在计算机内存中将只储存三角矩阵中的非零元素。因此，三角矩阵中非零元素的个数对内存的需要量及程序的计算速度将有直接的影响。
一般地讲，导纳矩阵非零元素的分布和分解后的三角矩阵是不同的。
例4.1：一个导纳矩阵其下三角部分数据为：
行号
列号
数据
2
1
-0.624024+j3.900156
3
1
-0.754717+j2.614509
3
2
-0.829876+j3.112033
4
2
0.00000+j63.49206
5
3
0.00000+j31.74603
求其经三角分解后的下三角元素为：
行号
列号
数据
2
1
-0.612227+j0.034973
3
1
-0.425687-j0.026671
3
2
-0.073971-j0.017193
4
2
-0.137743-j0.018393
4
3
-0.137743-j0.027718
5
3
-0.924654+j0.027559
5
4
-1.189287+j0.048151
在例4.1中导纳矩阵下三角部分共有五个非零元素，但把这个导纳矩阵分解以后的下三角矩阵却有七个非零元素，增加了非零元素l43和l54。这种元素称为消去过程或分解过程中的注入元素。
消去过程中产生注入元素的原因，可以直观地用电路的星网变换来解释。一般地讲，当消去节点k时，以k为中心的星形网络将变为一个以与节点K直接联系的节点为顶点的网型网络，如图4.1所示。
图（4.1）消去节点时的星网变换
如果与节点k相连的节点数为Jk，则网形网络的支路数为 Jk（Jk-1）。假设在节点k消去前，其周围Jk个节点间原有的支路数为Dk，则在消去节点k以后所增加的新支路（即注入元素的个数）为：
△b k= Jk（Jk-1）-Dk （4-1）
非对角元素lij的计算公式为：
lij =（Yij - ）/dij （4-2）
当Yij为零元素时，在这个位置上出现非零元素的条件为：
0 （4-3）
因为导纳矩阵三角分解过程中，一般对角元素dkk 0，所以上述条件相当于在下三角矩阵的i和j行中，行号k<j的任何一对 lik、ljk同时为非零元素时，就产生注入元素lij。
注入元素的多少与消去的顺序或节点编号有关，在图4-2中表示了一个简单电力网络的四种不同的节点编号方案和将导纳矩阵分解以后三角矩阵中出现注入元素的情况。显然，不同的节点编号方案所得到的注入元素的数目也不同。
节点编号图形
导纳矩阵
下三角矩阵
注入
1
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
*
。
。
*
*
。
。
*
*
*
。
6
2
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
*
。
。
*
*
。
3
3
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
*
。
1
4
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
。
0
。 ——非零元素
* ———非零注入元素
图（4-2） 节点编号对注入元素的影响
所谓节点编号的优化，就是要寻求一种使注入元素数目最少的节点编号方式。为此，可以比较四种不同的节点编号方案在三角分解中出现注入元素的情况，从中选取注入元素最少的节点编号方案。但这样做需要分析非常多的方案，例如对仅有5个节点的电力网络来说，其编号的可能方案就有5！=120个。一般，对有n个节点的电力网络来说节点编号的可能方案就有n！个，工作量非常大。因此，在实际计算工作中往往采用一些简化的方法，求出一个相对的节点编号优化方案，但不一定追求“最优”方案。
目前，节点编号优化的方法很多，大致可以分以下三类：
（1）静态按最少出线支路数编号
这种方法又称为静态优化法。在编号以前，首先统计电力网络各节点的出线支路数，然后，按出线支路数少的节点顺序编号，当有n个节点的出线支路数相同时，则可以按任意次序对这n个节点进行编号。
这种编号方法的依据是，在导纳矩阵中，出线支路数最少的节点所对应的行中非零元素也最少，因此在消去过程中产生注入元素的可能性也比较小；
（2）动态地按最少出线支路数编号
在上述方法中，各节点的出线支路数是按原始网络统计出来的，在编号过程中认为固定不变。事实上，在节点消去过程中，每消去一个节点以后，与该节点相连的各节点的出线支路数将发生变化（增加，减少或保持不变）。因此，如果在每消去一个节点后，立即修正尚末编号节点的出线支路数，然后选其中出线支路数最少的一个节点进行编号，就可以预期得到更好的效果。所谓动态地按最少出线支路数编号的方法和前一种方法不同之处就是在按最少出线最少原则编号时考虑了消去过程中各节点出线数目的变动情况。这种方法也称为半动态优化法。
（3）动态地按增加出线数最少编号
这种方法又称为动态优化法。用以上两种方法编号，只能使消去过程中出现新支路的可能性减少，但并不一定保证在消去这些节点时出现的新支路最少，比较严格的方法应该是按消去节点后增加出线数最少的原则编号。具体编号方案如下：
首先，根据星网变换的原理，按式（4-1）分别统计消去网络各节点时增加的出线数，选其中增加出线数最少的被消节点编为第1号节点。
确定了第1号节点以后，即从网络中消去该节点，相应地修改其余节点的出线数目。
然后，对网络中其余的节点重复上述过程，顺序编出第2号，第3号…一直到编完为止。
很明显，这种编号方法的工作量比以上两种方法大得多。所以这种编号方法一般不与采用。
考虑到编程的简单性、工作量大小以及是否节省内存和提高计算速度，在本程序中应用（1）静态按最少出线支路数编号对网络的各节点进行编号。
4.2 因子表
在实际计算中，常常遇到这种情况：对于方程组AX=F需要多次求解，每次仅改变其常数项F，而系数矩阵A通常是不变的。这时，为了提高计算速度，可以利用因子表对线形方程组求解。
因子表可以理解为高斯消去法解线性方程过程中对常数项F全部运算的一种记录表格。高斯消去法分为消去过程和回代过程。回代过程的运算对系数矩阵进行消去后得到的上三角矩阵元素确定。见式（4-4）、式（4-5）。
1
B=
1 （4-4）
1
（4-5）
为了对常数项进行消去运算（又称回代运算），还必须记录消去过程运算所需要的运算因子。消去过程又分为消去运算和规格化运算。由式（4-6）可知，消去过程对常数项F中的第i个元素fi的运算包括：
fi(k)= fi(k-1) - aik(k-1) fk(k) (k=1,2….i-1) （4-6）
fi(i)= fi(i-1)/ aik(i-1) （4-7）
将上式中的运算因子 及 琢行放在下三角部分和式（4-4）的上三角矩阵元素合在一起，就得到因子表：
（4-8）
其中下三角及对角线元素可用来对常数项F进行前代运算，上三角用来进行回代运算。因子表也可写成如下的形式：
D11 U12 U13 U14 U1n
L21 U22 U23 U24 U2n
L31 L32 U33 U34 U3n （4-9）
Ln1 Ln2 Ln3 Ln4 U1n
其中Dii=
Uij = （i<j）
Lij= （j<i）
不难看出，因子表式（4-9）中下三角部分的元素就是系数矩阵在消去过程中曾出现的元素，因此只要把它们保留在原来的位置，并把对角线取倒数就可以得到因子表。
