Understanding Power Flow

Joydeep Mitra

1 Introduction
The problem of power flow, also known as load flow, concerns the determination of how much power (both
real and reactive) flows through each element of the system under a given loading scenario. In other words,
it answers the following question: given that at a certain point of time the generators are generating known
amounts of real and reactive power, and the loads are consuming known amounts of real and reactive power,
how much power flows through each component of the system? The answer to this question would influence
decisions concerning system operation, or planning, or both. Power flow analysis is perhaps the most widely
used of all power system analyses.
The determination of power flows in various circuit elements is basically the same as determination of
the currents flowing through these elements, and multiplying the complex conjugate of each current with
the corresponding voltage phasor; hence it appears that the solution of the matrix equation

I = YV, (1)

where I is the vector of current phasors in the circuit elements, V the vector of nodal (or bus) voltage
phasors, and Y the bus admittance matrix [1], should be a linear problem. However, because of the manner
in which power system variables are defined in practice, the problem ends up being non-linear. Generation
and load variables are described in terms of complex powers injected into the system, and these values do
not immediately yield the corresponding values of injected currents, because the voltages are unknown. This
is what makes the problem non-linear.
In explaining power flow analysis methods, therefore, it is often a good idea to spend a little time
explaining methods used in solving sets (or systems) of non-linear equations simultaneously. Hence this
article contains a section that describes, by means of examples, methods of solving non-linear equations.
The remainder of the article describes how the power flow problem is formulated so that these methods may
be used, and the application of these methods is also described by means of examples.

2 Developing the Power Flow Problem

2.1 System Representation
The system is divided into two parts: the active and the passive. The passive part consists of the network,
comprising passive elements such as transmission lines, transformers, inductors and capacitors. The active
part consists of generation and load, which are simply called injections. Generation is considered a positive
injection since it transfers complex power into the network; loads are considered negative injections because
they draw complex power out of the network. Conceptually, the system can be represented as shown in
figure 1. The net injection into any bus is the difference between the generation and the load at the bus, as
shown.
While the analysis of this system can be performed either by node-based or mesh-based approaches,
it is common practice to used nodal methods in solving power flows for the bulk power system, while

1
© Joydeep Mitra, 2010
 Passive
Network Sinj = Sgen Sload
Sgen

to
network Sload

Figure 1: Passive network and injected power.

mesh analysis is sometimes used in the solution of distribution power flows. In this article, therefore, nodal
analysis will be described. Further, balanced operation will be assumed, and only the per unit phase a to
neutral equivalent will be analyzed.
The network is accordingly modeled using the bus admittance matrix. As stated before, the injections
are described in terms of complex power, rather than current phasors.

2.2 Basic Relationships

The relationship shown in (1) is really I = Ybus V, which can be rewritten as
N
Im = ∑ YmnVn (2)
n=1

where N is the number of buses in the system. Using (2), the complex power injection into bus m can be
expressed as [ ]∗
N
Sm = Vm Im∗ = Vm ∑ YmnVn . (3)
n=1

Now if we express the voltage phasors and the admittances in polar form, i.e., in the form

Vn = |Vn |̸ δn
Ymn = |Ymn |̸ θmn

then (3) can be rewritten as

N
Pm + jQm = |Vm | ∑ |Ymn ||Vn |̸ (δm − δn − θmn ) (4)
n=1

which separates out into

N
Pm = |Vm | ∑ |Ymn ||Vn | cos(δm − δn − θmn ) (5)
n=1
N
Qm = |Vm | ∑ |Ymn ||Vn | sin(δm − δn − θmn ). (6)
n=1

These equations are used in the various solution techniques, as we shall see later in this article.

2
3 Solution of Sets of Nonlinear Equations
It is evident from (3–6) that regardless of which formulation we use, we would have to simultaneously solve
a set of nonlinear equations in order to obtain a solution of the power flow problem. In this section therefore
we shall try to understand, through simple examples, two methods that are frequently used in power flow
analysis—the Gauss-Seidel method and the Newton-Raphson Method.
For our examples, we shall use the system of equations given below.

5x12 − x1 x2 = 14 (7)
3x22 − 2x1 x2 = 15. (8)

It is clear that obtaining a closed form solution, i.e., solving these equations in one step, is difficult.
We therefore attempt to obtain solutions using iterative methods, where we begin by making reasonable
estimates of the solutions and successively computing better estimates until these estimates converge, i.e.,
attain stable values, to what are presumably the solutions that we are seeking.

3.1 Gauss-Seidel Method

This method of solving nonlinear equations consists of first rewriting the equations in such a manner that
each variable is separately expressed by an implicit relationship involving some or all variables, including
itself. As an illustration, (7) may be rewritten as
( )
1 14
x1 = + x2 .
5 x1

Then the best estimates of the variables are applied on the right hand side to yield an improved estimate of
the variable on the left hand side. In the beginning, the ”best estimate” is the initial estimate; thereafter, it
is the most recent estimate. This is performed for every variable. Each round of improvements constitutes
an iteration. Iterations are continued until the difference in successive estimates of every variable becomes
negligible.
For this example, the Gauss-Seidel formulation for iterative solution is of the form given below, where
k denotes the iteration count.
( )
(k+1) 1 14 (k)
x1 = + x2 (9)
5 x(k)
( 1 )
(k+1) 1 15 (k+1)
x2 = + 2x1 . (10)
3 x(k)
2

Iterations are continued until

(k) (k−1)
|xi − xi |≤ε ∀i (11)
where ε is a prespecified tolerance.
(0) (0)
To further illustrate the method, consider the initial estimates x1 = 1, x2 = 1. Then, in the first
iteration, we use these values to compute
( )
(1) 1 14 (0)
x1 = + x2 = 3,
5 x(0)
1

3
 (1) (0)
and then use x1 = 3, x2 = 1 to compute
( )
(1) 1 15 (1)
x2 = (0)
+ 2x1 = 7.
3 x2

If we continue this process of successive improvement, then the estimates converge as follows.

Iteration count (k) x1 x2

1 3.000 7.000
2 2.333 2.270
3 1.654 3.305
4 2.354 3.082
5 1.806 2.826
6 2.116 3.180
7 1.959 2.879
8 2.005 3.073
9 2.011 2.968
10 1.986 3.009
(10) (9) (10) (9)
At this point we find that |x1 − x1 | = 0.025 and |x2 − x1 | = 0.041 and note that if we had set our
tolerance at ε = 0.1 then the solution x1 = 1.986, x2 = 3.009 has converged to our stipulated accuracy1 . If
we desired a higher level of accuracy, we would have to set a lower tolerance and continue the iterations
further.

3.2 Newton-Raphson Method

Suppose that we have a nonlinear function of x, given by f (x), and that we want to solve the equation

f (x) = c. (12)

This function is shown graphically in figure 2. Solving (12) consists of determining the value of x⋆ , such
that
c − f (x⋆ ) = 0. (13)
Define
∆ f = c − f (x). (14)
The solution process consists of driving ∆ f to zero. The Newton-Raphson method accomplished this as
follows. Suppose that our initial estimate of x⋆ is x(0) . Then the improved estimate of x⋆ is calculated by
using a linear approximation of the curve around x(0) (which is essentially the tangent to the curve at that
point), and determining where the approximate curve (i.e., the tangent) intersects the line f (x) = c. (See
figure 2.) Then the slope of the tangent to the curve at x = x(0) is given by
( ) d f 
′
f x(0) =  .
dx x(0)

It is clear from figure 2 that the improved estimate x(1) is given by

x(1) = x(0) + ∆x(0) .

1 It is easy to confirm that x1 = 2, x2 = 3 is the exact solution to the system of equations given by (7–8).

4
 f(x)

∆f (0)

f (0) = f(x(0))
∆x(0)

tangent to f(x) at x = x(0)

x
x(0) x* x(1)

Figure 2: Graphical illustration of Newton-Raphson method.

It is also evident from figure 2 that

∆ f (0) ( )
′
= f x(0) .
∆x (0)

Hence in the Newton-Raphson method successive improvements of x(k) are achieved using

∆ f (k)
x(k+1) = x(k) + ∆x(k) = x(k) + ( ) (15)
f ′ x(k)

where k is the iteration count.

Now let us see how this method can be used to solve the system of equations given by (7–8). Define the
vectors
[ ] [ ] [ ] [ ]
x1 5x12 − x1 x2 c1 − f1 (x1 , x2 ) 14 − 5x12 + x1 x2
X= , F(X) = and ∆F = = . (16)
x2 3x22 − 2x1 x2 c2 − f2 (x1 , x2 ) 15 − 3x22 + 2x1 x2

When working with vectors, the first derivative of F(X) takes the form of a matrix, called the Jacobian,
given by [ ]
∂ f1 ∂ f1
∂x1 ∂x2
J= ∂ f2 ∂ f2
. (17)
∂x1 ∂x2
For this example, then, [ ]
10x1 − x2 −x1
J= . (18)
−2x2 6x2 − 2x1
For such a system of nonlinear equations, then, the equation for iterative solution using the Newton-
Raphson method is
[ ]−1
X (k+1) = X (k) + J (k) ∆F (k) . (19)

5
Iterations are continued until
(k)
|∆Fi | ≤ ε ∀ i (20)
where ε is a prespecified tolerance.
To illustrate how the Newton-Raphson method applies to this example, assume once again that we begin
(0) (0)
with the initial estimates x1 = 1, x2 = 1. Then, in the first iteration, we use these values in (16) and (18)
to compute [ ] [ ] [ ]
4 10 9 −1
(0)
F = ; ∆F =
(0)
; J = (0)
.
1 14 −2 4
Then [ ]−1 [ ][ ] [ ]
0.1176 0.0294 10 1.5882
∆X (0)
= J (0)
∆F (0)
= = ,
0.0588 0.2647 14 4.2941
and [ ] [ ] [ ]
1 1.5882 2.5882
X (1)
=X (0)
+ ∆X (0)
= + = .
1 4.2941 5.2941
If we continue this process of successive improvement, then the estimates converge as follows.

Iteration count (k) x1 x2 ∆ f1 ∆ f2

1 2.5882 5.2941 −5.7924 −41.6782
2 2.0846 3.5260 −0.3777 −7.5976
3 2.0050 3.0457 6.5468 × 10−3 −0.6156
4 2.0000 3.0004 1.0158 × 10−4 −5.7008 × 10−3

At this point we note that if we had set our tolerance at ε = 0.01 then the solution x1 = 2.0000, x2 =
3.0004 has converged to our stipulated accuracy. If we desired a higher level of accuracy, we would have to
set a lower tolerance and continue the iterations further.
It is relevant to remark that although the Newton-Raphson method involves more effort per iteration than
does the Gauss-Seidel method, the former converges in much fewer iterations than the latter.

4 Solving the Power Flow Problem

Having understood the Gauss-Seidel and Newton-Raphson methods, we shall now see how these methods
are applied to solve the power flow problem. Here, too, we shall adopt an illustrative approach.

4.1 Gauss-Seidel Power Flow

In this method, we rewrite (3) as [ ]
N
∗
Sm = Vm∗ Im = Vm∗ ∑ YmnVn , (21)
n=1

and express the voltage phasor at bus m as

 
 
1  
∗ N
 Sm 
Vm =  − ∑ YmnVn  . (22)
Ymm  Vm∗ 
 n=1 
n ̸= m

6
We then use the following iterative form
[ ]
(k)∗ m−1 N
1 Sm
− ∑ YmnVn − ∑ YmnVn
(k+1) (k+1) (k)
Vm = (23)
Ymm Vm(k)∗ n=1 n=m+1

where k is the iteration number. The following example illustrates the application.
Problem Statement: In the system shown in figure 3, the bus types, bus generation and load, and the line
impedances are as given in the figure. (a) Build the bus admittance matrix. (b) Set up the load flow equations
for solution using the Gauss-Seidel method. (c) Beginning with flat-start estimates for the voltages at buses
2, 3 and 4, determine the power flow solution using the Gauss-Seidel method. Use 100 MVA as system base.

SG2 = 40+j30 MVA

V1 = 1.05 0 pu
Bus 1: slack j0.1 pu Bus 2: PQ

j0.25 pu SL2 = 10+j5 MVA

j0.2 pu
j0.2 pu
Bus 3: PQ j0.25 pu Bus 4: PQ
SL3 = 30+j12 MVA SL4 = 20+j15 MVA

Figure 3: Four-bus system.

Solution: We observe from figure 3 that V1 is completely specified, and we need to find only V2 , V3 and V4
in order to be able to determine all line flows. We also note that bus 1 is labeled “slack” and buses 2, 3 and
4 are all labeled “PQ”. A “PQ bus” is a bus at which both real and reactive power are specified. Now in
the system shown in figure 3, the total complex power generated at buses 1 and 2 must equal the complex
power loads at buses 2, 3 and 4, plus the complex power losses in the five transmission lines. However,
the losses remain unknown until the power flow has been solved. So while the complex power generated
at bus 2 is known, the generator at bus 1 must supply the difference in generated power and loads plus the
unknown losses. Hence bus 1 is known as the “slack bus”, and its complex power injection is unknown.
Having realized this, we proceed to solve the power flow problem as follows.
(a) First construct the Ybus :

1 2 3 4
j19 j10 j5 j4 1
j10 j15 0 j5 2
Ybus =
j5 0 j9 j4 3
j4 j5 j4 j13 4

(b) Gauss-Seidel power flow equations:

