Research Article  Open Access
Convergence Region of Newton Iterative Power Flow Method: Numerical Studies
Abstract
Power flow study plays a fundamental role in the process of power system operation and planning. Of the several methods in commercial power flow package, the NewtonRaphson (NR) method is the most popular one. In this paper, we numerically study the convergence region of each power flow solution under the NR method. This study of convergence region provides insights of the complexity of the NR method in finding power flow solutions. Our numerical studies confirm that the convergence region of NR method has a fractal boundary and find that this fractal boundary of convergence regions persists under different loading conditions. In addition, the convergence regions of NR method for power flow equations with different nonlinear load models are also fractal. This fractal property highlights the importance of choosing initial guesses since a small variation of an initial guess near the convergence boundary leads to two different power flow solutions. One vital variation of Newton method popular in power industry is the fast decoupled power flow method whose convergence region is also numerically studied on an IEEE 14bus test system which is of 22dimensional in state space.
1. Introduction
Since power flow calculation started in 1956 [1–3], a variety of numerical methods have been developed for power flow study, such as Gauss method [1], Newton method [4], and fast decoupled method [5]. Over the past 20 years, Newton method and its variations are widespread, utilized, and enhanced [6–11]. Power flow study (or load flow study) is the determination of steadystate conditions of a power system for a specified power generation and load demand. The power flow problem is the computation of voltage magnitude and phase angle at each bus in a power system under the following conditions:(i)balanced threephase steadystate conditions and the system is in sinusoidal steadystate;(ii)the transmission network is composed of constant, linear, and lumpedparameter branches;(iii)the specified real and reactive power demand at each (load) bus;(iv)the specified real power generation at each (generator) bus except one generator bus.
Power flow study is performed extensively both for system planning purposes to analyze alternative plans of future systems and for system operation purposes to evaluate different operating conditions of existing systems. Indeed, power flow study is used in transmission planning to check for branch overloads, bus voltage problems. In static contingency analysis, power flow study is used to assess the effect of branch and/or generator outages. In transfer capability analysis, (repetitive) power flow study is performed to calculate the power transfer limits.
In the literature, there was little work on the study of convergence regions of power flow solutions. The existing relevant papers that the authors were able to identify were mainly focused on the following aspects: power flow fractals and truncated fractals on a 3bus system [12–17], different convergence regions under polar and rectangular expresses Newton method [18]. This study of convergence regions provides insights to explain the reasons behind power flow divergence and may lead to improving the robustness of Newton power flow via some parameters adjustment.
In this paper, we numerically study the convergence region of each power flow solution under the NewtonRaphson (NR) method. One important variation of Newton method popular in power industry is the fast decoupled power flow method. We also numerically study the convergence region of fast decoupled power flow method.
The key results presented in this paper are summarized as follows.(1)Our numerical studies on a larger system, IEEE 14bus test system, as compared with the 3bus power system, confirm that the convergence region of NR method has a fractal boundary.(2)The convergence region shrinks with the increase in loading conditions. This property can be used to explain the difficulty of Newton method in finding a power flow solution under heavy loading conditions, in addition to the numerical illcondition.(3)The fractal boundary of convergence regions persist under different loading conditions (the lightly loaded condition, the medium loaded condition, and the heavily loaded condition).(4)Convergence regions of NR method for power flow equations with different load models all reveal fractal.(5)Convergence region of fast decoupled power flow method also has a fractal boundary.(6)The required number of iterations is large (i.e., greater than 10) when the initial guesses are located near the boundaries of convergence regions. It is interesting to note that the region of convergence with different required number of iterations is disconnected.
This fractal property highlights the importance of choosing initial guesses since a small variation of an initial guess near the convergence boundary leads to two different power flow solutions. This fractal property persists under different loading conditions. This study of convergence region provides insights of the complexity of the NR method as well as the fast decoupled power flow method in finding power flow solutions even for largescale power systems. This study also sheds lights on the difficulty of finding a power flow solution of heavilyloaded power systems.
2. Convergence Region of Newton Method
Newton method is a locally convergent iterative method for solving a set of nonlinear equations, which can be expressed as Here .
The Newton iteration can be expressed as where : state variables of the th iteration; : state variables of the ()th iteration; : error between and ; : Jacobian matrix of .
If the components of are differentiable, we define the Jacobian matrix in (2) as
It is well known that Newton method is locally convergent under the following assumptions [19].
Assumption 1. Equation (1) has a solution .
Assumption 2. is Lipschitz continuous near .
Assumption 3. is nonsingular.
The classic convergence theorem is stated as follows.
Theorem 4 ((local convergence) [19]). If Assumption 1 through Assumption 3 hold and the initial guess is sufficiently near , then the Newton sequence exists (i.e., is nonsingular for all ) and converges to and there is such that for sufficiently large, where the mismatch errors decrease quadratically.
Assumption 3 assures that the singular Jacobian matrix does not appear during iterations. By the same token, the initial point and the singularity of Jacobian matrix do make great effect on the convergence of Newton method. Starting from an initial guess , the algebraic equations (1) are iteratively solved. Generally, the Newton method converges within a reasonable number. However, the Newton method may diverge. To this end, the study of convergence region is relevant.
2.1. Fractal Convergence Boundaries
The convergence region of a power flow solution under the Newton method is the region in the state space such that initial points starting from this region can nicely converge, under the Newton method, to a power flow solution. It is known that the boundaries of convergence regions of power flow solutions obtained by Newton method are a set that is analogous to the “Julia set” [15–18]. These boundaries are known to have fractal geometric features. Using initial points close to any of the solutions would result in a rapid convergence. However, using initial points near the boundary of a convergence region, the Newton method becomes unpredictable. In this study, a region in state space in which the initial points do not converge (in a specified tolerance) to any power flow solution is termed the divergence region of Newton method.
We consider a bus power system with m generators. There exits three different types of buses: slack bus 1, bus (numbered from bus 2 to bus ), and bus (numbered from bus to bus ). The power flow equations defined in polar coordination can be expressed as [2] where and , are the active load powers of the buses; , are the real power output of the generators; , are the voltage phase angles of the buses; , are the voltage magnitude of the buses. Equation (7) characterizes the real power balance at buses and the real and reactive power at buses.
The Newton iteration for solving the power flow equation (7) proceeds as follows:
The Jacobian matrix in (8) is
To study the stability of solutions, a dynamic system is considered where . It is evident that the solutions of algebraic equation (7) correspond to the equilibrium points of the ordinary differential equation defined in (9). A solution of algebraic equation is exactly an equilibrium point of the ordinary differential equation. For a power flow solution, it is a stable equilibrium point of (9) if all of its eigenvalues lie in the open left half plane while a type1 equilibrium point of (9) has exactly one eigenvalue that lies in the open right half plane. Type2 equilibrium points of (9) are similarly defined.
We study the convergence region of Newton method on a modified 3bus system shown in Figure 1. The modified 3bus system contains 2 buses and 1 slack bus, with bus 1 being selected as the slack bus, and the remaining buses, bus 2 and bus 3, are both buses. Hence, there are only two real power equations to be solved with two unknowns (the voltage angles and ).
The power flow solutions obtained by using the Newton method are summarized in Table 1. The color column in Table 1 denotes the colors in the convergence regions and its fractal boundaries are highlighted in Figure 2, which is a grid of initial conditions in and , representing 160,000 power flow solutions. The scale in each voltage angle is 0 to in radians. The blue region is the convergence region of stable equilibrium point (SEP). The Newton power flow method starting from the initial conditions in this region all converges to the stable equilibrium point of (9). The green, purple, and skyblue regions correspond to those initial conditions that converge to the three type1 unstable equilibrium points (UEP), at which one of the eigenvalues of the Jacobian is positive and the other is negative. Finally, the red and yellow parts are the convergence regions of two type2 equilibrium points, at which both eigenvalues are negative. The numerical difference between pixels is 0.157 in radians.

