zantaiwending.asv

电力系统暂态稳定程序

function zantaiwending
clear;
clc;
n=14;%计算时间；
Pt=1;%机械功率；
w=zeros(n,1);
d=zeros(n,1);
e=0.00001;
t=0;
w(1)=1;%w的初值；
f0=50;%频率；
SN=300;SB=220;UB=209;
tj=6*SN./SB;
xd1=0.295;
xd2=0.432;
xt1=0.138;
xl1=0.243;
xl2=0.243;
xt2=0.122;
U=1;Q0=0.2;P0=1;
x1=xd1+xt1+xl1+xt2;%正常运行方式下的系统总电抗；
E1=sqrt((U+Q0*x1./U)^2+(P0*x1./U)^2)%发电机的暂态电动势；
d(1)=(atan(P0*x1./(U^2+Q0*x1)))*180./pi;%功率角的初始值；
h=0.05;%步长；
wzh=1;%故障位置！
%wzh=1故障发生在始端；不然发生在末端；
lx=3;%故障类型；
%lx=1单相短路；lx=2两相短路；lx=3两相接地短路；lx=4三相短路；
if wzh==1
    x22=(xd2+xt1)*(xl2+xt2)./(xd2+xt1+xl2+xt2);%负序电抗；
    x00=xt1*(xl1*4+xt2)./(xt1+xl1*4+xt2);%零序电抗；
else 
    x22=(xd2+xt1+xl2)*(xt2)./(xd2+xt1+xl2+xt2);
    x00=(xl1*4+xt1)*xt2./(xt1+xl1*4+xt2);
end
if wzh==1
    disp('线路始端发生故障！')
else
    disp('线路末端发生故障！')
end
if lx==1
    xa=x22+x00;%附加电抗；
    disp('发生单相接地故障！')
else if lx==2
        xa=x22;
        disp('发生两相短路故障！')
    else if lx==3
            xa=x22*x00./(x22+x00);
            disp('发生两相接地短路故障！')
        else 
            xa=0;
            disp('发生三相短路故障！')
        end
    end
end
x2=(xd1+xt1)+(xl1+xt2)+((xd1+xt1)*(xl1+xt2)./xa);%故障时候的电抗；
P2=E1*U./x2;%故障时发电机的最大功率；
x3=xd1+xt1+2*xl1+xt2;%故障切除后的电抗；
P3=E1*U./x3;%此时的最大功率；
dh=pi-asin(Pt./P3);
dcm=(180./pi)*acos((Pt*(dh-d(1)*pi./180)+P3*cos(dh)-P2*cos(d(1)*pi./180))./(P3-P2);%极限切除角；
disp('极限切除角:')
disp(dcm)
disp('各个时刻的切除角:')
for i=1:n
    j=1;
    if i>=4%故障时间；
        P=P3;
    else
        P=P2;
    end
    d1(i)=(w(i)-1)*360*f0;
    w1(i)=(1./tj)*(Pt-P*sin(d(i)*pi./180));
    d(i+1)=d(i)+h*d1(i);
    w(i+1)=w(i)+h*w1(i);
    d(i+1)=d(i)+h*((w(i)-1)*360*f0+(w(i+1)-1)*360*f0)./2;
    w(i+1)=w(i)+h*((Pt-P*sin(d(i)*pi./180))./tj+(Pt-P*sin(d(i+1)*pi./180))./tj)./2;
    while t==0%迭代过程
    d(i+1)=d(i)+h*((w(i)-1)*360*f0+(w(i+1)-1)*360*f0)./2;
    w(i+1)=w(i)+h*((Pt-P*sin(d(i)*pi./180))./tj+(Pt-P*sin(d(i+1)*pi./180))./tj)./2;
    a(j)=d(i+1);
    b(j)=w(i+1);
    if j==1
        j=j+1;
    else if abs(a(j)-a(j-1))<e&abs(b(j)-b(j-1))<e
            t=1;
        else
            j=j+1;
        end
    end
end
disp(d(i))
end
plot(d)