应用介绍

电力系统暂态稳定分析程序，

clear

clc

basemva = 100;  accuracy = 0.0001;  maxiter = 10;

%      Bus BusbVoltage Angle ---Load----  -----Generator-----Static Mvar

%       No code Mag. Degree  MW     Mvar      MW      Mvar   Qmin   Qmax Qc/-Ql

busdata=[1  1  1.06   0.0  00.00    00.00    0.00    00.00     0     0   0

         2  2  1.04   0.0  00.00    00.00  150.00    00.00     0   140   0

         3  2  1.03   0.0  00.00    00.00  100.00    00.00     0    90   0

         4  0  1.0    0.0 100.00    70.00   00.00    00.00     0    0    0

         5  0  1.0    0.0  90.00    30.00   00.00    00.00     0    0    0

         6  0  1.0    0.0 160.00   110.00   00.00    00.00     0    0    0];

%                                            Line code

%        Bus bus   R       X     1/2 B    = 1 for lines

%        nl  nr  p.u.     p.u.   p.u.     >1 or<1 tr. tap at bus nl

linedata=[1  4  0.035   0.225   0.0065   1.0

          1  5  0.025   0.105   0.0045   1.0

          1  6  0.040   0.215   0.0055   1.0

          2  4  0.000   0.035   0.0000   1.0

          3  5  0.000   0.042   0.0000   1.0

          4  6  0.028   0.125   0.0035   1.0

          5  6  0.026   0.175   0.0300   1.0];

 %       Gen.  Ra   Xd'    H

gendata=[ 1     0   0.20   20

          2     0   0.15    4

          3     0   0.25    5];

  %%%-------------------- Ene of input data --------------------------------

 %%------------------------Start creat Ybus----------------------------------

  j=sqrt(-1); i = sqrt(-1);

nl = linedata(:,1); nr = linedata(:,2); R = linedata(:,3);

X = linedata(:,4); Bc = j*linedata(:,5); a = linedata(:, 6);

nbr=length(linedata(:,1)); nbus = max(max(nl), max(nr));

Z = R + j*X; y= ones(nbr,1)./Z;        %branch admittance

for n = 1:nbr

if a(n) <= 0  a(n) = 1;

else end

Ybus=zeros(nbus,nbus);    

for k=1:nbr;

       Ybus(nl(k),nr(k))=Ybus(nl(k),nr(k))-y(k)/a(k);

       Ybus(nr(k),nl(k))=Ybus(nl(k),nr(k));

    end

end

for  n=1:nbus

     for k=1:nbr

         if nl(k)==n

         Ybus(n,n) = Ybus(n,n)+y(k)/(a(k)^2) + Bc(k);

         elseif nr(k)==n

         Ybus(n,n) = Ybus(n,n)+y(k) +Bc(k);

         else, end

     end

end

clear Pgg

  %%------------------------END creat Ybus----------------------------------

  %%------------------------Start Newton-Raphson load Flow ----------------------------------

  ns=0; ng=0; Vm=0; delta=0; yload=0; deltad=0;

nbus = length(busdata(:,1));

for k=1:nbus

n=busdata(k,1);

kb(n)=busdata(k,2); Vm(n)=busdata(k,3); delta(n)=busdata(k, 4);

Pd(n)=busdata(k,5); Qd(n)=busdata(k,6); Pg(n)=busdata(k,7); Qg(n) = busdata(k,8);

Qmin(n)=busdata(k, 9); Qmax(n)=busdata(k, 10);

Qsh(n)=busdata(k, 11);

    if Vm(n) <= 0  Vm(n) = 1.0; V(n) = 1 + j*0;

    else delta(n) = pi/180*delta(n);

         V(n) = Vm(n)*(cos(delta(n)) + j*sin(delta(n)));

         P(n)=(Pg(n)-Pd(n))/basemva;

         Q(n)=(Qg(n)-Qd(n)+ Qsh(n))/basemva;

         S(n) = P(n) + j*Q(n);

    end

end

for k=1:nbus

if kb(k) == 1, ns = ns+1; else, end

if kb(k) == 2 ng = ng+1; else, end

ngs(k) = ng;

nss(k) = ns;

end

Ym=abs(Ybus); t = angle(Ybus);

