PSP 253

Download as pdf or txt
Download as pdf or txt
You are on page 1of 109

Power System & Simulation Lab

1. LOAD FREQUENCY CONTROL OF SINGLE AREA POWER


SYSTEM
AIM:
To become familiar with modelling and analysis of the frequency and tie-line flow
dynamics of a power system without and with load frequency controllers (LFC) and to
design better controller for getting better responses using MATLAB/Simulink.

SOFTWARE REQUIRED:
MATLAB software package.

ALGORITHM:
STEP: 1 Click the Simulink icon.
STEP: 2 Open the untitled window and create the model for loan frequency control of
single area power system using various library functions available in the
Simulink browser.
STEP: 3 Save the file and simulate the model.
STEP: 4 Double click the scope and we can get the waveform as change in frequency
with respect to time.
STEP: 5 simulate the model for uncontrolled case and observe the waveform.
STEP: 6 Simulate the model for controlled case (i.e.), introduce one integral
control in the model and observe the waveform.
STEP: 7 Finally, compare the results obtained from simulation and comment the
result.

THEORY:
If the system is connected to numerous loads in a power system, then the system
frequency and speed change with the characteristics of the governor as the load changes. If
it’s not required to maintain the frequency constant in a system then the operator is not
required to change the setting of the generator. But if constant frequency is required the
operator can adjust the velocity of the turbine by changing the characteristics of the governor
when required. If a change in load is taken care by two generating stations running parallel
then the complex nature of the system increases. The ways of sharing the load by two
machines are as follow:

ASCET/EEE/B.TECH/20G25A0253 1
Power System & Simulation Lab

Block diagram:

Model graph:

ASCET/EEE/B.TECH/20G25A0253 2
Power System & Simulation Lab

1) Suppose there are two generating stations that are connected to each other by tie line. If the
change in load is either at A or at B and the generation of A is regulated so as to have
constant frequency then this kind of regulation is called as Flat Frequency Regulation.
2) The other way of sharing the load is that both A and B would regulate their generations to
maintain the frequency constant. This is called parallel frequency regulation.
3) The third possibility is that the change in the frequency of a particular area is taken care of
by the generator of that area thereby maintain the tie-line loading. This method is known as
flat tie line loading control.
4) In Selective Frequency control each system in a group is taken care of the load changes on
its own system and does not help the other systems, the group for changes outside its own
limits.
5) In Tie-line Load-bias control all the power systems in the interconnection aid in regulating
frequency regardless of where the frequency change originates.

SPEED GOVERNING SYSTEM:


MATHEMATICAL MODELLING OF A GENERATOR:
With the use of swing equation of a synchronous machine to small perturbation, we
have

Or in terms of small change in speed

Laplace Transformation gives,

Fig: Mathematical modelling Block Diagram of Load

ASCET/EEE/B.TECH/20G25A0253 3
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 4
Power System & Simulation Lab

MATHEMATICAL MODELLING FOR PRIME MOVER:


The source of power generation is the prime mover. It can be hydraulic turbines near
waterfalls, steam turbine whose energy come from burning of coal, gas and other fuels. The
model of turbine relates the changes in mechanical power output ∆Pm and the changes in the
steam valve position APV.

where the turbine constant is in the range of 0.2 -2.0s


MATHEMATICAL MODELLING FOR GOVERNOR:
When the electrical load is increased suddenly then the electrical power exceeds the
input mechanical power. This deficiency of power in the load side is compensated from the
kinetic energy of the turbine. Due to this reason the energy that is stored in the machine is
decreased and the governor sends signal for supplying more volumes of water, steam or gas
to increase the speed of the prime mover to compensate deficiency in speed.

In s-domain

The command AP, is transformed through amplifier to the steam valve position command
∆Pᵥ. We assume here a linear relationship and considering simple time constant we get this-
domain relation.

Combining all the above block diagrams, for a isolated area system we get the following:

Fig: complete block diagram of single area system


The closed loop transfer function that relates the load change to the frequency deviation

ASCET/EEE/B.TECH/20G25A0253 5
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 6
Power System & Simulation Lab

PROBLEM:
Single area system
An Isolated power station has the following parameters
Turbine time constant = 0.5 sec
Governor time constant = 0.2 sec
Generator inertia constant = 5 sec
Governor speed regulation = R per unit
The load varies by 0.8 percent for a 1 percent change in frequency i.e., D=0.8. The
governor speed regulation is set to R=0.05 pu. The turbine rated output is 250 MW at
nominal frequency of 60Hz. A sudden load change of 50 MW (APD=0.2 pu) occurs.
i). Find the steady state frequency deviation in Hz.
ii). Draw a MATLAB/Simulink model and also verify the output with the manual

ASCET/EEE/B.TECH/20G25A0253 7
Power System & Simulation Lab

SIMULINK Block Diagram for single area system:

ASCET/EEE/B.TECH/20G25A0253 8
Power System & Simulation Lab

Simulation Output:

RESULT:
Thus, the load frequency dynamics of single area power system is analysed using Simulink
and the outputs are verified.

ASCET/EEE/B.TECH/20G25A0253 9
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 10
Power System & Simulation Lab

2. LOAD FREQUENCY CONTROL OF TWO AREA POWER


SYSTEM
AIM:
To become familiar with modelling and analysis of the frequency and tie-line flow
dynamics of a power system without and with load frequency controllers (LFC) and to design
better controllers for getting better responses using MATLAB Simulink.

SOFTWARE REQUIRED:
MATLAB software package.

ALGORITHM:
STEP: 1 Click the Simulink icon.
STEP: 2 Open the untitled window and create the model for loan frequency control of
two area power system using various library functions available in the
Simulink browser.
STEP: 3 Save the file and simulate the model.
STEP: 4 Double click the scope and we can get the waveform as change in frequency
with respect to time.
STEP: 5 simulate the model for uncontrolled case and observe the waveform.
STEP: 6 Simulate the model for controlled case (i.e.), introduce one integral control in
the model and observe the waveform.
STEP: 7 Finally compare the results obtained from simulation and comment the result.

THEORY:
Modern day power systems are divided into various areas. For example, in India there
are five regional grids, e.g., Eastern Region, Western Region etc. Each of these areas is
generally interconnected to its neighbouring areas. The transmission lines that Connect an
area to its neighbouring area are called tie-lines. Power sharing between two areas occurs
through these tie-lines. Load frequency control, as the name signifies, regulates the power
flow between different areas while holding the frequency constant.
As we have that the system frequency rises when the load decreases if ∆Pref is kept at
zero. Similarly, the frequency may drop if the load increases. However, it is desirable to
maintain the frequency constant such that ∆f=0. The power flow through different tie-lines
are scheduled for example, area- i may export a pre-specified. Amount of power to area-j
while importing another pre-specified amount of power from area- k. However, it is expected
that to fulfil this obligation, area-i absorbs its own load change, i.e., increase generation to

ASCET/EEE/B.TECH/20G25A0253 11
Power System & Simulation Lab

Block Diagram:

Fig: Two are hydrothermal system

Model graph:

ASCET/EEE/B.TECH/20G25A0253 12
Power System & Simulation Lab

supply extra load in the area or decrease generation when the load demand in the area has
reduced. While doing this area- i must however maintain its obligation to areas j and k as far
as importing and exporting power is concerned. A conceptual diagram of the interconnected
areas is shown in Fig. 2.1.

Reasons for the Limits on Frequency:


