A Dynamized Power Flow Method based on Differential Transformation

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.


I. INTRODUCTION
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 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 T designs 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. PROBLEM STATEMENT
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 Ybus 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).
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= [e1,…, eN] T and f= [f1,…, fN] T are respectively the real and imaginary parts of the bus voltage phasor; N is the total number of buses; λ is the loading parameter. 0 ( , )   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.

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).
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:

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.
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).
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.
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 C1 and C2 are constants and vl(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 2: For Formulation 1, its purpose is to characterize how the power changes with time, i.e., the power increases with time when C1>0 and decreases with time when C1<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 C2>0 and decreases with time when C2<0. It can also be used to trace either monotonical or non-monotonical curve segments such as (a), (b) and (d) in Fig. 1.
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 pi sp and Δpi are constants and λ= λ(t) is a variable, therefore their transformations are: pi sp  pi sp (k) and Δpiλ  ΔpiΛ(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): Then, apply the rule (7-i) and (7-ii): Finally, using the rule (7-i), the RHS of (20) is obtained.

2) 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)=C1(k) holds. Then, it can be further written in (24) after replacing k by k-1 and using the definition of (k).
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.
where Y(k) 2 1 N   and Λ(k)   are variables representing the DT of y and λ respectively; aP,i, aQ,i, aV,i aE,i, aF,i From the Proposition, DTs (9c) of the nonlinear power flow equation satisfy formally linear equation (10a) with matrices Agy, Agλ, and Bg 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.
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 ( ), ( , , Calculate matrices Agy, Agλ, Ahy, Ahλ using (31) 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. CASE STUDIES
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 C1 and C2 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. 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. 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. 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 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 handle 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.

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

APPENDIX
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.
Proof of Lemma: Therefore, (32) holds with a, b and c given below.


According to the Lemma, the three terms are rewritten as: