Transient and stability analysis of heterogeneous micro-grid networks subject to uncertainties

: In this study, networks of interconnected heterogeneous micro-grids are studied. The transient dynamics is modelled as an averaging process whereby micro-grids are assimilated to dynamic agents in a network. An analysis of the convergence of the consensus dynamics is carried out under different assumptions on the damping and inertia parameters and the topology of the network. This study provides an insight into the relation between the network topology and the system's response. An analysis of the ways in which the heterogeneous inertial parameters affect the transient response of the network is also implemented. Additionally, the conditions that guarantee stability are identified when the system is under the influence of uncertain non-linear parameters. Finally, simulations are carried out based on a model calibrated on an existing network in the UK under parameter uncertainties.


Introduction
This paper presents an analysis of the transient dynamics of a network of micro-grids, mainly the influence of oscillations during the system response. A micro-grid is modelled using the swing dynamics and incorporating parameters for damping and inertia. The interrelation among micro-grids is modelled using the coupled oscillator archetype and the resulting dynamics is described by a graph-Laplacian matrix. The transient analysis is extended to a number of cases to gain insight on the role of the connectivity and the parameters.
Although the presented model is a simplified representation for individual micro-grids, it has been proven to be useful as a tool for the study of smart-grid related subjects, such as demand-side management in [1] and real-time pricing [2]. This paper provides a more in-depth stability analysis and yields some other highlights that might not be entirely practical in nature but are interesting nonetheless.
The parameter heterogeneity that is involved in this study can be useful in the design of modern power systems, since the impact of inertia is a present challenge [3,4], as more low-inertia systems such as renewables are being commonly integrated into current power networks [5].

Main contributions
As a first result, the relation between the transient stability and consensus dynamics is explored under the assumption that the micro-grids are homogeneous, namely every micro-grid in the network has equal parameters. This study sheds light on the ways in which different damping coefficients influence the frequency and power flow consensus values.
Secondly, stability analysis for the heterogeneous case is performed by estimating the system's eigenvalues based on the Gershgorin disc theorem. The conditions that guarantee absolute stability, namely the case in where the system's measurements are subject to non-linearities and uncertainties, are also explored.
Thirdly, simulations are performed using different topologies; reaffirming the ways in which the connectivity of the network affects the time constant of the transient response. Additionally, the present work also involves the adaptation of the proposed model to real instances in the UK electrical networks and the calibration of the nodes' parameters using data of the power capabilities of the micro-grids, simulation results also show the system's response when subject to parameter change over time.

Reviewed literature
Following a previous work of the authors in [6], this paper investigates the interplay between the transmission dynamics and the micro-grid dynamics by obtaining the conditions for stability when the system is subject to uncertainties. We are dealing with input/output systems interconnected through diffusive coupling as in [7,8], we differentiate our work by focusing on the heterogeneous case. This is also a continuation of the study in [9] about the effects of the parameters of homogeneous micro-grids on their transient stability. We now extend the approach to heterogeneous networks. A brief analysis of the influence of the parameters on a micro-grid, along with a simplified version of the model used in the present paper is found in [1]. The role of the Laplacian in the swing dynamics and the correspondence with the Kuramoto coupled oscillator is considered in [10]. The equivalent model for the connection between two micro-grids is based on the reduction shown in [11]. The link between the Laplacian and the swing dynamics for small phase angles is discussed in [12]; we use this approximation in the micro-grid model. The model for a micro-grid subject to uncertainties is explained in [2] based on a game-theoretic approach to describe the disturbances. A study of the power flow and demand response in a distributed system of micro-grids similar to the present one is carried out in [13]. As the main difference, we provide a stability analysis. Transient analysis and study of the impact of the damping to inertia ratio is in the same spirit as in [14], from which conditions obtained are also employed for the non-linearity sector calculation. The role of the damping parameters in electrical generators and the procedure to obtain them is discussed in [15]. We also borrow from [9] the idea of converting the one-line diagram into a dynamic network. This paper is structured as follows. In Section 2, we state the problem and introduce the model. In Sections 3 and 4, we present

