On the Applicability of Two Families of Cubic Techniques for Power Flow Analysis

: This work presents a comprehensive analysis of two cubic techniques for Power Flow (PF) studies. In this regard, the families of Weerakoon-like and Darvishi-like techniques are considered. Several theoretical ﬁndings are presented and posteriorly conﬁrmed by multiple numerical results. Based on the obtained results, the Weerakoon’s technique is considered more reliable than the Newton-Raphson and Darvishi’s methods. As counterpart, it presents a high computational burden. Regarding this point, the Darvishi’s technique has turned out to be quite efﬁcient and fully competitive with the Newton’s scheme.


Motivation
In power system analysis, PF is probably the most important computational tool, playing a vital role in a wide variety of applications such as power system planning and operation, economic dispatch or security analysis, among others. From a mathematical point of view, PF consists on solving the set of equations that model the nonlinear relationships between nodal voltages and power injections. Since the PF equations are nonlinear, this problem cannot be directly solved. In this sense, iterative solvers have been customary used for solving the PF equations. From a PF solver one expects two main characteristics: • Robustness: for effectively solving ill-conditioned systems, which are mainly encountered in presence of cascading failures, heavy loading conditions or badly initialization of the iterative procedure. • Efficiency: this aspect is crucial for online applications, since it determines the ability of a PF solver for quickly finding the solution of the PF equations.
Many efforts have been devoted during decades on developing different PF solvers (see Literature Review). Developed approaches have been mainly focused on separately enhancing one of the characteristics enumerated above. However, in a modern power system paradigm, it is necessary to use tools that are simultaneously efficient and robust enough to, for example, properly handle large-scale ill-conditioned cases. In this context, it is very difficult to find a PF solver which gathers both robust and efficient features. For instance, the conventional Newton-Raphson (NR) can be considered an efficient technique, however, it fails in ill-conditioned systems being so its universal applicability questionable. Recently, some studies have been conducted on studying the applicability of various highorder Newton-like (HONL) methods for PF analysis [1][2][3][4], manifesting good performance in well-conditioned systems. However, these works lack of a formal analysis of the mathematical properties of such solvers, more specifically, their stability and robustness properties were not properly studied. This work aims at being a first step on filling this gap.
Regularization-based methods have been recently applied for solving PF equations, especially in ill-conditioned cases [23,24]. This approach basically consists of avoiding the singularity of the Jacobian matrix when it is ill-conditioned, by adding the so-called regularizing operator. Due to the Jacobian matrix is modified, the solution of the iterative process may not be the actual solution of the PF problem, as commented in [24] and reported in [25]. This issue has limited the applicability of this kind of techniques. Homotopic and holomorphic principles have been also applied for solving PF problem [26,27]. This kind of methodologies solve the PF equations by posing a continuous version of them by which their solution is reaching by progressively increasing a well-known homotopic parameter. This kind of approaches are quite robust and have found certain interest on ill-conditioned cases. However, as in the case of regularization-based methods, they may reach a nonphysical solution which lacks of interest [28]. Recently, the authors have developed various self-developed robust solvers [28,29], which arise from combination of different approaches. Although this kind of techniques present good computational performance, they are not competitive with other conventional techniques in well-conditioned cases yet.
Recently, HONL methods hav gained popularity for solving PF equations. Basically, a HONL is a nonlinear solved based on NR with a convergence order higher than two. Thus, a HONL should converge employing less iterations than NR. If besides this feature the solver is computationally efficient, the obtained PF solution approach may be competitive with NR. In this context, Pourbagher and Derakhshandeh firstly compared various HONL in [1], showed very good results. The conclusions in this reference motivated the authors to further study this kind of methods in [2][3][4], confirming that HONL methods may suppose an attractive alternative to NR.

