Complex Couplings - A universal, adaptive and bilinear formulation of power grid dynamics

The energy transition introduces new classes of dynamical actors into the power grid. There is especially a growing need for so-called grid-forming inverters (GFIs) that can contribute to dynamic grid stability as the share of synchronous generators decreases. Understanding the collective behavior and stability of future grids, featuring a heterogeneous mix of dynamics, remains an urgent and challenging task. Two recent advances in describing such modern power grid dynamics have made this problem more tractable: First, the normal form for grid-forming actors provides a uniform, technology-neutral description of plausible grid dynamics, including grid-forming inverters and synchronous machines. Secondly, the notion of the complex frequency has been introduced to effortlessly describe how the nodal dynamics influence the power flows in the grid. The major contribution of this paper is to show how the normal form approach and the complex frequency dynamics of power grids combine, and how they relate naturally to adaptive dynamical networks and control affine systems. Using normal form and complex frequency, we derive a remarkably elementary and universal equation for the collective grid dynamics. Notably, we obtain an elegant equation entirely in terms of a matrix of complex couplings, in which the network topology does not explicitly appear. These complex couplings give rise to new adaptive network formulations of future power grid dynamics. We give a new formulation of the Kuramoto model with inertia as a special case. Starting from this formulation of the grid dynamics, the question of the optimal design of future grid-forming actors becomes treatable by methods from affine and bilinear control theory. We demonstrate the power of this perspective by deriving a quasi-local control dynamics that can stabilize arbitrary power flows, even if the effective network Laplacian is not positive definite.


Introduction
Due to the energy transition, conventional power plants are replaced by renewable energy sources (RES) such as wind and solar PV.RES are predominantly connected to the power grid via power electronic inverters.Furthermore, power electronic interfaces are also in use for many loads, and for connecting storage solutions such as batteries to power grids.Currently, most inverters are grid-following and have to rely on a stable grid to function.However, as we move towards a fully renewable grid, there is a growing need for so-called grid-forming inverters [2] (GFIs) which can contribute to the grid stability independently of conventional generation.Popella et al. have shown that following one of the scenarios in the German Network Development Plan [3], more than 80% of all inverters installed from 2021 need to be grid forming to assure a stable power grid in 2035 [4].This shift presents a significant challenge, as GFIs are a relatively new and complex technology, of which there is a limited practical and theoretical understanding.
In future grids, these GFI will coexist with a much reduced number of conventional generators.The collective dynamics of such inverter-dominated networks remain a challenging topic [5].For almost all proposed types of GFIs, there is an analysis of the stability of the individual machine.The stability of the entire network, however, is rarely addressed, with [6] and [7] as notable exceptions.Most research on the collective behavior of the power grid is still based on the paradigmatic Kuramoto model [8].For generators, we already know from the study of higher-order models that include voltage dynamics, that the Kuramoto Model can be misleading, for the individual dynamics [9] but also for the collective phenomena [10].As we move towards an inverter-dominated grid, we even lack the consensus on higher-order models of GFIs, as vastly different realizations of gridforming capabilities are possible.However, important phenomena occurring in inverter-dominated networks can not be captured by the Kuramoto model.To study collective phenomena in heterogeneous systems, a description of all grid-forming actors, conventional and future GFI-based, is needed.Kogler et al. [11] address this, by introducing the normal form, a dynamic technology-neutral formulation for the dynamics of grid-forming actors.This normal form approach provides a concise and meaningful parametrization of the space of plausible power grid actors, especially when they are near their desired operational state.Thus, it is well suited to the study of power systems with a heterogeneous mix of dynamical actors.Independently, Milano recently introduced the notion of the complex frequency [12] to elegantly describe how the nodal dynamics influence the power flows in the grid.In this paper, we show that the two formulations introduced above complement each other perfectly.Together, they provide an elementary and universal formulation of the space of plausible power grid dynamics.Further, by focusing on the evolution of the couplings between nodes, we can derive a new and surprisingly elegant formulation of the collective grid dynamics.This formulation has the remarkable property that the network adjacency matrix no longer explicitly appears in the equations, but is encoded in the initial condition.This immediately gives rise to an adaptive network formulation of the power grid dynamics that is different from that employed in [13].As a special case, we derive the Kuramoto model with inertia in this formulation, but the general formulations provided again make no simplifying assumptions such as lossless coupling or specific voltage dynamics.The power and voltage magnitude square, the most important control objectives of power grids, are linear combinations of the complex coupling.Thus, if we use the complex coupling dynamics as the dynamics to be controlled, we obtain an elegant quadratic-bilinear tracking problem.Solving this problem under the constraint of a fully decentral controller using only local information is challenging; however, the formulation provides a natural quasi-local Lyapunov-Function-Controller that is globally synchronizing while reducing the tracking error.A major contribution of this paper is to systematically collect and describe the various ways in which complex frequency dynamics of power grids, normal form descriptions of grid-forming actors, and concepts from control theory and dynamical systems intersect.To provide the reader with a background of the relevant fields we will begin with recalling the most relevant concepts from power grids, adaptive dynamical networks, and bilinear control systems in section 2. In the second part 3.1 we introduce the new formulations of the grid dynamics.Finally, we explore the control aspects in section 4. We provide an appendix with various useful related calculations.

