916 PIlOCEF.DINGS OF THE lEEE, VOL. 62, NO. 7, JULV 1974

Review oF Load-Flow Calculation MethodsBRIAN STOTT, HEMBER, IEEE

In!lÜed Paper

Af)stTact-A lI1UVeyia presented OU the cummt1y availablenumerical teclmiques for power-system load-fiow calc:'alation usingthe digital computer. The review deala with methoda that have re-ceived widespread prac:tical application, recent attractive develop-menta, and other methods that have interesting or useful 'Character-istics. The anaIytica1 bases, computational requirementa, and eem-parative numerica1 performancea of the methoda are discusaed. At-tention ia given to the problema and teclmiques of adjustments inload-fiow 8Olutiona, and the suitabllities of varlous methods formodem applications nch as securlty mOlÚtoring and optimal loadfiow are examined.

I. INTRODUCTION

]CAD FLOW (or power Bow) is the solution for thestatic operating conditian af an electric-power trans-mission system, and is the most frequently performed

aí routine digital-computer power network calculations. Overthe last 20 -years an enormous amount of effort has been ex-pended in research and development on the numerical calcu-lation processoFor instance, [1] gives a large but by no meansexhaustive bibliography on the subject, comprising more than200 "respectable" .papers in the English language alone. Outof these, 134 have appeared as publications of the IEEEPower Engineering Socíety and íts predecessor.

This review wiUc1earlybe unable to cover every aspect ofthe probiem and every proposed solution algorithm, nor is itintended to compete as a catalog of references with otherrecent sources [1], [2]. The aim is to present the underlyingprincipIes and techniques of the popularly accepted ap-proaches, those more recent methods that seem to offer par-ticular promise, and a selection of other methods that containideas of practical or theoretical interest.

Perhaps the most recurrent question arising in the load-Bow field is-which is the best method to choose for a givenapplication? The answer is rarely easy. The relative proper-ties and performances of different load-flow methods can beinfluenced substantiaHyby the types and sizes of problems tobe solved,by the computing facilities available, and by theprecise details of iroplementation. Any final choice is almostinvariably a compromise between the various criteria of good-ness by which load-fiow methods are to be compared witheach other. Every such criterion is directly or indirectJy as-socÍated with financial cost, in the actual calculation itselí, inthe engineering application, or in the computer hardware andsoftware requirements.

