A Three-Stage Algorithm Based on a Semi-Implicit Approach for Solving the Power-Flow in Realistic Large-Scale ill-Conditioned Systems

Solving the Power-Flow in realistic large-scale ill-conditioned systems supposes a challenging task for most of available solution methodologies. This paper tackles this issue by developing a novel efficient and robust Power-Flow method. It is mainly based on a Semi-Implicit approach but incorporates other numerical arrangements for enhancing its features. The resulting three-stage algorithm is validated using several realistic ill-conditioned systems ranging from 3012 to 70000-buses. Results show that the developed methodology constitutes an efficient and robust Power-Flow solution technique, outperforming the results obtained with other available approaches.


I. INTRODUCTION A. MOTIVATION
The Power-Flow (PF) problem has been customary solved using the Newton-Raphson technique (NR) [1]. This method is quite efficient and offers acceptable robustness properties for managing with the so-called well-conditioned systems. Oppositely, NR finds multiple convergence issues in ill-conditioned systems [2], [3].
A huge variety of robust methodologies have been developed in order to properly solve ill-conditioned power systems. However, these techniques are typically less efficient than NR, which limit its application in realistic large-scale systems.
Nowadays, ill-conditioned systems are becoming more frequent [4]. Consequently, the usage of robust methodologies in industrial tools is turning to be primordial. In that sense, developing robust PF solvers being able to solve realistic The associate editor coordinating the review of this manuscript and approving it for publication was Feiqi Deng . large-scale ill-conditioned systems efficiently, is currently an important topic.

B. LITERATURE REVIEW
The Iwamoto's method was proposed early in the 1980's [5]. It basically consists on calculating an optimal multiplier each iteration, which is devoted on modifying the increment vector for ensuring the convergence. The Iwamoto's method supposes a very popular robust PF solver and it is frequently adopted as benchmark methodology (see e.g. [6], [7]). It is very robust, however, it usually employs many iterations for converging due to the optimal multiplier tends to be very short as the solution is approached [8]. In [9], formulation of the Iwamoto's method was extended to polar coordinates.
The Continuation Power-Flow [10] has been profusely studied for voltage stability analysis. This tool prevents the singularity of the Jacobian matrix at turning points by using a parameterized formulation of the PF equations. Basically, the Continuation Power-Flow solves multiple PF problems for different loading levels in order to obtain the stability margin of the system. Due to numerous PF problems have VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see http://creativecommons.org/licenses/by/4.0/ to be solved, this tool is considered inefficient in comparison with NR when a PF problem is aimed to be solved [11]. Based on the Continuous Newton's method [12], an analogy between the Newton's method and the Forward-Euler methodology can be established. This idea was extended to the PF problem in [2], establishing a formal common framework between most of robust PF solvers and the Forward-Euler technique. On the basis of this analogy, any numerical integration method can be considered for developing robust PF solution techniques. The authors have exploited this idea in several references [13]- [15]. In addition, the Implicit counterpart of the Continuous Newton's method has been proposed for PF analysis in [16]. Although more robust than NR, these methodologies are still sensitive to badly initialized problems. Moreover, they might be very inefficient since several matrix factorizations have to be computed each iteration.
Using the Lyapunov theory, the PF problem can be posed as an artificial dynamic system [4], [17]. The resulting paradigm can be solved using integration routines like Matlab ode solvers. Although this solution framework is barely affected by the starting point, integration methods employed in its resolution are computationally inefficient in large-scale systems [4], which limit its application in realistic systems.
PF can be raised as an optimization problem. Consequently, it can be solved using minimization techniques. Applicability of the Levenberg's method for PF solution was explored in [18]. Theoretically, this technique should ensure a solution, at least, in the sense of sum of squares of residuals. However, solution of the established optimization problem does not always correspond with the actual solution of the PF equations [18]. In addition, convergence features of the Levenberg's method may lost if the starting point lie far away to the solution, resulting in an inefficient mapping. Finally, it involves denser factorizations due to the Jacobian matrix is multiplied by its transpose. An attempt for overcoming these drawbacks was made in [6], by proposing a high order Levenberg-like technique for PF analysis.