Balanced AC Power Grid Variables
Generally, power grids consist of several parallel circuits carrying phase-shifted AC voltages and currents.We will focus on the typical case of three phases.A three-phase signal 0 at all times.Balanced signals can be expressed in terms of a rotating 2D vector, often referred to as dq-coordinates [14].Defining the complex state v = v d + jv q , also called the Park vector, we can write the original signal as V l = Re(exp (jl2/3π)v) where l depicts the different phases.The complex voltage v can be written in terms of phase θ and amplitude V ≥ 0 as: Note that, with our definition, the amplitude V is the (peak) voltage magnitude of the three-phase signal.As we will assume balanced conditions throughout, the Park vector v(t) captures the full state space of the system.The same procedure is used to express the three-phase currents I h in terms of a single complex current i(t).
It is then a straightforward calculation to see that the real part of vi is 3 2 times the instantaneous power contained in a three-phase current i at voltage v.The imaginary part is called the instantaneous reactive power and is interpreted as the amount of energy that is injected into the transmission line without flowing through it.The full quantity S c = 2 3 vi is called the complex power.In what follows we will absorb the factor 2  3 into the coefficients and abuse terminology slightly to speak simply of as the complex power.

Complex Frequency
The frequency f is the most important variable for operating, controlling, and monitoring power grids [14].According to the accepted IEEE standard [15] the frequency f (t) of a signal x(t) with amplitude X and phase ϕ is defined as: This standard definition of frequency in AC power grids is only meaningful when a state with slowly varying amplitude is assumed.However, in power grids, fast amplitude changes in the voltage are possible.This standard definition can not distinguish between phase and amplitude variations in a meaningful way.The objective of the complex frequency, which has been introduced in [12], is to overcome this issue by providing an interpretation of the instantaneous frequency as a state space velocity.As long as it is non-zero, we can consider the time derivative of the logarithm of the voltage v(t) (1): here ω is the angular velocity of the signal and ρ is the rate of change of the log-voltage amplitude.
The complex number η is defined as the complex frequency [12].For slow-changing amplitudes, where both definitions are applicable, the standard definition of the frequency coincides with the imaginary part of the complex frequency.
The complex frequency has already found numerous applications ranging from the inertia estimation of virtual power plants [16] to analyzing the linear [17] and nonlinear stability [18] of networks which are composed of dispatchable virtual oscillators [19], a proposed type of GFI.Furthermore, [20] shows the connection between the complex frequency and power for several different converter controllers, including GFIs.The authors of [17] introduce complex frequency synchronization, a novel concept of synchronization.Equation (6) shows that to obtain full synchrony of the voltages, both the angular frequency and the rate-of-change-of-voltage have to match.During rate-of-change-of-voltage synchronization, the voltage amplitudes are still allowed to change, but at the same exponential rate.

Normal Form
The authors of [11] have introduced the normal form, a technology-neutral formulation for the dynamics of grid-forming actors.The main result of [11] is that the symmetry under global phase shifts of the nodal dynamics can be exploited to formulate the dynamics in terms of invariants of the symmetry.The entire interaction between grid-forming actors and the grid can be written in terms of the active P and reactive power Q, as well as the square of the voltage magnitude ν = vv, which form an elegant complete set of invariants for power grids.P , Q, and ν are exactly the quantities needed to describe the desired behavior of grid-forming actors, as these are responsible for balancing active power and providing enough reactive power to maintain the nominal voltage levels.Following [11] the dynamical behavior of grid-forming actors can be generally written in this form: where v is again the complex voltage, x c are the internal dynamics of the control of the grid forming actors and g v , and g xc are differentiable, potentially nonlinear, functions.Using the complex frequency (7) we see that we can rewrite the first equation simply as: which gives the dynamics a novel interpretation.The complex frequency η = v v determines, by definition, the voltage dynamics.As η is invariant under global phase shifts, it can only depend on other invariants if the entire system is to be invariant.In particular, η can not depend explicitly on v.

