Journal of Applied Mathematics

Volume 2013 (2013), Article ID 509496, 12 pages

http://dx.doi.org/10.1155/2013/509496

Research Article

Convergence Region of Newton Iterative Power Flow Method: Numerical Studies

Jiao-Jiao Deng¹ and Hsiao-Dong Chiang1,2

¹School of Electrical Engineering and Automation, Tianjin University, Tianjin 300072, China

²School of ECE, Cornell University, Ithaca, NY 14850, USA

Received 22 March 2013; Revised 10 October 2013; Accepted 14 October 2013

Academic Editor: Luigi Muglia

Copyright © 2013 Jiao-Jiao Deng and Hsiao-Dong 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.

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 Newton-Raphson (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 14-bus test system which is of 22-dimensional 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 steady-state 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 three-phase steady-state conditions and the system is in sinusoidal steady-state;(ii)the transmission network is composed of constant, linear, and lumped-parameter 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 3-bus 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 Newton-Raphson (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 14-bus test system, as compared with the 3-bus 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 ill-condition.(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 large-scale power systems. This study also sheds lights on the difficulty of finding a power flow solution of heavily-loaded 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 type-1 equilibrium point of (9) has exactly one eigenvalue that lies in the open right half plane. Type-2 equilibrium points of (9) are similarly defined.

We study the convergence region of Newton method on a modified 3-bus system shown in Figure 1. The modified 3-bus 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 ).

Figure 1: A modified 3-bus system.

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 sky-blue regions correspond to those initial conditions that converge to the three type-1 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 type-2 equilibrium points, at which both eigenvalues are negative. The numerical difference between pixels is 0.157 in radians.

Table 1: Multiple power flow solutions.

Equilibrium
Color

Initial numbers
Stability

1
Blue
0.1191
0.1195
33703
Stable

2
Green
−0.3151
3.0906
30430
Type-1

3
Red
−1.8206
2.3522
16463
Type-2

4
Purple
3.0919
−0.3137
30432
Type-1

5
Yellow
2.3525
−1.8220
16364
Type-2

6
Sky-blue
−3.0952
−3.0940
32587
Type-1

Figure 2: Convergence region of each power flow solution of a modified 3-bus test system.

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 3-bus 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 color-map 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.

Table 2: Required Iterations starting from different initials.

Convergence iteration number
Initial numbers

[, ]
235

(2, 4]
60197

(4, 6]
69038

(6, 8]
21945

(8, 10]
6162

(10, 12]
1776

(12, 14]
472

(14, 16]
140

(16, 18]
25

(18, 20]
9

>20
1

—
—

Figure 3: The convergence region and the number of iterations needed for the Newton method on () plane.

3. Convergence Regions under Different Conditions

In this section, the convergence regions on several conditions are investigated via a modified 3-bus system and the IEEE 14-bus system:

(i)base-case,(ii)various loading conditions,(iii)different load models (models representing load demands).

For the modified 3-bus system, there are 2 buses and 1 slack bus. Hence the state variable in power flow equations is a 2-dimensional vector, which means that the convergence region of the modified 3-bus system is 2-dimensional. For the IEEE 14-bus system, it contains 4 buses, 9 buses, and 1 slack bus. Hence the state variable is a 22-dimensional vector, which means that the convergence region of IEEE 14-bus system is 22-dimensional. For the sake of the visualization of the convergence region, the same two-dimensional cross-sections 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 Base-Case

Figure 4 shows the convergence region of IEEE 14-bus 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 type-1 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 3-bus system.

Figure 4: Convergence region of power flow solutions of the IEEE 14-bus system on () plane.

The number of power flow solutions can be large. We were able to find 20 power flow solutions for the base-case power system of the IEEE 14-bus power system in the region which is a on () plane and find 6 power flow solutions for the base-case power system of the modified 3-bus 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.

Table 3: Number of power flow solutions and fractal properties at different loading conditions.

Case name
Loading condition
Number of solutions
Number of initials converged to SEP
Fractal

3-bus
Loading factor
1
6
33703
Yes

4.02
4
34743
Yes

5
2
90051
Yes

9.1
0
0
—

14-bus
Loading factor
1
20
16274
Yes

2
8
16538
Yes

4
2
17125
Yes

5
0
0
—

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)

Figure 5: Convergence regions of modified 3-bus system on different loading conditions: (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 type-1 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)

Figure 6: Iteration numbers of Newton method for the modified 3-bus system with different load conditions on () plane: (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.

Table 4: Power flow solutions on different loading conditions.

Loading factor
Power flow solutions

SEP
UEP

theta2
theta3
theta2
theta3

6
0.7663
0.7689
2.5373
2.5433

7
0.9318
0.9350
2.3713
2.3777

8
1.1377
1.1416
2.1650
2.1714

9
1.5089
1.5140
1.7935
1.7994

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 base-case 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.

Table 5: Load parameters, the size of convergence region, number of power flow solutions, and fractal property for different load models.

Load model type
Load parameters for all buses
Number of solutions
Number of initials converged to SEP
Fractal

a
b
c

Constant impedance
1
0
0
70
16495
Yes

Constant current
0
1
0
56
16493
Yes

Constant power
0
0
1
20
16274
Yes

ZIP
0.33
0.33
0.33
34
16486
Yes

The convergence regions of IEEE 14-bus 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

Figure 7: Convergence region of IEEE 14-bus system with different load models on () plane: (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 3-bus system and the IEEE 14-bus system. We first consider the base-case, followed by various loading conditions and followed by different ratios.

4.1. Base-Case

Figure 8 shows the convergence region of IEEE 14-bus 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.

Figure 8: The convergence region of FDPF for IEEE 14-bus system on () plane.

The number of iterations needed for the FDPF method to converge is shown in Figure 9. The color-map 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.

Figure 9: The convergence region and the number of iterations needed for the FDPF method for IEEE 14-bus system on () plane.

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 3-bus 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)

Figure 10: The convergence region and the number of iterations needed for the FDPF method to converge at different loading conditions on () plane: (a) the lightly loaded condition, ; (b) ; (c) ; (d) the heavily loaded condition, . From an initial condition, the number of iterations required by the FDPF method increases with the increase of loading condition.

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 base-case () 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 heavily-loaded 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)

Figure 11: The convergence regions and the number of iterations required for the FDPF method to converge with different ratios on the () plane: (a) ; (b) ; (c) ; (d) . The convergence region decreases in size with the increase in the ratios; rendering the FDPF method does not work well in the distribution power networks which have high ratios.

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 lightly-loaded condition, the medium-loaded condition, and the heavily-loaded 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 large-scale power systems. This study also sheds lights on the difficulty of finding a power flow solution of heavily-loaded 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 (SGCC-MPLG028-2012).

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.

B. Stott, “Review of load-flow calculation methods,” Proceedings of the IEEE, vol. 62, no. 7, pp. 916–929, 1974.

S. Iwamoto and Y. Tamura, “Load flow calculation method for ill-conditioned power systems,” IEEE Transactions on Power Apparatus and Systems, vol. 100, no. 4, pp. 1736–1743, 1981.

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.

B. Stott and O. Alsac, “Fast decoupled load flow,” IEEE Transactions on Power Apparatus and Systems, vol. 93, no. 3, pp. 859–869, 1974.

R. Idema, D. J. P. Lahaye, C. Vuik, and L. van der Sluis, “Scalable Newton-Krylov solver for very large power flow problems,” IEEE Transactions on Power Systems, vol. 27, no. 1, pp. 390–396, 2012.

Y. Zhang and H. Chiang, “Fast newton-FGMRES solver for large-scale power flow study,” IEEE Transactions on Power Systems, vol. 25, no. 2, pp. 769–776, 2010.

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.

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.

P. J. Lagace, M. H. Vuong, and I. Kamwa, “Improving power flow convergence by Newton Raphson with a levenberg-marquardt method,” in IEEE Power and Energy Society General Meeting (PES '08), pp. 1–6, July 2008.

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.

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.

A. K. Bisoi, “Newton raphson method, scaling at fractal boundaries and mathematica,” Mathematical and Computer Modelling, vol. 21, no. 10, pp. 91–102, 1995.

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.

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.

M. Varghese and J. S. Thorp, “Analysis of truncated fractal growths in the stability boundaries of three-node swing equations,” IEEE Transactions on Circuits and Systems, vol. 35, no. 7, pp. 825–834, 1988.

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.

J. S. Thorp and S. A. Naqavi, “Load-flow fractals draw clues to erratic behavior,” IEEE Computer Applications in Power, vol. 10, no. 1, pp. 59–62, 1997.

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.

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.

H. Chiang, I. Dobson, R. J. Thomas, J. S. Thorp, and L. Fekih-Ahmed, “On voltage collapse in electric power systems,” IEEE Transactions on Power Systems, vol. 5, no. 2, pp. 601–611, 1990.