Comparative analysis of existing models for power-grid synchronization

The dynamics of power-grid networks is becoming an increasingly active area of research within the physics and network science communities. The results from such studies are typically insightful and illustrative, but are often based on simplifying assumptions that can be either difficult to assess or not fully justified for realistic applications. Here we perform a comprehensive comparative analysis of three leading models recently used to study synchronization dynamics in power-grid networks -- a fundamental problem of practical significance given that frequency synchronization of all power generators in the same interconnection is a necessary condition for a power grid to operate. We show that each of these models can be derived from first principles within a common framework based on the classical model of a generator, thereby clarifying all assumptions involved. This framework allows us to view power grids as complex networks of coupled second-order phase oscillators with both forcing and damping terms. Using simple illustrative examples, test systems, and real power-grid datasets, we study the inherent frequencies of the oscillators as well as their coupling structure, comparing across the different models. We demonstrate, in particular, that if the network structure is not homogeneous, generators with identical parameters need to be modeled as non-identical oscillators in general. We also discuss an approach to estimate the required (dynamical) parameters that are unavailable in typical power-grid datasets, their use for computing the constants of each of the three models, and an open-source MATLAB toolbox that we provide for these computations.

The dynamics of power-grid networks is becoming an increasingly active area of research within the physics and network science communities. The results from such studies are typically insightful and illustrative, but are often based on simplifying assumptions that can be either difficult to assess or not fully justified for realistic applications. Here we perform a comprehensive comparative analysis of three leading models recently used to study synchronization dynamics in power-grid networks-a fundamental problem of practical significance given that frequency synchronization of all power generators in the same interconnection is a necessary condition for a power grid to operate. We show that each of these models can be derived from first principles within a common framework based on the classical model of a generator, thereby clarifying all assumptions involved. This framework allows us to view power grids as complex networks of coupled second-order phase oscillators with both forcing and damping terms. Using simple illustrative examples, test systems, and real power-grid datasets, we study the inherent frequencies of the oscillators as well as their coupling structure, comparing across the different models. We demonstrate, in particular, that if the network structure is not homogeneous, generators with identical parameters need to be modeled as non-identical oscillators in general. We also discuss an approach to estimate the required (dynamical) parameters that are unavailable in typical power-grid datasets, their use for computing the constants of each of the three models, and an open-source MATLAB toolbox that we provide for these computations (available at https://sourceforge.net/projects/pg-sync-models). a) Electronic mail: t-nishikawa@northwestern.edu

I. INTRODUCTION
In the field of network dynamics, power grids have long served as one of the main model systems motivating the study of synchronization phenomena [1][2][3] . As the field became more mature, an increasing number of researchers have been seeking to apply ideas from past theoretical studies to power grid-specific problems [4][5][6][7][8][9][10] . The need for such applications comes from the fact that, despite the extensive engineering literature on power systems, there remains a largely under-explored problem of how the large-scale network structure influences the collective dynamics in power grids.
While previous studies have emphasized the detailed modeling of relatively small test systems, the increasing availability of data processing tools, substantial computing power, and theoretical developments in network synchronization are now making it possible to address large-scale properties of power-grid systems.
A major concern for power grids is the stability of desired states, in particular synchronization stability of the power generators, which is a condition required for their normal operation. A frequency-synchronous state of n g power generators is characterized bẏ where δ i = δ i (t) is the angle of rotation associated with the i-th generator at time t. Studying the stability of synchronous states of an alternating current interconnection against perturbations requires a network model capable of describing the coupled dynamics of power generators. Different models have been used in different publications in the physics literature, and there has been no comprehensive comparison to clarify how these models are related to each other. Providing such a comparative analysis is the main focus of this article.
Here, we discuss three leading models, which we refer to as the effective network (EN) model 8,11 , the structure-preserving (SP) model 12 , and the synchronous motor (SM) model 13 . Each model is described as a network of coupled phase oscillators whose dynamics is governed by equations of where ω R is a reference angular frequency for the system, and H i and D i are inertia and damping constants characterizing the oscillators, respectively. The differences in the three models are reflected in the definition of the parameters A i , K ij , and γ ij , as well as in the number of phase oscillators. The phase angle δ i of oscillator i may represent either a generator or load. While all three Here we used the example case consisting of two generators (nodes 1 and 2) and one load node (node 3), which is discussed in Section VI A. (B) Representation of the electrical properties of the components in the same network. There are three possible ways to represent the load node, which are used in three different models. (C) Representation of the system as a network of coupled oscillators for each choice of load representation in (B). For the system parameter values given in Section VI A, the network dynamics obeys Eq. (2) [or equivalently, Eqs. (28), (29), and (30) for the EN, SP, and SM model, respectively] with the indicated values of A i , K ij , and γ ij . Each of the three dynamical models has its own definition of nodes that are different from that used in (A). In (B), these nodes are shown as black dots and indicated by orange, blue, and green indices for the EN, SP, and SM models, respectively. The same coloring scheme is used for the node indices in (C). Note that in the SP model the generator terminals are treated as load nodes with zero power consumption (nodes 3 and 4), separately from the generator internal nodes (nodes 1 and 2), leading to the 5-node representation. models represent the n g generators as oscillators, they are distinguished mainly by their different modeling approaches for the loads, which are representations of individual or aggregated consumers that draw power from specific points in the (transmission) network. The EN model represent the loads as constant impedances rather than oscillators, thus focusing on the synchronization of generators as second-order oscillators. The SP model represents all load nodes as first-order oscillators (i.e., H i = 0), and each generator is represented by two oscillators, including one for its terminal (a point connecting the generator to the rest of the network). The SM model assumes that the loads are synchronous motors that are represented as second-order oscillators. Figure 1 shows a simple example network that illustrates how these differences lead to different network dynamics representations with different number and values of model parameters. Note that the parameters are denoted in the figure using appropriate superscripts (EN, SP, and SM), which is a convention we will use throughout the article.
The parameter A i , along with D i , determines the i-th oscillator's inherent frequency, ω * i := , which is the equilibrium frequency of oscillator i in the absence of the coupling term [the summation term in Eq. (2)]. For generators, this frequency is typically much higher than the system reference frequency ω R , as the two terms on the r.h.s. of Eq.
(2) balance each other in a steady state. In a realistic setting, however, the instantaneous frequency of a generator is unlikely to actually reach this inherent frequency, since the system operator would take control actions well before the frequency deviates too far from the designated system frequency ω R . In addition, the equation of motion for generators and motors we derive below assumes that the frequency remains close to ω R and thus would no longer be valid if the frequency deviates too far from ω R .
Nevertheless, this definition of inherent frequency can be useful in characterizing the nature of individual oscillators and, in fact, is analogous to the one usually used in studying networks of coupled phase oscillators. Note that "natural frequency" is a term commonly used to refer to ω * i in that context, but it has a very different meaning in the power systems literature; see Appendix A.
The parameter K ij ≥ 0 represents the strength of dynamical coupling between oscillators i and j, while γ ij represents a phase shift involved in this coupling. We note that Eq. (2) can be regarded as a second-order analog of the Kuramoto model 14 with arbitrary coupling structure. Similar forms of second-order equations for coupled oscillators have been used to study synchronization and phase transitions outside the context of power grids [15][16][17][18] .
The paper is organized as follows. In Sec. II, we discuss how to characterize and compute the steady state of a power-grid network. We then turn our attention in Sec . III to the problem   of describing the dynamics of individual generators, deriving the basic equation of motion and   describing the electric circuit representation of a generator under the classical model. In Sec IV, we derive the EN, SP, and SM models and discuss their differences and similarities. In Section V, we describe how to obtain model parameters required for numerical simulations. Finally, we discuss in Sec. VI an instructive example system and a selection of test systems and real power-grid datasets, focusing particularly on the heterogeneity of the constants A i and the inherent frequencies as well as on the coupling structure encoded in the constants K ij . We make concluding remarks in Sec. VII.

II. STEADY STATES OF POWER-GRID NETWORKS: THE POWER FLOW EQUATIONS
A power grid is a system of electrically coupled devices whose purpose is to deliver power from generators to consumers. To view the system as a network, we define a node to be a point in the system at which power is injected by a generator or extracted by power consumers, or a branching point through which power is redistributed among the transmission lines connected to the point.
A link is then defined as an electrical connection between a pair of such nodes and can represent a transmission line or a transformer. This network representation of the physical structure of a power grid is illustrated in Fig. 1A  The standard approach 19 , which we adopt here, is to represent the parameters of these models for transmission lines and transformers in terms of equivalent admittances (the inverse of impedances).
The structure of the entire physical network can then be represented 19 by the (complex-valued) where an off-diagonal element Y 0ij is the negative of the admittance between nodes i and j = i and a diagonal element Y 0ii is the sum of all admittances connected to node i (including the shunt admittances to the ground, which are parts of the models for transmission lines and transformers).
The operating condition of a power grid can be characterized by the distribution of power flow through the physical network of transmission lines and transformers. In an alternating current grid, power is represented as a complex number, whose real and imaginary components, P i and Q i , are called active and reactive power, respectively. Alternatively, the power flow state can be characterized uniquely by the complex voltage V i (with magnitude |V i | and phase angle φ i , so that V i = |V i |e jφ i ) and the complex power P i + jQ i injected into each node i in the network. Given the locations and parameters of power generators and loads, as well as the parameters of the network summarized in the admittance matrix Y 0 , fundamental laws of physics-the Kirchhoff's lawsdetermine the power flow state of the system. The laws can be used to derive the so-called power flow equations 19 , which provide the foundation for steady-state power system analysis: where the complex admittances are represented in polar form as Y 0ij = |Y 0ij |e jα 0ij , we define γ 0ij := α 0ij − π/2, and n is the number of nodes. To solve this set of 2n nonlinear equations for all 4n quantities that determine the power flow state of the system, we require that two real quantities per node are given as parameters. The usual assumptions are: • If node i is a generator node, the values of P i and |V i | are given as P i = P * g,i and This is reasonable because in real power grids the generators are scheduled to produce constant active power at a given time, and constant voltage magnitude is usually maintained by voltage regulators. Usually, one generator node is chosen to be a reference node, for which φ i is set to zero instead of specifying the value of P i . For the steady-state analysis, we need at least one node with unspecified P i because the total power generation cannot be known a priori even when the total power consumption is known, since the difference, which equals the power lost in the transmission lines and transformers, depends on the power flow state.
• If node i is a load node, the values of P i and Q i are given as For positive active and reactive power consumptions, P i and Q i take negative values. These values can be determined from historical data, load forecast, or measurement 11,19 .
Given the parameter values and the admittance matrix Y 0 , the nonlinear equations (3) and (4)  Most generators in today's power grids have a rotating mass (a rotor) that is driven by mechanical torque to generate electrical power. There exist models with various levels of sophistication and complexity, but in all such models, the rotational dynamics of the rotor is ultimately governed by a fundamental law of physics: the Newton's second law. The electrical output from such a generator is coupled to other generators in the grid through a network of transmission lines, transformers, and other devices, which serve to transport the generated power to consumers. The power demanded by the network is felt by the generator's rotor as an electrical torque, which is usually a decelerating torque. It is the balance between the mechanical input power and the electrical output power that determines the dynamics of the rotor.
To derive the equation of motion governing the dynamics of such a generator, we set the rate of change of the angular momentum of the rotor equal to the net torque acting on the rotor: where J is the moment of inertia in kg·m 2 , δ is the angle of the rotor relative to a frame rotating at the reference frequency ω R in rad/s,T m is the mechanical torque in N·m driving the rotor, D m is the damping coefficient in N·m·s for mechanical friction, ω is the angular frequency of the rotor in rad/s, R is the regulation parameter in rad/N·m·s characterizing the proportional frequency control by a governor, ∆ω := ω − ω R is the frequency deviation, D e is the damping coefficient in N·m·s for the electrical effect of the generator's damper windings, and T e is the typically decelerating torque in N·m due to electrical load in the network. Noting thatδ = ω − ω R = ∆ω, we can rewrite the equation as whereD := D m + D e + 1/R and T m :=T m − D m ω R is the net mechanical torque, accounting for frictional loss at the reference frequency. Multiplying both sides by ω and using the fact that torque in N·m multiplied by angular velocity in radians per second gives power in watts, the equation can be written in terms of power: where we defineP m := T m ω andP e := T e ω, and we assumed the factor ω R /ω to be nearly equal to one, i.e., that the generator's frequency ω is close to the reference frequency ω R .
This approximation, which can be formalized using singular perturbation analysis 22 , is valid and considered appropriate for studying power system stability, where the key question is whether the system desynchronizes after a perturbation from a steady-state operation near the reference frequency. If, for example, a disturbance leads to even 1% deviation of a generator's frequency from the 60 Hz reference frequency for a period of just one second, it would lead to an increase in angle difference equivalent to more than half a rotation. This much of angle deviation is highly likely to cause loss of synchronization in practical situations 11 . While the approximation ω ≈ ω R is widely accepted in this context, a more detailed electromechanical model that does not require this approximation is also available 23 .
We now divide both sides of the equation by the rated power P R (used as a reference) to makẽ P m andP e per unit (p.u.) quantities, which is a normalization procedure commonly used in power systems studies. The factor Jω R then becomes 2H/ω R , where we defined the inertia constant H := W/P R (in seconds) and the kinetic energy of the rotor W := 1 2 Jω 2 R (in joules). The factorDω R becomes D/ω R , where we defined the combined damping coefficient as D :=Dω 2 R /P R (in radians).
We then obtain what is known in the power systems literature as the swing equation 11,19,22 : which we use here as the fundamental equation of motion for a generator. The term P m represents the net mechanical power input to the rotor, while P e represents the electrical power demanded by the rest of the network and includes terms that depend explicitly on δ and the state variables of the other generators and loads in the network.

B. Electric circuit representation -The classical model
In general, both P m and P e have nonlinear dynamics, which may need to be accounted for in high-accuracy simulations that power system engineers require. There has been substantial effort in the power systems literature to model the dynamical behavior of the generator's internal magnetic flux, which affects the electrical power output P e 11,19,23-26 . One may also need to include the nonlinear dynamics of the governor that controls the generator frequency and the excitation system that controls the voltage magnitude at the generator terminal. Here we aim to study power grids as a complex dynamical system and focus on the way in which the network structure influences the synchronization dynamics of generators in simplest settings that are nevertheless realistic. For this purpose, we now describe a model of a generator used commonly in the engineering literature, particularly in theoretical studies, and forms a basis of all three network models we present below.
In this so-called classical model, a generator is represented as a voltage source with constant voltage magnitude |E| connected to the terminal node through a reactance x d > 0 (see Fig. 1B). The phase angle of the voltage source is assumed to be equal to the rotor's rotational angle and thus denoted by δ. In addition, the mechanical power input P m to the rotor in Eq. (8) is assumed to be constant.
Constant voltage magnitude and mechanical power are approximations that are valid for short-term dynamics. For analysis of transient stability after a disturbance, the approximation is considered valid for simulation of the "first swing" of the resulting oscillation of ω or δ, which typically covers a time period on the order of one second 11 . The classical model sometimes refers to the one in which damping is ignored [i.e., D = 0 in Eq. (8)], but here we consider damping effects explicitly, as they can have non-negligible effect on the stability of steady-state power grid operation and can potentially be used as tunable parameters for optimizing the stability 8 .
In a real power grid under a steady-state condition, the system frequency is continuously monitored and is in fact actively controlled (by adjusting P m of its generators) to stay close to the reference frequency of the system. Because of this, in power system stability analysis, it is usually assumed that all generators are initially synchronized at the reference frequency ω R in a steady state. Then, for a given steady state in which the generator is delivering active power P * g through its terminal, we haveδ =δ = 0 and P m = P e = P * g . This value of P m is then held constant when studying stability using the swing equation and the classical model of the generator.
Regardless of whether the grid is under steady-state condition or not, the expression for the electrical output power is given by the so-called power-angle equation, While H, D, and x d are constants that characterize physical and electrical properties of the specific generator, P * g and |E * | are parameters that depend on the steady-state distribution of power flow across the network. By eliminating the current I from the expression for the complex power, where the bar denotes the complex conjugate), and the Ohm's law, jx d I = E − V , applied to the transient reactance, the steady-state internal voltage magnitude |E * | can be calculated as where Q * g is the reactive power injected by the generator into the terminal and |V * | is the terminal voltage magnitude in the steady state. We thus see that the state of the terminal node given by P * g , Q * g , and |V * | determines some of the parameters required to simulate generator dynamics using Eq. (9). This is a crucial point that, as we will show below, has significant implications for the modeling of a power grid as a network of coupled oscillators. In a typical stability analysis using the classical model, P * g and |V * | are given as input data, and Q * g is obtained by solving the power flow equations (3) and (4) for the entire network, given necessary input data at other nodes. The power flow solution also provides values of the phase angles, δ * and φ * , in the steady state, which can be used as the initial condition for Eq. (9). Note, however, that the equation is not closed, as the dynamics of φ = φ(t) and V = V (t) are not yet specified. It is mainly the phase angle φ that serves as the medium through which this generator is coupled to the rest of the system.