Note that there are regions in the state space in which a small change in initial conditions results in a different convergence region. The number of iterations required for convergence is large in the fractal areas. Depending on the maximum iteration number used in the power flow program, these fractal areas may change. Note that due to the singularity of Jacobian, the Newton power flow will not converge from any initial condition at which the Jacobian is singular or from a point that maps on to a point where the Jacobian is singular.
2.2. Required Number of Iterations for Initials in Region of Convergence
The convergence of a Newton power flow has strong relationship with the maximum iteration number and the convergence threshold used in the computation. Newton power flow method usually converges within 10 iterations. For the study of required iterations starting from different initial points, we perform numerical studies on the modified 3bus system (see Figure 1) with the following parameters:(i)threshold for iterations: 20;(ii)threshold for convergence in power mismatch: .
Table 2 shows the required number of iterations starting from different initial conditions. The initial points that converge within 10 iterations constitute about 98.5%. Also, the required iterations from initials that lie in or near fractal boundaries are shown in Figure 3. The colormap in Figure 3 denotes the correspondence of colors with the required iterations. As can be seen from the figure, the initial points in convergence region can nicely converge with 10 iterations, which are always less than 10 and the initial points which require large iteration numbers all lie near the boundaries of convergence regions.

3. Convergence Regions under Different Conditions
In this section, the convergence regions on several conditions are investigated via a modified 3bus system and the IEEE 14bus system:(i)basecase,(ii)various loading conditions,(iii)different load models (models representing load demands).
For the modified 3bus system, there are 2 buses and 1 slack bus. Hence the state variable in power flow equations is a 2dimensional vector, which means that the convergence region of the modified 3bus system is 2dimensional. For the IEEE 14bus system, it contains 4 buses, 9 buses, and 1 slack bus. Hence the state variable is a 22dimensional vector, which means that the convergence region of IEEE 14bus system is 22dimensional. For the sake of the visualization of the convergence region, the same twodimensional crosssections are selected to display convergence regions. The basic data for study is listed as follows:(i)horizon axis is the voltage angle in bus 2, and the vertical axis is the voltage angle in bus 3;(ii)the scale in each coordinate is from 0 to ;(iii)each figure is grid of initial conditions on () plane, representing 160,000 power flow solutions.
3.1. Convergence Regions at BaseCase
Figure 4 shows the convergence region of IEEE 14bus system of the base case on the plane. The blue region is the convergence region of the SEP. The power flow starting from the initial conditions in this region all converges to the stable equilibrium point. The green region is the convergence region of one of the type1 equilibrium. The left colorized regions, like green, purple, yellow, and so forth, are regions corresponding to those initial conditions that converge to the corresponding equilibrium points. Finally, the white parts in grid are the divergence region, which occupies a great deal of area and is much larger than that in the modified 3bus system.
The number of power flow solutions can be large. We were able to find 20 power flow solutions for the basecase power system of the IEEE 14bus power system in the region which is a on () plane and find 6 power flow solutions for the basecase power system of the modified 3bus system. As the loading conditions increase, the number of power flow solutions decreases in pair, due to the saddle node bifurcations [20–22]. The total number of power flow solutions at each loading condition is summarized in Table 3.

