The University of Derby

Load Flow studies using Digital computing and Matlab

Power Flow Analysis Tools

Electrical Power Applications (6EJ503)

TUTOR: Dr Mahmud Shafik

STUDENT: Salau Oluwasesan Abiodun

Student number: 100277737

02/05/2013

1
Table of Contents
1.0 Executive Summary ………………………………………………………………………………… 3
2.0 Aims and Objectives ……………………………………………………………………………….. 4
3.0 Introduction …………………………………………………………………………………………… 4

4.0 Method No 1: Jacobi…………………………………………………………………………………….. 8

4.1 Calculation: .............................................................................................................. 8

.................................................................................................................... 9

.................................................................................................................. 10

.................................................................................................................. 11

................................................................................................................. 13

4.4 MATLAB Code: JACOBI (Without shunt capacitance) ............................................... 13

5.0 METHOD NO 2: Z- Matrix Approach using a block substitution method ............... 22

.............................................................................................................................. 25

5.1.1 ............................................................................................................... 26

5.1.2 ................................................................................................................ 27

6.0 CONCLUSIONS ............................................................................................................ 35

7.0 REFERENCES ............................................................................................................... 36

2
1.0 EXECUTIVE SUMMARY

Load flow study is a vital tool which involves numerical analysis of electrical Power systems in normal
steady-state operation by taking into account the voltage at every bus, and voltage drop in each
feeder and also the power flow in all branch and feeder circuits so as to determine whether the
system specifications are within stated limits under changing load conditions and to ensure the
integrity, safety and timely operations of electrical networks equipments such as the transformers
and conductors etc.

Load Flow computations provide active and reactive power flows and magnitude of the bus voltage
including its phase angle at all buses for a stated power system network and operating condition
subject to the generation and/or regulating abilities of generators, synchronous condensers, static
reactive power compensators, High Voltage D.C ( HVDC) controls, flexible AC Transmission Systems
(FACTS) controllers, tap changing under load transformers and stated net interchange amidst
individual operating systems (utilities).This knowledge is important for the constant evaluation of
the present performance of a power system and for the analyses of the efficiency of back-up plans
for system expansion to meet increasing load demands. These examinations require the estimation
of various power flow operating conditions ranging from normal and emergency conditions. The
output from load flow analysis usually provide the primary conditions required for other analysis like
transient stability, dynamic stability studies, short circuit fault analysis, economic load dispatch
analysis, Harmonic distortion profile, transformer compensation requirements for high tension
services and unit commitment etc.

Load flow studies are centred on a nodal voltage analysis of a power system. Load flow problems are
solved through some iterative methods such as: Gauss –Seidal method, Newton-Raphson, Jacobi
method and Z – matrix method (Block substitution and Forward substitution methods).

Nonetheless, the iterative methods mentioned above are becoming obsolete and are being replaced
by some stochastic search and simulated annealing approaches such as the Genetic algorithm,
Genetic tabulated search, Differential evolution, Evolutionary programming, Evolutionary strategy,
Ant colony optimization and Particle swarm optimization for solving power system network
problems. This is as a result of complexity that arises from load flow studies in distributed alternative
energy sources especially in remote geographical areas.

3
2.0 AIMS AND OBJECTIVES

The aim of this course work can be described as follows:

1. To develop the skill in solving power flow distribution problems in an electrical

power network
2. To develop some algorithms that describes the approach used in solving power flow
problems in an electrical network.
3. To be able to use MATLAB in analysing and solving power flow problems.
4. To understand that Power flow problems are non linear and therefore cannot be
solved by classical circuit analysis.
5. To understand how to obtain the bus admittance matrix of a simple power network.
6. To understand the term power flow analysis as it applies to interconnected power
network.
7. To understand that the flow of both real and reactive power account for real power
losses in a transmission system
8. To understand how to calculate voltages in a non linear power network.
9. To understand how to calculate both the real and the reactive powers on a generator
once the busbar voltages have been determined.

3.0 INTRODUCTION

The voltage at the slack bus bar is made reference at 1 p.u. The remaining bus bars supply
loads with various real power (P) and imaginary power (Q) as given in Fig 1A below.

The series impedance of interconnecting lines is given in table 1 below. Each line consists of
capacitive shunt admittance. The series impedance and capacitive shunt admittance of the
interconnecting lines are to be ignored and included at different stages of the analysis.
Consequently, the shunt capacitance effect in relation to the resulting nodal voltage
distribution of the system will be determined.

The network base quantities are MVA base = 100 MVA and V base = 132 kV.

4
 Figure 1: Five – Bus/Node, Seven – Line Network Power Network

From Bus bar To Bus bar R (p.u.) X (p.u.) B (p.u.)

Table 1 showing busbar to busbar resistance, reactance and susceptance in per unit

The line admittance can be calculated as shown in Table 2 below:

5
 Table 1: Line Admittance Calculation

Line Admittance Conductance G (S) Susceptance B (S)

Table 2: Line Admittances

Table 4 below shows the total network admittance without effect of the shunt line capacitance

Admittance Matrix – shunt capacitance

Table 3: Admittance Matrix without shunt capacitance

6
In table 5 below, the admittance was calculated taking into consideration the effect of shunt
admittance. Every shunt admittances connected to a specific node was added to the diagonal
element of that node. Incorporating the shunt capacitive elements into the admittance matrix
reduces the magnitude of the complex element in the diagonal matrix as shown below:

Admittance Matrix + shunt capacitance

Table 4: Admittance Matrix with shunt capacitance added

7
4.0 Method 1 – Jacobi Method

The Jacobi approach is the easiest and straightforward of all the iterative power flow approach to
program. This involves the calculation of successive sets of bus bar voltages using the formerly
obtained values. The total range of bus bar voltages are calculated before been used in the following
iteration.

