Mann-Iteration Process for Power Flow Calculation of Large-Scale Ill-Conditioned Systems: Theoretical Analysis and Numerical Results

Power Flow solution of realistic ill-conditioned systems has recently attracted huge attention. Nevertheless, there are still some gaps in this field. For example, most of available references do not provide exhaustive theoretical analysis about convergence properties of proposed approaches. In addition, efficient solution of large-scale ill-conditioned systems is still an open topic. This paper tackles these issues by comprehensively studying the suitability of the Mann Iteration Process for the solution of ill-conditioned systems. A comprehensive theoretical analysis is provided, from which is demonstrated that the Mann Iteration Process has with asymptotic stability, may achieve a high convergence rate and constitutes a robust methodology, improving the contractive properties of the Newton-Raphson method. Moreover, some interesting links with other Power-Flow approaches are obtained as by-product. Several numerical experiments serve to confirm the theoretical findings and to compare the performance of the Mann Iteration Process with other well-known PF solvers. In all cases, the results obtained with the Mann Iteration Process are superior to that obtained using other methodologies, being able to efficiently solve various large-scale ill-conditioned systems.


I. INTRODUCTION
A. MOTIVATION In power system analysis, Power-Flow (PF) supposes one of the most important computational tool, finding applications in control, operation and optimization of power systems, among others [1]. PF is a mathematical problem, which consists of solving the nonlinear equations that relate nodal voltages with power flows through the lines. As any other nonlinear problem, PF cases can be categorized as well or ill-conditioned [2]. While the formers are easily solvable The associate editor coordinating the review of this manuscript and approving it for publication was Youngjin Kim . using standard solvers like the Newton-Raphson method (NR) [3], the latter may suppose a challenge for conventional techniques [4].
During decades, PF ill-conditioned cases were very infrequent. However, this trend is changing nowadays and ill-conditioned PF problems are becoming more frequent in power system operation [5]. The main reason for that is an overall growing demand along the impossibility of upgrading infrastructures, which is leading the power systems to be operated close to their operability limits [6]. Nevertheless, other reasons are provoking ill-conditioning in PF cases, like the increasing deployment of electronic devices (FACTS and D-FACTS) [7], which contribute to modify the VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ line impedances so that they provoke ill-conditioning of the Jacobian matrix. Due to the reasons above, the PF solution of ill-conditioned systems is a hot topic nowadays, as evidenced by various recent works [4]. Normally, this kind of cases is addressed by using robust PF solvers, which present better numerical stability than NR, which can be considered the standard PF solution technique. Although many efforts have been made in addressing this issue (see Literature Review), most of the available literature lacks of an exhaustive theoretical analysis about the mathematical properties of the developed PF solvers. In this regard, it is necessary to check the robustness and convergence properties of the numerical techniques. As a sake of example, it is desirable that robust methodologies meet the Contraction Mapping Theorem (e.g. see [8]). However, this kind of analysis is infrequent in the related literature. Moreover, it has been well-reported in many references (e.g. see [9]), that most of available robust PF solvers are very inefficient, which has limited its application in industry tools. Finally, it is surprising that, despite the vast literature available, many promising nonlinear solvers have not been explored in PF studies yet. This paper aims to give a first step at filling the gaps above.