Contributions and Paper Organization
Motivated by the good results showed by HONL techniques in recent studies. This work aims at further analyzing this kind of techniques for PF analysis. So far, HONL methods have been only employed for solving well-conditioned systems, therefore, their robust properties have not been analyzed from a theoretical and empirical point of view. This work aims at supposing a first step to filling this gap. In this sense, we focus on two well-known family of cubic techniques. In this context, we comprehensively analyze the cubic methods proposed by Weerakoon [30] (3OW) and Darvishi [31] (3OD). The stability, efficient and robustness properties of these solvers are firstly theoretically deduced and analyzed and, posteriorly, these theoretical foundations are empirically confirmed by numerical experiments. The two cubic methods are compared with the most standard PF solver (i.e., NR) and the 4th order Runge-Kutta method (RK4) developed in [10] with the aim of providing a benchmark study for the applicability of these families in both well and ill-conditioned systems. It is worth mentioning that these two cubic methodologies suppose the background for the development of many HONL techniques (e.g., see [32]). Thereby, this work aims at supposing a valuable guidance for applying other HONL techniques that are in fact modifications of the studied cubic techniques.
Remainder of this paper is organized as follows. Section 2 outlines the necessary background. Stability of two families of cubic techniques is studied in Section 3. Section 4 is devoted on comparing two families of cubic methods with the Newton-Raphson approach using a well-known efficiency index. Section 5 presents several numerical experiments with results. The main conclusions of this work are duly drawn in Section 6.

PF Solution Using NR
In polar coordinates, the PF equations describe the nonlinear relationship between the nodal voltages and power injections. Thus, generically considering the i th system bus, the PF equations are given by [1][2][3][4]: g P , ∀i ∈ N l , N g g Q , ∀i ∈ N l (1) where, N l and N g are the sets of PQ and PV buses, respectively. P sch i ∈ R and Q sch i ∈ R are the active and reactive power injected at i th bus, respectively. V i ∠δ i ∈ C is the complex voltage at i th bus. Y ij ∠θ ij ∈ C is the ij th element of the admittance matrix, δ g ∈ R n g is the vector of voltage angles at PV buses, δ l ∈ R n l is the vector of voltage angles at PQ buses and V l ∈ R n l is the vector of voltage angles at PQ buses, n l ∈ N and n g ∈ N are the total number of PQ buses and PV buses, respectively.
Because (1) are strongly nonlinear, they cannot be explicitly inverted. In PF analysis, iterative techniques are customary employed for addressing this issue and solving (1). Thereby, NR is often considered as the standard approach for PF purposes. PF solution using NR is established as a discrete map given by [2][3][4]: where g = ∇ x g ∈ R n×n is the Jacobian matrix and [·] k is the iteration counter. It is wellknown that mapping (5) has quadratic order of convergence. Unlike to the mapping (5), HONL methods present order of convergence higher than two. In following Sections, we present two families of HONL techniques with third order of convergence.

Weerakoon-Like Methods for PF Analysis
Weerakoon and Fernando developed in [30] a nonlinear solver with cubic convergence rate (labelled 3OW in this work). Application of this technique to PF analysis leads to the following map: As observed, mapping (6) requires two Jacobian evaluations and two Lower-Upper (LU) decompositions. This approximately supposes doubling the computational cost of NR. Other nonlinear methodologies with cubic convergence rate based on mapping (6), have been developed and studied in [33][34][35][36].

Darvishi-Like Methods for PF Analysis
Darvishi and Barati studied a third order Newton-like technique in [31] (labelled by 3OD in this work). This approach was already applied to PF analysis in [1], leading to the following map: Alternatively, Babajee et al. [37] concluding in (7) by similitude with the Chebyshev's method. 3OD responds to the structure of the so-called multipoint iterative methods [38]. These approaches achieve (N + 1) th order of convergence, where N is the number of steps computed (two steps in (7)). In contrast to (5), the mapping (7) only requires a Lower-Upper (LU) factorization each iteration, however, the function g has to be computed twice. For the sake of self-sufficiency, Figure 1 is a flowchart for PF solution using NR, 3OW or 3OD. As seen, the flowchart is the same for the three considered solvers. This evidences that the methods 3OW and 3OD are easily implementable in standard PF codes, thus only requiring minor modifications.