m=2*nbus-ng-2*ns;

maxerror = 1; converge=1;

iter = 0;

% Start of iterations

clear A  DC   J  DX

while maxerror >= accuracy & iter <= maxiter 

for i=1:m

for k=1:m

   A(i,k)=0;      

end, end

iter = iter+1;

for n=1:nbus

nn=n-nss(n);

lm=nbus+n-ngs(n)-nss(n)-ns;

J11=0; J22=0; J33=0; J44=0;

   for i=1:nbr

     if nl(i) == n | nr(i) == n

        if nl(i) == n,  l = nr(i); end

        if nr(i) == n,  l = nl(i); end

        J11=J11+ Vm(n)*Vm(l)*Ym(n,l)*sin(t(n,l)- delta(n) + delta(l));

        J33=J33+ Vm(n)*Vm(l)*Ym(n,l)*cos(t(n,l)- delta(n) + delta(l));

        if kb(n)~=1

        J22=J22+ Vm(l)*Ym(n,l)*cos(t(n,l)- delta(n) + delta(l));

        J44=J44+ Vm(l)*Ym(n,l)*sin(t(n,l)- delta(n) + delta(l));

        else, end

        if kb(n) ~= 1  & kb(l) ~=1

        lk = nbus+l-ngs(l)-nss(l)-ns;

        ll = l -nss(l);

        A(nn, ll) =-Vm(n)*Vm(l)*Ym(n,l)*sin(t(n,l)- delta(n) + delta(l));

              if kb(l) == 0 

              A(nn, lk) =Vm(n)*Ym(n,l)*cos(t(n,l)- delta(n) + delta(l));end

              if kb(n) == 0  

              A(lm, ll) =-Vm(n)*Vm(l)*Ym(n,l)*cos(t(n,l)- delta(n)+delta(l)); end

              if kb(n) == 0 & kb(l) == 0 

              A(lm, lk) =-Vm(n)*Ym(n,l)*sin(t(n,l)- delta(n) + delta(l));end

        else end

     else , end

   end

   Pk = Vm(n)^2*Ym(n,n)*cos(t(n,n))+J33;

   Qk = -Vm(n)^2*Ym(n,n)*sin(t(n,n))-J11;

   if kb(n) == 1 P(n)=Pk; Q(n) = Qk; end   

     if kb(n) == 2  Q(n)=Qk;

         if Qmax(n) ~= 0

           Qgc = Q(n)*basemva + Qd(n) - Qsh(n);

           if iter <= 7                 

              if iter > 2               

                if Qgc  < Qmin(n),       

                Vm(n) = Vm(n) + 0.01;    

                elseif Qgc  > Qmax(n),   

                Vm(n) = Vm(n) - 0.01;end 

              else, end

           else,end

         else,end

     end

   if kb(n) ~= 1

     A(nn,nn) = J11;  

     DC(nn) = P(n)-Pk;

   end

   if kb(n) == 0

     A(nn,lm) = 2*Vm(n)*Ym(n,n)*cos(t(n,n))+J22;  

     A(lm,nn)= J33;        

     A(lm,lm) =-2*Vm(n)*Ym(n,n)*sin(t(n,n))-J44; 

     DC(lm) = Q(n)-Qk;

   end

end

DX=A\DC';

for n=1:nbus

  nn=n-nss(n);

  lm=nbus+n-ngs(n)-nss(n)-ns;

    if kb(n) ~= 1

    delta(n) = delta(n)+DX(nn); end

    if kb(n) == 0

    Vm(n)=Vm(n)+DX(lm); end

 end

  maxerror=max(abs(DC));

     if iter == maxiter & maxerror > accuracy 

   fprintf('\nWARNING: Iterative solution did not converged after ')

   fprintf('%g', iter), fprintf(' iterations.\n\n')

   fprintf('Press Enter to terminate the iterations and print the results \n')

   converge = 0; pause, else, end

end

if converge ~= 1

   tech= ('                      ITERATIVE SOLUTION DID NOT CONVERGE'); else, 

   tech=('                   Power Flow Solution by Newton-Raphson Method');

end   

V = Vm.*cos(delta)+j*Vm.*sin(delta);

deltad=180/pi*delta;

i=sqrt(-1);

k=0;

for n = 1:nbus

     if kb(n) == 1

     k=k+1;

     S(n)= P(n)+j*Q(n);

     Pg(n) = P(n)*basemva + Pd(n);

     Qg(n) = Q(n)*basemva + Qd(n) - Qsh(n);

     Pgg(k)=Pg(n);

     Qgg(k)=Qg(n);    

     elseif  kb(n) ==2

     k=k+1;

     S(n)=P(n)+j*Q(n);

     Qg(n) = Q(n)*basemva + Qd(n) - Qsh(n);

     Pgg(k)=Pg(n);

     Qgg(k)=Qg(n);  

  end

yload(n) = (Pd(n)- j*Qd(n)+j*Qsh(n))/(basemva*Vm(n)^2);

end

busdata(:,3)=Vm'; busdata(:,4)=deltad';

Pgt = sum(Pg);  Qgt = sum(Qg); Pdt = sum(Pd); Qdt = sum(Qd); Qsht = sum(Qsh);

  %%------------------------End Newton-Raphson load Flow ----------------------------------

%%------------------------Start transient Stability program ----------------------------------

global Pm f H E  Y th ngg

f=60;

ngr=gendata(:,1);

ngg=length(gendata(:,1));

%%

for k=1:ngg

zdd(ngr(k))=gendata(k, 2)+j*gendata(k,3);

H(k)=gendata(k,4);  

end

%%

for k=1:ngg

I=conj(S(ngr(k)))/conj(V(ngr(k)));

Ep(k) = V(ngr(k))+zdd(ngr(k))*I; 

Pm(k)=real(S(ngr(k)));            

end

E=abs(Ep); d0=angle(Ep);

for k=1:ngg

nl(nbr+k) = nbus+k;

nr(nbr+k) = gendata(k, 1);

R(nbr+k)  = real(zdd(ngr(k)));

X(nbr+k)  = imag(zdd(ngr(k)));

Bc(nbr+k)  = 0;

a(nbr+k) = 1.0;

yload(nbus+k)=0;

end

nbr1=nbr; nbus1=nbus;

nbrt=nbr+ngg;

nbust=nbus+ngg;

linedata=[nl, nr, R, X, -j*Bc, a];

%%%%%%%%---------- creat Ybus  before Fault -----------------------

%%-----------------Lowd flow YBUS -------------------------------

j=sqrt(-1); i = sqrt(-1);

nl = linedata(:,1); nr = linedata(:,2); R = linedata(:,3);

X = linedata(:,4); Bc = j*linedata(:,5); a = linedata(:, 6);

nbr=length(linedata(:,1)); nbus = max(max(nl), max(nr));

Z = R + j*X; y= ones(nbr,1)./Z;     

for n = 1:nbr

if a(n) <= 0  a(n) = 1; else end

Ybus=zeros(nbus,nbus);     

for k=1:nbr;

       Ybus(nl(k),nr(k))=Ybus(nl(k),nr(k))-y(k)/a(k);

       Ybus(nr(k),nl(k))=Ybus(nl(k),nr(k));

    end

end

for  n=1:nbus

     for k=1:nbr

         if nl(k)==n

         Ybus(n,n) = Ybus(n,n)+y(k)/(a(k)^2) + Bc(k);

         elseif nr(k)==n

         Ybus(n,n) = Ybus(n,n)+y(k) +Bc(k);

         else, end

     end

end

clear Pgg

%%-----------------End Lowd flow YBUS -------------------------------

for k=1:nbust

Ybus(k,k)=Ybus(k,k)+yload(k);

end

YLL=Ybus(1:nbus1, 1:nbus1);

YGG = Ybus(nbus1+1:nbust, nbus1+1:nbust);

YLG = Ybus(1:nbus1, nbus1+1:nbust);

Ybf=YGG-YLG.'*inv(YLL)*YLG;

clear YLL YGG YLG

%%%%%%%%---------- END Ybus  before Fault -----------------------

Y=abs(Ybf); th=angle(Ybf);

Pm=zeros(1, ngg);

for ii = 1:ngg

for jj = 1:ngg

Pm(ii) = Pm(ii) + E(ii)*E(jj)*Y(ii, jj)*cos(th(ii, jj)-d0(ii)+d0(jj));

end,

