Stability Analysis of Electrical Microgrids and Their Control Systems

The drive towards renewable energy generation is causing fundamental changes in both the structure and dynamics of power grids. Their topology is becoming increasingly decentralized due to distributed, embedded generation, and the emergence of microgrids. Grid dynamics are being impacted by decreasing inertia, as conventional generators with massive spinning cores are replaced by dc renewable sources. This leads to a risk of destabilization and places an upper limit on the volume of renewable power sources that can be installed. A wide variety of diﬀerent control schemes have been proposed to overcome this problem. Such schemes fall into two broad categories: so-called “grid-following” controllers that seek to match output ac power with grid frequency, and “grid-forming” systems that seek to boost grid stability. The latter frequently work by providing synthetic inertia, enabling dc renewable sources to emulate conventional generators. This paper uses the master stability function methodology to analyze the stability of synchrony in microgrids of arbitrary size and containing arbitrary control systems. This approach provides a powerful and computationally eﬃcient framework in which to benchmark the impact of any number of renewable sources on grid stability and thereby to support microgrid design strategies. The method is demonstrated by computing stability bounds for two diﬀerent grid-forming systems, providing bounds on the feasible number of generators that can be accommodated. In addition, we contrast our results with predictions from a simplistic but widely used phase-oscillator model, ﬁnding that such descriptions signiﬁcantly over estimate the grid stability properties.


I. INTRODUCTION
Many countries are rapidly shifting towards renewable forms of energy production, as part of global efforts to limit pollution and combat climate change.For example, the share of electricity generated by renewables in the UK has risen from 2.8% in the year 2000 to 39.3% in 2021 [1,2].These increases are set to continue in line with global initiatives such as the Paris Agreement [3], with some countries setting legally binding deadlines to reach a net zero-carbon grid within the next few decades [4].
Conventional power stations possess large amounts of rotational inertia owing to the spinning cores in their generators.This rotational inertia has a fundamental role in maintaining power-grid stability by damping out changes in frequency and balancing mismatches in supply and demand [5].However, renewable sources, such as solar and wind, possess no such rotational inertia.Increasing the proportion of renewable power therefore causes a drop in grid inertia, eventually leading to a loss of stability and thereby placing a limit on the volume of renewable sources that can be installed.This creates a critical engineering hurdle standing in the way of a fully decarbonized grid.In addition, renewable power sources tend to be intermittent, introducing fluctuations between renewable and conventional supply that have implications both for network stability and energy-market behavior [6].From a stability point of view, to mitigate the resulting supply-demand difficulties, there is an increasing drive to partition grids into so-called microgrids [7].These systems consist of a relatively small number of power consumers together with embedded renewable generators, connected to the external grid via a point of common coupling (PCC).Theoretically, these microgrids can decouple from the external grid during times of self-sufficiency [8].Modern grids are therefore undergoing drastic changes in both their structure and dynamics, posing significant theoretical and engineering problems that require urgent attention.
The dc electricity from solar or wind generation must go through a converter to output ac power.Presently, most of these converters tend to be managed by control systems of the so-called grid-following type.These systems usually comprise some simple control loops that match the frequency of the output ac power with the frequency of the grid.Unfortunately, this type of control system is unable to mitigate the loss of inertia caused by renewables [9].A host of novel control schemes that aim to boost grid stability have been proposed to replace the grid-following systems.Many of these new control systems aim to simulate the behavior of conventional generators, using batteries to create synthetic inertia.These are referred to as virtual synchronous machines (VSMs) [10] and fall into the broader category of grid-forming controllers [11].These control systems aim to maintain stable grid operation even in the absence of a strong connection to conventional generators.
Assessing the stability properties of these grid-forming systems is of vital importance.Most approaches in the electrical engineering literature rely on direct time simulation, usually of a single converter and control system connected to an infinite external grid.Typically these simulations employ industry standard software, such as Simulink [12] and PowerFactory [13].This timesimulation approach has also been used to model larger systems of converters up to country-scale grids, such as in Ref. [14].While being well suited to prototyping new converter designs and validating their performance under specific scenarios, this approach is computationally onerous: a new time simulation is required for each new grid topology, parameter value, and scenario.In addition, such an approach does not provide systematic insight into system features underpinning its (in)stability.An alternative approach used in the literature is to exploit the ordinary differential equation (ODE) description of a converter system directly, and perform a linear stability analysis (often referred to as "small signal analysis") in order to confirm the stability or otherwise of a given converter and grid architecture.However, as with direct time simulation, this method is particular to the given choice of grid topology and system parameters.Moreover, comprehensive presentation of such ODE systems to describe microgrid power dynamics and their various control systems, and their derivation, is not to our knowledge available in the literature.Numerical continuation was employed in Ref. [15] as an improved method to uncover the role of various parameters on the stability of a small system of converters, but again such an approach is specific to the particular choice of grid size and controller dynamics.Numerical continuation methods are also computationally expensive for larger systems.
More extensive study of power-grid systems has been undertaken in the context of high-power transmission grids.Examples include steady-state investigations of grid resilience to cascading failure [16][17][18], that reveal the importance of the network structure and triggering event on cascade severity; and optimization of grid structure for resilience [19][20][21].Power dynamics in such applications are typically modeled using the so-called swing equation [22].This description falls into the broader category of synchronous network oscillator models, for which there are well-developed numerical and analytical approaches to dynamical analysis; see Refs.[23,24] for an overview.This framework has been in used in Refs.[20,25,26] to demonstrate the impact of transmission network structure on grid synchrony and resilience.
Recent studies have adopted these network modeling approaches to understand how the proliferation of renewable power generation affects grid stability, a key feature of which being the variable and distributed nature of power generation.For example, employing steady-state dc approximation [27,28] and the swing equation in tandem with generation and usage data [29], the role of changing grid composition (associated with variations in daily and seasonal usage, and meteorological conditions) in significantly decreasing grid resilience was exposed.The impact of variability in wind-power generation on grid stability was studied in Ref. [30].However, such approaches either do not describe accurately low voltage distribution grids and microgrids, do not account for the dynamics of the converter systems necessary for renewable distributed power sources, and/or invoke network architectures at national (high voltage), rather than microgrid scale.
This paper first provides a comprehensive derivation of the dynamical system appropriate to describe the operation of microgrids of arbitrary size and under a given control system.A semianalytical framework to determine stability boundaries of such a system is then presented.This method is computationally efficient and provides stability bounds in terms of network size and power demand under variation of key control parameters, without the need for repeated simulations for different system sizes.Moreover, this approach is flexible with respect to converter control design, allowing different control prototypes to be rapidly benchmarked against each other.We demonstrate the efficacy of this method for grids under two different grid-forming control systems, and contrast these predictions with a widely used phase-oscillator model.