B. LITERATURE REVIEW
The first steps in PF solution of electrical networks were given at 60-70's, with the development of solution routines based on Gauss-Seidel [1] and NR [3]. Posteriorly, some authors focused on the inefficiency problems caused by Jacobian calculation in Newton-based solvers thus applying sparsity routines [10] and decoupled techniques [11], [12]. Since then, many advances have been made in the PF calculation from different perspectives, developing linear models [13]- [15], applying or developing high order Newton-like methods to the PF analysis [16], [17] or deriving alternative formulations [18]- [20].
Ill-conditioning in PF problems was firstly reported by Stott [21] and firstly addressing by Sasson [22] and Iwamoto [23]. Both authors aimed at improving the convergence characteristics of NR by posing minimization techniques. However, the algorithm developed by Sasson suffered of frequent local trapping issues. On the contrary, the method posed by Iwamoto resulted very effective for the first studied ill-conditioned cases, which motivated different developments based on this technique [24], [25]. Certainly, the Iwamoto's method worked well in the first studied ill-conditioned cases due to its small size (< 100 buses); however, this method has been reported to be quite inefficient in realistic large-scale ill-conditioned cases [2]. Moreover, it also suffers of local trapping issues [1]. Other authors tried to address PF calculation in ill-conditioned cases by developing second Taylor expansions of the PF equations [26], [23]. However, this idea was rapidly abandoned due to the high computational cost of calculating the Hessian matrix of the PF equations.
During 80-90's PF ill-conditioned cases were few studied because they were very infrequent in real life. However, this trending changed in early 2000's, due to the increment of ill-conditioned problems observed in real cases. This fact motivated Milano to apply the Continuous Newton's method to PF studies in [2]. The Continuous Newton's method establishes a common framework between systems of nonlinear equations and autonomous dynamic systems. On the basis of this analogy, it is determined that any methodology suitable for solving dynamic systems, can be successfully applied for solving PF problems. This was demonstrated in [2], where two solvers based on the Forward Euler method (FE) and the 4 th order Runge-Kutta technique (RK4) was developed. Results obtained with these solvers in realistic ill-conditioned cases were promising, outperforming other classical approaches. These good results motivated the authors to further apply the Continuous Newton's method by developing solvers based on the Adams-Bashforth's methods [27], Bulirsch-Stoer algorithm [28] and Runge-Kutta formulas [9]. Despite their apparent advantages, the solvers based on the Continuous Newton's method may be still inefficient in realistic cases and totally not competitive with NR in well-conditioned problems. Indeed, while NR only requires one Lower-Upper (LU) factorization each iteration, RK4 requires four. In addition, the convergence order of these methods is frequently linear, while NR converges with quadratic converge rate.
The other category of solvers is based on the Levenberg method [6], [7], [29]- [32]. By this technique, stability of NR is improved by intentionally well conditioning the Jacobian matrix. This is achieved by introducing the so-called damping factor and the regularization matrix. The methods based on the Levenberg procedure have a computational burden comparable to NR, requiring only an LU decomposition each iteration, and their ability to manage with ill-conditioned cases has been empirically proved [7], [31]. However, these techniques may converge to a non-physical solution, as reported in [33]. This is due to the damping factor modifies the PF equations when its value is high. However, a high value of the damping factor is necessary to ensure the stability of the algorithm and the full-rank of the Jacobian matrix. This issue is normally addressed by developing adaptive mechanisms for the damping factor, so that it is firstly fixed high to posteriorly be reduced as the algorithm evolves. Nevertheless, there is not still a unified a formal framework for the treatment of the damping factor in Levenberg-based solvers.
Instead of using conventional nonlinear iterative solvers, other authors have tried to address PF ill-conditioned problems from alternative point of views. Such is the case of the metaheuristic-like solvers [34]; Continuation approaches [6], [35], [36] and homotopy/holomorphic algorithms [37]- [39]. Nevertheless, these solvers still present serious problems. For example, continuation and metaheuristic solvers have been reported to be computationally inefficient [40], while the holomorphic approaches may lead to non-physical solutions [41]. In addition, this kind of methods normally presents a rigid coding structure, which hinders the simple integration of electronic devices, alternative formulations or control routines.
Some PF solvers avoid the calculation of the Jacobian matrix and thus its factorization. Such is the case of the Forward-Backward Sweep algorithm [42] and the implicit Z-bus formulation [43], [44]. These techniques have been successfully applied to radial (or weakly meshed) distribution networks. However, while the Forward-Backward Sweep method can be only applied to radial networks and frequently requires many iterations to converge (see results in [45]), convergence issues have been reported for the implicit Z-matrix approach in systems with many PV buses [46].

