EASTERN UNIVERSITY

DEPARTMENT OF ELECTRICAL AND ELECTRONIC ENGINEERING

Lab Manual

Course Title: Power System Analysis Lab

Course Code: EEE 392
Eastern University
Department of Electrical & Electronic Engineering
Course Name: Power System Analysis Laboratory
Course Code: EEE 392
Experiment No: 01
Name of the Experiments: Familiarization with MATLAB software

Objective: To understand the basics of MATLAB software and some basic operations.

What is MATLAB?
MATLAB is the language of technical computing. Millions of engineers and scientists worldwide use it
to analyze and design the systems and products transforming our world. MATLAB is in automobile
active safety systems, interplanetary spacecraft, health monitoring devices, smart power grids, and LTE
cellular networks. It is used for machine learning, signal processing, image processing, computer vision,
communications, computational finance, control design, robotics, and much more.
Key Feature:
• High-level language for scientific and engineering computing
• Desktop environment tuned for iterative exploration, design, and problem-solving
• Graphics for visualizing data and tools for creating custom plots
• Apps for curve fitting, data classification, signal analysis, control system tuning, and many other tasks
• Add-on toolboxes for a wide range of engineering and scientific applications
• Tools for building applications with custom user interfaces
• Interfaces to C/C++, Java®, .NET, Python, SQL, Hadoop, and Microsoft® Excel®
• Royalty-free deployment options for sharing MATLAB programs with end users
Desktop Basics
When you start MATLAB, the desktop appears in its default layout.
The desktop includes these panels:
• Current Folder — Access your files.
• Command Window — Enter commands at the command line, indicated by the prompt (>>).
• Workspace — Explore data that you create or import from files.
As you work in MATLAB, you issue commands that create variables and call functions.
For example, create a variable named a by typing this statement at the command line:
a=1
MATLAB adds variable a, to the workspace and displays the result in the Command Window.
a=1
Create a few more variables.
b=2
b=2
c=a+b
c=3
d = cos(a)
d = 0.5403
When you do not specify an output variable, MATLAB uses the variable ans, short for answer, to store
the results of your calculation.
sin(a)
ans = 0.8415
If you end a statement with a semicolon, MATLAB performs the computation, but suppresses the display
of output in the Command Window.
e = a*b;
You can recall previous commands by pressing the up- and down-arrow keys, ↑ and ↓. Press the arrow
keys either at an empty command line or after you type the first few characters of a command. For
example, to recall the command b = 2, type b, and then press the up-arrow key.

Matrices and Arrays

MATLAB is an abbreviation for "matrix laboratory." While other programming languages mostly work
with numbers one at a time, MATLAB® is designed to operate primarily on whole matrices and arrays.
All MATLAB variables are multidimensional arrays, no matter what type of data. A matrix is a two-
dimensional array often used for linear algebra.
Array Creation
To create an array with four elements in a single row, separate the elements with either a comma (,) or a
space.
a = [1 2 3 4]
a=1234
This type of array is a row vector.
To create a matrix that has multiple rows, separate the rows with semicolons.
a = [1 2 3; 4 5 6; 7 8 10]
a=
123
456
7 8 10
Another way to create a matrix is to use a function, such as ones, zeros, or rand. For example, create a 5-
by-1 column vector of zeros.
z = zeros(5,1)
z=
0
0
0
0
0

Matrix and Array Operations

