牛拉法潮流计算程序(附3机9节点结果对比).docx

摘要

电力系统潮流计算是研究电力系统稳态运行的一种重要方法,它根据给定的运行条件及系统接线情况确定整个电力系统各部分的运行状态,包括各母线的电压、线路的功率分布以及功率损耗等等。潮流计算主要用于电网规划和静态安全分析，它可为扩建电力网络，以达到规划周期内所需要的输电能力提供依据；也可以对预想事故进行模拟和分析，校核预想事故下的电力系统安全性。
本文简单介绍了牛顿-拉夫逊潮流计算的原理、模型与算法，然后用具体的实例，利用MATLAB对牛顿-拉夫逊法的算法进行了验证。

关键词：电力系统潮流计算 牛顿-拉夫逊法 MATLAB
一、牛拉法的数学模型
对一个N节点的电力网路，列写节点电压方程，即
I=Y V (1.1)
式中，I为节点注入电流列相量，Y为节点导纳矩阵，V为节点电压列相量。
由于异地测量的两个电流缺少时间同步信息，以注入功率替换注入电流作为已知量。即
(1.2)
其中,Yij=Gij+jBij，带入上式，得到有功功率和无功功率方程
Pi=Vij=1nVjGijcosθij+Bijsinθij (1.3)
Qi=Vij=1nVjGijsinθij-Bijcosθij (1.4)
大部分情况下，已知PQ，求解Vθ。考虑到电网的功率平衡，至少选择一台发电机来平衡全网有功功率，即至少有一个平衡节点，常选择调频或出线较多的发电机作为平衡节点。具有无功补偿的母线能保持电压幅值恒定，这类节点可作为PV节点。
潮流计算中节点分类总结如下：

表 1.1 节点分类

PQ节点
负荷母线、变电站母线
绝大部分是PQ节点
已知PQ，求Vθ

PV节点
有无功储备的发电机母线
很小部分是PV节点
已知PV（V=1），求Qθ

平衡节点
选择一个发电机母线
已知Vθ（V=1,θ=0）,求PQ
电压相量作为参考轴θ=0

已知电力系统有m个PQ节点，r个PV节点和1个平衡节点，则可以提取m+r个有功功率方程和m个无功功率方程，从而求解出m+r个θ和m个V，其余节点的有功和无功可通过式（1.3）、（1.4）求得，这样就完成了潮流计算。
二、潮流计算流程
用于潮流计算的有功和无功功率方程的阶数一般很高，常选择牛顿-拉夫逊迭代法进行求解高维的非线性方程组。
结束
maxF(Xk)<ε？

Xk

Xk+1=Xk+∆Xk
maxΔXk<ε？
ΔXk=-J-1XkFXk
牛顿-拉夫逊法的迭代过程如下：

图2.1 牛拉法迭代流程图
将功率方程改成Fx=0的形式
ΔPi=Pis-Vij=1nVjGijcosθij+Bijsinθij=0 (2.1)
ΔQi=Qis-Vij=1nVjGijsinθij-Bijcosθij=0 (2.2)
可以得到
ΔθjΔVj=-∂ΔPi∂θj∂ΔPi∂Vj∂ΔQi∂θj∂ΔQi∂Vj-1ΔPiΔQi (2.3)
为了便于计算，将上式改为
ΔθjΔVjVj=-∂ΔPi∂θj∂ΔPi∂VjVj∂ΔQi∂θj∂ΔQi∂VjVj-1ΔPiΔQi (2.4)
求取功率方程的雅克比矩阵
i≠j时，
∂ΔPi∂θj=-ViVjGijsinθij-Bijcosθij (2.5)
∂ΔPi∂VjVj=-ViVjGijcosθij+Bijsinθij (2.6)
∂ΔQi∂θj=ViVjGijcosθij+Bijsinθij (2.7)
∂ΔQi∂VjVj=-ViVjGijsinθij-Bijcosθij (2.8)
i=j时，
∂ΔPi∂θi=Vi2Bii+Qi (2.9)
∂ΔPi∂ViVi=-Vi2Gii-Pi (2.10)
∂ΔQi∂θi=Vi2Gii-Pi (2.11)
∂ΔQi∂ViVi=Vi2Bii-Qi (2.12)
牛拉法求解潮流的迭代过程如下

图2.2 牛拉法潮流迭代流程图
三、MATLAB编程实例
实例采用王锡凡主编的《现代电力系统分析》p326-327的例题。

图3.1 某小型电力网络图

例题中，1号发电机为平衡节点，2号和3号发电机为PV节点，其余为PQ节点。为了便于编程，可将PV节点和平衡节点放在节点描述图的最后，故要重新对节点进行编号。
。
图3.2 修改节点编号后的电力网络图

支路表如下：
%首节点 末节点 电阻 电抗 容纳之半 变比
1 2 0.010 0.085 0.088 1
1 3 0.017 0.092 0.079 1
2 4 0.032 0.161 0.153 1
3 6 0.039 0.170 0.179 1
4 5 0.0085 0.072 0.0745 1
5 6 0.0119 0.1008 0.1045 1
9 1 0 0.0576 0 1
7 4 0 0.0625 0 1
8 6 0 0.0586 0 1

节点表如下：
%节点 电压幅值 电压相角 节点有功 节点无功 节点类型 /“1”为PQ节点，“2”为PV节点，“3”为平衡节点
1 1.0 0 0 0 1
2 1.0 0 -1.2500 -0.5 1
3 1.0 0 -0.9 -0.3 1
4 1.0 0 0 0 1
5 1.0 0 -1 -0.35 1
6 1.0 0 0 0 1
7 1.025 0 1.6300 0 2
8 1.025 0 0.8500 0 2
9 1.040 0 0 0 3

程序清单：取ε=10-5
a=load('zhilu.txt');
b=load('jiedian.txt');
N=size(a,1);%支路数
Nbus=size(b,1);%节点数
pq=0;
for k=1:Nbus %PQ节点个数
if b(k,6)==1
pq=pq+1;
end
end
Y=zeros(Nbus);%节点导纳矩阵
for k=1:N
t1=a(k,1);t2=a(k,2);r=a(k,3);x=a(k,4);ban=a(k,5);K=a(k,6);
Y(t1,t1)=Y(t1,t1)+1/(r+j*x)+j*ban;
Y(t1,t2)=Y(t1,t2)-1/(K*(r+j*x));
Y(t2,t1)=Y(t2,t1)-1/(K*(r+j*x));
Y(t2,t2)=Y(t2,t2)+1/(K*K*(r+j*x))+j*ban;
end
G=real(Y);B=imag(Y);

precision=1;
t=0；%存储迭代次数
while precision>0.00001
P=zeros(Nbus,1);%存储所有节点的有功
Q=zeros(Nbus,1);%存储所有节点的无功
for m=1:Nbus %求pv、pq和平衡节点的有功
for n=1:Nbus
P(m,1)=P(m,1)+b(m,2)*b(n,2)*(G(m,n)*cos(b(m,3)-b(n,3))+B(m,n)*sin(b(m,3)-b(n,3)));
end
end
for m=1:Nbus %求pq、pv和平衡节点的无功
for n=1:Nbus
Q(m,1)=Q(m,1)+b(m,2)*b(n,2)*(G(m,n)*sin(b(m,3)-b(n,3))-B(m,n)*cos(b(m,3)-b(n,3)));
end
end
deltp=b(1:Nbus-1,4)-P(1:Nbus-1,1);%pq和pv节点有功差
deltq=b(1:pq,5)-Q(1:pq,1);%pq节点无功差
deltPQ=[deltp;deltq];
H=zeros(Nbus-1);
for m=1:Nbus-1 %求H矩阵
for n=1:Nbus-1
if m~=n
H(m,n)=-b(m,2)*b(n,2)*(G(m,n)*sin(b(m,3)-b(n,3))-B(m,n)*cos(b(m,3)-b(n,3)));
else
H(m,m)=b(m,2)*b(m,2)*B(m,m)+Q(m,1);
end
end
end
N=zeros(Nbus-1,pq);
for m=1:Nbus-1 %求N矩阵
for n=1:pq
if m~=n
N(m,n)=-b(m,2)*b(n,2)*(G(m,n)*cos(b(m,3)-b(n,3))+B(m,n)*sin(b(m,3)-b(n,3)));
else
N(m,m)=-b(m,2)*b(m,2)*G(m,m)-P(m,1);
end
end
end
J=zeros(pq,Nbus-1);
for m=1:pq %求J矩阵
for n=1:Nbus-1
if m~=n
J(m,n)=b(m,2)*b(n,2)*(G(m,n)*cos(b(m,3)-b(n,3))+B(m,n)*sin(b(m,3)-b(n,3)));
else
J(m,m)=b(m,2)*b(m,2)*G(m,m)-P(m,1);
end
end
end
L=zeros(pq,pq);
for m=1:pq %求L矩阵
for n=1:pq
if m~=n
L(m,n)=-b(m,2)*b(n,2)*(G(m,n)*sin(b(m,3)-b(n,3))-B(m,n)*cos(b(m,3)-b(n,3)));
else
L(m,m)=b(m,2)*b(m,2)*B(m,m)-Q(m,1);
end
end
end
Jacobi=[H N;J L];%雅克比矩阵
Correction=-Jacobi\deltPQ;%计算电压相角和幅值的修正量
for m=1:Nbus-1
b(m,3)=b(m,3)+Correction(m);
end
for m=1:pq
b(m,2)=b(m,2)+Correction(Nbus-1+m)*b(m,2);
end
precision=max(abs(deltPQ));
t=t+1
b
end
b(Nbus,4)=P(Nbus,1);
b(Nbus,5)=Q(Nbus,1);
for m=pq+1:Nbus-1
b(m,5)=Q(m,1);
end
b
运行结果：

四、结果分析与结论
运行结果与图3.2对应，例题给定的系统潮流与图3.1对应，对比之后发现结果完全一致，验证了牛顿-拉夫逊潮流计算的正确性。
程序迭代四次就完成，可见牛顿-拉夫逊算法的收敛性好，迭代次数少。但每一次迭代都要更新雅克比矩阵，而且PQ方程联立求解，方程维数很高。