C. CONTRIBUTIONS AND PAPER ORGANIZATION
A careful analysis of the most recent contributions confirms that efficient calculation of ill-conditioned systems is still an open topic, especially in realistic large-scale systems. Despite that, the authors have detected that scarce efforts have been made on applying some promising classical nonlinear iterative solvers in PF analysis. In this regard, the Fixed Point Iteration methods (FPI) may offer very good results. This supposition is founded on the enhanced convergence characteristics showed by this kind of techniques, which may be frequently encountered in a huge amount of studies available in the literature (see e.g. [47] and references therein).
To the best of our knowledge, only the reference [48] has covered the application of FPI in PF analysis. However, this work is only limited to pose NR and the classical Gauss-Seidel techniques in the form of FPI. Consequently, applicability of other FPI in PF studies has not been studied yet. With the aim to fill this gap, this paper deeply analyses the applicability of the Mann Iteration Process (MIP) in PF analysis. In order to be rigorous, some relevant theoretical results are found and posteriorly experimentally confirmed. As a by-product, some interesting links with other existing approaches are encountered. Suitability of MIP for PF analysis is also discerned by comparing its robustness and efficiency with other popular robust PF approaches in solving several practical large-scale ill-conditioned systems. From a methodological perspective, this paper contributes by developing a novel PF solver based on MIP, which is competitive with other well-known techniques, which is demonstrated by numerical results. For the sake of summarize, Fig. 1 shows the flowchart following for the development of this paper.
The remainder of this paper is organized as follows. Necessary background for properly tackling the targets of this paper is summarized in Section II. Convergence properties of MIP are analysed in Section III. Section IV presents a study about the stability of MIP. Section V presents some relevant comments about the so-called evolution parameter (µ). Some interesting links with other PF techniques are commented and illustrated in Section VI. Several numerical experiments with results are reported in Section VII. Finally, the paper is concluded with Section VIII.

A. PF PROBLEM AND ILL-CONDITIONED CASES
The PF problem in polar coordinates (see Appendix A for details), is defined as a set of n nonlinear equations as follows: where, g ∈R n −→ R n are the PF equations and x ∈ R n is the PF state vector. When the system (1) is well-conditioned, its solution can be easily found by using conventional iterative techniques like NR. However, when the problem is ill-conditioned, most of PF solvers might experience convergence difficulties. A formal definition of an ill-conditioned system is given below.

Definition 1 [2]
: Ill-conditioned system: solution of PF does exist, however, it is not reachable using NR and a flat start (e.g., all load voltage magnitudes equal to 1 and all voltage angles equal to 0).

B. THEORETICAL BACKGROUND
Some relevant definitions and Theorems are provided in this section.
Definition 2 [8]. Nonexpansive and Contractive Mappings: let C ⊂ R n be a nonempty convex subset of a real Banach space. A mapping T : C ⊂ R n → R n is called: Hyperbolic points: a fixed point x * of a map T : C ⊂ R n → R n is hyperbolic, if the Jacobian of the equationṪ = T (x) − x at x * has no eigenvalues with a zero real part. In addition, a hyperbolic point is asymptotically stable if all the eigenvalues ofṪ at x * have a modulus ≤ 1. A hyperbolic point can be: • Sink: if all the eigenvalues of the Jacobian ofṪ has a negative real part.
• Source: if all the eigenvalues of the Jacobian ofṪ has a positive real part. VOLUME 9, 2021 One expects that fixed points of an iterative map to be sinks, in order to ensure its stability. The following Theorem provides a sufficient condition about a sequence {x k } ∞ k=0 , defined by an iterative map x k+1 = T (x k ), for converging to x * Theorem 1 [8]. Contraction-Mapping Theorem: suppose that T : C ⊂ R n → R n maps a closed set C 0 ⊂ C into itself and that: where, α < 1 is the Lipschitz constant. Then, for any x 0 ∈ C 0 , the sequence x k+1 = T (x k ) converges to an unique fixed point x * of T in C 0 , and holds. By Theorem 1, it is established that contractive mappings present better convergence properties than expansive ones, since the latter do not meet the premises of Theorem 1. Moreover, it is worth mentioning that Theorem 1 establishes sufficient rather than necessary conditions, i.e. the sequence x k+1 = T (x k ) may converge to x * even if T is not contractive.

C. MIP FOR PF ANALYSIS
With the aim of finding fixed points of nonexpansive mappings, Picard introduced the following iterative formula [49]: where the subscript stands for the iteration counter. Clearly, most of PF solvers can be established in form of Picard iteration. Thus, as sake of example, a NR iterative procedure for solving the PF equations may be written as follows: where, g = ∇ x g ∈R n×n is the PF Jacobian matrix. It is well-known that Picard iterations of some nonexpansive mappings fail to converge even on a Banach space. In order to enhance the convergence properties of (4), Krasnoselskii introduced the following FPI [50]: Intuitively, a more general version of (6) can be established introducing a real parameter µ instead of taking it fixed. This idea was explored in [51] arising in the so-called Mann Iteration Process (MIP), which is given by: In this paper, µ will be called evolution parameter. As clearly seen, (6) and (7) are equivalent for µ = 1 2 . Now, it is worth questioning if the solution of PF (namely x * such that g x * = 0) is a fixed point of (7), which is proved by the following Theorem.
Then, x * is also a fixed point of the mapping G : C ⊂ R n → R n defined by (7).
Proof: Let us consider.
Since x * is a fixed point of T one can write.
Rearranging (9), the following equality is easily obtained: Which completes the proof. Trivially, it can be deduced that findings encountered in Theorem 1 are also applicable to (6), since it is a particular case of (7). It is worth mentioning that Theorem 2 establishes the necessary conditions. Thereby, one can also claim that MIP does not converge to x * if it is not a fixed point of T . Therefore, by Theorem 2, it is easily concluded that x * has to be a fixed point of T , since G would not constitute a PF solver otherwise.

III. CONVERGENCE ANALYSIS
In this Section, the convergence order and local convergence of MIP are studied.

A. ORDER OF CONVERGENCE OF MIP
The order of convergence is the characteristic of any FPI that defines how quickly the iterative procedure converges to a fixed point. For instance, it is well-known that NR (5) has a quadratic order of convergence while other methods present cubic or higher convergence rates [16], [17]. The following Theorem is devoted to studying the order of convergence of MIP, for which the Taylor expansion technique has been used (e.g. see [52] for further details).
Theorem 3. Order of Convergence of MIP: let T : C ⊂ R n → R n be a map with order of convergence p. Then, p is the maximum order of convergence attainable by the map G : C ⊂ R n → R n defined by (7), and it is achieved with µ = 1. Otherwise, its convergence remains linear.
Proof: Let us define the following error function: Applying (11) to MIP gives: Since T has an order of convergence p, one can write: The proof is completed. By Theorem 3, it is claimed that the convergence order of MIP is ruled by the mapping T . In this regard, MIP is not devoted to accelerating the convergence of T . Nevertheless, in contrast to most of available robust PF solvers, which present linear convergence (see e.g. [2] or [23]), MIP can take advantage of convergence features of the mapping T , in order to achieve high convergence speed. For example, let us suppose that the mapping T is defined by NR (5), then, MIP may achieve up to quadratic convergence. However, this finding does not suppose any advantage with respect to NR or other available solvers. Nevertheless, as shown in Section III.B, MIP makes wider convergent T . In other words, MIP is introduced for enhancing the robustness characteristics of T , while its convergence features may be preserved and achieved by just manipulating the evolution parameter. Immediately, one can guess the main advantage of MIP compared with (6), since the order of convergence of the latter always remain linear.

