电力系统课程设计极坐标表示的牛顿拉夫逊法潮流计算程序设计.docx

目录
1.任务书……………………………………（2）
2.模型简介及等值电路……………………（3）
3.设计原理…………………………………（5）
4.修正方程的建立…………………………（7）
5.程序流程图及MALAB程序编写……… （9）
6.结果分析…………………………………（16）
7.设计总结…………………………………（20）
8.参考文献…………………………………（20）

一.任务书
课程设计任务书
题 目
极坐标表示的牛顿拉夫逊法潮流计算程序设计
学生姓名

学号

专业班级
电气1001
设
计
内
容
与
要
求
1. 设计要求
掌握MATLAB语言编程方法；理解和掌握运用计算机进行潮流计算的基本算法原理；针对某一具体电网，进行潮流计算程序设计。
其目的在于加深学生对电力系统稳态分析中课程中基本概念和计算方法的理解，培养学生运用所学知识分析和解决问题的能力。
2. 内容
1）学习并掌握MATLAB语言。
2）掌握变压器非标准变比概念及非标准变比变压器的等值电路。掌握节点导纳矩阵的概念及导纳矩阵的形成和修改方法。
3）掌握电力系统功率方程、变量和节点分类。
4）掌握利用极坐标表示的牛-拉法进行潮流计算的方法和步骤。
5）选择一个某一具体电网，编制程序流程框图。
6）利用MATLAB语言编写该模型的潮流计算程序，并上机调试程序，对计算结果进行分析。
7）整理课程设计论文。
起止时间
2013年7月4日至2013 年 7月10日
指导教师签名
年 月 日
系（教研室）主任签名
年 月 日
学生签名
年 月 日

二.模型简介及解题思路
2.1课题模型及等值电路：
模型3
电力网络接线如下图所示，各支路阻抗标幺值参数如下：Z12=0.02+j0.06，Z13=0.08+j0.24， Z23=0.06+j0.18， Z24=0.06+j0.12， Z25=0.04+j0.12， Z34=0.01+j0.03， Z45=0.08+j0.24， k=1.1。该系统中，节点1为平衡节点，保持为定值；节点2、3、4都是PQ节点，节点5为PV节点，给定的注入功率分别为:
， ，，，。各节点电压（初值）标幺值参数如下：

节点

1

2

3

4

5

Ui（0）=ei（0）+jfi（0）

1.06+j0.0

1.0+j0.0

1.0+j0.0

1.0+j0.0

1.1+j0.0

计算该系统的潮流分布。计算精度要求各节点电压修正量不大于10-5。

图2.1电路图

节点1是平衡节点，节点2、3、4是PQ节点，节点5是PV节点。由题可得等值电路模型中各节点之间的导纳：y12=5.000-j15.000，y13=1.2500-j3.7500，y22=0.2750-j0.8250，y23=1.667-j5.000，y24=3.333-j6.667，y25=2.7500-j8.2500，y34=10.0000-j30.0000，y55=-0.25+j0.75

在图2-2中，将图2-1中的编号重新编排，节点2、3、4、5、1替换为1、2、3、4、5。则各节点之间的导纳变为y12=1.667-j5，y13=3.333-j6.667，y14=2.75-j8.25，y15=5-j15，y52=1.25-j3.75，y23=10-j30，y34=1.25-j3.75，y11=0.275-j0.825，y44=-0.25+j0.75。

2.2 求解思路
此电力系统是一个5节点，4支路的电力网络。其中包含3个PQ节点，一个PV节点，和一个平衡节点。综合比较牛顿拉夫逊法（直角坐标、极坐标）、PQ分解法等多种求解方法的特点，最后确定采用牛顿拉夫逊法（极坐标）。因为此方法所需解的方程组最少。
三.设计原理
本题采用了题目要求的牛顿－拉夫逊潮流计算的方法。牛顿法是数学中解决非线性方程式的典型方法，有较好的收敛性。解决电力系统潮流计算问题是以导纳距阵为基础的，因此，只要在迭代过程中尽可能保持方程式系数距阵的稀疏性，就可以大大提高牛顿法潮流程序的放率。
3.1潮流计算的定解条件
题中所给图表示一个五节点的简单电力系统，n个节点电力系统的潮流方程的一般形式是

或