[ ]
(k+1) 1 S2∗ (k)
V2 = −Y21V1 −Y24V4
Y22 V (k)∗
2

7
 [ ]
(k+1) 1 S3∗ (k)
V3 = −Y31V1 −Y34V4
Y33 V (k)∗
[ 3 ]
(k+1) 1 S4∗ (k+1) (k+1)
V4 = −Y41V1 −Y42V2 −Y43V3
Y44 V (k)∗
4

where V1 = 1.05̸ 0 pu, S2 = 0.3 + j0.25 pu, S3 = −0.3 − j0.12 pu, S4 = −0.20 − j0.15 pu and k is the
iteration number.
(c) The first Gauss-Seidel iteration:
(0) (0) (0)
Flat start ⇒ V2 = V3 = V4 = 1.0̸ 0 pu.
[ ]
1 0.3 − j0.25
− j10 × 1.05 − j5 × 1.0 = 1.0502̸ 1.091◦ pu
(1)
V2 =
− j15 1.0
[ ]
1 −0.3 + j0.12
− j5 × 1.05 − j4 × 1.0 = 1.0150̸ − 1.882◦ pu
(1)
V3 =
− j9 1.0
[ ]
(1) 1 −0.2 + j0.15 ◦ ◦
V4 = − j4 × 1.05 − j5 × 1.0502 1.091 − j4 × 1.0150 − 1.882
̸ ̸
− j13 1.0
= 1.0277̸ − 1.001◦ pu.

The second Gauss-Seidel iteration:

[ ]
1 0.3 − j0.25
− j10 × 1.05 − j5 × 1.0277 − 1.001 = 1.0581̸ 0.724◦ pu
◦
(2) ̸
V2 =
− j15 1.0502̸ − 1.091◦
[ ]
1 −0.3 + j0.12
− j5 × 1.05 − j4 × 1.0277̸ − 1.001 = 1.0266̸ − 2.254◦ pu
◦
(2)
V3 =
− j9 1.0150̸ 1.882◦
[ ]
(2) 1 −0.2 + j0.15 ◦ ◦
V4 = − j4 × 1.05 − j5 × 1.0581 0.724 − j4 × 1.0266 − 2.254
̸ ̸
− j13 1.0277̸ 1.001◦
= 1.0344̸ − 1.222◦ pu.

Convergence trend:

Iteration count V2 (pu) V3 (pu) V4 (pu)

1 1.0502̸ 1.091◦ 1.0150̸ − 1.882◦ 1.0277̸ − 1.001◦
2 1.0581̸ 0.724◦ 1.0266̸ − 2.254◦ 1.0344̸ − 1.222◦
3 1.0603̸ 0.635◦ 1.0295̸ − 2.323◦ 1.0362̸ − 1.270◦
4 1.0609̸ 0.615◦ 1.0303̸ − 2.338◦ 1.0366̸ − 1.280◦
5 1.0610̸ 0.610◦ 1.0306̸ − 2.341◦ 1.0368̸ − 1.282◦
6 1.0611̸ 0.609◦ 1.0306̸ − 2.341◦ 1.0368̸ − 1.282◦
7 1.0611̸ 0.609◦ 1.0306̸ − 2.341◦ 1.0368̸ − 1.282◦

After 7 iterations, the solution is said to have converged to a voltage tolerance of 0.0001 pu, i.e.,
(k) (k−1)
|Vi −Vi | ≤ 0.0001 pu ∀ i.

At this point, if the accuracy is deemed acceptable, we can stop the iterations and use the values obtained at
the 7th iteration to determine the following.

S1 = V1 I1∗ = V1 (Y11V1 +Y12V2 +Y13V3 +Y14V4 )∗ = 0.2000 + j0.0468 pu.

8
Notice that this really means that the generator at the slack bus must generate this amount of complex power,
i.e., 20 MW active power and 4.68 MVAr reactive power. If we tally the generated power and the load power,
we get:
Total generation = G1 + G2 = 20 + j4.68 + 40 + j30 = 60 + j34.68 MVA
and total load = L2 + L3 + L4 = 10 + j5 + 30 + j12 + 20 + j15 = 60 + j32 MVA.
This implies that system losses equal 0 + j2.68 MVA, which means that the reactive power losses in the lines
add up to 2.68 MVAr, while there is no active power loss. This is not surprising, since there is no resistance
in the given system, so there is no I 2 R loss. The power flow solution can be completed by determining the
line flows from the line data and the voltage solutions.
Inclusion of voltage controlled bus: The above example had a slack bus and PQ buses. Now we introduce a
third type of bus, the “PV bus” or voltage controlled bus. This is a bus where the real power injection and
the voltage magnitude are specified, and the reactive power injection can be controlled to hold the voltage
magnitude constant at the specified value. This can normally be achieved in a synchronous generator by
adjusting its field current.
It is appropriate at this point to point out the following. We have seen three bus types, with the following
properties.
bus type known variables unknown variables
slack voltage magnitude and angle real and reactive power injection
PQ real and reactive power injection voltage magnitude and angle
PV real power injection and voltage magnitude reactive power injection and voltage angle
Let us consider the case where bus 2 is a PV bus, with the following specification: PG2 = 40 MW and
|V2 | = 1.06 pu. The new system is shown in figure 4. The new power flow solution is determined as follows.

PG2 = 40 MW
V1 = 1.05 0 pu
|V2| = 1.06 pu
Bus 1: slack j0.1 pu Bus 2: PV
j0.25 pu SL2 = 10+j5 MVA
j0.2 pu
j0.2 pu
Bus 3: PQ j0.25 pu Bus 4: PQ
SL3 = 30+j12 MVA SL4 = 20+j15 MVA

Figure 4: Four-bus system including voltage-controlled bus.

(a) The bus admittance matrix remains unchanged.

(b) The Gauss-Seidel power flow equations are solved in the following manner. First determine updated
values of Q2 and V2 .
[ ( ) ]
(k+1) (k) (k) (k) ∗
Q2 = Im V2 Y21V1 +Y22V2 +Y24V4
[ (k+1)
]
(k+1)′ 1 P2 − jQ2 (k)
V2 = (k)∗
−Y21V1 −Y24V4
Y22 V2
( )
(k+1) (k+1)′ 1.06
V2 = V2 (k+1)′
.
|V2 |

9
Then update the remaining voltages.
[ ]
(k+1) 1 S3∗ (k)
V3 = −Y31V1 −Y34V4
Y33 V (k)∗
[ 3
]
(k+1) 1 S4∗ (k+1) (k+1)
V4 = −Y41V1 −Y42V2 −Y43V3
Y44 V (k)∗
4

where V1 = 1.05̸ 0 pu, |V2 | = 1.06 pu, P2 = 0.3 pu, S3 = −0.3 − j0.12 pu, S4 = −0.20 − j0.15 pu and k is
the iteration number.
(c) The first Gauss-Seidel iteration:
(0) (0) (0)
V2 = 1.06̸ 0 pu and V3 = V4 = 1.0̸ 0 pu.
[ ]
= Im 1.06 ( j10 × 1.05 − j15 × 1.06 + j5 × 1.0)∗ = 0.4240 pu
(1)
Q2
[ ]
(1)′ 1 0.3 − j0.4240
V2 = − j10 × 1.05 − j5 × 1.0 = 1.0602̸ 1.020◦ pu
− j15 1.06
= 1.06̸ 1.020◦ pu
(1)
V2
[ ]
1 −0.3 + j0.12
− j5 × 1.05 − j4 × 1.0 = 1.0150̸ − 1.882◦ pu
(1)
V3 =
− j9 1.0
[ ]
(1) 1 −0.2 + j0.15 ◦ ◦
V4 = − j4 × 1.05 − j5 × 1.06 1.020 − j4 × 1.0150 − 1.882
̸ ̸
− j13 1.0
= 1.0315̸ − 1.021◦ pu.

The second Gauss-Seidel iteration:

[ ]
= Im 1.06̸ 1.020 ( j10 × 1.05 − j15 × 1.06̸ 1.020 + j5 × 1.0315̸ − 1.021◦ )∗ = 0.2625 pu
(2)
Q2
[ ]
(2)′ 1 0.3 − j0.2625
V2 = − j10 × 1.05 − j5 × 1.0315 − 1.021 = 1.0600̸ 0.704◦ pu
̸ ◦
− j15 1.06̸ − 1.020
= 1.06̸ 0.704◦ pu
(2)
V2
[ ]
1 −0.3 + j0.12
− j5 × 1.05 − j4 × 1.0315̸ − 1.021 = 1.0283̸ − 2.261◦ pu
◦
(2)
V3 =
− j9 1.0150̸ 1.882◦
[ ]
(2) 1 −0.2 + j0.15 ◦ ◦
V4 = − j4 × 1.05 − j5 × 1.06 0.704 − j4 × 1.0283 − 2.261
̸ ̸
− j13 1.0315̸ 1.021◦
= 1.0357̸ − 1.227◦ pu.

Convergence trend:

Iteration count Q2 (pu) V2 (pu) V3 (pu) V4 (pu)

1 0.4240 1.06̸ 1.020◦ 1.0150̸ − 1.882◦ 1.0315̸ − 1.021◦
2 0.2625 1.06̸ 0.704◦ 1.0283̸ − 2.261◦ 1.0357̸ − 1.227◦
3 0.2389 1.06̸ 0.630◦ 1.0301̸ − 2.322◦ 1.0362̸ − 1.270◦
4 0.2356 1.06̸ 0.615◦ 1.0304̸ − 2.337◦ 1.0363̸ − 1.280◦
5 0.2351 1.06̸ 0.611◦ 1.0304̸ − 2.341◦ 1.0363̸ − 1.282◦
6 0.2351 1.06̸ 0.611◦ 1.0304̸ − 2.342◦ 1.0363̸ − 1.283◦
7 0.2351 1.06̸ 0.610◦ 1.0304̸ − 2.342◦ 1.0363̸ − 1.283◦

10
At this point, if the accuracy is deemed acceptable, we can stop the iterations and use the values obtained at
the 7th iteration to determine the following.

S1 = V1 I1∗ = V1 (Y11V1 +Y12V2 +Y13V3 +Y14V4 )∗ = 0.2000 + j0.0615 pu

S2 = P2 + jQ2 = 0.3000 + j0.2351 pu.

Notice that this really means that the generator at the slack bus must generate 20 + j6.15 MVA and the
generator at bus 2 must generate 40 + j28.51 MVA, so that the injection at bus 2 can be 30 + j23.51 MVA.
If we tally the generated power and the load power, we get:
Total generation = G1 + G2 = 20 + j6.15 + 40 + j28.51 = 60 + j34.66 MVA
and total load = L2 + L3 + L4 = 10 + j5 + 30 + j12 + 20 + j15 = 60 + j32 MVA.
This implies that system losses equal 0 + j2.66 MVA, which means that the reactive power losses in the lines
add up to 2.66 MVAr, which is not very different from before. The differences are explained as follows. In
the second case, the voltage at bus 2 is lower than in the first case, so the reactive power generated at bus 2
is lower. Bus 1 had to generate a little more reactive power than in the first case, so as to meet the reactive
power needs of the system.
The power flow solution can be completed by determining the line flows from the line data and the
voltage solutions:
( )
V1 −V2 ∗
S12 = V1 = −0.1185 − j0.1044 pu = −11.85 − j10.44 MVA.
z12

Similarly, S13 = 22.11 + j10.73 MVA, S14 = 9.75 + j5.85 MVA, S24 = 18.15 + j12.85 MVA and S34 =
−7.89 − j2.36 MVA.

4.2 Newton-Raphson Power Flow

In applying the Newton-Raphson method to the power flow problem, we separate out the variables P, Q,
|V | and θ, and usually formulate the problem entirely in terms of real numbers. While there are several
variants of the Newton-Raphson power flow formulation, we shall discuss the formulation described in the
text book [1], which expresses the real and reactive power in the form shown in equations (5) and (6), with
the voltages and admittances in polar form, i.e.,

S = P + jQ, V = |V |̸ δ, Y = |Y |̸ δ.

In commonly used implementations of Newton-Raphson power flow, all complex variables are separated
out into real and imaginary parts or into magnitude and angle, so that the entire solution is performed using
real variables. So for the above representation we would use
[ ] [ ]
[δ] [P(δ, |V |)]
X= and F(X) = ,
[|V |] [Q(δ, |V |)]

so the Jacobian is built in four blocks:

 [ ] [ ]
∂P ∂P
[ ] 
J1 J2  ∂δ ∂|V |  
J= = [ ] [ ] . (24)
J3 J4  ∂Q ∂Q 
∂δ ∂|V |

11
The entries in these blocks are described below in the context of the example, but table 6.5 in the text book
[1] shows the general forms of these terms.
The method is illustrated by solving the same problem as before, including the voltage-controlled bus.
(a) The bus admittance matrix remains unchanged.
(b) Newton-Raphson power flow equations:
Among the bus voltages, the following values are unknown: δ2 , V3 and V4 ; hence the vector of unknowns
is given by:  
δ2
 δ3 
 
X =  δ 4
.
 (25)
 |V3 | 
|V4 |
The following injections are known: P2 = 0.3 pu, S3 = −0.3 − j0.12 pu, and S4 = −0.20 − j0.15 pu. Hence
the following equality must hold:
   
P2 (X) 0.30
 P3 (X)   −0.30 
   
 P4 (X)  =  −0.20  pu (26)
   
 Q3 (X)   −0.12 
