A Novel Power Flow Solution Paradigm for Well and Ill-Conditioned Cases

This paper develops a novel four-stage power flow solver for ill-conditioned systems. Although the developed solver could be considered efficient, it is not competitive with the Newton-Raphson method in well-conditioned cases. With the aim of being fully competitive in a wide range of cases and scenarios, the developed algorithm is integrated within a novel efficient solution paradigm. As a result, a robust and efficient solution framework, competitive in both well and ill-conditioned cases, is obtained. The new proposals are tested in various well and ill-conditioned cases from 30-, to 13,659-buses. Results obtained with the developed solvers are promising.


I. INTRODUCTION A. MOTIVATION
Power flow calculation in ill-conditioned cases has been profusely studied for decades. As a result, a plethora of robust techniques is available in the literature. However, very few of these methods have found widespread usage in industry, and, currently, the conventional Newton-Raphson (NR) is still assumed as the most standard power flow solver. This is due to the fact that very few robust methods are available, which cannot compete with NR, especially in large-scale networks.
Well-conditioned systems are the most common situation in power system operation; thus, NR has offered a satisfying performance in most cases. However, ill-conditioned cases are nowadays more common due to the observed increasing energy demand in modern power systems and the difficulties in upgrading the network infrastructures. These facts currently are forcing the power components to operate close to their functional limits [1], [2]. Consequently, robust power flow solvers may presumably play a vital role in a future power system paradigm.
The associate editor coordinating the review of this manuscript and approving it for publication was Gokhan Apaydin .

