<<

. 5
( 17)



>>

for Use in 3-phase Load-¬‚ow Studies™, Proceedings of the IEE 121(10) 1155“1164.
Wedephol, L.M., 1963, ˜Application of Matrix Methods to the Solution of Travelling Wave Phenomena
in Poly-Phase Systems™, Proceedings of the IEE 110(12) 2200“2212.
Weedy, B.M., 1987, Electric Power Systems, John Wiley & Sons, Chichester.
4
Conventional Power Flow

4.1 INTRODUCTION

The main aim of a modern electrical power system is to satisfy continuously the electrical
power contracted by all customers. This is a problem of great engineering complexity where
the following operational policies must be observed: (1) nodal voltage magnitudes and
system frequency must be kept within narrow boundaries; (2) the alternating current (AC)
voltage and current waveforms must remain largely sinusoidal; (3) transmission lines must
be operated well below their thermal and stability limits; and (4) even short-term
interruptions must be kept to a minimum. Moreover, because of the very competitive nature
of the electricity supply business in an era of deregulation and open access, transmission
costs must be kept as low as possible.
To a large extent, several of these key issues in power system operation may be assessed
quite effectively by resorting to power ¬‚ow and derived studies (Arrillaga and Arnold, 1990;
Grainger and Stevenson, 1994; Stagg and El-Abiad, 1968; Wood and Wollenberg, 1984).
The main objective of a power ¬‚ow study is to determine the steady-state operating
condition of the electrical power network. The steady-state may be determined by ¬nding
out, for a given set of loading conditions, the ¬‚ow of active and reactive powers throughout
the network and the voltage magnitudes and phase angles at all buses of the network.
Expansion, planning and daily operation of power systems relies on extensive power ¬‚ow
studies (Kundur, 1994; Weedy, 1987). The information conveyed by such studies indicates
whether or not the nodal voltage magnitudes and active and reactive power ¬‚ows in
transmission lines and transformers are within prescribed operating limits. If voltage
magnitudes are outside bounds in one or more points of the network, then appropriate action
is taken in order to regulate such voltage magnitudes. Similarly, if the study predicts that the
power ¬‚ow in a given transmission line is beyond the power carrying capacity of the line,
then control action is taken.


4.2 GENERAL POWER FLOW CONCEPTS

The power ¬‚ow problem, is solved to determine the steady-state complex voltages at all
buses of the network, from which the active and reactive power ¬‚ows in every transmission


FACTS: Modelling and Simulation in Power Networks.
´ ´
Enrique Acha, Claudio R. Fuerte-Esquivel, Hugo Ambriz-Perez and Cesar Angeles-Camacho
# 2004 John Wiley & Sons, Ltd ISBN: 0-470-85271-2
94 CONVENTIONAL POWER FLOW

line and transformer are calculated (Stagg and El-Abiad, 1968). The set of equations
representing the power system are nonlinear. For most practical purposes, all power ¬‚ow
methods exploit the well-conformed nodal properties of the power network and equipment.
In its most basic form, these equations are derived by assuming that a perfect symmetry
exists between the phases of the three-phase power system (Arrillaga and Arnold, 1990).
Owing to the nonlinear nature of the power ¬‚ow equations, the numerical solution is reached
by iteration (Grainger and Stevenson, 1994).


4.2.1 Basic Formulation

A popular approach to assess the steady-state operation of a power system is to write
equations stipulating that at a given bus the generation, load, and powers exchanged through
the transmission elements connecting to the bus must add up to zero. This applies to both
active power and reactive power. These equations are termed ˜mismatch power equations™
and at bus k they take the following form:
ÁPk ¼ PGk À PLk À Pcal ¼ Psch À Pcal ¼ 0; °4:1Þ
k k k
ÁQk ¼ QGk À QLk À Qcal ¼ Qsch À Qcal ¼ 0: °4:2Þ
k k k

The terms ÁPk and ÁQk are the mismatch active and reactive powers at bus k, respectively.
PGk and QGk represent, respectively, the active and reactive powers injected by the generator
at bus k. For the purpose of the power ¬‚ow solutions it is assumed that these variables can be
controlled by the power plant operator. PLk and QLk represent the active and reactive powers
drawn by the load at bus k, respectively. Under normal operation the customer has control of
these variables, and in the power ¬‚ow formulation they are assumed to be known variables.
In principle, at least, the generation and the load at bus k may be measured by the electric
utility and, in the parlance of power system engineers, their net values are known as the
scheduled active and reactive powers:
Psch ¼ PGk À PLk ; °4:3Þ
k
Qsch ¼ QGk À QLk : °4:4Þ
k

The transmitted active and reactive powers, Pcal and Qcal , are functions of nodal voltages
k k
and network impedances and are computed using the power ¬‚ow equations. Provided the
nodal voltages throughout the power network are known to a good degree of accuracy then
the transmitted powers are easily and accurately calculated. In this situation, the
corresponding mismatch powers are zero for any practical purpose and the power balance
at each bus of the network is satis¬ed. However, if the nodal voltages are not known
precisely then the calculated transmitted powers will have only approximated values and the
corresponding mismatch powers are not zero. The power ¬‚ow solution takes the approach of
successively correcting the calculated nodal voltages and, hence, the calculated transmitted
powers until values accurate enough are arrived at, enabling the mismatch powers to be zero
or fairly close to zero. In modern power ¬‚ow computer programs, it is normal for all
mismatch equations to satisfy a tolerance as tight as 1e À 12 before the iterative solution can
be considered successful. Upon convergence, the nodal voltage magnitudes and angles yield
useful information about the steady-state operating conditions of the power system and are
known as state variables.
95
GENERAL POWER FLOW CONCEPTS

k m
zkm = zmk
Ik Im


Ek Em
Equivalent impedance
Figure 4.1



In order to develop suitable power ¬‚ow equations, it is necessary to ¬nd relationships
between injected bus currents and bus voltages. Based on Figure 4.1 the injected complex
current at bus k, denoted by Ik, may be expressed in terms of the complex bus voltages Ek
and Em as follows:

1
Ik ¼ °Ek À Em Þ ¼ ykm °Ek À Em Þ: °4:5Þ
zkm
Similarly for bus m,
1
Im ¼ °Em À Ek Þ ¼ ymk °Em À Ek Þ: °4:6Þ
zmk
The above equations can be written in matrix form as,
! ! !
ykm Àykm
Ik Ek
¼ ; °4:7Þ
Àymk ymk
Im Em
or
! ! !
Ik Ykk Ykm Ek
¼ ; °4:8Þ
Im Ymk Ymm Em

where the bus admittances and voltages can be expressed in more explicit form:

Yij ¼ Gij þ j Bij ; °4:9Þ
Ei ¼ Vi e ji ¼ Vi °cos i þ j sin i Þ; °4:10Þ
where i ¼ k; m; and j ¼ k; m.
The complex power injected at bus k consists of an active and a reactive component and
may be expressed as a function of the nodal voltage and the injected current at the bus:
Ã
Sk ¼ Pk þ j Qk ¼ Ek Ik
°4:11Þ
¼ Ek °Ykk Ek þ Ykm Em ÞÃ ;
Ã
where Ik is the complex conjugate of the current injected at bus k.
The expressions for Pcal and Qcal can be determined by substituting Equations (4.9) and
k k
(4.10) into Equation (4.11), and separating into real and imaginary parts:

Pcal ¼ Vk Gkk þ Vk Vm ½Gkm cos °k À m Þ þ Bkm sin °k À m ފ; °4:12Þ
2
k
Qcal ¼ ÀVk Bkk þ Vk Vm ½Gkm sin °k À m Þ À Bkm cos °k À m ފ: °4:13Þ
2
k
96 CONVENTIONAL POWER FLOW

For speci¬ed levels of power generation and power load at bus k, and according to
Equations (4.1) and (4.2), the mismatch equations may be written down as
È2 É
ÁPk ¼ PGk À PLk À Vk Gkk þ Vk Vm ½Gkm cos °k À m Þ þ Bkm sin °k À m ފ ¼ 0;
°4:14Þ
È2 É
ÁQk ¼ QGk À QLk À ÀVk Bkk þ Vk Vm ½Gkm sin °k À m Þ À Bkm cos °k À m ފ ¼ 0:
°4:15Þ

Similar equations may be obtained for bus m simply by exchanging subscripts k and m in
Equations (4.14) and (4.15).
It should be remarked that Equations (4.12) and (4.13) represent only the powers injected
at bus k through the ith transmission element, that is, Pik cal and Qik cal . However, a practical
power system will consist of many buses and many transmission elements. This calls for
Equations (4.12) and (4.13) to be expressed in more general terms, with the net power ¬‚ow
injected at bus k expressed as the summation of the powers ¬‚owing at each one of the
transmission elements terminating at this bus. This is illustrated in Figures 4.2(a) and 4.2(b)
for cases of active and reactive powers, respectively.


Qki cal
k k
m m
Pki cal
QGk
PGk
1 cal
P Qk cal
1
k


Pkn cal Qkn cal
QLk
PLk

(b)
(a)

Power balance at bus k: (a) active power, and (b) reactive power
Figure 4.2



The generic net active and reactive powers injected at bus k are:

X
n
¼ Pik cal ; °4:16Þ
Pcal
k
i¼1
Xn
Qcal ¼ Qik cal ; °4:17Þ
k
i¼1


where Pik cal and Qik cal are computed by using Equations (4.12) and (4.13), respectively.
As an extension, the generic power mismatch equations at bus k are:

X
n
ÁPk ¼ PGk À PLk À Pik cal ¼ 0; °4:18Þ
i¼1
X n
ÁQk ¼ QGk À QLk À Qik cal ¼ 0: °4:19Þ
i¼1
97
POWER FLOW SOLUTION METHODS

4.2.2 Variables and Bus Classi¬cation

In conventional power ¬‚ow theory each bus is described by four variables: net active power,
net reactive powers, voltage magnitude, and voltage phase angle.
Since there are only two equations per bus, two out of the four variables must be speci¬ed
in each bus in order to have a solvable problem. From a purely mathematical viewpoint, any
two variables can be speci¬ed; however, in engineering terms, the choice is based on which
variables at the bus can be physically controlled through the availability of a nearby
controller (Elgerd, 1982; Kundur, 1994; Weedy, 1987; Wood and Wollenberg, 1984). In the
broadest sense, one can think of voltage magnitudes and phase angles as state variables, and
active and reactive powers as control variables.
Buses are classi¬ed according to which two out of the four variables are speci¬ed:
 Load PQ bus: no generator is connected to the bus, hence the control variables PG and QG
are zero. Furthermore, the active and reactive powers drawn by the load PL and QL are
known from available measurements. In these types of buses the net active power and net
reactive power are speci¬ed, and V and  are computed.
 Generator PV bus: a generating source is connected to the bus; the nodal voltage
magnitude V is maintained at a constant value by adjusting the ¬eld current of the
generator and hence it generates or absorbes reactive power. Moreover, the generated
active power PG is also set at a speci¬ed value. The other two quantities  and QG are
computed. Constant voltage operation is possible only if the generator reactive power
design limits are not violated, that is, QG min < QG < QG max.
 Generator PQ bus: if the generator cannot provide the necessary reactive power support to
constrain the voltage magnitude at the speci¬ed value then the reactive power is ¬xed at
the violated limit and the voltage magnitude is freed. In this case, the generated active
power PG and reactive power QG are speci¬ed, and the nodal voltage magnitude V and
phase angle  are computed.
 Slack (swing) bus: one of the generator buses is chosen to be the slack bus where the
nodal voltage magnitude, Vslack, and phase angle, slack, are speci¬ed. There is only one
slack bus in the power system and the function of a slack generator is to produce suf¬cient
power to provide for any unmet system load and for system losses, that are not known in
advance of the power ¬‚ow calculation. The voltage phase angle at the slack bus slack is
chosen as the reference against which all other voltage phase angles in the system are
measured. It is normal to ¬x its value to zero.



4.3 POWER FLOW SOLUTION METHODS

4.3.1 Early Power Flow Algorithms

From the mathematical modelling point of view, a power ¬‚ow solution consists of solving
the set of nonlinear, algebraic equations that describe the electrical power network under
steady-state conditions. Over the years, several approaches have been put forward for the
solution of the power ¬‚ow equations. Early approaches were based on loop equations and
numerical methods using Gauss-type solutions. The method was laborious because the
network loops had to be speci¬ed beforehand by the systems engineer. Improved techniques
saw the introduction of nodal analysis in favour of loop analysis, leading to a considerable
98 CONVENTIONAL POWER FLOW

reduction in data preparation. Nevertheless, reliability towards convergence was still the
main concern. Further developments led to the introduction of the Gauss“Seidel method
with acceleration factors. The appeal of this generation of power ¬‚ow methods is their
minimum storage requirements and the fact that they are easy to comprehend and to code in
the form of computer programs. The drawback is that these algorithms exhibit poor
convergence characteristics when applied to the solution of networks of realistic size
(Elgerd, 1982). Power ¬‚ow solutions based on the nodal impedance matrix were brie¬‚y
experimented with (Brown, 1975), but problems with computer storage and speed became
insurmountable issues at the time. To overcome such limitations, the Newton“Raphson
method and derived formulations were developed in the early 1970s and have since become
¬rmly established throughout the power system industry (Peterson and Scott Meyer, 1971;
Stott, 1974; Stott and Alsac, 1978; Tinney and Hart, 1967).

4.3.2 The Newton“Raphson Algorithm

In large-scale power ¬‚ow studies the Newton“Raphson method has proved most successful
owing to its strong convergence characteristics (Peterson and Scott Meyer, 1971; Tinney and
Hart, 1967). This approach uses iteration to solve the following set of nonlinear algebraic
equations: 9
f 1 °x1 ; x2 ; Á Á Á ; xN Þ ¼ 0; >
>
>
f 2 °x1 ; x2 ; Á Á Á ; xN Þ ¼ 0; =
; or F°XÞ ¼ 0 °4:20Þ
.
. >
. >
>
;
f N °x1 ; x2 ; Á Á Á ; xN Þ ¼ 0;
where F represents the set of n nonlinear equations, and X is the vector of n unknown state
variables.
The essence of the method consists of determining the vector of state variables X by
performing a Taylor series expansion of F(X) about an initial estimate X(0):
    