end

nf=input('Please Enter fault location (Bus Number ?) ');

%%%------------- Ybus during fault ---------------------------------------

nbusf=nbust-1;

Ybus(:,nf:nbusf)=Ybus(:,nf+1:nbust);

Ybus(nf:nbusf,:)=Ybus(nf+1:nbust,:);

YLL=Ybus(1:nbus1-1, 1:nbus1-1);

YGG = Ybus(nbus1:nbusf, nbus1:nbusf);

YLG = Ybus(1:nbus1-1, nbus1:nbusf);

Ydf=YGG-YLG.'*inv(YLL)*YLG;

%%%%%--------------------------------------------------------------

%%%%%--------------- Ybus after Fault -------------------------------

nl=linedata(:, 1);  

nr=linedata(:, 2);

remove = 0;

while remove ~= 1

  fline=input('Enter removed line? for example [5,6]---->? ');

  nlf=fline(1); nrf=fline(2);

     for k=1:nbrt

       if nl(k)==nlf & nr(k)==nrf

       remove = 1;

       m=k;

       else, end

    end

    if remove ~= 1

       fprintf('\nThe line to be removed does not exist in the line data. try again.\n\n')

    end

end

linedat2(1:m-1,:)= linedata(1:m-1,:);

linedat2(m:nbrt-1,:)=linedata(m+1:nbrt,:);

linedat0=linedata;

linedata=linedat2;

%%%%%-----------------Load Flow  YBUS -------------------------------------

j=sqrt(-1); i = sqrt(-1);

nl = linedata(:,1); nr = linedata(:,2); R = linedata(:,3);

X = linedata(:,4); Bc = j*linedata(:,5); a = linedata(:, 6);

nbr=length(linedata(:,1)); nbus = max(max(nl), max(nr));

Z = R + j*X; y= ones(nbr,1)./Z;      

for n = 1:nbr

if a(n) <= 0  a(n) = 1; else end

Ybus=zeros(nbus,nbus);  

for k=1:nbr;

       Ybus(nl(k),nr(k))=Ybus(nl(k),nr(k))-y(k)/a(k);

       Ybus(nr(k),nl(k))=Ybus(nl(k),nr(k));

    end

end

for  n=1:nbus

     for k=1:nbr

         if nl(k)==n

         Ybus(n,n) = Ybus(n,n)+y(k)/(a(k)^2) + Bc(k);

         elseif nr(k)==n

         Ybus(n,n) = Ybus(n,n)+y(k) +Bc(k);

         else, end

     end

end

clear Pgg

%%%%%-----------------End of Load Flow  YBUS -------------------------------------

for k=1:nbust

Ybus(k,k)=Ybus(k,k)+yload(k);

end

YLL=Ybus(1:nbus1, 1:nbus1);

YGG = Ybus(nbus1+1:nbust, nbus1+1:nbust);

YLG = Ybus(1:nbus1, nbus1+1:nbust);

Yaf=YGG-YLG.'*inv(YLL)*YLG;

linedata=linedat0;

%%%%%--------------- End Ybus after Fault -------------------------------

tc=input('Enter  time of  removing fault in sec. tc = ');

tf=input('Enter  simulation time in sec.  tf = ');

clear t  x  del

t0 = 0;

w0=zeros(1, length(d0));

x0 = [d0,  w0];

tol=0.0001;

Y=abs(Ydf); th=angle(Ydf);

tspan=[t0, tc];                                     

[t1, xf] =ode23('trstability', tspan, x0); 

x0c =xf(length(xf), :);

Y=abs(Yaf); th=angle(Yaf);

tspan = [tc, tf];                         

[t2,xc] =ode23('trstability', tspan, x0c);      

t =[t1; t2]; x = [xf; xc];

for k=1:nbus1

    if kb(k)==1

    ms=k; 

        end

end

kk=0;

for k=1:ngg

    if k~=ms

    kk=kk+1;

    del(:,kk)=180/pi*(x(:,k)-x(:,ms));

    else, end

end

h=figure; figure(h)

plot(t, del)

title(['Phase angle difference (fault cleared at ', num2str(tc),'s)'])

xlabel('t, sec'), ylabel('Delta, degree'), grid

%%------------------------End transient Stability program ----------------------------------