Adaptive Dynamical Networks
Adaptive dynamical networks are a class of systems that change their coupling K hm (t) over time depending on the dynamical state y h of the nodes of the network.The review paper [21] gives an excellent overview of the applications, dynamical phenomena, and available mathematical methods for adaptive dynamical networks.Furthermore, [22] features a rundown of the concept of adaptivity and how it applies across different scientific disciplines.In this work, we will focus on weighted adaptive dynamical networks that are in general defined as: where f h describes the local dynamics of node h, γ is the coupling function and H the is the adaptation function.
Recently, it has been revealed that dynamical power grid models, based on phase oscillator models are connected to adaptive networks [13].In particular, it has been shown that the Kuramoto model with inertia and the Kuramoto model with inertia and voltage dynamics [23], in electrical engineering referred to as a third-order machine model [9], can be written as a system of adaptively coupled phase oscillators.Using this approach the voltage dynamics can be interpreted as an additional adaptivity term.It has been shown that phase oscillator models and oscillators on adaptive networks share many different dynamical phenomena such as solitary states [13,24].
In this work, we will introduce a type of complex coupling K hm that allows us to represent general dynamical power grid models as adaptive networks as long as the complex frequency η is well defined.

Bilinear and Control Affine Systems
Power grids are inherently non-linear systems.Therefore, to control them, non-linear control theory with all its associated difficulties is required.Instead, a linear approximation of the system could be used, but then non-linear effects cannot be taken into account and the controller will have limited effectiveness.
The complex frequency and normal form equations reveal that the power grid can be interpreted as a system in which the control variable is linear.Such systems are called control affine systems [25], and are of the form: where y and u are the state and control vector, respectively.The power grid equations we will find are bilinear systems [26], meaning that F(y) and G(y) are linear functions of y.This leads to the bilinear form of the control system: where y and u are again the state and control vector, respectively and N a are the system matrices.Bilinear systems are an area of intense research as they occur in various contexts [27], such as nuclear reactors.There is considerable work on bilinear control systems [26].Below, we will see that exploiting the control affine structure already allows us to derive interesting new control laws.

New Power Grid Equations
We begin by introducing the dynamics of the complex coupling, before deriving new adaptive formulations of power grids.

Complex Coupling Dynamics
In [12] Milano considers the dynamics of the complex power S h in terms of a dynamical coupling term.In Milano's paper, these coupling terms are referred to as s hm .To avoid confusion with the power flows on the lines, we will use K hm in this work.Take Y hm = 1 Z hm as the admittance on the line connecting nodes h and m.We then introduce the admittance weighted graph Laplacian L hm as: The dynamical coupling terms K hm , that will play the role of the adaptive coupling in the sense of equation ( 12), are defined as: where L hm is the complex conjugate of the Laplacian.These coupling terms K hm do not have a physical interpretation, but they have the important property that: holds, as: Using the voltage dynamics (7) the dynamics of the coupling are given by: Thus the dynamics for the nodal complex power S h is given by: which is equation (32) in [12].The diagonal coupling terms K mm describe the evolution of the local voltage magnitude square ν m up to a constant factor: From equation ( 21) we can see that the system has two possible sets of fixed points.The first set of fixed points coincides with K hm = 0 which can either be achieved when h and m are not connected by a transmission line or when there is a total voltage collapse at either node.The second type of fixed point is defined by η h = −η m throughout the connected component.If the system is connected at the fixed point, that is, for every two nodes there is a walk of non-zero K hm connecting them, then this implies that η h = jω global for all h.This shows an immediate advantage of considering a dynamical system given entirely in terms of invariants, the desirable operating states are fixed points rather than limit cycles, which simplifies system identification and control synthesis considerably.The stability of the synchronous fixed points can now be expressed using a master stability function [28], an approach that has recently been extended to a large class of adaptive networks [29].
Following [12] we assumed that the transmission line dynamics are faster than the nodal dynamics and that a quasi-steady state model can be employed for the transmission lines.In Appendix B we show briefly how to extend this description to more realistic line models.Furthermore, we point out that a similar, albeit less concise formulation can be achieved in power flow variables.The derivation can be found in the appendix as well C.The coupling dynamics ( 21) can be used to rewrite general power grid dynamics as adaptive dynamical networks, as introduced in section 2.4, as long as the complex frequency dynamics are well-defined.In the following, we will perform this reformulation for the classical Kuramoto model with inertia 3.2.Using the normal form 3.3, rather than the dynamics of the complex frequency, we derive an elegant equation entirely in terms of the matrix of complex couplings, in which the network topology does not explicitly appear.These reformulations generally come at the price of a larger phase space, which is typical for adaptive networks.In general, the phase space has 4N + 4E dimensions, where N is the number of nodes and E is the number of edges in the power grid.The normal form is a notable exception to this rule, as demonstrated in section 3.3, the phase space has a dimension of only 4E.