°0Þ °0Þ °0Þ
F° X Þ ¼ F X þJ X XÀX þ higher-order terms; °4:21Þ
where J(X (0)) is a matrix of ¬rst-order partial derivatives of F(X) with respect to X, termed
the Jacobian, evaluated at X ¼ X (0).
This expansion lends itself to a suitable formulation for calculating the vector of state
variables X by assuming that X(1) is the value computed by the algorithm at iteration 1 and
that this value is suf¬ciently close to the initial estimate X(0). Based on this premise, all
high-order derivative terms in Equation (4.21) may be neglected. Hence,
3
2 2 °1Þ °0Þ 3
qf 1 °XÞ 
2 °1Þ 3 2 °0Þ 3 qf 1 °XÞ qf 1 °XÞ  X1 À X1
f 1 °X Þ f 1 °X Þ ÁÁÁ
qxn 7
6 qx1 qx2 6 7
7
6 7 6 7 6 6 7
7
6 7 6 7 6 qf °XÞ qf °XÞ 6 °1Þ 7
qf 2 °XÞ 7
6 7 6 7 62 6 X À X °0Þ 7
2
ÁÁÁ 7
6 f 2 °X°1Þ Þ 7 6 f 2 °X°0Þ Þ 7 6 62 27
qxn 7
7 þ 6 qx1 qx2
6 7%6 7:
6
7
6 7 6 7 6 6 7
7
6 7 6 7 6. . .. .
. . 6 7
. . . 7
. .
6 7 6 7 6. . . . 7 .
. . 6 7
4 5 4 5 6 .
4 5
.
5
4 qf n °XÞ qf n °XÞ .
. . qf n °XÞ 
°1Þ °0Þ
f n °X Þ f n °X Þ °1Þ °0Þ
 Xn À Xn
qx1 qx2 qxn
|¬„¬„¬„¬„¬„¬„¬„¬„{z¬„¬„¬„¬„¬„¬„¬„¬„} |¬„¬„¬„¬„¬„¬„¬„¬„{z¬„¬„¬„¬„¬„¬„¬„¬„} X¼X°0Þ |¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„{z¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„}
|¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„{z¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„}
F°X°1Þ Þ F°X°0Þ Þ X°1Þ ÀX°0Þ
J°X°0Þ Þ

°4:22Þ
99
POWER FLOW SOLUTION METHODS

In compact form, and generalising the above expression for the case of iteration (i),
      
°iÞ °iÀ1Þ °iÀ1Þ °iÞ °iÀ1Þ
%F X þJ X X ÀX ; °4:23Þ
FX

where i ¼ 1, 2 . .Á. . Furthermore, if it is assumed that X°iÞ is suf¬ciently close to the solution
À À Á
X(*) then F X°iÞ % F X°ÃÞ ¼ 0: Hence, Equation (4.23) becomes
    
°iÀ1Þ °iÀ1Þ °iÞ °iÀ1Þ
þJ X X ÀX ¼ 0; °4:24Þ
FX

and, solving for X°iÞ ,
  
X°iÞ ¼ X°iÀ1Þ À JÀ1 X°iÀ1Þ F X°iÀ1Þ °4:25Þ

The iterative solution can be expressed as a function of the correction vector
ÁX°iÞ ¼ X°iÞ À X°iÀ1Þ ,
  
°iÞ À1 °iÀ1Þ °iÀ1Þ
ÁX ¼ À J ; °4:26Þ
X FX

and the initial estimates are updated using the following relation:
X°iÞ ¼ X°iÀ1Þ þ ÁX °iÞ : °4:27Þ
The calculations are repeated as many times as required using the most up-to-date values of
X in equation (4.26). This is done until the mismatches ÁX are within a prescribed small
tolerance (i.e. 1e À 12).
In order to apply the Newton“Raphson method to the power ¬‚ow problem, the relevant
equations must be expressed in the form of Equation (4.26), where X represents the set of
unknown nodal voltage magnitudes and phase angles. The power mismatch equations ÁP
and ÁQ are expanded around a base point (h(0),V(0)) and, hence, the power ¬‚ow Newton“
Raphson algorithm is expressed by the following relationship:
2 3
qP qP V °iÞ " Áh #°iÞ
!°iÞ
ÁP qh qV 5
¼ À4 ÁV : °4:28Þ
qQ qQ V
ÁQ
|¬„¬„¬„¬„¬„{z¬„¬„¬„¬„¬„} qh qV V
|¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„{z¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„¬„} |¬„¬„¬„¬„¬„{z¬„¬„¬„¬„¬„}
F°X Þ
°iÀ1Þ
ÁX°iÞ
J°X°iÀ1Þ Þ

The various matrices in the Jacobian may consists of up to (nb À 1Þ ‚ °nb À 1Þ elements
of the form:
9
q Pk q Pk
Vm ; >
>
; =
q m q Vm
°4:29Þ
>
q Qk q Qk
Vm ; >
;
;
q m q Vm
where k ¼ 1; . . . ; nb; and m ¼ 1; . . . ; nb but omitting the slack bus entries.
Also, the rows and columns corresponding to reactive power and voltage magnitude for
PV buses are discarded. Furthermore, when buses k and m are not directly linked by a
transmission element, the corresponding k“m entry in the Jacobian is null. Owing to the low
degree of connectivity that prevails in practical power systems, the Jacobians of power ¬‚ows
are highly sparse. An additional characteristic is that they are symmetric in structure but not
in value (Zollenkoff, 1970).
100 CONVENTIONAL POWER FLOW

It must be pointed out that the correction terms ÁVm are divided by Vm to compensate for
the fact that Jacobian terms °qPk =qVm ÞVm and °qQk =qVm ÞVm are multiplied by Vm. It is
shown in the derivative terms given below that this arti¬ce yields useful simplifying
calculations.
Consider the lth element connected between buses k and m in Figure 4.1, for which self
and mutual Jacobian terms are given below.
For k 6¼ m:
qPk; l
¼ Vk Vm ½Gkm sin °k À m Þ À Bkm cos °k À m ފ; °4:30Þ
qm; l
qPk; l
Vm; l ¼ Vk Vm ½Gkm cos °k À m Þ þ Bkm sin °k À m ފ; °4:31Þ
qVm; l
qQk; l qPk; l
¼À Vm; l ; °4:32Þ
qm; l qVm; l
qQk; l qPk; l
Vm; l ¼ : °4:33Þ
qVm; l qm; l

For k ¼ m:
qPk; l
¼ ÀQcal À Vk Bkk ; °4:34Þ
2
qk; l k

qPk; l
Vk; l ¼ Pcal þ Vk Gkk ; °4:35Þ
2
qVk; l k

qQk; l
¼ Pcal À Vk Gkk ; °4:36Þ
2
qk; l k

qQk; l
Vk; l ¼ Qcal À Vk Bkk : °4:37Þ
2
qVk; l k


In general, for a bus k containing n transmission elements l, the bus self-elements take the
following form:
X qPk; l
qPk n
¼ ; °4:38Þ
qk qk; l
l¼1

X qPk; l
qPk n
Vk ¼ Vk; l ; °4:39Þ
qVk qVk; l
l¼1

X qQk; l
qQk n
¼ ; °4:40Þ
qk qk; l
l¼1

X qQk; l
qQk n
Vk ¼ Vk; l : °4:41Þ
qVk qVk; l
l¼1

The mutual elements given by Equations (4.30)“(4.33) remain the same whether we have
one transmission element or n transmission elements terminating at bus k.
After the voltage magnitudes and phase angles have been calculated by iteration, active
and reactive power ¬‚ows throughout the transmission system are determined quite
straightforwardly.
101
POWER FLOW SOLUTION METHODS

An important point to bear in mind is that the mismatch power equations ÁP and ÁQ of
the slack bus are not included in Equation (4.28) and the unknown variables Pslack and Qslack
are computed once the system power ¬‚ows and power losses have been determined. Also,
QG in PV buses are calculated in each iteration in order to check if the generators are within
reactive power limits. However, the mismatch reactive power equations ÁQ of PV buses are
not included in Equation (4.28). Details of this computation are given in the next section.
One of the main strengths of the Newton“Raphson method is its reliability towards
convergence. For most practical situations, and provided the state variables, X°0Þ , are
suitably initialised, the method is said to exhibit a quadratic convergence characteristic; for
example,
f°X °1Þ Þ ¼ 1e À 1;
f°X °2Þ Þ ¼ 1e À 2;
f°X °3Þ Þ ¼ 1e À 4;
f°X °4Þ Þ ¼ 1e À 8;
for the maximum mismatch. Contrary to non-Newton“Raphson solutions, such a
characteristic is independent of the size of the network being solved and the number and
kinds of control equipment present in the power system. Aspects that may dent its quadratic
convergence performance are reactive power limit violations in generator PV buses and
extreme loading conditions.


4.3.3 State Variable Initialisation

The effectiveness of the Newton“Raphson method to achieve feasible iterative solutions is
dependent upon the selection of suitable initial values for all the state variables involved in
the study.
The power ¬‚ow solution of networks that contain only conventional plat components is
normally started with voltage magnitudes of 1 p.u. (per unit) at all PQ buses. The slack and
PV buses are given their speci¬ed values, which remain constant throughout the iterative
solution if no generator reactive power limits are violated. The initial voltage phase angles
are selected to be 0 at all buses.


4.3.4 Generator Reactive Power Limits

Even though the mismatch reactive power equation ÁQk of PV bus k is not required in
Equation (4.28), solution of Equation (4.17) for the PV bus is still carried out at each
iterative step to assess whether or not the calculated reactive power Qcal is within the
k
generator reactive power limits:
QG min k < QG k < QG max k : °4:42Þ
If either of the following conditions occur during the iterative process:
)
Qk ! QG max k ;
cal
°4:43Þ
QG min k ;
Qcal
k
102 CONVENTIONAL POWER FLOW

