Normal Form for Grid-Forming Power Grid Actors

Future power grids will be operating a large number of heterogeneous dynamical actors. Many of these will contribute to the fundamental dynamical stability of the system, and play a central role in establishing the self-organized synchronous state that underlies energy transport through the grid. We derive a normal form for grid-forming components in power grids. This allows analyzing the grids’ systemic properties in a technology neutral manner, without detailed component models. Our approach is based on the physics of the power ﬂow in the grid on the one hand, and on the common symmetry that is inherited from the control objectives grid-forming power grid components are trying to achieve. We provide an initial experimental validation that this normal form can capture the behavior of complex grid-forming inverters without any knowledge of the underlying technology, and show that it can be used to make technology-independent statements on the stability of future grids.


I. INTRODUCTION
The transport of energy through the power grid depends on a self-organized synchronous state of distributed dynamical actors at continental scale.These dynamical actors are called grid forming.They establish the dynamical state on which power flow is possible.In the current grid these are primarily conventional power plants, with control schemes designed around heavy rotating masses.The energy transition demands a shift from such conventional generators towards inverter-interfaced renewable energy sources, with their dynamics specified by power electronics.This poses a fundamental challenge to understand the collective phenomena of these novel grid-forming components, and ensure the existence and resilience of the self-organized synchronous state.Consequently, the collective dynamics of power grids has become a very active interdisciplinary area of research in recent years.
The design of conventional generators is well established and grounded in the physics of synchronous machines.There is a good understanding in the engineering community what level of model detail is required to study which questions.This is not the case for inverterinterfaced energy sources.Here the correct design of the control is an ongoing topic of research, e.g., see Refs.[1][2][3][4], especially for so-called grid-forming inverters, which, unlike grid-following inverters, do not rely on the preexistence of a stable grid.Consequently, there is, as of yet, no clear consensus on the proper dynamical modeling of such future sources of energy.
This situation is especially problematic as we expect the number of dynamical actors in a future power grid to increase by orders of magnitude, as energy generation is going to move more and more from the transmission to the distribution level.Furthermore, we expect a larger heterogeneity in dynamics, as the control is no longer structured around established principles, i.e., principles dictated by the physical components.Consequently, it cannot be expected that all the inverters participating in future power grids will use similar control designs.At the same time, the future power grid will face numerous dynamic stability issues, such as low inertia [5], (lack of) timescale separation [6], and greater fluctuations of power production [7,8] and consumption [9,10].The increasing number of dynamical actors and, in turn, heterogeneity mean that, in tackling these issues, we need to understand not only the intrinsic dynamics of individual actors, but more importantly their interactions in large complex networks [11].
This paper addresses these challenges by providing a complexity theoretic approach towards the modeling of future power grids.This means that rather than starting with detailed models of the individual actors and subsequently simplifying them, we focus on those features that are crucial for their interaction on the network.Our approach is based on the physical relationship between voltage and current, which provides the coupling on the network, on the one hand, and on the system desiderata and the symmetry implied by them on the other hand.We find that the most important desiderata are active and reactive power injection, and prescribed voltage levels.
Formulating the dynamics of the nodes in the power grid in terms of the natural invariants associated with their symmetry, we can give an order-by-order nonlinear approximation of their behavior in terms of the deviations from the local desiderata.The key insight is that, by choosing appropriate variables, the lowest nontrivial order in this approximation is able to capture the most fundamental nonlinearities of the power grid dynamics, while also being capable of expressing all local desiderata.The result is a normal form for the behavior of grid-forming actors that resembles controlled Stuart-Landau oscillators, analogous to the form that appears in bifurcation theory [12], parametrized by a latent linear input-output system.From a dynamical systems perspective, the crucial point here is that we fix not only the dynamics of the oscillators to be of Stuart-Landau type, but also provide a specific form for their coupling.
We show by means of comparison to experimental measurements from a sophisticated inverter implemented in the lab, as well as more broadly through numerical simulations, that our normal form indeed captures the quantitative and qualitative properties of a wide range of actors.We also find that it captures important features of the nonlinear interaction that occurs in highly heterogeneous power grids containing both conventional generators and grid-forming inverters.
This normal form provides a starting point for studying realistic models of future power grids from a truly transdisciplinary perspective.By providing a form that closely resembles the Stuart-Landau oscillator, it enables the application of a large range of dynamical system results in the context of power grids.The fact that the normal form is parametrized by a latent linear input-output system enables the use of tools like model reduction and system identification.As opposed to the phase models (e.g., the Swing equation) used in much of the control theoretic and complex systems literature on power grids, the normal form here is based on the instantaneous, physically relevant variables, and provides, order by order, all relevant dynamical aspects of the dynamical actors described.

A. Modeling power grids
We start by introducing our approach to the basic modeling of power grids as complex dynamical networks.This section serves a dual purpose.To make the paper selfcontained to researchers without a strong background in power system modeling, we briefly introduce the most salient physical properties of power grids.We then use this introduction to spell out a sequence of fairly general assumptions we make on our nodes in order to arrive at a highly general model of power grid components.For a more extensive background on overall power grid and inverter modeling, we refer the reader to Refs.[13,14] (and the references therein), which we largely follow in this regard.
Our overall approach is to consider nodes as specifying a voltage that in turn causes a current to flow on the power lines.In electrical engineering terms we thus think of nodes as capacitive elements, and consider our dynamical components as controlled voltage sources.While many of the concrete assumptions below are highly general and apply to all types of power grid components, the assumption of capacitive behavior is most natural for grid-forming voltage source inverters.
Notational aside.-Inaccordance with mathematical and dynamical systems literature, we use i = √ −1 to denote the imaginary unit and use j for the current (this is the opposite convention to the electrical engineering literature).