Adaptive Kuramoto Model with Inertia
The Kuramoto model with inertia is the classical model to describe a generator in a power grid [9].
It is used to describe the short-term behavior of the generator, also called the first swing.For this reason, it is also called the Swing equation in this context.It is given by: where θ h and P s h are the phase angle and active power set point of node h.D h is the damping coefficient and H h is the inertia constant.As the Kuramoto model with inertia has no voltage dynamics the complex frequency is purely imaginary η h = jω h .Using the relation between the coupling and the complex power (17) and the coupling dynamics (21) we can rewrite the Kuramoto model with inertia as the following adaptive network: It is important to note that this form of adaptive power grid dynamics is not the same as the formulation introduced in [13].As the complex coupling dynamics given above are formulated entirely in terms of invariants of the limit cycle, the phases do not appear explicitly.Thus synchronization at arbitrary frequencies is represented here by convergence to different fixed points, rather than limit cycles.For cluster synchronization, where several frequencies coexist, the clusters are internally approximately at fixed points, while the coupling between them oscillates rapidly.This formulation thus automatically realizes a fast-slow separation of variables for these common asymptotic states.

Self-contained complex coupling dynamics
By considering the equations of η, the normal form can be rewritten as an adaptive equation in a similar way to the Kuramoto model above.Instead, we can also directly combine the normal form equations for η and the complex couplings K km .These equations provide a self-contained dynamical system expressed entirely in terms of invariants of the limit cycle.For notational simplicity, we first neglect the internal variables.Then the entire power grid's dynamics are given by: Remarkably this system of differential equations depends only on the diagonal of the Laplacian L hm , i.e. the weighted degrees.The line admittances are only encoded in the initial conditions.
If there are internal degrees of freedom of the normal form, then there are dynamical variables at the node as well.Further, expressing g h in terms of S = P + jQ and S = P − jQ directly, we obtain the remarkable class of matrix dynamical systems: To illustrate the equivalence of the self-contained dynamical system, given by ( 29), to the dynamics of a full power grid we performed a simulation example.Using the algorithm introduced in [30] we generated a 3-bus power grid consisting of normal forms as the bus model.Figure 1 shows a comparison between the trajectories P, Q and ν of the test grid and the equivalent coupling system (29).Using equations ( 23) and ( 20) we calculate the voltage magnitude v and the active and reactive power from the coupling terms.From figure 1 we can find that the full power grid and the coupling dynamics give the same results for the entire trajectory.
Figure 1: Comparison between the power and voltage given by the dynamics of a full network and the coupling dynamics given by equation ( 21).It can be seen that both the coupling dynamics and the full network give the same results.
The possibility of expressing power networks via these elemental coupling dynamics provides many opportunities for future research, e.g. for stability analyses of the collective dynamics of future power grids.In the following, we will use the coupling terms to design complex frequency dynamics that stabilize the power grid.

GFI Design as a Bilinear Control Problem
In [11] the normal form is introduced as a way to parametrize the space of plausible grid actors.This allows us to obtain a unified description of inverter designs and existing machines.However, it is important to also ask which points in this space would provide stable dynamical laws for the power grid.We assume that we have a perfect voltage source that we can steer freely, an assumption that is typical in the study of grid-forming control [6,7].Below, we will see that the complex coupling variables allow us to cast this problem into the form of a bilinear-quadratic control problem.We also derive a quasi-local Lyapunov-function-based controller that is globally synchronizing.