bus k becomes a generator PQ bus with either of the following mismatch power equations
incorporated in Equation (4.28):

)
ÁQk ¼ QG max k À QL k À Qcal ;
k
°4:44Þ
ÁQk ¼ QG min k À QL k À Qcal ;
k



depending on the violated limit, together with the relevant Jacobian entries. The nodal
voltage magnitude at bus k is allowed to vary and Vk becomes a state variable.
It should be remarked that bus k may revert to being a generator PV bus at some point
during the iterative process if better estimates of Qcal , calculated with more accurate nodal
k
voltages, indicate that the reactive power requirements at bus k can, after all, be met by the
generator connected at bus k. Hence, reactive power limit checking is carried out at each
iteration. Programming wisdom indicates that limit checking should start after the ¬rst or
second iteration, since nodal voltage values computed at the beginning of the iterative
process may be quite inaccurate leading to misleading reactive power requirements. The
switching of buses from PV to PQ and vice versa impose additional numerical demands on
the iterative solution and retard convergence.



4.3.5 Linearised Frame of Reference

In order to illustrate how network components may be processed in the linearised frame of
reference afforded by the Newton“Raphson method (Fuerte-Esquivel et al., 1998) consider
the simple three-bus system shown in Figure 4.3. Bus 1 is selected to be the slack bus and
bus 2 is a generator bus. Bus 3 contains no generation and becomes a load bus. A
transformer and a transmission line link buses 1 and 2 and buses 2 and 3, respectively. One
shunt element and one load are connected at bus 3.
The concept of ˜power balance at a node™ may be used to great effect to account for bus
power injections in the Newton“Raphson solution. At a given bus, the power balance
is obtained by adding the contribution of each plant component connected to that bus.




Bus 2 Bus 3 PQ
Bus1
type

PV
Generator 1
Slack type Shunt Load
element
Generator 2

Figure 4.3 Three-bus network. Redrawn, with permission, from C.R. Fuerte-Esquivel, E. Acha,
S.G. Tan, and J.J. Rico, ˜Ef¬cient Object Oriented Power System Software for the Analysis of
Large-scale Networks Containing FACTS Controlled Branches™, IEEE Trans. Power Systems 3(2)
464“472, # 1998 IEEE
103
POWER FLOW SOLUTION METHODS

This is illustrated in Figure 4.4 with reference to Figure 4.3. The contribution of all three
buses is shown in this example for completeness, but it should be remembered that in
actual calculations active and reactive power mismatch entries are not required for the
slack bus. Likewise, the reactive power mismatch entry is not required for the generator
PV bus.


Mismatch
Mismatch
(specified power)
(calculated power)
P1 spec
P1 calc ’P 1
Ps gen 1 sen
P2 spec ’P 2
’Psen 2 P2 calc rec
Load
P3 spec
’Prec 3
Ps gen 2 ’Qsen1
’Ps load P3 calc
Generator Transformer
Q1 spec Line
’Qrec 2
Qs gen 1 ’Qsen 2 Q1 calc
Q2 spec
Qs gen 2 ’P 3
’ Q s load Q2 calc sen
’Qrec3
Q3 spec
Shunt
Q3 calc
’Qsen3

1 ∆P1
’Ptrans
1
Pgen 1
1
∆P2
2 ’Pline trans
’P
Pgen 2 2
2

=
3 ’Pline’P
+
3 ’Pload 3 ∆P3
shunt

∆Q1
1 ’Qtrans 1
Qgen 1
1
∆Q2
2 ’Qline’Qtrans
Qgen 2 2
2
3 ’Qline’Qshunt
’Qload ∆Q3
3
3
Mismatch
Mismatch Mismatch
(calculated power)
(specified power) vector

Figure 4.4 Power mismatch vector; subscripts ˜sen™ and ˜rec™ indicate the sending and receiving
ends. Reproduced, with permission, from C.R. Fuerte-Esquivel, E. Acha, S.G. Tan, and J.J. Rico,
˜Ef¬cient Object Oriented Power System Software for the Analysis of Large-scale Networks
Containing FACTS Controlled Branches™, IEEE Trans. Power Systems 3(2) 464“472, # 1998 IEEE




The construction of the Jacobian matrix is slightly more involved owing to the need to
evaluate self and mutual Jacobian terms, and ¬nding their location in the matrix.
Nevertheless, the basic procedure illustrated above, based on superposition, will also apply
to the formation of the Jacobian. For each plant component, relevant Jacobian equations are
chosen based on the type of buses to which the plant component is connected. These buses
determine the location of the individual Jacobian terms in the overall Jacobian structure.
The contributions of the line, transformer, and shunt components to the Jacobian are shown
in Figure 4.5. It should be noted that entries for the slack bus and the reactive power entry of
the generator bus are not considered in the Jacobian structure.
104 CONVENTIONAL POWER FLOW

COMPONENT
‚P2 ‚P2
‚θ 2 ‚θ 3
L L

‚P3 ‚P3
‚θ 2 ‚θ 3 JACOBIAN
L L


‚Q3
‚θ 3 3
2 3
L
Line
‚P2
‚V3 ‚P2 ‚P2
‚P2 ‚P2
L
+
2
‚P3 ‚θ3 ‚θ 3
‚θ 2 ‚θ 2 L L
L T
‚V3 L


‚Q3 ‚Q3
‚θ 2 ‚θ 3 ‚P3 ‚P3
‚P3 ‚P3 ‚P3
3 +
L L
‚θ 2 ‚θ 3 ‚θ 3 ‚V3 ‚V3
L S
L L S


‚P2
Transformer
‚θ 2 T

‚Q3 ‚Q3 ‚ Q3
‚Q3 ‚Q3
+
+
3
‚θ 2 ‚θ 3
‚θ 3 ‚θ 3 ‚V3
L L S
‚P3 L S
‚P3
‚θ3 ‚V3
S S
Shunt
‚Q3 ‚Q3
‚θ 3 ‚V3
S S



Figure 4.5 Jacobian structure. Reproduced, with permission, from C.R. Fuerte-Esquivel, E. Acha,
S.G. Tan, and J.J. Rico, ˜Ef¬cient Object Oriented Power System Software for the Analysis of
Large-scale Networks Containing FACTS Controlled Branches™, IEEE Trans. Power Systems 3(2)
464“472, # 1998 IEEE

4.3.6 Newton“Raphson Computer Program in Matlab1 Code

A computer program suitable for the power ¬‚ow solution of small and medium-sized power
systems is given in Program 4.1. The program is general, as far as the topology of the
network is concerned, and caters for any number of PV and PQ buses. Moreover, any bus in
the network may be designated to be the slack bus. Provisions are made for generator
reactive limit checking and to accommodate ¬x shunt compensation. No transformers are
represented in this base program and no sparsity techniques (Zollenkoff, 1970) are
incorporated.

PROGRAM 4.1 A program written in Matlab1 to calculate positive sequence power ¬‚ows
using the Newton“Raphson method
%***- - - Main Program

PowerFlowsData; %Read system data

[YR,YI] = YBus(tlsend,tlrec,tlresis,tlreac,tlsuscep,tlcond,shbus,...
shresis,shreac,ntl,nbb,nsh);
105
POWER FLOW SOLUTION METHODS

[VM,VA,it] = NewtonRaphson(nmax,tol,itmax,ngn,nld,nbb,bustype,...
genbus,loadbus,PGEN,QGEN,QMAX,QMIN,PLOAD,QLOAD,YR,YI,VM,VA);

[PQsend,PQrec,PQloss,PQbus] = PQ¬‚ows(nbb,ngn,ntl,nld,genbus,...
loadbus,tlsend,tlrec,tlresis,tlreac,tlcond,tlsuscep,PLOAD,...
QLOAD,VM,VA);

it %Iteration number
VM %Nodal voltage magnitude (p.u.)
VA = VA*180/pi %Nodal voltage phase angle(Deg)
PQsend %Sending active and reactive powers (p.u.)
PQrec %Receiving active and reactive powers (p.u.)

%End Main Program

%Build up admittance matrix
function [YR,YI] = YBus(tlsend,tlrec,tlresis,tlreac,tlsuscep,...
tlcond,shbus,shresis,shreac,ntl,nbb,nsh);
YR=zeros(nbb,nbb);
YI=zeros(nbb,nbb);
% Transmission lines contribution
for kk = 1: ntl
ii = tlsend(kk);
jj = tlrec(kk);
denom = tlresis(kk)^2+tlreac(kk)^2;
YR(ii,ii) = YR(ii,ii) + tlresis(kk)/denom + 0.5*tlcond(kk);
YI(ii,ii) = YI(ii,ii) - tlreac(kk)/denom + 0.5*tlsuscep(kk);
YR(ii,jj) = YR(ii,jj) - tlresis(kk)/denom;
YI(ii,jj) = YI(ii,jj) + tlreac(kk)/denom;
YR(jj,ii) = YR(jj,ii) - tlresis(kk)/denom;
YI(jj,ii) = YI(jj,ii) + tlreac(kk)/denom;
YR(jj,jj) = YR(jj,jj) + tlresis(kk)/denom + 0.5*tlcond(kk);
YI(jj,jj) = YI(jj,jj) - tlreac(kk)/denom + 0.5*tlsuscep(kk);
end
% Shunt elements contribution
for kk = 1: nsh
ii = shbus(kk);
denom = shresis(kk)^2+shreac(kk)^2;
YR(ii,ii) = YR(ii,ii) + shresis(kk)/denom;
YI(ii,ii) = YI(ii,ii) - shreac(kk)/denom;
end
% End of function YBus



%Carry out iterative solution using the Newton-Raphson method
function [VM,VA,it] = NewtonRaphson(nmax,tol,itmax,ngn,nld,nbb,...
bustype, genbus,loadbus,PGEN,QGEN,QMAX,QMIN,PLOAD,QLOAD,YR,YI,VM,VA)
106 CONVENTIONAL POWER FLOW

% GENERAL SETTINGS
D = zeros(1,nmax);
¬‚ag = 0;
it = 1;

% CALCULATE NET POWERS
[PNET,QNET] = NetPowers(nbb,ngn,nld,genbus,loadbus,PGEN,QGEN,...
PLOAD,QLOAD);

while ( it < itmax & ¬‚ag==0 )
% CALCULATED POWERS
[PCAL,QCAL] = CalculatedPowers(nbb,VM,VA,YR,YI);

% CHECK FOR POSSIBLE GENERATOR™S REACTIVE POWERS LIMITS VIOLATIONS
[QNET,bustype] = GeneratorsLimits(ngn,genbus,bustype,QGEN,QMAX,...
QMIN,QCAL,QNET, QLOAD, it, VM, nld, loadbus);

% POWER MISMATCHES
[DPQ,DP,DQ,¬‚ag] = PowerMismatches(nmax,nbb,tol,bustype,¬‚ag,PNET,...
QNET,PCAL,QCAL);

% JACOBIAN FORMATION
[JAC] = NewtonRaphsonJacobian(nmax,nbb,bustype,PCAL,QCAL,VM,VA,...
YR,YI);

% SOLVE FOR THE STATE VARIABLES VECTOR
D = JAC\DPQ™;

% UPDATE STATE VARIABLES
[VA,VM] = StateVariablesUpdates(nbb,D,VA,VM);


it = it + 1;
end
% End function Newton-Raphson


%Function to calculate the net scheduled powers
function [PNET,QNET] = NetPowers(nbb,ngn,nld,genbus,loadbus,PGEN,...
QGEN, PLOAD,QLOAD);
% CALCULATE NET POWERS
PNET = zeros(1,nbb);
QNET = zeros(1,nbb);
for ii = 1: ngn
PNET(genbus(ii)) = PNET(genbus(ii)) + PGEN(ii);
QNET(genbus(ii)) = QNET(genbus(ii)) + QGEN(ii);
end
for ii = 1: nld
PNET(loadbus(ii)) = PNET(loadbus(ii)) - PLOAD(ii);
107
POWER FLOW SOLUTION METHODS

QNET(loadbus(ii)) = QNET(loadbus(ii)) - QLOAD(ii);
end
%End function NetPowers

%Function to calculate injected bus powers
function [PCAL,QCAL] = CalculatedPowers(nbb,VM,VA,YR,YI)
% Include all entries
PCAL = zeros(1,nbb);
QCAL = zeros(1,nbb);
for ii = 1: nbb
PSUM = 0;
QSUM = 0;
for jj = 1: nbb
PSUM = PSUM + VM(ii)*VM(jj)*(YR(ii,jj)*cos(VA(ii)-VA(jj)) +...
YI(ii,jj)*sin(VA(ii)-VA(jj)));
QSUM = QSUM + VM(ii)*VM(jj)*(YR(ii,jj)*sin(VA(ii)-VA(jj)) “...
YI(ii,jj)*cos(VA(ii)-VA(jj)));
end
PCAL(ii) = PSUM;
QCAL(ii) = QSUM;
end
%End of functionCalculatePowers

%Function to check whether or not solution is within generators limits
function [QNET,bustype] = GeneratorsLimits(ngn,genbus,bustype,QGEN,...
QMAX,QMIN,QCAL,QNET, QLOAD, it, VM, nld, loadbus)
% CHECK FOR POSSIBLE GENERATOR™S REACTIVE POWERS LIMITS VIOLATIONS
if it > 2
¬‚ag2 = 0;
for ii = 1: ngn
jj = genbus(ii);
if (bustype(jj) == 2)
if ( QCAL(jj) > QMAX(ii) )
QNET(genbus(ii)) = QMAX(ii);
bustype(jj) = 3;
¬‚ag2 = 1;
elseif ( QCAL(jj) < QMIN(ii) )
QNET(genbus(ii)) = QMIN(ii);
bustype(jj) = 3;
¬‚ag2 = 1;
end
if ¬‚ag2 == 1
for ii = 1:nld
if loadbus(ii) == jj
QNET(loadbus(ii) = QNET(loadbus(ii)) - QLOAD(ii))
end
end
108 CONVENTIONAL POWER FLOW

