Afpek

% State variable representation of the multimachine system

% after fault. (for use with trstab)
%
% Copyright (C) 1998 H. Saadat

function xdot = afpek(t,x)
global Pm f H E Y th ngg

Pe=zeros(1, ngg);
for ii = 1:ngg
for jj = 1:ngg
Pe(ii) = Pe(ii) + E(ii)*E(jj)*Y(ii, jj)*cos(th(ii, jj)-
x(ii)+x(jj));
end, end
for k=1:ngg
xdot(k)=x(k+ngg);
xdot(k+ngg)=(pi*f)/H(k)*(Pm(k)-Pe(k));
end
xdot=xdot'; % use with MATLAB 5 (remove for MATLAB 4)


busout
% This program prints the power flow solution in a tabulated form
% on the screen.
%
% Copyright (C) 1998 by H. Saadat.

%clc
disp(tech)
fprintf(' Maximum Power Mismatch = %g \n',
maxerror)
fprintf(' No. of Iterations = %g
\n\n', iter)
head =[' Bus Voltage Angle ------Load------ ---
Generation--- Injected'
' No. Mag. Degree MW Mvar MW
Mvar Mvar '
'
'];
disp(head)
for n=1:nbus
fprintf(' %5g', n), fprintf(' %7.3f', Vm(n)),
fprintf(' %8.3f', deltad(n)), fprintf(' %9.3f', Pd(n)),
fprintf(' %9.3f', Qd(n)), fprintf(' %9.3f', Pg(n)),
fprintf(' %9.3f ', Qg(n)), fprintf(' %8.3f\n', Qsh(n))
end
fprintf(' \n'), fprintf(' Total ')
fprintf(' %9.3f', Pdt), fprintf(' %9.3f', Qdt),
fprintf(' %9.3f', Pgt), fprintf(' %9.3f', Qgt), fprintf('
%9.3f\n\n', Qsht)

defpek
% State variable representation of the multimachine system
% during fault. (for use with trstab)
% Copyright (c) 1998 by H. Saadat

function xdot = dfpek(t,x)
global Pm f H E Y th ngg
Pe=zeros(1, ngg);
for ii = 1:ngg
for jj = 1:ngg
Pe(ii) = Pe(ii) + E(ii)*E(jj)*Y(ii, jj)*cos(th(ii, jj)-
x(ii)+x(jj));
end, end
for k=1:ngg
xdot(k)=x(k+ngg);
xdot(k+ngg)=(pi*f)/H(k)*(Pm(k)-Pe(k));
end
xdot=xdot'; % use with MATLAB 5 (remove for MATLAB 4)


lf newton
% Power flow solution by Newton-Raphson method
% Copyright (c) 1998 by H. Saadat
% Revision 1 (Aug. 99) To include two or more parallel lines
ns=0; ng=0; Vm=0; delta=0; yload=0; deltad=0;
nbus = length(busdata(:,1));
kb=[];Vm=[]; delta=[]; Pd=[]; Qd=[]; Pg=[]; Qg=[]; Qmin=[];
Qmax=[]; % Added (6-8-00)
Pk=[]; P=[]; Qk=[]; Q=[]; S=[]; V=[]; % Added (6-8-00)
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;
%%%% added for parallel lines (Aug. 99)
mline=ones(nbr,1);
for k=1:nbr
for m=k+1:nbr
if((nl(k)==nl(m)) & (nr(k)==nr(m)));
mline(m)=2;
elseif ((nl(k)==nr(m)) & (nr(k)==nl(m)));
mline(m)=2;
else, end
end
end
%%% end of statements for parallel lines (Aug. 99)

% Start of iterations
clear A DC J DX
while maxerror >= accuracy & iter <= maxiter % Test for max. power
mismatch
for ii=1:m
for k=1:m
A(ii,k)=0; %Initializing Jacobian matrix
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 ii=1:nbr
if mline(ii)==1 % Added to include parallel lines (Aug. 99)
if nl(ii) == n | nr(ii) == n
if nl(ii) == n , l = nr(ii); end
if nr(ii) == n , l = nl(ii); 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);
% off diagonalelements of J1
A(nn, ll) =-Vm(n)*Vm(l)*Ym(n,l)*sin(t(n,l)-
delta(n) + delta(l));
if kb(l) == 0 % off diagonal elements of J2
A(nn, lk) =Vm(n)*Ym(n,l)*cos(t(n,l)- delta(n)
+ delta(l));end
if kb(n) == 0 % off diagonal elements of J3
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 % off diagonal
elements of J4
A(lm, lk) =-Vm(n)*Ym(n,l)*sin(t(n,l)- delta(n)
+ delta(l));end
else 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 % Swing bus P
if kb(n) == 2 Q(n)=Qk;
if Qmax(n) ~= 0
Qgc = Q(n)*basemva + Qd(n) - Qsh(n);
if iter <= 7 % Between the 2th & 6th
iterations
if iter > 2 % the Mvar of generator
buses are
if Qgc < Qmin(n), % tested. If not within
limits Vm(n)
Vm(n) = Vm(n) + 0.01; % is changed in steps of
0.01 pu to
elseif Qgc > Qmax(n), % bring the generator
Mvar within
Vm(n) = Vm(n) - 0.01;end % the specified limits.
else, end
else,end
else,end
end
if kb(n) ~= 1
A(nn,nn) = J11; %diagonal elements of J1
DC(nn) = P(n)-Pk;
end
if kb(n) == 0
A(nn,lm) = 2*Vm(n)*Ym(n,n)*cos(t(n,n))+J22; %diagonal
elements of J2
A(lm,nn)= J33; %diagonal elements of J3
A(lm,lm) =-2*Vm(n)*Ym(n,n)*sin(t(n,n))-J44; %diagonal of
elements of J4
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); %june 97
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); % June 1997
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);