Following are the reasons for keeping a strict limit on the system frequency variation:
1. The speed of the alternating current motors depends on the frequency of the power supply.
There are situations where speed consistency Is expected to be of high order.
2. The electric clocks are driven by the synchronous motors. The accuracy of the clocks are
not only dependent on the frequency but also is an integral of the frequency error,
3. If the normal frequency Id 50 Hertz and the system frequency falls below 47.5 Hertz or
goes up above 52.5 Hertz then the blades of the turbine are likely to get damaged so as to
prevent the stalling of the generator
4. The under-frequency operation of the power transformer is not desirable. For constant
system voltage if the frequency is below the desired level than the normal flux in the core
increases. This sustained under frequency operation of the power transformer results in low
efficiency and over-heating of the transformer windings.
5. The most serious effect of subnormal frequency operation is observed in the case of
Thermal Power Plants. Due to the subnormal frequency operation the blast of the ID and FD
fans in the power stations get reduced and thereby reduce the generation power in the thermal
plants. This phenomenon has got a cumulative effect and in turn is able to make complete
shutdown of the power plant if proper steps of load shedding technique is not engaged. It is
pertinent to mention that, in load shedding technique a sizable chunk of load from the power
system is disconnected from the generating units so as to restore the frequency to the desired
level.
We can therefore state that the load frequency control (LFC) has the following two objectives:
 Hold the frequency constant (∆f=0) against any load change. Each area must
contribute to absorb any load change such that frequency does not deviate.
 Each area must maintain the tie-line power flow to its pre-specified value.
The first step in the LFC is to form the area control error (ACE) that is defined as

ASCET/EEE/B.TECH/20G25A0253 13
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 14
Power System & Simulation Lab

Where Ptie and Psch are tie-line power and scheduled power through tie-line respectively
and the constant Bris called the frequency bias constant.
The change in the reference of the power setting ∆Pref,I, of the area- j is then obtained
by the feedback of the ACE through an integral controller of the form

Where Ki is the integral gain.


The ACE is negative if the net power flow out of an area is low or if the frequency
has dropped or both. In this case the generation must be increased. This can be achieved by
increasing ∆Pref,I, This negative sign accounts for this inverse relation between ∆Pref,I, and
ACE. The tie-line power flow and frequency of each area are monitored in its control centre.
Once the ACE is computed and ∆Pref,I, is obtained, commands are given to various turbine-
generator controls to adjust their reference power settings. The block diagram for two area
load frequency control is shown in below.

Fig: Two area load frequency control

ASCET/EEE/B.TECH/20G25A0253 15
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 16
Power System & Simulation Lab

PROBLEM:
Two area system
A two-area system connected by a tie line has the following parameters on a
1000MVA common base
Area 1 2
Speed regulation R1=0.05 R2=0.0625
Frequency sensitive load coefficient D1=0.6 D2=0.9
Inertia Constant H1=5 H2=4
Base power 1000 MVA 1000 MVA
Governor time constant
Turbine time constant

The units are operating in parallel at the nominal frequency of 50Hz. The synchronizing
power coefficient is computed from the initial operating condition and is given to be Ps=2 pu.
A load change of 187.5MW occurs in areal. Determine the new steady state frequency and
the change in tie-line flow.

ASCET/EEE/B.TECH/20G25A0253 17
Power System & Simulation Lab

SIMULINK Block Diagram for two area system:

ASCET/EEE/B.TECH/20G25A0253 18
Power System & Simulation Lab

Simulation Output:

RESULT:
Thus, the load frequency dynamics of two area Power system is analysed using
MATLAB Simulink and the outputs are verified.

ASCET/EEE/B.TECH/20G25A0253 19
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 20
Power System & Simulation Lab

3. COMPUTATION OF PARAMETERS AND MODELLING OF


TRANSMISSION LINES
AIM:
To develop a program in MATLAB for determining the transmission line parameters
and verify using MATLAB\ simulation.

SOFTWARE REQUIRED:
MATLAB software package.

THEORY:
A transmission line has *three constants R, L and C distributed uniformly along the
whole length of the line. The resistance and inductance form the series impedance. The
capacitance existing between conductors for 1-phase line or from a conductor to neutral for a
3-phase line forms a shunt path throughout the length of the line. Therefore, capacitance
effects introduce complications in transmission line calculations. Depending upon the manner
in which capacitance is taken into account; the overhead transmission lines are classified as:
Short transmission lines: When the length of an overhead transmission line is up to about
50 km and the line voltage is comparatively low (< 20 kV), it is usually considered as a short
transmission line. Due to smaller length and lower voltage, the capacitance effects are small
and hence can be neglected. Therefore, while studying the performance of a short
transmission line, only resistance and inductance of the. line are taken into account.
Medium transmission lines: When the length of an overhead transmission line is about 50-
150 km and the line voltage is moderately high (>20 kV < 100 kV), it is considered as a
medium transmission line. Due to sufficient length and voltage of the line, the capacitance
effects are taken into account. For purposes of calculations, the distributed capacitance of the
line is divided and lumped in the form of condensers shunted across the line at one or more
points.
Long transmission lines: When the length of an overhead transmission line is more than 150
km and line voltage is very high (> 100 kV), it is considered as a long transmission line. For
the treatment of such a line, the line constants are considered uniformly distributed over the
whole length of the line and rigorous methods are employed for solution
It may be emphasized here that exact solution of any transmission line must consider the fact
that the constants of the line are not lumped but are distributed uniformly throughout the
length of the line. However, reasonable accuracy can be obtained by considering these
constants as lumped for short and medium transmission lines.
Medium Transmission Lines (Nominal ∏- Method):
In this method, capacitance of each conductor (i.e., line to neutral) is divided into two halves;
one half being lumped at the sending end and the other half at the receiving end as shown in

ASCET/EEE/B.TECH/20G25A0253 21
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 22
Power System & Simulation Lab

below figure. It is obvious that capacitance at the sending end has no effect on the line drop.
However, it’s charging current must be added to fine current in order to obtain the total
sending end current.
Let ir= load current per phase
R= resistance per phase
XL = inductive reactance per phase
C= capacitance per phase
cos ØR= receiving end power factor (lagging)
VS = sending end voltage per phase

The phasor diagram for the circuit is shown in above figure. Taking the receiving end voltage
as the reference phasor, we have

Load current

Charging current at load end is

Line current

ASCET/EEE/B.TECH/20G25A0253 23
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 24
Power System & Simulation Lab

Sending end voltage

Charging current at sending end

Sending end current

Calculation of inductance for 3 phase transmission line system:


Inductance/phase/m = 10-7 [0.5+2log(gmd/r)]
Where gmd=geometric mean distance between conductors=(d12*d21*d31)^1/3
R=radius opf conductor
Calculation of capacitance for 3 phase transmission line system:
Capacitance/phase/m = 3.14*8.854*10-12/log(gmd/r)
Where gmd=geometric mean distance between conductors=(d12*d21*d31)^1/3
R=radius opf conductor

ABCD PARAMETERS:

Voltage regulation: When a transmission line is carrying current, there is a voltage drop in
the line due to resistance and Inductance of the line. The result is that receiving end voltage
(VR) of the line is generally less than the sending end voltage (VS). This voltage drop (VS -
VR) in the line is expressed as a percentage of receiving end voltage VR and is called voltage
regulation.

Transmission efficiency: The power obtained at the receiving end of a transmission line is
generally less than the sending end power due to losses in the line resistance.

ASCET/EEE/B.TECH/20G25A0253 25
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 26
Power System & Simulation Lab

FORMULAE:
Single phase Inductance = 10-7[1+4log(d/r)]
Capacitance = 3,14*8.854*10-12/log(d/r)
Three phase Inductance = 10-7 [0.5+2log(d/r)]
Capacitance = 2*3.14*8.854*10-12/log(d/r)
Dequivatent = [d1*d2*d3]1/3
were
d = spacing of conductors
r= radius of conductors