B. LITERATURE REVIEW
Power flow is likely the most important computational tool in power system analysis. It finds a wide variety of applications such as voltage stability analysis [3], optimal power flow [4] or state estimation [5], among others. Power flow calculation has been covered in multiple works since its first formulations in the mid-'50s [6]. Some references classified power flow problems into three different categories [7]: a. Well-conditioned cases: this is still the most common situation. In these cases, the conventional solvers, such as NR, converge to an accurate solution from a flat start. These cases do not entail convergence difficulties for most solvers. For this kind of system, apart from the conventional NR, a variety of decoupled techniques [8], linear models [9], hybrid AC/DC solvers [10] and high order Newton-like methods 11] have been proposed. b. Ill-conditioned cases: despite the fact that these systems are still solvable, conventional NR fails to reach the solution from a flat start [7]. Robust techniques have proved to be efficacious for solving ill-conditioned cases. c. Bifurcation points: on these points, the power flow solution does exist; however, the Jacobian matrix is singular or near singular, which provokes instability of 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/ the Newton-based methods. Continuation techniques 3] or Regularization schemes [12] prevent the singularity of the Jacobian matrix at bifurcation points. d. Unsolvable cases: in this case, the power flow solution does not exist. This is typically due to the loading level has surpassed the maximum loadability point of the system. Different approaches devoted to providing approximated solutions [13] or proposing countermeasures for restoring solvability [14] have been developed for unsolvable cases. This paper is focused on ill-conditioned cases. The power flow solution of this kind of case came to be studied at the end of '70s [15], [16] by the so-called optimal multiplier methods. These first works proposed a minimization approach derived from a second-order Taylor expansion of the power flow equations in rectangular coordinates. As a result, a so-called optimal multiplier is calculated, which modifies the Newton's increment vector to be conducted in a direction where the residuals are reduced. The optimal multiplier solution methods are, theoretically, no divergent. However, the calculated optimal multiplier tends to be very short or almost null when the solution is approached [17]. This undesirable behaviour nullifies the effect of the increment vector and provokes that the optimal multiplier solution methods might be trapped on a local minimum. The optimal multiplier techniques are not suitable for the polar form of the power flow equations due to the presence of transcendental functions [18].
The optimal multiplier techniques can be included within the family of line search techniques or trust-region methods [19]- [21]. These methods aim to be globally convergent; however, this important feature is typically achieved by deteriorating their convergence properties. Thus, these solvers frequently present a notably linear convergence, which supposes an important drawback w.r.t to NR or other conventional solvers. In addition, if the starting guess used for initializing the iterative procedure is far away from the solution, the convergence characteristics of the Newton's descent vector are lost, and the algorithm may be trapped on a local minimum [22].
Another category of robust techniques is based on the Continuous Newton's method [23]- [25]. This paradigm poses a formal analogy between the Forward Euler method and most of the available optimal multiplier approaches. On the basis of this analogy, Milano deduced that any integration routine (e.g. the family of the explicit Runge-Kutta formulas), can be considered for developing robust solution methods [7]. A solver based on the 4 th order Runge-Kutta formula was developed in this same reference, which was demonstrated to be more competitive than the optimal multiplier method in [16]. The Continuous Newton's philosophy has been exploited for developing a variety of robust solvers based on different Adams-Bashforth methods [26], Bulirsch-Stoer algorithm [27] and Runge-Kutta formulas [28]. As justified in [7], only the explicit formulation of the Continuous Newton's method is attractive for power flow analysis, due to the implicit formulation would require the factorization of the Hessian matrix. However, the same author recently proposed an implicit form in [20]. Nevertheless, in this latter reference, it has been proved that the effect of the Hessian matrix is immaterial and, consequently, the explicit formulation of the Continuous Newton's paradigm is by far more competitive than the implicit one. Those methods based on the Continuous Newton's paradigm are widely rather than globally convergent. In addition, they share some of the disadvantages of the Newton's based techniques, and for example, they fail at turning points or unsolvable cases. Also, these solvers are quite inefficient compared to NR; one should note that a power flow approach based on a Runge-Kutta formula would require as matrix factorizations as the order of the integration technique considered. Thus, the low order Runge-Kutta methods are usually more competitive than the high order ones [28].
The Levenberg-type methods have recently gained popularity for power flow analysis [2], [12], [29]- [31]. These techniques are, in fact, modifications of NR, which aim to avoid the vulnerability of the conventional technique at bifurcation points [32], [33]. It is achieved by adding some elements on the diagonal of the Jacobian matrix. The most common approach is based on adding the identity operator to the Jacobian matrix. The elements of this identity matrix are pondered by a series of damping factors, which ideally must approach zero as the algorithm evolves; thus, usually, an adaptive mechanism is preferred [2], [12], [30], [31]. The series of damping factors are the most critical aspect of the Levenberg-type technique. As commented in [29], the damping factor's value strongly influences the accuracy of the results obtained by the Levenberg method. Also, the series of the damping factors may lead to a very slow mapping if they are not carefully selected; however, a unified criterion for tuning it does not exist. It is worth mentioning that the Levenberg algorithm is more computationally expensive than NR, as an extra matrix-matrix product and an extra matrix-vector product have to be computed each iteration.
In some references, the power flow equations have been raised as an artificial dynamic system [1], [34], [35]. The equilibrium points of the proposed dynamic systems are, in fact, the solution of the power flow equations, so that solving the resulting set of ordinary differential equations is equivalent to solve the set of nonlinear algebraic power flow equations. By making use of the Lyapunov theory, it is demonstrated that the equilibrium points of the equivalent set of differential equations present an exponential asymptotic convergence. Thus, the resulting algorithm is barely affected by the starting guess; hence the obtained mapping is widely convergent. Nonetheless, the resulting algorithm has to be solved using costly calculations like the ode routines in Matlab, making this approach impractical for large realistic systems [1].
Apart from the solvers referenced above, a series of techniques primarily devoted to distribution systems have been studied for decades. Some of them are based on the Backward-Forward and ladder algorithms [36], [37], which iteratively ''sweep'' the network until the difference between two consecutive solutions becomes within a specific tolerance. These approaches are based on linear operations; therefore, they are barely affected by some network characteristics like the R/X ratio. Nevertheless, despite that they avoid heavy calculations like matrix factorizations, these techniques are usually time-consuming due to their strong linear convergence, which is especially relevant in heavy loading systems (see [38,Sec V.D]). In addition, this kind of method presents difficulties to manage with PV buses even with weakly meshed networks [39]. Another family of power flow solvers for distribution networks is, in fact, based on the well-known Gauss method. The most famous approaches of this family are the Z-implicit [40], [41], and the loop algorithms [42]. These kinds of methodologies normally are efficient than the Backward-Forward and ladder algorithms. In addition, they can easily accommodate meshed networks. Nonetheless, issues for incorporating PV buses are still found, as pointed out in [43].

