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

Abstract

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 our results. In Section 6, we provide numerical examples. Finally, in Section 7, we provide conclusions and discuss further directions.

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 micro-grid 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: 1${\stackrel{.}{P}}_i=\sum _{j\in N_i}\frac{T_{ij}}{{\sigma }_i}(f_j-f_i)-\frac{{\mu }_i}{{\sigma }_i}{P}_i,$where $f_i$ and $f_j$ are the frequencies of micro-grids i and j, respectively, $T_{ij}$ is the synchronising coefficient which represents the maximum power transfer between micro-grids i and j [10] in MVA, where $T_{ij}=|V_i\parallel V_j\parallel Y_{ij}|$ and $Y_{ij}$ is the inverse of the impedance $Z_{ij}$ of the transmission line $\{i,j\}$, Fig. 1 shows the equivalent circuit representation for the interconnection of two micro-grids [11], ${\sigma }_i$ and ${\mu }_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]: 2${\stackrel{.}{f}}_i=-\frac{D_i}{M_i}{f}_i+\frac{1}{M_i}{P}_i,$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).

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

The state-space representation of the system is obtained by introducing the state variables $P_i=x_1^{(i)}$, $f_i=x_2^{(i)}$ and taking $f_j=x_2^{(j)}$ as an external input. Model (1) and (2) is rewritten as 3$\left[\begin{array}{c}{\stackrel{.}{x}}_1^{(i)}\\ {\stackrel{.}{x}}_2^{(i)}\end{array}\right]=\left[\begin{array}{cc}-\frac{{\mu }_i}{{\sigma }_i}& -\frac{T_{ij}}{{\sigma }_i}\\ \frac{1}{M_i}& -\frac{D_i}{M_i}\end{array}\right]\left[\begin{array}{c}{x}_1^{(i)}\\ x_2^{(i)}\end{array}\right]+\left[\begin{array}{c}\frac{T_{ij}}{{\sigma }_i}\\ 0\end{array}\right]x_2^{(j)}.$A system of interconnected micro-grids can be modelled by a graph like the one shown in Fig. 3. Each node represents a micro-grid 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_{ij}$ is linked to the graph-Laplacian matrix, given by 4$\mathit{L}:=\left[\begin{array}{cccc}{T}_{11}& -T_{12}& \cdots & -T_{1n}\\ -T_{21}& T_{22}& \cdots & -T_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ -T_{n1}& -T_{n2}& \cdots & T_{nn}\end{array}\right],$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 $\mathit{A}$ of the network. Let us recall that the Laplacian of a graph is expressed as $\mathit{L}={\mathit{D}}_{\mathrm{out}}-\mathit{A}$, where ${\mathit{D}}_{\mathrm{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: 5$\left[\begin{array}{c}{\stackrel{.}{x}}_1^{(1)}\\ ⋮\\ {\stackrel{.}{x}}_1^{(n)}\\ {\stackrel{.}{x}}_2^{(1)}\\ ⋮\\ {\stackrel{.}{x}}_2^{(n)}\end{array}\right]=\left[\begin{array}{cccccc}-\frac{{\mu }_1}{{\sigma }_1}& \cdots & 0& -\frac{T_{11}}{{\sigma }_1}& \cdots & \frac{T_{1n}}{{\sigma }_1}\\ 0& \ddots & 0& & \ddots & \\ 0& \cdots & -\frac{{\mu }_n}{{\sigma }_n}& \frac{T_{n1}}{{\sigma }_n}& \cdots & -\frac{T_{nn}}{{\sigma }_n}\\ \frac{1}{M_1}& \cdots & 0& -\frac{D_1}{M_1}& \cdots & 0\\ 0& \ddots & 0& 0& \ddots & 0\\ 0& \cdots & \frac{1}{M_n}& 0& \cdots & -\frac{D_n}{M_n}\end{array}\right]\left[\begin{array}{c}{x}_1^{(1)}\\ ⋮\\ x_1^{(n)}\\ x_2^{(1)}\\ ⋮\\ x_2^{(n)}\end{array}\right].$Let ${\mathit{X}}_j=[x_j^{(i)}{]\phantom{a}}_{i=1,\dots ,n}$ then 6$\left[\begin{array}{c}{\stackrel{.}{\mathit{X}}}_1\\ {\stackrel{.}{\mathit{X}}}_2\end{array}\right]=\stackrel{\mathit{A}}{\overbrace{\left[\begin{array}{cc}-\mathrm{diag}\left(\frac{{\mu }_i}{{\sigma }_i}\right)& -\mathrm{diag}\left(\frac{1}{{\sigma }_i}\right)L\\ \mathrm{diag}\left(\frac{1}{M_i}\right)& -\mathrm{diag}\left(\frac{D_i}{M_i}\right)\end{array}\right]}}\left[\begin{array}{c}{\mathit{X}}_1\\ {\mathit{X}}_2\end{array}\right],$where $\mathrm{diag}(D_i/M_i)$ and $\mathrm{diag}({\mu }_i/{\sigma }_i)$ denote diagonal matrices with main diagonal entries equal to the damping to inertia ratio and $\mathrm{diag}(1/M_i)$ and $\mathrm{diag}(1/{\sigma }_i)$ are diagonal matrices with main diagonal entries equal to the inverse of the inertial constants $M_i$ and ${\sigma }_i$ of each micro-grid i. The state variables ${\mathit{X}}_1$ and ${\mathit{X}}_2$ are the vectors of power flows $P_i$ and frequencies $f_i$ of each micro-grid i for $i=1,\dots ,n$. Based on the micro-grid network model introduced above, we now leverage the two following derivations to help study its stability.

[IMAGE OMITTED. SEE PDF]

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 $\mathit{A}$. For an unweighted, undirected network of heterogeneous micro-grids with inertial coefficients $M_i,{\sigma }_i$ and damping coefficients $D_i,{\mu }_i$, to find the eigenvalues of system (6), the roots of $\det (\lambda I-\mathit{A})$ must be obtained. Taking $\mathit{A}$ as a square block matrix, its determinant is obtained as 7$\begin{array}{rl}\det (\lambda I-\mathit{A})& =\det \left(\left[\begin{array}{cc}\lambda I+\mathrm{diag}\left(\frac{{\mu }_i}{{\sigma }_i}\right)& \mathrm{diag}\left(\frac{1}{{\sigma }_i}\right)L\\ -\mathrm{diag}\left(\frac{1}{M_i}\right)& \lambda I+\mathrm{diag}\left(\frac{D_i}{M_i}\right)\end{array}\right]\right)\\ & =\det \left({\lambda }^{2}I+\lambda \left(\mathrm{diag}\left(\frac{D_i}{M_i}+\frac{{\mu }_i}{{\sigma }_i}\right)\right)+\mathrm{diag}\left(\frac{{\mu }_i{D}_i}{{\sigma }_i{M}_i}\right)+\mathrm{diag}\left(\frac{1}{{\sigma }_i{M}_i}\right)L\right).\end{array}$By denoting $\mathbf{\Psi }:=\mathrm{diag}(D_i/M_i)$, $\mathbf{\Phi }:=\mathrm{diag}(1/M_i)$, $\mathbf{\Gamma }:=\mathrm{diag}(1/{\sigma }_i)$ and $\kappa :=\mathrm{diag}({\mu }_i/{\sigma }_i)$, the determinant (7) is rewritten as 8$\det (\lambda I-\mathit{A})=\det ({\lambda }^{2}I+\lambda (\mathbf{\Psi }+\kappa )+\kappa \mathbf{\Psi }+\mathbf{\Gamma }\mathbf{\Phi }L).$System (6) is stable if all its eigenvalues ${\lambda }_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 ${\mathit{A}}_{nn}$ be an $n\times n$ matrix and let $a_{ij}$ be its ijth entry. For $i\in 1,\dots ,n$, let the radius $R_i={\sum }_{j\ne i}\left|a_{ij}\right|$ be the sum of the absolute values of the non-diagonal elements in the ith row. Let $\mathrm{\Delta }(a_{ii},R_i)$ be the closed disc centred at $a_{ii}$ with radius $R_i$. Such disc is called a Gershgorin disc.

$1$

Theorem

Every eigenvalue ${\lambda }_i$ of ${\mathit{A}}_{nn}$ lies within at least one of the Gershgorin discs $\mathrm{\Delta }(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 $\mathit{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,\dots ,m$. Then we obtain a disc ${\mathrm{\Delta }}_i^{(1)}$ in the first set which encloses the position of an eigenvalue ${\lambda }_i$ in the complex plane. Let $R^{(1)}=1/M_i$, then the disc ${\mathrm{\Delta }}_i^{(1)}$ is given by 9${\mathrm{\Delta }}_i^{(1)}\left(-\frac{D_i}{M_i},\frac{1}{M_i}\right)=\left\{\xi :\xi \in \mathbb{C}\parallel \xi +\frac{D_i}{M_i}|\le R^{(1)}\right\}.$Every disc ${\mathrm{\Delta }}_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 ${\mathrm{\Delta }}_i^{(2)}$, let us set ${\mu }_i>0$ for $i=m+1,$$m+2,\dots ,n$. Let $R^{(2)}={\sum }_j^N\left|l_{ij}\right|/{\sigma }_i$, then we obtain a disc ${\mathrm{\Delta }}_i^{(2)}$ in the second set given by 10${\mathrm{\Delta }}_i^{(2)}\left(-\frac{{\mu }_i}{{\sigma }_i},\sum _j^N\frac{1}{{\sigma }_i}\left|l_{ij}\right|\right)=\left\{\xi :\xi \in \mathbb{C}\parallel \xi +\frac{{\mu }_i}{{\sigma }_i}|\le R^{(2)}\right\}.$Every disc ${\mathrm{\Delta }}_i^{(2)}$ has a radius equal to $R={\sum }_j^N\left|l_{ij}\right|/{\sigma }_i$ and is centred in $-{\mu }_i/{\sigma }_i$ on the real axis of the complex plane. Here we denote by $\left|l_{ij}\right|=\left|T_{ij}\right|$ the absolute value of the ijth element of the Laplacian L. Let us recall that the spectrum of A is the set of eigenvalues $\{{\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n\}$.

$1$

Lemma

For the spectrum of matrix A, we have 11$\mathrm{spec}(\mathit{A})\in \bigcup _{i=1}^m{\mathrm{\Delta }}_i^{(1)}\left(-\frac{D_i}{M_i},\frac{1}{M_i}\right)\cup \bigcup _{i=m}^n{\mathrm{\Delta }}_i^{(2)}\left(-\frac{{\mu }_i}{{\sigma }_i},\sum _j^N\frac{1}{{\sigma }_i}\left|l_{ij}\right|\right).$

$$

Proof

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.

$1$

Assumption

The transmission damping coefficient ${\mu }_i$ is much larger than the swing damping coefficient $D_i$, ${\mu }_i\gg D_i$.

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: 12$\begin{array}{rl}& -\frac{{\mu }_n}{{\sigma }_n}<-\frac{{\mu }_{n-1}}{{\sigma }_{n-1}}\ll -\frac{{\mu }_{m+2}}{{\sigma }_{m+2}}<-\frac{{\mu }_{m+1}}{{\sigma }_{m+1}}\\ & \ll -\frac{D_m}{M_m}<-\frac{D_{m-1}}{M_{m-1}}\ll -\frac{D_2}{M_2}<-\frac{D_1}{M_1}.\end{array}$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 ${\lambda }_1$.

Before presenting the next result, let us define the closest point of each disc to the origin as the upper bound of ${\mathrm{\Delta }}_i^{(1)}$ and ${\mathrm{\Delta }}_i^{(2)}$, respectively, as 13${\overline{\lambda }}_i:=\frac{1-D_i}{M_i},{\hat{\lambda }}_i=\sum _j^N\frac{1}{{\sigma }_i}\left|l_{ij}\right|-\frac{{\mu }_i}{{\sigma }_i}.$Fig. 4 shows a possible configuration of the discs and their upper bounds in the complex plane.

[IMAGE OMITTED. SEE PDF]

$2$

Theorem

System (6) is stable if for the upper bound of its smallest eigenvalues it holds 14${\underset{\_}{\lambda }}_1:=\underset{i}{max}\{{\overline{\lambda }}_i,{\hat{\lambda }}_i\}<0.$Furthermore, the rate of convergence satisfies 15$|x(t)-x_{\mathrm{eq}}|\le {\mathbf{{\rm Y}}}_{2n}{\mathrm{e}}^{{\underset{\_}{\lambda }}_{1}t},$where $x_{\mathrm{eq}}$ is a generic equilibrium point and ${\mathrm{{\rm Y}}}_{2n}$ is an opportune $1\times 2n$ vector.

$$

Proof

Let any point ${\overline{p}}_i\in {\mathrm{\Delta }}_i^{(1)}$ and similarly, let ${\hat{p}}_i\in {\mathrm{\Delta }}_i^{(2)}$ be given, it holds that $\mathrm{\Re }[{\overline{p}}_i]<{\overline{\lambda }}_i;\mathrm{\Re }[{\hat{p}}_i]<{\hat{\lambda }}_i$. It follows that if condition (14) holds true, then $\mathrm{\Re }[{\overline{p}}_i],\mathrm{\Re }[{\hat{p}}_i]<0$, and therefore, the real part of the eigenvalues is negative.

By obtaining all eigenvalues $\{{\lambda }_1,\dots ,{\lambda }_n\}$ of $\mathit{A}$, an eigenvector matrix $\mathit{V}$ can be computed, as well as its inverse $\mathit{W}={\mathit{V}}^{-1}$. The modal transformation of $\mathit{A}$ is obtained from 16$\mathit{W}\mathit{A}\mathit{V}=\mathrm{diag}(\{{\lambda }_1,\dots ,{\lambda }_n\})=\mathbf{\Lambda },$which results in a diagonal matrix $\mathbf{\Lambda }$ where all elements contain an eigenvalue of $\mathit{A}$. The response of the system for a given initial state $x(0)$ and zero input is now expressed as 17$x(t)=\mathit{V}{\mathrm{e}}^{\mathbf{\Lambda }t}\mathit{W}x(0),$the rate of decay of the smallest eigenvalue ${\lambda }_1$ is dominant for the system's response. Since ${\underset{\_}{\lambda }}_1$ upper bounds, the smallest eigenvalue ${\lambda }_1$, every state of (17) is exponentially bounded by ${\underset{\_}{\lambda }}_1$. Namely, the system converges to an equilibrium $x_{\mathrm{eq}}$ as in (15). From the discs representation in Fig. 4, we can see that the radius of each disc ${\mathrm{\Delta }}_i^{(2)}$ depends proportionally on the topology by means of $\left|l_{ij}\right|$ 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 ${\sigma }_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.

$2$

Assumption

All damping and inertia coefficients in (6) are unitary, namely ${\mu }_i=D_i=1$, and ${\sigma }_i=M_i=1$ for all i, so that $\kappa =\mathbf{\Gamma }=\mathbf{\Phi }=\mathrm{\Psi }=I\in {\mathbb{R}}^n$.

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.

$3$

Theorem

Let Assumption 2 hold true. Then system (6) is stable. Furthermore, the maximal frequency of the oscillations is bounded by 18$\sqrt{2}\sqrt{d_{\max }}.$

$$

Proof

From Assumption 2, the determinant (8) is rewritten as 19$\det (({\lambda }^2+2\lambda )I+\mathit{L}+I)=\prod _i^n(({\lambda }^2+2\lambda )I+{\eta }_i),$where ${\eta }_i$ denotes the ith eigenvalue of $-(\mathit{L}+I)$. Taking (19) equal to zero, the eigenvalues of $\mathit{A}$, which we denote by ${\lambda }_i$ are then obtained as 20${\lambda }_{i,i+1}=\frac{-2\pm \sqrt{4+4{\eta }_i}}{2},\mathrm{for}i=1,\dots ,n.$From (20) and from the fact that by definition all ${\eta }_i\le -1$, it can be deduced that $\mathrm{\Re }({\lambda }_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 ${\eta }_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 $\sqrt{2}\sqrt{d_{\max }}$. □

For the second case, since matrices $\mathbf{\Phi }$ and $\mathbf{\Gamma }$ are diagonal, $\mathbf{\Gamma }\mathbf{\Phi }\mathit{L}$ is equal to the Laplacian $\mathit{L}$ scaled on each row by $1/{\sigma }_i{M}_i$: 21$\mathbf{\Gamma }\mathbf{\Phi }\mathit{L}=\left[\begin{array}{ccc}\frac{l_{11}}{{\sigma }_1{M}_1}& \cdots & \frac{l_{1n}}{{\sigma }_1{M}_1}\\ & \ddots & \\ \frac{l_{n1}}{{\sigma }_n{M}_n}& \cdots & \frac{l_{nn}}{{\sigma }_n{M}_n}\end{array}\right]=\left[\begin{array}{c}\frac{1}{{\sigma }_1{M}_1}{l}_1\\ ⋮\\ \frac{1}{{\sigma }_n{M}_n}{l}_n\end{array}\right],$where $l_i$ is the ith row of $\mathit{L}$. Note that by definition, the eigenvalues of a Laplacian $\mathit{L}$ are non-negative with its smallest eigenvalue equal to zero. Such eigenvalues are non-positive for $-\mathit{L}$. Scaling $\mathit{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 of different inertia parameters in the network is pivotal in the stability of the overall system.

$3$

Assumption

All damping to inertia ratios in (6) are unitary, namely ${\mu }_i/{\sigma }_i=D_i/M_i=1$ for all i, so that $\mathrm{\Psi }=\kappa =I\in {\mathbb{R}}^n$.

Let us define ${\sigma }_{\max }$ and $M_{\max }$ as the largest inertia coefficients in the system, namely ${\sigma }_{\max }:={max}_i\{{\sigma }_i\}$, $M_{\max }:={max}_i\{M_i\}$.

$4$

Theorem

Let Assumption 3 hold true. Then system (6) is stable. Furthermore, the maximum frequency of the oscillations is upper bounded by 22$\sqrt{2}\sqrt{d_{\max }/{\sigma }_{\max }{M}_{\max }}.$

$$

Proof

With Assumption 3 in mind, the same expression (20) can be obtained from the determinant (8), substituting ${\eta }_i$ with ${\hat{\mu }}_i$ which is the ith eigenvalue of $-(\mathbf{\Gamma }\mathbf{\Phi }\mathit{L}+I)$.

The scaling of the Laplacian shifts the discs closer to the origin. Furthermore, by adding the eigenvalues of - I we obtain that all ${\hat{\mu }}_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 ${\hat{\mu }}_i$ in (20) which in this case is $-2d_{\max }/{\sigma }_{\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.

$5$

Theorem

The number of clusters is obtained from 23$\sum _i\parallel \left[\bigcup _{j=1}^i{\mathrm{\Delta }}_j\cap \bigcup _{j=i+1}^n{\mathrm{\Delta }}_j=\mathrm{\varnothing }\right],$where $\mathbb{\parallel }[\cdot ]$ denotes the indicator function, and ${\mathrm{\Delta }}_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 24$-\frac{1+D_i}{M_i}<\frac{1-D_{i+1}}{M_{i+1}},\mathrm{for}\mathrm{any}i\in \{1,\dots ,m\}$is yielded, where the left-hand side describes the maximum distance of a point in ${\mathrm{\Delta }}_i^{(1)}$ from the origin, and the right-hand side is the minimum distance of any point in ${\mathrm{\Delta }}_{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 ${\mathrm{\Delta }}_i^{(2)}$. For the case where two or more discs overlap, suppose there exists a specific value for i denoted by $\tilde{i}$ such that satisfies 25$\bigcup _{j=1}^{\tilde{i}}{\mathrm{\Delta }}_j\cap \bigcup _{j=\tilde{i}+1}^n{\mathrm{\Delta }}_j=\mathrm{\varnothing },$which is the argument of (23). The above condition means that the union of the first $\tilde{i}$ discs $\{{\mathrm{\Delta }}_1,\dots ,{\mathrm{\Delta }}_{\tilde{i}}\}$ is disjoint from the union of the last $n-\tilde{i}$ discs $\{{\mathrm{\Delta }}_{i+1},\dots ,{\mathrm{\Delta }}_n\}$. Using the indicator function on (25) hints the separation of two clusters, the sum of the times $\mathbb{I}[\cdot ]$ 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 ${\psi }_i$. According to the proposed method we isolate the uncertainty in the feedback loop, as illustrated in Fig. 5, where ${\psi }_i(\cdot )$ denotes a sector non-linearity. An interpretation of the non-linearity in the feedback loop can be, for instance that the power and frequency measurements are subject to disturbances.

[IMAGE OMITTED. SEE PDF]

For the mentioned case, system (3) can then be rewritten as 26$\left[\begin{array}{c}{\stackrel{.}{x}}_1^{(i)}\\ {\stackrel{.}{x}}_2^{(i)}\end{array}\right]=\underset{\mathit{A}}{\underbrace{\left[\begin{array}{cc}-\frac{{\mu }_i}{{\sigma }_i}& -\frac{T_{ij}}{{\sigma }_i}\\ \frac{1}{M_i}& -\frac{D_i}{M_i}\end{array}\right]}}\left[\begin{array}{c}{x}_1^{(i)}\\ x_2^{(i)}\end{array}\right]+\underset{\mathit{B}}{\underbrace{\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]}}\left[\begin{array}{c}{\psi }_1(x_1^{(i)})\\ {\psi }_2(x_2^{(i)})\end{array}\right]+\underset{\mathit{U}}{\underbrace{\left[\begin{array}{c}\frac{T_{ij}}{{\sigma }_i}\\ 0\end{array}\right]}}{x}_2^{(j)},$where the non-linearities ${\psi }_i$ belong to the sector $[K_{\min },K_{\max }]$. This implies that the following inequality holds: 27$[\mathit{\psi }(\mathit{x})-{\mathit{K}}_{min}\mathit{x}){]\phantom{a}}^{\mathrm{T}}[\psi (x)-K_{\max }x)]\le 0,$where $K_{\min }=-{\gamma }_{2}I$ and $K_{\max }={\gamma }_{2}I$, ${\gamma }_2$ is a sufficiently small gain that determines the size of the non-linearity sector, and 28${\gamma }_1=\underset{\omega \in R}{\sup }{\varsigma }_{\max }[\mathit{G}(\mathrm{j}\omega )],$where ${\varsigma }_{\max }[\cdot ]$ denotes the maximum singular value of the system's transfer function $\mathit{G}(\mathrm{j}\omega )$ and ${\gamma }_1$ its upper bound; both ${\gamma }_1$ and ${\gamma }_2$ serve as tools to determine the size of the non-linearity that the system tolerates before becoming unstable. Conceptually, it is assumed that ${\psi }_i$ satisfies a sector condition, namely that ${\psi }_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 ${\psi }_i$ has to be checked to be Lipschitz.

[IMAGE OMITTED. SEE PDF]

Amplitude of uncertainty

Let us first point up the system's transfer function $G(s)$ derived from (26) as 29$\begin{array}{r}G(s)=\frac{1}{s^2+(({\mu }_i/{\sigma }_i)+(D_i/M_i))s+((T_{ij}+D_i{\mu }_i)/{\sigma }_i{M}_i)}\left[\begin{array}{cc}\frac{D_i}{M_i}+s& -\frac{T_{ij}}{{\sigma }_i}\\ \frac{1}{M_i}& \frac{{\mu }_i}{{\sigma }_i}+s\end{array}\right].\end{array}$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 $\mathit{G}(\mathrm{j}\omega )$ multiplied by its conjugate transpose, i.e. 30${\varsigma }_{max}=\sqrt{{\lambda }_{max}[{\mathit{G}}^{\mathrm{T}}(-\mathrm{j}\omega )\mathit{G}(\mathrm{j}\omega )]}.$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 ${\gamma }_1$, for illustrative purposes, the following assumption can be made:

$4$

Assumption

Both swing and transmission damping to inertia ratios, have similar values, such that $D_i/M_i\simeq {\mu }_i/{\sigma }_i$. Also, we assume $D_i>M_i$, namely, the dynamics are over-damped, as discussed in [14].

$6$

Theorem

Let Assumption 4 hold true. Then the maximum amplitude for the non-linearity sector in (26) is 31${\gamma }_2<1/\sqrt{M_i^2/(D_i^2+1)}.$

$$

Proof

Due to the complexity of the expression for ${\varsigma }_{\max }$ and for the sake of simplicity, taking $\delta$ for the shorthand of the polynomial in (29), $\mathit{G}(\mathrm{j}\omega )$ can be denoted as 32$\mathit{G}(\mathrm{j}\omega )=\frac{1}{\delta (\mathrm{j}\omega )}\left[\begin{array}{cc}{G}_{11}(\mathrm{j}\omega )& G_{12}(\mathrm{j}\omega )\\ G_{21}(\mathrm{j}\omega )& G_{22}(\mathrm{j}\omega )\end{array}\right],$ ${\mathit{G}}^{\mathrm{T}}(-\mathrm{j}\omega )$ can be defined similarly. ${\varsigma }_{\max }$ comes from the largest eigenvalue of ${\mathit{G}}^{\ast }\mathit{G}$, where ${\cdot }^{\ast }$ is the conjugate transpose, we can obtain it using the determinant $\mathrm{\Delta }$ and trace T. Expressions for both are simplified by accounting the following: $G_{12}(\mathrm{j}\omega )$ and $G_{21}(\mathrm{j}\omega )$ have no imaginary part, taking Assumption 4 as true, it holds that $\left|G_{11}\right|=\left|G_{22}\right|$. We can also assume the network to be undirected: $T_{ij}=1$, yielding $G_{12}=-G_{21};\left|G_{12}\right|=\left|G_{21}\right|$. Taking $f=\left|G_{11}\right|$ and $g=\left|G_{12}\right|$ and using $\overline{\cdot }$ to denote the conjugate we obtain 33$T=\frac{2}{\delta \overline{\delta }}(g^2+f^2),\mathrm{\Delta }=\frac{1}{{(\delta \overline{\delta })}^2}(g^2+f^2{)\phantom{a}}^2.$The eigenvalues of ${\mathit{G}}^{\ast }\mathit{G}$ are then obtained from (33) as follows: $\lambda =\frac{2(g^2+f^2)\pm \sqrt{(4/{(\delta \overline{\delta })}^2){(g^2+f^2)}^2-(4/{(\delta \overline{\delta })}^2){(g^2+f^2)}^2}}{2\delta \overline{\delta }}=\frac{g^2+f^2}{\delta \overline{\delta }}.$Then, ${\varsigma }_{\max }$ is obtained taking the square root as in (30) 34${\varsigma }_{\max }=\sqrt{(g^2+f^2)/\delta \overline{\delta }}.$Substituting all values and expressions in (34), we get 35${\varsigma }_{\max }=\sqrt{\frac{D_i^2}{M_i^2}+{\omega }^2+\frac{1}{M_i^2}/\frac{{(1+D_i^2)}^2}{M_i^4}+\frac{2(D_i^2-1){\omega }^2}{M_i^2}+{\omega }^4}$Substituting (35) into (28) we get ${\gamma }_1={\sup }_{\omega \in R}{\varsigma }_{\max }$, which refers to the largest value that (35) can get to as $\omega$ varies. If Assumption 4 holds, (35) is a monotonically decreasing function, with a least upper bound at $\omega =0$, obtaining ${\gamma }_1=\sqrt{M_i^2/(D_i^2+1)}$. Since we know from the small-gain theorem [19, p. 411] that ${\gamma }_1{\gamma }_2<1$, it can be inferred that inequality (31) holds. □

The above gives us a clearer definition of the size of the non-linearity 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 $\mathit{V}$ a candidate Lyapunov function $\mathit{V}={\mathit{x}}^{\mathrm{T}}\mathit{P}x$ for system (26).

$7$

Theorem

Given a small $\epsilon$, ${\gamma }_2$ and a symmetric positive definite matrix $\mathit{P}$ that satisfies the Riccati equation: 36$\mathit{P}\mathit{A}+{\mathit{A}}^{\mathrm{T}}\mathit{P}+\epsilon {\mathit{C}}^{\mathrm{T}}C+\frac{1}{{\gamma }_2^2}\mathit{P}{\mathit{B}\mathit{B}}^{\mathrm{T}}\mathit{P}\le 0,$then (26) is absolutely stable and $\mathit{V}$ is a Lyapunov function.

$$

Proof

The derivative of $\mathit{V}$ along the trajectories of system (26) is 37$\stackrel{.}{V}={\mathit{x}}^{\mathrm{T}}(\mathit{P}\mathit{A}+{\mathit{A}}^{\mathrm{T}}\mathit{P})\mathit{x}-2{\mathit{x}}^{\mathrm{T}}\mathit{P}B\mathit{\psi }.$ $\stackrel{.}{V}$ is strictly negative if the given $\stackrel{.}{V}$ plus a small quantity ${\gamma }_2^2{\mathit{\psi }}^{\mathrm{T}}$ are not larger than the small quadratic function $-\epsilon {\mathit{x}}^{\mathrm{T}}\mathit{L}\mathit{x}$, namely 38$\stackrel{.}{V}+{\gamma }_2^2{\mathit{\psi }}^{\mathrm{T}}\mathit{\psi }\le -\epsilon {\mathit{x}}^{\mathrm{T}}\mathit{L}\mathit{x}.$Then, inequality (38) can be expanded as 39${\mathit{x}}^{\mathrm{T}}(\mathit{P}\mathit{A}+{\mathit{A}}^T\mathit{P})\mathit{x}-2{\mathit{x}}^{\mathrm{T}}\mathit{P}\mathit{B}\mathit{\psi }+{\gamma }_2^2{\mathit{\psi }}^{\mathrm{T}}\mathit{\psi }\le -\epsilon {\mathit{x}}^{\mathrm{T}}\mathit{L}\mathit{x}.$To validate that (39) holds, let us rewrite it using $\mathit{L}={\mathit{C}}^{\mathrm{T}}\mathit{C}$ as 40$\left[\begin{array}{cc}{\mathit{x}}^{\mathrm{T}}& {\mathit{\psi }}^{\mathrm{T}}\end{array}\right]\stackrel{\mathit{M}}{\overbrace{\left[\begin{array}{cc}\mathit{P}\mathit{A}+{\mathit{A}}^{\mathrm{T}}\mathit{P}+\epsilon {\mathit{C}}^{\mathrm{T}}\mathit{C}& -\mathit{P}\mathit{B}\\ -{\mathit{B}}^{\mathrm{T}}\mathit{P}& -{\gamma }_2^{2}I\end{array}\right]}}\left[\begin{array}{c}\mathit{x}\\ \mathit{\psi }\end{array}\right]\le 0.$The negative definiteness of $\mathit{M}$ can be shown by imposing that its Schur complement is negative. Given that $-{\gamma }_2\le 0$, we take the Schur complement of the block $-{\gamma }_2^{2}I$: $\mathit{M}/[-{\gamma }_2^{2}I]:=\mathit{P}\mathit{A}+{\mathit{A}}^{\mathrm{T}}\mathit{P}+\epsilon {\mathit{C}}^{\mathrm{T}}C-(\mathit{P}B({\gamma }_2^{2}I{)\phantom{a}}^{-1}(-{\mathit{B}}^{\mathrm{T}}\mathit{P})),$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$, ${\sigma }_i=1\mathrm{MJ}/\mathrm{rad}$ synchronising coefficients $T_{ij}=1\mathrm{MVA}$, the number of iterations $N=1000$, the step size $\mathrm{d}t=0.01\mathrm{s}$. For illustrative reasons, different damping constants $D_i=1,2,4$ and ${\mu }_i=5,10,20\mathrm{MJ}-\mathrm{s}/\mathrm{rad}$ are used for different runs. The initial states of the frequency and power are obtained as random values in the interval [–0.5, 0.5]. The frequency and power variables are also reset every 3.3 s as a way to simulate periodic disturbances. The resulting plots have been scaled around 50 Hz and 30 MW for the frequency and power flow, respectively, to represent realistic values.

[IMAGE OMITTED. SEE PDF]

Fig. 8 shows the frequency response of each micro-grid. It can be seen that the response remains in the range between [49.5, 50.5] Hz and does not exceed in magnitude the desired frequency by >1 Hz. Fig. 9 displays the power flow of each micro-grid, the values remain in the range of [29.5, 30.5] MW as well. In both plots, different damping values are used from top to bottom. Observe that for larger values, both oscillations and settling times are reduced.

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

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.

[IMAGE OMITTED. SEE PDF]

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\mathrm{MJ}-\mathrm{s}/\mathrm{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.

[IMAGE OMITTED. SEE PDF]

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_{ij}$, the damping coefficients $D_i$, ${\mu }_i$ and the inertial coefficient $M_i$, ${\sigma }_i$ are different between micro-grids. First, based on the information in [17], a weighted and directed graph has been derived as shown in Fig. 12.

The synchronising coefficients $T_{ij}$ 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_{ij}$ 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 ${\sigma }_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/\pi 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 ${\mu }_i$ are chosen as ${\mu }_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.

Table 1 London City Road grid power capabilities

Number	Name	Capability $G_i$, MVA
1	City Road	1440
2	Devonshire Square	180
3	Beech Street	180
4	Mansell Street	190
5	Hoxton	60
6	Finsbury Market	198
7	Canal Street	132
8	City Road C	202
9	Seacoal Lane	76
10	City Road B	120

[IMAGE OMITTED. SEE PDF]

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.

[IMAGE OMITTED. SEE PDF]

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.

References

[1]

Genis Mendoza, F., Bauso, D., Konstantopoulos, G.: ‘Online pricing via Stackelberg and incentive games in a micro‐grid’. 2019 18th European Control Conf. (ECC), Naples, Italy, 2019, pp. 3520–3525

[2]

Namerikawa, T., Okubo, N., Sato, R. et al.: ‘Real‐time pricing mechanism for electricity market with built‐in incentive for participation’, IEEE Trans. Smart Grid, 2015, 6, pp. 2714–2724

[3]

Milano, F., Dorfler, F., Hug, G. et al.: ‘Foundations and challenges of low‐inertia systems (invited paper)’. 20th Power Systems Computation Conf., PSCC 2018, Dublin, Ireland, 2018

[4]

Dorfler, F., Bolognani, S., Simpson.Porco, J.W. et al.: ‘Distributed control and optimization for autonomous power grids’. 2019 18th European Control Conf., ECC 2019, Naples, Italy, 2019, pp. 2436–2453

[5]

Poolla, B.K., Bolognani, S., Dorfler, F.: ‘Optimal placement of virtual inertia in power grids’, IEEE Trans. Autom. Control, 2017, 62, (12), pp. 6209–6220

[6]

Mendoza, F.G., Bauso, D., Namerikawa, T.: ‘Transient dynamics of heterogeneous micro grids using second order consensus’. Int. Conf. on Wireless Networks and Mobile Communications (WINCOM), Rabat, Morocco, 2017

[7]

Scardovi, L., Sepulchre, R.: ‘Synchronization in networks of identical linear systems’, Automatica, 2009, 45, pp. 2557–2562

[8]

Li, Z., Duan, Z., Chen, G. et al.: ‘Consensus of multiagent systems and synchronization of complex networks: a unified viewpoint’, IEEE Trans. Circuits Syst. I, Regul.Pap., 2009, 57, (1), pp. 213–224

[9]

Bauso, D.: ‘Non‐linear network dynamics for interconnected micro‐grids’, Syst. Control Lett., 2018, 118, pp. 8–15

[10]

Bullo, F.: ‘Lectures on network systems’ (Kindle Direct Publishing, USA, 2016). Available at http://motion.me.ucsb.edu/book‐lns

[11]

Nabavi, S., Chakrabortty, A.: ‘Structured identification of reduced‐order models of power systems in a differential‐algebraic form’, IEEE Trans. Power Syst., 2017, 32, (1), pp. 198–207

[12]

Saber, R.O., Fax, J.A., Murray, R.M.: ‘Consensus and cooperation in multi‐agent networked systems’, Proc. IEEE, 2007, 95, (1), pp. 215–233

[13]

Okawa, Y., Namerikawa, T.: ‘Distributed optimal power management via megawatt trading in real‐time electricity market’, IEEE Trans. Smart Grid, 2017, 8, (6), pp. 3009–3019

[14]

Dörfler, F., Bullo, F.: ‘Synchronization and transient stability in power networks and nonuniform Kuramoto oscillators’, SIAM J. Control Optim., 2012, 50, (3), pp. 1616–1642

[15]

Hiskens, I.A.: ‘Introduction to power grid operation’. Ancillary Services from Flexible Loads CDC'13, Florence, Italy, 2013

[16]

Kothari, D.P., Nagrath, I.J.: ‘Modern power system analysis’ (Tata McGraw‐Hill Education, USA, 2003)

[17]

UK Power Networks: ‘Regional development plan city road city of London (excluding 33 kV)’, 2014. Available at https://library.ukpowernetworks.co.uk/library/en/RIIO/Asset_Management_Documents/Regional_Development_Plans/LPN/LPN_RDP_CityRoad_CityofLondon_33kV.pdf

[18]

Geršgorin, S.: ‘Über die abgrenzung der eigenwerte einer matrix’, Bull. de l'Acad. des Sci. de l'URSS Classe des Sci. Math. et na, 1931, 6, pp. 749–754

[19]

Khalil, H.K.: ‘Non‐linear systems’ (Pearson Education, Prentice‐Hall, USA, 2002)

[20]

Gao, Y.: ‘Existence and uniqueness theorem on uncertain differential equations with local Lipschitz condition’, J. Uncertain Syst., 2012, 6, pp. 223–232