Node model
Most ac power grids in operation today are three phase, i.e., at every node in the network there are three alternating voltages, V a (t), V b (t), V c (t) : R → R, which cause three separate alternating currents on each transmission element.The sum of the currents flowing from a node into the grid is then denoted I a (t), I b (t), I c (t) : R → R. Our first assumption is that the three phases of both voltage and current are balanced (symmetrical) at all times.That is, we assume that V a (t) + V b (t) + V c (t) = 0 and I a (t) + I b (t) + I c (t) = 0 for all t.This allows us to represent them each by a complex variable with the aid of what is called the Clarke or αβ transformation [15].We denote the complex nodal voltage by u(t) : R → C, with V a = √ 2/3Re(u), V b,c = √ 2/3Re(u exp(±2π i/3)), and the complex nodal current by j (t) : R → C, with I a = √ 2/3Re(j ), I b,c = √ 2/3Re(j exp(±2π i/3)).
Assumption 1 (Balanced phases).The nodal voltage and current are balanced at all times and can thus be described in terms of the complex variables u and j , respectively.
The main objective of ac power grid control is to reach and maintain an operating state where all nodal voltages and currents rotate at a uniform frequency s (usually 50 or 60 Hz), i.e., a quasisteady state described by a limit cycle where at each node we have a nodal voltage u s (t) ∼ exp(i s t) and current j s (t) ∼ exp(i s t).This operating state is further determined by specifying certain set points, i.e., desired values for the amplitude of the nodal voltage, ρ s := |u s |, and active as well as reactive power inputs, p s := Re(u s j * s ) and q s := (u s j * s ), respectively.These set points must be provided by some higher-level control in accordance with the desired power flow in the grid.Since we are primarily interested in the dynamics of the fast acting primary control that takes place at the subsecond scale, we assume that the control dynamics can be described as fully decentralized, i.e., we assume that the set points are given constants and, furthermore, we preclude any additional communication between the individual nodes such that information about the state of the power grid can only be inferred from local measurements.Since we consider controllable voltage sources, this means that the dynamical coupling between the nodes of the network can only be realized by measurements of the local nodal current.

Assumption 2 (Decentralized control). (a) The set points specifying the desired operating state are given constants and (b) the dynamics at a node depend on the other nodes only via the currents on the transmission lines.
While the collective network dynamics with this form of coupling can thus be described in terms of the nodal voltages and currents, an individual node may also have any number of internal dynamic variables.We assume that these are either balanced three-phase quantities (internal three-phase voltages and currents) or scalar quantities (e.g., frequency, dc voltages and currents, or auxiliary variables).For a particular component featuring L of the former and K of the latter, we denote them by z(t) : R → C L and x(t) : R → R K with their respective operating states z s (t) ∼ exp(i s t) and constant x s .This distinction is similar to that between ac and dc variables in Ref. [14] with the difference that we do not explicitly model phase variables.

Assumption 3 (Internal variables). The node dynamics may feature any number of internal variables comprising either balanced three-phase variables denoted by z or scalar variables denoted by x.
Next we assume that the system acts smoothly at the timescales considered.This is justified if typical timescales of switching events in the inverter are smoothed out by filters at the grid connection point, and the dynamics of stabilization of the synchronous component act faster than, e.g., intermittent renewable fluctuations [16,17] that act at the timescales of tens of seconds.If the latter is not the case, we would require modeling the system as including stochastic terms; we leave this challenge to future work.Assumption 4 (Smooth dynamics).The node dynamics is smooth and can be formulated as a system of ordinary differential equations in terms of u, z, and x, with an input given by j , i.e., there are no voltage jumps.Putting all of these assumptions together, the general form of the node dynamics that we are considering here is given by the system u = f u (x, z, z * , u, u * , j , j * ), ( 1 a ) where f u : R K × C L+2 → C, f z : R K × C L+2 → C L , and f x : R K × C L+2 → R K are assumed to be sufficiently smooth functions in the sense that their real and imaginary parts are differentiable analytic functions of the real and imaginary parts of u and j .From this, it follows that they can be written as holomorphic functions of u, u * , j , and j * ; see, e.g., Ref. [18].Equations (1) admit solutions congruent with the required operating state, i.e., for some operating state given by u s and j s , there exist z s and x s for which ) Our final assumption concerns the symmetry of the node dynamics (1).We assume that the node dynamics are homogeneous with respect to phase angles, i.e., there are no distinguished phase angles of u and z such that the dynamics can depend only on relative phase angles between the three-phase variables.This degree of freedom with respect to absolute phase angles translates to a symmetry under global phase shifts, extending the natural symmetry of the desired operating state (2) to the whole phase space.As a global phase shift is equivalent to a time shift for the uniformly rotating operating state, this last assumption is essentially the requirement of time invariance with respect to perturbations of the operating state, a highly desired property of control systems.Note, however, that this assumption is only valid when the transistor switching may be modeled as ideal, i.e., when the three-phase signals do not feature higher harmonics but are purely sinusoidal.This implies that the dynamics must be on a timescale where the switching may safely be neglected by averaging [19,20], which is a standard assumption in power grid control design, considering that the switching frequencies are typically in the range 2-20 kHz [13].

Networked model
When the general form ( 1) is combined with a model for the transmission lines providing the dynamics of the nodal currents, we obtain a connected model for the power grid.While the subsequent derivation of the normal form does not assume a particular model for the transmission lines, for the discussion, we make use of the very simple model of static currents, which we briefly introduce here.A more detailed discussion of dynamical current models is given in Appendix A.
We define the network as a graph G = (V, E) with the set of vertices V = {1, . . ., N }, the set of edges E = {1, . . ., M }, and its complex structure given by the incidence matrix B ∈ {−1, 0, 1} N ×M .Additionally, for each edge, we have its resistance r m ∈ R ≥0 and its inductance m ∈ R ≥0 .
Assuming that the line dynamics evolve much faster than the node dynamics, and that we are in the vicinity of the desired operating state with uniform frequency s , the nodal currents may be approximated by its quasisteady-state equations.In terms of the admittance matrix Y := B diag(r m + i s m ) −1 B T , these are simply given by It follows that we have the active and reactive power inputs or outputs, p n + iq n := u n j * n , given by Note that these equations, as well as the full dynamics from which they arise, are consistent with the symmetry of Assumption 5.This means that the full network model is invariant under global phase shifts.Quasisteady states are again given by an orbit of the symmetry.However a quasisteady state of the coupled model is not guaranteed to exist, unless the set points of the nodes are chosen in a manner compatible with Eq. ( 5).If this is not the case, for example, if there is a power imbalance in the system, the network's quasisteady state will deviate from the orbit specified by the set points, i.e., the desired operating state.