MATLAB allows you to process all of the values in a matrix using a single arithmetic operator or
function.
a + 10
ans =
11 12 13
14 15 16
17 18 20
sin(a)
ans =
0.8415 0.9093 0.1411
-0.7568 -0.9589 -0.2794
0.6570 0.9894 -0.5440
To transpose a matrix, use a single quote ('):
a'
ans =
147
258
3 6 10
You can perform standard matrix multiplication, which computes the inner products between rows and
columns, using the * operator. For example, confirm that a matrix times its inverse returns the identity
matrix:
p = a*inv(a)
p=
1.0000 0 -0.0000
0 1.0000 0
0 0 1.0000
Notice that p is not a matrix of integer values. MATLAB stores numbers as floating-point values, and
arithmetic operations are sensitive to small differences between the actual value and its floating-point
representation. You can display more decimal digits using the format command:
format long
p = a*inv(a)
p=
1.000000000000000 0 -0.000000000000000
0 1.000000000000000 0
0 0 0.999999999999998
Reset the display to the shorter format using
format short
format affects only the display of numbers, not the way MATLAB computes or saves them.
To perform element-wise multiplication rather than matrix multiplication, use the .* operator:
p = a.*a
p=
149
16 25 36
49 64 100
The matrix operators for multiplication, division, and power each have a corresponding array operator
that operates element-wise. For example, raise each element of a to the third power:
a.^3
ans =
1 8 27
64 125 216
343 512 1000
Complex Numbers
Complex numbers have both real and imaginary parts, where the imaginary unit is the square root of -1.
sqrt(-1)
ans =
0.0000 + 1.0000i
To represent the imaginary part of complex numbers, use either i or j.
c = [3+4i, 4+3j; -i, 10j]
c=
3.0000 + 4.0000i 4.0000 + 3.0000i
0.0000 - 1.0000i 0.0000 +10.0000i

Array Indexing
Every variable in MATLAB® is an array that can hold many numbers. When you want to access selected elements
of an array, use indexing.
For example, consider the 4-by-4 magic square A:
A = magic(4)
A=
16 2 3 13
5 11 10 8
9 7 6 12
4 14 15 1
There are two ways to refer to a particular element in an array. The most common way is to specify row and column
subscripts, such as
A(4,2)
ans =
14
If you try to refer to elements outside an array on the right side of an assignment statement, MATLAB throws an
error.
test = A(4,5)
Index exceeds matrix dimensions.
However, on the left side of an assignment statement, you can specify elements outside the current dimensions. The
size of the array increases to accommodate the newcomers.
A(4,5) = 17
A=
16 2 3 13 0
5 11 10 8 0
9 7 6 12 0
4 14 15 1 17
To refer to multiple elements of an array, use the colon operator, which allows you to specify a range of the form
start:end. For example, list the elements in the first three rows and the second column of A:
A(1:3,2)
ans =
2
11
7
The colon alone, without start or end values, specifies all of the elements in that dimension. For example, select all
the columns in the third row of A:
A(3,:)
ans =
9 7 6 12 0
The colon operator also allows you to create an equally spaced vector of values using the more general form
start:step:end.
B = 0:10:100
B=
0 10 20 30 40 50 60 70 80 90 100
If you omit the middle step, as in start:end, MATLAB uses the default step value of 1.

Text and Characters

To convert numeric values to characters, use functions, such as num2str or int2str.
f = 71;
c = (f-32)/1.8;
tempText = ['Temperature is ',num2str(c),'C']
tempText =
Temperature is 21.6667C

Calling Functions
MATLAB® provides a large number of functions that perform computational tasks. Functions are equivalent to
subroutines or methods in other programming languages. To call a function, such as max, enclose its input
arguments in parentheses:
A = [1 3 5];
max(A)
ans =
5
If there are multiple input arguments, separate them with commas:
B = [10 6 4];
max(A,B)
ans =
10 6 5
Return output from a function by assigning it to a variable:
maxA = max(A)
maxA =
5
When there are multiple output arguments, enclose them in square brackets:
[maxA,location] = max(A)
maxA =
5
location =
3
Enclose any character inputs in single quotes:
disp('hello world')
hello world
To call a function that does not require any inputs and does not return any outputs, type only the function name:
clc
The clc function clears the Command Window.
Line Plots
To create two-dimensional line plots, use the plot function. For example, plot the value of the sine function from 0 to
2π:
x = 0:pi/100:2*pi;
y = sin(x);
plot(x,y)

You can label the axes and add a title.

xlabel('x')
ylabel('sin(x)')
title('Plot of the Sine Function')
By adding a third input argument to the plot function, you can plot the same variables using a red dashed line.
plot(x,y,'r--')

To add plots to an existing figure, use hold.

x = 0:pi/100:2*pi;
y = sin(x);
plot(x,y)
hold on
y2 = cos(x);
plot(x,y2,':')
legend('sin','cos')

Until you use hold off or close the window, all plots appear in the current figure window.
Eastern University
Department of Electrical & Electronic Engineering
Course Name: Power System Analysis Laboratory
Course Code: EEE 392
Experiment No: 02
Name of the Experiments: Basic Power System Simulation with MATLAB software

Objective: The objective of this lab is to make the students familiar with some basics power system
simulation related operation using MATLAB software.

Sinusoidal Single phase Voltage, Current and Power

In this lab we consider only the steady state condition and the transient condition is not taken into
account. Suppose we have a 220V, 50Hz single phase AC supply which is fed to a 100Ω load. We would
like to plot the related signal. In general the ac voltage equation will be in the form, v = V msin(ωt), where
Vm=220V, ω=2πf is the angular frequency and t is time.
Sinusoidal voltage:
clear
clc
t=0:1/30000:1;
A=220;
f=60;
v=A*sin(2*pi*f*t);
R=100;
plot(t,v)
xlabel('time(sec)')
ylabel('voltage(V)')

250

200

150

100

50
voltage(V)

-50

-100

-150

-200

-250
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
time(sec)

Fig. 1 single phase sinusoidal voltage

Sinusoidal current:
The equation of current is
clear
clc
t=0:1/30000:1;
A=220;
f=60;
v=A*sin(2*pi*f*t);
R=100;
i=v/R;
plot(t,i)
xlabel('time(sec)')
ylabel('current(A)')

2.5

1.5

0.5
current(A)

-0.5

-1

-1.5

-2

-2.5
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
time(sec)

Fig. 2 single phase sinusoidal current for resistive load

Power:
The equation of power consumption by load is P=vi=i2R

clear
clc
t=0:1/30000:1;
A=220;
f=60;
v=A*sin(2*pi*f*t);
R=100;
i=v/R;
P=v.*i;
plot(t,P)
xlabel('time(sec)')
ylabel('Power(watt)')
500

450

400

350

300
Power(watt)

250

200

150

100

50

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
time(sec)

Fig. 3 single phase sinusoidal power

This figure is not clearly visible so, the figure below is given for power within 100ms

500

450

400

350

300
Power(watt)

250

200

150

100

50

0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
time(sec)

Fig. 4 single phase sinusoidal power

Voltage, Current and Power on the same plot:
As pure resistive load is used, the voltage and current of single phase circuit are in phase as seen from the figure
below:
clear
clc
t=0:1/30000:0.1;
A=220;
f=50;
va=A*sin(2*pi*f*t);
R=10;
ia=va/R;
P=va.*ia;
plot(t,va,t,ia)
xlabel('time(sec)')
grid on

250
voltage (V)
current (A)
200

150

100

50

-50

-100

-150

-200

-250
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
time(sec)

Fig. 5 single phase sinusoidal voltage and current for resistive load
If power is considered on the same plot then there is a tiny change in code in a line and the figure
becomes as follows:
plot(t,va,t,ia,t,P)

5000
voltage (V)
current (A)
Power (W )

4000

3000

2000

1000

-1000
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
time(sec)

Fig. 6 single phase sinusoidal voltage, current and power for resistive load
Now if we consider inductive load only of L=10mH, then the impedance becomes, Z=jωL=j2πfL.
So current equation will be i=v/Z or i=Imsin(ωt-90), where Im=Vm/Z.
clear
clc
t=0:1/30000:0.1;
A=220;
f=50;
va=A*sin(2*pi*f*t);
L=10*10^-3;
j=sqrt(-1);
Za=j*2*pi*f*L;
ia=(A/abs(Za))*sin(2*pi*f*t-angle(Za));
P=va.*ia;
plot(t,va,t,ia)
xlabel('time(sec)')
grid on

250
voltage(V)
current(A)
200

150

100

50

-50

-100

-150

-200

-250
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
time(sec)

Fig. 7 single phase sinusoidal voltage and current for inductive load
If power is considered on the same plot then it becomes a figure like below:

The change in the code will be only this

plot(t,va,t,ia,t,P)
8000
voltage(V)
current(A)
Power(W )
6000

4000

2000

-2000

-4000

-6000

-8000
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
time(sec)

Fig. 8 single phase sinusoidal voltage, current and power for inductive load
Up to this you will be able to plot the current, voltage and power for capacitive load of 6mF.

Sinusoidal Three Phase Voltage, Current & Power:

Suppose we have a three phase supply system which has three phase to neutral voltages va, vb and vc
who are 120 degree apart from each other and with a magnitude 220V. The supply system is 50Hz to a
three phase balanced Y connected load of 100Ω each.
Sinusoidal 3-Ҩ Voltage:
The three phase voltage equations are va = Vmsin(ωt), vb = Vmsin(ωt-120) and vc = Vmsin(ωt+120),
where the phase sequence is a-b-c.
clear
clc
t=0:1/30000:0.1;
A=220;
f=50;
va=A*sin(2*pi*f*t);
vb=A*sin(2*pi*f*t-120);
vc=A*sin(2*pi*f*t+120);
R=100;
plot(t,va,t,vb,t,vc)
xlabel('time(sec)')
ylabel('voltage(V)')
title('Voltage in three phase va, vb & vc')
Voltage in three phase va, vb & vc
250
va
vb
200 vc

150

100

50
voltage(V)

-50

-100

-150

-200

-250
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
time(sec)

Fig. 9 three phase sinusoidal voltages

Sinusoidal 3-Ҩ Current:
The equations of currents are , &

clear
clc
t=0:1/30000:0.1;
A=220;
f=50;
va=A*sin(2*pi*f*t);
vb=A*sin(2*pi*f*t-120);
vc=A*sin(2*pi*f*t+120);
R=100;
ia=va/R;
ib=vb/R;
ic=vc/R;
plot(t,ia,t,ib,t,ic)
xlabel('time(sec)')
ylabel('Current(A)')
title('Currents in three phase ia, ib & ic')
Currents in three phase ia, ib & ic
2.5
ia
ib
2 ic

1.5

0.5
Current(A)

-0.5

-1

-1.5

-2

-2.5
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
time(sec)

Fig. 10 three phase sinusoidal currents for balanced resistive load

Sinusoidal 3-Ҩ Power:
Practice it by yourself.
Load current with Resistive, Capacitive and Inductive Load:
Suppose we have a three phase balance Y connected load of 20Ω resistive, 30Ω capacitive reactance and
40Ω inductive resistance in series at each phase. So total impedance per phase= 20+j40-j30=20+j10. Now
if we plot the current of three phases we will find a figure below.
clear
clc
t=0:1/30000:0.1;
A=220;
f=50;
va=A*sin(2*pi*f*t);
vb=A*sin(2*pi*f*t-120);
vc=A*sin(2*pi*f*t+120);
R=20;
Xl=40;
Xc=30;
j=sqrt(-1);
Z=R+j*(Xl-Xc);
Zmag=abs(Z);
Zang=angle(Z);
ia=(A/Zmag)*sin(2*pi*f*t-Zang);
ib=(A/Zmag)*sin(2*pi*f*t-120-Zang);
ic=(A/Zmag)*sin(2*pi*f*t+120-Zang);
plot(t,va,t,vb,t,vc, 'linewidth',2)
hold on
plot(t,ia,t,ib,t,ic)
xlabel('time(sec)')
grid on
250
va
vb
200 vc
ia
ib
150 ic

100

50

-50

-100

-150

-200

-250
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
time(sec)

Fig. 11 three phase sinusoidal voltages & currents for balanced resistive, inductive & capacitive load in series
Now what will happen if the load is not balanced? Say we have unbalanced Y connected load of Za=15-
j12, Zb=20+j5 and Zc=18-j24. Now the following figure shows the effect of unbalanced load in current.
clear
clc
t=0:1/30000:0.1;
A=220;
f=50;
va=A*sin(2*pi*f*t);
vb=A*sin(2*pi*f*t-120);
vc=A*sin(2*pi*f*t+120);
j=sqrt(-1);
Za=15-j*12;
Zb=20j*5;
Zc=18-j*24;
Zmag_a=abs(Za);
Zang_a=angle(Za);
Zmag_b=abs(Zb);
Zang_b=angle(Zb);
Zmag_c=abs(Zc);
Zang_c=angle(Zc);
ia=(A/Zmag_a)*sin(2*pi*f*t-Zang_a);
ib=(A/Zmag_b)*sin(2*pi*f*t-120-Zang_b);
ic=(A/Zmag_c)*sin(2*pi*f*t+120-Zang_c);
plot(t,va,t,vb,t,vc)
hold on
plot(t,ia,t,ib,t,ic, 'linewidth',2)
xlabel('time(sec)')
grid on
250
va
vb
200 vc
ia
ib
150 ic

100

50

-50

-100

-150

-200

-250
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
time(sec)

Fig. 12 three phase sinusoidal currents for unbalanced resistive, inductive & capacitive load in series

Report:
1. Plot the single phase current for a supply of 220V, 60Hz system to a RLC series load of R=100Ω,
L=10mH and C=0.1µF and also show your code in the report.
2. Plot the Three line currents and three phase currents for a ∆ connected load of 100Ω each and
also show your code in the report. The phase sequence is a-c-b.
3. For an unbalanced ∆ connected load of R=100Ω, L=15mH & C=6µF at phase a, R=80Ω,
L=10mH & C=1µF at phase b and R=110Ω, L=5mH & C=3µF at phase c, plot the three phase
currents and also show your code in the report. The phase sequence is a-b-c.
Eastern University
Department of Electrical & Electronic Engineering
Course Name: Power System Analysis Laboratory
Course Code: EEE 392
Experiment No: 03
Name of the Experiments: Familiarization with MATLAB Simulink

Objective: To understand the basics of MATLAB Simulink software and simple design using the tools.

Create Simple Model

You can use Simulink® to model a system and then simulate the dynamic behavior of that
system. The basic techniques you use to create the simple model in this lab are the same
techniques that you use for more complex models. To create this simple model, you need four
Simulink blocks. Blocks are the model elements that define the mathematics of a system and
provide input signals:

 Sine Wave — Generate an input signal for the model.

 Integrator — Process the input signal.
 Bus Creator — Combine multiple signals into one signal.
 Scope — Visualize and compare the input signal with the output signal.

Simulating this model integrates a sine wave signal to a cosine signal and then displays the
result, along with the original signal, in a Scope window.
Open New Model in Simulink Editor
Use the Simulink Editor to build your models.
1. Start MATLAB®. From the MATLAB Toolstrip, click the Simulink button . A short
delay occurs the first time you open Simulink.
2. Click the Blank Model template, and then click the Create Model button. The Simulink
Editor opens with a new block diagram.

3. Select File > Save as. In the File name text box, enter a name for your model, For
example, simple_model.slx. Click Save.

Open Simulink Library Browser

From the Simulink Library Browser you can search for blocks to use in your model. Also, you
can create a new Simulink model, project, or Stateflow® chart.
1. From the Simulink Editor toolbar, click the Simulink Library button .
2. Set the Library Browser to stay on top of the other desktop windows. On the Library
Browser toolbar, select the Stay on top button .
Browse or Search for Specific Blocks
To browse through the block libraries, select a MathWorks® product and then a functional area in
the left pane. To search all of the available block libraries, enter a search term.
1. Search for a Sine Wave block. In the search box on the browser toolbar, enter sine, and
then press the Enter key. Simulink searches the libraries for blocks with sine in their name
or description, and then displays the blocks.
2. Get detailed information about a block. Right-click a block and then select Help for the
<block name>. The Help browser opens with the reference page for the block.
3. View block parameters. Right-click a block and then select Block Parameters. The
block parameters dialog box opens.
Add Blocks to Model
You build models by dragging blocks from the Simulink Library Browser window to the
Simulink Editor or single-clicking your model and entering a search term.
To build the simple model, begin by copying blocks from the Simulink Library Browser to the
Simulink Editor.
1. In the left pane of Simulink Library Browser, select the Sources library.
2. From the right pane, select the Sine Wave block.

3. Drag the Sine Wave block to the Simulink Editor. A copy of the Sine Wave block
appears in your model.

4. Add a Scope block using this alternative procedure:

a. Click within the block diagram.
b. After the search icon appears, type scope, and then from the list, select Scope.
5. Add the following blocks to your model using one of the approaches you used to add the
Sine Wave and Scope blocks.
Library Block
Continuous Integrator
Signal Routing Bus Creator

6. Your model should now have the blocks you need for the simple model.

Move and Resize Blocks

Before you connect the blocks in your model, arrange the block inputs and outputs to make
signal connections as straightforward as possible.
1. Move the Scope block after the Bus Creator block output. You can either:
 Click and drag the block.
 Select the block, and then press the arrow keys on your keyboard.
2. Move the blocks until your model looks similar to the following figure.
Block Connections
Most blocks have angle brackets on one or both sides. These angle brackets represent input and
output ports:
 The > symbol pointing into a block is an input port.
 The > symbol pointing out of a block is an output port.

You connect block output ports to input ports with lines. The lines represent signals with time-
varying values.
Draw Signal Lines between Blocks
Connect the blocks by drawing lines between output ports and input ports.
1. Position the cursor over the output port on the right side of the Sine Wave block.
The pointer changes to a cross hair (+).
2. Click, and then drag a line from the output port to the top input port of the Bus Creator
block.
While holding down the mouse button, the connecting line appears as a red dotted arrow.
3. Release the mouse button when the pointer is over the output port.
Simulink connects the blocks with a line and an arrow indicating the direction of signal
flow.

4. Connect the output port of the Integrator block to the bottom input port on the Bus
Creator block using a Ctrl key shortcut:
a. Select the Integrator block.
b. Press and hold the Ctrl key.
c. Click the Bus Creator block.
The Integrator block connects to the Bus Creator block with a signal line.
5. Connect the Bus Creator block to the Scope block by aligning ports:
a. Click and drag the Scope block until its input port is aligned with the Bus Creator
output port. A light blue line appears between the ports when the blocks line up
horizontally.

b. Release the mouse button. Ablue arrow appears as a suggested connection.

c. Click the end of the blue arrow. The arrow changes to a black signal line
connecting the blocks.
Draw Branched Signal Lines
Your simple model is almost complete. To finish the model, connect the Sine Wave block to the
Integrator block. This connection is different from the other connections, which all connect
output ports to input ports.
1. Hold down the Ctrl key.
2. Position the cursor where you want to start a branch line. Click, and then drag the cursor
away from the line to form a dotted-red line segment.

3. Drag the cursor to the Integrator input port, and then release the mouse button. The new
line, called a branch line, carries the same signal that passes from the Sine Wave block to
the Bus Creator block.
4. Drag line segments to straighten and align with blocks. Your model is now complete.
Define Configuration Parameters
Before you simulate the behavior of your model, you can modify the default values for the
configuration parameters. The configuration parameters include the type of numerical solver,
start time, stop time, and maximum step size.
1. From the Simulink Editor menu, select Simulation > Model Configuration
Parameters. The Configuration Parameters dialog box opens to the Solver pane.
Tip Alternatively, you can open the Model Configuration Parameters dialog box

by clicking the parameters button on the Simulink Editor toolbar.

2. In the Stop time field, enter 20. In the Max step size field, enter 0.2. With the Solver
parameter set to auto, Simulink determines the best numerical solver for simulating your
model.

3. Click OK.
Run Simulation
After you define the configuration parameters, you are ready to simulate your model.
1. From the Simulink Editor menu bar, select Simulation > Run.
The simulation runs, and then stops when it reaches the stop time specified in the Model
Configuration Parameters dialog box.
Tip Alternatively, you can control a simulation by clicking the Run simulation

button and Pause simulation button on the Simulink Editor toolbar or

Scope window toolbar.
Observe Simulation Results
After simulating a model, you can view the simulation results in a Scope window.
1. Double-click the Scope block.
The Scope window opens and displays the simulation results. The plot shows a sine wave
signal with the resulting cosine wave signal.
2. On the Scope window toolbar, click the Style button.

A Style dialog box opens with display options.

3. Change the appearance of the display. For example, select white for the display color and
axes background color (icons with a pitcher).
4. Select black for the ticks, labels, and grid colors (icon with a paintbrush).
5. Change signal line colors for the Sine Wave to blue and the Integrator to red. To see your
changes, click OK or Apply.
Report:
1. Design a simple block diagram of sinusoidal signal using the differential operator.
Eastern University
Department of Electrical & Electronic Engineering
Course Name: Power System Analysis Laboratory
Course Code: EEE 392
Experiment No: 04
Name of the Experiments: Demonstration on Bus Admittance and Bus Impedance Matrices

Objective: To develop a computer program to form the bus admittance matrix, Ybus of a power
system.

Introduction:
Consider a Power system with n-nodes or bus bars. The voltage-current relationship of the
system shown in Figure 4.2 is given by the following equations

1 3
z12

z14 z34 z23

z24
2
4

Fig. 4.1. A four node system

3 y23 2

y13 y34
y24
1
y14 4

y11
y33 y44 y22

Ground

Here, y11, y22 ,y33, and y44 are self load of buses.
I1 = y11E1 + y13(E1 – E3) + y14(E1 – E4) .....................(4.1)
I2 = y22E2 + y23(E2 – E3) + y24(E2 – E4)
0 = y33E3 + y31(E3 – E1) + y32(E3 – E2) + y34(E3 – E4)
0 = y44E4 + y41(E4 – E1) + y42(E4 – E2) + y43(E4 – E3)
After rearranging,
I1 = (y11 + y13 + y14 )E1 – y13E3 – y14E4 .....................(4.2)
I2 = (y22 + y23 + y24 )E2 – y23E3 – y24E4
0 = -y31E1 – y32E2 + (y31 + y32 + y33 + y34)E3 – y34E4
0 = -y41E1 – y42E2 – y43E3 + (y41 + y42 + y43 + y44)E4

In general form
I1 = Y11E1 + Y12E2 + Y13E3 + Y14 E4 .....................(4.2)
I2 = Y21E1 + Y22E2 + Y23E3 + Y24 E4
I3 = Y31E1 + Y32E2 + Y33E3 + Y34 E4
I4 = Y41E1 + Y42E2 + Y43E3 + Y44 E4
Here,
y12 = y21 = 0. Because, there is no transmission line between bus 1 and 2. Therefore,
transmission line impedance = , which means that the admittance y12 = 1/z12 = 0.
Also diagonal elements
Y11 = y11 + y13 + y14 = y11 + y12 + y13 + y14 [ y21 = y12 = 0 ]
Y22 = y22 + y23 + y24 = y21 + y22 + y23 + y24 .....................(4.3)
Y33 = y31 + y32 + y33 + y34
Y44 = y11 + y12 + y13 + y14
Off diagonal elements of the matrix can be can be calculated as
In matrix form,
 I 1  Y11 Y12 Y13 Y14   E1 
 I  Y Y22 Y23 Y24   E 2 
 2    21
 I 3  Y31 Y32 Y33 Y34   E 3 
    
 I 4  Y41 Y42 Y43 Y44   E 4 

In Fig, 4.1., it is clear that there are only two injected currents supplied by the two generators.
So, I3 and I4 are equal to zero. Also Y12 = Y21 = -y21 = 0.

 I 1  Y11 0 Y13 Y14   E1 

I   0 Y Y23 Y24   E 2 
 2   22

 0  Y31 Y32 Y33 Y34   E 3 

    
 0  Y41 Y42 Y43 Y44   E 4 
Problem 4.1: Determine the YBUS of the following power system. All impedance values are in
per unit.

1 3
j0.2

j0.5 j0.125 j0.25

j0.1 2
4

Given Self load for different buses are as follows:

z11 = j0.25; z22 = j0.5; z33 = j0.4 and z44 = j0.125.

Solution: At first we have to calculate all admittances.

y11 = 1/z11 = 1/j0.25 = -j4. Similarly,

y22 = - j2;
y33 = - j2.5
and y44 = - j8.
Also we can calculate transmission line impedances.

y13 = - j5
y14 = - j2
y23 = -j4
y24 = - j10
y34 = -- j8

Y13 = - y13 = j5
Y14 = - y14 = j2
Y23 = - y23 = j4
Y24 = - y24 = j10
Y13 = - Y34 = -- j8

Y11 = y11 + y13 + y14 = -j(4+5+2) = - j11

Admittance Matrix Formation (Y bus formation):

MATLAB Code for Y bus formation of given power system network

clear all;
clc;
display('----------Ybus formation code--------');
nbranch=input('enter the number of branches in system = ');
display('enter line data');
for n=1:1:nbranch
fb=input('enter from bus = ');
tb=input('enter to bus = ');
r=input('enter value of resistance = ');
x=input('enter value of reactance = ');
B=input('enter the value of line charging admittance(b/2) =
');
z=r+i*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=i*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 linked =
');
x0=input('enter value of ground reactor = ');
if x0~=0
y0=1./(i*x0);
else
y0=0;
end
Yreact(k,:)=[bs y0];
end

bs=Yreact(:,1);
y0=Yreact(:,2);
% off diagonal element

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

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

Yb=Y

Result: The Y bus matrix was formed for the given system by direct inspection method and the
results were verified using MATLAB program.

Example:

Q.4(a) Calculate the YBUS of the following power system

1 3 Self Loads:
0.3+j1.2
Gen1: 0.1+ j3.4
j5 0.2+j1.125 j4.25

0.3+j2.1 2 All impedances are

4
given in per unit.
Eastern University
Department of Electrical & Electronic Engineering
Course Name: Power System Analysis Laboratory
Course Code: EEE 392
Experiment No: 05
Name of the Experiments: Demonstration on Complex Power Flows [Load Flow Study]

5.1. Objective:
1. The load flow solution gives the nodal voltages and phase angles and hence the power
injection at all the buses and power flows through interconnecting power channels.
2. System transmission loss minimizes.

5.2. Introduction:
In order to determine the characteristics of power system under steady state operating mode
complex power flow(load flow) techniques are normally used. This is known as load flow study.

In load flow operations, active power generations are normally specified and generator bus
voltage magnitudes are normally maintained at a specified level (by AVR and machine
excitation). The loads are usually specified as constant power type with real and reactive
components being given. The loads are assumed to be unaffected by small variations of voltage
and frequency expected during normal steady state operation.

Why are load flow studies essential?

1. In planning the future development of the system as it brings out the effect of
a. interconnection with other power system
b. new loads
c. new generating stations
d. new transmission lines
Even before they are installed
2. In determining the best operation of the existing systems as it provides information on
a. line flows (active and reactive power)
b. line losses
c. voltage throughout the system
3. To provide data for stability studies on the system.
What is load flow study?
By load flow study we mean the determination of the following four quantities at each
and every node of any power system.
1) The voltage magnitude, V
2) The phase angle of the voltage, 
3) The real power injection, P
4) The reactive power injection, Q

