Analytical solution to swing equations in power grids

Summary

Objective

To derive a closed-form analytical solution to the swing equation describing the power system dynamics, which is a nonlinear second order differential equation.

Existing challenges

No analytical solution to the swing equation has been identified, due to the complex nature of power systems. Two major approaches are pursued for stability assessments on systems: (1) computationally simple models based on physically unacceptable assumptions, and (2) digital simulations with high computational costs.

Motivation

The motion of the rotor angle that the swing equation describes is a vector function. Often, a simple form of the physical laws is revealed by coordinate transformation.

Methods

The study included the formulation of the swing equation in the Cartesian coordinate system, which is different from conventional approaches that describe the equation in the polar coordinate system. Based on the properties and operational conditions of electric power grids referred to in the literature, we identified the swing equation in the Cartesian coordinate system and derived an analytical solution within a validity region.

Results

The estimated results from the analytical solution derived in this study agree with the results using conventional methods, which indicates the derived analytical solution is correct.

Conclusion

An analytical solution to the swing equation is derived without unphysical assumptions, and the closed-form solution correctly estimates the dynamics after a fault occurs.

Keywords:

Simulation and modeling – Eigenvalues – Inertia – Rotors – Differential equations – Electrical faults – System instability – Mechanical energy

Introduction

Electric power loads are expected to be fulfilled continuously in modern society, and when a load is not satisfied it is termed an “event”. The infrastructure to generate and transport electricity to end-consumers is called a power system. In the United States, this infrastructure comprises 19,023 individual, commercial generators (6,997 power plants) [1], 70,000 substations [2], and 360,000 miles of lines [3]. The number of power electronic devicess is more than a million, and non-anticipated losses of system components inevitably occur. The Federal Energy Regulatory Commission in the United States regulates the interstate transmission and wholesale of electricity. According to a report submitted to them [4], the 1-in-10 standard is a widely used reliability standard across North America. To meet this standard in a large-scale network with many system components, the power system must be able to withstand sudden disturbances (such as electric short circuits or non-anticipated loss of system components). It should be noted that most disturbances (including the failure of components) do not lead to an event. When a disturbance occurs, a governor regulates the speed of a machine to adjust the output power of a generator according to the network conditions. In general, the timeframe of governor action is approximately 0.1–10 s [5]. Therefore, it is necessary to assess if the power system is stable approximately 10 s after a disturbance occurs, which is the subject of the transient stability assessment.

Viewed overall, power systems consist of mechanical and electrical systems that obey energy conservation and Kirchhoff’s laws, which are integrated as the so-called swing Eq [5],

This paper is structured as follows: the first section lists three approaches, assumptions, and limitations; the second section formulates the problem in a different coordinate system than the polar coordinate system (PCS), and discusses the differences between the two; the third section lists the solution process to solve the reconstructed problem; the fourth section presents the analytical solution; the fifth section shows a set of examples; and the sixth section lists the conclusions and future studies.

Approaches for solving the swing equation

Coupled oscillator model

If the complexity associated with power systems is ignored, there might be a problem such as the swing equation. The first approach is to ignore the complexity of power systems, and to apply the knowledge from a different domain of science. It was found that in an oscillatory motion, if the frequency is spread more than the coupling between the oscillators, each oscillator runs at its own frequency. Otherwise, the system spontaneously maintains synchronization [6]. The field of synchronization in networks is reviewed in [7]. The Kuramoto model [8] provides an analytic function as follows: dδiKdt=ωiK+ςKN∑j=1Nsin(δiK−δjK).

For the coupled oscillator model (COM), a system of particles aims to minimize the energy function E(δK)=∑{i,j}aij[1−cos(δiK−δjK)]−∑kωkKδkK≅12∑{i,j}aij(δiK−δjK)2−∑kωkKδkK. Term E(δ^K) features a phase-cohesive minimum with interacting particles no further than a certain angle γ^K. A solution has cohesive phases if there is an angle γ^K between 0 and π/2 that is the maximum phase distance among all the pairs of connected oscillators, i.e., |δ_i−δ_j| ≤ γ^K for a line connecting two nodes i and j. In [9], it is demonstrated that the mechanical analogy applies to the electric power grid to yield the phase cohesiveness. This approach involves a very light computation cost, but the applicability may be limited due to the difference between a generic network and the electric power grids.

Lyapunov stability

The second approach is to solve the equation for a simple problem without ignoring the nonlinear feature of the equation. To construct a simple problem, we consider a system with a single machine and infinite bus where the loss is ignored (γ_ik = 0). In physics, we often define zero potential at an infinitely remote location, so that the existing field is not affected. Similarly, we define an infinite bus so that the terminal voltages do not vary with the internal voltage angle (and therefore with time). However, the terminal voltages may change with the network condition. Then, the nonlinear component |y˜ik|EiEksin(−θk+δi−γik) is rewritten as p_maxsinδ_i, where pmax(=|y˜ik|EiEk) is a constant dependent on the network condition. There are three different conditions regarding the condition for a disturbance: pre-fault, post-fault, and on-fault (See Fig 1).

Fig 1 shows that the corresponding p_max decreases in the order of pre-fault, post-fault, and on-fault conditions (pmaxpre−fault>pmaxpost−fault>pmaxon−fault). In the pre-fault condition, the operation point is determined where the sine curve intersects with a horizontal line pimech−giiEi2 (Point a in Fig 1(B), note that losses are ignored, i.e., g_ii = 0). At a, the generator is synchronized with the system frequency. During a fault, the sine function drops down the on-fault sine curve, but the internal voltage angle remains unchanged (δ_pre), because the angle must be continuous over time (a → b). Therefore, the electric power output is less than the mechanical power input, and the difference must be stored in the rotor in terms of mechanical energy to meet the energy balance. The internal voltage angle increases accordingly (b → c), which makes the rate of the internal voltage angle deviate from the system frequency. The fault will shortly be monitored and cleared, and the sine curve follows the post-fault curve that the internal voltage angle follows (c → d). Similarly, the internal voltage moves to a new curve at δ_clear in which the mechanical power input is less than the electric power output. The angle will increase further until the stored energy is exhausted at δ_p (d → e). The integration of the left-hand side between two synchronized points of the swing equation (after ignoring the damping term) yields zero (Area 1 = Area 2). This is called the equal-area criterion. At δ_p, the generator recovers the synchronization with the system frequency, but it cannot maintain the operation point because the electric output is greater than the mechanical power input and the stored energy is already exhausted. The rotor angle decreases until δ_n, which is set by the equal-area curve (Area 3 = Area 4) and the rotor, swings between the two angles [δ_n, δ_p] if no damping exists. With the damping, the swing range decreases and settles down to δ_post, where the mechanical power input and electric output are balanced.