Grid Forming Control as a Normal Form Tracking Problem
As noted above, the variables used in the normal form to express the coupling of the grid into the grid-forming actor are exactly those that a grid-forming actor seeks to control.We can assume that the grid-forming actor has (possibly time-dependent) set-points ν s , P s , Q s for the desired active and reactive power and voltage amplitude square.Then we can write the normal form in terms of the error coordinates e(t), the difference between the desired and actual quantities.The control problem of grid-forming actors is thus reduced to a tracking problem with a feedback loop.
The set-points act as the reference r(t) which should be tracked, and the complex frequency is the control input u(t) by which the power flow on the grid should be controlled.This is illustrated in 2.
Figure 2: A block diagram that represents a single grid forming actor coupled to the power grid.
The GFI prescribes a control input u, the complex frequency, to the power grid, which then leads to an observed output y(t), the power flow and voltage, which is used as control feedback by the normal form.
If we require the feedback loop to be a linear time-invariant system, we arrive exactly at the internally linearized normal form equations already discussed in [11], that are valid near the desired state.

Complex Voltage
Viewed this way, the challenge of designing a good grid forming control is a decentral bilinear feedback control problem.If we use the complex voltages v(t) as the states y(t) and the complex frequency as the control input u(t), the underlying bilinear system is reads as: where each system matrix F a has only one entry.The rows and columns of F a are the indices of the nodes in N .The elements of F a are given by: where h, m stands for the nodes in N .This system has the advantage that it is very simple and elegant and is directly formulated in physically meaningful variables.However, the observed power and voltage amplitude mismatches that are required for the feedback are quadratic functions of the state variables, and the asymptotic state that should be achieved is not a fixed point but a limit cycle.The complex coupling formulation solves both of these issues.

Complex Couplings
The power grid dynamics are fully described by the coupling terms K hm .We consider only those lines l whose admittances are non-zero, corresponding to the set of links L (including self-loops) in the underlying power grid.We will treat K hm , K hm , K mh and K mh as separate variables, and the two directions of the edge between h and m ̸ = h as two distinct links l = (h, m) ̸ = (m, h) = l ′ .Thus the state vector y is given by: The nodal complex frequencies u(t) = (η 1 , η 1 , .., η n , η n ) T act as the control input into the system.We define the system matrices R a and Ra to bring the system into a bilinear form 14. We have two controls per node h: η h , η h .The system matrices can be decomposed into the following blocks: We introduce the shorthand notation o(l) and t(l) for the origin and target of a link: Then the origin matrix O a projects onto those links whose origin is a i.e. o(e) = a.The rows and column indices of O a correspond to links of the graph, and we can write explicitly: where l, l ′ are the links in L.
The target matrix T a projects onto those links whose target is h: With these matrices, we have the following bilinear form of the control problem: Taking these equations as the underlying system to be controlled by the complex frequency means that the tracking errors are now just linear combinations of the system state and reference signal.Thus, we have cast the challenge of finding stabilizing grid-forming controllers for the power grid, in terms of a decentral linear feedback control synthesis of a bilinear tracking control problem (14).We should note that this formulation is not without challenges.The fact that we have one coupling per edge, while we have one voltage per node, implies that the complex frequency can not achieve arbitrary couplings from any starting position.Our linear state space decomposes into reachable layers.As we will see further below, this complicates the synthesis of elegant dynamical laws in this formulation.Furthermore, we want to note that the dynamics for the power flow (81) and the voltage square (82) have a self-contained bilinear structure as well.The derivation can be found in the appendix C.1.