The type of buses and corresponding parameters being specified are as under:
1. Load bus(or PQ bus): The total injected power (Pi + jQi ) is specified and normally the
loads are constant power type and remain unaffected by small variations in bus voltage.
2. Slack bus(or swing or reference bus): This type of bus arises because the system losses
are not known in advance before the load flow calculation. Any generator bus may be
chosen as slack bus, which is assumed to supply the line losses. The slack bus voltage is
assigned as system reference angle and its complex bus voltage (V + j0) is specified.
3. The voltage controlled bus ( or PV bus): The total injected power (active) P i is specified
while the voltage magnitude is maintained at |V| (using reactive power injection). This
type of bus usually corresponds to a generator bus or even a load bus where the bus
voltage is fixed by supplying reactive power from reactive power compensators.
The load flow solution is expected to give the final result about voltage magnitudes with angles,
active and reactive power flows in the individual line, losses and reactive power absorbed or
generated at PV bus.
A complex power flow(or load flow) solution have the following solution properties:
 High computational speed
 Simplicity of the program
 Flexibility of the program
 Low computer storage
 Reliability of the solution

5.3. ANAYTICAL FORMULATION OF COMPLEX POWER FLOW SOLUTION

In any power flow problem, it is required to have four variables at each bus i of the system, e.g.
Pi (the network active bus power where i stands for bus no.), Qi (the network reactive bus power),
|Vi| (the bus voltage magnitude), δi(the voltage phase angle). Only two of these variables are
known a priori and the load flow solution provides the solution of the remaining two variables at
any bus.

