如何用matlab实现无限大功率电源供电系统的三相短路分析

暂态上机计算教学指导书(适用于Matlab)
一、基础知识
1.
为完成本次设计，需要掌握的Matlab的基础知识有:数组的创建和计算，循环语句，条件选择语句，条件判断语句，绘制二维图的语句;
2.
提供了Matlab的学习参考《精通Matlab6.5》光盘版;
3.
Matlab是区分大小写的;
4.
本次设计需要学习和掌握的专业知识:《电力系统暂态分析》相关知识;
5.
matlab
在三角函数计算时默认是弧度，而不是度;
6.
Matlab
数组(矩阵)的下标从1开始;
7.
Matlab
默认小数点后保留4位;
8.
如果不想显示变量的数值，则在该行末尾加“;”符号。
三、程序设计指导
提供了两个程序，(1)不对称短路，名字是buduichen.m;(2)稳定计算的程序，名字是wending.m。所有程序里“%”表示该行为注释行。
1.
不对称短路程序
function
buduichen
%数据来源于教材《电力系统暂态分析》P129例(5-1)
%给学生提供的不对称短路故障的部分程序。
%要求的工作(1)自动形成导纳阵;(2)完善短路类型，两相相间短路和两相短路接地;(3)求短路后发电机的电压
clear
%清除所有变量的值
%输入支路各序阻抗，ZZ1表示支路正序阻抗，ZZ2表示支路负序阻抗，ZZ0表示支路零序阻抗
ZZ1(1,2)=j*0.1;
ZZ1(1,3)=j*0.1;
ZZ1(2,3)=j*0.1;
ZZ2(1,2)=j*0.1;
ZZ2(2,3)=j*0.1;
ZZ2(1,3)=j*0.1;
ZZ0(1,2)=j*0.2;
ZZ0(2,3)=j*0.2;
ZZ0(1,3)=j*0.2;
%输入节点导纳矩阵，Y1表示不计发电机节点的正序网络节点导纳阵，Y2表示不计发电机节点的负序网络节点导纳阵，Y0表示不计发电机节点的零序网络节点导纳阵，
Y1=[-j*26.6266
j*10
j*10
j*10
-j*33.2933
j*10
j*10
j*10
-j*19.96];
Y2=[-j*26.6266
j*10
j*10
j*10
-j*33.2933
j*10
j*10
j*10
-j*19.96];
Y0=[-j*30
j*5
j*5
j*5
-j*50
j*5
j*5
j*5
-j*10];
YY1=[-j*39.96
j*10
j*10
j*20
0
%YY1表示计及发电机节点的正序网络节点导纳阵，YY2表示计及发电机节点的负序网络节点导纳阵
j*10
-j*59.96
j*10
0
j*40
j*10
j*10
-j*19.96
0
0
j*20
0
0
-j*30
0
0
j*40
0
0
-j*60];
YY2=YY1;
N1=3;
%输入网络的节点数
N2=5;
%输入包括所有发电机节点的网络的节点数
k=3;%输入短路点的节点号
fault=1;
%输入短路类型
;f(1)=1单相短路接地;f(2)=2两相相间短路;f(1,1)=3两相短路接地
%第一部分:计算所有节点的a,b,c三相电压
for
p=1:N1
if
p==k
I(p)=1;
else
I(p)=0;
end
end
%正序
Z1=eye(N1)/Y1;%
eye(N1)是3阶的单位矩阵,求节点导纳矩阵Y1的逆阵,即节点阻抗矩阵Z
Z1(:,k)=Z1*I';%短路点k的自阻抗和互阻抗
Zk1=Z1(:,k);
%负序
Z2=eye(N1)/Y2;
Z2(:,k)=Z2*I';
Zk2=Z2(:,k);
%零序
Z0=eye(N1)/Y0;
Z0(:,k)=Z0*I';
Zk0=Z0(:,k);
%计根据故障类型选择不同的计算公式，算故障点各序电流
if
fault==1
Ik1=1/(Z1(k,k)+Z2(k,k)+Z0(k,k));
Ik2=Ik1;Ik0=Ik1;
else
if
fault==2
else
if
fault==3
end
end
end
end
%计算各序电流的故障分量
for
p=1:N1
p
if
p==k
I1(p)=-1*Ik1;
I2(p)=-1*Ik2;
I0(p)=-1*Ik0;
else
I1(p)=0;
I2(p)=0;
I0(p)=0;
end
end
%计算故障点各序电压的故障分量
P96公式(4-13)
uu1(:,k)=Y1\I1.';
%“\”表示左除
A\B=inv(A)*B
inv是求逆阵。
即AX=B
uu2(:,k)=Y2\I2.';
uu0(:,k)=Y0\I0.';
%假设故障前空载，各节点电压正序分量为1
for
p=1:N1
U1(p)=1;
end
%计算故障点发生故障后的电压值:正常运行的分量+故障分量，P37公式(5-39)
u1=U1'+uu1(:,k);
%计算所有节点正序电压
u2=uu2(:,k);
%计算所有节点负序电压
u0=uu0(:,k);
%计算所有节点零序电压
a=-0.5+j*sqrt(3)/2;
T=[1
1
1
%
T为对称分量法的变换矩阵，见P94公式(4-4)
a^2
a
1
a
a^2
1];
for
p=1:N1
U=[u1(p)
u2(p)
u0(p)];
p
Uabc=T*U.'
%Uabc表示p节点的a,b,c三相电压向量,见P94公式(4-4)
UUabc=abs(Uabc)
%UUabc表示p节点的a,b,c三相电压有效值
end
%第二部分:计算支路电流
for
m=1:N1
for
n=1:N1
if
mmn(1)=m;mn(2)=n;
mn
I1(m,n)=(u1(m)-u1(n))/ZZ1(m,n);%支路正序电流
I2(m,n)=(u2(m)-u2(n))/ZZ2(m,n);%支路负序电流
I0(m,n)=(u0(m)-u0(n))/ZZ0(m,n);%支路零序电流
Iabc=T*[I1(m,n)
I2(m,n)
I0(m,n)].';%由序分量合成相分量
Iabc
abs(Iabc)
end
end
end
%第三部分:计算发电机的端电压，只写了部分程序，余下的工作由同学们完成。要考虑变压器绕组是11点接线的
for
p=1:N2
if
p==k
II(p)=-Ik1;
else
II(p)=0;
end
end
vv1(:,k)=YY1\II.';
vv2(:,k)=YY2\II.';
for
p=1:N2
V1(p)=1;
end
2.稳定计算的程序
function
wending
%下面是《电力系统暂态分析》210页例题8-1的计算过程
clear
%清除所有的变量值
PGN=240;
%发电机额定有功功率
SB=220;
%
选取的基准功率
f0=50;
%
系统额定频率
U=1;
%无穷大系统电压标幺值为1
PT=220/SB;
%正常运行时发电机向无穷大系统传输的有功功率标幺值
P=PT;
Q=P*tan(acos(0.98))
Tj=6*(PGN/0.8)/SB
%归算后的发电机惯性时间常数
X1=0.798;
%正常运行时网络等值电抗
X2=2.8;
%故障存在时网络的等值电抗
X3=1.041;
%故障切除时网络的等值电抗
E=sqrt((U+Q*X1/U)^2+(P*X1/U)^2)
%发电机暂态电动势
D0=atan((P*X1/U)/(U+Q*X1/U))*180/pi
%
求正常运行时的功角。
matlab
在三角函数计算时默认是弧度，而不是度。
P2M=E*U/X2
%故障存在时发电机的最大功率
P3M=E*U/X3
%故障存在时发电机的最大功率
DH=180-asin(1/1.35)*180/pi
%
《电力系统暂态分析》P207图8-2中
h点所对应的功角
DCM=acos((PT*(pi/180)*(DH-D0)+P3M*cos(DH*pi/180)-P2M*cos(D0*pi/180))/(P3M-P2M))*180/pi
%求出了极限切除角
%下面是《电力系统暂态分析》214页例题8-2的计算过程,可以用逐段计算法,也可以用改进欧拉法.
%用逐段计算法时,注意(1)在第一个时段求功角增量的公式与其它时段是不同的;(2)在0.15S时功率有跃变
%(1)计算极限切除时间
dt=0.05;
K=360*f0*dt^2/Tj;
DD1=0;
D1=D0;
for
i=1:5
D(i)=D1;
t(i)=dt*(i-1);
P2=P2M*sin(D1*pi/180);
DP2=PT-P2;
if
i==1
DD1=DD1+K*DP2/2;
else
DD1=DD1+K*DP2;
end
D1=D1+DD1;
end
t,D
plot(t,D,'-')
%(1)当故障切除时间为0.15S时计算功角随时间变化的曲线.下面的程序自己写