B. Normal form of the node dynamics
System (1) is highly general and captures most models for grid-forming components, as long as they do not explicitly feature higher-level control layers, asymmetric phases, nonsmooth features (like current limitation, as in Ref. [21] for example), algebraic equations without a closed-form solution, or model-free control (e.g., data-based approaches [22]).The central insight is that by exploiting the symmetry we can cast the model in a form that is given by an explicit dependence on the voltage u that is tightly constraint by the symmetry, and a set of meaningful quadratic invariants.The latter uniquely specify different orbits of the symmetry ( 3).An orbit of the symmetry is a set of states related by phase shifts.Thus, the quasisteady states correspond to such orbits.Therefore, the invariants provide us with a sensible notion of the distance to the desired operating state (or any other specified quasisteady state).We can then develop the dynamics order by order in these quadratic invariants.This provides us with a normal form for grid components with dynamics of the form given by Eqs. ( 1) that is valid in the vicinity of the desired operating state.

Derivation
To derive the normal form, we make a temporary change of variables based on quadratic invariants of symmetry (3).Recall the equations for the node dynamics Since the symmetry is defined by a one-dimensional U(1) action on a (2L + 4)-dimensional space (disregarding the scalar variables x for that matter), there are 2L + 3 functionally independent invariants associated with it (see, e.g., Theorem 2.17 of Ref. [23]).These invariants are sufficient to specify the orbits of the group action for an individual node and thus also its desired operating state.Therefore, we want to choose physically meaningful invariants, including, in particular, the quantities that are typically used to define the operation point of the power grid.While there is a certain freedom of choice involved, we propose the set of invariants comprising uu * =: ρ 2 , the voltage amplitude squared, (uj * + u * j )/2 =: p and (uj * − u * j )/(2i) =: q, active and reactive power inputs and outputs, as well as (uz * + u * z)/2 =: ψ and (uz * − u * z)/(2i) =: χ , which may be interpreted as internal active and reactive power flows with respect to the terminal voltage.As the remaining variable that is needed to fully describe the node dynamics, we choose to keep the complex voltage u.
In this combination all invariants are real valued, polynomial, and, together with the complex voltage u, do not require division by the current or any of the internal variables when used to express the original set of variables, i.e., the coordinate transformation is never singular away from the origin u = u * = 0. Employing this new set of variables and defining ξ : The explicit transformation is spelled out in Appendix E.
Note that instead of depending on u and u * we have u and ρ 2 = uu * .Effectively, we have expressed a nonholomorphic function of one complex variable through a function of one complex and one real variable that is holomorphic in the former.Also, note that we did not explicitly write down the remaining dynamical equations for ρ 2 , p, and q, as we are ultimately interested in a normal form expressed in terms of u and j directly and will not thus need these equations in the end.Now symmetry conditions (3) can be stated in terms of the infinitesimal generator of the group action [23] as These equations can be readily solved, e.g., by using the fact that the f are analytic in u and thus we can compare terms in the Taylor expansion of the left-and righthand sides.This provides the explicit dependence of the dynamics on u: with some continuous nonlinear functions g u : R K+2L+3 → C and g ξ : R K+2L+3 → R K+2L .What we have thus achieved is a unified description of all possible node dynamics congruent with our modeling assumptions that is based on the common symmetry of its quasisteady state and preserves the physical relationship between current, voltage, and power.As it turns out, we may describe the node dynamics as a complex oscillator that is augmented by some internal dynamics and is coupled to the other nodes in the network via active and reactive power.Expressing the system in this way now allows for the aforementioned Taylor expansion with a well-defined notion of closeness to some desired operating state.
For notational convenience, we define the vector of instantaneous invariants y(t) := (ξ(t) T , ρ(t) 2 , p(t), q(t)) T : R → R K+2L+3 and a given constant vector y 0 := (ξ T 0 , ρ 2 0 , p 0 , q 0 ) T ∈ R K+2L+3 with ρ 0 > 0 denoting the point around which we want to expand the node dynamics.While the latter can in principle be chosen arbitrarily, of most practical interest are the cases when y 0 represents either (i) a valid operating state considering the entire network, that is, the active and reactive power and the voltage amplitude are consistent with a power flow solution, or (ii) the set points (or also educated guesses thereof if not explicitly available), which may or may not be consistent with a power flow solution.However, in both cases we would have g u (y 0 ) = i s and g ξ (y 0 ) = 0.
Denoting the deviation of y from y 0 by δy := y − y 0 , up to first order we have with the respective coefficients K+2L) , and A u,ξ := g u,ξ (y 0 ).

Normal form of the node dynamics
Having carried out the expansion, we may now return to using the dynamical variables u, u * and j , j * , so that our model may be easily connected by providing equations for the currents flowing through the network according to some model for the transmission lines.We thus arrive at the normal form Note that the back transformation does not change the quality of the approximation; this normal form is still accurate up to terms of order δy 2 .Moreover, since we are dealing with asymptotically stable systems, we expect the error to be bounded and, in most practical cases, quickly decrease over time as long as the trajectories remain within the basins of attraction of both the original system and its normal form.A streamplot example for the phase space of the normal form is depicted in Fig. 1.Voltages are given in per-unit (pu) values.
The rest of this paper will explore the implications and properties of this normal form.