The complex power injected at any bus is given by

S i  Vi I i*
Here, Vi  Vi  i
n
I i   Yik V k  Yi1V1  Yi 2V2  Yi 3V3  .................  YinVn
k 1
....................................................(2.1)
n
Ii = Net current injected into the network at bus i = Yi1V1 + Yi2V2 + …… + YinVn = Y
k 1
ik Vk

Yik = Element of YBUS [ith row and kth column]

In complex form Si can be written as

Si = Pi + jQi
*
n  n n n
= Vi  Yik Vk   Vi  Yik*Vk*   ViVk*Yik*   Vi Vk Yik  i   k   ik
 k 1  k 1 k 1 k 1

................(2.2)

Separating into real and imaginary parts, we obtain

n n
Pi   ViVk Yik cos( i   k   ik )  VV Y i k ik cos( ik   k   i )
k 1 k 1
n n
.................(2.3)
Qi   ViVk Yik sin( i   k   ik )    ViVk Yik sin( ik   k   i )
k 1 k 1

The above two equations represents the polar form of the power flow also and is known as
steady state load flow equations (SLFE). It provides the calculated values of the network real and
reactive power injected at any bus i.

Let us consider that there are n no. of buses. Out of them m number of generators are there.
Therefore no. of load buses will be n-m.
Summary of power flow variables
Type of buses No. of buses Quantity Specified
 Slack (or swing or reference bus) one |Vi| and δi
 Voltage controlled (or PV) m-1 |Vi| and Pi
