A Common Framework for Developing Robust Power-Flow Methods with High Convergence Rate

Featured Application: This work proposes a common framework for developing robust Power-Flow solvers that may ﬁnd multiple applications in power system analysis tools. Abstract: This paper presents a novel Power-Flow solution paradigm based on the structure of the members of the Runge–Kutta family. Solution approaches based on the introduced solution paradigm are intrinsically robust and can achieve high-order convergences rates. It is demonstrated that some well-known Power-Flow solution methods are in fact special cases of the developed framework. Explicit and embedded formulations are discussed, and two novel solution methodologies based on the Explicit Heun and Embedded Heun–Euler’s methods are developed. The introduced solution techniques are validated in the EU PEGASE systems, considering different starting points and loading levels. Results show that the developed methods are quite reliable and efﬁcient, outperforming other robust and standard methodologies. On the basis of the results obtained, we can afﬁrm that the introduced solution paradigm constitutes a promising framework for developing novel Power-Flow solution techniques.


Motivation
The Power-Flow (PF) is one of the main computational tools in power system analysis. It solves the nonlinear model of a power system for a given power generation/consumption profile. Nowadays, the Newton-Raphson (NR) is still the most widely used technique for solving the nonlinear equations involved in the PF problem. This is because of its simplicity, easy codification in standard programming languages, acceptable robustness, and quadratic convergence. These features typically bring superior characteristics with respect to (w.r.t.) other conventional PF solvers.
The PF problems are customary categorized as well or ill-conditioned. Conventionally, a PF problem is said to be ill-conditioned when its solution by using standard techniques such as NR is hardly reachable [1]. This has supposed an important barrier for the wide application of Newton-like techniques; nevertheless, this issue has not been very problematic due to most PF cases being well-conditioned. However, ill-conditioned cases are becoming more frequent [2]. Solving these kinds of problems is still challenging due to the difficulties showed by standard methods or the low degree of efficiency manifested by most of the available robust PF solution techniques. In this context, developing novel efficient and robust PF solution approaches is still interesting, relevant, and demanding.

