zhouxichao.asv

这个是电力系统的很经典的潮流算法

a1=fopen('date1.txt'); 
b1=fgetl(a1);
fclose(a1);
c=str2num(b1);
S0=c(1);           
V=c(2);
Tj=c(3);
N=c(4);
Ni=c(5);
dt=c(6);
A2=fopen('date2.txt');
C2=fgetl(A2);
fclose(A2);
B=str2num(C2);    
disp('初始功率S0=');
disp(S0);
disp('无穷大母线电压V=');
disp(V);
disp('惯性时间常数Tj=');
disp(Tj);
disp('读入时段数N=');
disp(N);
disp('哪个时段发生故障 Ni=');
disp(Ni);
disp('每段时间间隔 dt=');
disp(dt);
disp('系统等值电抗B=');
disp(B);
for i=1:Ni
    C(i,1)=1;
    C(i,2)=B(1);
end
for i=Ni:N
    C(i,1)=0;
    C(i,2)=B(2);
end
E0=sqrt((V+imag(S0)*B(1)./V)^2+(real(S0)*B(1)./V)^2);
dtj0=atan(real(S0)*B(1)./(V*(V+imag(S0)*B(1)./V)));
Q(1)=dtj0*dtj0*180./pi;
K1=360*50*dt^2./Tj;
P2m=E0*V./B(2);
P3m=E0*V./B(3);
dtk=pi-asin(real(S0)./P3m);
dtjx=acos((real(S0)*(dtk-dtj0)+P3m*cos(dtk)-P2m*cos(dtj0))./(P3m-P2m));
dtjx=dtjx*180./pi;
for K=2:N+1
    P(K-1)=E0*V./C(K-1,2);
    if C(K-1,1)==1
        Peq(K-1)=P(K-1)*sin(Q(K-1)*pi./180);
        DP(K-1)=real(S0)-Peq(K-1);
        if K==2
            DQ(K)=K1*DP(K-1)./2;
        else
            DQ(K)=DQ(K-1)+K1*DP(K-1)./2;
        end
        Q(K)=Q(K-1)+DQ(K);
    elseif C(K-1,1)==0
        P1(K-1)=P(K-2);
        P2(K-1)=P(K-1);
        Peq(K-1)=P1(K-1)*sin(Q(K-1));
        DP(K-1)=real(S0)-Peq(K-1);
        DP2(K-1)=real(S0)-P2(K-1)*sin(Q(K-1));
        DQ(K)=DQ(K-1)+K1*(DP(K-1)+DP2(K-1))./2;
        Q(K)=Q(K-1)+DQ(K);
    end
end
disp('程序计算结果如下：');
disp('不平衡功率DP=');
disp(DP);
disp('功角增量DQ=');
disp(DQ);
disp('功角 Q=');
disp(Q);
disp('电磁功率Peq=');
disp(Peq);
disp('发电机转子摇摆曲线如图所示：');
t=0:dt:dt*N;
plot(t,Q,t,dtjx);
xlabel('时间（秒）');
ylabel('功角（度数）');
title('发电机摇摆曲线');