Stability of the Considered Cubic Methodologies
In order to study the stability of 3OW and 3OD, let us introduce the following definition:  In contrast, a fixed point is called unstable if it is not stable.
The Definition 1 can be used to study the stability of an iterative map by assimilating it to a dynamic paradigm [10]. This way, one requires for a PF solution (namely x * such that g(x * ) = 0) to be asymptotically stable. In addition, it is desirable that x * be a sink.

Stability of 3OW
In the case of mapping (6) one can write: By differentiating (8) with respect x one obtains: On the other hand, by evaluating the first step of (6) at x * one obtains: Then, evaluation of (9) at x * yields: where I ∈ R n×n is the identity operator. By (11), we can easily deduce that x * is asymptotically stable for the mapping (6). In addition, since all the eigenvalues of (11) have a negative real part, x * is also a sink for 3OW.

Stability of 3OD
For (7) one has that: By differentiating (12) with respect x one obtains: As observed, equation (10) is also valid for the mapping (7). Consequently, one could write: As checked by (14), x * is a sink for the mapping (7), however, it is not asymptotically stable. From the results obtained in this Section, it is concluded that the Weerakoon-like techniques are more stable than the Darvishi-like approaches. Table 1 summarizes the main computations involved in an iteration of the mappings (5)-(7). From Table 1 may look that the Newton-Raphson technique is the most efficient mapping. However, the order of convergence should be also taken into account. In order to properly consider this aspect, let us introduce the following efficiency index [39]:

A Comparison of the Considered Cubic Techniques
where, p ∈ R + is the order of convergence and CO stands for the total computational cost of an iteration in terms of flops. In this regards, the cost of a LU decomposition may be estimated using the following Theorem. Theorem 1. The number of products and quotients required for solving q linear systems of equations with the same matrix of coefficients, by using LU factorization, is: The proof of the Theorem 1 can be found in [40].
Finally, by considering the cost of each function o(g) = n and Jacobian evaluation o g = n 2 , the total computational cost of an iteration might be expressed as follows: where K 1 ∈ N is the number of function evaluations, K 2 ∈ N is the number of Jacobian evaluations, n LU is the number of LU factorizations and q m indicates how many linear systems are solved using the m th LU factorization. The value of these parameters for each studied solver were summarized in Table 1. With this information, the following Theorems are devoted on estimating the total computational cost of the maps (6) and (7), respectively.  Table 1 we have that K 1 = 1, K 2 = 2 and K 3 = 2. On the other hand, each LU decomposition is devoted on solving a linear system, therefore, q = 1 for both LU decompositions. Finally, it is well-known that p = 3 for mapping (6) (a proof is provided in [30]). Thus:

Proof of the Theorem 2. From
Therefore: which completes the proof.
Proof of the Theorem 3. From Table 1 we have that K 1 = 2, K 2 = 1 and K 3 = 1. On the other hand, the LU decomposition is devoted on solving two linear systems, therefore, q = 2. Finally, it is well-known that p = 3 for mapping (7) (a proof is provided in [31] and [37]). Thus: Therefore: which completes the proof.
The results of the Theorems 2 and 3 allow to compare the methods 3OW and 3OD with NR. In this sense, Figure 2 plots the value of (12) for the mappings (5)- (7). Clearly, 3OD is the most efficient solver, which was expected since this technique only requires one LU factorization but achieves cubic convergence. In contrast, 3OW lower efficiency index than NR. This low degree of efficiency is undoubtedly due to the extra LU factorization in comparison with (5) and (7). One should note that, especially in large systems, the factorization of the Jacobian matrix is the heaviest part of any PF calculation [10].

Numerical Experiments
This Section is devoted on presenting some numerical results, which aim to illustrate the advantages and disadvantages of the considered cubic methods in PF studies. For this purpose, the considered families of cubic methods are compared with the conventional NR and the PF solver based on the 4th order Runge-Kutta (RK4) developed in [10]. Throughout this Section, the following cases will be considered:

•
The IEEE 30-, and 300-bus test systems [41,42]. For the sake of self-sufficiency, Table 2 summarizes the main characteristics of the studied systems. Along NR and RK4 in polar coordinates, the iterative procedures defined by (6) and (7) have been coded in Matpower [46]. In all cases, g(x) ∞ < 10 −6 has been imposed as convergence criterion. The accuracy of the different techniques has been compared with the so-established correct solution of each system (calculated by using NR and the initial point defined in Matpower), concluding that solution achieved by all the considered solvers in the case studies are in essence the same operational point. In this sense, all considered techniques have been concluded to be accurate enough and a further comparative analysis have not been included as it lacks of importance.

