电力系统分析大作业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+m*2,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=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
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软件的使用有了更深层次的了解。我们已经将如此复杂的问题经过矩阵这样的方式得以简化, 但依然有庞大的计算过程是难以经过手算方式进行解决的。同时也深深感叹于电力系统的庞大而精细, 为自己能在以后为之付出感到期待。