Interpretation
To connect Eqs.(10) to more familiar models for power grids, it is insightful to consider them from the point of view of phase-amplitude coupling in power grids.To this , and H u = 1.For this parametrization, the trajectories converge to a stable limit cycle.The dynamical behavior changes for deviations of the power: no power deviation δp = δq = 0 (left), an active power deviation δp = −0.9,δq = 0 (middle), and a reactive power deviation δp = 0, δq = 0.7 (right).It can be seen that, for this choice of parameters, a change in active power changes the angular velocity (i.e., the frequency), whereas a change in reactive power changes the amplitude of the limit cycle (i.e., the voltage amplitude).
end, consider the complex voltage in terms of phase and the logarithm of the amplitude σ , i.e., u = e σ +iφ and As δξ , δρ 2 , δp, and δq are all real valued, the real and imaginary parts of A u , B u , C u , G u , and H u respectively control the influence of the corresponding terms on amplitude and phase dynamics.For example, in the absence of internal variables the explicit impact of amplitude deviations on phase dynamics can be read off immediately from the imaginary part of C u .However, internal dynamics can mix phase and amplitude reactions in more subtle ways.Form (10) expands the internal dynamics to linear order in the quadratic invariants.The approximate internal dynamics is a linear multiple-input-multiple-output system with two outputs, given by the real and imaginary parts of the right-hand side of u/u, and three nonconstant inputs, given by δρ 2 , δp, and δq.The explicit form is given in Appendix F.
This opens the door to introducing methods from the analysis of linear time-invariant systems, such as model order reduction, to the study of power grid models in a systematic fashion.One way in which we touch upon this possibility later in the paper is by using system identification techniques to fit a model of a fixed complexity to measurement data of a real inverter.This provides semi-black-box models for components (see Sec. II D 1).
We also want to remark that, while we chose to work with the complex voltages (as it makes the expression for active and reactive power particularly simple, is in line with the recently developed concepts of virtual oscillator control [4,24,25], and bears resemblance to an already established dynamics as discussed in Sec.II C 3), it is not mandatory to work with these coordinates.If desired, it is also possible to carry out our approach in terms of u α := Re(u) and u β := (u) or amplitude and phase directly, i.e., u =: ρe iφ .

C. Relation to other models
By approximating the fully nonlinear system (1) with our normal form (10), we have, in principle, reduced all qualitative differences between various concrete models to quantitative differences in their respective coefficients.Thereby, we can distinguish certain classes of models by the coefficients being nonzero.In this section, we introduce a few simple examples of different classes of power grid models and their corresponding normal form, provide working parameter values, and discuss how the normal form relates to the well-known Stuart-Landau oscillator.

Examples
To understand the relationship of the normal form to established low-dimensional models of grid-forming power grid components, we provide the normal form approximation of a variety of them, and demonstrate which coefficients occur in various classes of models.These examples will also be used in the numerical experiments for a heterogeneous network in Sec.II D 3.
We begin by giving a concrete implementation of the abstract derivation of the preceding section for the droopcontrolled inverter model introduced in Ref. [13] [to be self-contained, we give the equations in Eqs.(H2) of Appendix H].In this model the dynamical equation for the frequency, the only internal variable, is already linear with respect to the invariants.As the dynamical equations for the voltage are formulated in terms of amplitude and phase, i.e., ρ = |u| and φ = arg u, we can make use of the relationship to translate them into complex form.This yields Now carrying out the expansion for y 0 = (ω d , (ρ d ) 2 , p d , q d ) T , as prescribed by Eqs.(9), we arrive at The table of normal form coefficients in terms of the original parameters is given as follows: In a similar fashion, we can consider other classes that result from well-known models.
Example 1 (Pure phase oscillators.).Pure phase oscillators are oscillators without any amplitude or internal dynamics.As discussed above in light of Eq. ( 11), the imaginary part of the coefficients provides the phase dynamics, the real part the amplitude dynamics.Pure phase oscillators thus yield the normal form The canonical example from this class is the well-known Kuramoto model [26,27], which further has H u = 0 and assumes that the nodes are coupled by purely inductive transmission lines, i.e., Re(Y) = 0.
Example 2 (Phase-frequency oscillators.).The ubiquitous swing equation [28] (or its nonlinear variant [29]) falls into the class of phase-frequency oscillators.We still have no amplitude dynamics, but now allow for an internal variable: the frequency of the oscillator.The normal form of this type of oscillator is given by As for the Kuramoto model, the standard swing equation further has H ω = 0 and typically assumes coupling with Re(Y) = 0.
Example 3 (Phase-amplitude oscillators.).The quite recent development of virtual oscillator control [4,24,25] represents phase-amplitude oscillators, i.e., there is now amplitude dynamics but no internal dynamics.The corresponding normal form is given by Example 4 (Phase-amplitude-frequency oscillators.).
Synchronous machines can also be cast into normal form.There are both phase and amplitude dynamics as well as internal dynamics.Considering third-order models [30], there is only the frequency as the internal variable, and the normal form for this type of model reads with C u , G u , H u ∈ R if ω should be the true frequency of the nodal voltage.As we see from Eqs. ( 12), the droopcontrolled inverter of Refs.[3,25] also falls into this class.

Normalized parameters values
If one wants to work with our normal form given by Eqs.(10), as an abstract model for grid-forming components detached from any concrete models, there is the question of reasonable parameter values from which to start exploring the parameter space and its dynamic features.To spare the reader the tedious work of trial and error, we thus want to give a starting point here.As the dynamics of the system will highly depend on the given network structure and node types, we here consider the simple case of two identical grid-forming components connected by a RL transmission line.We work in dimensionless units such that |r + i s | = 1 and define tan κ := s /r to account for the ratio of inductivity and resistivity of the transmission line (more details on this are given in Appendix C).Prescribing an operating state with ρ 0,1 = ρ 0,2 = 1 (in dimensionless units) and a relative phase angle of | φ 0 | < cos −1 (1/2), the following sets of parameters yield convergence to the operating state on a timescale of 5-10 s.
(a) Phase-amplitude oscillators [Eq.(13)] For more realistic sets of parameters, we refer the reader to the simulations of Sec.II D, where we numerically compare the normal form to concrete models.The parameters used there are of similar order of magnitude however.