按变量的不同，一般将节点分为三种类型。
1 节点
这类节点的有功功率和无功功率是给定的，节点是待求量。通常变电所都是这一类型节点。由于没有发电设备，故其发电功率为零。有些情况下，系统中某些发电厂输出的功率在一段时间内是固定时，该发电厂母线也作为节点。因此，电力系统中绝大多数节点属于这一类型。
2 节点
这类节点有功功率P和电压幅值V是给定的，节点的无功功率Q和电压的相位是待求量。这类节点必须有足够的可调无功容量，用以维持给定的电压幅值，因此又称为电压控制节点。一般选择有一定无功储备的发电厂和具有可调无功电源设备的变电所作为节点。
3 平衡节点
在潮流分布算出以前，网络中的功率损耗是未知的，因此，网络中至少有一个节点的有功功率P不能给定，这个节点承担了系统系统的有功功率平衡，故称之为平衡节点。

3.2潮流计算的约束条件
1所有节点电压必须满足

2所有电源节点的有功功率和无功功率必须满足

3某些节点之间的电压应满足

3.3牛顿-拉夫逊的基本原理
   设欲求解的非线性代数方程为

设方程的真实解为，则必有。用牛顿-拉夫逊法求方程真实解的步骤如下： 
 首先选取余割合适的初始估值作为方程的解，若恰巧有，则方程的真实解即为若，则做下一步。取 则
其中为初始估值的增量，即。设函数具有任意阶导数   
若所取的足够小，则含的项及其余的一切高阶项均可略去，并使其等于零，即：　　　        

故得                          

可见，只要≠0，即可根据上式求出第一次的修正估值，若恰巧有 =0，则方程的真实解即为。若，则用上述方法由再确定第二次的修正估值。如此反复叠代下去，直到求得真实解为止。设第K次的估值为第（K+1）次的修正估值，则有

迭代过程的收敛数据为
或
其中，，为预先给定的小正数。

四.修正方程的建立
（1） 网络中共有n个节点，编号为1,2,3…，n，其中包含一个平衡节点，编号为s。
（2） 网络中有（m-1）个PQ节点，编号为1,2,3，…，m，其中包含编号为s的平衡节点。
（3） 网络中有（n-m）个PV节点，编号为m+1,m+2,…,n。
极坐标表示的牛拉法修正方程
= （4-44）
式中留出了（n-m）行空格和（n-m）列空列。式中的有功、无功功率不平衡量、分别由式

得出
（4-45a）
（4-45b）
而式中雅克比矩阵的各元素则分别为
；；; （4-46）

式（4-44）中将改为只是为使公式（4-46）中个偏导数的表示形式上更相似，为求取这些偏导数，可将、分别展开如下
(4-47a、b)
计及
（4-48）
时，由于对特定的j，只有该特定节点的，从而特定的是变量，由式（4-46）到式（4-48）可得
(4-49a)
相似的，由于对特定的j，只有该特定节点的是变量，可得
(4-49b)
j=i时，由于是变量，所有都是变量，可得
（4-49c）
相似的，由于是变量，可得
（4-49d）

五.程序流程图及MATLAB程序编写
5.1程序流程图（见下页）

5.2 MATLAB程序编写
MATLAB语言来编写N-L极坐标潮流计算程序,程序的清单如下:
％程序中的各个变量的说明：
m、n 代表系统中的节点总数，把平衡节点标为最大号；
N1 代表P-Q节点的总数；
pp 代表有功功率的不平衡量ΔPi
qq代表无功功率的不平衡量ΔQi
u 代表各节点电压的大小
delt代表各节点电压的相角
S 代表线路的功率
g(1,1)=0.275;b(1,1)=0.825;
g(1,2)=1.667;b(1,2)=-5;
g(1,3)=3.333;b(1,3)=-6.667;
g(1,4)=2.75;b(1,4)=-8.25;
g(1,5)=5;b(1,5)=-15;
g(2,1)=1.667;b(2,1)=-5;
g(2,2)=-0;b(2,2)=0;
g(2,3)=10;b(2,3)=-30;
g(2,4)=0;b(2,4)=0;
g(2,5)=1.25;b(2,5)=-3.75;
g(3,1)=3.33;b(3,1)=-6.667;
g(3,2)=10;b(3,2)=-30;
g(3,3)=0;b(3,3)=0;
g(3,4)=1.25;b(3,4)=-3.75;
g(3,5)=0;b(3,5)=0;
g(4,1)=2.75;b(4,1)=-8.25;
g(4,2)=0;b(4,2)=0;
g(4,3)=1.24;b(4,3)=-3.75;
g(4,4)=-0.25;b(4,4)=0.75;
g(4,5)=0;b(4,5)=0;
g(5,1)=5;b(5,1)=-15;
g(5,2)=1.25;b(5,2)=-3.75;
g(5,3)=0;b(5,3)=0;
g(5,4)=0;b(5,4)=0;
g(5,5)=0;b(5,5)=0;
d(1,4)=0.275-0.825*j;
d(4,1)=-0.25+0.75*j
d(1,5)=0;d(5,1)=0;
%求取节点导纳矩阵。
for m=1:5
for n=1:5
if m==n
G(m,m)=g(m,1)+g(m,2)+g(m,3)+g(m,4)+g(m,5);
B(m,m)=b(m,1)+b(m,2)+b(m,3)+b(m,4)+b(m,5);
else
G(m,n)=-g(m,n);
B(m,n)=-b(m,n);
end
end
end
Y=G+j*B
%//////////////////////////////////////////////////////////////
%//设定节点起始节点电压，并给出已知功率值。
%//下面将题中节点 2、3、4、5、1分别替换为节点1、2、3、4、5即节点4为PV节点，节点5为平衡节点。
delt(1)=0;delt(2)=0;delt(3)=0;delt(4)=0;
u(1)=1.0;u(2)=1.0;u(3)=1.0;
p(1)=0.2;q(1)=0.2;
p(2)=-0.45;q(2)=-0.15;p(3)=-0.4;q(3)=-0.05;p(4)=-0.5;q(4)=0.0;
%置迭代次数k=0，并计算节点功率的不平衡量 ， 。
k=0;precision=1;
k,delt,u
N1=4;
while precision>0.00001
u(4)=1.1; delt(5)=0;u(5)=1.06;
for m=1:N1
for n=1:N1+1
pt(n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));
qt(n)=u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));
end
pp(m)=p(m)-sum(pt);qq(m)=q(m)-sum(qt);
end
pp,qq
%///////////////////////////////////////////////////////
%求取雅可比矩阵元素 （m=n 时）。
for m=1:N1
for n=1:N1+1
h0(n)=u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));
n0(n)=-u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));
j0(n)=-u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));
l0(n)=-u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));
end
H(m,m)=sum(h0)-u(m)^2*(G(m,m)*sin(delt(m)-delt(m))-B(m,m)*cos(delt(m)-delt(m)));
N(m,m)=sum(n0)-2*u(m)^2*G(m,m)+u(m)^2*(G(m,m)*cos(delt(m)-delt(m)) +B(m,m)*sin(delt(m)-delt(m)));
J(m,m)=sum(j0)+u(m)^2*(G(m,m)*cos(delt(m)-delt(m))+B(m,m)*sin(delt(m)-delt(m)));
L(m,m)=sum(l0)+2*u(m)^2*B(m,m)+u(m)^2*(G(m,m)*sin(delt(m)-delt(m)) -B(m,m)*cos(delt(m)-delt(m)));
end

for m=1:N1-1
JJ(2*m-1,2*m-1)=H(m,m);JJ(2*m-1,2*m)=N(m,m);
JJ(2*m,2*m-1)=J(m,m);JJ(2*m,2*m)=L(m,m);

end
for m=N1:N1
JJ(2*m-1,2*m-1)=H(m,m);
end

%///////////////////////////////////////////////////////
%求取雅可比矩阵元素 （m不等于n时）
for m=1:N1
for n=1:N1
if m==n
else
H(m,n)=-u(m)*u(n)*(G(m,n)*sin(delt(m)-delt(n))-B(m,n)*cos(delt(m)-delt(n)));
J(m,n)=u(m)*u(n)*(G(m,n)*cos(delt(m)-delt(n))+B(m,n)*sin(delt(m)-delt(n)));
N(m,n)=-J(m,n);L(m,n)=H(m,n);
end
end
end
for m=1:N1-1 %求前六行六列元素(非对角)
for n=1:N1-1
if m==n
else
JJ(2*m-1,2*n-1)=H(m,n);JJ(2*m-1,2*n)=N(m,n);
JJ(2*m,2*n-1)=J(m,n);JJ(2*m,2*n)=L(m,n);
end
end
end
for m=N1
for n=1:N1-1 %求取第七行的元素
JJ(2*m-1,2*n-1)=H(m,n);JJ(2*m-1,2*n)=N(m,n);
end
end
for n=N1
for m=1:N1-1 %求的第七列元素
JJ(2*m-1,2*n-1)=H(m,n);JJ(2*m,2*n-1)=J(m,n);
end
end
%///////////////////////////////////////////////////////
%解修正方程式，由 ， 和 计算电压修正量 和 。
for m=1:N1-1
PP(2*m-1)=pp(m);PP(2*m)=qq(m);
end
for m=N1
PP(2*m-1)=pp(m);
end
uu=-inv(JJ)*PP';precision=max(abs(uu));uu
%///////////////////////////////////////////////////////
%若结果不收敛，执行下列语句。
for n=1:N1-1
delt(n)=delt(n)+uu(2*n-1);
u(n)=u(n)+uu(2*n);
end
for n=N1
delt(n)=delt(n)+uu(2*n-1);