IV. COLLECTIVE DYNAMICS OF GENERATORS IN POWER-GRID NETWORKS
The dynamical state of the network is characterized by four functions of time for each node: P i , The power flow equations (3) and (4), being valid at each instant of time for timevarying variables, provide two out of the four required equations per node for uniquely determining the dynamical state of the system. We thus need two additional equations/assumptions at each node. For a generator node, the swing equation effectively provides one equation, and the usual practice is to make the additional assumption that the reactive power Q g injected by the generator at its terminal node is constant over short time scales and equals its steady-state value Q * g , which is computed from the (static) power flow equations. We would then have a complete set of equations that allows us to determine all four state variables of the generator, |E|, δ, P g , and Q g , as a function of time, given the initial condition and the state variables (for all t) at all the other nodes. Note that E = |E|e jδ is related to the terminal voltage V = |V |e jφ through the transient reactance at all times. Of the four generator variables, δ directly relates to Eq. (1) and thus is the most relevant for the problem of synchronization stability.
Supplying the full set of equations for load nodes requires modeling of static and/or dynamic behavior of the loads. Load modeling is a hard problem because the power consumption at a node in a transmission network is typically an aggregate of a large number of loads of a wide variety of types and sizes, each of which can be time dependent and influenced by human activity. As a consequence, a number of different models have been proposed and used in the power system literature. As mentioned above, it is desirable to have a simple but realistic model that can clarify the role of network structure in the problem of power-grid synchronization stability. We thus provide in this article a comparative analysis of three approaches that can be used to recast the problem as that of synchronization in complex networks of coupled oscillators. These approaches lead to the three different models, the EN, SP, and SM models, which are all described by Eq. (2) but with different values and interpretations for the parameters A i , K ij , and γ ij .