Connection to the Stuart-Landau oscillator
The normal form derived above is closely related to the classical model known as the Stuart-Landau oscillator [31].Key to this connection is that the invariants we chose to represent the deviation from the limit cycle, in particular the voltage amplitude squared, are polynomial with respect to u and u * .
Using the quasi-steady-state approximation for the current equation ( 4), and the normal form for amplitude-phase dynamics given by Eqs.(10) yields the networked system where we have absorbed the expansion points into the coefficients, i.e., It can be seen that the dynamics of the complex voltages u n is that of Stuart-Landau oscillators with a particular nonlinear coupling.Since the admittance matrix Y is given as a linear combination of Laplacian matrices, the coupling between the individual oscillators may be interpreted as diffusive, albeit with a state-dependent diffusion matrix and also involving the complex conjugate voltages of the connected nodes.
This resemblance is of course no coincidence as the Hopf bifurcation (from which the Stuart-Landau oscillator results as the corresponding normal form [12]) prescribes the same U(1) symmetry for the emerging limit cycle in the vicinity of the bifurcation point.Furthermore, the nonlinearity of the amplitude squared arises naturally in the context of a Taylor expansion.This close relationship opens the door to applying methods from the study of coupled Stuart-Landau oscillators, as, for example, in Refs.[32][33][34], to the dynamics of power grids.

D. Validation and probabilistic stability
While the normal form is guaranteed to be a valid approximation in a small neighborhood of the desired operating state, it is not a priori clear whether it can successfully approximate the behavior of real systems under realistic perturbation.Furthermore, if the quasisteady state of the network starts to deviate from the desired operating state of the individual units, we might be stretching the validity of the normal form even further.This section will show initial evidence that real systems (Sec.II D 1) and large perturbations (Sec.II D 2) can be accurately captured by the normal form.Furthermore, we show in Sec.II D 3 that, for an adapted standard IEEE test network, the normal form approximation correctly captures both persistent deviations from the desired operating state as well as probabilistic stability properties.

Lab experiment
We begin with an empirical test of our model using measurement data of a grid-forming inverter with an elaborate control scheme, devised and built at TECNALIA labs [35].The data had originally been gathered to validate numerical simulation tools [36].
We use the normal form as a semi-black-box model that we fit to this measurement data.Note that while a detailed model of the inverter may include many inner control loops, all of which need to be modeled correctly to reproduce the measurements, the semi-black-box model can be chosen to be much simpler.In fact, we use a single internal variable, thus obtaining an effective model of reduced order (dimensionality) for the inverter.Details on the measurement setup and the fitting procedure can be found in Sec.IV.
Figure 2 depicts the results.On the left, we see the measurement data against the trajectories of the fitted normal form with the input that has been used for the fitting.On the right, we show the trajectories for an input from a different measurement while using the same set of parameters, thereby validating the fit on an independent dataset.We see that the normal form is able to capture most of the dynamical behavior of the grid-forming inverter very accurately, only showing overshoots during the sudden shifts in frequency.This is to be expected, however, as we have used a model of quite low dimensionality.

Simulations-infinite bus
We now turn to simulation studies of the normal form approximation in an infinite bus bar setting, as considered in Appendix C. We consider droop-controlled inverters of Schiffer et al. [3] [see Eqs.(H2)] and third-order models of synchronous machines as used by Schmietendorf et al. [30] [see Eqs.(H1)].
Figures 3(a) and 3(d) show the resulting trajectories for a large power perturbation during which the desired power input at the node is doubled for two seconds.We see that the qualitative agreement of the trajectories is excellent in both cases.In fact, the trajectories almost completely match for the droop-controlled inverter.For the synchronous machine, for which the trajectory drops to extremely low voltage levels, deviations are more noticeable.
To explore the reaction to a large perturbation more systematically, we consider various slices of the system's phase space.That is, given the operating state specified in amplitude-angle coordinates by (φ o , ω o , ρ o ), we consider trajectories starting from coordinates of the forms (φ, ω, ρ o ) and (φ o , ω, ρ).We then ask whether the system returns to the operating state from these initial conditions, i.e., we consider slices of its basin of attraction.The model's basin is depicted in yellow, the basin for the normal form in red, and their overlap in orange.
Figure 3(b) shows that, for the droop-controlled inverter, the basins show full agreement.In Fig. 3(c) we see some deviations at low voltage amplitudes very far from the operating state.This is to be expected as the normal form approximation involved expands the voltage dynamics around the desired operating state.In this case the normal form underestimates the stability region of the system.In Appendix II D 3 we see that this is not always the case.
For the third-order model, we have similar results in Figs.

IEEE-14 bus system
To address the question of whether the normal form can also capture the complex interactions between different components, as well as persistent deviations in a realistic power grid, we turn to the IEEE 14-bus test system [37].We adapt this test system by placing various grid-forming components at the nodes, using synchronous machines, different types of droop control, and dispatchable virtual oscillator controls.A detailed setup of the simulation is given in Appendix D.
We then consider the normal form for all of them given their design set points, ignoring deviations that result in the actual networks' operating state due to imbalances and losses.We study the system at a range of power flows by scaling the active and reactive power at the inverter nodes by a common factor f s between 0 and 2. This leads to a variety of operating states that include some deviation from the desired operating points of the inverters.Figure 4 shows in red the minimal voltage amplitude that occurs in the network for these operating states, for both the original and the normal form models.We see that the normal form is capable of describing this behavior very accurately.In black we show the single-node basin stability [38,39] of node 1, i.e., the probability that a random large perturbation at node 1 destabilizes the system.Again, the results for the normal form and the full model agree within the uncertainty bounds.
This demonstrates that the normal form, by taking a "network first" approach to modeling and keeping the power flow equations fully accurate, is capable of capturing sophisticated properties of complex highly heterogeneous power grids.