end

k=k+1;
k,delt,u;
end
%///////////////////////////////////////////////////////
%若结果收敛，计算各节点电压，平衡节点功率，PV节点功率和线路功率 。
for n=1:N1+1
U(n)=u(n)*(cos(delt(n))+j*sin(delt(n)));
end
for m=1:N1+1
I(m)=Y(5,m)*U(m);
end
S5=U(5)*sum(conj(I)); %平衡节点功率
%PV节点功率
for n=1:N1+1
q4(n)=u(4)*u(n)*(G(4,n)*sin(delt(4)-delt(n))-B(4,n)*cos(delt(4)-delt(n)));
end
Q4=sum(q4)
%线路功率 a
for m=1:N1+1
for n=1:N1+1

S(m,n)=U(m)*(conj(U(m))*conj(d(m,n))+(conj(U(m))-conj(U(n)))*conj(-Y(m,n)))

end
end
%///////////////////////////////////////////////////////
%显示运行结果（至结束）。
Y
JJ
S
B
pp
qq
uu
U
k
Q4
S5
六.计算结果及分析
6.1、计算结果：
求得导纳矩阵：
Y =

13.0250 -34.0920i -1.6670 + 5.0000i -3.3330 + 6.6670i -2.7500 + 8.2500i -5.0000 +15.0000i
-1.6670 + 5.0000i 12.9170 -38.7500i -10.0000 +30.0000i 0 -1.2500 + 3.7500i
-3.3300 + 6.6670i -10.0000 +30.0000i 14.5800 -40.4170i -1.2500 + 3.7500i 0
-2.7500 + 8.2500i 0 -1.2400 + 3.7500i 3.7400 -11.2500i 0
-5.0000 +15.0000i -1.2500 + 3.7500i 0 0 6.2500 -18.7500i
经迭代后的电压和相角：
U =

1.0966 - 0.0677i 1.0624 - 0.0881i 1.0663 - 0.0920i 1.0959 - 0.0945i 1.0600

k =

6
平衡节点功率：
Q4 =

-0.6895

S5 =

1.2297 - 1.0678i
雅可比矩阵：
JJ =

-40.9565 -15.9241 5.8963 1.8286 7.9334 3.7268 10.0491
15.5240 -41.3566 -1.8286 5.8963 -3.7268 7.9334 -3.0796
5.8139 2.0756 -44.1859 -14.2290 34.2658 11.2955 0
-2.0756 5.8139 15.1290 -43.8859 -11.2955 34.2658 0
7.7421 4.1061 34.1899 11.5231 -46.3468 -16.3011 4.4148
-4.1061 7.7421 -11.5231 34.1899 17.1011 -46.2468 -1.4718
9.8870 3.5657 0 0 4.4149 1.4597 -14.3020
功率分布：
S =

0 0.1839 + 0.1398i 0.2968 + 0.1152i 0.5722 + 0.9064i -0.8529 + 1.0305i
-0.1812 - 0.1319i 0 0.0687 - 0.1735i 0 -0.3375 + 0.1554i
-0.2917 - 0.1052i -0.0683 + 0.1744i 0 -0.0400 - 0.1193i 0
-0.5407 - 0.8120i 0 0.0407 + 0.1226i 0 0
0.8825 - 0.9416i 0.3472 - 0.1262i 0 0 0

B =

-34.0920 5.0000 6.6670 8.2500 15.0000
5.0000 -38.7500 30.0000 0 3.7500
6.6670 30.0000 -40.4170 3.7500 0
8.2500 0 3.7500 -11.2500 0
15.0000 3.7500 0 0 -18.7500

pp =

1.0e-004 *

-0.1416 0.0200 0.0231 0.0377

qq =

-0.0000 0.0000 0.0000 0.6895

uu =

1.0e-006 *

-0.0000
-0.9640
-0.0000
-0.1828
-0.0000
-0.2303
-0.0001
迭代过程中各节点功率的不平衡量

K
△P1
△Q1
△P2
△Q2
△P3
△Q3
△P4
△Q4

0

0.5000

2.7500

-3.750

0.0750

-0.2750

0.3250

-0.6375

-0.4125

1

-0.1013

-0.3711

0.0261

-0.0269

0.0501

-0.0454

0.0208

0.8565

2

-0.1620

-0.0428

0.0013

0.0018

0.0017