暂态上机计算教学指导书（适用于Matlab）
一、基础知识
1． 为完成本次设计，需要掌握的Matlab的基础知识有：数组的创建和计算，循环语句，条件选择语句，条件判断语句，绘制二维图的语句；
2． 提供了Matlab的学习参考《精通Matlab6.5》光盘版；
3． Matlab是区分大小写的；
4． 本次设计需要学习和掌握的专业知识：《电力系统暂态分析》相关知识；
5． matlab 在三角函数计算时默认是弧度，而不是度；
6． Matlab 数组（矩阵）的下标从1开始；
7． Matlab 默认小数点后保留4位；
8． 如果不想显示变量的数值，则在该行末尾加“;”符号。

三、程序设计指导
提供了两个程序，（1）不对称短路，名字是buduichen.m；（2）稳定计算的程序，名字是wending.m。所有程序里“%”表示该行为注释行。
1． 不对称短路程序

function buduichen
%数据来源于教材《电力系统暂态分析》P129例（5-1）
%给学生提供的不对称短路故障的部分程序。
%要求的工作（1）自动形成导纳阵；（2）完善短路类型，两相相间短路和两相短路接地；（3）求短路后发电机的电压
clear %清除所有变量的值

%输入支路各序阻抗，ZZ1表示支路正序阻抗，ZZ2表示支路负序阻抗，ZZ0表示支路零序阻抗
ZZ1(1,2)=j*0.1;
ZZ1(1,3)=j*0.1;
ZZ1(2,3)=j*0.1;
ZZ2(1,2)=j*0.1;
ZZ2(2,3)=j*0.1;
ZZ2(1,3)=j*0.1;
ZZ0(1,2)=j*0.2;
ZZ0(2,3)=j*0.2;
ZZ0(1,3)=j*0.2;
%输入节点导纳矩阵，Y1表示不计发电机节点的正序网络节点导纳阵，Y2表示不计发电机节点的负序网络节点导纳阵，Y0表示不计发电机节点的零序网络节点导纳阵，
Y1=[-j*26.6266 j*10 j*10
j*10 -j*33.2933 j*10
j*10 j*10 -j*19.96];
Y2=[-j*26.6266 j*10 j*10
j*10 -j*33.2933 j*10
j*10 j*10 -j*19.96];
Y0=[-j*30 j*5 j*5
j*5 -j*50 j*5
j*5 j*5 -j*10];
YY1=[-j*39.96 j*10 j*10 j*20 0 %YY1表示计及发电机节点的正序网络节点导纳阵，YY2表示计及发电机节点的负序网络节点导纳阵
j*10 -j*59.96 j*10 0 j*40
j*10 j*10 -j*19.96 0 0
j*20 0 0 -j*30 0
0 j*40 0 0 -j*60];
YY2=YY1;