III. DISCUSSION
This paper introduces a normal form for grid-forming power grid components.This form is arrived at by using symmetry arguments to restrict the functional form FIG. 4. IEEE 14-bus test system [37].The system is studied for a range of power demands, and all active and reactive power set points at the inverter nodes are scaled by a common factor f s between 0 and 2. In red we plot the minimal voltage amplitude ρ min = min n ρ o n in the resulting operating state, the solid line is the original model, and the dashed line is the normal form.In black we give the single-node basin stability of node 1.
and expanding order by order in physically meaningful quadratic nonlinearities, concretely the power and voltage mismatch at the nodes.At lowest order the normal form is parametrized by a linear time-invariant system with three time varying inputs and two outputs.The main factor in the complexity of the normal form is the number of internal states of this linear system.
The normal form can be derived from more detailed analytic models, for which we gave detailed examples, but it is also possible to directly infer it from experimental measurements.In the latter case, we can fix the number of internal states a priori to obtain an empirically best model of the system at a given complexity.We give a proof of concept of this approach by fitting the data of a gridforming inverter built at TECNALIA to a low complexity normal form with one internal variable.
A more systematic exploration of this approach will require adapting tools from system identification to this context, especially for dealing with noise in the measurements.
We saw in numerical experiments that the normal form is capable of describing the nonlinear behavior of the power grid in the vicinity of the desired operating states.This was explored for both single machines at an infinite bus and a highly challenging heterogeneous network of diverse grid-forming actors.While the quality of the lowest-order approximation for a single node is relatively clear from the derivation and the numerical experiments, a more thorough understanding of the limits of the approximation when we consider a whole network of oscillators will require more work.
Besides being of interest in itself, the normal form presented here also provides a starting point for the transdisciplinary study of realistic models of future power grids [40].
The form closely resembles Stuart-Landau oscillators, thus opening the door to adapting a large body of dynamical systems research to the study of future power grids.As it is based on very general principles and physically meaningful variables, all relevant dynamical aspects of the dynamical actors can be described by it.This opens up the possibility to transport results from theoretical research on control and complex system aspects of power grids, often based on highly conceptual phase models [41], to models that are accurate with regard to the real power grid.
The normal form also opens up further novel research avenues.For example, it can serve not just as a model for concrete systems but as a specification for the behavior of future designs.The study of the linear stability of the normal form can be considered a proof of concept in this direction.
While this work focuses on grid-forming nodes, we expect the approach to be fruitful more broadly.By choosing different invariants, and different variables, it is possible to arrive at normal forms that will be suited to other classes of grid actors.Nonsmooth behavior might also be modeled as switching between different normal forms.
Finally, the mathematical approach taken is highly general.Whereas ordinary phase reduction approaches require a small coupling assumption, here we were able to work from the assumption that coupling and oscillation are well adapted to each other.It is rare for oscillating systems to exist and develop in isolation, and we posit that such an approach is likely to be fruitful in other fields of complex systems science, synchronization, and oscillator networks.

IV. MATERIALS AND METHODS
For the empirical validation of our approach, we use measurement data of a grid-forming inverter devised and built at TECNALIA labs.The inverter control design basically consists of conventional droop control with a low-pass filter for the measured power output that emulates inertia.This basic design is similar to Eqs. (H2), but further includes additional filters for voltage and frequency measurements, as well as a virtual impedance [35].
Besides the inverter, the lab setup contains an ac power source and an Ohmic load.The load is connected to the power source by an emulated line containing a series of inductances and resistances.The inverter is connected to the load by a transformer.More details on the setup and the parameterization of the components can be found in technical report [42].In this particular test case, we vary the voltage angle frequency at the power source and measure the voltages and currents for two of the three phases directly at the inverter.The measured time PRX ENERGY 1, 013008 (2022) Re(u) lm(u) FIG. 5. TECNALIA inverter measurement.The depicted lines are the measured time series for current and voltage at the gridforming inverter described in Ref. [35].
series are depicted in Fig. 5. Assuming the three phases to be balanced, we can thus directly calculate the complex nodal voltage and current (u and j ), and further the active and reactive power outputs (p and q), as well as the nodal voltage amplitude ρ.The frequency can be determined by numerical differentiation of the voltage phase angle.
We take the normal form with the voltage angle frequency as the only internal variable and fit to this measurement data.To reliably fit such models, it will be necessary to properly adapt system identification techniques to this setting.To obtain a proof of concept, we instead opt for a straightforward two-step approach with generic tools.First, we perform a linear regression for obtaining rough parameter estimates, then we fine-tune these using scientific machine learning tools.For the linear regression, we first numerically calculate the derivative of the complex voltage and the frequency to obtain the left-hand side of differential equation ( 14) and subsequently get an estimate for the parameters using the method of least squares.For the fine-tuning, we use the current signal as a data-driven input for a dynamical simulation of the normal form model using the DifferentialEquations.jlpackage [43] and optimize the least-squares fit of the trajectories with stochastic gradient decent using the DiffEqFlux.jlpackage [44].
The code to reproduce the results and figures of this paper is available from a Zenodo archive [45] or a GitHub repository [46].

ACKNOWLEDGMENTS
We would like to thank Meng Zhan and Sebastian Liemann for detailed comments on a draft of this manuscript.The authors acknowledge support from the BMBF, CoNDyNet2 FK. 03EK3055A.This work is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Grants No. KU 837/39-1, No. RA 516/13-1, and No. HE 6698/4-1.All authors gratefully acknowledge the European Regional Development Fund (ERDF), the German Federal Ministry of Education and Research, and the Land Brandenburg for supporting this project by providing resources on the highperformance computer system at the Potsdam Institute for Climate Impact Research.

