13BK1D0XX

SPEC

EXPERIMENT-1
AIM: To Develop MATLAB program for YBUS formation.

SOFTWARE-REQUIRED:
MATLAB:7.8
THEORY:
Bus admittance is often used in power system studies. In most of the power
system studies it is required to form y- bus matrix of the system by considering
certain power system parameters depending upon the type of analysis. Y-bus may be
formed by inspection method only if there is no mutual coupling between the lines. Every
transmission line should be represented by p- equivalent. Shunt impedances are added to
diagonal element corresponding to the buses at which these are connected. The off
diagonal elements are unaffected. The equivalent circuit of Tap changing transformers is
included while forming Y-bus matrix.
FORMATION OF Y-BUS MATRIX:
Generalized Y-bus = yii .. yid
ydi ydd
where, Yii = Self admittance
Ydi = Transfer admittance
PROCEDURE:
1. Enter the command window of the MATLAB.
2. Create a new M file by selecting File - New M File
3. Type and save the program in the editor window.
4. Execute the program by either pressing tools Run.
5. View the results.

1
M.TECH

EPS

13BK1D0XX

SPEC

PROGRAM:
% loine data
%
bus
%
nl
linedata=[1
1
2

bus
nr
2
3
3

R
pu
0.02
0.01
0.0125

X
pu
0.04
0.03
0.025

1/2B
pu
0
0
0

Tap
1
1
1];

This program obtains the Bus Admittance Matrix

for power flow solution
%

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
disp(Ybus);
clear Pgg

2
M.TECH

EPS

13BK1D0XX

SPEC

RESULTS
Ybus =[20.0000 -50.0000i -10.0000 +20.0000i -10.0000 +30.0000i
-10.0000 +20.0000i 26.0000 -52.0000i -16.0000 +32.0000i
-10.0000 +30.0000i -16.0000 +32.0000i 26.0000 -62.0000i]

3
M.TECH

EPS

13BK1D0XX

SPEC

EXPERIMENT-2
AIM : To Develop MATLAB program for G-S Load Flow Analysis .
SOFTWARE REQUIRED:
MATLAB -7.8
THEORY:
The GAUSS SEIDEL method is an iterative algorithm for solving a set of non-linear
load flow equations.
PROCEDURE:
1. Enter the command window of the MATLAB.
2. Create a new M file by selecting File - New M File
3. Type and save the program in the
editor Window
4. Execute the program by either pressing tool Run.
5. View the results.
PROGRAM:

%PF studies for the Gauss Siedel Method

clc;
clear all;
basemva=100;
accuracy=0.001;
maxiter=100;
%bus data
%
Bus
Inj
%
No
busdata=[1
0
2
0
3
0];

Bus

Vmag

Angle

-Load-

-Gen-

Qmin

Qmax

Code
1

1.05

deg
0

MW
0

MVAr
0

MW
0

MVAr
0

MVAr
0

MVAr
0

1.00

400

250

1.04

200

-100

R
pu
0.02
0.01
0.0125

X
pu
0.04
0.03
0.025

1/2B
pu
0
0
0

200

% line data
%
bus
%
nl
linedata=[1
1
2
Lfybus;

bus
nr
2
3
3

Tap
1
1
1];

4
M.TECH

EPS

13BK1D0XX

SPEC

Lfgauss;
Busout;
Lineflow;

MATLAB PROGRAM FOR GUASS SEIDAL METHOD

%

Power flow solution by Gauss-Seidel method

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
DV(n)=0;
end
num = 0; AcurBus = 0; converge = 1;
Vc = zeros(nbus,1)+j*zeros(nbus,1); Sc = zeros(nbus,1)+j*zeros(nbus,1);
while exist('accel')~=1
accel = 1.3;
end
while exist('accuracy')~=1
accuracy = 0.001;
end
while exist('basemva')~=1
basemva= 100;
end
while exist('maxiter')~=1
maxiter = 100;
end
iter=0;
maxerror=10;
while maxerror >= accuracy & iter <= maxiter
iter=iter+1;
for n = 1:nbus;
YV = 0+j*0;
for L = 1:nbr;
if nl(L) == n, k=nr(L);
YV = YV + Ybus(n,k)*V(k);
elseif nr(L) == n, k=nl(L);
YV = YV + Ybus(n,k)*V(k);
end
end
Sc = conj(V(n))*(Ybus(n,n)*V(n) + YV) ;

5
M.TECH

EPS

13BK1D0XX

SPEC

Sc = conj(Sc);
DP(n) = P(n) - real(Sc);
DQ(n) = Q(n) - imag(Sc);
if kb(n) == 1
S(n) =Sc; P(n) = real(Sc); Q(n) = imag(Sc); DP(n) =0; DQ(n)=0;
Vc(n) = V(n);
elseif kb(n) == 2
Q(n) = imag(Sc); S(n) = P(n) + j*Q(n);

buses are

if Qmax(n) ~= 0
Qgc = Q(n)*basemva + Qd(n) - Qsh(n);
if abs(DQ(n)) <= .005 & iter >= 10 % After 10 iterations
if DV(n) <= 0.045
% the Mvar of generator

limits Vm(n)
0.005 pu
to bring
within the

if Qgc < Qmin(n),

% tested. If not within

Vm(n) = Vm(n) + 0.005;

% is changed in steps of

DV(n) = DV(n)+.005;

% up to .05

elseif Qgc > Qmax(n),

% the generator Mvar

pu in order

Vm(n) = Vm(n) - 0.005;

% specified limits.
DV(n)=DV(n)+.005; end
else, end
else,end
else,end
end
if kb(n) ~= 1
Vc(n) = (conj(S(n))/conj(V(n)) - YV )/ Ybus(n,n);
else, end
if kb(n) == 0
V(n) = V(n) + accel*(Vc(n)-V(n));
elseif kb(n) == 2
VcI = imag(Vc(n));
VcR = sqrt(Vm(n)^2 - VcI^2);
Vc(n) = VcR + j*VcI;
V(n) = V(n) + accel*(Vc(n) -V(n));
end