C. CONTRIBUTIONS AND PAPER ORGANIZATION
Motivated by a constant effort to develop efficient and robust power flow solvers, this paper proposes a novel robust four-stage algorithm (4SA). Although it is quite efficient and competitive with the most available robust solvers, it is not competitive with NR in well-conditioned systems. For expanding the suitability of 4SA, it is also integrated within a developed solution paradigm.
The developed techniques find some similitudes with the optimal multiplier or line search techniques. For example, our proposals are also based on the truncation of the Newton's increment vector (i.e. modifying the increment vector by multiplying it by a step size value). However, the introduced paradigms aim at overcoming some of the difficulties posed by this kind of technique. Specifically, they are developed with the aim to avoid the slow convergence and locally trapping phenomena that some robust algorithms (e.g. [16], [19] or [28]), may experience when the starting guess is not suitable or lies far away from the solution. Unlike other line search techniques based on finding the largest step-size that satisfies the Armijo-Goldstein conditions, our proposals aim to find the most suitable truncation (step-size) of the Newton's increment vector using a novel search strategy over the state vector space. Thus, the a priori truncation of the Newton's increment vector is determined by finding the socalled saddle point, which is defined as that point on the Newton's direction beyond which the residuals grow; in other words, the Newton's increment vector is not a descent direction beyond the saddle-point. This saddle-point is calculated through various stages, devoted to progressively refining the value of the calculated truncation in order to ensure a good trade-off between robustness and fast convergence. The truncation which produces the saddle-point is assumed to be the most suitable approach for updating the vector of variables as the quadratic convergence features of the NR technique does not deteriorate excessively. Thus, by this approach, the global convergence of the iterative procedure is not ensured to promote a fast convergence.
As commented, the developed 4SA is not competitive with NR in a well-conditioned system. Therefore, the developed algorithm is integrated within a suitable solution framework which, by making proper use of some results obtained at early stages of 4SA, is able to heuristically determine if the studied system is well or ill-conditioned; in the former case, our method proceeds as NR, while 4SA solves illconditioned cases. Thus, the proposed paradigm can efficiently solve both well and ill-conditioned systems, which is not typically encountered on the available robust solvers. This outstanding characteristic makes the developed solution paradigm interesting and suitable for industrial applications, where the system's conditioning is frequently unknown, and a high degree of efficiency in the calculations is typically required.
The computational burden of the developed paradigm is comparable to NR in both well and ill-conditioned cases, making our method scalable to large and even complex-scale networks. This paper is limited to present a simple illustrative version of our paradigm in polar coordinates; however, it is sufficiently versatile to be adapted to different formulations, functional limits and scenarios. Various numerical results are provided in order to validate and show the performance of the developed paradigm.
The remainder of this paper is organized as follows. Due to our proposal is described in this paper using the power flow formulation in polar coordinates, the power flow equations in this form and its solution using NR are described in Section II. The developed 4SA for ill-conditioned cases is introduced in Section III. In Section IV, 4SA is fully integrated within a developed paradigm for effectively solving well and illconditioned cases. Various numerical results on various well and ill-conditioned systems are presented in Section V. The paper is concluded with Section VI.

II. PRELIMINARIES
The power flow equations in polar coordinates are defined as a set of n nonlinear algebraic equations as follows: where: where, P sp i ∈ R and Q sp 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. n l ∈ N is the total number of PQ buses, n g ∈ N is the total number of PV buses, δ g ∈ R n g is the VOLUME 9, 2021 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.
As (1) are nonlinear equations, their solution cannot be obtained directly. The power flow problem has been customarily handled using a plethora of iterative techniques. The most popular is NR method, which is defined by the following map: where, g = ∇ x g ∈ R n×n is the Jacobian matrix and the subindexes denote the iteration counter. Indeed, equation (4a) corresponds with the Newton's increment vector, which is desirable to be a descent direction in order to reduce the residuals as the mapping (4) evolves. Nonetheless, it is well-known that the mapping (4) may be unstable in ill-conditioned systems [7].

III. THE DEVELOPED FOUR STAGE ALGORITHM
The developed 4SA aims to improve the convergence abilities of NR in ill-conditioned systems. Fig. 1 shows a schematic summary of the developed 4SA. Each stage is further explained in the following sections.
where, z j ∈ R + are different truncations, which are collected in the so-called exploration vector z. The exploration vector can be built using Algorithm 1. In this case, the infinity norm of the residuals is assumed to be a good indicator of the proximity of the solution [29].
As observed, the largest truncation is equal to 1, i.e. the largest truncation in fact does not modify φ. For those problems where NR successfully converge, the lowest residual is Algorithm 1 Construction of the Exploration Vector 1: Let ξ k be given 2: Initialize j ← 1 and flag ← 0 is a descent direction on the y-space. However, in ill-conditioned problems, this behaviour is not generally observed and, beyond some determined truncation, the descending trending is lost and the residual grows (see Fig. 2), it means that the Newton's increment vector is not a descent direction over the whole y-space. The furthest point (4a) is still a descent direction for g, called the saddle point. The developed 4SA assumes that the saddle point offers the best trade-off between robustness and convergence rate; consequently, its corresponding truncation should be considered for modifying the Newton's increment vector. The main objective of the Exploration stage is to locate an approximation of the saddle point and its corresponding truncationĥ , which can be achieved by using Algorithm 2. The elements of the exploration vector are calculated in increasing order using the so-called jump parameter ξ , which is given by: As observed, beyond this point the residual grows. In this example, the approximated saddle point is located for z 4 . One should note that the accuracy withŷ is estimated, strongly depends on the value of the jump parameter.
As observed in Fig. 2, the larger ξ , the less accurate estimation of the saddle point. However, this latter aspect is not especially relevant at this stage since only a first approximation to the saddle point aims to be calculated. The Exploration stage usually ends when an approximation of the saddle point Algorithm 2 Exploration Stage 1: Let x k , z k and φ k be given 2: Initialize j ← 1 and saddle ← 0 3: Take y 0 k ← x k 4: while saddle == 0 do 5: if z j k == 1 then 6:ĥ k ← 1 7: break# No saddle point 8: end if 9: Calculate if j == 1 then 12: break # Exploration fails 13: end if 17: end if 18: j ← j + 1 19: end do 20: return approximated saddle truncationĥ k namelyŷ is located. One should observe in Algorithm 2 that the Exploration stage may fail if the saddle point is located at y 0 , i.e. (4a) is not a descent direction. It may also occur that y does not exist, this is typical of well-conditioned cases. In the former case, it is determined that the system is not solvable for the given loading conditions and the algorithm stops prematurely. In the latter case, the largest truncation (i.e.ĥ = 1) is taken.