In Jacobi approach, a solution is normally obtained in almost around 29 iterations, which makes it
the longest among all the iterative techniques. The solution is obtained when there is convergence
of the bus bar voltages as will be shown in the later part of this report and the procedure for this
iterative solution is based on the general equation below:

4.1 Calculation:
The solution procedure is centred on the total network admittances and bus bar conditions and a
primary voltage referenced to each bus bar. A flow diagram of the iterative procedure is shown
below.

The number of nodes “n” = 5, and as the voltage at node (1) is stated, the real range of K = 2, 3, 4, 5.

8
4.2 Flowchart for JACOBI Method:

Figure 2: Flowchart for Jacobi Method

The step by step procedure is determined by numerically solving the first set of iterations for the
network shown in Figure 1.

4.3 Procedure:

9
A negative sign is assigned to indicate that the parameters being calculated are loads, and is
therefore taking supply from the source; otherwise, a positive signed is assigned.

10
11
12
Note: Shunt capacitance is included in the calculations above.

This step includes the replacement of all the values of with , and then repeating steps 1
through 3 until a convergence of two voltages is seen within a preselected tolerance.

4.4 MATLAB Code: JACOBI (Without shunt capacitance)

% Matlab Program.

% BEng(Hons) Electrical and Electronic Engineering

% Module 6EJ035

% Assignment No2

% MODELLING LOAD FLOW STUDIES USING MATLAB

% USING THE JACOBI METHOD

% YEAR 2012/13

% STUDENT: SALAU OLUWASESAN ABIODUN (100277737)

% NETWORK (see Fig 1)

% Input (1)
% Insert the Network Admittance Matrix
% (see page 2 of supplied assignment paperwork ignoring shunt capacitive
admittance)

format short
Y(1,1) = 10.958905 - 26.027398i;
Y(1,2) = -3.424658 + 7.534247i;
Y(1,3) = -3.424658 + 7.534247i;
Y(1,4) = 0.00;
Y(1,5) = -4.109589 + 10.958904i;

Y(2,1) = -3.424658 + 7.534247i;

Y(2,2) = 11.672080 - 26.090949i;
Y(2,3) = 0.00;
Y(2,4) = -4.123711 + 9.278351i;
Y(2,5) = -4.123711 + 9.278351i;

Y(3,1) = -3.424658 + 7.534247i;

Y(3,2) = 0.00;
Y(3,3) = 6.351487 - 13.87571i;
Y(3,4) = -2.926829 + 6.341463i;
Y(3,5) = 0.00;

Y(4,1) = 0.00;

13
Y(4,2) = - 4.123711 + 9.278351i;
Y(4,3) = - 2.926829 + 6.341463i;
Y(4,4) = 11.174251 - 24.898165i;
Y(4,5) = - 4.123711 + 9.278351i;

Y(5,1) = -4.109589 + 10.958904i;

Y(5,2) = -4.123711 + 9.278351i;
Y(5,3) = 0.00;
Y(5,4) = -4.123711 + 9.278351i;
Y(5,5) = 12.357012 - 29.515606i;

% Input (2)
% Input the given node loadings .....see diagram on page 1 ....
% .... of the supplied assignment paperwork
% A positive (+ sign) is for generated power,
% The negative (- sign) indicates load taken from the network.
% Per unit loading is used with a base loading taken at 100MVA

P2 = - 0.4; % p.u. active power loading

P3 = - 0.25; % p.u. active power loading
P4 = - 0.4; % p.u. active power loading
P5 = - 0.5; % p.u. active power loading

Q2 = - 0.2i; % p.u. reactive power loading

Q3 = - 0.15i; % p.u. reactive power loading
Q4 = - 0.2i; % p.u. reactive power loading
Q5 = - 0.2i; % p.u. reactive power loading

% The voltage at the generation busbar, node (1) is 1.0 p.u.

% This voltage is fixed and therefore determines the reference ....
% ... or SLACK BUSBAR must be Node(1)

% JACOBI SOLUTION
% The number of iteration will need to be set
% Set the value of 'l' the chosen number of iterations
% Typically l=40 should be sufficient for this problem.

l = 40; % sets the number of iterations

m = l + 1; % allocates the space used to store each iteration
vector = [1:m]; % assigns a row with the appropriate spaces

row = ones(size(vector));

Vnode1 = row; % each node voltage has the required storage spaces
Vnode2 = row;
Vnode3 = row;
Vnode4 = row;
Vnode5 = row;

for n = 1:m

Vnode1(n) = 1.0 + 0.00i; % inserts the assumed initial node voltages

Vnode2(n) = 1.0 + 0.00i;
Vnode3(n) = 1.0 + 0.00i;
Vnode4(n) = 1.0 + 0.00i;
Vnode5(n) = 1.0 + 0.00i;

end

S2star = P2 - Q2; % congugates of the specified node loadings

S3star = P3 - Q3;
S4star = P4 - Q4;
S5star = P5 - Q5;

% Now Proceed with the JACOBI SOLUTION

% Following the same method as in the worked example

14
% Page 5 of the assignment notes

for n = 1:l

% NOW UPDATE THE NODE VOLTAGES (l=40) TIMES

% SIMILAR TO THE CALCULATION PROCEDURES
% AS DETAILED ON PAGE 9 OF THE ASSIGNMENT NOTES

V2star(n) = conj(Vnode2(n)); % conjugate of V2

V3star(n) = conj(Vnode3(n)); % conjugate of V3
V4star(n) = conj(Vnode4(n)); % conjugate of V4
V5star(n) = conj(Vnode5(n)); % conjugate of V5

% (Step 1) Solving for the updated values of the currents

I2(n) = S2star/V2star(n);
I3(n) = S3star/V3star(n);
I4(n) = S4star/V4star(n);
I5(n) = S5star/V5star(n);