end
maxerror=max( max(abs(real(DP))), max(abs(imag(DQ))) );
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= ('
else,
tech=('
Method');
end
k=0;
for n = 1:nbus

ITERATIVE SOLUTION DID NOT CONVERGE');

Power Flow Solution by Gauss-Seidel

6
M.TECH

EPS

13BK1D0XX

SPEC

Vm(n) = abs(V(n)); deltad(n) = angle(V(n))*180/pi;

if kb(n) == 1
S(n)=P(n)+j*Q(n);
Pg(n) = P(n)*basemva + Pd(n);
Qg(n) = Q(n)*basemva + Qd(n) - Qsh(n);
k=k+1;
Pgg(k)=Pg(n);
elseif kb(n) ==2
k=k+1;
Pgg(k)=Pg(n);
S(n)=P(n)+j*Q(n);
Qg(n) = Q(n)*basemva + Qd(n) - Qsh(n);
end
yload(n) = (Pd(n)- j*Qd(n)+j*Qsh(n))/(basemva*Vm(n)^2);
end
Pgt = sum(Pg); Qgt = sum(Qg); Pdt = sum(Pd); Qdt = sum(Qd); Qsht =
sum(Qsh);
busdata(:,3)=Vm'; busdata(:,4)=deltad';
clear AcurBus DP DQ DV L Sc Vc VcI VcR YV converge delta

7
M.TECH

EPS

13BK1D0XX

SPEC

RESULTS
Ybus=

20.0000 -50.0000i -10.0000 +20.0000i -10.0000

+30.0000i
-10.0000 +20.0000i 26.0000 -52.0000i -16.0000
+32.0000i
-10.0000 +30.0000i -16.0000 +32.0000i 26.0000
-62.0000i
Power Flow Solution by Gauss-Seidel
Method
Maximum Power Mismatch = 7.88577e-005
No. of Iterations = 6
Bus Voltage
---Generation--No. Mag.
Mvar
Mvar
1
140.871
2
0.000
3
146.171

Angle
------Load-----Injected
Degree
MW
Mvar

1.050
0.000
0.000
0.972
-2.697
0.000
1.040
-0.499
0.000

Total
287.042

MW

0.000

0.000

218.436

400.000

250.000

0.000

0.000

0.000

200.000

400.000

250.000

418.436

0.000
Line Flow and Losses

--Line-Transformer
from to
tap

--Line loss--

MW

Mvar

MVA

218.436
179.366
39.064

140.871
118.726
22.117

259.921
215.100
44.890

8.393
0.183

16.787
0.548

-400.000 -250.000
1 -170.972 -101.940

471.699
199.056

8.393

16.787

1
2
3
2

Power at bus & line flow

MW

Mvar

8
M.TECH

EPS

13BK1D0XX

SPEC
3 -229.035 -148.043

3
1
2

200.000
-38.881
238.881

146.171
-21.568
167.736

272.715

9.846

19.693

247.722
44.463
291.890

0.183
9.846

0.548
19.693

18.423

37.028

Total loss

EXAMPLE
%PF studies for the Gauss Siedel Method 5-bus
system
clear all
clc
basemva=100;
accuracy = .0001;
maxiter = 10;
% IEEE 5-bus system
% bus bus
% no code
busdata=[1
1
2
2
3
2
4
0
5
0

vaolt angle ----load-- ---generator-injected

mag
degree MW MVAR
MW MVAR Qmin Qmax MVAR
1.06
0
0
0
0 0
10 50
0
1.045
0
20
10
40 30 10 50
0
1.03
0
20
15
30 10 10 40
0
1
0
50
30
0
0
0 0
0
1
0
60
40
0
0
0 0
0];

%line dadta
%
busnl busnr R X 1/2B tapsetting
linedata=[ 1
2
0.02
0.06
0.03
1
3
0.08
0.24
0.025
2
3
0.06
0.18
0.020
2
4
0.06
0.018
0.023
2
5
0.04
0.12
0.015
3
4
0.01
0.03
0.010
4
5
0.08
0.24
0.025
Lfybus;
Decouple;
Busout;
Lineflow;

1
1
1
1
1
1
1];

9
M.TECH

EPS

13BK1D0XX

SPEC

EXPERIMENT-3
AIM: To Develop MATLAB program for N-R Load Flow Analysis.
SOFTWARE REQUIRED:
MATLAB:7.8
THEORY:
The N-R method of load flow analysis is an iterative method which
approximates the set of nonlinear simultaneous equations to a set of linear
simultaneous equations
using Taylors series and the terms are limited to first order approximation.
The load flow equations for Newton Raphson method are non-linear equations in terms
of real and imaginary part of bus voltages.

PROGRAM:

%PF studies for the Newton raphson

Method

clc;
clear all;
basemva=100;
accuracy=0.001;
maxiter=100;
%bus data
%
Bus
Inj
%
No
busdata=[1
0

Bus

Vmag

Angle

Code
1

1.05

deg
0

-LoadMW
0

MVAr
0

-GenMW
0

MVAr
0

Qmin

Qmax

MVAr
0

MVAr
0

10
M.TECH

EPS

13BK1D0XX

0
0];

SPEC

1.00

400

1.04

R
pu
0.02
0.01
0.0125

X
pu
0.04
0.03
0.025

1/2B
pu
0
0
0

250

200

-100

0
200

% line data
%
bus
%
nl
linedata=[1
1
2
Lfybus;
Lfnewton;
Busout;
Lineflow;

bus
nr
2
3
3

Tap
1
1
1];

MATLAB PROGRAM FOR Newton-Raphson METHOD

%

Power flow solution by Newton-Raphson method

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

11
M.TECH

EPS

13BK1D0XX

SPEC

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

12
M.TECH

EPS

13BK1D0XX

SPEC

J2

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

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

13
M.TECH

EPS

13BK1D0XX

SPEC

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
%clear A DC DX

Ybus=

J11 J22 J33 J44 Qk delta lk ll lm

J11 J22 J33 Qk delta lk ll lm

RESULTS
0.0000 -50.0000i -10.0000 +20.0000i -10.0000 +30.0000i

14
M.TECH

EPS

13BK1D0XX

SPEC

-10.0000 +20.0000i 26.0000 -52.0000i -16.0000 +32.0000i

-10.0000 +30.0000i -16.0000 +32.0000i 26.0000 -62.0000i
Power Flow Solution by Newton-Raphson Method
Maximum Power Mismatch = 0.000216612
No. of Iterations = 3
Bus Voltage Angle
No. Mag.
1 1.050

Degree
0.000

------Load-----MW
0.000

---Generation--- Injected

Mvar

MW

Total

Mvar

0.000 218.403 140.848

2 0.972 -2.696 400.000 250.000

3 1.040 -0.499

Mvar

0.000

0.000

0.000

0.000

0.000

0.000 200.000 146.161

0.000

400.000 250.000 418.403 287.010

0.000

Line Flow and Losses

--Line-- Power at bus & line flow
from to
1

MW

Mvar

--Line loss-- Transformer

MVA

MW

Mvar

tap

218.403 140.848 259.881

2 179.362 118.734 215.101
3 39.061 22.118 44.888

8.393 16.787
0.183

0.548

-400.000 -250.000 471.699

1 -170.968 -101.947 199.056

8.393 16.787

3 -229.032 -148.053 272.718

9.847 19.693

200.000 146.161 247.716

1 -38.878 -21.569 44.461
2 238.878 167.746 291.893

Total loss

0.183

0.548

9.847 19.693

18.423 37.028

EXPERIMENT-4
AIM: To Develop MATLAB program for FDLF Load Flow Analysis.

15
M.TECH

EPS

13BK1D0XX

SPEC

SOFTWARE REQUIRED:
MATLAB -7.8
THEORY:
The GAUSS SEIDEL method is an iterative algorithm for solving a set of non-linear
load flow equations.
PROCEDURE:
1. Enter the command window of the MATLAB.
2. Create a new M file by selecting File - New M File
3. Type and save the program in the
editor Window
4. Execute the program by either pressing tool Run.
5. View the results.
PROGRAM:

PF studies for the fast decouple load

flow Method
%

clc;
clear all;
basemva=100;
accuracy=0.001;
maxiter=100;
%bus data
%
Bus
----Inj--%
No
busdata=[1
2
3

Bus

Volt

Angle

Code
1
0
2

mag
1.05
1.00
1.04

deg
0
0
0

-LoadMW MVAr
0
0
400 250
0
0

-----Gen----MW
MVAr Qmin Qmax MVAr
0
0
0
0
0
0
0
0
0
0
200 0
-100
200 0];

% loine data
%
bus
%
nl
linedata=[1
1
2
Lfybus;
decouple;
Busout;
Lineflow;

bus
nr
2
3
3

R
pu
0.02
0.01
0.0125

X
pu
0.04
0.03
0.025

1/2B
pu
0
0
0

Tap
1
1
1];

MATLAB PROGRAM FOR Fast Decoupled METHOD

%

Fast Decoupled Power Flow Solution

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

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

16
M.TECH

EPS

13BK1D0XX

SPEC

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
if kb(n) == 1, ns = ns+1; else, end
nss(n) = ns;
end
Ym = abs(Ybus); t = angle(Ybus);
ii=0;
for ib=1:nbus
if kb(ib) == 0 | kb(ib) == 2
ii = ii+1;
jj=0;
for jb=1:nbus
if kb(jb) == 0 | kb(jb) == 2
jj = jj+1;
B1(ii,jj)=imag(Ybus(ib,jb));
else,end
end
else, end
end
ii=0;
for ib=1:nbus
if kb(ib) == 0
ii = ii+1;
jj=0;
for jb=1:nbus
if kb(jb) == 0
jj = jj+1;
B2(ii,jj)=imag(Ybus(ib,jb));
else,end
end
else, end
end
B1inv=inv(B1); B2inv = inv(B2);
maxerror = 1; converge = 1;
iter = 0;
% Start of iterations
while maxerror >= accuracy & iter <= maxiter % Test for max. power
mismatch
iter = iter+1;
id=0; iv=0;
for n=1:nbus
nn=n-nss(n);
J11=0; J33=0;

17
M.TECH

EPS

13BK1D0XX

SPEC

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));
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;
Qgc = Q(n)*basemva + Qd(n) - Qsh(n);
if Qmax(n) ~= 0
if iter <= 20
% Between the 1th & 6th
iterations
if iter >= 10
% the Mvar of generator buses
are
if Qgc < Qmin(n),
% tested. If not within limits
Vm(n)
Vm(n) = Vm(n) + 0.005;
% is changed in steps of 0.05
pu to
elseif Qgc > Qmax(n),
% bring the generator Mvar
within
Vm(n) = Vm(n) - 0.005;end % the specified limits.
else, end
else,end
else,end
end
if kb(n) ~= 1
id = id+1;
DP(id) = P(n)-Pk;
DPV(id) = (P(n)-Pk)/Vm(n);
end
if kb(n) == 0
iv=iv+1;
DQ(iv) = Q(n)-Qk;
DQV(iv) = (Q(n)-Qk)/Vm(n);
end
end
Dd=-B1\DPV';
DV=-B2\DQV';
id=0;iv=0;
for n=1:nbus
if kb(n) ~= 1
id = id+1;
delta(n) = delta(n)+Dd(id); end
if kb(n) == 0
iv = iv+1;
Vm(n)=Vm(n)+DV(iv); end
end
maxerror=max(max(abs(DP)),max(abs(DQ)));
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')