%clear A DC DX J11 J22 J33 J44 Qk delta lk ll lm
%clear A DC DX J11 J22 J33 Qk delta lk ll lm

Lfybus
% This program obtains th Bus Admittance Matrix for power flow
solution
% Copyright (c) 1998 by H. Saadat

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); % initialize Ybus to zero
% formation of the off diagonal elements
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
% formation of the diagonal elements
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


trstab
% TRANSIENT STABILITY ANALYSIS OF A MULTIMACHINE POWER SYSTEM
NETWORK
%
% Program ``trstab'' must be used in conjuntion with the power
flow
% program. Any of the power flow programs ``lfgauss'',
``lfnewton'',
% ``decoupled'' or ``perturb'' can be used prior to ``trstab''
program.
% Power flow program provide the power, voltage magnitude and
phase
% angle for each bus. Also, the load admittances are returned in a
% matrix named ``yload''. In addition to the required power flow
data,
% Transient reactance, and inertia constant of each machine must
be
% specified. This is defined in a matrix named ``gendata'', Each
row
% contain the bus No. to which a generator is connected, armature
% resistance, transient reactance, and the machine inertia
constant.
% Program ``trstab'' obtains the prefault bus admittance matrix
including
% the load admittances. Voltage behind transient reactances are
obtained.
% The reduced admittance matrix before, during and after fault are
found.
% Machine equations are expressed in state variable form and the
MATLAB
% ode23 is used to solve the multimachine equations. The phase
angle
% difference of each machine with respect to the slack bus are
plotted.
% The simulation can be repeated for a different fault clearing
time, or
% a different fault location.
%
% Copyright (c) 1998 by H. Saadat
%
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); % new
end
for k=1:ngg
I=conj(S(ngr(k)))/conj(V(ngr(k)));
Ep(k) = V(ngr(k))+zdd(ngr(k))*I; % new
Pm(k)=real(S(ngr(k))); % new
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];
[Ybus, Ybf]=ybusbf(linedata, yload, nbus1,nbust);
fprintf('\nPrefault reduced bus admittance matrix \n')
Ybf
Y=abs(Ybf); th=angle(Ybf);
Pm=zeros(1, ngg);
disp([' G(i) E''(i) d0(i) Pm(i)'])
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,
fprintf(' %g', ngr(ii)), fprintf(' %8.4f',E(ii)),
fprintf(' %8.4f', 180/pi*d0(ii))
fprintf(' %8.4f \n',Pm(ii))
end
respfl='y';
while respfl =='y' | respfl=='Y'
nf=input('Enter faulted bus No. -> ');
rtn=isempty(nf);
if rtn==1; nf=-1; end
while nf <= 0 | nf > nbus
fprintf('Faulted bus No. must be between 1 & %g \n', nbus)
nf = input('Enter Faulted Bus No. -> ');
rtn=isempty(nf);
if rtn==1; nf=-1; end
end

