Manual For A Power Flow Programme: (Version 1.0)

Manual for
A Power Flow Programme

(Version 1.0)
Prepared by
Dr. K.N. Shubhanga

Mr. Chandra S.R.
Department of Electrical Engineering

NITK, Surathkal
Srinivasnagar, Mangalore - 575025
KARNATAKA, INDIA
Version-1.0
Contents
List of Figures iii
List of Tables iv
1 Load Flow Studies: An Introduction 1

1.1 Introduction: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Load flow Problem Formulation: . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2.1 Types of Buses: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.2.2 Handling of Bus Types in Load Flow Problem: . . . . . . . . . . . . 6
1.2.2.1 Handling of PV-Buses: . . . . . . . . . . . . . . . . . . . . 6
1.2.2.2 Handling of Slack Bus: . . . . . . . . . . . . . . . . . . . . 7
1.2.3 Load Flow Solution Methods: . . . . . . . . . . . . . . . . . . . . . 8
2 Newton-Raphson Method 9
2.1 Introduction To Newton-Raphson Method: . . . . . . . . . . . . . . . . . . 9
2.2 Application of the Newton-Raphson Method to Load Flow Solution: . . . . 11
2.2.1 Calculation of Elements of Submatrices: . . . . . . . . . . . . . . . 14
2.3 Newton-Raphson Algorithm: . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.3.1 PV-buses Without Q-limits Specifications: . . . . . . . . . . . . . . 17
2.3.1.1 Load Flow Study Example: . . . . . . . . . . . . . . . . . 19
2.3.2 PV-buses With Q-limits Specifications: . . . . . . . . . . . . . . . . 22
2.3.3 A Case Study to Demonstrate the Effect of Point of Enforcement
of Q-limits Based on Iteration Count: . . . . . . . . . . . . . . . . . 25
2.3.4 Convergence Results for the IEEE Systems: . . . . . . . . . . . . . 25
3 Case Studies with A Test System: NR method 27

3.1 The IEEE 14-bus System . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.2 Newton-Raphson method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.2.1 Format of Data Files . . . . . . . . . . . . . . . . . . . . . . . . . . 28
NITK Surathkal i Electrical Dept.

Version-1.0
4 Fast Decoupled Load Flow Method 33

4.1 Introduction to FDLF Method: . . . . . . . . . . . . . . . . . . . . . . . . 33
4.2 Development of FDLF Equations: . . . . . . . . . . . . . . . . . . . . . . . 33
4.2.1 Formation of B ′ and B ′′ Matrices: . . . . . . . . . . . . . . . . . . . 38
4.3 FDLF Implementation Algorithm: . . . . . . . . . . . . . . . . . . . . . . . 38
4.3.1 PV-buses Without Q-limits Specifications: . . . . . . . . . . . . . . 39
4.3.2 PV-buses With Q-limits Specifications: . . . . . . . . . . . . . . . . 42
4.3.3 Convergence Results for the IEEE Systems: . . . . . . . . . . . . . 45
4.4 Handling of Static Load Models in Load Flows . . . . . . . . . . . . . . . . 46
4.5 Representation HVDC Systems For Power Flow Solution: . . . . . . . . . . 49
4.5.1 DC Link Specification Schemes . . . . . . . . . . . . . . . . . . . . 51
4.5.2 Solution of DC Link Equations . . . . . . . . . . . . . . . . . . . . 52
4.5.3 Incorporation of Control Variable Limits: . . . . . . . . . . . . . . . 54
4.5.3.1 Rectifier and Inverter Transformers’ Tap Adjustments in
Mode-1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.5.3.2 Rectifier and Inverter Transformers’ Tap Adjustments in
Mode-2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4.5.4 AC-DC Load Flow Algorithm: . . . . . . . . . . . . . . . . . . . . . 57
5 Case Studies with A Test System: FDLF method 59

5.1 Format of Data Files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.2 IEEE 14-bus Example for Inclusion of HVDC Link . . . . . . . . . . . . . 62
A Important Expressions Used in Load Flow Programme 69
Bibliography 73
NITK Surathkal ii Electrical Dept.

Version-1.0
List of Figures
1.1 Representation of a typical bus p in load flow studies. . . . . . . . . . . . . 3
2.1 Flow chart for the implementation of Newton-Raphson method . . . . . . . 18

2.2 One line diagram for load flow example. . . . . . . . . . . . . . . . . . . . 19
2.3 Flow chart for the standard bus-type switching logic. . . . . . . . . . . . . 23
3.1 The IEEE 14-bus system . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
4.1 Flow chart for the modified FDLF selective iteration scheme. . . . . . . . . 40
4.2 Flow chart for modified bus-type switching logic for FDLF. . . . . . . . . . 43
4.3 Implementation of load modelling module in load flow. . . . . . . . . . . . 49
4.4 Block diagram to illustrate the interfacing of DC system with AC system
in load flow. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4.5 Flowchart of DC module. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.6 Flowchart for Mode-1 calculations. . . . . . . . . . . . . . . . . . . . . . . 55
4.7 Flowchart for Mode-2 calculations. . . . . . . . . . . . . . . . . . . . . . . 56
4.8 Flowchart for the inclusion of DC calculations in AC load flow. . . . . . . . 58
5.1 The modified IEEE 14-bus system. . . . . . . . . . . . . . . . . . . . . . . 62
A.1 Nominal π Model of transmission lines. . . . . . . . . . . . . . . . . . . . . 69

A.2 Transformer with off-nominal tap. . . . . . . . . . . . . . . . . . . . . . . . 70
A.3 Equivalent circuit of transformer. . . . . . . . . . . . . . . . . . . . . . . . 70
NITK Surathkal iii Electrical Dept.

Version-1.0
List of Tables
2.1 Line data for the example. . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2.2 Bus data for the example. . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.3 Converged load flow results: without Q-limits . . . . . . . . . . . . . . . . 22
2.4 Converged Load flow results: accounting Q-limits at bus 4. . . . . . . . . . 24
2.5 Effect of iteration count as the starting criteria (NR-method). . . . . . . . 25
2.6 IEEE test system results without accounting Q-limits at PV-buses (NR). . 26
2.7 IEEE test system results accounting Q-limits at PV-buses (NR). . . . . . . 26
4.1 Converged load flow results without Q-limits (Selective iteration scheme). 42
4.2 Converged load flow results with Q-limits (Selective iteration scheme) . . 45
4.3 IEEE test system results without accounting Q-limits at PV-buses (FDLF). 46
4.4 IEEE test system results accounting Q-limits at PV-buses (FDLF). . . . . 46
4.5 Study results with voltage- and frequency-dependent load models. . . . . . 48
4.6 DC link data for the IEEE 14-bus system. . . . . . . . . . . . . . . . . . . 58
NITK Surathkal iv Electrical Dept.

Version-1.0
Chapter 1
Load Flow Studies: An Introduction
1.1 Introduction:
One of the most common computational procedures used in power system analysis is the
load flow calculation. Based on the assumed load demands, generation levels and transmis-
sion network configuration, the load flow calculations provides the complete steady-state
information about the system such as the overall voltage profile (bus voltage magnitudes
and bus angles), power flows in the branches, line currents, line losses, and other related
variables for a given set of specified conditions. Such a steady-state snap-shot analysis
provides the basic data for various studies involving planning, design, and operation of
power systems. Application areas of load flow studies are listed below [1, 2]:
1. System planning: This includes some of following tasks:
(a) Transmission planning: It involves design, analysis, sizing of conductors, trans-

formers, reactors, shunt capacitors, and future planning of transmission cir-
cuits.
(b) Interchange studies to identify line over-loadings, violation of bus voltages lim-
its, inappropriately large bus phase angles etc.
(c) Generation adequacy studies and generation planning.
(d) Cost to benefit analysis of system additions.
2. System operation and control: Some of the important studies carried out include
(a) Economic dispatch of the generating stations (i.e.,calculating the power levels
at each generating unit such that the system is operated most economically).
(b) Contingency analysis (analysis of outages and other forced operating condi-
tions).
NITK Surathkal 1 Electrical Dept.

Load Flow Studies: An Introduction Version-1.0
(c) Studies that ensure power pool coordination.

(d) Generation dispatch and load demand control.
(e) Operation tasks related to day to day and week to week operation of the system.
(f) Studies pertaining to scheduling of phase shifters, tap changing transformers,
reactive power controlling devices etc.
(g) Studies related to settings of protective gears.
3. The load flow results provides the basic data for various other studies, for exam-
ple: short circuit analysis, small-signal and transient stability studies, security con-
strained stability analysis: as optimal power flow routine, long-term small-signal
voltage stability studies: continuation power flow programme, etc.
1.2 Load flow Problem Formulation:

The load flow analysis of a power system involves the following steps:
1. Formulation of a suitable mathematical model for the network.
2. Specifications of power and voltage constraints at various buses of the network.
3. Numerical solution of the power flow equations subject to the above constraints.
The load flow calculation is a network solution problem. The necessary equations can
be established by using either the bus or loop frame-of-reference [3]. Due to the inherent
advantages such as simplicity of data preparation and ease of network representation
and computation (exploitation of sparsity), the bus frame-of-reference is normally used
for formulation of load flow problem. The equation describing the performance of the
network of a power system using the bus frame-of-reference in admittance form is given
by
I = YBU S V (1.1)
where
I is the vector of total positive sequence currents injected into the network nodes.
V is the vector of total positive sequence voltages at the network nodes.
YBU S is the network admittance matrix.
If either I or V were known, the solution for the unknown quantities could be obtained
by the solution of linear algebraic equations represented by (1.1). Partly because of
tradition and partly because of the physical characteristics of generation and load, the

terminal conditions at each bus are normally described in terms of a vector of net injected
active and reactive power (P and Q). The injected bus current I is calculated as:
(P + jQ)∗
I= (1.2)
V∗
Combining (1.1) and (1.2) yields
(P + jQ)∗

= YBU S V (1.3)
V∗
Note that the above equation is nonlinear and can not be readily solved by closed
form matrix techniques. Because of this, load flow solutions are obtained by employing
numerical techniques involving iterative procedures.
Let Pp and Qp denote the net injected real and reactive power at bus p as shown in
Figure 1.1.
Ip
Sp Vp
Figure 1.1: Representation of a typical bus p in load flow studies.
Then, the complex power injected at bus p is given by
Sp = Pp + jQp = Vp Ip∗
where Qp represents inductive reactive power and a negative Qp represents capacitive

reactive power.
From (1.1), the injected current at bus p, Ip can be calculated as
n
X
Ip = Ypq Vq , p = 1, 2, ...., n (1.4)
q=1
where n represents the number of buses in a system and Vq = |Vq |6 δq

Using (1.4) we can write an expression for apparent power Sp at bus p as

n
X
Sp = V p Ypq∗ Vq∗ (1.5)
q=1
In the above equation, equating the real and imaginary parts on both sides we get the
active and reactive power balance equations as follows, which constitute the power-flow
equations.
" n
#
X
Pp = Re Vp Ypq∗ Vq∗ (1.6)
p=1
" n
#
X
Qp = Im Vp Ypq∗ Vq∗ (1.7)
q=1
for p = 1, 2, . . . , n
Note that the above set constitutes 2n number of nonlinear equations in 6n variables
with 6 variables at each bus p: bus voltage magnitude |Vp |, bus voltage angle δp , active
power generation PGp , reactive power generation QGp , active power load demand PLp ,
reactive power load demand QLp . Normally, the generation and load power at a bus
are combined to give net injected bus powers, (Pp and Qp ), resulting in 4 variables at
each bus. Thus, in load flow problem 2n nonlinear equations in 4n variables need to be
solved. From this, it is clear that it is not possible to obtain a solution for any unknown
variables unless we reduce the number of “unknowns” to 2n which is equal to the number
of connecting equations. In other words, at each bus two of the four variables needs to be
a priori specified and the remaining two variables are required to be calculated using the
two power balance equations available at that bus. This has lead to the general practice
of identification of different types of buses in the network.
1.2.1 Types of Buses:

To enable a priori fixation of some of the variables, the 6 variables at each bus are classified
as follows [1, 4]:
1. Demand variables : The demand variables PLp and QLp , can be specified from the
knowledge of customer load demand.
2. Generation (or control) variables: The generation variables PGp and QGp , can be
specified using the scheduled generation data.
3. State variables: The variables such as voltage magnitude |Vp | and bus angle δp
completely define the state of a system. These variables are normally treated as

state (or unknown) variables.
Based on this classification, the buses in a system are categorized into 3 types as follows:
1. Load buses (PQ-buses): A load bus is any bus that does not have generation. Here,
real power and reactive powers are specified, because for loads, one generally specifies
the real power demand, PLp and the reactive power demand, QLp . In some cases,
such a bus need not have a load, it may simply be an interconnection point for two
or more lines. Those buses are also treated as load buses with PLp = QLp = 0.
2. Generator buses (PV-buses): A generator bus is usually called a voltage controlled

or PV-bus. Here, real power generation, PGp and voltage magnitude, |Vp | are nor-
mally specified. At such buses voltage magnitude is held constant at a specified value
by adjusting the reactive power output of the generator. Based on its capability,
maximum and minimum limits on reactive power generation (QGmax and QGmin )
are specified at PV-buses. Certain PV-buses without real power generation may
have voltage control capability, such buses are also designated as voltage-controlled
buses at which the real power generation is set to zero. A generator bus could also
have load connected to it. Based on whether load is present or not at a PV-bus, for
the purpose of ease of handling, a PV-bus is further categorized as follows:
(a) Simple PV-bus : These are buses in which there is no load connected.
(b) PV-bus with loads : In addition to normal PV- bus specification, loads are also
connected at these buses.
3. Swing bus: This is also known as a reference bus or slack bus. Swing bus is a special
type of generator bus that is needed by the solution processes. This requirement
comes from the fact that in load flow analysis, the losses in the system cannot
be determined unless the converged solution is obtained. This slack in the power
balance is accounted in the reference bus generation. Therefore in a power system,
! ! !
Sum of complex Sum of complex Total (complex) power
= +
power of generators power of loads loss in transmission lines
This model assumes that only the slack bus generates the necessary real and reactive
power required for overcoming all system losses. Hence, for the slack bus, the
magnitude and phase angle of bus voltage are specified and real and reactive power
generations are obtained at the end of the load flow solution. The voltage angle of
the slack bus serves as reference for the angles of all other bus voltages.
A swing bus could also have load connected to it. Based on whether load is present
or not swing bus is further categorized into:

(a) Swing bus without load

(b) Swing bus with load
Thus, in a n bus system, there will be one slack bus and there may be m PV-buses and
rest of the buses are designated as PQ-buses. The load flow problem formulation results
in 2n nonlinear algebraic equations in 4n variables. Of these 4n variables, 2 variables
are specified for the slack bus, 2m variables are specified for PV-buses and 2(n − m − 1)
variables are specified for PQ-buses. These specifications make 4n − (2 + 2m + 2(n −
m − 1)) = 2n variables to be determined from the load analysis. However, in load
flow problem, the nonlinear equations are solved only for 2n − m − 2 unknown variables
(voltage magnitudes and bus angles). The remaining (m + 2) unknown variables (such as
PGs , QGs at slack bus and QGp at m PV-buses) are computed. In the solution of nonlinear
equations, normally the unknown voltage magnitudes are initialized to 1.0 p.u. and angles
are initialized to 0o . In some cases, the unknown bus voltage vector is initialized to the
specified slack bus voltage, which is referred to as flat voltage start [4].
NOTE:
1. The following convention is followed with respect to apparent powers. The apparent
power drawn by the load is SLp = PLp + jQLp . Where QLp denotes inductive load
and a negative QLp represents capacitive load.
2. The apparent power injected by the generator is SGp = PGp + jQGp . Where QGp
represents inductive reactive power and a negative QGp represents capacitive reac-
tive power.
1.2.2 Handling of Bus Types in Load Flow Problem:

In load flow study, a separate strategy is followed to handle different types of buses. In
the following lines the procedure of handling PV-buses and slack bus is explained:
1.2.2.1 Handling of PV-Buses:
Ii is known that by controlling the the reactive power generation at PV-bus, the voltage
magnitude at the PV-bus is maintained at the specified value. If limits on reactive power
generation is specified at PV-buses, the checking of Q-limits are carried out in the following
manner:
1. While checking for Q-limit specifications at a Simple PV-bus, the computed reactive
power (= Qp ) is considered.

2. While checking for Q-limits at a PV-bus with load, the net reactive power generation
(QGp ) is considered. This QGp is equal to the sum of reactive power load QLp at
that bus and the computed reactive power Qp , (i.e, QGp = QLp + Qp ).
If any Q-limit violation is detected, a PV-bus is handled in any one of the following ways
[5]-[7]:
1. Bus-type switching : In any iteration, if Q-limits are violated at a PV-bus, then it

is switched to PQ-type with the reactive generation set at the limiting value. In any
consequent iteration if Q-limits are satisfied at such a bus, then that bus is reverted
back to PV-type with the original bus voltage specification. This is referred to as
back-off. Thus, handling of PV-buses in load flow analysis by a way of bus switching
calls for repeated adjustments of number of unknowns.
2. Adjusting the specified voltage at PV-buses: In any iteration, if Q-limits are violated
at a PV-bus, then the specified voltage at the violated PV-bus is adjusted in such a
way that it remains as PV-bus and meets the Q constraints. This procedure avoids
the readjustment of number unknowns.
In bus-type switching, it is important to consider the time at which the Q-limits are
enforced at PV-buses in the solution process as it may effect the convergence charac-
teristics of an algorithm. If not handled properly, such an enforcement may slow down
the convergence or even cause the solution to oscillate or even diverge [9]. The main
parameters used for the starting criteria are:
1. Iteration count [8]: In this case, the bus-type switching at each PV-buses is delayed
until the second or later iterations.
2. Bus power mismatches [9]: In this case, the bus-type switching at each PV-buses
is delayed until the reactive power mismatch is sufficiently small (i.e., load flow
solution is moderately converged).
1.2.2.2 Handling of Slack Bus:
The real and reactive power generations at slack bus are calculated from the converged
load flow results as follows: n
X
SGs = Ss = Vs Ysq∗ Vq∗ (1.8)
q=1
If any load is present at slack bus, then the real and reactive components of load are
simply added to the computed values to obtain the net power generation at slack bus.
SGs = (PGs + jQGs ) = Ss + (PLs + jQLs ) (1.9)

1.2.3 Load Flow Solution Methods:

Load flow studies for power system analysis are based on static network models assuming
sinusoidal steady state. As a result the equations are algebraic in form rather than dif-
ferential. The usual power flow problem formulation does not consider time variation of
loads, generation, or network configuration. The load flow problem formulation is nonlin-
ear in nature even for a linear transmission network since power specification is employed
which is a product of voltage and current. Additional nonlinearities arise due to specifi-
cation and use of complex voltages and currents. Also, there are transmission component
nonlinearities, which may be considered (such as tap changing transformers in which the
tap is adjusted to hold a given bus voltage magnitude fixed). Thus, load flow problem
generally involves the simultaneous solution of many nonlinear algebraic equations. Since
it is not feasible to obtain solution in a closed form directly for nonlinear system of equa-
tions, numerical techniques are normally employed involving iterative procedure. Some of
the popularly used methods are fixed-point (Gauss) methods and Newton based methods.
The earliest algorithms were based on the Gauss-Seidel methods, which are very simple
and reliable. However, it takes extremely large computation times when applied to sys-
tems of increasing size due to its poor convergence characteristics [1, 3, 10]. To overcome
the drawbacks in the these methods, a new method called ”Newton’s method” was intro-
duced in early sixties [8]. This method basically comprises the repeated solution of a large
set of linear equations which are obtained by linearizing the load flow equations. New-
ton’s methods have very powerful convergence properties, but with poor computational
efficiency and also it requires large storage memory. To improve the computational speed
and storage requirements of the Newton Raphson method, B. Stoot and O. Alasc have
developed ‘fast decoupled load flow method’ in early 1970’s [7] which is a modification of
the Newton method.

Version-1.0
Chapter 2
Newton-Raphson Method
2.1 Introduction To Newton-Raphson Method:

The application of Newton’s method to load flow problem is essentially the n-dimensional
generalization of the well known Newton-Raphson method [6, 11] for the solution of a
nonlinear equation in one variable. Taylor’s series expansion for a function of two or more
variables is the basis for the Newton-Raphson method of solving the load-flow problem
[12]. Let the following represent n nonlinear equations in n unknowns :
f1 (x1 , x2 , · · · , xn ) = b1
f2 (x1 , x2 , · · · , xn ) = b2 (2.1)
··············· ··· ···
fn (x1 , x2 , · · · , xn ) = bn
If the iterations start with an initial estimate of x01 , x02 , · · · , x0n for n unknowns and if
∆x1 , ∆x2 , · · · , ∆xn are the corrections necessary to the estimates so that the equations
are exactly satisfied, we have
f1 (x01 + ∆x1 , x02 + ∆x2 , · · · , x0n + ∆xn ) = b1 (2.2)

f2 (x01 + ∆x1 , x02 + ∆x2 , · · · , x0n + ∆xn ) = b2
···································· ··· ···
fn (x01 + ∆x1 , x02 + ∆x2 , · · · , x0n + ∆xn ) = bn
Each of the above equations can be expanded using Taylor’s theorem. The expanded

Newton-Raphson Method Version-1.0
form of the ith equation is

∂fi ∂fi
 
0 0 0

 f1 (x1 , x2 , · · · , xn ) + ∆x1 + ∆x2 



 ∂x1 ∂x2 


fi (x01 + ∆x1 , . . . , x0n + ∆xn ) = ∂fi = bi
+...+ ∆xn + terms with



 ∂xn 



higher powers of ∆x1 , . . . , ∆xn
 
The terms of higher powers can be neglected, if our initial estimate is close to the true
solution. The resulting linear set of equations in matrix form is
 
∂f1 ∂f1 ∂f1
···

 ∂x10 0 ∂x2 ∂xn
 0
 
     o
b1 − f1 (x01 , x02 , · · · , x0n )  ∂f ∂f
2 2 ∂f2  ∆x1
···
 
 b2 − f1 (x01 , x02 , · · · , x0n )   ∂x1 0 ∂x2 0 ∂xn 0  
    
 =  ∆x2 
 (2.3)
 ·····················    ··· 
    
bn − f1 (x01 , x02 , · · · , x0n )

 ··· ··· ··· ··· 
 ∆xn
 
 
 
 ∂fn ∂fn ∂fn 
···
∂x1 0 ∂x2 0 ∂xn 0
or in general,
∆F = J ∆x
where J is referred to as the Jacobian. If the estimates x01 , · · · , x0n were exact, then ∆F
and ∆x would be zero. However, as x01 , · · · , x0n are only estimates, the errors ∆F are
finite. The above equation provides a linearized relationship between the errors ∆F and
the corrections ∆x through the Jacobian of the simultaneous equations. A solution for
∆x can be obtained by applying any suitable method for the solution of a set of linear
equations. An updated value of x can be calculated as
x1i = x0i + ∆xi
The process is repeated until the errors ∆Fi are lower than a specified tolerance. The
Jacobian has to be recalculated at each step.
Newton’s method possesses quadratic convergence characteristics, a feature not shared
by other methods. Its drawback is that the initial guess used should be close to the solution
for the process to converge. For power systems this is not a serious drawback since a good
initial guess is generally available from experience [8].

2.2 Application of the Newton-Raphson Method to

Load Flow Solution:
In general, there are three types of approaches available in literature for solving load flow
equations using Newton’s method [10]. Those are:
1. Rectangular power-mismatch version: This version uses the real and imaginary parts
of the bus voltages (i.e., rectangular co-ordinate system ei + jfi ) as variables.
2. Polar power-mismatch version: This version uses the magnitudes and angles of the
bus voltages (i.e., polar co-ordinate system (|V |6 δ) as variables.
3. Current-mismatch version: In this version, rectangular and polar coordinate version

are constructed using the real and imaginary parts of the complex current expression.
Polar power-mismatch version is widely used version of the Newton’s method, because of
its more reliability and accuracy over rectangular power-mismatch version [10]. So in this
report, polar power-mismatch version is considered [6, 11]. In this version we denote
Ypq = Gpq + jBpq

Vp = |Vp |6 δp = |Vp |(cos 6 δp + j sin 6 δp ) (2.4)
Vq = |Vq |6 δq = |Vq |(cos 6 δq + j sin 6 δq )
Using (2.4) in (1.5), the complex power injected at bus p can be represented as
n
X
Pp + jQp = |Vp |6 δp ((Gpq + jBpq )|Vq |6 δq )∗ (2.5)
q=1
Expanding (2.5) and equating real and imaginary parts on both sides, we get
n
X
Pp = |Vp | ((Gpq cos δpq + Bpq sin δpq) )|Vq |) (2.6)
q=1
Xn
Qp = |Vp | ((Gpq sin δpq − Bpq cos δpq) )|Vq |) (2.7)
q=1
where p = 1, 2, . . . , n and δpq = δp − δq

From the above, it can be seen that Pp and Qp at each bus are functions of voltage

magnitude |V | and angle δ at all buses and the load flow equations are formulated as:
P1 (δ1 , · · · , δn , |V1|, · · · , |Vn |) = P1sp (2.8)

······························
Pn (δ1 , · · · , δn , |V1|, · · · , |Vn |) = Pnsp
Q1 (δ1 , · · · , δn , |V1 |, · · · , |Vn |) = Qsp
1
······························
Qn (δ1 , · · · , δn , |V1 |, · · · , |Vn |) = Qsp
n
Note that the above equation is similar to (2.1). The RHS terms denote the power
specifications at each bus.
Computation of bus power mismatches:
Having specified real power generation PGp , and real power load PLp at bus p, the net
specified real power injected into the network at bus p is computed as
Ppsp = PGp − PLp
Similarly, having specified reactive power generation QGp , and reactive power load QLp at
bus p, the net specified reactive power injected into the network at bus p is computed as
Qsp
p = QGp − QLp
NOTE:
1. The following convention is followed with respect to apparent powers. The apparent
power drawn by the load is SLp = PLp + jQLp . Where QLp denotes inductive load.
And a negative QLp represents capacitive load.
2. The apparent power injected by the generator is SGp = PGp +jQGp . Where QGp rep-
resents inductive reactive power. And a negative QGp represents capacitive reactive
power.
Following the general procedure described earlier for the application of the Newton’s
method, we have

0
∂P1 ∂P1 ∂P1 ∂P1

··· ···
 ∂δ1 ∂δn ∂|V1 | ∂|Vn | 
 
 
 
 ··· ··· ··· ··· ··· ··· 
 
0 0
P1sp
   
− P1 (δ10 , · · · , δn0 , |V10 |, · · · , |Vn0 |) 


 ∆δ1
 ∂Pn ∂Pn ∂Pn ∂Pn 
······························ ···
   
  
 ∂δ ··· ···   

 P sp − P (δ0 , · · · , δ0 , |V 0 |, · · · , |V 0 |)

1 ∂δn ∂|V1 | ∂|Vn |  
 
 n n 1 n 1 n
   ∆δn 
=  (2.9)
    
 sp
 Q1 − Q1 (δ10 , · · · , δn0 , |V10 |, · · · , |Vn0 |)
  

 ∂Q1

∂Q1 ∂Q1 ∂Q1  
  ∆|V1 | 
··· ···
  
······························ ··· 
     
   ∂δ1 ∂δn ∂|V1 | ∂|Vn |  
sp
Qn − Qn (δ10 , · · · , δn0 , |V10 |, · · · , |Vn0 |)
 



 ∆|Vn |
 ··· ··· ··· ··· ··· ··· 
 
 
 
 
 ∂Qn ∂Qn ∂Qn ∂Qn 
··· ···
∂δ1 ∂δn ∂|V1 | ∂|Vn |
The above equation is similar to (2.3). The LHS vector represents the power mismatch
at each buses [∆P , ∆Q]T . The real power mismatch at bus p, pertaining to the first
iteration (k=0), ∆Pp0 is calculated as
∆Pp0 = Ppsp − Pp0 (2.10)
Using (2.6), the above equation is re-written as

n
X
∆Pp0 = Ppsp − |Vp0 | 0
((Gpq cos δpq 0
+ Bpq sin δpq) )|Vq0 |) (2.11)
q=1
p = 1, 2, . . . , n and p 6= slack bus
Similarly, reactive power mismatch at bus p, pertaining to the first iteration (k=0),
∆Q0p is calculated as
∆Q0p = Qsp 0
p − Qp (2.12)
Using (2.7), the above equation is re-written as

n
X
∆Q0p = Qsp 0
p − |Vp |
0
((Gpq sin δpq 0
− Bpq cos δpq) )|Vq0 |) (2.13)
q=1
p = 1, 2, . . . , n; p 6= slack bus, p 6= PV-bus

To bring-in symmetry in the elements of the Jacobian, the linearized load flow equation
given in (2.9) is modified (for the k th iteration) as follows:
" #(k−1) " #(k−1) " #(k−1)
∆P H N ∆δ
= ∆|V | (2.14)
∆Q M L |V |
where the submatrices elements are defined as
∂P ∂P
[H](n−1)×(n−1) = [N](n−1)×(n−1−m) = ∂|V |
∂δ
|V |
∂Q ∂Q
[M](n−1−m)×(n−1) = [L](n−1−m)×(n−1−m) = ∂|V |
∂δ
|V |
NOTE:
1. In H submatrix the row and column corresponding to slack bus are not present.
2. In L submatrix the rows and columns corresponding to slack bus and PV-buses are
not present.
3. In N submatrix the row and column corresponding to slack bus and columns corre-
sponding to PV-buses are not present.
4. In M submatrix the row and column corresponding to slack bus and rows corre-
sponding to PV-buses are not present.
2.2.1 Calculation of Elements of Submatrices:

Each submatrices (H, N, M, L) will have its own diagonal and off-diagonal elements.
When p = q: it represents the diagonal elements and p 6= q represents the off-diagonal
elements. The elements H, N, M, and L matrices are calculated as follows.
1. Calculation of H submatrix:
∂Pp
(a) Off-diagonal element Hpq is given by, Hpq =
∂δq
Rewriting (2.6) we,have
n
X
Pp = |Vp | ((Gpq cos δpq + Bpq sin δpq) )|Vq |) (2.15)
q=1
where p = 1,2,. . . ,n

And differentiating Pp with respect to δq , we have

 
∂Pp ∂  (Gp1 cos(δp − δ1 ) + Bp1 sin(δp − δ1 ))|V1 | + . . .
 

= |Vp | (Gpq cos(δp − δq ) + Bpq sin(δp − δq ))|Vq | + . . .
∂δq ∂δq  
(Gpn cos(δp − δn ) + Bpn sin(δp − δn ))|Vn |
 
where p = 1,2,. . . ,n
Simplifying above equation, we get
Hpq = |Vp ||Vq |(Gpq sin δpq − Bpq cos δpq) ) (2.16)
∂Pp
(b) Diagonal elements Hpp is given by, Hpp =
∂δp
Differentiating (2.6) with respect to δp , we,have

 
∂Pp ∂  (Gp1 cos(δp − δ1 ) + Bp1 sin(δp − δ1 ))|V1 | + . . .
 

= |Vp | (Gpq cos(δp − δq ) + Bpq sin(δp − δq ))|Vq | + . . .
∂δp ∂δp  
(Gpn cos(δp − δn ) + Bpn sin(δp − δn ))|Vn |
 
where p = 1,2,. . . ,n
Simplifying the above equation, we get
n
X
Hpp = |Vp | −Gpq sin(δp − δq ) + Bpq cos(δp − δq ))|Vq | (2.17)
q=1
q 6= p
Using (2.7), the above expression can be simplified as
Hpp = −Qp − Bpp |Vp |2 (2.18)
Following the procedure detailed above we can obtain the expressions for the ele-
ments of other submatrices, M, N and L as follows:
2. Calculation of M submatrix:
(a) Off diagonal element Mpq is given by:
∂Qp
Mpq = = −|Vp ||Vq |(Gpq cos δpq + Bpq sin δpq) ) (2.19)
∂δq

(b) Diagonal element Mpp is given by:
∂Qp
Mpp = = Pp − Gpp |Vp |2 (2.20)
∂δp
3. Calculation of N submatrix:
(a) Off-diagonal element Npq is given by:
∂Pp
Npq = ∂|Vq |
= |Vp ||Vq |(Gpq cos θpq + Bpq sin θpq) ) (2.21)
|Vq |
(b) Diagonal elements Npp is given by:
∂Pp
Npp = ∂|Vp |
= Pp + Gpp |Vp |2 (2.22)
|Vp |
4. Calculation of L submatrix:
(a) Off diagonal element Lpq is given by:
∂Qp
Lpq = ∂|Vq |
= |Vp ||Vq |(Gpq sin δpq − Bpq cos δpq) ) (2.23)
|Vq |
(b) Diagonal element Lpp is given by:
∂Qp
Lpp = ∂|Vp |
= Qp − Bpp |Vp |2 (2.24)
|Vp |
Using (2.18), we can write that
Lpp = −Hpp − 2Bpp |Vp |2 (2.25)
From above expressions it can be seen that Hpq = Lpq and Npq = −Mpq . This is due to
the fact that the voltage magnitude increment ∆|V | at a PQ-bus is divided by |V |.
NOTE:
1. The real and reactive power mismatch vectors are updated along with the Jacobian
elements at the beginning of solving (2.14). The maximum value of mismatch vector,
[∆P 0 , ∆Q0 ]T will be considerable at the beginning of power flow solution depending
on the closeness of the initial guess (|V |0 , δ0 ) to the final solution. If the iterative

solution process tends to converge, then the maximum value of the mismatch vector
tends to zero. Normally, the solution process is terminated if the maximum value of
the power mismatch falls below a tolerance factor ǫ (which is set to 0.001 or 0.0001),
i.e. if
max|[∆P , ∆Q]T | ≤ ǫ
then the load flow solution is said to be converged.
2. In each of the (k − 1)th iterative step, a solution of the linearized system of equation
is carried out. Here, k denotes the iterative count and is set equal to 1.
3. In (k−1)th iterative step, the mismatch vector [∆P (k−1) , ∆Q(k−1) ]T and the Jacobian
elements are calculated using |V |(k−1) and δ (k−1) values.
(k−1)
4. In (k − 1)th iterative step, the correction vector [∆δ (k−1) , ∆|V |
|V |(k−1)
] is obtained by
solving the linear system of equations. Using this correction vector the voltage
magnitude and angle are updated at bus p as follows:
δpk = δp(k−1) + ∆δp(k−1) (2.26)
∆|Vp |(k−1)
|Vp |k = |Vp |(k−1) + |Vp |(k−1)
|Vp |(k−1)
or
∆|Vp |(k−1)

k (k−1)
|Vp | = |Vp | 1+ (2.27)
|Vp |(k−1)
2.3 Newton-Raphson Algorithm:

In this section the steps involved in the N-R method of load flow analysis is presented.
The algorithm has been implemented for the following two different cases:
1. PV-buses without Q-limits specifications.
2. PV-buses with Q-limits specifications.
2.3.1 PV-buses Without Q-limits Specifications:

The implementation steps are depicted in flow chart shown in Figure 2.1. The ‘Q-bit’ is
set equal to 0 for for this case.

START
Read system data
Construct Y−bus matrix
Assume initial bus voltages as ,

|V p | = 1.0 p.u, for all PQ−buses ,
sp
|Vp | = Vp for all PV−buses & at slack bus
o
bus angles δ = 0 for all buses
Set k = 1
Compute Pp and Q p
Is
Yes
Q−bit = 1 A B
?
No
Compute P and Q
k = k+1
Test
for convergence Yes
max [ P , Q ] <ε
? compute line flows
No & slack bus power
Form the Jacobian matrix (J)
END
Compute voltage magnitudes and angle
corrections
Update bus voltages
Figure 2.1: Flow chart for the implementation of Newton-Raphson method

2.3.1.1 Load Flow Study Example:
The single line diagram of a 4 bus, 2 machine power system is shown in Figure 2.2. The
system details have been adopted from Example 9.2 [11]. Base used is 100 MVA.
1
1
0
0
1 2
0
1
3 4
Figure 2.2: One line diagram for load flow example.
The bus data and line data are as follows:
From To r x y(total)
(pu) (pu) (pu)
1 2 0.01008 0.05040 0.1025
1 3 0.00744 0.03720 0.0775
2 4 0.00744 0.03720 0.0775
3 4 0.01272 0.06360 0.1275
Table 2.1: Line data for the example.
Bus Generation Load Type

No.p PGp QGp PLp QLp |Vp |
(MW) (MVAr) (MW) (MVAr) (pu)
1 - - 50 30.99 1.00 Slack-bus
2 0 0 170 105.35 - PQ-bus
3 0 0 200 123.94 - PQ-bus
4 318 - 80 49.58 1.02 PV-bus
Table 2.2: Bus data for the example.
NOTE: The data files have been provided in the folder 4_bus_example.

1. Constructing the bus admittance matrix YBU S , using bus admittance algorithm. we
have,
 
8.9852 − j44.8360 −3.8156 + j19.0781 −5.1696 + j25.8478 0
 −3.8156 + j19.0781 8.9852 − j44.8360 0 −5.1696 + j25.8478
 

 
 −5.1696 + j25.8478 0 8.1933 − j40.8638 −3.0237 + j15.1185 
 
0 −5.1696 + j25.8478 −3.0237 + j15.1185 8.1933 − j40.8638
 
 
2. Assume |V1 | = 1.0 pu, |V2 | = 1.0 pu and |V3 | = 1.0 pu. And δ2 , δ3 , δ4 = 0. And
the convergence factor ǫ = 0.001
Iteration 1:
1. The specified bus powers and calculated bus powers (using the initial bus voltages)
are :
Bus No.p Ppsp Pp Qsp
p Qp
2 -1.7 -0.1034 -1.0535 -0.6070
3 -2.0 -0.0605 -1.2394 -0.4049
4 2.38 0.1671 - 0.7291
2. Using the computed bus power mismatches and the Jacobian elements we have (as
per (2.14)),
2 32 3 2 3
2 45.4429 0 −26.3648 8.8818 0 ∆δ2 −1.5966
6 76 7 6 7
6 76 7 6 7
6 76 7 6 7
3 6
6 0 41.2687 −15.4209 0 8.1328 7
76
6 ∆δ3 7 6 −1.9395 7
7 6 7
6 76 7 6 7
6 76 7 6 7
∆δ4
6 76 7 6 7
7 = 6 2.2129 7 (2.28)
6 −26.3648
4 6 −15.4209 41.7857 −5.2730 −3.0842 7 6
76 7 6 7
6 76 7 6 7
6 76 7 6 7
6 76 ∆|V2 | 7 6 7
2 6
6 −9.0886 0 5.2730 44.2290 0 76
76 |V2 |
7 6 −0.4465 7
7 6 7
6 76 7 6 7
4 54 5 4 5
∆|V3 |
3 0 −8.2537 3.0842 0 40.4590 |V3 |
−0.8345
3. Solving the above equation, we get the corrections to the voltage magnitudes and
angles as:
Bus No.p ∆δp ∆|Vp |

(deg) (pu)
2 -0.9309 -0.0166
3 -1.7879 -0.0290
4 1.5438

4. At the end of the first iteration, the set of updated bus voltages are:
Bus No.p 1 2 3 4
δp (deg) 0 -0.9309 -1.7879 1.5438
|Vp | (per unit) 1.00 0.9834 0.9710 1.02
Iteration 2:
1. The specified bus powers and the bus powers calculated (using the updated bus
voltages at the end of 1st iteration) are:

p Qp
2 -1.7 -1.6677 -1.0535 -1.0193
3 -2.0 -1.9355 -1.2394 -1.1774
4 2.38 2.3441
per (2.14)),
    
2 44.37499 0 −25.6778 7.0208 0 ∆δ2 −0.0323
    
 
 
 
   
3  0 39.7018 −14.7736 0 5.7887   ∆δ3   −0.0645 

    
    
    
   
4  −26.1256 −15.1217 41.2473 −4.0609 −2.1193   ∆δ4  =  0.0359 

    
    
    
   ∆|V2 |   
2 
 −10.3562 0 6.2998 42.3363 0 
  |V2 |  
  −0.0342 

    
    
∆|V3 |
3 0 −9.6597 3.8597 0 37.3470 |V3 |
−0.0620
angles as:

(deg) (pu)
2 -0.0451 -0.0009
3 -0.00841 -0.0019
4 0.0207
At the end of 2nd iteration: The specified and the bus powers calculated (using the
updated bus voltages at the end of 2nd iteration) are:

Bus No.p Ppsp Pp ∆Pp Qsp

p Qp ∆Qp
2 -1.7 -1.7000 0.00 -1.0535 -1.0535 0.00
3 -2.0 -1.9999 0.0001 -1.2394 -1.2392 0.0002
4 2.38 2.38 0.0 - - -
From above table, it is observed that the calculated bus power mismatches falls below ǫ.
The converged results are shown in Table 2.3 :
Bus No. p |Vp | δp PGp QGp PLp QLp

(pu) (deg) (pu) (pu) (pu) (pu)
1 1.000000 0.000000 1.867948 1.144877 0.500000 0.309900
2 0.982422 -0.976035 0.000000 0.000000 1.700000 1.053500
3 0.969010 -1.871975 0.000000 0.000000 2.000000 1.239400
4 1.020000 1.523147 3.180000 1.814193 0.800000 0.495800
Table 2.3: Converged load flow results: without Q-limits
2.3.2 PV-buses With Q-limits Specifications:

In N-R method, PV-buses are handled by a way of bus-type switching and iteration count
is used as the starting criteria since it is easy to implement and is well suited for the N-R
method because of its unified nature of its solution for [∆δ, ∆|V |]T . The implementation
logic is shown in Figure 2.3. The ‘Q-bit’ set to 1 in Figure 2.1.
For the test system shown in Figure 2.2, the maximum reactive power limit at PV-bus
4 is QGmax = 0.4 and the minimum reactive power limit is QGmin = 0. The load flow
results are given as follows. Note that in the results given, the starting criteria is not
accounted just to demonstrate the intermediate steps. This done by setting start_cri =
0 in jacob_form.m
Iteration 1:
1. The specified bus powers and the bus powers calculated (using the initial bus volt-
ages) are :

p Qp
2 -1.7 -0.1034 -1.0535 -0.6070
3 -2.0 -0.0605 -1.2394 -0.4049
4 2.38 0.1671 - 0.7291

Is
No
Iteration count
k>2
?
Yes
Compute Q G = Q + Q L for all PV−buses
Is Are
Yes there any Q−limit No there any previously No
violation at any violated
PV−buses PV−buses
? ?
Yes
Treat violated PV−buses as PQ−buses
sp
with Q G = Q Gmax (or) QGmin , Set | Vp | = | V p | at back−off buses
whichever limit is violated
Readjust the number of PV and PQ buses
sp
Recompute bus powers (Pp , Qp ) and Qp
Figure 2.3: Flow chart for the standard bus-type switching logic.
2. Here, at bus 4 a real load of PL4 = 0.8 pu. and a reactive load of QL4 = 0.4958
pu. is present. While checking the Q-limits at PV-bus, the required reactive power
generation is calculated as: QG4 = Q4 + QL4 = 0.7291 + 0.4958 = 1.2249 pu. Note
that QG4 > QGmax . Therefore, this bus is treated as a PQ-bus with QG4 = 0.4 pu.
Further, Qsp sp
4 = QG4 − QL4 = 0.4 - 0.4958 = -0.0958 pu. and ∆Q4 = Q4 − Q4 =
-0.0958 - 0.7291 = -0.8249 pu. The initial bus voltage magnitude is 1.02 pu.
per (2.14))
∆δ2
2 32 3 2 3
2 45.4429 0 −26.3648 8.8818 0 −5.2730 −1.5966
6 76 7 6 7
6 76 7 6 7
∆δ3
6 76 7 6 7
3 6
6 0 41.2687 −15.4209 0 8.1328 −3.0842 776
6 7 6 −1.9395 7
7 6 7
6 76 7 6 7
6 76 7 6 7
∆δ4
6 76 7 6 7
4 6 −26.3648 −15.4209 41.7857 −5.2730 −3.0842 8.6914 7 6 7 6 2.2129 7
6 76 7 6 7
7=6
6 76 7 6 7
6 76 7
6 76 ∆|V2 | 7 6 7
2 6 −9.0886
6 0 5.2730 44.2290 0 −26.3648 7
76
6 |V2 |
7 6 −0.4465 7
7 6 7
6 76 7 6 7
6 76 7 6 7
6 76 ∆|V3 | 7 6 7
3 6
6 0 −8.2537 3.0842 0 40.4590 −15.4209 7
76
6
|V3 |
7 6 −0.8345 7
7 6 7
6 76 7 6 7
4 54 5 4 5
4 ∆|V3 |
5.2730 3.0842 −8.3571 −26.3648 −15.4209 43.2438 −0.8249
|V3 |

The above equation is identical to (2.28) except for an additional row and column
pertaining to bus 4.
angles as:

(deg) (pu)
2 -5216 -0.0512
3 -1.5240 -0.0511
4 2.2467 -0.0591
5. At the end of the first iteration the updated bus voltages are:
Bus No. p 1 2 3 4
δp (deg) 0 -0.5216 -1.5240 2.2467
|Vp | (per unit) 1.00 0.9488 0.9489 0.9609
At the end of 3rd iteration: The specified bus powers and the bus powers calculated
(using the updated bus voltages at the end of 3rd iteration) are:

p Qp ∆Qp
2 -1.7 -1.7000 0.000 -1.0535 -1.0535 0.000
3 -2.0 -1.9999 0.000 -1.2394 -1.2394 0.000
4 2.38 2.38 0.000 -0.0958 -0.0958 0.000
From above table, it is observed that the calculated bus power mismatches fall below ǫ.
The converged results are shown in Table 2.4 :
Bus No.p |Vp | δp PGp QGp PLp QLp

