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

TL;DR

This paper introduces a novel holomorphic embedding continuation method that efficiently finds multiple power flow solutions by combining large step sizes with strategic switching around singularities, validated on IEEE test cases.

Contribution

The paper presents a new holomorphic embedding-based continuation method that improves efficiency in locating multiple power flow solutions by integrating large step sizes and singularity handling.

Findings

Method reduces computational steps compared to traditional routines.

Effective in traversing non-singular and singular solution curve segments.

Validated on standard IEEE test cases with favorable results.

Abstract

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.

Tables

Table 1. TABLE I: Numerical Results by Predictor-Corrector Method

Method	Predictor-Corrector
Case	overall steps	overall time (s)	# Solutions
3TS	2962	1.257	6
3	3349	1.143	6
4gs	4281	1.405	6
4BBc	8838	2.660	12
4BB0	13791	3.719	14
5Salam	11465	3.626	10
5loop	17049	4.568	10
6ww	9421	3.209	6
7Salam	5195	1.978	4
9	22264	10.945	8
14	151423	102.401	30
30	5358518	6054.987	472
33bw	311957	249.736	16
39	3009935	3758.195	176
57	14647351	23864.005	606

Table 2. TABLE II: Numerical Results by HEBC Method

Method	HEBC
Routine	Holomorphic		Predictor-Corrector		overall time (s)
Case	# steps	time (s)	# step	time (s)
3TS	253	0.373	228	0.133	0.859
3	290	0.345	387	0.155	0.798
4gs	481	0.625	766	0.197	1.248
4BBc	850	1.174	1358	0.362	2.102
4BB0	1212	1.489	2805	0.773	2.890
5Salam	1128	1.754	1224	0.290	2.661
5loop	1695	2.448	1627	0.358	3.537
6ww	995	1.402	1143	0.288	2.249
7Salam	564	1.068	362	0.132	1.676
9	1668	3.158	4026	1.668	5.572
14	13350	34.443	27238	12.784	50.013
30	403181	2077.966	910664	674.828	2813.249
33bw	15904	81.896	65351	43.519	129.323
39	184458	1247.794	1044166	930.796	2204.543
57	835550	10565.59	3078609	3598.361	14304.691

Equations

$$\mathbf{V}:={\mathbf{V}}_{\mathrm{d}}+j{\mathbf{V}}_{\mathrm{q}}$$

$$V_k^{\star }\sum _{n=1}^{N_{\mathrm{bus}}}{Y}_{n,k}{V}_n=S_k^{\star }$$

$$P_k$$

$$\displaystyle+V_{\rm{q},\it{k}}\sum_{n=1}^{N_{\rm{bus}}}\big{(}G_{n,k}V_{\rm{q},\it{n}}+B_{n,k}V_{\rm{d},\it{n}}\big{)}$$

$$Q_k$$

$$\displaystyle-V_{\rm{d},\it{k}}\sum_{n=1}^{N_{\rm{bus}}}\big{(}G_{n,k}V_{\rm{q},\it{n}}+B_{n,k}V_{\rm{d},\it{n}}\big{)}$$

$$P_k$$

$$\displaystyle+V_{\rm{q},\it{k}}\sum_{n=1}^{N_{\rm{bus}}}\big{(}G_{n,k}V_{\rm{q},\it{n}}+B_{n,k}V_{\rm{d},\it{n}}\big{)}$$

$$V_{\mathrm{m},k}^2$$

$$V_{\mathrm{m},s}^2$$

$$0$$

$$PF(\mathbf{U}):=\{f_i(\mathbf{U})={\mathbf{U}}^T{\mathbf{M}}_i\mathbf{U}-r_i,i=1,\cdots ,2N_{\mathrm{bus}}\}$$

$$\mathcal{E}:PF(\mathbf{U})\to EF(\mathbf{U})$$

$$EF(\mathbf{U}):=\{g_i(\mathbf{U})={\mathbf{U}}^T{\mathbf{H}}_i\mathbf{U}-{\gamma }_i,i=1,\cdots ,2N_{\mathrm{bus}}\}$$

$$E{F}_{l-}(\mathbf{U}):=EF(\mathbf{U})-\{g_l(\mathbf{U})\}$$

$$E{F}_{l,\alpha }(\mathbf{U},\alpha ):=E{F}_{l-}(\mathbf{U})\cup \{g_l(\mathbf{U})-\alpha ,\alpha \in \mathbb{R}\}.$$

$$\mathcal{Z}(E{F}_{l-};\mathbf{U})=\mathcal{Z}(E{F}_{l,\alpha };\mathbf{U})$$

$$\mathscr{Z}(EF_{l-};\mathbf{U})=\mathscr{Z}\big{(}\{PF(\mathbf{U})-\alpha\mathscr{E}^{-1}\mathbf{e}_{l}\};\mathbf{U}\big{)}$$

$$\displaystyle\mathscr{E}^{-1}\big{(}EF(\mathbf{U})-\alpha\mathbf{e}_{l}\big{)}$$

$$=PF(\mathbf{U})-\alpha {\mathcal{E}}^{-1}{\mathbf{e}}_l.$$

$$\mathscr{Z}(EF_{l,\alpha};\mathbf{U})=\mathscr{Z}\big{(}\{PF(\mathbf{U})-\alpha\mathscr{E}^{-1}\mathbf{e}_{l}\};\mathbf{U}\big{)}.$$

$$\mathscr{Z}(EF_{l-};\mathbf{U})=\mathscr{Z}\big{(}\{PF(\mathbf{U})-\alpha\mathscr{E}^{-1}\mathbf{e}_{l}\};\mathbf{U}\big{)}.$$

$$P_k(\alpha )$$

$$Q_k(\alpha )$$

$$\sum _{n=1}^{N_{\mathrm{bus}}}{Y}_{n,k}{V}_n(\alpha )$$

$$\displaystyle+K_{\rm{q},\it{k}}\alpha)Q_{k,0}\bigg{)}W_{k}^{\star}(\alpha^{\star})$$

$$V_k(\alpha )W_k(\alpha )$$

$$\displaystyle\sum_{n=1}^{N_{\rm{bus}}}\bigg{(}Y_{n,k}\sum_{i=0}^{\infty}\mathit{v}_{n,i}\alpha^{i}\bigg{)}$$

$$\displaystyle+K_{\rm{q},\it{k}}\alpha)Q_{k,0}\bigg{)}\sum_{i=0}^{\infty}\mathit{w}_{n,i}^{\star}\alpha^{i}$$

$$\sum _{i=0}^{\infty }{v}_{n,i}{\alpha }^i\sum _{i=0}^{\infty }{w}_{n,i}{\alpha }^i$$

$$\sum _{n=1}^{N_{\mathrm{bus}}}{Y}_{n,k}{V}_n(\alpha )$$

$$V_k(\alpha )V_k^{\star }({\alpha }^{\star })$$

$$V_k(\alpha )W_k(\alpha )$$

$$\displaystyle\sum_{n=1}^{N_{\rm{bus}}}\bigg{(}Y_{n,k}\sum_{i=0}^{\infty}\mathit{v}_{n,i}\alpha^{i}\bigg{)}$$

$$\displaystyle-j\sum_{i=0}^{\infty}q_{k,i}\alpha^{i}\bigg{)}\sum_{i=0}^{\infty}\mathit{w}_{n,i}^{\star}\alpha^{i}$$

$$\sum _{i=0}^{\infty }{v}_{n,i}{\alpha }^i\sum _{i=0}^{\infty }{v}_{n,i}^{\star }{\alpha }^i$$

$$\sum _{i=0}^{\infty }{v}_{n,i}{\alpha }^i\sum _{i=0}^{\infty }{w}_{n,i}{\alpha }^i$$