A. Basic model of converters
This section introduces the mathematical model of the inverter-filter system used to convert dc output power from distributed generators (DGs) into usable ac power for the grid.This inverter-filter system will form the base model for each DG.Subsequent sections will show how grid dynamics and control system dynamics can be built on top of this base model in such a way that allows for convenient stability analysis of the whole network.
Figure 1 shows a diagram of a dc source coupled to an inverter and filter.The inverter consists of banks of semiconductors, which effectively chop the constant dc input signal up into a noisy square-wave approximation of a three-phase signal with frequency ω.In other words, it produces three channels of approximate sinusoids each shifted by 2π/3 relative to each other, each with a current amplitude i 0 , and voltage amplitude v 0 .These values, together with ω, are in general set by some control system, which we leave unspecified for now.These approximate ac signals comprise a spectrum of frequencies, including high-frequency noise, in addition to the desired frequency ω.To remove these additional modes, the signals are sent through a low-pass filter.These filters typically comprise inductors with inductance L and capacitors with capacitance C arranged in series and parallel respectively, as shown in Fig. 1.The values of L and C are chosen to provide maximal attenuation of high-frequency modes, but low attenuation of the desired mode ω.Smooth three-phase power with frequency ω, voltage v, and current i then flows out from the filter into the external network.
In order to build up the mathematical model of the inverter-filter system shown in Fig. 1, we first note that each of the three-phase quantities are described by timevarying three vectors.For example, the output voltage is , where a, b, and c refer to each of the three channels.The internal behavior of the inverter will not be modeled because the timescale of its semiconductor dynamics, referred to as the switching time, is orders of magnitude shorter than the timescale of grid dynamics.For instance, the switching time of the commonly used IGBT inverter is typically on the order of nanoseconds to microseconds, while the shortest grid dynamics typically occur over a few milliseconds [31,32].The remainder of the system can be modeled with six linear ODEs: three describing the filter inductors in each phase, and three for the filter capacitors in each phase; these are omitted here for brevity but recorded in Appendix A for completeness.For convenience, the system is transformed into a rotating reference frame with frequency ω g , where ω g is the nominal grid frequency (typically 2π × 50 Hz for Europe, Africa, and Asia, and 2π × 60 Hz for the Americas).This standard transformation is referred to as the dq transform, and allows the sinusoidal three-vector quantities to be expressed as two vectors with components d and q.Performing the dq transform on the six filter ODEs yields the following four linear ODEs (full details are given in Appendix A) determining the currents (i 0,d , i 0,q ) and voltages (v d , v q ): where the d and q superscripts denote the d and q channels, and time dependence of currents and voltages has been suppressed for notational clarity.The output current (i d , i q ) will be determined by the external grid dynamics, while the output frequency ω will be specified by a control system.Equations ( 1)-( 5) define the base model of DG dynamics.
Sections II B and II C will examine how grid and control dynamics, respectively, can be incorporated into this base model to study grid stability.

B. Microgrids and dynamically equivalent networks
Each DG, as described by Eqs. ( 1)- (5), must be coupled to some external loads via a network of connections.This section describes this network model for the case of microgrids: small scale collections of DGs and loads coupled to an external transmission grid via a PCC.
Various microgrid models exist [33,34]; in the following, we concentrate on a treelike architecture as illustrated in Fig. 2 both as it is the most commonly employed and since it provides the most extreme case from the point of view of network instability.Each DG is connected to a local load, which could be the power demand of a household, for instance.This load is represented as a complex impedance.The DG and local load are then connected upstream to the PCC via a transmission connection with impedance Z T , which may again be complex.The PCC then couples to the external grid, which is represented as another potentially complex impedance Z grid .Kirchoff's laws of conservation can then be used to model the system, resulting in coupled sets of ODEs for each impedance 2. Schematic of a microgrid.Distributed generators, consisting of a source and filter as shown in Fig. 1, are coupled into a network.There may be a local load, perhaps representing household demand, and a transmission impedance Z T .All houses and DG units are then connected via a PCC to the external network.
The external network is modeled as an impedance Z grid .
in the network.However, this results in a potentially highdimensional set of equations wherein the role of network topology and coupling is ambiguous.To recast this system in a form more amenable to analysis, and thereby to exploit some of the methods of network science and dynamical systems, a natural approach is to adopt a description as a network of nodes, each of identical type, whose coupling can be described with an adjacency or network Laplacian matrix.For high-voltage long-range transmission grids, this separation is typically achieved by assuming that the voltage and frequency are constant across the grid and that the resistance of transmission lines is negligible.This leads to the swing-equation network formulation commonly used in the complex systems literature, such as in Ref. [26].Unfortunately, these assumptions do not, in general, hold for lower voltage microgrids.In the following, we therefore develop the following method for analyzing microgrids governed by Eqs.(1)-( 5), together with a control system to be described in Sec.II C, without approximation as complex networks.The resulting formulation, separating as it does the nodal dynamics from network structure, provides a convenient framework for the interrogation of different network structures or control systems.
First, consider a microgrid of n DGs and local loads, indexed by k = 1, . . ., n, arranged as in Fig. 2. The local loads can each be described by two linear ODEs, representing inductances in both the d and q channels, and absorbed into the description of each DG by appending them to the set of Eqs. ( 1)- (5).For simplicity, let us assume that local loads are zero, so that all power must flow to the PCC; this limiting case represents the maximal case of supply-demand mismatch and hence our stability analysis that follows provides the least upper bound on microgrid size.The d and q currents on the kth transmission edge are given by where L T and R T are the edge inductance and resistance, respectively, and v d PCC and v q PCC are the d and q channels of the voltage at the PCC.The edge currents, i d k and i q k , are the output currents for the kth DG as featured in Eqs.(3) and (4).Equations ( 6) and (7) are subsumed into the dynamics for the kth node, so that the dynamics of edge k are appended to the dynamics of DG k.This is illustrated in Figs.3(a Illustration of the process of converting the microgrid model, shown in Fig. 2 into a dynamically equivalent network that is more amenable to analysis.In (a), all DG units are described by (1)- (5).In (b), the additional ODE governing the transmission line currents on the incident edges has been absorbed into the nodes' dynamics.In (c), the equations for the PCC load have also been absorbed, resulting in effective all-to-all coupling between nodes.
The PCC voltage v PCC in these equations depends upon the external grid load.In the case of a purely Ohmic grid load, then by Kirchoff's law, this voltage is proportional to the sum of the input currents from each of the k nodes: where R d is the grid resistance.This term is then incorporated into the equations in each of the k nodes, so that Eqs. ( 6) and ( 7) become where A is the adjacency matrix of a complete graph on n nodes.Effectively all nodes are now coupled directly to each other via the PCC, meaning that the microgrid is dynamically equivalent to a complete graph with the dynamics of each described by Eqs. ( 1)-( 5), (9), and (10).
The coupling between nodes is linear and occurs through the i d k , i q k variables.This process of converting a microgrid to a dynamically equivalent complete graph is illustrated in Fig. 3.We reiterate that this methodology preserves all physical quantities; in other words, the voltages and currents at every point in the system are modeled.In the above presentation we have assumed the PCC load impedance to be purely resistive; however, as shown in Appendix F, the same process can be conducted for the case of a complex load leading to an equivalent description albeit with more complicated terms.A full stability analysis for this case included therein demonstrates that the inclusion of the reactive load component produces only a very small divergence from the results of the purely resistive case: a substantial inductive load of 50 mH causes only 3% shift in the location of the bifurcation point.Therefore, for pedagogical purposes, we shall stick with the purely resistive case in what follows.
Denoting the state vector by 6 for the kth node, the dynamics for this base model, given by Eqs. ( 1)-( 5), (9), and (10), can be compactly written as A kl H (x l ), (11) for each node k.Here, σ = R d /L T is the coupling strength; f : R m → R m describes the intrinsic, uncoupled dynamics of the kth node; A is the adjacency matrix; and H : R m → R m is the coupling function.In this case H is given by 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 since coupling occurs linearly though only the i d k and i q k components of x k as described above.In what follows, we will exploit the master stability methodology as a convenient and powerful framework with which to analyze this system.To do so, the coupling matrix in Eq. ( 11) must have rows that sum to zero [35].This is not true of the adjacency matrix A introduced above.It is, however, true of the Laplacian matrix.We can therefore exploit the relationship where L is the Laplacian matrix and D the degree matrix, whose kth diagonal element is the degree d k of node k to rewrite Eq. (11) as where In the case of the complete graph, this is simply The network model, described by Eqs. ( 1)-( 5), (9), and (10), for the microgrid is now encapuslated by Eq. ( 14), in which the function F : R m → R m for this model is defined by Any additional control systems, VSM or otherwise, are accommodated by suitable addition to this specification of the intrinsic nodal dynamics F(x k ), as will be demonstrated in Sec.II C.

C. Control systems
The previous section introduced a model formulation in Eq. ( 14) for distributed generators in which the contributions from internal nodal dynamics F and interaction terms H are separated.The derivation presented assumed no control system governing input voltage and frequency.In this section we show how an arbitrary control systems for each generator can be incorporated while retaining this convenient model structure; in other words, in such a way that changes only F, leaving the interaction terms unchanged.
There are many different control schemes for distributed generators.Our objective here is not exhaustively to survey them, but to demonstrate how any scheme can be incorporated into the formalism of Eq. ( 14).Subsequent sections will then demonstrate the use of the master stability function to benchmark and compare different control systems.We will focus here on two different grid-forming control schemes: a VSM and a simple voltage-frequency droop controller.The VSM simulates the dynamics of a conventional rotating generator, with the aim of providing synthetic inertia to stabilize the grid.Using a model based on Refs.[14,36], the VSM control system can be described by the following seven ODEs: dε dγ Equations ( 16) and ( 17) control the output frequency ω k by emulating the swing dynamics with inertia J and damping D p (see Sec. III B for more details on the swing equation and a comparison with the microgrid formulation), while Eqs.( 18)-( 22) represent proportional and integral current control.The input voltages in Eqs. ( 1) and (2) are then computed via where i 0,d, * k and i 0,d, * k are reference currents given by Equations ( 1)-( 4), together with Eqs. ( 9), (10), and ( 16)-( 22) now describe the dynamics of the kth node as a set of 13 ODEs [and the algebraic relations Eqs. ( 23)-( 26); note that Eq. ( 5) no longer applies].A microgrid containing n DGs controlled by VSMs can then be described as an n × 13 set of ODEs, with the dynamics of each node written in the form where x k ∈ R 13 .Refer to Appendix E 1 for further details on the derivation of the VSM equations and for the specific structure of F and H above. Appendix E 2 provides a derivation for an alternative voltage and frequency droop controller.

D. Microgrid stability analysis
The dynamics of an n-node microgrid is now described by the n × m system given in Eq. ( 27); in the case of the VSM model developed above, m = 13.By utilizing the so-called master stability function (MSF) method, first developed by Carroll and Pecora [35], the stability of the synchronous oscillatory network state can be determined using a system of only m equations.This effectively allows the dynamics of the whole network to be interrogated using the equations of only a single node.This subsection gives a brief overview of the MSF method and its application to the microgrid model Eq.(27).
Let y(t) ∈ R m be a time-periodic (with period T) solution to the uncoupled dynamics of Eq. ( 27), such that ẏ = F(y).Linearizing the microgrid equations via 1 and Taylor expanding provides the variational equation L kj z j , for k = 1, . . ., n, (28) where DF(y) and DH (y) are the Jacobian matrices of F and H , respectively, both evaluated along the synchronous solution y(t).Equation ( 28) can then be simplified by projecting z into the eigenspace of the Laplacian matrix L, yielding (29) where γ k = σ λ k and λ k is the kth eigenvalue of L. The block structure of Eq. ( 29) leads us to consider the generic The MSF is then defined as the function mapping γ to the Floquet exponent of Eq. ( 30) with the largest real part.
The synchronous solution y(t) is then stable if the MSF is negative for all γ = σ λ k , from which we hence obtain the critical point of instability γ c .Recall from Sec. II B that the microgrid is dynamically equivalent to a complete graph on n nodes whose Laplacian hence has one eigenvalue equal to zero and all others equal to n.Therefore, we need only to consider the case λ k = n since the zero eigenvalue simply corresponds to perturbations along the synchronous manifold and can be ignored for the purposes of a stability analysis.The coupling strength is σ = R d /L T and so The critical point of instability γ c therefore occurs at α c /L T , assuming the transmission impedences L T remain constant.Once this critical point has been found, the maximum number of nodes that can be supported on the network as a function of load impedance is given by Evaluating the MSF for the system in Eq. ( 30) and finding γ c therefore gives a stability bound in terms of network size and network load, without having to do many different experiments for increasing n.In other words, the computational complexity is constant in n.

III. RESULTS
To demonstrate the dynamics of the VSM model Eq. ( 27), and to motivate the results that follow, Fig. 4 shows a direct simulation of a test microgrid system containing n = 2 DGs.The parameters used are summarized in Table I.The results in Fig. 4(a) show the frequencies of each DG.A perturbation to the frequency of DG 1 is added at time t 0 = 0.5 s.This perturbation can be seen in the discontinuous jump in the frequency ω 1 of node 1.This increase rapidly decays over approximately 0.1 s back to the original frequency of 2π × 50 Hz; the frequency of the second node is not visibly affected.Panel (b) shows the resulting three-phase voltage output of node 1.The perturbation can be seen at t 0 = 0.5 s.The voltage quickly returns to regular output.This system response reflects a stable operating regime, in which the network quickly returns to synchrony on experiencing a disturbance.
Figures 5 and 6 summarize the stability properties of microgrids of DGs as a function of key model parameters (coupling strength α, grid resistance R d , and synthetic inertia J ) and network size n.Default parameter values are listed in Table I.
The maximum Floquet exponent μ max computed as a function of α = γ L T for a microgrid under VSM control is shown in Fig. 5(a), indicating the critical coupling strength α c = γ c L T = 5778 below which a given network is in a regime of stable operation.Recall that this computation is independent of network size, with subsequent inferences of permissible network size and load obtained via Eq.(32).Panels (b) and (c) highlight the dependence of this critical threshold on synthetic inertia J , this being a natural control parameter in the VSM system.In (b) the locus of the stability threshold is presented in the (n, R d ) plane for two TABLE I. Default parameters for the VSM model defined in Eq. ( 27), or equivalently: Eqs. ( 1)-( 4), ( 9), (10), and ( 16)- (22) 27), or equivalently Eqs. ( 1)-( 4), ( 9), (10), and ( 16)- (22).Parameters are given in Table I. different values of J .The region of stable operation lies below the curve, and so these results highlight an inverse relationship between feasible microgrid size n and load resistance R d .As expected, the window of stable grid operation is decreased with decreases in inertia J .We remark that for values of load resistance in (on average) physically reasonable ranges (R d ∼ 60 , as listed in Table I), the inverse relationship between n and R d facilitates strong control over feasible grid size via the grid inertia [see also Fig. 6(a)].An alternative interpretation is that for fixed J , in this load-resistance regime, small deviations in R d can have dramatic consequences for grid stability.As the load resistance increases to larger values, control over this stable window requires higher inertia.Figure 6 compares grid-stability predictions obtained under VSM and for a voltage frequency droop controller, as presented in Ref. [37] and treated widely in the literature (see Appendix E 2 for more details), and highlights the consequences of these grid-stability results in terms of power demand P. Panel (a) illustrates the significantly larger stable region in the (n, R d ) plane accessible under VSM control [results here are shown for J = 1, D p = 50, corresponding to the green line in Fig. 5(b)].Figure 6(b) recasts this in terms of power demand, this being computed from the (n, R d ) values in panel (a) via the method presented in Appendix C. Microgrids under VSM control that lie above the green line in Fig. 6(b) [which corresponds to the green stable boundary identified in panel (a)] exhibit destabilization of the synchronous state by becoming too large for the available power supply.The orange curve is the equivalent stability boundary for microgrids under droop control.The additional dashed line demarks a 6 A (occurring at R d = 36.9 ) safety cutoff beyond which the associated currents become too large for feasible operation.This illustrative threshold is chosen to correspond to suitable domestic ratings; however, alternative such thresholds are straightforward to calculate (see Appendix C for detail).The shaded region hence describes grid regimes for each control system which are both stable (according to previously described analyses) and within operating tolerances, further highlighting the significantly more restricted regime of operation accessible for grids under droop control.
Direct simulation of the full set of m × n equations representing given microgrid configurations under VSM control highlighted by points (i) and (ii) in Fig. 6(b) confirm our stability predictions (see Appendix H for numerical details).Figure 6(c) provides schematic representations of these microgrids, which comprise n = 6 VSMs operating with power demand within the stable region [marked by point (i) in panel (b); P = 1297 W, for which R d = 60] and on the boundary [point (ii); P = 85.4W; R d = 952.2],with node colors indicating the voltage frequency produced by each node's control system after 0.25 s and emphasizing that the former exhibits stable synchrony at ω = 2π × 50 Hz while the latter has desynchronized.The destructive effect of this destabilization is illustrated in panel (d), which shows an example timeseries of the three-phase voltage recorded from a VSM within our six-node microgrid that experiences a drop in power demand from point (i) to (ii) in panel (b) at t = 0.125 s.

A. Microgrid heterogeneity
Previous results correspond to constant power demand; however, within microgrid operation significant variation over short timescales will naturally occur and as can be observed from Fig. 6(b), depending on microgrid size, such deviations need not be large to place the system in an unstable regime.Moreover, particular network architectures may deviate from the treelike organization considered so far.Figure 7 provides exemplar results for microgrid networks under VSM control indicating the consequences for stability of considering some simple forms of heterogeneity in network structure and power demand.
Panels (a)-(c) compare the stability threshold value of load resistance R d in three small networks of varying topology with n = 4 identical nodes.The ODE structure of these nontree networks is presented in Appendix G.The results highlight that the tree structure provides the strictest stability constraint, having the lowest critical value of R d .In panel (d), the impact of variability in power demand on In each case the local load of each node is drawn from a normal distribution with variance ν, and mean R dl .Each system is then time stepped to determine if a stable state can be found.The results are then plotted in the (P, ν) plane, where P is the power delivered to the grid and is calculated from R dl using the equations in Appendix C. The shaded region indicates the stable regime for the homogeneous system, computed using the MSF (Sec.II D).
the stability of networks of n = 6 nodes with tree structure is evaluated via direct simulation.In each network, nodes are allocated a (constant) local load R dl from a normal distribution with given mean R dl and variance ν.The grid load is set at R d and acts in parallel with the local load, as in Fig. 2. The networks are then time stepped to determine if a stable state can be reached.The results are then plotted in the (P, ν) plane, where P is the mean power each network delivers to the grid, and is computed from R d and R dl using the equations in Appendix C. The values R d and R dl are chosen so as to place all of the networks close to the stability boundary calculated for the homogenous case using the MSF.We see that even in the presence of significant fluctuation in local demand, networks lying in the stable region previously identified for homogeneous nodes are able to function stably; indeed, some stable configurations are observed even for node demand that places the network outside of this region.In summary, the results presented in Figs.5-7 highlight how the MSF method presented in Sec.II allows analysis of microgrids under different forms of control and with different network architectures without resorting to direct simulation.Thereby, efficient interrogation of grid function under variation of key control parameters (such as synthetic inertia, in the case of VSM control) can be undertaken to determine bounds on stability that dictate safe operating regimes.The stability analysis considers networks of identical nodes, and provides a form of extreme bound for network stability; however, significant deviations in power demand (or equivalently R d ) will naturally occur within such networks.Exemplar simulation results presented in Fig. 7(d) demonstrate that these stability bounds remain informative even under introduction of significant noise.Together, these results provide an indication of how the MSF method could support the design of dynamic control schemes incorporating manipulation of VSM damping or inertia, or battery-charging schemes.

