Analytical solution to swing equations in power grids

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.


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], The swing equation is a heterogeneous nonlinear second-order differential equation with multi-variables. There is no known method to solve the differential equation in an analytical fashion. Instead, several approaches to analyze the problem are suggested: 1) simplify the problem by ignoring the difficult components; 2) solve the problem for a "simple" system and extend the knowledge to a complex system; and 3) adopt a numerical approach. While these three approaches provide practical assessment for some cases, their applicability is limited due to their assumptions. 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.

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: , injection currents at Bus j when no voltages are applied; p i , real power injection from a generator i; p elec i!k , electrical power output from a generator i to Bus k; p mech i , mechanical power input to a generator i; p max , maximum power injection, p max jỹ ik jE i E k ; q i , reactive power injection from a generator I to Bus k; rad, radian; s(λ), condition of λ, sð�Þ ¼ ju H L u R j; t 0 , time when a disturbance (a physical anomaly and/or a control action) occurs; u L , left eigenvector of Mat, Mat ¼ �u H L ; u R , right eigenvector of Mat, Matu R = λu R ; v j , voltages at Bus j; v x , real part of the voltages v; v y , imaginary part of the voltages v; x I , real part of the loss-reflected voltage at Ibuses; y I , imaginary part of the loss-reflected voltage at Ibuses; jỹ ik j, magnitude of line admittance connecting a generator i and Bus k.
For the coupled oscillator model (COM), a system of particles aims to minimize the energy function Eðd K Þ ¼ 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 jỹ ik jE i E k sinðÀ y k þ d i À g ik Þ is rewritten as p max sinδ i , where p max ð¼ jỹ ik jE i E k Þ 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). ). In the pre-fault condition, the operation point is determined where the sine curve intersects with a horizontal line p mech 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 postfault 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][11][12][13][14][15][16][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 ϑ(kΔxk 2 ) [20], and that of the Runge-Kutta method is ϑ(kΔxk 4 ) [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].

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 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 righthand 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 À ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi part of admittance; the superscript 0+ represents "immediately after contingency"; and t > t 0 , 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) 2 , the square of the rate to deviate from the synchronization, is in general smaller than ðÀ q i þ b ii E 2 i Þ=M i while maintaining synchronization. According to O2, while maintaining synchronization, |O i (t)| is a slowly varying function that follows the Karamata representation theorem [25]: 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 . 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 E i E k sinðÀ y k þ d i À g ik Þ ¼ 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 x 2 i þ y 2 i ¼ E 2 i . 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 d j ¼ d 0 j þ m j ðdd j =dtÞ [8], where the first component d 0 j 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 d j ¼ m j ðdd j =dtÞ þ d 0 j [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 i l ¼ i cc l þ y ci l v l , 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 i L ¼ diagðy ci L Þv L þ i cc L . 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: Combined with the load modeling constant current and constant impendence, the current is expressed in terms of voltage 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., i KM Note that induction, frequency-dependent, and time-dependent loads are not included in the currentvoltage relationship. With no direct connection between Ibus and Mbus, Ohm's law over the network modeled in this work (Fig 2 (B)) is 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: The parameters i 0 Kbus and i 0 Mbus represent the constant current and constant impedance components, respectively. Note that at Kbuses and Mbuses, only loads (including zero loads) exists because the loads in this study excludes the induction, frequency-dependent, and time- dependent loads. Eq (1) yields the following relationship in terms of modified voltages at Ibus Using (2), the real and the imaginary components of the bus voltages are 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: The definitions of x i and y i lead to With (4) and (5), the swing equation in CCS becomes In comparison to the swing equation in the PCS, (6) is a heterogeneous nonlinear secondorder 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 Note that (8) involves no approximation. It is interesting that the reactive power injection appears in (8), while the swing equation considers only the real power balance. The constraints are conditions that the differential equations hold under. From O2, ½ðÀ Furthermore, the voltages at Kbus are expressed in terms of w I as shown in (2) Note thatR I associated withw I is an identity matrix. To distinguish the modified voltages derived by the load from a physical internal voltage,w I is introduced in (11). The buses where the frequency-or time-dependent loads are located belong to the Ibus. For a system with NI of Ibuses and NK of Kbuses, (10) is further simplified with the constraints (O1 and O2 are small) as follows: are the modified voltages derived from the frequency-or time-dependent load at Bus i. Let z be [w I ; dw I /dt;w I ], and then with the constraints (O1 and O2 are small) (12) and (13) becomes (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., d 0þ i . When a disturbance occurs, the mechanical power input does not change instantaneously, but the electric power output p elec i!k exhibits a sudden but finite change. The change is represented by a step function f ðx 2 R 1 Þ ¼ X i m i w A i ðxÞ where m i is a finite real number and χ A is an i . 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 dz G dt ¼ Tz G is as follows: where Ф is the coefficient matrix to be determined. The base solution u may contain terms associated with complex eigenvalues, because T is generally an asymmetric real-valued matrix. For any complex eigenvalue, its complex conjugate is also a complex eigenvalue. Let such a pair be l l þ l lþ1 ffi ffi ffi ffi ffi ffiffi À 1 p and l l À l lþ1 ffi ffi ffi ffi ffi ffiffi À 1 p . Then, the l th corresponding pair in the base solution is It is possible to represent u in terms of a real-valued function as follows: . .
where n and l are the indices of real and complex eigenvalues, respectively. Differentiating u with respect to time introduces a block-diagonal matrix D l such that du/dt = D l u, where , λ Re is the real eigenvalues, l Re Co and l Im Co are real and imaginary components of the complex eigenvalues of T, respectively. It is noteworthy that D l is obtained from the block Schur decomposition [31] T ¼ U À 1 T D l U T , so that D l has a blockdiagonal structure. Therefore, T and D l are related through the similarity transformation. In (15), Ф is the time-invariant coefficient matrix, and the differential equation yields Because u is not a null vector, (18) introduces a Sylvester equation; that is, FD l −TF = 0. The Bartels-Stewart algorithm is an efficient way to solve a Sylvester equation [31], and its computational cost is ϑ (N 3 ). It involves the orthogonal reduction of D l and T matrices into triangular form using the QR factorization, and then solving the resulting triangular system via back-substitution. Because D l is a block-diagonal matrix comprised of the eigenvalues of T, their eigenspaces overlap. Therefore, the Bartels-Stewart algorithm is not applicable to solve (18). Consideration of the null space yields vecðFÞ 2 null½ðI � T À D T l � IÞ T �, where � represents the Kronecker product. Therefore, the QR-factorization of I � T À D T l � I reveals the space of Ф. Let ψ l and Ѱ l be the l th vector spanning the null space of I � T À D T l � I and the matrix form of the vector, respectively. Then, Ф is a linear combination of all Ѱs.
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: Using o 0þ With a set of vectors ξ, ξ l = C l u(t = 0), β is determined from the following least square process: Analytical solution to swing equations The analytical solution to the swing equation that satisfies Kirchhoff's laws is If all the real parts of the eigenvalues are negative, u converges to zero as time increases. Therefore, under the circumstance and when O1 and O2 remain sufficiently small, 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) 3 ] [31], because the dimension of T is (4NI+2ND) × (4NI+2ND). Therefore, the first process is predominant if Nb > (NI+ND) 2 . 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 2 +y 2 −1 = 0. One may eliminate λ, and find a differential equation as follows: u 0 x þ u 2 þv 2 À y x 2 þy 2 ¼ 0; v 0 þ1 y þ u 2 þv 2 À y x 2 þy 2 ¼ 0 with the following constraints: x 2 +y 2 = 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 2 +y 2 = 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Ê i ¼ 1 À ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi  (10) is exact except for one element each row by an approxima- 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], Eq (25) indicates that the impact of the change in O i (t) is a linear change in the eigenvalues of T, the block-diagonal elements in D l . LetT andD l be the true matrices; δT and δD l be the difference matrices between the true and the approximated matrices (i.e.,T ¼ T þ dT andD l ¼ D l þ dD l ); andF and δΦ be the true solution matrix corresponding toT andD l and the difference matrix, respectively. Further, define ρ and σ: Eq (26) can be rewritten in terms of the Kronecker product, Using the Kronecker product, one finds the upper/lower bounds of the Frobenius norm of a product between two matrices A and B as follows: Eqs (27) and (28), and the triangle inequality theorem yield and using the triangular inequality and the fact that the Frobenius norm is no less than 2-norm [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., ε = kδTk F /kTk F �Δ T for all t�τ ε . Let y z be the first derivative of z with respect to time, then (14) yields: y z = Tz 0 +b where y z = dz/dt| t = 0 and z = z 0 +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 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: |ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl {zffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl ffl } 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]: 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, q c i in q c i þ b ii E 2 i : 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 b ii E 2 i term; therefore, its impact on the transient stability is highly limited. Unsettled mechanical power, p mech i in p mech i À g ii E 2 i : 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, p mech i 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, frequencyor 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 |� m2[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 kL † pk ε,1 �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 8 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 ε (= kδTk F /kTk F ), Δ T , is set to 1% for all the numerical evaluations. The value of kTk 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 kδTk F is 1.67 × 10 −2 (= 1% × 16.71). The     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 kTk 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 kδTk 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 kδTk F is 9.88 × 10 −2 (kδTk F /kTk F = 5.61 × 10 −3 ). 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. The trajectories of (A) δ, of (B) E, of (C) O i (t), and of (D) the estimated magnitudes of internal voltages due to the entire loss of the load at Bus 8 at t = 1 sec (red arrow in Fig 5). https://doi.org/10.1371/journal.pone.0225097.g007 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. kTk F is 7.88, and the threshold for kδTk F is 7.88 × 10 −2 (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 kTk F immediately after clearing the fault is 11.43; hence, the tolerance to the impact of O2 toward kδTk F is 0.114 (= 1% × 11.43). The impact of O2 is 1.79 (kδTk F /(kTk 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: max j jl j Àl j j � ε½min j sðl j Þ� À 1 . It was found that the largest change in the eigenvalue was l Re max ¼ 0:3879 before the update, and the corresponding updated eigenvalue waŝ l Re max ¼ 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 matricesM w andM 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 salientpole 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.
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 ¼ ½ d T o TwT � T , and rewrite the swing equation: where A ¼ 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.