fprintf('\nFaulted reduced bus admittance matrix\n')
Ydf=ybusdf(Ybus, nbus1, nbust, nf)
%Fault cleared
[Yaf]=ybusaf(linedata, yload, nbus1,nbust, nbrt);
fprintf('\nPostfault reduced bus admittance matrix\n')
Yaf
resptc='y';
while resptc =='y' | resptc=='Y'
tc=input('Enter clearing time of fault in sec. tc = ');
tf=input('Enter final 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);
%[t1, xf] =ode23('dfpek', t0, tc, x0, tol); % Solution during
fault (use with MATLAB 4)
tspan=[t0, tc]; %use with
MATAB 5
[t1, xf] =ode23('dfpek', tspan, x0); % Solution during fault (use
with MATLAB 5)
x0c =xf(length(xf), :);
Y=abs(Yaf); th=angle(Yaf);
%[t2,xc] =ode23('afpek', tc, tf, x0c, tol); % Postfault solution
(use with MATLAB 4)
tspan = [tc, tf]; % use with MATLAB 5
[t2,xc] =ode23('afpek', tspan, x0c); % Postfault solution
(use with MATLAB 5)
t =[t1; t2]; x = [xf; xc];
fprintf('\nFault is cleared at %4.3f Sec. \n', tc)
for k=1:nbus
if kb(k)==1
ms=k; else, end
end
fprintf('\nPhase angle difference of each machine \n')
fprintf('with respect to the slack in degree.\n')
fprintf(' t - sec')
kk=0;
for k=1:ngg
if k~=ms
kk=kk+1;
del(:,kk)=180/pi*(x(:,k)-x(:,ms));
fprintf(' d(%g,',ngr(k)), fprintf('%g)', ngr(ms))
else, end
end
fprintf(' \n')
disp([t, del])
h=figure; figure(h)
plot(t, del)
title(['Phase angle difference (fault cleared at ',
num2str(tc),'s)'])
xlabel('t, sec'), ylabel('Delta, degree'), grid
resp=0;
while strcmp(resp, 'n')~=1 & strcmp(resp, 'N')~=1 &
strcmp(resp, 'y')~=1 & strcmp(resp, 'Y')~=1
resp=input('Another clearing time of fault? Enter ''y'' or
''n'' within quotes -> ');

if strcmp(resp, 'n')~=1 & strcmp(resp, 'N')~=1 & strcmp(resp,
'y')~=1 & strcmp(resp, 'Y')~=1
fprintf('\n Incorrect reply, try again \n\n'), end
end
resptc=resp;

end
resp2=0;
while strcmp(resp2, 'n')~=1 & strcmp(resp2, 'N')~=1 &
strcmp(resp2, 'y')~=1 & strcmp(resp2, 'Y')~=1
resp2=input('Another fault location: Enter ''y'' or ''n''
within quotes -> ');
if strcmp(resp2, 'n')~=1 & strcmp(resp2, 'N')~=1 &
strcmp(resp2, 'y')~=1 & strcmp(resp2, 'Y')~=1
fprintf('\n Incorrect reply, try again \n\n'), end
respf1=resp2;
end
if respf1=='n' | respf1=='N', return, else, end
end

ybusaf
% This function forms the bus admittance matrix including load
% admittances after removal of faulted line. The corresponding
reduced
% bus admittance matrix is obtained for transient stability study.
%
% Copyright (c) 1998 by H. Saadat
function [Yaf]=ybusaf(linedata, yload, nbus1,nbust, nbrt);
global Pm f H E Y th ngg

nl=linedata(:, 1); nr=linedata(:, 2);
remove = 0;
while remove ~= 1
fprintf('\nFault is cleared by opening a line. The bus to bus
nos. of the\n')
fprintf('line to be removed must be entered within brackets,
e.g. [5, 7]\n')
fline=input('Enter the bus to bus Nos. of line to be removed ->
');
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;
lfybus
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;

ybusbf
% This function forms the bus admittance matrix including load
% admittances before fault. The corresponding reduced bus
% admittance matrix is obtained for transient stability study.
%
% Copyright (c) 1998 by H. Saadat

function [Ybus, Ybf] = ybusbf(linedata, yload, nbus1, nbust)
global Pm f H E Y th ngg

lfybus
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;

Ybusdf
% This function forms the bus admittance matrix including load
% admittances during fault. The corresponding reduced bus
% admittance matrix is obtained for transient stability study.
%
% Copyright (c) 1998 by H. Saadat

function Ypf=ybusdf(Ybus, nbus1, nbust, nf)
global Pm f H E Y th ngg
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);
Ypf=YGG-YLG.'*inv(YLL)*YLG;