•
Methodologies developed in this work are intrinsically robust, which means, unlike those in [1,[7][8][9][10][11][12], it is not necessary to include a discrete step size in the formulation to improve their convergence properties. • Methodologies based on the introduced solution paradigm achieve, at least, a quadratic convergence rate, which makes them competitive with standard solvers such as NR.
The latter point is especially relevant since PF solution techniques based on the proposed solution framework might be also high order Newton-like methods, which is reflected in the superior convergence characteristics in comparison with NR and other robust PF methods. In this paper, Explicit and Embedded formulations of the developed solution paradigm are discussed, and two novel robust and efficient PF solution techniques with cubic order of convergence are developed. In order to validate the developed methodologies, they have been tested in several realistic systems based on the European Transmission system.
The remainder of this paper is organized as follows. Section 2 introduces the developed common framework for developing high-order robust PF solution methods. Two novel solution techniques based on the introduced solution paradigm are developed in Section 3. Several numerical experiments and results are presented in Section 4. Finally, the paper is concluded in Section 5.

Analogy between NRJ and the Explicit Midpoint Method
Let us consider the PF problem in polar coordinates as a set of nonlinear equations as follows: where g : R n −→ R n are the PF equations and x ∈ R n is the PF state vector. Due to strong nonlinearities in (1), its solution must be achieved by using iterative techniques. The most popular one in PF analysis is NR, which is performed as follows: where [ * ] −1 : R n → R n is the inverse operator, the superscript denotes the k th iteration of the iterative procedure, and g ∈ R n×n is the Jacobian matrix of the PF equations, which is formed by the first partial derivatives of g w.r.t. x. In contrast to (2), NRJ poses a two-step iterative scheme for solving the set in Equation (1) described by: In [10], the Explicit Midpoint rule [19] was exploited for solving the PF equations on the basis of the Continuous Newton's principle. This technique poses a two-stage iterative solver, which is given by the mapping in (4).
where ∆t ∈ R + is the step size (or truncation) and f : R n −→ R n is given by: Further insight on mappings (3) and (4) reveals a clear analogy between these two solvers. In fact, the iterative map (3) could be rewritten as follows: where the step size has vanished and the auxiliary function h : R n −→ R n is defined by: The analogy described above reveals that NRJ and the Explicit Midpoint rule are actually equivalent mappings. Keeping this in mind, NRJ may be seen as the application of the Explicit Midpoint method for solving the function h, whereas the different methodologies proposed in [1,10] are devoted to solving f. In the same way, NR could be assimilated as the direct application of the Explicit Euler methodology for solving (7).

The Novel Family of Robust PF Techniques
In the previous section, we showed that NR and NRJ can be seen as the application of the Explicit Euler and Midpoint techniques for solving (7). This analogy suggests that any other member of the family of Explicit Runge-Kutta methods can be adapted to be applied to the PF problem. Similar to the generic form of the Runge-Kutta methods (e.g., see [10] (Equation (11)), we could develop a generic form for the application of this family of numerical methods to the PF problem, as follows: where the a and b values ∈ R. In essence, these parameters are introduced to copy the standard structure of the Runge-Kutta methods, which are broadly defined by the total number of stages involved (p), the value of the elements of the Runge-Kutta matrix (a values) and weights (b values) [19]. Conventionally, the Runge-Kutta methods are presented by their so-called Butcher tableaus [19] (p. 95), which succinctly collects all the necessary information to completely define a Runge-Kutta technique, as follows: It is worth mentioning that (9) is lower triangular, since we are dealing with explicit Runge-Kutta formulas (the implicit formulation of (8), although possible, may be computationally expensive, since it would require the factorization of the inverse of (7) [1]). In addition, it is worth noting that nodes [19] are not included in the formulation developed in this paper. This is because the application of the Runge-Kutta formulas to PF analysis is developed from the assimilation of the set of equations (1) as an autonomous set of differential equations. In such formulations, the nodes of a numerical method lack importance [19].
One can note that (8) can be adapted to any Explicit Runge-Kutta formula. At this point, readers can easily check that both NR and NRJ are actually special cases of (8). Thus, for NR, one has: While NRJ is fully described by: Interested readers can easily check in the vast literature that (11) and (12) correspond with the Runge-Kutta matrix and weight vectors of the Explicit Euler and Midpoint techniques, respectively.

Embedded Formulation
The solution paradigm established in the previous section allows considering any explicit numerical method for developing robust PF solution techniques. However, not only explicit numerical integration techniques are available in the literature. The Embedded Runge-Kutta techniques arise as an alternative of Explicit methods. Such methodologies allow accurately estimating the local truncation error of a single step and properly adapting the step size to obtain high accuracy [19]. In this case, due to the step size having vanished in the developing formulation, its ability to estimate the local truncation error lacks importance; nonetheless, it is still relevant to rewrite (8) into its Embedded counterpart as follows: where b * values ∈ R. Unlike paradigm (8), each weight in (13) is described by a pair of parameters (b and b * ). In this case, two numerical methods are nested, each contributing their own characteristic weights. In numerical analysis, this arrangement is very useful for achieving higher accuracy, while in our case, these parameters have been introduced just for copying the standard structure of the Embedded Runge-Kutta formulas [19] (p. 204). Thus, instead of directly calculating the value of the state vector for the k + 1 iteration through a single step arrangement, formulation (13) firstly calculates an approximation of the state vectorx that is posteriorly refined by means of other numerical methods defined by the b * values (see [20] for further details). This way, any PF solver based on the developed paradigm (13) could be defined by two matrixes as follows: It is worth noting that unlike the techniques developed in [1,[7][8][9][10][11][12], any truncation of the Newton's increment vector based on the influence of the step size is considered in the developed solution frameworks. It means that the developed solution techniques are (as shown in Section 4) intrinsically robust. This point is relevant, since the influence of the step size vanishes, facilitating its codification and application in industry tools.

Handling Equipment Limits
Handling equipment limits (such as those imposed by e.g., generators, on-load tap changers, etc.) is an important topic in PF analysis. The developed solution frameworks (8) and (13) do not directly take into account these limits; however, they could be easily incorporated. In this case, we have considered the same strategy used in [21]. Particularizing to the case of reactive limits in generation units, this approach consists of checking the reactive limits of all generators after a PF solution is calculated. Then, if any limit has been violated, those buses are converted to PQ, and the violated limit is taken as the reactive power injection of this new PQ bus. Then, taking the obtained solution as a starting point, the process is repeated until a feasible solution is reached. This strategy is described in the flowchart of Figure 1. One should note that other control parameters such as transformer tap changers can be taken into account using similar strategies. It is worth mentioning that the developed solution framework is versatile enough to easily incorporate other strategies, such as that described in [2].

Stability of the Developed Solution Paradigms
In order to study the stability of the mappings (8) and (13), let us introduce the following definition: Contrarily, a fixed point is said to be unstable if it is not stable.
By Definition 1, one requires a PF solution (i.e., x * such that g(x * ) = 0) to be asymptotically stable. In addition, we desire from x * to be a sink. By using definition 1, we can study the numerical stability of the developed solution frameworks, which is performed in the following lemma.
Lemma 1: a solution of the PF equations (i.e., x * such that g(x * ) = 0) is a sink for the nonlinear solution paradigms (8) and (13) Proof of Lemma 1. By differentiating (7) with respect x one obtains: One can easily check that at an equilibrium point (i.e., x * such that g(x * ) = 0), the developed solution paradigms (8) and (13) take the form: Substituting (17) into (16), one obtains: where I ∈ R n×n is the identity matrix. Then, one can replace (18) into the differential equation associated with the mapping (8) obtaining: While for the paradigm (13), one has: Equations (19) and (20) show that at the equilibrium point, all the eigenvalues of the Jacobian of the differential equations associated with the developed solution paradigms do not have imaginary parts. In addition, their real parts are equal to − ∑ m i = 1 {b i } in the case of (8) and − ∑ m for the mapping (13). Therefore, the equilibrium point x * is asymptotically stable for the paradigm (8) if: While for the mapping (13), the following inequality must hold:

Developed PF Solvers
For the sake of exemplifying, we have developed two novel robust techniques based on the introduced solution paradigms defined by (8) and (13).
Firstly, we have taken the Explicit Heun's method for solving (7), which brings the following iterative procedure for solving the PF equations: On the other hand, taking the developed formulation (13) and considering the Embedded Heun-Euler method, the following iterative procedure for solving the PF equations is introduced: The developed PF solution techniques require two matrix evaluations (and factorizations) along with a function evaluation each iteration. In addition, they achieve cubic convergence (see Appendix A). As a result of these features, one can guess that the developed techniques are more efficient than those presented in [1,9,10]. It is worth remarking on the outstanding features of the developed methods since, besides their intrinsic robustness characteristics that are empirically proved by results (see Section 4), they present a high convergence rate. These two characteristics are very difficult to be found simultaneously in a PF solver.
For the sake of summarizing, Table 1 shows the A matrix and b vectors (or B matrix), of NR, NRJ, and the developed PF solution methods (23) and (24). To complete the section, the PF solution procedure using any of the developed techniques is summarized in Algorithm 1 using pseudocode. This algorithm aims to serve as a guidance for easily implementing the developed solvers in standard programming languages such as MATLAB or C. It is worth mentioning that this algorithm can be adapted to NR and NRJ by just replacing Equation (23) or (24) by (2) or (6), respectively. In addition, we have considered the failure of our algorithms when they do not converge in a predetermined number of iterations k max . This situation usually indicates divergence. In the pseudocode of Algorithm (1), the parameter tolerance makes mention of a predefined value beyond which it is assumed that the calculated solution is accurate enough and the algorithm may be determined successfully convergent. It is worth noting that the iterative nonlinear algorithms are unable to calculate the exact solution as in the case of linear systems; therefore, it is customary to establish a threshold for determining when the correct solution is sufficiently approached. Typical values for convergence tolerances in PF studies are fixed below 10 −4 [7][8][9][10][11].

Numerical Experiments
A variety of numerical experiments have been carried out in order to compare the performance of the developed techniques with other well-known PF solvers. (23) and (24), respectively) are compared with the standard NR, the Reverse Bulirsch-Stoer solver (RBS) [9], the Newton-Raphson with Jacobian adjustments (NRJ) [6], and the implicit implementation of the Continuous Newton's method using the Backward-Euler technique (BEM) [12]. In this latter case, we have used the implementation labelled 'Implicit-Continuous Newton Method with approximated Hessian and the solution of the inner loop is approximated by its first iteration (ICNM-J1)' in reference [12].

The developed Explicit Heun and Embedded Heun-Euler solution approaches (Equations
All simulations have been obtained using Matpower [22] and running on a 3 GHz Intel Core i5-8500 personal laptop (8.00 GB RAM). In all experiments, a convergence tolerance of 10 −5 and k max = 2000 have been considered. The maximum number of iterations was tuned high enough to only discriminate the divergence cases. In order to avoid the impact of parallel computational activities, the reported execution times have been calculated taken as the average values of 1000 simulations. To determine the goodness of the solutions calculated, they have been compared with those obtained using NR and the default starting point provided in Matpower. Thereby, if both solutions differ more than 0.1 per unit (p.u.), the obtained solution is considered inaccurate.

Studied Systems
Five realistic systems based on the European Transmission system from the EU PE-GASE project have been considered [23,24]. These meshed networks range from 89 to 13,569 buses, thus comprising a wide variety of scenarios. For the sake of self-sufficiency, Table 2 reports the main characteristics of the studied systems. A starting point, built by modifying the correct solution of the studied cases using a Gaussian distribution and a determined standard deviation, has been considered. These starting guesses aim at emulating a real situation in which an available system solution is used for initializing the PF analysis. In our opinion, this approach is more realistic than the conventional flat start, because the iterative algorithm starts from an explicit solution of the system; this way, the starting point could be considered more "natural" than a flat guess. We call the reader's attention toward a similar approach that has been also considered in [12]. In order to enable the reproducibility of the results reported, the studied systems with the considered starting points can be found in [25].

Convergence Rates for Base Cases
Our first test is devoted to comparing the convergence rates for base cases. Table 3 provides the total iterations for the studied systems. As observed, NR failed in all studied cases, which denotes a high degree of ill-conditioning [1]. RBS only successfully solved the 89-and 1354-bus cases, while NRJ failed in the 9241-and 13659-bus cases. On the other hand, BEM converged to the low voltage solution in the 13659-bus system. The developed PF solvers successfully solved all the studied cases. In addition, our proposals usually required less iterations than the remainder of the methodologies. These results serve to empirically demonstrate the robustness characteristic and efficiency salient features of the developed PF solvers.

Convergence Rates with Reactive Limits
Next, we explore the ability of the tested methods to solve the studied systems with reactive limits. To do that, the strategy described in Section 2.4. has been implemented within the considered algorithms. Table 4 reports the total iterations for the studied systems. Similar conclusions as for the base case can be extracted for this scenario.

Solution Times
In this case, we are interested in comparing the computational efficiency of the tested PF solvers. This is done by comparing their solution times. Tables 5 and 6 provide the  solution times for the cases reported in Tables 2 and 3, respectively. As observed, those solvers based on the developed solution paradigm were the most competitive techniques, outperforming RBS and BEM in all studied systems. Competitive results obtained with NRJ and the solution techniques based on Equations (23) and (24) are undoubtedly due to them typically requiring less factorizations than other solvers. In this regard, it is worth commenting that the factorization of the Jacobian matrix is, by far, the heaviest calculation in PF analysis [1]. This fact is studied in Tables 7 and 8, where the total factorizations required for the cases analysed in Tables 2  and 3, respectively, are provided.

Influence of the Loading Level
It is well known that the loading level may affect the overall performance of a PF solver. Typically, PF solvers require more iterations as the loading level grows. In some cases, this factor may even provoke instability. Figure 2 shows the total iterations of the studied solvers for different loading levels without reactive limits. The loading levels for the studied cases have been modified by multiplying the injected active and reactive power at PQ buses along with the injected reactive power at PV buses by a real positive factor (called loading level). As observed, our solvers and NRJ normally required less iterations than BEM and RBS. Nevertheless, it is worth noting that the developed techniques required fewer iterations than NRJ for high loading levels. To provide a clearer interpretation of the Figure 2, Table 9 reports the total iterations required by the different solvers only at the maximum loadability point.

Conclusions and Future Works
Inspired by the similitude observed between the methodology described in [6] and the Explicit Midpoint method, we have developed a novel framework for developing efficient and robust PF solution approaches based on Runge-Kutta formulas. This solution paradigm allows adapting any member of the Runge-Kutta family to the PF problem. Although this idea is similar to that described by Milano in [1], the developed framework presents two important advantages:

•
The influence of the step size vanishes from the formulation, which means that the developed solution paradigm is intrinsically robust.

•
High convergence rates can be achieved, even higher than two.
Explicit and Embedded formulations of the developed solution paradigm have been discussed, and two PF solution techniques based on the Explicit Heun and Embedded Heun-Euler's methods with cubic order of convergence have been developed.
The developed techniques have been validated using the EU PEGASE systems considering different demanding scenarios. Results have been compared with those obtained with other well-known robust and conventional PF solutions approaches. Our proposals turned out to be more stable and efficient than the other tested techniques, which occasionally showed an unstable behaviour or a high computational cost.
On the basis of the results obtained, it is concluded that the introduced solution paradigm constitutes a promising framework for developing novel PF solution techniques. This opens the door for further applying the developed solution paradigms. Future works should be focused on developing and validating other PF solution techniques, applying other members of the Runge-Kutta family to the developed solution framework. Similarly, analogies between the introduced solution paradigm and other families of numerical integration techniques such as Adams-Bashforth's methods [8] or Bulirsch-Stoer's scheme [9] can be found and, consequently, other high-order robust PF solution approaches can be developed. We expect the emergence of this kind of techniques.
The developed solution techniques may be also adapted to other Power System tools such as the Continuation Power-Flow [26]. Here, their robustness along its high convergence order may bring very good results. Future works should study the application of the introduced PF solution techniques to this problem.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A Proof of Cubic Convergence of the Developed Solvers
In this section, we prove that the developed PF solution techniques present a cubic convergence rate. To this end, we use the Taylor series technique (e.g., see [27]).
Theorem A1. Let g be sufficiently differentiable at each point of an open neighborhood D of r ∈ R n ; this is a solution of the system g(x) = 0. Let us suppose that g(x) is continuous and nonsingular in x. Then, the developed method (23) converges to r with order three.
Proof of Theorem A1. Taylor expansion of g(x) and g (x) about r yields: g x (k) = g (r) e (k) + C 2 e (k) 2 + C 3 e (k) 3 + O e (k) 4 (A1) g x (k) = g (r) I + 2C 2 e (k) + 3C 3 e (k) 2 where e (k) = x (k) − r ∈ R n , e i = (e, e, . . . , e) i−times , C j = (1/j!) g (r) −1 g (j) (r) ∈ L i (R n , R n ), g (j) ∈ L(R n × . . . × R n , R n ) and g (r) −1 ∈ L(R n ). Now, let us assume that: where c s ∈ R. Considering the following inverse definition: and solving the resulting linear system, one obtains: Finally, we calculate the error function at k + 1 as follows: (A9) As observed in (A9), the error function decreases cubically for each iteration, and therefore, the proof is completed. Theorem A2. Let g be sufficiently differentiable at each point of an open neighborhood D of r ∈ R n ; this is a solution of the system g(x) = 0. Let us suppose that g(x) is continuous and nonsingular in x. Then, the developed method (24) converges to r with order three.
Proof of Theorem A2. Taking the Taylor expansions of g x (k) and g x (k) −1 about r from (A1) and (A2), respectively, let us define e (k) x =x (k) − r. This way, we can write: At this point, we can easily calculate: