电力系统潮流及短路电流计算程序.doc

班级：
姓名：
学号：
一、 作业要求
编写程序计算图1所示算例系统的潮流及三相短路电流。
潮流计算：方法不限，计算系统的节点电压和相角。
短路电流：4号母线发生金属性三相短路时(zf=0）,分别按照精确算法和近似算法计算短路电流、系统中各节点电压以及网络中各支路的电流分布，并对两种情况下的计算结果进行比较.
二、 电路图及参数
图1 3机9节点系统
表1 9节点系统支路参数
支路
R（p.u。）
X（p。u。)
B/2（TK） （p.u。）
1~4
0
0.0576
1.0
2~7
0
0.0625
1。0
3～9
0
0.0586
1.0
4～5
0.01
0.085
0。088
4～6
0。017
0.092
0.079
5～7
0.032
0.161
0.153
6~9
0。039
0。17
0.179
7~8
0.0085
0.072
0.0745
8~9
0。0119
0。1008
0。1045
表2 9节点系统发电机参数
发电机编号
节点类型
PG（p.u。)
VG（p。u。)
（p。u.）
（p.u.)
1
平衡

1.04
0.3
1.137
2
PV
1.63
1.025
0.3
1。211
3
PV
0。85
1。025
0。3
1。043
表3 9节点系统负荷参数
节点编号
节点类型
Pi(p。u。）
Qi(p.u.）
4
PQ
0
0
5
PQ
1。25
0。5
6
PQ
0.9
0.3
7
PQ
0
0
8
PQ
1
0.35
9
PQ
0
0
三、 计算步骤
(1) 进行系统正常运行状态的潮流计算，求得
(2) 形成不含发电机和负荷的节点导纳矩阵YN；
(3) 将发电机表示为电流源（)和导纳(）的并联组合;节点负荷用恒阻抗的接地支路表示;形成包括所有发电机支路和负荷支路的节点导纳矩阵Y,即在YN中的发电机节点和负荷节点的自导纳上分别增加发电机导纳和负荷导纳（）；
(4) 利用，计算节点阻抗矩阵,从而得到阻抗矩阵中的第f列；
(5) 利用公式(6-7）或(6—10)计算短路电流；
(6) 利用公式（6-8)或(6—11）计算系统中各节点电压；
(7) 利用公式（6—9）计算变压器支路的电流;对输电线路利用П型等值电路计算支路电流.
四、 计算结果
节点导纳矩阵Yn：
Columns 1 through 5
0 —17.3611i 0 0 0 +17.3611i 0
0 0 -16。0000i 0 0 0
0 0 0 —17.0648i 0 0
0 +17。3611i 0 0 3。3074 -39。3089i —1.3652 +11.6041i
0 0 0 -1。3652 +11.6041i 2.5528 -17。3382i
0 0 0 -1。9422 +10.5107i 0
0 0 +16。0000i 0 0 -1.1876 + 5.9751i
0 0 0 0 0
0 0 0 +17.0648i 0 0
Columns 6 through 9
0 0 0 0
0 0 +16.0000i 0 0
0 0 0 0 +17.0648i
-1.9422 +10.5107i 0 0 0
0 -1.1876 + 5。9751i 0 0
3.2242 —15.8409i 0 0 —1.2820 + 5。5882i
0 2.8047 —35.4456i -1.6171 +13。6980i 0
0 -1.6171 +13.6980i 2。7722 —23.3032i —1。1551 + 9.7843i
—1.2820 + 5.5882i 0 —1.1551 + 9。7843i 2.4371 -32。1539i
电压幅值：
1。0400 1.0250 1。0250 1。0258 0。9956 1。0127 1.0258 1。0159 1.0324
电压相角：
0 0。1620 0。0814 —0。0387 —0.0696 —0.0644 0。0649 0。0127 0.0343
节点有功：
0.7164 1。6300 0。8500 0.0000 —1.2500 -0。9000 —0。0000 -1.0000 -0。0000
节点无功：
0.2705 0.0665 -0.1086 0。0000 -0.5000 -0。3000 —0.0000 —0。3500 -0.0000
修正后的节点导纳矩阵Y：
Columns 1 through 5
0 —20.6944i 0 0 0 +17.3611i 0
0 0 -19。3333i 0 0 0
0 0 0 —20。3982i 0 0
0 +17.3611i 0 0 3.3074 —39.3089i -1.3652 +11。6041i
0 0 0 -1。3652 +11.6041i 3.8716 -17.6627i
0 0 0 -1.9422 +10。5107i 0
0 0 +16。0000i 0 0 —1。1876 + 5.9751i
0 0 0 0 0
0 0 0 +17。0648i 0 0
Columns 6 through 9
0 0 0 0
0 0 +16。0000i 0 0
0 0 0 0 +17。0648i
-1.9422 +10.5107i 0 0 0
0 -1。1876 + 5.9751i 0 0
4。1321 -16.0184i 0 0 —1。2820 + 5。5882i
0 2。8047 -35。4456i -1.6171 +13.6980i 0
0 -1.6171 +13.6980i 3。7323 -23.6669i -1。1551 + 9。7843i
—1。2820 + 5。5882i 0 —1.1551 + 9。7843i 2.4371 -32.1539i
节点阻抗矩阵Z的第4列：
0。0463 + 0.1252i
0。0329 + 0。0693i
0.0316 + 0。0707i
0。0552 + 0.1493i
0。0589 + 0.1204i
0.0562 + 0.1226i
0.0397 + 0。0838i
0.0416 + 0.0814i
0。0378 + 0。0845i
精确计算结果:
短路电流：
模值：
6.4459
相角：
—71。9365
节点电压模值：
0.1831 0。5687 0.5427 0.0000 0。1466 0。1506 0。4537 0。4463 0.4495
支路电流：
i j Iij
1 4 0.5779—3.1264i
2 7 1.3702-1.4433i
3 9 0.64294—1.4808i
4 5 —0.77968+1。5248i
4 6 -0.6411+1.477i
5 7 —0.89528+1.6436i
6 9 —0。73353+1。5487i
7 8 0。50734+0。10234i
8 9 0.062766+0.056451i
近似计算结果：
短路电流：
模值:
6.2838
相角：
—69.7198
节点电压模值：
0。1611 0.5214 0.5157 0.0000 0。1827 0.1675 0。4227 0.4348 0。4217
五、程序流程图
Y
N
N
Y
结 束
输出结果
计算潮流分布
输出潮流不收敛信息
k←k+1
k>Kmax?
Max|fi（x[k]）|<ɛ?
求解修正方程获得∆x[k]
x[k+1] ←x[k]+∆x[k]

计算节点导纳矩阵参数
设置节点电压初值x[0]
设置k←0及最大迭代次数Kmax
计算f（x[k]）及雅克比矩阵
开 始
输入电网及节点注入数据

六、程序及输入文件
input_data.xls 文件：
Sheet1
9节点系统支路参数
首节点编号i
尾节点编号j
R（p。u。）
X（p.u.)
B/2（p。u.）
1
4
0
0。0576
0
2
7
0
0.0625
0
3
9
0
0.0586
0
4
5
0。01
0.085
0.088
4
6
0。017
0.092
0.079
5
7
0.032
0.161
0。153
6
9
0。039
0.17
0。179
7
8
0。0085
0。072
0.0745
8
9
0.0119
0。1008
0.1045
Sheet2
9节点系统发电机及负荷参数
节点编号
节点类型
P（p。u。）
Q(p。u。）
U（p。u.）
Xd（p.u.）
E（p。u.）
1
平衡

1.04
0.3
1.137
2
PV
1.63

1。025
0。3
1。211
3
PV
0。85

1.025
0.3
1。043
4
PQ
0
0

5
PQ
-1。25
—0。5

6
PQ
—0.9
—0.3

7
PQ
0
0

8
PQ
—1
-0.35

9
PQ
0
0

powerflow_cal.m 文件：
l=9；%支路数
n=9；%节点数
m=6；%PQ节点数
Yn=zeros（n）；％初始化节点导纳矩阵Y
DATA1=xlsread（'input_data.xls',1）；
％计算节点导纳矩阵Y
for k=1：l
i=DATA1(k，1);
j=DATA1（k，2)；
R=DATA1（k,3）；
X=DATA1（k，4）；
B2=DATA1（k，5);
Yn(i,i）=Yn（i,i)+1i＊B2+1/（R+1i*X）;
Yn(j,j）=Yn(j，j)+1i*B2+1/（R+1i＊X）；
Yn（i，j）=Yn(i,j)—1/（R+1i＊X);
Yn(j，i）=Yn（j，i)—1/（R+1i＊X）；
end
disp('节点导纳矩阵Yn：'）；
disp（Yn);
G=real（Yn);
B=imag(Yn);
DATA2=xlsread（’input_data。xls’,2）;
P=zeros(1，n);
Q=zeros(1，n)；
U=ones（1，n)；
P(2:n）=DATA2(2：n，3）；
Q（4：n）=DATA2（4：n,4)；
U(1：3）=DATA2（1:3，5)；
%设置节点电压初值
e(1)=DATA2(1,5）;
e(2:n)=1。0;
f(1：n)=0.0；
％设置迭代次数
t=0;
tmax=10;
while t〈=tmax
%计算f(x）
a(1：n)=0。0;
c(1：n)=0.0；
for i=2：n
for j=1：n
a(i）=a（i）+G（i，j)＊e(j）—B(i，j)*f（j）;
c（i)=c(i)+G(i，j）＊f（j)+B（i,j)＊e(j）;
end
end
for i=2:n
deltaP（i）=P（i)-e（i）＊a（i）—f(i）*c(i）;
end
for j=4：n
deltaQ(j)=Q（j）-f(j）＊a(j）+e（j）＊c(j);
end
for k=2：3
deltaU2（k）=U(k)＊U（k)—e(k)＊e（k）-f（k)*f(k);
end
fx=[deltaP（2:n) deltaQ(4：n） deltaU2（2:3)］';
％计算雅克比矩阵J
for i=2：n
for j=2:n
if i~=j
H(i,j）=—（G(i,j）*e(i）+B（i,j）*f（i））；
N（i，j)=B(i,j)*e（i)-G（i，j)*f(i)；
else
H（i，j)=—a（i）-（G(i,i）＊e(i）+B（i,i）*f(i))；
N(i，j）=-c（i）+(B(i，i)*e(i)-G（i，i）*f（i));
end
end
end
for i=4：n
for j=2:n
if i~=j
M(i，j）=B（i，j)*e（i)-G（i,j）＊f（i）；
L（i，j)=G(i，j)*e(i)+B(i，j）*f（i)；
else
M（i，j)=c(i）+（B(i，i）＊e(i)—G（i，i)＊f(i）)；
L(i,j)=—a（i）+（G（i,i）＊e(i）+B（i,i）*f(i));
end
end
end
for i=2：3
for j=2：n
if i～=j
R(i,j)=0；
S(i,j)=0;
else
R(i，j)=—2＊e(i);
S(i,j）=-2＊f(i）;
end
end
end
J=［H（2：n，2:n) N（2:n，2：n）;M（4：n,2：n） L（4：n，2：n）;R(2：3，2：n） S(2:3,2:n)］；
if max（abs(fx））〈0.0001
％输出结果
break;
else
%求解修正方程获得dx
dx=-J^（-1)*fx;
dx=dx’;
e（2:n)=e（2:n）+dx(1：n—1）；
f（2：n）=f(2:n)+dx(n:2*(n—1）);
t=t+1；
end
end
if t〉tmax
str='潮流计算不收敛’;
disp(str);
else
a（1：n）=0。0；
c(1:n)=0.0；
for i=1：n
for j=1：n
a(i）=a（i)+G(i，j)＊e(j）-B（i，j)*f(j）;
c(i）=c(i)+G(i，j）＊f(j)+B(i，j）＊e（j）；
end
end
for i=1:n
U（i)=e(i)+1i＊f(i)；
amp(i）=abs（U（i)）;
arg(i）=angle（U（i））;
P(i)=e(i)＊a（i）+f（i）*c(i）;
Q（i)=f(i）*a(i）-e(i）*c（i);
end
disp（'电压幅值：’）；
disp（amp)；
disp(’电压相角：')；
disp(arg）;
disp（'节点有功：’）;
disp(P）；
disp（’节点无功：'）;
disp（Q）;
end
％计算短路电流
f=4；
zf=0.0;
％修正节点导纳矩阵
Xd=DATA2(1：3,6）；
E=DATA2(1：3,7）;
for i=1：3
Ii(i）=E（i）/(1i*Xd(i)）；
end
Y=Yn；
for i=1：3
Y(i，i)=Y(i,i)+1/（1i＊Xd（i));
end
for j=4:n
Y(j,j）=Y（j,j)+（—P(j)+1i*Q（j））/(U（j)＊U(j));
end
disp(’修正后的节点导纳矩阵Y:'）；
disp（Y）;
Z=Y^(—1）;
disp（’节点阻抗矩阵Z的第4列：’)；
disp(Z(:，4）);
%精确计算
disp('精确计算结果:’);
U0=U；
If=U0（f)/(Z(f，f）+zf)；
amp=abs(If)；
arg=atand（imag(If)/real（If））;
disp（'短路电流：’);
disp(’模值:’)；
disp（amp）;
disp（'相角：’）;
disp（arg）;
for i=1:n
U(i）=U0(i）-Z（i,f）*If；
amp=abs(U);
end
disp（’节点电压模值：’);
disp（amp);
disp（’支路电流： '）；
str=[’i ' ’j ' ' Iij'］;
disp（str）;
for k=1：l
i=DATA1(k,1);
j=DATA1（k，2)；
r=DATA1(k,3);
x=DATA1（k，4);
z=r+1i＊x;
I=(U(i）—U（j))/z;
str=[num2str(i) ' ' num2str(j） ’ ’ num2str(I）］；
disp（str）；
end
%近似计算
disp（’近似计算结果:');
U0(1：n）=1.0;
If=U0（f)/（Z(f,f)+zf）;
amp=abs(If)；
arg=atand（imag(If）/real(If））；
disp（'短路电流:’）；
disp('模值：’）；
disp（amp）;
disp(’相角：')；
disp(arg);
for i=1：n
U(i）=U0(i）-Z（i，f)*If;
amp=abs(U）；
end
disp(’节点电压模值：')；
disp（amp）；
10