B. ON THE CONTRACTIVE PROPERTIES OF MIP
Next, let us study the contractive properties of MIP.
Theorem 4. Contractive properties of MIP: let T : C ⊂ R n → R n be an inexpansive map. Let G : C ⊂ R n → R n be a mapping defined by (7). Suppose that both T and G map a closed set C 0 ⊂ C into itself. Then, the following inequality holds: Proof: By observing (7), one can write: By triangular inequality, the following condition holds: Since T is nonexpansive, the worst case would be T (x) − T (y) = x − y . Keeping this in mind, the equality (17) holds Consequently, one can write: Which completes the proof.
In the light of the result of Theorem 4, the following ideas can be deduced: • If T converges, G also converges, since contractivity on T necessarily entails contractivity on G.
• If T diverges, G may still converge. One can observe that inequality (18) may be strict even when T (x) − T (y) = x − y . In other words, G may be contractive when T is strictly inexpansive. Consequently, the mapping defined by (7) could meet the premises of Theorem 1 even when T does not.
• A direct deduction from the preceding point leads to assert that mapping G is always wider convergent than mapping T .
• It is worth remarking that MIP is not globally convergent. Theorem 4 establishes sufficient conditions for convergence, if and only if x 0 ∈ C 0 (see Theorem 1). From the points above is deduced that MIP constitutes a robust map, at least, more robust than T , which makes it suitable for PF calculation of ill-conditioned systems. Now, let us enunciate the following conjecture, which will be empirically validated in Section VII.
Conjecture 1: Lipschitz constant of the contractive mapping (7): let G : C ⊂ R n → R n be a map defined by (7). Consider a closed set C 0 ⊂ C If the following inequality holds, then G converges to an unique fixed point x * ∈ C 0 . Equation (19) is an indirect way to calculate the Lipschitz constant α of MIP. It is worth mentioning that inequality (19) holds just in C 0 where MIP is contractive. In addition, (19) does not suppose a necessary condition for convergence of MIP.

IV. STABILITY OF MIP
The stability of MIP at x * such that g (x * ) = 0 is studied using the Definition 3, by assimilating it to a dynamic system (see [2] for further details). Without loss of generality, the map T defined by (5) is considered.
Theorem 5. Hyperbolic points of MPI: let G : C ⊂ R n → R n be a map defined by (7). Let T : C ⊂ R n → R n be a map defined by (5). Suppose that x * is a fixed point of G such that g (x * ) = 0. Then, all the eigenvalues of the Jacobian of the functionĠ = G (x) − x at x * have a real part equal to µ Proof: the functionĠ of MPI is given by: By differentiating (20) with respect to x, one obtains: where, I ∈ R n×n is the identity matrix. In the case of T is defined by (5), one can write: Substituting (22) into (21), one obtains: Evaluating (23) at x * ∇ xĠ x * = µ −∇ x g x * −1 g x * − I = −µI (24) The proof is completed. Result of Theorem 5 deserves various comments: VOLUME 9, 2021 • If the procedure described in (20)-(24) is applied to the map defined by (5), it is obtained that all eigenvalues at equilibrium point has a real part equal to −1 (see [2, Sec. III.A]). Therefore, if µ k < 1, then it can be asserted that all eigenvalues of MIP lie closer to the origin. From this assertion, it is deduced that MIP is generally more stable than NR.
• It is worth noting that (24) only holds at x * . For the remainder iterative procedure, one should consider (23), from which is deduced that all the eigenvalues of its associated Jacobian have a real part proportional to g (x) −1 g (x). Here, the evolution parameter plays a crucial role in order to counteract the effect of these factors and keep the map G as stable as possible.

V. CONDITIONS ON THE EVOLUTION PARAMETER
Now, some conditions about the evolution parameter are established. Typically, the following two conditions are imposed to µ (see e.g. [53]): i.
k < ∞ Basically, conditions (i) and (ii) force µ ∈ (0, 1). However, in PF analysis, other conditions can be deduced from the results obtained in the provided Theorems: a) In order to achieve high convergence speed, µ should be increased to be equal to 1 (Theorem 3). b) In order to keep the equilibrium point x * asymptotically stable, the evolution parameter has to be positive. In addition, µ k ≤ 1, ∀k(Theorem 5). c) In order to keep the mapping (7) as stable as possible, the evolution parameter, at least for initial iterations, should be inversely proportional to g (x) −1 g (x) and, in all cases, µ 0 < 1 (Theorem 5). From a-c, the following three conditions are imposed on the evolution parameter: iii. lim It is worth noting that conditions (iii)-(v) are more specific than (i)-(ii). The formers have been imposed in order to obtain the best performance of MIP in PF analysis. For example, although by Theorem 5 it is recommend that µ k < 1, this point may be considered too much restrictive, since MIP reaches its maximum order of convergence for µ = 1 (see Theorem 3). Attending to deductions (a)-(c) and conditions (iii)-(v), the following simple rule is proposed for updating the evolution parameter.
It can be easily observed that (25) addresses the conditions (i)-(v).

