International Journal of Scientific & Engineering Research, Volume 5, Issue 5, May-2014 1412

ISSN 2229-5518

Comparison between Load Flow Analysis Methods in Power System using MATLAB

Kriti Singhal

Abstract— Now these days load flow is a very important and fundamental tool for the analysis of any power systems and in the operations as well as planning stages. Certain applications, particularly in distribution automation and optimization of a power system, require repeated load flow solutions. In these applications it is very important to solve the load flow problem as efficiently as possible. Since the invention and widespread use of digital computers and many methods for solving the load flow problem have been developed. Most of the methods have “grown up” around transmission systems and, over the years, variations of the Newton method such as the fast decoupled method; have become the most widely used. Some of the methods based on the general meshed topology of a typical transmission system are also ap- plicable to distribution systems which typically have a radial or tree structure. Specifically, we will compare the proposed method to the standard Newton Method, and the implicit Z bus Gauss method.

Our goal is to develop a formulation and solutions algorithm for solving load flow in large 3-phase unbalanced systems which exploits the radial topological structure to reduce the number of equations and unknowns and the numerical structure by using Gauss-Seidel Method, Newton Raphson Method and Fast Decoupled Method.

Index Terms— Fast Decoupled, Gauss-Seidel, Load Flow, Newton Raphson, Numerical Analysis, Power System, Solutions Algorithm, Z- bus.

1 INTRODUCTION

——————————  ——————————

OAD flow study also known as power flow study, is an important tool involving numerical analysis applied to a power system. A power-flow study usually uses simpli-
fied notation such as a one-line diagram and per-unit system, and focuses on various forms of AC power (i.e.: voltages, volt- age angles, real power and reactive power). Load-flow studies are performed to determine the steady-state operation of an electric power system. It calculates the voltage drop on each feeder, the voltage at each bus, and the power flow in all branch and feeder circuits. Determine if system voltages re- main within specified limits under various contingency condi- tions, and whether equipment such as transformers and con- ductors are overloaded. It is used to identify the need for addi- tional generation, capacitive, or inductive support, or the placement of capacitors and/or reactors to maintain system voltages within specified limits. Losses in each branch and total system power losses are also calculated. It is Necessary for planning, economic scheduling, and control of an existing system as well as planning its future expansion.
There are two popular numerical methods for solving the power-flow equations. These are the Gauss-Seidel (G-S) and the Newton-Raphson (N-R) Methods (Grainger and Steven- son, 1994; Elgend, 1982; Glover and Sharma, 1994). The N-R method is superior to the G-S method because it exhibits faster convergence characteristics. However, the N-R method suffers from the disadvantages that a “flat start” is not always possi- ble since the solution at the beginning can oscillate without converging towards the solution. In order to avoid this prob- lem, the load-flow solution is often started with a G-S algo-

————————————————

• Kriti Singhal is pursuing B-tech in electrical & electronics engineering in

Merrut Institute of Technology.E-mail: singhal.kriti23@gmail.com
rithm followed by the N-R algorithm after a few iterations. There is also an approximate but faster method for the load- flow solution. It is a variation of the N-R method, called the faster-decoupled method, which was introduced.

2 POWER FLOW PROBLEM FORMULATION

The goal of a power flow study is to obtain complete voltage angle and magnitude information for each bus in a power sys- tem for specified load and generator real power and voltage conditions. Once this information is known, real and reactive power flow on each branch as well as generator reactive pow- er output can be analytically determined[1].

The solution to the power flow problem begins with identify- ing the known and unknown variables in the system. The known and unknown variables are dependent on the type of bus. A bus without any generators connected to it is called a Load Bus. A bus with at least one generator connected to it is called a Generator Bus. The exception is one arbitrarily- selected bus that has a generator. This bus is referred to as the slack bus.

In the power flow problem,if the real power and reactive pow- er at each Load Bus are known. For this reason, Load Buses are also known as PQ Buses. For Generator Buses, it is as- sumed that the real power generated PG and the voltage mag- nitude V is known. For the Slack Bus, it is assumed that the voltage magnitude |V| and voltage phase Θ are known. Therefore, for each Load Bus, the voltage magnitude and an- gle are unknown and must be solved for; for each Generator Bus, the voltage angle must be solved for; there are no varia- bles that must be solved for the Slack Bus. In a system with N buses and R generators, there are then