A. The effective network model
The dynamical interaction between a pair of generators is mediated by the paths of transmission lines connecting the generators and is expressed as the sinusoidal dependence on the terminal voltage phase in Eq. (9). While the nature of the coupling between generators is ultimately shaped by the structure of the transmission network and the distribution of loads, it would be insightful and useful if one could represent the coupling in a single term that depends on the state variables of the generators. This can indeed be achieved by modeling loads as constant impedances and reducing the physical network to what we call the effective network of interactions between generators. Since the model is derived through this reduction process, it is also called the network-reduction model or network-reduced model in the literature.
In order to see how the reduction process leads to the effective network, which then determines the coupling constants K ij governing the network dynamics in Eq. (2), let n g and n denote the number of generator terminal nodes and load nodes, respectively, so that n = n g +n . In addition to these n nodes, we define an additional node to be a point between the internal transient reactance and the constant voltage source for each generator (see the two generators depicted in Fig. 1B), making the total number of nodes N := 2n g +n = n g +n. By suitable re-indexing, we may assume that i = 1, . . . , n g correspond to the generator internal nodes, i = n g + 1, . . . , 2n g to generator terminal nodes, and i = 2n g + 1, . . . , N to load nodes. Define Y d as the n g × n g diagonal matrix whose diagonal elements are the admittances (jx d,1 ) −1 , . . . , (jx d,ng ) −1 of the generator transient reactances. We write the admittance matrix Y 0 for the physical network as where Y gg 0 , Y g 0 , Y g 0 , and Y 0 are the four blocks resulting from separating the first n g rows (columns) from the last n rows (columns).
The active power P * ,i and reactive power Q * ,i consumed at load node i (and possibly also at a generator terminal node) in a steady state is represented by an equivalent impedance to the ground 11 whose admittance is where |V * i | is the voltage magnitude at the node. The constants P * ,i and Q * ,i are part of the input data for Eqs. (3) and (4) to solve for a power flow solution, which provides |V * i |. The constant admittance Y ,i is computed from these steady-state values and added to the corresponding diagonal components of Y gg 0 and Y 0 to obtain Ỹ gg 0 andỸ 0 . This representation requires the assumption that the power demand is constant. We thus consider the dynamics of the system on time scales that are short enough for the validity of the assumption, but long enough to address the problem of synchronization. We can now define an ) that includes the links representing the transient reactances connecting the generators' internal and terminal nodes, as well as the equivalent impedances for the loads: where 0 denotes matrices (of appropriate sizes) whose elements are all zeros.
Now let V g , V t , and V be the vectors of node voltages for generator internal nodes, generator terminal nodes, and load nodes, respectively, and stack them vertically in that order to form the voltage vector V. Also let I g be the vector of currents injected into the system at the generator internal nodes. Because the loads are modeled as constant impedances to the ground, all nodes have zero injection currents except for the generator internal nodes. The Kirchhoff's current law can then be written in the form I = Y 0 V, or equivalently The system is then converted to I g = Y V g by eliminating V t and V , a procedure known as Kron where 1 denotes the n g × n g identity matrix. The inverse of Y d always exists, while the matrix Using the effective admittance matrix Y EN , the single sinusoidal coupling term in Eq. (9) can be replaced by an expression for P e that comes from a power balance equation equivalent to Eq. (3) with |V i | replaced by |E * i |, Y 0 by Y EN , and φ i by δ i : where Y EN ij = |Y EN ij | exp(jα EN ij ) and all quantities must be expressed in p.u. with respect to the system base. This can be used to show that writing the swing equation (8) for each generator leads to an equation of the same form as Eq. (2), and the EN model is thus given by where G EN ii is the real part of the complex admittance Y The equation for each generator is also of the same form but with H i > 0, and is in fact exactly Eq. (9) because the only link from the generator's internal node is to its terminal node through a purely imaginary admittance, Y SP i,i+ng = −(jx d,i ) −1 , i = 1, . . . , n g . Putting together, the SP model is given by for generator internal nodes and for load nodes (including the generator terminal nodes), where G SP ii denotes the real part of Y SP ii .
Note that these equations have the same form as Eq. (2), and in this case the number of phase oscillators is N = n g + n, since each of the generator internal and terminal nodes, as well as load nodes, is represented by a single variable δ i . In contrast to the EN model, the SP model has the advantage that the physical structure of the transmission network represented by Y 0 is preserved as part of the admittance matrix Y SP and reflected in the coupling constants K SP ij . This, however, comes at the expense of additional complexity associated with the larger set of equations (real power grids can have many more nodes than the number of generators) and increased uncertainty associated with the additional parameters D i that must be estimated. Note that in the SP model and mentioned in Sec. IV A for the EN model (also assuming lossless transmission lines) has been shown in the same publication to be applicable also to the SP model. Variations of the SP model that represent droop inverters as first-order oscillators 32 and incorporate stochastic deviations of frequencies 33 have been used to study synchronization stability in networks with renewable energy sources.

