4.第四章 潮流计算中的特殊问题
第四章 潮流计算中的特殊问题
第一节 负荷的静态特性
负荷的功率是系统频率和电压的函数。在潮流计算中可以认为频率变化不大。但由于发电机或输电设备的开断会引起电压较大的变化，在潮流计算中计及负荷的静态电压特性是合理的。
负荷的电压静态特性就是负荷的有功和无功功率与电压大小的关系，一般表达如下：
????
?????????????+???? ??+???? ??=??
?
?????+???? ??+???? ??=Qi is i Qi is i
Qi Di Di Pi is i Pi is i Pi Di Di c V V b V V a Q Q c V V b V V a P P 2
)0(2)0( （4-1） 式中系数满足
1
1=++=++Qi Qi Qi Pi Pi Pi c b a c b a
)0(Di P 、)
0(Di Q 是在设定电压is V 下的负荷值。组成负荷的三部分被分别看做恒定阻
抗部分、恒定电流部分和恒定功率部分，所以（4-1）称为负荷的ZIP 模型。
当0=Pi a 、0=Qi a 时，忽略电压的二次项。
'
潮流计算中计及负荷的静态电压特性的方法：
1、节点功率的不平衡量计算:
????
?????-????????+???? ??+???? ??-=--=?-???
?????+???? ??+???? ??-=--=?)
,(),()
,(),(2
)0(2)0(θθθθV Q c V V b V V a Q Q V Q Q Q Q V P c V V b V V a P P V P P P P i Qi is i Qi is i
Qi Di Gi i Di Gi i i Pi is i Pi is i
Pi Di Gi i Di Gi i (4-2) 2、牛顿法雅可比矩阵子矩阵N 和L 的对角线元素要增加
i i V P ???和i
i
V Q ???
3、P-Q 分解法，Q-V 迭代的系数矩阵B ''的对角线元素也应增加i
i
V Q ???，这样B ''不再是常数了。为了节省计算量，
i
i
V Q ???也可取为常数，如忽略二次项取0=Qi a ，或不改变B ''，但功率不平衡量要按（4-2）计算。
负荷电压静态特性模型的指数形式
???
????
?
??? ??=?
??
? ??=β
α
is i Di Di is i Di
Di V V Q Q V V P P )0()0( （4-3） 8.1~5.0=α、6~5.1=β
在动态潮流计算中，不能不考虑频率的变化。考虑频率变化时式（4-1）、（4-3）变为。
????
????????? ??-+????????+???? ??+???? ??=???? ??-+???
?????+???? ??+???? ??=002
)0(002)0(11f f f k c V V b V V a Q Q f f f k c V V b V V a P P Qi Qi is i Qi is i
Qi Di Di Pi Pi is i Pi is i
Pi Di Di ,
???
???????? ??-+???? ??=???
? ??-+???? ??=00)0(00)0(11f f f k V V Q Q f f f k V V P P Qi is i Di Di Pi is i
Di Di β
α
当考虑频率变化时，频率也是待求的未知量，应出现在潮流方程中。 模型中系数的选取属于负荷建模的问题，仍未得到很好的解决。
第二节 节点类型的相互转换
一、PV 节点转换为PQ 节点
当在迭代过程中出现PV 节点无功功率越限时，可以再迭代几次，如果无功
仍越限，说明PV 节点电压设置不合理，应进行调整：
如果无功功率越下限，检查是否电压设置过低如是可适当提高电压设定值，或转换为PQ 节点，无功定值置下限值。
如果无功功率越上限，说明节点无功功率不能支持设定的电压，可适当调低电压设定值，或转换为PQ 节点，无功定值取上限值。
PV 节点转换为PQ 节点的处理方法： —
1、直角坐标方式的节点不平衡量由2i V ?变为i Q ?；
2、牛顿法极坐标方式的修正方程加1个Q ?方程；
3、P-Q 分解法，θ-P 迭代不变，V Q -迭代的系数矩阵有两种处理方法： （1）B ''增加一行一列，如增加到最后：
??
?
???''=''ii T
i
B B B B B i ~ （4-4） 新的矩阵的因子表可由右下角加边的因子表修正法求出。 （2）B ''的对角元加大数
在形成B ''时包含PV 节点对应的导纳，但PV 节点的对角元加一个很大的数。这样在正常Q-V 迭代时，PV 节点的电压修正零接近于0，不会影响其他节点的电压修正量。当PV 转换为PQ 节点时，将加的大数去掉。
ΔB B B -''=''~
（4-5） 采用因子表秩1修正法得到新的因子表
?
二、PQ 节点转换为PV 节点
当在迭代过程中出现PQ 节点电压越限时，可以再迭代几次，如果电压仍越限，说明PQ 节点无功设置不合理，应进行调整：
如果电压越下限，说明无功设置较低，可适当提高无功设定值，或转换为
PV 节点，电压定值取下限值。
如果电压越上限，说明节点无功设定偏高，可适当调低无功设定值，或转换为PV 节点，电压定值取上限值。
PQ 节点转换为PV 节点的处理方法：
1、直角坐标方式的节点不平衡量由i Q ?变为2i V ?；
2、牛顿法极坐标方式的修正方程减1个Q ?方程；
3、P-Q 分解法，θ-P 迭代不变，V Q -迭代的系数矩阵有两种处理方法： （1）在B ''中划去将要转换为PV 节点的节点所在的行和列，重新形成因子表。
（2）在B ''中将要转换为PV 节点的节点对应的对角元加一个很大的数，用因子表秩1修正法得到新的因子表
,
三、因子表修正方法
1、因子表秩1修正法
设系数矩阵A 已因子化为如下的形式
LDU A = （4-6）
由于某种原因，A 变化为：
A A N M A A ?+=+=T a ~
（4-7） 其中M 和N 为1?n 的列矢量，a 为标量。新矩阵A ~
的因子表为： U D L A ~
~~~= （4-8） 将（4-8）、（4-6）代入（4-7）有：
T M aN LDU U D L +=~
~~ （4-9）
&
为了求出U D L ~
~~中的各元素，将U D L ~
~~和LDU 各矩阵的第一行和第一列单独列
出，并写成分块矩阵的形式：
????
??=11L l L 1
??????=1D D 1d ???
??
?=1U u U 11 (4-10) 和
??????=11
L l L ~~
1
~ ?????
?=1D D ~~
~1d ??
?
???=1U u U ~~1~1 (4-11) 及
??????=11M M m ???
???=11N N n （4-12）
将（4-10）代入（4-6），A 矩阵可写为：
?
?
?
???+=111111111
11U D L u l l u A d d d d
（4-13） 将（4-11）代入（4-8），A ~
矩阵可写为：
??
????+=1111111
11
11
~~~~~~~
~~
~~
~U D L u l l u A d d d d （4-14） #
将（4-12）代入T
MaN A =?
??
??
??=?T T a an a m an m 11111
11
1N M M N A （4-15） 将（4-13）、（4-14）、（4-15）代入（4-7）有
??
??
??+??????+=??????+T T
a an a m an m d d d d d d d d 111
11
11
11111111
11111111111
11
11
~~~~~~~
~~
~~
N M M N U D L u l l u U D L u l l u （4-16） 根据等号两端矩阵对应元素相等，可得：
（1） 1111~
an m d d += （4-17）
（2） T a m d d 111111~~
N u u +=将（4-17）变为1
111~
an m d d -=代入，有 T a m d 111111~~~N u u -+= （4-18）
其中
1111~
u N N n T
T -= （4-19）
*
（3） 111111~~an d d M l l +=将（4-17）变为1111~
an m d d -=代入，有
111111~
~
~
-+=d an M l l （4-20） 其中
1111~
m l M M -= （4-21） 由上（1）、（2）、（3）可计算出新矩阵因子表上三角矩阵第一行元素、下三角矩阵第一列元素和对角线矩阵第一个元素。
（4） T a d d 1
1111111111111~
~~~~
~N M U D L u l U D L u l ++=+ 重写为
111111111111111111~~~~~~A U D L N M U D L u l u l U D L ?+=++-=T
a d d （4-22）
其中T a d d 1
11111111~~
~N M u l u l A +-=?将（4-18）、（4-20）代入得 T T
T T
T T
T T
T T T
T T T
T T
T T T T T T
T T
a a m d an a a m d an n a m a m d an an m a m a m d an a an a m an m a m d an a an m an an m a m an m a m d an a an m n a m an m a a m d an an a m an m a a m d d d an d d an a m d d d an m d a a m d d d an an m d a d d 1
1
1
111111
111111111111
111111*********
111111111111111111
11111111111111111111111111
1111111111111111111111
111111111111111111
1111111111111111111111111111111
111111111111111111
11111111~~~~)~(~~~~)()(~~~)()(~~~~~~~~~)()(~~~~~~~~~~~~~~~~~~)~~(~)~~()~(~~~N
M N M N M u N l M N M u l M N l M N M N M u M N l u l N M N M u l u M u l N l u l N M N M u l M u N l u l N M N M u M N l u l N M N M u M N l u l u l u l N M N u M l u l N M u l u l A =-=---=----=-+--=-++-+--=-+-----=+----=+-----=+++--=+-=?-------------
（4-23）
>
其中
)~(~1111a m d an a a --= （4-24）
因此，（4-22）可写为如下的形式
T a 1
1111111~
~~
~
~~N M U D L U D L += （4-25） （4-25）与（4-9）有同样的形式，可用（1）、（2）、（3）的方法分别求出1
11~
~~U D L
矩阵的第一行第一列元素。
因子表的秩1修正过程就是递归使用式（4-17）、（4-18）和（4-20）逐次求出新矩阵的因子表的过程。程序流程为：
n
n n n i i T i i i i i T
i i i i i i i i i i i i i i i an m d d end loop
a
m d a-an a a m d n d an m l an
m d d ,...,n-,loop i i
i i +=????
????
???????=+=-=+=-=+==---111121N u u u N N M l l M M （4-26） 由于M 和N 是非零元素非常少的稀疏矢量，使用稀疏矢量的排零运算后，上述程序流程的计算量很小。
2、原矩阵右下角加边的因子表修正法 已知原矩阵A 的因子表
LDU A = （4-27）
}
在原矩阵A 右下角增加m 行m 列后形成的加边矩阵A ~
及因子表表示为：
??
?
????????
???????=???
???=a M a N u u U d D
l l L
a N M A A ~
（4-28） 其中M 是m n ?阶矩阵，N 是n m ?阶矩阵。由因子分解方法可知，原矩阵因子表的元素不变化。右边展开
??
?
???+=??????ααdu l Du l DU l LDu LDU a N M A M N N M
（4-29） 对照、比较得
M L D u 1
1--=M （4-30） 1
1--=D NU l N
（4-31）
M N Du l a du l -=αα （4-32）
对矩阵M N Du l a -进行三角分解可求出a l 、d 和a u 矩阵各元素。当1=m 时，
1=a l 、1=a u ，（4-32）变为
M N Du l a d -= （4-33）
…
其结果为一标量。
回忆因子表法线性方程组的求解过程的前代、规格化运算 b LDUx =
b L y b Ly 1-=?=（前代）
b L D y D z y Dz 111---==?=（规格化）
可见（4-30）就是用L 和D 对M 中的各列进行的前代、规格化运算。将（4-31）
两边取转置T
T T T T N N U D D NU l 1111)()()(----==，所以（4-31）就是用U 和D 对N 中
的各行进行的前代、规格化运算。不别进行矩阵的求逆运算。
第三节 多θV 节点的潮流计算
有时我们只关心一个大电网中的一部分网络的运行状态，为了简化计算将电网分解为内部网、边界和外部网三部分。将外部网的所有节点用高斯消去法消掉，用边界节点间的等值支路来等值替代。如图所示。
要对等值后的网络进行计算，需要知道边界节点处来自外部等值网络的注入
功率。内部网和边界节点的运行参数可由SCADA系统得到，进而用状态估计方法可求出边界节点处的电压后，这样等值网络就是一个具有多θV节点的网络。
设系统节点数为N，其中有S个节点θV给定，R个节点PV给定，其余N-S-R 个节点为PQ给定。这样可写出N-S个有功约束方程和N-S-R个无功约束方程，待求的为N-S-R个节点的电压相角θ和N-S个节点的电压幅值V。待求变量数量与方程数相等，方程可解。牛顿法或P-Q分解法均可使用。
上述多θV的潮流方程收敛后，可求出边界节点处由外部网注入的等值功率。将此注入功率作为已知量，可进一步研究内部网在各种运行方式下的潮流分布。即认为在以后的分析中外部网提供的注入功率不变。
第四节潮流方程解的存在性、多值性以及病态潮流解法
一、潮流方程解的存在性、多值性
(
（1）、潮流方程是一组非线性方程组，理论上是多解的。如不同的给定值、不同的初始条件得到的解可能不同；
（2）、潮流计算得到的解可能是有实际意义的，即与实际运行状态相符或实际能运行的解；也可能是无实际意义的，如实际无法实现或无法接受的解；
（3）、潮流方程无解，或无实数解，不能收敛；
（4）、潮流方程有解，但算法不完善，不能收敛，如修正方程病态；
（5）、采用“平值启动”或状态估计的参数作为初始条件得到的解一般是有意义的。
（6）、潮流方程解的存在性、多值性是仍未解决的困难课题。
二、病态潮流及其解法
主要是指修正方程的病态。
1、病态方程
一个线性方程组 AX b =，若右端向量b 或系数矩阵A 的微小变化就会引起方程组的解发生很大的变化，则称AX b =为病态方程组。方程组的系数矩阵A 的条件数()1Cond A A A -=刻画了方程组的性态，若()1Cond A ≥，则称AX b =为“病态”方程组；若()Cond A 相对较小，则称AX b =为“良态”方程组。良态方程组用GAUSS 消去法和JACOBI 等简单的迭代法就可以得到比较好的计算解，而对于病态方程组，一般的直接法和迭代法会有较大的误差，甚至严重失真。所以，在解方程组时，有必要先对方程组的性态进行研究，采用相应的算法，才能得到比较精确的计算解。利用方程组的条件数来判断就是一个很好的办法。
}
下面的一些直观的现象可作为判别病态矩阵的参考：
（1）在主元消去法的过程中出现小主元，则A 有可能是病态矩阵，但病态矩阵未必一定有这种小主元；
（2）若解方程组时出现很大的解，则A 有可能是病态矩阵，但病态矩阵也可能有一个小解；
（3）从矩阵本身来看，若元素间数量级相差很大且无一定规律；或矩阵的某些行（列）近似线性相关，即矩阵的行列式接近于0，这样的矩阵就有可能是
病态的。
当然，这些现象只能帮助我们做初步的判断，并且很多病态矩阵也不一定会出现这些现象。最可靠的判别方法是求出矩阵的条件数。
2、病态潮流方程的求解 （1）最优乘子法
潮流方程解不收敛可能是由于修真方程病态使修正量偏大，造成解的振荡，为此在第k 次修正时采用如下的公式。
(k)(k))(k μΔx x x -=+1 （4-34）
μ是标量乘子，应满足
'
2
)()()(m in k k f x x ?-μμ （4-35）
上式表明μ的值应使第k 次修正后潮流方程0f(x)=的适配量的平方和为最小值。所以叫最优乘子法。式（4-35）是最小二乘目标函数，展开为
∑=?-=?-n
i k k i
k k f f 1
2)()
(2
)
()
())((min
)
(min x x
x
x
μμμ
μ
（4-36）
（4-35）取最小值的条件是 0)
(2
)()(=?-du
f d k k x x μ
（4-37）
为了能求出满足（4-37）的μ，做如下处理： （1）将)()()(k k f x x ?-μ用泰勒展开的前两项表示
)()()()()()()()(k k k k k f f x x f x x x ?'-=?-μμ （4-38）
（2）（4-36）中最小二乘函数变为
∑∑∑===?+?-=?-=?-n i k k i k k i k i k i n
i k k i k i n
i k k i
f f f x
x
f 1
2)()('2)()(')()(21
2
)()(')(1
2
)
()
()
))(()()(2)(())()(())((x x f x x f x x x x f x μμμμ （4-39）
@
（3）对（4-39）式求对μ的导数
∑∑∑===?+?-=?+?-=?-n
i k k i n
i k k i
k i n
i k k i k k i k i k k f f du
f d 1
2
)()('1
)
()
('
)
(12)()(')()(')(2
)()())((2))()((2)
))((2)()(2()
(x x f x x f x x x f x x f x x x μμμ（4-40）
（4）令（4-40）等于0，得
∑∑==??=
n i k k i n
i k k i k i
f 1
2
)()('1
)()(')
())(()
)()((x x f x x f x
μ
（4-41）
如潮流有解，目标函数会逐渐减小为零；如果无解，目标函数会逐渐减小到一个定值上，μ越来越小，最后为零，此时得到的解为潮流方程的最小二乘解。
（2）、非线性规划法
解潮流方程是寻求满足0f(x)=的x 的值。如果x 满足
0)(21)()(21)(1
2
===∑=m i i T f F x x f x f x ，则x 是0f(x)=的解。潮流方程的求解转化为非线
性规划问题。
)()(2
1)(min x f x f x T
x
F =
（4-42） 如果上述非线性规划问题的解使目标函数为0，则非线性规划问题的解就是潮流方程的解，潮流方程有解；如果上述非线性规划问题的解使目标函数非0，表明潮流方程无解，非线性规划问题的解是潮流方程的最小二乘解。
式（4-42）非线性规划目标函数取极小值的条件是
？
),...,2,1( 0)
(n i x F i
==??x （4-43）
即
),...,2,1( 0)()())(21()
(112
n i x f f x f x F m
j i
j j i m j j i ==??=??=??∑∑==x x x x （4-44）
上式为n 个方程组成的非线性方程组。n 为潮流方程变量数，m 为潮流方程数。可见这种方法不要求潮流方程数与变量数相等。
方程组（4-44）求解计算量太大，实用的是最优乘子法。
第五节 潮流方程中的二次型
一、二次型函数的泰勒展开
考虑只包含二次项的代数函数为如下形式
22),(Zy Bxy Ax y x f ++= （4-45） 、
该函数三阶以上的导数为零，其各一、二导数为
Z f B f A f B f Bx
Zy f By Ax f yy
yx xx
xy
y x 2222=''=''=''=''+='+='
（4-46）
该函数在0x ，0y 处泰勒展开为
)
,(),(),(),()
222(2
1
),(),(),()),(),(2),((21
),(),(),(),(000000220000002000020000000000y x f y y x f x y x f y x f y Z y x B x A y y x f x y x f y x f y y x f y x y x f x y x f y
y x f x y x f y x f y y x x f y x y x yy xy xx
y x ??+?'+?'+=?+??+?+?'+?'+=?''+??''+?''+?'+?'+=?++?+（4-47）
即二次型函数的泰勒展开二次项部分与原函数有同样的形式，只是变量换成了变量的增量。由于三阶以上的导数为零，所以（4-47）是函数的精确表达。
二、潮流方程中的二次型及求解
直角坐标形式的潮流方程为仅包含二次项的非线性代数方程组
11 )(121
0)()(121 0)()(222??
??
???-+--=+===+--=-==++-=∑∑∑∑∈∈∈∈),...,n r r,n n (i f e V ),...,n-r-,(i e B f G e f B e G f Q ),...,n ,(i e B f G f f B e G e P i i is
i j j ij j ij i i j j ij j ij i is i j j ij j ij i i j j ij j ij i is （4-48）
)(x y y =sp （4-49）
>
其中sp y 为12?n 列矢量，是已知量；)(x y 为12?n 潮流方程的函数矢量；x 为待求变量，是12?n 维矢量。
按照二次型函数的特点，将（4-49）在解的初值0x 处展开式为
)()(0x y x J x y y ?+?+=sp （4-50） J 是)(x y 在0x 的雅可比矩阵。
（4-50）是（4-49）的精确表达，没任何近似。对于潮流方程取平启动电压为初值，由（4-50）解出x ?并用x x x ?+=0修正可得原方程的解。但（4-50）也是一个非线性方程组，其形式较原方程一致，所以不是一个可取的解法。
由（4-50）可以得到一个迭代式
))()(()(01)1(k sp k x y x y y J x ?--=?-+ （4-51） 上式中J 可保持平启动电压时的值不变，所以是一种定雅可比的算法。x ?的初值可取0。式（4-51）为高斯迭代法，可以使用高斯-赛德尔迭代法改善收敛性。特点是避免了牛顿法每次迭代需要重新计算雅可比矩阵的问题。
（4-51）能得到可靠稳定解的条件是J 非病态，如果J 病态可采用最优乘子法。将（4-51）修改为
))()(()(01)1(k sp k x y x y y J x ?--=?-+μ （4-52）
μ值的选取是使（4-50）两端的偏差量的平方和为最小，先写出偏差量的表
达式
2
2
00)()()
()()(μμμμμμμc b a x y x J x y y x y x J x y y x ++=?+?--=?+?--=?sp sp f （4-53）
??
?
???-=?-=-=)()(0x y c x J b x y y a sp （4-54）
目标函数
)()(2
1)(min 22μμμμμμc b a c b a ++++=T F
对目标函数取μ的的一阶导数，并令其为零，得
0332210=+++μμμg g g g （4-55） 其中
???????==+==c c c b c
a b b b a T T
T T T g g g g 2324
210 （4-56） 解（4-55）一元三次方程得最优乘子μ。