International Journal of Scientific & Engineering Research, Volume 5, Issue 5, May-2014 1413

ISSN 2229-5518

2(N − 1) − (R − 1)unknowns

I i = yioVi + yi1 (Vi − V1 ) + yi 2 (Vi − V2

)+ ... + y (V − V )

In order to solve for the above equation, the possible equations to use are power balance equations, which can be written for real and reactive power for each bus. The real power balance

= (yio + yi1 + yi 2 + ... + yin )Vi − yi1Vi − yi 2V2 − ... − yinV

Representing eq. (1) in summation form;

(1)

equation is:

I i = Vi ∑ yij

− ∑ yijV j

j ≠ i

(2)

0 = − Pi + ∑ Vi Vk (Gik cosθ ik + Bik sin θ ik )

The complex power at its ith bus is;

j = 0

j =1

k =1

P + jQ

= V I ∗

(3)
Where, Pi is the net power injected at bus i,

i i i i

This is for reactive power.

Gik

is the real part of the element in the bus admit-

P − jQ

tance matrix,

I = i i

i V

(4)
YBUS corresponding to the ith row and kth column, i

Bik is the imaginary part of the element. The reactive power balance equation is:

This is for active power.
Substituting for Ii in (2) yields;

N P − jQ n n

0 = −Qi + ∑ Vi Vk (Gik sin θ ik − Bik cosθ ik )

V = Vi ∑ yij − ∑ yijV j

j ≠ i

(5)

k =1

i∗ j =0

j =1

Where, Qi is the net reactive power injected at bus i,
Equation (5) is an algebraic non linear equation which must be
and θ ik

= δ i − δ k .

solved by iterative techniques. We use Methods like Gauss- Seidel, Newton Raphson and Fast Decoupled Method for fur-
Equations included are the real and reactive power balance
equations for each Load Bus and the real power balance equa-
tion for each Generator Bus. Only the real power balance
equation is written for a Generator Bus because the net reac-
tive power injected is not assumed to be known and therefore
including the reactive power balance equation would result in

an additional unknown variable. For similar reasons, there are
no equations written for the Slack Bus.
ther solution[2].

4 GAUSS-SEIDEL METHOD

The Gauss–Seidel method, also known as the Liebmann meth- od or the method of successive displacement, is an iterative method used to solve a linear system of equations. It is named after the German mathematicians Carl Friedrich Gauss and Philipp Ludwig von Seidel, and is similar to the Jacobi method.

4.1 For P-Q bus

Equation (5) is solved for Vi solved iteratively;

sch i

− jQ sch

(k )

V ∗(k )

+ ∑ yijV j

V (k +1) = i

i y

j ≠ i

(6)

∑ ij

Fig. 1- Power System Analysis.

3 POWER FLOW EQUATION

For Generator buses or P-V bus (where real and reactive pow- ers are injected), Pisch and Q isch have positive values. Load bus- es or P-Q bus (real and reactive powers flow away from the bus), P isch and Qi sch have negative values.
Eq (5) can be solved for Pi and Q i;

P (k +1) =



RVi



∗(k ) 

Vi



(k )

∑ yij

j =0

− ∑ yijV

j =1

(k ) 

j 



j ≠ i

(7)

Q (k +1) = −



ImVi



∗(k ) 

Vi



(k )

∑ yij

j =0

− ∑ yijV

j =1

(k ) 

j 



j ≠ i

(8)
Fig 2- A typical bus of the power system
Applying KCL to this bus results in;
The power flow equation is usually expressed in terms of the elements of the bus admittance matrix, Ybus , shown by upper case letters, are Yij = -yij , and the diagonal elements are Yii = ∑ yij .
Hence eq (6), (7) and (8) can be written as;

International Journal of Scientific & Engineering Research, Volume 5, Issue 5, May-2014 1414

ISSN 2229-5518

sch i

− jQ sch

− ∑

Y V (k )

V (k +1) =

∗(k )

j ≠i ij j

(9)

(k +1)

 

 ∗(k ) 
(k )



n (k ) 

Q (k +1)

= RVi







= − ImV





∗(k ) 

V

Yii +