Q4 (X) −0.15

where

P2 (X) = |V2 |2 |Y22 | cos θ22 + |V2 | {|Y21 ||V1 | cos(δ2 − δ1 − θ21 ) + |Y24 ||V4 | cos(δ2 − δ4 − θ24 )}
P3 (X) = |V3 |2 |Y33 | cos θ33 + |V3 | {|Y31 ||V1 | cos(δ3 − δ1 − θ31 ) + |Y34 ||V4 | cos(δ3 − δ4 − θ34 )}
P4 (X) = |V4 |2 |Y44 | cos θ44 + |V4 | {|Y41 ||V1 | cos(δ4 − δ1 − θ41 ) + |Y42 ||V2 | cos(δ4 − δ2 − θ42 )
+|Y43 ||V3 | cos(δ4 − δ3 − θ43 )}
Q3 (X) = −|V3 |2 |Y33 | sin θ33 + |V3 | {|Y31 ||V1 | sin(δ3 − δ1 − θ31 ) + |Y34 ||V4 | sin(δ3 − δ4 − θ34 )}
Q4 (X) = −|V4 |2 |Y44 | sin θ44 + |V4 | {|Y41 ||V1 | sin(δ4 − δ1 − θ41 ) + |Y42 ||V2 | sin(δ4 − δ2 − θ42 )
+|Y43 ||V3 | sin(δ4 − δ3 − θ43 )} .

Representing the left hand side of (26) as F(X), we can represent the deviations in F as
 
0.30 − P2 (X)
 −0.30 − P3 (X) 
 
∆F =  −0.20 − P 4 (X) .
 (27)
 −0.12 − Q3 (X) 
−0.15 − Q4 (X)

Then the equation representing the Newton-Raphson iterations would be given by

[ ]−1
X (k+1) = X (k) + J (k) ∆F (k) (28)

where k represents the iteration number and the elements of the Jacobian matrix [J] are given by

∂P2
J11 = = −|V2 | {|Y21 ||V1 | sin(δ2 − δ1 − θ21 ) + |Y24 ||V4 | sin(δ2 − δ4 − θ24 )}
∂δ2

12
 ∂P2
J12 = =0
∂δ3
∂P2
J13 = = |V2 ||Y24 ||V4 | sin(δ2 − δ4 − θ24 )
∂δ4
∂P2
J14 = =0
∂|V3 |
∂P2
J15 = = |V2 ||Y24 | cos(δ2 − δ4 − θ24 )
∂|V4 |
∂P3
J21 = =0
∂δ2
∂P3
J22 = = −|V3 | {|Y31 ||V1 | sin(δ3 − δ1 − θ31 ) + |Y34 ||V4 | sin(δ3 − δ4 − θ34 )}
∂δ3
∂P3
J23 = = |V3 ||Y34 ||V4 | sin(δ3 − δ4 − θ34 )
∂δ4
∂P3
J24 = = 2|V3 ||Y33 | cos θ33 + |Y31 ||V1 | cos(δ3 − δ1 − θ31 ) + |Y34 ||V4 | cos(δ3 − δ4 − θ34 )
∂|V3 |
∂P3
J25 = = |V3 ||Y34 | cos(δ3 − δ4 − θ34 )
∂|V4 |
∂P4
J31 = = |V4 ||Y42 ||V2 | sin(δ4 − δ2 − θ42 )
∂δ2
∂P4
J32 = = |V4 ||Y43 ||V3 | sin(δ4 − δ3 − θ43 )
∂δ3
∂P4
J33 = = −|V4 | {|Y41 ||V1 | sin(δ4 − δ1 − θ41 ) + |Y42 ||V2 | sin(δ4 − δ2 − θ42 ) + |Y43 ||V3 | sin(δ4 − δ3 − θ43 )}
∂δ4
∂P4
J34 = = |V4 ||Y43 | cos(δ4 − δ3 − θ43 )
∂|V3 |
∂P4
J35 = = 2|V4 ||Y44 | cos θ44 + |Y41 ||V1 | cos(δ4 − δ1 − θ41 ) + |Y42 ||V2 | cos(δ4 − δ2 − θ42 )
∂|V4 |
+ |Y43 ||V3 | cos(δ4 − δ3 − θ43 )
∂Q3
J41 = =0
∂δ2
∂Q3
J42 = = |V3 | {|Y31 ||V1 | cos(δ3 − δ1 − θ31 ) + |Y34 ||V4 | cos(δ3 − δ4 − θ34 )}
∂δ3
∂Q3
J43 = = −|V3 ||Y34 ||V4 | cos(δ3 − δ4 − θ34 )
∂δ4
∂Q3
J44 = = −2|V3 ||Y33 | sin θ33 + |Y31 ||V1 | sin(δ3 − δ1 − θ31 ) + |Y34 ||V4 | sin(δ3 − δ4 − θ34 )
∂|V3 |
∂Q3
J45 = = |V3 ||Y34 | sin(δ3 − δ4 − θ34 )
∂|V4 |
∂Q4
J51 = = −|V4 ||Y42 ||V2 | cos(δ4 − δ2 − θ42 )
∂δ2
∂Q4
J52 = = −|V4 ||Y43 ||V3 | cos(δ4 − δ3 − θ43 )
∂δ3
∂Q4
J53 = = |V4 | {|Y41 ||V1 | cos(δ4 − δ1 − θ41 ) + |Y42 ||V2 | cos(δ4 − δ2 − θ42 ) + |Y43 ||V3 | cos(δ4 − δ3 − θ43 )}
∂δ4

13
 ∂Q4
J54 = = |V4 ||Y43 | sin(δ4 − δ3 − θ43 )
∂|V3 |
∂Q4
J55 = = −2|V4 ||Y44 | sin θ44 + |Y41 ||V1 | sin(δ4 − δ1 − θ41 ) + |Y42 ||V2 | sin(δ4 − δ2 − θ42 )
∂|V4 |
+ |Y43 ||V3 | sin(δ4 − δ3 − θ43 ).
(c) Solution using Newton-Raphson method:
The first Newton-Raphson iteration: Using V1 = 1.05̸ 0 pu, |V2 | = 1.06 pu, and flat start estimates, we have:
     
0.0 0.3 16.43 0 −5.3 0 0
 0.0   −0.3   0 9.25 −4 0 0 
     
X =
(0)    
 0.0  pu, ∆F =  −0.2  pu, and J =  −5.3 −4 13.5
(0) (0)
0 0  pu.
 1.0   0.13   0 0 0 8.75 −4 
1.0 0.35 0 0 0 −4 12.5
This gives us
 
0.0108
[ ]−1  −0.0424 
 
∆X (1) = J (0)
∆F = 
(0)
 −0.0232  pu

 0.0324 
