潮流计算的快速分解法(可编辑修改word版)

潮流计算的快速分解法
摘要：本文采用快速分解法进行潮流计算，分析其基本理论，并使用 MATLAB 软件进行编程设计。

最后运用实例进行验证。

结果表明快速分解法具有较好的迭代速度。

关键词：潮流计算快速分解法 MATLAB 编程，实例验证
1引言
潮流计算是电力系统分析最基本、最重要的计算，是电力系统运行、规划以及安全性、可靠性分析和优化的基础，也是各种电磁暂态和机电暂态分析的基础和出发点。

潮流计算要求具有可靠的收敛性，占用内存少，计算速度快，调整和修改容易，使用灵活方便。

各种算法的改进以及新算法的提出，很多都是为了使潮流计算能更好地满足计算要求。

本文应用快速分解法进行潮流计算，并给出算例分析。

2潮流计算的快速分解法
研究表明，用牛顿-拉夫逊法计算潮流时，每次迭代都要重新形成雅可比矩阵，然后重新对它进行因子表分解并求解修正方程。

为避免每次迭代重新形成雅可比矩阵及其因子表，人们研究用定雅可比矩阵取代随迭代过程不断变化的雅可比矩阵，这种方法叫定雅可比法。

此外，人们还结合电力系统的物理特点，发展了各种版本的解耦潮流算法，20 世纪 70 年代初提出的快速分解法是这一阶段的主要研究成果。

关于快速分解潮流算法，有三项里程碑意义的研究成果。

其一是 Stott 在1974年发现的 XB 型算法；其二是 Van Amerongen 在1989 年发现的 BX 型算法；
其三是Monticelli 等人在1990 年所作的关于快速分解潮流算法收敛机理的理论阐述。

这些研究工作不仅是电力系统计算方面的典范，也揭示了这样一个事实：工程上有效的方法一定有其深刻的理论来支持。

2.1快速分解法的修正方程及迭代格式
将极坐标型定雅可比法的修正公式重写如下：
H L
- ⎡ B H - G N ⎤⎡V ∆⎤ = ⎡∆P V ⎤ ⎢G B ⎥⎢ ∆V ⎥ ⎢∆Q V ⎥
（2.1）
⎣ M L ⎦⎣ ⎦ ⎣
⎦ 经验表明，电力系统中有功功率主要受电压相角的影响，而无功功率主要受电压幅值的影响，同时由于高压电网大部分线路的电阻比电抗小，因此在牛顿- 拉夫逊迭代中可以忽略雅可比矩阵的非对角块，即将G N ， G M 设为零，从而实 现有功和无功潮流修正方程的解耦。

Stott 通过大量的计算实践发现，为了获得最好的收敛性，还要对雅可比矩阵的对角块作特殊的常数化处理：对系数矩阵 B H
，忽略支路电阻和接地支路的影响，即用- 1 x 为支路电纳建立的节点电纳矩阵 B '
代替 B ；对系数矩阵 B ，用节点导纳矩阵中不包含 PV 节点的虚部 B '' 代替；
V ∆
前的电压幅值用标幺值 1 代替。

于是可得简化的修正方程式如下：
- B '∆
= ∆P V
（2.2）
- B ''∆V = ∆Q V
（2.3）
在潮流计算中，上述两个修正方程式依次交替迭代，Stott 把在此基础发展起来的潮流算法称为快速分解法（fast decoupled load flow ）。

假定当前点为 (
(k ) ,V (k ) ) ，则求解((k +1) ,V (k +1) ) 的连续迭代格式如下：
⎧∆V (k ) = - B ''-1∆Q ((k ) ,V (k ) ) V (k )
（2.4）
⎨ (k +1)
⎩V
= V (k ) + ∆V (k ) ⎧∆(k ) = - B '-1
∆P ((k )
,V (k +1) ) V (k +1)
⎨ (k +1) (k ) (k )
（2.5）
⎩ = + ∆
快速分解法公式的特点是：① P -
和Q - V 迭代分别交替进行；②功率偏
差计算时使用最近修正过得电压值，且有功无功偏差都用电压幅值去除；③ B ''
和B ' 的构成不同，B '应用-1x 建立，并忽略所有接地支路（对非标准变比变压
器支路，变比可取为 1），而B ''就是导纳矩阵的虚部，不包括PV 节点。

在快速
分解法的实施中，这些技术细节缺一不可，否则程序的收敛性将受到影响。

1989 年，荷兰学者 Van Amerongen 通过大量仿真计算发现了另一版本的快速分解潮流算法，他把该算法称为 BX 型算法，而把 Stott 的算法称为 XB 型算法，用以区分二者。

