A Dynamized Power Flow Method based on Differential Transformation
>> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) < 1 Abstract — This paper proposes a novel method for solving and tracing power flow solutions with changes of a loading parameter. Different from the conventional continuation power flow method, which repeatedly solves static AC power flow equations, the proposed method extends the power flow model into a fictitious dynamic system by adding a differential equation on the loading parameter. As a result, the original solution curve tracing problem is converted to solving the time domain trajectories of the reformulated dynamic system. A non-iterative algorithm based on differential transformation is proposed to analytically solve the aforementioned dynamized model in form of power series of time. This paper proves that the nonlinear power flow equations in the time domain are converted to formally linear equations in the domain of the power series order after the differential transformation, thus avoiding numerical iterations. Case studies on several test systems including a 2383-bus system show the merits of the proposed method.
Index Terms — Continuation power flow; dynamized power flow; differential transformation; power flow; power-voltage curve; voltage stability; voltage collapse I. I NTRODUCTION
RACING solution curves of power flow equations with the changes of a trending parameter such as the loading level is usually computation-intensive task in power system operations and planning to prevent steady-state voltage instability and other insecurities [1]-[4]. An example is computation of the power-voltage (P-V) curves for critical buses with the increase of load. Traditionally, the continuation power flow (CPF) method [5]-[9] is widely used to solve P-V curves, which adopts a prediction-correction scheme to identify a series of power flow solutions along the solution curve where each prediction step gives an initial guess and the following corrector step performs numerical iterations to find the converged solution [10]-[12]. However, the CPF method may suffer from huge computation burdens since it requires solving nonlinear AC power flow equations for multiple times using numerical iteration methods [13]-[14]. Moreover, the computation burden of the CPF method can further grow and become unacceptable with modern power grids being integrated with high renewable energy resources and demand response, where such a solution process with numerical
This work was supported in part by the ERC Program of the NSF and DOE under NSF Grant EEC-1041877 and in part by NSF Grant ECCS-1610025. Y. Liu, K. Sun and J. Dong are with the Department of EECS, University of Tennessee, Knoxville, TN 37996 USA (e-mail: [email protected], [email protected], [email protected]). iterations needs to be repeated many times for multiple contingencies. In the literature, some techniques are proposed to reduce the computation burden of the CPF method [15]-[17], falling into two categories. The first category of techniques aim to design a more effective predictor than a standard tangent predictor [10]-[12] as adopted by many commercial CPF programs [13]-[14]. For example, paper [15] proposes three types of nonlinear predictors to predict a new solution from more than one previous solutions using interpolation and polynomial approximation including the Lagrange’s polynomial interpolation formula, Newton’s forward and backward divided difference formula, and cubic spline interpolation method. Besides, a secant predictor is used in [16] and a holomorphic embedding-based predictor is proposed in [17]. Methods in the first category are able to generate a better initial guess for the Newton Raphson (NR) method so as to reduce the total number of iterations; however, they still require many numerical iterations for solving nonlinear power flow equations. The second category of techniques focus on more efficient correctors than the standard NR method-based corrector. For example, authors in [15] propose a hybrid corrector allowing the switches between a NR method (taking its merit of robustness) and a fast decoupled NR method (taking its merit of fast speed) until a pre-defined maximum total number of iterations. Methods in this category are mainly used to improve the convergence performance of numerical iterations, but still require solving nonlinear AC power flow equations repeatedly. Overall, a major bottleneck for the above two categories of methods lies in their inherent solution mechanism that the power-flow equations are essentially solved as algebraic equations in an iterative manner. To more efficiently trace solution curves of power flow equations, this paper proposes a novel dynamized power flow (DPF) method that extends the power flow model into a fictitious dynamic system, called a “dynamized” power flow model, by adding a differential equation about a fictitious time, and then solve the complete time-domain trajectory of the dynamic system instead of repeatedly solving power flow equations for a series of conditions. A differential transformation (DT) method, which is proved effective for solving power system transient stability simulation in our recent works [18]-[20], is applied to solve the dynamized model, named as dynamized power flow (DPF) method. This paper proves that the nonlinear AC power flow equations are converted to formally linear equations after DT, and further
A Dynamized Power Flow Method based on Differential Transformation
Yang Liu,
Student Member , IEEE , Kai Sun,
Senior Member , IEEE , Jiaojiao Dong,
Member , IEEE T REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) < 2designs an efficient algorithm to solve the time domain trajectory without numerical iterations. Case studies on several test systems including a 2383-bus system demonstrate the accuracy, computational complexity and time performance of the proposed approach compared with a CPF solver. The rest of the paper is organized as follows. Section II gives the problem description. Section III presents the proposed method. Section IV is case study and Section V draws conclusions. II. P ROBLEM S TATEMENT
The conventional power flow equations are given in (1a) where S is a vector of the complex power injections, V is a vector of bus voltage phasor, and Y bus is the bus admittance matrix. By adding the product of a loading parameter λ and a constant vector b to the left-hand side, a general continuum of power flow equations is formulated in (1b). *bus *bus S=V(Y V) (a)S+ b=V(Y V) (b) (1) Note that the vector b is defined to reflect an arbitrarily direction of load changes, for example, uniform increases of all generation and load, or increases of generation and load at certain buses or zones. Meanwhile, practical operating constraints such as the reactive power limit of generators can be considered during the load change. Equation (1b) is further written as the general form in (2) where g is a nonlinear vector field; y is the bus voltage vector under rectangular coordinates defined as y =[ e T , f T ] T , where e = [ e ,…, e N ] T and f = [ f ,…, f N ] T are respectively the real and imaginary parts of the bus voltage phasor; N is the total number of buses; λ is the loading parameter. g y (2) The goal is to determine how power flow solution y changes with loading parameter λ, shown in (3). After (3) is obtained, the other system variables (such as voltage magnitude and power injections) are easily calculated to draw P-V curves. ( ) y y (3) Generally, analytical expression of (3) is unavailable due to the nonlinearity of g in (2). Therefore, a prediction-correction scheme and numerical iterations are needed in conventional CPF method. III. P ROPOSED D YNAMIZED P OWER F LOW M ETHOD A. Introduction of the Differential Transformation
A smooth nonlinear function of time x ( t ) can be approximated by a K -th order polynomial function of time as shown in (4), where X ( k ) is the k th order power series coefficient and can be calculated by (5). ( ) ( ) K kk x t X k t (4) ( )1( ) ! k k t d x tX k k dt (5) Generally, these power series coefficients are calculated in a recursive manner from k =0 to k = K , and many mathematical methods can be used such as the Adomain decomposition method [21]-[22] and the power series-based method in [23]-[24]. However, the applications of the above methods are limited by their huge computational burdens in deriving power series coefficient X ( k ) using the complicated high order derivative operations. As an emerging mathematical tool, DT [25]-[28] considers power series coefficient X ( k ) as a transformation of x ( t ) at the k th order as shown by (6). When multiple functions like x ( t ) are to be calculated and analyzed, their high order power series coefficients can directly be operated and calculated based on transformation rules introduced by DT. Thus, there is no need to calculate complicated high order derivatives of each function. ( ) ( ) x t X k (6) Our recent paper [18]-[19] introduces DT to the power system field to effectively solve power system nonlinear differential-algebraic equations (DAEs) for transient stability simulation. New transformation rules for nonlinear functions in power system models are proved in [18]-[19]. Five rules are given in (7) and will be utilized in this paper. Here, X ( k ) and Y ( k ) are DTs of functions x ( t ) and y ( t ), c is a constant, and is the Kronecker delta function: i) ( ) ( ) ( ) ( ); ii) ( ) ( )iii) ( ) ( ) ( ) ( ) ( ) ( )1, 0( )iv) ( 1) ( 1); v) ( )= 0, 0 km x t y t X k Y k cx t cX kx t y t X k Y k X m Y k mkdx t k X k c c k kdt (7) B. Conceptual Description of the Proposed Method
The proposed method has following four steps, where each step is first briefed below and then described in detail in Section III-C to Section III-F respectively. First, the algebraic equation (2) is extended to a set of DAEs by introducing a fictitious time t and adding two new equations, i.e., (8a) and (8b). Differential equation (8a) is a dynamic system to trace the changes of system variables such as power or voltages, where x ( t ) is a state variable and f (·) is a vector field. Algebraic equation (8b) is an ancillary equation that builds the relationship between the newly introduced variable x ( t ) and the original variables y ( t ) and λ( t ). Note that the ancillary function h may not be needed if x ( t ) is selected from one of the variables in y ( t ) and λ( t ). The details of designing (8a) and (8b) are described in Section III-C. ( ) ( ( ), ( ), ( )) ( )0 ( ( ), ( ), ( )) ( )0 ( ( ), ( )), . ., . (2) ( ) x t f x t t t ah x t t t bt t i e Equ c yyg y (8) Second, the DTs of (8a)-(8c) are derived in (9a)-(9c) respectively, using the transformation rules in (7). Specifically, the nonlinear power flow equation (8c) is converted to a new set of equations (9c) that couples the power series coefficients of y ( t ) and λ( t ) in all orders, i.e., Y (0)… Y ( k ), Λ(0)… Λ( k ). REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) < 3 ( 1) ( 1) ( (0 : ), (0 : ), (0 : )) ( )0 ( (0 : ), (0 : ), (0 : )) ( )0 ( (0 : ), (0 : )) ( ) k X k F X k k k aH X k k k bk k c YYG Y (9) Third, we prove that both (9c) and (9b) satisfy formally linear equations about the k th order power series coefficients Y ( k ) and Λ( k ), as shown in (10a) and (10b), respectively, where A matrices are functions of Y (0) and Λ(0) and B matrices are functions of Y (0: k -1) and Λ(0: k -1). As a result, Y ( k ), i.e., the k th order power series coefficient of bus voltage vector, is analytical solved from Y (0: k -1) and Λ(0: k -1), either by (11) or by (12), depending on if the ancillary function h is needed when designing the differential equation in (8). gy g ghy h h k k ak k b A Y A BA Y A B (10) ( ) ( ( ) ) gy g g k k Y A A B (11) ( )( ) gy g ghy h h kk Y A A BA A B (12) Finally, we design a non-iterative algorithm based on (9a) and (11) or (12) to solve power series coefficients
X(k ), Y ( k ) and Λ( k ) from k =0 to any order K in a recursively manner, and approximate variables x ( t ), y ( t ) and λ( t ) as power series of time, shown in (13). After y ( t ) and λ( t ) are solved, the solution curves of power flow equations are directly obtained, as illustrated in Section IV-A. ( ) (0) (1) (2) ... ( )( ) (0) (1) (2) ... ( )( ) (0) (1) (2) ... ( ) KKK x t X X t X t X K tt t t K tt t t K t y Y Y Y Y (13) Among the above four steps, only the last step needs to be performed online, while the first three steps can be conducted in the offline stage because they are mainly used to derive expressions of matrices A and B in (10) and function F in (9a), which is a one-time effort. Remarks : there are two important observations: 1) from (10a) that the nonlinear power flow equation (2) about y ( t ) are converted to a formally linear equation about power series coefficients Y ( k ) after DT; 2) coefficients on bus voltages are explicitly solved by (11) or (12) and then used to calculate bus voltages by (13) in a straightforward manner, which is different from using a conventional power flow solver to calculate bus voltages by numerical iterations. The proposed DT based method for solving solution curves of power flow equations differentiates itself from the traditional continuation power flow method that relies on iterative numerical methods such as the family of Newton Raphson methods. C. Step 1: Dynamization of the Power Flow Equation
Two formulations of (8) are proposed to dynamize the power flow equation (2), shown in (14) and (15) respectively, where C and C are constants and v l ( t ) is the voltage magnitude of a load bus l . In (14), there is no ancillary equation (8b) because the selected state variable λ ( t ) has existed in (2). In (15), the ancillary equation gives the relationship between bus voltage magnitude and the rectangular coordinate components. Formulation 1: ( )0 ( ( ), ( )), . ., .(2) t C t t i e Equ g y (14) Formulation 2: ( )0 ( ) ( ) ( )0 ( ( ), ( )), . ., .(2) l l l l v t Cv t e t f tt t i e Equ g y (15) For Formulation 1, its purpose is to characterize how the power changes with time, i.e., the power increases with time when C >0 and decreases with time when C <0. It can be used to trace curve segment in various shapes, either monotonically or non-monotonically, such as the curves (a), (b) and (c) in Fig. 1. For Formulation 2, its purpose is to characterize how the voltage magnitude changes with time, i.e., the voltage magnitude increases with time when C >0 and decreases with time when C <0. It can also be used to trace either monotonical or non-monotonical curve segments such as (a), (b) and (d) in Fig. 1. Fig. 1. Illustration of the two dynamized formulations for tracing curve segments of power flow equations
The above two formulations can be flexibly used to trace the full solution curve of a power flow equation. For example, the high voltage solutions in a P-V curve can be traced by Formulation 1 with C >0, the low voltage solutions can be traced by Formulation 1 with C <0, and the solution curves near the nose point can be traced by Formulation 2 with C <0. D. Step 2: Deriving DTs of the Dynamized Power Flow Model 1)
DTs of Nonlinear Power Flow Equation
The nonlinear power flow equation (2) is written into (16)-(19) under rectangular coordinates, where Ω PQ , Ω PV , Ω REF are the set of PQ buses, PV buses and reference bus respectively, p and q are active and reactive power, e and f are the real and imaginary parts of bus voltages, g and b are real and imaginary parts of the admittance, v is the voltage magnitude, superscript sp means the value is specified, subscript i and j are the index of buses. ( , ) if spi p N Ni ij i j i j ij i j i jj jPQ PV p g p g e e f f b f e e fi y (16) REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) < 4 ( , ) if spi q N Ni ij i j i j ij i j i jj jPQ q g q b e e f f g f e e fi y (17) ( ) ( ) , if spi v i i PV v g e f i y (18) ( ) , if ( ) , if spi e i REFspi f i REF e g e if g f i yy (19) The DTs of (16)-(19) are derived in (20)-(23), respectively. ( ) ( , ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) if Y spi p Ni ij i j i jjN ij i j i jj PQ PV p k G p k g E k E k F k F kb F k E k E k F ki (20) ( ) ( , ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) if Y spi q Ni ij i j i jjN ij i j i jj PQ q k G q k b E k E k F k F kg F k E k E k F ki (21) ( ) ( ) ( ) ( ) ( ) ( ) ( ), if Y spi vi i i i PV v k GE k E k F k F k i (22) ( )= ( ) ( ), if ( )= ( ) ( ), if YY spi e i REFspi f i REF e k G E k if k G F k i (23) For details, the derivation of (20) is elaborated below as an example. The remaining equations (21)-(23) are obtained in a similar procedure. The left-hand-side (LHS) and the first term in the right-hand-side (RHS) of (20) are obtained by applying the rule (7-i), (7-ii) and (7-v) to the corresponding terms of (16). Note that p isp and Δ p i are constants and λ = λ ( t ) is a variable, therefore their transformations are: p i sp p i sp ( k ) and Δ p i λ Δ p i Λ ( k ). The remaining terms in the RHS of (20) are obtained from the corresponding terms of (16) by following steps: First, apply the rule (7-iii): ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) i j i j i j i ji j i j i j i j e e E k E k f f F k F kf e F k E k e f E k F k Then, apply the rule (7-i) and (7-ii): ( ) ( ) ( ) ( )( ) ( ) ( ) ( ) ij i j i j ij i j i jij i j i j ij i j i j g e e f f g E k E k F k F kb f e e f b F k E k E k F k
Finally, using the rule (7-i), the RHS of (20) is obtained. DTs of the Designed Differential Equations
For the differential equation in
Formulation
1, i.e., (14), its DT is derived as follows. After applying the rule (7-iv) for LHS and (7-v) for RHS, ( k +1) Λ ( k +1)= C ( k ) holds. Then, it can be further written in (24) after replacing k by k -1 and using the definition of ( k ). ( ) ( 1) k C k (24) For Formulation
2 in (15), the DT of the differential equation is in (25) where the derivation is similar as (24) and is omitted here; the DT of the ancillary equation is in (26) after applying the rule (7-iii) to both sides. ( ) ( 1) l V k C k (25) ( ) ( ) ( ) ( ) ( ) ( ) l l l l l l V k V k E k E k F k F k (26) E. Step 3: Proving the Nonlinear Power Flow Equations Satisfy Formally Linear Equations after DT
Proposition 1 : The transformed power flow equations (20)-(23) respectively satisfy formally linear equations (27)-(30). P, i i i PQ PV k p k i a Y (27) Q, i i i PQ k q k i a Y (28) V, i i PV k k i a Y (29) E,F, a Ya Y spi i REFspi i REF k k e k ik k f k i (30) where Y ( k ) N and Λ ( k ) are variables representing the DT of y and λ respectively; a P, i , a Q, i , a V, i a E, i , a F, i N and ε i , μ i , ζ i are parameters given in (44)-(49) respectively. The detailed proof of Proposition 1 is in Appendix. From the Proposition , DTs (9c) of the nonlinear power flow equation satisfy formally linear equation (10a) with matrices A gy , A g λ , and B g given by (31). For notation simplicity, here we let buses 1 to M be PQ buses, buses M +1 to N -1 be PV buses and bus N be the reference bus. ,PQ ,PQ PQ,PV ,PV PV,REF ,REF REF , , yygy g gy A A BA A BA A BA A B (31) y M M MM M M pqpq aaA A Baa
0 0
M MM MMy N N NN N pp aaA A Baa REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) < 5
E,,REF ,REF REFF, ( ) ( ) a 0A A Ba 0 spN Ny spN N e kf k Besides, the DT (26) of the ancillary equation in Formulation II also satisfies a formally linear equation in (10b) with proof in the Appendix. F. Step 4: Non-iterative Algorithm for Solving Variables as Power Series of Time
Following the basic idea in Section III-B, an algorithm is designed to solve power series coefficients
X(k ), Λ ( k ), Y ( k ) up to any desired order. Note that these coefficients are explicitly calculated with no numerical iteration. Algorithm: Solve Coefficients
Input : initial values (0), (0), (0) x y Output : any order coefficients ( ), (, ,) ) 0(
X Kk kk k Y Steps : Initialization:
X k x tkk tt yY Calculate matrices A gy, A g λ , A hy, A h λ using (31) Calculate gy A if Formulation I or hy hgy g A AA A if Formulation II for k K Calculate matrix B g, B h using (31) if Formulation I Solve k using ( ) ( 1) k C k Solve ( ) k Y using ( ) ( ( ) ) gy g g k k Y A A B if Formulation II Solve l V k using ( ) ( 1) l V k C k Solve ( ), kk Y using ( )( ) gy g ghy h h kk Y A A BA A B end
After the power series coefficients are calculated, y ( t ) and λ( t ) are calculated by evaluating the power series of time in (13) and the solution curves are directly obtained. In practical, the multi-time window strategy [18]-[19] can be used to extend the convergence region of power series of time and ensure the accuracy. The time step length as well as the order K of the power series of time are usually selected from trial simulations [19], and the impact of K and time step length are also studied in [18]-[19]. The proposed method can also be applied to more complicated power system models such as 1) considering reactive power limit of generators, 2) ZIP load model. First, to consider the reactive power limit of generators, the proposed method can be slightly modified as follows: if a generator meets the reactive power limit, then it is changed from a PV bus to a PQ bus; correspondingly, the matrices A and B need to be re-calculated using (31). Later in Section IV-A, we demonstrate the proposed method for tracing PV curves considering reactive power limits. Second, for a nonlinear power flow equation with ZIP loads, we proved in [29] that its DTs still satisfy formally linear equations. Therefore, the proposed method can be directly applied with slight modification on matrices A and B . Regarding the computational complexity, the proposed method has two unique features: First, it shifts most of the computation burden to the offline stage, i.e., deriving the equation for calculating matrices A and B , which is a one-time effort (the matrices A and B derived in this paper can be directly used by others without deriving them again); and the online stage only involves explicit matrix operation and evaluation of analytical solutions, which do not require any numerical iterations. Second, the proposed method can reduce the frequency of solving linear equations compared with the CPF method, thus having better computational efficiency. This is because the CPF method needs to solve a linear equation in every iteration and every prediction-correction step, while the proposed method only needs to solve linear equations once in each time step, and the total number of time steps are greatly reduced benefiting from the high order approximation. IV. C ASE S TUDIES
The proposed DPF method is first tested on the IEEE 9-bus system [13] to demonstrate the basic idea, the impact of load change directions, and the impact of reactive power limit of generators. Then, the accuracy, computational complexity, and computation time are compared with the CPF method in MATPOWER using several large systems including the IEEE 39-bus system, IEEE 57-bus system and a Polish 2383-bus test system [13]. At last, the proposed approach is applied to N -1 contingency analysis. Simulations are conducted in MATLAB R2017a on a personal computer with i5-8250U CPU. Without specification, generations and loads of all buses are uniformly increased. For the CPF method, various simulation control parameters are adjusted for the best time performance, including using the pseudo arc-length for parameterization, enabling adaptive step size, increasing the maximum allowed step size and disabling the incremental curve plotting in each iteration, etc. For the DPF method, parameters C and C are set as 1, K is set as 6 from trail simulations, and the time step length is fixed at 0.05s for 2383-bus system and 0.1s for other systems. A. Demonstration on the 9-bus Power System
To demonstrate the idea of the proposed method, Fig. 2 and Fig. 3 respectively give the time domain trajectories of the solved dynamized power flow model and the obtained PV curve. In the first 1.63s, the loading parameter λ increases with time in a constant rate while the voltage magnitude of bus 9 drops from 0.9956 p.u. to 0.6268 p.u., indicating high voltage solutions. During the time period between t =1.63s and t =1.68s, the voltage is decreased from 0.6268 p.u. to 0.5439 p.u., while the loading parameter is first increased from 1.63 to reach the maximum value 1.64 and then decreased to 1.63, indicating the dynamic process of passing the nose point. Finally, both the loading parameter and the bus voltage are decreased after t =1.68s, indicating the low voltage solutions. The obtained loading limit 1.64 is the same as the limit from the CPF method. REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) < 6 Fig. 2. Time domain trajectory of the dynamized power flow model B u s vo lt a g e m a gn it ud e ( p . u . ) Fig. 3. Solution curve of load bus 9 on 9-bus system
Two scenarios are designed to demonstrate the capability of the proposed method on handling load changes with 1) non-uniform directions and 2) reactive power limits. Fig. 4a shows the PV curve of load bus 9 when increasing generation at bus 3 and load at bus 7 by 50 MW in active power and 10 MVar in reactive power. Fig. 4b further shows the PV curve of the same bus when reactive power limit of generators is considered. It shows the calculated maximum loading limits are reduced from 8.17 to 7.79 due to the reactive power limit. These results demonstrate the performance of the proposed method on practical power system models and applications. (a) (b) Fig. 4. Solution curves under non-uniform load change direction, a) not consider reactive power limit, b) consider reactive power limit B. Performance on Large Systems: Accuracy, Computation Complexity, and Time Performance
Respectively for the 39-bus system, the 57-bus system and the 2383-bus system, the proposed DPF method is compared with the CPF method. In all following studies, the CPF method is tested using the commercial MATPOWER package while the proposed DPF method is tested using our research code. Fig. 5 to Fig. 7 show the PV curves of three load buses, obtained by both the proposed method and the CPF method. Respectively for the three test systems, the calculated loading limits are 1.12, 0.88 and 0.89 for DPF method, and 1.13, 0.89 and 0.89 for the CPF method. These results demonstrate the accuracy of the proposed method.
Fig. 5. Solution curves on 39-bus system Fig. 6. Solution curves on 57-bus system Fig. 7. Solution curves on 2383-bus system
For both the CPF method and the DPF method, a major computation burden is in solving linear equations. Table I gives how many times linear equations are solved for both methods. It shows that the proposed approach is 10 times fewer than the CPF method for all the three test systems. This is because the CPF method needs to solve a linear equation in each iteration and for every prediction-correction step while the proposed method only solves a linear equation once in each time step.
TABLE
I N
UMBERS OF T IMES OF S OLVING L INEAR E QUATIONS
Test Systems CPF DPF 39-bus system 174 11 57-bus system 108 10 2383-bus system 424 18
Table II further gives the computation times of both methods. It shows the proposed DPF method is around 9 times, 12 times, and 2 times faster than the CPF method, respectively, for the three test systems. The speed up on the 2383-bus system is less than speedups on the other two smaller systems because our current academic research code that implements the DPF method in MATLAB has not been optimized to as efficiently REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) < 7handle large-scale matrix operations as the commercial CPF solver in the MATPOWER. However, these test results do demonstrate the potential of the proposed DPF method for online power flow solution tracing and voltage stability assessment.
TABLE
II C
OMPARISON OF T IME P ERFORMANCE (U NIT : SECOND ) Test Systems CPF DPF 39-bus system 0.26 0.03 57-bus system 0.50 0.04 2383-bus system 24.45 10.13 C. Application to N-1 Contingency Analysis
The proposed approach is further applied to screen N -1 contingencies. For the 39-bus system, 46 contingencies are created each with the loss of each single branch. Fig. 8 shows the maximum loading condition identified by both methods. Using the CPF results as benchmarks, the DPF method is accurate and reliable for all the contingencies. Regarding the computation time, the CPF method and the DPF method respectively takes 12.0 s and 1.4 s, showing that the DPF method can identify insecure contingencies much faster than the CPF method, and thus can scan more contingencies than the CPF method within limited time in the real-time environment. Fig. 8. Identified maximum loading condition by CPF and DPF method V. C ONCLUSION
In this paper, a novel dynamized power flow method has been proposed to efficiently trace solution curves of power flow equations. The original curve tracing problem for steady-state power flow solutions is converted to a time domain simulation problem about a dynamized model after adding a differential equation on changes of the operating condition. An DT-based approach is proposed for efficiently solving the dynamized model without numerical iterations. Simulation results have shown high accuracy, reduced computational complexity and improved time performance of the proposed DPF method compared with a CPF solver in MATPOWER. Besides, the proposed method can deal with practical engineering constraints such as the non-uniform load change directions and reactive power limits of generators. A
PPENDIX
To make the proofs more compact, the following Lemma is first proved. In the Lemma, the transformation of multiplication operation from time domain to the convolution operation in the domain of power series orders is well-known in many DT literatures, however, the resulted linear relationship in (33)-(34), despite their simplicity and being straightforward, are rarely noticed and exploited as far as the authors know.
Lemma:
The DT of z ( t )= x ( t ) y ( t ), satisfies a formally linear equation in (32). Especially, when x ( t )= y ( t ), (33) holds. ( ) ( ) ( ) ( ) ( ) Z k X k Y k aX k bY k c (32) ( ) ( ) ( ) 2 ( )
Z k X k X k aX k c (33)
Proof of Lemma: ( ) ( ) ( ) ( ) ( )(0) ( ) ( ) (0) ( ) ( ) km km
Z k X k Y k X m Y k mX Y k X k Y X m Y k m Therefore, (32) holds with a, b and c given below. (0), (0), ( ) ( ) km a Y b X c X m Y k m Proof of Proposition 1:
Use (27) as an example. The RHS of (20) is rewritten as
Term 11, Term 21 Term
RHS ( ) ( ) ( ) ( ) ( )( ) ( ) ( ) ( )( ) ( ) ( ) ( ) i ii i i i iN ij i j i jj j iN ij i j i jj p k g E k E k F k F kg E k E k F k F kb F k E k E k F k According to the Lemma, the three terms are rewritten as:
Term 1 2 (0) ( ) 2 (0) ( )( ) ( ) ( ) ( ) ii i i ii i ik kii i i ii i im m g E E k g F F kg E m E k m g F m F k m
Term 2 (0) ( ) (0) ( )(0) ( ) (0) ( )( ) ( ) ( ) ( )
N ij j i i jjj iN ij j i i jjj iN k kij i j i jj m mj i g E E k E E kg F F k F F kg E m E k m F m F k m
11 1 11 1 1
Term 3 (0) ( ) (0) ( )(0) ( ) (0) ( )( ) ( ) ( ) ( )
N ij j i i jjN ij j i i jjN k kij i j i jj m m b E F k F E kb F E k E F kb F m E k m E m F k m
Finally, (27) is obtained by summating the above three terms, with vector a P, i and parameter ε i in (34) and (39). Similarly, (28)-(30) can be proved with vectors a P, i , a Q, i , a V, i , a E, i , a F, I and parameters ε i , μ i , ζ i in (34)-(41). REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) < 8 P, 1 111 , where(0) (0), (0) (0), if (0) (0 (0) (0)(0) (0) (0) (0) i i i ij ijij ij i ij i ij ij i ij iNii ij j ij j ii i ii ijNii ij j ij j ii i ii ij g E b F g F b E j ig E b F g E b Fb E g F b E g F a (34) Q, 1 111 , where(0) (0), (0) (0), if (0) (0 (0) (0)(0) (0) (0) (0) i i i ij ijij ij i ij i ij ij i ij iNii ij j ij j ii i ii ijNii ij j ij j ii i ii ij b E g F b F g E j ib E g F b E g Fg E b F g E b F a (35) V, i i i E F a (36) E, a i (37) F, a i (38) ( ), where: ( ) ( ) ( ) ( ): ( ) ( ) ( ) ( ) N Ni ij ij ij ij ij jk kij i j i jm mk kij i j i jm m g c b d p kc E m E k m F m F k md F m E k m E m F k m (39) ( )
N Ni ij ij ij ij ij j b c g d q k (40) ( ) i ii i c v k (41) Proof of (10b) from (26):
From the
Lemma , there are ( ) ( ) 2 (0) ( ) ( ) ( )( ) ( ) 2 (0) ( ) ( ) ( ) kl l l l l lmkl l l l l lm
E k E k E E k E m E k mF k F k F F k F m F k m Then, (26) is rewritten as: ( ) ( ) 2 (0) ( ) 2 (0) ( )+ ( ) ( ) ( ) ( ) l l l l l lk kl l l lm m
V k V k E E k F F kE m E k m F m F k m
Finally, (10b) is obtained with a l and ξ l given in (42). l l lk kl l l l l l lm m E FE m E k m F m F k m V k V k a (42) R EFERENCES [1]
J. W. Simpson-Porco, F. Dörfler, F. Bullo, "Voltage collapse in complex power grids,"
Nature communications , vol. 7, p. 10790, 2016. [2]
I. Dobson, H.D. Chiang, "Towards a theory of voltage collapse in electric power systems,"
Systems & Control Letters , vol. 13, no. 3, pp. 253-262, 1989. [3]
H.D. Chiang, et al, "On voltage collapse in electric power systems,"
IEEE Tran. Power Syst. , vol. 5, no. 2, pp. 601-611, May 1990. [4]
I. Dobson, "Observations on the geometry of saddle node bifurcation and voltage collapse in electrical power systems,"
IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications , vol. 39, no. 3, pp. 240-243, March 1992. [5]
V. Ajjarapu, C. Christy, "The continuation power flow: a tool for steady state voltage stability analysis,"
IEEE Tran. Power Syst. , vol. 7, no. 1, pp. 416-423, Feb. 1992. [6]
H. Song, S. Kim, B. Lee, S. H. Kwon, V. Ajjarapu, "Determination of interface flow margin using the modified continuation power flow in voltage stability analysis,"
IEE Proceedings-Generation, Transmission and Distribution , vol. 148, no. 2, pp. 128-132, 2001. [7]
A. Dukpa, B. Venkatesh, M. El-Hawary, "Application of continuation power flow method in radial distribution systems,"
Electric Power Systems Research , vol. 79, no. 11, pp. 1503-1510, 2009. [8]
X.P. Zhang, P. Ju, E. Handschin, "Continuation three-phase power flow: A tool for voltage stability analysis of unbalanced three-phase power systems,"
IEEE Tran. Power Syst. , vol. 20, no. 3, pp. 1320-1329, 2005. [9]
H.D. Chiang, et al, "CPFLOW: A practical tool for tracing power system steady-state stationary behavior due to load and generation variations,"
IEEE Tran. Power Syst. , vol. 10, no. 2, pp. 623-634, 1995. [10]
P. Kundur, N. J. Balu, M. G. Lauby,
Power system stability and control . McGraw-hill New York, 1994. [11]
F. Milano,
Power system modelling and scripting . Springer Science & Business Media, 2010. [12]
T. Van Cutsem, C. Vournas,
Voltage stability of electric power systems . Springer Science & Business Media, 2007. [13]
R. D. Zimmerman, et al, “Matpower: Steady-state operations, planning and analysis tools for power systems research and education,”
IEEE Tran. Power Syst. , vol. 26, no. 1, pp. 12–19, 2011. [14]
F. Milano, "An open source power system analysis toolbox,"
IEEE Tran. Power Syst. , vol. 20, no. 3, pp. 1199-1206, 2005. [15]
S. H. Li, H.D. Chiang, "Nonlinear predictors and hybrid corrector for fast continuation power flow,"
IET generation, transmission & distribution , vol. 2, no. 3, pp. 341-354, 2008. [16]
D. A. Alves, L. C. P. da Silva, C. A. Castro, V. F. da Costa, "Continuation fast decoupled power flow with secant predictor,"
IEEE Tran. Power Syst. , vol. 18, no. 3, pp. 1078-1085, 2003. [17]
B. Wang, C. Liu, K. Sun, "Multi-stage holomorphic embedding method for calculating the power-voltage curve,"
IEEE Tran. Power Syst. , vol. 33, no. 1, pp. 1127-1129, 2018. [18]
Y. Liu, K. Sun, R. Yao, B. Wang, "Power system time domain simulation using a differential transformation method,"
IEEE Trans. Power Syst. , vol. 34, no. 5, pp. 3739-3748, Sept. 2019. [19]
Y. Liu, K. Sun, "Solving Power System Differential Algebraic Equations Using Differential Transformation,"
IEEE Trans. Power Syst. , in press. [20]
Y. Liu, K. Sun, " Differential Transformation of a Motor Load Model for Time-Domain Simulation," arXiv preprint arXiv: , 2019. [21]
N. Duan, K. Sun, "Power system simulation using the multi-stage Adomian Decomposition Method,"
IEEE Trans. on Power Syst ., vol. 32, no. 1, pp 430-441, Jan. 2017. [22]
G. Gurrala, D. L. Dinesha, A. Dimitrovski, et al, "Large multi-machine power system simulations using multi-Stage Adomian Decomposition,"
IEEE Trans. on Power Syst ., in press. [23]
E. Abreut, B. Wang, K. Sun, "Semi-analytical fault-on trajectory simulation and its application in direct methods," in Proc.
Power and Energy Society General Meeting , 2017, pp. 1-5. [24]
B. Wang, N. Duan, K. Sun, “A Time-Power Series Based Semi-Analytical Approach for Power System Simulation,”
IEEE Trans. Power Syst. , vol. 34, No. 2, pp. 841-851, March 2019 [25]
E. Pukhov, G. Georgii, "Differential transforms and circuit theory,"
International Journal of Circuit Theory and Applications , vol. 10, no. 3, pp. 265-276, Jun. 1982. [26]
I. H. A. Hassan, "Application to differential transformation method for solving systems of differential equations,"
Applied Mathematical Modelling , vol. 32, no. 12, pp. 2552-2559, Oct. 2007. [27]
V. S. Erturk, Z. M. Odibat, S. Momani, "The multi-step differential transform method and its application to determine the solutions of non-linear oscillators,"
Advances in Applied Mathematics and Mechanics , vol. 4, no. 4, pp. 422-438, Aug. 2012. [28]
Z. Odibat, S. Momani, "A generalized differential transform method for linear partial differential equations of fractional order,"
Applied Mathematics Letters , vol. 21, no. 2, pp. 194-199, Feb. 2008. [29]
Y. Liu, K. Sun, "Differential Transformation of Nonlinear Power Flow Equations," arXiv preprint arXiv:2004.08017arXiv preprint arXiv:2004.08017