Problem statement and model
This study mainly addresses the analysis of the transient dynamics of a network system comprised of interconnected micro-grids and the influence of the parameters and topology of the network on the stability, more specifically, the ways in which the heterogeneous parameters in the network influence the eigenvalues of the overall system. Furthermore, we investigate conditions that guarantee absolute stability, that is, the maximum magnitude of uncertain non-linear parameters that the system can withstand when subject to such uncertainties.
Our approach allows to link the transient response to the connectivity of the equivalent network graph, approximate the eigenvalues' position and bounds depending on the different microgrid parameters, derive conditions for the stability and determine the maximum amplitude of uncertainties in terms of the parameters of each micro-grid.
The model of a single micro-grid i in a network describes the dynamics of the power flow P i in the micro-grid, which follows the first-order differential equation: where f i and f j are the frequencies of micro-grids i and j, respectively, T i j is the synchronising coefficient which represents the maximum power transfer between micro-grids i and j [10] in Fig. 1 shows the equivalent circuit representation for the interconnection of two micro-grids [11], σ i and μ i are the transmission inertia and damping coefficients, respectively. If micro-grid i is connected to multiple other micro-grids, then the first term is a sum of the adjacent micro-grids N i to micro-grid i as we shall explain later. The second term helps describe the characteristic response of the power transmission. From (1), it can be seen that the power flow depends on the frequency error f j − f i . A physical interpretation of this is that if f i < f j then the power flows from micro-grid j to micro-grid i; in contrast if f i > f j then the power flows from micro-grid i to micro-grid j. The model of micro-grid i also involves the dynamics of frequency f i represented by the swing equation [1,2]: where D i and M i denote the micro-grid's swing damping and inertia coefficients, respectively. In electrical systems the damping D i is obtained as a result of changing loads and control loops [15] and is measured in MJ/rad. M i is the moment of inertia caused by the rotors of the electric generators in the micro-grid [16, p. 438] and is measured in MJ-s/rad. Fig. 2 shows the system block representation of the system (1) and (2). The state-space representation of the system is obtained by introducing the state variables as an external input. Model (1) and (2) is rewritten as A system of interconnected micro-grids can be modelled by a graph like the one shown in Fig. 3. Each node represents a microgrid and each edge the power line that connects two of them; the connectivity of a micro-grid is captured by the degree d i of the corresponding node, which is equal to the its number of connections. In the unweighted and undirected case, d i is equal to the number of edges that are incident to node i. By extending (3) to the case of a system of n micro-grids, we obtain the state-space model (5). The block matrix that contains the synchronisation parameters T i j is linked to the graph-Laplacian matrix, given by where its diagonal entries correspond to the sum of the weights of the outgoing edges, while the off-diagonal entries are the weights of the adjacency matrix A of the network. Let us recall that the Laplacian of a graph is expressed as L = D out − A, where D out is a diagonal matrix whose elements are the out-degree of the nodes. The Laplacian matrix is then used to represent the system dynamics in matrix form as follows: where diag(D i /M i ) and diag(μ i /σ i ) denote diagonal matrices with main diagonal entries equal to the damping to inertia ratio and diag(1/M i ) and diag(1/σ i ) are diagonal matrices with main diagonal entries equal to the inverse of the inertial constants M i and σ i of each micro-grid i. The state variables X 1 and X 2 are the vectors of power flows P i and frequencies f i of each micro-grid i for i = 1, …, n. Based on the micro-grid network model introduced above, we now leverage the two following derivations to help study its stability.

Preliminary derivations
In this section, we review a couple of preliminary results on the determinant of the micro-grid network system and the Gershgorin disc theorem which will be used in the following sections to show the ways in which the eigenvalues are obtained and subsequently the conditions for the system's stability.