C. The synchronous motor model
Another way to model the dynamics of the loads while preserving the physical network structure is to use synchronous motors to represent the loads. A synchronous motor is essentially the same type of machine as a generator, except that the flow of power is in the opposite direction, converting electrical power into mechanical power, and thus can be modeled by the same swing equation, Eq. (8) with P m < 0 and P e < 0. In addition to the n g internal nodes we defined for the generators in the EN model, in this case one defines n internal nodes for the synchronous motors representing the loads, which makes the total number of nodes 2n. Since the motors are modeled in exactly the same way as the generators, this network model is mathematically equivalent to the EN model in which all n nodes of the physical network has "generators," n of which have negative mechanical power P m,i < 0. The matrix Y 0 will be replaced by in this case, where Y d is the n × n diagonal matrix whose diagonal elements are the admittances

V. OBTAINING MODEL PARAMETERS REQUIRED FOR SIMULATIONS
The parameters required to run a simulation of the synchronization dynamics using Eq. (2) are the system reference frequency ω R , the parameters of individual generators/motors, H i , D i , and x d,i , as well as the constants for the network dynamics equation, A i , K ij , and γ ij . The parameters ω R , H i , D i , and x d,i must be given, while A i , K ij , and γ ij are computed by solving the power flow equations (3) and (4) and performing Kron reduction. As described in Section II, a power flow solution is usually calculated given the parameters P * g,i and |V * i | at each non-reference generator, P * ,i and Q * ,i at each load node, and |V * i | at the reference generator, as well as the admittance matrix Y 0 representing the physical network structure. The type of power system dataset most common in engineering is suitable for this calculation and thus contains the required parameters at each node, as well as sufficient information to compute Y 0 . Given such a dataset, standard power systems software, such as PST 20 11,19 . This observation, combined with the assumption that the rated power (which is unavailable in many datasets) is proportional to P * g,i , implies that x d,i is inversely proportional to P * g,i when x d,i is expressed in p.u. with respect to a reference power common to the entire network. In our previous work 8 we found that the available data for the first three systems in Table 1 follows our formula more closely than the standard approach. For D i , one may use