B. STAGE 2-EXPLOITATION
The objective of the Exploitation stage is to calculate a more accurate approximation of the saddle point. To do that, the vicinity ofĥ is further explored. The exploitation search space is defined by two parameters η 1 < 1 and η 2 > 1 as observed in Fig. 3. The parameters η 1 and η 2 determine the amplitude of the search space. The search space should not be too wide in order to facilitate its efficient exploration. Indeed, the search space is in practice split into N intervals (see Fig. 4). So that, the wider search space, the higher N in order to ensure an accurate estimation of the saddle point  (h ). This fact is relevant since this stage has a computational burden proportional to N . Hence, if N has to be high, the computational cost of this stage would also be high, affecting the performance of 4SA. Nevertheless, the parameter N might be selected high at first iterations, when the solution is presumably still far and the algorithm may potentially fail, and progressively decrease it as the algorithm satisfactory evolves (see Section IV). It is found that η 1 ∼ 0.5 and η 2 ∼ 1.2 work satisfactorily in most cases.
The entire process for defining the Exploitation search space is shown in Algorithm 3 (where the operator (a : : b) builds a vector whose first element is a, the last element is b, and the step is ). As a result, the exploitation vector is built and collects N different truncations j (see Fig. 4).

Algorithm 3 Construction of the Exploitation
The Exploitation stage is carried out as follows. Up to N elements as (7) are calculated using the different truncations collected in the exploitation vector. Then, the accurate saddle truncationh is determined as that element of which produces the vector (7) with the lowest residual. This process is summarized in Algorithm 4 using pseudocode.