0.0384
 (1) 
δ2    
 (1)  0.0108 0.6182◦
 δ3   −0.0424   −2.4318◦ 
 (1)     
X (1)
= X (0) + ∆X (1) =   
 δ4  =  −0.0232
 pu = 
  −1.3267◦ .

 (1)   1.0324   1.0324 pu 
 |V3 | 
(1) 1.0384 1.0384 pu
|V4 |
The second Newton-Raphson iteration: Using values from the first iteration, the second iteration yields:
   
−6.8463 × 10−3 16.6295 0 −5.5002 0 0.1799
 0.0127   0 9.7204 −4.2872 −0.3029 −0.0797 
   
∆F = 
(1) −3  
 5.0380 × 10  pu, and J =  −5.5002 −4.2872 14.1474
(1)
0.0801 −0.1975   pu.
 −0.0102   0 −0.3127 0.0827 9.1852 −4.1288 
−0.0193 0.1868 −0.0827 −0.2050 −4.1527 13.3729
This gives us
 
  (2)
δ2  
−1.3967 × 10−4   0.6102◦
 1.5663 × 10−3   (2)
δ3   −2.3422◦ 
     
∆X = 
(2)
 7.5925 × 10−4  pu, and X = 

(2)

(2)
δ4 =
  −1.2832◦ .

 −1.9726 × 10−3   (2)   1.0304 pu 
 |V3 | 
−2.0326 × 10−3 (2) 1.0363 pu
|V4 |
The third Newton-Raphson iteration: Using values from the second iteration, the third iteration yields:
   
−9.6026 × 10−6 16.6189 0 −5.4896 0 0.1751
 2.9582 × 10−5   0 9.6759 −4.2707 −0.2912 −0.0762 
   
∆F (2) =  −6  
 2.4013 × 10  pu, and J =  −5.4896 −4.2707 14.1118
(2)
0.0766 −0.1930  pu.

 −2.7952 × 10−5   0 −0.3000 0.0789 9.1574 −4.1210 
−4.3274 × 10−5 0.1815 −0.0789 −0.2000 −4.1446 13.3276

14
This gives us
 
  (2)
δ2  
−1.7164 × 10−7   0.6102◦
 3.3408 × 10−6   (2)
δ3   −2.3420◦ 
     
∆X = 
(3)
 1.0766 × 10−6  pu, and X = 

(3)

(2)
δ4 =
  −1.2831◦ .

 −5.1115 × 10−6   (2)   1.0304 pu 
 |V3 | 
−4.7983 × 10−6 (2) 1.0363 pu
|V4 |

At this point we can see that the solution has stabilized, since X (3) and X (2) are almost identical, and
(k) (k−1)
|Vi −Vi | ≤ 0.0001 pu ∀ i,

so we can stop the iterations. We could even have stopped after the second iteration, considering how small
∆F (2) is, meaning the condition given by (49) is practically satisfied at the end of the second iteration. In
fact, the latter criterion is more commonly used in Newton-Raphson methods.
Now we have the voltage solution
   
V1 1.05̸ 0
 V2   1.06̸ 0.610◦ 
   
 V3  =  1.0304̸ − 2.342◦  pu
V4 1.0363̸ − 1.283◦
from which we can calculate the unknown real and reactive power injections:

P1 = |V1 |2 |Y11 | cos θ11 + |V1 | {|Y12 ||V2 | cos(δ1 − δ2 − θ12 ) + |Y13 ||V3 | cos(δ1 − δ3 − θ13 )
+|Y14 ||V4 | cos(δ1 − δ4 − θ14 )} = 0.2000 pu
Q1 = −|V1 |2 |Y11 | sin θ11 + |V1 | {|Y12 ||V2 | sin(δ1 − δ2 − θ12 ) + |Y13 ||V3 | sin(δ1 − δ3 − θ13 )
+|Y14 ||V4 | sin(δ1 − δ4 − θ14 )} = 0.0615 pu
Q2 = −|V2 |2 |Y22 | sin θ22 + |V2 | {|Y21 ||V1 | sin(δ2 − δ1 − θ21 ) + |Y24 ||V4 | sin(δ2 − δ4 − θ24 )} = 0.2351 pu.

These are the same results as those we obtained using the Gauss-Seidel method. The power flow solution
can be completed by determining the line flows.

4.3 Fast Decoupled Power Flow

In the Newton-Raphson power flow solution above, observe that the off-diagonal blocks of the Jacobian,
i.e., the J2 and J3 blocks in 24, are small in magnitude, compared to the J1 and J4 blocks. In other words, the
coupling between real power and voltage magnitudes is weak, as is that between reactive power and voltage
angles. The method of fast decoupled power flow takes advantage of this, and neglects the J2 and J3 blocks
in the Jacobian. It also uses some approximations that enable easy computation of the J1 and J4 blocks, and
does not update them over the iterations. Hence considerable computational effort is avoided by precluding
the computation and factorization of the Jacobian every iteration.
The method of fast decoupled power flow can be understood as follows. First, the equation
[ ] [ ][ ]
∆P(δ, |V |) J1 J2 ∆δ
= (29)
∆Q(δ, |V |) J3 J4 ∆|V |
is simplified to the form [ ] [ ][ ]
∆P(δ) J1 0 ∆δ
≈ (30)
∆Q(|V |) 0 J4 ∆|V |

15
which reduces to

[∆P] ≈ [J1 ] [∆δ] (31)

[∆Q] ≈ [J4 ] [∆|V |] . (32)

Then the terms of J1 and J4 are simplified as follows. Consider (see table 6.5 in the text book [1])
∂Pi
= |Vi ||Yi j ||V j | sin(δi − δ j − θi j ) (33)
∂δ j
∂Pi N
= −|Vi |2 |Yii | sin θii − |Vi | ∑ |Yi j ||V j | sin(δi − δ j − θi j ) = −|Vi |2 |Yii | sin θii − Qi (34)
∂δi j=1
∂Qi
|V j | = |Vi ||Yi j ||V j | sin(δi − δ j − θi j ) (35)
∂|V j |
∂Qi N
|Vi | = −|Vi |2 |Yii | sin θii + |Vi | ∑ |Yi j ||V j | sin(δi − δ j − θi j ) = −|Vi |2 |Yii | sin θii + Qi . (36)
∂|Vi | j=1

To understand how we introduced Qi in (34) and (36) above, see equation (6.6.3) in the text [1]. Assuming
that δi ≈ δ j ∀ i, j and Qi ≪ |Vi |2 |Yii | sin θii ∀ i, these equations can be simplified to