Buses 2,3,4,…………m
 Load bus(or PQ bus) n-m Pi and Qi
Buses m+1, m+2, ……….n

5.4 Gauss-Seidal method for the solution of non-linear equations

2x1 + x1x2 =1 .....................................................................................(5.4)

2x2 – x1x2 = -1 .....................................................................................(5.5)

Though there are two equations and two unknowns but they can not be solved directly because
they are non-linear equations. Iterative methods are to be used to solve them.
x1 = - 0.5 – 0.5x1x2
x2 = 0.5 –0.5 x1x2

Let us consider initial values x1 and x2 are [ 0 , 0].

Iteration 1
Putting initial values in (5.5), we have,
x11= 0.5, and putting x1 and x2 are [ 0.5 , 0] in (5.6), we have
x21= - 0.5.
Iteration 2
Putting x1 and x2 are [ 0.5 , -0.5] in (2.5), we have
x12= 0.625, and putting x1 and x2 are [ 0.625 , -0.5] in (5.6), we have
x22= - 0.6563.

Iteration 3
Putting x1 and x2 are [ 0.625 , -0.6563] in (5.5), we have
x13= 0.7051, and putting x1 and x2 are [ 0.7051 , -0.6563] in (5.6), we have
x23= - 0.7314.

Iteration 4
Putting x1 and x2 are [0.7051 , - 0.7314] in (5.5), we have
x14= 0.7578, and putting x1 and x2 are [0.7578, - 0.7314] in (5.6), we have
x24= - 0.7771. This is the result after 4th iteration.
.
.
.
After few iterations, the result will be [0.99 , - 0.99]