Convergence Rates
The first test is focused on comparing the convergence rates showed by the considered methods. A priori, one can guess that 3OW and 3OD solvers will exhibit higher convergence rate than NR and RK4 due to their cubic convergence features. However, a deeper comparison among 3OW and 3OD techniques has to be empirically carried out. Thereby, Figure 3 plots the convergence profiles for the studied systems. As expected, the cubic techniques required less iterations than NR and RK4. On the other hand, 3OW clearly exhibited the highest convergence rate. It is worth mentioning that 3OD failed in the 2736-, and 9241-bus cases, which confirms, to some extent, the conclusions drawn in Section 3.

Solution Times
Now, the efficiency of the different studied solvers is analyzed. To do that, solution times employed by the different techniques for solving the cases studies are compared. Table 3 reports the solution times in milliseconds, required by the studied techniques. These results have been obtained on a 64-bit i5-9400F Intel Core personal computer (2.90 GHz, 8 GB of RAM). In order to avoid the influence of other computational activities, the reported results have been calculated as the average values of 100 simulations. As observed, 3OD was clearly the most efficient technique (even more than NR). Nevertheless, we have to remark that this algorithm failed in the 2736-, and 9241-bus systems. On the other hand, 3OW resulted not competitive w.r.t. NR as it normally takes longer solution times. To get a better overview, Figure 4 provides a comparison of the computational saving offered by 3OW and 3OD w.r.t. NR and RK4. As seen, while 3OD was able to improve by~30% the results obtained with NR, 3OW was usually inefficient in comparison with the mapping (5). As expected, computational superior of cubic techniques is more noticeable in comparison with RK4, because this solver requires 4 LU factorizations per iteration and presents a very low convergence rate. These results empirically confirm the conclusions drawn in the efficiency analysis performed in Section 4.  The reason behind the low degree of efficiency of the mapping (6), lies in the fact that it requires two LU factorizations each iteration. As commented, the factorization of the Jacobian matrix may be considered as the heaviest part of any PF calculation, due to its computational cost is o n 3 . As checked in Table 4, 3OW often required more LU decompositions than the remainder techniques. At the same time, one can check that 3OD typically computed less factorizations than NR and RK4, which explains why this technique resulted competitive compared with NR.

Influence of the Loading Level
Next, les us examine the influence of the loading level in the overall performance of the studied algorithms. To do that, we have modified the power profiles of the studied systems as follows: where λ ∈ R + is the loading level. Figure 5 shows the total number of iterations required by the different studied solvers for different loading levels. For the sake of simplicity, only results in in the 30-, and 2869-bus systems have been reported, since similar conclusions were observed by the authors in the other cases. As expected, 3OW and 3OD outperformed NR in all studied systems. It is worth observing that 3OW is less affected by the loading level than 3OD. Results for RK4 are not shown in Figure 5 for facilitating its visualization (note that this technique employed many iterations even for solving the base case).

Contractive Properties
Next, we aim to check the contractive properties of the analyzed mappings. Firstly, it is suitable to present the following Theorem: Theorem 4. Contraction map Theorem: Let x * be a root of the mapping G. If there exists a neighbourhood S(x * , δ) and a positive constant η ∈ (0, 1) such that G(x) − G(x * ) ≤ η x − x * , then the sequence {x k } ∞ k=0 defined by x k+1 = G(x k ) converges to x * for any x 0 ⊂ S.
Proof of the Theorem 4 can be found in [47]. Theorem 4 provides sufficient (but not necessary) conditions for convergence. Here, two aspects deserve to be remarked:

•
Convergence of a mapping is ensured if the following holds. • Intuitively, one can guess that the lowest η the highest degree of robustness and stability.
This way, one can analyse the contractive properties (and indirectly the stability and robustness features) of an iterative mapping by observing the value of η. This analysis has been performed for the studied solvers in various cases studies and the results are shown in Figure 6 plots the value of η for NR, RK4, 3OW and 3OD in the 30-, 2736-, and 9241-bus cases. In this case, the point x * has been calculated using the Newton-Raphson technique with the initial guess provided in Matpower with convergence tolerance equals to 10 −9 . From Figure 6 it is clearly observed that 3OW presents better contractive properties than the other solvers. Indeed, in the 2736-, and 9241-bus systems the ratio η is occasionally greater than one for NR and 3OD, while the relation (24) holds in all cases when 3OW is used. In the case of RK4, it still present good contractive properties due to the condition (24) holds in all studied cases, however, the ratio η decreased very slowly due to this solver presents a very low convergence rate.

Influence of the Initial Guess
It is well-known that NR has local convergence [47]. This entails that successful convergence of this algorithm is strongly influenced by the initial guess x 0 . Generally, one can affirm that the mapping (5) converges if x 0 lies sufficiently close to x * . However, since x * is normally unknown, choosing a suitable x 0 might be a critical task. In this Section we aim to discern if the considered cubic techniques are wider or narrower convergent that NR technique. To do that, we have performed a statistical analysis, which a set of different initial guesses is built as follows: where N (x * , σ) is a Gaussian distribution with mean x * and standard deviation σ. We have considered up to 100 different initial guesses for different standard deviations. Figure 7 plots the total number of scenarios successfully solved by different methods, in the 30-, and 300-bus cases. As observed, NR and RK4 presented similar performance and were outperformed by 3OW, which resulted to be the less sensitive method. In contrast, 3OD has turned out to be scarcely reliable. To complete the Section, we have considered the ill-posed versions of the 1345-, 2869and 9241-bus cases, which can be found in [48]. These cases present an initial point for which NR fails to converge. Thereby, considering the definition provided in [10], these cases can be considered ill-conditioned. In order to provide a better overview in such systems, we have included RK4 (10) in the analysis, since this solver could be considered robust. Figure 8 shows the convergence profiles in these cases. As observed, while NR typically showed an oscillatory pattern 3OD diverged in all the studied cases, 3OW was successful in all the studied cases. These results strengthen the ideas outlined in Section 3. It is worth remarking that, despite RK4 could be considered a robust solver, it failed in the 2869-and 9241-bus cases. Finally, comparison of the computational burden between RK4 and 3OW is addressed in Table 5 for the ill-posed version of the 1354-bus case. As seen, RK4 consumed much more time than 3OW. A further insight allows to compare the average time consumed in each iteration of both techniques. As observed, the computational burden of RK4 is much higher than that observed for 3OW, which is expected as the RK4 requires up to four LU factorizations each iteration.

Concluding Remarks
This work has presented a comprehensive analysis of two families of HONL techniques, in order to validate them for PF studies. This analysis has been supported on theoretical findings and numerical results.
This way, The Weerakoon (3OW) and Darvishi (3OD) families of cubic techniques have been considered, since they suppose the main background for developing multiple HONL approaches. This way, this work aims at serving as valuable benchmark for developing and applying other multiple HONL techniques for PF analysis in the future.
The first block of this manuscript has developed various theoretical findings on the basis of well-known theorems. This analysis has yielded as main conclusion that 3OW is more stable but less efficient than 3OD.
In a second part, various numerical experiments have been performed to confirm the theoretical findings encountered in the previous study. In these experiments, the considered cubic methods have been compared with the well-known NR and RK4 techniques. The results reported allow to assert that 3OW generally result to be more robust and presents better convergence rate compared with NR and 3OD. In addition, this technique seems to be less affected by the loading level and the initial guess. The main drawback of 3OW is its low degree of efficiency. In this aspect, it has been usually outperformed by the 3OD, which has turned out to be very efficient (even more than NR). However, 3OD has resulted very weak and totally not suitable for ill-conditioned systems.
On the light of the results obtained, we believe that the Weerakoon-like techniques may be fully competitive with the Newton's scheme, due to their intrinsic robust features. In this regard, future works should be focused on reducing its computational burden by proposing alternative schemes.

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