Transient dynamics of the micro-grid network system
The first preliminary derivation deals with the transient dynamics of system (6). To this purpose, we need to obtain the eigenvalues of matrix A. For an unweighted, undirected network of heterogeneous micro-grids with inertial coefficients M i , σ i and damping coefficients D i , μ i , to find the eigenvalues of system (6), the roots of det(λI − A) must be obtained. Taking A as a square block matrix, its determinant is obtained as (see (7)) By denoting (7) is rewritten as System (6) is stable if all its eigenvalues λ i lie in the left-half side of the complex plane. The following theorem illustrates the ways in which an estimation of the eigenvalues of matrix A can be obtained.

Gershgorin disc theorem
This theorem is a well-known tool used to bound the possible values of the eigenvalues of a square matrix in the complex plane. Let A nn be an n × n matrix and let a i j be its ijth entry. For i ∈ 1, …, n, let the radius R i = ∑ j ≠ i a i j be the sum of the absolute values of the non-diagonal elements in the ith row. Let Δ(a ii , R i ) be the closed disc centred at a ii with radius R i . Such disc is called a Gershgorin disc. Theorem 1: Every eigenvalue λ i of A nn lies within at least one of the Gershgorin discs Δ(a ii , R i ).
Proof: We refer the reader to the original paper by Geršgorin [18] for full details of the proof. □

Stability and response analysis
In this section, we present results regarding the transient of the system. Firstly, an estimation of the eigenvalues of system (6) is provided using the Gershgorin disc theorem. Furthermore, a discussion of the ways in which the eigenvalue that is closest to the origin affects the system's response is presented. Secondly, the eigenvalues are obtained for the case when the damping to inertia ratios are normalised to one and the inertia is either equal to one or has different values for each micro-grid. Thirdly, a procedure to identify regions containing the eigenvalues of the system, is shown.

Eigenvalue location and response bounding
In the following, we focus on obtaining an estimation for the eigenvalues of A taking mainly into account the different damping to inertia ratios of the micro-grids. For the analysis utilising the Gershgorin disc theorem, two sets of discs are obtained. For the first one, we take D i > 0 and M i > 0 for i = 1, 2, …, m. Then we obtain a disc Δ i (1) in the first set which encloses the position of an eigenvalue λ i in the complex plane. Let R (1) = 1/M i , then the disc Every disc Δ i (1) has a radius equal to R = 1/M i and is centred in −D i /M i on the real axis of the complex plane. For the second set of discs denoted by Δ i (2) , let us set μ i > 0 for i = m + 1, m + 2, …, n.
Let R (2) = ∑ j N l i j /σ i , then we obtain a disc Δ i (2) in the second set given by Every disc Δ i (2) has a radius equal to R = ∑ j N l i j /σ i and is centred in −μ i /σ i on the real axis of the complex plane. Here we denote by l i j = T i j the absolute value of the ijth element of the Laplacian L. Let us recall that the spectrum of A is the set of eigenvalues {λ 1 , λ 2 , …, λ n }.
Lemma 1: For the spectrum of matrix A, we have Recalling the Gershgorin disc theorem, all eigenvalues of the system are contained within the union of all areas of the discs. The centre of each disc is situated on each of the diagonal elements of A in (6), the radius of each disc is equal to the sum of the rest of the elements in the matching row. □ In the following, we present some results in the case where the transmission dynamics is much faster than the swing dynamics. This is an assumption that is commonly found in the literature, since it yields the standard swing equation [1,2,6,9] because of the difference in the parameter's magnitude.
Assumption 1: The transmission damping coefficient μ i is much larger than the swing damping coefficient If Assumption 1 holds, let us then assume without loss of generality that the nodes are ordered decreasingly in the damping to inertia ratio as follows: In other words, the ratio −D 1 /M 1 corresponds to the micro-grid with the smallest damping to inertia ratio and is the centre to the disc that encircles the smallest eigenvalue λ 1 . Before presenting the next result, let us define the closest point of each disc to the origin as the upper bound of Δ i (1) and Δ i (2) , respectively, as  (6) is stable if for the upper bound of its smallest eigenvalues it holds Furthermore, the rate of convergence satisfies where x eq is a generic equilibrium point and Υ 2n is an opportune 1 × 2n vector. Proof: Let any point p i ∈ Δ i (1) and similarly, let p i ∈ Δ i (2) be given, it holds that ℜ[p i ] < λ¯i; ℜ[p i ] < λ^i. It follows that if condition (14) holds true, then ℜ[p i ], ℜ[p i ] < 0, and therefore, the real part of the eigenvalues is negative.
By obtaining all eigenvalues {λ 1 , …, λ n } of A, an eigenvector matrix V can be computed, as well as its inverse W = V −1 . The modal transformation of A is obtained from which results in a diagonal matrix Λ where all elements contain an eigenvalue of A. The response of the system for a given initial state x(0) and zero input is now expressed as the rate of decay of the smallest eigenvalue λ 1 is dominant for the system's response. Since λ 1 upper bounds, the smallest eigenvalue λ 1 , every state of (17) is exponentially bounded by λ 1 . Namely, the system converges to an equilibrium x eq as in (15). From the discs representation in Fig. 4, we can see that the radius of each disc Δ i (2) depends proportionally on the topology by means of l i j of the Laplacian. Therefore, the eigenvalues are tied to the network's connectivity. □ From Fig. 4, it can be seen that the imaginary part of the eigenvalues is bounded by the radii of the discs. Namely, the maximal amplitude of the frequency of the oscillations is bounded by the radius. Moreover, without altering the topology, if the inertia coefficients σ i and M i increase, the discs are shifted to the left with a reduced radius. Conversely, if decreased, the discs shift to the right and the radii expand, leading to larger and faster oscillations in the transient.

Effect of inertia on eigenvalues
In this section, two cases are analysed. In the first one, the eigenvalues of the system are obtained for the ideal case of when all of its parameters are normalised to one. In the second, all inertia coefficients are considered different in order to emphasise the ways in which such parameters affect the transient. A result for each case is shown below. Let us now state the first assumption.
Assumption 2: All damping and inertia coefficients in (6) are unitary, namely μ i = D i = 1, and σ i = M i = 1 for all i, so that This is a strong assumption, however, it is useful for the purpose of isolating the effect of the topology in the network's eigenvalues and subsequently illustrate the rate of convergence towards the consensus value. The above will be relaxed by Assumption 3. Let us denote d max as the maximum degree of all the nodes in the network, namely d max := max i {d i } which identifies the node with most connections and its quantity.
Theorem 3: Let Assumption 2 hold true. Then system (6) is stable. Furthermore, the maximal frequency of the oscillations is bounded by where η i denotes the ith eigenvalue of −(L + I). Taking (19) equal to zero, the eigenvalues of A, which we denote by λ i are then obtained as From (20) and from the fact that by definition all η i ≤ − 1, it can be deduced that ℜ(λ i ) is negative for all eigenvalues, hence the network system is stable. As discussed in [10], the smallest eigenvalue of -L is lower bounded by −2d max . By extension, the lowest bound for the smallest η i is equal to −2d max − 1. From this, we can infer bounds on the argument of the square root of (20) which is the imaginary part of the eigenvalues. This, in turn, establishes that the maximal oscillation frequency of the system's response depends directly on the topology, substituting said bound in (20) we get 2 d max . □ For the second case, since matrices Φ and Γ are diagonal, ΓΦL is equal to the Laplacian L scaled on each row by 1/σ i M i : where l i is the ith row of L. Note that by definition, the eigenvalues of a Laplacian L are non-negative with its smallest eigenvalue equal to zero. Such eigenvalues are non-positive for −L. Scaling L by any positive scalar, will not affect the sign of the eigenvalues but these will be compressed or expanded depending on the scalar. The next result shows that stability is guaranteed under weaker conditions than the ones in Theorem 3 and serves to highlight the effect of the inertia coefficients on the eigenvalues. This can be justified in the sense that, as emphasised by the authors in [3,4] and references therein, a current problem in power systems is the tendency towards low-inertia micro-grids; and knowing the effect The scaling of the Laplacian shifts the discs closer to the origin. Furthermore, by adding the eigenvalues of -I we obtain that all μ i < − 1 and therefore, all eigenvalues are negative. The bound of the imaginary part of the eigenvalues is obtained by substituting the lowest bound for the smallest μ i in (20) which in this case is −2d max /σ max M max − 1. □