APPENDIX A: LINE DYNAMICS
As an addendum to Sec.II A, we want to briefly discuss the applicable line models for our normal form.As we take nodes to be a voltage that reacts to a current, that is, it behaves like a capacitor, our lines need to be modeled as providing a current in reaction to the terminal voltages, thus behaving like inductances.
Using the same notation as in the main section, we only need to include further dynamical variables corresponding to the currents flowing on each transmission line denoted by which would replace the algebraic relationship between voltage and current, Eq. ( 4), that we used in the main section.The latter is actually derived by considering the quasisteady state where j e,m (t) ∼ exp i s t, i.e., setting d dt j s e,m = i s j s e,m .
More sophisticated models of lines, which include line capacitances, can be naturally coupled to our nodal ordinary differential equations (ODEs) as long as they provide ODEs for the terminal currents.This is the case, for example, for τ models and iterated τ models of transmission lines.Models such as the π model that have ODEs for the terminal voltages lead to algebraic constraints on the system.We also want to note that in the special case of a uniform ratio between resistance and inductance across the whole network, i.e., m /r m =: τ for all m ∈ E, the current dynamics may be directly expressed in terms of the admittance matrix Y and also greatly simplified with respect to dimensionality by eliminating the line currents and writing the nodal current dynamics as

APPENDIX B: FITTED PARAMETERS
We get the following parameters by fitting the model as described in Sec.II D 1: All parameters have the units 1/s in the per-unit system with P base = 10 kW and V base = 393.4V.

APPENDIX C: LINEAR STABILITY (INFINITE BUS BAR)
A crucial aspect of obtaining a normal form is that it allows us to make highly general analytic statements that apply directly (if approximately) to a wide range of potential power grid components.To demonstrate this point with a proof of concept, and to further improve our understanding of the coefficients in Eqs.(10), we consider the linear stability when connected to an infinite bus (or slack node) via Eq.( 5).We limit ourselves to the normal form with a single internal frequency variable, i.e., Eqs.(14).
We work in the reference frame corotating with the infinite bus and fix its phase angle at zero such that the nodal voltage at the infinite bus is given by a constant V s ∈ R >0 .Our goal is to derive conditions for the parameters that ensure local asymptotic stability for some valid equilibrium point with synchronized frequency, i.e., some y 0 = (0, ρ 2 0 , p 0 , q 0 ) T for which A u (y 0 ) = A ω (y 0 ) = 0, and such that there exists ϕ ∈ [0, 2π) : p 0 + iq 0 = Y * (ρ 2 0 − ρ 0 V s e iϕ ).Note that, with a slight abuse of notation, Y denotes the admittance of the single transmission line here.For convenience, we make the change of coordinates σ R + iσ I := ln u, with the subscripts R, I denoting real and imaginary parts in the following, and set B u I = 1 without loss of generality.The system we are considering here is thus given by with the Jacobian and the constants p ± 0 , q ± 0 defined as Invoking the Routh-Hurwitz criterion [47] we can ensure that all three eigenvalues lie in the left complex half-plane if the inequalities trJ (y 0 ) < 0, (C3a) trJ are satisfied [see Eqs.(G1) in Appendix G for these inequalities in terms of the parameters].While these conditions are necessary and sufficient, they are too intricate to yield any qualitative insights, so we consider the special case in which the response of the node to active and reactive power is adapted to the behavior of the power line: for some k u > 0, k ω > 0, and tan κ := −Y I /Y R .This makes the coupling behave like conventional droop control [3] for a purely inductive network (Y R = 0) even when Y R ≥ 0 (an idea which that been used, e.g., for dispatchable virtual oscillator control [4]).To state the sufficient stability conditions, we first define the short forms encoding the ratios between the coefficients specifying the system's reaction on voltage amplitude and power deviations, as well as the ratio between the infinite bus voltage and the desired voltage amplitude at the node (although in practical cases we usually have R V ≈ 1/2).Additionally, we define the angle With these definitions and assumptions (C4), we can state that system (C1) is asymptotically stable if the following conditions are satisfied: The derivations of these inequalities can be found in Appendix G.The conditions show that the chief determinants for the stability of a certain equilibrium (with relative phase angle ϕ and voltage amplitude ρ 0 ) are the ratios between the coefficients specifying the system's reaction on voltage amplitude and power deviations.A geometric representation of inequalities (C5) and (C6) is depicted in Fig. 6.Assuming that V s ≈ ρ 0 (as is usually the case for power grids), it can be seen that, for relative angles |ϕ| < π/3, it is enough to ensure the correct sign for the coefficients C u and C ω , i.e., C u < 0 and sign(ϕ)C ω ≥ 0, while respecting the bound on R ω .Here R ω parametrizes the impact of the amplitude on the frequency of the system, relative to the power droop.Thus we find that in a moderately loaded scenario, the crucial factor for the linear stability of the system is the amplitude-phase coupling.If a greater load, and thus a greater relative phase angle, needs to be guaranteed stable, we must further ensure that the influence of voltage amplitude deviations dominates that of power deviations, as given by inequality (C6).Lastly, we note that inequality (C6) is actually a necessary bound under assumptions (C4).Inequalities (C5) are only sufficient, i.e., they may be relaxed by invoking stricter bounds on B ω and B u R , which are rather technical however [see Eq. (G4)].Now, given a concrete model from the class of amplitude-frequency oscillators, we only have to substitute the coefficients with the corresponding partial derivatives as given by Eqs. ( 9) in order to translate the stability conditions to the specific model parameters.

APPENDIX D: NUMERICAL SIMULATIONS
This appendix describes the model and simulation setup for Sec.II D, as well as some further results for the networked case.All simulations are performed using Dif-ferentialEquations.jl [43].The network model is built with PowerDynamics.jl [36].The simulations are performed with a RADAU solver [48] with relative tolerances set to 10 − 6; for more details, we refer the reader to the code accompanying the paper.