A quasi-local, error-minimizing and globally synchronizing grid control
So far, we have formulated a bilinear tracking control problem with linear feedback.The next problem to solve is to design a controller that effectively stabilizes the system.It is natural to consider a quadratic cost function in the error coordinates: This cost function can serve as a Control-Lyapunov function (CLF) for our system.CLFs are an extension of the Lyapunov function from general dynamical systems to systems with control inputs u [31].According to Artstein's theorem [32] if a CLF for a system exists then the system is asymptotically stabilizable, meaning that for any state y a control u(y) can be constructed that asymptotically guides the system back to a fixed point.
Recall that a control-affine system is given by For such systems, a stabilizing control can be explicitly constructed from any suitable V for which the gradient V does not vanish except at the fixed point.The most well-known formula for constructing a stabilizing controller is given by Sontag's formula [32].However, in our context we have F(y) = 0, and a stabilizing control is already given by: where u c is the control input, leading to the dynamics V = −u c uc .
Using the square norm of the error coordinates of the tracking problem 39 in this formula leads to the quasi-local control law: The full derivation of this control law, as well as a more in-depth introduction to CLF, can be found in the appendix D. The control law (40) has the property that: Thus, the mismatch between the realized grid state and the desired grid state can not grow under this dynamical law.As the cost function is positive definite and vanishes at the desired grid state, we immediately obtain that this control law also stabilizes the system's operating state.However, we can not conclude that the system will always find the global minimum of the cost function.
While the cost function is quadratic and thus convex on the full space of couplings, its projection on the reachable layers is not.Therefore, we have to contend with local minima of the cost function.However, these local minima are still globally synchronized because they still imply that V and thus δη go to zero.We show in Appendix D.2 that the local minima of V , subject to the constraint that the h arise from some v, exactly require η c h = 0.It is also noteworthy that this control law seems to stabilize fair power-sharing at the global minimum of the cost function, without coupling the power imbalance to the frequency (see Section 4.4.1).The system always synchronizes exactly at the design frequency.All points that satisfy the power flow equations have V = 0 and thus are stable.This is in direct contrast to the requirements typically imposed on the power flow solution for guaranteeing stability.For example, [17] require that the difference between the phase angles between voltages is smaller, and in particular, not larger than π/2.Thus v i v j has to have a positive real part.Stability of the Kuramoto model is likewise only guaranteed if the effective network Laplacian, weighted with the cosine of the phase angle differences, is positive definite [8].This follows easily from the condition that the differences are smaller than π/2 if the system is lossless.However, the condition is complicated in general.The disadvantage of this controller is that it is not decentral but only quasi-local as η c h requires information about the state of the lines and nodes adjacent to h.Thus, the resulting controller is not in the space of behaviors parametrized by the normal form [11].It is, however, similar in its communication needs to distributed averaging controllers proposed to tackle secondary control objectives in power grids [33,34].It should also be noted that this controller minimizes the cost function, but is not an optimal control in the sense that it minimizes the integral of the cost over time.We consider this controller and the intriguing dynamics it provides as a proof of concept that the formulations given here are highly promising for further analytical investigation.It will be especially worthwhile to investigate how to incorporate a decentrality constraint.Furthermore, it has to be noted that the controller highly depends on the choice of the cost function to use as CLF.Hence, it is valuable to find and analyze other possible CLFs and their corresponding controllers.This is beyond the scope of this work and will be addressed in subsequent papers.

Examples
Using the same algorithm as in section 3.1, we generate a 10-bus synthetic power grid to validate the control law in a simulation example.We implement the complex frequency η c given by the control law (42).We calculate the voltage transients according to equation 7 and then define the currents and powers using Ohm's law and Kirchhoff's law.As an example, we decrease the voltage at a bus to 0.5 [p.u.] and study the CLF V and its derivative dV as shown in figure 3. We can see that the CLF is rapidly decreasing in the first milliseconds of the simulation.Furthermore, we also analyze which of the error coordinates ∆P, ∆Q, ∆ν, gives the largest share to the CLF.From the relative shares, we find that this rapid decrease is due to the reduction of the error in the active and reactive power.After this first initial drop, the CLF decreases further, which corresponds to the reduction in the voltage magnitude errors.Then the CLF saturates at 10 −17 after approximately 4 seconds as also reflected by its derivative which reaches 10 −8 after approximately 3 seconds.Figure 3: The CLF V and its derivative dV after a voltage perturbation to a node.It can be seen that the CLF reduces rapidly in the first milliseconds by compensating the active and reactive power errors, which can be seen in the right-hand figure which shows the relative shares of the CLF to its three components ∆P, ∆Q, ∆ν.After the initial drop, the CLF decreases further and then saturates after approximately 4 seconds.
Additionally, we use the control law for a proper tracking problem.In this example, we simulate a random superposition of sinusoidal fluctuations of the active power set-points at all nodes.Each fluctuation has a different period and amplitude and runs for one respective period.As these are drawn independently the resulting set points are inconsistent during the fluctuation.Note that due to the non-zero derivatives of the set points, we do not have the condition that V < 0. As soon as the variation of the set-points stops, V has to decrease.Nevertheless, the dynamics is constantly working to minimize the cost function and successfully tracks the desired set points. Figure 4 provides the transients of the voltage magnitude and the active and reactive power.We find that the active and reactive powers are adjusted to the new set points.The system does not become unstable, although the set points are highly unbalanced during the transient period.
Remarkably, the voltage magnitudes stay very close to their nominal values, although the system is under severe stress.This is a further success of the designed control law. Figure 5 shows the errors of the rated power and the set points.As the active power set points are highly unbalanced during the fluctuations, it is impossible for the system to perfectly follow the set points.Instead, our control law adjusts the power outputs such that all nodes share the mismatch, which is also referred to as equal sharing.After all fluctuations have vanished, the control law smoothly brings the system back into the previously synchronized and stable state.Further examples, such as the so-called black start capabilities and simulations for larger networks, can be found in the appendix D.1.