Clusterisation
In this subsection, we introduce a tool that helps to discern the area(s) in the complex plane where the eigenvalues can be located. The union of the area of a number of overlapping discs derived from system (6) can be referred to as a cluster.
Theorem 5: The number of clusters is obtained from where ∥[ ⋅ ] denotes the indicator function, and Δ j is any disc in the complex plane. Proof: Depending on the values of D i and M i and ordering the discs as in (12), there can be an instance where the equality is yielded, where the left-hand side describes the maximum distance of a point in Δ i (1) from the origin, and the right-hand side is the minimum distance of any point in Δ i + 1 (1) from the origin, (24) means that there is at least a partial overlap between both discs. A similar inequality can be derived for all Δ i (2) . For the case where two or more discs overlap, suppose there exists a specific value for i denoted by i~ such that satisfies which is the argument of (23). The above condition means that the union of the first i~ discs {Δ 1 , …, Δ i~} is disjoint from the union of the last n − i~ discs {Δ i + 1 , …, Δ n }. Using the indicator function on (25) hints the separation of two clusters, the sum of the times I[ ⋅ ] yields a positive result equals to the number of clusters in which the eigenvalues lie. □

Stability under uncertainties
In this section, we extend the analysis to the case where both frequency and power values in each micro-grid are subject to uncertainties ψ i . According to the proposed method we isolate the uncertainty in the feedback loop, as illustrated in Fig. 5, where ψ i ( ⋅ ) denotes a sector non-linearity. An interpretation of the nonlinearity in the feedback loop can be, for instance that the power and frequency measurements are subject to disturbances. For the mentioned case, system (3) can then be rewritten as where the non-linearities ψ i belong to the sector [K min , K max ]. This implies that the following inequality holds: where K min = − γ 2 I and K max = γ 2 I, γ 2 is a sufficiently small gain that determines the size of the non-linearity sector, and where ς max [ ⋅ ] denotes the maximum singular value of the system's transfer function G(jω) and γ 1 its upper bound; both γ 1 and γ 2 serve as tools to determine the size of the non-linearity that the system tolerates before becoming unstable. Conceptually, it is assumed that ψ i satisfies a sector condition, namely that ψ i is at equilibrium at the origin and is locally Lipschitz in the system output's domain. Fig. 6 shows an example of the sector non-linearity and its bounds. The utility of this method is to determine if the origin is asymptotically stable for all non-linearities in the sector, yielding absolute stability to the system, this is also referred as Lure's problem [19]. Although the existence and uniqueness of a solution to the system can potentially be verified through the Lipschitz condition; the presence of uncertainties complicate the analysis, calling for a substantially different approach as touched in [20] and references therein. An advantage of the sector non-linearity method is that to determine stability, as mentioned above, given a positive real system only ψ i has to be checked to be Lipschitz.