BX 型算法与 XB 型算法的主要不同在于雅可比矩阵对角块的形成
进行简化时，保留了支路电阻的上。

BX 型算法的处理方式是：在对系数矩阵B
H
影响，但忽略了接地支路项。

BX 型算法的迭代格式与 XB 型算法是相同的。

计算
经验表明，BX 型和 XB 型两种快速分解潮流算法在大部分情况下性能接近，在某
些情况下 BX 型算法收敛性略好。

快速分解法只对雅可比矩阵作了简化，但节点功率偏差量的计算及收敛条
件仍是严格的，因此收敛后的潮流结果仍然是准确的。

由于方程的维数减小了，
且B' 和B '' 是常数矩阵，只需在迭代计算之前形成一次，然后分解成因子表，并
一直在迭代过程中使用，所以计算效率大幅提高。

快速分解法是一种定雅可比法，虽然只具有线性收敛速度，但由于其鲁棒性好，适应性强，在电力工业界被广泛采用，特别适合在线计算。

2.2快速分解法的理论基础
Stott 的快速分解法提出时并没有任何理论解释，它是计算实践的产物。

多年来，人们普遍认为在满足r <<x 的系统中，快速分解法才能有较好的收敛性。

但
在许多实际应用中，当r >x 时，快速分解法也能很好收敛。

因此，从理论上解
释快速分解法的收敛机理，便成为一个有趣的研究课题。

20 世纪 80 年代末，Monticelli 等人的研究工作对这一问题做了比较完整的解释，在一定程度上阐明
了 XB 型和BX 型快速分解潮流算法的收敛机理。

Monticelli 等人的分析工作是以定雅可比牛顿-拉夫逊迭代方程为出发点的。

具体过程如下：①通过高斯消去法，把牛顿-拉夫逊法的每一次迭代等价地
细分为三步计算；②对每一步计算作详细分析，证明了在连续的两次牛顿-拉夫
逊迭代中，上一次迭代的第三步和下一次迭代的第一步可以合并，从而导出等效的
两步式分解算法；③论证了该两步式分解算法的系数矩阵与快速分解法的系数
L
M
H ⎩ P P L 矩阵是一致的。

推导过程并未引用任何解耦的假设。

为以后书写方便，将式（2.1）中的∆P V 用∆P 代替， ∆Q V 用∆Q 代替， 而V ∆用∆代替，则给出的定雅可比法的修正公式改写如下
- ⎡ H N ⎤⎡∆⎤ = ⎡∆P ⎤ ⎢M L ⎥⎢∆V ⎥ ⎢∆Q ⎥
（2.6）
式中
⎣ ⎦⎣ ⎦ ⎣ ⎦ H = B H ≈
∂∆P
， N = -G ∂T N
≈ ∂∆P
∂V T M = G M
≈ ∂∆Q ， L = B ∂T L
≈ ∂∆Q ∂V T 整个推导分为三步。

1) 将原问题分解成 P ， Q 子问题
首先，对式（2.6）用高斯消去法消去子块 N ，有
- ⎡H - NL -1M 0⎤⎡∆⎤ = ⎡∆P - NL -1∆Q ⎤ ⎢ M L ⎥⎢∆V ⎥ ⎢ ∆Q ⎥
（2.7）
记
⎣ ⎦⎣ ⎦ ⎣ ⎦
~ -1 ~ -1 H = H - NL M ， ∆P = ∆P - NL ∆Q
并定义
∆V = -L -1∆Q ， ∆V
= -L -1M ∆
则式（2.7）的解可以表示为
⎧∆= - ~
-1∆P
⎨
∆V = ∆V + ∆V
M
上式中对∆~
的计算可以采用较简单的方法。