end
end
end
%End function Generatorslimits


%Function to compute power mismatches
function [DPQ,DP,DQ,¬‚ag] = PowerMismatches(nmax,nbb,tol,bustype,...
¬‚ag,PNET,QNET,PCAL,QCAL);
% POWER MISMATCHES
DPQ = zeros(1,nmax);
DP = zeros(1,nbb);
DQ = zeros(1,nbb);
DP = PNET - PCAL;
DQ = QNET - QCAL;
% To remove the active and reactive powers contributions of the slack
% bus and reactive power of all PV buses
for ii = 1: nbb
if (bustype(ii) == 1 )
DP(ii) = 0;
DQ(ii) = 0;
elseif (bustype(ii) == 2 )
DQ(ii) = 0;
end
end
% Re-arrange mismatch entries
kk = 1;
for ii = 1: nbb
DPQ(kk) = DP(ii);
DPQ(kk+1) = DQ(ii);
kk = kk + 2;
end
% Check for convergence
for ii = 1: nbb*2
if ( abs(DPQ) < tol)
¬‚ag = 1;
end
end
%End function PowerMismatches


%Function to built the Jacobian matrix
function [JAC] = NewtonRaphsonJacobian(nmax,nbb,bustype,PCAL,QCAL,...
VM,VA,YR,YI);
% JACOBIAN FORMATION
% Include all entries
JAC = zeros(nmax,nmax);
iii = 1;
109
POWER FLOW SOLUTION METHODS

for ii = 1: nbb
jjj = 1;
for jj = 1: nbb
if ii == jj
JAC(iii,jjj) = -QCAL(ii) - VM(ii)^2*YI(ii,ii);
JAC(iii,jjj+1) = PCAL(ii) + VM(ii)^2*YR(ii,ii);
JAC(iii+1,jjj) = PCAL(ii) - VM(ii)^2*YR(ii,ii);
JAC(iii+1,jjj+1) = QCAL(ii) - VM(ii)^2*YI(ii,ii);
else
JAC(iii,jjj) = VM(ii)*VM(jj)*(YR(ii,jj)*sin(VA(ii)-VA(jj))...
-YI(ii,jj)*cos(VA(ii)-VA(jj)));
JAC(iii+1,jjj) = -VM(ii)*VM(jj)*(YI(ii,jj)*sin(VA(ii)...
-VA(jj))+YR(ii,jj)*cos(VA(ii)-VA(jj)));
JAC(iii,jjj+1) = -JAC(iii+1,jjj);
JAC(iii+1,jjj+1) = JAC(iii,jjj);
end
jjj = jjj + 2;
end
iii = iii + 2;
end
% Delete the voltage magnitude and phase angle equations of the slack
% bus and voltage magnitude equations corresponding to PV buses
for kk = 1: nbb
if (bustype(kk) == 1)
ii = kk*2-1;
for jj = 1: 2*nbb
if ii == jj
JAC(ii,ii) = 1;
else
JAC(ii,jj) = 0;
JAC(jj,ii) = 0;
end
end
end
if (bustype(kk) == 1) | (bustype(kk) == 2)
ii = kk*2;
for jj = 1: 2*nbb
if ii == jj
JAC(ii,ii) = 1;
else
JAC(ii,jj) = 0;
JAC(jj,ii) = 0;
end
end
end
end
%End of function NewtonRaphsonJacobian
110 CONVENTIONAL POWER FLOW

%Function to update state variables
function [VA,VM] = StateVariablesUpdates(nbb,D,VA,VM)
iii = 1;
for ii = 1: nbb
VA(ii) = VA(ii) + D(iii);
VM(ii) = VM(ii) + D(iii+1)*VM(ii);
iii = iii + 2;
end
%End function StateVariableUpdating



%Function to calculate the power ¬‚ows
function [PQsend,PQrec,PQloss,PQbus] = PQ¬‚ows(nbb,ngn,ntl,nld,...
genbus,loadbus,tlsend,tlrec,tlresis,tlreac,tlcond,tlsuscep,PLOAD,...
QLOAD,VM,VA);
PQsend = zeros(1,ntl);
PQrec = zeros(1,ntl);
% Calculate active and reactive powers at the sending and receiving
% ends of tranmsission lines
for ii = 1: ntl
Vsend = ( VM(tlsend(ii))*cos(VA(tlsend(ii))) + ...
VM(tlsend(ii))*sin(VA(tlsend(ii)))*i );
Vrec = ( VM(tlrec(ii))*cos(VA(tlrec(ii))) + ...
VM(tlrec(ii))*sin(VA(tlrec(ii)))*i );
tlimped = tlresis(ii) + tlreac(ii)*i;
current =(Vsend - Vrec) / tlimped + Vsend*( tlcond(ii) + ...
tlsuscep(ii)*i )*0.5 ;
PQsend(ii) = Vsend*conj(current);
current =(Vrec - Vsend) / tlimped + Vrec*( tlcond(ii) + ...
tlsuscep(ii)*i )*0.5 ;
PQrec(ii) = Vrec*conj(current);
PQloss(ii) = PQsend(ii) + PQrec(ii);
end
% Calculate active and reactive powers injections at buses
PQbus = zeros(1,nbb);
for ii = 1: ntl
PQbus(tlsend(ii)) = PQbus(tlsend(ii)) + PQsend(ii);
PQbus(tlrec(ii)) = PQbus(tlrec(ii)) + PQrec(ii);
end



% Make corrections at generator buses, where there is load, in order to
% get correct generators contributions
for ii = 1: nld
jj = loadbus(ii);
for kk = 1: ngn
ll = genbus(kk);
111
POWER FLOW SOLUTION METHODS

if jj == ll
PQbus(jj) = PQbus(jj) + ( PLOAD(ii) + QLOAD(ii)*i );
end
end
end
%End function PQ¬‚ows




4.3.7 The Fast Decoupled Algorithm

It was demonstrated in the late 1970s that the storage and computing requirements of the
Newton“Raphson method could be reduced very signi¬cantly by introducing a series of
well-substantiated, simplifying, assumptions in Equation (4.28). These assumptions are
based on physical properties exhibited by electrical power systems, in particular in high-
voltage transmission systems.
The resulting formulation is no longer a Newton“Raphson method but a derived
formulation described as ˜fast decoupled™ (Stott, 1974; Stott and Alsac, 1978). The power
mismatch equations of both methods are identical but their Jacobians are quite different; the
Jacobian elements of the Newton“Raphson method are voltage-dependent whereas those of
the fast decoupled method are voltage-independent (i.e. constant parameters). Moreover, the
number of Jacobian entries used in the fast decoupled method is only half of those used in
the Newton“Raphson method. The trade-off lies in the weakening of the strong convergence
characteristic exhibited by the Newton“Raphson method; the convergence characteristics of
the fast decoupled method are linear as opposed to quadratic.
For a typical power ¬‚ow problem, where convergence to a tight tolerance is mandatory,
the 6 iterations normally taken by the Newton“Raphson method to converge will increase to
25 iterations and above when the fast decoupled method is employed. In fact, the number of
iterations taken by the fast decoupled method may be in¬‚uenced by the size of the system
being solved, how loaded the system is, the number of power system controllers, and the
ratio of resistance to reactance in the transmission elements “ although there are simple
programming arti¬ces to circumvent this problem. However, an asset of the fast decoupled
method is that one of its iterations only takes a fraction of the time required by one of
the Newton“Raphson method iterations. Hence, in power ¬‚ow studies of high-voltage
networks with a small number of system controllers, the use of the fast decoupled method
may be advantageous.
With a view to developing the fast decoupled formulation, simpli¬cations will be
introduced into the Jacobian of Equation (4.28). It has been observed that, during normal
operation, incremental changes in voltage magnitude produce almost no change in active
power ¬‚ow and that, likewise, incremental changes in voltage phase angle produce almost
no change in reactive power ¬‚ow. Hence, the following Jacobian elements may be assumed
to be zero:


! !
qP qQ
V; : °4:45Þ
and
qV qh
112 CONVENTIONAL POWER FLOW

Accordingly, the overall problem stated in Equation (4.28) reduces to the following two
subproblems:
!
qP °iÞ
°iÞ
½ÁhŠ°iÞ ;
½ÁPŠ ¼ À °4:46Þ
|¬„¬„¬„{z¬„¬„¬„} |¬„¬„¬„{z¬„¬„¬„}
qh
|¬„¬„¬„{z¬„¬„¬„} ÁX°iÞ
F1 °X1 Þ
°iÀ1Þ

J1 °X1 Þ
°iÀ1Þ 1



! !
qQ °iÞ ÁV °iÞ
°iÞ
½ÁQŠ ¼ À : °4:47Þ
V
|¬„¬„¬„{z¬„¬„¬„} qV V
|¬„¬„¬„¬„¬„¬„{z¬„¬„¬„¬„¬„¬„} |¬„¬„¬„¬„{z¬„¬„¬„¬„}
F2 °X2 Þ
°iÀ1Þ

J2 °X2 Þ
°iÀ1Þ °iÞ
ÁX2

Further simplifying assumptions pertaining to high-voltage transmission networks, which
are relevant to the problem at hand, are as follows:
 X ) R in all transmission lines and transformers of the network.
 The difference in voltage phase angles between two adjacent buses is small and hence the
following relations apply: sin°k À m Þ ¼ k À m , and cos°k À m Þ ¼ 1:
 The nodal voltage magnitudes are close to 1 p.u. at every bus.
 Current ¬‚ows in shunt-connected elements may be grouped together with the equivalent
loads and generator currents.
Incorporating these assumptions in the Jacobian elements of Equations (4.46) and (4.47)
we obtain the following set of equations:
½ÁPŠ°iÞ ¼ À ½B0 Š ½ÁhŠ°iÞ ; °4:48Þ
|¬„¬„¬„{z¬„¬„¬„} |¬„¬„¬„{z¬„¬„¬„}
F1 °X1 Þ
°iÀ1Þ °iÞ
ÁX1

½ÁQŠ°iÞ ¼ À ½B00 Š ½ÁVŠ°iÞ ; °4:49Þ
|¬„¬„¬„{z¬„¬„¬„} |¬„¬„¬„{z¬„¬„¬„}
F 2 °X 2 Þ
°iÀ1Þ °iÞ
ÁX2

where B0 corresponds almost exactly to the negative of the imaginary part of the nodal
admittance matrix. Owing to the requirements of the power ¬‚ow problem, the row and
column corresponding to the slack bus is not included in B0 . Matrices B0 and B00 are identical
if no generator buses exist in the system. However, in the more general case, when generator
buses do exist in the system then the row and column corresponding to each generator bus
are removed from matrix B00 .
Equations (4.48) and (4.49) are very simple compared with that of the full Newton“
Raphson method given by Equation (4.28). Matrices B0 and B00 are symmetric in structure
and, provided no phase-shifting transformers are present in the system, they are also
symmetric in value. These matrices are inverted only once, during the ¬rst iteration, and
then remain constant throughout the iterative process. This is in contrast to with the
Newton“Raphson method, where the Jacobian is evaluated and inverted (factorised;
Zollenkoff, 1970) at each iteration.

4.3.8 Fast Decoupled Computer Program in Matlab1 Code

Program 4.2 is fully equivalent to the Newton“Raphson power ¬‚ow program given in
Section 4.3.6 (Programe 4.1). The functions PowerFlowsData, YBus and PQ¬‚ows are also
113
POWER FLOW SOLUTION METHODS

used here. The function FastDecoupled replaces NewtonRaphson, with the new function
using all the functions called by NewtonRaphson except for NewtonRaphsonJacobian,
which is replaced with FastDecoupledJacobian.

PROGRAM 4.2 Program written in Matlab1 to calculate positive sequence power ¬‚ows
by means of the fast decoupled method.

%- - - Main Program


PowerFlowsData; %Function to read data


% Form the bus admittance matrix
[YR,YI] = YBus(tlsend,tlrec,tlresis,tlreac,tlsuscep,tlcond,shbus,...
shresis,shreac,ntl,nbb,nsh);


[VM,VA,it] = FastDecoupled(nmax,tol,itmax,ngn,nld,nbb,bustype,...
genbus,loadbus,PGEN,QGEN,QMAX,QMIN,PLOAD,QLOAD,YR,YI,VM,VA);

[PQsend,PQrec,PQloss,PQbus] = PQ¬‚ows(nbb,ngn,ntl,nld,genbus,...
loadbus,tlsend,tlrec,tlresis,tlreac,tlcond,tlsuscep,PLOAD,QLOAD,...
VM,VA);

it %Iteration number
VM %Nodal voltage magnitude (p.u.)
VA = VA*180/pi %Nodal voltage phase angle(Deg)
PQsend %Sending active and reactive powers (p.u.)
PQrec %Receiving active and reactive powers (p.u.)

% End of Main Program

% Fast Decoupled function
function [VM,VA,it] = FastDecoupled(nmax,tol,itmax,ngn,nld,nbb,...
bustype,genbus,loadbus,PGEN,QGEN,QMAX,QMIN,PLOAD,QLOAD,YR,YI,VM,VA);

% GENERAL SETTINGS
¬‚ag = 0;
B1 = zeros(nbb,nbb);
B2 = zeros(nbb,nbb);
% CALCULATE NET POWERS
[PNET,QNET] = NetPowers(nbb,ngn,nld,genbus,loadbus,PGEN,QGEN,PLOAD,...
QLOAD);


% BEGINNING OF ITERATIVE LOOP


it = 1;
while ( it < itmax & ¬‚ag==0 )
114 CONVENTIONAL POWER FLOW

% CALCULATED POWERS
[PCAL,QCAL] = CalculatedPowers(nbb,VM,VA,YR,YI);

% CHECK FOR POSSIBLE GENERATOR™S REACTIVE POWERS LIMITS VIOLATIONS
[QNET,bustype] = GeneratorsLimits(ngn,genbus,bustype,QGEN,QMAX,...
QMIN,QCAL,QNET, QLOAD, it, VM, nld, loadbus);

% POWER MISMATCHES
[DPQ,DP,DQ,¬‚ag] = PowerMismatches(nmax,nbb,tol,bustype,¬‚ag,PNET,...
QNET,PCAL,QCAL);

% OBTAIN INVERTED JACOBIANS DURING THE FIRST ITERATION
[B1,B2] = FastDecoupledJacobian(nbb,bustype,DP,DQ,YI,B1,B2,it);
% SOLVE FOR THE STATE VARIABLES VECTOR
DVA = B1*DP™;
DVM = B2*DQ™;
% Re-arrange state variables entries
kk = 1;
for ii = 1: nbb
D(kk) = DVA(ii);
D(kk+1) = DVM(ii);
kk = kk + 2;
end

% UPDATE THE STATE VARIABLES VALUES
[VA,VM] = StateVariablesUpdating(nbb,D,VA,VM,it);

it = it + 1;
end
% End of function FastDecoupled