Amplitude of uncertainty
Let us first point up the system's transfer function G(s) derived from (26) as (see (29)) . In the case where the transmission and swing dynamics have similar parameters, we can obtain a sufficient condition for the maximum size of the non-linearity sector for which absolute stability holds. Such value is equal to the square root of the maximum eigenvalue of the transfer function matrix G(jω) multiplied by its conjugate transpose, i.e.
With the purpose of obtaining such value and for the rest of the analysis, the external input x 2 (i) in (26) is assumed to be equal to zero. While obtaining an expression of γ 1 , for illustrative purposes, the following assumption can be made: Assumption 4: Both swing and transmission damping to inertia ratios, have similar values, such that D i /M i ≃ μ i /σ i . Also, we assume D i > M i , namely, the dynamics are over-damped, as discussed in [14].
Theorem 6: Let Assumption 4 hold true. Then the maximum amplitude for the non-linearity sector in (26) is (31) Proof: Due to the complexity of the expression for ς max and for the sake of simplicity, taking δ for the shorthand of the polynomial in (29), G(jω) can be denoted as G T ( − jω) can be defined similarly. ς max comes from the largest eigenvalue of G * G, where ⋅ * is the conjugate transpose, we can obtain it using the determinant Δ and trace T. Expressions for both are simplified by accounting the following: G 12 (jω) and G 21 (jω) have no imaginary part, taking Assumption 4 as true, it holds that G 11 = G 22 . We can also assume the network to be undirected: T i j = 1, yielding G 12 = − G 21 ; G 12 = G 21 . Taking f = G 11 and g = G 12 and using ⋅ to denote the conjugate we obtain The eigenvalues of G * G are then obtained from (33) as follows: λ = 2(g 2 + f 2 ) ± (4/(δδ¯) 2 )(g 2 + f 2 ) 2 − (4/(δδ¯) 2 )(g 2 + f 2 ) 2 2δδ= g 2 + f 2 δδ¯.
Then, ς max is obtained taking the square root as in (30) Substituting all values and expressions in (34), we get Substituting (35) into (28) we get γ 1 = sup ω ∈ R ς max , which refers to the largest value that (35) can get to as ω varies. If Assumption 4 holds, (35) is a monotonically decreasing function, with a least upper bound at ω = 0, obtaining γ 1 = M i 2 /(D i 2 + 1). Since we know from the small-gain theorem [19, p. 411] that γ 1 γ 2 < 1, it can be inferred that inequality (31) holds. □ The above gives us a clearer definition of the size of the nonlinearity sector as a function of the parameters.

Lyapunov approach
In the general case, where there is no simple expression for the maximum singular value, a numerical Lyapunov stability approach can be carried out. Some other alternatives to corroborate stability include the application of a loop transformation of the system into feedback-connected passive dynamical systems, and the utilisation of either the Popov or the circle criterion when applicable [9,19].
Let us denote by V a candidate Lyapunov function V = x T Px for system (26).
Theorem 7: Given a small ε, γ 2 and a symmetric positive definite matrix P that satisfies the Riccati equation: then (26) is absolutely stable and V is a Lyapunov function.
Proof: The derivative of V along the trajectories of system (26) is V˙ is strictly negative if the given V˙ plus a small quantity γ 2 2 ψ T are not larger than the small quadratic function −εx T Lx, namely Then, inequality (38) can be expanded as To validate that (39) holds, let us rewrite it using L = C T C as The negative definiteness of M can be shown by imposing that its Schur complement is negative. Given that −γ 2 ≤ 0, we take the Schur complement of the block −γ 2 2 I: and by setting it to less than or equal to zero, we obtain the Riccati equation (36) itself. □ Other ways in which the linear matrix inequality (40) can be solved are via a graphical method or by recurring to a system of algebraic Riccati equations (AREs).

Simulations
In this section, real instances of power network topologies are simulated. The first one covers the case when the network is considered homogeneous, unweighted and undirected. The second example touches on a different network with different topology and shows the influence of the connectivity on the response. The third set of simulations contains parameter uncertainties, in such case the network is heterogeneous, weighted and directed.