Unfortunately, this equal-area criterion cannot generally be extended to a large-scale network with multiple machines. According to Lyapunov stability in nonlinear dynamic theory, one can tell if a system will be stable in the future if a Lyapunov function is identified to meet some conditions [10]. While there is no rule to construct a Lyapunov function, DM is developed under various assumptions [10–17] that often includes a lossless network (zero damping), no consideration of reactive power, and constant voltage magnitudes that are difficult to justify physically. Under these physically unjustifiable assumptions, the system is stable if one can find a Lyapunov function. On the other hand, the failure to find such a function (or the nonexistence of the function) does not necessarily mean the system will be unstable. Rather, the stability would be undetermined. Even though the advantage of this approach involves relatively inexpensive computation costs, it yields a stability region subset in addition to the physically infeasible assumptions.

Both the COM and the DM ignore the losses and variations of the voltage magnitudes. Further, because both assumptions are not physically justifiable, the application of the approaches is either highly limited or a conservative application is performed.

Numerical approach

The third approach is to accommodate all the details of the swing equation and Kirchhoff’s laws, and conduct a numerical computation. Often, a numerical integration or differentiation approach is applied to solve nonlinear differential equations and recent advancements in fast computation makes the numerical approach an attractive option [18], [19]. For the swing equation, a numerical method (such as the Euler or Runge–Kutta method) is used to find the temporal changes of internal and terminal voltages [20], [21]. The advantage of this approach is the capability to solve any complex function. However, in addition to the high computation costs, the numerical method involves a truncation error. For example, the errors associated with the Euler and modified Euler methods are in the order of ϑ(‖Δx‖²) [20], and that of the Runge–Kutta method is ϑ(‖Δx‖⁴) [22]. If there is a subset of modes that diverge exponentially, but in which the magnitudes are initially very small, the modes may be underestimated. The numerical approach also requires precise estimation of the model parameters. It should be noted that numerical simulation did not capture the impact of the 1996 WSCC (Western System Coordinating Council) system outage [23].

Write down the problem

Feynman’s algorithm

According to Dr. Thomas (an Emeritus Professor at Cornell University), Feynman’s algorithm for solving hard problems is a process of thinking about problems in a different way. The swing equation, Mid2δidt2+Didδidt=pimech−giiEi2−|y˜ik|EiEksin(−θk+δi−γik), is formulated in the PCS because the equation describes the angular motion of a rotor inside a generator. The PCS provides a concise form of the equation of the motion; the differential terms are multiplied with scalar constants that result in a linear form for the left-hand side; and on the right-hand side, the magnitudes of internal voltage and terminal voltage are linear. However, the nonlinearity involved in voltage angles makes the swing equation difficult to solve. In this paper, instead of the conventional approach using the PCS, we construct the problem in the Cartesian coordinate system (CCS) and aim to solve it without physically unjustifiable assumptions.

Observations in transient studies

We begin with a list of facts that are observed in transient studies.

O1. E (internal voltage of a generator) at Ibus (see Section “Network flows” for its definition) is constant [24], which is a general assumption applied to the transient studies in power systems, O1=‖[1−xi2+yi2/Ei⋯1Ei2xi(dxidt)+2yi(dyidt)⋯]T‖F.