在给定的电压幅值和相角初值附近，保持电压相角不变，考虑只有电压幅值的变化∆V L 时，有功功率的偏差量 为
∆P (,V + ∆V L
) ≈ ∆P (,V ) + ∂∆P
∆V ∂V T
L
= ∆P (,V ) - NL -1∆Q = ∆~
⎩
H ∆P ( ,V （2.8）
综合上述结果，如果当前的迭代点为(
(k ) ,V (k ) ) ，则第k 次迭代对式（2.6）
的计算可以分解为以下三步。

①
⎧∆V (k ) = -L -1∆Q ((k ) ,V (k ) ) ⎨ ~ L V (k +1) = V (k ) + ∆V (k )
（2.9）
②
⎩ L
⎧⎪∆(k ) = - ~ -1 (k )
~(k +1)
（2.10）
③
⎨⎪(k +1) = (k ) + ∆(k )
⎧∆V (k ) = -L -1 M ∆(k )
⎨ M ~ V (k +1) = V (k +1) + ∆V (k ) ⎩ M
（2.11）
2) 简化无功迭代步骤
按①～③完成第k 次迭代后，下面再考察第k + 1次迭代的①，有
⎧∆V (k +1) = -L -1∆Q ((k +1) ,V (k +1) ) ⎨ ~ L
V (k +2) = V (k +1) + ∆V (k +1) ⎩ L
（2.12）
利用式（2.11），上式中的无功功率偏差为
∆Q ((k +1) ,V (k +1) ) = ∆Q ((k +1)
~(k +1)
+ ∆V (k ) ) ,V M
≈ ∆Q ((k +1)
~(k +1) ) + ∂∆Q ∆V (k )
,V ∂V T
M = ∆Q ((k +1)
~(k +1)
) + L ∆V (k ) ,V M
（2.13）
代入式（2.12），经整理得
∆V (k +1) + ∆V (k ) = -L -1
∆Q (
(k +1)
~
(k +1)
)
L M
,V
（2.14）
)
V H ⎩
H H
H H 式（2.14）说明，如果将第k 次迭代的①计算出的 ~
(k +1) 和②计算出的
(k +1) ，
用于计算第k + 1次迭代的无功偏差量，即式（2.14）中的∆Q ，则所求得的第k + 1
次迭代的电压修正量将自动包含第k 次迭代的③的式（2.11）与第k + 1次迭代的①的式（2.12）合并，只需保留式（2.9）和式（2.10）。

因此，第k 次迭代对式（2.6）的计算可以用以下两步计算完成：
⎧⎪∆V (k ) = -L -1∆Q ((k ) ,V (k ) )
（2.15）
⎨ ⎪⎩
V (k +1)
= V (k )
+ ∆V (k ) ⎧⎪∆(k )
= - ~ -1 ∆P ((k )
,V (k +1) )
（2.16）
⎨⎪(k +1) = (k ) + ∆(k )
在式（2.6）处已说明， ∆P 实际是∆P V ， ∆Q 实际是∆Q V ， ∆实际是 V ∆，式（2.15）和式（2.16）和快速分解法迭代格式相同。

显然，这种迭代算法是否与快速分解法等效，取决于系数矩阵 L 和 ~。

与 XB 型快速分解法的修正 方程相比，系数矩阵 L 是导纳矩阵的虚部，这与 B
'' 相同，所以关键要看
~
是否
与 B ' 有相同或相似的关系。

3) 简化有功迭代矩阵 ~
~
H 的定义，有
~ -1 -1
H = H - NL M = B H +G N B L G M
（2.17）
对于一般的电网， ~ 可能有较复杂的结构。

为了对 ~
有直观的认识，假定
H H
网络中无 PV 接点，则式（2.17）中各矩阵的维数相等，并且接点导纳矩阵可用节点支路关联矩阵 A 和支路导纳对角矩阵（分别用的b 和 g 表示电纳和电导）表 示。

下面将证明，对于树形电网或所有支路的r x 比值都相同的环形网络，
~
与
B ' 相等。

由
~
H
~
如果网络是树状的，其关联矩阵 A 是方阵且非奇异，此时对式（2.17）有
H = AbA T + ( A gA T )( A bA T )-1 ( AgA T )
= A (b + gA T A -T b -1 A -1 Ag ) A T
= A (b + gb -1b ) A T
= Ab ' A T
= B '
（2.18）
式中， b
' 为以- 1 x 为支路电纳组成的对角线矩阵； B ' 为以- 1 x 为支路电纳建立的节点电纳矩阵。

这说明对树形电网， ~
就是 XB 型快速分解
法中的 B ' 阵。

对于环形网络，如果电网是均一网，即对任一支路l 有r l x l = ，则得
g l = r l
r 2 + x 2 = x l r 2 + x 2 = -b l
l
l
l
l
并有
- r 2
x 1 b + g b 1
g = (1 +2 )b = (1 + l
)(- l ) == - l l l l l x 2 r 2 + x 2 x
l l l l
所以
g = -b ， (1 +
2
)b = b '
故有
H = AbA T + ( A gA T
)( A bA T
)-1 ( AgA T )
如果电
元素处，相应 H H
= AbA T +2
( AbA T )( AbA T )-1 ( AbA T )
= (1 +
2 ) AbA T
= Ab ' A T = B '
（2.19）
网不是均一网，上述结论不再严格成立。

但 ~
和 B ' 相比，在 B ' 的零
H
~ 的元素近似等于零；在 B ' 的非零元素处，相应 ~
的元素近似和 B '
H ⎣
2 ⎣ 2 H
H 1⎦⎣ ⎦ ⎦
的非零元素相等。

这可以用下面的例子来说明。