3.2. Convergence Regions at Various Loading Conditions
We next study the convergence region of each power flow solution under different loading conditions. Since the number of power flow solutions decreases with the increase of loading conditions, we made efforts to trace each power flow solution. The loading effect on the convergence region is simulated via a loading factor . Although the convergence regions differ with different loading conditions, the fractals persist for different loading conditions, which are shown in Figure 5.
(a)
(b)
(c)
However, the number of power flow solutions decreases while the convergence region of SEP increases in size with the increase of loading factor , which is shown in Figure 5. In other words, the number of power flow solutions decreases when the loading conditions become heavier. During this process of load increases, saddle node bifurcations occur which brings about the solutions that vanished in pair, and then at the end there are only two power flow solutions left (one SEP and one type1 UEP); finally these two power flow solutions come together when the loading factor reaches its bifurcation value at which these two power flow solutions disappear due to the saddle node bifurcation (SNB).
Figure 6 shows the numbers of iterations required by the Newton power flow method to converge to a power flow solution under different loading conditions. The required number of iterations for those initial points in convergence region increases when the loading factor increases. It implies that it is more difficult to converge when the system is in heavy loading conditions. Moreover, the two solutions (corresponding to SEP and UEP of dynamical system) become closer as the loading conditions get heavier, and they coincide when the loading factor is 9 at which the saddle node bifurcation occurs.
(a)
(b)
(c)
(d)
The power flow solutions are illustrated in Table 4 to give a contrast. It is interesting to note that the coordinate of stable power flow solutions is quite different from that of the unstable power flow solutions.