Conclusion
The stability of future power grids, particularly in the context of the wide-scale integration of renewable energy sources via power-electronic inverters, remains a challenging topic.While there have been numerous stability analyses of power-electronic inverters connected to an infinite bus bar [6,19], and stability analysis of networks consisting purely of conventional generation [8,10], the stability of inverter-dominated networks is not well understood and is an area of active research.To address these challenges, the normal form [11], a dynamic technology-neutral model for gridforming actors, has recently been introduced.The normal form can describe the dynamics of all grid-forming actors, including grid-forming inverters and conventional synchronous generators, and can thus capture the transitional periods when a reduced number of synchronous generators are still connected to the grid, as well as future highly heterogeneous inverter-dominated systems.In this paper, we have shown that the concept of the complex frequency [12] combined with newly introduced dynamic complex couplings (21), allows for an elegant adaptive network formulation of the Kuramoto model with inertia.When combined with the normal form description of gridforming actors, we obtain a self-contained matrix dynamics for general power grid dynamics in which the grid topology does not explicitly appear.Beyond what we have shown here, the complex coupling dynamics allow for rewriting a large class of dynamical equations as adaptive.As considerable progress has been made recently in understanding the properties of adaptive networks [21,22], this opens up a new avenue for understanding power grids beyond the currently used generator models.We expect that our formulations will be useful in deriving widely applicable analytic stability results in the future.Furthermore, we have seen that the design of stable grid-forming actors can be cast as a bilinear control problem.The control problem simplifies due to the control-affine nature of the bilinear system and allows us to define a stabilizing control input.By defining a suitable Control-Lyapunov function, a control law has been derived that leads to global synchronization.However, it should be noted that the resulting controller is not completely decentralized, as it requires information from adjacent nodes.Thus it resembles distributed averaging control [33] which is usually considered for secondary control objectives, such as restoring frequency and power balance.Both of these are achieved by our controller as well.The formulation of power grids using the combination of complex frequency and normal form shows surprising elegance.This paper has mapped out some, though by no means all, of the relevant results from control theory, electrical engineering, and statistical physics to demonstrate that there is considerable synergy in this approach.Excitingly, the methods include the full physics of balanced three-phase power grids, including losses and generic voltage dynamics.These are known to be important to obtain a full understanding of the collective dynamics of the system [10,24].The proposed models and formulations thus offer a pathway for further research and development in this area, to ensure the stability and reliability of future renewable energy systems.
The right-hand side is well-defined without complexification.It can be evaluated using only the real functions f x and f y and their derivatives.The derivative operators ∂ z/z = ∂ x ∓ j∂ y applied to a complex-valued function on the real domain of x and y is called the Wirtinger derivative.Note that the Cauchy-Riemann equation is simply ∂ z f = 0. Note that we have as might be expected.

A.2 Extrema with Wirtinger derivatives
To illustrate the use of these equations, consider the optimization of a real function V (x, y).The first order condition for extrema is If we write the function in terms of V (z, z) we can write this condition as Take for example V = xy + x 2 + y 2 , then V = zz + 1 4j (z + z)(z − z).The complex derivative conditions then give Which can readily be solved by inserting the second equation into the first.Note that the second equation can be obtained from complex conjugating the first, which can also be verified directly from the definition of the Wirtinger derivatives applied to a real function.Depending on the functions at hand, it can be considerably simpler or harder to work in complex coordinates.

B Line Dynamics
The main body of the paper is concerned with the dynamics of power systems and power lines in the quasi-steady state approximation, that is, the current on the lines is given by equation ( 18): In reality, a change in voltage at one end of the line does not lead to an instantaneous change in current at the other.A dynamical model for lines that takes this into account better is the RL-line model.This describes lines as a series circuit of an inductance and a resistor R. The inductance is typically called L. To avoid confusion with the Laplacian matrix we call it Λ.With complex voltage and current and in a frame rotating with Ω this leads to the differential equation: introducing the complex impedance Z hm = R hm −jΩΛ hm , which is the reciprocal of the admittance Y hm we can write this as For the power on the line, we immediately obtain And S qss hm follows the equation of Section C.1.A remarkable simplification arises if R hm Λ hm is assumed to be homogeneous throughout the network.This is plausible as both R and Λ are proportional to the length of the line.Call their ratio α rl = R Λ , then the dynamics of the power can be given in terms of the K hm directly: Thus, this model of dynamic lines leads to the same adaptive coupling formulation in terms of K hm , the only difference is that the nodal dynamics get augmented by what is essentially a low-pass filter for the part of the nodal power variation that arises due to the coupling.

