A Holomorphic Embedding Based Continuation Method for Identifying Multiple Power Flow Solutions

In this paper, we propose an efficient continuation method for locating multiple power flow solutions. We adopt the holomorphic embedding technique to represent solution curves as holomorphic functions in the complex plane. The holomorphicity, which provides global information of the curve at any regular point, enables large step sizes in the path-following procedure such that non-singular curve segments can be traversed with very few steps. When approaching singular points, we switch to the traditional predictor-corrector routine to pass through them and switch back afterward to the holomorphic embedding routine. We also propose a warm starter when switching to the predictor-corrector routine, i.e. a large initial step size based on the poles of the Pad\'{e} approximation of the derived holomorphic function, since these poles reveal the locations of singularities on the curve. Numerical analysis and experiments on many standard IEEE test cases are presented, along with the comparison to the full predictor-corrector routine, confirming the efficiency of the method.


I. INTRODUCTION
The electric power grid is a critical energy infrastructure for power generation, transmission, and distribution in modern society. The inherent nonlinearity of power grid introduces a great challenge to analyze its dynamical behaviors when subject to disturbances, especially when penetrated with a large amount of intermittent renewable energies. Identifying the region of attraction about the operating condition, i.e. a stable equilibrium point (SEP) of the underlying dynamical system, can significantly improve the situational awareness and, therefore, will be of great importance to avoid blackouts. Characterizing this region requires the knowledge of a special type of unstable equilibrium point (UEP) which is called the type-1 UEP [1], [2]. Determining them usually requires locating all nearby equilibria. In classical model [3], equilibria are the solutions to the power flow equations [4], [5], [6].
A high-voltage solution in the range of [0.9, 1.1] 1 p.u. represents a steady state under which the system can be welloperated. This solution is usually the SEP in transient stability analysis, while other solutions are UEPs. For tree structured networks, the high-voltage solution is unique [7]. However, it †: Laboratory for Information and Decision Systems, Massachusetts Institute of Technology, Cambridge, MA, danwumit@mit.edu; ‡: Department of Electrical and Computer Engineering, Taxes A&M University, College Station, TX, binwang@tamu.edu 1 A more restricted range may be assumed to be [0.95, 1.05] for transmission systems.
is possible that mesh networks can have multiple high-voltage solutions with either circulated flow [8] or reversed power flow [9]. Although being avoided in the normal operations, circulated flow can happen during fault transients. Meanwhile, the reversal of power flow can become very common in future power grids as distributed energy resources (DER) keeps penetrating to distribution networks. To better characterize stability region and to examine other high-voltage operating points, finding multiple power flow solutions plays a key role.
Nowadays a single high-voltage solution can be solved very efficiently. For example, systems with about 10,000 buses can be solved within a second [10]. However, the largest system that can be provably solved for all solutions is the 14 bus system [11]. The efficiency of solving a single high-voltage solution comes from the knowledge of a good initial guess for local solvers to converge. But it is rather hard to acquire appropriate initial guesses for other solutions. If using random seeds, the complexity increases exponentially as the system size increases. Therefore, a systematic method is required to find these solutions.
Early attempts to find multiple power flow solutions dates back to 1970s when [12] examined a 3-node system which admits 0, 2, 4 or 6 solutions. In 1989, [13], [14] introduced the probability-one homotopy continuation method to find all the complex-valued solutions to the power flow problem. The homotopy continuation method requires estimating the total number of solutions to the power flow problem, which is still an ongoing research. In 1982, [15] sharpened the solution number bound from the classic Bezout's bound, 2 2N bus −2 , to a combinatorial bound, C N bus −1 2N bus −2 , where N bus is the number of nodes in a power grid. Recently, [11] applied a polyhedral homotopy continuation method to completely solve the IEEE standard 14-bus system by the Bernstein-Khovanskii-Kushnirenko (BKK) bound which is sharper than the Bezout's bound. However, evaluating the BKK bound is very expensive. To further explore a simpler bound, [16] introduced the adjacent polytope bound, which is sharper than the BKK bound and more computable.
While progressive, the homotopy method usually ends up with a huge amount of complex-valued solutions which are fictitious power flow solutions. To only identify actual power flow solutions, [17] introduced the idea of curve design which connects different real solutions by some 1-dimensional curves. Following these curves power flow solutions can be reached one by one 2 . Though efficient, [19] provided a counter-example for [17]. To rectify their method, an elliptical formulation of the power flow problem is used in [20] to restrict the curve design on high dimensional ellipses. It helps solve all the standard IEEE test cases which can be verified by the homotopy method in a reasonable time 3 , including the counter-example in [19]. The existence and construction of elliptical formulation were provided in [21] and extended to the optimal power flow problem to find multiple local extrema for hard problems in [22].
The curve tracing routine performed in [20], [21], [22] is a traditional predictor-corrector algorithm which adopted a quadratic predictor [23], Newton's method for corrector, and an adaptive step-length control [24]. Many variations of the predictor-corrector algorithm exist, however, most of them depend only on the local information or previously solved points. To accelerate the curve tracing, in this paper we design a new hybrid algorithm called the holomorphic embedding based continuation (HEBC) method to replace the traditional predictor-corrector algorithm during most of the curve tracing periods. It applies the holomorphic embedding technique to quickly pass through the non-singular curve segments by utilizing the global information of that curve, and uses a predictor-corrector routine to travel across singularities.
The holomorphic embedding method (HEM) was introduced by Trias [25] in 2012 as a new power flow solver. The basic idea is to parameterize a polynomial system by an extra free variable and acquires the solution curve information by power series. Early attempts to use parameterization and power series for solving power flows started with [26] and followed by [27], [28]. Recently, HEM was extended to some applications with different modelings [29], [30], [31], [32]. Two features of HEM are particularly useful in our circumstance to improve searching efficiency. First, HEM can release us from local predictor-corrector scheme and provide with very long arc steps on the solution curve. This can largely reduce the burden of repeatedly solving linear systems in the corrector part. Moreover, the smallest real-valued pole of Padé approximation can be used to design an appropriate step length when passing through singular point. It avoids overly large step sizes to improve numerical stability, and keeps step sizes progressive to maintain efficiency.
The contributions of this paper are summarized below.
1) Showed an equivalent curve design for the elliptical formulation of the power flow problem; 2) Proposed a hybrid numerical continuation method HEBC for finding multiple power flow solutions; 3) Proposed a warm starter to quickly initiate the predictorcorrector routine for passing through singularities; 4) Showed that HEBC outperforms the traditional predictor-corrector algorithm [20] for all the tested cases; 5) Computed solution sets 4 for several large test cases which currently are intractable by homotopy continuation method or the similar.