B. Limitations of phase-oscillator microgrid models
As highlighted in Sec.I, the so-called swing equation has been widely employed in studies of power-grid dynamics, stability, and resilience, including more recently in studies of smaller-scale distributed generation.However, despite its historical importance such a model does not describe accurately the low-voltage and converter dynamics of key importance to microgrid operation of interest here.In this section, we contrast the predictions obtained from employing such a model with those obtained from the detailed microgrid description described above.
For clarity, we first state the swing-equation model equations, but refer the reader to [29,[38][39][40] for a full derivation, discussion, and historical context.In brief, the swing equation describes each generator in a grid of size n as a rotating machine, balancing inertial, dissipative, and transmitted electrical power.The deviation θ k (t) (k = 1, . . ., n) of the phase angle of each from the grid reference frequency ω g = 2π × 50 Hz is governed by in which γ k = 2D k /J k describes the balance between effective damping D k and inertia J k at each node k.Denoting the power demand by P m k , P k describes the effective power on each node defined as K kl is the coupling strength between nodes k and l, given by K kl = V 2 0 /(X kl J k ), where V 0 is the common grid voltage and X kl = 2π × 50 × L is the reactance of the edge connecting k and l, with L being the edge inductance.We emphasize that the voltage V 0 is held constant in the swing-equation description.Self-evidently, therefore, voltage dynamics are not captured in the swing-equation description.Another simplification inherent in this description is the assumption of reactive-only edge impedances.
Figure 8 shows direct simulation and numerical continuation results corresponding to a microgrid of n = 6 generators, governed by Eq. ( 33); numerical details are given in Appendix H.The grid voltage is chosen as V 0 = 240 V and edge inductance L = 0.02 H; values for inertia and damping parameters are chosen to reflect the low-inertia nature of the renewable power sources of interest here: J k = D k = 0.01.Panel (a) shows the frequency response of each generator when power demand for each node is chosen as P k =: P = 1 W, at which the results shown in Fig. 5(d) for the microgrid converter model indicate instability.In contrast, however, we observe relaxation to a synchronous state within a few seconds, indicating stable operation; this is confirmed in panel (b) where results from numerical continuation of Eq. ( 33) indicates stable operation for all power demands between 0 and approximately equal to 30 kW.
These results highlight clearly that the swing-equation model significantly overestimates the stability properties of microgrid operation, in comparison to those obtained from realistic descriptions of relevant power-converter systems, as illustrated in Figs.4-6, and Appendices E-G.

IV. DISCUSSION
This paper has provided a framework to analyze the stability characteristics of electrical microgrids, a theoretical and engineering problem of increasing importance, as the drive towards decentralized and renewable power generation continues.Typically in such grids, dc electricity from, e.g., solar or wind is converted to output ac power via a converter, managed by a control system that aims to match output ac frequency with the grid.Such "grid-following" control systems unfortunately have negative consequences for stable grid operation.Alternative "grid-forming" control systems attempt to boost stability, for example, by simulating the behavior of traditional generators via the creation of synthetic inertia.
Understanding the destabilizing effects of renewable incorporation into the grid is crucial to determine potential upper limits of these power sources that can safely be accommodated, and to design new control strategies or grid specifications for the amelioration of such effects.The majority of existing analyses either (i) follow extensive previous work in the context of high-power transmission grids, and thereby concern models that do not describe accurately low-voltage distribution grids and the dynamics of the converter systems necessary for renewable distributed power sources, and/or invoke inappropriate network architectures, or (ii) employ direct time simulation to study the dynamics of specific converter and/or grid systems.
Here, we address these issues to provide a semianalytical framework for the analysis of microgrid stability.Crucially, the method is flexible to the incorporation of different control systems and network structures, and for the study of networks of arbitrary size, without the need for comprehensive time simulation as is typically employed in industry-standard approaches.We illustrate the utility of our approach by computing stability bounds on the volume of renewable generation that can be accommodated in a microgrid for two different grid-forming control systems.Specifically, we contrast results for VSM and droop control, demonstrating how power demand, grid size, and control-system parameters (in particular, we consider synthetic inertia and damping in the case of VSM control) combine to determine regimes of feasible operation.We additionally contrast our predictions with those obtained from an archetypal phase-oscillator model, whose relative simplicity makes it amenable to analysis and is hence very widely studied.While such a description is appropriate for high-power transmission grids, we show that it significantly overestimates the stability properties of microgrids and is hence of limited utility for studies aimed at aiding microgrid design or control.The method demonstrated herein shows that the same powerful insights on grid stability, as are straightforwardly obtained from such simplified models, can be derived for faithful microgrid controller descriptions with relatively little increase in mathematical or computational cost.
Our analysis provides a form of extreme bound on stability, considering, in the main, a treelike network connectivity architecture and assuming as it does that local loads are absent so that the supply to the grid is maximal.Realistic microgrid dynamics will feature heterogeneity in both generation and usage, and will employ a range of network designs.Comprehensive study of the former necessitates, for example, incorporation of a suitable model for local (household) loads into our dynamical model; while the MSF extension to weak heterogeneity caused by, e.g., PV shading or small variations in radiation intensity is relatively straightforward [41], methods for treating such strongly heterogeneous networks are not available.Such a study necessitates comprehensive direct simulation and so lies outside the scope of the current work whose focus lies in the exposition of the MSF method for efficient benchmarking of grid-control strategies.Nevertheless, results considering alternative networks and the introduction of a simple model of noise indicate that stability bounds provided by the MSF method remain informative.We remark in passing that in line with our focus on electrical grid operation, the analysis here concerns stability of the synchronous state; the MSF method presented is flexible to consider other more exotic states and we direct the interested reader to Ref. [42] and references therein.
Lastly, we highlight that as part of this work, we provide a comprehensive derivation of the dynamical system suitable to describe microgrid dynamics in a form convenient for simulation and analysis, thereby providing a valuable resource for future mathematical studies.This formulation will enable complex network-based analyses, such as those in Refs.[16,19,20], to be extended to incorporate higher fidelity low-voltage dynamics and control systems.For example, natural and important future work includes comprehensive consideration of the heterogeneous and dynamic features of generation or usage described above (in the spirit of Ref. [29], for example), and their implications for usage-specific control strategies.The current and voltage over the inductor and capacitor is given by The so-called dq transform recasts this system in a rotating reference frame; for a three-phase signal x = (x a , x b , x c ), this is defined or more compactly where θ = ω g t, with ω g being a reference frequency, say 2π × 50 Hz.The resulting dq signal is x = (x d , x q , x null ).In all cases where the network is balanced, the x null component is 0, meaning only the d and q components need to be considered.An example of this transform is shown in Fig. 9.In practice, networks will have small deviations from the nominal frequency, resulting in oscillatory behavior in this dq frame.Applying Eq. (A4) to Eqs. (A1) and (A2) yields the base model given in Sec.II A. Considering Eq. (A2), for concreteness, we obtain It can be shown that and hence we obtain

APPENDIX B: PHASOR NOTATION
Phasors are a commonly used notational shorthand in electrical engineering to represent sinusoidal signals.They allow the steady-state amplitude and phase shift of ac circuits to be found in a similar vein to Ohms law for dc circuits.Phasors are used in Appendix C to compute the steady states and power-flow levels of the inverter, and we include a brief explanation of this notation here for completeness.
In essence, phasors represent time-varying currents in the form of complex exponentials.For example, consider a time-varying voltage V(t) = V 0 e j (ωt+δ) , where j is the imaginary unit (to avoid ambiguity with currents) and δ is a phase shift.The phasor of V is V := V 0 e j δ ; in other words, V(t) = V e j ωt .It is common to employ the notation V = V 0 δ, where the • notation is shorthand for exp(j •).
Recovering the real sinusoidal signal from V requires multiplying the phasor by the time-varying component and taking the real part: The utility of the phasor notation is that it allows for reactive components to be conveniently added in series and parallel in the same way as a resistive element.10.Single-phase model of a decoupled distributed generator with no control system.Input voltage and current i 0 and v 0 flow through a low-pass filter.A complex load with impedance Z grid together with a transmission-line inductance L T model the grid and transmission line.
For clarity of presentation, in the following example calculation we compute the steady state of a single-phase system, shown in Fig. 10.Power values for three phase power can then simply be obtained by multiplication by a factor of 3.
The single-phase system is described by Eqs.(A1) and (A2) (in single-phase form) and where i and L T are the current and inductance of the load, respectively, and R d is the grid load resistance.We note in passing that the characteristic polynomial of the Jacobian of Eqs.(A1)-(C1) satisfies the Routh-Hurwitz criterion, ensuring all three eigenvalues have negative real part and hence the system is stable for any physical value of L, L T , C, and R d .To solve this system, we employ phasor notation described in Appendix B, above.Considering Eq. (A2), we obtain where V 0 , V, and I denote the phasors of v 0 , v, and i.The load current phasor is defined by Combining Eqs.(C2) and (C3) yields the phasor of the output voltage as or, equivalently where z ∈ C denotes the denominator of Eq. (C4).
The amplitude I of the current I flowing through the inverter can now be found from V by dividing through by the impedance value: The complex power S of a three-phase system can also be found from V: wherein the factor of 3 comes from the fact we want three-phase power and the 1/2 arises from employing RMS values.The real and reactive power are then P = Re(S) and Q = Im(S), respectively.Figure 11 shows the power, voltage, and current as a function of R d computed using the equations above.All parameter values are as given in Table I.For low R d , there is a power spike, as shown in (a).Correspondingly, (b) and (c) show a collapse in voltage and a spike in current for these low values of R d .This is dangerous short-circuiting behavior and would not be allowed in reality.The dashed line in each panel (and in the corresponding panels in Fig. 5 in the main text) indicates an illustrative safety cutoff of 6A, which occurs at R d = 36.9, computed via Eq.(C6).This threshold is chosen to be within the range of fuse ratings for typical domestic premises, but alternative cutoff points would be straightforward to calculate by solving Eq. (C6) for R d .
For the case of an inverter connected to both a grid load and a local load R dl in parallel, the total power delivered can be calculated by substituting R dl R d /(R dl + R d ) for R d in Eqs.(C4) and (C7).To compute only the power delivered to the grid when a local load is present, one can use Eq.(C7) with a modified z.This modified z can be obtained by considering the system illustrated in Fig. 10 but with a load R dl added in parallel.After some algebraic manipulation, one obtains the real and imaginary components of this modified z as and

APPENDIX D: COMPUTATION OF FLOQUET EXPONENTS
In order to compute the master stability function and determine the stability of the microgrid, we must compute the Floquet exponents μ i of the linearized Eq. (30).That is, we need the Floquet exponents of the m-dimensional where the matrix J ∈ R m×m is defined as in Eq. (30) as J (y) = DF(y) + γ DH (y), and y(t) is the synchronous orbit whose dynamics obey ẏ = F(y).Now let M (t) ∈ R m×m be a fundamental matrix solution of Eq. (D1), so that Ṁ = J (y)M , (D2) and let the synchronous orbit have period T. The Floquet multipliers ρ i = exp(μ i T) are defined as the eigenvalues of the monodromy matrix B = M −1 (0)M (T).Setting M (0) equal to the identity matrix I m×m then gives the monodromy matrix as B = M (T).If T is known, then B can be found numerically by timestepping the (m 2 + m)-dimensional system ẏ = F(y), (D3) from t = 0 to t = T. Once M (T) and its eigenvalues ρ i are found, the Floquet exponents of the system can be computed as In the case of the microgrid model of interest, the period T of any synchronous state is not known a priori.To deal with this, the timestepping procedure described above is augmented with a QR decomposition-based reorthogonalization scheme, as described in Ref. [44].This entails timestepping the system for some time t max much larger than any expected T.Then, at regular intervals of t the current state M (t) is decomposed with a QR decomposition.The Floquet multipliers are the diagonal elements of R. The timestepping is then resumed with Q being set as the new initial state of the matrix M .This continues until t max is reached.This process is summarized in Algorithm 1.The computation also requires an initial point p on the stable attractor of a single node.This is found by timestepping Eq. (D3) for 100 s with an arbitrary initial condition.
There is no rigorous way to choose t max and t.In general, the larger t max and the smaller t, the more accurate will be the final result.Here, this balance between accuracy and computational cost is chosen as t = 0.1 s and t max = 100 s by trial and error.

APPENDIX E: DERIVATION OF DG CONTROL SYSTEM GOVERNING EQUATIONS
In this Appendix we demonstrate how the equations governing DG control systems are derived and appended to the base microgrid ODE system, thereby to provide a convenient and flexible description of microgrid dynamics for simulation and analysis.We re-emphasize that such a complete derivation is not available as standard in the electrical engineering literature, the concentration there being more on direct simulation in specific industry-standard software, and so this section potentially provides a valuable resource for future analyses of microgrid dynamical properties.

VSM control system
The underlying idea of the VSM control method is to allow the frequency ω k and the voltage v 0 k of the inverter to Algorithm 1. Compute Floquet exponents.
react dynamically to grid dynamics in a way that emulates the classical swing dynamics of a conventional generator.This is achieved by calculating the frequency according to where P is the measured instantaneous real power of the inverter, ω g is the nominal grid frequency, and P ref is the reference real power (see Appendix C).J and D p are the virtual inertia and damping coefficients, respectively.The phase angle of the inverter then evolves according to The dynamical behavior of the voltage is emulated by computing a virtual version of the classical droop dynamics [45] according to where E k is a reference voltage to be fed to the inverter, K and D q are virtual droop parameters, V is the magnitude of the measured voltage that is output by the inverter system, and Q is the measured instantaneous reactive power.Q ref and V ref are the set reference reactive power and voltage magnitude, respectively (see Appendix C).The measured power and voltage quantities are calculated via where i d/q k and v d/q k are the d/q components of the current and voltage at the output of the inverter system, as measured in the ω k dq-reference frame.
Once the dynamical reference voltage E k and the frequency ω k have been computed, they are fed through the current and voltage loop system illustrated in Fig. 13.These systems are common to many control systems and consist of a rack of four interlocked proportional-integral (PI) controllers; see Fig. 12.Each PI controller takes some input signal x(t) and computes an output y(t) according to The output signal is thereby a sum of the reference signal and its integral, weighted by proportional and integral coefficients K p and K i , respectively.The presence of four PI units in the control system results in four additional ODEs.
The current and voltage loop system (blue area of Fig. 13) then proceeds as follows.First, a so-called outer reference voltage for the d, q components (as indicated by the d/q shorthand) is computed: 12. Schematic of a PI controller.An output signal y(t) is computed as the sum of an input signal x(t) and its time integral, weighted by coefficients K p and K i , respectively.Four PI controllers are built into the current and voltage control loops, as illustrated in Fig. 13.  13.Schematic of the full control scheme.This system provides a reference signal that is fed to the igbt inverter illustrated in Fig. 1.The inverter then outputs three-phase sinusoidal voltage and current slaved to this signal.Using measurements from the inverter and low-pass filter (see also Fig. 1), the control scheme dynamically varies this reference signal, emulating the classical swing and droop dynamics of a conventional generator, in order to keep the inverter's voltage and current output constant and stable.The system works by first computing a phase angle φ k according to the virtual swing equation given in Eq. (E1).This component is highlighted in the orange block.A reference voltage amplitude E k is then computed according to the virtual droop equation given in Eq. (E3).The dq components of this voltage are then computed according to E k cos φ k and E k sin φ k , respectively.This droop control process is illustrated by the green block.Once the reference voltage and phase angle have been computed by the droop and virtual swing equations, they are fed into the current and voltage control loops, highlighted in the blue block.This block comprises four interlocked PI controllers (each as illustrated in Fig. 12) and modulates the reference voltage using measurements from the inverter and lowpass filter.The required measurements are first transformed into the φ k -rotating reference frame using dq transforms.The equations describing the blue block are given compactly in Eqs.(E8)-(E12).After passing through the PI controllers, the resulting dq signal is sent through an inverse dq transformer to convert it back to a three-phase sinusoid.This is then passed into the igbt inverter in Fig. 1, via a pulse-width-modulator (not illustrated).Note that the current and control loops are ubiquitous to most control schemes.Only the orange and green blocks will change between different schemes.
The rate of change of the drift of this reference voltage from the actual voltage is then tracked in each case: This is then used within the first two PI controllers to compute a reference internal current: The form of Eq. (E10) comes from the operation of the PI controllers.In this case the output i 0,d/q, * k is computed via a weighted sum of an input v and the input's integral ε d/q k .These are weighted by a proportional K pv and integral K iv coefficient, respectively.Finally, the trailing terms in Eqs.(E10) come from the cross wiring of the PI controllers as illustrated in Fig. 13.
The rate of change of the drift of this reference internal current from the actual internal current is then tracked: and finally, these values are used by the second pair of PI controllers to compute the inner reference voltage: Note that the forms of Eqs.(E12) have a similar structure to Eqs. (E10), since they too describe PI controllers.The inverter will then output three-phase ac voltage at a frequency ω k /2π , with dq amplitudes v 0,d, * k and v 0,q, * k .To re-express Eqs.(E1)-(E12) as a coupled system in (ω k , φ k , E k , ε d/q k , γ d/q k ) as presented in Sec.II C, we combine as follows.First Eqs.(E5) in (E1) gives and then Eqs.(E4) and (E6) into (E3) gives Equation (E8) into (E9): Equations (E8) into (E10): Eqs. (E16) into (E11): Eqs. (E16) into (E12): Equations ( E2), (E13)-(E15), and (E17) are the seven additional ODEs described in Sec.II C, while input voltages and reference currents are computed from these variables via Eqs.(E16) and (E18).We now write these dynamics, together with those of the usual LC filter, in the form: where x k ∈ R 13 .Let us first split the internal dynamics F into linear terms A ∈ R 13×13 and nonlinear and constant terms B : R 13 → R 13 , so that F(x k ) = Ax k + B(x k ).For the VSM, the matrix A has nonzero entries: L is the graph Laplacian, and σ and H are as defined in Sec.II B.
In line with its value as a resource for future mathematical studies, we remark that the majority of the system encoded in Eq. ( E19) is common to a variety of microgrid controller systems.Only the specific dynamics in Eqs.(E1) and (E3) will differ.This is demonstrated in the following subsection.

Droop controller
So far the results presented have been specific to the VSM control scheme.However, most of the modeled components are common to any general inverter control system.All systems will feature a low-pass filter and some form of voltage-current control loop.This generality enables other control systems to be modeled as a dynamical system in the form of Eq. (E19), allowing the stability analysis presented Sec.II D be extended to other control systems.In general, adapting the framework for other control systems will require the modification of Eqs.(E13) and (E14) only.
Here we provide a worked example for a voltagefrequency droop controller, as described in Ref. [37].This droop controller system and variants of it are prolific in the literature.Whereas the VSM employs a second-order swing equation to calculate the inverter phase angle, the voltage-frequency controller uses the first-order ODE: where m p is a droop coefficient.This gives a 12dimensional system for each node.When written in the form of Eq. (E19), the entries of the matrix A and vector B are where ω = ω g + m p (P ref − P).Note the similarity to the equivalent entities in Appendix E 1.
Using the master stability framework as before, the maximum number of nodes before instability as a function of R d are plotted in Fig. 6.

APPENDIX F: COMPLEX LOAD CASE
Section II B demonstrated how a microgrid with a purely Ohmic load can be reduced to a dynamically equivalent complete graph where all nodes are of the same type.This allows easy application of the MSF method.Here, we describe how the same process can be achieved in the case of a more general complex load.
First, recall that each generator system k is coupled to the rest of the grid via its output currents i d k and i q k .This gives output current ODEs for the kth generator of the form where L T and R T are the transmission inductance and voltage, respectively, and i d/q k and v d/q k are the generator's internal current and voltage for each channel.The coupling depends upon the voltage at the PCC v d/q PCC , which in turn depends upon the type of grid load.In Sec.II B, the grid load was assumed to be purely Ohmic.This meant the PCC voltage was simply a linear sum over the currents of each generator v Section II B exploited this to show how this load term could then be distributed amongst all of the generators, producing an all-to-all network with Laplacian coupled output currents.Here, we consider the inductive load case.For simplicity, we will work in the single-phase nonrotating reference frame.
For a single-phase microgrid with a complex load, the PCC voltage is where R d and L d are the load resistance and inductance, respectively.The output current for the kth generator system is given by For the purposes of clarity of exposition, we consider graph reduction for the case of a network of two DGs (n = 2).Let x := (di 1 /dt) and y := di 2 /dt.Then we can write Eq. (F3) for the two nodes as where Eliminating y in favour of x and returning to original variables provides describing the output current of DG1; the corresponding equation for DG2 is obtained by transposing the indices 1 and 2. This process has eliminated the PCC node and incorporated its dynamics into the two DGs.The system is now a network with two nodes of identical type.An identical method can be applied for networks with n ≥ 3, although the algebra becomes increasingly lengthy.We obtain the general form of the output current ODE for the kth DG in a network of size n as where Exploiting the relationship n l=1 i l = i k + n l=1 A lk i l , where A = D − L is the n by n adjacency matrix of a complete graph on n nodes with Laplacian L and degree matrix D, and after some manipulation, we obtain where This is the form of single-phase output current equation for each node k in a microgrid with n generators.For a three-phase system we have one ODE of the form Eq. (F9) for each phase a, b, and c.Applying the dq transform, as described in Appendix A one obtains an ODE each for the d and q channels of the current: (F11) These two ODEs, combined with the ODEs in Eqs. ( 1)-( 4), which describe the voltage and current evolution for each DG, provide a set of ODEs for each node k describing the full dynamics of the network.This can be written in the form of Eq. ( 14): where the coupling function H is a linear function: The master stability function for this system is shown in Fig. 14(a) for two different values of L d .The MSF was computed using the algorithm outlined in Appendix D. Figure 14(a) demonstrates that the critical point at which the MSF crosses the axis is not substantially different in the case of a purely Ohmic case (L d = 0 H) and for a highly inductive load (L d = 20 mH).

APPENDIX G: EQUATIONS FOR PERTURBED GRID TOPOLOGIES
Section III A considered nontree grid arrangements.Here, we develop the ODEs necessary to timestep such systems.Recall that in the tree case, with n DGs connected directly to a point of common coupling bus PCC with voltage v PCC , the k th node's output current i k evolves according to where L T and R T are the transmission inductance and resistance, respectively and v k is the node's voltage at the point of its low-pass filter.The current and voltage variables here are all three-phase vectors.After a dq transform, Eq. (G1) yields Eqs. ( 6) and (7).Now consider the case where some subset of the nodes, labeled S, are connected to an intermediate bus before connecting to the PCC.This is illustrated in Figs.15(a) and 15(b).In this case each node connected directly to the PCC will have an output current evolving according to Eq. (G1), while the nodes in the set S will have output currents evolving according to where R ∞ is a virtual resistance whose value is typically chosen to be around 10 8 or higher and is used to increase the rank of the circuit equations without altering current flow.The current i S of the combined DGs in S and evolves according to where the PCC voltage in this case is After applying a dq transform to Eq. (G3), one obtains an extra two ODEs that must be included in the system to account for the current on the line between the two buses.This can be easily generalized to the case where there are multiple distinct sets, each grouped into a distinct intermediate bus before the PCC.In this case there would be two additional ODEs following the form of Eq. (G3) for each grouping.Deeper levels of branching can be accommodated by including two additional equations for the form Eq. (G3) at each bus in order to model the dynamics of the current formed of the sum of the currents of all of the child nodes of the bus.Next, consider the case illustrated in Figure 15(c) where two adjacent DGs, with output currents i 1 and i 2 , respectively, are linked together by a bridge carrying a current i B before joining the the PCC.This creates a looped network.Figure 15(d) shows a detailed view of the looped area.The output currents i k of any other DGs which are linked directly to the PCC evolve according to Eq. (G1), while i 1 and i 2 evolve according to where v 1 and v 2 are the voltages of the two paired nodes measured at each of their low-pass filters.The currents i s 1 , i s 2 , and i B evolve according to and The PCC voltage in this looped case is proportional to the sum of all currents incident on the PCC bus; that is the two currents associated with the loop, i s 1 and i s 2 , together with the remaining n − 2 currents from the DGs that are connected directly: After dq transforming Eqs.(G7)-(G9), one obtains six additional ODEs that must be included in the system to model the loop.This can be extended to multiple loops, where each loop would be described by six such equations.... ...

Direct simulation
All models are implemented in Python, with Python's integrate library [46] used for numerical integration.This uses either a Runge-Kutta implementation or, if the system is stiff, the LSODA algorithm [47].For computational speed up, the system's Hessian matrix is supplied during direct simulation.

Continuation methods for the swing equation
The system must be cured of its zero eigenvalue.To do this, we form a new variable θk := θ k − θ n , which is the kth node's phase angle minus the phase angle of the n th node.Let θ be the n − 1 dimensional vector of each nodes reduced phase, excluding that of the nth node, which is zero anyway.Then the reduced dynamics can be written in n − 1 dimensional matrix vector form as where Ẽ is the reduced network's incidence matrix and Ãn is an (n − 1) × (n − 1) matrix whose rows are duplicates of the nth row of the adjacency matrix of the full network.This produces a system whose Jacobian is of full rank and can then be used to find steady states via Newton's method.These steady states are then continued using arclength continuation via Julia's BifurcationKit suite [48].

FIG. 1 .
FIG.1.General structure of a DG.A dc source on the left is converted into three approximate sinusoidal signals each with a frequency of ω and a voltage amplitude of v 0 .These three signals then pass through a low-pass filter, comprising inductors L and capacitors C in series and parallel, respectively.These scrub high-frequency oscillations from the signals, producing smooth three-phase sinusoidal voltage v and current i to be transmitted to the rest of the network.

FIG. 5 .FIG. 6 .
FIG. 5. Microgrid stability results under VSM control.(a) The maximum real part of the Floquet exponent μ max as a function of α = γ L T .The critical value is α = 5778 .(b)The maximum number of nodes before instability as a function of R d for increasing inertia for different values of J .These curves are computed using Eq.(32).(c) The value of α c as a function of J .Increasing inertia J leads to higher values of α c , corresponding to a larger window of stability.Parameter values are given in Table I except as stated.

FIG. 7 .
FIG.7.The effects of heterogeneity.Diagrams (a)-(c) illustrate three different grid topologies with varying levels of heterogeneity.(a) shows four generators coupled directly to a PCC bus in a tree configuration.(b) shows the same system but with a branch added between two of the generators creating a looped network.(c) shows an asymmetric nested-tree arrangement.The bar chart shows the value of R d at which each network becomes unstable.The worst performing is the tree network in (a).The ODEs governing the dynamics of the network structures (b) and (c) are given in Appendix G. Panel (d) shows simulation results for an ensemble of six-node VSM networks.In each case the local load of each node is drawn from a normal distribution with variance ν, and mean R dl .Each system is then time stepped to determine if a stable state can be found.The results are then plotted in the (P, ν) plane, where P is the power delivered to the grid and is calculated from R dl using the equations in Appendix C. The shaded region indicates the stable regime for the homogeneous system, computed using the MSF (Sec.II D).

PFIG. 8 .
FIG. 8. (a) Exemplar dynamics of a microgrid network governed by the swing Eq. (33), with n = 6 generators and a uniform power demand on each node P k =: P = 1 W; other model parameters specified in the main text.The phase deviation θ k (t) of each generator is shown in different colors.(b) Numerical continuation of Eq. (33), with P as the bifurcation parameter.

FIG. 9 .
FIG. 9. Demonstration of the dq transform.Three-phase sinusoidal signals in (a) are transformed into two stationary signals in (b).

FIG. 11
FIG. 11.(a) Real (P) and reactive (Q) power (b) voltage, and (c) current, as a function of load resistance R d computed from Eqs. (C7), (C4), and (C6), respectively.The dashed line indicates a nominal safety cutoff of 6A, below which there is a voltage collapse and current spike.This would be unsafe power flow in a domestic setting.

FIG. 14 .
FIG. 14.(a) Master stability as a function of PCC resistance R d for the complex load case with load impedance L d = 0 (red) and L d = 2 mH (blue).(b) The locus of the critical point at which the MSF crosses the axis as a function of load impedance.Across a realistic range of L d spanning 2 orders of magnitude the critical point moves by approximately 10 , corresponding to a power change of only 0.1 W.

FIG. 15 .
FIG. 15.(a) Illustration of a distribution network where a subset S of nodes (highlighted by the dashed box) are connected to an intermediate bus before the PCC bus.Thick black bars indicate the buses.For simplicity in this example there are only two nodes in S. (b) An enlargement of the dashed boxed area in (a) to show circuit detail.All shown inductors and resistors have inductance and resistance values L T and R T , respectively, except for the virtual resistors R ∞ .The i variables denote the various current values of interest and their (arbitrarily) assigned direction.(c) A distribution network where two adjacent nodes are connected forming a loop (highlighted by the dashed box).(d) A detailed view of the dashed boxed area in (c).