% (Step 2) Solving the summation of Y(kj)V(j) for all busbars

Sum2(n) = (Y(2,1) * Vnode1(n)) + (Y(2,4) * Vnode4(n)) + (Y(2,5) *

Vnode5(n));
Sum3(n) = (Y(3,1) * Vnode1(n)) + (Y(3,4) * Vnode4(n));
Sum4(n) = (Y(4,2) * Vnode2(n)) + (Y(4,3) * Vnode3(n)) + (Y(4,5) *
Vnode5(n));
Sum5(n) = (Y(5,1) * Vnode1(n)) + (Y(5,2) * Vnode2(n)) + (Y(5,4) *
Vnode4(n));

% (Step 3) Solving for the new value of node voltage

V2(n) = ( I2(n) - Sum2(n) ) / Y(2,2);

V3(n) = ( I3(n) - Sum3(n) ) / Y(3,3);
V4(n) = ( I4(n) - Sum4(n) ) / Y(4,4);
V5(n) = ( I5(n) - Sum5(n) ) / Y(5,5);

u = 1 + n;

Vnode2(u) = V2(n);
Vnode3(u) = V3(n);
Vnode4(u) = V4(n);
Vnode5(u) = V5(n);

end

disp (':')
disp (':')
disp (':')
disp (':')
disp (' Vnode2 Vnode3 Vnode4 Vnode5
')
disp(':')

d = 0;
for c = 1:l
fprintf(' %1.4f %1.4fi\t, %1.4f %1.4fi\t, %1.4f %1.4fi\t, %1.4f
%1.4fi\n',...
real(Vnode2(1+d)), imag(Vnode2(1+d)), real(Vnode3(1+d)),...
imag(Vnode3(1+d)), real(Vnode4(1+d)), imag(Vnode4(1+d)),
real(Vnode5(1+d)), imag(Vnode5(1+d)));
d=d+1;
end

plot(vector,Vnode2,vector,Vnode3,vector,Vnode4,vector,Vnode5);

15
 Vnode2 Vnode3 Vnode4 Vnode5
1 1.0000 0.0000i , 1.0000 0.0000i , 1.0000 0.0000i , 1.0000 0.0000i
2 0.9879 -0.0099i , 0.9842 -0.0108i , 0.9873 -0.0104i , 0.9882 -0.0120i
3 0.9789 -0.0178i , 0.9781 -0.0156i , 0.9741 -0.0213i , 0.9802 -0.0186i
4 0.9712 -0.0240i , 0.9719 -0.0206i , 0.9659 -0.0279i , 0.9731 -0.0247i
5 0.9656 -0.0285i , 0.9680 -0.0236i , 0.9586 -0.0338i , 0.9680 -0.0289i
6 0.9611 -0.0321i , 0.9646 -0.0263i , 0.9535 -0.0378i , 0.9638 -0.0323i
7 0.9576 -0.0347i , 0.9621 -0.0281i , 0.9493 -0.0410i , 0.9607 -0.0348i
8 0.9550 -0.0367i , 0.9601 -0.0296i , 0.9461 -0.0434i , 0.9583 -0.0367i
9 0.9529 -0.0382i , 0.9587 -0.0307i , 0.9436 -0.0452i , 0.9564 -0.0381i
10 0.9513 -0.0394i , 0.9575 -0.0315i , 0.9417 -0.0466i , 0.9549 -0.0392i
11 0.9501 -0.0403i , 0.9566 -0.0321i , 0.9402 -0.0476i , 0.9538 -0.0400i
12 0.9491 -0.0409i , 0.9559 -0.0326i , 0.9391 -0.0484i , 0.9529 -0.0407i
13 0.9484 -0.0414i , 0.9553 -0.0330i , 0.9382 -0.0490i , 0.9523 -0.0412i
14 0.9478 -0.0418i , 0.9549 -0.0332i , 0.9375 -0.0495i , 0.9517 -0.0415i
15 0.9474 -0.0421i , 0.9546 -0.0335i , 0.9370 -0.0498i , 0.9513 -0.0418i
16 0.9470 -0.0423i , 0.9544 -0.0336i , 0.9366 -0.0501i , 0.9510 -0.0420i
17 0.9468 -0.0425i , 0.9542 -0.0337i , 0.9363 -0.0503i , 0.9508 -0.0422i
18 0.9466 -0.0426i , 0.9540 -0.0338i , 0.9360 -0.0504i , 0.9506 -0.0423i
19 0.9464 -0.0427i , 0.9539 -0.0339i , 0.9358 -0.0506i , 0.9504 -0.0424i
20 0.9463 -0.0428i , 0.9538 -0.0340i , 0.9357 -0.0506i , 0.9503 -0.0425i
21 0.9462 -0.0429i , 0.9537 -0.0340i , 0.9356 -0.0507i , 0.9502 -0.0425i
22 0.9461 -0.0429i , 0.9537 -0.0340i , 0.9355 -0.0508i , 0.9502 -0.0425i
23 0.9461 -0.0429i , 0.9536 -0.0341i , 0.9354 -0.0508i , 0.9501 -0.0426i
24 0.9460 -0.0430i , 0.9536 -0.0341i , 0.9354 -0.0508i , 0.9501 -0.0426i
25 0.9460 -0.0430i , 0.9536 -0.0341i , 0.9353 -0.0509i , 0.9501 -0.0426i
26 0.9460 -0.0430i , 0.9536 -0.0341i , 0.9353 -0.0509i , 0.9500 -0.0426i
27 0.9459 -0.0430i , 0.9536 -0.0341i , 0.9353 -0.0509i , 0.9500 -0.0426i
28 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9353 -0.0509i , 0.9500 -0.0427i
29 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
30 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
31 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
32 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
33 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
34 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
35 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
36 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
37 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
38 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
39 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
40 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
Table 5: Command Window Output showing successive iteration values for JACOBI method (without shunt capacitance)

From the table above, it can be seen that the convergence occurs after 28 iteration, which
characterize the Jacobi method (without shunt capacitance) of solving load flow problems. A
graphical representation of the voltage at each of the nodes is shown below:

16
 Figure 3: Graph of solution to the voltages at each of the nodes

4.5 MATLAB CODE: JACOBI METHOD (WITH SHUNT CAPACITANCE)

%Matlab Program.

% BEng(Hons) Electrical and Electronic Engineering

% Module 6EJ035

% Assignment No2

% MODELLING LOAD FLOW STUDIES USING MATLAB

% USING THE JACOBI METHOD

% YEAR 2012/13

% STUDENT: SALAU OLUWASESAN ABIODUN (100277737)

% NETWORK (see Fig 1)

% Input (1)
% Insert the Network Admittance Matrix
% (see page 2 of supplied assignment paperwork with shunt capacitive admittance)

17
format short
Y(1,1) = 10.958905 – 25.997397i;
Y(1,2) = -3.424658 + 7.534247i;
Y(1,3) = -3.424658 + 7.534247i;
Y(1,4) = 0.00;
Y(1,5) = -4.109589 + 10.958904i;

Y(2,1) = -3.424658 + 7.534247i;

Y(2,2) = 11.672080 - 26.060948i;
Y(2,3) = 0.00;
Y(2,4) = -4.123711 + 9.278351i;
Y(2,5) = -4.123711 + 9.278351i;

Y(3,1) = -3.424658 + 7.534247i;

Y(3,2) = 0.00;
Y(3,3) = 6.351487 - 13.875707i;
Y(3,4) = -2.926829 + 6.341463i;
Y(3,5) = 0.00;

Y(4,1) = 0.00;
Y(4,2) = - 4.123711 + 9.278351i;
Y(4,3) = - 2.926829 + 6.341463i;
Y(4,4) = 11.174251 - 24.898165i;
Y(4,5) = - 4.123711 + 9.278351i;

Y(5,1) = -4.109589 + 10.958904i;

Y(5,2) = -4.123711 + 9.278351i;
Y(5,3) = 0.00;
Y(5,4) = -4.123711 + 9.278351i;
Y(5,5) = 12.357012 - 29.485605i;

% Input (2)
% Input the given node loadings .....see diagram on page 1 ....
% .... of the supplied assignment paperwork
% A positive (+ sign) is for generated power,
% The negative (- sign) indicates load taken from the network.
% Per unit loading is used with a base loading taken at 100MVA

P2 = - 0.4; % p.u. active power loading

P3 = - 0.25; % p.u. active power loading
P4 = - 0.4; % p.u. active power loading
P5 = - 0.5; % p.u. active power loading

Q2 = - 0.2i; % p.u. reactive power loading

Q3 = - 0.15i; % p.u. reactive power loading
Q4 = - 0.2i; % p.u. reactive power loading
Q5 = - 0.2i; % p.u. reactive power loading

% The voltage at the generation busbar, node (1) is 1.0 p.u.

% This voltage is fixed and therefore determines the reference ....
% ... or SLACK BUSBAR must be Node(1)

% JACOBI SOLUTION
% The number of iteration will need to be set
% Set the value of 'l' the chosen number of iterations
% Typically l=40 should be sufficient for this problem.

l = 40; % sets the number of iterations

m = l + 1; % allocates the space used to store each iteration
vector = [1:m]; % assigns a row with the appropriate spaces

row = ones(size(vector));

Vnode1 = row; % each node voltage has the required storage spaces
Vnode2 = row;
Vnode3 = row;

18
Vnode4 = row;
Vnode5 = row;

for n = 1:m

Vnode1(n) = 1.0 + 0.00i; % inserts the assumed initial node voltages

Vnode2(n) = 1.0 + 0.00i;
Vnode3(n) = 1.0 + 0.00i;
Vnode4(n) = 1.0 + 0.00i;
Vnode5(n) = 1.0 + 0.00i;

end

S2star = P2 - Q2; % congugates of the specified node loadings

S3star = P3 - Q3;
S4star = P4 - Q4;
S5star = P5 - Q5;

% Now Proceed with the JACOBI SOLUTION

% Following the same method as in the worked example
% Page 5 of the assignment notes

for n = 1:l

% NOW UPDATE THE NODE VOLTAGES (l=40) TIMES

% SIMILAR TO THE CALCULATION PROCEDURES
% AS DETAILED ON PAGE 9 OF THE ASSIGNMENT NOTES

V2star(n) = conj(Vnode2(n)); % conjugate of V2

V3star(n) = conj(Vnode3(n)); % conjugate of V3
V4star(n) = conj(Vnode4(n)); % conjugate of V4
V5star(n) = conj(Vnode5(n)); % conjugate of V5

% (Step 1) Solving for the updated values of the currents

I2(n) = S2star/V2star(n);
I3(n) = S3star/V3star(n);
I4(n) = S4star/V4star(n);
I5(n) = S5star/V5star(n);

% (Step 2) Solving the summation of Y(kj)V(j) for all busbars

Sum2(n) = (Y(2,1) * Vnode1(n)) + (Y(2,4) * Vnode4(n)) + (Y(2,5) *

Vnode5(n));
Sum3(n) = (Y(3,1) * Vnode1(n)) + (Y(3,4) * Vnode4(n));
Sum4(n) = (Y(4,2) * Vnode2(n)) + (Y(4,3) * Vnode3(n)) + (Y(4,5) *
Vnode5(n));
Sum5(n) = (Y(5,1) * Vnode1(n)) + (Y(5,2) * Vnode2(n)) + (Y(5,4) *
Vnode4(n));

% (Step 3) Solving for the new value of node voltage

V2(n) = ( I2(n) - Sum2(n) ) / Y(2,2);

V3(n) = ( I3(n) - Sum3(n) ) / Y(3,3);
V4(n) = ( I4(n) - Sum4(n) ) / Y(4,4);
V5(n) = ( I5(n) - Sum5(n) ) / Y(5,5);

u = 1 + n;

Vnode2(u) = V2(n);
Vnode3(u) = V3(n);
Vnode4(u) = V4(n);
Vnode5(u) = V5(n);

19
end

disp (':')
disp (':')
disp (':')
disp (':')
disp (' Vnode2 Vnode3 Vnode4 Vnode5 ')
disp(':')

d = 0;
for c = 1:l
fprintf(' %1.4f %1.4fi\t, %1.4f %1.4fi\t, %1.4f %1.4fi\t, %1.4f
%1.4fi\n',...
real(Vnode2(1+d)), imag(Vnode2(1+d)), real(Vnode3(1+d)),...
imag(Vnode3(1+d)), real(Vnode4(1+d)), imag(Vnode4(1+d)), real(Vnode5(1+d)),
imag(Vnode5(1+d)));
d=d+1;
end
plot(vector,Vnode2,vector,Vnode3,vector,Vnode4,vector,Vnode5);

Vnode2 Vnode3 Vnode4 Vnode5

1 1.0000 0.0000i , 1.0000 0.0000i , 1.0000 0.0000i , 1.0000 0.0000i
2 0.9888 -0.0104i , 0.9842 -0.0108i , 0.9885 -0.0109i , 0.9891 -0.0124i
3 0.9806 -0.0186i , 0.9786 -0.0158i , 0.9759 -0.0221i , 0.9817 -0.0193i
4 0.9733 -0.0250i , 0.9727 -0.0209i , 0.9684 -0.0209i , 0.9751 -0.0256i
5 0.9681 -0.0296i , 0.9691 -0.0241i , 0.9615 -0.0350i , 0.9703 -0.0299i
6 0.9638 -0.0333i , 0.9659 -0.0268i , 0.9567 -0.0391i , 0.9664 -0.0334i
7 0.9606 -0.0360i , 0.9636 -0.0287i , 0.9527 -0.0425i , 0.9664 -0.0334i
8 0.9581 -0.0381i , 0.9617 -0.0302i , 0.9497 -0.0449i , 0.9635 -0.0359i
9 0.9562 -0.0396i , 0.9603 -0.0313i , 0.9474 -0.0468i , 0.9612 -0.0379i
10 0.9546 -0.0408i , 0.9592 -0.0322i , 0.9456 -0.0482i , 0.9549 -0.0393i
11 0.9535 -0.0417i , 0.9584 -0.0328i , 0.9442 -0.0493i , 0.9580 -0.0404i
12 0.9526 -0.0424i , 0.9577 -0.0333i , 0.9431 -0.0501i , 0.9570 -0.0413i
13 0.9519 -0.0429i , 0.9572 -0.0337i , 0.9423 -0.0507i , 0.9561 -0.0419i
14 0.9513 -0.0433i , 0.9568 -0.0340i , 0.9416 -0.0511i , 0.9555 -0.0424i
15 0.9509 -0.0436i , 0.9565 -0.0342i , 0.9411 -0.0515i , 0.9550 -0.0428i
16 0.9506 -0.0438i , 0.9563 -0.0343i , 0.9407 -0.0518i , 0.9546 -0.0431i
17 0.9503 -0.0439i , 0.9561 -0.0345i , 0.9404 -0.0520i , 0.9543 -0.0433i
18 0.9501 -0.0441i , 0.9559 -0.0346i , 0.9402 -0.0521i , 0.9541 -0.0434i
19 0.9500 -0.0442i , 0.9558 -0.0346i , 0.9400 -0.0522i , 0.9539 -0.0436i
20 0.9499 -0.0442i , 0.9558 -0.0346i , 0.9399 -0.0523i , 0.9538 -0.0436i
21 0.9498 -0.0443i , 0.9558 -0.0347i , 0.9398 -0.0524i , 0.9537 -0.0437i
22 0.9497 -0.0443i , 0.9557 -0.0347i , 0.9397 -0.0525i , 0.9536 -0.0438i
23 0.9497 -0.0444i , 0.9556 -0.0348i , 0.9396 -0.0525i , 0.9535 -0.0438i
24 0.9496 -0.0444i , 0.9556 -0.0348i , 0.9396 -0.0525i , 0.9534 -0.0439i
25 0.9496 -0.0444i , 0.9556 -0.0348i , 0.9396 -0.0525i , 0.9534 -0.0439i
26 0.9496 -0.0444i , 0.9555 -0.0348i , 0.9395 -0.0526i , 0.9534 -0.0439i
27 0.9495 -0.0444i , 0.9555 -0.0348i , 0.9395 -0.0526i , 0.9534 -0.0439i
28 0.9495 -0.0445i , 0.9555 -0.0348i , 0.9395 -0.0526i , 0.9533 -0.0439i
29 0.9495 -0.0445i , 0.9555 -0.0348i ,0.9394 -0.0526i , 0.9533 -0.0439i
30 0.9495 -0.0445i , 0.9555 -0.0348i ,0.9394 -0.0526i , 0.9533 -0.0439i
31 0.9495 -0.0445i , 0.9555 -0.0349i ,0.9394 -0.0526i , 0.9533 -0.0439i
32 0.9495 -0.0445i , 0.9555 -0.0349i ,0.9394 -0.0526i , 0.9533 -0.0439i
33 0.9495 -0.0445i , 0.9555 -0.0349i ,0.9394 -0.0526i , 0.9533 -0.0439i

20
 34 0.9495 -0.0445i , 0.9555 -0.0349i ,0.9394 -0.0526i , 0.9533 -0.0439i
35 0.9459 -0.0445i , 0.9555 -0.0349i ,0.9394 -0.0526i , 0.9533 -0.0439i
36 0.9459 -0.0445i , 0.9555 -0.0349i ,0.9394 -0.0526i , 0.9533 -0.0439i
37 0.9459 -0.0445i , 0.9555 -0.0349i ,0.9394 -0.0526i , 0.9533 -0.0439i
38 0.9495 -0.0445i , 0.9555 -0.0349i ,0.9394 -0.0526i , 0.9533 -0.0439i
39 0.9495 -0.0445i , 0.9555 -0.0349i ,0.9394 -0.0526i , 0.9533 -0.0439i
40 0.9495 -0.0445i , 0.9555 -0.0349i ,0.9394 -0.0526i , 0.9533 -0.0439i
Table 6: Command Window Output showing successive iteration values for JACOBI method (with shunt capacitance)

It can be observed from the table above that the convergences occur after 29 iterations when the
effect of the shunt admittance is been considered. Same graph as that of the one obtained when the
effect of the shunt capacitance is ignored is also obtained on the command window of the matlab.

21
5.0 METHOD NO 2: Z- Matrix Approach using a block substitution method
The Y –matrix used with the Jacobi and Gauss Seidal approaches are typically known for their poor
convergence and in some cases, divergence occurs in successive iterations values. Much effort has
been initiated to overcome this problem with various conflicting results. A better approach to
solving power flow problems is the Z – matrix method. In Y – matrix approach, the network
relationship are expressed in an admittance matrix format.

Inverting the equation gives the Z – matrix as seen below:

After this has been done, an iterative strategy is adopted to produce a convergence to the required
set of system voltages.

5.1 Calculations:
Since bus bar 1 is taken as the slack bus bar, for which the voltage conditions are completely stated,
the remaining nodal admittance equations will have the form:

In this form, the Y-matrix is not square and therefore cannot be inverted. Rearranging:

This can be written as

Where is square and non-singular and can be inverted.

After inversion, the following equation was obtained:

22
Where the second summation is constant since is fixed and:

The necessary calculation strategy is one that uses these last two equations iteratively until a
convergence of the voltage solution is attained. The flowchart illustrating the form of the computer
program is shown below:

Figure 4: Flowchart of Z - Matrix Block Substitution Method

The successive approach to solving power flow problems using the Z – matrix approach is
mathematically shown below:

The – below indicates the matrix elements , where and

Note that the first row and column removed. Note that a resultant 4x4 matrix is obtained,
which is presented below:

23
The equivalent can be obtained by using the relation: . But, this process can be
very lengthy and cumbersome. A quicker way to obtain the is by using . The
screen shot below shows the procedure for obtaining a Z – matrix from a Y- matrix.

5.2 MATLAB Code for the Inverse of Y – Matrix

format short
Y(1,1) = 10.958905 - 26.027398i;
Y(1,2) = -3.424658 + 7.534247i;
Y(1,3) = -3.424658 + 7.534247i;
Y(1,4) = 0.00;
Y(1,5) = -4.109589 + 10.958904i;

Y(2,1) = -3.424658 + 7.534247i;

Y(2,2) = 11.672080 - 26.090949i;
Y(2,3) = 0.00;
Y(2,4) = -4.123711 + 9.278351i;
Y(2,5) = -4.123711 + 9.278351i;

Y(3,1) = -3.424658 + 7.534247i;

Y(3,2) = 0.00;
Y(3,3) = 6.351487 - 13.87571i;
Y(3,4) = -2.926829 + 6.341463i;
Y(3,5) = 0.00;

Y(4,1) = 0.00;
Y(4,2) = - 4.123711 + 9.278351i;
Y(4,3) = - 2.926829 + 6.341463i;
Y(4,4) = 11.174251 - 24.898165i;
Y(4,5) = - 4.123711 + 9.278351i;

Y(5,1) = -4.109589 + 10.958904i;

Y(5,2) = -4.123711 + 9.278351i;
Y(5,3) = 0.00;
Y(5,4) = -4.123711 + 9.278351i;
Y(5,5) = 12.357012 - 29.515606i;

Ysub = Y(2:5,2:5);

Z = inv(Ysub)

24
The command window for the resulting 4x4 matrix is shown below:

Figure 5: Command window Screenshot for Z-matrix

The first step involves calculating the constant , for

25
5.1.1

26
5.1.2
Calculating all busbar voltages

27
5.3 MATLAB CODE FOR Z- MATRIX
% BEng(Hons) Electrical Electronic Engineering

% Module 6EJ503

% Assignment No2

28
% MODELLING LOAD FLOW STUDIES USING MATLAB

% YEAR 2012/13

% STUDENT: SALAU OLUWASESAN ABIODUN (100277737)

% NETWORK (see Fig 1)

% Input (1)
% Insert the Network Admittance Matrix
% (see page 2 of supplied assignment paperwork ignoring shunt capacitive
admittance)

format short
Y(1,1) = 10.958905 - 26.027398i;
Y(1,2) = -3.424658 + 7.534247i;
Y(1,3) = -3.424658 + 7.534247i;
Y(1,4) = 0.00;
Y(1,5) = -4.109589 + 10.958904i;

Y(2,1) = -3.424658 + 7.534247i;

Y(2,2) = 11.672080 - 26.090949i;
Y(2,3) = 0.00;
Y(2,4) = -4.123711 + 9.278351i;
Y(2,5) = -4.123711 + 9.278351i;

Y(3,1) = -3.424658 + 7.534247i;

Y(3,2) = 0.00;
Y(3,3) = 6.351487 - 13.87571i;
Y(3,4) = -2.926829 + 6.341463i;
Y(3,5) = 0.00;

Y(4,1) = 0.00;
Y(4,2) = - 4.123711 + 9.278351i;
Y(4,3) = - 2.926829 + 6.341463i;
Y(4,4) = 11.174251 - 24.898165i;
Y(4,5) = - 4.123711 + 9.278351i;

Y(5,1) = -4.109589 + 10.958904i;

Y(5,2) = -4.123711 + 9.278351i;
Y(5,3) = 0.00;
Y(5,4) = -4.123711 + 9.278351i;
Y(5,5) = 12.357012 - 29.515606i;

% Input (2)
% Input the given node loadings .....see diagram on page 1 ....
% .... of the supplied assignment paperwork
% A positive (+ sign) is for generated power,
% The negative (- sign) indicates load taken from the network.
% Per unit loading is used with a base loading taken at 100MVA

P2 = - 0.4; % p.u. active power loading

P3 = - 0.25; % p.u. active power loading
P4 = - 0.4; % p.u. active power loading
P5 = - 0.5; % p.u. active power loading

Q2 = - 0.2i; % p.u. reactive power loading

Q3 = - 0.15i; % p.u. reactive power loading
Q4 = - 0.2i; % p.u. reactive power loading

29
Q5 = - 0.2i; % p.u. reactive power loading

% The voltage at the generation busbar, node (1) is 1.0 p.u.

% This voltage is fixed and therefore determines the reference ....
% ... or SLACK BUSBAR must be Node(1)

% Z-MATRX SOLUTION
% The sub-matrix of the actual 5x5 network admittances (Y) ....
% but with the first row and column removed now needs to be obtained
% This sub-matrix is called Ysub and is easily derived from ....
% the network admittance matrix by using the following command

Ysub = Y(2:5,2:5);

% The inverse of this sub-matrix gives the required Z-matrix

Z = inv(Ysub);

% Some care needs to be taken in respect of the Z-matrix suffixes

% This can now be rectified by adding a further ....
% row(1) and column(1) of zeros to the Z-matrix.
% Note the prime (') after the matrix is a transpose command
% Rows are changed columns and columns become rows

A = [0;0;0;0];
Z = [A Z];
B = [0;0;0;0;0];
Zt = [B Z'];

% Zt is the transpose of the Z-matrix with an added set of zeros

% for the first row and column
% A further transposition brings back the modified Z-matrix

Z = Zt';

% The number of iteration will need to be set

% Set the value of 'l' the chosen number of iterations
% l=40 will be used since this is the number in the JACOBI solution.

l = 40; % sets the number of iterations

m = l + 1; % allocates the space used to store each iteration
vector = [1:m]; % assigns a row with the appropriate spaces

row = ones(size(vector));

Vnode1 = row; % each node voltage has the required storage spaces
Vnode2 = row;
Vnode3 = row;
Vnode4 = row;
Vnode5 = row;

for n = 1:m

Vnode1(n) = 1.0 + 0.00i; % inserts the assumed initial node voltages

Vnode2(n) = 1.0 + 0.00i;
Vnode3(n) = 1.0 + 0.00i;
Vnode4(n) = 1.0 + 0.00i;
Vnode5(n) = 1.0 + 0.00i;

30
end

S2star = P2 - Q2; % conjugates of the specified node loadings

S3star = P3 - Q3;
S4star = P4 - Q4;
S5star = P5 - Q5;

% Now Proceed with the Z-MATRIX SOLUTION

% Following the same method as in the worked example
% Pages 28 and 29 of the assignment notes

% STEP 1.
% Obtain the values of the constant terms ....
% from the summation of Z(kj)Y(j1)V(1)

for n = 1 : m
Const2 = Z(2,2) * Y(2,1) * Vnode1(n) + Z(2,3) * Y(3,1) * Vnode1(n) + Z(2,4)
* Y(4,1) * Vnode1(n) + Z(2,5) * Y(5,1) * Vnode1(n);
Const3 = Z(3,2) * Y(2,1) * Vnode1(n) + Z(3,3) * Y(3,1) * Vnode1(n) + Z(3,4)
* Y(4,1) * Vnode1(n) + Z(3,5) * Y(5,1) * Vnode1(n);
Const4 = Z(4,2) * Y(2,1) * Vnode1(n) + Z(4,3) * Y(3,1) * Vnode1(n) + Z(4,4)
* Y(4,1) * Vnode1(n) + Z(4,5) * Y(5,1) * Vnode1(n);
Const5 = Z(5,2) * Y(2,1) * Vnode1(n) + Z(5,3) * Y(3,1) * Vnode1(n) + Z(5,4)
* Y(4,1) * Vnode1(n) + Z(5,5) * Y(5,1) * Vnode1(n);

% STEP 2.
% Obtain the values of the node current input terms ....
% from Skstar/Vkstar

V2star(n) = conj(Vnode2(n));
V3star(n) = conj(Vnode3(n));
V4star(n) = conj(Vnode4(n));
V5star(n) = conj(Vnode5(n));

I2(n) = S2star / V2star(n);

I3(n) = S3star / V3star(n);
I4(n) = S4star / V4star(n);
I5(n) = S5star / V5star(n);

% STEP 3
% Calculate the new (update) value of node Voltage

V2(n) = Z(2,2) * I2(n) + Z(2,3) * I3(n) + Z(2,4) * I4(n) + Z(2,5) * I5(n) -

Const2;
V3(n) = Z(3,2) * I2(n) + Z(3,3) * I3(n) + Z(3,4) * I4(n) + Z(3,5) * I5(n) -
Const3;
V4(n) = Z(4,2) * I2(n) + Z(4,3) * I3(n) + Z(4,4) * I4(n) + Z(4,5) * I5(n) -
Const4;
V5(n) = Z(5,2) * I2(n) + Z(5,3) * I3(n) + Z(5,4) * I4(n) + Z(5,5) * I5(n) -
Const5;

u = 1 + n;

Vnode2(u) = V2(n);
Vnode3(u) = V3(n);
Vnode4(u) = V4(n);
Vnode5(u) = V5(n);

end

31
disp (':')
disp (' Electrical Power Application 6EJ503 ')
disp (':')
disp (' Z-MATRIX Solution of Load Flow Assignment(2) 2012/13')
disp (':')
disp (' Successive Iteration Values ')
disp (':')
disp (' Vnode2 Vnode3 Vnode4
Vnode5 ')
disp(':')

d = 0;
for c = 1:m
fprintf(' %1.4f %1.4fi\t, %1.4f %1.4fi\t, %1.4f %1.4fi\t, %1.4f
%1.4fi\n',...
real(Vnode2(1+d)), imag(Vnode2(1+d)), real(Vnode3(1+d)),...
imag(Vnode3(1+d)), real(Vnode4(1+d)), imag(Vnode4(1+d)),
real(Vnode5(1+d)), imag(Vnode5(1+d)));
d=d+1;
end
plot (vector,V2, 'r-', vector, V3, 'b-', vector,V4, 'g-', vector, V5, 'y-');

The successive iteration values for the Z – Matrix method is displayed in the table below:

Vnode2 Vnode3 Vnode4 Vnode5

1 1.0000 0.0000i , 1.0000 0.0000i , 1.0000 0.0000i , 1.0000 0.0000i
2 0.9508 -0.0431i , 0.9574 -0.0342i , 0.9412 -0.0510i , 0.9546 -0.0426i
3 0.9462 -0.0428i , 0.9537 -0.0340i , 0.9356 -0.0507i , 0.9502 -0.0425i
4 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
5 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
6 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
7 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
8 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
9 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
10 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
11 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
12 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
13 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
14 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
15 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
16 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
17 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
18 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
19 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
20 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
21 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
22 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
23 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
24 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
25 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
26 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
27 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i

32
 28 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
29 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
30 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
31 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
32 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
33 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
34 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
35 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
36 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
37 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
38 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
39 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
40 0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
0.9459 -0.0430i , 0.9535 -0.0341i , 0.9352 -0.0509i , 0.9500 -0.0427i
Table 9: Z – Matrix Iteration table

In contrary to the JACOBI method, the Z – Matrix method reached convergence after just 4
iterations as shown in the table above. Therefore, it can be inferred that the Z – Matrix method is
faster than the JACOBI method. The graph of the Z matrix method is presented below:

33
Figure 6: Graph of solution to the voltages at each of the nodes

34
6.0 CONCLUSIONS
The Power flow study was done using two iterative approaches – the Y- matrix approach applying
the Jacobi algorithm and the Z – matrix approach using the block substitution algorithm. It can be
deduced from the results obtained that the Z- matrix method reaches convergence faster than the Y
– matrix approach (Jacobi). As it is apparent from the two methods used, computer based analysis is
important for achieving accurate power flow solutions of any practically sized power network. A
computer based analysis usually applies divers’ numerical techniques such as sparsity and optimal
ordering techniques so as to reduce memory and storage needs. There are various excellent power
flow programs available that are vastly used by engineers in the industry (power) for periodic system
studies. While industry – grade power flow software is expensive, CAD software like the MATLAB is
an inexpensive and adequate alternative in solving and studying load flow problems in a small scale
networks.

35
7.0 REFERENCES
< https://www.courseresources.derby.ac.uk/bbcswebdav/pid-783447-dt-content-rid-
3525226_1/courses/2012-
6EJ503/EPA_Power%20Flow%20Studies%20using%20Advanced%20Computing%20Practical_2_2012
-2013.pdf > [Accessed on 25th March 2013].

< http://www.scribd.com/doc/36050693/PSCh9LoadFlow > [Accessed on 29th March 2013].

< http://www.scribd.com/doc/15853061/LoadFlow-Analysis-in-Power-Systems > [Accessed on 1st

April 2013]

< http://www.en.wikipedia.org/wiki/power_flow_study > [Accessed on 2nd April 2013]

< http://www.emersonnetworkpower.com/en-
us/services/ERS/PLCservices/Design/Engineeringservices/pages/LoadFlowStudies.aspx > [Accessed
on 2nd April 2013]

< http://www.powerapps.org/PAES_LoadFlow.aspx > [Accessed on 3rd April 2013]

< http://www.intechopen.com/books/simulated-annealing-advances-applications-and-
hybridizations/applications-of-simulated-annealing-based-approaches-to-electric-power-systems >
[Accessed on 5th April 2013]

Alsaadi, and B. Gholami, [2009], An Effective Approach for Distribution System Power Flow Solution,
World Academy of Science, Engineering and Technology, [Internet] Available from: <
http://www.waset.org/journals/waset/v49/v49-45.pdf> [Accessed 8th April 2013]

Trani, [2001], Introduction to MATLAB MATLAB Matrix Manipulations and Graphics, Civil and
Environmental Engineering Virginia Polytechnic Institute and State University [Internet] Available
from: < http://128.173.204.63/courses/matlab/matlab_matrices.pdf> [Accessed 8th April 2013]

Badrul .H. Chowdhury, [2004], Load – Flow Analysis In Power Systems, Electrical & Computer
Engineering Department, University of Missouri-Rolla [Internet] Available from:
<http://www.accessengineeringlibrary.com/mghpdf/007144940x_ar011.pdf >[Accessed on 10th April
2013]

36
37
38