The final answer is [ 1, -1].

5.5 Gauss-Seidal method for load flow solution

We know ,
Si = Pi + jQi = Vi Ii*
n
Si* = Pi – jQi = Vi*Ii = Vi*  Yik Vk  Vi* (Yi1V1  Yi 2V2  Yi 3V3  .................  YinVn )
k 1
.................(5.7)
This equation may be expressed in the folllowing form also

n
Si* = Pi – jQi = Vi*YiiVi  Vi*  Yik Vk for i = 2, 3,4,...................,n .......................(5.8)
k 1
k i
n
Vi*YiiVi = Si* - Vi*  Yik Vk , for i = 2, 3,4,...................,n
k 1
k i

 
1  S i* n
Vi =   Yik Vk  , for i = 2, 3,4,...................,n .....................(5.9)
Yii  Vi* k 1 
 k i 

In Gauss-Seidal algorithm, the above equation(2.9) is utilized to find the final bus voltages using
successive steps of iterations, where
 
1  S *i n
Vi new
=   Yik Vk  , for i = 2, 3,4,...................,n ..........(5.10)


Yii  Viold 
*
k 1
k i



5.6 Gauss-Seidal method for load flow solution – An Example

Example: For the following 3-bus system

i) Calculate all bus voltages
ii) Total Line loss
iii) Swing Bus active and reactive power
iv) PV bus reactive power