Infinite bus bar
In the infinite bus simulations of Sec.II D 2 we consider a model of a droop-controlled inverter, Eqs.(H2), and a third-order approximation of a synchronous machine, Eqs.(H1), connected to an infinite bus.We use pu units with voltages 1 at the nodes and a per-unit power chosen such that the admittance of the line, which we take to be purely inductive, is Y = −1i.The line model also includes shunt capacitance of Y s = 0.2i.For the inverter model (H2), we choose the time constant τ p = 2.5 and droop gains k p = 5, k q = 0.1, and expand around ρ 0 = V d = 1, ω 0 = ω d = 0, p 0 = p d = 0.5, and q 0 = q d = 0.2.For the synchronous machine model (H1), we choose the damping constant γ = 0.2, time constant α = 2, and internal reactance X = 1, and expand around ρ 0 = E f = 1, ω 0 = 0 (the model is given in the corotating reference frame), p 0 = p m = 0.5, and q 0 = 0 as no explicit set point for reactive power is given in the model.

Network
The network model is based on the IEEE-14 bus dynamical test system, augmented with grid-forming components.The distribution of components is shown in Fig. 7.We use synchronous machines at nodes 1, 3, 6, and 8 [Eqs.(H1)] and inverters with different types of control at most other nodes, namely at 4, 5, 9, and 12 we use the droop controlled inverter (H2), at nodes 11 and 14 we attach the inverter model (H3) (which adds inertia to the amplitude dynamics), and at buses 10 and 13 we place inverters with dispatchable virtual oscillator control (H4).For the detailed parameter choices, we refer the reader to the code.
To stress the system, we scale the active and reactive power demands at all inverter controlled nodes by a common factor f s varying between 0 and 2. Note that the set points do not provide a solution of the power flow, and we have persistent deviations from the desired operating points.
The main part of the paper discusses the single-node basin stability for a variety of f s values.Here we also briefly present more phase space slices in the style of Fig. 3.Each location in the slice corresponds to an initial condition for two variables of the system.All other variables in the system are initialized at the quasisteady operating state of the network.Thus, the plots show very large, instantaneous perturbations affecting only a single node.The results are shown in Fig. 8.
We see that in this highly challenging scenario with large perturbations, the agreement with respect to the shape of the stability region becomes significantly worse the further we stray from the operating state.The most problematic deviations appear in the voltage amplitude.This should be kept in mind, when analyzing fault scenarios that feature such deviations.However, many important qualitative features are still captured by the normal form.

APPENDIX F: LINEAR TIME-INVARIANT INPUT-OUTPUT FORM OF THE INTERNAL DYNAMICS
We can write Eqs.(10)  This form is the most suitable for using tools from the study of linear time-invariant systems for power grid models.

APPENDIX G: LOCAL ASYMPTOTIC STABILITY
Full inequalities (C3) in terms of the normal form coefficients and the expansion point are given by the rather FIG. 8. Various basin slices for the IEEE 14-bus test case.Color code as in Fig. 3. Top row: bus 1 (left), bus 4 (middle), bus 6 (right).Bottom row: bus 8 (left), bus 9 (middle), bus 12 (right).lengthy expressions 0 > 2ρ 2 0 C u + G u p + 0 +H u q − 0 +B ω , (G1a) Employing the additional assumptions in Eqs.(C4) and the definitions of R u , R ω , and R V yields the more compact inequalities with From these inequalities we can immediately deduce the conditions given in Appendix C. By requiring B ω < 0, inequality (G2) is satisfied if W 1 < 0, which yields the first inequality of inequalities (C5).By further requiring |ϕ| ≤ π/2, inequality (G4) is satisfied if 0 ≥ W 2 sin ϕ, which is equivalent to sign(ϕ)B u R ≥ 0 and the second inequality of inequalities (C5).For inequality (G3), we make use of the trigonometric identity Since W 1 < 0 implies that R u > 0, this yields inequality (C6).

APPENDIX H: MODELS USED IN SEC. II D
For completeness, we give the models we base the heterogeneous network in Sec.II D on.Here the parameters are kept in line with the notation used in the original papers that they were taken from.We also provide the normal form coefficients when expanded around the design set points.
(a) Third-order approximation of synchronous machines [30]: q n E n , (H1) FIG.1.Streamplot example for the normal form model without any internal variables.The parameters are A u = 1 + i, C u = −1, G u = i, and H u = 1.For this parametrization, the trajectories converge to a stable limit cycle.The dynamical behavior changes for deviations of the power: no power deviation δp = δq = 0 (left), an active power deviation δp = −0.9,δq = 0 (middle), and a reactive power deviation δp = 0, δq = 0.7 (right).It can be seen that, for this choice of parameters, a change in active power changes the angular velocity (i.e., the frequency), whereas a change in reactive power changes the amplitude of the limit cycle (i.e., the voltage amplitude).

FIG. 2 .
FIG. 2. Left: model fit to data, with one internal variable.Right: validation of the fitted parameters against different test runs.Blue denotes measurement data.Orange denotes the output of the normal form.
3(e) and 3(f), with a slight overestimation of stability in some areas.The qualitative features of the original model are fully reproduced by the normal form though.The details of the numerical setup are given in Appendix D.

FIG. 3 .
FIG. 3. Results for the infinite bus scenario.Trajectories for the (a) inverter, (d) synchronous machine (the vertical bars signify the beginning and end of the perturbation at t = 10 s and t = 12 s).Basin cross sections for the (b),(c) inverter, (e),(f) synchronous machine with four possible cases for each state: points inside the basins of both the original model and its normal form (orange), only inside the original model's basin (yellow), only inside the normal form's basin (red), and outside both basins (white).The crosses depict the actual operating state (φ o , ω o , ρ o ) of the full models.
j e,m (t) : R → C. From the standard laws of electrical circuit elements, a line modeled as a resistor and an inductance in series satisfies m d dt j e,m = −r m j e,m + N n=1 B nm u n ,