VI. LINKS WITH OTHER APPROACHES A. LINKS WITH EULER METHOD AND ZHANG DYNAMICS
Equation (20) can be conceived as a set of autonomous differential equations as follows: Again, if map T is defined by (5), the function (26) can be written as follows: Rewriting (27) in the form of discrete map, one obtains: Equation (28) corresponds to the integration of function (27) using the Forward Euler methodology taking µ k as step size. Immediately, one can find a link with the Zhang dynamics [54]. In this case, (28) corresponds with [55, eq. (13)] using linear activation function and a step size equal to µ k .

B. DIFFERENCES WITH HOMOTOPY
A further insight on the equation (7) may lead to deducing that MIP is actually a homotopy method, due to its striking similitude with the linear homotopy (see [56,Sec. I]). Nevertheless, if the conditions marked in [57,Th. 3] are observed, it can be deduced that MIP is not a homotopy method, since some conditions of this Theorem are not met by (7). An even clearer difference can be found using the Davidenko's method [57], which computes the variation of x as a function of the evolution parameter as follows: By evaluating (29) at x * for equation (7) one obtains: Equation (30) says that, at the equilibrium point x * , the state vector x does not vary. In other words, x * is a fixed point of (7), which is a distinctive feature of FPI maps. On the other hand, if one evaluates (29) at x * for the linear homotopy using, for example, the Newton formula for the easier system [55], the following result is obtained: Equation (31) indicates that x * is not a fixed point of the linear homotopy, in other words, the linear homotopy is not a FPI. Indeed, homotopy methods find the solution of g (x) following a connected curve between two different systems (easy and difficult ones). Therefore, unlike FPI which undoubtedly describe an iterative map, homotopy techniques does not actually iterate. In that sense, MIP clearly defines an iterative map rather than a homotopy. Algorithm 1 PF Solution Procedure Using MIP 1: Let x 0 , ε, k max and µ 0 < 1 be given 2: Initialize iteration counter k ← 0 3: while g (x k ) ∞ ≥ ε do 4: x k+1 ← Solve (7) 5: if g (x k+1 ) ∞ < ε 6: break #Convergence 7: endif 8: µ k+1 ← solve (25) 9: k ← k + 1 10: if k > k max 11: break #Fail 12: endif 13: end 14: return solution x k

VII. NUMERICAL EXPERIMENTS
This Section presents various numerical results which aim at validating MIP for PF calculation, along to check their mathematical characteristics empirically. To this end, the PF calculation using MIP is summarized in Algorithm 1 using pseudocode. Here, g ∞ has been imposed as convergence criterion. Moreover, the iterative procedure is considered failed if the iteration counter k surpasses a predefined threshold namely k max , which normally indicates divergence. Similar procedures have been used for other solvers considered in the simulations, simply replacing the MIP mapping (7) by the corresponding numerical iterative scheme.

A. EMPIRICAL VALIDATION OF CONJECTURE 1
This Section is devoted on empirically validating the Conjecture 1. To do that, the ratio η is evaluated for different values of the evolution parameter. In this case, µ has been taken constant during the iterative procedure, hence line #8 of Algorithm 1 has been omitted. Fig. 2 plots the value of η manifested when MIP with the mapping T defined by (5) is used for solving the IEEE 9-bus system [58], from a flat start during the first 10 iterations. In all cases, MIP successfully converged. Clearly, it is observed that in all cases η ≈ 1 − µ after 5 ÷ 6 iterations. Next, the performance of MIP is analysed in the 3012-bus snapshot of the Polish Transmission system at winter 2007-08 evening peak [59]. This system has been documented as ill-conditioned in some references (see e.g. [28] or [33]). Fig. 3 is analogue to Fig. 2 for this case. As observed, for µ = 0.7 frequently η > 0.3, which has supposed divergence. In remainder cases, η often traced a random pattern for initial iterations. Nevertheless, it tends to be stable and equal to 1 − µ, which has led to successful convergence.
From the results reported above, Conjecture 1 may be empirically validated. Basically, this conjecture establishes that the Lipschitz constant of MIP is equal to α = 1 − µ (see Theorem 1). It is also remarkable that the Theorem 1 and Conjecture 1 establish sufficient rather than necessary conditions for convergence. As observed in Fig. 3, MIP is able to converge even when eventually η > 1. In this latter case, one can simply deduce that MIP is not contractive during the whole iterative procedure.
NR is typically considered as the most standard PF solution solver, while the Iwamoto's method is often considered as a benchmark robust technique. On the other hand, both RK4 and RBS might be considered representative of modern robust PF techniques. It is worth mentioning that both RK4 and RBS involve several parameters within their iterative procedures. In this work, these parameters have been tuned following the guidelines provided in their respective references. Furthermore, all considered solvers have been coded in Matpower v7.0 [60].
The 3012-bus snapshot of the Polish Transmission system at winter 2007-08 evening peak [59], and the 13659-bus portion of the European Transmission system from the PEGASE project [61], [62], are used as benchmark cases. These systems have been reported as ill-conditioned in some references (see e.g. [28] or [33]). In addition, they are clearly large-scale and can be considered realistic enough for covering the targets of this work. Table 1 reports the total iterations required by the considered PF techniques for solving the studied systems from a flat start. As expected, NR failed in all cases due to the ill-conditioning character of the studied systems. RK4 failed in solving the 3012-bus case. On the other hand, the Iwamoto's method converged to the low voltage solution instead of the high voltage one in the 13659-bus system. Oppositely, RBS, Kranoselskii and MIP converged to the stable (high voltage) solution in all the studied systems. In all cases, MIP required less iterations than remainder techniques. It is worth mentioning that the results obtained in the studied systems may be considered a confirmation of Theorem 4, since divergence of NR (mapping T ) has not necessarily implied divergence of MIP.  Fig. 4 plots the convergence profile and the value of the evolution parameter in the 3012-bus system using MIP with ε = 10 −6 . It is clearly seen that MIP achieves quadratic convergence when µ k = 1, which may be considered as an empirical validation of Theorem 3. To get a better overview on the convergence features of the different PF solvers, Fig. 5 plots the convergence curves of the studied methods in the 13659-bus case with ε = 10 −3 . As seen, NR failed due to the residual grew during a large number of iterations. On the other hand, RBS showed the best convergence rate before the 8 th iteration, when the convergence of MIP was accelerated because the increased value of the evolution parameter, as appreciated in Fig. 3, thus leading to quadratic convergence, outperforming RBS and converging in just 8 iterations. The performance of the studied solvers have been also validated under stressing conditions. To that end, Table 2 reports the total number of iterations for limit load cases. As in other references (e.g. see [28] or [63]), the loading level in the studied systems has been progressively increased until achieving the operability limit. In this case, NR failed in all the studied systems, RK4 successfully solved the 3012-bus case, however, it converged to the low voltage solution in the 13659-bus case while RBS failed to solve the 13659-bus case. In contrast, the Iwamoto's method along the Krasnoselskii's solver and MIP were able to correctly solve all the studied systems. Convergence features of MIP were deteriorated because the stressing operating conditions, employing more iterations compared to the base case. Nevertheless, it is less affected by the convergence tolerance than other methods. This is the case of RBS, which increases its iteration counter by three when the convergence tolerance was reduced to 10 −6 , while the iteration counter of was only increased by one. Moreover, MIP has turned out to be more reliable correctly solving all the studied cases while RBS failed in the 13659-bus system.
It is also interesting how the studied PF solvers perform in the presence of control limits. On typical example are the reactive limits of generators (PV buses). To this end, a simple PV-PQ switching mechanism [64] can be incorporated to the solution routines of the studied techniques, which basically consists on converting to PQ those PV buses in which a reactive limit has been violated, taking the hit limit as reactive power injection of the new PQ bus, thus repeating the PF solution procedure until achieve a feasible solution. Table 3 reports the total iterations for this case. As expected, NR failed in all the cases. RK4 failed in the 3012-bus case while the Iwamoto's solver converged to the low voltage solution in the 13659-bus case. The remainder solvers successfully converge in all the considered cases. In this scenario, MIP showed superior convergence features thus notably outperforming the remainder solvers.
Along with the iteration number, the execution time is widely used as indicator for comparing PF solvers. In this regard, all considered methods have been run under Windows 10 on a 64-bit i5-9400F Intel Core personal computer (2.90 GHz, 8 GB of RAM). Table 4 provides the    execution times (in milliseconds) for base cases. These times have been calculated as the average value of 100 simulations, in order to avoid the influence of other computational activities. As observed, MIP was the most efficient technique. The reason behind these good results, is the low computational cost required by MIP. In this sense, the computational burden of a numerical solver is quite proportional to the total number of LU factorizations computed in the solution procedure, since it is by far the heaviest calculation [2]. In this sense, MIP can be considered quite efficient, only requiring a LU factorization each iteration being so comparable to NR. Table 5 reports the total LU factorizations required by the different solvers for solving the base cases of the studied systems. As seen, MIP required much less computations than the others methods. Same conclusions can be extracted for the other scenarios analysed, whose results are reported in Tables 6 and 7.

VIII. CONCLUSION AND FUTURE WORKS
This work aims to fill some gaps found in the available literature about PF calculation of realistic large-scale ill-conditioned systems. Firstly, it was identified that efficient calculation of this kind of systems is still an open topic. Secondly, in most of the available literature, the robustness properties of different PF approaches are empirically demonstrated on the basis of the results obtained, without providing sufficient theoretical analysis, which it is believed to be fundamental. Finally, despite that a vast literature about solution of systems of nonlinear equations is available, the applicability of many promising classical nonlinear solvers in PF analysis has not been studied yet.
In this paper, first step at filling this gap is given by profusely studying the suitability of MIP for PF calculation of ill-conditioned systems. Firstly, a rigorous theoretical analysis has been provided by which has been proved that MIP presents the same order of convergence of T , its contractive properties have been derived and asymptotical stability proved. Secondly, some interesting links and differences with other approaches have been commented. In this sense, it has been demonstrated that MIP has some similtudes with homotopy approaches, however, in contrast to homotopy methods, MIP supposes an iterative mapping. Finally, several numerical experiments served for empirically confirming the theoretical findings and to compare the performance of MIP with other well-known PF solvers. The results obtained have confirmed that MIP constitutes a robust methodology with superior contractive properties compared with other benchmark robust methodologies. In addition, its efficiency has been clearly manifested, outperforming other conventional and modern PF approaches.
Due to the promising results obtained, MIP should be taken into account and deeply explored for PF calculation of realistic large-scale ill-conditioned systems. In that sense, its suitability in other related problems like the Continuation Power Flow [35] or contingency analysis [65] should be analysed in future works. Maybe, the main issue relative to the PF calculation using MIP is its local convergence, which entails sensitivity with respect the initial guess x 0 . Some globally convergent mechanisms might be applied within the MIP procedure in future studies in order to overcome this difficulty. Finally, this work has been limited to study the mapping T defined by (5), however, any other mapping which the PF solution is a fixed point may be used. This idea opens the door to consider other alternative mappings like the high order Newton-like techniques [16], [17], [45].

APPENDIX A. PF EQUATION IN POLAR COORDINATES
In polar form, voltage angles at PV and PQ buses, along voltage magnitudes at PQ buses constitute the PF variables, so that, the PF state vector is given by: where, δ PV ∈R n g is the vector of voltage angles at PV buses, δ PQ ∈ R n g is the vector of voltage angles at PQ buses and V PQ ∈ R n l is the vector of voltage angles at PQ buses. On the other hand, n g ∈ N and n l ∈ N represent the total number of PV and PQ buses, respectively.
The nonlinear relationships between the power nodal injections and nodal voltages are given by: g (x) = g P , For all buses g Q , For PQ buses (33) where, 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.