图 2.1 四节点电力系统
以图 2.1 所示的四节点系统为例，图中给出了支路阻抗。

该例中 H ， ~
和
B ' 分别为
⎡ ⎢
H = ⎢
⎢
1.5 - 1 0
- 1 1.2 - 0.2
0 - 0.2 0.7
- 0.5⎤
⎥ ⎥ - 0.5⎥ ⎢- 0.5 0 - 0.5 ⎥ ⎡ 1.9 - 0.8 - 0.1 - 1⎤ ⎡ 2 - 1 0 - 1⎤ ~ ⎢- 0.8 1.6 - 0.8 0⎥
⎢- 1 2 - 1 0⎥
H = ⎢ ⎢- 0.1 - 0.8 1.9 ⎥ - 1⎥ B ' = -⎢
⎢ 0 - 1 ⎥ 2 - 1⎥ ⎢ - 1 0 - 1 ⎥ ⎢- 1 0 - 1 ⎥
可见， B ' 比 H 更接近于
~
，而用 B ' 代替
~
即得到 XB 型快速分解法。

3 程序流程
计算时，只要有电路图及其各支路阻抗初始值、节点类型，即可进行迭代计算得出结果。

利用快速分解法编程流程图如下。

max ∆P(k)≤
否
max ∆Q (k)≤
i
i
i
i
i
计算相角修正量∆(k ) ，求得(k +1) =∆(k ) +(k )
4程序
function [Y,U,P,Q,deltaSij,a,S,Sij,Sji,sumdeltaS]=PQ()
%电力系统潮流计算程序；
%输出：U——节点电压，P--节点有功，Q--节点无功，deltaSij--支路功率损耗，
%Sij--从节点 i 流向节点 j 的功率,S--节点复功率,sumdeltaS--网络总损耗
%输入参数：point 为节点信息矩阵，branch 为支路信息矩阵；
[x]=3;%节点数 x
[y]=3;% 支路数 y
e=0.00005;%误差要求
point= [1 1 0 -2 -1;2 1.01 0 0.5 0.25;3 1 0 0 0];%从exel 中读取节点信息矩阵
branch=[1 2 0.01 0.2 0.02 0 0 0;1 3 0 0 0 0.01 0.1 1.05;2 3 0.02 0.2
0.04 0 0 0];%从 exel 中读取支路信息矩阵
TYPE=zeros(x,1);%TYPE 为节点类型矩阵
U=zeros(x,1);%U 为节点电压矩阵
a=zeros(x,1);%a 为节点电压相角矩阵
P=zeros(x,1);%P 为节点有功功率
Q=zeros(x,1);%Q 为节点无功功率
I=zeros(y,1);%I 为起始节点编号矩阵
J=zeros(y,1);%J 为终止节点编号矩阵
Rij=zeros(y,1);%R 为线路电阻
Xij=zeros(y,1);%X 为线路电抗
Zij=Rij+j*Xij;%Yij 为线路阻抗Y=zeros(x);%Y
为n 阶节点导纳方阵G=zeros(x);%G 为n 阶节
点电导方阵B=zeros(x);%B 为n 阶节点电纳方
阵B0=zeros(y,1);%B0 为n*1 阶线路对地电纳
值
RT=zeros(y,1);%RT 为ij 支路y( 矩阵 branch 的行数)*1 阶变压器电阻
XT=zeros(y,1);%XT 为ij 支路y*1 阶变压器电抗
ZT=RT+j*XT;%求变压器阻抗
KT=zeros(y,1); %K 为ij 支路y*1 阶变压器变比,若k=0 表示无变压器,K=1 则为标准变比,k 不等于 1 为非标准变比
% ----------------------------- 矩阵赋初值：
TYPE=point(:,1);%将 point 矩阵的第一列赋给 TYPE,以下类似
U=point(:,2);
a=point(:,3);
P=point(:,4);
Q=point(:,5);
I=branch(:,1);
J=branch(:,2);
Rij=branch(:,3);
Xij=branch(:,4);
Zij=Rij+j*Xij;
B0=branch(:,5);
RT=branch(:,6);
XT=branch(:,7);
ZT=RT+j*XT;
KT=branch(:,8);
% ----------------------------- 求节点导纳矩阵 Y
for m=1:y %求 Y 中非对角元元素 Yij
if KT(m)==0%若无变压器,则 Yij 直接为线路阻抗分之一取负值． Y(I(m),J(m))=-1/Zij(m);
Y(J(m),I(m))=-1/Zij(m);
else %有变压器时，Ｙij 为线路阻抗乘以ＫＴ后分之一再取负值Y(I(m),J(m))=-1/(KT(m)*ZT(m));
Y(J(m),I(m))=-1/(KT(m)*ZT(m));
end
end
for m=1:x %求 Y 中的Ｙii
for n=1:y
if KT(n)==0 %无变压器时Ｙii 为Ｙij 加上线路对地电导的一半乘 j
if(I(n)==m|J(n)==m)
Y(m,m)=Y(m,m)-Y(I(n),J(n))+j*B0(n)/2;
end
elseif I(n)==m %有变压器时，若支路起始节点为 m,则变压器等值模型的对地支路的(1-KT)/KT^2 算到Ｉ(m)节点
Y(m,m)=Y(m,m)-Y(I(n),J(n))+(1-
KT(n))/(KT(n)^2)*(1/ZT(n));
elseif J(n)==m %有变压器时，若支路终止节点为 m,则变压器等值模型的对地支路的(KT-1)/KT 算到J(m)节点
end end
Y
Y(m,m)=Y(m,m)-Y(I(n),J(n))+(KT(n)-1)/KT(n)*(1/ZT(n)); else Y(m,m)=Y(m,m);
end
G=real(Y);
% ---------------------- 求B'矩阵及其逆矩阵 B1
B=imag(Y);%求Y 的虚部,节点电纳矩阵
[y1]=[y-1];
[x1]=[x-1];
B11=zeros(x1,x1);
for m=1:y1
B11(I(m),J(m))=1/(Xij(m)+XT(m));
B11(J(m),I(m))=1/(Xij(m)+XT(m));
end
for m=1:x1
for n=1:y
if(I(n)==m|J(n)==m)
B11(m,m)=B11(m,m)-1/(Xij(n)+XT(n));
end
end
end
ph=find(TYPE(:,1)==3);%找出平衡节点编号
B11(:,ph)=[];%平衡节点编号对应行置空
B11(ph,:)=[];%平衡节点编号对应列置空
B11
B1=B11;
B1=inv(B1);%B1 求逆后得 B1 矩阵
% ---------------------- %求 B''及其逆矩阵 B2 phpv=find(TYPE(:,1)>1);%找出非PQ 节点的编号
B22=B; %BB 矩阵为中间变量
B22(:,phpv)=[];%非PQ 节点编号对应行置空
B22(phpv,:)=[];%非PQ 节点编号对应行置空
B22
B2=B22;
B2=inv(B2);%求得 B2 矩阵
% ------------ 计算各节点有功功率不平衡量 deltaPi
k=0; %k 为迭代次数
kp=0; %计算P 不平衡量 deltaPi 的收敛标志(0 表示不收敛,1 表示收敛) kq=0; %计算U 不平衡量 deltaQi 的收敛标志(0 表示不收敛,1 表示收敛) notph=find(TYPE(:,1)<3);%找出非平衡节点编号
deltaPi=zeros(x-1,1);%deltaPi 为x*1 阶矩阵,x 即为节点数
pq=find(TYPE(:,1)==1);%找出 PQ 节点编号
pqnum=size(B2);
pqnum=pqnum(1);%求PQ 节点的个数(因 B1 矩阵的维数等于 PQ 节点数) deltaQi=zeros(pqnum,1);%deltaQi 为pqnum*1 阶矩阵
while((~kq|~kp)&(k<100))
k=k+1;
for m=1:(x-1)%求 deltaPi
sum1=0;
for n=1:x
sum1=sum1+U(notph(m))*U(n)*(G(notph(m),n)*cos(a(notph(m))- a(n))+B(notph(m),n)*sin(a(notph(m))-a(n)));
end
deltaPi(m)=P(notph(m))-sum1;
end
Up=U; %Up 为中间变量Up(ph)=[];%
将平衡节点所在行置空
Unotph=Up;%求得除平衡节点外的电压列向量
deltaa=(-B1*(deltaPi./Unotph))./Unotph;%求相角a 的不平衡量
for m=1:(x-1) %求相角a 的新迭代值矩阵
a(notph(m))=a(notph(m))+deltaa(m);
end
max1=abs(deltaPi(1)/U(notph(1)));%求deltaP/U 绝对值的最大值for m=1:(x-2)
if
abs(deltaPi(m)/U(notph(m)))<abs(deltaPi(m+1)/U(notph(m+1))) max1=abs(deltaPi(m+1)/U(notph(m+1)));
end
end
if max1<=e %如果最大值满足要求,则 kp 置为"1",表示收敛kp=1;
end
for m=1:pqnum %求 deltaQi
sum2=0;
for n=1:x
sum2=sum2+U(pq(m))*U(n)*(G(pq(m),n)*sin(a(pq(m))- a(n))-B(pq(m),n)*cos(a(pq(m))-a(n)));
end
deltaQi(m)=Q(pq(m))-sum2;
end
Uq=U;%Uq 为中间变量
Uq(phpv)=[];%将非 PQ 节点所在行置空
Upq=Uq;%求得包括 PQ 节点电压的电压列向量
deltaU=-B2*(deltaQi./Upq);% 求 U 的不平衡量 deltaU
max2=max(abs(deltaQi./Upq)); %求deltaQ/U 绝对值的最大值
if max2<=e %如果最大值满足要求,则kq 置为"1",表示收敛kq=1;
end
for m=1:pqnum %求 U 的迭代新值
U(pq(m))=U(pq(m))+deltaU(m);
end
end
sum3=0+j*0;%求平衡节点功率 Sph
for m=1:x
sum3=sum3+conj(Y(ph,m))*(U(m)*cos(a(m))-i*U(m)*sin(a(m))); end
Sph=(U(ph)*cos(a(ph))+j*U(ph)*sin(a(ph)))*sum3;%求平衡节点功率fid=fopen('PQj.txt','wt'); %打开输出文件,PQ.txt,向其写入结果数据fprintf(fid,'潮流计算输出结果\n');
fprintf(fid,'迭代次数k 为: %d \n',k);
fprintf(fid,'节点电压: \n');
for m=1:x
fprintf(fid,'第%d 个节点: %f\n',m,U(m));
end
fprintf(fid,'节点相角: \n');
for m=1:x
fprintf(fid,'第%d 个节点: %f\n',m,a(m));
end
fprintf(fid,'节点有功功率: \n');
for m=1:x
fprintf(fid,'第%d 个节点: %f\n',m,P(m));
end
fprintf(fid,'节点无功功率: \n');
for m=1:x
fprintf(fid,'第%d 个节点: %f\n',m,Q(m));
end
fclose(fid); %关闭文件,程序结束
5算例验证
5.1求解系统的节点导纳矩阵
如图5.1（a）所示的一个三母线电力系统，在母线①和母线③之间的输电线的母线③端连接有一个纵向串联加压器，可在同一电压等级改变电压幅值。

该系统的网络原件用图（b）所示的等值电路表示，串联支路用电阻和电抗表示，并联支路用电纳表示。

支路（1,3）用一个变比可调的等值变压器支路表示，非标准变比t 1.05 ，在节点①侧。

求解该网络的节点导纳矩阵。

（a）
（b）
⎢ ⎢ = ⎥ ⎥⎢ ⎥⎢ ⎥
1 ⎥
⎥ t
T
⎢ ⎢
（c ）
图 5.1 三母线电力系统
解 首先对支路编号并规定串联支路的正方向如图 5.1（c ）所示，则广义节点支路关联矩阵是
⎡- 1 - 1 1 ⎤
t A = ⎢ 1 - 1 ⎥
⎢ 1 1 1⎥
⎢⎣ ⎥⎦
A 中行与节点对应，列与支路对应。

支路导纳是
y 12 y 13 y 23 = 1
0.01+ j 0.2
= 1
0.01+ j 0.1
= 1
0.02 + j 0.2 = 0.2494 - j 4.9875
= 0.9901- j 9.9010
= 0.49505 - j 4.9505
y 10 = y 20 = y 30 = j 0.01 j 0.03 j 0.02
建立节点导纳矩阵如下：
⎡ y 12
⎤⎡ - 1 1
⎤
⎡
1
⎤ ⎢
y
⎥⎢- 1 ⎥
⎢- 1 - t 1 ⎢ 13
y
⎥⎢ 1⎥ ⎥⎢ ⎥
Y = Ay b A ⎢ - 1 1 ⎥ ⨯ ⎢ 23
y
⎥⎢ - 1 1⎥
⎥ ⎢
1 1 ⎢⎣ 1⎥ ⎢ ⎦⎥ ⎢ ⎣
10
y 20
⎥⎢ 1 ⎥ y 30 ⎦⎣
⎥ 1 ⎥ 1⎥⎦
1
⎢ ⎥
于是，Y 的各元素是
Y = y +
y 13 + y
11
12
t
2
10
= 0.2494 - j 4.9875 + = 1.1474 - j 13.9580 Y 22 = y 12 + y 23 + y 20
1
1.052
(0.9901 - j 9.9010) + j 0.01
= 0.2494 - j 4.9875 + 0.49505 - j 4.9505 + j 0.03 = 0.74445 - j 9.908 Y 33 = y 13 + y 23 + y 30
= 0.9901 - j 9.901 + 0.49505 - j 4.9505 + j 0.02 = 1.48515 - j 14.8315
Y 12 = Y 21 = - y 12 = -0.2494 + j 4.9875
Y 13 = Y 31 = - y 13 t =
- 0.9901 + j 9.901
= -0.9430 + j 9.43 1.05
Y 23 = Y 32 = - y 23 = -0.49505 + j 4.9505
导纳矩阵是对称矩阵，只需求出上三角部分。

最后求得导纳矩阵为
⎡ 1.1474 - j 13.9580 Y = ⎢- 0.2494 + j 4.9875 - 0.2494 + j 4.9875
0.74445 - j 9.908 - 0.9430 + j 9.430⎤ - 0.49505 + j 4.9505⎥ ⎣⎢
- 0.9430 + j 9.430 - 0.49505 + j 4.9505 1.48515 - j 14.8315⎥⎦ 4.2 用快速分解法计算潮流
如图 5.2 所示的电力系统，其导纳矩阵在上面已经求出。

如果给定各节点的发电和负荷功率以及节点电压，写出极坐标形式的潮流方程，并用快速分解法计算潮流。

已知 P D 1 + jQ D 1 = 2 + j 1， P D 2 + jQ D 2 = 0.5 + j 0.25 ， P G 2 = 1，V 2 = 1.01 ， V 3∠3 = 1∠0 。

⎩ ⎦
D 1 1 11 1 2 12 12 12 D 1 5 ⎣ 1
图 5.2 三母线电力系统
解 如图 5.2 所示节点①是 PQ 节点，节点②是 PV 节点，节点③是V 节点。

待求的状态变量是 x = [ 1
2
V ] ，共有两个有功潮流方程和一个无功潮流方
T
程：
⎧∆P 1 = -P D 1 - P 1 (V ,) = 0 ⎪∆P = P - P - P (V ,) = 0 ⎨ 2 G 2 D 2 2
⎪∆Q = -Q - Q (V ,) = 0
1
其具体表达式如下：
⎧∆P = -P - V 2G - V V (G cos + B sin
) ⎪
1 D 1 1 11 1
2 12 12 12 12 ⎪
- V 1V 3 (G 13 cos 13 + B 13 sin 13 ) = 0 ⎪ ⎪
∆P = P - P - V 2G - V V (G cos + B sin
)
⎪ 2
G 2 D 2 2 22 2 1 21 21 21 21
⎨
- V 2V 3 (G 23 cos 23 + B 23 sin 23 ) = 0 ⎪ ⎪
⎪∆Q = -Q + V 2
B - V V (G sin - B cos ) ⎪ - V V (G sin - B cos ) = 0
⎪
1 3 13 ⎪⎩
13 13 13
这个系统的节点导纳矩阵在 5.1 中已求出，为
⎡ 1.1474 - 0.2494 - 0.9430⎤ ⎡- 13.9580 4.9875 9.430⎤ Y = G + jB = ⎢- 0.2494 0.74445 - 0.49505⎥ + j ⎢ 4.9875 - 9.908
4.9505⎥ ⎢ ⎢⎣- 0.9430 - 0.49505 ⎥ ⎢ 1.48515⎥⎦ ⎢⎣
9.430 4.9505 ⎥ - 14.8315⎥⎦
快速分解法中的 B ' 和 B '' 分别为
B ' =
⎡- 15
5⎤ , B ''
= [- 13.9580] ⎢ - 10⎥
其中 B ' 是由- 1 x 建立所得，即通过
⎧
B (i , i ) = ∑ 1 ⎪
j ∈i , j ≠i x ij ⎨ 1
求出。