C. STAGE 3-PONDERATION
As a result of the Exploitation stage, a second approximation to the saddle truncationh is obtained, which is presumably more accurate thanĥ . However, we cannot assert that this truncation exactly yields the saddle point since it is only VOLUME 9, 2021 obtained with very accurate Exploration and Exploitation stages (ξ ↓↓ or N ↑↑), which is unaffordable in practice. Therefore, the truncation obtained at Stage 2 should not be considered for modifying the state vector (2d) and directly updating it. This stage aims to determine which truncation is more secure for updating the state vector. If one considers all truncations up to the exact saddle point as safe and all evolution beyond it as unsafe, we can divide the exploitation vector into two zones. However, as commented, the exact saddle point is unknown. Therefore, it is not suitable to reduce all the deductions to only two zones; instead, consider three ones. Thus, a caution space is added between the safe and unsafe zones. This new zone lies aroundh , and it can be defined as these points that one cannot assert if the evolution is safe or not. As commented, each truncation of can be categorized as safe, caution or unsafe. Nevertheless, one should also put into value how each truncation affects the convergence rate of the developed 4SA. Thus, the lower truncations are associated with a low convergence rate, while the higher ones normally produce a faster evolution of the state vector (close to quadratic convergence). A priori, one shall consider the fastest evolution of the state vector for converging employing less iterations; however, as observed in Fig. 5, those elements of the exploitation vector with higher convergence rate usually lie on the unsafe zone. As shown in this figure, the exploitation vector can also be divided as a function of the speed of convergence; thus, the lower truncations are associated with a low convergence rate while the higher ones normally produce a faster evolution of the state vector. To sum up, the different j can be categorized attending to different criteria. The main issue is that some ponderations offer good and bad properties at the same time. To solve this problem, the Ponderation stage aims at assigning a weight to each element of following the rules collected in Table 1. Although these rules are based on heuristic foundations and, consequently, they are not universal, the adopted approaches are typically well-behaved and work quite well, as demonstrated in Section V. As observed, the assigned weights aim to balance between a safe and fast evolution of the state vector. Following the guidelines of Table 1, the different weights are assigned using the following function: (8) where, h j is the j th element of the so-called ponderation vector h, σ ∈ R + is the standard deviation, e is the natural exponential and µ = η µh with η µ slightly lower than 1. The function (8) is used in this work since it naturally weights the different elements of the exploitation vector according to the rules presented in Table 1. Also, the function (8) compensates the heuristic foundations adopted in Table 1; hence those elements of the exploitation vector with presumably bad characteristics (for instance, points within the unsafe zone), are not nullified. Instead, the function (8) automatically assigns them a low weight to reduce their influence on updating the state vector. This approach also compensates for the intrinsic not-totally accurate estimation of the saddle point after carrying out the Exploration and Exploitation stages. Thus, even the impact on the selection of some arbitrary parameters (η 1 , η 2 , N and η µ ) is reduced in order to facilitate the implementation of the developed 4SA.
The main steps of the Ponderation stage, along the process for building the ponderation vector, are summarized in Algorithm 5 using pseudocode. Fig. 6 plots the weights assigned by (8) as a function of h for different values of the parameter σ . As observed, the different weights are assigned according to the criteria collected in Table 1. It is worth noting that the highest weight is assigned to µ, which is a few lower thanh . Thus, the highest weight is always assigned to µ. As commented, η µ should be slightly lower than 1. Indeed, since the saddle point is never exactly determined, it is preferable to assign the highest weights to those points lower thanh , since these values are presumably close to the safe zone and, consequently, algorithm's robustness is more guaranteed. The different zones are determined based on heuristic criteria; thus, the lower truncations lie in the safe, low convergence zone, while most of the higher truncations are in the unsafe high convergence area. It is worth noting that it is unsuitable to assign high weights toh since, as commented, one cannot know if this truncation has lied in the safe or unsafe zone. It is also remarkable how the parameter σ affects to (8). Low values of σ provoke that (8) gives a very high weight to µ, while the other truncations of the exploitation vector are nullify. In other words, if σ is fixed very low, it is assumed that µ is the most suitable truncation for the Newton's increment vector, so this approach is not generally recommended since, as commented, the saddle point is never determined with total accuracy. Instead, it is preferable to assign a high value of σ . This approach assumes that the exploitation search space is not too wide, allowing a correct exploration of the Exploitation vector. In practice, it is observed that σ ∼0.5 − 1 works satisfactorily in most cases.

D. STAGE 4-UPDATING
Finally, the state vector (2d) is updated, taking into account the weights calculated at the previous stage. Thus, the considered most suitable truncation is given by: And the state vector is updated as: The main steps of Stage 4 are summarized in Algorithm 6 using pseudocode.

Algorithm 6 Updating Stage
1: Let x k , φ k , w k and h k be given 2: Calculate The main calculations of the developed 4SA are summarized in Table 2. The new proposal only requires a Jacobian evaluation and a LU factorization. In that sense, its computational burden is similar to NR. However, 4SA needs multiple function evaluations in the stages 1 and 2. Although this computation is generally cheap, it undoubtedly supposes an extra computational burden with respect to NR. In the following section, we propose a paradigm for avoiding these calculations.