II. DESCRIPTION OF POWER FLOW PROBLEM
Throughout this paper we adopt the power flow formulation in rectangular coordinates.

A. Power Flow Equations in Rectangular Coordinates
Consider a connected power grid with N bus nodes. Let the node voltage vector be where V ∈ C N bus ; V d ∈ R N bus and V q ∈ R N bus are the real and imaginary parts of V, respectively . For the PQ bus we have where V k and V n are the corresponding entries of V; Y n,k is the (n, k)-th entry of the bus admittance matrix Y ∈ C N bus ×N bus ; S k ∈ C is the complex power load at bus k; superscript star ⋆ represents the conjugate operator. Separating the real and imaginary parts of Equation (2) gives the two equations about a PQ bus where P k ≤ 0 and Q k ≤ 0 5 are the fixed active and reactive power loads at bus k; G n,k and B n,k are the (n, k)-th entries of the bus conductance matrix G and the bus susceptance matrix B 6 ; V d,k , V d,n , V q,k and V q,n are the corresponding entries of V d and V q , which are unknown variables that should be determined.
For the PV bus we have 4 Solution sets will be available online soon. 5 Usually a load absorbs reactive power, but it can possibly generate reactive power. In that case Q k ≥ 0.
where P k is a fixed active power injection at bus k which is usually positive but can be negative; V m,k is the fixed voltage magnitude at bus k. For the slack bus with an angle reference we have where subscript s is the slack bus number; V m,s is the slack bus voltage magnitude.
One can further substitute (5b) in (5a), (4) and (3) to eliminate V q,s . Finally, (3), (4), and (5) together are the power flow equations we will investigate in this paper. Note that they are in quadratic form, thus can be written succinctly as bus is a symmetric constant matrix for the quadratic part; r i ∈ R is the constant scalar part.

B. Equivalent Curve Design of Elliptical Formulation of Power Flow Equations
As introduced in Section I, [19] presented a counterexample that fails the proposed algorithm in [17] for finding all the real-valued power flow solutions. Then, [20] introduced the concept of elliptical formulation of power flow equations which substantially changes the topology of path following curves and succeeded for that example. Later, [21] showed the existence of elliptical formulation under mild conditions and constructed it in a systematical way.
We start our discussion with a given invertible linear map E ∈ R 2N bus ×2N bus that sends Equation (6) to a set of high dimensional ellipses EF (U). The construction of E can be found in [21], [20]. Consider Since EF (U) defines a determined algebraic system, its algebraic set is generically 0-dimensional in R 2N bus . By removing one equation from EF (U), EF l− (U) acquires one degree of freedom and defines a 1-dimensional algebraic set in R 2N bus . On the other hand, adding one extra degree of freedom to EF (U) makes the algebraic set of EF l,α (U, α) 1-dimensional in R 2N bus +1 . The following Lemma 1 shows an equivalence between these two 1-dimensional algebraic sets.
The proof is trivial and omitted here. Next, we state the equivalent curve design of elliptical formulation in Theorem 1.
where e l ∈ R 2N bus is a unit column vector with the j-th entry being 1.
Proof: By definition, EF l,α (U, α) can also be expressed as {EF (U) − αe l }. Then we have Since E is an invertible linear map, it is a homeomorphism. Hence, Finally, by Lemma 1 we conclude that

III. HOLOMORPHIC EMBEDDING TECHNIQUE
Theorem 1 states that the 1-dimensional curves derived from the elliptical formulation EF l− can be acquired alternatively from a particular parameterized power flow problem P F (U) − αE −1 e l . In Section II this α is restricted to a realvalued scalar to support one extra degree of freedom. If we allow α to be a complex number, the parameterized curve resides in the complex plane and becomes a 2-dimensional surface in the real space. If this complex-value parameterized curve happens to be governed by holomorphic functions, it is called the holomorphic embedding. The advantage of being holomorphic is that the global information of the embedded curve is determined and singularities on the curve can be predicted by analytic continuation techniques.

A. Holomorphic Embedding of Power Flow Equations 1) PQ Bus Embedding:
We start with the basic complex power balance equation for PQ bus in Equation (2). Note that where α ∈ C; K p,k and K q,k are obtained from E −1 e l for some l; P k,0 and Q k,0 are the fixed starting active and reactive power which admit a known solution.
If we define a new variable W k := V −1 k for V k = 0, and restrict parameterized W k (α) to be reflective such that W k (α) = W k (α ⋆ ), then Equation (2) can be written as Note that on the right hand side of (9a) we use W ⋆ k (α ⋆ ) instead of W ⋆ k (α) since they are equal by the reflective property 7 . Since V k (α) and W k (α) are holomorphic [33], we can use power series to represent them. Then, (9) can be re-written as where v n,i and w n,i are the power series coefficients. Matching up coefficients for every monomial of α in (10a) and (10b) we can solve (v k,1 , v k,2 , · · · ) and (w k,1 , w k,2 , · · · ) recursively as long as v k,0 and w k,0 are provided.
2) PV Bus Embedding: Next, we consider the holomorphic embedding for PV bus equations. To retain holomorphicity, we need to bring back the reactive power balance equation (3b) to (4) and consider reactive power input as a new variable. Again, by defining W k := V −1 k for V k = 0 and restricting parameterized W k (α) to be reflective we have the holomorphic embedded equations where V k,m ∈ R is the fixed voltage magnitude at bus k; K v,k is obtained from the corresponding entry of E −1 e l . By the holomorphic structure, we represent parameterized unknowns V n (α), W k (α), and Q k (α) through their power series. Then, (11) are re-written as where q k,i 's are the power series coefficients of Q k (α).
Matching up coefficients for every monomial of α in (12a), (12b) and (12c) we can solve u i , w i , and q i as well.
3) Slack Bus Embedding: Consider the slack bus voltage magnitude equation (5a). Its holomorphic embedded equation 7 A more detailed discussion on the reflective requirement can be found in [33].
where V s,m is the slack bus voltage magnitude, K s is the corresponding entry from E −1 e l .
Substituting the power series of V s (α) into Equation (13) and matching up each monomial of α we have Combining the corresponding equations from the PQ bus, PV bus and slack bus equations we finally solve the power series coefficients for each degree-i. In practice, every degree requires solving a real-valued linear system (sparse) with its size (4N bus + N gen − 3) × (4N bus + N gen − 3) where N gen is the number of PV nodes. As i goes to infinity, the power series converges to the actual curve in the convergence range. To compromise accuracy and speed, we usually stop at a given maximum degree i max 8 .