⎪B (i , j ) = - ⎪⎩
x ij
而 B '' 则直接从导纳矩阵虚部中取。

将Y 矩阵的具体数值代入功率偏差方程
有
12
⎣
⎤ [
1 2 1 1 2 12 12 ⎧∆P = -2 - 1.1474V 2 - V V (-0.2494 cos + 4.9875sin )
⎪
1 1 1
2 12 12 ⎪
⎪ ⎪∆P - V 1V 3 (-0.9430 cos 13 + 9.430 s in 13 ) = 0 = 1 - 0.5 - 0.74445V 2 - V V (-0.2494 c os + 4.9875sin )
⎪ 2
2 2 1 21 21 ⎨
- V 2V 3 (-0.4951cos 23 + 4.951sin 23 ) = 0 ⎪ ⎪
⎪∆Q = -1 + 13.9580V 2
- V V (-0.2494 sin - 4.9875 cos ) ⎪ - V V (-0.9430 s in - 9.4300 c os ) = 0
⎪ 1 3 13 13
⎪⎩
从平启动开始，所有 PQ 节点电压为1∠0 ，这时电压初值为
⎡ • •
⎢V 1 V 2 •
V 3 ⎥⎦ = 1∠0 1.01∠0 1∠0 ]
代入有功潮流方程中计算第一次迭代时的有功不平衡量
∆P (0)
= -2 -
1.1474 - 1.01⨯ (-0.2494) + 0.9430 = -1.9525 ∆P (0) = 0.5 - 0.74445 ⨯1.012 - 1.01⨯ (-0.1494) - 1.01⨯ (-0.4951) = 0.49248
所以
⎡∆P ⎤
(0)
⎡ ∆P V ⎤ (0) ⎡- 1.9525⎤ ⎢ ⎥ = ⎢ 1 ⎥ = ⎢ ⎥ ⎣ V ⎦
计算相角修正量 ⎣∆P 2 V 2 ⎦ ⎣ 0.48761⎦
⎡∆P ⎤
(0)
⎡- 15
5⎤-1
⎡- 1.9525⎤ ⎡- 0.1367 ⎤ ∆(0) = -B '-1 ⎢ ⎥
= -⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎣ V ⎦
⎣ 5 - 10⎦ ⎣ 0.48761⎦ ⎣- 0.01959⎦
(1)
=
(0)
+ ∆
(0)
= -⎡0⎤ + ⎡ - 0.1367⎤ = ⎡- 0.1367 ⎤
⎢0⎥ ⎢- 0.01959⎥ ⎢- 0.01959⎥ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦
式中， ∆前得电压幅值项取值为 1。