$$V_{\mathrm{s}}(\alpha )V_{\mathrm{s}}^{\star }({\alpha }^{\star })=V_{\mathrm{s},\mathrm{m}}^2+K_{\mathrm{s}}\alpha$$

$$v_{\mathrm{s},i}$$

Full text

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]$111A more restricted range may be assumed to be $[0.95,~{}1.05]$ for transmission systems. p.u. represents a steady state under which the system can be well-operated. 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 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 $1970$s when [12] examined a 3-node system which admits [math], $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_{\rm{bus}}-2}$, to a combinatorial bound, $C_{2N_{\rm{bus}}-2}^{N_{\rm{bus}}-1}$, where $N_{\rm{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 one222A very special type of test cases can be solved much more efficiently by some techniques from algebraic geometry. Interested readers are referred to [18]. However, there is no such efficient algebraic geometry method for solving a general power flow case at present.. 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 time333Currently, there is no rigorous theoretical guarantee to show that the elliptical formulation can always connect all the real solutions. It is an ongoing research., 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.

Showed an equivalent curve design for the elliptical formulation of the power flow problem; 2. 2.

Proposed a hybrid numerical continuation method HEBC for finding multiple power flow solutions; 3. 3.

Proposed a warm starter to quickly initiate the predictor-corrector routine for passing through singularities; 4. 4.

Showed that HEBC outperforms the traditional predictor-corrector algorithm [20] for all the tested cases; 5. 5.

Computed solution sets444Solution sets will be available online soon. 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.

II-A Power Flow Equations in Rectangular Coordinates

Consider a connected power grid with $N_{\rm{bus}}$ nodes. Let the node voltage vector be

[TABLE]

where $\mathbf{V}\in\mathbb{C}^{N_{\rm{bus}}}$; $\mathbf{V}_{\rm{d}}\in\mathbb{R}^{N_{\rm{bus}}}$ and $\mathbf{V}_{\rm{q}}\in\mathbb{R}^{N_{\rm{bus}}}$ are the real and imaginary parts of $\mathbf{V}$, respectively .

For the PQ bus we have

[TABLE]

where $V_{k}$ and $V_{n}$ are the corresponding entries of $\mathbf{V}$; $Y_{n,k}$ is the $(n,k)$-th entry of the bus admittance matrix $\mathbf{Y}\in\mathbb{C}^{N_{\rm{bus}}\times N_{\rm{bus}}}$; $S_{k}\in\mathbb{C}$ is the complex power load at bus $k$; superscript star $\star$ represents the conjugate operator.

Separating the real and imaginary parts of Equation (2) gives the two equations about a PQ bus

[TABLE]

where $P_{k}\leq 0$ and $Q_{k}\leq 0$555Usually a load absorbs reactive power, but it can possibly generate reactive power. In that case $Q_{k}\geq 0$. 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 $\mathbf{G}$ and the bus susceptance matrix $\mathbf{B}$666$\mathbf{Y}=\mathbf{G}+j\mathbf{B}$; $V_{\rm{d},\it{k}}$, $V_{\rm{d},\it{n}}$, $V_{\rm{q},\it{k}}$ and $V_{\rm{q},\it{n}}$ are the corresponding entries of $\mathbf{V}_{\rm{d}}$ and $\mathbf{V}_{\rm{q}}$, which are unknown variables that should be determined.

For the PV bus we have

[TABLE]

where $P_{k}$ is a fixed active power injection at bus $k$ which is usually positive but can be negative; $V_{\rm{m},\it{k}}$ is the fixed voltage magnitude at bus $k$.

For the slack bus with an angle reference we have

[TABLE]

where subscript $s$ is the slack bus number; $V_{\rm{m},\it{s}}$ is the slack bus voltage magnitude.

One can further substitute (5b) in (5a), (4) and (3) to eliminate $V_{\rm{q},\it{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

[TABLE]

where $\mathbf{U}:=[\mathbf{V}_{\rm{d}}^{T}~{}~{}\mathbf{V}_{\rm{q}}^{T}]^{T}$ is the unknown variable vector; $\mathbf{M}_{i}\in\mathbb{SR}^{2N_{\rm{bus}}\times 2N_{\rm{bus}}}$ is a symmetric constant matrix for the quadratic part; $r_{i}\in\mathbb{R}$ is the constant scalar part.

II-B Equivalent Curve Design of Elliptical Formulation of Power Flow Equations

As introduced in Section I, [19] presented a counter-example 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 $\mathscr{E}\in\mathbb{R}^{2N_{\rm{bus}}\times 2N_{\rm{bus}}}$ that sends Equation (6) to a set of high dimensional ellipses $EF(\mathbf{U})$. The construction of $\mathscr{E}$ can be found in [21, 20]. Consider

[TABLE]

with

[TABLE]

where $\mathbf{H}_{i}\in\mathbb{SR}^{2N_{\rm{bus}}\times 2N_{\rm{bus}}}$ and $\mathbf{H}_{i}\succ 0$; $\gamma_{i}>0$.

Let $\mathscr{Z}(h;x)$ be the operator that takes the projection of $\{(x,y)|h(x,y)=0\}$ onto $x$; define

[TABLE]

[TABLE]

Since $EF(\mathbf{U})$ defines a determined algebraic system, its algebraic set is generically 0-dimensional in $\mathbb{R}^{2N_{\rm{bus}}}$. By removing one equation from $EF(\mathbf{U})$, $EF_{l-}(\mathbf{U})$ acquires one degree of freedom and defines a 1-dimensional algebraic set in $\mathbb{R}^{2N_{\rm{bus}}}$. On the other hand, adding one extra degree of freedom to $EF(\mathbf{U})$ makes the algebraic set of $EF_{l,\alpha}(\mathbf{U},\alpha)$ 1-dimensional in $\mathbb{R}^{2N_{\rm{bus}}+1}$. The following Lemma 1 shows an equivalence between these two 1-dimensional algebraic sets.

Lemma 1**.**

[TABLE]

The proof is trivial and omitted here. Next, we state the equivalent curve design of elliptical formulation in Theorem 1.

Theorem 1**.**

[TABLE]

where $\mathbf{e}_{l}\in\mathbb{R}^{2N_{\rm{bus}}}$ is a unit column vector with the $j$-th entry being 1.

Proof.

By definition, $EF_{l,\alpha}(\mathbf{U},\alpha)$ can also be expressed as {$EF(\mathbf{U})-\alpha\mathbf{e}_{l}$}. Then we have

[TABLE]

Since $\mathscr{E}$ is an invertible linear map, it is a homeomorphism. Hence,

[TABLE]

Finally, by Lemma 1 we conclude that

[TABLE]

∎

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 $PF(\mathbf{U})-\alpha\mathscr{E}^{-1}\mathbf{e}_{l}$. In Section II this $\alpha$ is restricted to a real-valued scalar to support one extra degree of freedom. If we allow $\alpha$ 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.

III-A Holomorphic Embedding of Power Flow Equations

III-A1 PQ Bus Embedding

We start with the basic complex power balance equation for PQ bus in Equation (2). Note that $S_{k}=P_{k}+jQ_{k}$ we define

[TABLE]

where $\alpha\in\mathbb{C}$; $K_{\rm{p},\it{k}}$ and $K_{\rm{q},\it{k}}$ are obtained from $\mathscr{E}^{-1}\mathbf{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_{k}^{-1}$ for $V_{k}\neq 0$, and restrict parameterized $W_{k}(\alpha)$ to be reflective such that $W_{k}(\alpha)=W_{k}(\alpha^{\star})$, then Equation (2) can be written as

[TABLE]

Note that on the right hand side of (9a) we use $W_{k}^{\star}(\alpha^{\star})$ instead of $W_{k}^{\star}(\alpha)$ since they are equal by the reflective property777A more detailed discussion on the reflective requirement can be found in [33]..

Since $V_{k}(\alpha)$ and $W_{k}(\alpha)$ are holomorphic [33], we can use power series to represent them. Then, (9) can be re-written as

[TABLE]

where $\mathit{v}_{n,i}$ and $\mathit{w}_{n,i}$ are the power series coefficients.

Matching up coefficients for every monomial of $\alpha$ in (10a) and (10b) we can solve $(\mathit{v}_{k,1},\mathit{v}_{k,2},\cdots)$ and $(\mathit{w}_{k,1},\mathit{w}_{k,2},\cdots)$ recursively as long as $\mathit{v}_{k,0}$ and $\mathit{w}_{k,0}$ are provided.

III-A2 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_{k}^{-1}$ for $V_{k}\neq 0$ and restricting parameterized $W_{k}(\alpha)$ to be reflective we have the holomorphic embedded equations

[TABLE]

where $V_{k,\rm{m}}\in\mathbb{R}$ is the fixed voltage magnitude at bus $k$; $K_{\rm{v},\it{k}}$ is obtained from the corresponding entry of $\mathscr{E}^{-1}\mathbf{e}_{l}$.

By the holomorphic structure, we represent parameterized unknowns $V_{n}(\alpha)$, $W_{k}(\alpha)$, and $Q_{k}(\alpha)$ through their power series. Then, (11) are re-written as

[TABLE]

where $q_{k,i}$’s are the power series coefficients of $Q_{k}(\alpha)$.

Matching up coefficients for every monomial of $\alpha$ in (12a), (12b) and (12c) we can solve $u_{i}$, $w_{i}$, and $q_{i}$ as well.

III-A3 Slack Bus Embedding

Consider the slack bus voltage magnitude equation (5a). Its holomorphic embedded equation is

[TABLE]

where $V_{\rm{s},m}$ is the slack bus voltage magnitude, $K_{\rm{s}}$ is the corresponding entry from $\mathscr{E}^{-1}\mathbf{e}_{l}$.

Substituting the power series of $V_{\rm{s}}(\alpha)$ into Equation (13) and matching up each monomial of $\alpha$ we have

[TABLE]

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_{\rm{bus}}+N_{\rm{gen}}-3)\times(4N_{\rm{bus}}+N_{\rm{gen}}-3)$ where $N_{\rm{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_{\rm{max}}$888[25] claims that degree $i$ will deplete double precision digits after $60$. How to choose an appropriate $i_{\rm{max}}$ is beyond the scope of this paper. We choose $i_{\rm{max}}=15$ in our numerical experiments by empirical experience considering speed and accuracy..

III-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 $\mathit{v}_{k}(\alpha)$ for some $k$. Suppose its first $N$ coefficients are known.

[TABLE]

Let its Padé approximation be

[TABLE]

where we specify $N_{\rm{n}}+N_{\rm{d}}=N$, $N_{\rm{n}}\geq N_{\rm{d}}$, and $N_{\rm{n}}-N_{\rm{d}}\leq 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)\times(N+1)$ complex-valued sparse linear system. If we compute the power series to the maximum degree $i_{\rm{max}}$, the system size in the real space is $2(i_{\rm{max}}+1)\times 2(i_{\rm{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 threshold999In our numerical experiments, this threshold is set at $10^{-3}$ p.u. 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.

IV-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].

We start Algorithm 1 with a known solution $x_{1}$ which can be solved by Newton’s method or other techniques101010This step relies on the past extensive research of solving a high voltage solution to the power flow problem. Many mature solvers are able to do this job for very large systems.. After several initialization steps, designing the curve $\{PF-\alpha\mathscr{E}^{-1}\mathbf{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.

IV-B Inner Loop for Curve Tracing

Instead of tracing a curve by the traditional predictor-corrector 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.

IV-B1 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 $|\alpha_{k_{h}+1}-\alpha_{k_{h}}|$ is smaller than a threshold value $d\alpha_{h,\rm{min}}$. Both suggest that current holomorphic step is close to singular (or at least badly scaled with respect to $\alpha$).

IV-B2 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_{\rm{hp}}$ from the singular point to the last holomorphic point. We specifically choose the initial step interval $S_{\rm{pc}}$ to be $1/5$ of the estimated distance $d_{\rm{hp}}$ and to be no greater than $0.45$ of the last holomorphic step size. Then, using $S_{\rm{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_{\rm{pc}}$. It makes the predictor-corrector routine quickly pass through the singular point as shown by the rest green triangles.

IV-B3 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 $(\alpha_{m_{c}}-\alpha_{m_{c}-1})(\alpha_{m_{c}+1}-\alpha_{m_{c}})<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_{\rm{m}}$.

[TABLE]

As long as $R_{\rm{m}}$ drops to a threshold $R_{\rm{max}}$, say, $2\times 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_{\rm{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_{\rm{max}}$ times to generate coefficients for all degrees. Therefore, the computational complexity for the power series is

[TABLE]

In Padé approximation, the complexity is

[TABLE]

The total complexity of a holomorphic prediction is $C_{\rm{Holo}}=C_{\rm{Tl}}+C_{\rm{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

[TABLE]

Suppose $i_{max}$ is fixed, $N_{gen}=0.2N_{bus}$111111The number of PV buses usually occupies a small fraction of the total number of buses., and each corrector takes 3 Newton’s iterations to converge for both the holomorphic step and the traditional predictor-corrector step, we have

[TABLE]

It suggests that one holomorphic step takes about four predictor-corrector steps computations asymptotically with the 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.

VI Numerical Experiments

This section presents a comprehensive numerical evaluation of the proposed HEBC method on several standard power system test cases including “case3TS”, “case3”, “case4gs”,“case4BBc”, “case4BB0”, “case5Salam”, “case6ww”, “case7Salam”, “case9”, “case14”, “case30”, “case33bw”, “case39”, “case57”121212Tap ratios are removed in this case to reduce the number of solutions., which can be found in the Matpower libary [36], and “case5loop” [19]. To avoid numerical instability and structurally unstable solutions, small resistance at $10^{-4}$ p.u. is added to lossless lines. The HEBC method and the full predictor-corrector method are coded in Matlab R2017b and executed on a PC with $2.8$GHz Intel i7-7700HQ CPU and $16$GB RAM.

VI-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 case14131313No existing literature claims complete solution sets for larger IEEE test cases.. 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 comparison141414A more recent progress in [16] successfully reduced the computational time of case14 to 5 minutes, however, the proposed HEBC is still much faster..

VI-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. Both methods provide the same solution sets for all cases, but the HEBC method is more efficient than the traditional predictor-corrector method. Some hard151515A curve is hard to follow in the sense that it contains too many singularities or some singular points are very sharp when turning directions. 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.

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.

To reveal the efficiency of HEBC, we compute the equivalent number of predictor-corrector steps $N_{\rm{eqv}}$

[TABLE]

where $N_{\rm{pc}}$ is the number of full predictor-corrector steps; $N_{\rm{he,pc}}$ is the number of predictor-corrector routine steps in HEBC; and $N_{\rm{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$ predictor-corrector 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 predictor-corrector 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.

VI-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

[TABLE]

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

[TABLE]

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.

Acknowledgments

We would like to acknowledge the helpful discussions with Dr. Wenqiang Feng from DST Systems, Dr. Honghao Zheng from Siemens, Prof. Konstantin Turitsyn at MIT and Prof. Bernard Lesieutre at UW-Madison. The authors gratefully acknowledge support from the National Science Foundation under grant CRISP 1735513.

Appendix A Sample Curves

Sample curves from simulations.