∂Pi
≈ −|Vi ||Yi j ||V j | sin θi j = −|Vi ||V j |Bi j (37)
∂δ j
∂Pi
≈ −|Vi |2 |Yii | sin θii = −|Vi |2 Bii (38)
∂δi
∂Qi
|V j | ≈ −|Vi ||Yi j ||V j | sin θi j = −|Vi ||V j |Bi j (39)
∂|V j |
∂Qi
|Vi | ≈ −|Vi |2 |Yii | sin θii = −|Vi |2 Bii (40)
∂|Vi |
where Bi j is the imaginary part of Yi j .
Now we can express the terms of (31) and (32) as
N
∂Pi N
∆Pi N
∆Pi = ∑ ∂δ j ∆δ j ≈ −|Vi | ∑ |V j |Bi j ∆δ j ⇒ |Vi | ≈ − ∑ |V j |Bi j ∆δ j (41)
j=1 j=1 j=1

and
N
∂Qi N
∆Qi N
∆Qi = ∑ ∂|V j | ∆|V j | ≈ −|Vi | ∑ Bi j ∆|V j | ⇒ |Vi |
≈ − ∑ Bi j ∆|V j |. (42)
j=1 j=1 j=1

In (41), a further approximation is made: |V j | ≈ 1.0 pu ∀ j, so that

∆Pi N
≈ − ∑ Bi j ∆δ j . (43)
|Vi | j=1

Now (43) and (42) can be written in matrix form as

[ ]
∆P [ ]
≈ B′ [∆δ] (44)
|V |
[ ]
∆Q [ ]
≈ B′′ [∆|V |] (45)
|V |

16
where the B′ and B′′ matrices are composed of the negative values of the imaginary parts of the corresponding
terms from the bus admittance matrix. The difference between B′ and B′′ is that while the rows and columns
of the former correspond to all buses except the slack, those of the latter correspond to PQ buses only. The
form shown in (44) and (45) represents the approximations used in fast decoupled power flow as originally
proposed by Stott and Alsac in [2]. This reference also spells out some other approximations that will not
be repeated here.
The reason the values of |V j | are approximated to 1.0 pu in (41), but this approximation is not extended
to the voltage magnitudes in the denominators of (44) and (45), is that the former affect the coupling, i.e.,
the effect of the reactive power flows on ∆δ while the latter affect the behavior of the defining functions
∆P(δ, |V |) and ∆Q(δ, |V |), as explained in [2].
The text book [1], however, states that all voltage magnitudes may be approximated to 1.0 pu, and some
references accept this, while others use the form proposed by Stott and Alsac [2]. If we were to accept
the assertion that all voltage magnitudes may be approximated to 1.0 pu, then (31) and (32) may simply be
expressed as
[ ]
[∆P] ≈ [J1 ] [∆δ] ≈ B′ [∆δ] (46)
[ ′′ ]
[∆Q] ≈ [J4 ] [∆|V |] ≈ B [∆|V |] . (47)
The application of (46) and (47) in a fast decoupled power flow solution is illustrated by the same
problem as before, including the voltage control bus.
(a) The bus admittance matrix remains unchanged.
(b) Fast decoupled power flow equations:
In solving this problem, we will represent the complex power in rectangular coordinates, and the voltages
and admittances in polar coordinates (as in the text book [1]), i.e., in the form
S = P + jQ, V = |V |̸ δ, Y = |Y |̸ θ.
Among the bus voltages, the following values are unknown: δ2 , V3 and V4 ; hence the vectors of un-
knowns are given by:  
δ2 [ ]
  |V3 |
δ = δ3 and |V | = . (48)
|V4 |
δ4
The following injections are known: P2 = 0.3 pu, S3 = −0.3 − j0.12 pu, and S4 = −0.20 − j0.15 pu. Hence
the following equality must hold:
   
P2 (δ, |V |) 0.30
 P3 (δ, |V |)   −0.30 
   
 P4 (δ, |V |)  =  −0.20  pu. (49)
   
 Q3 (δ, |V |)   −0.12 
Q4 (δ, |V |) −0.15
where
P2 (δ, |V |) = |V2 |2 |Y22 | cos θ22 + |V2 | {|Y21 ||V1 | cos(δ2 − δ1 − θ21 ) + |Y24 ||V4 | cos(δ2 − δ4 − θ24 )}
P3 (δ, |V |) = |V3 |2 |Y33 | cos θ33 + |V3 | {|Y31 ||V1 | cos(δ3 − δ1 − θ31 ) + |Y34 ||V4 | cos(δ3 − δ4 − θ34 )}
P4 (δ, |V |) = |V4 |2 |Y44 | cos θ44 + |V4 | {|Y41 ||V1 | cos(δ4 − δ1 − θ41 ) + |Y42 ||V2 | cos(δ4 − δ2 − θ42 )
+|Y43 ||V3 | cos(δ4 − δ3 − θ43 )}
Q3 (δ, |V |) = −|V3 |2 |Y33 | sin θ33 + |V3 | {|Y31 ||V1 | sin(δ3 − δ1 − θ31 ) + |Y34 ||V4 | sin(δ3 − δ4 − θ34 )}
Q4 (δ, |V |) = −|V4 |2 |Y44 | sin θ44 + |V4 | {|Y41 ||V1 | sin(δ4 − δ1 − θ41 ) + |Y42 ||V2 | sin(δ4 − δ2 − θ42 )
+|Y43 ||V3 | sin(δ4 − δ3 − θ43 )} .

17
Now we represent ∆P and ∆Q as
 
0.30 − P2 (δ, |V |) [ ]
  −0.12 − Q3 (δ, |V |)
∆P = −0.30 − P3 (δ, |V |) and ∆Q = . (50)
−0.15 − Q4 (δ, |V |)
−0.20 − P4 (δ, |V |)

Having defined all the variables, we can now extract B′ and B′′ directly from the bus admittance matrix:
 
15 0 −5 [ ]
′   ′′ 9 −4
B = 0 9 −4 and B = . (51)
−4 13
−5 −4 13

Then the equations representing the fast decoupled iterations would be given by
[ ]−1 (k)
δ(k+1) = δ(k) + B′ ∆P (52)
[ ]−1
|V |(k+1) = |V |(k) + B′′ ∆Q(k+0.5) (53)

where k represents the iteration number. The value of ∆Q in (53) is calculated from |V |(k) from the previous
iteration and δ(k+1) in the present iteration; hence its iteration count is incremented by 0.5.
(c) Solution using fast decoupled method:
The first iteration: Using V1 = 1.05̸ 0 pu, |V2 | = 1.06 pu, and flat start estimates, we have:
 (0)   
δ2 0.0 [ ] [ ]
(0)
 (0)   |V | 1.0
δ =  δ3  = 0.0  pu, and |V | =
(0) (0) 3
(0) = pu.
(0) 0.0 |V | 1.0
δ4 4

This gives us
   
0.3 [ ]−1 (0) 0.0119
∆P(0) =  −0.3  pu, ∆δ(1) = B′ ∆P =  −0.0442  pu
−0.2 −0.0244
 (1)     
