{'fraction_non_alphanumeric': 0.039092998145020624, 'fraction_numerical': 0.01686093191948836, 'mean_word_length': 4.297741273100616, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the Network Topology Dependent Solution Count of the Algebraic Load Flow Equations

A large amount of research activity in power systems areas has focused on developing computational methods to solve load flow equations where a key question is the maximum number of isolated solutions.Though several concrete upper bounds exist, recent studies have hinted that much sharper upper bounds that depend the topology of underlying power networks may exist. This paper establishes such a topology dependent solution bound which is actually the best possible bound in the sense that it is always attainable. We also develop a geometric construction called adjacency polytope which accurately captures the topology of the underlying power network and is immensely useful in the computation of the solution bound. Finally we highlight the significant implications of the development of such solution bound in solving load flow equations.

1. Introduction. Engineers are regularly required to perform power flow computations for designing, operating, and controlling power systems [36]. The solutions of the power flow or load flow equations, a system of multivariate nonlinear equations, are used to ensure the normal operating conditions as well as to perform contingency, stability, and bifurcations analysis. In addition, these solutions, which we refer to as solutions of the load flow equations from now on, are also important in system security assessment and optimal dispatching. In general, load flow equations may have more than one solutions [34]. There are indeed quite a few existing methods for finding one or many load flow solutions [1, 5, 8-10, 16, 18-20, 23, 30, 31, 38, 41, 43-45, 47, 49, 52, 54, 57-60, 62, 63] (see [48] for a recent review). Out of the few methods that guarantee to find all load flow solutions, i.e., the interval based approach [57], Gröbner bases technique [16,54,58,59] and the numerical polynomial homotopy continuation (NPHC) method [8,9,41,47,49,60], the NPHC method appears most promising in scalability with increasing system sizes in that it has already found all load flow solutions of up to IEEE 14 bus systems [49] (and 18 oscillators case for the Kuramoto model [47]) and is inherently parallel: formulating load flow equations as system of polynomial equations, the NPHC method, whose roots are in complex algebraic geometry, finds all isolated complex systems which obviously include all the isolated real solutions. In all these computational methods, the knowledge of the number of solutions play a crucially important role.
While several concrete upper bounds exist, several studies [25,53] has hinted that much sharper upper bounds that depend the topology of underlying power networks may exist. The main goal of this paper is to establish such a topology dependent solution bound. The remainder of the paper is organized as follows: In Sec. 2, we set up a formulation of the algebraic load flow equations, define the problems precisely, and provide a brief historical account on the available results. Sec. 3 describes our rigorous results: in Sec. 3.1 we describe the tight bound, called the BKK bound, on the number of isolated complex solutions for our systems. We also refine the BKK bound in the context of algebraic load flow systems. In Sec. 3.2 we introduce a novel upper bound called adjacency polytope (AP) bound. Sec. 3.3 discusses the computational issues, and in Sec. 3.4 we highlight the significant implications of the development of these solution bounds in homotopy methods for solving load flow equations. In Sec. 4, we compute the BKK and AP bounds and the actual number of complex solutions for numerous power network examples from literature and compare with the previously known bounds. In Sec. 6, we discussion our results in a larger context and conclude.
2. Algebraic formulation. In this paper, we focus on the mathematical abstraction of a power network which is captured by a graph G = (B, E) together with a complex matrix Y = (Y ij ). Here B is the finite set of nodes representing the "buses", E is the set of edges (also called lines or branches) representing the connections between buses, and the matrix Y is the nodal admittance matrix 1 which assigns a nonzero complex value Y ij (mutual admittances) to each edge in (i, j) ∈ E. For any (i, j) / ∈ E (i.e., nodes i and j are not directly connected), Y ij = Y ji = 0. Here, Y is not assumed to be symmetric, but we require Y ij and Y ji to be either both zero or both nonzero (i.e., the underlying graph is "undirected"). As a convention, we further require all nodes to be connected with itself via a "loop" to reflect the nonzero diagonal entries Y ii known as the self-admittances. That is, we require (i, i) ∈ E for any i ∈ B.
For brevity, we define n to be the number of non-reference buses (i.e., |B| = n + 1) and label the nodes by integers 0, 1, 2, . . . , n. The corresponding complex valued voltages will be denoted by v 0 , v 1 , . . . , v n . Here we fix node 0 to be the designated reference bus for which the voltage is fixed to a nonzero real number, that is, v 0 ∈ R is a nonzero constant. In this setup, the load flow equation takes the form of which is a system of n equations in the n variables v 1 , . . . , v n since v 0 , corresponding to the reference bus, is fixed (and hence not a variable). Here v * i and Y * ij denotes the complex conjugates of v i and Y ij respectively, and S i ∈ C are the injected power. The equations 1 may represent either a transmission or a distribution network, with PQ buses. It is the network topology along with other features that can distinguish these cases: a mesh topology would usually correspond to transmission networks, whereas radial (tree-like) topology would correspond to distribution networks.
A solution to (1) is said to be isolated if it is the only solution in a sufficiently small neighborhood. i.e., an isolated solution has no degree of freedom. Solutions with some v k = 0 are said to be deficient, and non-deficient otherwise. By a standard application of Sard's Theorem [60], it can be verified that under a generic perturbation of S 1 , . . . , S n , the system (1) has no deficient solutions (though not a completely unlikely case in some load flow systems). We therefore focus only on the non-deficient solutions. The problem central to this paper is a counting problem for the isolated non-deficient load flow solutions: Problem Statement 1. For a power network, what is the maximum number of isolated non-deficient solutions to the system (1)?
Following the fruitful algebraic approach taken by works such as [3,42], we shall "embed" this problem into a bigger algebraic root counting problem: That is, we consider a polynomial system whose solution set captures all the solutions of the above (non-algebraic) system. To that end, we introduce new variables Substituting them into (1), we obtain the system of algebraic equations This is a system of n equations in 2n variables. However, a "square" system where the number of variables and equations match is often much more convenient from an algebraic point of view. We therefore extract n hidden equations by taking the complex conjugate of both sides of each of the above equations and obtain We now sever the tie between u = (u 1 , . . . , u n ) and v = (v 1 , . . . , v n ) and consider them to be variables independent from one another. That is, we ignore (2). Then (3) and (4) combine into a system of 2n from Y .
2 polynomial equations in the 2n variables: Note that here the values of v 0 and u 0 are fixed, as they correspond to the reference node and are hence constants in the above system. For brevity, this system will be referred to as the algebraic load flow equations in the following discussions. This formulation first appeared in [42] to the best of our knowledge.
Other polynomial formulations of the load flow equations have been also known (see, e.g., [2,3,8,9,46,49]). It is worth noting that in P G,Y , the topology (i.e., the edges) of the underlying power network, in a sense, is encoded in the set of monomials while entries of Y appear as the coefficients. Developing a solution count that exploits network topology via the monomial structure is the main goal of the current article.
Clearly, for every solution v = (v 1 , . . . , v n ) of the original (non-algebraic) system (1), we naturally have P G,Y (v, v * ) = 0. That is, the solution set of P G,Y = 0 captures all solutions of the original (nonalgebraic) load flow system (1). In the following we shall focus on the following root counting problem: Problem Statement 2. Among the power networks of a given topology provided by a graph G, what is the maximum number of isolated solutions of P G,Y in (C \ {0}) 2n for all possible choices of Y ?
Here, the phrase "maximum number" means the lowest upper bound that is also attainable and shall be distinguished from a mere "upper bound". Of course, the existence of such a "maximum number" is not a priorily guaranteed, after all the lowest upper bound taken over an infinite family (choices of Y 's), even if exists, may not be attainable. One of the goals of this paper is to establish the validity of the above question. Clearly, any answer to Problem 2 provides an upper bound for the answer to Problem 1.
Various upper bounds for Problem 2 have been proposed in the past (see [53] for a recent review). The classical Bézout number (CB number) provides a simple upper bound: it is the product of the degrees of each polynomial equation in (5) (hence also known as the total degree). It is a basic fact in algebraic geometry that the number of isolated complex solutions of a polynomial system is bounded above by the CB number. Therefore, for a power network of n (non-reference) buses, and one reference bus, the CB bound is 2 2n , since there are 2n equations in (5) each of degree 2. A much tighter upper bound on the number of isolated complex solutions, 2n n , was derived for the special case of completely interconnected lossless networks with all nodes being power-voltage nodes in Ref. [3] and then for the general case in Ref. [42] (see [46] for a recent alternative derivation of this bound). We shall refer to this bound as the Baillieul-Byrne-Li-Sauer-Yorke (BBLSY) bound. However, neither of these bounds exploit the network topology of a given power system. The link between network topology and complex solution count was first hinted in [25], however, a concrete and computable answer remains elusive.
In a recent study [53], with extensive numerical experiments via the NPHC methods, it was observed that the number of isolated complex solutions is generally significantly lower than both the CB and BBLSY bound for sparsely connected graphs. Based on these observations, it was anticipated that the key to exploiting the network structure of the power system may be to exploit the underlying topology of the power system. In the present work, we show that this maximum number exists and it is essentially the Bernshtein-Kushnirenko-Khovanskii (or BKK) bound. We then develop a novel approximation of this maximum number, the "adjacency polytope" (or AP) bound, which has tremendous computational advantage over the BKK bound. Moreover, we show, via experiments, that the AP bound is exactly the BKK bound for many concrete cases, making it an ideal surrogate for the BKK bound.
3. The Maximum Number of Solutions. We shall focus on Problem 2 and show that the BKK bound provides the best possible bound for any given power network topology which is always attainable by some choice of the Y -matrix. Then in §3.2 we propose a simplified formulation that is easier to compute and analyze. Computational issues are explored in §3.3.
3.1. The BKK bound. Problem 2 is a special case of the root counting problem for systems of polynomial equations which is an important problem in algebraic geometry and mathematics in general that has wide range of applications [55,56,60]. Two basic root counts are provided by the CB and BBLSY bounds described above. One common weakness of these two upper bounds is that they only utilize the rather incomplete information about the polynomial system -the degree (or "multi-degree" of the polynomials). In the context of the algebraic load flow equation, this means that they do not take into consideration the topology of the underlying power network. Following up from the observations in Ref. [53] we refine these bounds using the theory of BKK bound [4] which accurately captures the network topology of the power systems. Here we state the powerful, albeit abstract, Bernshtein's theorem in the context of the algebraic load flow equation. The technical terms are explained in Appendix A for completeness. Here, it is sufficient to take the following interpretation: Part (A) establishes a computable upper bound for the number of isolated solutions that depends on the geometric configuration of the monomial structure, and part (B) shows this upper bound is generically exact. 3 The original proof was given in [4]. An alternative proof that gives rise to the development of polyhedral homotopy is developed in [28]. More detail can be found in references such as [15,29,41]. In [7], the root count in the above theorem was nicknamed the BKK bound 4 after the works of Bernshtein [4], Kushnirenko [37], and Khovanskii [33]. In general, it provides a much tighter bound on the number of isolated zeros of a polynomial system compared to variants of Bézout bounds. More importantly, it is sensitive to the monomial structure and hence the topology of the underlying power network.
It is important to note that the "generic exactness" expressed in part (B) of the above theorem only holds when one ignores the tie between Y ij and Y * ij as well that between S i and S * i . That is, one must allow Y ij and Y * ij to vary independently in interpreting the above theorem. We shall now bring back the restriction that all the (Y ij , Y * ij ) and (S i , S * i ) must be conjugate pairs and investigate the sharpness of the BKK bound under these restriction. Since the BKK bound is sensitive to the monomial structure in the equation (5) and hence topology of the underlying power network, we shall fix the network topology in the following discussion. That is, we shall fix the sparsity pattern of the Y matrix but allow its entries to vary among the set of nonzero complex numbers. In the following, we shall establish that within the set of admittance matrices of the same sparsity pattern: The BKK bound of the algebraic load flow equations is always exact for some choices of the admittance matrix. In other words, we have the following assertions: Theorem 2. Given a graph G, there exists an admittance matrix Y on G for which the number of isolated solutions of the corresponding algebraic load flow equation Proof. For convenience, let Z = (Z 1 , . . . , Z ) which collect all the nonzero entries of Y jk and S j . That is, Z contains all the nonzero coefficients in (5). By the interpretation given in Appendix B, we simply have to show that there exists a choice of Z ∈ (C \ {0}) such that D(Z, Z * ) = 0.
Suppose no such choice of Z exist, then D(Z, Z * ) = 0 for all Z ∈ (C\{0}) . By Lemma 7 in Appendix C, D(Z, W ) = 0 for all (Z, W ) ∈ (C\{0}) 2 . This means D must be a zero polynomial, which is a contradiction. Therefore, we must conclude that there is always a choice of Z (and hence Y and S) such that D(Z, Z * ) = 0.
3.2. Solution bound via adjacency polytope. In this section, we develop a significantly simplified version of the BKK bound for the algebraic load flow equation (5) that can be analyzed and computed more easily. Here we borrow a construction from coding theory and encode the given graph into a polytope (which, roughly speaking, is a geometrical object whose all sides are flat) which we shall call the "symmetric AP". The definition requires the following notations: We define e 0 := 0 ∈ R n . For i = 1, . . . , n, e i denotes the vector in R n that has an entry 1 on the i-th position and zero everywhere else. (e i , e j ) is simply the concatenation of the two vectors e i ∈ R n and e j ∈ R n that forms a vector in R 2n . Finally, "conv" denotes the convex hull operator which produces the smallest convex set containing a given set.
With this, we define the symmetric adjacency polytope of G to be

Implicit in this definition is the assumption that the graph
The polytope ∇ G is a geometric encoding of the connectivity of the underlying power network with connections manifested as points. It is only dependent on the connectivity among the nodes (and hence the sparsity patterns in Y ) but independent from the actual entries in Y .

Remark 1.
It is evident in the construction of (5) that equations in the system always contain many common monomials. Indeed, if (i, j) is an edge in the graph, then the monomial v i u j appear in both the i-th and the (n + j)-th equation. In other words, the monomial structure of (5) has certain level of built-in redundancy. Such redundancy is removed in the construction of ∇ G : In Definition 3, one takes the union of the set of points representing the edges, common monomials in (5) will therefore coalesce into the same points. Consequently, the polytope ∇ G , in a sense, contains much less information than the monomial structure in (5). Therefore the encoding ∇ G is advantageous from a computational point of view.
With the constructions above, we propose a new upper bound on the isolated non-deficient complex solutions to the algebraic load flow equation that takes into consideration the connectivity of the underlying power network: Theorem 4. The number of isolated solutions the algebraic load flow system (5) has in (C \ {0}) 2n is bounded above by which we shall define as the the adjacency polytope bound (or simply AP bound) for any P G,Y (with any choice of admittance matrix Y ).
Here "NVol 2n " denotes the normalized volume in R 2n . It is a volume measurement commonly used in the study of "lattice polytopes", and it is defined so that the standard "corner simplex" (the corner of a unit hypercube) has volume 1.
Proof. For a nonsingular 2n × 2n matrix M , we can form the new system M · P G,Y by forming the formal matrix-vector product where P G,Y is considered as a column vector. This technique is known as randomization. Clearly, M · P G,Y (v, u) = 0 if and only if P G,Y (v, u) = 0 and the number of isolated solutions (in (C \ {0}) 2n ) remains the same under this transformation. It is easy to verify that the support of the randomized system M · P G,Y is unmixed of type 2n, and the Newton polytope is precisely the symmetric adjacency polytope ∇ G defined in Definition 3. Then by the unmixed form of Bernshtein's Theorem [28], the BKK bound of this randomized system is precisely the normalized volume NVol 2n (∇ G ). Example 1. Consider, for example, the simple "path graph" G = ({0, 1, 2, 3}, E) of 4 nodes as shown in Figure 2 where each node i is connected the next node i + 1. Recall that we also require each node to have a loop to itself (to reflect the nonzero diagonal entries of the nodal admittance matrix Y ), so the edges in the graph are Since 0 is already included in Γ G , the symmetric adjacency polytope ∇ G is precisely the convex hull of Γ G . With programs for computing volume of convex polytopes to be listed in §3.3, one can easily obtain that µ G = NVol 6 (∇ G ) = NVol 6 (conv Γ G ) = 8. That is, the algebraic load flow equation for such a path graph has at most 8 isolated non-deficient complex solutions.
It is quite easy to understand the situation under which the introduction of a new connection in the power network will change the AP bound: Since the AP bound is formulated in terms of the volume of a polytope which is nondecreasing (i.e., it will either increase or remain unchanged when new points are added), this upper bound must also be nondecreasing with respect to the additions of new connections.
Theorem 5. For a graph G = (B, E) and two of its nodes i and j that are not directly connected (i.e., be the new graph constructed by adding the edge between i and j to G. Then Proof. Recall that each edge in a graph contributes certain points (which may or may not be vertices) in the construction of the symmetric adjacency polytope. Since the edges of G is a subset of the edges of G , we can see that ∇ G ⊆ ∇ G with the equality hold precisely when the points contributed by (i, j) are already inside ∇ G . With these observations in mind, both parts of this theorem are direct implications of normalized volume being nondecreasing.
Based on this observation, it can be shown that the AP bound is never more than the BBLSY bound. This is essentially our alternative proof of the BBLSY bound: Theorem 6. For a graph G = (B, E), let |B| = n + 1. Then, µ G ≤ 2n n Proof. Fixing the set of buses, Theorem 5 states that the AP bound is nondecreasing as new edges are added to the graph. Consequently, the AP bound for any network constructed from this set of buses is bounded above by the AP bound for the graph with most edges, that is, a complete graph. It is easy to verify that for a complete graph G = (B, E) (with loops), ∇ G ⊆ (conv A) + (conv B) where  (5) in (C \ {0}) 2n , the BKK bound provided by Theorem 1, the AP bound given in Theorem 4, the BBLSY bound given in [3,42], and the CB number are related as follows: 3.3. Computing BKK and AP bounds. The BKK bound, can be computed by using efficient software packages such as DEMiCs [51], Gfan [32], MixedVol [22], MixedVol-2.0 [39]. For larger power networks involving many buses, the induced algebraic load flow equation may contain a large number of terms. Parallel computing technology will be essential in the computation of the BKK bounds for such large scale polynomial systems. MixedVol-3 [13,14] (with an improved version integrated in Hom4PS-3 [12]) is capable of computing the BKK bound for large polynomial systems in parallel on a wide range of hardware architectures including multi-core systems, NUMA systems, and computer clusters. As noted in Remark 1, however, there is built-in level of redundancy in the Newton polytopes (see Appendix A) of the algebraic load flow equations. The formulation of the AP bound takes advantage of this natural redundancy and can generally be computed much more easily than the BKK bound for larger power networks. Since the AP bound is formulated in terms of the volume of a convex polytope (the symmetric AP), any software that can compute such volume exactly can be used to provide this bound. A survey on the various algorithms for exact volume computation can be found in [6]. A particularly versatile software package that utilizes a wide range of algorithms is the Vinci program [21].
3.4. Homotopy methods for solving load flow equations. The previous sections described two upper bounds -the BKK bound and the AP bound -for the number of isolated non-deficient complex solutions of the algebraic load flow equation induced by a given network topology. It is worth reiterating that the BKK bound is more than just an upper bound : As shown in Theorem 2, it is actually the maximum for the given network topology in the sense that there always exits some choice of the admittance matrix Y (and S) for which the total number of isolated non-deficient complex solutions is exactly the BKK bound. While the AP bound, in general, may be larger than the BKK bound, as we shall show in §4, the two coincide for all the networks we have investigated in the present work. The family of numerical methods known as homotopy methods have been proved to be a robust and efficient approach for solving algebraic load flow equations. One great strength of these methods lies in the pleasantly parallel nature: in principle, each solution can be computed independently. This feature is of particular importance in dealing with larger power networks containing many buses and connections (hence more complicated equations). It is therefore a natural question to ask: is there a homotopy method that can solve the algebraic load flow equation by tracking BKK bound number of homotopy paths? This section establishes the answer in the affirmative.
This homotopy method is the polyhedral homotopy method developed in [28]. Here we briefly state the construction. We choose a pair of random (rational) numbers ω ij and ω ij for each i, j = 1, . . . , n. With these we define the homotopy function An apparent difficulty is the identifications of the "starting points". After all, at t = 0, H G,Y (u, v, t) either becomes constant or undefined since each of its nonconstant term contains a power of t. This, fortunately, is surmounted via a device known as mixed cells which are themselves the by-product of the computation of the BKK bound. The technical detail is outside the scope of the present article, we therefore refer to standard references such as [28,41]. The application of the polyhedral homotopy to load flow equations will be explored in future works, here we simply emphasize that with the polyhedral homotopy method, the number of solution paths one needs to track is precisely the BKK bound of the particular polynomial formulation of the load flow equations that we have used (eq. (5)) -the best possible network topology-based complex solution bound.

Solution bound for certain power networks.
In this section we provide concrete computation results for the BKK bound and the AP bound described above when applied to specific graphs. We should reiterate that we require all the graphs to have self-loops for each node to reflect the nonzero diagonal entries of the nodal admittance matrix Y . In all cases, the BKK and the AP bounds are computed via MixedVol-3 [13,14]. Complex solutions count of specific load flow systems are computed by actually solving the systems via Hom4PS-3 [11,12] which implements the polyhedral homotopy method described above.
For comparison, in each case we show the "5-way comparison" among the complex solution count, BKK bound, AP bound, the BBLSY Bound, and the CB number. Here, "complex solution count" refers to the number of isolated complex solutions of (5) in (C \ {0}) 2n where n is number of non-reference nodes. Since the actual complex solution count may depend on the coefficients (Y ij and S i in (5)), they are computed based on a sample of randomly chosen set of coefficients for each graph.

Paths.
We first consider an extremely sparsely connected family of graphs -path graphs. In a path graph G = (B, E) with nodes 0, 1, . . . , n where n = |B| − 1 (the number of non-reference nodes), a node i is connected to two of its neighbors i + 1 and i − 1, with the exceptions that node 0 is only connected to node 1 and the last node, node n, is only connected to node n − 1. Table 1 shows the "5-way comparison" described above. It is important to note that in all cases computed (100 in total, with 10 random admittance matrices chosen for each |B|), the actual complex solution count, the BKK bound, and the AP bound are exactly the same.

Rings.
By joining the two ends of a path graph studied above, a "ring graph" is obtained. More precisely, in a ring G = (B, E) with nodes 0, 1, . . . , n where n = |B| − 1, a node i is connected to two of its neighbors i + 1 and i − 1, and node reference node, node 0 is connected to node n. See, for example, the graph shown in Figure 3a. Table 2 shows the similar "5-way comparison". Again, in all cases computed (100 in total, with 10 random admittance matrices chosen for each |B|), the actual complex solution count, the BKK bound, and the AP bound are all exactly the same.
(a) A ring.

Clusters.
Real power networks generally exhibit certain level of "clustering", that is, certain subset of buses are densely connected while on a larger scale, the connections among such subsets are sparse. Here for simplicity, we focus on the most extreme cases where a larger network is created by joining completely connected subnetworks. This section lists concrete results of the BKK and AP bounds for certain simple types of clusters. For comparison, in each case we only show the "3-way comparison" among the BKK bound, AP bound, and the BBLSY bound (due to the large amount of data).  4.4.1. Subnetworks sharing nodes. First, these subnetworks can be jointed by sharing common nodes. See, for example, the networks shown in Figure 3. Table 4 shows the "3-way comparison" for cases where two completely connected subnetworks share a single (non-reference) common bus. These cases have been studied in [26]. Our computational results agree with their assertion. Table 5 shows the similar "3-way comparison" for cases where two completely connected subnetworks share two (non-reference) common bus. These cases have been extensively studied in [53] via numerical methods. The results and conjectures in that work are precisely reproduced by our computation. For larger networks, the AP bounds are generally much easier to compute than the BKK bound using existing implementations. In Table 6 Table 9 shows the AP bounds of the more general cases where the networks consist of multiple completely connected subnetworks of the same sizes connected via edges to form chains-like structure. See, for example, the network shown in Figure 4b where five cliques each of size three are connected via edges that, on a macro level, resembles a chain.   8  1  2  4  8  16  32  64  128  2  2  8  32  128  512  2048  8192  32768  3  6  72  864  10368  124416  1492992  17915904  4  20  800  32000  1280000  5  70  9800  1372000  6  252  127008  7  924  1707552  8 3432 Table 9: The AP bound for graphs consisting of n cliques each of size c. Here we show that the isolated complex solution count, BKK bound, and AP bound are much smaller than previously studied solution bounds. In particular, the BKK bound in our formulation of the load flow equations is 427680. This means the polyhedral homotopy method described in Sec. 3.4 need to trace at most 427680 paths to obtain all isolated nondeficient complex solutions, which compared with the BKK bound in a previous formulation of the load flow equations [49], 49283072, is around a factor 115 of reduction. This is well within the reach of modern polyhedral homotopy implementations. In particular, with a random choice of the Y -matrix, Hom4PS-3 [11,12] was able to find all solutions in less than 5 minutes (around 297 seconds) on a single workstation (with two Intel Xeon processor).  6. Conclusion and Outlook. Solutions of the load flow equations play an important role in design and operations of power systems and a huge amount of research activities have been focused on developing efficient computational methods to solve the load flow equations [48]. In this paper we focus on computing a tight upper bound on the number of load flow solutions for the following reasons: the NPHC method [41,60] introduced in power systems areas in [8,9,47,49,53] requires an upper bound on the number of isolated (complex) solutions of the system. The CB bound or most of the other existing upper bounds, and in particular the BBLSY bound [3,42], do not capture the network topology [53] of the power system into account. To make the computation of the NPHC method efficient, and to provide a definitive stopping criterion for other iterative methods, one has to come up with a tighter upper bound on the number of isolated complex solutions.
In this paper, for technical reasons, we ignore the physically less likely solutions at which at least one of the voltages is zero. We described a specific version of the algebraic formulation of the load flow equations and a corresponding tighter upper bound known as the BKK bound. We showed that there exists at least some generic parameter values for which the BKK bound is attainable.
Another novel feature in this paper is the introduction of a novel bound, called an AP, for sparse systems. With this bound, we proved one of our main theorems, Theorem 4, which gives a more explicit bound on the number of complex solutions. The new bound is computationally much more tractable for large systems compared to the BKK bound.
Our bounds not only give the precise number of complex solutions for the given power system, but also provide a particular specialization of the NPHC method that in turn achieves all the solutions for the system but more efficiently than the previous ones that are based on the crude bounds. We examined our results on upper bounds on the number of complex solutions with the actual number of complex solutions obtained with the help of the polyhedral NPHC method for various cases including complete graphs; path graphs; ring graphs; wheel graphs; random tree graphs; graphs having two cliques sharing exactly one, two and three nodes; graphs having two cliques both connected by a single edge; graphs with multiple cliques of equal size connected via edges, and the IEEE standard 14-bus system. In all the instances we have tested, we find that (1) the BKK bound is always equal to the adjacency bound, and (2) the BKK bound is always equal to the actual number of complex solutions, i.e., the BKK or the AP bound are the actual number of isolated complex solutions for generic coefficients and completely capture the specific network topology of the given power systems. A proof supporting the numerical results will be published elsewhere. Our numerical investigations reproduced all the previously known results [3,27,53]. In addition, we provide further novel results for aforementioned list of graphs.
We emphasis that the BKK bound is dependent on the polynomial formulation used to compute it, though the number of complex solutions and number of physical solutions are obviously independent of the polynomial formulations. The phrase "BKK bound" throughout the paper, unless specified otherwise, is computed from the polynomial formulation we have described in this paper (given in eq. (5)). In other words, our polynomial formulation seems optimal since the BKK bound computed using it is generically equal to the number of isolated complex solutions. To give an example, for the IEEE standard 14-bus system, the BKK using a previous formulation was 49, 283, 072, whereas the BKK in the present formulation is 427, 680, i.e., a factor of around 115 reduction with our present formulation.
Our results open up several future directions: though a less likely case, extending our results to the situtation where at least one of the voltages is zero would further complete the discussion on number of complex solutions. Further refinement will be accomplished by computing or estimating the number of physical solutions, and analysing their stability properties for different network topologies. Investigating for which parameter values the number of physical solutions is equal or at least of the same order of magnitude as the number of complex solutions of the systems will also provide a better tool to analyse efficiencies of different computational methods: For example, for the 3 nodes complete graph case, it is known that the number of real solutions does attain the tight-most upper bound on the number of complex solutions, 6 [2,35,61]. However, for 4 nodes complete graph case, so far parameter points with only 14 physical solutions are known as opposed to 20 complex solutions [3]. This has an important consequences on the tracing-methods for solving load flow equations such as proposed in Refs. [40,44]. Leaving aside the lack of proof for guaranteeing all the physical solutions of the load flow equations in the tracing-methods, for such parameter points at which the number of physical solutions is equal or of the same order of magnitude as the BKK bound, the number of traces the method has to go through will be equivalent to the polyhedral NPHC method will also have to go through. Hence, in such cases, there will be no gain even if one ignores the guarantee to find all solutions that the NPHC method provides. On the other hand, our bounds can provide