C Power Flow Variables
The complex power flow S hm on a line connecting node h and m is defined as: The dynamics for the power flow then lead to: If we augment these with the dynamics for the voltage square This again provides a self-contained set of equations for the power flow and voltage amplitude as a function of the complex frequency at the nodes.

C.1 Power Flow Bilinear Structure
For the power flow dynamics, the state vector y is given by: S l = S mh denotes the power flow on the link l as seen from t(l) = m and S l ′ denotes the flow in the opposite direction, as seen from h.To bring the system into a bilinear form, we have to define system matrices.We have two controls per node h, η h , η h , and we denote the corresponding system matrices as N a and Ña respectively.The matrices can then be decomposed into the following block matrix form: where the blocks F a , O a and T a have been defined in sections 4.2-3.1.M a is a mixed block handling the interaction between the voltages at the nodes and the power flow on the links and incorporates the line admittances.The column indices of M a are links and the row indices are vertices.The elements of the mixed block are defined as: The block M a is the element-wise complex conjugate of M a .

D Control Lyapunov Function
A function V (y) is a CLF for a controlled dynamical system ẏ = f (y, u) and the origin 0, if: It is continuously differentiable; It is positive-definite, meaning that V (y) > 0 for all y ̸ = 0 and V (0) = 0; For all y ̸ = 0 there exists a u such that and there exists a u * such that V (0, u * ) = 0. Control Lyapunov functions are particularly useful for control affine systems: The derivative of a potential CLF V for such a system is given by: As in our case a(y) = 0 a stabilizing controller is simply given by: which then results in the following derivative for V : and thus shows that V is a CLF up to the condition that V must not become zero away from the origin.We present the following CLF V for our system: Using equation ( 20) and (82) to get the derivatives of S h and ν h respectively we get the following expression for V : V (96) After switching the indices m and h in the second sum we obtain V in the form of equation ( 88) which allows us to directly define the control: Using the control law (89) we define the control for η c as: then V is given by: which shows that V is a CLF and η c results in global asymptotic stability.For both, the power flow and the coupling dynamics, the state space is larger than the space of physically realizable flows.There is no requirement that a voltage v h exists that can realize the K hm given.However, if we start on the physical manifold, then by construction the dynamics driven by the complex frequency will stay on that manifold.The calculation showing this can be found in section D.2.

D.1 Example Trajectories
To study a more severe control task we also consider a black start of the 10-bus system, which means that we decrease all voltages to 0.01[p.u.].We used this voltage magnitude as a voltage of 0.0[p.u.] is a fixed point of the dynamical system as discussed in section 3.1.From figure 6 we can see that our control can perform a black start for the test system in 35 seconds.However, the black start capabilities are not given for arbitrary systems.Figure 7 shows an attempted black start of a 100-node system generated by the algorithm given in [30].It can be seen that the system voltages are not returning to their nominal values.The test system returns to another fixed point that is not the operating point of the grid.

D.2 Existence of the voltage
For both, the power flow and the coupling dynamics, the state space is larger than the space of physically realizable flows.There is no requirement that a v h exists that can realize the K hm given.
However, if we start on the physical manifold, then by construction the dynamics driven by the complex frequency will stay on that manifold.The calculation showing this can be found in the appendix.To understand where the choice of η c we consider the local minima of V on the physical manifold.
min S,S,ν,λ S ,λ S ,λ ν ,v,v The variation are then given by: Thus with this quasi-local control (the control depends on the power imbalance at the neighbors and the state of the coupling on the edges), V is a Lyapunov function.If the set points satisfy the power flow equations, then the system at the stable power flow has V = 0 and is Lyapunov stable.The control law η c always drives the system to the local minima of V on the physical manifold, and these minima are fixed points as the local minima satisfy η c = 0.

Figure 4 :
Figure 4: Voltage magnitude, active and reactive power transients during the set-point fluctuations of the 10-node test system.

Figure 5 :
Figure 5: Errors of the active and reactive power transients during the set-point fluctuations of the 10-node test system.

Figure 6 :
Figure 6: Voltage magnitude, active and reactive power transients during the black start of the 10 node test system.

Figure 7 :
Figure 7: Attempted black start of a system with 100 nodes.