δ2 0.0119 0.6796◦
  
δ(1) = δ(0) + ∆δ(1) =  δ(1)
3 
= −0.0442  pu =  −2.5317◦ 
(1)
δ4 −0.0244 −1.3991◦
[ ] [ ]
0.1241 [ ′′ ]−1 (0.5) 0.0296
∆Q(0.5) = pu, ∆|V | = B
(1)
∆Q = pu
0.3445 0.0356
[ ] [ ]
(1)
|V | 1.0296
|V |(1) = |V | + ∆|V | =
(0) (1) 3
(1) = pu.
|V4 | 1.0356

Notice how within the iteration the updated value of δ(1) is used to determine the value of ∆Q; hence ∆Q
has been assigned the iteration number 0.5. This results in better values of ∆|V |(1) and |V |(1) than we could
have obtained using δ(0) .
The second iteration: Using values from the first iteration, the second iteration yields:
   
−0.0311 [ ] −1.3444 × 10−3
−1
∆P(1) =  0.0231  pu, ∆δ(2) = B′ ∆P(1) =  3.5345 × 10−3  pu
0.0210 2.1841 × 10−3

18
  (2)
    
δ2 0.0105 0.6025◦
  
δ(2) = δ(1) + ∆δ(2) =  δ(2)
3 
= −0.0407  pu =  −2.3292◦ 
(2)
δ4 −0.0222 −1.2739◦
[ ] [ ]
4.4453 × 10−3 [ ′′ ]−1 (1.5) 8.2237 × 10−4
∆Q(1.5) = pu, ∆|V | = B
(2)
∆Q = pu
6.3175 × 10−3 7.3900 × 10−4
[ ] [ ]
(2)
|V3 | 1.0304
|V |(2) = |V | + ∆|V | =
(1) (2)
(2) = pu.
|V |
4
1.0363

The third iteration: Using values from the second iteration, the third iteration yields:
   
3.1029 × 10−3 [ ]−1 (2) 1.4891 × 10 −4

∆P(2) =  −1.4787 × 10−3  pu, ∆δ(3) = B′ ∆P =  −2.4157 × 10−4  pu

−2.0385 × 10 −3 −1.7386 × 10−4
 (3)     
δ2 0.0107 0.6111◦
  
δ(3) = δ(2) + ∆δ(3) =  δ(3) 3 
= −0.0409  pu =  −2.3430◦ 
(3)
δ4 −0.0224 −1.2839◦
[ ] [ ]
−9.5029 × 10−5 [ ′′ ]−1 (2.5) −2.0122 × 10−5
∆Q (2.5)
= pu, ∆|V | = B
(3)
∆Q = pu
−1.9924 × 10−4 −2.1517 × 10−5
[ ] [ ]
(3)
|V3 | 1.0304
|V |(3)
= |V | + ∆|V | =
(2) (3)
(3) = pu.
|V4 | 1.0363

The fourth iteration: Using values from the third iteration, the fourth iteration yields:
   
−3.2246 × 10−4 [ ] −1.6646 × 10−5
−1
∆P(3) =  1.0875 × 10−4  pu, ∆δ(4) = B′ ∆P(3) =  1.8551 × 10−5  pu
1.9823 × 10−4 1.4554 × 10−5
 (4)     
δ2 0.0106 0.6101◦
  
δ(4) = δ(3) + ∆δ(4) =  δ(4) 3 
= −0.0409  pu =  −2.3419◦  .
δ
(4) −0.0224 −1.2831◦
4

At this point we find that using the values of δ(4) and |V |(3) , we get
 
[ −6
] 3.3858 × 10−5
4.9816 × 10
∆Q(3.5) = pu, and ∆P(3.5) =  −8.1853 × 10−6  pu
1.1540 × 10−5
−1.9155 × 10−5

where all mismatches have dropped below 10−4 pu, and hence, if this tolerance is deemed acceptable, we
can stop the iterations. But we did not need to complete the fourth iteration, so in this case the solution is
said to have converged to a tolerance of 10−4 pu in 3.5 iterations.
Now we have the voltage solution
   
V1 1.05̸ 0
 V2   1.06̸ 0.610◦ 
   
 V3  =  1.0304̸ − 2.342◦  pu
V4 1.0363̸ − 1.283◦

19
from which we can calculate the unknown real and reactive power injections:

P1 = |V1 |2 |Y11 | cos θ11 + |V1 | {|Y12 ||V2 | cos(δ1 − δ2 − θ12 ) + |Y13 ||V3 | cos(δ1 − δ3 − θ13 )
+|Y14 ||V4 | cos(δ1 − δ4 − θ14 )} = 0.2000 pu
Q1 = −|V1 |2 |Y11 | sin θ11 + |V1 | {|Y12 ||V2 | sin(δ1 − δ2 − θ12 ) + |Y13 ||V3 | sin(δ1 − δ3 − θ13 )
+|Y14 ||V4 | sin(δ1 − δ4 − θ14 )} = 0.0615 pu
Q2 = −|V2 |2 |Y22 | sin θ22 + |V2 | {|Y21 ||V1 | sin(δ2 − δ1 − θ21 ) + |Y24 ||V4 | sin(δ2 − δ4 − θ24 )} = 0.2351 pu.

These are the same results as those we obtained using the Gauss-Seidel and Newton-Raphson methods. The
power flow solution can be completed by determining the line flows.
It is noteworthy that although this method took half an iteration more than the Newton-Raphson method
took to converge to the same tolerance, the overall computational effort was much less, since we did not
have to compute and factorize the Jacobian at every iteration.

5 Calculation of Power Flows and Losses

So far we have seen how to determine the voltages at all the nodes of the system. With this information the
final step in the power flow solution can be completed, i.e., the flows in the various network elements can be
calculated. The calculation of complex power flow in a shunt element using the shunt admittance and bus
voltage is straightforward, as is the calculation of complex power flow in a series element using the series
admittance and the voltages at the terminal buses. The calculation of flows in an element with a π model,
however, merits comment.
Consider a transmission line with a π model, as shown in figure 5. The complex power Smn entering
the line from bus m must be calculated as the sum of the complex power flowing into the shunt and series
branches. The complex power Snm entering the line from bus n must be similarly calculated. The complex
power loss in the line is given by
Sloss,mn = Smn + Snm . (54)

Smn Snm

Vm Vn

Figure 5: Transmission line with π model.

References
[1] J. D. Glover, M. S. Sarma, T. J. Overbye, Power System Analysis and Design. Fourth Edition: Thomp-
son, 2007.

[2] B. Stott, O. Alsac, “Fast Decoupled Load Flow,” IEEE Transactions on Power Apparatus and Systems,
vol. 93, no. 3, pp. 859–869, May 1974.

20