%Form the Jacobian for the Fast Decoupled Method
function [B1,B2] = FastDecoupledJacobian(nbb,bustype,DP,DQ,YI,...
B1,B2,it);
DVA = zeros(nbb);
DVM = zeros(nbb);
if ( it == 1 )
% Include all entries
B1 = zeros(nbb,nbb);
B2 = zeros(nbb,nbb);
B1 = -YI;
B2 = -YI;
% Delete the voltage magnitude and phase angle equations of the slack
% bus and voltage magnitude equations corresponding to PV buses
for ii = 1: nbb
if (bustype(ii) == 1)
for jj = 1: nbb
115
POWER FLOW SOLUTION METHODS

if ii == jj
B1(ii,ii) = 1;
B2(ii,ii) = 1;
else
B1(ii,jj) = 0;
B1(jj,ii) = 0;
B2(ii,jj) = 0;
B2(jj,ii) = 0;
end
end
end
if (bustype(ii) == 1) | (bustype(ii) == 2)
for jj = 1: nbb
if ii == jj
B2(ii,ii) = 1;
else
B2(ii,jj) = 0;
B2(jj,ii) = 0;
end
end
end
end
B1 = inv(B1);
B2 = inv(B2);
end
% End of FastdecoupledJacobian Function



4.3.9 A Benchmark Numerical Example

A small network (Stagg and El-Abiad, 1968) is used to illustrate the power ¬‚ow solutions
given by the Newton“Raphson and the fast decoupled methods. As shown in Figure 4.6, this
is a ¬ve-bus network containing two generators and seven transmission lines. The data are
given in function PowerFlowsData, suitable for use with either the Newton“Raphson or
the fast decoupled Matlab1 programs. The power ¬‚ow results are superimposed on the one-
line diagram of the network, and the bus voltages are given in Table 4.1.
Function PowerFlowsData, to read data for the ¬ve-bus test network, is as follows:
%The following convention is used for the four types of buses available
%in conventional power ¬‚ow studies:
%bustype = 1 is slack or swing bus
%bustype = 2 is generator PV bus
%bustype = 3 is load PQ bus
%bustype = 4 is generator PQ bus
%
%The ¬ve buses in the network shown in Figure 4.6 are numbered for the
% purpose of the power ¬‚ow solution, as follows:
116 CONVENTIONAL POWER FLOW

%North = 1
%South = 2
%Lake = 3
%Main = 4
%Elm = 5
%
%Bus data
%nbb = number of buses
%bustype = type of bus
%VM = nodal voltage magnitude
%VA = nodal voltage phase angle
nbb = 5 ;
bustype(1) = 1 ; VM(1) = 1.06 ; VA(1) =0 ;
bustype(2) = 2 ; VM(2) = 1 ; VA(2) =0 ;
bustype(3) = 3 ; VM(3) = 1 ; VA(3) =0 ;
bustype(4) = 3 ; VM(4) = 1 ; VA(4) =0 ;
bustype(5) = 3 ; VM(5) = 1 ; VA(5) =0 ;
%
%Generator data
%ngn = number of generators
%genbus = generator bus number
%PGEN = scheduled active power contributed by the generator
%QGEN = scheduled reactive power contributed by the generator
%QMAX = generator reactive power upper limit
%QMIN = generator reactive power lower limit
ngn = 2 ;
genbus(1) = 1 ; PGEN(1) = 0 ; QGEN(1) = 0 ; QMAX(1) = 5 ; QMIN(1) = -5 ;
genbus(2) = 2 ; PGEN(2) = 0.4 ; QGEN(2) = 0 ; QMAX(2) = 3 ; QMIN(2) = -3 ;
%
%Transmission line data
%ntl = number of transmission lines
%tlsend = sending end of transmission line
%tlrec = receiving end of transmission line
%tlresis = series resistance of transmission line
%tlreac = series reactance of transmission line
%tlcond = shunt conductance of transmission line
%tlsuscep = shunt susceptance of transmission line
ntl = 7 ;
tlsend(1) = 1 ; tlrec(1) = 2 ; tlresis(1) = 0.02 ; tlreac(1) = 0.06 ;
tlcond(1) = 0 ; tlsuscep(1) = 0.06 ;
tlsend(2) = 1 ; tlrec(2) = 3 ; tlresis(2) = 0.08 ; tlreac(2) = 0.24 ;
tlcond(2) = 0 ; tlsuscep(2) = 0.05 ;
tlsend(3) = 2 ; tlrec(3) = 3 ; tlresis(3) = 0.06 ; tlreac(3) = 0.18 ;
tlcond(3) = 0 ; tlsuscep(3) = 0.04 ;
tlsend(4) = 2 ; tlrec(4) = 4 ; tlresis(4) = 0.06 ; tlreac(4) = 0.18 ;
tlcond(4) = 0 ; tlsuscep(4) = 0.04 ;
117
POWER FLOW SOLUTION METHODS

tlsend(5) = 2 ; tlrec(5) = 5 ; tlresis(5) = 0.04 ; tlreac(5) = 0.12 ;
tlcond(5) = 0 ; tlsuscep(5) = 0.03 ;
tlsend(6) = 3 ; tlrec(6) = 4 ; tlresis(6) = 0.01 ; tlreac(6) = 0.03 ;
tlcond(6) = 0 ; tlsuscep(6) = 0.02 ;
tlsend(7) = 4 ; tlrec(7) = 5 ; tlresis(7) = 0.08 ; tlreac(7) = 0.24 ;
tlcond(7) = 0 ; tlsuscep(7) = 0.05 ;
%
%Shunt data
%nsh = number of shunt elements
%shbus = shunt element bus number
%shresis = resistance of shunt element
%shreac = reactance of shunt element:
%+ve for inductive reactance and “ve for capacitive reactance
nsh = 0 ;
shbus(1) = 0 ; shresis(1) = 0 ; shreac(1) = 0 ;
%
%Load data
%nld = number of load elements
%loadbus = load element bus number
%PLOAD = scheduled active power consumed at the bus
%QLOAD = scheduled reactive power consumed at the bus
nld = 4 ;
loadbus(1) = 2 ; PLOAD(1) = 0.2 ; QLOAD(1) = 0.1 ;
loadbus(2) = 3 ; PLOAD(2) = 0.45 ; QLOAD(2) = 0.15 ;
loadbus(3) = 4 ; PLOAD(3) = 0.4 ; QLOAD(3) = 0.05 ;
loadbus(4) = 5 ; PLOAD(4) = 0.6 ; QLOAD(4) = 0.1 ;

%General parameters
%itmax = maximum number of iterations permitted before the iterative
%process is terminated “ protection against in¬nite iterative loops
%tol = criterion tolerance to be met before the iterative solution is
%successfully brought to an end
itmax = 100;
tol = 1e-12;
nmax = 2*nbb;
%End of function PowerFlowsData



As expected, the Newton“Raphson and the fast decoupled methods yield practically the
same results when the power ¬‚ow problem is solved to a prescribed tolerance of 1e À 12.
The former method takes 6 iterations to converge whereas the latter takes 27 iterations.
However, it should be mentioned that one iteration of the fast decoupled method executes
much faster than one iteration of the Newton“Raphson method; no inversion (refactorisa-
tion) of the Jacobian is required at each iterative step of the fast decoupled.
It can be observed from the results presented in Table 4.1 that all nodal voltages are
within accepted voltage magnitude limits (i.e. 100 Æ 6% in the UK). The largest power ¬‚ow
118 CONVENTIONAL POWER FLOW

40 + j5
45 + j15
90.82
131.12


Main
North Lake 19.4
41.8 40.3 19.3

89.3 74.02
16.8 4.7
17.5 2.9 6.6 0.55

<<

. 5
( 17)



>>