O2. O2 = |O_i(t)—O₊⁰⁺| where O_i(t) = (−qi+biiEi2)/Mi+(dδi/dt)2; b_ii is the imaginary part of admittance; the superscript 0+ represents “immediately after contingency”; and t > t₀, i.e., the sum of the reactive power injection and the rotor speed changes slowly over time strictly after a disturbance occurs before losing synchronization. The square of the rotor speed (dδ_i/dt)², the square of the rate to deviate from the synchronization, is in general smaller than (−qi+biiEi2)/Mi while maintaining synchronization. According to O2, while maintaining synchronization, |O_i(t)| is a slowly varying function that follows the Karamata representation theorem [25]: |O_i(t)| =|−qi+biiEi2Mi+(dδidt)2|=exp[η(t)+∫t1tμ(y)ydy] where {limt→∞η(t)=η∞limt→∞μ(t)=0 for t ≥ t₁ (> t₀). The variation of O_i(t) over a short time period is negligible so that O_i(t) = (−qi+biiEi2)/Mi+(dδi/dt)2≈(−qi0+biiEi2)/Mi+(dδi/dt)2|t=0 = O_i(t₀). O2 is defined as |O_i(t)–O_i(t₀)|.

When a disturbance occurs it is possible for some generators to stray slightly off perfect synchronization. However, the rate of deviation from synchronization is always kept small before losing synchronization. The impacts of these approximations are discussed in a later section.

Write down the equation in CCS

Definition of variables

A natural choice of variables in CCS would be the real and the imaginary components of voltages. The choice leads to E_iE_ksin(−θ_k+δ_i−γ_ik) = −cosγ_ik(E_icosδ_i)(E_ksinθ_k)+cosγ_ik(E_isinδ_i)(E_kcosθ_k) −sinγ_ik(E_icosδ_i)(E_kcosθ_k)+sinγ_ik(E_isinδ_i)(E_ksinθ_k). A better choice is to rotate the voltage angles by the phase angle of the line between the internal and terminal buses (δ_i−γ_ik), and the choice yields only two terms by EiEksin(−θk+δi−γik)=−xivyk+yivxk where x_i = E_i cos(δ_i−γ_ik), y_i = E_i sin(δ_i−γ_ik), vxk=Ekcosθk, and vyk=Eksinθk. For CCS, two variables are necessary (x and y) instead of one (δ) to describe the angular motion of internal voltages, and an additional equation is imposed to preserve the constant internal voltage magnitudes xi2+yi2=Ei2.

Load modeling

The development of a new load model is beyond the scope of this study. Instead, we attempt to integrate the existing load models for our proposed formulation. A widely used load model is outlined in [26]. Loads located at a same bus are integrated into a single load, and this load is separated and modeled into four subgroups: (Category I) an induction load, (Category II) a frequency-dependent or a time-dependent load, (Category III) a load with constant impedance or a load with constant current, and (Category IV) the remaining loads. While Categories I and II are integrated in the swing equation, the admittance matrix Y is modified to incorporate Categories III and IV in terms of the voltage–current relationship.

Category I: An induction load. This has nonzero inertia, meaning it appears in the swing equation as a synchronous machine that remains in the swing equation. The load is modeled similarly to that of a generator, or more specifically, to that of a negative generator.

Category II: A frequency dependent or a time dependent load. A frequency dependent load is dj=dj0+mj(dδj/dt) [8], where the first component dj0 represents a load with a fixed impedance, and the second term m_j(dδ_j/dt) is the frequency-dependent load. If a synchronous machine is located at the same bus as the load, the coefficient of the first time-derivative term is the sum of both the damping and frequency-dependent terms. Similarly, a time dependent load is dj=mj(ddj/dt)+dj0 [27]. For this type of load, the constant term is integrated into Category IV.

Category III: A load with constant impedance or with constant current. These are integrated in il=ilcc+ylcivl, where i_l^cc and y_l^ci model the constant current and the constant impedance at Bus l, respectively. This expression makes it possible to convert these loads into the diagonal shunt element in the admittance matrix Y_bus (y_l) or a constant in iL=diag(yLci)vL+iLcc.

Category IV: a remaining load. The remaining load introduces nonlinear characteristics, and it cannot characterize load characteristics on a constant voltage node in the system due to non-dependency on the voltage angle at the node. The BIG model in [28] and [29] is an attempt to integrate the load as a linear model (similar to loads with constant impedance or current), and show a good match for static loads. Interested readers may wish to read a literature survey on modeling loads in [30]. The type of load may be interpreted with a Taylor series expansion in terms of voltages near the operation point: Il=sl*vl*=(il0vl+yl¯|vl|2)*vl*=il0+y¯i*vl. Combined with the load modeling constant current and constant impendence, the current is expressed in terms of voltage iL=iL0+diag(yL0)vL.

Network flows. To establish a set of equations, an electric power grid is redefined, and a set of buses is introduced to model the internal voltage of a synchronous machine, Ibus; in other words, a Kbus is directly connected to an Ibus. Note that synchronous machines include induction motors, frequency- and time-dependent loads, and synchronous generators (but not asynchronous renewable generators). No Ibus is connected to multiple Kbuses, while a Kbus can be connected to multiple Ibuses. If two Kbuses are directly connected, an Mbus is inserted between them. Fig 2 illustrates the definitions of the Ibus, Kbus, and Mbus. Fig 2(A) shows the one-line diagram of a modified IEEE 3-bus system with three generators, and Fig 2(B) is the one-line diagram of the system modeled in this study. The generators are modeled as Ibuses (in red); the top two buses where the generators are directly connected are modeled as Kbuses (in green); and the bottom bus that is not connected to a generator is modeled as an Mbus (in blue), respectively. To prevent two Kbuses from being directly connected, another Mbus (vertical bus in blue) is inserted. Because all the buses in the original network remain in the model, and because NI of Ibuses and additional Mbus are introduced, the number of buses for this network model is always greater than that for the original network. It should be emphasized that unlike DM, in which the network reduction is applied to make the problem simple, this study preserves the topology of the network.

Suppose a disturbance occurs between t = 0 and t = 0⁺ (immediately after the disturbance). During the disturbance, the angular motion of the physical rotor must be continuous; that is, instantaneous change in the angular motion is zero. Unlike rotor angles, voltages do change abruptly, while Kirchhoff’s laws are obeyed both before and after the disturbance. Using Kirchhoff’s laws, the voltages are computed to meet the real and the reactive power balances at Kbus and Mbus immediately after the disturbance. After the disturbance, the angular motion of the rotors and voltages at all the buses correspondingly adjusts to meet both the swing equations and Kirchhoff’s laws. The admittance matrix of the redefined grid that accommodates Categories III and IV loads is in a block structure associated with both Ibuses and Mbuses because there is no direct connection between them. For the Kbuses and the Mbuses in the network model shown in Fig 2(B), Ohm’s law finds i = Yv where the currents i and the voltages v at the buses are also expressed in terms of loads, i.e., iKM=iKM0+diag(yKM)vKM. Note that induction, frequency-dependent, and time-dependent loads are not included in the current–voltage relationship. With no direct connection between Ibus and Mbus, Ohm’s law over the network modeled in this work (Fig 2 (B)) is (iIbusiKbusiMbus)=[YIIYIK0YKIYKKYKM0YMKYMM](vIvKvM). Because an Ibus is only connected to one Kbus, the inclusion of the buses is a radial extension of the original network. Therefore, the currents on the Kbus and the Mbus at a given voltage must be invariant to the choice of network model:

Revisiting swing equation

Typically, swing equations are written in the PCS, as this is convenient for describing angular motion. The PCS inevitably introduces the exponential function (i.e., sinusoidal function) to express the real power injection, which makes it difficult to solve analytically. In a DC power flow model [26], we often linearize the sinusoidal function by taking sinθ ≈ θ and cosθ ≈ 1 as θ ≈ 0. However, in the swing equation, δ to describe the angular motion is not small enough to apply the approximation. The power injection from Ibus i to Kbus k is as follows:

With (4) and (5), the swing equation in CCS becomes

In comparison to the swing equation in the PCS, (6) is a heterogeneous nonlinear second-order differential equation with multi-variables; the time derivatives are multiplied with the variables; and the last term is nonlinear. Additionally, the conditions related to the constant internal voltage magnitude and to its derivatives are also nonlinear. If one combines the conditions with (6), the resulting equation becomes highly complex. Instead, we derive

Eqs (6) and (7) lead to

Furthermore, the voltages at Kbus are expressed in terms of w_I as shown in (2)

(14) is a constrained homogeneous linear first-order differential equation with multi-variables. There is no known solution to the constrained differential equation. Fig 3 illustrates our approach to solving the swing equation. The arrows with solid lines refer equivalent, and those with dashed lines represent similar in the proximity proportional to the relaxation applied. S1 is the solution to the conventional swing equation formulated in PCS, which is identical to S2 in CCS and S3 because (8) is equivalent to (10) with the constraint of O2 = 0. Equality constraints are identical to the inequality constraints if the upper bounds are zeros, i.e., S3 = S4. S5 is the solution to the problem relaxed by ϑ(Δ_E). Note that t_E + δt is the first time when any of the relaxed constraints is binding, and that δt is strictly positive. In the range of [0, t_E], no constraints are binding, which makes S6 involve nonbinding constraints. In the optimization theory, the solution to a constrained problem equals that of the same problem without the constraints that are not binding. Therefore, S7 is the same as that to S6 even though no constraints are taken into consideration in solving S7. The process finds S7 is in the proximity of S1 by ϑ(Δ_E) if δt≪t_E and Δ_E is kept small enough to ensure the validity of the solution. The justification is outlined in the section describing the validity region with projections so that the error remains small. This concludes the first step of Feynman’s algorithm, write down the problem.

Think hard: Solve the swing equation

Initial condition

To solve a first-order differential equation, a set of initial conditions are required. δ represents a physical angular motion that is continuous in time; i.e., δi0+=δi0. When a disturbance occurs, the mechanical power input does not change instantaneously, but the electric power output pi→kelec exhibits a sudden but finite change. The change is represented by a step function f(x∈R1)=∑imiχAi(x) where m_i is a finite real number and χ_A is an indicator function over an interval A, χA(x)={1ifx∈A0otherwise. Because the time interval between 0 and 0⁺ has a finite length, the Lebesgue integral of pi→kelec is ∫pi→kelecdt=ml(A), where l(A) is the length of the time interval between 0 and 0⁺, ∫t=0t=0+pi→kelecdt≃0. The value of ω_i (= dδ_i/dt) immediately after the disturbance occurs is computed from the integration of the swing equation over the time interval ωi0+=ωi0. Similarly, for the frequency- or time-dependent loads, ω˜i0+=ω˜i0. The rotor angle δ (in terms of time) is differentiable, and its first derivative ω is bounded; hence, the rotor angle δ is Lipschitz continuous over time.

General solution to the swing equation

The general solution to (14) that satisfies dzGdt=TzG is as follows:

It is possible to represent u in terms of a real-valued function as follows:

Solution to meet initial conditions

We claim that the following function is the solution to (14):

It is straightforward to prove this claim as follows:

With a set of vectors ξ, ξ_l = Ψ_lu(t = 0), β is determined from the following least square process:

The analytical solution to the swing equation that satisfies Kirchhoff’s laws is

Write down the solution

Computational complexity

There are several steps to determining the complexity: finding voltages immediately after the disturbance, decomposing eigenvalues of T, and computing β. The first process has the same complexity as solving power flow equations ϑ(Nb^1.5) [32]; the second and last processes involve ϑ[(NI+ND)³] [31], because the dimension of T is (4NI+2ND) × (4NI+2ND). Therefore, the first process is predominant if Nb > (NI+ND)². Even though additional computations are necessary at τ, the first process does not need to be solved because the voltages are readily available from (2). Therefore, the additional computation may not necessarily increase the computational cost significantly.

Stability assessment

The conventional stability assessment based on COM or DM certifies when a state is eventually stable. We propose another stability assessment—that the eigenvalues of T play a key role.

• Type I: stable if the real parts of all the eigenvalues are non-positive or Θ = 0 (no disturbance).

• Type II: operationally stable if the largest positive real eigenvalue λ_m is small enough such that 1/λ_m ≤ T_op, where T_op is a time scale where transient stability is concerned. Typical time scales of transient stability studies are between sub-seconds to tens of seconds [33].

• Type III: operationally stable if the coefficients corresponding to the terms in the base solution with positive real eigenvalues are small enough for the system to remain stable within T_op.

• Type IV: unstable if the system divulges rapidly.

In this assessment scheme, the positive real parts of eigenvalues and the corresponding coefficients play a key role. Apart from Type III (in the presence of a disturbance), assessment is possible with the eigenvalues of the T matrix. If no positive eigenvalues exist, the state is eventually stable.

Validity region

The general solution of (20) to (10) is obtained without considering the constraints. It would be ideal to solve the ordinary differential Eq (10) with the constraints–DAE. DAE integrates the differential variables and the algebraic variables. While DAE is complete, to the best of our knowledge, no analytical solution has currently been identified. Therefore, the approach in this study is to identify the validity region where the constraints hold. To discuss the change of the constraints (drift-off phenomena), if we do not explicitly consider them while solving the differential equation (such as P7 in Fig 3), suppose we want to solve the following pendulum problem, dx/dt = u;du/dt = −λx;dy/dt=v;dv/dt = −λy−1;x²+y²−1 = 0. One may eliminate λ, and find a differential equation as follows: u′x+u2+v2−yx2+y2=0;v′+1y+u2+v2−yx2+y2=0 with the following constraints: x²+y² = 1 and xu+yv = 0. The numerical solution to the differential equation is evaluated without considering the constraints explicitly. Fig 4(A) shows the drift-off phenomena of the constraints, and the errors associated with x²+y² = 1 and xu+yv = 0 grow quadratically and linearly, respectively [34].

Eq (10) is analogous to constrained mechanical systems that can incorporate the repeated projection of a numerical solution onto the solution manifold; projections on position constraints and on velocity constraints for improving the stability (see Fig 4(B)). The error remains negligible with well-developed projections; hence, the solution (20) to (10) is valid without consideration of the constraints. However, the projections are not a linear process, which makes it difficult to integrate them into (10). Instead, we define a validity region where solution (20) is valid, and the projections are made at the boundary of the validity region to compensate for the drift-off. The induced errors and the approach to compensate the errors are discussed in the following section.

Internal voltages, E_I. Suppose, in a time range of [0, Δt], the internal voltages deviate from the nominal values significantly, which affects z in (20), meaning the constraints do not hold. Because (20) yields the analytical expression of x_I and y_I in terms of t, one can evaluate a projection vector Δ_con that includes the error in E_I (position projection, p⌣Ei=1−xi2+yi2/Ei) and the error in the first derivative (velocity projection, v⌣Ei=2xi(dxi/dt)+2yi(dyi/dt)/Ei), i.e., Δcon=‖(p⌣ETv⌣ET)T‖2. The constraints hold if the 2-norm of the projection vector does not exceed a threshold Δ_E.

Error in T. The matrix T in (10) is exact except for one element each row by an approximation that O_i(t) = (−qi+biiEi2)/Mi+(dδi/dt)2=(−qi0+biiEi2)/Mi+(dδi/dt)2|t=0 = O_i⁰⁺. There might be a discrepancy between the true z and the estimated z due to the nonnegligible O2. The impact will affect two components: (1) the inaccuracy in estimating T and the eigenvalues of T, and (2) the error in the coefficients according to Φ, δΦ.

The impact of O2 in the inaccuracy is as follows: When the change in O_i(t) is in the order of ε_T, the impact propagates to T, D_l, and correspondingly Φ and Ψ. The change in T due to O_i(t) affects the eigenvalues of T that are in the block-diagonal in D_l. Suppose λ is a simple eigenvalue of T and u_L and u_R are the corresponding left- and right- eigenvectors. ϑ(ε_T) changes in T can induce ε_T/s(λ) in the eigenvalues [31],

With the definition of ρ and σ, (30) becomes:

Because a time-varying O_i(t) determines ε, one can compute the τ_ε that the relative error in T remain a threshold Δ_T, i.e., ε = ‖δT‖_F/‖T‖_F≤Δ_T for all t≤τ_ε. Let y_z be the first derivative of z with respect to time, then (14) yields: y_z = Tz₀+b where y_z = dz/dt|_{t = 0} and z = z₀+ydt. In the time range, there might be an error in T involving the error that propagates in y_z and z as follows:

Beyond the validity region. The boundary of the validity region is defined either when the constraints do not hold, or when the error in T exceeds the threshold. At the boundary of the validity region, the values for z and T are updated as described. With the updated values, (14) still holds because the errors in (14) and the constraints remain less than the threshold. Therefore, the problem in (14) will be solved with the initial condition that is the updated values for z. It is necessary to ensure consistency among the initial values; hence, a consistent initialization problem. This problem has been studied widely; (1) a small artificial step with the backward Euler method [35], [36], (2) Taylor series expansion [37], [38], and (3) graph theoretic algorithm to obtain the minimal set of equations for differentiation to solve for consistent initial values [39], [40].

Eqs (12) and (13) and the constraints are rewritten as follows:

The cardinalities of f and g are 2NI+2ND and 2NI, respectively. By the definition of the variables in (33), u and s are classified as differential variables and algebraic variables. Accordingly (33) is called a DAE. At a given w_I⁰ and u_I⁰, the Newton-Raphson method finds:

The term in (34) is compensated to yield a consistent initial point immediately beyond the validity region. With the consistent initial point and updated T, the analytical solution in (20) is identified. Most components inside the curly bracket in (34) are constant. Therefore, only a few visits to the boundaries of the validity limits does not increase the computation time significantly if efficient rank-update techniques are employed.

Insight on the dynamics

As in (12), T has a block structure, and it can be broken into two submatrices as follows:

The first term is invariant with the operation point of specific contingency, but the second term varies. We consider the eigenvalue decomposition of T as consecutive applications of the eigenvalue update and down-date of T_op from the eigenvalues of T_sys; that is, the sum of multiple rank-1 [41] or rank-2 [42] updates (+) and multiple down-dates (-). The impacts of T_op are to shift the eigenvalues of T_sys, which is bounded by the Bauer-Fike theorem [31]: minλ∈λ(Tsys)|λ−λ^|≤(‖UΦ‖p‖UΦ−1‖p)‖Top‖p where UΦ−1TsysUΦ=diag(λ^1,⋯,λ^N). If the diagonal elements in T_op are all zeros (or negligible), the Gershgorin circle theorem can be applied [31]. The Gershgorin disk is isolated from the other disks so that a disk contains precisely one eigenvalue of T. Motivated readers may find the process in [43]. Note that T_sys only depends on the system, and T_op varies with the dynamics. The following section outlines the terms affecting the stability of the system.

Reactive power, qicinqic+biiEi2: This term negatively affects L_con, shifting the eigenvalues to the right. Reactive power supports voltages, which makes it difficult for a system to settle into new voltages. However, the reactive power injection is very small in comparison to the biiEi2 term; therefore, its impact on the transient stability is highly limited.

Unsettled mechanical power, pimech in pimech−giiEi2: These terms are associated with the mechanical power in L, and have the same magnitudes with opposite signs; hence, their impacts on shifting the real eigenvalues of T cancel out. However, their impacts are on shifting the imaginary components of the eigenvalues. It is intuitively correct that unsettled mechanical power enhances the oscillating motion. In many cases, the resistance of an internal generator model is zero, meaning g_ii is zero. Unlike the reactive power injection, pimech plays a key role in stability assessment. A synchronization condition based purely upon the power injections is peoposed in [9].

Inter-voltage sensitivity matrix, H_KI: As the sensitivity of H_KI increases in a positive direction, the voltages at Kbus swing tightly together with the generators. The negative sensitivities in H_KI imply that the voltages at Kbus swing against the angular motion of rotors. They act as a dragging force to stabilize the system. The impact of the inter-voltage sensitivities is to shift in the same direction as the sign of H_KI.

Damping coefficient. The damping term appears in the lower diagonal in T. Because the damping terms are non-negative, the impact is to shift the eigenvalues to the left. Therefore, the damping terms help to stabilize the system against disturbance.

Inertia. The terms reactive power, unsettled mechanical power, voltage sensitivity H_KI, and damping are all scaled by the inertia constant M_i. Regardless of their signs, the amplitudes are normalized in terms of inertia. The impact of inertia becomes clear in TDS in a way that agrees with this observation.

Loads. Loads are classified into three categories: synchronized induction load, frequency- or time-dependent load, and remaining loads. The remaining loads affect the sensitivity matrix H_KI, which is a part of the L, L˜, and T_op matrices. The frequency- and time-dependent loads, and the synchronized induction loads, are taken into consideration in the swing Eq (14). It is noteworthy that the T matrix has a block diagonal structure between the two loads; therefore, the Eigen-space of each block is independent so that their subspaces are orthogonal.

Illustrative examples and discussion

Simulations are performed for various IEEE model systems (IEEE 9, 14, 30, and 118-bus systems). To compare the results with both a coupled oscillator model [9] and DM, which assume a lossless system (γ = 0 for all the lines) and constant voltage magnitudes, the resistance components are all ignored. A phase cohesiveness for synchronization is also introduced: |θ_i−θ_j|≤m∈[0,π/2]. Electric power is generated at Ibus and injected into the grids; therefore, the angles at Ibus (δ_I) are greater than those at Kbus. (δ_i≥θ_k). The power flow from Ibus i to Kbus k (power injection at PV buses of the original network) is proportional to sin(δ_i≥θ_k) [26]. Therefore, the phase cohesiveness is equivalent to the constraints imposed on the injection between two directly connected oscillators. We found that as the system becomes unstable, the maximum angle differences after the disturbance are significantly higher than those in the stable case. Even though the values for the dynamics are not listed in [9], the phase cohesiveness for synchronization finds the trend correctly for the disturbances we tested. Similarly, the synchronization condition is also checked, and it was found that there is a γ^K appearing in [9] to satisfy ‖L^†p‖_ε,∞≤sinγ^K for the chosen equilibrium points. Because there can be many equilibrium points, it can be difficult to analyze the stability region of each point [13]. For the sake of visual presentation, the simulation results on the IEEE 9 bus system are discussed in this paper. Fig 5 shows the one-line diagram of the system, and Table 1 lists the data relevant to this study. In the proposed network modeling in this study, there are three Ibuses, three Kbuses, and six Mbuses. The pre-fault power flows are computed using the unified method based on the Kronecker product [44], and the threshold values for Δ_T and Δ_E are 1% and 10%, respectively. The scaling factor h in this study is 0.1. All the numerical computations are performed using a Mac pro with two 2.93 GHz 6-core Intel Xeon processors.

The results are summarized in Table 2. In the table, 10% d₈ in the first column refers to 10% loss of loads at Bus 8; Fig in the second column is the figures associated with the event at the first column; Δδ_max under coupled oscillator column is the synchronization condition proposed in [9] where the threshold for the system is 0.129; ΔV_st under the DM column is the stability margin that is defined in the Appendix. For some cases, the stability assessments are undetermined and shown as “-“, because the certificate is not issued. Time under the time domain simulation (TDS) and proposed model columns represent the computation time for numerical computations.

Loss of loads

Three sets of cases are performed to simulate the loss of load at Bus 8; no loss, 10% loss, and entire loss. Prior to the loss of load, the system was in the steady state condition that was identified using the power flow study. When no disturbance occurs, the system should stay in the same steady state at t = 0. Immediately after the loss of load, a new operating point is found by solving (2) with the updated loads and w_I at the steady state, because the rotor angles do not change instantaneously.

The stability of the post-fault operation point is estimated using DM and the synchronization condition proposed in [9]. The trajectory during the transient state is numerically calculated in terms of TDS. The details of TDS and of DM based on an energy function are outlined in the Appendix. For the case with no disturbance, even though the eigenvalues of T are non-zeros, the coefficient Θ is zero (Type I) because the constant term in (20) is zero, meaning the system stays in the steady state at the pre-fault condition.

Slight change in the load at Bus 8

Fig 6 illustrates the trajectories of (A) rotor angles δ, of (B) voltage magnitudes E, of (C) O_i(t), and of (D) computed magnitudes of the internal voltages when a 10% loss of the load at Bus 8 occurs at t = 1 s.

Stability of the system. In the post-fault state, the maximum angle differences in voltage angles immediately after the disturbance is 0.116 rad, while the synchronization condition reported from [9] is 0.129 rad. As shown Fig 6(B), the variations of the voltage magnitudes are not significant, which indicates that a condition of both COM and DM holds. Because the maximum angle difference is less than the threshold, the stability assessment predicts the convergence to a stable state. The energy functions V(δ, ω) for DM are evaluated at all equilibrium points to check the stability of the system. Based on to the stability margin of 8.96 (= V_cr−V_cl > 0), a stability certificate is issued by DM. TDS is also performed with a time step Δt of 0.01 s. All three stability assessment approaches yield the same estimate—converging to a new equilibrium point.

The positive eigenvalues of T are small, and the corresponding Θ is numerically negligible (Type II). As shown in Fig 6(A), the rotor angles are all synchronous, and the voltage magnitudes quickly settle down to a new equilibrium point. The analytical solution and the numerical results from TDS are visually indistinguishable, as shown in Fig 6(A). For clear presentation, the voltage magnitudes and the rotor speeds are generated by the proposed analytical approach in Fig 6(B).

Validity region. In this study, the tolerance for ε (= ‖δT‖_F/‖T‖_F), Δ_T, is set to 1% for all the numerical evaluations. The value of ‖T‖_F in the case of 10% loss of the load at Bus 8 is 16.71; hence, the tolerance of the impact of O2 toward ‖δT‖_F is 1.67 × 10⁻² (= 1% × 16.71). The variations of O_i(t) are negligible, as shown in Fig 6(C)—the value of ‖δT‖_F is 1.21 × 10⁻⁴ (‖δT‖_F/‖T‖_F = 7.24 × 10⁻⁶) is ϑ(10⁻⁴). This means the impact of O2 is negligible, and the validity region extends to the entire time domain of the transient in this case. Fig 6(C) indicates that η(t) for the generators are constant after assuming μ(t) = 0. Therefore, |O_i(t)| follows the Karamata representation theorem (i.e., a slowly varying function), which also leads to the extended validity region. Fig 6(D) shows that the estimated voltage magnitudes vary within the threshold (Δ_E = 10%) of the nominal voltages, and that O1 is less than the threshold most times. This implies only a few times of crossing the validity regions in 10 s. It is noteworthy to mention that in most simulations, O1 and O2 are small.

Entire loss of the load at Bus 8

A significant change in the load at Bus 8 occurs at t = 1 s, and the generator dynamics are simulated. Fig 7 illustrates the trajectories of (A) δ, (B) the voltage magnitudes at the terminal voltages, (C) O_i(t), and (D) the estimated magnitudes of the internal voltages.

Stability of the system. After the loss of load, the stability of the system is tested; (1) the maximum angle difference across the lines is 0.113 (<0.129), (2) the stability margin is 9.21 (>0), and (3) TDS with Δt = 0.01 s. All three methods find the stability of the system, and converge to a new equilibrium point after the disturbance. Due to the loss of loads at Bus 8, the reactive power injected to the grid exceeds the demands, and the uncompensated reactive power makes the voltages at the terminal buses increase (Fig 7(B)). A condition of COM and DM that the terminal voltage magnitudes are constant holds marginally.

Fig 7(A) shows the dynamics of rotor angle for the entire loss of the load at Bus 8, and exhibits no discrepancies between the analytic approach and TDS. The positive eigenvalues of T are small, but the corresponding Θ is numerically negligible (Type II). Therefore, the proposed analytical solution also estimates the stability of the system correctly.

Validity region. Because the value of ‖T‖_F in the case of the entire loss of the load at Bus 8 is 17.61, the tolerance to the impact of O2 toward ‖δT‖_F is 0.18 (= 1% × 17.61). The variation in O_i(t) of three generators over 10 s is approximately 8, and the corresponding ‖δT‖_F is 9.88 × 10⁻² (‖δT‖_F/‖T‖_F = 5.61 × 10⁻³). Therefore, the validity region extends to the entire time domain of the transient if O1 is small. As shown in Fig 7(C), η(t) (after assuming μ(t) = 0) does not converge as t goes to infinity; therefore, the Karamata representation theorem may not be applicable. Consistent with this observation, O_i(t) is not slowly varying, but the impacts on T are within the threshold. O1 often (5 peaks) reaches Δ_E (= 10%), as shown in Fig 7(D). The projection method identifies the peaks (the boundaries of the validity regions), and the consistent initial points are evaluated using (34) for the application of (20), which increases the computation time.

Line fault: On-fault trajectory

At t = 1 s, a three-phase fault occurs on the line connecting Buses 7 and 8 near Bus 7 (red cross mark shown in Fig 4). Using the pre-fault power flow solution by applying the unified power flow analysis approach in [44], the parameters to formulate the swing equations are identified. Due to the continuation of the rotational movement, the rotor angle δ remain continuous in time regardless of the change in the network. The element in Y_bus corresponding to Bus 7 is increased to represent a high admittance to ground, and the voltage at Bus 7 collapses. With this modification, the on-fault voltages and the parameters for the swing equations in the on-fault trajectory are computed. Fig 8 illustrates the on-fault trajectories of (A) δ, of (B) E, of (C) O_i(t), and of (D) the estimated magnitudes of the internal voltages when the fault is not cleared.

Stability of the system. After the disturbance, Bus 2 is isolated from the system that is directly connected to Generator 2, and its voltage magnitude reduces to zero (Fig 8(B)). Because no load is located at Bus 2, the electric power injection becomes zero, and the electrical power input from the generator is stored in the rotor in the form of mechanical power (increased rotor speed). For the rest of the system, the change in the power supply by the isolated Bus 2 (effectively the loss of Generator 2) is compensated by the other generators. This leads to a change in rotor angle, as shown in Fig 8(A). Clearly, the reactive power injection changes abruptly (see Fig 8(C)).

The maximum angle difference for COM is 0.137 rad, while the value for γ in [9] is 0.129, which does not meet the synchronization condition. If the method fails to issue a certificate, the stability of the system is undetermined. However, the short-circuit makes Bus 2 disconnected from the rest of the network, and the remaining network different from the original network. Therefore, the prediction based on the synchronization condition may not be exact. The closest unstable equilibrium for DM is found, and the stability margin is found to be -34.87 (<0), which means the stability of the system is not certified. Similar to COM, the stability of the system is undetermined if a certificate is not issued.

Fig 8(A) shows large deviations in the results for Generator 2 between the proposed analytic solution and TDS after t = 1.5 s, but both indicate system instability. They both indicate that Generator 2 will lose synchronization quickly after the line fault, and be disconnected from the network (Type IV, unstable). In the on-fault condition, Generator 2 is isolated from the rest of the system. Therefore, as shown in Fig 8(A), the rotor angular motion becomes faster than motion of the other generators.

Validity region. Fig 8(C) shows the fault-on trajectories of O_i(t) if the fault is not cleared. O_i(t) for Generator 2 is not a slowly varying function, and it does not follow the Karamata representation theorem because η(t) divulges for Generator 2, meaning η(t) is not bound. ‖T‖_F is 7.88, and the threshold for ‖δT‖_F is 7.88 × 10⁻² (1% × 7.88), but the impact of O2 is 1.00, which is greater than its threshold. Fig 8(D) lists the estimated magnitudes of the internal voltages. Shortly after the line fault, the internal voltages at Bus 2 increase too quickly, and O1 becomes large. Therefore, the validity region is limited quickly after 1 s, and neither O1 nor O2 are small during the on-fault trajectory.

Line fault: Post-fault trajectory

To prevent the loss of synchronization, the fault should be cleared. At t = 1.1 s, the fault is cleared by opening up the circuit breakers of the line between Buses 7 and 8. With the updated topology, the Y_bus is updated; accordingly, a new operation point is identified. However, the rotor angular motions are continuous. Fig 9 show the post-fault trajectories of (A) rotor angle, and (B) estimated magnitudes of the internal voltages, respectively. Fig 9(A) exhibits that the analytical approach does not predict the rotor dynamics properly after 2.5 s, and the discrepancy increases with time. The TDS has a time step of 0.01 s.

Stability of the system. The maximum angle difference in voltage angles immediately after clearing the fault is 0.155 rad, which COM concludes is the undetermined stability of the system. The closest unstable equilibrium point is observed, and the stability margin based on DM is -1.98 (<0), which also renders DM unable to estimate the stability of the system. There is a discrepancy in the results of the rotor angles between the proposed analytic method and TDS. This discrepancy increases after t = 3 s as the system evolves over time.

Validity region. The value of ‖T‖_F immediately after clearing the fault is 11.43; hence, the tolerance to the impact of O2 toward ‖δT‖_F is 0.114 (= 1% × 11.43). The impact of O2 is 1.79 (‖δT‖_F/(‖T‖_F = 0.156), which is beyond the threshold of 1% after clearing the fault. Fig 9(D) indicates that the voltage magnitudes suddenly increase after 2 s, meaning O1 is not small after 2 s. The validity region boundaries are identified with the projections, and the consistent initial points are identified to correct the errors.

Post-fault trajectory: Beyond the validity region

At the boundary of the validity region, T and z are modified according to (32) and (34). Because they are still close to the true ones at the boundary, the consistent initial points are evaluated at the boundaries. With the updated values for T and z, one can compute the coefficient to construct the analytical solution under the constraints that O1 and O2 are small with the updated values. Fig 10 show the post-fault trajectories beyond the validity region of (A) router, (B) the magnitudes of the bus voltages, (C) O_i(t), and (D) the estimated magnitudes of the internal voltages, respectively.

Stability of the system. Fig 10(A) shows an improved similarity between two models. Fig 10(B) illustrates the post-fault trajectories of the terminal voltages that indicate non-negligible changes of voltage magnitudes. The left and the right eigenvectors are identified corresponding to each eigenvalue λ of T to compute the condition of the eigenvalue, s(λ) from (25). The largest deviation of the eigenvalue corresponds to the smallest condition of the eigenvalue: maxj|λj−λ^j|≈ε[minjs(λj)]−1. It was found that the largest change in the eigenvalue was λmaxRe=0.3879 before the update, and the corresponding updated eigenvalue was λ^maxRe=0.3734, and s(λ) was approximately 1.8. It turns out that the eigenvalue has the largest real part in the positive, but the corresponding coefficients are not large enough to make the system divulge before t = 10 s (Type III stable). As shown Fig 10(A) and 10(B), both the analytical approach and TDS expect the stability of the system.

Validity region. As discussed with Fig 10(D), the validity region boundaries are frequently revisited (15 times), and consistent initializations are performed by the projections using (34). Note that the codes are not currently optimized to utilize the sparse structure of the matrices and to explain the partial update of the Jacobian matrices M^wandM^u in (34), which makes the computation process inefficient. The voltage magnitudes of Generator 2 swings widely because the generator is closest to the location of the line fault among all the generators.

Future works

In this work, for the sake of simplicity, the dynamics of the generators are represented using the classical model. It is important to model the generators correctly to discuss the power system dynamics [23]. We plan to accommodate generator models that have been widely used for decades in numerous commercial simulation programs for modeling round-rotor and salient-pole synchronous generators [26] (the GENROU/GENSAL models) and for the detailed treatment of saturation (the GENTPF/GENTPJ models) [45]. The immediate adjustment in modeling is on Ibus where the voltage magnitude E_i does not stay constant, which requires a modification in O1.

Conclusions

In this paper, we derived an analytical solution to the swing equations to assess the transient stability of power grids. To the best of our knowledge, our solution is a unique analytical solution to the swing equations without physically unacceptable assumptions, while obeying Kirchhoff’s laws. The solution indicates the factors affecting the stability after a disturbance occurs. Based on the solution, a new stability assessment approach is proposed. The assessment tool is different from the conventional assessment tools (COM, DM, and TDS approach) in that the derivation of our solution does not require the unphysical assumptions required for both COM and DM (such as lossless grids, constant voltages at all buses, and no consideration of reactive power). Moreover, its computational complexity is manageable. In addition to the low computational complexity, the approach proposed in this study explores the components affecting power system dynamics by examining the structure of the T matrix in terms of system dependent (T_sys) and operation point dependent (T_op) submatrices. The simulation results show that O1 and O2 are small in most cases. However, even in a case when O1 and O2 are large, it is possible to maintain O1 and O2 small by introducing the validity region based on the projection methods. The consistent initializations make it possible to identify the trajectories reliably.

Appendix

Time domain simulation, TDS

For comparison, TDS is performed using the Runge–Kutta method [46]. At first, the loads are converted to equivalent admittance to construct a linear model (i.e., the current is linear with the internal voltages). From the linear current and voltage relationship, the parameters in (2) are evaluated at the post-fault state. Because the magnitudes of internal voltages are assumed constant, the voltage at the Ibuses can be evaluated in terms of δ. Once the voltages at the Ibuses are computed, one can compute all the terminal voltages and real power injection at the Ibuses. Using the terminal voltages and the real power injection, it is possible to update δ.

We first define z=[δTωTw˜T]T, and rewrite the swing equation: y=dzdt=Az+b where A=[0I00−diag(D/M)000−L˜],b=(0f(z)−l˜), and f(z) = diag(1/M)[p_mech−p_elec(z)]. At the first sub step at iteration, with a step size Δt > 0, one finds

The update of z allows finding a power flow solution corresponding to the value of z for updating f(zm{1}). (A1) is generalized at the k^th sub step for the m^th iteration as follows:

The update of bm−1{k} follows (A2). Then the last update at the m^th iteration is

In averaging the four increments, a larger weight is given to the middle increments. Multiple simulations are performed with various time steps Δt in a decreasing order until the simulation results are invariant.

Direct method based on an energy function

The Lyapunov stability theorem gives sufficient conditions to determine the stability of the system. Numerical computation of the underlying ordinary differential equations is not necessary to derive the stability properties. The Lyapunov stability theorem gives two conditions for stability for a Lyapunov function V on an open set U. (1) V(x)={=0ifx=0>0otherwise and (2) dV(x)dt=∑j∂V∂xjdxjdt=∇V⋅dxdt≤0∀x≠0.

For a power system with multiple machines, the traditional DM utilizes the energy function under the assumptions of losses (γ_ik = 0) and of the fixed voltage magnitudes (E_k is fixed). The simplest energy function would be the mechanical interpretation [8] as follows:

The rotor angle δ^ref is the measure of the internal voltage angle with respect to the terminal voltage angle of the reference bus. For convenience, the reference frame for the rotor angle is redefined with respect to the rotating COI, i.e., δi=δiref−δCOI,ωi=ωiref−ωCOI where δCOI=1MT∑iMiδi, ωCOI=1MT∑iMiωi,MT=∑iMi. Equilibrium points are found at a set of the rotor angles that satisfies dδ_i/dt = 0 and dω_i/dt = 0, which leads to δik=sin−1(pimech|y˜ik|EiEk). Beside the stable equilibrium point, there will be the NI neighboring unstable equilibrium points to the stable equilibrium point. At the unstable equilibrium points, the “kinetic energy” (E_KE) is zero, and the energy equals the “potential energy”, E_PF. Along the unstable equilibrium points, one can find the closest unstable equilibrium point. At the point, the critical energy V_cr is defined so that any trajectory starting from a point with a lower energy than V_cr is guaranteed to converge into the stable equilibrium point if no other equilibrium points are contained in the set [10]. As a result, one can certify the system stability based on the stability margin, V_margin = V_cr−V_cl where Vcl is the current energy at the clearing time.