电力系统分析大作业matlab三机九节点潮流计算报告.doc

电力系统分析大作业

一、设计题目
本次设计题目选自课本第五章例5-8,美国西部联合电网WSCC系统简化三机九节点系统，例题中已经给出了潮流成果，计算成果可以与之对照。取ε=0.00001 。

二、计算环节
第一步，为了以便编程，修改节点序号，将平衡节点放在最后。如下图：
9
2
1
3
2
7
4
5
6
8
3

第二步，这样得出系统参数如下表所示：

第三步，形成节点导纳矩阵。

第四步，设定初值：
；
，。
第五步，计算失配功率
=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.
Δ=P2
-1.25
0.0379
0.0005
0.
Δ=P3
-0.9
0.0439
0.0005
0.
Δ=P4
0
-0.0421
-0.0012
-0.
Δ=P5
-1
0.061
0.0009
0.
Δ=P6
0
-0.0269
-0.0007
-0.
Δ=P7
1.63
-0.0579
-0.0004
-0.
Δ=P8
0.85
-0.0336
-0.0002
-0.
Δ=Q1
0.8614
-0.0501
-0.0004
-0.
Δ=Q2
-0.259
-0.0714
-0.0012
-0.
Δ=Q3
-0.042
-0.0424
-0.0006
-0.
Δ=Q4
0.6275
-0.1875
-0.0021
-0.
Δ=Q5
-0.171
-0.0241
-0.0004
-0.
Δ=Q6
0.7101
-0.0828
-0.0007
-0.
max
1.63
0.061
0.0009
0.
迭代过程中节点电压变化状况如下表：
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);
format short
for j=1:n
T1(j,1)=Data_B1(j,3)+Data_B1(j,4)*1i;
T1(j,2)=Data_B1(j,5)*1i;
end
for j=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);
for j=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
for j=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);
for i=1:n
if B1(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;
for i=1:n1
if B2(i,5)==2
m=m+1;
end
end
m %PQ节点个数
l=0;
for i=1:n1
if B2(i,5)==1
l=l+1;
end
end
l %PV节点个数
Mismatch_power=zeros(l+m2,1);
for i=1:n1-1
Pj=0;
for j=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
for k=n1:(l+m*2)
Qj=0;
for j=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)
for i=1:(n1-1)
Pj=0;
for j=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
for k=n1:(l+m*2)
Qj=0;
for j=1n1
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
for i=1:(n1-1)
for j=1:(n1-1)
if i==j
P_H=0;
for k=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
for i=1:(n1-1)
for j=1:m
if i==j
P_N=0;
for k=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
for i=1:m
for j=1:(n1-1)
if i==j
P_K=0;
for k=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
for i=1:m
for j=1:m
if i==j
P_L=0;
for k=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)；%求解修正方程
for i=1:(n1-1) %角度初值加变化量
B2(i,4)=B2(i,4)+S(i,1);
end
for i=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)；%初始化电压向量
for i=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)；%将电压角度弧度值转为角度值
for i=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);
for i=1:n;
if B1(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软件使用有了更深层次理解。咱们已经将如此复杂问题通过矩阵这样方式得以简化，但依然有庞大计算过程是难以通过手算方式进行解决。同步也深深感叹于电力系统庞大而精细，为自己能在后来为之付出感到期待。