N1=3; %输入网络的节点数
N2=5; %输入包括所有发电机节点的网络的节点数
k=3;%输入短路点的节点号
fault=1; %输入短路类型 ;f(1)=1单相短路接地;f(2)=2两相相间短路;f(1,1)=3两相短路接地

%第一部分：计算所有节点的a,b,c三相电压

for p=1:N1
if p==k
I(p)=1;
else
I(p)=0;
end
end

%正序
Z1=eye(N1)/Y1;% eye(N1)是3阶的单位矩阵,求节点导纳矩阵Y1的逆阵,即节点阻抗矩阵Z
Z1(:,k)=Z1*I';%短路点k的自阻抗和互阻抗
Zk1=Z1(:,k);

%负序
Z2=eye(N1)/Y2;
Z2(:,k)=Z2*I';
Zk2=Z2(:,k);

%零序
Z0=eye(N1)/Y0;
Z0(:,k)=Z0*I';
Zk0=Z0(:,k);

%计根据故障类型选择不同的计算公式，算故障点各序电流
if fault==1
Ik1=1/(Z1(k,k)+Z2(k,k)+Z0(k,k));
Ik2=Ik1;Ik0=Ik1;
else
if fault==2

else
if fault==3

end
end
end
end

%计算各序电流的故障分量
for p=1:N1
p
if p==k
I1(p)=-1*Ik1;
I2(p)=-1*Ik2;
I0(p)=-1*Ik0;
else
I1(p)=0;
I2(p)=0;
I0(p)=0;
end
end