C. CONTRIBUTIONS
From the Literature Review is deduced that most of available robust PF solvers are, in fact, totally inefficient to be applied in realistic cases. This paper tackles this topic by developing a robust PF solution method suitable for large-scale ill-conditioned systems. This is achieved by proposing a robust but efficient yet algorithm based on the Semi-Implicit approach (SIA) developed in [19]. The resulting iterative scheme achieves outstanding features by combining several numerical arrangements: • Lavrentiev's regularization [20], for enhancing the robustness characteristics at turning points.
• Chebyshev-like method with cubic convergence [21], for speeding up the convergence features.
• Heun's method [22], for obtaining a good balance between robustness and efficiency.
By combining the aforementioned numerical techniques, a three-stage algorithm is developed (3S-SIA). The developed methodology is validated using various realistic large-scale systems under different demand conditions like heavy load, cascading failures, high R/X ratio or equipment limits enforcement. Proposed algorithm is compared with other well-known standard and robust PF solution approaches.

D. PAPER ORGANIZATION
Remainder of this paper is organized as follows. Section II briefly explains SIA. Section III outlines the PF problem and ill-conditioned power systems. Developed PF solution algorithm is explained in Section IV. Several numerical experiments are carried out in Section V. Finally, the paper is concluded with Section VI.

II. BACKGROUND
SIA was proposed in [19] for solving systems of nonlinear equations. Developed 3S-SIA is based on this methodology, therefore, a brief explanation about SIA is here included for the sake of self-sufficiency.

A. SIA FOR SOLVING NONLINEAR EQUATIONS
Let us begin with the simplest 1-dimensional case given by: where, y ∈ R, f : R → R and ϕ : R → R is a primitive iterative method as follows: where the superscript denotes the k th iteration of the iterative solution. Eq. (2) can be modified to become a new expression with the same roots as: where, ρ ∈ R. The new degree of freedom introduced by parameter ρ in (8) opens the opportunity to optimize the convergence rate.

B. SIA FOR SOLVING SYSTEMS OF NONLINEAR EQUATIONS
Now, let us consider the more general n-dimensional case.
To do that, equation (3) is extended to its multivariable counterpart as follows: where y ∈ R n , I ∈ R n×n is the identity matrix, A ∈ R n×n and ψ : R n −→ R n is defined as follows: where, ϕ : R n −→ R n is an old iteration method. In this case, it is given by: where, f : R n −→ R n . Elements of matrix A are determined by differentiating (5): where f = ∇ T y f (y) ∈ R n×n is the Jacobian matrix formed by the first partial derivatives of f with respect y and: Elements of the matrix R are taken as free parameters which should be adapted each iteration in order to obtain quasi-monotonous and superlinear convergence. One should note that the matrix A is only defined if f is invertible. Grouping (5)-(8), a generic k th iteration of SIA for solving systems of nonlinear equations may be written as follows:

III. PF PROBLEM AND ILL-CONDITIONED SYSTEMS
Let us consider the PF in polar coordinates as a set of nonlinear equations as follows: where g :R n −→ R n are the PF equations and x ∈ R n is the PF state vector. Solution of (1) using NR is carried out as follows: where g = ∇ T x g (x) ∈ R n×n is the Jacobian matrix of PF equations. As it was commented, PF cases can be broadly categorized as well or ill-conditioned. A widely employed measure for determining the degree of ill-conditioning of a power system is the condition number defined by Neumann and Goldstine [23]. One can affirm that, if the condition number grows, the degree of ill conditioning grows as well.
However, an ill-conditioning system can be still solved using standard techniques. In fact, instability experienced by the map produced by NR is due to a poor starting point rather than ill conditioning of equations. Following this idea, we consider a system ill-conditioned if, despite its solution does exists, it is not reachable using NR and a flat start [2].