B. Padé Approximation
The above subsection shows that each node voltage (as well as reactive power at PV bus) can be embedded as a holomorphic function, and demonstrates a recursive way to obtain the coefficients. In practice the holomorphic function can only be evaluated by a finite sequence of power series. Thus, the accuracy of the sequence deteriorates when approaching the singularities of the holomorphic function. To achieve a better convergence performance and to predict the location of singular point, we further compute the Padé approximation. It approximates the holomorphic function by a rational function in which the numerator and denominator are polynomials. According to [34], [35], the Padé approximation has the maximum convergent domain if the degrees of its numerator and denominator have the minimum difference. It provides a criterion for determining the best degree(s) that should be chosen.
Consider an embedded voltage variable v k (α) for some k. Suppose its first N coefficients are known.
Let its Padé approximation be where we specify N n +N d = N , N n ≥ N d , and N n −N d ≤ 1.
To reach a unique coefficient set, let l k,0 = 1. Matching up the coefficients for each monomial we can solve u k,n 's and l k,n 's in a (N + 1) × (N + 1) complex-valued sparse linear system. If we compute the power series to the maximum degree i max , the system size in the real space is 2(i max + 1) × 2(i max + 1).
Once the Padé approximation has been calculated, we can move along the parameterized curve by evaluating Padé approximated values until a power mismatch threshold 9 has been reached. We can also compute the real-valued zeros to the denominator function of Padé. These zeros reveal the locations of singularities on the parameterized curve, which can further assist us designing appropriate arc length for passing through these singular points by the traditional predictor-corrector algorithm. Next section will discuss these designs in detail.

IV. HOLOMORPHIC EMBEDDING BASED CONTINUATION METHOD
The proposed HEBC method can be divided into an outer loop part and an inner loop part. The outer loop focuses on new solution updates and sequential curve designs; while the inner loop primarily follows the curve fed by the outer loop and returns the solution set found on that curve.

A. Outer Loop for Solution Search
To make this article self-sustained, we briefly explain the search strategies in the outer loop and summarize it in Algorithm 1. Interested readers can refer to [20]. k ← k + 1 ⊲ Update counting number 8: x 0 ← x k ⊲ Update starting solution 9: for l = 1, 2, · · · , N eqn do 10: We start Algorithm 1 with a known solution x 1 which can be solved by Newton's method or other techniques 10 . After several initialization steps, designing the curve {P F − αE −1 e l } which is equivalent to {EF l− } for l. Following the curve from l = 1 to the last one by Algorithm 2 (which will be discussed shortly below) and collect new solutions. When finished tracing curves, assigning the starting point x 0 to a newly found solution, say, x 2 , and repeating the procedure.
The whole loop terminates upon every solution having been assigned to a starting point.
Algorithm 1 presents a procedure to follow each curve sequentially. However, the curve designs at the same starting solution are independent with each other, suggesting a parallel computing framework to simultaneously trace these curves. The parallel computing is not performed in this article, but can be done with ease and increase speed drastically.

B. Inner Loop for Curve Tracing
Instead of tracing a curve by the traditional predictorcorrector algorithm, we apply the holomorphic embedding technique to quickly pass through the regular curve segments. The predictor-corrector algorithm is only executed for traveling across singularities. It is switched back to the holomorphic embedding as soon as current steps leave a singular point. Figure 1(a) shows four holomorphic steps on a selected curve from a 5-bus case [13]. They reach the singular point very quickly. On the other hand, the blue curve in Figure 1(a) was generated by the traditional predictor-corrector algorithm. It took dozens of steps to reach the same singularity.

1) Criterion to Enter Predictor-Corrector Routine:
Two indicators are considered to trigger the switch of the algorithms in Algorithm 2. The first indicator appears when the corrector steps fail to make the holomorphic prediction converge within a certain number of iterations. Another indicator comes when |α k h +1 − α k h | is smaller than a threshold value dα h,min . Both suggest that current holomorphic step is close to singular (or at least badly scaled with respect to α).