VI. HETEROGENEITY OF MODEL PARAMETERS
Given the input data for the system, the constants A i , K ij , and γ ij for the network dynamics governed by Eq. (2) can generally be distributed heterogeneously across the network. Since A i are related to the inherent frequencies ω * i through the same formula ω * i = ω R (1 + A i /D i ) for all three models, the A i -heterogeneity quantifies the heterogeneity of the intrinsic properties of the generator dynamics. While the heterogeneity in the generator parameters is naturally a source of heterogeneity in A i , it is actually possible to have heterogeneous A i even in the complete absence of heterogeneity in the generator parameters H i , D i , x d,i , P * g,i , and |V * i |. The first two parameters actually do not enter into the calculation of A i (as described in the previous section) and thus do not affect the A i -heterogeneity. The other three affect the values of A i , but do not fully determine A i , since the calculation also involves the admittance matrix Y 0 as well as P * ,i and Q * ,i for the loads. In the following subsections, we first show that A i -heterogeneity can indeed arise without heterogeneity of generator parameters using a small network example, and generalize it to larger networks. We then discuss the A i -heterogeneity as well as the coupling matrices K ij and γ ij for examples of realistic systems.

A. System with two generators and one load
We consider the simplest case with which the heterogeneity of generator parameters can be discussed: a system with two generators, each connected to a single load by a transmission line, as shown in Fig. 1. For this system we have n g = 2, n = 1, n = n g + n = 3, and N = 2n g + n = 5.
Assuming lossless and inductive transmission lines (i.e., r 13 = r 23 = 0, x 13 > 0, and x 23 > 0), the power flow equations for this network read where B 011 = Im(Y 011 ) = b 13 2 − 1 x 13 > 0, B 022 = b 23 2 − 1 x 23 > 0, and B 033 = B 011 + B 022 . Setting φ 1 = 0 and using it as the reference angle, we solve this set of five equations for five unknowns, φ 2 , where SP model: where SM model: where Note that the phase shifts γ SP ij and γ SM ij do not appear in the equations, since the assumption of lossless and inductive transmission lines implies that they are both zero, as discussed in Sec. IV B and Sec. IV C.
Assuming that P * g,1 , P * g,2 , |V * 1 |, and |V * 2 | are all nonzero, Eqs. (22) and (23) imply  The effect of network asymmetry demonstrated here for this small example system when using the EN and SM models suggests generalization to larger systems, as well as to other forms of network asymmetry, such as heterogeneities in the distribution of power demand, generator locations, and the connectivity of generators/loads in terms of the number of transmission lines. This will be discussed in the next section.