ALGORITHM:
STEP 1: Find that the given transmission line is single phase or three phase
STEP 2: If it is single phase, get the value of distance between the conductors.
STEP 3: Get the radius of the conductor
STEP 4: Using the appropriate formula, find inductance and capacitance. __
STEP 5: If the given system is three phase, classify whether it is symmetrical or
unsymmetrical
STEP 6: If symmetrical, get the distance between the conductors and radius of
the conductor.
STEP 7: Using the appropriate formula, find inductance and capacitance
STEP 8: If unsymmetrical, get the distance between the conductors and radius of the
conductor. Using the appropriate formula, find inductance and capacitance.
PROBLEM:
Determine the sending end voltage, current, power & power factor for a 160km
section of 3phase line delivering 5OMW at 132kV and P.F 0.8 lag. Also find the efficiency
and regulation of the line. Resistance per line 0.1557o0hm per km, spacing 3.7m, 6.475m,
7.4m transposed. Evaluate the A, B, C, D parameters also. Diameter is 1.956cm. Write and
execute a MATLAB program and also verify the output with the manual calculation results.

ASCET/EEE/B.TECH/20G25A0253 27
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 28
Power System & Simulation Lab

PROGRAM:

%analysis of transmission tine


clc;
display(‘………………………………………………………………....’)
display(‘output for medium transmission line parameters-19G21A0235’)
display(‘…………………………………………………………………’)
I=input(‘enter length of medium 3 phase transmission line In km:');
f=input(‘enter operating frequency In Hz:');
d=input(‘enter diameter of conductor in meters:');
r=d/2;
d12=input(‘enter distance between lines 1&2:');
d23=input(‘enter distance between lines 2&3:');
d31=input(‘enter distance between lines 3&1:');
R=input(‘enter resistance per km of a line:');
RT=R*I;
pr=input(‘enter receiving power in MW:');
vr=input(‘enter receiving end line voltage in kv:');
pf=input(‘enter power factor at output:');
w=2*pi*f;
L=2*(10^-7)*I*(10^3)*log((nthroot((d12*d23*d31),3))/(0.7788*r))
C=(2*pi*I*(10^3)*8.854*(10^-12))/log((nthroot((d12*d23*d31),3))/r)
Z=complex(RT,(w*L))
Y=complex(0,(w*C))
A=(1+((Y*Z)/2))
B=Z
C=Y*(1+((Y*Z)/4))
D=A
irold=(pr*(10^3))/((nthroot(3,2)*vr*pf))
ir=irold*(pf-1i*sin(acos(pf)))

ASCET/EEE/B.TECH/20G25A0253 29
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 30
Power System & Simulation Lab

vp=(vr*10^3)/(nthroot(3,2))
vs=(A*vp)+(B*ir)
is=(C*vp)+(D*ir)
anglevs=angle(vs)*(180/pi)
angleis=angle(is)*(180/pi)
inputpf=cos(angle(vs)-angle(is))
ps=3*abs(vs)*abs(is)*inputpf
efficiency=((pr*10^6)/ps)*100
regulation=((abs(vs)-abs(vp))/abs(vp))*100

simulation output:
………………………………………………………………....
output for medium transmission line parameters-19G21A0235
………………………………………………………………....
enter length of medium 3 phase transmission line in km:160
enter operating frequency in Hz:50
enter diameter of conductor in meters:1.956*10^-2
enter distance between lines 1&2:3.7
enter distance between lines 2&3:6.475
enter distance between lines 3&1:7.4
enter resistance per km of a line:0.1557
enter receiving power in MW:50
enter receiving end line voltage in kv:132
enter power factor at output:0.8
L=0.2113
C =1.4010e-06
Z=24.9120 +66.3840i
Y=0 +4.4014e-004i
A=0.9854 + 0.0055i
B =24.9120 +66.3840i

ASCET/EEE/B.TECH/20G25A0253 31
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 32
Power System & Simulation Lab

C = -1.2065e-06 +4.3692e-04i
D =0.9854 + 0.0055i
irold =273.3666
ir =2.1869e+02 -1.6402e+02i
vp =7.6210e+04
vs =9.1433e+004 +1.0849e+004i
is =2.1631e+02 -1.2713e+02i
anglevs =6.7671
angleis =-30.4436
inputpf =0.7964
ps =5.5195e+007
efficiency =90.5882
regulation =20.8167

RESULT:
Thus, the computation of line parameters of a given power transmission system is
done by using MATLAB and the output is verified.

ASCET/EEE/B.TECH/20G25A0253 33
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 34
Power System & Simulation Lab

4. FORMATION OF BUS ADMITTANCE MATRIX

AIM:
To form the bus admittance matrix (Y-bus) for a given power system using MATLAB
program.

SOFTWARE REQUIRED:
MATLAB software package.

THEORY:
The Ybus / Zbus matrix constitutes the models of the passive portions of the power
network. Ybus matrix is often used in solving load flow problems. It has gained widespread
applications owing to its simplicity of data preparation and the ease with which the bus
admittance matrix can be formed and modified for network changes. Of course, sparsity is
one of its greatest advantages as it heavily reduces computer memory and time requirements.
In short circuit analysis, the generator and transformer impedances must also be taken into
account. In contingency analysis, the shunt elements are neglected, while forming the Z-bus
matrix, which is used to compute the outage distribution factors.
This can be easily obtained by inverting the Y-bus matrix formed by inspection
method or by analytical method. The impedance matrix is a full matrix and is most useful for
short circuit studies. Initially, the Y-bus matrix is formed by inspection method by
considering line data only. After forming the Y-bus matrix, the modified Y-bus matrix is
formed by adding the generator and transformer admittances to the respective diagonal
elements and is inverted to form the Zbus matrix.
The performance equation for a n-bus system in terms of admittance matrix can be
written as,

The admittances Y11, Y12,... Y1n are called the self-admittances at the nodes and all other
admittances are called the mutual admittances of the nodes.

ASCET/EEE/B.TECH/20G25A0253 35
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 36
Power System & Simulation Lab

Formulae Used:

Main diagonal element in Y-bus matrix =


Where, Bij is the half line shunt admittance in mho, Yij is the series admittance in mho.
Off-diagonal element in Y-bus matrix, Yij= -Yij
Where, Yij is the series admittance in mho.

ALGORITHM:
STEP1: Read all the data namely R and X for the system
STEP2: Calculate the mutual or transfer reactance for the reactance between i and jand
i=j=1, 2, 3, 4...
STEP3: Calculate the self- admittance or point admittance bus i=1, 2, 3, 4...
STEP4: Output the Y-bus matrix.
STEP5: Print the result

PROBLEM:
Form the Ybus matrix for the power system connected with bus numbers and
impedance and self-admittance given for buses as
Bus number impedance Bus Self-admittance
1-2 0.06+i0.18 1 10.05
1-3 0.02+i0.06 2 10.06
2-3 0.04+i0.12 3 10.05

PROGRAM:
%formation of bus admittance matrix
clc;
display(‘Formation of bus admittance matrix (19G21A0235)’);
nbranch=input(‘enter the number of branches in the system =');
display(‘enter line data’);
for n=1:1:nbranch
fb=input(‘enter from bus=');
tb=input(‘enter to bus=');

ASCET/EEE/B.TECH/20G25A0253 37
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 38
Power System & Simulation Lab

r=input(‘enter the value of resistance’);


x=input(‘enter the value of reactance=’);
b=input(‘enter the value of line charging admittance(b/2)=’);
z=r+ 1i*x;
y=1/z;
Ldata(n,:)=[fb tb r x b y];
end
fb=Ldata(:,1);
tb=Ldata(:,2);
r=Ldata(:,3);
x=Ldata(:,4);
b=Ldata(:,5);
y=Ldata(:,6);
b=1i*b;
nbus=max(max(fb),max(tb));
Y=zeros(nbus,nbus);
for k=1:nbus
bs=input(‘enter the bus at which ground reactor is connected=');
x0=input(‘enter value of ground reactor=');
if x0~=0
y0=1/(1i*x0);
else
y0=0;
end
Yreact(k,:)=[bs y0];
end
bs=Yreact(:, 1);
y0=Yreact(:,2);
%off diagnal element
for k=1:nbranch

ASCET/EEE/B.TECH/20G25A0253 39
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 40
Power System & Simulation Lab

Y (fb(k),tb(k))=Y(fb(k),tb(k))-y(k);
Y(tb(k),f(k))=Y(fb(k),tb(k));
end
% diagonal element
for m=1:nbus
for n=1:nbranch
iffb(n)==m
Y(m,m)=Y(m,m)+y(n)+b(n);
elseiftb(n)==m
Y(m,m)=Y(m,m)+y(n)+b(n);
end
end

for k=1:nbus
Y(k,k)=Y(k,k)+y0(k);
end
Yb=Y

OUTPUT SIMULATION:
Formation of bus admittance matrix(19G21A0235)
enter the number of branches in the system =3
enter line data
enter from bus=1
enter to bus=2
enter the value of resistance=0.06
enter the value of reactance=0.18
enter the value of line charging admittance(b/2)=0.05
enter from bus=1
enter to bus=3
enter the value of resistance=0.02

ASCET/EEE/B.TECH/20G25A0253 41
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 42
Power System & Simulation Lab

enter the value of reactance=0.06


enter the value of line charging admittance(b/2)=0.06
enter from bus=2
enter to bus=3
enter the value of resistance=0.04
enter the value of reactance=0.12
enter the value of line charging admittance(b/2)=0.05
enter the bus at which ground reactor is connected=1
enter value of ground reactor=0
enter the bus at which ground reactor is connected=2
enter value of ground reactor=0
enter the bus at which ground reactor is connected=3
enter value of ground reactor=0
Yb =
6.6667 -19.8900) -1.6667 + 5.0000i -5.0000 +15.0000i
-1.6667 + 5.0000; 4.1667 -12.4000i -2.5000 + 7.5000i
-5.0000 +15.0000i -2.5000+ 7.50001 7.5000 -22.3900i

RESULT:
Thus, the formation of bus admittance matrix is done by using MATLAB and the
output is verified for the given power system.

ASCET/EEE/B.TECH/20G25A0253 43
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 44
Power System & Simulation Lab

5. TRANSIENT STABILITY ANALYSIS OF A SINGLE MACHINE


INFINITE BUS SYSTEM
AIM:
To analyse the transient stability of a single machine infinite bus system by point-by-
point method using MATLAB.

SOFTWARE REQUIRED:
MATLAB software package.

ALGORITHM:
STEP 1: Evaluate the accelerating power Pa.
STEP 2: From the swing equation d2δ/dt 2 = α (0) = Pa(Qt)/M Where δ is the
acceleration factor. Evaluation δ.
STEP 3: The change in angular velocity for the first interval is calculated.
STEP 4: The change in rotor angle for the first interval is also calculated.
STEP 5: If the discontinuity occurs due to removal of the fault or due to switching
operation there are three possibilities.
a. The discontinuity occurs at the beginning of i" interval.
b. The discontinuity occurs at the middle of the i interval.
STEP 6: To evaluate Pa when under first situation, one should use the value
corresponding to average value of Q-accelerating power i.e., power before and
after clearing the fault.
STEP 7: To evaluate Pa under the situation, a weighted average value of Pa before and
after the discontinuity may be used.
STEP 8: Thus, equating for nth interval can be written as
Pa(n-1) = Ps — Pe (n-1)
δn = 5n-1 + del δn
There are used for plotting the curve.
STEP 9: To evaluate Pa under second situation no procedure is required. It is taken as
the value at the beginning of the interval.

ASCET/EEE/B.TECH/20G25A0253 45
Power System & Simulation Lab

Model graph:

Formulae:
Reactive power Qe = sin(cos-1(p.f))

Stator Current It =

Voltage behind transient condition

Angular separation between E1 and EB

Prefault Operation:

ASCET/EEE/B.TECH/20G25A0253 46
Power System & Simulation Lab

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

PROBLEM:
Consider a system which consists of generator having 2 rating of 50 MMVA & H=2.7 -
MJ/MVA at rated speed, E=1.05, V=1, Xd’ =0.2, X1=%2=0.4 pu. The generator supplies 50
MW to the infinite bus when a 3phase fault occurs at middle of line 2.
i) Plot swing curve for a sustained fault up to 0.5 SEC.
ii) Plot the swing curve if the fault is cleared in 0.1 sec by simultaneous opening of
breakers at both ends of line 2.
iii) Find the critical clearing angle & clearing time.
iv) Write and execute a MATLAB program and also verify the output with the
manual calculation results.

ASCET/EEE/B.TECH/20G25A0253 47
Power System & Simulation Lab

During Fault Condition:


Pe= Psi = 0
Find out X from the equivalent circuit during fault condition

Post fault Condition:


Find out X from the equivalent circuit during post fault condition

ASCET/EEE/B.TECH/20G25A0253 48
Power System & Simulation Lab

PROGRAM:

clc;
display(‘….Transient stability output <19G21A0235>….’);
t=0;
tf=0;
f=input('Enter the input frequency’);
s= input(‘Enter the machine rating’);
ang(1)=input(‘Enter the initial angle’);
h=input(‘Enter the moment of inertia constant’);
tfinal=input('Enter the end time’);
tstep=input('Enter the change in time’);
tc=input('enter the clearing time’);
pm=input('enter the power transfer’);
pmaxbf=input(‘enter the prefault power’);
pmaxdf=input(‘enter the power during fault’);
pmaxaf=input(‘enter the postfault power’);
m= (s*h)/(180*f);
delta=ang(1)*pi/180;
i=2;
ddelta=0;
time(1)=0;
while t<tfinal
if t==tf
paminus=pm-pmaxbf*sin(delta);
paplus=pm-pmaxdf*sin(delta);
paav=(paminus+tpaplus)/2;
pa=paav;
end
if(t==tc)

ASCET/EEE/B.TECH/20G25A0253 49
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 50
Power System & Simulation Lab

paminus=pm-pmaxdf*sin(delta);
paplus=pm-pmaxaf*sin(delta);
paav=(paminus+paplus)/2;
pa=paav;
end
if(t>tf&t<tc)
pa=pm-pmaxdf*sin(delta);
end
if(t>tc)
pa=pm-pmaxaf*sin(delta);
end
t
pa
ddelta=ddelta+(tstep*tstep*pa/m);
delta=((delta*180/pi)+ddelta)*(pi/180);
deltadeg=delta*180/pi;
t=t+tstep;
pause
time(i)=t;
ang(i)=deltadeg;
i=i+1;
end
axis([0 0.6 0 180])
plot(time,ang,’ko-‘)

EXECUTION:
….Transient stability output <19G21A0235>….
Enter the input frequency 50
Enter the machine rating 50
Enter the initial angle 4/11

ASCET/EEE/B.TECH/20G25A0253 51
Power System & Simulation Lab

OUTPUT GRAPH:

ASCET/EEE/B.TECH/20G25A0253 52
Power System & Simulation Lab

Enter the moment of inertia constant 2.7


Enter the end line 0.5
Enter the change in time 0.01
Enter the clearing time 0.1
Enter the power transfer 1
Enter the prefault power 11/4
Enter the power during fault 11/8
Enter the postfault power 11/6
t=
0
pa=
0.9869
t=
0.1000
pa=
0.9911
t=
0.0200
pa=
0.9908
t=
0.0300
pa=
0.9903

RESULT:
Thus, the transient stability analysis for single machine connected to infinite bus and
fault clearance in different time was studied.

ASCET/EEE/B.TECH/20G25A0253 53
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 54
Power System & Simulation Lab

6. LOAD FLOW SOLUTION BY USING GAUSS-SEIDEL METHOD

AIM:
To find load flow analysis using Gauss-Seidel method In MATLAB

SOFTWARE REQUIRED:
MATLAB software package.

ALGORITHM:
STEP1: The slack bus voltage magnitude and angle are measured usually V1=1 p.u
with the load profile known at each bus, we allocate Pl and Qi to all generating
Stations with this step, bus Injections (Pi+Qi) are known at all buses other than
the slack bus.
STEP2: Assembly of bus admittance matrix: with the ling and shunt admittance data
stored in the computer, Y bus Is assembled by using the algorithm developed
earlier, Alternatively bus ls assembled using Y bus=ATYA where the Input is
in the form of primitive admittance matrix Y and singular connection bus
incidence matrix A
STEP3: Iterative computation of bus voltages (Vi, 1=1,2...n) to start Iteration a set of
initial values is assumed, since in a power system the voltage spread is not too
wide, it is normal practice to use a flat voltage start, le. Initially all voltages
are set equal to (1+j0) expect slack bus voltage which is fixed . this reduced
the n equations in complex number which are to solved iteratively for finding
complex voltages V2,V3,...Vn. IF complex no options are net available in a
computer, the equation is real unknown js

We also defined,

Now for (r+1) the iteration the voltage becomes,


Vi(r+1) =[Ai /(Vir) * -Σi-1 (Bik Vk(r+1) -ΣbikVk r]
The iterative process is continued till the change in magnitude of bus voltage
|∆Vi(r+1) Between two consecutive is less than a certain for all bus voltages
i.e.,
∆Vi (r+1) = | Vi(r+1)-Vi (r)| ≤εi=2......,n

ASCET/EEE/B.TECH/20G25A0253 55
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 56
Power System & Simulation Lab

Also we see if
|∆Vi | min $ [AVE [max I=2 …………. n
If not, we fix |∆Vi |at one of extreme values i.e.
|∆Vi| min if |∆ Vi | ≤ |∆ Vi| min or |∆Vi | max if |∆Vi| ≤ |∆Vi|
STEP4: Computation of slack bus power; substitution of all bus voltages computed in
step3 with Vi and I=1 yield real and reactive power at slack bus i.e.,
S1=P1+jQ1
STEPS: Computation of line flows; this is the last step in the load flow analysis where
in the power flows on the various lines of the network are computed. This also
enables us to check whether any line overloaded. Consider the line connecting
buses | and k. The line and transformer at each end can be represented by a
circuit with series admittance Yik and to shunt admittances Yiko. As the
current fed by bus line to the line can be expressed as
lik = lik1+lik0= (Vi-Vk)Yik+ Viyik0

THEORY:
Load flow analysis is the most frequently performed system study by electric utilities.
This analysis is performed on a symmetrical steady-state operating condition of a power
system under ‘normal’ mode of operation and aims at obtaining bus voltages and
line/transformer flows for a given load condition. This information is essential both for long
term planning and next day operational planning. In long term planning, load flow analysis
helps in investigating the effectiveness of alternative plans and choosing the ‘best’ plan for
system expansion to meet the projected operating state. In operational planning, it helps in
choosing the ‘best’ unit commitment plan and generation schedules to run the system
efficiently for them next day’s load condition without violating the bus voltage and line flow
operating limits.
The Gauss seidel method is an iterative algorithm for solving a set of non- linear
algebraic equations. The relationship between network bus voltages and currents may be
represented by either loop equations or node equations. Node equations are normally
preferred because the number of independent node equation is smaller than the number of
independent loop equations.
The network equations in terms of the bus admittance matrix can be written as,

For an bus system, the above performance equation can be expanded as,

ASCET/EEE/B.TECH/20G25A0253 57
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 58
Power System & Simulation Lab

where n is the total number of nodes.


Vp is the phasor voltage to ground at node p.
Ip is the phasor current flowing into the network at node p.

At the pth bus, current injection:

At bus p , we can write

Hence, the current at any node p is related to P, Q and V as follows:

(For any bus p except slack bus s)


Substituting for Ip, in Equation Vp,

Ip has been substituted by the real and reactive powers because normally in a power system
these quantities are specified.

ASCET/EEE/B.TECH/20G25A0253 59
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 60
Power System & Simulation Lab

PROBLEM:
Figure shows the one-line diagram of a three-bus power system with generation at
bus1.The magnitude of voltage at bus1 is adjusted to 1.05 per unit. The scheduled loads at
buses 2 & 3 are as marked on the diagram. Line impedances are marked in per unit on a
100MVA base and the line charging susceptance are neglected.
i) Using Gauss-Seidel method, determine the phasor values of the voltage at the
load buses 2 & 3 (PQ buses) accurate to 4 decimal places. .
ii) Write and execute a MATLAB program and also verify the output with the
manual calculation results.

PROGRAM:

clc;
display (‘----- Formation of Ybus By G-S Method<19G21A0235>-----‘);
n= Input ('Enter the number of buses ');
fprintf (‘Enter your choice’);
p= input ('1. impedance, 2. admittance’);
if (p==1)
for q= 1:n
for r=q+1:n
fprintf ('Enter the impedance value between %d-%d',q,r);
z(q,r)=input(':');

ASCET/EEE/B.TECH/20G25A0253 61
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 62
Power System & Simulation Lab

if (z(q,r)==0)
y(q,r)=0;
else
y (q,r)=inv(z(q,r));
end
y(r,q)= y(q,r);
end
end
elseif (p==2)
for a= 1:n
for b=a+1:n
fprintf (‘Enter the admittance value between %d-%d',a,b);
y(a,b)=input(":');
y(b,a)= y(a,b);
end
end
else
fprintf ('enter the correct choice’);
end
ybus=zeros(n,n);
for a=1:n
for b=1:n
if (a==b)
for c=1:n
ybus(a,a)= ybus(a,a)+ y(a,c);
end
else
ybus(a,b)=-y(b,a);
end
end

ASCET/EEE/B.TECH/20G25A0253 63
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 64
Power System & Simulation Lab

end
ybus
%......BUSDATA……%
busdata =[1 1 1.05 0 0 0 0 0 0 0;
2 3 1.00 0 0 0 256.6 110.2 0 0;
3 3 1.00 0 0 0 138.6 45.2 0 0];
bus = busdata(:,1);
type = busdata(:,2);
V = busdata(:,3);
th = busdata(:,4);
GenMW = busdata(:,5);
GenMVAR = busdata(:,6);
LoadMW = busdata(:,7);
LoadMVAR = busdata(:,8);
Qmin = busdata(:,9);
Qmax = busdata(:,10);
nbus = max(bus);
P = GenMW - LoadMW;
Q = GenMVAR - LoadMVAR;
Vprev =V;
toler = 1;
iteration = 1;
disp(' Bus number | 1.Slack 2.PQ 3.PV | V| angle|Pg | Qg | PL | QL | Qmin | Qmax ');
busdata
%----- VOLTAGE CALCULATION -—---%
while (toler> 0.1)
for i = 2:nbus
sumyv = 0;
for k = 1:nbus
if i~=k

ASCET/EEE/B.TECH/20G25A0253 65
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 66
Power System & Simulation Lab

sumyv = sumyv + ybus(i,k)* V(k);


end
end
if type(i) == 2
Q(i) = -imag(conj(V(i))*(sumyv + ybus(i,i)*V(i)));
if (Q(i) >Qmax(i)) || (Q(i) <Qmin(i))
if Q(i) <Qmin(i)
Q(i) = Qmin(i);
else
Q(i) = Qmax(i);
end
type(i) = 3;
end
end
V(i) = (1/ybus(i,i))* ((P(i)-j* Q(i))/conj(V(i)) - sumyv);
if type(i) == 2
V(i) = pol2rect(abs(Vprev(i)), angle(V(i)));
end
end
iteration = iteration + 1;
toler = max(abs(abs(V) - abs(Vprev)));
Vprev=V;
end
iteration
Vmag = abs(V)
Ang = 180/pi*angle(V)
sum=0;
%----- REAL AND REACTIVE POWER CALCULATION -----%
For i=1:nbus
If i==1

ASCET/EEE/B.TECH/20G25A0253 67
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 68
Power System & Simulation Lab

for f=1:nbus
sum=sum+(ybus(i,f)*V(f));
real_power(i)=(real(V(i)*sum))*100;
reactive_power(i)=-(imag(V(i)*sum))*100;
end
else
end
end
real_power
reactive_power

SIMULATION OUTPUT:
----- Formation Of Ybus By G-S Method<19G21A0235>-----
Enter the number of buses 3
Enter your choice1. impedance, 2. Admittance 1
Enter the impedance value between 1-2:0.02+j*0.04
Enter the impedance value between 1-3:0.01+j*0.03
Enter the impedance value between 2-3:0.0125+1j*0.025
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

Bus number | 1.Slack 2.PQ 3.PV | V| angle|Pg | Qg | PL | QL | Qmin | Qmax


busdata =
1.0000 1.0000 1.0500 0 0 0 0 0 0 0
2.0000 3.0000 1.0000 0 0 0 256.6000 110.2000 0 0
3.0000 3.0000 1.0000 0 0 0 138.6000 45.2000 0 0
iteration =
22
Vmag=
1.0500

ASCET/EEE/B.TECH/20G25A0253 69
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 70
Power System & Simulation Lab

2.0926
1.9439
Ang=
0
-71.8342
-83.9538
real_power =
1.1570e+04
reactive_power =
-619.9159

RESULT:
Thus the load flow analysis using is performed by Gauss-Seidal method and a
program is developed using MATLAB to find the solution of load flow for given power
system and the output is verified.

ASCET/EEE/B.TECH/20G25A0253 71
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 72
Power System & Simulation Lab

7. LOAD FLOW SOLUTION BY USING NEWTON-RAPHSON


METHOD
AIM:
To develop in MATLAB to find the solution of power flows using Newton Raphson
method.

SOFTWARE REQUIRED:
MATLAB software package

ALGORITHM:
STEP 1: Assume a flat profile 1+j0 for all buses except the slack bus in the Specified
voltage and it is not modified in any iteration.
STEP 2: Assume a suitable value of ε called convergence criterion. Hence € is a
specified change in the residue that is used to compare the actual residue at the
end of each iteration.
STEP 3: Set the iteration count K=0 and assumed voltage profile of the buses are
denoted as V10,V20......vn0.
STEP 4: Set the bus count p=1
STEP 5: Check for slack bus. If it is a slack bus then go to step 13. Otherwise go to
next step.
STEP 6: Calculate the real & reactive power of bus p using the following equation,

STEP 7: Calculate the change In real power, ∆Pk=Pp spec -Ppk Where, Pp
spec=specified real power of bus p.
STEP 8: Check for generator bus. If it Is a generator bus gob to next step otherwise go
to step 12.
STEP 9: Check for generator bus. If it is a generator reactive power limits Violation of
generator buses. For this compare the calculated reactive power Q with
specified limits. If the limits are violated go to step11. Otherwise go to next
step.
STEP 10: If the calculated reactive power is within the specified limit then consider this
bus as generator bus. Now calculate the voltage residue using the equation

ASCET/EEE/B.TECH/20G25A0253 73
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 74
Power System & Simulation Lab

Where |Vp spec| = specified voltage magnitude for generation bus. Then go to
step 13.
STEP 11: If the reactive power limit violated the treat this bus as a load bus. Now the
specified reactive power for this bus will correspond to Limit violated.
STEP 12: Calculate the change in reactive power for load bus change in reactive power,

STEP 13: Repeat the step 5 to 12 until all residues are calculated for increment the bus
count n. by 1 to 5 steps until the bus count is n.
STEP 14: Determine the largest of the absolute value of the residue (i.e.) Find the largest
value among ∆Ppk , ∆Qpk, or |∆Vpk|2
STEP 15: Compare ∆E and E, if ∆E<E then go to step 20. If ∆E>E go to next Step.
STEP 16: Determine the element, the load flow equation using th iteration Value.
STEP 17: Calculate the increment in real and reactive part of voltage ∆epk and ∆fpk by
solving the matrix equation B=JC.
STEP 18: Calculate the new bus voltage.
STEP 19: Advance the iteration count k=k+1 and go to step 4.
STEP 20: Calculate the line flows.

Theory:
Load flow study in power system parlance is the steady state solution of the power
system network. The main information obtained from this study comprises the magnitudes
and phase angles of load bus voltages, reactive powers at generator buses, real and reactive
power flow on transmission lines, other variables being specified. This information is
essential for the continuous monitoring of current state of the system and for analysing the
effectiveness of alternative plans for future system expansion to meet increased load demand.
Newton-Raphson method is an iterative method that approximates the set of non
linear ay simultaneous equations to a set of linear simultaneous equations using Taylor’s
series expansion and the terms are limited to first approximation. The rate of convergence is
fast as compared to the FDLF program and also it is suitable for large size system. So, we go
for N-R method. The non-linear equations governing the power system network are,

where Ip is the current injected into bus p.


The complex power in pe bus is given by

ASCET/EEE/B.TECH/20G25A0253 75
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 76
Power System & Simulation Lab

In polar co-ordinates, the power on pth bus Is given as,

Separating the Real and Imaginary parts we get,

The Newton —Raphson method requires that a set of linear equations be formed expressing
the relationship between the changes in real and reactive powers and the components of the
bus voltages as follows:

where, the coefficient matrix is known as Jacobian matrix.


In the above equation, bus 1 is assumed to be the slack bus. The Jacobian matrix gives the
linearized relationship between small changes in voltage angle ∆δi(r) and voltage magnitude
|∆Vi(r)| the small changes in real and reactive power ∆Pi(r) and ∆Qi(r) Elements of the Jacobian
matrices are the partial derivatives of (2) and (3) evaluated at ∆δi(r) and |∆Vi(r)|
The above relationship can be written in a compact form as,

ASCET/EEE/B.TECH/20G25A0253 77
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 78
Power System & Simulation Lab

The elements of Jacobian matrix are defined as,

on when solved for ∆δ, ∆V gives the correction to be applied to |V| and δ, i.e.

Next, we get a new set of linear equations evaluated at (r+1)th iteration and the process is
repeated. Convergence is tested by the power mismatch criteria. This method converges to
high accuracy nearly always in 2 to 5 iterations from a flat start (|V| = 1 p.u. and θ=0) for all
buses
where |V|, θ are unknown, independent of system size.
At PV bus at the end of an iteration and if it violates the limits, the PV bus is switched to a
PQ bus. When Q is within limits, then it is switched back to PV bus.

PROBLEM:
Figure shows the one-line diagram of a three-bus power system with generators at buses 1& 3.
The magnitude of voltage at bus1 is adjusted to 1.05 per unit. The magnitude of voltage at
bus 3 is fixed at 1.04pu with a real power generation of 200 MW. A load consists of 400 MW
and 250 MVAR is taken from bus2. Line impedances are marked in per unit on a 100MVA
base and the line charging susceptance are neglected.

i) Obtain the load flow solution by Newton-Raphson method.


ii) Write and execute a MATLAB program and also verify the output with the
manual calculation results.

ASCET/EEE/B.TECH/20G25A0253 79
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 80
Power System & Simulation Lab

PROGRAM:

clc;
display(‘……………………………………………………………………………..’);
display(‘Load Flow Solution By Using Newton-Raphson Method<19G21A0235>’);
display(‘……………………………………………………………………………..’);
v=[1.05;1.0;1.04];
d=[0;0;0];
ps=[-4;2.0);
qs=-2.5;
n= input('Enter the number of buses ‘);
fprintf('Enter your choice’);
p= input ('1, impedance, 2. admittance’);
if (p==1)
for q= 1:n
for r=q+1:n
fprintf('Enter the impedance value between %d-%d',q,r);
z(q,r)=input(':');
if (z(q,r)==0)
y(q,r)=0;
else
y(q,r)=inv(z(q,r));
end
y(r,q)= y(q,r);
end
end
elseif (p==2)
for a= 1:n
for b=a+1:n
fprintf('Enter the admittance value between %d-%d',a,b);

ASCET/EEE/B.TECH/20G25A0253 81
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 82
Power System & Simulation Lab

y(a,b)=input(':');
y(b,a)= y(a,b);
end
end
else
fprintf(‘enter the correct choice');
end
ybus=zeros(n,n);
for a=1:n
for b=1:n
if (a==b)
for c=1:n
ybus(a,a)= ybus(a,a)+ y(a,c);
end
else
ybus(a,b)=-y(b,a);
end
end
end
ybus
y=abs(ybus);
t=angle(ybus);
iter=0;
pwracur=0.00025; % Power accuracy
dc=10; % Set the maximum power residual to a high value
while max(abs(dc))>pwracur
iter=iter+1;
p=[v(2)*v(1)*y(2, 1) *cos(t(2,1)-d(2)+d(1))+v(2)^2 *y(2,2) *cos(t(2,2))+v(2)*v(3)
*y(2,3)*cos(t(2,3)-d(2)+d(3));
v(3)*v(1)* y(3,1)* cos(t(3,1)-d(3)+d(1))+v(3)^2*y(3,3)* cos(t(3,3))+v(3)*v(2)*y(3,2)*
cos(t(3,2)-d(3)+d(2))];

ASCET/EEE/B.TECH/20G25A0253 83
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 84
Power System & Simulation Lab

q=-v(2)*v(1)*y(2,1)*sin(t(2,1)-d(2)+d(1))-v(2)^2*y(2,2)*sin(t(2,2))-
v(2)*v(3)*y(2,3)*sin(t(2,3)-d(2)+d(3));
j(1,1)=v(2)*v(1)*y(2,1)*sin(t(2,1)-d(2)+d(1))+v(2)*v(3)*y(2,3)*sin(t(2,3)-d(2)+d(3));
j(1,2)=-v(2)*v(3)*y(2,3)*sin(t(2,3)-d(2)+d(3));
j(1,3)=v(1)*y(2,1)*cos(t(2,1)-d(2)+d(1))+2*v(2)*y(2,2)*cos(t(2,2))+v(3)*y(2,3)*cos(t(2,3)-
d(2)+d(3));
j(2,1)=-v(3)*v(2)*y(3,2)*sin(t(3,2)-d(3)+d(2));
j(2,2)=v(3)*v(1)*y(3,1)*sin(t(3,2)-d(3)+d(1))+v(3)*v(2)*y(3,2)*sin(t(3,2)-d(3)+d(2));
j(2,3)=v(3)*y(2,3)*cos(t(3,2)-d(3)+d(2));
j(3,1)=v(2)*v(1)*y(2,1)*cos(t(2,1)-d(2)+d(1))+v(2)*v(3)*y(2,3)*cos(t(2,3)-d(2)+d(3));
j(3,2)=-v(2)*v(3)*y(2,3)*cos(t(3,2)-d(2)+d(3));
j(3,3)=-v(1)*y(2,1)*sin(t(2,1)-d(2)+d(1))-2*v(2)*y(2,2)*sin(t(2,2))-v(3)*y(2,3)*sin(t(2,3)-
d(2)+d(3));
dp=ps-p;
dq=qs-q;
dc=[dp;dq]
j
dx=j\dc
d(2)=d(2)+dx(1);
d(3)=d(3)+dx(2);
v(2)=v(2)+dx(3);
v,d,delta=180/pi*d;
end
p1=v(1)^2*y(1,1)*cos(t(1,1))+v(1)*v(2)*y(1,2)*cos(t(1,2)-
d(1)+d(2))+v(1)*v(3)*y(1,3)*cos(t(1,3)-d(1)+d(3));

q1=-v(1)^2*y(1,1)*sin(t(1,1))-v(1)*v(2)*y(1,2)*sin(t(1,2)-d(1)+d(2))-
v(1)*v(3)*y(1,3)*sin(t(1,3)-d(1)+d(3));

q3=-v(3)*v(1)*y(3,1)*sin(t(3,1)-d(3)+d(1))-v(3)*v(2)*y(3,2)*sin(t(3,2)-d(3)+d(2))-
v(3)^2*y(3,3)*sin(t(3,3));

ASCET/EEE/B.TECH/20G25A0253 85
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 86
Power System & Simulation Lab

EXECUTION:
………………………………………………………………………………
Load Flow Solution By Using Newton-Raphson Method<19G21A0235>
………………………………………………………………………………
Enter the number of buses 3
Enter your choice1. impedance, 2. admittance1
Enter the impedance value between 1-2:0.02+0.04j
Enter the impedance value between 1-3:0.01+0,03)
Enter the impedance value between 2-3:0.0125+0.025}
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
iter =
1
dc=
-2.8600
1.4384
-0.2200
j=
54.2800 -33.2800 24.8600
-33.2800 64.1664 -16.6400
-27.1400 16.6400 49.7200
dx=
-0.0455
-0.0080
-0.0265
v=
1.0500
0.9735

ASCET/EEE/B.TECH/20G25A0253 87
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 88
Power System & Simulation Lab

1.0400
d=
0
-0.0455
0.0080
iter =
2
dc =
-0.0992
0.0367
-0.0509
j=
51.7246 -31.7678 21.3026
-32.9797 63.7409 -15.3834
-28,5386 17.3988 48.1036
dx =
-0,0016
-0.0007
-0.0018
v=
1.0500
0.9717
1.0400
d=
0
-0.0471
-0.0087
iter =
3
dc =

ASCET/EEE/B.TECH/20G25A0253 89
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 90
Power System & Simulation Lab

-0,0002
0.0014
-0.0001
j=
51.5967 -31.6941 21.1474
-32.9337 63.6841 -15.3520
-28,.5482 17.3966 47.9549
dx=
1.0e-004 *
0.1465
0.2778
-0,0428
Y=
1,0500
0.9717
1.0400
d=
0
-0.0471
0.0087
iter =
4
dc=
1.0e-004 *
-0.0000
-0.5314
-0.0001
j=
51.5964 -31.6937 21.1471
-32.9337 63.6846 -15.3516

ASCET/EEE/B.TECH/20G25A0253 91
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 92
Power System & Simulation Lab

-28.5482 17,3969 47.9545


dx =
1.0e-005 *
-0.0750
-0.1223
-0.0003
v=
1.0500
0.9717
1.0400
d=
0
-0.0471
-0.0087

RESULT:
Thus the load flow analysis using is performed by Newton-Raphson method and a
program is developed using MATLAB to find the solution of load flow for given power
system and the output is verified.

ASCET/EEE/B.TECH/20G25A0253 93
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 94
Power System & Simulation Lab

8. SHORT CIRCUIT FAULT ANALYSIS

AIM:
To develop a MATLAB program to perform the short circuit fault analysis for a given
power system.

SOFTWARE REQUIRED:
MATLAB software package

ALGORITHM:
STEP 1: Get the data for positive and negative sequence impedance materials.
STEP 2: Enter the bus code number, base kV, base MVA number,
STEP 3: Implementation of Z-Bus using BUILDING ALGORITHM.
STEP 4: Enter the choices.
STEP 5: Enter the impedance value to be included according to the corresponding
choices.
STEP 6: Enter the fault impendence value, fault bus code, and fault bus voltage.
STEP 7: Print the result.

Short circuit analysis:


A three phase symmetrical fault is caused by application of three equal fault
impedances 2; to the three phases. If Z; = 0 the fault is called a solid or a bolted fault. These
faults an be of two types: (a) line to line to line to ground fault (LLLG fault) or (b) line to line
to line fault LLL fault). Since the three phases are equally affected, the system remains
balanced. That is why, this fault is called a symmetrical or a balanced fault and the fault
analysis is done on per phase basis. The behaviour of LLLG fault and LLL fault is identical
due to the balanced nature of the fault. This is a very severe fault that can occur in a system
and if 2f = 0, this is usually the most severe fault that can occur in a system. Fortunately, such
faults occur infrequently and only about 5% of the system faults are three phase faults.
Symmetrical or Balanced three phase fault analysis:
In this type of faults all three phases are simultaneously short circuited. Since the
network remains balanced, it is analysed on per phase basis. The other two phases carry
identical currents but with a phase shift of 120°. A fault in the network is simulated by
connecting impedances in the network at the fault location. The faulted network is then
solved using Thevenin’s equivalent network as seen from the fault point. The bus impedance

ASCET/EEE/B.TECH/20G25A0253 95
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 96
Power System & Simulation Lab

matrix is convenient to use for fault studies as its diagonal elements are the Thevenin’s
impedance of the network as seen from different buses. Prior to the occurrence of fault, the
system is assumed to be in a balanced steady state and hence per phase network model is
used. The generators are represented by a constant voltage source behind a suitable reactance
which may be sub-transient, transient or normal d-axis reactance. The transmission lines are
represented by their n-models with all impedances referred to a common base. Further, a
balanced three phase fault, through fault impedance Zf is assumed to occur at kth bus as
shown in the figure. A pre-fault load flow provides the information about the pre-fault bus
voltage.

Let [ VBus(0)] be the prefault bus voltage vector [V2(0).... VK(0)....Vn(0)]T p.u. The fault at kth
bus through an impedance Zf will cause a change in the voltage of all the buses [∆VBUS] due
to the flow of heavy currents through the transmission lines. This change can be calculated by
applying a voltage Vk (0) at kth bus and short circuiting all other voltage sources. The sources
and loads are replaced by their equivalent impedances. This is shown in below figure

In above figure, Zj and Zk are the equivalent load impedances as bus i and k respectively, Zik
is the impedance of line between ith and kth buses. Xdl is the appropriate generator reactance,
Zf is the fault impedance,Ik(F) is the fault current and Vk(0) is the prefault voltage at kth bus,
From the superposition theorem, the bus voltages due to a fault can be obtained as the sum of
prefault bus voltages and the change in bus voltages due to fault, i.e.,

ASCET/EEE/B.TECH/20G25A0253 97
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 98
Power System & Simulation Lab

Where,
[ Vbus(F)]= Vector of bus voltages during fault = [V1(F)…..Vk(F)…...Vn(F)]T
[ Vbus(0)]= Vector of pre-fault bus voltages= [V1(0).... Vk (0)....Vn(0)]T
[∆Vbus) = Vector of change in bus voltages due to fault = [∆V1.... ∆Vk... ∆Vn]T
Also the bus injected current [IBUS] can be expressed as,

Where, [Vbus] is the bus voltage vector and [Ybus] is the bus admittance matrix. With all the
bus currents, except of the faulted bus k, equal to zero, the node equation for the network of
above figure can be written as

As the fault current Ik(F) is leaving the bus it is taken as a negative current entering the bus.
Hence,

[Vbus] can be calculated as:

where, [Zbus] is the bus impedance matrix = [Ybus]-1


Substituting the expression of [∆Vbus]

Expanding the above equation one can write,

ASCET/EEE/B.TECH/20G25A0253 99
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 100
Power System & Simulation Lab

The bus voltage of kth, bus can be expressed as:

Also from Fig

For a bolted fault Zf=0 and hence, Vk(F) = 0. Thus the fault current Ik(F) for bolted fault can
be expressed

For faulty with non-zero fault impedance Zf, the fault current can be calculated as:

The quantity Zkk in equations is the Thevenin’s impedance or open circuit impedance of the
network as seen from the faulted bus k. From equation [Vbus(F)] the bus voltage after fault for
the healthy buses can be written as:

Substituting Ik(F)

Where Zij is the impedance of line connecting buses i and j.

PROBLEM:
Consider the 3bus system as shown in figure. The generators are 100 MVA with
transient reactance 10% each. Both the transformers are 100 MVA with a leakage reactance
of 5%. The reactance of each of the lines to a base of 100 MVA, 110kV is 10%. Obtain the
short circuit calculation for a 3phase solid short circuit on bus 3. Assume pre-fault voltages to
be 1pu and pre-fault currents to be zero.

ASCET/EEE/B.TECH/20G25A0253 101
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 102
Power System & Simulation Lab

Write and execute a MATLAB program and also verify the output with the manual
calculation results.

PROGRAM:
clc;
display('........short circuit analysis.......');
display('...............................................’);
display(‘………19G21A0235……….’);
display('...............................................’);
n=input (‘enter number of busses in the power system’);
x=input ('enter 1.impedence given between buses 2.admittance given between buses:’);
y=zeros(n,n);
%reading power System line data
for(i=1:n)
for(j=1:n)
if (j>=i)
if (i==j)
fprintf('Enter the half line charging admittance between buses %d-%d',i,j)
y(i,j)=1i*input(‘:’);
elseif(x==1)

ASCET/EEE/B.TECH/20G25A0253 103
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 104
Power System & Simulation Lab

fprintf('Enter the impedance value between buses %d-%d',i,j)


z=input (‘:’);
y(i,j)=1/z;
else
fprintf('Enter the admittance between buses %d-%d',i,j)
y(i,j)=input(':');
end
end
else
y(I,j)=y(j,i);
end
end
end
%Ybus calculation
Y=zeros(n,n);
for (i=1:n)
for(j=1:n)
if (i==j)
%Diagonals of Ybus matrix
for (k=l:n)
Y(i,i)=Y(i,i)t=+y (i,k);
end
else,
%off diagonal elements of Ybus matrix
Y(i,j)=-y(j,i);
end
end
end
%Display Ybus
Y

ASCET/EEE/B.TECH/20G25A0253 105
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 106
Power System & Simulation Lab

Z=zeros(n,n);
Z=inv(Y)
k=input('fault occurred at bus’);
vpf=input('enter prefault voltage');
zff=input('fault impedance');
ifault=vpf/(Z(k,k)+zff)
for (i=l:n)
fprintf(‘Voltage at bus %d under fault’,i);
v-vpf=(g(t,k)*ifault)
fprintf(‘/n’);
end

OUTPUT:
………….short circuit analysis……………
……………………………………………..
………………19G21A0235………...…….
…………………………………………….
Enter number of buses In the power system:3
enter 1. impedance given between buses 2.admittance given between buses:1
Enter the half line charging admittance between buses 1-1:-1/0.051
Enter the Impedance value between buses 1-2:0.11
Enter the impedance value between buses 1-3:0.11
Enter the half line charging admittance between buses 2-2:-1/0.051
Enter the Impedance value between buses 2-3:0.11
Enter the half line charging admittance between buses 3-3:0
Y=
0-39.60781 0+10.00001 0 +10.0000i
0+10.00001 0-59.21571 0 +10.0000i
0+10.00001 0+10.00001 0-20.0000i

ASCET/EEE/B.TECH/20G25A0253 107
Power System & Simulation Lab

ASCET/EEE/B.TECH/20G25A0253 108
Power System & Simulation Lab

Z=
0+0.03281 0+0.0091i 0+0.0210i
0+0.0091i 0+0.0210i 0+0.0150i
0+0.0210i 0+0.0150i 0+0.0680i
Fault occurred at bus 3
Enter prefault voltage
fault impedance 0
ifault =
0 -14.7081i
voltage at bus 1 under fault
v=
0.6917
/n voltage at bus 2 under fault
v=
0.7791
/n voltage at bus 3 under fault
v=0

RESULT:
Thus, a program is developed using MATLAB to perform the short circuit fault
analysis for a given power system network and the output is verified.

ASCET/EEE/B.TECH/20G25A0253 109

You might also like