PROJECT REPORT ON
“ Mini Project on Selected topic ”
Course title : POWER SYSTEM-1
Course Code : EEE-305
Submission Date: 24-12-2024
Submitted by: Submitted to:
1. ALI HOSSAIN
ID: 231901013 DR.MD. AHSANUL ALAM
2. TASNIA TABASSUM ASSOCIATE PROFESSOR,
ID: 231901014 DEPARTMENT OF EEE
3. NAHID AHMED NISAT GREEN UNIVERSITY OF
BANGLADESH
ID: 231901023
4. ABU SALEH NASIM
ID: 231901003
5. DIBBO CH. ROY
ID: 231901002
Table of Contents
1. Input Data for Analysis ……………………………………………………………………….. 03
1.1 Global Variables
1.2 Bus Data
1.3 Line Data
1.4 Generator Data
2. Forming Bus Admittance Matrix …………………………………………………………. 04
2.1 Program for lfybus
3. Load Flow Solution ……………………………………………………………………………… 04- 07
3.1 Program for lfnewton
4. Forming Zbus Including the Load …………………………………………..…………… 07- 09
4.1 Program for zbuildpi
5.Symmetrical Fault Analysis ……………………………………………………………………. 10- 12
5.1 Program for symfault
6.Three-Phase Fault Results………………………………………………………………………. 12- 14
6.1 Example: Three-Phase Fault at Bus 6
2
Input data for the analysis:
clc, clear global basemva basemva = 100; accuracy = 0.001; accel =
1.8; maxiter = 100;
% IEEE 6-BUS TEST SYSTEM
% Bus Bus Voltage Angle ---Load---- -------Generator----- Static Mvar
% No code Mag. Degree MW Mvar MW Mvar Qmin Qmax +Qc/-Ql busdata=[1 1 1.06 0.0 0.0
0.0 0.0 0.0 0 0 0
2 0 1.04 0.0 0.0 10.0 150.0 30 0 0 0
3 0 1.03 0.0 0.0 15.0 100.0 0.0 0 0 0
4 0 1.0 0.0 100.0 70.0 0.0 0.0 0 0 05 0
1.0 0.0 90.0 30.0 0.0 0.0 0 0 0
6 0 1.0 0.0 160.0 110 0.0 0.0 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
1 5 0.025 0.105 0.0045 1 1 6 0.040 0.215 0.0055 1
2 4 0.000 0.035 0.0000 1
3 5 0.000 0.042 0.0000 1
4 6 0.028 0.125 0.0035 1
5 6 0.026 0.175 0.0300 1];
% Gen. Ra Xd'
gendata=[1 0 0.20
2 0 0.15
3 0 0.25];
lfybus % form the bus admittance matrix
lfnewton % Load flow solution by Newton Raphson method
Zbus=zbuildpi(linedata, gendata,yload); %Forms Zbus including the Load symfault(linedata,
Zbus, V) %3-phase fault including load current
Program For lfybus:
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));
3
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
Program for lfnewton:
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 % Test for max. power mismatch for i=1 m
for k=1 m
A(i,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
i=1 nbr
if nl(i) == n | nr(i) == n
4
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);
% 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
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
5
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
Program for zbuildpi:
function [Zbus, linedata] = zbuildpi(linedata, gendata, yload)
ng=length(gendata(:,1));
nlg=gendata(:,1); nrg=gendata(:,2);
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
6
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;
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
7
% 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
Program for Symfault & SCC:
8
function symfaul(zdata, Zbus, V) global
basemva
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('\Three-phase balanced fault analysis \n')
ff = 999; while ff > 0 nf = input('Enter
Faulted Bus No. -> ');
while nf <= 0 | nf > nbus
fprintf('Faulted bus No. must be between 1 & %g \n', nbus)
nf = input('Enter Faulted Bus No. -> '); end
fprintf('\nEnter Fault Impedance Zf = R + j*X in ') Zf =
input('complex form (for bolted fault enter 0). Zf = '); 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;
%SSC1=sqrt(3)*Ifm*abs(V0(nf)*le-3
SCC=Ifm*basemva; %Considering VL=VB; here Ifm in pu fprintf('Total fault
current = %8.4f per unit \n\n', Ifm) fprintf('SCC = %8.4f MVA \n\n',SCC)
%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')
9
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 &
10
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
Result :
1)Three phase fault at bus 6 Enter
Faulted Bus No. -> 6
Enter Fault Impedance Zf = R + j*X in complex form (for bolted fault enter 0). Zf = 0
Balanced three-phase fault at bus No. 6
Total fault current = 8.8058 per unit
SCC = 880.5796 MVA
Bus Voltages during fault in per unit
Bus Voltage Angle
No. Magnitude degrees
1 0.5071 -2.5947
2 0.5015 2.9322
3 0.4969 2.1292
4 0.4064 -1.3000
5 0.4225 -0.4474
6 0.0000 0.0000
Line currents for fault at bus No. 6
From To Current Angle
Bus Bus Magnitude degrees
1 4 0.4444 -88.5384
1 5 0.7998 -89.6085
1 6 2.3192 -81.9878
2 4 2.8768 -69.7362
3 5 1.8379 -73.6258
4 6 3.1733 -78.6491
5 6 2.3901 -81.6961 6 F 8.8058 -75.3527
2.
i)Fault at bus 4
Enter Faulted Bus No. -> 4
Enter Fault Impedance Zf = R + j*X in complex form (for bolted fault enter 0). Zf = 0
11
Balanced three-phase fault at bus No. 4
Total fault current = 10.4091 per unit
SCC = 1040.9071 MVA
Bus Voltages during fault in per unit
Bus Voltage Angle
No. Magnitude degrees
1 0.5509 1.0136
2 0.1709 3.7151
3 0.6063 4.0775
4 0.0000 0.0000
5 0.5523 2.0596
6 0.3330 3.3196
Line currents for fault at bus No. 4
From To Current Angle
Bus Bus Magnitude degrees
1 4 2.4200 -80.0608
1 6 1.0001 -81.7883
2 4 4.8833 -86.2849
3 5 1.3745 -66.2338
4 F 10.4091 -76.4007
5 1 0.0967 6.8873
5 6 1.2429 -80.6421
6 4 2.5997 -74.0295
ii)Fault at Bus 5
Enter Faulted Bus No. -> 5
Enter Fault Impedance Zf = R + j*X in complex form (for bolted fault enter 0). Zf = 0
Balanced three-phase fault at bus No. 5
Total fault current = 9.8606 per unit
SCC = 986.0582 MVA
Bus Voltages during fault in per unit
Bus Voltage Angle
No. Magnitude degrees
1 0.4452 -2.0027
2 0.6243 5.3388
3 0.1403 1.8307
4 0.5576 2.1232
5 0.0000 0.0000
12
6 0.3838 2.7534
Line currents for fault at bus No. 5
From To Current Angle
Bus Bus Magnitude degrees
1 5 4.1255 -78.5832
2 4 2.1280 -59.8286 3 5 3.3402 -
88.1693
4 1 0.5197 -62.9169
4 6 1.3582 -76.5608
5 F 9.8606 -76.9089
6 1 0.3226 71.3080 6 5 2.1709 -78.4953
3.Comparing SCC:
SCC at bus 4 1040.9071 MVA
SCC at bus 5 986.0582 MVA
SCC at bus 6 880.5796 MVA
THEEND
13