(k )Y

∑YijV j

j =1 j ≠ i

+ ∑Y V





( k ) 



j ≠ i

(10)
(11)

i i i

 

ii ij j

J =1 

Equation (5) is solved untill Vik+1 -Vi k <0.0001.

4.2 For P-V bus

Firstly Qi k+1 is calculated from eq (11);

(k +1) = −

 ∗(k ) 

(k ) +

(k ) 

Qi ImVi



Vi



Yii

∑YijV j

J =1J ≠ i





j ≠ i

For a P-V bus the upper limit Qmax and lower limit Qmin of Q to hold the generations vars within limits is also given. That is
Qi(min) < Qi < Qi(ma x)
Therefore, the calculated Qi(k+1) is checked for limits of Qi , that is
Qi(min) < Qi (k+1) < Qi(max)
If Qi(k+1) is within its limits then calculate new value of Ck and Vi(k+1) from eq (6) and treat the ith bus as P-Q bus. The compu- tation is then continued as in case of a P-Q bus.[2]

4.3 Flow Diagram

Fig 3- Flow chart of Gauss-Seidel Method[3].

International Journal of Scientific & Engineering Research, Volume 5, Issue 5, May-2014 1415

ISSN 2229-5518

5 NEWTON RAPHSON METHOD

∆Q(k )

= Q sch

− Q(k )

(21)

If you've ever tried to find a root of a complicated function algebraically, you may have had some difficulty. Using some basic concepts of calculus, we have ways of numerically eval- uating roots of complicated functions. Commonly, we use the Newton-Raphson method. This iterative process follows a set guideline to approximate one root, considering the function, its derivative, and an initial x-value.
Power flow equations formulated in polar form. For the sys- tem in Fig 2, Eqn (2) can be written in terms of bus admittance matrix as;[2]

Continue until scheduled errors ΔPi(k) and ΔQi(k) for all load buses are within a specified tolerance.

5.1 Figures and Tables

I i = ∑ YijV j

j =1

(12)
Expressing in polar form;

I i = ∑ Vij j =1

V j ∠θ ij

+ δ j

(13)
Substituting for I i from Eqn (21) in Eqn (4);

Pi − jQi

= Vi ∠ − δ i

∑ Vij V j ∠θ ij + δ j

j =1

(14)
Separating the real and imaginary parts;

Pi =

∑ Vi V j

j =1

Vij

cos(θ

− δ i

+ δ j )

(15)

Qi = −

∑ Vi V j

j =1

V sin θ

− δ i

+ δ j )

(16)

δ i is phase angle.

Expanding Eqn (15) & (16) in Taylor's series about the initial estimate neglecting high order terms we get;

The Jacobian matrix gives the linearized relationship between
small changes in Δδi (k) and voltage magnitude Δ[Vi k] with the
small changes in real and reactive power ΔPi (k) and ΔQi(k) .
Fig 4- Flow chart of Newton Raphson Method[3].

 ∆P 

 J 1

J 2   ∆δ 

(17)

6 FAST DECOUPLED METHOD

  = 

  

∆Q

 J 3

J 4  ∆ V 

Fast decoupled load flow method is a variation on Newton-
The diagonal and the off-diagonal elements of J1 are,
Raphson that exploits the approximate decoupling of active

∂Pi

∂δ i

= ∑ Vi V j

j ≠i

Yij

cos(θ

− δ i

+ δ j )

(18)
and reactive flows in well-behaved power networks, and addi- tionally fixes the value of the Jacobian during the iteration in

∂Pi

∂δ j

= − Vi V j

Yij

sin (θ

− δ i

+ δ j )

(19)
order to avoid costly matrix decompositions.
In this power transmission lines have high X/R ratio. Real
Similarly we can find the diagonal and off-diagonal elements of J2,J3 and J4;
The terms ΔPi(k) and ΔQi(k) are the difference between the scheduled and calculated values, known as the power resid- uals.
power changes are less sensitive to voltage magnitude chang-
es and are most sensitive to changes in phase angle Δδ. Simi-
larly, reactive power changes are less sensitive to changes in
angle and are mainly dependent on changes in voltage magni-
tude. Therefore the Jacobian matrix in Equation (17) can be

∆Pi

(k )

= Pi