%计算故障点各序电压的故障分量 P96公式（4-13）
uu1(:,k)=Y1\I1.'; %“\”表示左除 A\B=inv(A)*B inv是求逆阵。 即AX=B
uu2(:,k)=Y2\I2.';
uu0(:,k)=Y0\I0.';

%假设故障前空载，各节点电压正序分量为1
for p=1:N1
U1(p)=1;
end

%计算故障点发生故障后的电压值:正常运行的分量+故障分量，P37公式（5-39）
u1=U1'+uu1(:,k); %计算所有节点正序电压
u2=uu2(:,k); %计算所有节点负序电压
u0=uu0(:,k); %计算所有节点零序电压

a=-0.5+j*sqrt(3)/2;
T=[1 1 1 % T为对称分量法的变换矩阵，见P94公式（4-4）
a^2 a 1
a a^2 1];

for p=1:N1
U=[u1(p) u2(p) u0(p)];
p
Uabc=T*U.' %Uabc表示p节点的a,b,c三相电压向量,见P94公式（4-4）
UUabc=abs(Uabc) %UUabc表示p节点的a,b,c三相电压有效值
end

%第二部分：计算支路电流
for m=1:N1
for n=1:N1
if m mn(1)=m;mn(2)=n;
mn
I1(m,n)=(u1(m)-u1(n))/ZZ1(m,n);%支路正序电流
I2(m,n)=(u2(m)-u2(n))/ZZ2(m,n);%支路负序电流
I0(m,n)=(u0(m)-u0(n))/ZZ0(m,n);%支路零序电流
Iabc=T*[I1(m,n) I2(m,n) I0(m,n)].';%由序分量合成相分量
Iabc
abs(Iabc)
end
end
end

%第三部分：计算发电机的端电压，只写了部分程序，余下的工作由同学们完成。要考虑变压器绕组是11点接线的
for p=1:N2
if p==k
II(p)=-Ik1;
else
II(p)=0;
end
end
vv1(:,k)=YY1\II.';
vv2(:,k)=YY2\II.';
for p=1:N2
V1(p)=1;
end

2．稳定计算的程序
function wending
%下面是《电力系统暂态分析》210页例题8-1的计算过程
clear %清除所有的变量值
PGN=240; %发电机额定有功功率
SB=220; % 选取的基准功率
f0=50; % 系统额定频率
U=1; %无穷大系统电压标幺值为1
PT=220/SB; %正常运行时发电机向无穷大系统传输的有功功率标幺值
P=PT;
Q=P*tan(acos(0.98))
Tj=6*(PGN/0.8)/SB %归算后的发电机惯性时间常数
X1=0.798; %正常运行时网络等值电抗
X2=2.8; %故障存在时网络的等值电抗
X3=1.041; %故障切除时网络的等值电抗
E=sqrt((U+Q*X1/U)^2+(P*X1/U)^2) %发电机暂态电动势
D0=atan((P*X1/U)/(U+Q*X1/U))*180/pi % 求正常运行时的功角。 matlab 在三角函数计算时默认是弧度，而不是度。
P2M=E*U/X2 %故障存在时发电机的最大功率
P3M=E*U/X3 %故障存在时发电机的最大功率
DH=180-asin(1/1.35)*180/pi % 《电力系统暂态分析》P207图8-2中 h点所对应的功角
DCM=acos((PT*(pi/180)*(DH-D0)+P3M*cos(DH*pi/180)-P2M*cos(D0*pi/180))/(P3M-P2M))*180/pi %求出了极限切除角

%下面是《电力系统暂态分析》214页例题8-2的计算过程,可以用逐段计算法,也可以用改进欧拉法.
%用逐段计算法时,注意(1)在第一个时段求功角增量的公式与其它时段是不同的;(2)在0.15S时功率有跃变

%(1)计算极限切除时间
dt=0.05;
K=360*f0*dt^2/Tj;
DD1=0;
D1=D0;
for i=1:5
D(i)=D1;
t(i)=dt*(i-1);
P2=P2M*sin(D1*pi/180);
DP2=PT-P2;
if i==1
DD1=DD1+K*DP2/2;
else
DD1=DD1+K*DP2;
end
D1=D1+DD1;
end
t,D
plot(t,D,'-')

%(1)当故障切除时间为0.15S时计算功角随时间变化的曲线.下面的程序自己写