3648 IEEE Transactions on Power Apparatus and Systems, Vol. PAS-101, No. 10 October 1982 LOAD-FLOW SOLUTIONS FOR ILL-CONDITIONED POWER SYSTEMS BY A NEWTON-LIKE METHOD S.C. Tripathy G. Durga Prasad O.P. Malik G.S. Hope Dept. of Electrical Engineering Dept. of Electrical Engineering Indian Institute of Technology University of Calgary New Delhi -110016 India Calgary, Alta. T2N lN4, Canada ABSTRACr (v) choice of acceleration factors. In this paper mathematician K.M. Brown's method is Networks which have the above features are desusse to solve load-flow problems. The method is par-cribed by a system of ill-conditioned nonlinear algeticuularl effective for solving of ill-conditioned non-braic equations. In other words, a small change in linear algebraic equations. It is a variation of New-parameter produces a large'change in the solution. ton's method incorporating Gaussian elimination in such a way that the most recent infonnation is always used In this paper a new algorithm using K.M. Brown's at each step of the algorithm; similar to what is done method [7,8] is presented. It is useful when formlaii the Gauss-Seidel process. The iteration converges tion of the load-flow problem results in ill-conditionloccall and the convergence is quadratic in nature. A ed equation and the nodal admttance matrix. This is general discussion of ill-conditioning of a system of a quadratically convergent Newton-like method based upalgeebrai equations is given, and 'it is also show by on Galussian elimination. In Brown's method, each equatth fixed-point formulation that the.-proposed method tion is expanded in the approximate Taylor's series; falls in the general category of sucessive approxima-however, the most recent information available is imtiio methods. Digital computer solutions by the pro-mediately used in the construction of the next fuction posed method are given for cases for which the standard argment This procedure is similar to the procedure load-flow methods failed to converge, namely 11-, 13-in the Gauss-Seidel process for the solution of nonlinaan 43-bus ill-conditioned test systems. A comparison ear sets of equations. This contrasts sharply with of this metlhod with the stanidard load-flow methods is Newton's method in which all equations are treated simalls presented for the well-conditioned AEP 30-' and ultaneously It has been proven that Brown's method 57-bus systemns. [7] is not equivalent to Newton's method. An optional feature of Brown's method is that the partial derivatiive of the equations are replaced by their first differrenc quotient approximations similar to the discrete INTRODUJCIION Newton's method. It is shown by the fixed-point formulattio [13] that the proposed load-flow method falls Load-flow calculations are performed in system into the general category of successive approximation methods such as Gauss-Seidel and Newton's method. Digi-planning. An early approach was the Gauss-Seidel iter-tal comuter results of the proposed Brown's load-fl ative method [1] using the nodal-admittance-method and algorithm are given for 11-, 13-and 43-bus ill-condi-this was further improved by using the nodal-impedance tioned power systems [6,12], forwhichstandardload-matrix method [2,9].-Newton's method [3] using a nodal flw mtos 1,4] naey Gus-id,Nwoadmitttanc matrix has gained-widespread popularity becaaus of its quadratic convergence characteristics. Raphson and fast decoupled algorithms failed to conii has limitations in small-core-computer ap-verge. Although the Brameller and Denmead method [2] However, ~~~~~~~~~~~and Fletcher-Powell method [5] work for ill-conditione'd plications where the weakly convergent Gauss-Seidel proble they shw poor cn rg A comaisoneo methd isgenrall resuitble Extnsie meory problems, they show poor convergence. A,co, arlson of method is generally more suitable. Extensive meiry the proposed Brown's load-flow algorithm with the raequiremt in larepow system calculations mor ti-standard load-flow methods [1,3,4,5] is presented for vatedthe xplotatio of parsty wih orered 11-, 13-and 43-bus ill-conditioned test systems, -as elimination and skilful progranuning in the Newt,on n~th well as for1 AEP 14-, 30-and 57-bus well-conditioned od f3]. Recently, advantage has been taken of loose test physical interaction between MW and MVAr flows by math-systems [9]. ematically decoupling the MW-e and MVAr-V-calculations LIST OF SYMBOLS [4]. Despite the substantial progress there are still k = condition number some difficulties with som'eof the above methods [5,6]. [J] = Jacobian matrix Features which cause instability and divergence in load-flow solution methods are: ' . .: ~~~~~~~~~N =total number of ukown variables (i) position of the reference slack bus ....-.x=vectorofukownvariables (ii) existence of negative line reactance (iii) certain types of radial systems x° = initial assumption to the vector of knowns (iv) high ratio of long-to-short line reactance for x* = solution vector lines terminating on the same bus f= ith fuction of system of equations f gj = ith function of system of equations , F = iteration function of Brown's method fi = af /ax. partial derivative of function fi w.r.t. x. 82 WM 021-4 A paper recomended and approved by the3 IEEE Power System Engineering Committee of the IEEE gi = ag./ax. partial deriv5ative of fuction gi Power Engineering Society for presentation at the IEEE 3 w.r.t. x. PES 1982 Winter Meeting, New York, Nlew York, January 31-D February 5, 1982. Man^uscript s-ubmitted January 21, 1981; T = transpose of matrix (superscript) made available for printing Ootober 23, 1981. n = iteration cout (superscript) -~~~~~0018-9510/82/1000-3648$00.75 ©B 1982 IEEE3649 P. = ith bus real power mismatch stable. K.M. Brown's method is described in detail in 1. = ith bus reactive power mismatch Appendix 1. The application of this method to a system Q. = ith bus reactive power mismatch of nonlinear equations is also demonstrated in Appendix Pi(cal) = ith bus calculated real power 1. Brown's method can be formalized by writing terms Q = ith bus calculated reactive power of the iteration function F = (fl, f**.N) used, be-Q(cal) = ginning with a starting guess x°, to form successively: P = ith bus scheduled real power i(sch) xn+1 =F( ) n 0 p2 3 Q = ith bus scheduled reactive power x = P(xn) n = 0, 1, 2, .( i(sch) P. = ith bus net real power The iteration fuction, F, in this method is given by 1 Qi = ith bus net reactive power i-l V .,V = voltage magnitude at bus i and j F.(x1, * x1) = Xi -, /g x 1 i i XN) Xi (gN-i+l~~~~~~~~~~~~~~~,x. N-i+l,x. e., O.= voltage phase angles at bus i and j J ' 1 Y..= G. +jB. = (i,j)th element of bus admittance (F x i 1, 2, .,N (4) matrix n [H] = transverse upper triangular matrix formulated where I = 0 whenever m > n (5) in forward part of Brown's method j=m E = specified tolerance and gi is defined as follows: = machine tolerance A= eigenvalue. g= f1(x1, x2, ..., xN) ILL-CONDITIONED SYSTEMS AND EIGENVAL[E ANALYSIS gi =fi(Xl, X2 .. XN-i+l LN-i+2 * LN) (6) A computer formulation of a problem-is defined to be ill-conditioned if computed values are very -sensitiiv to small changes in input value. A matrix may N fN(Xl1 L2 L3 N**, LN) have some eigenvalues which are very sensitive to small changes in its elements while others are comparatively The linear terms L. *L** LN, i 1, * N-l1 insensitive. It is convenient to have some number are themselves functlons oi x. and are obtained recurcallle a 'Condition Number' which,defines the condition sively by successive substitution in the system: of a matrix with respect to the computing problem. Ideally, it should give. some overall assessment of the rate of change of the solution with respect to changes L = x. -y n ) (L -Xi) in coefficients. It may be recalled that the size of 1 1 J ) the condition number k([J]) of a Jacobian matrix [J] is defined to be IIJII IIJ-11 . It'gives a good indica--n /g i = N, N-1, ..., 2 (7) tion of the sensitivity of [JI-to small perturbation 1. N4+1,X in [JM. It is also called the spectral condition nunbbe because of its dependence on the spectral norm with' g = f1 and L1 = xl, so that gN is just a function u 1 1*. of the single variable xl. Now expanding, linearizing Thle condition number of a symmetric positive defi-and solving for x we obtain nite matrix [JTJ], whose eigenvalues are all real and x n n (8) positive can be computed by, xl x gN,x ) Ix I k([JTJ]) max (1) The point xl thus obtained is used as the next approximmi mation xkhl to the first component xl of the solution -vec'tor x* li6nce, is renam'ed as-.L, and the t ys The condition nunber of [J] is the square root of o H e v 1 i a sy. k([JTJ]). For a well-conditioned system, the value of tem of eqn (7) is back-solved to get improved approxi- is 1. A very high value of k indicates that the sys-mations to the other compoents of x*. Hre is tem is ill-conditioned. If k exceeds 0a, where a taken as the value obtained.for L. when back-s6lving equals decimal precision of the digital computer used, eqn., (7). TIhe "successive substitu%ion!' nature of the it is not possible to obtain a solution [11]. Solutions algorithm allows the mst recent information available to eqn. (1) are. given in Table I for the AEP 30-bus to. be used in the construction of next function arguweellconditioned and 11-, , and 43-bus ill-condi-ment, similar to the Gauss-Seidel process for linear tioned test systems. and nonlinear systems of equations. K.M. BROWN'S METHOD M pe a -. -----, .~ M-triM reetatim . Consider the following real continuously differentiaabl system of N nonlinear equations in N unknowns, For the sake of -definiteness, it may be stated x. (i = 1, 2, ..., N). In vector notation that the variables are eliminated in the.order xN, f. (x) = 0 ; i-= 1, 2, .., N (2) xN, T, x2. Using thechain rule for differentia-,-1 -* * ~~~~~~tionto expand each derivative g. , gives the follow-K.M. Brown proposed a. local method [7,8] which ing matrix representation for tlieXforward part of the handles the functions of eqn. (2) one at a time so that method: information obtained from working with f1 can be. incor--n+ n porated. when working with f2, etc. A successive sub-[H] (x -x ) = -[g] (9) stitution scheme is used rather thanl the. simTultaneous treatment of the f. which is the characteristic of New-where the muatrix [H] = (h..) is given by ton's method. It is a quadratically convergent local technique (in the vicinity of a root). It is fast and3650 TABLE I Maximum and Minimum Eigenvalues and Condition Number Type of No. of Maximum Minimum Condition Re ]rs System Buses Eigenvalue Eigenvalue Number (k) * emar WVell-conditioned 3 0 X13 Moeael el syteml-conditione | 30 .1087xl1 .2322x10 .4681x10 deratelywll-system ~~~~~~~~~~~~~~~conditio'ned (fair k) 11 .1222xlO3 .1126x100 .1086xl04 Ill-conditioned (bad k) 2 14 Ill-conditioned 13 .2905x10 .1442xlO-1 .2014x10 Ill-conditioned system (bad k) 43 .2426xl04 .9476xl0_1 . 2560x105 Ill-conditioned (bad k) * Ideal value of condition number: k = 1 f ~~~~~~j=1,..,N flj f f lInitialize all the variablesll 1,N-i+2 1N f ff -Statieainc-_t> 2j 2,N-i+2 2N S h.. = i+l Ifij f ' f | i = 2, i, N | Create an array which permits (-l) ~~~~~~~~~~~~~~a partial pivoting ef'ect witi 1,N-i J 1 N out having to physically inf11N-2 fi Lerchange rows and columns 1-l1,N-i+2 i-1,N D , . ..(10) | Jwhere K -runction numberl where the argunents f.j are progressive arguments gen-+ ve erated successively anA [g] is a vector-valued function Kgiive by eqn. (10). It is observed that hi; = 0 for j > N-i+l, that is, the matrix [H] is transverse upper = KMIN . K-I triangular. The matrix [H] can be obtained by transforrmin the Jacobian matrix [J] into transverse upper Solue the first KMIN rows of triangular form using Gaussian elimination with partial a triangularized linear systpivooting em for improved values ofe x in terms of' previous one*s When the condition number of the matrix [H] << the I condition number of the Jacobian [J] convergence occuurs Thus Brown's method gives convergence in cases where Newton's method fails. vet up approximate partial derviatiee of the kth function and find the Computational Efficiency of Brown's Method largest absolute alue or thp a A count of the number of function values of the fi YES i e l~~argest absolut YS needed per iteration of Brown's method, is given in Ap-apro.artial der pendix 1. The first step requires N+l evaluations of vprox-partial deO fl, the second step N evaluations of f2, the third N-l in Try a different i evaluations of f3, etc., so that the total is NO tial .approximation N+1 1 2 Set up coefricients for kth ro i2 = 12~N+3)(1 or' triangular linear system usii2 ed to back solve for the f'irst Newton's. method requires N partial derivative evaluations and N function component evaluations per B xcksubstitute to obtain iterative step. Thus there is a corresponding savings in storage locations required: from (N2+N) for Newton's method to -I(N +3N) locations for Brown's method. How-YES__j ever, it must be stressed that here the savings in function values applies only to saving function values of the fi of the original system (2), because Brown's method adds a number of other fuctions to be evaluated Fig. 1. Simnplified flow-chart for Brow's method namely the linear functions, L}{. The evaluation of the Lkr o included in the cout above, 1(N2+N show that Brow's method also falls into the general L ar no '-N+3) category of successive approximation methods like the The equnceof omptatonsis how ina smpl Gass-Seidel method and Newton-Raphson method, differfllowchart for Brown's method in Fig. 1 ing only slightly in their iteration fuction. The -~~~~~~~solution of a vector set of equations f. (x) = 0 iS FIXED-POINT FORMUATION sought by means of a suitable recursive Tormla [14] put into the form: The fixed-point formulation of Ref. 13 is used to3651 xn+l = p(xn) (12) same. The stopping criterion is when all the function values are less than the prespecified where the mapping 1p(xn) is given by tolerance. n(x) = xn [Rn] f(xn) (13) 5. Calculate all bus powers and line flows, and --print. By a successive approximation procedure, solution to the fixed-point formulation is given by xn+l = xn -[R'] h(xn, xn+l) n = 0 1 2,... (14) Nodes are numbered so that at each step of the ---' ' ' ' ' Gaussian elimination the next node to be eliminated is where [Rn] is a square matrix and h(.,.) is a continu-the one having the least number of nonzero elements. ous vector-valued fumction. This is generally preferred for Newton-Raphson method [3,16]. The same strategy is used,for Brown's method The Gauss-Seidel method uses a simple diagonal ma-also. Henice, the above procedure can be easily incortrri [R] having scalar acceleration factors for the porated into the existing Newton-Raphson load-flow proreea and imaginary parts of bus voltages. For Newton's grams. method matrix [R] corresponds to the inverse of the Jacobian matrix, whereas in Brown's method it corres-c) Results and CIarison: Line diagrams of 11-, 13-, ponds to the inverse of Brown's Jacobian matrix [H] and 3-us ill-conditioned test systems are shown in given in eqn. (9). Figs. 2, 3 and 4 respectively. The system data of these ill-conditioned systems are given in Appendix II. Load LOAD-FLOW BY BROWN'S METHOD and generation data are shown only for 13-bus system, whereas, for 11-, and 43-bus systems only injection inaa Problem formulation: The load-flow problem involves formation is available. The system in Fig. 3 is diffitth solution of a system of nonlinear algebraic equa-cult to solve [6] because of the two series capacitors tions, that is and the position of the slack-generator. The systems in Figs. 2 and 4 are also difficult to solve [4,12] beff (x) = 0 * i = 1, 2, ..., N (2) cause of low X/R ratios and some negative line reactannces This is also clear from Table I which gives the For the solution to be unique, one equation must be eigenvalue analysis and condition number of the Jacob-Forcifiedsolution for each unknwni e eqaThe following ian matrix [J] of the load-flow equations. The results specified for each unknown variable, h oloig in Table I show that the condition number of the ill-equations define more clearly the load-flow problem on i oned s t a conditi smuch highe in the basis of the bus admittance matrix. Corresponding conditioned systems Jacobian matrix is much higher in to fi in eqn. (2), the load-flow mismatch equations comparison to the well-conditioned systems. equal to the number of unknown voltage phase-angles and Slack-node magnitudes are given by 1 2 Real Power Mismatch: APi. =P. S i(Cal)=0 (15) ( I i(Sch) icl Reactive Power 7 Mismatch: AQi = Qi(Sch) Qi(Cal) (16) The equations for the calculated real and reactive bus powers at bus i are given by n ( P. = I V.V.(G.. cos(6. -0.) + B.. sin(e. -O.)) 1j= 1 1 3 13 1 3 1 J 1 * * (17) 11 n Qi = l V.V.(G.. sin(0. -0.) -Bi; cos(0. -0.)) Fig. 2. 11 Bus ill-conditioned power system P a 500 Slack Pa O (18) 0a ? generotor 0.? The total number of nonlinear algebraic load-flow equations (15) and (16) that are to be solved by Paso 650 Brown's method is equal to: (n -1) real power mismatch equations and (n -1 -number of PV buses) reactive Pao power mismatch equations. b) Solution algorithm: Equations (15) and (16) are to be solved alternately as in [3], for each bus by K.M. Brown's method as follows: 1. Read network line data, generation and load data, and form the Y-bus matrix elements according to Piaso _2 the optimal ordering described below. .2 2. Formulate the system of equations for real and re-P'O Pa 500 active bus power mismatch equations (16 and 17) by sa? using optimal ordering. 3. Initialize the unknown vector xn with alternate Q a30' voltage phase-angles and voltage miagnitudes, i.e., voltage magnitudes to 1.0 p.u. and angles to 0°. Fig. 3. 13 Bus ill-conditioned power system 4. Solve the bus power mismatch equations for the new estimates xl+± by Brown's method and update the3652 Slack-.node A comparison of the proposed Brown's load-flow method with the standard methods [1,3,4,5] in terms of convergence characteristics is given in Table II. The proposed Brown's method converges in fewer iterations for both ill-and well-conditioned systems. A comparison for the well-conditioned systems betwwee Brown's load-flow method and the standard Newtoon' method in terms of ratio of computing time per 12 23 16D1 (9 (!) ( iteration is given in Table III. Since the Newton-Raphson method does not work for any of the 11-, 13-, and 43-bus ill-conditioned test systems, comparison resuult are not given for this method. Digital computer 22 28 37 31 @results were obtained by using non-optimized Fortran compilation on ICL 2960 computer with George 2+ operatiin systems. Even though Brown's method takes 15% more 21) 30 38 33) ( g time per iteration than Newton's method, for large systeem the total number of iterations is less. This makes Brown's load-flow method comparable to Newton's 36 -\9method in total computing time. -The storage requirement of Brown's method is slightly more (not exceeding 10%) than that of Newton's method. However, this disadvantage is offset by the superiority of Brown's method for ill-conditioned systeems CONCLUSIONS The main objective of this paper is to show that the K.M. Brown technique is suitable for solution of 3 13 18 ill-conditioned power systems equations. The technique has quadratic convergence characteristics. The pro-Fig. 4. 43 Bus ill-conditioned power system cedure adopted in this method is relatively simple and can be easily incorporated in the existing Newton-Raph-The proposed Brown's load-flow method has been ap-son algorithm. The optimal ordering strategy and plied to the ill-conditioned test systems and AEP 14-, sparsity programming technique with compact storage 30-, and 57-bus well-conditioned test systems. The re-scheme can be easily implemented. Also transformers sults are shown in Fig. 5 for the convergence rate of with no-load tap changers and phase shifters can be the largest absolute mismatch for ill-conditioned as taken into account while developing production type well as for well-conditioned test systems. It can be load-flow programs. It is not necessary to include observed that the convergence rate of the absolute max-transformer resistance in the algorithm to help the imum mismatch is quite rapid. The proposed method con-load-flow to converge to a solution. Brown's load-flow verges close to the solution on the first iteration and algorithm converges to the solution of ill-conditioned requires few iterations to obtain the final solution. system in a few iterations, whereas the standard loadfllo methods either show poor convergence or diverge. 11 bus system ACKNOWLEDGEMNTS 13 bus system I ll -cond; t ioned -3 bus system The authors would like to acknowledge Dr. S. 1-bs-Iwamoto of the Department of Electrical Engineering, \1I. bus system Tokai University, Japan for providing data for 11-and 30 bus syslem Well-conditioned 43-bus ill-conditioned systems. -101-\ *--57 bus system t REFERENCES X o°X0 [1] A.F. Glimn, and G.W. Stagg, "Automatic calculation of load-flows", AIEE Trans., 1957, PAS-76, pp. 817-825. [2] A. Brameller, and J.K. Demnead, "Some improved methods for digital network analysis", Proc. IEE, l0.2 1962, 109, A, pp. 109-116. w '<>--\ \ \[3] W.F. Tinney, and C.E. Hart, "Power flow solution A. iJ3L by Newton's method"v, I EEE Trans., 1967, PAS-86, : \\ X ~~~~~~pp. 1449-1460. 10.4l 2 3 4 5 [4] B. Stott, and 0. Alsac, "'Fast decoupled load-flow" ITERATIONS IEEE Trans., 1974, PAS-93, pp. 859-869. [5] A.M. Sasson, "tNonlinear programming solutions for Fig. 5. Gonvergence rate of Brown' s load-flow method load-flow, minimum loss and economic dispatch for well-conditioned and ill-conditioned problems", IEEE Trans., 1969, PAS-88, pp. 399-409. systems3653 TABLE II Convergence Characteristics of Different Methods (Number of iterations) System No. of Gauss-Newton-Fast Fletcher Brown's Type Buses Seidel Raphson Decoupled (Steps) Method Ill-11 divergent divergent divergent divergent 5 conditioned 13 divergent divergent divergent 47 4 systems 43 divergent divergent 22 divergent 5 Wlell-14 24 3 4 32 2 conditioned 30 33 3 4 68 3 systems 57 59 4 41 135 3 TABLE III tems of equations", (Academic Press, New York, 1966). Comparison of Computing Time per Iteration [15] R. Fletcher, and M.J.D. Powell, "A rapidly convergeen descent method for minimization", Computer One iteration time in terms Journal, Vol. 6, pp. 163-168, 1963. System No. of ~~~~~of ratios System obuss (Newton's method as 1 unit) [16] W.F. Tinney and J.W. Walker, "Direct solutions of Type buses sparse network equations by optimally ordered tri-Newton's method Brown's method angular factorization", Proc. IEE, Vol. 55, pp. 1801-1809, Nov. 1967. Well-14 1.0 1.17 conditioned 30 1.0 1.12 [17] S. Iwamoto, Personal communication with the ausysstem 57 1.0 1.16 APPENDIX I K.M. Brown's Method [6] K. Zollenkopf, "Load-flow calculation using loss minimization techniques", Proc. IEE, 1968, 115(1), A step by step synthesis of Brown's method [7] of pp. 121-127. solving the system of nonlinear algebraic equations is as follows: [7] K.M. Brown, "A quadratically convergent Newtonliik method based upon Gaussian elimination", SIAM Let the system be in vector notation as J. Numer. Anal., 1969, pp. 560-569. f.(x)=0*i=1,.,N(19) [8] K.M. Brown, "Computer oriented algorithms for f(x) 0 ;i-1 2, solving systems of simultaneous nonlinear alge-Step 1 braic equations", in G.D. BRYNE and C.A. HALL, (eds.), "Nunerical solutions of systems of nonlin-Let x' denote an approximation to the solution x* ear algebraic equations", (Academic Press, 1973), of (9). Eipand the first function in the Taylor series pp. 281-348. expansion about point _,. Alternatively f, can be expannde in an approximate Taylor series with the actual [9] L.L. Freris, and A.M. Sasson, "Investigation of partial derivatives replaced by first difference quotth load-flow problem", Proc. IEE, 1968, 115, (10) tient approximations. In the later case the. user does pp. 1459-1470. not have to provide the partial derivative expressions in his program. Retaining only first order terms gives [10] G.W. Stagg, and A.E. El-Abiad, "Computer methods a first order approximation in power system a;nalysis", (McGraw-Hill, 1968). 1 (x) e f1 (x ) + f (xn) (x -xn) + f ( n) rx-x) [11] J.H. Wlilkinson, "The algebraic eigenvalue prob-2 2 led'l, (Clarendon P-ress, Oxford, 1965). [12] S. Iwamoto, and Y. Tamura, "A load-flow calcula-l +fx (x%(x -xE) (20) tion method for ill-conditioned power systems"', Npressente at the IEEE PES Summer Meeting, Vanicouuver Canada, 15-20th July, 1979, paper A79 441-af. f.(x'n + -n) _ i = , ,.., ~~~~~~~~f [13] J. Meisel, and R.D. Barnard,. "Application of ) ) h~ . .. (21) fixed-point techniques to load-flow studies", IEEE Trans., 1970, PAS-89, pp. 136-140. where ej denotes the jth unit vector and the scalar hn is normally chosen such that hn = 0(1 If(xn) |). With [14] A.M. Ostrowski, "Solutions of equations and sys-this choice it can be proven 1181 that the discrete3654 Newton's method has second order convergence. This is the Ws are obtained by back-substitution in the N-1 the same rate of convergence as the ordinary Newton's rowed Itriangular linear system which now has the form method. If xn is close enough to x*, fl(x ) Owe can L= -I /gf Nl )(Lg -xj, equate (70) to zero and solve for the ffrst variable. J=1 i This is the variable, say XN, whose corresponding par--g tial derivative, f (xD, is largest in absolute val-/g-NN+l,x.2 5 ue. This gives lxN J, N-1 n with g-f1 and L1 = xl so that gn is just a function XN . (fXNflx )(Xj-X. X f/f (22) of the single variable xl. Now expanding, linearizing j=1 j N N and solving for xl, we obtain The constants f'1 /f'n1 are saved for later use. Brown = -'n/a (26) lxi lxN X, i1 glg, 2) [7] has shown that under the usual hypotheses for New-1 ton's method there is always at least one non-zero par-Thepoint x thus obtained is used as the next approxi-tial derivative, and, of course, the corresponding ap-T proximate partial (21) is non-vanishing. Thus the sol-mation, x'1+ to the first coponent, xl, of the soluuttio procedure (22) is well defined. By choosing the tion vector x*. Renane xl as L1 and back solve the Li approximate partial derivative of largest absolute val-system (25) to get an improved approximation of' comonuu as the division, a partial pivoting effect is a-ents of x*. Here x4+1 is taken as the value obtained chieved similar to what is often done when using the for L. wheii back solving (25). Gaussian elimination process for solving linear sys-n tems. This enhances the numerical stability of the Choice of h in the Computer Implementation method. We observe from (22) that xN is a linear functiio of the N-1 variables xl, x2, ..., XN X . For pur-In order to guarantee quadratic convergence, the poses of clarity the right-hand side of (221 is renamed following strategy is used for choosing hM, by which to as L.(x1, x2, *.., xNl) increment x' when working with the i finction fi in eqn. (21): Step 2 hn7 = maxCn. 5 x 10-6+2) Define a new function g2 of the N-1 variables xl, 1) .., xN_1 which is related to the second function f2 Of where (19) as follows { fln n2^***^|gn| ) ; *001 x xn |} 2(X **..., XN-1) f2(X1i, ..., XN-1 Ln(xi, 1' ' n fN-1m ..(27) **(23) where S is the machine tolerance, i.e. the nunber of Now expand g2 in an approximate Taylor series expansion significant digits carried by the machine. about the point (x, ..., N_). Linearize (ignore higher order terms) and solve for that variable, say EXAMPLE: Application of K.M. Brown's method xN_, whose corresponding partial derivative, x is largest in magnitude: N-1 Let us consider a nonlinear system with two var-N-2 iables: xN-= ~N-11 j -j=1 22X/2x ) -) 92/2x f (X)=X -2x + 1 = ° j1 ~~~~~~~~~~~~~~~~~~~~~~~~~~~(... (24) f (x) = x + 2x -3 = O X_l is a linear function of the remaining N-2 vari-2( 1 2 ables. Remember the right-hand side of (24) is denoted which has the solution at [1,1]T. For x° = [0,0]T as by Li (xl, ...`, xN2). Again this forms the ratio, the initial approximation, the application procedure is N2NI/gteration 1, ..., N-2, and gn/gn is stored demonstrated here using exact partial derivatives. 92x 2N-1 922N-1Ieato for future use in any camputer inplementation of this algorithm. Each step of the algorithm adds one more linear expression to a linear system. During the 0 Expand f1(x) i the Taylor series about the poit (k+l)st step of the algorithm, it is necessary to eval-x; retain only first-order terms and thus obtain the uate gk+l, i.e fk1 for various arguments. The values Iinear approximation: of the last k components of the argument of fk+1 are 0 + x + x obtained by back-substitution in the linear system LN, f1(X) f1(X + f1X (Xi X0 + f1X (X2 XP LN-1, .. ., 'N-k+1 which has been built up. These argu-2 ... (2) ments are required to determine the quantities gktl and where ,n+l , j = 1, ..., N-k, needed for the elimination of the (k+l)st variable, say xN..k, by the basic process of 1-' lx x2 expansion, linearization and solution of the resulting teRSo q.(9 ozr n ov o h expression. The process results, for each k, in the Equ.ate te oeq..2)t eoad ov o h (k+l)st variable, say ._,i expressed as a linear variable x2, whose corresponding partial derivative is combination, L-k of th~ekremaining N-k+l variables, largest in absolute value: Step N xl = 0.5 At this stage gN ENf(x1 L2, L3, ... LN) where Substituting this value into function f2 (x), we have3655 Substitute this value in function f2(x) in (28) and redeffin as g2(x): Thus2-xl = 2.0 g2(x) = 2.5x2 -1.95 Expanding the function g2(x) into the Taylor series 2 9g2(x) = g2( + g2x (x2 -x2)=j (31) 2 Iteration 211 where g(x ) = -1.0 gl = 2.5 By repeating the samne procedure with new estimates 2 1(X 4.0 f1 4.0 f1 2.0 Solving for the variable x2 from (31) gives the new f1(x ) 4.0 ; f = 4.0 ; f -2.0lx2 estimate after the second iteration 2 Substituting in x1 =1.20 f (x) = f(x1) + f (x -xl + f x2 2 X 22 f1x 1 f1(x f 1-fx -2x20)9 1;nd equatingho2 en(30)x( = and equating the RHS of eqn. (30) to ze-ro gives By repeating this same procedure the problem converged x= 0.5x2 + 0.75 in five iterations to xi = 1.0 and x2 = 1.0. APPENDIX II A2.2 Description of the 13 bus ill-conditioned system System Data for Ill-conditioned Systems The line and bus data has been taken from [6]. A2.1 Description of the 11 bus ill-conditioned system TABLE A2.3 The line and bus data has been taken from [17] and Line data for 13 bus ill-conditioned system is available only in the form of Y-bus matrix elements and net bus powers. Suscep-Branch From To Resistance Reactance tance TABLE A2.1 number node node (p.u.) (p.u.) (pu) Y-bus matrix elements of 11 bus system 1 0 1 0.0040 0.0850 0 2 0 2 0.0040 0.0947 0 k -mGi B k -m G B 3 4 3 0.0040 0.0947 0 km km km km 4 3 2 0.0074 0.1430 0.436 1 1 0.0 -14.939 5 5 2.581 -5.889 55 1 0.0481 0.4590 0.246 1 2 0.0 14.148 6 6 0.0 -55.556 6 5 6 0.,0090 0.1080 0.016 2 2 12.051 -33.089 7 7 3.226 -4.304 7 2 0.0121 0.2330 0.712 4 7 -0.592 0.786 11 11 0.283 -2.785 Base 6 7-01000 MVA 2 3 0.0 6.494 7 8 -2.213T 2.959870 0 2 4 -12.051 13.197 8 8 2.893 -5.468 9 8 9 0.0105 0.2020 0.620 3 3 2.581 -10.282 8 9 -0.138 1.379 10 9 10 0 -0.1500 0 3 5 -2.581 3.789 8 10 -0.851 1.163 11 10 11 0.,0086 0.1665 0.508 Operating~~ ~ ~ ~ ~ ~ ~ ~~ ~1codiio 12 11070bus6 systm4-8 4 4 12.642 -74.081 9 9 0.104 -1.042 12 11 12 0075 .16 048 4 5 0.0 2.177 10 10 1.346 -6.110 10 0 4 6 0.0 56.689 10 11 -0.374 3.742 6 -0. 090 -0 .068~~Bas 100 W 4 7 -0.592 0.786 11 11 0.283 -2.785 TABLE A2.4 TABLE A2.2 Transformer data for 13 bus system Operating condition of 11 bus system__________________ Branch From TO Tap Bus Voltage Phase-Net real. Net reactive number node node setting No. magnitude angle power powe-r V (p.u.) e (deg.) P (p.u.) Q (p.u.) 1 0 1 + 5%90 2 0 2 +10% 1 1.024* 0.0* 003 4 +10% 3 -0.1028 -0.062 4 0.0 0.0 3 -0.128 -0.080 6 -0.090 -0.068 *Slack-bus input data3656 TABLE A2. 5 TABLE A2. 7 Bus data for 13 bus system Operating condition of 43 b-us system Assumed Assumed Generation Load Busbus phasue-Generation Lea Bus Voltage Phase-Net real Net reactive Bus voltage agle pG QG PL QL No. magnitude angle power power V (p.u.) e (deg.) 4) (4VAr)) (M) (WAr) * V (p.u.) e (deg.) P (p.u.) (P.U.) 0 1. 0* 0. 0* --1650 560 1 1. 136* 0. 0* --1 1.0 0.0 0 0 0 0 2 0.0 0.0 2 '1.0 0.0 0 0 0 0 3 -0.160 -0.120 3 1.0 0.0 0 0 0 0 4 0.0 0.0 4 1.0* 0.0 0 -0 0 5 -0.530 -0.400 5 1.037* 0.0 500 -50 30 6 0.0 0.0 6 1.063 0.0 0 0 0 0 7 -1.600 -1.200 7 1.100* 0.0 0 -0 0 8 0.0 0.0 8 0.943* 0.0 500 -0 0 9 0.0 0.0 9 1.100* 0.0 0 -0 0 10 0.0 0.0 10 1.0 0.0 0 0 50 30 11 0.0 0.0 11 1.0 0.0 0 0 50 32 12 -0.800 -0.600 12 1.0 0.0 0 0 0 0 13 0.0 0.0 *Input data 14 -0.800 -0.600 15 0.0 0.0O A2.3 Description of the 43 bus ill-conditioned system 16 -0.640 -0.480 The line and bus data has been taken from [17] and 17 0.0 0.0 is available only in the form of Y-bus matrix elements 18 -0.240 -0.180 and net bus powers. 19 0.0 0.0 20 -0.880 -0.660 TABLE A2.6 21 0.0 0.0 Y-bus matrix elements of 43 bus system 22 0.0 0.0 :_____________________________________________ _ . -23 0.0 0.0 k -m Gk Bk k-G B 24 -0.640 -0.480 .km~ -ki 25 0.0 0.0 1 1 0.0 -30.609 1 2 0.0 30.609 26 -0.800 -0.600 2 2 481.288 -1545.194 2 5 -277.195 873.583 27 -0.320 -0.240 2 6 -34.368 108.124 2 15 -169.726 534.322 28 0.0 0.0 3 3 0.0 -5.714 3 4 0.0-6.015 29 0.0 0.0 4 4 61.331 -69.160 4 13 -61.331 62.874 30 0.0. 0.0 5 5 277.195 -916.892 5 7 0.0 21.277 31 1.160 0.520 5 8 0.0 20.513 6 6 34.368 -118.699 32 2.900 2.570 6 12 0.0 10.638 7 7 0.0 -20.000 33 0.285 0.300 8 8 452.840 -482.861 8 9 -288.938 295.777 34 0.0 0.0 8 23 -163.902 167.191 9 9 300.983 -317.044 35 0.580 0.560 9 10 -12.045 12.342 9 16 0.0 8.796 36 -0.005 0.030 10 10 12.045 -20.855 10 11 0.0 2.857 37 0.0 0.0 10 17 0.0 5.714 11 11 0.0 -2.857 38 -1.440 -1.020 12 12 0.0 -10.000 13 13 92.381 -100.709 39 0.0 0.0 13 18 0.0 6.015 13 25 -31.050 31.640 40 0.0 0.0 14 14 0.0 -15.015 14 43 0.0 15.400 41 -0.800 -0.300 15 15 340.398 -916.783 15 19 0.0 8.649 42 -2.240 -1.680 15 20 0.0 15.791 15 28 -170.673 357.003 43 0.0 0.0 16 16 0.0 -8.576 17 17 0.0 -5.714 18 18 0.0a -5.714 19 19 164.292 -280.783 *Slack-bus input data 19 22 -164.292 272.805 20 20 0.0 -15.002 21 21 104.312 -143.609 21 24 0.0 9.267 21 29 -104.312 133.623 22 22 164.292 -282.281 Dr. S.C. Tripathy, Professor of Electrical Engin-22 26 0.0 9.023 23 23 321.579 -328.810 eering at I.I.T. Delhi, India, is presently a visiting 23 29 -157.677 161.760 24 24 0.0 -8.572 professor at The University of Calgary. His research 25 25 87.150 -106.814 25 27 0.0 9.023 interests are in power system analysis and control. 25 29 -56.100 65.824 26 26 0.0 -8.572 27 27 0.0 -8.572 28 28 373.447 -612.837 r. G. Dlurga Prasad is a research scholar at I.I.T. 28 39 -202.775 256.136 29 29 318.089 -372.311 DelTli, with interests in power system analysis. 29 30 0.0 3.766 29 37 0.0 7.895 30 30 125.789 -524.464 30 32 0.0 30.769 'Dr. O.P. Malik, SMIEEE, Professor of Electrical 30 38 0.0 4.131 30 40 -125.789 485.547 Engineering afid Associate Dean -Academic, The Univer-31 31 0.0 -13.038 31 37 0.0 13.038 sity of Calgary is interested in the real-time digital 32 32 0.0 -30.769 33 33 0.0 -3.320 control of electrical machines and power systems. 33 38 0.0 3.320 34-34 0.0 -7.365 34 38 0.0 6.852 35 35 0.0 -6.180 Dr. G.S. Hope, SIEBE, Professor of Electrical En-35 38 0.0 6.180 36 36 0.0 -2.703 gineering, Thie University of Calgary, has research in-36 38 0.0 2.703 37 37 0.0 -21.348 terests in the area of digital systems, anld real-time 38 38 0.0 -22.398 39 39 512.581 -663.260 control and protection of power systems. 39 41 0.0 15.015 39 43 -309.806 392.255 40 40 125.789 -508.837 40 42 0.0 21.622 41 41 0.0 -15.015 42 42 0.0 -20.000 43 43 309.806 -408.0293657 Discussions As regards 11 bus system results. The discussers have the following observations to make: J. Nanda, D. P. Kothari and D. L. Shenoy (Indian Institute of 1) 11 bus test system as reported by the authors pertains to an Technology, New Delhi,. India): The authors are to be commended for Australian Distribution System originally reported by Mr. B. Stott. establishing their stand that K. M. Brown's method is reliable for solv-2) For the same 11 bus system the discussers made several tests which ing ill-conditioned power systems. revealed that this system has converged by Newton-Raphson (rec-Authors' definition and consequently the condition they impose for tangular co-ordinate version) with a flat voltage start (1 +jo) for all "condition number" for ill-conditioning of computer formulation of a load buses in 6 iterations within an accuracy of 0.01 MW/MVAR, while problem are quite interesting indeed. It is stated in the text of the paper B. Stott's FDLF failed (1]. The results-for the 11 bus system (obtained that for a well conditioned system, the value of the condition number K by Newton-Raphson method) are given in Table-1. Could the authors is unity. However, in Table I which shows the maximum and minimum please clarify why they could not get convergence? eigen values and the condition number K for well-conditioned 30 bus and ill-conditioned 11.13 and 43 bus systems. The condition number for the 30 bus test system is given as 0.468 x 103, which is quite large. Authors have indicated in the remarks column this number K as fair. Would the authors indicate the ranges in which the condition number could be said to be good (well-conditioned), fair (moderately well-REFERENCE conditioned) and bad (ill-conditioned systems)? [1] P. G. Murthy, D. L. Shenoy, J. Nanda and D. P. Kothari, "Perfor-Table I mance of Typical Power Flow Algorithms with reference to Indian 11 Bus System Solution as obtained by Newton-Raphson (rectangular Power Systems", presented at Second Symposium on Power Piant co-ordinate version) method. Dynamics and Control, Bharat Heavy Electricals Ltd., Hyderabad, India, February, 1979. Load in Voltage obtained Manuscript received February 18, 1982. MW MVAR Mag Angle in in p.u. degree i --1.024 0.000(slack bus) 2 0.0 0.0 1.0552 1.215 S. C. Tripathy, G. D. Prasad, O. P. Malik, and G. S. Hope: The 3 12.8 6.2 1.0444 4.608 authors are thankful for a stimulating discussion on the present paper. 4 0.0 0.0 1.0306 2.830 The condition number k is a pointer to the extent of separation of 5 16.5 8.0 1.0329 4.830 system eigenvalues and ideally its value is one if the minimum and max-6 9.0 6.8 1.0501 2.916 imum values are -equal. It gives a qualitative assessment of the ill-7 0.0 0.0 0.8310 12.066 conditioning effect, so it is improper to define the range of its values 8 0.0 0.0 0.9131 14.857 quantitatively. Regarding the second question on the solution of the 9 2.6 0.9 1.1264 15.911 1 1-bus system, the results of the angle in degrees given in Table 1 by the 10 0.0 0.0 0.8266 21.033 discussers are in doubt because of their positive signs. 11 15.8 5.7 1.0186 23.996 Manuscript received April 15, 1982.