B. Larger power-grid networks
Consider a system with identical generator parameters: P * g,i = P * g , |V * i | = |V * |, and and 38, respectively. The data for the Guatemala and Northern Italy systems were provided by F.
Milano (University of Castilla -La Mancha), and the data for the Poland system were provided as part of the MATPOWER software package 21 . The required dynamic data that are not available from the respective sources cited above were computed as described in Section V. Figure 3 shows the computed inherent frequencies ω * i = ω R (1 + A i /D i ) under the three models for each of the systems in Table I. Note that those realistic systems have heterogeneous generator parameters, which is a source of heterogeneity in A i in addition to the effect described above that results from the heterogeneity and asymmetry of the physical network structure. We see that the inherent frequencies can be significantly different from the system reference frequency ω R , with mostly higher frequencies for the generators and lower frequencies for the load nodes (except for the EN model, in which the load nodes are eliminated and thus have no well-defined inherent frequencies). This is consistent with the general picture that the voltage phase angles of the generators tend to rotate at higher frequencies and pull the frequency of the rest of the grid upward, while the load phase angles tend to rotate at lower frequencies and pull the rest of the grid downward. The variation of the inherent frequencies among generators and load nodes can be extremely large. Indeed, we find a generator with ω * i > 7 × 10 2 Hz for the 50-generator system, and a load node with ω * i < −1.1 × 10 3 Hz (corresponding to a rotation in the opposite direction with frequency > 1.1 × 10 3 Hz) for the Northern Italy system. As noted in Sec. I, frequencies far from ω R are not likely to be observed in reality, but these values of inherent frequency do reflect the inherent dynamical property of the oscillators representing generators or loads. The observed large variation in inherent frequency across the network may indicate a high level of frustration in the system in the sense that oscillators with vastly different frequencies are forced to rotate at  Table I  grid. Figure 4 shows a visualization of the coupling matrix K = (K ij ) and the phase shift matrix Γ = (γ ij ) for the 50-generator system. Recall first that the size of the coupling and phase shift matrices is different for different models: n g ×n g for the EN model, N ×N for the SP model (where N = 2n g +n ), and n×n for the SM model (where n = n g +n ). Despite the fundamental differences in the dimensionality and in the type of equations used in the three models, we can identify some similarity within the coupling structure. We see in Fig. 4A that there is a non-trivial structure in the dynamical coupling among the generator internal nodes under the EN model. Although this coupling matrix is fully-populated, the strength of coupling varies widely and spans across many orders of magnitude, from 10 −9 to order one. Notice that a similar coupling structure can be identified in the middle box in Fig. 4B, corresponding to the sparse physical network connecting the generator terminal nodes in the SP model. We can also identify a similar structure in the top left box in Fig. 4C for the SM model. Most of the γ ij values are close to zero (red in Fig. 4D-F), which indicates that the coupling between oscillators i and j tends to make the oscillators synchronized with a small phase angle difference. A zero phase shift corresponds to an admittance (physical or effective, depending on the model) that has positive imaginary part and zero real part, which means that it represents a pure capacitive reactance. There are, however, some pairs of oscillators for which γ ij is significantly different from zero (green to blue in Fig. 4D-F). Such large phase shifts correspond to resistive or inductive components.