sch

− Pi

(k )

(20)
written as;

International Journal of Scientific & Engineering Research, Volume 5, Issue 5, May-2014 1416

ISSN 2229-5518

 ∆P 

 J 1

0   ∆δ 

(22)
Y(1,5)=X(2)+1i*0.020;

  = 

  

∆Q  0

J 4  ∆ V 

 ∂P 

Y(2,3)=X(3)+1i*0.025;
Y(2,5)=X(4)+1i*0.020;

∆P = J 1 ∆δ = 

 ∆δ

(23)
Y(3,4)=X(5)+1i*0.020;

 ∂δ 

Y(3,5)=X(6)+1i*0.010;

  Y(4,5)=X(7)+1i*0.075;

∆Q = J

4 ∆ V

∂Q

=  ∆ V

 ∂ V 

(24)
Y(1,1)=Y(1,2)+Y(1,5); Y(1,3)=0;
The diagonal elements of J1 given by Eqn.(18) is written as;

Y(1,4)=0;

∂Pi =

V V Y

cos(θ
− δ + δ

)− V 2 Y

sin θ
Y(2,1)=Y(1,2);
∂δ i

∑ i

j =1

j ij

ij i

j i ii ii

Y(2,2)=Y(1,2)+Y(2,3);
Y(2,4)=0;
Replacing the first term with –Qi from (19)
Y(3,1)=0;

∂Pi

∂δ i

= −Qi

− Vi

Yii

sin θ ii

= −Qi

− Vi

Bii

Y(3,2)=Y(2,3); Y(3,3)=Y(2,3)+Y(3,4)+Y(3,5);
Bii = sum of susceptances of the entire elements incident to bus i. In a typical power system, B ii » Q i therefore we may neglect Qi .
Furthermore, [Vi]2 ≈ [Vi] Ultimately;
Y(4,1)=0;
Y(4,2)=0;
Y(4,3)=Y(3,4);
Y(4,4)=Y(3,4)+Y(4,5);
Y(5,1)=Y(1,5);

∂Pi

∂δ i

= − Vi

Bii

(25)
Y(5,2)=Y(2,5); Y(5,3)=Y(3,5);
In equation (19) assuming θii-δi+δj ≈ θii, the off diagonal ele- ments of J1 becomes,
Y(5,4)=Y(4,5); Y(5,5)=Y(1,5)+Y(2,5)+Y(3,5)+Y(4,5);

∂Pi

∂δ j

= − Vi V j

Bij

[r,c]=size(Y);
Y1=Y;

for i=1:r

Assuming [Vi] ≈ 1 we get;

for j=1:c

∂Pi

∂δ j

= − V j

Bij

(26)

if i~=j

Y(i,j)=-Y(i,j);

end

Similarly we can simplify the diagonal and off-diagonal ele-
ments of J4 as,

end end

∂Qi
∂ Vi

= − Vi

Bii

∂Q

; i
∂ V j

= − Vi

Bij

theta=angle(Y);
del_tolerance=1;
With these assumptions, equations (25) & (26) can be written in the following form,
V=[1.05 1 1 1 1.02];
del=[0 0 0 0 0];
sp=0;

∆P = − B′∆δ

∆Q

= − B ′ ∆ Vi

(27)
(28)
sq=0;

for k=1:5

for n=1:5

sp=sp+abs(V(k)*V(n)*Y(k,n))*cos(theta(k,n)+del(n)-del(k));
B’ and B’’ are the imaginary part of the bus admittance matrix Ybus . Since the elements of the matrix are constant, need to be triangularized and inverted only once at the beginning of the iteration[2].

7 COADING IMPIMENTATION

clc
clear all
close all
x=[0.02+1i*0.10 0.05+1i*0.25 0.04+1i*0.20 0.05+1i*0.25
0.05+1i*0.25 0.08+1i*0.40 0.10+1i*0.50];
X=1./x;
Y=zeros(5);
Y(1,2)=X(1)+1i*0.030;

end

P(k)=sp;
P(5)=24;

end

figure('color',[0 .5 1],'menubar','none')
plot(1:5,P,'color','r','marker','.','markeredgecolor','k','linewidth',
2)
hold on

for k=1:5

for n=1:5