Load-fiow calculations are performed in power-systemplanning, operational planning, and operation(control. Theyare increasingly being used to solve very large systems, tosolve multiple cases for purposessuch as outage securityassessment, and withiri more complicated calculations such asoptimization and stability. Table I gives a brief summary of

ManU1!Crlptreceived Deeember 14, 1973; revised January 21, 1974.The author i8 with the Univeraity of Waterloo, Waterloo, Ont.,

Canada, on leave from the Uníversity of Manchester Institute of Sclenceand Technology, Manchester, U. K.

TABLE I

Loxn-Ft.ow CALCULATlON5- TvpltS AND REQUIIUDlENTS

Types of Solution

Accurate ApproximateUnadjusted AdjuatedOff-Line On-LineSingIe-Case Multiple Cases

Propertíes Requíred of Load-Flow Solution Method

Hlgh speed espedally for: Iarge systemsreal-time applicationsmultiple casesinteractive applications

especially for: 1arge sy&temscomputers with small core-store avail-

abilityespedally for: ill-conditioned problema

outage studiesreal-time applications

abllity to handle conventional anel special features(adjustments, representation of power-aystem appa-ratus) ; auItability for incorporation into more eomplí-cated processes

ease (and costj of coding, maintaining, and enhancingthe algorithm anel computer prograrn based on it.

Lowstorage

Re1iability

Veraatilit y

Simpücity

some of the main types of 10acl-fIowsolutions currently inapplication, and the requirements imposed on the numericalprocesses. Each application requires a combination aí thetypes shown, e.g., some forms of security assessment call forapproximate, unadjusted, on-Iíne, rnultíple-case solutions.

A. Brief Histc1'y of Load Flow

Prior to, and for some time after, the advent of digitalcomputers, load-fíow solutions were obtained using networkanalyzers. The first really practical automatic digital solutionmethods appeared in the literature in 1956 and subsequently[10J-[12J. These Y-matrix iterative methods were weIl-suite to the early generations of computers since they re-quire minimal computer storage. Although they performsatisfactorily on many problems, they converge sIowly, andtoo often not at ali. The incentive to overcome this deficiencyled to the development of Z-matrix methods [19 J-[21],which converge more reliably but sacrifice some of the ad-vantages of Y-matrix iterative methods, notably storage andspeed when applied to large systems. Around the same time,the Newton(-Raphson) method was shown to have verypowerful convergence properties [24J, [25], but was computa-tionally uncompetitive. Major breakthroughs in power-systernnetwork computation carne in the mid-19óO's, with the de-velopment by Tinney and others of very efficient sparsity-programmed ordered elimination (31. One of its earliest sue-cesses was in dramatically improving the computing speedand storage requirements of Newton's method, which has nowcome to be widely regarded as the preeminent general-purposeload-flow approach [26J, and has been adopted by much of

r<

'·"--'7./I'" WAD-now CALCULATm. xsraons

industry. Currently, with the stimulus of increasing problemsizes, on-line applications, and system optimization, newermethods are emerging which are also expected to find wideapplication.

11. NOTATIONS

A. Comj>kx Q-uantities at Bus is, .•.ViI8,=e;+jh Nadal voltage.I, - Nadal injected current.S, -=Pi+iQ, Net nodal injected power./lS. =/lP.+j/lQ, Power mismatch,/lI, Current mismatch./lE,=/l V,IM, Correction to nodal voltage.

B. Matrices

Y=G+jBZJ

Nadal admittance matrix,Nadal impedance matrix.Jacobian matrix.

C. General

Number of buses in system.Reactance of branch between buses i

and k.

Slack bus index.

Superscripts IP,••1,and * denote specified value, calculatedvalue, and cornplex conjugare, respectively.

Subscripts G and L denote generation and load, respec-tively.

AlI unsubscripted upper-case ,quantities denote vectors arrnatrices.

kw-i denotes a bus k directly connected to bus i.kJli denotes a PQ-type bus directly connected to bus s.kEiis thesameas kwi, except that it íncludes the case k=i.

111. BASIC ANALYTICAL FORMULATION OF

LOAD-FLOW PROBLEMThe basie fo~muIation o~ the Ioad-flow problem is now

weH-known [4]-[7}, and is presented only briefly in this sec-tion.>Virtually ali load-flow ca1culation methods have beendevised initially for the solution of this basic problem. Exten-sions to 'it are usually "grafted on" afterw.ards, and Section IXdiscusses some of them.

A. Basic Equaticns and Bus. Types

A balaneed three-phase power system is assumed, and thetransmission system is rep~esented bt its positive-phase-sequence network of linear lumped serid and shunt branches,Nodal analysis is almost univ1ersallypre~erred in this applica-tion, and leads to the network nod j admittance matrixequation

1= Y·E

where matrix Y is square, sparse, and symmetrical (in theabsence of phase shifters or mutual coul?lings represented bynonbilateral network branches). I I

The static operating state lof the system is specified by theeonstraints on power andj or voltage at the network buses.Buses are categorized into three types, as follows:

A PQ bus i is one at which the total injected power isspecified:

917

S,"" = P,IP +jQ,"P

= PG,IP - PdP +j(QG,1P - QdP) = EJ,*. (2)

A PV bus i is one at which the total injected active poweris specified, and the voltage magnitude is maintained at aspecífied value by reactíve-power injection. The bus con-straints are

PtP = PGilP - PdP = Re (EJ,*) (3)

• V,SP = (e,' +/0')1/1 = V. = 'E,I. (4)

The system slack (M swing) bus is a fictitious concept,created by the load-flow analyst. It arisesbecause the systemPR.losses are not known precísely in advance of the load-flowcalculation. Therefore, the total injected active power cannotbe specified at every single bus. I tis eustomary to choose oneof the available voltage-controlled buses as slaek, and to re-gard its active power as unknown. The difference between itsexpected and solved megawatt output then represents theerrar in the prior estimate of system li R losses. li the slaekbus is a point of)arge generation, this error will be a smallproportion of the net bus generation. The slack-bus voltageangle is usually assigned as the' system complex phase refer-ence, and its complex voltage E. is therefore known.

B. Bus MismatcMs and Solution Acçuracy CriteriaAn exact load-flow solution is achieved when the bus con-

straints (2)-(4) are satisfied. Nearly alI load-fíow algorithmsare devised to satisfy (4) exactly, and in these cases, the de-gree to which the specified and ca1culated bus powers in (2)and (3) equate with each other constitutes ao accuracy test forthe solution.

For each PQ bus i the eomplex power mismatch is derivedfrom (2) and (1) as

AS, = S,"P - EJ,* = P,"P + iQ.1P - E, L Ya.*E,,*. (5).\Ei

Equatíon (5) can be separated into its real and imaginaryparts, giving active- and reactive-power mismatch expressionsrelevant to both (2) and (3). Rectangular- and polar-coordi-nate versíons are presented as follows.

For each PQ or PV bus i:

t:.Pi=P,IP- Re {Cei+H,) ~ (Gi.l:-jBa.)(e.l:-jj.l:)} (6a)

=P.IP_ Vi L (G •..•cos 0•..•+B •..•sin Oa.)V". (6b).\E,

For each PQ bus i:

/lQ,=Q,"p-Im {Ce,+jf,) ~ (Ga.-jBi.i;)(ek-jjt)} (7a)

=Q,"P-V,:E (Gi.i;sin9a.-Ba.eos9i.i;)V,I;. (7b).\E,

(1)Current-mismatch expressions are also available. From (5) wecan define

M, = AS,*/E,* = (P,"P - jQ,"P)/E,* -:E Ya.E,I;. (8).I:E'

The most common convergence eriterion used in practice ís

/lP, S Cp,

t:.Q, :::;Cq,

for ali PQ and PV buses ifor ali PQ buses i

where Cp and Cq are tolerances chosen typically in the range

r--------7 ../!/'18

I 0.01 to 10 MW/mvar. These tolerances put an approximateupper bound on theerrors in any calculated line flows, andtherefore have a direct engineering interpretation. AIternativemore complicated criteria, e.g., using the sum of the absolutevalues or squares of the mismatches, are probably superfluous.Bus voltage-change tests are often used for load-flow algo-rithms in which mismatches are not readily .available. Suchtests are sensitive to the convergence rate of the solutionprocess, and are usually used as inicial stopping criteria, afterwhich the mismatches are computed and tested.

IV. Y-MATlUX ITEUTIVE METBODS

A. Rela%aJÍtmaM lhe GafUs-Seidel MetTsodThe Y-matrix iterative methods of load-flow calculation

are based on the iterative solutíon of the linear equation (1)for the bus voltages, using the relaxation algorithm [8]. As-suming temporarily that the bus currents I are known, thecurrent mismatch (residual) at each bus is, from (1):

AI, = I, - ~ Yss:iE,

Each application of the relaxation algorithm liquidates tJ.I,byreevaluating E,. The usualcomputatíonal form of the algo-rithm is

An a1ternative form, expressed in terms of the change in E;, is

(10b)

By far the most popular way of applying algorithm (10) isto use a systematic, single-sweep, sucçessíve-dísplacementsmode of iteration, i.e., the Gauss-Seidel method, since thishas the attractions of simplicity, comparatively good per-formance, and no need to store previous values. The nameGauss-Seidel refers specifically to the above numericalmethod, but is sometimes used erroneously in the load-flowliterature to be synonymous with "successive displacements."

The matrix theory shows that the method converges if thelargest eigenvalue-modulus p of its iteration matrix (a matrixderived from Y) is less than unity [9]~A more useful sufficient,though over-stringent, condition fOf'Oconvergence is that Yshould possess strict diagonal dominance. This condition maynot be satisfied for practical power networks, but. the self-admittances of the .buses are usually large reiative to themutual admittances, and convergence is obtained. The "spec-trai radius" p is associated with the least diagonally dominantrow (or colurnn} of Y, which gives the means for identifyingnetwork characteristics that affect convergence. Junctions oívery high and low series impedances, and large capacitances,encountered with cable circuits, long EHV lines, series andshunt compensation, are detrimental to convergence becausethey weaken diagonal dominance.

B. The Slack Bus

For a Ioad-flow calculation, the bus constraints have to beincorporated into the relaxation algorithm (10). As far as theslack bus ís concerned, its complex voltage is known, andtherefore it ís simply not reevaluated. This is equivalent toreducing the order of (1) by deleting the slack bus row andcolumn from Y, and thus p is changed. The choice of slack buscan affect convergence considerably in difficult cases-by re-

PROCEEDINGS OF THE IEEE, JULY 1974

moving the least diagonally dominant row and column of Y itis possible to reduce p.

C. cue« aM Stagg Met1wdIn this, the most straightforward and standard load-Bow

method [111, the nonlinear bus .constraints are incorporatedby substituting for the unknown I, in (10) from (2). For eachPQ bus, the algorithms (10a) and (lOb), respectively, become

E, = {(Pi"P - jQ,~)/E,· - E YüE.} / Y" (lla)

and

tJ.E,= AI,/Y" "'"llS"/E,*Y"

(l1b)

(9)

The problem is now "weakly linear," with the degree of non-linearity determined by the variation during the solution ofthe denominator term E,*. The factors affectíng convergencein the purely linear problem (10) remain generally valid unlessthis variation is large.

For a PV bus, V, is specified and Q.,8P is not available. lu(l1a), Q.JJP can be replaced by Q,••.I, the calculated reactivepower, which is the same as setting.1Q, in (l1b) to zero. Afterreevaluating E" its magnitude V, is reset to V,"P while main-taining its new angle. In this strategy, the angular change t:JJ,in E, is being strongly associated with 11P, (since 11Q,=O), asseen from (l1b). A minor variant employs the correspondingQ- V association principie. Starting with the value Q,cal, Q,in (11a) is adjusted by linear extrapolation until E, takes onthe specified magnitude [13J.

D. Zero-Mismatch MethodsIn the nonlinear load-flow problem, each application of

(11) does not entirely liquidate 11f, or 11Si. Several variantsmake more effort to satisfy the bus power constraints at eachiteration. In (l1a) an inner iteration loop for E; can be createdby substituting the newly reevaluated E, frorn the left-handside back into the denominator term.E,* on the right-handside, If íterated to convergence, this inner-loop process sue-ceeds in liquidating tJ.St. An alternative way of doing the saroething is to construct and solve a paír of real quadratic equa-tions in e, and f, from the lef't-hand E, and ríght-hand E,·terms [14]. For a PV bus one of these equations isreplaced bythe magnitude constraint

(12)

For certain difficult systems, ·this gives better convergence,but otherwise the extra computation is not justified.

E. Ward and Hale Met1wdThis roethod [10] solves Iinearized versions of the previ-

ously mentioned quadratic equations at each iteration. Thesolution, for 11e, and l1f"constitutes one iteration of a localNewton solution. The method performs very similarly to theGlimn and Stagg version, from which it differs only slightly.It was derived in the original paper from the incrementalexpression, for PQ buses,

Pi'P + jQ;.P = Ei(I, + y ••11E.)*+ l1E;l,· (13)

which differs from (11) only in the last term tJ.E;l,*.

-I1srcrr: LOAD-FLOW CALCULATION METHODS

F. Polar FormulationY-matrix iterative methods are coded most effidently with

the complex voltages in rectangular formo A polar version,which can be shown to approximate to (l1b) under the as-sumptions that G"«B,, and that AO, and A Vi are small, is

AO. = - AP./B"V.2

AV. = - AQ./B"V •.

(14a)

(14b)

This approach, either in the Hubert and Hayes implementa-tion[lS], or in a more conventional successive-dísplacementsiteration mode, does not appear to offer any advantages overthe rectangular forms [1].

G. Secondary-Adjustment Method [16], [17]Whenever algorithm (11) is applied to bus i, the resulting

voltage change AEi will in part be responsible for the voltagechange at each connected bus k within the next iteratióncycle. If bus k is a PQ type, an approximate relation betweenAE, and the change AEk due to AE, is

This correction transmittal can be anticipated in advance bysubstituting for Ek on the right-hand side of (lla) the voltageE.+AEk, where AE. is given by (15).

The resulting algorithm, in the forms (l1a) or (l1b), turnsout to be unchanged, except that the normal admittanceterm Y" used in it is replaced by

Y"efl = Yi, - .L: Y"Yk'/YU.k",

Note that P V buses connected to bus i are omitted from (16).The final iteration scheme is to apply (11) but with y"efl andthen to adjust the voltage of each connected PQ bus using(15). The method gíves faster and more reliable convergencefor most systerns [1], except those containing many PVbuses.For well-conditioned systems, the reduction in the number ofiterations compared with the normal algorithm (li) is offsetby the extra computation. Its added reliability can be a greatasset, however.

H. A cceleration

Y-matrix iterative methods converge slowly, because ofthe loose mathernatical coupling between the buses-at eachiteration cycle, an improvement in each bus voltage can onlyaffect the voltage improvements of the buses directly con-nected to it. Acceleration techniques are invariably used inpractice to speed up the convergence.

The most popular acceleration method is successive over-relaxation (SOR), which can be explained using (l1b). A fixedempirically determined acceleration factor a (1<a <2) is ap-plied to each voltage-change reevaluation thus:

and can be interpreted as a smaIl linear extrapolation of Ei.With a good choice of a the convergence rate is usual!y irn-proved by a factor of two or more, compared with no accelera-tion (a = 1), and an otherwise divergent case can sometimesbe made to converge. For a given system, it is often foundthat a near-optirnal choice of a remains valid over a range ofoperating conditions. Although a cornplex value of a can beused, it is more normal to use a real value, or to accelerate the

919

real and irnaginary parts of AE, with slightly different valuesof a (11].

Several attempts have been made to devise variable ratherthan fixed srnall-extrapolation methods [1S], [18]. Introduc-íng the iteration index p, one such method [18] uses insteadof (17) the additive factor {3 rather than the multiplicativefactor a:

(18)

where

(1S)

{3,(p) = sss»:» + b2t:.S,<P-2) + b3t:.S,(p-3) + ...and b is a fixed "boost" factor usually in the range 0.75 to0.92. Hence a bus whose previous convergence is poor (Iargemismatches) wiII receive a correspondingly large extrapola-tion. Compared with SOR, this method gives faster con-vergence on some systems but not on others, and rarelyjustifies the extra computation at each iteration [1].

Large nonlinear extrapolation can be used if at some stagethe iterative process is detected to be converging sufficientlymonotonically [12], [13J. Aitken's A2 process is based on theassumption of linear or geometric convergence, i.e., that

E,ool - E;<p)

Ei&Ol - E';<p-l)

E,ooI - Ei<P-l)

E.&01 - Ei<P-2)(19)

(16)

Provided that the necessary previous values have been storedand the test for monotonic convergence of ali voltages issatisfied, (19) can be solved for the "exact" solution Eiso1• Inpractice, of course, convergence is not strictly linear, and thisis not the true solution, and the iterations must be continued.This form of acceleration.can be used concurrently with SOR,and it is usually applied automatically several times duringthe calculation. The computational overheads are high, andmost often not justified by the reduction in the number ofiterations.

(17)

I. Com putational Efficiency of Y-Matrix Iteratioe Meihods

Computationally, the salient feature of these methods isthat the nurnber of elements in the summation term in (10a)is small, averaging about three for well-developed power-system networks. The off-diagonal elements of the Y-matrix,for use in the summation, are therefore stored and addressedcompactly, often taking advantage of Y-matrix symmetry.Both the storage requirements for the network and the com-putation per iteration are then small, and roughly propor-tional to the number of buses n, Since the number of iterationsfor a large well-conditioned system is of order n, the totaliterative computing time varies approximately with n'. As thesize of problems to be solved increases, Y-matrix iterativemethods becorne Iess and less competitive with newer methodssuch as Newton's. However, their storage requirements arevery low. The Y-matrix compacting takes Iittle computation,and efficient coding of the methods is relatively very easy.

V. Z-MATRIX METHODS

A. Z-Matrix Algoritkms [19J-[22]

The major difference in principIe between the Y-matrixiterative methods and the Z-matrix methods is that in thelatter, the linear equatíon (1) is solved directly for E in termsof I using the inverse of Y:

E = y-I.[ = Z·I. (20a)

·/1./1------;7IaI

An alternative version eliminates the equation for the slackbus from (1) before inverting Y. This leads to the (n-l)thorder equation

E = Z·I+ C·E. (20b)

where rnatrix Z is different from the one in (20a), E and I donot now contain E. and 1., respectively, and C is a eonstantvector.

The basie scheme in either version incorporates the busconstraints by substituting for currents, just as in SectionIV-C. Each bus current li (i;;oés) is evaluated from

where Q, •••Q.IIP for a PQ bus, and Q.", Q.cal for a P V bus(equivalent to setting AQ. = O). The bus voltages are reevalu-ated by the iterative application of (20) and (21), taking theangular voltage change only for a PV bus. Using version(20a), it is best to update I. at every stage to reflect thechanges in the bus voltages [22]. Versions (20a) and (20b)are then numerically identical except for rounding errors.

The conventíonal Z-matrix approach uses the more power-fui successive-displacernents iteration mode, and thereíore, bydefinition, the nonsparse explicit impedance matrix Z asshown in (20). A technique of great general value is to repre-sent a portion of each bus load as a fixed shunt admittance,which in effect reduces the value of the nonlinear bus powerconstraint [22]. A number of variants on the Z-matrix algo-rithm are available [22J. As in Section IV-D, zero-rnismatchtechniques can be used, either by reiterating 'each voltage inan inner loop between (20) and (21) or by solving a pair ofreal quadratic equations [20]. This usually improves thespeed and reliability of convergence.

B. Th« Impeda.nce MatricesThe driving-point and transfer impedance matrices Z in

(20a) and (20b) are referred to ground and the slack bus, re-spectively. In the forrner case, it is usually necessary to createat least one strong artificial tie to ground to avoid numericaldifficulties when obtaining Z, because in the 'abseuce of this,Yis ill-conditioned or even singular. A large shunt admittanceinserted at the slack bus most simply achieves the desiredeffect.

The impedance matrix is conventionally computed by anordered bus-by-bus assembly process,reflecting the shunt andseries branches into the formative matrix as soon as possibleusing the inverse-matrix modification algorithm [21]. Moremodern sparsity-oriented inversion algorithms are probably

.faster [9].C. ConvergenceProperties

By virtue of the direct network solution, each bus voltageis coupled with ali bus currents. Convergence is thereforerapid and reliable, compared with Y-matrix iterative methods,and typically takes 8 to 20 iterations for medium or largesystems, but sometimes much more.

Systems with no P V buses tend to converge best. To ex-plain this phenomenon, consider in (20) the reevaluation ofany bus voltage Ei, by multiplying the relevan:t row of Z withthe vector I. This process is unable to distinguish whethereach of the currents used belongs to a P V or a PQ bus. Inother words, it fails to transmit to B, the information that themagnitude of each PV-bus voltage E., (k;;oéi) is to be held fixed,i.e., that tbe equation in (20) for Ele is in reality an equation in

PROCEEDINGS OF THE IEEE. JULY 1974

angle only. The símultaneous-displacements version of theZ-matrix method is particularly badly affected by the presenceof P V buses, although it works well without them. This isunfortunate, because the sparse triangular factors of matrix Ycan be used instead of Z in this version, with a great saving instorage and computation per iteration (see Section VIII-C).

As with all load-flow methods, a' rigorous convergenceanalysis is not easy. The incremental form of the aIgorithm is

(22)

(21)

In the símultaneous-dísplacements version, an approximateexpression for AE at iteration p is

where [J'SP-jQ] is a diagonal matrix. For convergence, AE(l')must tend to zero as p tends to infinity. Tberefore, con-vergence is best if the bus powers are small (hence the successof representing loads as shunt admittances), and if the spectralradius p of Z is small. The Z-matrix method is not usually toosensi tive to the choice of slack bus. N evertheless, a good choiceis the bus whose row surn in Z is largest, since this in effectreduces p [23].

D. Cotrlputational Efficuncy oi Z-Matrix Met1wdsThe dominating feature of Z-matrix methods is the need

to obtain, store, and iterate with the nonsparse Z matrix.Even utilizing symmetry, more than n' real words of storageare required for it, which can make entirely in-core solutionsunattractive or impossible, depending on the sizes of theproblem and the computer. Block transfers of rows of Z toand from backing store may have to be programmed.

The initial construction of Z by the assembly method takesan amount of computing of orderbetween nl and nl, depend-ing on the network configuration. The Z matrix of a gívensystem can be stored for use from study to study. Each indi-vidual network change by inverse-rnatrix modification takes2n' real operations (utilizing symmetry). More than 4n' realoperations are performed at each iteration cycle,

For problems that can be solved by both, Z-matrix meth-ods are rarely competitive with Y-matrix iterative methods,and therefore certainly not with modero methods. Computa-tion and storage becorne astronomical for véry large problems.

VI. NEWTON METHODs

A. The Newron AlgorithmThe generalized Newton(-Raphson) method is an iterative

algorithm for solving a set of simultaneous nonlinear equationsin an equal number of unknowns F(X) =O. At a given itera-tion point, each functionf.(X) is approximated by its tangenthyperplane. This linearized problem is constructed as theJacobian-rnatrix equation

F(X) = - J·AX (24)

which is then solved for the correetion AX. The squareJacobian matrix J is defined by Jü.;=df;jdx", and representsthe slopes of the tangent hyperplanes. Matrix J is highlysparse in the load-flow application, and (24) is solved directlyand rapidly by sparsity-programmed ordered triangulationand back-substitution [3].

Tbe Newton method's quadratic convergence is fasterthan that of any other load-fiow approach, and the process"homes in" very rapidly when the solution point is dose. Its

J!--7i'r STOTT: LOAD-FLOW CALCULATION METHODS

performance is sensitíve to the behaviors of the functionsF{X) and hence to their formulation. The more linear theyare, the more rapidly and reliably Newton's method con-verges. Nonsmoothness, i.e., humps, in any function j;(X) inthe regíon of interest can cause convergence delays, totalfailure, or misdirection to a nonuseful solution,

Since the chosen load-tlow functions F(X) tend not to betoo nonlinear, and reasonably good initial estimates areavailable, these difficulties are encountered infrequent1y. Infact, applied to the vast majority of practical load-flow prob-lems, Newton's method is very reliable and extremely fast inconvergence. Its sensitivity is minimal to factors causing poorconvergence with many other load-Row methods, súch aschoice of slack bus and series capacitors.

The Newton load-fIow formulations adopted to date usefor F(X) the bus power or current mismatch expressions, anddesignate the unknown bus voltages as the problem variablesX. Mathematically speaking, the complex load-flow equationsare nonanalytic, and cannot be differentiated in complexformo In order to apply Newton's method, the problem isseparated into real equations and variables. Rectangular orpolar coordinates may be used for the bus voltages.

B. The Polar Power-Mismatch Version [24]-[26]

This is the most widely used of all formulations, whoseJ acobian-rnatrix equatíon (24) can be written for convenienceof presentation in the partitioned form:

[AP] fliJIil ~§] ~[~]~~J@]

Slack-bus mismatches and voltage corrections are not includedin (25), and likewise ilQi and ilVi for each P V bus are absent.The order of the equation is therefore (npo + 2npq). The su b-matrices H, N, M, and L represent the negated partial deriva-tives of (ôb) and (7b) with respect to the relevant8's and V's,e.g., Hik= -ailPi/a8k. li buses i and k are not directly con-nected, their "mutual" terms in the J matrix are zero, and Jis thus highly sparse, with positional but not numerical syrn-metry. I n (25), each right-hand correction ilVi is usually di-vided by Vi, to gain computationally by making the terms inM and L similar to those in N and H, respectively. This doesnot alter the nurnerical perforrnance of the method.

The polar power-mismatch version converges to high ac-curacy nearly always in 2 to 5 iterations from a Elat start(V = 1 per-unit, 8 = O), independent of systern size. The ac-cepted formulation (25) can be improved by a minor modifica-tion, which very often reduces the number of iterations byone, and can avoid divergence in some extreme cases. Notingthat the performance of Newtan's method is closely associatedwith the degree af problem nanlinearity, the best left-handdefining functions are the most linear ones. If (7b) is dívidedthroughaut by Vi, only one term Qi"P/ Vi on the right-handside of this equation is not linear in V. Moreover, for practicalvalues of Qi'P and Vi, this nonlinear term is numerically rei a-tively small. I t is therefore preferable to use a problern-defining function ilQ/ V on the left-hand side of (25) in placeof ilQ. Dividing ilP by V can also be helpful, but has a rela-tively small effect, since the active-power component of theproblem is not strongly coupled with voltage magnitudes.

C. Rectangular Pouier=Mismatck Version [25]

T'his version uses the real and irnagiriary parts of the busvoltages as variables, and hence derives frorn (6a) and (7a).

921

The number of equations and variables is greater than thatfor (25), by the number of P V buses. A voltage-magnitude-squared mismatch equation from (12) is needed for each" P Vbus. With sparsity programming, this increase in order is ofIittle significance. Indeed, each iteration is rnargínally fasterthan for (25) because there are no time-consuming trigono-metrical functions (noting though that even the polar versíonavoids these as far as possible by using rectangular arithmeticin constructing (25» [26]. Frorn a series of tests carried outby the author, the rectangular version seems to be slightly lessreliable and rapid in convergence than the polar version.

By neglecting alI the "mutual" terms in the Jacobianmatrix, a .simultaneaus-displacements version oi the Wardand Hale method is obtained.

(25)

D. Current-Mismatcb VersionsRectangular or polar coordinate Newton versions can be

constructed from the real and irnaginary parts of the complexcurrent-mismatch equation (8). In both versions, it is con-venient to revert to the power-mismatch equation (6) for aP V bus, adding a voltage-squared mismatch equation if (õa)is chosen. These versions generally perform less satisfactorilythan the power-mismatch methods. The main interest lieswith the rectangular version. For PQ buses, its J-matrix"mutual" terms are simply network branch admittances, andon no-Ioad the "self" terms become likewise. The version canthus perforrn relatively reliably for difficult light-load prob-lems [27], [28]. N ate that for systerns with no P V buses andon no-load, the Jacobian-matrix equation is exactly equivalentto the Z-matrix equation (22).

E. Algorithmic Enhancements

A number of schemes are available for attempting to irn-prove the performance of Newton's method [27]. One of thesimplest oI these is to impose limits on the perrnissible sizes ofthe voltage correctíons at each iteration, thereby helping tonegotiate humps in the. defining functions. If the lirnits areperrnanently applied they should not be toa small, otherwiseconvergence will be slowed down for well-behaved problems.A better approach is to back-track as soon as divergence isseen to have started, and then apply small lirnits.

With its quadratic cqnvergence, Newton's method takesmaximal advantage of good initial voltage estimates. As withali other load-flow methods, previously stored solutions cansometimes be used as starting values. Some prograrns performone or two Gauss-Seidel íterations beIore the Newton pro-cess [26]. This is beneficial provided that the relatively weakGauss-Seidel method does not diverge when faced with adifficult problem. A most rapid and reliable Newton programcan be created by calculating good initial angular estirnatesusing the dc load flow (see Section X-A) and also good voltage-magnitude estimates by a similar technique [29].

I teration time can be saved by using the same triangulatedJacobian matrix for two or more iterations [26]. However, thelower triangular factor of J now needs to be stored, and thealgorithm loses some speed and reliability of convergence.

F. Construction and Solution of Jacobion-Mctrix EquaiionFor nonsmall systems (say, larger than 20 buses) Newton's

method relies for its computational efficiency on skilIful pro-gramming in the construction and solution of the Jacobian-rnatrix equation. The remarks on procedure in this section arebroadly applicable to alI N ewton versions, but specific detailsrefer to the standard version (25).

",-'--·'---7

f'Prior to iteration, the system buses are reordered to mini-

mize fill-in during the triangulation process later. The twomismatches and the two voltage corrections for each PQ busare ordered consecutívely to make computational use of com-mon terms in constructing (25), to ensure large rnatrix di-agonal elements, and to preserve positional symmetry of J[26J. Row (Crout) elimination for the upper triangular reduc-tioo is most efficient in sparsity programmiog. A choice of ap-proaches is then available: a) the rows (or row pairs) of (25)are computed and eliminated one at a time, or b) the completemismatch vector and compact matrix J are computed to-gether, enabling further computational use to be made ofcommon terms in the mismatch and J-matrix expressions, andenabling convergence testing to be performed on the mís-matches before deciding whether another solution of (25) isnecessary.

Approach b) requires one less solution of (25) than ap-proach a), but needs more storage if the number of elementsin the J matrix is larger than the number in its upper triangu-lar reduction. Unless the power-system network is very largeor complicated, and if a standard dynamic ordering method isused [3], this is usualIy the case.

G. Computational Effic~y of Neuuon Met1wds

Although programming technique is important in ali load-flow methods for obtaining fast execution and storage econ-omy, it is the cornerstone of methods such as Newton's thatuse ordered elimination for solving the large sparse matrixequations [3], [26]. These methods realize their full computa-tional potential only if the programming is of the very higheststandard. Packages for solving a given sparse network-matrixequation such as the Newton equation (25) are becominggenerally available, but the need to construct the equationcompactly and efficiently remains.

If these programming requirements are fully satisfied, thenthe computing time per iteration of Newton's method rises\ittle more than linear1y with the number of buses in the sys-tem, on average. This is ao empirical íunctioo of the networktopology and the success of the bus-ordering scherne. Sincethe number of iterations is size-ínvariant, the superíority ofNewton's method speedwise over previous methods increasesrapidly as the size of the system to be solved increases. Fortypical large systems, the computing time-tor one Newtoniteration is roughly equivalent to seven Gauss-Seidel itera-tions [26]. The usually satisíactory simple dynamic bus-ordering process takes less time than one Newton iteration. Iftherefore we take the example of a 5OO-bus problem that takes500 iterations using Gauss-Seidel, and four iterations (to ahigher accuracy) using Newton, the speed advantage of thelatter for the iterative calculation is about 15: 1.

Like speed, Newton's storage requirements are not exactlypredictable, because the fillup in the sparse upper triangularreduction of the J matrix is not known in advance. They alsodepend on the speed-storage tradeoff decisions made by theprogrammer. However, in the approaches most commonlyadopted by industry, a practical upper Iimit for the storagerequiredfor the J-matrix equation and its solution is about18n real words plus their addressing integers. Approach a) ofSection VI-F will usually require less. With large moderocomputers, this extra storage cornpared with, say, the Y-matrix iterative methods, does not prevent very large systernsfro m being solved in core.

PROCEEDINGS OF THE IEEE, JUL Y 1974

VII. DECOUPLED METHODS

A. The Decoupling Principie

An inherent characteristic of any practícal electric-powertransmission systern operatíng in the steady state is thestrong interdependence between active powers and bus volt-age angles, and between reactive powers and voltage magni-tudes. Correspondingly, the coupling between these "P~"and "Q- V" components of the problem is relatively weak. Ap-plied numerical methods are generally at their most efficientwhen they are able to take advantage of the physical proper-ties of the system being solved. ln the load •.flow problern,there has been a recent trend towards this objective by "de-coupling" (solving separate1y) the P~ and Q- V problems.

B. Decoupled Neunon Met1wds

ln any formal Newton method, half of the elements of theJacobian matrix represent the weak coupling referred toabove, and therefore have relatively small numerical values.These elements may be neglected. Any such approximationsto J inevitably sacrifice the true quadratic convergence prop-erty, but compensating cornputational benefits can ac-crue (33].

Many algorithrnic possibilities emerge frorn this. Forbrevity, only one recommended "good" decoupled Newtonversion is presented here. In (25) the e1ements to be neglectedare those contained in submatrices N and M. Equation (25) isthen separated into two smaller matrix equations, viz., theP~ and Q-V problems. Adopting the minor formulationmodifications of Section VI-B that can improve the polarpower-mismatch method, the decoupled Newton equationsbecome

t..P/V = A ·t..8 (26a)

and

AQ/V = C·t..V (26b)

where A and C are negated J acobian matrices. In this method,dividing t:.Q by Vis important, since it substantial1y reducesthe nonlinearity of the Q- V problem. An alternative but notsuperior version achieves a similar effect by replacíng t:.Q/Vby Irn (AI) [34].

Equations (26a) and (26b) can be constructed and solvedsirnultaneously witheach other at each iteration. A muchbetter approach is to conduct each iteration cycle by firstsolving (26a) for MJ, and use the updated 8 in constructing andthen solving (26b) for AV. Each of these construction/solu-tions can be perforrned alternately in the same storage area.The advantages of this block-suceessive iteration scheme areapparent frorn the beginning of the solution processo The firstcalculated values of e are accurate to within a few degrees,even when starting frorn e = O. The first solution of (26b) thenusually gíves remarkably good values for V, within say 0.3percent of the final solution. Subsequent iterations refine thesevalues, with the 8's taking longer to converge than the V's.

In its form (26) the decoupled method converges at leastas reliably as the formal Newton version (25) from which itwas derived. For convergence to practical accuracíes, it usu-ally takes a similar number of iterations. For very high ac-curacies ít takes more iterations, because overall quadraticconvergence is lost. However, it should be observed that theaccuracy "overkill" produced by the final iteration of a formal

;,~-#

/ srorr: LO.u>-"OW=C~T<ON ~T'OD'

f Newton solution is a bonus that has no practical value formost applications.

The decoupled Newton version saves by a factor of fouron the storage for the J matrix and its triangulation. This is ofcourse only a proportion of the total core required by a com-plete program, and the overaIl saving is more in the region of35 to 40 percent. The computation per iteration is 10 to 20percent less than for the formal Newton method.

C. Fast Decoupleâ Method

A variety of contributions have led to the evolution of thismethod [29], [33]-[37J, In its fully developed form [38], it isderived from (26) by making constant approximations tomatrices A and C. Equation (26) is then transformed, in per-unit, to

f.P/V = B'·f.(J (27a)

and

t.Q/V = B"· f.V (27b)

where

Bik' = - l/Xü (i ~ k),

B'/ = ~ l/Xii, and Bü" = - Bis.l:wi

Matrices B' and B" represent constant approximations tothe slopes of the tangent hyperplanes of the functions t.P I Vand f.Q/V, respectively, In effect, they are very elose to beingthe Jacobian matrices A and C, evaluated at system no-load.The fixed-tangent iteration method [8] embodied in (27a) and(27b) is not vulnerable to hurnps under loaded system condi-tions in the left-hand problem-defining functions.

Pursuing the decoupling principie to its logical conclusion,network elements that primarily affect the Q- V problem (e.g.,shunt susceptances and transformer off-nominal taps) shouldnot be represented in B'. Similarly, phase shifts should not berepresented in BI/. ln fact, the method has been found to workjust as well if phase shifts are left out of B' also. Consequently,both B' and B" are always symmetrical, and their constantsparse upper triangular factors are calculated and stored onceonly at the beginning of the solution. To solve (27a) and(27b), forward and backward substitutions are performedusing these factors.

The method converges very reliably, usually in 2 to 5iterations for practical accuracy on large systerns. The methodhas the decoupled property of giving a very good approximatesolution after the first one or two iterations. Provided thatthe f.PI V and t.Q/ V functions are calculated efficient1y, thespeed per iteration is roughly 5 times that of the formalNewton method, and two-thirds that of the Gauss-Seide!methods, The storage requirements for the B' and B" factorsare in between the corresponding J-matrix requirements ofthe formal and decoupled Newton methods, Using a standardtriangulation package, programming is easier than the New-ton methods, whose compact J matrices have to be cori-structed at eaeh iteration.

A further algorithmic development may sometimes finduseful application. For systerns with only one infeed point (noPV buses), rnatrices B' and B" have exactly the same struc-ture. Provided that the system does not contain significantlylarge shunt susceptances (which for best convergence should

923

be exeluded from B'), B' and B" have almost the same values.Therefore only one B matrix need be used, thus saving storageand triangulation time. Even P V buses can be eatered forusing this single-B approach [39]. For each such bus, f.Q, canbe set to zero whensolving (27b), and the calculated f. Vi canbe ignored. This is highly analogous to the treatrnent of P Vbuses in the Z-matrix approaeh, and suffers from the sameweakness. However, where storage is at a premium, and whenthe method converges, it will be attractive.

VIII. MISCELLANEOUS LOAD-FLOw METHODS

This section deals brietly with a selection of proposedmethods, not fitting into the popular categories, that haveattracted some interest within the last few years.

A. Minimization Methoâs [40], [41]In these methods, the load-fíow problem is conyerted into

the minimization of an unconstrained scalar objectíve func-tion fo, usually equal to the sum of the squares of the currentor power mismatches. The global minimum point fo = O coin-cides with the load-fíow solution, since each mismatch mustthen be zero. The nonlinear programming field provides manyminimum-seeking processes from which to choose. The mini-mization approach does not appear to be well-suited to thebasic load-fiow calculation, either from a theoretical stand-point or from practical results obtained. Most numericalmethods for solving a set of equations are at their most sue-cessfuJ when the problem formulation retains as much linear-ity as possible. In Iormulatíon j's, ali the nonlinearities in themismatch functions beco me more nonlinear by squaring thern.

B. Hybrid Newton/ Minimizaiion [42]This ingeníous idea attempts to accelerate Newton's

method. After one Newton jteration, a straight line throughthe first and second iteration points in the multidimensionalcoordinate space should pass reasonably close to the solutionpoint, if the process ís converging. This solution point is alsothe solution fo = O of the minimization forrnulation of SectionVIII-A based on the sum of the squares of the N ewton mis-match functíons. A search along the defined line is thereforemade for the minimum value of fo, and once found, this pointis accepted as the accelerated Newton iteration point. Thescheme can be applied at each Newton iteration. The gradientcomponents of fo needed during the search can be obtained

. from the evaluation of the Newton Jacobian matrix, Never-theless, the search is time-consuming, and it is probably onlyworth performing after the first iteration, if at aIl, since closerto the solution the quadratic Newton proeess homes in rap-idly. However, if the aforementioned straight line fortu-itously passes very close to the solution point, the scherne isvery successful. The method also guaran tees nondivergence,but the interpretation of a located minirnum fo~O (local orglobal?) is difficult.

C. Hybrid Newton/Z-Matrix Metkod [43]

As noted in Section V-C, the attraction of the simulta-neous-displacements version oí the Z-matrix method is thatthe sparse triangular factors of the Y matrix (utilizing sym-metry) can be used. However, the method's convergenee isweak if the system con tains P V buses. The hybrid methodthereíore applies UZ-matrix" iteration only to PQ buses, i.e,removing the rows and columns of voltage-controlled buses

frorn Y before factorizing it. The polar power-mismatchNewton version is used to solve separately for the voltageangles of P V buses. Block-successive iteration is performedbetween the Newton and "Z-matrix" iterations. If the numberof P V buses is small, the storage required for the matrix solu-tions is a Iittle more than half of that for a standard N ewtonsolution, and the computing time per iteration is about halfas much. On well-conditioned systems, the hybrid methodtakes about twice as many iterations as Newton's method.

D. Transient-Response Analog [44]

If a digital step-by-step transient-response calculation iscontinued until the system reaches its steady state, this is thesolution to the system load-flow problern. This idea has beenadapted by expressing each P V-bus active-power constraintas the differential equation

(28a)

where H, is an artificial "inertia." In one version of themethod, PQ buses are eliminated at each step by insertingshunt admittances in the network equivalent to their loads.ln another version, each PQ bus is represented by (28a) andthe reactive-power constraint

(28b)

The system is solved by step-by-step integration in theusual transient-response manner, where the initial erroneousbus-voltage estimates provide the disturbance to the system.A simple but effective discontinuous darnping device is intro-duced to encourage the dynamic system to settIe down to thesteady state rapidly.

The authors claim very high reliability of convergence.However, the method is unlikely to be competi tive with othermodern methods.

IX. AUTOllATIC CONTROLS AND AnJUSTllENTSIN LOAD-FLow SOLUTIONS

Most practical general-purpose load-flow programa con-tain adjustment facilities to simulate features of real-lifepower systerns that are not included within the basic problemformulation as in Section Ill. These features are important,but because they tend, in principie at least, to be simple anddo not contribute to the computational elegance of a proposedload-flow aIgorithm, they are frequently not discussed in theliterature. This has tended to obscure their influence on thecomparative merits (especially speeds) of the various load-flow methods: the adjustments nearly always increase thenurnber of iterations for a solution, compared with no adjust-ments.

Many of the features to be represented are automaticpower-systern devices operating under single-criterion control.The conventional approach to incorporating them in the load-flow solution process is to adjust the reIevant pararnetersand/or variables from íteration to iteration. The precise de-tails of implementing these adjustments can have a majoreffect on the overall speed of the solution. This section givessome general guides, and describes some of the features mostcommonly encountered.

A. Principies ojError-Feedback Adjusimen:We shall consider a single-criterion contrai in which a

parameter x varies in order to maintain a system quantity y

PROCEEDINGS OF TRE IEEE, JULY 1974

at a scheduled value yap, within some tolerance. The adjust-ment algorithm takes the form of a closed-loop feedbackmechanism, where x is corrected (within lirnits, if applicable)at successive load-flow iterations, in an attempt to reduce theerror t:.y = (y"P - )""'1) below the given tolerance. The adjust-ment process should not be initiated until the load-flow calcu-lation has converged sufficiently to give reliable, though notnecessarily very accurate, values for y"""I,

The sizes of the corrections t:.x must be caordinated care-íully with the convergence speed of the particular load-flowmethod used. If the corrections are toa small, overall con-vergence wilI be dictated by the slow convergence of theadjustment loop, If the corrections are toa big, they willintroduce successive perturbations into the load-flow solutionthat will slow down convergence, or more seriously, producedivergence.

In the physical apparatus, the parameter x is either a con-tinuously variable quantity, or is limited to discrete steps.If in the Jatter case x is restricted to these discrete valuesduring the adjustments, it is not always easy to avoid thepreviously mentioned convergence problems. A typical diffi-culty is hunting between the variousadjustments, A commonapproach is to ignore discreteness in the adjusted parametersuntil initial overall convergence has been achieved. Only thenare the parameters converted to their nearest discrete values,after which the load fiow is reconverged with no furtheradjustments. This can produce minor anomalies, since it maythen be found that some of the errors t:.y are slightly abovetheir tolerances.

The general "continuous" adjustment formula is simply

t:.x = a·t:.y (29)

where a is the "feedback gain," whose choice is important foreach type of control, each load-flow method, and in some caseseach system. The objective in choosing a is to minimize thetotal number of iterations while preserving reliable con-vergence.

The slowly converging load-flow methods (especially Y-matrix iterative) tend to suffer much less than the rapidlyconverging ones (Newton and decoupled) Irorn the effects ofadjustments. With the former, a smaIl empiricaI value of a isfound, which enables the adjustment algorithm to convergeat a similar slow rate to that of the load-Bow algorithm. In therapid load-fíow methods, and Newton's in particular, theiteration following a set of adjustments tends to give a fairlyaccurate load-flow solution. This means that the convergencerates of the various adjustments dictate the overall conver-gence rate much more than that of the load-flow algorithm.Consequently, the solution takes far more iterations than thecorresponding unadjusted solution. For exarnple, a typicaladjusted Newton solution will require 10 iterations.

For these rapidly converging load-flow methods, a in (29)should be approximately the sensitivity between x and y atthe operating point. For a given systern, a suitable fixedestimate of this can be calculated or found ernpirieally. Alter-natively, when this is not convenient, the estimation can bebuilt into the adjustment processo When the adjustment pro-cess is ínítiated, a trial correction t:.x(l) (not toa small) is madeon the basis of an error t:.y(O). One or more load-flow iterationsare then performed until moderate convergence accuracy isachieved, and the new error is .6.yw. An estimate of a can nowbe calculated thus:

STO'tT: LOAD-FLOWCALCUUnON llETHODS

a = ~(1) /(.6y(1) - .6y(O».

Care should be taken that the denominator in (30) is not verysmal!, otherwise a might be too inaccurate. It is useful to im-pose an upper Iimit on the change áe at any one step, in case.6y is initially very Iarge, to avoid perturbing the load-flowa1gorithm too much.

In some cases the parameter x affects the values of thenetwork admittances. Y-matrix changes can easily be accom-rnodated in Ioad-flow methods that use the admittance ele-ments directly in the network calculation at each iteration,e.g., Y·matrix iterative methods and Newton's method. Withother methods, it is more efficientto reflect the x-changes inthe nodal injections of the relevant buses.

Since most adjustment processes converge Iinearly, it is adistinct advantage to choose the tolerance on each .6y to beas large as possible, commensurate with the practical accuracyrequirements of the solution.

B. Controlled Transformer In-Phase TapsAn automatic on-Ioad tap-changing transforrner usually

controls a local bus voltage, within the sensitivity of thevoltage relay [6]. Discrete tap-by-tap adjustment ís some-times used, especially in slowly converging load-fíow methodswhere an adjustment is perforrned only every few iterations.More usually, the continuous approach is used, in which thealgorithm for an off-nominal tapping ti controlling the voltageof bus k becomes

where .6 Vi= (V,l:IP-Vi""l). For the rapidly converging load-fIow methods, a fixed ai = ± 1per-unít is frequently satisfac-tory, particularly if bus k is a transformer terminal.

Transformer tapping must give precedence to a reactive-power source, if both are designated to control the same busvoltage, i.e., the transformer only taps if the source Q limit isopera tive. The effect of line-drop compensation can easily berepresented by modifying the voltage signal Vl:cal.Instead ofcontrolling a bus voltage, the transformer may control thereactive flow through itself or one or more neighboring lines.In this case Q-flow error feedback adjustment is perforrned.

C. PV-Bus Reactioe-Pouier LimitsWhen the reactive-power source at a P V bus reaches its

upper or lower limit of output and its ability to maintain thescheduled voltage V,IP is Iost, either bf two Iimit-enforcernentmethods are available: a) the bus is converted to a PQ typewith QG,IP set to the limit, or b) an error-feedback process ofthe form (29) is entered, such that at each subsequentiteration, ViS]) is adjusted according to the present limitviolation tlQGi:

(32)

Method a) can be inserted very easily into the Y-matrixiterative and Z-matrix load-fíow methods. With Newtonmethods, bus-type conversion changes the structure of theJ-matrix equation, but with most programming schemes theextra work involved is slight, since the equation is constructedand solved at each iteration. In other load-flow methods, bus-type conversion is unattractive where it requires the refac-torization of otherwise constant sparse matrices, especialIy ifa number of Q Iimits are likely to be exceeded at different

925

(30) stages of the ca1culation. Method b) is then preferred. BothQ-limit methods should provide the opportunity for subse-quent lirnit 'backoff due to the possible effects of other adjust-ments.

D. PQ-Bus Voltage LimitsIf the voltage magnitude of a PQ bus violates specified

lirnits, reactive-power support may be switched in. Thismight represent the automatic action of physical apparatusconnected to the system. I tis also used in planning studies todetermine reactive compensation, and is sometimes used as aprogram device to prevent divergence on difficuIt systems (thearnount of reactive power is gradually reduced as the ítera-tions progress), The adjustments are the reverse of those usedfor P V-bus Q limits,

E. A utomatic P hase-Shift ControlAn automatically controlled phase-shifting transformer

regulates the fiow of active power through itself, or throughone or more neighboring lines, by varying its angle of shift <Pi.Where the load-flow method takes advantage of matrix syrn-metry, this shift is best represented by bus-injection tech-niques. The usual error-feedback adjustment approach is ap-plicable. With the rapidly converging load-flow methods, thescheme (30) has been found to perform very rapidly andreliably. Depending on the type of physical device, the ín-phase voltage-ratio variation with <P, should be represented.

(31)F. Area Interchange Control

Each area in a multiarea load-flow probIem is likely tohave a scheduled net megawatt import/export. This ínter-change constraint can be satistied for each area by sumrningits tie-line flows at each load-fíow íteration, and adjusting thescheduled generations. Transfer control is often assigned toone generator bus per area, although it can be shared betweenany number of generators in the area in some defined propor-tions. The usual error-feedback adjustment is used. For eacharea, the change in its generation .6Pa for an error .6Pr in thescheduled interchange is

tlPó = CZ' tlPz. (33)

Note that L: .6Pz over alI areas should be zero. For rapidlyconverging load-flow methods, a = 1 is suitable.

G. Functional Loads

In the basic load-flow formulation, constant PQ-bus loadsare normally assumed on the understanding that voltagevariation on the PQ bus is compensated for by tap changes atthe load side. In their absence, or if the load fiow simulates atime period in the system over which the taps do not opera te,constant-power loads do not represent the system correctly,ln fact, they can sometimes give totaIly wrong results, e.g.,that the system appears to be unstable. One of the main diffi-culties is in defining the loadjvoltage characteristics, butwhere these are available, they are often expressed in theform, for active power,

and símilarly for reactive power.If the function (34) is substituted for PdP in the load-flow

equations, it must be adjusted at each iteration as Vi changes.This causes no difficulty-rather, it tends to assist conver-

gence, because in problem terms the function is less nonlinearthan a total constant-PQ load of the same nominal size. Theterm in V' can be represented as a fixed shuntadmittance. InNewton's rnethod, an addition to the J-matrix diagonal canbe made to account for the term in V (note that this is notnecessary in {ormulations using tlP / V and tlQ/ V).

H. AuWmaJie Adjustments 11$Ne'Wt<m's Met1wdThe mathematically general methods such as Newton's

and the minimization approach are able to absorb contínuousadjustments of the type (29) into the problem formu\ation,rather than relying on the usual "parasitic" adjustment loops.The adjusted parameters beeome variables of the main load-fiow problem, and are solved for directIy by the load-llowa1gorithm. The possibilities of this approach have been ex-amined in some depth for Newton's method [27J, [30J-[32J,with the object of benefiting from Newton's powerful con-vergence while removing the obstacles to obtaining solutionsín the "normal" number of iterations. Practical success hasbeen demonstrated with transformers and phase shifters [31],and interarea transfers {30].

Details of the approach cannot be given here. However,for general applications it has two main drawbacks. The firstis that it cannot overcorne the problem oí discontinuities(limits), which inevitably require extra iterations when active.More seriously, the sparsity structure of the Jacobian matrixis usually a1tered. Controlled devices regulating quantitiesflowing through themselves or at their terminals do not causeproblems, apart from sometimes slightly weakening Newton'sconvergence properties, but other controls create zero-valueddiagonal J-matrix elements. To overcome this, the structureof the J matrix can be altered, at the expense of additionalprogramming complexity, and loss of efficiency in the triangu-lation processo

X. MULTIPLE-CASE SOLUTIONS

Many power-system calculations involve the repeated so-lution of the load-flow problem, with modifications to the net-work and/or bus constraints in each case. This section givesan introductíon to three of the more comrnon applications andtechniques, only in so far as the load-âow solutions are con-cerned.

A. Outage AssessmentsOutage assessments are made off-line in system planning

and operational planning, and on-line or semi-on-line inoperation/control (45J and interactive work [46J. From agiven base-case load-fíow solutíon, a sequence of load-flowcalculations is carried out, each simulating the loss of one ormore transmission elements, generating units or even loads.Since the purpose is to identify potential operating difficulties,high solution accuracy is frequently not required. However,computing speed has a high príority, because of the very largenumber oí cases that may have to be examined. This ís espe-cially true in real-time applications, where the reliability ofthe solutions is also very important. Startíng from the base-case solution voltages is nearly always a good approach. Theinclusion of adjustments is dependent on the time period afterthe outage which it is desired to sirnulate, and the degree ofaccuracy needed.

The simple de load flow has been ignored so far in this re-view, since it is not capable of giving more than a very ap-proximate megawatt fiow solution. However, it is found to be

PROCEEDINGS OF THE IEEE, JULY 1974

adequate in some power systems for assessing megawatt over-loads, and has the advantage of being extremely fast. Theoutage-simulation techniques used in it will also serve as auseful introduction to the similar teehniques applied to the acproblem.

The base-case dc load-Bow matríx equatíon, omitting theslack bus, is

OB = (B')-l.1'-r> = D· 1'-r> (35)

where the de admittance matrix B' is the same as that de-fined for (27a). Outages involving bus powers affect only thevector of bus powers P-P, and are easy to compute. For theoutage of a series branch of admíttance b between buses i andk, the base-case matrix B' changes. The new matrix -is ob-tained by adding to B' the highly sparse square matrixõ' M'· M, where t signifies transpose and M is a null row vectorexcept for M,==l and Mi:=-1 (i or k~s). Applying the ín-verse-matrix modification formula, the solution (JiII for thisoutage case ís

where DU: is the subtraction of columns i and k of D, and c ísthe scalar (l/b+M·DU:)-l. Given the base-case matrix D, theoutage solution (36) is performed in Iittle more than 1$arith-metie operations. The method can easily handle multipleoutages.

MatrÍx D has the usual disadvantages of nonsparsity re-ferred to in Section V-O. An alternative approach is to usethe sparse upper triangular factor of B', and use it to generateeach vector DU:as required by a repeat solution with M' as theindependent vector, While saving a great deal of storage andeomputation for D, this approach takes at least three times aslong for each outage case (programming to take the sparsity ofM' into account in the solution for Dilr). The overall speedcomparison between the twocases must be carefully evaluatedfor a particular application [9], [47]. For approximate acoutage assessments, the same teehniques can be applied to theeomplex equation (20) [48].

For greater accuracy, a true ac load-flow method must beused. The rapidly converging methods sueh as Newton's andthe decoupled versions are generally the most suitable. Sincethe accuracy requirements of most outage-assessment studiesare not too high, however, it is frequently possible to avoidreíormíng and refactorizing the network matrix or matricesfor each branch-outage case, by retaining the matrix faetors ofthe base-case solution throughout. The assumption is madethat the outage of a branchchanges only its terrninal-bus selfand mutual elements in the base-case matrix. Compensationfor these changes then has to be included in eaeh repeat solu-tion for the unknown voltages using the base-case factors. Inthe formal Newtan method, one approach [49 J is to calculateand iterate with the outaged-branch terminal-bus sensitivities,which can be obtained efficiently for ali branches at one time[50]. The same technique can be used with decoupled Newtonrnethods.

Those deeoupled methods whose network matrices aresymmetrical enable the compensation technique as in (36) tobe used [37J, [38J. For unadjusted moderate-accuracy solu-tions, sue h as in securíty monitoring, a single iteration, startedírorn the base-case voltage solution, is often sufficient, al-though in general it is desirable to provide the opportunity ofreiterating if the accuracy is inadequate (írequently the most

77srorr : r.oxo-rt.ow '''''ULATION ."ao"important outage cases). Using the method of (27) for line-overload monitoring, for instance, the following scheme maybe used. One outage iteration of (27a) is performed, and onthe basis of the new f)'s, the branch megawatt flows are ob-tained. If no lines are approaching overload, this outage caseis discontinued; otherwise (27b) is solved, and the solution isiterated to within acceptable mismatches [38].

On the bases of speed, storage, and accuracy in these ap-plications, there seems to be no point in using the formal ordecoupled Newton methods in preference to the symmetricdecoupled versions. Note that the symrnetric decoupled for-mulations of [35] and [37] give direct soIutions for 8 and Vrather than /j,f) and /j, V. These are simply alternative arrange-ments of the equivalent "Newton-type" equation forms suchas that of (27).

B. Optimal Load Flow [51]This section discusses optimalload flow only in the context

of methods that use conventional load-flow solutions withinthe processoAn optimal load-fíow calculation optimizes theactive- and reactive-power dispatch of a system, incloding ascontrol variables those single-criterion-control parametersthat are adjusted during an ordinary load-fiow solution. Allthe relevant static limit constraints on the systern operationare enforced.

Some practical approaches to this problem, including thenow classical Dommel and Tinney method [52], adjust thecontrol variables in some optimurn-seeking manner in betweenconventional load-fíow solutions. Polar-coordinate Ioad-flowformuJations are the most natural choice for the optimizationapplication, since voltage magnitudes in particular are thenexplicitly available as variables for control and lirnít-enforce-ment purposes.

Apart from the effects of transrnission-apparatus param-eter adjustments, the network remains fixed in admittancevalues and configuration (although outage constraints canalso be incorporated [53]). Each resolution of the load flow tocalculate the efíects of the control-variable adjustments usesthe previous point to start it. If the adjustments are not large,a rapidly converging method such as Newton's wiII give highaccuracy in one or two iterations. Newton's method is alsovaluable because its Jacobian matrix contains the informationneeded to determine the next incremental adjustments of thecontrol variables [52]. If the latest J matrix frorn the load- ,flow solution is to be used, its lower triangular factor must bestored. The fast decoupled method of Section VII-C has alsobeen used for the load-fíow solutions and in calculating theincremental costs (Lagrange multípliers) [54]. It requires lessstorage, and, being much faster than Newton's method forpractical accuracies, it gives an overall speed advantage inthe optimal load flow of more than 2: 1.

C. Multimachine Dynamic Response CalculationsThe conventional approach to the calculation of power-

system stability uses a step-by-step process in which themachine differential equations are integrated alternately witha load-fiow solution of the transmission system. The solutionmethods and degrees of system representation vary widely,and it is only possible here to make some general commentson the Ioad-flow portion of the overalI calculation [55].

The load-flow problem to be solved after or during eachintegratíon step ís characterízed by multiple slack buses. Eachof these represents the temporarily fixed complex voltage be-

927

hind the relevant machine impedance(s), transformed into thenetwork reíerence frame. The convergence powers of ali load-flow methods are enhanced by having many slack buses,which are easy to handle in the respective Ioad-flow programs.The load-flow solutions are also assísted by representing loadsor parts of loads by fixed shunt admittances, by netwark re-duction (where appropriate), and by the good starting volt-ages available from the last solution (except after switchingoperations) .

Saliency in the machine impedances introduces time-varying nonbilateral network branches, usually Nortonshunts [56]. The variations of these impedances can be repre-sented by bus-injection techniques if the load-flow methodrelies on a constant and/or symmetrical network matrix.

Often, some or ali bus power constraints are representedas fixed shunt admittances. Network reduction can then beundertaken, but the degradation of possible network-sparsityexploitation must be weighed carefully against the expectedimprovement in load-flow convergence. Nonadrnittance loadsmust be corrected at each iteration (see Section IX-G), andcan gíve convergence problerns, particularly constant-powerloads, if their bus voltages go very low or high.

Final\y, the integration step sizes used can affect thechoice of load-flow method. Small step sizes, where the varia-tion of system conditíons over the step is small, do not imposesuch stringent requirements on the convergence power of theload-fíow method.

The Y·matrix iterative methods can very convenientlyhandle the whole range of system representational require-ments, and need minimum storage, depending on whether anyform of network reduction is performed. They are best suitedto small integration step lengths. The method of Section IV-Gseems to be advantageous [or its reliability, since there are noPV buses. Z-matrix methods have also been widely used, andprofit from the same facto The sparse-factor form will usuallybe preferable to the explicit inverse, although the latter'sability to iterate with successive displacements gives it anadvantage for difficult cases.

Newton's method is suitable for larger integration steplengths, and can handle saliency (the J matrix is nonsy rn-metrical and is constructed repeatedly) and nonadmittanceloads well. The rectangular current-rnismatch version is proboably best [S6J, since there are no PV buses, and loads repre-sented as shunt admittances present it with a no-load situa-tion (see Section VI-D). Under these conditions, the methodis the same as (22), apart from the tirne-varying rnachine-impedance terms in it. If the J matrix is formed and triangu-lated at each iteration, one iteration per solution is usualIysufficient, depending on the step size and integration ap-proach. However, this is still computationally expensive, andconstant J-matrix upper and lower triangular factors, con-taining mean values of the machine impedances, can be usedover many integration steps without updating.

It should be noted that for large-disturbance studies, thedecoupled methods are not appropriate, since the decouplingprincipie will be violated.

XI. CONCLUDING REMARKS

This review has tried to make the necessary compromisebetween presenting a very large number of the available load-flow techniques, and giving due attention to the relevantprinciples and computational implications. For the most part,therefore, it has been possible to consider only the mainstream

1ideas and algorithms. Scores of other methods and variants,some of them highly original, have had to be omitted. Inaddition, a number of other specific aspects of the load-fiowproblem have not been touched on.

One of the topics not covered is that of system tearing.This subject lost a lot of impetus in the Ioad-fíow applicationafter the introduction of sparse rnatrix methods (with whichit is not incompatible, however). The sizes of problems to besolved have increased at a faster rate than the developrnent ofcheaper computer core memory elements. With this, and theprospects of multicomputer applications in power systems,diakoptical [57] and block-iteration [58] techniques for loadflow have been receiving fresh attention. Tearing is also oneway of trying toavoid excessive fillup in the sparse-rnatrixtriangulation processoThe subject of load-flow equivalents iscurrentIy of great importance, especially for real-time controland monitoring. Another topic not mentioned in the review isthe efficient simulation of controlled HV dc Iinks in the load-flow calculation. The effect of HV dc links on the ac-systerncalculatíon is similar to the adjustment problern described inSection IX. A more fundamental issue is the nonuniqueness ofload-flow solutions [59]. This is when two or more differentsolutions, each equally feasible from an engineering view-point, can be obtained. Although, fortunately, this phenorne-non arises inírequentIy in practice, there seems to be no re-liable and inexpensive way of identifying cases where severalsuch solutions exist, nor is there a criterion for deterrniningwhich solution most properly represents the Iikely operationof the system. Finally, the question of data uncertainty andits effects on the credibility of the calculated load flows isimportant, and is currentIy being examined [60]-[62].

Algorithmic improvements in load-flow calculation overthe years have tended to be measured by speed, rather morethan by the other criteria of goodness listed in Table I. Forthe traditional off-line single-case solution, algorithmic speedhas ceased to be relevant. Using well-programmed versions ofthe modern methods with sparse-rnatrix triangulation, inputjoutput is much more the limiting time facto r than theiterative solution process, and should therefore be given, moreattention in the design and coding of programs (NB, thefastest algorithms and computers can now perforrn accurateiterative solutions for 1000-bus load fíows in under one CPU-second). The need for speed remains, however, for repeatedsolutioos such as in outage-security assessment, especialIy inreal-time applications, and where tbe dc load flow does notgive an acceptable minimum accuracy.

Storage is a perennial problem for many users, for twomain reasons. Small computers are now fast enough to solvemedium-large systems. Those organizations with large com-puters tend to want to solve very large systems, up to severalthousand buses. Moreover, computer managers are usual1yreluctant to allocate most or all of the core area to a single(and noncommercial) job. This problem might be combattedby developing even more effective ordering methods insparsity exploitation (for very large systerns), andjor systemtearing applications. Unless speed is of the essence, a Iimitednumber of block data transfers to and from fast backing storemay be tolerable. Of the modern algorithms, the formalNewton method is the most storage-consumíng, especiallywhen the lower triangular factor of the Jacobian matrix isneeded. The efficient row-elimination approach to triangula-tion is not well-suited to block transferring. The fast de-coupled method offers some scope, since with about 25 trans-

PROCEEDINGS OF THE IEEE, JULY 1974

fers per iterative solution, its storage requirements are similarto those of an in-core Y-matrix iterative solution.

Some of the modern load-flow algorithms converge com-paratively very reliably, but no known method that is fastenough for routine application can guarantee to solve anyfeasible problem. One of the main values of absolute reliabilityis that, in the event of failure to converge, the engineer canthen assume with confidence that the data or the systernmodeled are inadequate. For operational purposes, it is prob-ably enough that the method should not break down for anyphysically operable systern condition. Newton's method, withenhancernents, and the fast decoupled method seem to beapproaching this requirement. For planning, even greater re-liability properties are needed, since the engineer will some-times want to solve cases that do not represent practicaloperation, e.g., to investigate stability margins,

The Newton methods are certainly the most versatilernethods available. They have been used in a wide variety ofsystem optirnization calculations, they provide sensitivityanalyses, and they can be used in modero dynamic-responseand outage-assessment calculations. As far as adjustments areccncerned, the case for pursuing the formal incorporation ofsingle-criterion controls into Newton's method does not ap-pear to be strong, with the present state of art in sparsematrix techniques. Only a restricted set of controls can bethus incIuded efficient1y,and without considerable extra com-plications.As soon as one conventional adjustment has to beintroduced into a program, the benefits of the formal1yincorporated controls are largely lost , A better approach mightbe to retain the simple conventional adjustment philosophythroughout, with careful implementation.Where speed is stillimportant, the fast decoupled load-fíow method can be used.An adjusted solution that takes say 12 iterations with thismethod is equivalent in time to 3 iterations of the standardNewton method.

This review is concluded by emphasizing the importanceof sparse-matrix techniques in the exploitation and furtherdevelopment of load-fiow algorithms. Modero methods usingsparse-matrix factorization need great attention toprogram-ming efficiency if tbey are to ful1y realize their considerableadvantages over previous methods. The onus is thereforeplaced on the load-flow analyst to gain sufficient understand-ing of the details and implications of these techniques.

REFERENCES.Ge1ural

[1] M. R. G. Al-Shakarchl, "Nodal iteratíve load flow," M.Sc. thesis,Uníversity of Manche8ter Institute of Science and Tedlnology,Maochester, U. K., 1973. A copy of the bibliography 00 load-fiowcalculation contained in this disaertation ia available frorn B. Stott.

[2] D. A. Conner, "Repreeentatíve bibliography ou load-fiow analysiaand related topícs," preeeoted at the IEEE PES Winter Meet.,NewYork, Ian. 28--Feb. 2,1973. Paper C73-104-7.

[3] W. F. Tinneyand I.W. Walker, "Direct solutions of sparee networkequatíons by optimally ordered triangular iactorizatioo, n Proc,IEEE, vol. 55, pp. 1801-1809, Nov. 1967.

[41 M. A: Laughtoo and M. W. Humphrey Davies, "Numerical teca-níquel! in the solutíon of power system load fiow problema, n Proc,Inst. Elec. E1I1. voL, 111, pp..1575--1588,Sept. 1964.

[51 A. M. Sassoo and F. J. Jaimes, "Digital methods applied to powertlaw atudíes, " IEEE Trens. Pouur App. Svst., vol. PA&86, pp. 860-867, Iu1y 1967.

[6] G. W. Stagg and A. H. EI-Abiad, Comp..úr MeIJwds ;" POfIler Syslm<A1JOlysis. NewYock: McGraw-Hill, 1968.

[7] O. EIgerd, E!edric EurlY Systems Tlw>ry. New York: McGraw-HUl. 1971.

[8] A. Ra1aton, A First Course ;11 Numeric<ú A1JOlysis. New Yock:McGraw-Hill,1965.

[9] W. F. Tinney and W. L. Powell, "Comparisoo af matrix iover&Íon

------fi1/

/lj!I STOTT: LOAD-FLOW CALCULATION METHODSf

and sparse triangular factorization for solution of power networkproblema. " presented at Romania/U.S.A. Research Serninar, Bucha-rest, Romania, June J-6, 1974.

Y -M atri» I terative Load-Floui M ethods[IOJ J. B. Ward and H. W. Hale, "Digital cornputer solution of power-

flow problema," AIEE Trens. (Poioer App. Syst.), vol. 75, pp. 398-404, June 1956. .

{lI] A. F. Glimn and G. W. Stagg, «Automatic calculation of load fiows,"AIEE Trons. (PO'UJt:l' App. Syst.), vol. 76, pp. 817-828, Oct. 1957.

[12J R. J. Brown and W. F. Tinney, "Digital solutions for large powernetworks," AIEE Trens. (Pawt:l' App. Syst.), voI. 76, pp. 347~355,June 1957.

(13] D. G. Taylor and J. A. Treece, "Load fiow analysis by the Gaues-Seidel method," presented at the Symp, on Power Systems LoadFlow Analysís, University of Manchester Institute of Scíence andTechoology, Maochester, U. K., 1967.

[14] C. Trevino, "Cases of diflicult convergence ln load flow problema,"presented at IEEE Winter Power Meet., New York, Jao. 31-Feb. 5, 1971. Paper 71 CP 62-PWR.

(15J F. J. Hubert and D. R. Hayes, "A rapíd digital cornputer solution forpower system network load-flow," IEEE Trens, PO'IPer App. Syst.,voI. PAS-90, pp. 934-940, May/June 1971.

(16J W. W. Maslin, S. T. Matraszelc, C. H. Rush, and J. G. Irwin, "Apower system planning cornputer program package emphasiziogfiexibility and compatíbílity," presented at lhe IEEE SummerPower Meet., Los Aageles, July 1970. Paper 70 CP 684-PWR.

(17] R. Podmore and J. M. Undrill, "Modified nodal iterative load flowalgorithm to handle eeries capacitive branches," IEEE Trens. PfYf1JerApp, Syst., voI. PAS-92, pp. 1379-1387, July/ Aug. 1973.

[18) J. A. Treece, "Bootstrap Gauss-Seidel load fíow," Proc. Inst, Elec,Eal., voI. 116, pp. 866-870, May 1969.

Z-Matrix Load-FICIIIJ Met/wd.s[19J P. P. Gupta and M. W. Humphrey Davies, "Digital computers ín

power eystem analysís," Proc, Inst, Elec, E1I,., voI. 108A, pn, 383-404, Jao. 1961.

[20] A. Brameller and J. K. Deomead, "Some improved methods ofdigital network analvsía," Proc, I ••st. Elec, E1II., voI. 100A, pp. 109-116, Feb. 1962. . .

[21] H. E. Brown, G. K. Carter, H. H. Happ, and C. E. Person, "Powerflow solution by impedance rnatrix iterative method," IEEE TrlJ1IS.PO'UJt:l'App. Syst., vol. PAS-82, pp. 1-10, Apr. 1963.

[22J --, "Z-rnatrlx a1gorithms in load-flow programe," IEEE Trens,Power App. Syst., voI. PAS-87, pp. 807-814, Mar. 1968.

(23] L. L. Frerls and A. M. Sasson, "Investigationa on the load-flowproblern," Proc, Inst, Elec, E1Ig., vol. 114, p. 1960, Dec, 1967.

N etDto1I Load-PlcwJ M et/wd.s[24] J. E. Van Ness, "Iteration methods for digital load flow studies,"

AIEE Trens. (PC1fDerApp. Syst.), voI. 78, pp. 583-588; Aug. 1959.[25] J. E. Van Ness and J. H. Griffin, "Elirnination methods for load-

Bdw studies," AIEE Trens, (PlYrDer App. Syst.), vol. 50, pp, 299-304, June 1961.

[26J W. F. Tinney and C. E. Hart, ·Power fiow 8OIutlon by Newton'smethod," IEEE TrlJ1IS. Poeer App. SysI., vol. PAS-86, pp. 1«9-1456, Nov. 1967.

[27] H. W. Domme1, W. F. Tlnney, and W. L. Powell, "Further develop-ments in Newton's method for power system applicatlons," pre-eented at IEEE Winter Power Meet., New York, Jan, 25-30, 1970.Paper 70 CP 161-PWR.

[28J H. W. Dommel, disc:ussion of [26], IEEE Trcfts. PotDU App. Syst.,vol. PAS-86, pp. 1456-1458, Nov. 1967.

[29J B. Stott, "EfIective starting procees for Newton-Raphson loadfiows," Proc. bul. EUc. EfJ,., vaI. 118, pp. 983-987, Aug. 1971.

[30] J. P. Brittoo, "Improved area interchange control foe Newton'smethod load flows," IEEE Tra1lS. P07ID' App. Syst., vol. PAS-88,pp. 1577-1581, Oct. 1969.

[31] N. M. Petersoo and W. S. Meyer, "Automatic adjustment of traoa-former and phase-shifter tape ln the Newton power fiow," IEEETra1lS. PO'IIJer App. Sysl., vol. PAS-90, pp. 103-108, Jan./Feb. 1971.

[32] J. P. Britton, "Improved load flow performance through a mOl'egeneral equatioo form," IEEE TrlJff,S. PClllJer App. Syst., vol. PAS-90, pp. 109-116, Jao./Feb. 1971.

Deccvpled Load-FIO'Ill M el/wd.s[33J G. W. Stagg and A. G. Phadlee, "Real-time evaluation Df power-

system contingeodeB-Üetection of steady l!tate overloads," pre-sented at the IEEE Summer Power Meet., Loe Angelee, Ju1y 1970.Paper 70 CP 692-PWR.

[34J B. Stott, "Decoupled Newton load flow," IEEE TrlJ1IS. PO'It1t:l'App.Syst., vol. PAS-91, pp. 1955-1959, Sept./Oct. 1972.

[35J S. T. Despotovic, B. S. Babic, and V. P. Maatilovic, "A rapld andrellable method for 80Iving load flow problemB,» IEEE Tra1lS.P07Ier App. Sysl., vol. PAS-90, pp. 123-130, Jan./Feb. 1971.

929

{36] K. Uemura, "Power flow solution by a Z-matrix type method and itsapplication to contingency evaluation,' presented at IEE1i: PICAConf., Boston, 1971.

[37J N. M. Peterson, W. F. Tinney, and D. W. Bree, Jr., "Iterative linearac power fíow solution for fast approximate outage studies." IEEETra?ts. Power App. Syst., vol. PAS-91 , pp. 2048-2056, Sept./Oct.1972.

(38] B. Stott and O. Alsac, "Fast decoupled load flow,' presented atIEEE PES Summer Meet., Vancouver, July 15-20, 1973. PaperT 73 463-7. Also, IEEE Trens. Pouier App. Syst., vol. PAS-93, pp.859-869, May/June 1974.

[39J E. Hobson, discussion of {38J, IEEE Trans, PO'IPer App. Syst ; vol.PAS-93, p. 868, May/June 1974..

MisceJlatuOUS Load-Floui MetMds[40] Y. Wallach, "Gradient methods for load fíow-problerns," IEEE

Trens. Pouier APp. Syst., voI. PAS-87, pp. 1314-1318, May 1968.[41] A. M. Sasson, "Nonlinear programmiog solutiona for the load-flow,

mlnímum-lose, and economic diapatcb.iog problema," IEEE TrlJ1IS.PO'UJerApp. Svst., vol. PAS-88, pp. 399-409, Apr. 1969.

[42] A. M. Saason, C. Trevíno, and F. Aboytes, "Irnproved Newton's .load flow through a minlrnizatioo technique," presented at IEEEPICA Conf., Boeton, 1971.

[43] Y. P. Dusonchet, S. N. Taluledar, H. E. Sinnott, and A. H. EI-Abiad,"Load f1oW8usiog a combination ai point Jacobi and Newton'smethcds," IEEE Trans, PO'It1er App. Syst., vol. PAS-90, pp. 941-949,May/June 1971.

(44] R. H. Galloway, ]. Taylor, W. D. Hogg, and M. Scott, "New ap-proach to power-system load-flow analysis in a digital computer,"Proc, Inst, Elec. Eng., voI. 117, pp. 165-169, Jan. 1970.

M rdtiple-C ase Load-Floui Soluticms[45] T. E. Dy Liacco, "Real-time cornputer control Df power systems, "

thís Issue, oo. 884-891.[46] J. M. Undrill, F. P. deMello, T. E. Koetyniak, and R. J. Mills,

"Interactive computation io power system analysis,' this issue, PIlo1009-1018.

[47] B. Fox aod A. M. Revíagton, "Network calculatíons for on-Iínecontrol of a power system," in Proc, 1EE C01If. 011 Comlntlers ;11Pouier Syslem OperlÚio1I and C01itrol'(Bournemouth, U. K.), pp. 261-275, 1972.

[48] H. E. Brown, "Contiogencies evaluated by a Z-rnatrix rnethod,"IEEE Trens, PO'IIJerApp. Syst., vol, PAS-88, DP. 409--U2, Apr, 1969.

[49] M. S. Sachdev and S. A. Ibrahím, "A fast approximate technique foroutage studies in power system planníng and operatíon," presentedat IEEE PES Summer Meet., Vancouver, July 15-20, 1973. PaperT 73469-4.

[56) K. Takahashi, J. Fagan, and M. Chen, "Formatíon of a sparse busimpedance matrix and its application to short circuit atudy," pre-sented at IEEE PICA Conf., Mioneapolis, 1973.

[51J A. M. Saseon and H. M. Merrill, "Power system optitriization tech-niques, " th.is íssue, pp. 959-972.

[52] H. W. Dornrnel and W..F. Tinney, "Optimal power flow solutíons,"IEEE Trens, POUJt:l' App. Syst., voI. PAS-87, pp. 1866-1876, Oct.1968.

[53] O. Alsac aod '8. Stott, "Optimal load fiow with steady-state se-curity," presented at IEEE PES Summer Meet., Vancouver, July15-20, 1973. Paper T 73 484-3.

(54J B. Stott and O. Alsac, cJoeure of [38], IEEE Trans, Poeer App.. Syst., vol. PAS-93, p. 869, May/June 1974.

[55] B. Stott, "Fonnulation and solution of matrix equationa 10 mult!-machine stability calcuíatíons," preeented at Symp. on Power Sys-tem Dynarníca, University of Manchester Instítute of Science andTechoology, Manchester, U. K., Sept, 1973.

[56] H. W. Dommel and N. Sato, "Fast transieot stability soIUtlOll8,"IEEE TraftS. PO'IlJt:l' App. SysI., vol. PAS-91, pp. 1643-1650, July/Aug.1972.

Ftutller Aspects af Loa.d FlcwJ

[57] H. H. Happ, "DialcoptiClt-The solution of system problema bytearing, " this ÍS8ue,pp. 930-940.

[58J W. E. Boearge, J. A. Jordan, and W. A. Murray, "A non-linear block:SOR-Newú>n \oad f/ow algorithm,' presented at IEEE PES Summer·Meet., Vancouver, JUly 15-20, 1973. Paper C 73 644-5.

[59J A. J. Kocsak, "On the queetion of uniqueness Df stable load-flowsoIutions,' IEEE TraftS. PlYrDer App. Syst., vol. PAS-91, pp. 1093-1100, May/June 1972.

[6OJL. S. Van SIyck. and J. F. Dopazo, "Conventional load flow notauited for real-time power system monitoriog," preaented at IEEEPICA Conf., Minneapolis, Minn., 1973.

[61] B. Borleowska, "Probabilistic load fiow,' preeeoted at IEEE PESSummer Meet., Vancouver, July 15-20, 1973. Paper T 73485-0.

[62J J. F. Dopazo, O. A. Klitio and A. M. Sasson, "Stochastic loadflOW8,».presented at IEEE PES Summer Meet., Anabeim, Calif.,July 14-19,1974. Paper T 74308-3.