VII. CONCLUSIONS
We have presented first-principle derivations for three different network-level models of powergrid synchronization dynamics. These models are all based on the same fundamental equation of motion for mechanical rotation (the swing equation) and the same electric circuit representation (the classical model) of power generators and synchronous motors. The models are, however, clearly distinguished by the different sets of assumptions adopted for the loads. All three models are formulated in their full generality, which helps clarify how different choices of assumptions lead to different special cases of these models that have been used in previous studies. We have shown that even when we assume all the generators to be completely identical, the parameters of the phase oscillators representing them in the EN and SM models may be quite different. We have also shown that the oscillator parameters are vastly heterogeneous in a selection of test systems and real power-grid datasets.
While the validity and appropriateness of the models depend on the purpose, the choice of additional assumptions, and the system under consideration, these three models lay a solid foundation for the study of synchronization dynamics in complex power-grid networks. If the purpose is to isolate the effect of the network topology on synchronization, making simplifying assumptions (as has been done in many previous studies) to homogenize all other aspects of the system may be appropriate. If the purpose is to investigate the interplay between heterogeneities in the parameters of the generators, loads, and transmission lines, it may be helpful to use a simple, random, and/or regular network topology. If the purpose is to understand the network structure and the component parameters of realistic power systems, then keeping the details in a given dataset and comparing the synchronization properties with null models will provide the most useful insights.
The models we considered here can serve as bases for stability analysis against small or large perturbations, for which a variety of analytical and numerical techniques are available. As shown in Ref. 8, stability against small perturbations can be fully addressed using the powerful conceptual framework of master stability functions 29 based on linear stability analysis and network spectral theory. In contrast, results on stability against large perturbations tend to be more limited in scope.
If the relevant admittance matrix has zero real components, the problem can be tackled by energy function-based approaches 24 , which can also be used in combination with bi-stable representation of transmission lines for the modeling of the dynamics of cascading failure propagation in power grids 39 . Other approaches include estimating the region of attraction of a synchronous state 40 , quantifying the size of the region using a recently proposed measure called the basin stability 41 , and analyzing nonlinear modes using Koopman operators 42,43 . Such analyses of synchronization stability against small and large perturbations may help develop strategies for improving the existing power grids by elucidating the relation between network structure and stability, which we have recently demonstrated to be highly non-intuitive in general models of coupled oscillator networks 44 .
Understanding the dynamics of the models we considered here can potentially serve as a stepping stone for studying other types of instabilities relevant for power grids. For example, studying the so-called dynamic bifurcation problem 45 for these models, in which slow variation of the parameters A i , K ij , and/or γ ij drives the system toward a bifurcation point, can provide insights into the issue of voltage instability in power grids 46 . To capture the fast dynamics that follows the bifurcation, one may supplement the models we studied with additional equations that describe time-varying voltage magnitudes [47][48][49] . Such model-based studies of voltage instability complement model-free approaches for predicting the bifurcation point based on time series measurement 50 . Other causes for instability include fluctuations in power demand and failure of system components, and modeling them as stochastic processes can help optimize response strategies for power grid operators 51 .
With the anticipated transition to smart grids, there will be additional fluctuations in the system from new components such as intermittent energy sources and increased use of plug-in electric vehicles 52 . Such fluctuations will have impact on the synchronization stability of power grids, and its implications will be difficult to understand unless we understand the short-term dynamics of the models considered here. By providing comprehensive comparative analysis of these models at the intersection of the physics of complex systems, network science, and power system engineering, we hope to contribute to a better understanding and control of the dynamics of existing power grids as well as to a better design of future power grids. Appendix: Single generator connected to a large power grid

ACKNOWLEDGMENTS
Consider a generator connected through a single transmission line to a large power grid. Since the power grid is large, the dynamics of the generator does not affect the behavior of the rest of the grid much, and we thus approximate the voltage magnitude and the phase angular frequency at the other end of the line to be constant. This approximation corresponds to having the generator connected to a sort of "heat bath," which can be modeled as an artificial generator that has infinite inertia and can absorb an arbitrary amount of active and reactive power without changing its terminal voltage magnitude and phase frequency (and thus is often called an infinite bus, where the term "bus" is used by power system engineers to refer to a network node). We measure the generator's phase angle δ relative to the phase of the voltage at the power-receiving end of the transmission line, which has magnitude |V ∞ | and phase rotating at the reference frequency ω R . Combining the transient reactance and the transmission line impedance Z = R + jX to obtain a total admittance of Y = [R + j(X + x d )] −1 = |Y |e jα = G + jB (which defines α), the active power output P e from the generator at its internal node can be calculated as P e = |E * | 2 G + |E * V ∞ Y | sin(δ − γ), where γ := α − π/2. Equation (8) then becomes where A = P m − |E * | 2 G and K = |E * V ∞ Y |. Upon changing the variable to θ := δ − γ and defining M := 2H/ω R andD := D/ω R , the equation becomes Mθ +Dθ = A−K sin θ, which is the equation of motion for a forced damped pendulum, where θ is the angle the pendulum arm makes with the vertical, M is the mass,D is the damping coefficient, and A is the forcing torque. We thus see that a pair of stable and unstable equilibria exist if and only if A < K, and they are characterized byθ * = 0 and A = K sin θ * and correspond to the states of the generator synchronized to that of the rest of the grid at the frequency ω R .
IfD 2 > 4M √ K 2 − A 2 , the effect of damping is strong enough, andδ will converge exponentially to zero (corresponding to the generator's instantaneous frequency converging to ω R ) after a small perturbation is applied when the system is at the stable equilibrium. Otherwise, the damping is too weak, and the frequency will oscillate around ω R with decaying amplitude with a characteristic frequency determined by the parameters. WhenD = 0, the equilibrium is neutrally stable and the (angular) frequency of oscillation around the equilibrium is (K 2 − A 2 ) 1/4 / √ M , which is referred to as the "natural frequency" of the generator in the power systems engineering literature 11 (3) and (4) (with n = 2) give where φ is the voltage phase angle at the terminal relative to that at the other end of the transmission line and Q * g is the reactive power injected into the line. This can be used to express Q * g in terms of the given parameters as Q * g = |V * V ∞ | 2 /|Z| 2 − (P * g ) 2 . This, together with P m = P * g and Eq. (10), determines the dependence of the parameters of Eq. (A.1) on the steady-state power flow of the system, specifically as a function of P * g , |V * |, and |V ∞ |.