Graph modelling from the existing network
Simulations were carried out using data of the London City Road power network, as found in [17]. Fig. 7 shows a simplified diagram extracted from the one-line diagrams in [17], this contains the names of the generators and their respective load buses; from this, a network graph was derived which is the example in Fig. 3. The graph was modelled as unweighted and undirected, assuming that the influence between any two nodes is bidirectional. The objective of the first set of simulations is to analyse the transient dynamics and investigate the convergence of the frequency and power to the desired reference. When this occurs, the network achieves synchronisation. In the present simulations, all micro-grids are considered homogeneous. The corresponding parameters were selected as follows: the number of nodes n = 10, inertial constants M i = 1, σ i = 1 MJ/rad synchronising coefficients To show that the previous results are scalable, a different section of the London City Road network was selected. The derived undirected unweighted graph is shown in Fig. 10. It is worth mentioning that on average, this topology has 2.75 connections per node, in contrast to the 2.5 of the previous example. The rest of the parameters are unchanged.
The frequency response is shown in Fig. 11. Comparing these results against the previous simulations, it can be seen that under the second topology the system converges about half a second faster; this is more evident in the top plots, where D i = 1 MJ − s/rad, reaffirming our justification for Assumption 2. This implies that a larger connectivity yields a smaller time constant in the overall system. Furthermore, all frequency responses have less oscillations with a smaller magnitude during the transient. We omit showing the power flow plot since it has no significant differences from the previous one.

Parameter varying and heterogeneity
To account for heterogeneity, we now consider the system shown in Fig. 3, where all nodes contain different parameters and the influence from nodes i to j differs to the one from j to i. The following simulations shed light on the transient response when the synchronising coefficient T i j , the damping coefficients D i , μ i and  First, based on the information in [17], a weighted and directed graph has been derived as shown in Fig. 12.
The synchronising coefficients T i j have been selected depending on the power in MVA that flows in and out of each micro-grid as found in [17], i.e. if micro-grid i outputs 60 MVA to micro-grid j, its T i j will vary within the range [59, 61] MVA. This range has been introduced with the aim of including uncertainties in the system. Due to the unavailability of data on the exact parameters of the network, approximations were done in accordance to typical data from Westinghouse in [16, p. 436]. The inertial coefficients M i and σ i depend on the capacity G i of each micro-grid as shown in Table 1. The constant H i is assigned randomly from a range of values in [6,9]. For the swing inertia, we take M i = G i H i /π f i , where f i is the nominal frequency, which in this case is 50 Hz. For the damping constant D i , a random value in the interval [4.5M i , 5.5M i ] is assigned to each micro-grid for the simulation. For illustrative purposes the values of μ i are chosen as μ i = 15D i . As a way to subject the system to non-linearities, the parameters change their value randomly within their assigned range every 0.1 s during the simulations. Also, the states are reset every 3.3 s as in previous examples.
It can be seen in Fig. 13 that the power flow is contained within the acceptable tolerance of 1 MW for the entirety of the simulation. Let us finally mention that the random change of topology every 0.1 s produces barely noticeable oscillations that do not modify the behaviour or the consensus value. For the sake of brevity, we omit to show simulations on the effect of uncertainties that result in an unstable response. However, it is straightforward to choose an uncertainty amplitude value that causes larger measurement variations, namely an amplitude under which condition (31) does not hold for the given inertia and damping parameters.

Conclusions
As a progression of previous works which are focused on networks of homogeneous micro-grids, we have now extended the analysis to the case where heterogeneity is involved in the form of the different parameters for each micro-grid.
We have investigated the transient stability, and shown the ways in which the heterogeneity of the parameters between micro-grids in the network affects both the response and eigenvalues of the overall system. We mainly focused on the inertial parameters since studying their effect is a current issue in the design of modern power systems.
We obtained a few interesting observations regarding the displacement of the eigenvalues, which depends on the multiple heterogeneous parameters in the network; plus the way of deriving the clusterisation of the areas where the eigenvalues might reside in.
We have studied the maximal magnitude of the non-linearities that the system can accept while remaining stable and expressed it as a function of the parameters.
Finally, we have illustrated the scalability of the model by simulating different topologies and shown the role of the connectivity in the network's response.
The future direction of this work involves the analysis of the impact of stochastic disturbances due to renewable generation and demand response under on-line dynamic pricing. Fig. 12 Weighted directed graph of the London City Road network [17]