-0.0013

0.0047

0.7059

3

-0.0016

-0.0042

-0.0002

0.0004

0.0002

0.0001

0.0005

0.6907

4

-0.152.3
*e-3

-0.0004
*e-3

0.0199
*e-3

0.0001
*e-3

0.0209
*e-3

0.0000
*e-3

0.0433
*e-3

0.6894
*e-3

5

-0.1416
*e-4

-0.0000
*e-4

0.0200
*e-4

0.0000
*e-4

0.0023
*e-4

0.0000
*e-4

0.0378
*e-4

0.6893
*e-4

迭代过程中的各节点电压的修正量

k
△δ1
△U1
△δ2
△U2
△δ3
△U3
△δ4

0

-0.0648

0.1115

-0.0882

0.0757

-0.0918

0.0801

-0.0876

1

0.0031

-0.0114

0.0054

-0.0089

0.0056

-0.0090

0.0016

2

-0.0000

-0.0013

0.0000

0.0007

-0.0000

-0.0008

-0.0001

3

-0.0004
*e-3

-0.1190
*e-3

-0.0004
*e-3

-0.0475
*e-3

-0.0005
*e-3

-0.5320
*e-3

-0.0010
*e-3

4

-0.0000
*e-4

-0.1073
*e-3

-0.0000
*e-4

-0.2950
*e-4

-0.0000
*e-4

-0.0351
*e-4

-0.0001
*e-4

5

-0.0000
*e-6

-0.9640
*e-6

-0.0000
*e-6

-0.1828
*e-6

-0.0000
*e-6

-0.2303
*e-6

-0.0001
*e-6

迭代过程中各节点电压

K
δ1
U1
δ2
U2
δ3
U3
δ4
U4
δ5
U5
0
0
1.0000
0
1.0000
0
1.0000
0
1.1000
0
1.0600
1
-0.0648
1.1115
-0.0882
1.0757
-0.0918
1.0801
-0.0876
1.1000
0
1.0600
2
-0.0617
1.1002
-0.0828
1.0668
-0.0861
1.0711
-0.0861
1.1000
0
1.0600
3
-0.0617
1.0989
-0.0828
1.0661
-0.0861
1.0703
-0.0860
1.1000
0
1.0600
4
-0.0617
1.0987
-0.0828
1.0660
-0.0861
1.0703
-0.0861
1.1000
0
1.0600
5
-0.0617
1.0987
-0.0828
1.0660
-0.0861
1.0703
-0.0861
1.1000
0
1.0600
6
-0.0617
1.0987
-0.0828
1.0660
-0.0861
1.0703
-0.0861
1.1000
0
1.0600

各线路的功率
j
i
1
2
3
4
5
1
0
0.1839 + 0.1398i
0.2968 + 0.1152i
0.5722 + 0.9064i
-0.8529 + 1.0305i
2
-0.1812 - 0.1319i
0
0.0687 - 0.1735i
0
-0.3375 + 0.1554i
3
-0.2917 - 0.1052i
-0.0683 + 0.1744i
0
-0.0400 - 0.1193i
0
4
-0.5407 - 0.8120i
0
0.0407 + 0.1226i
0
0
5
0.8825 - 0.9416i
0.3472 - 0.1262i
0
0
0

6.2结果分析：
由上述程序运行结果可知：该电力系统模型在极坐标表示的牛顿拉夫逊法潮流计算方式下经过6次的迭代后，其精确度达到了要求的各节点电压的修正量不大于10-5。
第三部分 设计总结及参考文献
七.设计总结
此次课程设计首先让我明白了要使电力系统运行的稳定，必须经过精密的设计和计算。在进行课题设计的过程中，加深了我对潮流计算的认识，尤其是对牛顿拉夫逊潮流计算的求解思路有了比较透彻的理解。同时由于求解过程中用到求节点导钠矩阵，求矩阵的逆等等，又使我对以前所学的知识有了一次很好的温习。同时也看到了研究性学习的效果，从研究中去学习，理论结合实际，将理论运用到实际，同时在实践中发现问题，然后解决问题，
而且在此次课程设计中，我发现了自己的基础知识有很多的不足。这些基础的缺乏给我的设计计划造成了不小的障碍。在这个过程中，我明白了，只要用心去做，认真去做，持之以恒，就会有新的发现，有意外的收获。

八.参考文献
《电力系统分析》 何仰赞 华中理工大学出版社
《电力系统稳态分析》 陈珩 中国电力出版社
《MATLAB程序设计与应用》 张智星 清华大学出版社
《MATLAB语言实用教程》 马莉 清华大学出版社