IV. DEVELOPED PF SOLVER A. FOUNDATIONS
As it was commented, mapping defined by (9) is only stable if the Jacobian matrix is continuously invertible. This is assumed in PF solution by NR, however, the Jacobian matrix may become near singular close to a turning point. In PF analysis, this situation is not infrequent since the maximum loadability point of a system is a well-documented turning point [10]. Regularization schemes are very popular methodologies for improving the stability of a nonlinear solver in the vicinity of a turning point. Basically, regularization consists on adding some elements to the Jacobian matrix in order to avoid its singularity.
Most of available robust PF solvers are very inefficient, at least, in comparison with NR. A reason of this low degree of efficiency is their linear convergence rate, which is reflected in a huge amount of iterations for achieving a solution (one should note that it is the main disadvantage of the Iwamoto's method). This low convergence rate is because some robust methodologies are, in fact, based on NR, which has quadratic convergence. Alternatively, other nonlinear solvers with high convergence rate can be considered. However, this high-order methodology should be efficient. In that sense, those techniques that require more than a matrix factorization per iteration should be avoided. This motivates us to use the Chebyshev-like method with cubic convergence defined in [21], which only require a matrix factorization each iteration, resulting in a similar computational burden in comparison with NR.
It is worth mentioning that cubic methodologies are typically weaker than NR [24]. In order to preserve an acceptable robustness but maintaining the resulting algorithm efficient enough, it is compulsory to establish some mechanism for balancing efficient and robust properties of the considered numerical arrangements. These characteristics were found on the Heun's method [22], which presents a particular structure which allows to balance the characteristics of two iterative procedures. As result of combining the aforementioned numerical arrangements, a 3-stage algorithm is obtained. Thus, at first step, the Lavrentiev's regularization is introduced while the Chebyshev-like [21] and Heun's methods are respectively used at second and third stage. Finally, authors have found in SIA an estimable framework for gathering all mentioned methodologies in an elegant way. Conceptual idea of the developed 3S-SIA is sketched in Fig. 1.
B. 3S-SIA FORMULATION SIA described in Section II is widely convergent, however, its numerical stability is not always guaranteed. Robustness of SIA can be notably enhanced using regularization. Basically, two regularization schemes namely Tikhonov [25] and Lavrentiev [20], are available in the literature. While the former is more stable and produces better condition numbers [26], the latter is simpler and computationally cheaper. For those reasons, we have preferred to consider the Lavrentiev's regularization mechanism. However, one should know that the Tikhonov's regularization could be also introduced in the same manner. In order to properly employ the Lavrentiev's regularization mechanism within SIA procedure, let us define the following increment vector: where α ∈ R + is the so-called damping factor which also appears in the Levenberg-like techniques [18] and h :R n −→ R n is given by: where H : R n −→ R n is defined by: where (k) = g x (k) . In this case, (12) is used for calculating an intermediate value of the PF state vector as follows: where R was already defined in (8). At this point, (15) can be used for calculating other increment vector following, in this case, the same structure of the cubic Chebyshev-like method [21]. Thus, this increment vector is calculating as follows: One can check that (12) represents a robust evolution of x in the sense of it uses the Lavrentiev's regularization. On the other hand, (16) is devoted to speedup the convergence rate in the sense of it shares structure with the cubic methodology [21]. Higher convergence rate is in this case provided by evaluating (10) in two points namely x (k) and x (k) . This is the basis of the multipoint iterative methods [27], which achieve higher convergence rate without extra matrix factorizations. Since (12) and (16) bring robust and efficient properties, respectively, it would be interesting provide some mechanism for balancing the influence of each one. Thus, the PF state vector is finally updated using the Heun's method as follows: By (17), x (k+1) is influenced by both (12) and (16).

C. PARAMETERS TUNING
3S-SIA have two degrees of freedom namely R and α.
For promoting an easy implementation of the developed PF solution paradigm, these parameters should be updated each iteration using adaptive mechanisms. In this case, SIA provides easy updating schemes as it will be showed during this section. Firstly, matrix R can be initialized as follows: where δ is the Kroenecker delta and ρ 0 ∈ (0, 1). A logic strategy is fixed ρ 0 close to 1 and reduces it as solution is successfully approached [19]. Hence, ∞ seems a good indicator for determining if the PF is being correctly solved. Thus, R is updated for k ≥ 1 as follows: It is worth mentioning that R evolves according to the value of α, this is because both parameters should follow similar patterns. Let us observe equation (12), one can easily check that, a priori, this increment vector does not lead to the actual solution of PF equations (10) for α = 0 [18]. However, if α is fixed too short, robustness properties of the Lavrentiev's scheme are lost due to elements added to the diagonal of the Jacobian matrix are not large enough to ensure its nonsingularity. Thus, it is easily deduced that α, as elements of R, should be large at first iterations while they should be progressively reduce as the iterative procedure progresses (this behaviour is also encountered in other updating mechanisms used in some Levenberg-like methods [28], [29]). Following this guideline, we propose updating α as follows: To complete the section, developed PF solution procedure based on 3S-SIA is summarized in Algorithm 1 using pseudo code. Here, ∞ has been considered as convergence criterion. Algorithm also stops if the iteration counter surpasses a predefined threshold (k max ), this situation normally indicates divergence. Solve (12)

D. COMPUTATIONAL BURDEN
Developed PF solver described in Algorithm 1, only requires a matrix factorization each iteration. If we assume the factorization of the Jacobian matrix as the heaviest computational part of a PF calculation [2], we can affirm that the developed PF solution procedure is quite efficient and its computational burden is comparable to NR.

E. HANDLING EQUIPMENT LIMITS
Developed formulation does not directly take into account the equipment limits. However, they can be easily incorporated. For example, the generators' reactive limits can be considered using the subroutine described by the following steps: Step 1: Solve the PF problem using the developed method. If there are not PV buses and PF has not found a solution, the problem is declared unfeasible.
Step 2: Check if any reactive limit has been violated.
Step 3: Create a Q + set with these generators' indices whose higher limit is violated.
Step 4: Create a Q − set with these generators' indices whose lower limit is violated.
Step 6: For each j ∈ Q + set Q j = Q max j , where Q j is the injected reactive power at bus j and Q max j denotes the upper bound of the injected reactive power at bus j.
Step 7: For each j ∈ Q − set Q j = Q min j , where Q min j denotes the lower bound of the injected reactive power at bus j.
Step 8: Set x (0) = x (k) and return to Step 1. Although the indicated steps define a simple mechanism for incorporating the generators' reactive limits, one can easily check that other equipment control or limits might be easily incorporated using similar strategies.

F. ROBUSTNESS OF 3S-SIA
The developed 3S-SIA is mainly devoted on solving ill-conditioned systems. As it was commented, this kind of systems are, in fact, badly initialized. Consequently, the developed algorithm is expected to be, at least, more robust than NR with respect the initialization of the iterative procedure. In order to validate this feature, we have carried out a statistical analysis which a set of starting points is given by: where x * ∈ R n denotes the solution of the system and N (x * , σ ) is a Gaussian distribution with mean x * and standard deviation σ . Thus, different starting points are built for different values of σ , then, these initial guesses are used for solving the PF problem in the IEEE 30-bus test system [30]. Fig. 2 shows the total number of PF successfully solved by NR and 3S-SIA for different standard deviations. Clearly, the developed 3S-SIA is able to solve more scenarios than NR. Therefore, we can deduce that the developed algorithm is more robust than NR.

V. NUMERICAL EXPERIMENTS
Several numerical experiments are carried out in order to validate the developed 3S-SIA. In that sense, the following PF solvers are considered for comparison: • Standard NR in polar coordinates. • Iwamoto's method in polar coordinates [9]. • High-order Levenberg-Marquardt's method (HOLM) [6]. In this case, the damping factor is updated according to 1.3 ∞ • 4 th order Runge-Kutta PF solver (RK4) developed in [2]. • Reverse-Bulirsch-Stoer (RBS) defined in [14]. • Implicit solution of the Continuous Newton's method based on the Backward-Euler method (BEM). On the basis of the results obtained in [16], we have considered the implementation namely ICNM-J1 in this reference.
• Developed 3S-SIA algorithm. Parameters involved in HOLM, RK4, RBS and BEM have been adjusted according to the guidelines provided in their respective references. Regarding solution procedure based on 3S-SIA, ρ 0 = 0.9 has been taken.
All considered methodologies have been coded in Matpower [31]. On the other hand, all simulations have been run on a 3.4 GHz Intel Core i7-8750H CPU 2.2 GHz personal laptop (16.00 GB RAM), under Windows 10 and Matlab R2014b. In addition, one of the following situations may be reported in simulations: • Convergence: in this case, a PF solver is able to achieve the correct solution within a predefined convergence tolerance before surpassing k max . In that sense, solution found by NR using the starting point provided in Matpower is considered the correct one, therefore, if the solution found does not differ from the correct solution greater than 0.01, this is considered accurate enough. On the other hand, ε = 10 −6 and k max = 2000 have been taken in simulations.
• Fail: in this case, residual grows or oscillates during a considerable number of iterations. Thus, a PF solver is not able to achieve a solution within the required convergence tolerance before surpassing k max .
• Low voltage solution (labelled † ): due to the quadratic form of PF equations, they have two solutions namely high and low voltage solutions. The former has been considered the correct one, however, we have also reported when a PF solver converged to the low voltage one.
• Inaccurate solution (labelled * ): in this case, a PF solver successfully converged, however, the achieved solution differs from the correct one greater than 0.01. All iterative solutions have been initialized using a flat start (e.g. all voltage magnitudes equal to 1 pu at PQ buses and all voltage angles equal to 0 deg). Finally, it is worth mentioning that reported execution times have been obtained as the average value of 100 simulations.

A. STUDIED SYSTEMS
Four realistic large-scale ill-conditioned systems have been considered: • System #1: the 3012-bus snapshot of the Polish transmission system at 2007-08 winter evening peak [32].
• System #4: the 70000-bus synthetic grid [30], [35]. All considered systems are available at Matpower's database [32] and their main characteristics are reported in Table 1. Observing the condition numbers of studied systems and compared them with those reported in [6], one can deduce that the studied systems are ill-conditioned. Nevertheless, it is confirmed later since NR failed in all cases using a flat start (see Section III).

B. BASE CASE SCENARIO
Firstly, the studied systems are solved at base case scenario. Table 2 reports the total number of iterations required by each considered PF solver at this scenario.  The Iwamoto's method typically required many iterations for converging. HOLM employed a reasonable number of iterations in the Systems #1 and #2, however, its convergence turned slow in the Systems #3 and #4. RK4 converged in the System #3 employed 25 iterations. RBS and BEM performed similarly employing 13-15 iterations for converging. Finally, 3S-SIA outperformed remainder methodologies solving all studied systems employing 4-5 iterations. Superior convergence features of 3S-SIA are also observed in Fig. 3, where the convergence curves are plotted. In this figure, it can be appreciated as 3S-SIA produced the lowest residual after third iteration, which is also characteristic of good robustness properties. Table 3 reports the execution time consumed by each considered PF solution technique for solving the base case scenario in the studied systems. As is clearly observed, the developed 3S-SIA was considerably faster than the remainder methods. This is undoubtedly due to its low computational burden. Table 4 reports the total number of factorizations required by each PF solver at base case scenario. As it was commented, factorization is the heaviest part of any PF calculation, in that sense, 3S-SIA was the most efficient  method due to it requires less factorizations than the remainder methodologies.

C. GENERATORS' REACTIVE LIMITS ENFORCEMENT
In this case, let us consider that the generators' reactive limits are enabled. To do that, strategy described in Section IV.E has been used. Table 5 summarizes the results obtained at this scenario. As it can be seen, more iterations are required in comparison with the base case scenario, which was expected since several PF problems have to be solved for achieving a feasible solution. Nevertheless, similar conclusions as for the base case scenario can be drawn in this case.

D. LIMIT LOAD CASES
The loading level of the system typically affects the convergence properties of PF solution techniques or provoking instability [36]. In that sense, it is interesting to study the performance of the developed 3S-SIA close to the maximum loadability point, where the condition of the system is notably deteriorated. To do that, the injected active and reactive power at PQ buses along the injected reactive power at PV buses, are multiplied by a loading factor namely λ ∈ R + . The loading level is increased in steps of 0.0001 pu until all considered methodologies failed, then, the immediately lower λ is taken. For example, for the System #2 we consider λ = 1.1586pu since λ = 1.1587pu gives rise divergence. Table 6 reports the results obtained at this scenario.
NR failed in all studied systems. The Iwamoto's method and BEM failed in the System #4. RK4 failed in the Systems #1, #2 and #4 and converged to the low voltage solution in the System #3. RBS failed in the Systems #3 and #4. HOLM and 3S-SIA solved all studied systems, however, HOLM obtained an inaccurate solution in the system #3. As for preceding experiments, 3S-SIA outperformed remainder methodologies in terms of iteration numbers and computation time.

E. INFLUENCE OF THE R/X RATIO
The overall R/X ratio of the system also affects the computational performance of PF solution methodologies in a similar way of the loading level [36]. We have modified the R/X ratio of the studied systems decreasing the branch impedances by a constant factor namely η, thus, if branch impedances are reduced 0.4, the overall R/X ratio is increased ∼1.67. In Fig. 4, the total number of iterations required by 3S-SIA for solving the studied systems for different R/X ratios is showed. As it can be seen, this topic does not considerably affect the performance of the developed algorithm, in fact, it simply needed a larger iteration in the System #2 when the R/X was increased, preserving its convergence properties in remainder cases.

F. CASCADING FAILURES
When an element of a power system fails (branch, transformers. . . ), the condition of the system is deteriorated and PF becomes harder to be solved [36]. One typical situation VOLUME 8, 2020 in power system operation is the cascading failure of various elements. To simulate this episode, we have generated 10 different scenarios which a predefined number of branches/transformers are out of service. Outage elements are determined randomly. Fig. 5 shows the average number of iterations, required for different methodologies for solving the PF of the generated failure scenarios in the System #1. The total number of scenarios successfully solved by each technique is also showed. One should note that RK4 and NR are not reported in this figure, since they even failed for solving the base case scenario.
As it can be seen, all methodologies solved the same number of scenarios. However, the developed algorithm was considerably more efficient, since it always required less iterations than the remainder considered methodologies.

VI. CONCLUSIONS AND FUTURE WORKS
A novel three-stage algorithm (3S-SIA) for efficiently solving realistic large-scale ill-conditioned systems has been developed. It is mainly based on a Semi-Implicit approach (SIA) but incorporates other numerical arrangements for enhancing either the robustness or efficiency features.
The developed 3S-SIA is quite efficient in the sense that it just requires a matrix factorization each iteration, resulting in similar computational burden in comparison NR. In addition, it was proved that 3S-SIA is notably more robust than NR.
Several numerical experiments in different ill-conditioned systems have been carried out. Various demand scenarios have been considered. Performance of the developed 3S-SIA has been compared with other well-known PF solution methods.
Results have showed that 3S-SIA constitutes a robust and efficient methodology. Among all studied methodologies, only 3S-SIA successfully solved all the considered scenarios outperforming the remainder considered methodologies. Finally, it is barely affected by the loading level or R/X ratios and it has turned out to be quite reliable and efficient in presence of cascading failures.
Future works should study the applicability of the developed algorithm in other related problems.