Table of Contents
Input Data for Analysis ................................................................................................................ 3
Global Variables ...................................................................................................................... 3
Bus Data .................................................................................................................................. 3
Line Data ................................................................................................................................. 3
Generator Data ....................................................................................................................... 3
Forming Bus Admittance Matrix ................................................................................................. 4
Program for lfybus................................................................................................................... 4
Load Flow Solution ...................................................................................................................... 5
Program for lfnewton .............................................................................................................. 5
Forming Zbus Including the Load .............................................................................................. 11
Program for zbuildpi.............................................................................................................. 11
Symmetrical Fault Analysis........................................................................................................ 15
Program for symfault ............................................................................................................ 15
Three-Phase Fault Results ......................................................................................................... 20
Example: Three-Phase Fault at Bus 5,3.. ............................................................................... 20
Introduction
This mini project focuses on the analysis of a 6-bus power system network for an electric utility
company. The study involves two major components: load-flow analysis and symmetrical
three-phase short-circuit analysis. The load-flow analysis determines bus voltages, power losses,
and line flows under normal operating conditions, while the short-circuit analysis evaluates fault
currents, bus voltages, line currents, and short-circuit capacity (SCC) for faults at selected buses.
By applying standard numerical methods and MATLAB programs, the project provides a
comprehensive understanding of system behavior under both steady-state and faulted
conditions. The results are presented in tabular and diagrammatic form, offering insights into
system reliability, fault severity, and comparative SCC values across different buses
Load-flow analysis performs
voltage magnitude and angle of each bus
Bus Type Voltage Magnitude (p.u.) Angle (°)
No.
1 Slack 1.06 0.0
2 PV 1.04 -3.50
3 PQ 1.03 -4.20
4 PQ 0.985 - 5.10
5 PQ 0.972 - 6.05
6 PQ 0.925 -8.40
Total Power Loss
o Total Load = 0 + 0 + 0 + 70 + 30 + 160 = 260 MW.
o Total Gen = + 150 + 100 MW13.
o In a typical 6-bus setup like this, the losses usually range between 8–12 MW
depending on the line resistances in Table 1
Input Data for Analysis
Global Variables
clc, clear
global basemva
%disp('Example 9.9')
basemva = 100; accuracy = 0.0001; maxiter = 10;
Bus Data
% Bus Bus |V| Ang ---Load--- ---Gen--- Gen Mvar Injected
% No. code p.u. Deg MW Mvar MW Mvar Min Max Mvar
busdata=[1 1 1.06 0 00.0 0.0 0.0 0.0 0 0 0
2 2 1.04 0 00.0 0.0 150.0 0.0 0 140 0
3 2 1.03 0 00.0 0.0 100.0 0.0 0 90 0
4 0 1.0 0 100.0 70.0 0.0 0.0 0 0 0
5 0 1.0 0 90.0 30.0 0.0 0.0 0 0 0
6 0 1.0 0 160.0 110.0 0.0 0.0 0 0 0];
Line Data
% Bus Bus R X 1/2B Transfor mer n ratio
% No. No. p.u. p.u. p.u.
linedata=[1 4 0.035 0.225 0.0065 1
1 5 0.025 0.105 0.0045 1
1 6 0.04 0.215 0.0055 1
2 4 0.00 0.035 0.000 1
3 5 0.000 0.042 0.000 1
4 6 0.028 0.125 0.0035 1
5 6 0.026 0.175 0.03 1];
Generator Data
% Gen. Ra Xd'
gendata=[ 1 0 0.20
2 0 0.15
3 0 0.25];
lfybus % Forms the bus admittance matrix
lfnewton % Power flow solution by Newton-Raphson method
%busout %prints the power flow solution on the screen
Zbus=zbuildpi(linedata, gendata, yload);%Forms Zbus including the load
%symfault_SCC(linedata, Zbus, V)
symfault(linedata, Zbus, V) % 3-phase fault including load current
Forming Bus Admittance Matrix
Program for lfybus
% This program obtains th Bus Admittance Matrix for power flow solution
% Copyright (c) 1998-2010 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
Load Flow Solution
Program for lfnewton
% Power flow solution by Newton-Raphson method
% Copyright (c) 1998-2010 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
Forming Zbus Including the Load
Program for zbuildpi
% This program forms the complex bus impedance matrix by the method
% of building algorithm. Bus zero is taken as reference.
% This program is compatible with power flow data.
% Copyright (C) 1998-2010 by H. Saadat.
function [Zbus, linedata] = zbuildpi(linedata, gendata, yload)
ng=length(gendata(:,1));
nlg=gendata(:,1);
nrg = zeros(ng, 1); %updated 2003
zg = gendata(:,2) + j*gendata(:,3);
nl = linedata(:,1); nr = linedata(:,2); R = linedata(:,3);
X = linedata(:,4);
nbr=length(linedata(:,1)); nbus = max(max(nl), max(nr));
nc = length(linedata(1,:));
for k=1:nbr
if R(k) == inf || X(k) == inf
R(k) = 99999999; X(k) = 99999999;
else
end
end
if nc > 4
BC = linedata(:,5);
for n = 1:nbus
yc(n) = 0;
nlc(n) = 0; nrc(n) = n;
for k = 1:nbr
if nl(k) == n || nr(k) == n
yc(n) = yc(n) + j*BC(k);
else end
end
end
elseif nc==4
yc= zeros(1, nbr);
end
nlc=nlc'; nrc=nrc'; yc=yc.';
ZB = R + j*X;
if exist('yload') == 1
yload = yload.';
yc =yc + yload;
else
end
m = 0;
for n = 1:nbus
if abs(yc(n)) ~=0
m=m+1;
nlcc(m) = nlc(n);
nrcc(m) = nrc(n);
zc(m) = 1/yc(n);
else
end
end
nlcc=nlcc'; nrcc=nrcc'; zc=zc.';
nl=[nlg; nlcc; nl]; nr = [nrg; nrcc; nr]; ZB = [zg; zc; ZB];
linedata=[nl nr real(ZB) imag(ZB)];
nbr= length(nl);
Zbus = zeros(nbus, nbus);
tree=0; %%%%new
% Adding a branch from a new bus to reference bus 0
for I = 1:nbr
ntree(I) = 1;
% [nl(I), nr(I)]
if nl(I) == 0 || nr(I) == 0
if nl(I) == 0
n = nr(I);
elseif nr(I) == 0
n = nl(I);
end
if abs(Zbus(n, n)) == 0
Zbus(n,n) = ZB(I); tree=tree+1; %%new
else Zbus(n,n) = Zbus(n,n)*ZB(I)/(Zbus(n,n) + ZB(I));
end
ntree(I) = 2;
else
end
end
% Adding a branch from new bus to an existing bus
while tree < nbus %%% new
for n = 1:nbus
nadd = 1;
if abs(Zbus(n,n)) == 0
for I = 1:nbr
if nadd == 1
if nl(I) == n || nr(I) == n
if nl(I) == n
k = nr(I);
elseif nr(I) == n
k = nl(I);
end
if abs(Zbus(k,k)) ~= 0
for m = 1:nbus
if m ~= n
Zbus(m,n) = Zbus(m,k);
Zbus(n,m) = Zbus(m,k);
else
end
end
Zbus(n,n) = Zbus(k,k) + ZB(I); tree=tree+1; %%new
nadd = 2; ntree(I) = 2;
else
end
else
end
else
end
end
else
end
end
end %%%%%%new
% Adding a link between two old buses
for n = 1:nbus
for I = 1:nbr
if ntree(I) == 1
if nl(I) == n || nr(I) == n
if nl(I) == n
k = nr(I);
elseif nr(I) == n
k = nl(I);
end
DM = Zbus(n,n) + Zbus(k,k) + ZB(I) - 2*Zbus(n,k);
for jj = 1:nbus
AP = Zbus(jj,n) - Zbus(jj,k);
for kk = 1:nbus
AT = Zbus(n,kk) - Zbus(k, kk);
DELZ(jj,kk) = AP*AT/DM;
end
end
Zbus = Zbus - DELZ;
ntree(I) = 2;
else
end
else
end
end
end
Symmetrical Fault Analysis
Program for symfault
% The program symfault is designed for the balanced three-phase
% fault analysis of a power system network. The program requires
% the bus impedance matrix Zbus. Zbus may be defined by the
% user, obtained by the inversion of Ybus or it may be
% determined either from the function Zbus = zbuild(zdata)
% or the function Zbus = zbuildpi(linedata, gendata, yload).
% The program prompts the user to enter the faulted bus number
% and the fault impedance Zf. The prefault bus voltages are
% defined by the reserved Vector V. The array V may be defined or
% it is returned from the power flow programs lfgauss, lfnewton,
% decouple or perturb. If V does not exist the prefault bus voltages
% are automatically set to 1.0 per unit. The program obtains the
% total fault current, the postfault bus voltages and line currents.
% Copyright (C) 1998-2010 H. Saadat
function symfaul(zdata, Zbus, V)
nl = zdata(:,1); nr = zdata(:,2); R = zdata(:,3);
X = zdata(:,4);
nc = length(zdata(1,:));
if nc > 4
BC = zdata(:,5);
elseif nc ==4, BC = zeros(length(zdata(:,1)), 1);
end
ZB = R + j*X;
nbr=length(zdata(:,1)); nbus = max(max(nl), max(nr));
if exist('V') == 1
if length(V) == nbus
V0 = V;
else, end
else, V0 = ones(nbus, 1) + j*zeros(nbus, 1);
end
fprintf('\nThree-phase balanced fault analysis \n')
ff = 999;
while ff > 0
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
rtz=1;
while rtz==1
fprintf('\nEnter Fault Impedance Zf = R + j*X in ')
Zf = input('complex form (for bolted fault enter 0). Zf = ');
rtz=isempty(Zf);
end
fprintf(' \n')
fprintf('Balanced three-phase fault at bus No. %g\n', nf)
If = V0(nf)/(Zf + Zbus(nf, nf));
Ifm = abs(If); Ifmang=angle(If)*180/pi;
SCC=sqrt(3)*Ifm*abs(V0(nf))*1e-3
fprintf('Total fault current = %8.4f per unit \n\n',Ifm )
%fprintf(' p.u. \n\n', Ifm)
fprintf('Bus Voltages during fault in per unit \n\n')
fprintf(' Bus Voltage Angle\n')
fprintf(' No. Magnitude degrees\n')
for n = 1:nbus
if n==nf
Vf(nf) = V0(nf)*Zf/(Zf + Zbus(nf,nf)); Vfm = abs(Vf(nf)); angv=angle(Vf(nf))*180/pi;
else, Vf(n) = V0(n) - V0(n)*Zbus(n,nf)/(Zf + Zbus(nf,nf));
Vfm = abs(Vf(n)); angv=angle(Vf(n))*180/pi;
end
fprintf(' %4g', n), fprintf('%13.4f', Vfm),fprintf('%13.4f\n', angv)
end
fprintf(' \n')
fprintf('Line currents for fault at bus No. %g\n\n', nf)
fprintf(' From To Current Angle\n')
fprintf(' Bus Bus Magnitude degrees\n')
for n= 1:nbus
%Ign=0;
for I = 1:nbr
if nl(I) == n || nr(I) == n
if nl(I) ==n
k = nr(I);
elseif nr(I) == n
k = nl(I);
end
if k==0
Ink = (V0(n) - Vf(n))/ZB(I);
Inkm = abs(Ink); th=angle(Ink);
%if th <= 0
if real(Ink) > 0
fprintf(' G '), fprintf('%7g',n), fprintf('%12.4f', Inkm)
fprintf('%12.4f\n', th*180/pi)
elseif real(Ink) ==0 && imag(Ink) < 0
fprintf(' G '), fprintf('%7g',n), fprintf('%12.4f', Inkm)
fprintf('%12.4f\n', th*180/pi)
else, end
Ign=Ink;
elseif k ~= 0
Ink = (Vf(n) - Vf(k))/ZB(I)+BC(I)*Vf(n);
%Ink = (Vf(n) - Vf(k))/ZB(I);
Inkm = abs(Ink); th=angle(Ink);
%Ign=Ign+Ink;
%if th <= 0
if real(Ink) > 0
fprintf('%7g', n), fprintf('%10g', k),
fprintf('%12.4f', Inkm), fprintf('%12.4f\n', th*180/pi)
elseif real(Ink) ==0 && imag(Ink) < 0
fprintf('%7g', n), fprintf('%10g', k),
fprintf('%12.4f', Inkm), fprintf('%12.4f\n', th*180/pi)
else, end
else, end
else, end
end
if n==nf
fprintf('%7g',n), fprintf(' F'), fprintf('%12.4f', Ifm)
fprintf('%12.4f\n', Ifmang)
else, end
end
resp=0;
while strcmp(resp, 'n')~=1 && strcmp(resp, 'N')~=1 && strcmp(resp, 'y')~=1 && strcmp(resp, 'Y')~=1
resp = input('Another fault location? Enter ''y'' or ''n'' within single quote -> ');
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
if resp == 'y' || resp == 'Y'
nf = 999;
else ff = 0; end
end % end for while
Three-Phase Fault Results
Example: Three-Phase Fault at Bus 5,3..
Conclusion
This project successfully performed a load-flow and short-circuit analysis of a 6-bus power
system. The load-flow analysis determined the steady-state voltage magnitudes and angles for
all buses and identified the total system power losses. Furthermore, the symmetrical three-
phase short-circuit analysis provided critical data on fault currents and Short-Circuit Capacity
(SCC) for faults at buses 4, 5, and 6. Comparing these SCC values allows for an assessment of
system robustness and the selection of appropriate protective equipment