V1 = 1.030 pu
1 2
z12 =0.01+j0.2
V2 = |1.04| pu

P2 = 80MW

SL1 =50+j30 MVA z13 =0.02+j0.3 z23 =0.03+j0.4

SL2 =40+j30 MVA

3

Step 1: Calculation of Ybus from line data

Primitive admittances: y12 = 1/z12 = 0.2494 -j 4.9875 pu = y21

y13 = 1/z13 = 0.2212 - j3.3186 pu = y31
y23 = 1/z23 = 0.1865 - j2.4860 pu= y32
Also given, B/2= j0.2 [ For Transmission line, Capacitive admittance, Y = G +jB]
Y11 = y12 + y13 + j0.02 + j0.02 = 0.4706 – j8.2661 pu
Y22 = y21 + y23 + j0.02 + j0.02 = 0.4358 - j 7.4335 pu
Y33 = y31 + y32 + j0.02 + j0.02 = 0.4077 - j5.7646 pu
Y12 = Y21 = -y12 = - 0.2494 + j 4.9875 pu
Y13 = Y31 = -y13 = - 0.2212 + j3.3186 pu
Y23 = Y32 = -y23 = - 0.1865 + j2.4860 pu
Y11 Y12 Y13   0.4706 - j8.2661 - 0.2494 + j4.9875 - 0.2212 + j3.3186 
YBus = Y21 Y22 Y23  = - 0.2494 + j4.9875 0.4358 - j7.4335 - 0.1865 + j2.4860 
Y31 Y32 Y33  - 0.2212 + j3.3186 - 0.1865 + j2.4860 0.4077 - j5.7646 

Step 2: Assumption of initial values of voltages and corresponding angles

Consider initial values are: V1 = 1.03 0 Slack bus. Magnitude of V and angle are given.
V2 = 1.04 0 PV bus. Magnitude of V given but angle is assumed to 0 .
V3 = 1.000 PQ or load bus. Both the magnitude of V and angle are unknown.

Step 2: Caculation of current and hence Q for PV bus only

To calculate new value of voltage using equation (2.10), it is necessary to know the value of S =
P + j Q. But in the case of PV bus, though P is given but Q is unknown. To calculate Q for PV
bus, we have to calculate S = VI* , where I = YV.
Here, Bus 2 is PV bus.
I2 = Y21V1 + Y22V2 + Y23V3 = 0.0100 -j 0.1077 pu
Q2 = imag(V2I2*) = 0.1120.

Step 3: Caculation of complex power S and hence V

S2 = P2 + jQ2 = (PG2 – PL2) + Q2 = 0.4 + 0.112

V2 = 1.0429 + j 0.0502, here |V2| = 1.0442 pu. This magnitude change is not accepted because
it is a PV bus. Its magnitude is fixed and equal to 1.04 pu. However, angle change is acceptable.
It is 2.7573 degree.
 V2 = 1.042.7573 pu.

Bus 3 is load or PQ bus. Complex power S is given.

S3 = P3 + jQ3 = (PG3 – PL3) + j (QG3 – QL3) = - 0.4000 - j0.3000
V3 = 0.9841 - j 0.0443, here |V3| = 0.9851 pu. and angle equal to – 2.5766 degree. Both
magnitude and angle change are accepted because it is a PQ bus where all were unknown and
have to be determined.

So, at the end of 1st iteration the result is

V1 = 1.030 pu

V2 = 1.042.7573 pu

For 2nd iteration

I2 = Y21V1 + Y22V2 + Y23V3 = 0.4944 - j 0.1081pu