sq=-(sq+abs(V(k)*V(n)*Y(k,n))*sin(theta(k,n)+del(n)-del(k)));

end

Q(k)=sq;
Q(5)=11;

end

plot(1:5,Q,'color','k','marker','.','markeredgecolor','r','linewidth',
2)

International Journal of Scientific & Engineering Research, Volume 5, Issue 5, May-2014 1417

ISSN 2229-5518

grid on
xlabel('Number of buses');
ylabel('Powers')
legend('Real Power','Reactive power','location','best')

%Load and Generator

PsL=[0 96 35 16 24];
del(k));

for N=[1 k+1:5];

t1=t1+abs(V(k)*V(N)*Y(k,N))*cos(theta(k,N)+del(N)-

end

JQD(k-2,n-1)=t1;
t1=0;
QsL=[0 62 14 8 11];
PsG=[NaN 0 0 0 48];
QsG=[NaN 0 0 0 NaN];
P2sch=(PsG(2)-PsL(2))/100;
P3sch=(PsG(3)-PsL(3))/100;
P4sch=(PsG(4)-PsL(4))/100;
P5sch=(PsG(5)-PsL(5))/100;
Q2sch=(QsG(2)-QsL(2))/100;
Q3sch=(QsG(3)-QsL(3))/100;
Q4sch=(QsG(4)-QsL(4))/100;

%mismatching in power calculations

DelP2=P2sch-P(2);
DelP3=P3sch-P(3);
DelP4=P4sch-P(4);
DelP5=P5sch-P(5);
DelQ2=Q2sch-Q(2);
DelQ3=Q3sch-Q(3);
DelQ4=Q4sch-Q(4);
final_vector=[DelP2 DelP3 DelP4 DelP5 DelQ2 DelQ3 DelQ4];

%% formation of Jacobian

s=0;
idx=0;
t=0;
t1=0;
t2=0;

for k=2:5

for n=2:5

if k==n

for N=[1 k+1:5]

s=s+abs(V(k)*V(N)*Y(k,N))*sin(theta(k,N)+del(N)-
del(k));

end

JPD(k-1,n-1)=s;
s=0;

elseif k~=n

JPD(k-1,n-1)=-abs(V(k)*V(n)*Y(k,n))*sin(theta(k,n)-
del(k)+del(n));

end

if n>=3

if k==n

for N=[1 k+1:5]

t=t+abs(V(N)*Y(k,N))*cos(theta(k,N)+del(N)-del(k));

end

t1=2*abs(V(k)*Y(k,n))*cos(theta(k,n))+t;
JPV(k-1,n-2)=t1;
t=0;

elseif k~=n

JPV(k-1,n-2)=abs(V(k)*Y(k,n))*cos(theta(k,n)-
del(k)+del(n));

end

if k>=3

if k==n

elseif k~=n

JQD(k-2,n-1)=-abs(V(k)*V(n)*Y(k,n))*cos(theta(k,n)-
del(k)+del(n));

end

if k>=3 && n>=3

if k==n

for N=[1 k+1:5]

t2=t2+abs(V(N)*Y(k,N))*sin(theta(k,N)-
del(k)+del(N));

end

t3=-2*abs(V(k)*Y(k,n))*sin(theta(k,n))-t2;
t2=0;
JQV(k-2,n-2)=t3;

elseif k~=n

JQV(k-2,n-2)=-abs(V(k)*Y(k,n))*sin(theta(k,n)+del(n)-
del(k));

end

final_jacobian=[JPD JPV;JQD JQV];
jac_inv=inv(final_jacobian)*final_vector';
count=2;

for kk=1:length(jac_inv)

if kk<=4

Del(kk+1)=jac_inv(kk);

elseif kk>3

count=count+1;
Delv(count)=jac_inv(kk);

end

del_tolerance=max(abs(Del-del));
del=Del+del;
V1=Delv+V;
f=find(V1>2);
f1=find(V1<1);
V1(f )=rand(size(f ));
V1(f1)=rand(size(f1));

%% Gauss seidel method

P2sch=(PsG(2)-PsL(2));
P3sch=(PsG(3)-PsL(3));
P4sch=(PsG(4)-PsL(4));
Q2sch=(QsG(2)-QsL(2));
Q3sch=(QsG(3)-QsL(3));
Q4sch=(QsG(4)-QsL(4));
V2=(1/Y(2,2))*((P2sch-1i*Q2sch)/V(2)-(Y(2,1)*V(1))-
(Y(2,3)*V(3))-(Y(2,4)*V(4))-(Y(2,5)*V(5)));
V3=(1/Y(3,3))*((P3sch-1i*Q3sch)/V(3)-(Y(3,1)*V(1))-
(Y(3,2)*V(2))-(Y(3,4)*V(4))-(Y(3,5)*V(5)));
V4=(1/Y(4,4))*((P4sch-1i*Q4sch)/V(4)-(Y(4,1)*V(1))-
(Y(4,2)*V(2))-(Y(4,3)*V(3))-(Y(4,5)*V(5)));

%% Updating reactive power on PV bus

International Journal of Scientific & Engineering Research, Volume 5, Issue 5, May-2014 1418

ISSN 2229-5518

Q=-
1i*[V(5)*(Y(5,1)*V(1)+Y(5,2)*V2+Y(5,3)*V3+Y(5,4)*V4+Y(5,5)*V(
5))];
V(2:4)=[V2 V3 V4];
V=abs(V);
k=1:5;
figure
plot(k,V1,'r','linewidth',2,'marker','o','markersize',10)
hold on
grid on
plot(k,V,'k','linewidth',2,'marker','o','markersize',10)
xlabel('Number of Buses','fontsize',10,'fontweight','bold')
ylabel('Voltage','fontsize',10,'fontweight','bold')
title('Graph between Voltage vs Number of bus-

es','fontsize',10,'fontweight','bold')

legend('Newton Raphson','Gauss seidel')

8 RESULT

8.1 Comparison between Real power and Reactive

Power.

Fig 5- Graph between Number of buses vs Power.

8.2 Comparison between Newton Raphson and Gauss- Seidel.

Fig 6- Graph between Voltage vs Number of buses.

9 CONCLUSION

Analyses, designing and comparison between different load flow system solving techniques i.e. Gauss-Seidel Method, Newton-Raphson Method in Power System using MATLAB has been successfully done and observed the desired result. In Gauss-Seidel, rate of convergence is slow. It can be easily pro- gram and the number of iterations increases directly with the number of buses in the system and in Newton-Raphson, the convergence is very fast and the number of iterations is inde- pendent of the size of the system, solution to a high accuracy is obtained. The NR Method convergence is not sensitive to the choice of slack bus. Although a large number of load flow methods are available in literature it has been observed that only the Newton-Raphson and the Fast Decoupled load-flow methods are most popular. The fast decoupled load flow is definitely superior to the Newton-Raphson Method from the point of view of speed and storage.

10 END SECTIONS

10.1 Acknowledgment

It is my proud privilege to express my sincere thanks to my supervisor, Mr. Sanjiv Kumar, Assistant Professor, Depart- ment of Electrical & Electronics Engineering, for valuable guidance and constant encouragement throughout the dura- tion of the dissertation.
My sincere thanks goes to my another guide Mr. Sajid Ali, Incharge Project & Seminar, EN, for his valuable suggestions and encouragement throughout the duration of my work and also giving an opportunity to work on MATLAB.
I am also thankful to my fellow classmates and friends for their cooperation and moral support during my stay.
Finally, my heartiest gratitude goes to my Parents, Brother,

International Journal of Scientific & Engineering Research, Volume 5, Issue 5, May-2014 1419

ISSN 2229-5518

Sister,..etc.. for their encouraging support as well as financial support throughout my studies which allowed me to concen- trate on my work. They have always been a great role model and lead me to greater success in pursuit of higher knowledge.

10.2 References

[1] John Grainger; William Stevenson's classic, Elements of Power System Analysis; Science/Engineering/Math; 1

ed.,pp.783., January 1, 1994.
[2] Hadi Saadat, McGraw Hill WCB;Power System Analysis,
[3] Ashfaq Husain;CBS Publishers & Distributors Pvt. Ltd.,5th
ed.,2007,pp.364-369.
[4] Sanjiv Kumar- Narendra Kumar;Kamal Jagasia publish-
er,Power System Analysis;1st ed.,2010,pp.150-202