电力系统分析大作业matlab三机九节点潮流计算报告.docx
《电力系统分析大作业matlab三机九节点潮流计算报告.docx》由会员分享，可在线阅读，更多相关《电力系统分析大作业matlab三机九节点潮流计算报告.docx（17页珍藏版）》请在冰豆网上搜索。
电力系统分析大作业matlab三机九节点潮流计算报告
电力系统分析大作业
一、设计题目
本次设计题目选自课本第五章例5-8,美国西部联合电网WSCC系统的简化三机九节点系统，例题中已经给出了潮流结果，计算结果可以与之对照。
取ε=0.00001。
二、计算步骤
第一步，为了方便编程，修改节点的序号，将平衡节点放在最后。
如下图：
第二步，这样得出的系统参数如下表所示：
第三步，形成节点导纳矩阵。
第四步，设定初值：
；
，
。
第五步，计算失配功率
=0，
=-1.25，
=-0.9，
=0，
=-1，
=0，
=1.63，
=0.85；
=0.8614，
=-0.2590，
=-0.0420，
=0.6275，
=-0.1710，
=0.7101。
显然，
。
第六步，形成雅克比矩阵（阶数为14×14）
第七步，解修正方程，得到：
-0.0371，
-0.0668，
-0.0628，
0.0732，
0.0191，
0.0422，
0.1726，
0.0908；
0.0334，
0.0084，
0.0223，
0.0372，
0.0266，
0.0400。
从而
-0.0371，
-0.0668，
-0.0628，
0.0732，
0.0191，
0.0422，
0.1726，
0.0908；
1.0334，
1.0084，
1.0223，
1.0372，
1.0266，
1.0400。
然后转入下一次迭代。
经三次迭代后
。
迭代过程中失配功率的变化情况如下表：
k
0
1
2
3
Δ=P1
0
-0.0106
0.0001
0.0000000845
Δ=P2
-1.25
0.0379
0.0005
0.0000000896
Δ=P3
-0.9
0.0439
0.0005
0.0000000778
Δ=P4
0
-0.0421
-0.0012
-0.0000003421
Δ=P5
-1
0.061
0.0009
0.0000001845
Δ=P6
0
-0.0269
-0.0007
-0.0000001631
Δ=P7
1.63
-0.0579
-0.0004
-0.0000000278
Δ=P8
0.85
-0.0336
-0.0002
-0.0000000103
Δ=Q1
0.8614
-0.0501
-0.0004
-0.0000000561
Δ=Q2
-0.259
-0.0714
-0.0012
-0.0000002774
Δ=Q3
-0.042
-0.0424
-0.0006
-0.0000001236
Δ=Q4
0.6275
-0.1875
-0.0021
-0.0000003279
Δ=Q5
-0.171
-0.0241
-0.0004
-0.0000000805
Δ=Q6
0.7101
-0.0828
-0.0007
-0.0000000799
max
1.63
0.061
0.0009
0.0000001845
迭代过程中节点电压变化情况如下表：
k
U1
U2
U3
U4
U5
U6
0
1
1
1
1
1
1
1
1.0334
1.0084
1.0223
1.0372
1.0266
1.0400
2
1.0259
0.9958
1.0128
1.0259
1.0160
1.0324
3
1.0258
0.9956
1.0127
1.0258
1.0159
1.0324
迭代收敛后各节点的电压和功率：
k
U
θ
P
Q
1
1.0258
-2.2168
0.0000
0.0000
2
0.9956
-3.9888
-1.2500
-0.5000
3
1.0127
-3.6874
-0.9000
-0.3000
4
1.0258
3.7197
0.0000
0.0000
5
1.0159
0.7275
-1.0000
-0.3500
6
1.0324
1.9667
0.0000
0.0000
7
1.0250
9.2800
1.6300
0.0665
8
1.0250
4.6648
0.8500
-0.1086
9
1.0400
0.0000
0.7164
0.2705
同课本上给出的潮流相比较，结果完全一致，证明计算过程与程序编写正确。
最后得出迭代收敛后各支路的功率和功率损耗：
i
j
Pij
Qij
Iij
Pji
Qji
Iji
PL
QL
1
2
0.4094
0.2289
0.4572
-0.4068
-0.3869
0.5639
0.0026
-0.1579
1
3
0.3070
0.0103
0.2995
-0.3054
-0.1654
0.3430
0.0017
-0.1551
2
4
-0.8432
-0.1131
0.8545
0.8662
-0.0838
0.8484
0.0230
-0.1969
3
6
-0.5946
-0.1346
0.6020
0.6082
-0.1807
0.6146
0.0135
-0.3153
4
5
0.7638
-0.0080
0.7447
-0.7590
-0.1070
0.7546
0.0048
-0.1150
5
6
-0.2410
-0.2430
0.3368
0.2418
0.0312
0.2362
0.0009
-0.2118
9
1
0.7164
0.2705
0.7363
-0.7164
-0.2392
0.7363
0.0000
0.0312
7
4
1.6300
0.0665
1.5916
-1.6300
0.0918
1.5916
0.0000
0.1583
8
6
0.8500
-0.1086
0.8360
-0.8500
0.1496
0.8360
0.0000
0.0410
三、源程序及注释
由于计算流程比较简单，所以编写程序过程中没有采用模块化的形式，直接按顺序一步步进行。
disp（'【节点数:
】'）;
[n1]=xlsread（'input.xls','A3:
A3'）%节点数
disp（'【支路数:
】'）;
[n]=xlsread（'input.xls','B3:
B3'）%支路数
disp（'【精度:
】'）;
Accuracy=xlsread（'input.xls','B4:
B4'）%精度
[branch]=xlsread（'input.xls','E4:
K12'）;
[node]=xlsread（'input.xls','M4:
S12'）;
Data_B1=branch;%支路参数
Data_B2=node;%节点参数
T1=zeros（n,2）;
T2=zeros（n1,3）;
i=sqrt（-1）;
formatshort
forj=1:
n
T1（j,1）=Data_B1（j,3）+Data_B1（j,4）*1i;
T1（j,2）=Data_B1（j,5）*1i;
end
forj=1:
n1
T2（j,1）=Data_B2（j,1）+Data_B2（j,2）*1i;
T2（j,2）=Data_B2（j,3）+Data_B2（j,4）*1i;
end
B1=zeros（n,6）;
B2=zeros（n1,5）;
forj=1:
n
B1（j,1）=Data_B1（j,1）;
B1（j,2）=Data_B1（j,2）;
B1（j,3）=T1（j,1）;
B1（j,4）=T1（j,2）;
B1（j,5）=Data_B1（j,6）;
B1（j,6）=Data_B1（j,7）;
end
forj=1:
n1
B2（j,1）=T2（j,1）;
B2（j,2）=T2（j,2）;
B2（j,3）=Data_B2（j,5）;
B2（j,4）=Data_B2（j,6）;
B2（j,5）=Data_B2（j,7）;
end
disp（'【支路参数矩阵:
】'）;
B1%显示支路参数矩阵
disp（'【节点参数矩阵:
】'）;
B2%显示节点参数矩阵
%以上为从excel中导入初值的程序
Y=zeros（n1）;
fori=1:
n
ifB1（i,6）==0%不含变压器的支路
p=B1（i,1）;
q=B1（i,2）;
Y（p,q）=Y（p,q）-1/B1（i,3）;
Y（q,p）=Y（p,q）;
Y（p,p）=Y（p,p）+1/B1（i,3）+0.5*B1（i,4）;
Y（q,q）=Y（q,q）+1/B1（i,3）+0.5*B1（i,4）;
else%含有变压器的支路
p=B1（i,1）;
q=B1（i,2）;
Y（p,q）=Y（p,q）-1/（B1（i,3）*B1（i,5））;
Y（q,p）=Y（p,q）;
Y（p,p）=Y（p,p）+1/B1（i,3）;
Y（q,q）=Y（q,q）+1/（B1（i,5）^2*B1（i,3））;
end
end
disp（'【导纳矩阵:
】'）;
Y%显示导纳矩阵
m=0;
fori=1:
n1
ifB2（i,5）==2
m=m+1;
end
end
m%PQ节点个数
l=0;
fori=1:
n1
ifB2（i,5）==1
l=l+1;
end
end
l%PV节点个数
Mismatch_power=zeros（l+m*2,1）;
fori=1:
n1-1
Pj=0;
forj=1:
n1
Pj=Pj+（B2（i,3）*B2（j,3）*（real（Y（i,j））*cos（B2（i,4）-B2（j,4））+imag（Y（i,j））*sin（B2（i,4）-B2（j,4））））;
end
Mismatch_power（i,1）=real（B2（i,1））-real（B2（i,2））-Pj;
end
fork=n1:
（l+m*2）
Qj=0;
forj=1:
n1
Qj=Qj+B2（（k-n1+1）,3）*B2（j,3）*（real（Y（（k-n1+1）,j））*sin（B2（（k-n1+1）,4）-B2（j,4））-imag（Y（（k-n1+1）,j））*cos（B2（（k-n1+1）,4）-B2（j,4）））;
end
Mismatch_power（k,1）=imag（B2（（k-n1+1）,1））-imag（B2（（k-n1+1）,2））-Qj;
end
%Mismatch_power%计算失配功率
times=0;
while（max（Mismatch_power）>Accuracy）
fori=1:
（n1-1）
Pj=0;
forj=1:
n1
Pj=Pj+B2（i,3）*B2（j,3）*（real（Y（i,j））*cos（B2（i,4）-B2（j,4））+imag（Y（i,j））*sin（B2（i,4）-B2（j,4）））;
end
Mismatch_power（i,1）=real（B2（i,1））-real（B2（i,2））-Pj;
end
fork=n1:
（l+m*2）
Qj=0;
forj=1:
n1
Qj=Qj+B2（（k-n1+1）,3）*B2（j,3）*（real（Y（（k-n1+1）,j））*sin（B2（（k-n1+1）,4）-B2（j,4））-imag（Y（（k-n1+1）,j））*cos（B2（（k-n1+1）,4）-B2（j,4）））;
end
Mismatch_power（k,1）=imag（B2（（k-n1+1）,1））-imag（B2（（k-n1+1）,2））-Qj;
end
disp（'【当前迭代次数:
】'）;
times
disp（'【失配功率:
】'）;
Mismatch_power
Jacobian=zeros（l+m*2）;%雅克比矩阵7*7
%————————————————————————————————————H
fori=1:
（n1-1）
forj=1:
（n1-1）
ifi==j
P_H=0;
fork=1:
n1
P_H=P_H+B2（i,3）*B2（k,3）*（real（Y（i,k））*sin（B2（i,4）-B2（k,4））-imag（Y（i,k））*cos（B2（i,4）-B2（k,4）））;
end
Jacobian（i,i）=P_H-B2（i,3）*B2（i,3）*（0-imag（Y（i,i）））;
else
Jacobian（i,j）=0-B2（i,3）*B2（j,3）*（real（Y（i,j））*sin（B2（i,4）-B2（j,4））-imag（Y（i,j））*cos（B2（i,4）-B2（j,4）））;
end
end
end
%————————————————————————————————————N
fori=1:
（n1-1）
forj=1:
m
ifi==j
P_N=0;
fork=1:
n1
P_N=P_N+B2（k,3）*（real（Y（i,k））*cos（B2（i,4）-B2（k,4））+imag（Y（i,k））*sin（B2（i,4）-B2（k,4）））;
end
Jacobian（i,n1-1+i）=0-B2（i,3）*real（Y（i,i））-P_N;
else
Jacobian（i,n1-1+j）=0-B2（i,3）*（real（Y（i,j））*cos（B2（i,4）-B2（j,4））+imag（Y（i,j））*sin（B2（i,4）-B2（j,4）））;
end
end
end
%————————————————————————————————————K
fori=1:
m
forj=1:
（n1-1）
ifi==j
P_K=0;
fork=1:
n1
P_K=P_K+B2（i,3）*B2（k,3）*（real（Y（i,k））*cos（B2（i,4）-B2（k,4））+imag（Y（i,k））*sin（B2（i,4）-B2（k,4）））;
end
Jacobian（n1-1+i,i）=0+B2（i,3）*B2（i,3）*real（Y（i,i））-P_K;
else
Jacobian（n1-1+i,j）=B2（i,3）*B2（j,3）*（real（Y（i,j））*cos（B2（i,4）-B2（j,4））+imag（Y（i,j））*sin（B2（i,4）-B2（j,4）））;
end
end
end
%————————————————————————————————————L
fori=1:
m
forj=1:
m
ifi==j
P_L=0;
fork=1:
n1
P_L=P_L+B2（k,3）*（real（Y（i,k））*sin（B2（i,4）-B2（k,4））-imag（Y（i,k））*cos（B2（i,4）-B2（k,4）））;
end
Jacobian（n1-1+i,n1-1+i）=0-P_L+B2（i,3）*imag（Y（i,i））;
else
Jacobian（n1-1+i,n1-1+j）=0-B2（i,3）*（real（Y（i,j））*sin（B2（i,4）-B2（j,4））-imag（Y（i,j））*cos（B2（i,4）-B2（j,4）））;
end
end
end
S=zeros（l+m*2,1）;%初始化电压角度变化量
S=inv（Jacobian）*（0-Mismatch_power）;%求解修正方程
S=（Jacobian）\（0-Mismatch_power）;%求解修正方程
fori=1:
（n1-1）%角度初值加变化量
B2（i,4）=B2（i,4）+S（i,1）;
end
fori=1:
m%电压初值加变化量
B2（i,3）=B2（i,3）+S（n1-1+i,1）;
end
disp（'【雅克比矩阵:
】'）;
Jacobian%显示雅克比矩阵
%S=inv（Jacobian）
times=times+1;
end
times=times-1;
disp（'【共计迭代次数:
】'）;
times%显示迭代次数
U_It=zeros（n1,1）;%初始化电压向量
fori=1:
n1
U_It（i,1）=B2（i,3）*cos（B2（i,4））+B2（i,3）*sin（B2（i,4））*1j;
end
angle_It=zeros（n1,1）;%将电压角度的弧度值转为角度值
fori=1:
n1
angle_It（i,1）=B2（i,4）*180/pi;
end
Node_S_It=U_It.*（conj（Y）*conj（U_It））;%求解节点功率
disp（'【迭代收敛后各节点的电压幅值:
】'）;
Node_U_It=abs（U_It）%显示迭代收敛后各节点的电压幅值
disp（'【迭代收敛后各节点的电压角度:
】'）;
angle_It%显示迭代收敛后各节点的电压角度
disp（'【迭代收敛后各节点的功率:
】'）;
Node_S_It%显示迭代收敛后各节点的功率
Branch_It=zeros（n,10）;
fori=1:
n;
ifB1（i,6）==0;%不带变压器支路
m=B1（i,1）;%得到支路号
n=B1（i,2）;
Branch_It（i,1）=m;%显示支路号
Branch_It（i,2）=n;
a=U_It（m,1）*（conj（U_It（m,1））*conj（B1（i,4））*0.5+（conj（U_It（m,1））-conj（U_It（n,1）））/conj（B1（i,3）））;
Branch_It（i,3）=real（a）;%显示Pij
Branch_It（i,4）=imag（a）;%显示Qij
b=U_It（m,1）*B1（i,4）*0.5+（U_It（m,1）-U_It（n,1））/B1（i,3）;
Branch_It（i,5）=sqrt（real（b）^2+imag（b）^2）;%显示Iij
c=U_It（n,1）*（conj（U_It（n,1））*conj（B1（i,4））*0.5+（conj（U_It（n,1））-conj（U_It（m,1）））/conj（B1（i,3）））;
Branch_It（i,6）=real（c）;%显示Pji
Branch_It（i,7）=imag（c）;%显示Qji
d=U_It（n,1）*B1（i,4）*0.5+（U_It（n,1）-U_It（m,1））/B1（i,3）;
Branch_It（i,8）=sqrt（real（d）^2+imag（d）^2）;%显示Iji
e=a+c;
Branch_It（i,9）=real（e）;%显示线路损耗有功分量
Branch_It（i,10）=imag（e）;%显示线路损耗无功分量
else%带变压器支路（同以上内容）
m=B1（i,1）;
n=B1（i,2）;
Branch_It（i,1）=m;
Branch_It（i,2）=n;
a=U_It（m,1）*（conj（U_It（m,1））/conj（B1（i,3））-conj（U_It（n,1））*conj（1/（B1（i,5）*B1（i,3））））;
Branch_It（i,3）=real（a）;
Branch_It（i,4）=imag（a）;
b=U_It（m,1）*（B1（i,5）-1）/B1（i,3）/B1（i,5）+（U_It（m,1）-U_It（n,1））/（B1（i,5）*B1（i,3））;
Branch_It（i,5）=sqrt（real（b）^2+imag（b）^2）;
c=U_It（n,1）*（conj（U_It（n,1））/（conj（B1（i,5）*B1（i,5）*B1（i,3）））-conj（U_It（m,1））*conj（1/（B1（i,5）*B1（i,3））））;
Branch_It（i,6）=real（c）;
Branch_It（i,7）=imag（c）;
d=U_It（n,1）*（1-B1（i,5））/B1（i,5）/B1（i,5）/B1（i,3）+（U_It（n,1）-U_It（m,1））/B1（i,5）/B1（i,3）;
Branch_It（i,8）=sqrt（real（d）^2+imag（d）^2）;
e=a+c;
Branch_It（i,9）=real（e）;
Branch_It（i,10）=imag（e）;
end
end
disp（'【迭代收敛后各支路的功率和功率损耗:
】'）;
Branch_It%显示迭代收敛后各支路的功率和功率损耗
%%—————————————————————————————向Excel表中输出数据
%Node_S_It_Real=real（Node_S_It）;
%Node_S_It_imag=imag（Node_S_It）;
%xlswrite（'output.xls',Node_U_It,1,'B3'）;
%xlswrite（'output.xls',angle_It,1,'C3'）;
%xlswrite（'output.xls',Node_S_It_Real,1,'D3'）;
%xlswrite（'output.xls',Node_S_It_imag,1,'E3'）;
%xlswrite（'output.xls',Branch_It,1,'G3'）;
程序中还有将数据从Excel表格中读入输出的xlsread和xlswrite功能，直接将数据输入到Excel表格中，可以省略将数据写在程序中或者一一输入的步骤，适用于任何节点的电力系统潮流计算。
四、程序运行结果
五、手算结果（第一次迭代）
六、总结
通过采用计算机和手算进行潮流计算，我对潮流计算的计算过程和MATLAB软件的使用有了更深层次的了解。
我们已经将如此复杂的问题通过矩阵这样的方式得以简化，但仍然有庞大的计算过程是难以通过手算方式进行解决的。
同时也深深感叹于电力系统的庞大而精细，为自己能在以后为之付出感到期待。
升级会员 