IV. IMPLEMENTATION OF THE DEVELOPED ALGORITHM
As commented in the previous section, the developed 4SA is inevitably less efficient than NR. This disadvantage makes our proposal less attractive for solving well-conditioned systems, where NR undoubtedly offers better performance. To tackle this issue, the developed 4SA is integrated within a novel solution paradigm. This novel solution framework avoids the usage of 4SA in those cases easily solvable with NR. Thus, the developed 4SA is only used in ill-conditioned systems or when the convergence does not seem ensured. Now, the main question arises: how is it determined if 4SA is suitable or not? To answer this question, one can take advantage of the information brought by 4SA for employing the following rules: • The developed 4SA can heuristically determine if the case is well or ill-conditioned using the jump parameter. It is worth noting that this parameter is, in fact, inversely proportional to the infinity norm of the Newton's increment vector (see eq. (6)). When a nonlinear system is ill-conditioned, some elements of the Newton's vector tend to be abnormally large due to the ill-conditioning of the Jacobian matrix. Therefore, it may be reasonable to assume that the degree of ill-conditioning of a system is inversely proportional to the value of ξ . Of course, the well-known condition number of the Jacobian matrix may also be used for this purpose [28]; however, the threshold considered for determining if a nonlinear system is ill-conditioned based on the condition number is strongly case-dependent. Indeed, it is observed that the condition number not only depends on the conditioning of the Jacobian matrix, and it is also affected by other characteristics; for instance, largescale cases frequently present large condition numbers even if they are not ill-conditioned [28]. In addition, calculating the condition number may be computationally expensive. For these reasons, the jump parameter is preferred to estimate the degree of ill-conditioning of the problem and act accordingly. Thus, at the beginning of the algorithm, it can be determined if NR is suitable for solving the studied case or not. Based on our experience, most well-conditioned systems yield ξ 0 > 0.5; hence this threshold can be considered for determining the usage of NR or 4SA.
• As 4SA evolves, the size of the exploitation and ponderation vectors can be reduced using the following rule: (11) where, round (·) rounds to the nearest integer. As observed, the rule (11) is only activated if the residual is reduced at k + 1, so that the computational burden is reduced as less function evaluations are required during the process. When N = 1, the Exploration and Ponderation stages are immaterial and, consequently, NR can be used instead.
• If the result of both the Exploration and the Exploitation stages is equal to 1, one can guess that the Newton's evolution vector is safe and, consequently, NR can be reliably employed. In conclusion, the developed paradigm consists of carrying out 4SA and switching to NR when some of the above conditions are matched. Furthermore, the flowchart of the developed power flow solution paradigm is shown in Fig. 7.
It is worth mentioning that a more straightforward approach may consist of simply running NR and 4SA simultaneously and, after convergence, take the more reliable solution. However, this procedure is not functional at all and the developed algorithm would outperform this procedure. For example, in the case of ill-conditioned systems, NR just diverges so it is unnecessary to run this algorithm. In addition, simultaneously running both methods suppose an important computational burden that can be avoided by implementing the developed solution framework. On the other hand, it is difficult to know a priori if a power system is naturally ill or wellconditioned. In this sense, it is necessary an indicator that determines the degree of ill-conditioning of a system without needing human assistance. Thus, the algorithm automatically determines which method is more suitable for solving the network. In the developed framework, this issue is addressed by the parameter ξ .

V. NUMERICAL EXPERIMENTS
The developed 4SA and the solution paradigm, described in Section IV and summarized in the flowchart in Fig. 7, have been tested along NR, the conventional Backtracking line search technique (BT) (e.g. see [44]), the High order Levenberg-Marquardt method (HLM) [30], the Backward-Euler technique (BE) [22], and the Reverse Bulirsch-Stoer power flow solution method (RBS) [27]. For HLM, the socalled damping factor is updated each iteration by g 1.3 ∞ . Regarding BE, the implementation denoted by ICNM-J1 [22] has been tested. In the case of RBS, the guidelines provided in [27] have been followed. For our methods, N 0 = 20, η 1 = 0.50, η 2 = 1.20, η µ = 0.75 and σ = 1 have been taken. All tested techniques have been implemented within the Matpower environment [45]. Matpower provides an extensive Matlab-based library for power system analysis. This library fully exploits the Matlab capabilities for handling large sparse systems as those typically encountered in power system analysis. Hence, factorizations and linear systems solutions are carried out using sparse facilities such as the minimum degree ordering algorithm or partial pivoting. In addition, all the tested cases are available in Matpower database [46]- [48]. They have been denoted as they are labelled in this software. Two cases are contemplated; thus, the studied cases are solved from a flat start or the default starting point provided in Matpower.
All reported times have been obtained on a 64-bit i5-9400F Intel Core personal computer (2.90 GHz, 8 GB of RAM), as the average value of 100 simulations; on the other hand, g ∞ < 10 −6 has been taken as a convergence criterion. Table 3 provides the results obtained for base cases. As seen, both 4SA and the developed paradigm successfully solved all studied systems. Both techniques were well-behaved generally for different cases. However, some methods outperformed 4SA in a few cases. This issue is solved by the developed solution paradigm, which offered the best trade-off between robustness and efficiency. The good convergence properties observed when using the developed techniques, are undoubtedly due to the accurate estimation of the most suitable truncation of the Newton's increment vector after carrying out the Exploration, Exploitation and Ponderation stages. Thus, the developed algorithms do not excessively deteriorate the quadratic convergence features of NR, as other algorithms do. For instance, RBS and BE truncate the Newton's increment vector using a so-called step size, which follows an algorithm-based update rule rather than observing the residuals' behaviour as 4SA and the developed solution paradigm do. Therefore, the quadratic convergence behaviour of NR is usually preserved and the solution is typically reached in a reasonable number of iterations.

