Locating Power Flow Solution Space Boundaries: A Numerical Polynomial Homotopy Approach
11 Locating Power Flow Solution Space Boundaries: ANumerical Polynomial Homotopy Approach
Souvik Chandra, Dhagash Mehta, Aranya Chakrabortty
Abstract
The solution space of any set of power flow equations may contain different number of real-valued solutions. The boundariesthat separate these regions are referred to as power flow solution space boundaries. Knowledge of these boundaries is importantas they provide a measure for voltage stability. Traditionally, continuation based methods have been employed to compute theseboundaries on the basis of initial guesses for the solution. However, with rapid growth of renewable energy sources these boundarieswill be increasingly affected by variable parameters such as penetration levels, locations of the renewable sources, and voltageset-points, making it difficult to generate an initial guess that can guarantee all feasible solutions for the power flow problem. Inthis paper we solve this problem by applying a numerical polynomial homotopy based continuation method. The proposed methodguarantees to find all solution boundaries within a given parameter space up to a chosen level of discretization, independent ofany initial guess. Power system operators can use this computational tool conveniently to plan the penetration levels of renewablesources at different buses. We illustrate the proposed method through simulations on 3-bus and 10-bus power system exampleswith renewable generation.
Index Terms
Power flow solution boundary, renewable energy sources, numerical homotopy, voltage stability
I. I
NTRODUCTION
Over the past few decades, power system networks in different parts of the world have been going through transformationalchanges in structure and operation due to integration of large amounts of renewable resources with associated new devicessuch as power electronic converters, loads, smart transformers, and new transmission lines [1]. Addition of all of these newdevices has tremendously facilitated power system operations and markets, but at the same time has also provoked threats onpower system stability [2], especially voltage stability [3]–[5]. Stimulated by several major voltage collapses, the frameworkof voltage stability as defined in [6] has been extensively studied since the 1990s. Voltage stability can be classified intosmall disturbance and large disturbance categories. Small disturbance voltage stability refers to the ability of a power systemto maintain steady voltage levels following small disturbances experienced through continuous changes in load [7]. But withexponentially increasing renewable generation, operational uncertainties and variability not only arise from loads but fromgenerations as well, directly affecting the dynamic performance and voltage fluctuations in the grid. Typically, the transientvoltage stability limit of the grid is characterized by a power flow solution space boundary or a loadability boundary at whichthe Jacobian matrix of the power flow equations is singular [2]. System operators choose operating points that are well withinthe loadability boundary to reduce risk of disruption of service by balancing off deficit powers through reserves. But withhigh penetration of renewables and their associated intermittency this type of balancing will become more difficult from bothoperational standpoint and economic viability [8]. Therefore, it will be critical for operators to have complete information ofpower flow solution boundaries in terms of the injections levels of the different renewable energy sources.In this paper we present a new numerical approach for computing all possible solution boundaries of a power flow problem asa function of a parameter set that models the amount of power injected into the grid from extraneous sources. We do not makeany specific distinction between the source of these extraneous injections, meaning that these injections can happen throughnew renewable installations such as wind and solar generators as well as through new conventional synchronous generators.Several papers have been written in recent past to address similar problems of finding the power flow solution boundaryin a certain parameter space [9]–[11]. These methods identify the trajectory of the solution boundaries by continuation-basedtracking, starting from an initial guess. The work that is most relevant to our problem is by Hiskens et al. in [11] that proposes agradient based predictor-corrector method to identify the power flow solution boundary in a given multi-dimensional parameterspace. However, the success of this method is critically dependent on two factors, namely, a good initial guess and continuity ofthe solution boundaries, both of which may be difficult to achieve in face of variability due to intermittent power injection from
S. Chandra was with the Department of Electrical and Computer Engineering, North Carolina State University, Raleigh, NC, USA, when the work in thispaper was carried out. He is with Schweitzer Engineering Laboratories, Pullman, WA, USA, e-mail: [email protected]. Mehta was with the Department of Applied and Computational Mathematics and Statistics, University of Notre Dame, Notre Dame, IN, USA, andCentre for the Subatomic Structure of Matter, Department of Physics, School of Physical Sciences, University of Adelaide, Adelaide, South Australia 5005,Australia, when the work in this paper was carried out. He is with the Systems Department at United Technologies Research Center, East Hartford, CT, USA.A. Chakrabortty is with the Department of Electrical and Computer Engineering, North Carolina State University, Raleigh, NC, USA, e-mail:[email protected] was supported by the NSF through grant NSF-ECCS-1509036 and an Australian Research Council DECRA fellowship no. DE140100867.Chakrabortty was supported partially through NSF grant ECCS 1509137. a r X i v : . [ c s . S Y ] A p r renewables. For example, the solution boundary may be non-smooth due to inequality constraints in the power flow problemoccurring due to generator over-excitation limits or reactive power limits. [12], [13]. Moreover, depending on different levelsof penetration and their relative locations with respect to the existing generators and loads, the structure of these boundariesmay potentially change. Therefore, continuity-based methods such as those in [11] may not be able to find all solution spaceboundaries over a given parameter space.In contrast, our approach in this paper is based on a numerical polynomial homotopy-based continuation method (NPHC).This approach guarantees to find all solution boundaries within a parameter space, and can, therefore, help an operator to makea more thorough evaluation of voltage stability limits. The principal distinction of our algorithm from [11] is that our methoddoes not look for the singularity of the Jacobian of the power flow equations. Rather we compute all the isolated real solutionsof the power flow equations directly at different points on the parameter space. At the points on the solution boundary, whichare also denoted as bifurcation points, the number of real-valued solutions of the power flow problem changes. Therefore,one only needs to identify those points in the parameter space. The advantage is that one can obtain all possible equilibriafor a given parameter set, not just a local equilibrium that is closest to the initial guess [14]–[17]. NPHC has recently beendemonstrated on IEEE prototype power system models in [18]. The algorithm in this paper is inspired by a variant of NPHCknown as parameter-coefficient NPHC, reported in [19]–[21].The remainder of this paper is organized as follows. In Section II, we introduce the problem statement with the formulationof the power system equilibrium model and a discussion on solving the power flow boundaries via the traditional continuation-based methods as in [11]. In Section III we solve the power flow equations using the NPHC approach, and highlight on theuse of power system network topology to reduce computation. In Section V we present a discussion on a novel choice of theupper bound for the number of solutions of NPHC utilizing the structure of the power system model. Section V also presentscase studies on a 3-bus power system with variable active power injection at one bus, and a 10-bus example with two variablerenewable power generators. Section VI concludes the paper with discussions and future directions of research.II. P OWER FLOW SOLUTION BOUNDARIES IN POWER SYSTEMS WITH RENEWABLE GENERATION
A. Problem formulation
We consider a power system with N buses and n generators. These generators are either classified as synchronous generatorswith indices belonging to a set S =: { , . . . , n s } , or generators with renewable sources with indices belonging to a set R =: { , . . . , n r } , where n s + n r = n . Without loss of generality we classify the buses into 3 sets as follows: Bus to n s arethe synchronous generator buses whose indices belong to a set N s =: { , . . . , n s } , bus n s + 1 to n s + n r are the renewablegeneration buses whose indices belong to a set N r =: { n s + 1 , . . . , n } and the rest of the variables are for load buses with theirindices belonging to the set N l =: { n + 1 , . . . , N } . The equilibrium for each of the different buses are obtained by consideringpower balance with the neighboring buses [22]. The superscript e for any variable is used to indicate its equilibrium value(s).Usually in a power flow problem, one of the synchronous generator bus is considered to be the slack bus, which in this caseis assumed to be bus without any loss of generality. Thus the reference voltage level | V e | and angle θ e are respectively equalto 1 and 0. Correspondingly, the steady-state active and reactive power flow, P e and Q e from slack bus can be expressedas, P e =Re (cid:40) N (cid:88) k =2 (cid:18) (1 − V ek ) Z k (cid:19) ∗ (cid:41) + P eL , (1a) Q e =Im (cid:40) N (cid:88) k =2 (cid:18) (1 − V ek ) Z k (cid:19) ∗ (cid:41) + Q eL , (1b)where k denotes the index of neighboring buses, V k is the voltage of the bus k and Z k is the line impedance connecting buses and k . P L and Q L are the active and reactive power of the loads connected to bus . All other synchronous generatorbuses, i ∈ N s , such that i (cid:54) = 1 are considered to be PV buses for which the active power P ei and | V ei | are specified. The powerbalance for bus i with the neighboring buses is shown as below, P ei =Re N (cid:88) k =1 ,k (cid:54) = i V ei (cid:18) V ei − V ek Z ik (cid:19) ∗ + P eLi (2a) Q ei =Im N (cid:88) k =1 ,k (cid:54) = i V ei (cid:18) V ei − V ek Z ik (cid:19) ∗ + Q eLi . (2b)The bus j ∈ N r , connected to a renewable source, may be either PV or PQ depending on the type of the renewable generators.But in this work, for the sake of simplicity we limit our model to only those types of renewable generators for which internal control of both active and reactive power at the output bus is available. Thus, we consider the renewable bus to be a PV buswhich maintains a given reference voltage (cid:12)(cid:12) V ej (cid:12)(cid:12) and injects active power P ej . The power balance is shown as below, P ej =Re N (cid:88) k =1 ,k (cid:54) = j V ej (cid:18) V ej − V ek Z jk (cid:19) ∗ + P eLj (3a) Q ej =Im N (cid:88) k =1 ,k (cid:54) = j V ej (cid:18) V ej − V ek Z jk (cid:19) ∗ + Q eLj . (3b)The total active and reactive power injections at these buses are assumed to be a product of a factor λ j and constant unitpowers P er and Q er respectively, where j ∈ N r . Here, λ j , for example, may represent the total number of wind turbines behindthe point of common coupling with the grid, assuming every turbine produces the same amount of power in equilibrium.In general, λ j is an indicator of the total amount of injection at bus j , and will be treated as a variable in the power flowequations. The variable parameter λ j can also be incorporated in (2) for an increase in synchronous generation resulting intransmission expansion. For convenience of analysis we limit the use of λ j to (3) to introduce variability only through therenewable sources. Correspondingly, equation (3) can be rewritten as, λ j P er =Re N (cid:88) k =1 ,k (cid:54) = j V ej (cid:18) V ej − V ek Z jk (cid:19) ∗ + P eLj (4a) λ j P er =Im N (cid:88) k =1 ,k (cid:54) = j V ej (cid:18) V ej − V ek Z jk (cid:19) ∗ + Q eLj , (4b)For load bus, j ∈ N l , usually the active and reactive load power P eLj , Q eLj are specified. The power balance is shown as below, N (cid:88) k =1 ,k (cid:54) = j V ej (cid:18) V ej − − V ek Z jk (cid:19) ∗ + P eLj (5a) N (cid:88) k =1 ,k (cid:54) = j V ej (cid:18) V ej − V ek Z jk (cid:19) ∗ + Q eLj . (5b)The complete set of power flow equations for a power system with renewable injections is given by the set of nonlinearalgebraic equations shown in (1)-(5). These equations are parametrized by λ ∈ R n r , the vector of all renewable penetrationlevels. (1)-(5) can be represented as a set of equations P as shown below, P ( x, λ ) = 0 . (6)where x ∈ R N is a vector of unknown power flow variables. Equation (6) for the slack bus is solved for x which includes, P e and Q e for slack bus, Q ei and θ ei where i ∈ ( N s ∪ N r ) , i (cid:54) = 1 for the PV buses, θ ej and (cid:12)(cid:12) V ej (cid:12)(cid:12) where j ∈ N l for the PQbuses. B. Problem to be solved
The power flow solution boundary represents those solutions of the power flow equations (1)-(5) for which the Jacobian ofthe equations is singular. From (6), it is clear that the Jacobian is a function of renewable penetration levels λ . The objectiveof this work is to find out all boundary solutions of (1)-(5) for a given set of points in the parameter space of λ , and obtainthe corresponding solution boundaries. C. The traditional method
In [11] Hiskens et al. propose a technique based on a predictor-corrector method to identify the power flow solution boundaryof P defined in (6) over a parameter space formed by λ . The solution method consists of two steps. Initially the solution spaceboundaries are represented as a solution of the set of equations as shown below, P ( x, λ ) = 0 (7a) g ( x, λ, v ) = P x ( x, λ ) .v = 0 (7b) h ( v ) = v T v = 1 . (7c) Here the Jacobian matrix of P with respect to x is given by P x . v ∈ R N is a right eigenvector corresponding to a singulareigenvalue of P x . As a first step of the method, one of the parameters, say λ is varied while all others are kept constant.Consequently, (7) turns into a set of n + 1 equations with n + 1 unknowns which can be solved via computational methodssuch as Newton-Raphson with an initial guess. The solution of this first step, say x , will provide an initial point on thesolution space boundary. Starting from x , a gradient based predictor-corrector tracking is implemented to obtain the otherpoints on the boundary. First, the algebraic equations in (7) are represented as, φ ( z ) = P ( z ) g ( z ) h ( z ) , (8)where z = (cid:2) x v λ (cid:3) (cid:48) . An initial guess z is obtained which satisfies the condition φ ( z ) = 0 . Starting from z , thealgorithm iteratively solves for each z i by solving the following equations at each step i , φ ( z i ) = 0 , (9a) ( z i − z i − ) (cid:48) ν = (cid:15). (9b) (cid:15) is a small constant that determines the accuracy of the numerical approximation. ν , on the other hand, is a vector that istangential to the solution boundary at z i − . Repeated solutions of z i can provide closed boundary contours of the power flowequations. Since this method tracks the solution boundary iteratively, the choice of the initial point z is particularly importantfor the success of this method. As specified earlier, the boundaries can be often disjoint and non-smooth. In such scenariosthis method does not guarantee to obtain all solution boundaries. To solve this problem, we introduce a numerical homotopycontinuation method, detailed in the next section.III. P ARAMETER H OMOTOPY C ONTINUATION A LGORITHM
In this section we solve for all real equilibria of the power system model shown in (1)-(5) using the parameter-coefficientNPHC method [23], [24]. As a preconditioning, we represent the algebraic equations shown in (1)-(5) as multivariate polyno-mials. To achieve this, the voltage phasor V ej at any bus j ∈ N is expanded in terms of its magnitude (cid:12)(cid:12) V ej (cid:12)(cid:12) and angle θ ej , suchthat V ej = (cid:12)(cid:12) V ej (cid:12)(cid:12) cos θ ej + J (cid:12)(cid:12) V ej (cid:12)(cid:12) sin θ ej , J being the complex operator. Using such expansion, energy balance in (3) for instancecan be expressed as, P ej =Re N (cid:88) k =1 ,k (cid:54) = j S jk + P eLj (10a) Q ej =Im N (cid:88) k =1 ,k (cid:54) = j S jk + Q eLj , (10b)where power flow between bus j and k can be represented as, S jk = (cid:12)(cid:12)(cid:12) V ej (cid:12)(cid:12)(cid:12) (cid:16)(cid:12)(cid:12)(cid:12) V ej (cid:12)(cid:12)(cid:12) cos θ ej − (cid:12)(cid:12) V ek (cid:12)(cid:12) (cid:16) cos θ ej sin θ ek − J sin (cid:16) θ ej − θ ek (cid:17)(cid:17)(cid:17) arg (cid:16) Z ∗ jk (cid:17) . Z ∗ jk represents the complex conjugate of the line impedance connecting bus j and bus k . In the modified system of equations,power balance of slack bus 1 has P e and Q e as variables yielding two first-order equations. The power balance of PQ bus j have (cid:12)(cid:12) V ej (cid:12)(cid:12) cos θ ej and (cid:12)(cid:12) V ej (cid:12)(cid:12) sin θ ej as the variables with two second-order equations. Both equations are in quadratic polynomialforms, as indicated in (10). For a PV bus, however, the variables are Q ej and θ ej , due to which (10) no longer retains itspolynomial form for this type of a bus. To solve this problem, the unknowns are alternatively defined as Q ej , cos θ ej , and sin θ ej ,all of which are independent of each other. An extra constraint equation now needs to be added as cos θ ej + sin θ ej = 1 , (11)where j ∈ N s ∪ N r . As a consequence, n − more equations are added for the n − PV buses considered in our system.The resulting system turns into a set of N + n − quadratic equations in x and λ as, P ( x, λ ) = 0 , (12)where, x ∈ R (2 N + n − is the vector of unknowns, and λ ∈ R n r is a constant parameter vector. It is noted that the symbol forthe function P in (6) has been changed to P to indicate that the latter represents a set of quadratic equations in the variable x . In essence, both equations (6) and (12) are equivalent. Gen 1 Bus 1 Bus 2Bus 3 Gen 2
Gen 3 j j j Load 1 Load 2 (cid:14) j (cid:14) j (cid:74) (cid:74) (cid:68)(cid:32) (cid:145) V (cid:32) (cid:145) V (cid:68) (cid:32) (cid:145) V (a) Power system model (b) Number of real solutions of the power flow problem Fig. 1: 3-bus power system with two parametersAs a first step of the NPHC method, an upper bound of the number of complex isolated solution of (12) is determined.Next, a homotopy H ( x, t ) is defined as shown below, H ( x, t ) = η h (1 − t ) Q ( x ) + t P ( x ) , (13)where Q ( x ) is an arbitrary start system which is easily solvable. Q ( x ) is chosen in such a way that the number of isolatedsolutions of Q ( x ) = 0 is equal to the estimated upper bound of isolated solutions of (12). η h is a generic complex number and t is a continuous parameter varying from 0 to 1. Therefore, the solution set of H ( x, t ) = 0 for ≤ t ≤ actually consists ofa finite number of smooth paths parametrized by t ∈ [0 , . For a generic η h ∈ C , it is proven in [24] that each of the pathswill be well-behaved, i.e., either they will converge to H ( x,
1) = 0 , or will diverge to infinity. Hence, for a generic value of η h , the NPHC method guarantees to find all isolated complex solutions [25] of P ( x ) = 0 . The crux of the algorithm is to trackeach solution of H ( x, t ) = 0 for t ∈ [0 , using an efficient predictor-corrector method [26] to obtain all complex solutionsfor P ( x ) = 0 . Next we discuss the issue of maximum number of paths which should be tracked for a system of polynomials P ( x ) to guarantee all complex solutions of P ( x ) = 0 .Now, a system of m polynomials can have a maximum of (cid:81) mi =1 d i number of isolated complex solutions, where d i is thedegree of the i th polynomial, as specified by the classical B´ezout theorem. This poses an upper bound on the number ofsolutions to be tracked for NPHC method known as classical B´ezout bound (CBB). For solving our power flow problem of(12) there will be (2 N + n − algebraic equations, each of degree 2. Thus the number of paths to be tracked is (2 N + n − .This translates to the fact that in our model the number of paths to be tracked grows exponentially with the number of buses.Converting the power flow equations to their quadratic forms as in (12) actually pays off here by allowing a tighter upperbound on the number of complex isolated solutions or the CBB. However, this crude upper bound does not capture the specificcomplex algebraic structure of the polynomials of (12). Moreover, we are interested in solving (12) over a parameter spacegiven by specific penetration levels λ . Solving such parametric systems for every parameter-point from scratch using the NPHCmethod can be computationally very expensive. We, therefore, use a more sophisticated approach based on a variant of NPHC,parameter-coefficient homotopy [20], an earlier version of which was called Cheater’s homotopy [19]. This method uses the factthat for a parametric system of polynomial equations, the maximum number of isolated complex solutions over all parameter-points is same for a generic complex parameter-point. Hence, we can solve P ( x, λ ) = 0 at a generic complex parameter-point λ ∗ ∈ C m , using the NPHC method with the help of some crude upper bound on the number of complex solutions such as theCBB. Although such a complex parameter-point is physically not meaningful as λ represents the level of renewable penetration,solving the system at such a point reduces the computation for all other physically relevant parameter-points. Next, we choose P ( x, λ ∗ ) = 0 as the start system for all other parameter-points λ ∈ C m − { λ ∗ } . Each solution of this start system needs to betracked with the following homotopy: H ( x, λ, t ) = t P ( x, λ ∗ ) + (1 − t ) P ( x, λ ) = 0 , (14)from t = 1 to t = 0 . This procedure again guarantees all isolated complex solutions at each of the chosen parameter-points,independent of the upper bound chosen to solve the system at λ ∗ in the first step. In our simulations we will show thatdepending on operating conditions and structure of the power system, the number of paths to be tracked in the first step (e.g.the CBB) for power system models can dramatically reduce to a very small integer in the second step. The method thus becomes computationally very cheap once the first step is solved. Moreover, the process of solving the second step for aparameter vector λ is independent of a different parameter vector λ . Hence, if one intends to solve (12) for multiple λ atthe same time the process can be parallelized. For our simulations, we used a novel computational package called Paramotopy [21] which efficiently implements the above mentioned procedure with appropriate parallelization. Once all real solutions areobtained at each parameter point one can accurately estimate the power flow solution boundaries where the number of realsolutions changes, given the parameter space is densely represented.
A. On the Network Topology and Upper Bound on the Number of Equilibria
An actual system of polynomials may have fewer complex solutions as compared to its CBB. Thus, in the parameter-coefficient homotopy, it would be computationally wasteful to track the paths which would eventually diverge. Thus, for solvinglarge sets of power flow equations one should exploit the underlying structure or the sparsity of the network connectivity ofany given power system model, and try to compute a tighter upper bound. Recent reviews on the existing results on upperbounds are provided in [27] and [28]. An upper bound of (cid:0) N − N − (cid:1) was computed in [29]–[31] for a generic power flow problemwith N buses, although it did not still exploit the network topologies. In [32], the number of complex solutions for networkswith cliques with exactly one common node was shown to be equal to the product of number of complex solutions for theindividual cliques as independent networks. In [27], this result is extended to other related network topologies, though severalof the patterns for the particular topologies are still not well understood.Sparsity of network connectivity, for example, may result in lesser number of complex solutions of (12) than the usual CBBwhich is (2 N + n − . This is because the coefficients of the polynomial are related to each other according to the networktopology, and hence certain sparsity patterns may be such that a large subset of all the possible monomials of degree up tothe highest degree of the polynomials do not appear in those polynomials. For example, consider the power balance of bus 1in the 3-bus system in Figure 1a which can be written as, P e =Re (cid:40) (cid:88) k =2 S k (cid:41) + P eL (15a) Q e =Im (cid:40) (cid:88) k =2 S k (cid:41) + Q eL . (15b)where the power flow from bus 1 to any bus k is given as, S k = | V e | (cid:0) | V e | cos θ e − | V ek | (cos θ e sin θ ek − J sin ( θ e − θ ek )) (cid:1) arg ( Z ∗ k ) . The quadratic equations in (15) are not the densest polynomials of degree as the interaction of the unknown variablesoccur in only a specified structural form. For example, if the physical connection or the transmission line between bus 1 andbus k does not exist, then / arg ( Z ∗ ) = 0 . This eliminates the interaction terms of the variables associated with bus 1 andbus k . Thus the number of complex solutions at a generic complex parameter point can always be expected to be lesser thanthe CBB, particularly for a power flow problem. IV. E XAMPLES
A. 3 bus system
First we explore the solution space boundary for a 3-bus power system as shown in Figure 1a, which is a modification of anexample in [11] with added loads. The active power input at bus 1 and 2 are the variable parameters λ and λ respectively.Bus 3 is assumed to be the swing bus whose voltage equals ∠ . The unknown variables of the power flow problem whichconstitute the vector x in (7) are the active and reactive power input at bus 3, the reactive powers and angles of bus 1 and2. However as shown in Section III we represent the angles in rectangular form to limit the order of the algebraic equationsto two. Thus we have the sine and cosine of the angle of bus 1 and bus 2 as the unknown variables. The problem, therefore,reduces to the computation of the unknown vector, x = (cid:2) P Q Q Q sin δ cos δ sin δ cos δ (cid:3) (cid:48) , over a set of parameter values λ and λ . As mentioned earlier in Section III, we solve the system of equations P ( x, λ ∗ ) = 0 at a generic complex parameter point λ ∗ . The generic complex vector has two elements which are chosen from uniformdistributions such as { a + ib : − ≤ a, b ≤ } , and are normalized to ensure that they are inside unit circle. For the startproblem P ( x, λ ∗ ) = 0 , we followed the CBB and tracked = 64 paths to obtain all the isolated complex solutions. However,the start problem for this case yielded only complex solutions. Correspondingly, in the second stage of the algorithm,following the parameter-coefficient NPHC method we to track the 6 paths for each parameter point starting from the startsolution. Once all the complex solutions have been obtained the real solutions are identified by doing a numerical check. As (cid:42)(cid:72)(cid:81)(cid:3)(cid:20)(cid:42)(cid:72)(cid:81)(cid:3)(cid:21) (cid:42)(cid:72)(cid:81)(cid:3)(cid:22)(cid:42)(cid:72)(cid:81)(cid:3)(cid:23)(cid:37)(cid:88)(cid:86)(cid:3)(cid:20)(cid:37)(cid:88)(cid:86)(cid:3)(cid:21) (cid:47)(cid:82)(cid:68)(cid:71)(cid:3)(cid:20) (cid:47)(cid:82)(cid:68)(cid:71)(cid:3)(cid:21)(cid:58)(cid:76)(cid:81)(cid:71)(cid:3)(cid:83)(cid:79)(cid:68)(cid:81)(cid:87)(cid:3)(cid:20) (cid:58)(cid:76)(cid:81)(cid:71)(cid:3)(cid:83)(cid:79)(cid:68)(cid:81)(cid:87)(cid:3)(cid:21)(cid:37)(cid:88)(cid:86)(cid:3)(cid:28)(cid:37)(cid:88)(cid:86)(cid:3)(cid:24) (cid:37)(cid:88)(cid:86)(cid:3)(cid:25) (cid:37)(cid:88)(cid:86)(cid:3)(cid:20)(cid:19)(cid:37)(cid:88)(cid:86)(cid:3)(cid:26) (cid:37)(cid:88)(cid:86)(cid:3)(cid:27) (cid:37)(cid:88)(cid:86)(cid:3)(cid:22)(cid:37)(cid:88)(cid:86)(cid:3)(cid:23) (a) power system model (b) Number of real solutions of the power flow problem Fig. 2: 10-bus power system with two wind power plants with penetration levels λ and λ respectively (cid:21) (cid:21)(cid:17)(cid:24) (cid:22) (cid:22)(cid:17)(cid:24) (cid:23) (cid:23)(cid:17)(cid:24) (cid:24) (cid:24)(cid:17)(cid:24) (cid:25)(cid:19)(cid:19)(cid:17)(cid:24)(cid:20)(cid:20)(cid:17)(cid:24)(cid:21)(cid:21)(cid:17)(cid:24)(cid:22)(cid:22)(cid:17)(cid:24)(cid:23) (cid:79) (cid:20) (cid:79) (cid:21) (cid:3) (cid:3) (cid:76)(cid:81)(cid:76)(cid:87)(cid:76)(cid:68)(cid:79)(cid:29)(cid:3) (cid:79) (cid:20) (cid:32)(cid:21)(cid:17)(cid:23)(cid:24)(cid:21)(cid:23)(cid:15)(cid:3) (cid:79) (cid:21) (cid:32)(cid:21)(cid:76)(cid:81)(cid:76)(cid:87)(cid:76)(cid:68)(cid:79)(cid:29)(cid:3) (cid:79) (cid:20) (cid:32)(cid:24)(cid:17)(cid:24)(cid:23)(cid:26)(cid:25)(cid:15)(cid:3) (cid:79) (cid:21) (cid:32)(cid:21)(cid:76)(cid:81)(cid:76)(cid:87)(cid:76)(cid:68)(cid:79)(cid:29)(cid:3) (cid:79) (cid:20) (cid:32)(cid:21)(cid:17)(cid:19)(cid:26)(cid:27)(cid:28)(cid:15)(cid:3) (cid:79) (cid:21) (cid:32)(cid:22)(cid:76)(cid:81)(cid:76)(cid:87)(cid:76)(cid:68)(cid:79)(cid:29)(cid:3) (cid:79) (cid:20) (cid:32)(cid:23)(cid:17)(cid:19)(cid:22)(cid:27)(cid:28)(cid:15)(cid:3) (cid:79) (cid:21) (cid:32)(cid:22) Fig. 3: Power flow solution boundary tracking with various initial pointsseen in Figure 1b the number of real solutions vary from 0 to 6. The boundaries can be identified by the change in the numberof real solutions where the Jacobian becomes singular. ( λ , λ ) parameter space, as shown in Figure 1b is discretized into agrid of equispaced parameter points of dimension × . All solutions for each of the discrete points on the parameterspace are obtained henceforth by the application of numerical homotopy.Following Hiskens et al. [11], as a first step we keep λ fixed, and find an initial point on the solution boundary. Essentiallywe solve (7) by a Newton-Raphson method with different initial points for λ . When λ = 2 , we find initial points with λ = 2 . and . . When λ = 3 , the initial points are λ = 2 . and . . As seen in Figure 1b, all of thesepoints are located on the outer elliptical solution boundary. Correspondingly, in the second step when we solve (9) in aniterative form starting from these points they only track the outer boundary as shown in Figure 3. Thus following the methodin [11], it is difficult to identify all the power flow solution boundaries unless one has a sound knowledge of the solution spacefor a given set of parameters. On the contrary our alternate numerical homotopy based method can guarantee to identify all solution boundaries for a given power flow problem. In Figure 1b, it can be seen that via homotopy continuation method, wecould identify regions on the parameter space with different number of real solutions identified by the different colors. Theboundaries between these regions are the power flow solution boundaries. It can be noted that our method identifies boundariesbetween , , and 6 real solutions while the conventional method identifies the boundary between 0 and 2 solutions only. B. 10 bus system
In the next example we use a 10-bus, 4-synchronous machine power system as shown in Figure 2a. Two wind power plantsare connected at bus 9 and bus 10 whose penetration levels are represented as λ and λ respectively. The parameters of thesimulation are given in Appendix A. We vary the active power output of each of the wind plants between 0.1 and 0.7 p.u. ona 100 MVA base and find the solution of the power flow problem for each point on the plane defined by λ and λ . We first solve the problem at a generic complex parameter point using the homotopy based algorithm of section III. The problem has25 unknown equations which require paths to be tracked by continuation in the start problem to ensure all the isolatedroots. However, it turns out that the start system has only isolated roots. Correspondingly, we look for only paths inthe subsequent steps saving a lot of computational effort in finding roots for actual parameter values. Figure 2b shows thenumber of real solutions of the power flow problem for given values of λ and λ . The different colored regions of Figure 2bdemonstrate varying number of real solutions of the power flow. Thus the boundaries between the regions constitute the powerflow solution boundaries. It can be observed that the geometry of the solution boundary for the 10-bus case with varying windpenetration levels is strikingly different as compared to the 3-bus case. Also the number of isolated boundaries are more thanthat of the 3-bus case. All these observations point to the fact that identifying these boundaries by an initial guess and localapproximation is totally intractable for systems with large dimension. Thus, if the system has multiple wind power plants, thenour algorithm can provide the power system operator to choose an optimal set of power injections at the renewable buses.The operating points can be post processed for different robustness criteria and placed at a suitable distance away from theloadability boundary of the system. In this case as well, the parameter space ( λ , λ ) as shown in Figure 2b is discretizedinto a grid of equispaced parameter points of dimension × .Although finding a novel upper bound on the number of power flow solutions is not our goal, we still observed from theabove simulations that the number of complex solutions at a generic complex parameter-point is dramatically small comparedto the CBB of the system. It is evident from the simulation results that this number is a new and tighter upper bound on thenumber of complex isolated solutions compared to the previously known upper bounds. For example, for the 3 bus case, theCBB was whereas the number of complex solutions at a generic point was only , which is the same as the binomial boundmentioned above. However, for the 10 bus case, the CBB is while the number of complex solutions is = 512 , whichis much smaller compared to the binomial bound .V. C ONCLUSION
We developed a numerical power flow solution method that guarantees to identify all power flow solution boundaries ina power system in presence of variable generation parameters. Determination of all solution boundaries will be critical inthe foreseeable future due to the variability imposed by the rapid intrusion of renewable energy penetration. The essenceof our algorithm is based on a homotopy continuation concept, which also has the potential for accommodating topologicalinformation of the system. Our future research direction would include deriving an explicit relationship between the numberof solutions and the structure of the power system, that may lead to real-time usage of this tool during different contingencies.VI. A
CKNOWLEDGEMENT
The authors would like to thank Dr. Daniel K. Molzahn and Dr. Konstantin Turitsyn for their helpful suggestions anddiscussions on this topic. A
PPENDIX AM ODEL P ARAMETERS FOR THE
BUS SYSTEM
The conventional generators have the following power output on 100-MVA base: P e = 35 . p.u., P e = 17 . p.u., P e = 10 . p.u., P e = 40 . p.u.Line parameters in per unit on 100-MVA base: Z = (0 .
25 + j . e − , Z = (0 .
25 + j . e − , Z = (0 .
25 + j . e − , Z = (0 .
25 + j . e − , Z = (0 .
10 + j . e − , Z = (0 .
10 + j . e − , Z = (2 .
20 + j . e − , Z =(0 .
25 + j . e − , Z = (0 .
25 + j . e − Load parameters in per unit on 100-MVA base: S L = 30 + j . and S L = 70 + j .The rated power output of the wind plants at bus 9 and bus 10 are P er = 25 p.u. and P er = 15 p.u..R EFERENCES[1] National Renewable Energy Laboratory. 20% wind energy by 2030: Increasing wind energy’s contribution to U.S. electricity supply. Technical ReportDOE/GO-102008-2567, US Dept. of Energy, July 2008.[2] Prabha Kundur, Neal J Balu, and Mark G Lauby.
Power system stability and control , volume 7. McGraw-hill New York, 1994.[3] E Muljadi and YC Zhang. Wind power plant voltage stability evaluation. In
International Conference on Wind Energy Grid-Adaptive Technologies ,2014.[4] Yongning Chi, Yanhua Liu, Weisheng Wang, and Huizhu Dai. Voltage stability analysis of wind farm integration into transmission network. In , pages 1–7. IEEE, 2006.[5] Ce Zheng and Mladen Kezunovic. Impact of wind generation uncertainty on power system small disturbance voltage stability: A pcm-based approach.
Electric Power Systems Research , 84(1):10–19, 2012.[6] Thierry Van Cutsem and Costas Vournas.
Voltage stability of electric power systems , volume 441. Springer Science & Business Media, 1998.[7] Walmir Freitas, Luiz CP Da Silva, and Andre Morelato. Small-disturbance voltage stability of distribution systems with induction generators.
IEEETransactions on Power Systems , 20(3):1653–1654, 2005.[8] Eilyan Y Bitar, Ram Rajagopal, Pramod P Khargonekar, Kameshwar Poolla, and Pravin Varaiya. Bringing wind energy to market.
IEEE Transactionson Power Systems , 27(3):1225–1235, 2012. [9] A Jepson and A Spence. Folds in solutions of two parameter systems and their calculation. part i.
SIAM journal on numerical analysis , 22(2):347–368,1985.[10] Werner C Rheinboldt. Computation of critical boundaries on equilibrium manifolds.
SIAM Journal on Numerical Analysis , 19(3):653–669, 1982.[11] Ian Hiskens, Robert J Davy, et al. Exploring the power flow solution space boundary.
Power Systems, IEEE Transactions on , 16(3):389–395, 2001.[12] Michael E Karystianos, Nicholas G Maratos, and Costas D Vournas. Maximizing power-system loadability in the presence of multiple bindingcomplementarity constraints.
Circuits and Systems I: Regular Papers, IEEE Transactions on , 54(8):1775–1787, 2007.[13] Magnus Perninge and Lennart S¨oder. On the validity of local approximations of the power system loadability surface.
Power Systems, IEEE Transactionson , 26(4):2143–2153, 2011.[14] FMA Salam, L Ni, S Guo, and X Sun. Parallel processing for the load flow of power systems: the approach and applications. In
Decision and Control,1989., Proceedings of the 28th IEEE Conference on , pages 2173–2178. IEEE, 1989.[15] V. Ajjarapu and C. Christy. The continuation power flow: a tool for steady state voltage stability analysis.
IEEE Transactions on Power Systems ,7(1):416–423, 1992.[16] Weimin Ma and James S Thorp. An efficient algorithm to locate all the load flow solutions.
Power Systems, IEEE Transactions on , 8(3):1077–1083,1993.[17] Chih-Wen Liu, Chen-Sung Chang, Joe-Air Jiang, and Guey-Haw Yeh. Toward a cpflow-based algorithm to compute all the type-1 load-flow solutionsin electric power systems.
Circuits and Systems I: Regular Papers, IEEE Transactions on , 52(3):625–630, 2005.[18] D. Mehta, H. Nguyen, and K. Turitsyn. Numerical Polynomial Homotopy Continuation Method to Locate All The Power Flow Solutions.
Preprint: http://arxiv.org/abs/1408.2732, 2014.[19] TY Li, Tim Sauer, and JA Yorke. The cheater’s homotopy: an efficient procedure for solving systems of polynomial equations.
SIAM Journal onNumerical Analysis , 26(5):1241–1251, 1989.[20] Alexander P Morgan and Andrew J Sommese. Coefficient-parameter polynomial continuation.
Applied Mathematics and Computation , 29(2):123–160,1989.[21] DANIEL J Bates, DANIEL A Brake, and MATTHEW E Niemerg. Paramotopy: Parameter homotopies in parallel, 2012.[22] Souvik Chandra, Dhagash Mehta, and Aranya Chakrabortty. Equilibria analysis of power systems using a numerical homotopy method. In
Power &Energy Society General Meeting , pages 1–5. IEEE, 2015.[23] Tien Yien Li. Solving polynomial systems by the homotopy continuation method.
Handbook of numerical analysis , 11:209–304, 2003.[24] Andrew John Sommese and Charles Weldon Wampler.
The Numerical solution of systems of polynomials arising in engineering and science , volume 99.World Scientific, 2005.[25] Alexander Morgan and Andrew Sommese. Computing all solutions to polynomial systems using homotopy continuation.
Applied Mathematics andComputation , 24(2):115–138, 1987.[26] Daniel J Bates, Jonathan D Hauenstein, Andrew J Sommese, and Charles W Wampler.
Numerically solving polynomial systems with Bertini , volume 25.SIAM, 2013.[27] Daniel K Molzahn, Dhagash Mehta, and Matthew Niemerg. Toward topologically based upper bounds on the number of power flow solutions. arXivpreprint arXiv:1509.09227 , 2015.[28] Dhagash Mehta, Daniel K Molzahn, and Konstantin Turitsyn. Recent advances in computational methods for the power flow equations. arXiv preprintarXiv:1510.00073 , 2015.[29] J. Baillieul and C.I. Byrne. Geometric Critical Point Analysis of Lossless Power System Models.
IEEE Trans. Circuits Syst. , 29(11), 1982.[30] Tien-Yien Li, Tim Sauer, and James A Yorke. Numerical solution of a class of deficient polynomial systems.
SIAM journal on numerical analysis ,24(2):435–451, 1987.[31] Jakub Marecek, Timothy McCoy, and Martin Mevissen. Power flow as an algebraic system. arXiv preprint arXiv:1412.8054 , 2014.[32] S.X. Guo and F.M.A. Salam. Determining the solutions of the load flow of power systems: Theoretical results and computer implementation. In