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