A. GENERAL COMPARISON FOR BASE CASES
It is worth observing that the developed paradigm offered the same results compared with NR in well-conditioned cases. This result is logic since as commented in Section IV and described in Fig. 7, the developed solution algorithm proceeds as NR in well-conditioned systems. On the other hand, the developed solution framework is able to detect ill-conditioned systems and perform as 4SA in such cases until the convergence is ensured, then, the algorithm switches to NR to speed up the convergence.
The conventional NR showed good performance with the default initialization; however, it failed in some cases when a flat start was used. When converged, NR was along the developed solution paradigm, the most competitive technique. This is due to the proposed solution framework was able to identify ill-conditioned cases accurately. According to the definition given in the Introduction and [7], those systems where NR failed from a flat start can be categorized as ill-conditioned. Such is the case of 'case3012wp', 'case3375wp' and 'case13659pegase', as seem in Table 3. As shown in Table 4, the value of ξ 0 allows to discriminate between well and ill-conditioned cases; thereby, the value of this variable is quite short in these three cases. Thus, the developed solution paradigm performed as NR in those wellconditioned cases, saving iterations and solution time. It is also worth noting that the value of ξ 0 abruptly varies from well to ill-conditioned cases; hence, the threshold adopted in this paper (ξ 0 > 0.5), which was taken based on personal experience, has not to be defined finely (note that one could take ξ 0 > 0.8 instead and the result would be the same). In other words, this threshold has not to be selected carefully; actually, it could take different values, and the degree of illconditioning would still be correctly estimated. In the light of these results, it is demonstrated that the value of ξ 0 supposes a good indicator to determine the degree of ill-conditioning of a system; in addition, unlike the condition number, the jump parameter seems case-independent, which enables its broad utilization.
The considered line search method (BT), successfully solved all the studied cases. However, one can see a clear deterioration of the convergence properties in some cases with respect to NR. Typically, this technique required a large number of iterations to find the solution; this is especially remarkable when a flat start is used, which makes this technique few practical in ill-conditioned cases.
HLM was quite robust. However, its efficiency was deteriorated as the size of the system grows. This is due to it requires a matrix-matrix product for each iteration. This operation is much costly in large-scale systems. In addition, it shows a low convergence rate when the starting point lies far away from the actual solution. Thus, it generally required many iterations for converging from a flat start.
RBS is a reliable method but inefficient due to the enormous amounts of LU factorizations required. Following the guidelines presented in [27], this technique requires six factorizations at the first iteration, four at the second iteration and 2 for the remaining iterations. In addition, linear convergence is observed during the whole iterative procedure due to the influence of the step, which barely reaches its optimal value for the last iterations.
BE under its implementation, namely ICNM-J1 in [22], is quite efficient as only a LU factorization is required, and the calculation of the Hessian matrix is avoided. As RBS, this technique frequently shows a very linear convergence pattern. In addition, it occasionally converged to the low voltage solution while the developed techniques always converged to the stable operating point. Figure 8 provides a better overview of the convergence properties of the different studied solvers. More precisely, this figure plots the convergence profiles for the 'case3012wp' with flat initialization. As seen in this figure, while NR shows an oscillatory pattern which indicates divergence, remainder solvers successfully converged. The developed solution paradigm and 4SA performs similarly during the first iterations. This is due to the developed paradigm is able to identify this case as ill-conditioned and, consequently, 4SA is run at the beginning of the iterative procedure. After a few iterations, the developed solution framework determined that the convergence is ensured and switched to the conventional NR, thus speeding up the convergence.   in all cases, only the latter has been reported in this figure). Here, the loading level has been modified according to:

B. INFLUENCE OF THE LOADING LEVEL
where λ ∈ R + is the loading level. As seen from Fig. 8, the performance of HLM gets worse as the loading level grows, employing too many iterations close to the maximum loadability point. BT, RBS and BE are barely affected by the loading level and their convergence rates stay almost constant. However, the developed methods were the most competitive, employing fewer iterations than the remainder techniques.