1 1.000000 0.000000 1.884566 2.667856 0.500000 0.309900
2 0.941383 -0.533208 0.000000 0.000000 1.700000 1.053500
3 0.942312 -1.631389 0.000000 0.000000 2.000000 1.239400
4 0.951521 2.636509 3.180000 0.400001 0.800000 0.495800
Table 2.4: Converged Load flow results: accounting Q-limits at bus 4.
NOTE:
• It is observed that with the enforcement of Q-limits specifications at PV-buses causes

the Newton’s method to take more number of iterations.

• When Q-limits at PV-buses are handled by a way of bus-type switching it requires

resizing of jacobian whenever there is any new Q-limit violations and back-offs in
each iterative step. In MATLAB, this involves the solution of a linear system of
equations in the usual manner (see Appendix).
• Since Q-limit has been imposed at bus 4, the bus voltage is no longer held at 1.02
pu. which is the specified value.
2.3.3 A Case Study to Demonstrate the Effect of Point of En-

forcement of Q-limits Based on Iteration Count:
Case studies are carried out on IEEE 30 and 300 bus systems [13], to demonstrate the
effect of point of enforcement of Q-limits using ‘iteration count’ (k) as the starting criteria.
The results are presented in the following lines:
No of iterations
Type of system At starting k > 1 k > 2 k > 3
IEEE 30 bus system 6 3 3 2
IEEE 300 bus system NC NC 8 4
Table 2.5: Effect of iteration count as the starting criteria (NR-method).
From Table 2.5, it is observed that the enforcement of Q-limits at starting causes
divergence in the IEEE 300-bus system and also it makes the solution to take more
number of iterations in the IEEE 30-bus system. For value of k > 2, the IEEE 300-bus
system converges. For k > 3, no Q-limits are enforced and hence the solution converges
without accounting Q-limits at PV-buses for the IEEE 30-bus system. However, for k > 3,
the solution converges within 4 iterations for the IEEE 300-bus system. In this report a
value of k > 2 has been used as the starting criteria for all systems. This done by setting
start_cri = 2 in jacob_form.m.
2.3.4 Convergence Results for the IEEE Systems:

The results obtained using the IEEE test systems [13, 14] are presented in Tables 2.6 and
2.7. In the tables NC denotes failure to converge. Base used is 100 MVA. The approximate
execution time is obtained in a P4 computer (3 GHz, 512 MB DDR 400 MHz memory,
910 GL chipset) and it also includes the time for report generation.

IEEE system Iterations Execution time

(s)
14-bus 2 0.032
30-bus 2 0.032
57-bus 3 0.047
118-bus 3 0.094
145-bus NC -
162-bus 3 0.156
300-bus 4 0.281
Table 2.6: IEEE test system results without accounting Q-limits at PV-buses (NR).
IEEE system Iterations Execution time (s)
14-bus 2 0.047
30-bus 3 0.062
57-bus 3 0.062
118-bus NC -
300-bus 8 0.481
Table 2.7: IEEE test system results accounting Q-limits at PV-buses (NR).

Version-1.0
Chapter 3
Case Studies with A Test System:

NR method
3.1 The IEEE 14-bus System

The single line diagram of the IEEE 14-bus system is shown in Figure 3.1. The system
details are adopted from [13].
G GENERATORS
C SYNCHRONOUS
CONDENSERS
13
14
12
11
10
G
9
8
6
1
C C
7
5
4
Figure 3.1: The IEEE 14-bus system

Case Studies with A Test System: NR method Version-1.0
3.2 Newton-Raphson method

To run the load flow programme using Newton-Raphson method, the file required is:
loadflow.m. This file in turn calls the following .m files:
1. ybus_form.m: It constructs the YBU S
2. jacob_form.m: It constructs the Jacobian matrix and performs the solution of load
flow linearized equation till the solution is converged.
3. powerflow.m: It calculates the line flows and line losses.
4. lfl_result.m: It constructs the result files: lfl.dat and report.dat.
The above files require the following data files:
1. busno.dat : System details- number of lines, buses, transformers, etc.
2. nt.dat : Transmission line and transformer data.
3. pvpq.dat : Generation data and load data.
4. shunt.dat : Shunt data.
5. Qlim_data.dat : Reactive power generation limits data at PV-buses.
On successful run, it generates two output files: lfl.dat and report.dat. The
converged loadflow results are available in lfl.dat.
3.2.1 Format of Data Files

In the following lines the format of each of the data file has been given using the IEEE
14-bus system data:
System details:
File name: busno.dat

--------------------------------------------------------------------------
1 -----> Slack bus number.
0.001 -----> Loadflow convergence tolerance.
14 -----> Number of buses in the system.
17 -----> Number of lines.
3 -----> Number of transformers.
4 -----> Number of PV buses = (Number of generators - 1).
1 -----> To account Q-limits set this bit 1, otherwise 0. (Q-bit)

11 -----> Number of load buses (including loads at PV and slack buses).

1 -----> Number of shunts.
1.06 -----> Slack bus voltage magnitude.
---------------------------------------------------------------------------
Network data:
File name: nt.dat

-----------------------------------------------------------------------
From To R X y (total)/Tap ratio Remarks
-----------------------------------------------------------------------
1 2 0.01938 0.05917 0.0528 ---> Line 1
1 5 0.05403 0.22304 0.0492 ---> Line 2
2 3 0.04699 0.19797 0.0438 ---> Line 3
2 4 0.05811 0.17632 0.034 ---> Line 4
2 5 0.05695 0.17388 0.034 ---> Line 5
3 4 0.06701 0.17103 0.0128 ---> Line 6
4 5 0.01335 0.04211 0.00 ---> Line 7
7 8 0.0 0.17615 0.0 ---> Line 8
7 9 0.0 0.11001 0.0 ---> Line 9
9 10 0.03181 0.08450 0.0 ---> Line 11
6 11 0.09498 0.19890 0.0 ---> Line 11
6 12 0.12291 0.25581 0.0 ---> Line 12
6 13 0.06615 0.13027 0.0 ---> Line 13
9 14 0.12711 0.27038 0.0 ---> Line 14
10 11 0.08205 0.19207 0.0 ---> Line 15
12 13 0.22092 0.19988 0.0 ---> Line 16
13 14 0.17093 0.34802 0.0 ---> Line 17
4 7 0.0 0.20912 0.978 ---> Transformer data starts here.
4 9 0.0 0.55618 0.969
5 6 0.0 0.25202 0.932
------------------------------------------------------------------------

Generation and load data:
File name: pvpq.dat

--------------------------------------------------------------------------
Bus No. Vg/PL0 Pg0/QL0 Remarks
--------------------------------------------------------------------------
2 1.045 0.4 ---->Generator buses other than the slack bus
3 1.010 0 are specified as PV buses.
6 1.070 0
8 1.090 0
2 0.217 0.127 ----> Load data starts here (including loads at
3 0.942 0.19 and slack buses).
6 0.112 0.075
4 0.478 -0.039
5 0.076 0.016
9 0.295 0.166
10 0.090 0.058
11 0.035 0.018
12 0.061 0.016
13 0.135 0.058
14 0.149 0.050
--------------------------------------------------------------------------
Shunt admittances:
File name: shunt.dat

-----------------------
Bus No. G B
-----------------------
9 0 0.19
-----------------------

NOTE: The file Qlim_data.dat will be used if ‘Q-bit’= 1.

Q-Limits data at PV-buses.
File name: Qlim_data.dat

------------------------------------------
PV Bus No. Qmin Qmax
------------------------------------------
2 -0.4 0.5
3 0.0 0.4
6 -0.60 0.24
8 -0.60 0.24
------------------------------------------
The system has taken 2 iterations to reach convergence and the converged load flow
results accounting Q-limits at PV-buses are as follows:
File name: lfl.dat

----------------------------------------------------------------------------
Bus No VbO thetaO PGO QGO PLO QLO
----------------------------------------------------------------------------
1 1.060000 0.000000 2.323125 -0.165537 0.000000 0.000000
2 1.045000 -4.980701 0.400000 0.435235 0.217000 0.127000
3 1.010000 -12.721706 0.000000 0.250445 0.942000 0.190000
4 1.017716 -10.310331 0.000000 0.000000 0.478000 -0.039000
5 1.019556 -8.771159 0.000000 0.000000 0.076000 0.016000
6 1.070000 -14.216144 0.000000 0.126609 0.112000 0.075000
7 1.061581 -13.357956 0.000000 0.000000 0.000000 0.000000
8 1.090000 -13.357956 0.000000 0.175853 0.000000 0.000000
9 1.055992 -14.938240 0.000000 0.000000 0.295000 0.166000
10 1.051039 -15.096276 0.000000 0.000000 0.090000 0.058000
11 1.056948 -14.787257 0.000000 0.000000 0.035000 0.018000
12 1.055207 -15.071036 0.000000 0.000000 0.061000 0.016000
13 1.050404 -15.152100 0.000000 0.000000 0.135000 0.058000
14 1.035576 -16.032367 0.000000 0.000000 0.149000 0.050000
----------------------------------------------------------------------------

The partial content of file report.dat is shown below. The expressions used for the
calculation of line flows and the system losses are given in Appendix.
Line flows:
-----------------------------------------------------------------------------
Line flows Line flows
_____________________ _______________________
From To P-flow Q-flow From To P-flow Q-flow
-----------------------------------------------------------------------------
1 2 1.5683 -0.2039 2 1 -1.5253 0.2765
1 5 0.7549 0.0384 5 1 -0.7273 0.0224
2 3 0.7322 0.0356 3 2 -0.7090 0.0160
2 4 0.5612 -0.0157 4 2 -0.5444 0.0304
2 5 0.4150 0.0118 5 2 -0.4060 -0.0205
3 4 -0.2329 0.0445 4 3 0.2366 -0.0481
4 5 -0.6117 0.1583 5 4 0.6168 -0.1421
7 8 0.0000 -0.1713 8 7 0.0000 0.1759
7 9 0.2810 0.0578 9 7 -0.2810 -0.0498
9 10 0.0522 0.0423 10 9 -0.0520 -0.0420
6 11 0.0736 0.0354 11 6 -0.0730 -0.0342
6 12 0.0778 0.0250 12 6 -0.0771 -0.0235
6 13 0.1775 0.0720 13 6 -0.1754 -0.0678
9 14 0.0942 0.0362 14 9 -0.0931 -0.0337
10 11 -0.0380 -0.0160 11 10 0.0382 0.0163
12 13 0.0161 0.0075 13 12 -0.0161 -0.0075
13 14 0.0565 0.0174 14 13 -0.0560 -0.0162
4 7 0.2809 -0.0969 7 4 -0.2809 0.1139
4 9 0.1609 -0.0043 9 4 -0.1609 0.0174
5 6 0.4407 0.1249 6 5 -0.4407 -0.0807
------------------------------------------------------------------------------
Total real power losses in the system = 0.133859
Total reactive power losses in the system = 0.089733

Version-1.0
Chapter 4
Fast Decoupled Load Flow Method
4.1 Introduction to FDLF Method:

FDLF method is characterized by a decoupling of the linearized loadflow equations with
constant load flow matrices. This method is an alternative to the Newton method provided
the following three conditions are met:
1. The voltages are around their nominal values.
2. The angle difference across the lines is small.

X
3. The R
ratios are relatively large for all branches.
Fast decoupled load flow method is an approximate version of the Newton-Raphson

Method. The principle underlying the decoupled approach is based on the following
two standard observations in normal situations corresponding to a stable operating state
of a power system:
1. Change in the bus angle δ at a bus predominantly affects the real power flow in the
transmission lines and leaves the reactive power flow relatively unchanged.
2. Change in the voltage magnitude |Vp | at a bus primarily affects the flow of reac-
tive power in the transmission lines and leaves the flow of real power relatively
unchanged.
The development of FDLF method [7] from the Newton’s method has been explained
in the following sections.
4.2 Development of FDLF Equations:

Following the general procedure described earlier for the application of the Newton-
Raphson method in section 2.2, we have

Fast Decoupled Load Flow Method Version-1.0
∂P ∂P
 
" #  ∂δ ∂|V |  " #
△P   △δ
=  (4.1)
△Q 
 ∂Q ∂Q  △|V |

∂δ ∂|V |
∂P ∂Q
The first observation indicated above, essentially states that ∂δ
is much larger than ∂δ
∂Q ∂P
and the second observation states that ∂|V |
is much larger than ∂|V |
. Employing these
approximations in the Jacobian of (4.1), we get,
∂P
 
" #  ∂δ 0
" #
△P   △δ
=  (4.2)
△Q 
 ∂Q  △|V |

0
∂|V |
Expanding (4.2) we have
∂P1 ∂P1
    
△δ1 △P1
 ∂δ1 · · · ∂δn     
    
    
 .. .. ..  
    
= (4.3)
  
 . . . 
  ...   ... 
   

    
    
 ∂Pn ∂Pn     
··· △δn △Pn
∂δ1 ∂δn
 ∂Q1 ∂Q1   △V1   △Q 

|V1 | · · · |Vn | 1
 ∂|V1 | ∂|Vn |   |V1 |   
    
    
.. .. .. ..  
    
= .  (4.4)
  

 . . . 
 .   .. 
    
    
 ∂Qn ∂Qn   △Vn   
|V1 | · · · |Vn | △Qn
∂|V1 | ∂|Vn | |Vn |
The above equations are decoupled in the sense that the voltage-angle corrections △δ are
calculated using only real power mismatches △P , while the voltage magnitude correc-
tions are calculated using only △Q mismatches. However the coefficient matrices of (4.3)
and (4.4) are still interdependent. Because of this interdependency, it would still require
the evaluation and factoring of the two coefficient matrices at each iteration. To avoid
such computations, the elements of the coefficient matrices are simplified as follows [11]:

Modified off-diagonal elements: Rewriting the expression for off-diagonal elements of

H and L submatrices, from (2.16) and (2.23), we have
∂Pp ∂Qp
= = |Vp ||Vq |(Gpq sin δpq − Bpq cos δpq ) (4.5)
∂δq ∂|Vq |
|Vq |
Using the following assumptions in the above expression,
1. The angular differences (δp − δq ) between typical buses of the system are usually so
small that
cos(δp − δq ) = 1; sin(δp − δq ) = (δp − δq ) (4.6)
2. The susceptance Bpq is many times larger than the conductance Gpq so that
Gpq sin(δp − δq ) << Bpq cos(δp − δq ) (4.7)
we have,
∂Pp ∂Qp
= ∂|Vq | = −|Vp ||Vq |Bpq (4.8)
∂δq
|Vq |
Modified diagonal elements: Rewriting the expression for diagonal element of H

submatrix, from (2.18), we have
∂Pp
= −Qp − Bpp |Vp |2 (4.9)
∂δp
Assuming Qp << |V 2 |Bpp in (4.9) yields
∂Pp
= −Bpp |Vp |2 (4.10)
∂δp
From the expression for diagonal elements of L submatrix, (2.24), we have
∂Pp ∂Qp
= ∂|Vp | = −Bpp |Vp |2 (4.11)
∂δp
|Vp |
Using (4.8) and (4.11) in (4.3) and (4.4), and expanding, we get

    
−|V1 V1 |B11 −|V1 V2 |B12 ··· −|V1 Vn |B1n △δ1 △P1
    
    
    
 −|V V |B −|V2 V2 |B22 · · · −|V2 Vn |B2n   △δ   △P 
 2 1 21  2   2 
    
  =  (4.12)
.. .. .. ..   ..   .. 
    


 . . . .  .   . 
   
    
    
−|Vn V1 |Bn1 −|Vn V2 |Bn2 · · · −|Vn Vn |Bnn △δn △Pn
and
   △|V |   
−|V1 V1 |B11 −|V1 V2 |B12 ··· −V1 Vn |B1n 1 △Q1
  |V1 |    
    
    
 −|V V |B −|V2 V2 |B22 · · · −|V2 Vn |B2n   △Q 
 2 1 21   △|V2 |  
  2 
  
   |V | 
2 =
  (4.13)
.. .. .. .. ..  .. 
   
 

 . . . . 
 .   . 
 

   △|Vn |    
   
|Vn |

−|Vn V1 |Bn1 −|Vn V2 |Bn2 · · · −|Vn Vn |Bnn △Qn
(4.14)
Now, the first row of the coefficient matrix in (4.13) is multiplied by the correction
vector and then it is divided by |V1 | to obtain,
△Q1
−B11 △|V1| −B12 △|V2 | · · · −B1n △|Vn | = (4.15)
|V1 |
Each row of (4.13) can be similarly treated by representing the reactive power mismatch
△QP
at bus p by the quantity |VP |
, where, p = 1, 2, . . . , n. With these, (4.13) is modified as
    △Q 
−B11 −B12 · · · −B1n △|V1 | 1
    |V1 | 
    
    
 −B
21 −B22 · · · −B2n △|V2 | 
  
    △Q2 
    
   =  |V2 |  (4.16)
.. .. .. .. .. ..
    
    

 . . . . 
 .  
  . 

    △Qn 
    
|Vn |
−Bn1 −Bn2 · · · −Bnn △|Vn |

Using the first row of the coefficient matrix in (4.12) we have,
△P1
−|V1 |B11 △δ1 −|V2 |B12 △δ2 · · · −|Vn |B1n △δn = (4.17)
|V1 |
Applying the above steps to the remaining rows of the coefficient matrix, we get,
   △P
 
−|V1 |B11 −|V2 |B12 · · · −|Vn |B1n △δ1 1
    |V1 | 
    
    
 −|V |B
1 21 −|V2 |B22 · · · −|Vn |B2n   △δ   
  △P2
  2   
  
   =  |V2 |  (4.18)
.. .. .. ..   ..   ..
    
 
 . . . .  .   . 
    
    △Pn 
    
|Vn |
−|V1 |Bn1 −|V2 |Bn2 · · · −|Vn |Bnn △δn
The coefficient matrix in (4.18) is made the same as that in (4.16) by setting |V2 |,|V2 |, · · ·
|Vn | equal to 1.0 pu in the left hand side expression of (4.18). With these modification,
we get
   △P 
−B11 −B12 ··· −B1n △δ1 1
    |V1 | 
    
    
 −B
21 −B22 · · · −B2n   △δ   
  △P2
  2   
  
   =  |V2 |  (4.19)
.. .. .. ..   ..   ..
    
 
 . . . .  .   . 
    
    △Pn 
    
|Vn |
−Bn1 −Bn2 · · · −Bnn △δn
In FDLF method, the P − δ equation, (4.19) and Q − |V | equation, (4.16) are written
in a compact form as:
△P
= [B ′ ][△δ] (4.20)
|V |
△Q

= [B ′′ ][△|V |] (4.21)
|V |
Where B ′ and B ′′ are susceptances matrices and are obtained as negative imaginary part
of bus admittance matrix. The correction vectors △δ and △|V | are obtained by solving
(4.20) and (4.21) respectively. These corrections are used to update the angle and voltage
magnitude at bus p in k th iteration as follows:
δpK = δp(K−1) + △δp(K−1) (4.22)

|VpK | = |Vp(K−1) | + |△Vp(K−1) | (4.23)
4.2.1 Formation of B ′ and B ′′ Matrices:

Some of the typical formulations of B ′ − B ′′ matrices, in order, are B-B, B-X, X-B and
X-X type [16]. The letter B represents the inclusion of both the branch resistances and
reactances in the formation of coefficient matrix and the letter X indicates construction
of coefficient matrix by neglecting the branch resistances. In this paper, the standard X-
B-version is considered for the formulation of B ′ and B ′′ matrices because of its excellent
convergence properties under normal operating conditions [17].
Formation of B ′ Matrix:
Form Yd with the following assumptions:
1. Omit the representation of those network elements that predominantly affect MVAr
flows, i.e., shunt susceptances (inductive or capacitive) if any, line charging, and
off-nominal in-phase transformer taps.
2. Neglect series resistances.
Then
B ′ = −Im[Yd ]
Note that B ′ does not contain the row and column corresponding to slack bus and its
size is (n − 1) × (n − 1).
Formation of B ′′ Matrix:
Construct YBU S after neglecting the effect of phase shifters, if any and then
B ′′ = −Im[YBU S ]
Note that B ′′ does not contain rows and columns corresponding to slack and PV-buses
and its size is (n − m − 1) × (n − m − 1).
4.3 FDLF Implementation Algorithm:

There are two types of schemes available for the implementation of FDLF method, these
are as follows:
1. Selective iteration scheme [7]: In this approach, △P and △Q mismatches are verified
individually for convergence. This leaves the possibility to skip one or more P − δ
and/or Q − |V | iterations as soon as the related power mismatches are converged.

2. Successive iteration scheme [16, 18]: In this approach, both △P and △Q are verified
at a time for convergence, i.e., if △P convergence has reached early compared to
△Q convergence, even then P − δ cycle is repeated until the △Q convergence has
reached.
In the following sections, the FDLF schemes have been demonstrated with and without
accounting Q-limits at PV-buses, only for selective iteration scheme.
4.3.1 PV-buses Without Q-limits Specifications:

The implementation steps are depicted in flow chart shown in Figure 4.1. The flow chart
also shows the provision made for the inclusion of Q-limits at PV-buses and voltage- and
frequency-dependent load models in the normal FDLF method [19]. In this case set ‘Q-
bit’ = 0, ‘VM-bit’ = 0 and ‘FM-bit’= 0. This setting by-passes the load model routine
[(C)-(D)], and Q-limit imposing routine [(X)-(Y)].
The sample 4 bus system shown in Figure 2.2 is used to demonstrate the FDLF method
without accounting Q-limits at PV-buses.
1. Construct the bus admittance matrix YBU S as usual.

 
8.9852 − j44.8360 −3.8156 + j19.0781 −5.1696 + j25.8478 0
 −3.8156 + j19.0781 8.9852 − j44.8360 0 −5.1696 + j25.8478
 

 
 −5.1696 + j25.8478 0 8.1933 − j40.8638 −3.0237 + j15.1185 
 
0 −5.1696 + j25.8478 −3.0237 + j15.1185 8.1933 − j40.8638
 
 
Using this bus admittance matrix, B ′′ matrix is obtained by removing rows and
columns corresponding to PV-bus and slack bus as:
" #
44.8360 0
B ′′ =
0 40.8638
2. Construct the bus admittance matrix Yd , following the assumptions explained in

section (4.2.1).

0 0
Initialize V, δ and f=0
Set k =1
Compute Pp and Q p
Is
VM−bit =1 Yes
( or) C D
FM−bit = 1
?
No
Compute P and Q
P Yes Q Yes
k=k+1 convergence convergence
? ?
No
No
Solve (P −δ ) cycle and update δ
Compute Pp and Q p
Set the updated voltages as Output results

Is Yes
the new specified voltages Q−bit = 1 X Y
at the violated PV− buses ?
No
Compute P and Q
Q Yes P Yes
convergence convergence
Yes ? ?
No Is No No
Q−bit = 1
Solve (Q− | V | ) cycle and update | V |
?
Figure 4.1: Flow chart for the modified FDLF selective iteration scheme.
 
0 − j46.723 0 + j19.8413 0 + j26.8817 0
 0 + j19.8413 0 − j46.7230 0 0 + j26.8817
 

 
Yd = 
 0 + j26.8817 0 0 − j42.6050 0 + j15.7233 

0 0 + j26.8817 0 + j15.7233 0 − j42.6050
 
 
Using this bus admittance matrix, B ′ matrix is obtained by removing the row and
column corresponding to the slack bus as:
 
46.7230 0 −26.8817
B′ =  0 42.6050 −15.7233 
 
−26.8817 −15.7233 42.6050

Iteration 1:
1. The specified bus powers and the bus powers calculated (using the initial bus volt-
ages) are :
△Pp
Bus No. p Ppsp Pp |Vp |
2 -1.7 -0.1034 -1.5966
3 -2.0 -0.0605 -1.9395
4 2.38 0.1671 2.2129
2. Using the computed real power mismatches and the elements of B ′ matrix, we solve
(4.20), to get the angle correction vector as:
Bus No. p 2 3 4
△δp (deg) -1.1312 -2.0780 1.4369
3. Using the updated angles, the computed reactive power mismatches are:
△Qp
Bus No. p Qsp
p Qp |Vp |
2 1.0535 -0.2652 -0.7883
3 1.2394 0.0177 -1.2571
4. Using the computed reactive power mismatches and the elements of B ′′ matrix,
we solve (4.21), to get the voltage magnitude correction vector and the voltage
magnitude vector as:
Bus No. p 2 3
|△Vp | (pu) -0.0176 -0.0308
|Vp |(pu) 0.9824 0.9692
At the end of 3rd iteration (both P − δ and Q − |V | cycle):

The specified bus powers and the bus powers calculated (using the updated bus volt-
ages at the end of 3rd iteration) are:

p Qp ∆Qp
2 -1.70 -1.7000 0 -1.0535 -1.0535 0.000
3 -2.0 -2.0001 0.0001 -1.2394 -1.2394 0.000
4 2.38 2.3801 0.0001 - - -
From the above table, it is observed that the calculated bus power mismatches fall
below ǫ. The converged results are shown in Table 4.1.

Bus No. p |Vp | δp PGp QGp PLp QLp

No (pu) (deg) (pu) (pu) (pu) (pu)
1 1.000000 0.000000 1.868101 1.145057 0.500000 0.309900
2 0.982421 -0.976037 0.000000 0.000000 1.700000 1.053500
3 0.969003 -1.872244 0.000000 0.000000 2.000000 1.239400
4 1.020000 1.523196 3.180000 1.814320 0.800000 0.495800
Table 4.1: Converged load flow results without Q-limits (Selective iteration scheme).
4.3.2 PV-buses With Q-limits Specifications:

In FDLF method, the Q-limits at PV-buses are handled by a way of bus-type switching
(see section 1.2.2) as is implemented in the N-R method. Here, instead of iteration count,
the reactive power mismatch at PQ-buses (max|∆Q| ≤ 0.1) is used as the starting criteria
since it provides a direct way to control the degree of convergence [9]. It is observed that
the enforcement of Q-limit specifications at PV-buses by a way of bus-type switching,
causes an oscillatory convergence. To overcome these convergence difficulties, a procedure
which is a combination of the bus-type switching and adjusting the specified voltage at
PV-buses, is suggested to the standard FDLF. The adjustment procedure is as follows:
• Following Q-limit violations at PV-buses in a certain iterative step, the buses are
switched to PQ-type and the voltage corrections at those buses are calculated by
solving Q − |V | cycle.
• The obtained voltage corrections are used to adjust the bus voltages at the violated
PV-buses.
• The adjusted bus voltages are set as the new specified voltages at the PV-buses
for the next iteration. This implies that the back-offs at PV-buses are neglected.
This is justified based on the observation made in [9] that the number of back-
offs is normally small when Q-limits are enforced after the solution has moderately
converged.
With the above cited steps, Figure 2.3 is modified as shown in Figure 4.2 and the
modified scheme is enabled in Figure 4.1 by setting ‘Q-bit’ = 1. The variable, Qstart_cri
is set to 0.1 in the programme fdlf_jacob_form.m
The sample 4 bus system shown in Figure 2.2 is used to demonstrate the FDLF method
accounting Q-limits at PV-buses using selective iteration scheme. The ‘Q-bit’ is set to 1,

Is
power mismatch No
max Q < 0.1
?
Yes
Compute Q G = Q + Q L for all PV−buses
Is
there any Q − limit No
violation at any
PV− buses
?
Yes
Treat violated PV−buses as PQ−buses
with Q G = Q Gmax (or) Q Gmin , whichever limit is violated
Readjust the number of PV and PQ buses and accordingly modify[ B ]
sp
Recompute Q p
Figure 4.2: Flow chart for modified bus-type switching logic for FDLF.
and ‘VM-bit’ and ‘FM-bit’ are set to 0.
Iteration 1:
The bus voltage magnitudes and angles at the end of 1st P- and Q-iteration are as
follows:
Bus No. p 1 2 3 4
|Vp | 1.00 0.9824 0.9692 1.02
δp 0 -1.1312 -2.0780 1.4369
It is observed that the enforcement of Q-limits takes place at the beginning of 2nd Q-
iteration. Therefore, the steps involved in the above calculation are identical to that of
1st iteration in a case where Q-limits are not accounted at PV-buses (see section 4.3.1.1).
Iteration 2
1. The specified bus powers and the bus powers calculated (using the updated bus
voltages at the end of 1st P- and Q-iteration) are:

△Pp
Bus No. p Ppsp Pp |Vp |
2 -1.7 -1.7814 0.0829
3 -2.0 -2.1188 0.1225
4 2.38 2.4425 -0.0612
2. Using the computed bus power mismatches and the elements of B ′ we have (as per
(4.20),
    
0.0829 46.7230 0 −26.8817 △δ2
    
    
    
 0.1225  =  0 42.6050 −15.7233  (4.24)
  △δ3
 
   
    
    
−0.0612 −26.8817 −15.7233 42.6050 △δ4
3. Solving the above equation, we get the angle corrections and updated bus angles as:
Bus No p 2 3 4
△δp (deg) 0.1505 0.1961 0.0850
δp (deg ) -0.9807 -1.8819 1.5219
4. Recomputed reactive power mismatches △Q (using updated angles) are

△Qp
Bus No. p Qsp
p Qp |Vp |
2 -1.0536 -1.0530 -0.0006
3 -1.2397 -1.2289 -0.0108
4 - 1.3144 -
From the above tabulated values of △Qp , it is observed that the max|△Q| < 0.1.
Therefore, enforcement of Q-limits takes place now as per the flow chart given in
Figure 4.2. At bus 4, the real power load is PL4 = 0.8 pu. and reactive power load
is QL4 = 0.4958 pu. While checking the Q-limits at PV bus, the required reactive
power generation is calculated as: QG4 = Q4 + QL4 = 1.3144 + 0.4958 = 1.8102 pu.
Note that QG4 > QGmax (= 0.4 pu.). Therefore, this bus is treated as a PQ-bus with
QG4 = 0.4 pu. with the initial bus voltage magnitude as |V4 | = 1.02 pu. Due to this,
the number of PQ-buses are 3. Further, Qsp 4 = QG4 − QL4 = 0.4 - 0.4958 = -0.0958.
sp
∆Q4 = Q4 − Q4 = -0.0958 -1.3144 = 1.4102. And ∆Q 4
|V4 |
= 1.3826. Accordingly, the
B ′′ matrix is modified as:
 
44.8360 0 −25.8478
B ′′ =  0 40.8638 −15.1185 
 
−25.8478 −15.1185 40.8638

5. Using the computed bus power mismatches and the elements of B ′′ we have (as per
(4.21),
    
−0.0006 44.8360 0 −25.8478 △V2
 −0.0108  =  0 40.8638 −15.1185   △V3  (4.25)
    
−1.3826 −25.8478 −15.1185 40.8638 △V4
6. Solving the above equation, we get
Bus No. p 2 3 4
|△Vp | -0.0393 -0.0255 -0.0681
|Vp | 0.9432 0.9438 0.9519
7. Now, set |V4sp | = 0.9519 pu. for the next iterations.

At the end of 5th iteration: The specified bus powers and the bus powers calculated
(using the updated bus voltages at the end of 5th iteration) are:
p Qp ∆Qp
2 -1.7 -1.6998 0.0002 -1.0535 -1.0528 -0.0007
3 -2.0 -1.9998 0.0002 -1.2394 -1.2388 -0.0006
4 2.38 2.3796 0.0004 - -0.0963 -
From above table, it is observed that the bus power mismatches fall below ǫ. Further, it
is found that the Q-limit constraints at bus 4 are almost satisfied at this point (QG4 =
−0.0963 + 0.4958 = 0.3995). In this example, △Q mismatch reaches convergence at
the end of 4th iteration. And △P convergence has reached at the end of 5th iteration.
Converged load flow results are as follows:
Bus No.p |Vp | δp PGp QGp PLp QLp
1 1.000000 0.000000 1.884557 2.666903 0.500000 0.309900
2 0.941405 -0.533504 0.000000 0.000000 1.700000 1.053500
3 0.942331 -1.631538 0.000000 0.000000 2.000000 1.239400
4 0.951528 2.635645 3.180000 0.399517 0.800000 0.495800
Table 4.2: Converged load flow results with Q-limits (Selective iteration scheme)
4.3.3 Convergence Results for the IEEE Systems:

The results obtained using the IEEE test systems [13, 14] are presented in Tables 4.3 and
4.4. Base used is 100 MVA. The approximate execution time is obtained in a P4 computer
(3 GHz, 512 MB DDR 400 MHz memory, 910 GL chipset) and it also includes the time
for report generation.

IEEE system P-iter Q-iter. Execution time

(s)
14-bus 3 3 0.032
30-bus 3 3 0.032
57-bus 4 3 0.047
118-bus 4 3 0.078
145-bus 9 7 0.156
162-bus 4 3 0.109
300-bus 6 5 0.157
Table 4.3: IEEE test system results without accounting Q-limits at PV-buses (FDLF).
IEEE system P-iter Q-iter. Execution time

(s)
14-bus 3 3 0.032
30-bus 3 3 0.032
57-bus 4 3 0.047
118-bus 4 3 0.063
300-bus 6 5 0.156
Table 4.4: IEEE test system results accounting Q-limits at PV-buses (FDLF).
4.4 Handling of Static Load Models in Load Flows

To account load modelling in load flow programmes, polynomial expression-based static
load model has been employed, where voltage- and frequency-dependency of loads have
been considered [20]-[22]. The procedure involved has been discussed in the following
lines.
1. Voltage-dependent load models: Here the term PLp (i.e., real power load at bus p)
in the specified real power expression given by Ppsp = PGp − PLp , is replaced by the
polynomial expression for real power as:
( 2 )
|Vp | |Vp |
PLp = PLop a1p + a2p + a3p (4.26)
|Vop | |Vop |

Similarly, the reactive power load at bus p, QLp is given by

( 2 )
|Vp | |Vp |
QLp = QLop b1p + b2p + b3p (4.27)
|Vop | |Vop |
where, PLop and QLop are nominal values of active and reactive components of load
powers at nominal voltage |Vop |. a1 to a3 and b1 to b3 , are coefficients of voltage-
dependent real and reactive power loads at bus p with
a1p + a2p + a3p = 1

b1p + b2p + b3p = 1
2. Frequency-dependent load models: Here the term PLp (i.e., real power load at bus
p) is replaced by PLf p and is given as:

PLf p = PLp 1 + Kpf △f (4.28)
Similarly, the term QLp is replaced by QLf p and is given as:

QLf p = QLp 1 + Kqf △f (4.29)
where, Kpf and Kqf are coefficients of frequency-dependent real and reactive power
loads at bus p.
Incorporation of frequency-dependent load models in load flow studies can be im-
sp
plemented by knowing the specified real power generation at slack bus, PGs . The
implementation procedure is as follows:
(a) Compute the slack bus real power generation PGs = Ps + PLs , where, Ps repre-
sents the computed real power at slack bus, PLs represents the real power load
at slack bus if any.
(b) Initialize the per unit frequency deviation △f = 0
(c) Compute the frequency-deviation correction factor (△Fc ) using,
sp
△Fc = (PGs − PGs ) × Cf (4.30)
where Cf is referred to as a frequency convergence factor and is selected ap-

propriately for a given system.
(d) Update the frequency-deviation △f , in k th iterative step as
△f k = △f (k−1) + △Fck

(e) Use the updated △f, in (4.28 and 4.29).

(f) Compute the system frequency f = fo (1 + ∆f ), where fo is the nominal fre-
quency.
(g) While checking for convergence, the real power generation mismatch at slack
bus is also considered.
NOTE:
• The value of Cf is chosen so as to provide the minimum possible number of iterations.

This value is system dependent and is to be decided by trail and error method. To
get an approximate initial value of Cf , the following rule-of-thumb may be used:
 
sp
1  PGs
Cf = P 
10 PLop
all loads
• The governor characteristics are not modelled.
Thus, the specified bus power itself becomes a function of voltage magnitude and frequency
deviation and it varies from iteration to iteration. This in turn is used to calculate the
bus power mismatches in the solution process. The implementation flowchart for handling
the load models in load flows is given in Figure 4.3. ‘VM-bit’ and ‘FM-bit’ are set to 1 or
0, if voltage- and/or frequency-dependent load models are to be considered. Refer Figure
4.1 for the complete flowchart.
Test results are obtained are tabulated in Table 4.5 considering Q-limits at PV-buses
with voltage- and frequency-dependent load models for all loads. This is done by setting
‘Q-bit’ = 1, ‘VM-bit’ = 1 and ‘FM-bit’ = 1. The load model parameters used are a1 = 0.4,
a2 = 0.3 and a3 = 0.3 for real power loads and b1 = 0.4, b2 = 0.3, and b3 = 0.3 for reactive
power loads. In addition, Kpf and Kqf are set to 2 and -1.5, respectively. The nominal
frequency is 50 Hz.
sp
IEEE PGs Cf f P-itr Q-itr Time
system pu (Hz) s
14-bus 2.3715 0.1 50.0194 7 6 0.063
30-bus 2.6389 0.15 49.9700 8 7 0.078
118-bus 4.9048 0.01 49.9994 5 4 0.093
300-bus 4.5337 0.001 50.0003 14 11 0.281
Table 4.5: Study results with voltage- and frequency-dependent load models.
sp
The specified nominal real power generation at the slack bus (PGs ) is obtained for the
base case loads accounting voltage-dependent loads and Q-limits at PV-buses.

Is No
VM−bit =1
?
Yes
Read a1 , a2 , a3 and b1 , b2 , b3
sp sp
Compute the voltage dependent loads PL , Q L and P , Q
Is
No
FM−bit =1
?
Yes
sp
Read Cf , K pf , K qf and PGs at slack bus
Compute Fc and update f
sp sp
Compute the frequency dependent loads PLf , Q Lf and P , Q
Figure 4.3: Implementation of load modelling module in load flow.
NOTE:
The execution time indicated are not exact and is used for comparison purpose only.
These values are obtained using an in-built MATLAB function, etime and clock.
4.5 Representation HVDC Systems For Power Flow

Solution:
Power flow analysis requires joint solution of the DC and AC system of equations. The
sequential approach where solution is alternated between single AC and DC iterations is
employed [23]-[27]. In this method, the AC and DC link equations are solved separately
and thus the integration into the standard load flow program is carried out without
significant modifications of the AC load flow algorithm. The flow of various variables
between DC and AC solutions are illustrated in Figure 4.4.
Here, Eacr and Eaci are considered to be the input quantities for the solution of DC
system equations. They are known from the previous AC iteration. Using these bus

Interface
D.C system
Pdr
α Q dr
desired
Rectifier
α Eacr
min
A.C system
Transformer
Idc D.C line
bus
Eaci
γ
min Q di
Inverter
I − Pdi
m
Figure 4.4: Block diagram to illustrate the interfacing of DC system with AC system in
load flow.
voltages the converter bus powers, Pdr , Qdr , Pdi and Qdi are computed and are treated
as fictitious loads at the converter buses in the next iteration for solving the AC system
equations. This process is repeated until convergence is reached in AC iterations.
The basic converter equations, for both rectifier and inverter operations, describing
the relationship between the AC and DC variables are given below [12]:
√
3 2
Vdo = B T Eac
π
3
Vdr = Vdor cos α − Xcr Id Br
π
3
Vdi = Vdoi cos γ − Xci Id Bi
π
Φ = cos−1 (Vd /Vdo ) (4.31)
Pd = Vd Id
Qd = Pd tan Φ
where, Vdo is the ideal no-load direct voltage, α is the rectifier ignition delay angle, γ is the
inverter extinction advance angle, φ is the phase angle between the AC voltage Eac and the
fundamental AC current (i.e. power factor angle), Rc = π3 Xc is the commutating resistance
per bridge/phase, Xc is the leakage reactance of the converter transformer/phase, B is
the number of series connected bridges in each converter, T is the converter transformer
tap-ratio, Tn represents nominal tap setting. N denotes the number of tap positions, Rdc
represents resistance of the DC line and Im is the specified current margin.

4.5.1 DC Link Specification Schemes

The basic variables per converter are Vd , α(γ), T , φ, Pd and Qd , and the variable Id is
common to both converters. Thus in all, there are 13 variables available for any two-
terminal DC link. Two major specification schemes that are generally employed in DC
systems to solve for unknown variables, are as follows:
sp
1. DC link-current specification : In this case, the link-current, Idc is specified.
2. DC link-power specification : In this case, either rectifier-end DC power Pdr or

inverter-end DC power Pdi is specified.
In each of the above specifications, other variables such as α, Vdr , Tr at rectifier, and
γ, Vdi , Ti at inverter may be specified. Thus, it may lead to many combinations of the
specified variable set. The solution of DC link equations will be different for each of the
specification set [28]. In this paper, the following specification sets are used:
sp
1. Idc and Vdisp are specified.
sp
2. Pdr and Vdisp are specified.
The dependent and independent variables in the solution of DC equations depend on

the mode of operation of converters. For a given power direction, generally, 3 modes of
operations are defined: Mode-1: rectifier is on constant current (CC) control and inverter
is on constant extinction angle (CEA) control, Mode-2: rectifier is on constant ignition
angle (CIA) control and inverter is on constant current (CC) control, and Mode-3: rectifier
is on constant ignition angle (CIA) control and inverter is on constant extinction angle
(CEA) control. It is found that in power flow studies it is usually sufficient to consider
Mode-1 and Mode-2 [12]. For any given system conditions, the rectifier and inverter modes
of operation may not be known prior to the solutions of system equations. The following
procedure has been adopted to identify the mode of operation of a converter (see Figure
4.5):
1. Knowing AC side bus voltages, compute α using the following expression.

−1 Vdr + Rcr Br Id
α = cos (4.32)
Vodr
2. If α > αmin , Mode-1 condition is satisfied.
3. If α ≤ αmin , Mode-2 condition is satisfied.

Read converter bus voltages from AC load flow
sp
Set Id = Idc for current specified links.
C C
Compute I dc for power specified links and set I d = Idc
Compute α
α > α min Check for α < α min

mode of
operation
E ? R
Mode 1 calculations Mode 2 calculations

F T
Compute converter bus power

Pdc = [ Pdr , − Pdi ] and Q dc= [ Q dr , Q di ].
sp sp
Use Pdc and Qdc to modifiy P and Q
Figure 4.5: Flowchart of DC module.
4.5.2 Solution of DC Link Equations

DC link-current specification:
sp
Having specified Idc and Vdisp , it is possible to solve for the remaining unknown variables
with the knowledge of the AC voltages. The steps involved in Mode-1 and Mode-2 are as
follows:
1. Initially set Vdr = Vdrn , where Vdrn is the nominal rectifier-end DC voltage.
sp
2. Compute α using (4.32) with Id = Idc and identify the mode of operation.
3. Perform Mode-1 and/or Mode-2 calculations using (4.39) or (4.40) for each link.
Rectifier-end DC power specification:

sp
Having specified Pdr and Vdisp , it is possible to solve for the remaining unknown vari-
ables with the knowledge of the AC voltages. The steps involved are as follows:
sp
C
1. Compute the DC link current order Idc knowing Pdr and Vdisp . The equations used
C
in the computation of DC current Idc are :
Pdi = Vdisp Idc

C
(4.33)
sp C 2
Pdr = Pdi + (Idc ) Rdc (4.34)

C
Solving (4.33) and (4.34), we get a feasible expression for Idc as,

1
q
C sp sp 2 sp
Idc = −Vdi + (Vdi ) + 4 Pdr Rdc (4.35)
2Rdc
C
2. Initially set Vdr = Vdrn , and compute α using (4.32) with Id = Idc and identify the
mode of operation.
3. Perform Mode-1 or Mode-2 calculations using (4.39) or (4.40) for each link.
The expression that is used to calculate DC link-current order in Mode-2 is obtained as

follows:
We know that in Mode-2,
sp
Pdr = Vdr Id (4.36)
Vdr = Vdor cos αmin − Rcr Br Id (4.37)
Using (4.36) and (4.37), and solving for a feasible Id , we get
1 h p i
Id = Vdor cos αmin − (Vdor cos αmin )2 − 4Pdr Rcr (4.38)
2Rdc
1. Mode-1 calculation:
sp
Id = Idc for current specification.
C
Id = Idc for power specification, see (4.35).
V sp + Rci Bi Id
Ti = √ di
3 2
π Bi Eaci cos γmin
√
3 2
Vdoi = Bi Ti Eaci
π
φi = cos−1 (Vdisp /Vdoi ) (4.39)
Pdi = Vdisp Id
Qdi = Pdi tan φi
Vdr = Vdisp + Rdc Id
√
3 2
Vdor = Br Tr Eacr
π
−1 Vdr + Rcr Br Id
α = cos
Vodr
φr = cos−1 (Vdr /Vdor )
Pdr = Vdr Id
Qdr = Pdr tan φr

2. Mode-2 calculations:
sp
Id = Idc − Im for current specification.
Id is obtained using (4.38), for power specification.
√
3 2
Vdor = Br Tr Eacr
π
Vdr = Vdor cos αmin − Rcr Br Id
φr = cos−1 (Vdr /Vdor )
Pdr = Vdr Id for current specification
Qdr = Pdr tan φr
Vdi = Vdr − Rdc Id (4.40)
√
3 2
Vdoi = Bi Ti Eaci
π
−1 Vdi + Rci Bi Id
γ = cos
Vodi
φi = cos−1 (Vdi /Vdoi )
Pdi = Vdi Id
Qdi = Pdi tan φi
4.5.3 Incorporation of Control Variable Limits:

At each converter, the angle (α or γ) and the transformer tap (T ) can be controlled
directly to achieve (i) current control (ii) DC voltage control (iii) power control or (iv)
control of reactive power. The control strategies of tap change control are explained in
the following sections.
4.5.3.1 Rectifier and Inverter Transformers’ Tap Adjustments in Mode-1
The tap changer control at the rectifier is designed to maintain the delay angle, α within
the nominal range (say 10◦ to 20◦ ) in order to achieve certain voltage margin for the
purpose of current control. This has been implemented as follows:
1. If α violates the upper or lower limit of the nominal range, then Tr is adjusted as
Tr(k+1) = Tr(k) ± △Tr (4.41)
(Trmax − Trmin )
where, △Tr =
Nr
Trmax = specified rectifier transformer upper limit.

Trmin = specified rectifier transformer lower limit.

Nr = specified rectifier transformer number of tap positions.
2. If transformer tap-limits are reached, then the taps are set at the limiting values.
3. Compute Vdor and α.
When converters are operating in Mode-1, normally Vdisp is specified. This specified
value must correspond to the inverter tap setting which is within the specified range of
tap positions. However, if Vdisp specified leads to a tap setting out of the specified tap
range, then the transformer tap is set at the limiting values and the Vdisp is recomputed
as per (4.42).
√
3 2
Vdoi = Bi Ti Eaci
π
Vdi = Vdoi cos γmin − Rci Bi Id (4.42)
The above procedures are depicted in Figure 4.6.

E
Set γ = γ min and

sp
Vdi = Vdi for all mode 1 links
Compute Ti
Is
there any T i No
limit
violations
? Yes
Set Ti = Timax (or) Timin ,

Compute Vdi and Vdoi
compute α
Check α Yes
for the nominal
range
No ?
Re compute α Adjust T r
Is
No there any Tr
limit
violations
?
Yes
Set Tr= Trmax (or) Trmin ,
Re compute α
Compute Pdr, Qdrand Pdi , Qdi
Figure 4.6: Flowchart for Mode-1 calculations.

4.5.3.2 Rectifier and Inverter Transformers’ Tap Adjustments in Mode-2
In this mode, the rectifier transformer tap is adjusted to maximize DC voltage and the
inverter transformer tap is adjusted so that γ > γmin and var consumption is minimized.
In Figure 4.7, only inverter transformer tap-setting objective is met. The adjustment
procedure is illustrated below.
1. If γ < γmin then adjust the Ti as follows:
Ti(k+1) = Ti(k) − △Ti (4.43)
(Timax − Timin )
where, △Ti =
Ni
Timax = specified inverter transformer upper limit.
Timin = specified inverter transformer lower limit.
Ni = specified inverter transformer number of tap positions.
2. If transformer-tap-limits are reached, then the taps are set at the limiting values.
3. Compute γ and Vdoi .
sp
Set Id = I dc − Im for current
C
specified links. Compute Idc for C
power specified links and set I d= Idc .
Set α = α min for all mode 2 links.
Compute Vdr and Vdor
compute γ
Is
γ > γ min Yes
?
No
Re compute γ Adjust T i
Is
No there any T i
limit
violations
?
Yes
Set Ti = Timax (or) Timin ,
Re compute γ
Compute Pdr , Q drand Pdi, Qdi
Figure 4.7: Flowchart for Mode-2 calculations.
NOTE: If γ < γmin at the limiting values of Ti , it may be required to re-adjust γmin
and/or Timax /Timin specifications.

4.5.4 AC-DC Load Flow Algorithm:

After obtaining the rectifier and inverter bus powers Pdr , Qdr , Pdi and Qdi , from (4.39)
and (4.40), Pdc and Qdc are calculated as
Pdc = [Pdr , −Pdi ]T (4.44)

Qdc = [Qdr , Qdi ]T . (4.45)
Treating the above powers as fictitious loads, these quantities are interfaced with the
existing AC load flow as:
PL = PLo + Pdc (4.46)

QL = QLo + Qdc (4.47)
Where, PLo and QLo represent the vector of nominal loads in the AC system.
NOTE:
1. The DC system filter admittances Bsh are considered as shunts and are accounted
in a usual manner in the bus admittance matrix.
2. If voltage-dependent and/or frequency-dependent load modelling is to be considered,

then PLo and QLo are accordingly modified (see section 4.4).
These fictitious load powers are used to compute the bus power specifications P sp =
P G − P L and Qsp = QG − QL in the AC routine (see Figure 4.5). If bus power mismatches
are not below the prescribed tolerance value, i.e., if AC solutions are not converged, the
AC load flow routine computes a new set of AC bus voltages and are used in the next
DC calculations. The complete flowchart for the FDLF programme with the inclusion of
HVDC links and load models is given in Figure 4.8. In the figure, the module ‘AC Load
Flow’ represents the standard FDLF scheme (with a provision for inclusion of Q-limits at
PV-buses using an index ‘Q-bit’ as discussed in section 4.3.2). A flag, ‘DC-bit’ facilitates
the inclusion of DC links in the AC load flow.
A sample run for the IEEE 14-bus system with 2 two-terminal HVDC links and with
the following specification (refer Table 4.6), is shown in section 5.2.

Start
Read AC and DC system data.

Initialize AC and DC variables.
Set k = 1
Compute P and Q
Is
Yes
DC− bit = 1 V W
?
DC calculations
No
Load model depending on
k = k+1 ‘VM−bit’ and ‘FM−bit’ settings
Set the updated voltages as

the new specified voltage Compute P and Q
at the violated PV − buses
Check
for Yes
convergence
Is
( Q−bit = 1 &
? Compute line flows
DC−bit = 1) No & slack bus power
?
AC Load flow
End
Figure 4.8: Flowchart for the inclusion of DC calculations in AC load flow.
Parameter Link (5-4) Link (1-2)

Rectifier Inverter Rectifier Inverter
Bus No. 5 4 1 2
B 1 1 1 1
Xc 0.126 0.07275 0.1512 0.0873
αmin , γmin 5o 15o 5o 15o
T ±15% ±15% ±15% ±15%
N 27 19 25 20
Tn 1 1 1 1
Bsh 0.4902 0.6301 1.9608 1.0
Rdc 0.0034 0.0035
sp
Pdr 0.5865 1.5858
Vdisp 1.284 1.237
sp
Idc 0.4562 1.2774
Table 4.6: DC link data for the IEEE 14-bus system.

Version-1.0
Chapter 5
Case Studies with A Test System:

FDLF method
The single line diagram of the IEEE 14-bus system is shown in Figure 3.1. The system
details are adopted from [13]. To run the load flow programme using FDLF method, the
file required is: fdlf_loadflow.m. This file in turn calls the following .m files:
1. B_bus_form.m: It constructs B ′ and B ′′ matrices.
2. fdlf_jacob_form.m: It performs the solution of load flow equations till the solution
is converged.
3. dcflow.m: It performs the mode checks and calculates bus powers at the HVDC
link converter buses.
4. dcdata.m: It prepares the necessary data for the DC power flow using the input
data files.
5. powerflow.m: It calculates the line flows and line losses.
6. lfl_result.m: It constructs the result files: lfl.dat and report.dat.
The above files require the following data files:
1. busno.dat : System details- number of lines, buses, transformers, etc.
2. nt.dat : Transmission line and transformer data.
3. pvpq.dat : Generation data and load data.
4. shunt.dat : Shunt data.
5. Qlim_data.dat : Reactive power generation limits data at PV-buses.

Case Studies with A Test System: FDLF method Version-1.0
6. load_model.dat : Voltage-dependent load models data.
7. freq_model.dat : Frequency-dependent load models data.
8. hvdc.dat : DC link data -transformer taps, number of bridges per pole, α etc.
9. spdata.dat : DC link Specifications- current/power specification for each link.
On successful run, it generates three output files: lfl.dat, report.dat and hvdc_res.dat.
The converged loadflow results are available in lfl.dat. The final converged DC link re-
sults are given in hvdc_res.dat.
5.1 Format of Data Files

In the following lines the format of each of the data file has been given using the IEEE
14-bus system data:
System details:

--------------------------------------------------------------------------
17 -----> Number of lines. (existing ac lines are replaced by DC links)
0 -----> To account volt-dep. loads set this bit 1, otherwise 0. (VM-bit)
0 -----> To account freq-dep. loads set this bit 1, otherwise 0. (FM-bit)
0 -----> No of HVDC links (DC-bit).
----------------------------------------------------------------------------
NOTE: The formats of files: nt.dat, pvpq.dat, shunt.dat , lfl.dat and Qlim_data.dat
are identical to that given for NR-method -see section 3.2.1.

Voltage-dependent load models:
File name: load_model.dat

------------------------------------------------------
Load Bus No. a1 a2 a3 b1 b2 b3
------------------------------------------------------
2 0.4 0.3 0.3 0.4 0.3 0.3
3 0.4 0.3 0.3 0.4 0.3 0.3
6 0.4 0.3 0.3 0.4 0.3 0.3
2 0.4 0.3 0.3 0.4 0.3 0.3
3 0.4 0.3 0.3 0.4 0.3 0.3
6 0.4 0.3 0.3 0.4 0.3 0.3
4 0.4 0.3 0.3 0.4 0.3 0.3
5 0.4 0.3 0.3 0.4 0.3 0.3
9 0.4 0.3 0.3 0.4 0.3 0.3
10 0.4 0.3 0.3 0.4 0.3 0.3
11 0.4 0.3 0.3 0.4 0.3 0.3
12 0.4 0.3 0.3 0.4 0.3 0.3
13 0.4 0.3 0.3 0.4 0.3 0.3
14 0.4 0.3 0.3 0.4 0.3 0.3
------------------------------------------------------
Frequency-dependent load models:
File name: freq_model.dat

------------------------------------------------------------------
0.1 ----> Frequency convergence factor (Cf).
2.0 ----> The frequency coefficient for real power (Kpf).
-1.5 ----> The frequency coefficient for reactive power (Kqf).
2.371545 ----> The real power generation at slack bus. P^{sp}_{Gs}
------------------------------------------------------------------
NOTE:
1. The load_model.dat is used if ‘VM-bit’= 1. Constant power type load models are
assumed if there are no entries for loads (both real and reactive power) at a load
bus.
2. The freq_model.dat is used if ‘FM-bit’= 1.

5.2 IEEE 14-bus Example for Inclusion of HVDC

Link
In the original IEEE 14-bus system, two AC lines are replaced by HVDC links: one DC
sp
link between buses 5 and 4 (with Idc specification), and the second DC link between buses
sp
1 and 2 (with Pdr specification) -see Figure 5.1. In addition, Q-limits at PV-buses have
been accounted. The DC link details (in pu.) are given in Table 4.6 (adopted from [28]).
The base values used are: 100 MVA and 100 kV. See Appendix for base selection for DC
system.
G GENERATORS
C SYNCHRONOUS
CONDENSERS
13
14
12
11
10
G
9
8
6
1
C C
7
5
4
Figure 5.1: The modified IEEE 14-bus system.

The format of various input files pertaining to DC links is as follows:

HVDC Link data:
File name: hvdc.dat

------------------------------------------------------------------------------
Link No. Rec/Inv bus Br/Bi Tr_max/Ti_max Tr_min/Ti_min Tn_r/Tn_i
------------------------------------------------------------------------------
1 5 1 1.15 0.85 1.0
1 4 1 1.15 0.85 1.0
2 1 1 1.15 0.85 1.0
2 2 1 1.15 0.85 1.0
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Xcr/Xci Alpha_min/Gama_min Nr/Ni Bshr/Bshi
------------------------------------------------------------------------------
0.126 5 27 0.4902 ----> link 1 Rectifier data
0.07275 15 19 0.6301 ----> link 1 Inverter data
0.1512 5 25 1.9608 ----> link 2 Rectifier data
0.0873 15 20 1.0 ----> link 2 Inverter data
-----------------------------------------------------------------------------
HVDC Link Specifications-related data:
File name: spdata.dat

--------------------------------------------------------------------------
Link No. Rec-bus Inv-bus Vdis Idc Pdr Im Rdc
--------------------------------------------------------------------------
1 5 4 1.284 0.4562 0.5865 0.04 0.0034
2 1 2 1.237 1.2774 1.5858 0.12 0.0035
--------------------------------------------------------------------------
--------------------------------------------------------------------------
Alp_rmin Alp_rmax cn pn Vbase (kV) Vdrs(initial)
--------------------------------------------------------------------------
10 20 1 0 100 1 ----->Link 1 specified data
10 20 0 1 100 1 ----->Link 2 specified data
--------------------------------------------------------------------------
cn -set 1 if the dc link is on link current specification.

pn -set 1 if the dc link is on rectifier-end power specification.

For this case, where Q-limits are accounted, the files busno.dat and nt.dat are mod-
ified as follows:

--------------------------------------------------------------------------
15 * -----> Number of lines. (existing ac lines are replaced by DC links)
0 -----> To account volt-dep. loads set this bit 1, otherwise 0. (VM-bit)
0 -----> To account freq-dep. loads set this bit 1, otherwise 0. (FM-bit)
2 * -----> No of HVDC links (DC-bit).
----------------------------------------------------------------------------
* -indicates the changes made with respect to the original busno.dat.

Network data:
File name: nt.dat

-----------------------------------------------------------------------
From To R X y (total)/Tap ratio Remarks
-----------------------------------------------------------------------
%1 2 0.01938 0.05917 0.0528 ---> Line 1
1 5 0.05403 0.22304 0.0492 ---> Line 2
2 3 0.04699 0.19797 0.0438 ---> Line 3
2 4 0.05811 0.17632 0.034 ---> Line 4
2 5 0.05695 0.17388 0.034 ---> Line 5
3 4 0.06701 0.17103 0.0128 ---> Line 6
%4 5 0.01335 0.04211 0.00 ---> Line 7
7 8 0.0 0.17615 0.0 ---> Line 8
7 9 0.0 0.11001 0.0 ---> Line 9
9 10 0.03181 0.08450 0.0 ---> Line 11
6 11 0.09498 0.19890 0.0 ---> Line 11
6 12 0.12291 0.25581 0.0 ---> Line 12

6 13 0.06615 0.13027 0.0 ---> Line 13

9 14 0.12711 0.27038 0.0 ---> Line 14
10 11 0.08205 0.19207 0.0 ---> Line 15
12 13 0.22092 0.19988 0.0 ---> Line 16
13 14 0.17093 0.34802 0.0 ---> Line 17
4 7 0.0 0.20912 0.978 ---> Transformer data starts here.
4 9 0.0 0.55618 0.969
5 6 0.0 0.25202 0.932
--------------------------------------------------------------------------
NOTE:
Since two AC lines: 1-2 and 4-5 are replaced by two-terminal HVDC links, these two
lines must be deleted from the file nt.dat. Instead, one may comment out these two rows
as shown above.
With a tolerance factor, ǫ = 0.001 and without considering load models, the solution
has taken 4 P -iterations and 3 Q-iterations to reach convergence. The partial content of
file report.dat is shown below.
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Loadflow results:
----------------------------------------------------------------------------
Bus No VbO thetaO PGO QGO PLO QLO

----------------------------------------------------------------------------
1 1.060000 0.000000 2.286187 -1.250493 0.000000 0.000000
2 1.047204 -4.081198 0.400000 -0.397088 0.217000 0.127000
3 1.010581 -11.736848 0.000000 0.000133 0.942000 0.190000
4 1.058625 -10.249079 0.000000 0.000000 0.478000 -0.039000
5 1.029111 -8.156934 0.000000 0.000000 0.076000 0.016000
6 1.070000 -13.530646 0.000000 0.018110 0.112000 0.075000
7 1.079937 -13.131210 0.000000 0.000000 0.000000 0.000000
8 1.090000 -13.131210 0.000000 0.062267 0.000000 0.000000
9 1.073443 -14.659714 0.000000 0.000000 0.295000 0.166000
10 1.065463 -14.750492 0.000000 0.000000 0.090000 0.058000
11 1.064285 -14.285811 0.000000 0.000000 0.035000 0.018000
12 1.056551 -14.401321 0.000000 0.000000 0.061000 0.016000
13 1.052942 -14.527435 0.000000 0.000000 0.135000 0.058000
14 1.046708 -15.585947 0.000000 0.000000 0.149000 0.050000
-----------------------------------------------------------------------------

Line flows:
----------------------------------------------------------------------------
Line flows Line flows
_____________________ _______________________
From To P-flow Q-flow From To P-flow Q-flow
------------------------------------------------------------------------------
1 5 0.7004 -0.0010 5 1 -0.6768 0.0448
2 3 0.7284 0.0445 3 2 -0.7055 0.0058
2 4 0.6000 -0.2478 4 2 -0.5781 0.2765
2 5 0.4347 -0.0364 5 2 -0.4249 0.0297
3 4 -0.2365 -0.1956 4 3 0.2425 0.1973
7 8 0.0000 -0.0617 8 7 -0.0000 0.0623
7 9 0.2811 0.0675 9 7 -0.2811 -0.0596
9 10 0.0522 0.0817 10 9 -0.0520 -0.0810
6 11 0.0736 -0.0039 11 6 -0.0731 0.0048
6 12 0.0767 0.0199 12 6 -0.0760 -0.0185
6 13 0.1767 0.0517 13 6 -0.1747 -0.0478
9 14 0.0961 0.0615 14 9 -0.0947 -0.0585
10 11 -0.0381 0.0230 11 10 0.0382 -0.0227
12 13 0.0150 0.0025 13 12 -0.0150 -0.0025
13 14 0.0548 -0.0075 14 13 -0.0543 0.0085
4 7 0.2811 0.0200 7 4 -0.2811 -0.0058
4 9 0.1622 0.0437 9 4 -0.1622 -0.0305
5 6 0.4390 0.1704 6 5 -0.4390 -0.1246
-------------------------------------------------------------------------------
Total real power losses in the system = 0.096122
Total reactive power losses in the system = -4.533410
-------------------------------------------------------------------------------
NOTE:
The expressions used for the calculation of line flows and the system losses with HVDC
links, are given in Appendix.

Results of The D.C.Link: 1

----------------------------------------------------------
Parameter Rectifier Inverter

----------------------------------------------------------
Bus no 5 4
D.C.Voltage(pu) 1.285551 1.284000
Transformer tap Position(pu) 1.011111 0.952757
Control Angles(Deg) 17.466154 15.000000
Commutation overlap Angles(Deg) 11.397456 8.162801
Real Power flow(pu) 0.586468 0.585761
Reactive power consumption(pu) 0.258883 0.207395
Power factor 0.914833 0.942658
Current in the D.C.Link(pu) 0.456200
----------------------------------------------------------
Voltage Base in kV = 100.00
----------------------------------------------------------
Results of The D.C.Link: 2
----------------------------------------------------------
Parameter Rectifier Inverter

----------------------------------------------------------
Bus no 1 2
D.C.Voltage(pu) 1.241471 1.237000
Transformer tap Position(pu) 1.012000 0.983495
Control Angles(Deg) 10.173883 15.000000
Commutation overlap Angles(Deg) 32.968526 20.629222
Real Power flow(pu) 1.585800 1.580089
Reactive power consumption(pu) 0.953687 0.812281
Power factor 0.856966 0.889365
Current in the D.C.Link(pu) 1.277356
----------------------------------------------------------
Voltage Base in kV = 100.00
----------------------------------------------------------

converged D.C Load flow Results:
File name: hvdc_res.dat

---------------------------------------------------------------------------
Rec/ Vdr/ Tr/ Pdr/ Qdr/ Alpha/ P.f Vbase
Inv Vdi Ti Pdi Qdi Gama
---------------------------------------------------------------------------
5 1.285551 1.011111 0.586468 0.258883 17.466154 0.914833 100.00
1 1.241471 1.012000 1.585800 0.953687 10.173883 0.856966 100.00
4 1.284000 0.952757 0.585761 0.207395 15.000000 0.942658 100.00
2 1.237000 0.983495 1.580089 0.812281 15.000000 0.889365 100.00
----------------------------------------------------------------------------

Version-1.0
Appendix A
Important Expressions Used in Load

Flow Programme
Computation of Line Flows

Transmission network mainly consists of transmission lines and transformers. The mod-
elling of these components are briefly discussed in the following sections:
Transmission Lines:
Transmission Lines are modelled as a nominal π circuit [6] as shown in Figure A.1.
I I qp
pq
From To
Node Z
Node
(p) (q)
y y
2 2
Figure A.1: Nominal π Model of transmission lines.
where,
Z: represents the series impedance of the line.
y
2
: represents half of the total line charging y, at each node.

Version-1.0
Transformers:
The transformers are generally used as inter-connecting (IC) transformers and generator
transformers. These transformers are usually with off-nominal turns ratio and notations
used are shown in Figure A.2. The π equivalent circuit [6] of transformer is given in
Figure A.3.
a:1 yt
Ipq Iqp
p q
Figure A.2: Transformer with off-nominal tap.
I y t/a No tap
Tap pq I qp
(p) side
side
(q)
(1−a ) (a−1)
yt yt
a2 a
Figure A.3: Equivalent circuit of transformer.
where,
yt = z1t
zt : represents the series impedance at nominal-turns-ratio.
a: represents per unit off-nominal tap position.
1. For transmission lines: For the above model, the line currents Ipq and Iqp are com-
puted as
(Vp − Vq ) y
Ipq = + Vp
Z 2
(Vq − Vp ) y
Iqp = + Vq (A.1)
Z 2

Version-1.0
2. For transformers : In Figure A.3 assume that,
yt
Y1 =
a
(1 − a)
Y2 = yt
a2
(a − 1)
Y3 = yt
a
The currents Ipq and Iqp are computed as
Ipq = (Vp − Vq )Y1 + Vp Y2

Iqp = (Vq − Vp )Y1 + Vq Y3 (A.2)
The complex apparent powers Spq and Sqp are calculated as
∗
Spq = Ppq + jQpq = Vp Ipq
∗
Sqp = Pqp + jQqp = Vq Iqp
The net power loss in the component connected between nodes p and q is given by
Spq−loss = Spq + Sqp
If any shunt element (Yshp = Gp + jBp ) is present at node p, the shunt current and
hence the complex apparent power at node p are given by
Ishp = Vp Yshp
∗
Sshp = Vp Ishp
The total apparent power loss in the system is given by

X X X
Sloss = Spq−loss + Sshp = Id2 Rdc
all lines and transf ormers all shunts all dc links
(A.3)
n
X
= [PGp − PLp ] + j [QGp − (QLp + Qdcp )]
p=1
or in terms of bus power injection Sp at bus p, the Sloss is given by

n
X X
Sloss = Sp − (PLs + jQLs ) − (Pdcs + jQdcs ) + Id2 Rdc (A.4)
p=1 all dc links

Version-1.0
Computation of Overlap Angles:

The commutation overlap angles at rectifier and inverter bus are computed as,

−1 Vdr
µr = cos 2 − cos α − α (A.5)
Vdor

−1 Vdi
µi = cos 2 − cos γ − γ (A.6)
Vdoi
Base Selection for DC System

Base quantities for the DC system are chosen as follows:
SBDC = SBAC
VBDC = VBAC
(A.7)
SBDC
IBDC =
VBDC
VBDC
ZBDC =
IBDC
MATLAB Implementation:
1. The linearized load flow equation given by (2.14) for N-R method, and P − δ and
Q − |V | equations given by (4.20) and (4.21), respectively, for FDLF method, are
of the form Ax = b. In MATLAB, the solution for the correction vector is obtained
by using ‘backslash’ command as follows:
x = A\b
2. In the programming, all variables have been vectorized and declared as sparse to
exploit the advantage of sparsity solution techniques which is inherent to MATLAB.

Version-1.0
Bibliography
[1] IEEE Recommended practice for industrial and commercial power Systems Analy-
sis.IEEE Standards, 399 -1997
[2] G.T.Heydt, Computer Methods in Power System Analysis, Macmillan Publishing

Company, New York, 1986.
[3] Glenn W. Stagg and Ahmed H. El-Abiad, Computer Methods in Power System Anal-
ysis, McGraw-Hill Inc, New York, 1968.
[4] Elgerd O.I, An Introduction Electric Energy Systems Theory, Tata Mc Graw-Hill
Publishing Company Limited, Bombay, 1983.
[5] Saadat. H, Power System Analysis, McGraw-Hill Publishers, Singapore, 1999.
[6] M.A. Pai, Computer Techniques in Power System Analysis, Tata McGraw-Hill Pub-
lishing Company Ltd, New Delhi, 1979.
[7] Stott B and O. Alsac, “Fast Decoupled Load Flow method”, IEEE Trans on Power
Apparatus and Systems, Vol.PAS-93, pp. 859-869, 1974.
[8] W.F. Tinny and C.E. Hart, “Power flow solution by Newton’s method”, IEEE Trans.
on Power Apparatus and Systems, Vol.86, No. 8, pp. 1449-56, November, 1967.
[9] Show-kang Chang, and Vladimir Brandwajn, “Adjusted solutions in fast decoupled
load flow”, IEEE Trans. on Power Systems, Vol. 3, No. 2, pp. 726-733, May 1988.
[10] Stott. B, “Review of loadflow calculation methods”, Proceedings of the IEEE, Vol.62,
pp. 916-929, 1974.
[11] John J. Grainger and William D. Stevenson Jr, Power System Analysis, McGraw-Hill
Inc, New York, 1994.
[12] P.Kundur, Power System Stability and Contol, McGraw-Hill Inc, New York, 1994.
[13] Power system test cases archive, UWEE, www.ee.washington.edu/research/pstca.

Bibliography Version-1.0
[14] IEEE committe report “Transient stability test systems for direct stability methods”,
IEEE Trans on power systems Vol.7, pp. 37-43, Feb., 1992.
[15] Using MATLAB, Version 5.3, Release 11, The Math Works Inc.
[16] Robert A.M. Van Amerongen, “A general-purpose version of the Fast Decoupled
Load flow”, IEEE Trans. on Power Systems, Vol. 4,No.2, pp. 760-770, May, 1989.
[17] Hubbi,W., “Effects of of neglecting resistances in XB and BX load flow methods”,

IEE Proceedings, Vol.C, No.138, pp. 452-456, 1991.
[18] S.A Soman, S.A Khaparde and Shubha Pandit, Computational Methods for Large
Sparse Power Systems Analysis: An Object Oriented Approach, Kluwer Academic
Publishers., Boston/London, 2002.
[19] Kusic G.L, Computer Aided Power System Analysis, Prentice, New Delhi, 1989.
[20] K.R. Padiyar, Power System Dynamics - Stability and Control, BS Publications,
Hyderabad, India, 2002.
[21] P S R Murty, “Load modelling for power flow solutions”, J Inst Eng (India)
Electr.Eng Div., 58 pt EL3., pp. 162-165, December 1977.
[22] IEEE Committee Report,“Standard load models for power flow and dynamic perfor-
mance simulations IEEE Trans. on Power Systems, Vol. 10, No. 3, pp. 1302-1313,
August, 1995.
[23] Kimbark,E.W, Direct Current Transmission , Vol. 1, Wiley Limited, New york, 1971.
[24] K.R.Padiyar, HVDC Power Transmission Systems Technology and System Interac-
tions, Wiley Eastern Limited, New Delhi, 1990.
[25] Dusan Povh, Fellow, IEEE, “Use of HVDC and FACTS”, Proceedings of the IEEE.,
Vol. 88, No. 2, February 2000.
[26] T.J. Hammons, D. Woodford, J. Loughtan, M. Chamia, J. Donahoe, D. Povh, B.

Bisewski, and W. Imng, “Role of HVDC transmission in future energy development”,
IEEE Power Engineering Review, pp. 10-25, February, 2000.
[27] Willis Long and Stig Nilson, “HVDC transmission: yesterday and today”, IEEE
Power & Energy System magazine (Special issue on HVDC systems (pp. 32-61)),
Vol.5, No.2, pp. 22-31, march/ April 2007.
[28] Arrillaga.J, and Arnold C.P, Computer Analysis of Power Systems, John Wiley Lim-
ited, New Delhi, 1990.

Version-1.0
Acknowledgments
The authors thank Prof. A.M.Kulkarni at IIT Bombay for his valuable suggestions and
discussions regarding the programme implementation.
Please report bugs to: knsa1234@yahoo.com or knsa@nitk.ac.in

Manual For A Power Flow Programme: (Version 1.0)

Uploaded by

Document Informationclick to expand document information

Copyright:

Available Formats

Manual For A Power Flow Programme: (Version 1.0)

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Manual For A Power Flow Programme: (Version 1.0)

Uploaded by

Copyright:

Available Formats

Manual for

A Power Flow Programme

Dr. K.N. Shubhanga

Department of Electrical Engineering

List of Figures iii

1 Load Flow Studies: An Introduction 1

3 Case Studies with A Test System: NR method 27

NITK Surathkal i Electrical Dept.

4 Fast Decoupled Load Flow Method 33

5 Case Studies with A Test System: FDLF method 59

A Important Expressions Used in Load Flow Programme 69

NITK Surathkal ii Electrical Dept.

1.1 Representation of a typical bus p in load flow studies. . . . . . . . . . . . . 3

2.1 Flow chart for the implementation of Newton-Raphson method . . . . . . . 18

3.1 The IEEE 14-bus system . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

5.1 The modified IEEE 14-bus system. . . . . . . . . . . . . . . . . . . . . . . 62

A.1 Nominal π Model of transmission lines. . . . . . . . . . . . . . . . . . . . . 69

NITK Surathkal iii Electrical Dept.

2.1 Line data for the example. . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

NITK Surathkal iv Electrical Dept.

Load Flow Studies: An Introduction

1. System planning: This includes some of following tasks:

(a) Transmission planning: It involves design, analysis, sizing of conductors, trans-

NITK Surathkal 1 Electrical Dept.

(c) Studies that ensure power pool coordination.

1.2 Load flow Problem Formulation:

1. Formulation of a suitable mathematical model for the network.

2. Specifications of power and voltage constraints at various buses of the network.

NITK Surathkal 2 Electrical Dept.

Combining (1.1) and (1.2) yields

Figure 1.1: Representation of a typical bus p in load flow studies.

Then, the complex power injected at bus p is given by

where Qp represents inductive reactive power and a negative Qp represents capacitive

where n represents the number of buses in a system and Vq = |Vq |6 δq

NITK Surathkal 3 Electrical Dept.

Using (1.4) we can write an expression for apparent power Sp at bus p as

1.2.1 Types of Buses:

NITK Surathkal 4 Electrical Dept.

state (or unknown) variables.

2. Generator buses (PV-buses): A generator bus is usually called a voltage controlled

NITK Surathkal 5 Electrical Dept.

(a) Swing bus without load

1.2.2 Handling of Bus Types in Load Flow Problem:

1.2.2.1 Handling of PV-Buses:

NITK Surathkal 6 Electrical Dept.

1. Bus-type switching : In any iteration, if Q-limits are violated at a PV-bus, then it

1.2.2.2 Handling of Slack Bus:

SGs = (PGs + jQGs ) = Ss + (PLs + jQLs ) (1.9)

NITK Surathkal 7 Electrical Dept.

1.2.3 Load Flow Solution Methods:

NITK Surathkal 8 Electrical Dept.

2.1 Introduction To Newton-Raphson Method:

f1 (x01 + ∆x1 , x02 + ∆x2 , · · · , x0n + ∆xn ) = b1 (2.2)

NITK Surathkal 9 Electrical Dept.

form of the ith equation is

x1i = x0i + ∆xi

NITK Surathkal 10 Electrical Dept.

2.2 Application of the Newton-Raphson Method to