2) Using A Warm Starter to Accelerate Predictor-Corrector
Steps: One can initiate the predictor-corrector routine from a minimum step size, and increase it gradually. We refer it to a cold starter. To avoid slow "warming up" steps, a warm starter is proposed and implemented. It relies on an estimated distance d hp from the singular point to the last holomorphic point. We specifically choose the initial step interval S pc to be 1/5 of the estimated distance d hp and to be no greater than 0.45 of the last holomorphic step size. Then, using S pc to compute two backward steps to initiate a quadratic predictor. For example, the first two green triangles on the upper curve segment in Figure 1(b) are the backward points evaluated by Padé approximation at the step length S pc . It makes the predictor-corrector routine quickly pass through the singular point as shown by the rest green triangles. Initialize the 1st holomorphic step size δ h . 6: Prepare parameters for holomorphic embedding. 7: Compute power series of holomorphic embedding. 8: Compute Padé approximation. 9: Evaluate voltage values from Padé and update α k h +1 . 10: Evaluate power mismatch dP mis from computed voltages. 11: while minimum pole p min is not determined do 12: Compute roots {ζ i } from Padé denominator. 13: p min ← ζ min if the minimum real root ζ min has correct sign. 14: end while 15: Increase δ h while dP mis < dP max and |current point| < |p min |. 16: Decrease δ h while dP mis ≥ dP max or |current point| ≥ |p min |.

17:
Correct current holomorphic predicted point by Newton's method. 18: if Correction succeeds then 19: Record current point. 20: else 21: Delete current point and compute a starter for switching algorithm. 22: Break. 23: end if 24: if α k h +1 α k h < 0 then 25: Find a solution nearby. 26: if Fail to locate the solution then 27: Delete current point and compute a cold starter for switching algorithm. 28

3) Criterion to Exit Predictor-Corrector Routine:
When travelling across a singular point, the direction of curve changes. Numerically, there exists a particular step m c such that (α mc − α mc−1 )(α mc+1 − α mc ) < 0. After this moment, we continue the predictor-corrector routine for a while until the curve's slope value returns from infinity back to a tractable value. Instead of evaluating the actual slope of the curve, we monitor the maximum variable secant slope R m . (17) As long as R m drops to a threshold R max , say, 2 × 10 4 , we jump out of the predictor-corrector routine and start a new sequence of holomorphic steps.

V. COMPUTATIONAL COMPLEXITY COMPARISON
The holomorphic prediction consists of two sub-routines: 1) construct the power series; 2) compute Padé approximation based on the power series. Both sub-routines require solving sparse linear systems. The sparsity reduces computational efforts in practice but makes analysis hard. To get a rough idea of the complexity, we simply assume the matrices are dense in the analysis, but solve them in sparse form practically.
In Section III computing the power series coefficients requires solving a sequence of linear systems up to the highest degree i max . A favorable observation is that all these linear systems share the same constant matrix. Thus, the LU factorization only needs to be performed once, while forward and backward substitutions need to be performed i max times to generate coefficients for all degrees. Therefore, the computational complexity for the power series is In Padé approximation, the complexity is The total complexity of a holomorphic prediction is C Holo = C Tl + C Pd . On the other hand, in the traditional predictor-corrector algorithm the Newton's iterations in correctors are the most computational complex part. Again, suppose a dense Jacobian matrix (sparse in practice) the complexity of solving one Newton's iteration is Suppose i max is fixed, N gen = 0.2N bus 11 , and each corrector takes 3 Newton's iterations to converge for both the holomorphic step and the traditional predictor-corrector step, we have It suggests that one holomorphic step takes about four predictor-corrector steps computations asymptotically with the 11 The number of PV buses usually occupies a small fraction of the total number of buses. dense matrix LU factorization. So an average holomorphic step size which is greater than 4 times the average step size of the predictor-corrector algorithm can potentially reduce the computational time under the same assumptions.

A. Comparison To Homotopy Continuation Method
To demonstrate the superiority of computational efficiency in finding multiple power flow solutions, we begin with a comparison of the proposed HEBC method to the homotopy continuation method. The homotopy continuation is performed by the PHCpack [37].
The HEBC method finds all the actual power flow solutions in this comparison as well as case14 13 . Figure 2 shows execution time (in logarithmic scale) comparison between two methods. For test cases smaller than 5 buses, the PHCpack runs faster than the proposed HEBC method. However, for cases more than 5 buses, the HEBC outperforms the homotopy continuation method substantially. Considering the HEBC method is coded in Matlab and is not optimized to reach the most computational performance, the time reductions from HEBC are impressive. Test cases larger than 9 buses cannot be solved by PHCpack within 24 hours, thus are not considered in this comparison 14 .

B. Comparison To Full Predictor-Corrector Algorithm
In this part, we testify the traditional full predictor-corrector method from [20] and the proposed HEBC method on the same set of test cases, and compare their numerical performances. 12 Tap ratios are removed in this case to reduce the number of solutions. 13 No existing literature claims complete solution sets for larger IEEE test cases.
14 A more recent progress in [16] successfully reduced the computational time of case14 to 5 minutes, however, the proposed HEBC is still much faster.
Both methods provide the same solution sets for all cases, but the HEBC method is more efficient than the traditional predictor-corrector method. Some hard 15 sample curves are presented in Appendix A. One can see from the left plots of Figure 6 that the traditional full predictor-corrector method, though with quadratic predictor and automatic step length adaption, takes very dense points to trace curves. On the other hand, the right plots of Figure 6 are primarily sparse. Small dense point periods only occur around singularities when HEBC switches to the predictor-corrector routine for passing through those singularities. Summaries of the numerical results are collected in Table I and II.  HEBC provides the same solution sets for all the cases as in Table  I.
Comparing the results in Figure 3, the total number of steps for HEBC is about 1/6 to 1/3 of the total number of steps for the full predictor-corrector method. This ratio, not surprisingly, should depend on the problem structure. In general, fewer singularities and longer horizontal curve segments favor the HEBC more.

Fig. 4: Execution Times with Full Predictor-Corrector and HEBC
To reveal the efficiency of HEBC, we compute the equivalent number of predictor-corrector steps N eqv N eqv := (N pc − N he,pc )/N he,holo (22) where N pc is the number of full predictor-corrector steps; N he,pc is the number of predictor-corrector routine steps in HEBC; and N he,holo is the number of holomorphic routine steps in HEBC. From Table I and II we calculate that one holomorphic step on average can represent 8.5 predictorcorrector steps, with the worst case of 7 steps and the best case of 15 steps. In Figure 4 the first 9 small cases up to case7Salam show a limited time saving by HEBC. However, starting at case9 the HEBC method outperforms the full predictorcorrector method by up to 50% of the execution time. Larger cases also exhibit at least 30% time saving in the lower plots of Figure 4.

C. Average Number of Steps on Each Dimension
Recall that the HEBC method calls the Newton's method at each step to correct the predicted point. These predicted points are sequentially determined over the curve tracing process. Thus, the HEBC method can be regarded as a systematic way to choose initial points for solving the power flow equations, where the number of initial points equals the number of steps in Table II, i.e. the sum of entries in the second and forth columns for each case. From this point of view, one can assess the efficiency of HEBC by computing the average number of initial points (steps) allocated in each dimension R eq := N 1/d (23) where N is the total number of initial points, d is the dimensionality of the problem. R eq represents the number of points required in each single dimension such that the total number of initial points composed by their direct combinations achieves the same amount of initial points N for the whole d-dimensional problem. Specifically for our problem, R eq is computed as R eq = (N he,pc + N he,holo ) 1/(2N bus −1)  Figure 5 depicts the trend of R eq as system size increases. One can clearly see that the average number of steps distributed on each dimension decreases to nearly 1. Hence, despite the increase of total number of steps, the average number of steps on each dimension seems to decrease in an asymptotic sense.

VII. CONCLUSIONS
In this paper, we proposed an efficient hybrid method to solve multiple power flow solutions. We derived an equivalent curve design to the elliptical formulation of the power flow equations. Based on this design, a holomorphic embedding continuation method was introduced to replace the traditional predictor-corrector algorithm for regular curve tracing. Singular points were passed by the predictor-corrector routine. The complexity of one holomorphic step is around four times the complexity of a predictor-corrector step under certain assumptions. Numerical simulations showed that one holomorphic step size is equivalent to over eight predictor-corrector step size on average, and saved up to half of the computational time for some large test cases.
A possible future direction of research can use the proposed method to find multiple power flow solutions for dynamic stability analysis, especially in characterizing the stability boundary of an equilibrium point. Another interesting topic would be using this method for solving optimal power flow problems.