C. VALIDATION WITH REACTIVE LIMITS ENABLED AT PV BUSES
Although the developed formulation of 4SA does not take into account the equipment limits, they can be easily incorporated using any mechanism available for NR. Let us here consider that switch logic strategy described in [49], for incorporating the generators' reactive power limits within the power flow calculation using NR. This strategy basically consists of converting those PV buses with any limit violation to PQ buses, taking the hit limit as injected reactive power. In these converted buses, the injected reactive power is now taken as an independent variable. The nodal voltage freely varies; then, starting from the calculated solution without reactive violations, the process is repeatedly carried out until a feasible solution is reached. Otherwise, the case is declared unfeasible if a feasible solution has not been found with all the PV buses converted to PQ. It is worth noting that this strategy is quite universal, and other references and software have adopted similar approaches (e.g. see [7], [45] or [50]). Nevertheless, the developed solution algorithms are versatile enough to incorporate other more sophisticated approaches such as those described in [1] and [51]. The aim of this section is to simply validate the developed power flow solution techniques when some functional limits are taken into account; hence the strategy described in [49] and above described has been implemented within 4SA and the developed paradigm. Fig. 10 shows the nodal voltage angles and magnitudes for the 'case30' with reactive limits obtained by 4SA and NR from a flat start (since the developed paradigm performs in this case as NR, their results are not shown). As observed, 4SA obtained very accurate results. Fig. 11 illustrates the implemented bus switching mechanism. As shown, the injected reactive power and voltage magnitude at bus #2 during execution of 4SA. As observed, the reactive power grows to keep the nodal voltage magnitude constant until its limit is reached. After that, this bus is converted to PQ type. Then, the reactive power is kept constant, and the nodal voltage magnitude varies until a feasible solution is reached. It is worth noting that this same strategy can be considered for taking into account other controls or functional limits like transformers tap changers.
It is well known that incorporating an outer switching mechanism to consider the equipment limits may deteriorate the convergence properties of a PF solver; this is because more PQ buses and, consequently, variables and equations are incorporated into the system each iteration of the outer loop. This is illustrated in Fig. 12, where  the convergence profiles of NR, 4SA and the developed paradigm for the 'case3012wp' from a flat start are plotted. As observed, NR failed to solve the base case because it oscillates during a large number of iterations. In this case, the developed paradigm and 4SA performed well in case of considering the generators reactive limits. As seen, the proposed methods successfully follow the switching mechanism to a feasible solution after the solution of 3 power flow problems. VOLUME 9, 2021 In order to fully show how the developed techniques perform in the presence of reactive limits, we have compared them with the other considered solvers. In all cases, the switching strategy previously described has been incorporated within the different studied PF techniques. Table 5 compares the different studied solvers to consider the generators' reactive limits. As seen for the base cases, the developed techniques were well-behaved and showed very competitive results.

VI. CONCLUSION AND FUTURE WORKS
In this paper, a novel power solution method has been proposed to solve well and ill-conditioned cases. It is carried out through four stages, devoted to determining the optimal evolution of the Newton's increment vector. Although the developed method is more efficient than most available robust solvers, the developed 4SA is still not competitive with NR in well-conditioned cases. This disadvantage would limit the potential application of the developed technique in commercial packages. To overcome this drawback, 4SA has been integrated within an efficient paradigm, which can guess if the studied case is well or ill-conditioned. In the former case, NR is used while 4SA is employed for solving the illconditioned cases. The developed paradigm also switches to NR when the convergence seems ensured.
Various numerical results have been provided in a variety of cases and scenarios. The results showed that 4SA is more competitive than most robust, tested techniques; however, it is still less efficient than NR in most cases. The integration of 4SA within the developed framework brought superior characteristics. Thus, the developed solution paradigm was generally the most competitive method, offering the best trade-off between robustness and efficiency.
Despite the promising results obtained, the developed solution framework still has some limitations as it fails at turning points and unsolvable cases. Future works should be devoted to expanding the applicability of our proposal to these cases. TALAL ALHARBI (Member, IEEE) received the B.Sc. degree in electrical engineering from Qassim University, Saudi Arabia, in 2010, and the M.A.Sc. and Ph.D. degrees in electrical and computer engineering from the University of Waterloo, ON, Canada, in 2015 and 2020, respectively. In 2020, he rejoined the Department of Electrical Engineering, Qassim University, where he is currently an Assistant Professor. His research interests include electric vehicle, power system operation and optimization, energy management, energy storage systems, power electronic and artificial intelligence applications in power systems, and smart grids.
OMAR ALRUMAYH (Member, IEEE) received the B.Sc. degree in electrical engineering from Qassim University, Saudi Arabia, in 2011, and the M.A.Sc. and Ph.D. degrees in electrical and computer engineering from the University of Waterloo, ON, Canada, in 2016 and 2021, respectively. In 2021, he joined the Department of Electrical Engineering, Qassim University, where he is currently an Assistant Professor. His research interests include power system operations, energy management, energy storage systems, and smart grids.
SALAH KAMEL received the dual Ph.D. degree (international) from the University of Jaen, Spain (Main), and Aalborg University, Denmark (Host), in January 2014. He is currently an Associate Professor with the Electrical Engineering Department, Aswan University. He is also a Leader of the Power Systems Research Group, Advanced Power Systems Research Laboratory (APSR), Aswan, Egypt. His research interests include power system analysis and optimization, smart grid, and renewable energy systems.
FRANCISCO JURADO (Senior Member, IEEE) was born in Linares, Jaén, Spain. He received the M.Sc. and Dr.Ing. degrees from the National University of Distance Education, Madrid, Spain, in 1995 and 1999, respectively. Since 1985, he has been a Professor with the Department of Electrical Engineering, University of Jaén, Jaén. His current research interests include power systems, modeling, and renewable energy.