Q2 = imag(V2I2*) = 0.1370.
S2 = P2 + jQ2 = (PG2 – PL2) + Q2 = 0.4 + 0.137
V2 = 1.0387 + j 0.0361, here |V2| = 1.0393 pu. This magnitude change is not accepted because
it is a PV bus. Its magnitude is fixed and equal to 1.04 pu. However, angle change is acceptable.
It is 1.9887 degree.
 V2 = 1.041.9887 pu.

Bus 3 is load or PQ bus. S data is given.

S3 = P3 + jQ3 = (PG3 – PL3) + j (QG3 – QL3) = - 0.4000 - j0.3000
V3 = 0.9806 - j 0.0486, here |V3| = 0.9818 pu. and angle equal to – 2.839 degree.

So, after 2nd iteration the result is

V1 = 1.030 pu

V2 = 1.041.9887 pu

Step 3: Caculation of Line flows i.e. active and reactive power flows (All are in pu)
I12=(V(1)-V(2))*y(1,2) + V(1)*j yself(1,2) = -0.1720 + 0.0386i
I21=(V(2)-V(1))*y(1,2) + V(2)*j yself(2,1) = 0.1720 - 0.0386i
I13=(V(1)-V(3))*y(1,3) + V(1)*j yself(1,3) = 0.1756 - 0.1784i
I31=(V(3)-V(1))*y(1,3) + V(3)*j yself(3,1) = -0.1756 + 0.1784i
I32=(V(3)-V(2))*y(3,2) + V(3)*j yself(3,2) = -0.2190 + 0.1498i
I23=(V(2)-V(3))*y(3,2) + V(2)*j yself(2,3) = 0.2190 - 0.1498i

S12= V(1)*conj(I12) = -0.1770 - 0.0610i

S21= V(2)*conj(I21) = 0.1774 + 0.0244i
S13= V(1)*conj(I13) = 0.1804 + 0.1381i
S31= V(3)*conj(I31) = -0.1793 - 0.1622i
S23= V(2)*conj(I23) = 0.2226 + 0.1230i
S32= V(3)*conj(I32) = -0.2207 - 0.1378i

Loss = 0.0034 - 0.0756i pu

The complex power generation from the buses are

S1 = 0.5034 + 0.3771i
S2 = 0.8000 + 0.4473i
S3 = -0.0000 + 0.0000i = 0 That means complex power injection = complex power
received. or S31 + S32 = PL3 + jQL3

5.7 Gauss-Seidal method for load flow solution – Computer program – Line capacitance is
given – pi model for transmission line is used.

clc
clear
%From to R X B/2
A=[1 2 .01 .2 .02;
1 3 .02 .3 .02;
2 3 .03 .4 .02];
P1=max(A(:,1));
P2=max(A(:,2));
n=max(P1,P2);
d1=size(A);
rt=d1(1); % Total no of rows
y=zeros(n,n); % Initialize primitive admittance matrix to zero.
% y matrix caculation(primitive admittance)
for k=1:rt
D=A(k,:);
r=D(1); c=D(2); R=D(3); X=D(4); B=D(5);
z(r,c)=R+j*X;
y(r,c)=1/z(r,c);
y(c,r)=y(r,c);
y(r,r)=y(r,r)+j*B;
y(c,c)=y(c,c)+j*B;
yself(r,c)=B;
yself(c,r)=B;
end
for k=1:n
for m=1:n
if m==k % It is diagonal element of YBUS
P=y(k,:);
Y(k,k)=sum(P);
else
Y(k,m)=-y(k,m);
end
end
end

% All data in pu 100 MVA base

%Busno Pg Qg PL QL |V| Buscode 1-swing, 2-PV, 3-load
PQ=[1 0 0 .5 .3 1.03 1;
2 .8 0 .4 .3 1.04 2;
3 0 0 .4 .3 1 3];
% Let the base MVA = 100
PG=PQ(:,2); % MW
QG=PQ(:,3); % MVar
PL=PQ(:,4); % MW
QL=PQ(:,5); % MVar
V0=PQ(:,6); % pu
BC=PQ(:,7); % To identify swing(1),PV(2) and Load buses(3).

P=(PG-PL); % Line Flows pu

Q=(QG-QL); % pu
S=P+j*Q;
V=V0;
V00=V0; % for Convergence checking
Vtol=.000001; Dtol=.000001; % Tolerance checking Vtol(pu), Dtol(rad)
for it=1:100
for k=2:n
if BC(k)==2 %PV Bus code is 2
sumI=0;
for m=1:n

sumI=sumI+Y(k,m)*V(m);
end
I(k)=sumI;

Q(k)=imag(V(k)*conj(I(k)));
end
S(k)=P(k)+j*Q(k);

sumV=0;

for m=1:n
if m ~=k
sumV=sumV+Y(k,m)*V(m);
end
end

V(k)=(conj(S(k)/V(k))-sumV)/Y(k,k);

D=angle(V(k))

if BC(k)==2 % PV bus
V(k)=V0(k)*(cos(D)+j*sin(D));
end
end

% Checking convergence Vtol=.0001 pu Dtol=.0001 rad

DV=abs(abs(V)-abs(V00));
DVV=floor(DV/Vtol);
SDVV=sum(DVV);
DD=abs(abs(angle(V))-abs(angle(V00)));
DDD=floor(DD/Dtol);
SDDD=sum(DDD);

if SDVV==0 & SDDD ==0

break,end;
V00=V;
end

% Line Flow calculation

Loss=0;
for a=1:n
for b=1:n
if a~=b

I(a,b)=(V(a)-V(b))*y(a,b)+V(a)*j*yself(a,b);
SS(a,b)=V(a)*conj(I(a,b));
Loss=Loss+SS(a,b);
end
end
end

%For Generation Checking

for k=1:n
SSS(k)=sum(SS(k,:))+PL(k)+j*QL(k); % Total Generation
end