http://www.apollocode.net/a/1086.html
电力系统暂态稳定分析程序，
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 ----------------------------------