3.3. Convergence Regions with Different Load Models
Several popular load models are used in the power flow analysis to better reflect static load behaviors. They are (1) constant power model (the power does not change with the movement of voltage magnitude); (2) constant impedance (the power is direct proportion of the square of voltage magnitude); (3) constant current (the power is direct proportion of voltage magnitude); ZIP model (the hybrid model of the above three models).
The polynomial equations for representing the ZIP load model can be expressed as where and . The parameters , , and represent the proportions of constant impedance, constant current, and constant power load, respectively. At the basecase condition, all loads are modeled as constant power. Table 5 lists ZIP load parameters, number of power flow solutions, and sizes of convergence regions for different load models.

The convergence regions of IEEE 14bus system with constant impedance, constant current, constant power, and ZIP load models are plotted in Figure 7. One interesting observation is that these convergence regions are similar in shape with fractal boundaries. Nevertheless, the number of power flow solutions is different among different load models. In accordance with the number of solutions with various load models, the descending order is as follows: constant impedance, constant current, ZIP and constant power, as is shown in Table 5.
(a) Constant impedance load
(b) Constant current load
(c) Constant power load
(d) ZIP load
4. Convergence Regions of Fast Decoupled Newton Method
Fast Decoupled Power Flow (FDPF) method is another powerful power flow technique which is based on the property that the variation of active power is mainly related to the variation of voltage phase angles in buses and the variation of reactive power with the variation of voltage magnitudes in buses. FDPF was proposed on the basis of the following numerical properties of practical power networks:(i)usually ; therefore ;(ii)usually is small, so .
The Jacobian matrix (see (9)) of Newton method is approximated by the following:
The fast decoupled method solves the following equations:
Compared with Newton method, the FDPF method brings about some merits: a smaller memory space and computational speed. However, FDPF method has the following issues:(1)convergence rate is linear;(2)it may diverge when one of the situations occurs:(a)the line parameters have a high ;(b)the system is heavily loaded.
Based on the above considerations, we numerically investigate the convergence regions of FDPF on the modified 3bus system and the IEEE 14bus system. We first consider the basecase, followed by various loading conditions and followed by different ratios.
4.1. BaseCase
Figure 8 shows the convergence region of IEEE 14bus system at the base case under the FDPF method on the () plane. The colors blue and white present the initial points which converge to the power flow solution and the initial points that diverge, respectively. Similar with the convergence region of the Newton method, the convergence region of the FDPF also has fractal boundary. While only one power flow solution is obtained starting from different initial points which is different from the Newton method under which multiple solutions are reached. This implies that the FDPF method is not appropriate for the computation of multiple solutions.
The number of iterations needed for the FDPF method to converge is shown in Figure 9. The colormap in Figure 9 denotes the required iteration numbers. As can be seen from the figure, the initial points which lie near the power flow solution will converge with small number of iterations, while the initial points which lie near the boundary of the convergence region require a large number of iterations or even diverge.
4.2. Various Loading Conditions
Normally, the voltage angle difference between two buses during a light loaded condition is small. However, the increase in loading condition causes a large voltage angle difference between two buses, which may lead to a violation of an assumption needed in the FDPF method. We next study the impact of loading conditions on the convergence regions of FDPF method on the modified 3bus system. The loading effect on the convergence region is studied with the help of a loading factor . Figure 10 shows the iteration numbers required by the FDPF method to converge to a power flow solution.
(a)
(b)
(c)
(d)
It is obvious that starting from an initial condition the number of iterations needed to converge to a solution increases with the increase of the loading factor. All initial points that converge to a power flow solution under the FDPF within a maximum number of 30 of the basecase () are shown in Figure 10. As expected, the convergence region shrinks significantly when the loading factor is 8, which is near the load margin of the system in which many initial points need the maximum number of iterations to converge or diverge. This fact can serve to explain why the FDPF method performs poorly during heavilyloaded conditions.
4.3. Different Ratio
One condition used to derive the FDPF method is that the resistance is much smaller than the reactance of the transmission lines of the power system. With the increment of ratio, the numbers of iterations required by the FDPF method will increase. When the ratio is high, the FDPF method tends to diverge. We next study the impact of ratio on the convergence region of the FDPF method.
Figure 11 shows the number of iterations needed for the FDPF method to converge to a solution with different ratios. In the base case, the ratio is ; the number of iterations required from different initial conditions for the FDPF method is shown in Figure 11(a). Figures 11(b)–11(d) show the number of iterations needed for the FDPF method to converge with 5, 10, and 13 times of the ratio of the base case. It is manifest that the ratio has great effect on the convergence of FDPF method. These numerical studies agree with the consensus that the FDPF method tends to diverge when the ratio is high.
(a)
(b)
(c)
(d)
5. Conclusion
This paper numerically studies the convergence region of power flow solutions under the Newton method and the fast decoupled power flow method which is a vital variation of the Newton method. The simulation results of the Newton method indicate that (1) the convergence regions of Newton method have fractal boundary, and the initial points which lie near the boundary of convergence regions require large iteration numbers; (2) these fractal boundary of convergence regions persists under different loading conditions (the lightlyloaded condition, the mediumloaded condition, and the heavilyloaded condition) and different load models (constant power, constant current, constant impedance, and ZIP); (3) the convergence region shrinks with the increase in loading conditions, and the required number of iterations for those initial points in convergence region increases when the loading factor increases. The simulation results of the FDPF method show that (1) the convergence region of the FDPF also has fractal boundary, but multiple solutions cannot be detected under the FDPF method; (2) the convergence region shrinks with the increase of the loading condition; (3) the numbers of iterations required by the FDPF method will increase with the increment of the ratio. These indicate that the FDPF method does not have a good convergence performance in heavy loading condition and high ratio condition.
These properties highlight the importance of choosing initial guesses for the Newton method and the FDPF method. This study of convergence region provides insights of the complexity of the NR method as well as the fast decoupled Newton method in finding power flow solutions even for largescale power systems. This study also sheds lights on the difficulty of finding a power flow solution of heavilyloaded power systems. These properties may give some reference to the application and enhancements of Newton power flow method and fast decoupled power flow method.
Acknowledgments
This work was supported by National Basic Research Program of China (2012CB215102) and State Grid Corporation of China, Major Projects on Planning and Operation Control of Large Scale Grid (SGCCMPLG0282012).
References
 J. B. Ward and H. W. Hale, “Digital computer solution of power flow problems,” AIEE Transactions, vol. 75, no. 3, pp. 398–404, 1956. View at: Google Scholar
 B. Stott, “Review of loadflow calculation methods,” Proceedings of the IEEE, vol. 62, no. 7, pp. 916–929, 1974. View at: Publisher Site  Google Scholar
 S. Iwamoto and Y. Tamura, “Load flow calculation method for illconditioned power systems,” IEEE Transactions on Power Apparatus and Systems, vol. 100, no. 4, pp. 1736–1743, 1981. View at: Google Scholar
 W. F. Tinney and C. E. Hart, “Power flow solution by Newton's method,” IEEE Transactions on Power Apparatus and Systems, vol. 86, no. 11, pp. 1449–1460, 1967. View at: Google Scholar
 B. Stott and O. Alsac, “Fast decoupled load flow,” IEEE Transactions on Power Apparatus and Systems, vol. 93, no. 3, pp. 859–869, 1974. View at: Google Scholar
 R. Idema, D. J. P. Lahaye, C. Vuik, and L. van der Sluis, “Scalable NewtonKrylov solver for very large power flow problems,” IEEE Transactions on Power Systems, vol. 27, no. 1, pp. 390–396, 2012. View at: Publisher Site  Google Scholar
 Y. Zhang and H. Chiang, “Fast newtonFGMRES solver for largescale power flow study,” IEEE Transactions on Power Systems, vol. 25, no. 2, pp. 769–776, 2010. View at: Publisher Site  Google Scholar
 T. Ochi, D. Yamashita, K. Koyanagi, and R. Yokoyama, “The development and the application of fast decoupled load flow method for distribution systems with high R/X ratios lines,” in IEEE PES Innovative Smart Grid Technologies (ISGT '13), pp. 1–6, February 2013. View at: Google Scholar
 P. R. Bijwe, B. Abhijith, and G. K. V. Raju, “Robust three phase fast decoupled power flow,” in IEEE/PES Power Systems Conference and Exposition (PSCE '09), pp. 1–5, March 2009. View at: Publisher Site  Google Scholar
 P. J. Lagace, M. H. Vuong, and I. Kamwa, “Improving power flow convergence by Newton Raphson with a levenbergmarquardt method,” in IEEE Power and Energy Society General Meeting (PES '08), pp. 1–6, July 2008. View at: Publisher Site  Google Scholar
 S. Abe, N. Hamada, A. Isono, and K. Okuda, “Load flow convergence in the vicinity of a voltage stability limit,” IEEE Transactions on Power Apparatus and Systems, vol. 97, no. 6, pp. 1983–1993, 1978. View at: Google Scholar
 H. Susanto and N. Karjanto, “Newton's method's basins of attraction revisited,” Applied Mathematics and Computation, vol. 215, no. 3, pp. 1084–1090, 2009. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 A. K. Bisoi, “Newton raphson method, scaling at fractal boundaries and mathematica,” Mathematical and Computer Modelling, vol. 21, no. 10, pp. 91–102, 1995. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 J. H. E. Cartwright, “Newton maps: fractals from Newton's method for the circle map,” Computers and Graphics, vol. 23, no. 4, pp. 607–612, 1999. View at: Publisher Site  Google Scholar
 N. Varghese, J. S. Thorp, and H. D. Chiang, “Fractally deformed basin boundaries of pendulum systems: new approaches in the study of swing dynamics,” in Proceedings of the 26th IEEE Conference on Decision and Control, vol. 26, pp. 885–888, December 1987. View at: Google Scholar
 M. Varghese and J. S. Thorp, “Analysis of truncated fractal growths in the stability boundaries of threenode swing equations,” IEEE Transactions on Circuits and Systems, vol. 35, no. 7, pp. 825–834, 1988. View at: Publisher Site  Google Scholar
 J. S. Thorp, S. A. Naqavi, and H.D. Chiang, “More load flow fractals,” in Proceedings of the 29th IEEE Conference on Decision and Control, vol. 6, pp. 3028–3030, December 1990. View at: Google Scholar
 J. S. Thorp and S. A. Naqavi, “Loadflow fractals draw clues to erratic behavior,” IEEE Computer Applications in Power, vol. 10, no. 1, pp. 59–62, 1997. View at: Google Scholar
 C. T. Kelly, Solving Nonlinear Equations with Newton's Method, vol. 1 of Fundamentals of Algorithms, SIAM, 2003.
 I. Dobson and H. Chiang, “Towards a theory of voltage collapse in electric power systems,” Systems and Control Letters, vol. 13, no. 3, pp. 253–262, 1989. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 H. Chiang, S. Rovnyak, J. S. Thorp, and R. J. Thomas, “Toward a general power system voltage collapse model,” in Proceedings of the 28th IEEE Conference on Decision and Control, vol. 1, pp. 338–339, December 1989. View at: Google Scholar
 H. Chiang, I. Dobson, R. J. Thomas, J. S. Thorp, and L. FekihAhmed, “On voltage collapse in electric power systems,” IEEE Transactions on Power Systems, vol. 5, no. 2, pp. 601–611, 1990. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2013 JiaoJiao Deng and HsiaoDong Chiang. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.