然后进行无功迭代。

先计算节点无功不平衡量。

这时节点电压相角用最新修正后的值
⎡1 ⎤ ⎡ - 0.1367 ⎤
V (0)
= ⎢1.01⎥ ，(1)
= ⎢- 0.01959⎥ ⎢ ⎥ ⎢⎣1 ⎥⎦
⎢ ⎥
⎣⎢
0 ⎥⎦
只有一个 PQ 节点，故只有一个无功平衡方程
1
1 ∆Q (0) = -1 - 13.9580 - 1.01⨯[-0.2494 ⨯ sin(-0.1171) - 4.9875 ⨯ cos(-0.1171)]
-[-0.9430 ⨯ sin(-0.1367) - 9.430 ⨯ cos(-0.1367)]
= -14.9580 - 1.01⨯ (-4.92421) + 9.21352 = -0.7715
(
∆Q 1 )(0) = -0.7715
V 1
再计算电压修正值：
∆V (0) = -B ''-1
(
∆Q 1 )(0) =
1
⨯ (-0.7715) = -0.055273 V 1
13.958
V (1) = V (0) + ∆V (0) = 1 - 0.055273 = 0.94473
1
1
1
再利用新计算出的电压转入第二次有功迭代。

整个迭代过程中功率不平衡量、电压幅值和角度的修正量以及电压幅值和角度的变化量列在表 5.1 中，其中∆P 1 、∆P 2 、∆Q 1 的列用于判断收敛。

表 5.1 快速分解法潮流计算迭代过程（= 0.00005 ）
迭代 次
数
1 -1.95260
0.49248 0.13670 0.01959 -0.13670 -0.01959 -0.77151 0.05527 0.94473
2 -0.13518
0.01872 0,.01071 0.00350 -0.14741 -0.02309 -0.08405 0.00637 0.93835
3 -0.01486
0.00431 0.00110 0.00012 -0.14851 -0.02321 -0.00904 0.00076 0.93759
4 -0.00167
0.00060 0.00012 0.00000 -0.014862 -0.02321 -0.00176 0.00009 0.93750
5 -0.00019
0.00007 0.00001 0.00000 -0.14854 -0.02321 -0.00014 0.00001 0.93749
6 -0.00002
0.00001
-0.00001
由表 5.1 可知第六次迭代时用功功率和无功功率不平衡量已小于预先制定的收敛阙值= 0.00005 ，因此不用继续做电压相角和幅值的修正。

潮流计算程序结果截图如图 5.3 所示：
1
∆P 1 ∆P 2 ∆1 ∆2
1 2
∆Q 1 ∆V 1 V 1
（a）
（b）
图5.3 P-Q 潮流计算程序输出结果截图
6参考文献
[1]张伯明，陈寿孙，严正.高等电力网络分析[M].北京：清华大学出版社，2007.
[2]韩平.潮流计算方法的实用化改进[D].北京：华北电力大学,2009.
[3]刘国栋.电力网络潮流计算的模块开发[D].吉林：吉林大学,2009.
[4]王忠礼，段慧达，高玉峰.MATLAB 应用技术——在电气工程与自动化专业中
的应用[M].北京：清华大学出版社，2007.
[5]龚纯，王正林.MATLAB 语言常用算法程序集[M].北京：电子工业出版社，
2008.。