18
M.TECH

EPS

13BK1D0XX

SPEC

converge = 0; pause, else, end

end
if converge ~= 1
tech= ('
ITERATIVE SOLUTION DID NOT CONVERGE');
else,
tech=('
Power Flow Solution by Fast Decoupled
Method');
end
k=0;
V = Vm.*cos(delta)+j*Vm.*sin(delta);
deltad=180/pi*delta;
clear A; clear DC; clear DX
i=sqrt(-1);
for n = 1:nbus
if kb(n) == 1
S(n)=P(n)+j*Q(n);
Pg(n) = P(n)*basemva + Pd(n);
Qg(n) = Q(n)*basemva + Qd(n) - Qsh(n);
k=k+1;
Pgg(k)=Pg(n);
elseif kb(n) ==2
S(n)=P(n)+j*Q(n);
Qg(n) = Q(n)*basemva + Qd(n) - Qsh(n);
k=k+1;
Pgg(k)=Pg(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);
clear Pk Qk DP DQ J11 J33 B1 B1inv B2 B2inv DPV DQV Dd delta ib id ii
iv jb jj

RESULTS

Ybus=

20.0000 -50.0000i -10.0000 +20.0000i -10.0000 +30.0000i

-10.0000 +20.0000i 26.0000 -52.0000i -16.0000 +32.0000i
-10.0000 +30.0000i -16.0000 +32.0000i 26.0000 -62.0000i
Power Flow Solution by Fast Decoupled Method

19
M.TECH

EPS

13BK1D0XX

SPEC
Maximum Power Mismatch = 0.00045882
No. of Iterations = 13

Bus Voltage Angle

No. Mag.
1 1.050

------Load------

Degree
0.000

MW
0.000

---Generation--- Injected

Mvar

MW

0.000 218.408 140.854

2 0.972 -2.697 400.000 250.000

3 1.040 -0.499
Total

Mvar

0.000

0.000

Mvar
0.000

0.000

0.000

0.000 200.000 146.185

0.000

400.000 250.000 418.408 287.039

0.000

Line Flow and Losses

--Line-- Power at bus & line flow

--Line loss-- Transformer

from to

MW

MW

Mvar

MVA

Mvar

tap

218.408 140.854 259.888

2 179.359 118.723 215.093
3 39.061 22.117 44.888

8.393 16.785
0.183

0.548

-400.000 -250.000 471.699

1 -170.967 -101.937 199.050

8.393 16.785

3 -229.028 -148.037 272.707

9.846 19.692

200.000 146.185 247.730

1 -38.879 -21.569 44.461
2 238.874 167.729 291.880

Total loss

0.183

0.548

9.846 19.692

18.421 37.025

EXPERIMENT-5
AIM:
To become familiar with modelling and analysis of power systems under faulted
condition and to compute the fault level, post-fault voltages and currents for different
types of faults, both symmetric and unsymmetric.
SOFTWARE REQUIRED:
MATLAB-7.8
20
M.TECH

EPS

13BK1D0XX

SPEC

EXERCISE:

The one line diagram of a simple power system is shown in figure. The neutral of each
generator is grounded through a current limiting reactor of 0.25/3 per unit on a 100MVA
base.The system data expressed in per unit on a common 100 MVA base is tabulated
below. The generaters are running on no load at their rated voltage and rated frequency
with their emfs inphase.
Determine the fault current for the following faults.
(a) A balanced three phase fault at bus 3 through a fault impedance Zf = j0.1 per unit.
(b) A single line to ground fualt at bus3 through a fault impedance Zf = j0.1 per unit.
(c) A line to line fault at bus3 through a fault impedance Zf = j0.1 per unit.
(d) A double line to ground fault at bus3 through a fault impedance Zf = j0.1 per unit.

21
M.TECH

EPS

13BK1D0XX

SPEC

\
22
M.TECH

EPS

13BK1D0XX

SPEC

PROGRAM:
zdata1 = [0 1 0 0.25
0 2 0 0.25
1 2 0 0.125
1 3 0 0.15
2 3 0 0.25];
zdata0 = [0 1 0 0.40
0 2 0 0.10
1 2 0 0.30
1 3 0 0.35
2 3 0 0.7125];
zdata2 = zdata1;
Zbus1 = zbuild(zdata1)
Zbus0 = zbuild(zdata0)
Zbus2 = Zbus1;
symfault(zdata1,Zbus1)
lgfault(zdata0, Zbus0, zdata1, Zbus1, zdata2, Zbus2)
llfault(zdata1, Zbus1, zdata2, Zbus2)
dlgfault(zdata0, Zbus0, zdata1, Zbus1, zdata2, Zbus2)

23
M.TECH

EPS

13BK1D0XX

SPEC

Programme for symmetrical fault calculations

%
%
%
%
%
%
%
%
%
%
%
%
%
%

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.

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('\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;
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

24
M.TECH

EPS

13BK1D0XX

SPEC

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

25
M.TECH

EPS

13BK1D0XX

SPEC

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

26
M.TECH

EPS

13BK1D0XX

SPEC

Programme for L-G fault calculations

%
%
%
%
%
%
%
%
%
%
%
%
%
%
%

The program lgfault is designed for the single line-to-ground

fault analysis of a power system network. The program requires
the positive-, negative- and zero-sequence bus impedance matrices,
Zbus1 Zbus2,and Zbus0.The bus impedances matrices 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, bus voltages and line currents during the fault.

function lgfault(zdata0, Zbus0, zdata1, Zbus1, zdata2, Zbus2, V)

if exist('zdata2') ~= 1
zdata2=zdata1;
else, end
if exist('Zbus2') ~= 1
Zbus2=Zbus1;
else, end
nl = zdata1(:,1); nr = zdata1(:,2);
nl0 = zdata0(:,1); nr0 = zdata0(:,2);
nbr=length(zdata1(:,1)); nbus = max(max(nl), max(nr));
nbr0=length(zdata0(:,1));
R0 = zdata0(:,3); X0 = zdata0(:,4);
R1 = zdata1(:,3); X1 = zdata1(:,4);
R2 = zdata1(:,3); X2 = zdata1(:,4);
for k=1:nbr0
if R0(k)==inf | X0(k) ==inf
R0(k) = 99999999; X0(k) = 99999999;
else, end
end
ZB1 = R1 + j*X1; ZB0 = R0 + j*X0;
ZB2 = R2 + j*X2;
if exist('V') == 1
if length(V) == nbus
V0 = V;
else, end
else, V0 = ones(nbus, 1) + j*zeros(nbus, 1);
end
fprintf('\nLine-to-ground 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 ')

27
M.TECH

EPS

13BK1D0XX

SPEC

Zf = input('complex form (for bolted fault enter 0). Zf = ');

fprintf(' \n')
fprintf('Single line to-ground fault at bus No. %g\n', nf)
a =cos(2*pi/3)+j*sin(2*pi/3);
sctm = [1
1
1; 1 a^2 a; 1 a a^2];
Ia0 = V0(nf)/(Zbus1(nf,nf)+Zbus2(nf, nf)+ Zbus0(nf, nf)+3*Zf); Ia1=Ia0;
Ia2=Ia0;
I012=[Ia0; Ia1; Ia2];
Ifabc = sctm*I012;
Ifabcm = abs(Ifabc);
fprintf('Total fault current = %9.4f per unit\n\n', Ifabcm(1))
fprintf('Bus Voltages during the fault in per unit \n\n')
fprintf('
Bus
-------Voltage Magnitude------- \n')
fprintf('
No.
Phase a
Phase b
Phase c \n')
for n = 1:nbus
Vf0(n)= 0 - Zbus0(n, nf)*Ia0;
Vf1(n)= V0(n) - Zbus1(n, nf)*Ia1;
Vf2(n)= 0 - Zbus2(n, nf)*Ia2;
Vabc = sctm*[Vf0(n); Vf1(n); Vf2(n)];
Va(n)=Vabc(1); Vb(n)=Vabc(2); Vc(n)=Vabc(3);
fprintf(' %5g',n)
fprintf(' %11.4f', abs(Va(n))),fprintf(' %11.4f', abs(Vb(n)))
fprintf(' %11.4f\n', abs(Vc(n)))
end
fprintf(' \n')
fprintf('Line currents for fault at bus No. %g\n\n', nf)
fprintf('
From
To
-----Line Current Magnitude---fprintf('
Bus
Bus
Phase a
Phase b
Phase c
for n= 1:nbus
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
Ink1(n, k) = (Vf1(n) - Vf1(k))/ZB1(I);
Ink2(n, k) = (Vf2(n) - Vf2(k))/ZB2(I);
else, end
else, end
end
for I = 1:nbr0
if nl0(I) == n | nr0(I) == n
if nl0(I) ==n
k = nr0(I);
elseif nr0(I) == n k = nl0(I);
end
if k ~= 0
Ink0(n, k) = (Vf0(n) - Vf0(k))/ZB0(I);
else, end
else, end
end
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

\n')
\n')

28
M.TECH

EPS

13BK1D0XX

SPEC

Inkabc = sctm*[Ink0(n, k); Ink1(n, k); Ink2(n, k)];

Inkabcm = abs(Inkabc); th=angle(Inkabc);
if real(Inkabc(1)) > 0
fprintf('%7g', n), fprintf('%10g', k),
fprintf(' %11.4f', abs(Inkabc(1))),fprintf(' %11.4f',
abs(Inkabc(2)))
fprintf(' %11.4f\n', abs(Inkabc(3)))
elseif real(Inkabc(1)) ==0 & imag(Inkabc(1)) < 0
fprintf('%7g', n), fprintf('%10g', k),
fprintf(' %11.4f', abs(Inkabc(1))),fprintf(' %11.4f',
abs(Inkabc(2)))
fprintf(' %11.4f\n', abs(Inkabc(3)))
else, end
else, end
else, end
end
if n==nf
fprintf('%7g',n), fprintf('
F'),
fprintf(' %11.4f', Ifabcm(1)),fprintf(' %11.4f', Ifabcm(2))
fprintf(' %11.4f\n', Ifabcm(3))
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
%Ink0
%Ink1
%Ink2

29
M.TECH

EPS

13BK1D0XX

SPEC

Programme for L-L fault calculations

%
%
%
%
%
%
%
%
%
%
%
%
%
%
%

The program llfault is designed for the line-to-line

fault analysis of a power system network. The program requires
the positive- and negative-sequence bus impedance matrices,
Zbus1, and Zbus2.The bus impedances matrices 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, bus voltages and line currents during the fault.

function llfault(zdata1, Zbus1, zdata2, Zbus2, V)

if exist('zdata2') ~= 1
zdata2=zdata1;
else, end
if exist('Zbus2') ~= 1
Zbus2=Zbus1;
else, end
nl = zdata1(:,1); nr = zdata1(:,2);
R1 = zdata1(:,3); X1 = zdata1(:,4);
R2 = zdata2(:,3); X2 = zdata2(:,4);
ZB1 = R1 + j*X1; ZB2 = R2 + j*X2;
nbr=length(zdata1(:,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('\nLine-to-line 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('Line-to-line fault at bus No. %g\n', nf)
a =cos(2*pi/3)+j*sin(2*pi/3);
sctm = [1
1
1; 1 a^2 a; 1 a a^2];
Ia0=0;
Ia1 = V0(nf)/(Zbus1(nf,nf)+Zbus2(nf, nf)+Zf); Ia2=-Ia1;
I012=[Ia0; Ia1; Ia2];
Ifabc = sctm*I012;
Ifabcm = abs(Ifabc);

30
M.TECH

EPS

13BK1D0XX

SPEC

fprintf('Total fault current = %9.4f per unit\n\n', Ifabcm(2))

fprintf('Bus Voltages during the fault in per unit \n\n')
fprintf('
Bus
-------Voltage Magnitude------- \n')
fprintf('
No.
Phase a
Phase b
Phase c \n')
for n = 1:nbus
Vf0(n)= 0;
Vf1(n)= V0(n) - Zbus1(n, nf)*Ia1;
Vf2(n)= 0 - Zbus2(n, nf)*Ia2;
Vabc = sctm*[Vf0(n); Vf1(n); Vf2(n)];
Va(n)=Vabc(1); Vb(n)=Vabc(2); Vc(n)=Vabc(3);
fprintf(' %5g',n)
fprintf(' %11.4f', abs(Va(n))),fprintf(' %11.4f', abs(Vb(n)))
fprintf(' %11.4f\n', abs(Vc(n)))
end
fprintf(' \n')
fprintf('Line currents for fault at bus No. %g\n\n', nf)
fprintf('
From
To
-----Line Current Magnitude---fprintf('
Bus
Bus
Phase a
Phase b
Phase c

\n')
\n')

for n= 1:nbus
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
Ink0(n, k) = 0;
Ink1(n, k) = (Vf1(n) - Vf1(k))/ZB1(I);
Ink2(n, k) = (Vf2(n) - Vf2(k))/ZB2(I);
Inkabc = sctm*[Ink0(n, k); Ink1(n, k); Ink2(n, k)];
Inkabcm = abs(Inkabc); th=angle(Inkabc);
if real(Inkabc(2)) < 0
fprintf('%7g', n), fprintf('%10g', k),
fprintf(' %11.4f', abs(Inkabc(1))),fprintf(' %11.4f',

abs(Inkabc(2)))

fprintf(' %11.4f\n', abs(Inkabc(3)))

elseif real(Inkabc(2)) ==0 & imag(Inkabc(2)) > 0
fprintf('%7g', n), fprintf('%10g', k),
fprintf(' %11.4f', abs(Inkabc(1))),fprintf(' %11.4f',
abs(Inkabc(2)))

fprintf(' %11.4f\n', abs(Inkabc(3)))

else, end
else, end
else, end

end
if n==nf
fprintf('%7g',n), fprintf('
F'),
fprintf(' %11.4f', Ifabcm(1)),fprintf(' %11.4f', Ifabcm(2))
fprintf(' %11.4f\n', Ifabcm(3))
else, end

end
resp=0;
while strcmp(resp, 'n')~=1 & strcmp(resp, 'N')~=1 & strcmp(resp,
'y')~=1 & strcmp(resp, 'Y')~=1

31
M.TECH

EPS

13BK1D0XX

SPEC

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

Programme for LL-G fault calculations

32
M.TECH

EPS

13BK1D0XX

%
%
%
%
%
%
%
%
%
%
%
%
%
%
%

SPEC

The program dlgfault is designed for the double line-to-ground

fault analysis of a power system network. The program requires
the positive-, negative- or zero-sequence bus impedance matrices,
Zbus1 Zbus2,and Zbus0. The bus impedances matrices 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, bus voltages and line currents during the fault.

function dlgfault(zdata0, Zbus0, zdata1, Zbus1, zdata2, Zbus2, V)

if exist('zdata2') ~= 1
zdata2=zdata1;
else, end
if exist('Zbus2') ~= 1
Zbus2=Zbus1;
else, end
nl = zdata1(:,1); nr = zdata1(:,2);
nl0 = zdata0(:,1); nr0 = zdata0(:,2);
nbr=length(zdata1(:,1)); nbus = max(max(nl), max(nr));
nbr0=length(zdata0(:,1));
R0 = zdata0(:,3); X0 = zdata0(:,4);
R1 = zdata1(:,3); X1 = zdata1(:,4);
R2 = zdata2(:,3); X2 = zdata2(:,4);
for k = 1:nbr0
if R0(k) == inf | X0(k) == inf
R0(k) = 99999999; X0(k) = 999999999;
else, end
end
ZB1 = R1 + j*X1; ZB0 = R0 + j*X0;
ZB2 = R2 + j*X2;
if exist('V') == 1
if length(V) == nbus
V0 = V;
else, end
else, V0 = ones(nbus, 1) + j*zeros(nbus, 1);
end
fprintf('\nDouble line-to-ground fault analysis \n')
ff = 999;
while ff > 0
nf = input('Enter Faulted Bus No. -> ');
while nf <= 0 | nf > nbus

33
M.TECH

EPS

13BK1D0XX

SPEC

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('Double line-to-ground fault at bus No. %g\n', nf)
a =cos(2*pi/3)+j*sin(2*pi/3);
sctm = [1
1
1; 1 a^2 a; 1 a a^2];
Z11 = Zbus2(nf, nf)*(Zbus0(nf, nf)+ 3*Zf)/(Zbus2(nf, nf)+Zbus0(nf, nf)
+3*Zf);
Ia1 = V0(nf)/(Zbus1(nf,nf)+Z11);
Ia2 =-(V0(nf) - Zbus1(nf, nf)*Ia1)/Zbus2(nf,nf);
Ia0 =-(V0(nf) - Zbus1(nf, nf)*Ia1)/(Zbus0(nf,nf)+3*Zf);
I012=[Ia0; Ia1; Ia2];
Ifabc = sctm*I012; Ifabcm=abs(Ifabc);
Ift = Ifabc(2)+Ifabc(3);
Iftm = abs(Ift);
fprintf('Total fault current = %9.4f per unit\n\n', Iftm)
fprintf('Bus Voltages during the fault in per unit \n\n')
fprintf('
Bus
-------Voltage Magnitude------- \n')
fprintf('
No.
Phase a
Phase b
Phase c \n')
for n = 1:nbus
Vf0(n)= 0 - Zbus0(n, nf)*Ia0;
Vf1(n)= V0(n) - Zbus1(n, nf)*Ia1;
Vf2(n)= 0 - Zbus2(n, nf)*Ia2;
Vabc = sctm*[Vf0(n); Vf1(n); Vf2(n)];
Va(n)=Vabc(1); Vb(n)=Vabc(2); Vc(n)=Vabc(3);
fprintf(' %5g',n)
fprintf(' %11.4f', abs(Va(n))),fprintf(' %11.4f', abs(Vb(n)))
fprintf(' %11.4f\n', abs(Vc(n)))
end
fprintf(' \n')
fprintf('Line currents for fault at bus No. %g\n\n', nf)
fprintf('
From
To
-----Line Current Magnitude---fprintf('
Bus
Bus
Phase a
Phase b
Phase c
for n= 1:nbus
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
Ink1(n, k) = (Vf1(n) - Vf1(k))/ZB1(I);
Ink2(n, k) = (Vf2(n) - Vf2(k))/ZB2(I);
else, end
else, end
end
for I = 1:nbr0
if nl0(I) == n | nr0(I) == n
if nl0(I) ==n
k = nr0(I);
elseif nr0(I) == n k = nl0(I);
end
if k ~= 0

\n')
\n')

34
M.TECH

EPS

13BK1D0XX
Ink0(n, k) =
else, end
else, end

SPEC
(Vf0(n) - Vf0(k))/ZB0(I);

end
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
Inkabc = sctm*[Ink0(n, k); Ink1(n, k); Ink2(n, k)];
Inkabcm = abs(Inkabc); th=angle(Inkabc);
if real(Inkabc(2)) < 0
fprintf('%7g', n), fprintf('%10g', k),
fprintf(' %11.4f', abs(Inkabc(1))),fprintf(' %11.4f',
abs(Inkabc(2)))
fprintf(' %11.4f\n', abs(Inkabc(3)))
elseif real(Inkabc(2)) ==0 & imag(Inkabc(2)) > 0
fprintf('%7g', n), fprintf('%10g', k),
fprintf(' %11.4f', abs(Inkabc(1))),fprintf(' %11.4f',
abs(Inkabc(2)))
fprintf(' %11.4f\n', abs(Inkabc(3)))
else, end
else, end
else, end
end
if n==nf
fprintf('%7g',n), fprintf('
F'),
fprintf(' %11.4f', Ifabcm(1)),fprintf(' %11.4f', Ifabcm(2))
fprintf(' %11.4f\n', Ifabcm(3))
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

RESULT:
Finally, became familiar with modelling and analysis of power systems under faulted condition and

35
M.TECH

EPS

13BK1D0XX

SPEC

to compute the fault level, post-fault voltages and currents for different types of faults, both
symmetric and unsymmetric.

Zbus1 =
0 + 0.1450i
0 + 0.1050i
0 + 0.1300i

0 + 0.1050i
0 + 0.1450i
0 + 0.1200i

0 + 0.1300i
0 + 0.1200i
0 + 0.2200i

Zbus0 =
0 + 0.1820i
0 + 0.0545i
0 + 0.1400i

0 + 0.0545i
0 + 0.0864i
0 + 0.0650i

0 + 0.1400i
0 + 0.0650i
0 + 0.3500i

Enter Faulted Bus No. -> 1

Enter Fault Impedance Zf = R + j*X in complex form (for bolted fault enter 0).
Zf = 0.1
Balanced three-phase fault at bus No. 1
Total fault current = 5.6773 per unit
Bus Voltages during fault in per unit
Bus
No.
1
2
3

Voltage
Angle
Magnitude degrees
0.5677 -55.4077
0.6115 -33.6063
0.5741 -46.8769

Line currents for fault at bus No. 1

From
To Current Angle
Bus
Bus Magnitude degrees
G
1
3.2929 -55.4077
1
F 5.6773 -55.4077
G
2
2.3845 -55.4077
2
1
1.8167 -55.4077
2
3
0.5677 -55.4077
3
1
0.5677 -55.4077
Another fault location? Enter 'y' or 'n' within single quote -> 'y'
Enter Faulted Bus No. -> 2
Enter Fault Impedance Zf = R + j*X in complex form (for bolted fault enter 0). Zf =
0.1

36
M.TECH

EPS

13BK1D0XX

SPEC
Balanced three-phase fault at bus No. 2
Total fault current = 5.6773 per unit
Bus Voltages during fault in per unit
Bus
No.
1
2
3

Voltage
Angle
Magnitude degrees
0.6115 -33.6063
0.5677 -55.4077
0.5852 -41.3715

Line currents for fault at bus No. 2

From
To
Bus
Bus
G
1
1
2
1
3
G
2
2
F
3
2

Current Angle
Magnitude degrees
2.3845 -55.4077
1.8167 -55.4077
0.5677 -55.4077
3.2929 -55.4077
5.6773 -55.4077
0.5677 -55.4077

Another fault location? Enter 'y' or 'n' within single quote -> 'y'
Enter Faulted Bus No. -> 3
Enter Fault Impedance Zf = R + j*X in complex form (for bolted fault enter 0). Zf =
0.1
Balanced three-phase fault at bus No. 3
Total fault current = 4.1380 per unit
Bus Voltages during fault in per unit
Bus
No.
1
2
3

Voltage
Angle
Magnitude degrees
0.5567 -23.5688
0.5852 -20.5560
0.4138 -65.5560

Line currents for fault at bus No. 3

37
M.TECH

EPS

13BK1D0XX

SPEC
From
To
Bus
Bus
G
1
1
3
G
2
2
1
2
3
3
F

Current Angle
Magnitude degrees
2.1518 -65.5560
2.4828 -65.5560
1.9863 -65.5560
0.3310 -65.5560
1.6552 -65.5560
4.1380 -65.5560

EXPERIMENT-6
AIM:
To become familiar with various aspects of the transient and small signal stability
38
M.TECH

EPS

13BK1D0XX

SPEC

analysis of Single-Machine infinite Bus (SMIB) system.

PROGRAM-REQUIRED:
MATLAB:7.8
THEORY:
Stability: Stability problem is concerned with the behaviour of power system when it is
subjected to disturbance and is classified into small signal stability problem if the
disturbances are small and transient stability problem when the disturbances are large.
Transient stability: When a power system is under steady state, the load plus
transmission loss equals to the generation in the system. The generating units run a
synchronous speed and system frequency, voltage, current and power flows are steady.
When a large disturbance such as three phase fault, loss of load, loss of generation etc.,
occurs the power balance is upset and the generating units rotors experience either
acceleration or deceleration. The system may come back to a steady state condition
maintaining synchronism or it may break into subsystems or one or more machines may
pull out of synchronism. In the former case the system is said to be stable and in the later
case
it
is
said
to
be
unstable.
Small signal stability: When a power system is under steady state, normal operating
condition, the system may be subjected to small disturbances such as variation in load
and generation, change in field voltage, change in mechanical toque etc., The nature
of system response to small disturbance depends on the operating conditions, the
transmission system strength, types of controllers etc. Instability that may result
from small disturbance may
be
of
two
forms,
(i) Steady increase in rotor angle due to lack of synchronising torque.
(ii) Rotor oscillations of increasing magnitude due to lack of sufficient damping torque.

FORMULA:

39
M.TECH

EPS

13BK1D0XX

SPEC

PROCEDURE:
1. Enter the command window of the MATLAB.
40
M.TECH

EPS

13BK1D0XX

SPEC

2. Create a new M file by selecting File - New M File

3. Type and save the program.
4. Execute the program by either pressing tools Run
5. View the results.
EXERCISE:
1. A 60Hz synchronous generator having inertia H = 5 constant J/MVA and a direct axis
transient reactance Xd
1 = 0.3 per unit is connected to an infinite bus through a purely
reactive circuit as shown in figure. Reactances are marked on the diagram on a common
system base. The generator is delivering real power Pe = 0.8 per unit and Q = 0.074
per unit to the infinite bus at a voltage of V = 1 per unit.

a) A temporary three-phase fault occurs at the sending end of the line at point F.When
the fault is cleared, both lines are intact. Determine the critical clearing angle and
the critical fault clearing time.
b) Verify the result using MATLAB program.
PROGRAM:
% This program obtains the power angle curves for a one-machine system
% before fault, during fault and after the fault clearance.
% The equal area criterion is applied to find the critical clearing
angle
% for the machine to stay synchronized to the infinite bus bar
Pm = 0.8; E = 1.17; V = 1.0;
X1 = 0.65; X2 = inf; X3 = 0.65;
function eacfault(Pm, E, V, X1, X2, X3)
if exist('Pm')~=1
Pm = input('Generator output power in p.u. Pm = '); else, end
if exist('E')~=1
E = input('Generator e.m.f. in p.u. E = '); else, end
if exist('V')~=1
V = input('Infinite bus-bar voltage in p.u. V = '); else, end
if exist('X1')~=1
X1 = input('Reactance before Fault in p.u. X1 = '); else, end
if exist('X2')~=1
X2 = input('Reactance during Fault in p.u. X2 = '); else, end
if exist('X3')~=1
X3 = input('Reactance aftere Fault in p.u. X3 = '); else, end
Pe1max = E*V/X1; Pe2max=E*V/X2; Pe3max=E*V/X3;
delta = 0:.01:pi;

41
M.TECH

EPS

13BK1D0XX

SPEC

Pe1 = Pe1max*sin(delta); Pe2 = Pe2max*sin(delta); Pe3 =

Pe3max*sin(delta);
d0 =asin(Pm/Pe1max); dmax = pi-asin(Pm/Pe3max);
cosdc = (Pm*(dmax-d0)+Pe3max*cos(dmax)-Pe2max*cos(d0))/(Pe3max-Pe2max);
if abs(cosdc) > 1
fprintf('No critical clearing angle could be found.\n')
fprintf('system can remain stable during this disturbance.\n\n')
return
else, end
dc=acos(cosdc);
if dc > dmax
fprintf('No critical clearing angle could be found.\n')
fprintf('System can remain stable during this disturbance.\n\n')
return
else, end
Pmx=[0 pi-d0]*180/pi; Pmy=[Pm Pm];
x0=[d0 d0]*180/pi; y0=[0 Pm]; xc=[dc dc]*180/pi; yc=[0 Pe3max*sin(dc)];
xm=[dmax dmax]*180/pi; ym=[0 Pe3max*sin(dmax)];
d0=d0*180/pi; dmax=dmax*180/pi; dc=dc*180/pi;
x=(d0:.1:dc);
y=Pe2max*sin(x*pi/180);
y1=Pe2max*sin(d0*pi/180);
y2=Pe2max*sin(dc*pi/180);
x=[d0 x dc];
y=[Pm y Pm];
xx=dc:.1:dmax;
h=Pe3max*sin(xx*pi/180);
xx=[dc xx dmax];
hh=[Pm h Pm];
delta=delta*180/pi;
if X2 == inf
fprintf('\nFor this case tc can be found from analytical formula. \n')
H=input('To find tc enter Inertia Constant H, (or 0 to skip) H = ');
if H ~= 0
d0r=d0*pi/180; dcr=dc*pi/180;
tc = sqrt(2*H*(dcr-d0r)/(pi*60*Pm));
else, end
else, end
%clc
fprintf('\nInitial power angle
= %7.3f \n', d0)
fprintf('Maximum angle swing
= %7.3f \n', dmax)
fprintf('Critical clearing angle = %7.3f \n\n', dc)
if X2==inf & H~=0
fprintf('Critical clearing time = %7.3f sec. \n\n', tc)
else, end
h = figure; figure(h);
fill(x,y,'m')
hold;
fill(xx,hh,'c')
plot(delta, Pe1,'-', delta, Pe2,'r-', delta, Pe3,'g-', Pmx, Pmy,'b-',
x0,y0, xc,yc, xm,ym), grid
Title('Application of equal area criterion to a critically cleared
system')
xlabel('Power angle, degree'), ylabel(' Power, per unit')
text(5, 1.07*Pm, 'Pm')
text(50, 1.05*Pe1max,['Critical clearing angle = ',num2str(dc)])
axis([0 180 0 1.1*Pe1max])

42
M.TECH

EPS

13BK1D0XX

SPEC

hold off;

CASE-II
Pm = 0.8; E = 1.17; V = 1.0;
X1 = 0.65; X2 = 1.8; X3 = 0.8;
eacfault(Pm, E, V, X1, X2, X3)
RESULT:
Finally,became familiar with various aspects of the transient and small signal stability
analysis of Single-Machine-Infinite Bus (SMIB) system.
Application of equal area criterion to a critically cleared system
Critical clearing angle = 84.7745

1.8
1.6

Power, per unit

1.4
1.2
1
Pm

0.8
0.6
0.4
0.2
0

20

40

60

80
100
120
Power angle, degree

140

160

180

EXPERIMENT-7

43
M.TECH

EPS

13BK1D0XX

SPEC

AIM:

To perform PSPICE simulation on the following converter circuits.

a) Single phase full converter with R load
b) Single phase full converter with R-L load
c) Single phase full converter with R-L-E load
APPARATUS:

PSPICE software
CIRCUIT DIAGRAM FOR R LOAD:

SUB CIRCUIT:

44
M.TECH

EPS

13BK1D0XX

SPEC

PROGRAM FOR R LOAD:

SINGLE PHASE FULL CONVERTER WITH R LOAD

VS 1 2 SIN(0 169.7 50)
VG1 3 7 PULSE(0 10 1666.6US 1NS 1NS 100US 20MS)
VG2 6 2 PULSE(0 10 1666.6US 1NS 1NS 100US 20MS)
VG3 4 7 PULSE(0 10 11666.6US 1NS 1NS 100US 20MS)
VG4 5 1 PULSE(0 10 11666.6US 1NS 1NS 100US 20MS)
XT1 1 7 3 7 SCR
XT2 0 2 6 2 SCR
XT3 2 7 4 7 SCR
XT4 0 1 5 1 SCR
R 7 0 10OHM
.SUBCKT SCR 1 2 3 2
S1 1 5 6 2 SMOD
RG 3 4 50OHM
VX 4 2 DC 0V
VY 5 7 DC 0V
DT 7 2 DMOD
RT 6 2 1OHM
CT 6 2 10UF
F1 2 6 POLY(2) VX VY 0 50 11
.MODEL SMOD VSWITCH (RON=0.0125 ROFF=10E+5 VON=0.5V VOFF=0V)
.MODEL DMOD D(IS=2.2F-15 BV=1800V TT=0)
.ENDS SCR
.TRAN 1US 60MS
.PROBE
.END

CIRCUIT DIAGRAM FOR R-L LOAD:

45
M.TECH

EPS

13BK1D0XX

SPEC

SUB CIRCUIT:

PROGRAM FOR R-L LOAD:

46
M.TECH

EPS

13BK1D0XX

SPEC

SINGLE PHASE FULL CONVERTER WITH R-L LOAD

VS 1 2 SIN(0 169.7 50)
VG1 3 7 PULSE(0 10 1666.6US 1NS 1NS 100US 20MS)
VG2 6 2 PULSE(0 10 1666.6US 1NS 1NS 100US 20MS)
VG3 4 7 PULSE(0 10 11666.6US 1NS 1NS 100US 20MS)
VG4 5 1 PULSE(0 10 11666.6US 1NS 1NS 100US 20MS)
XT1 1 7 3 7 SCR
XT2 0 2 6 2 SCR
XT3 2 7 4 7 SCR
XT4 0 1 5 1 SCR
R 7 8 10OHM
L 8 0 10MH
.SUBCKT SCR 1 2 3 2
S1 1 5 6 2 SMOD
RG 3 4 50OHM
VX 4 2 DC 0V
VY 5 7 DC 0V
DT 7 2 DMOD
RT 6 2 1OHM
CT 6 2 10UF
F1 2 6 POLY(2) VX VY 0 50 11
.MODEL SMOD VSWITCH (RON=0.0125 ROFF=10E+5 VON=0.5V VOFF=0V)
.MODEL DMOD D(IS=2.2F-15 BV=1800V TT=0)
.ENDS SCR
.TRAN 1US 60MS
.PROBE
.END

CIRCUIT DIAGRAM FOR R-L-E LOAD:

47
M.TECH

EPS

13BK1D0XX

SPEC

SUB CIRCUIT:

48
M.TECH

EPS

13BK1D0XX

SPEC

PROGREM FOR R-L-E LOAD:

SINGLE PHASEFULL CONVERTER WITH RLE-LOAD
VS 1 2 SIN(0 169.7 50)
VG1 3 7 PULSE(0 10 1666.6US 1NS 1NS 100US 20MS)
VG2 6 2 PULSE(0 10 1666.6US 1NS 1NS 100US 20MS)
VG3 4 7 PULSE(0 10 11666.6US 1NS 1NS 100US 20MS)
VG4 5 1 PULSE(0 10 11666.6US 1NS 1NS 100US 20MS)
XT1 1 7 3 7 SCR
XT2 0 2 6 2 SCR
XT3 2 7 4 7 SCR
XT4 0 1 5 1 SCR
R 7 8 10OHM
L 8 9 10MH
VE 9 0 DC 5V
.SUBCKT SCR 1 2 3 2
S1 1 5 6 2 SMOD
RG 3 4 50OHM
VX 4 2 DC 0V
VY 5 7 DC 0V
DT 7 2 DMOD
RT 6 2 1OHM
CT 6 2 10UF
F1 2 6 POLY(2) VX VY 0 50 11
.MODEL SMOD VSWITCH (RON=0.0125 ROFF=10E+5 VON=0.5V VOFF=0V)
.MODEL DMOD D(IS=2.2F-15 BV=1800V TT=0)
.ENDS SCR
.TRAN 1US 60MS
.PROBE
.END

RESULT FOR R LOAD:

49
M.TECH

EPS

13BK1D0XX

SPEC

RESULT FOR R-L LOAD:

50
M.TECH

EPS

13BK1D0XX

SPEC

RESULT FOR R-L-E LOAD:

51
M.TECH

EPS

13BK1D0XX

SPEC

EXPERIMENT-8
AIM: To Develop MATLAB model for Sinusoidal Pulse Width Modulation.

52
M.TECH

EPS

13BK1D0XX

SPEC

SOFTWARE REQUIRED:
MATLAB-(SIMULINK)-7.8
Theory:
Bipolar PWM Modulation Approach
Bipolar PWM switching is a classical switching scheme for single-phase inverter. The
switch pairs (T1, T4) and (T2, T3) in Figure 6.2 on the different legs are switched on and
off simultaneously. This results a bipolar voltage output Vout because there is no zero
output voltage state exists. The output waveform is the same as the point voltage V in
Figure 1.2 but the amplitude doubles. The principle of bipolar PWM can be summarized
in equation 2.1:

vout vd

When

vcontrol vcarrier

vout vd

When

vcontrol vcarrier

Figure

2.0 gives a demonstration of

bipolar modulation approach with parameters mi=0.8 and mf=15. From the figure we
know that VBO is exactly the reverse of VAO at any time, so there is no zero state in the
53
M.TECH

EPS

13BK1D0XX

SPEC

output voltage Vout=VAO-VBO, which makes output signal bipolar. In this bipolar PWM
modulation scheme, if the frequency modulation ratio is selected to be odd, the output
waveform Vout will be odd and half-wave symmetry with origin. This results a
disappearance of all the even harmonics in Vout.
OPERATION OF UNIPOLAR SINUSOIDAL PWM TECHNIQUE
In a unipolar PWM modulation, the switch pairs (T1, T4) and (T2, T3) in inverter
H-bridge are not switched simultaneously as in the bipolar approach. Instead, the leg A
and B in Figure 6.3 are controlled independently by the comparison of two control
signals Vcontrol1 and Vcontrol2 with the same carrier signal. Figure 2.2 illustrates the unipolar
PWM switching scheme and output voltage at mi=0.8 and mf=12
.

Figure 6.3 Unipolar PWM scheme and output

The switching algorithm of unipolar PWM is given in equation below:

54
M.TECH

EPS

13BK1D0XX

SPEC

Vout=Vd

when T1,T4 is on

Vout=-Vd

when T2,T4 is on

Vout=0

when T1,T3 or T2,T4

From the above Figure 2.1 it is quite clear that the output voltage changes between 0 and
+Vd or between 0 and - Vd in each half fundamental period, so this PWM scheme is
called unipolar PWM. The voltage change at each switching transition is Vd in
comparison with the 2Vd voltage change in bipolar PWM approach. Also in this PWM
approach, it is possible to choose the frequency modulation ratio mf as even number to
cancel the harmonic components at the switching frequency (fs ).
OPERATION OF MODIFIED UNIPOLAR SINUSOIDAL PWM TECHNIQUE
A modified unipolar PWM approach, which is quite different from the traditional bipolar
or unipolar PWM schemes, is presented in this section. The two legs in the H- Bridge of
the single-phase inverter are switched at different frequencies: leg A is switched at the
fundamental frequency f1 while leg B is switched at the carrier frequency fs , which is
much higher than f1. This PWM switching scheme establishes an unbalanced switching
speed between two legs of the inverter. It also generates a unipolar output voltage that
changes between 0 and +Vdor between 0 and -Vd. Figure 6.4 is a demonstration of this
modified unipolar PWM approach with parameters m i=0.8 and

mf =15. To obtain odd

and half-wave symmetry of the output voltage to suppress the even number harmonics,
the frequency modulation ratio should be chosen as an odd integer in this PWM
approach.

55
M.TECH

EPS

13BK1D0XX

SPEC

Modified unipolar PWM scheme and output

The modified unipolar approach has the advantage of reducing high frequency electrical
magnetic interference (EMI).As shown in Figure 6.4, unlike the previous two PWM
approaches, the modulation signal Vcontrol in this modified unipolar PWM approach is a
discontinuous signal.

56
M.TECH

EPS

13BK1D0XX

SPEC

SIMULINK DIAGRAM:

57
M.TECH

EPS

13BK1D0XX

SPEC

RESULT

58
M.TECH

EPS

13BK1D0XX

SPEC

EXPERIMENT-9
Aim: Developing MATLAB model for Closed Loop Speed Control of Separately Excited D.C
Motor.

Theory:
The principle of speed control for dc motors is developed from the basic emf equation of
the motors. Torque, flux current, induced emf, and speed are normalized to present the
motor characteristics. Two types of control arc available: armature control and field
control. These methods arc combined to yield a wide range of speed control. The torquespeed characteristics of the motor are discussed for both directions of rotation and
delivering both motoring and regenerating torques in any direction of rotation. Such an
operation, known as four quadrant operation, has a unique set of requirements on the
input voltage and current to the dc motor armature and field. These requirements are
identified for specifying the power stage [8]. Modern power converters constitute the
power stage for variable-speed dc drives. These power converters are chosen for a
particular application depending on a number of factors such as cost, input power source,
harmonics, power factor, noise, and speed of response Controlled-bridge rectifiers fed
from single-phase and three-phase ac supply are considered in this chapter. The theory,
operation and control of the three-phase controlled-bridge rectifier is considered in detail
because of its wide-spared use. A model for the power Converter is derived for use in
simulation and controller design. Two and four quadrant dc motor drives and their
control are developed.

Speed control techniques of separately excited dc motor:

The speed of a separately excited dc motor could be varied from zero to rated speed
mainly by varying armature voltage in the constant torque region. Whereas in the
constant power region, field flux should be reduced to achieve speed above the rated
speed. The motor drives a mechanical load characterized by inertia J, friction coefficient
Dm, and load torque TL. The specifications of the dc motor are detailed as follows:
Specifications
Shaft power 5 hp
Rated voltage 240 V
Armature resistance 0.6 _
Armature inductance 0.012 H
Field resistance 240 _
Field inductance 120 H
Total inertia (J) 1 kgm2
Viscous friction coefficient (B) 0.02 Nms
Coulomb friction torque (Tf) 0 Nm
a) Modeling and control of SEDM using MATLAB/SimPowerSystems: Fig. 1
showsthespeedcontrolcircuitofanarmaturecontrolledseparatelyexciteddcmotor
using chopper circuit, and in Fig. its MATLLAB/SimPowerSystems model [31] is
59
M.TECH

EPS

13BK1D0XX

SPEC

shown.ItconsistsofaseparatelyexciteddcmotorfedbyaDCsourcethrougha
choppercircuit.AsingleGTOthyristorwithitscontrolcircuitandafreewheeling
diodeformthechoppercircuit.Themotordrivesamechanicalloadcharacterizedby
inertiaJ,frictioncoefficientB,andloadtorqueTL.Thecontrolcircuitconsistsofa
speedcontrolloopandacurrentcontrolloop.Aproportionalintegral(PI)controlled
speedcontrolloopsensestheactualspeedofthemotorandcompares itwiththe
referencespeedtodeterminethereferencearmaturecurrentrequiredbythemotor.
Onemaynotethatanyvariationintheactualspeedisameasureofthearmature
currentrequiredbythemotor.Thecurrentcontrolloopconsistsofahysteresiscurrent
controller(HCC).Theblockdiagramofahysteresiscurrentcontrollerisshownin
Fig.HCCisusedtogenerateswitchingpatternsrequiredforthechoppercircuitby
comparingtheactualcurrentbeingdrawnbythemotorwiththereferencecurrent.A
positivepulseisgeneratediftheactualcurrentislessthanreferencearmaturecurrent,
whereasanegativepulseisproducediftheactualcurrentexceedsreferencecurrent.
Hysteresis current control is a method of controlling a power electronic converter
so that an output current is generated which follows a reference current waveform. A
hysteresis current controller is implemented with a closed loop control. The difference
between the desired current, and the current being injected is used to control the
switching of the chopper circuit. When the error reaches an upper limit namely upper
hysteresis limit, GTO is switched to force the current the current down. On the other hand
when the error reaches the lower hysteresis limit, a positive pulse is produced to increase
the current. The minimum and maximum values of the error signal are e min and emax. The
range of the error signal, emax-emin, directly controls the amount of ripple in the output
current and is called the hysteresis band. Thus the armature current is forced to stay
within the hysteresis band determined by the upper and lower hysteresis limits.
b) Modeling and control of SEDM using Simulink model: The speed control circuit of
a SEDM using Simulink is shown in Fig. , where
Vt
Eb
Ra
La
Rf
Lf
If
Ia
wm

Dm

Supply voltage (V)

Back emf (V)
Armature resistance ( )
Armature inductance (H)
Field resistance ( )
Field inductance (H)
Field current (A)
Armature current (A)
Speed (rad/s)
Rotor inertia of motot
(kgm2)
Viscous friction of motor
(Nms)

In Fig.4, the GTO is modeled using a switch. The switch block has three inputs: the
middle input controls which of the two other inputs is routed to the output. If the control
input is one, 240 V is routed to the output, on the other hand if the control input is zero, a
zero will be routed to the output.
60
M.TECH

EPS

13BK1D0XX

SPEC

61
M.TECH

EPS

13BK1D0XX

SPEC

Results:

SpeedcontrollersystembasedonNARMAL2controllerhasbeensuccessfully
developedusingMATLABtocontrolthespeedofaseparatelyexciteddcmotor.

62
M.TECH

EPS

13BK1D0XX

SPEC

EXPERIMENT-10
AIM: Simulation of RLC Circuit using PSPICE

APPARATUS:
PSPICE software
Theory:
Linear Inductors in PSpice

The next passive element we add to our parts list is the linear inductor. This part name
begins with the letter, L, in column 1 of the source listing. Be aware that PSpice enables
this part to access a nonlinear model description. That will be explained in a later
tutorial. For now, our inductor model is a linear device incapable of saturation.The
inductor stores energy in its magnetic field. This makes it necessary to be able to specify
its initial current in a simulation. Although we can include inductors in DC circuit
simulations, there is usually little advantage in doing so because the inductor behaves as
a short circuit under steady-state DC excitation. In steady-state AC simulations the
inductor behaves as an imaginary impedance. We do not specify initial current in an
inductor in either of those steady-state conditions. However, when imulating transient
operations, we often need to specify this initial current

The figure shown above shows the circuit symbol for an inductor with node designations
of "1" and "2," n initial current of 2.5 A, and a value of 50 mH. An appropriate code
listing for entering this element nto a PSpice circuit file is:

63
M.TECH

EPS

13BK1D0XX

SPEC

Note hat the initial current is assumed to flow from the first node in the node list through
the inductor owards the second node in the node list. If there is a need to change the
direction of this initial current, either reverse the order of the nodes in the node list or
place a minus sign in front of the value of the nitial current. For better readability, the
above line could be written as:

"H" for henrys and the "A" for amps will be ignored by PSpice.

Linear Capacitors in PSpice

The capacitor is the second energy storing circuit component we add to our parts list. We
will assume hat the capacitor is ideal in the sense of being linear and lossless. Since it
can store energy, PSpice povides a method for specifying the initial voltage across the
capacitor. This is useful for simulations of tansient behavior of circuits with capacitors.
The figure shown below illustrates a capacitor with node designations and an initial
voltage of 20 from node 4 to node 5. The part name for a capacitor must start with the
letter C.

An appropriate code listing to represent this capacitor in a PSpice listing is:

The capacitance of the above element is 50F. This can be represented as "50u" in
PSpice. (See PSpice utorial No. 1 for a description of the system for metric prefixes in
PSpice.) Note that the polarity of the iitial voltage (as shown) is such that the positive
side is the first node in the list with the negative side on the second node in the list. To
reverse the polarity of the initial voltage for the simulation, either reverse the order of the
nodes in the node list or place a minus sign in front of the value in the "IC=" phrase. For
better clarity, the above capacitor could be coded as:

PSpice would ignore the "F" for farads and the "V" for volts.

64
M.TECH

EPS

13BK1D0XX

SPEC

Circuit diagram:

Programe:
60 Hz AC Circuit
Vs 1 0 AC 120V 0
R 1 2 0.5
L 2 3 3.183mH
C 3 0 132.8uF
.AC LIN 1 60 60
.PRINT AC VM(3) VP(3) IM(Rm) IP(Rm) IM(Cx) IP(Cx)
.END

65
M.TECH

EPS