Supermodal Decomposition of the Linear Swing Equation for Multilayer Networks.

We study the swing equation in the case of a multilayer network in which generators and motors are modeled differently; namely, the model for each generator is given by second order dynamics and the model for each motor is given by first order dynamics. We also remove the commonly used assumption of equal damping coefficients in the second order dynamics. Under these general conditions, we are able to obtain a decomposition of the linear swing equation into independent modes describing the propagation of small perturbations. In the process, we identify symmetries affecting the structure and dynamics of the multilayer network and derive an essential model based on a 'quotient network.' We then compare the dynamics of the full network and that of the quotient network and obtain a modal decomposition of the error dynamics. We also provide a method to quantify the steady-state error and the maximum overshoot error. Two case studies are presented to illustrate application of our method.

INTRODUCTION

Modeling and analysis of power grid dynamics have been the subject of many papers [1], [9], [12], [13], [15], [16], [20], [21], [26], [29], [34], [35], [38], [39], [41], [48]. A fundamental tool to describe the dynamics of power grids is the swing equation. In the presence of small disturbances, this equation can be approximated by the linear swing equation, which can be decomposed into independent modes [2], [5], [19], [24], [37], [46]. However, most papers in the literature introduce two unrealistic assumptions in order to derive a modal decomposition: 1- the same model is used for both generators and motors, while they are intrinsically different and 2- all the individual systems are characterized by the same effective damping. An exception to 2 is provided by Tyloo et al. [46] who relax the constant inertia to damping ratio assumption and show that their derivation is still valid with heterogeneous dynamical parameters. Here we remove both assumptions 1 and 2. We introduce intrinsically different mathematical models for generators and motors, where motors have typically negligible inertia compared to generators and we allow for the damping terms to be arbitrarily chosen. After removing both assumptions 1 and 2 we show that a modal decomposition of the linear swing equation can still be obtained.

While performing this analysis we also look for the effects of symmetries in the network topology. Symmetries play a significant role in the study of networked systems. References [3], [6], [11], [14], [17], [31], [33], [36], [43]–[45], [49] have proposed tools based on graph theory and group theory to analyze the dynamics of complex networks with symmetries. A recent paper [10] has proposed indices for the characterization of symmetries in complex networks. Reference [7] has discussed how the analysis of symmetries may be used to explore associations in image features to infer similarities/dissimilarities among features.

The presence of symmetries in power grid networks has been documented in [14], [45]. Reference [22] has analyzed how network symmetries may affect the synchronization modes of power grid networks, while [18] has studied structure and dynamics preserving reduction methods for power grids using a graph theoretical approach. Reference [5] has analyzed the effect of symmetries in networks using the linear swing equation and developed techniques for easy calculation of maximum flow and error using a modal decomposition. Reference [27] has found that certain heterogeneous parameters enhance the stability of power grid networks, which provides a motivation for us to remove the assumption of homogeneous parameters. Many real systems do not possess exact symmetries but approximate symmetries. Recent work has investigated the presence of such approximate symmetries in both natural and engineering systems and found that approximate symmetries are widespread in these systems, see e.g., [28], [30]. Reference [42] has shown a connection between the stability of symmetric states and approximately symmetric states.

Different from large part of the literature, in this paper we model powergrids as multilayer networks, where one layer is entirely made up of generator nodes and the other layer of motor nodes (‘loads’). We first obtain a simplified description of the multilayer network dynamics in terms of its quotient network. This description is useful as it provides an essential model for the dynamics in which redundancies due to symmetry are ‘removed’ [44], under the assumption that symmetric nodes produce/absorb the same amount of power. In what follows, we withdraw this assumption and quantify the error/deviation dynamics between the quotient network model and the full network model. Further, we develop a method to compute the maximum overshoot of the aforementioned error dynamics.

The are several novel contributions of this paper, which we briefly review next. They are: (i) a study of the linear swing equation applied to a multilayer network in which loads and generators are modeled differently, (ii) a dynamical reduction of the dynamics in the case of heterogeneous damping coefficient based on the solution of the quadratic eigenvalue problem, (iii) an expansion of the solution in terms of first order supermodes and second-order supermodes, and (iv) a characterization of the error dynamics of the quotient network when compared to the full-network in terms of supermodes.

The rest of the paper is organized as follows. Section II covers the modeling of a powergrid as a multilayer network using the linear swing equation and block diagonalization of the network equation into its quotient block and transverse block using the IRR transformation. Section III describes the block-diagonalization of the error dynamics. Section IV deals with the maximum error computation using the concept of a ‘supermode’ and Section V provides validation using one example network and one real-life power grid network. Finally, conclusions are given in Section VI.

NETWORK MODEL

Several alternative models of the swing equation have been proposed [25], [32]. For this paper, we specialized Bergen and Hill’s Structure Preserving Model [4], [32] to multilayer networks. The power grid is formed of a collection of rotating machines (‘nodes’) connected by transmission lines (‘edges’). Nodes can be of either one of two types: generator nodes which produce power, and motor nodes (or loads) which consume power and have typically a lower inertia than generator nodes.

We thus consider a multilayer network with two layers: the generator layer denoted with n nodes and the load layer denoted with n nodes. The nodes in are modeled as second order oscillators, with dynamics described by the following equation, where θ(t) is the nodal displacement of generator node i, γ > 0 is the damping coefficient of generator i and is the power produced at generator node i. The nodes in are modeled as first order oscillators which obey the following equation, where ϕ(t) is the nodal displacement of motor node j, and is the power consumed at load j.

Definition 1: The network is said to be balanced if .

We assume the network is balanced throughout the paper unless specified otherwise. The multilayer network can be represented using the ‘Supra-Adjacency matrix’ which is an n × n matrix, n = n + n, where the n × n-dimensional symmetric matrix (the n × n-dimensional symmetric matrix ) describes the intra-layer connectivity of layer Namely, if generator nodes i and j are connected and otherwise. Analogously, if motor nodes i and j are connected and otherwise. The interlayer connectivity is given by the n × n-dimensional matrix and by the the n × n-dimensional matrix , . If generator node i and motor node j are connected (not connected), then

Like in the case of single-layered networks, multilayer networks may be affected by the presence of symmetries [14]. In general, a symmetry is a permutation of the network nodes that results in a network that is isomorphic to the original one. The symmetry group is the set of all symmetries with the operation composition. The set of all symmetries in the group will only permute certain subsets of nodes (the orbits or clusters) among each other. The nodes in each subset are mapped into each other by application of one or more symmetries in ; however, there is no symmetry in that will map into each other nodes in different subsets. We refer to such subsets of nodes as ‘clusters’ or ‘orbits’ of the symmetry group [5].

For the case of multilayer networks, with layers composed of nodes of different types, the above definition of symmetry needs to be specialized to account for the fact that symmetries may only move nodes (may not move nodes) from the same layer (from different layers.)

We first introduce the group of symmetries for each layer: for layer which satisfies, and for layer , which satisfies,

Definition 2: A symmetry for the set of (1) and (2) is defined as a block diagonal permutation matrix of the form, with , , and such that and s = , where the n-dimensional vector .

In addition to the equalities (3) and (4), the equation also requires satisfaction of the following two conjugacy relationships [14],

By using (5a) and (5b), we can define the following two subgroups of and [14], Following [14] we can prove that the two sets and are subgroups of and , respectively.

The symmetry group partitions the set of nodes of the generator layer into a set of l clusters, , k = 1, …, l and the set of nodes of the load layer into a set of l clusters, , k = 1, …, l. For each layer, each cluster is formed of nodes that are mapped into each other by application of all the symmetries in ; however there is no symmetry in that maps into each other nodes in different clusters. We call l = l + l the total number of clusters. In what follows we call i* the cluster to which generator i belongs, and j* the cluster to which motor j belongs.

Lemma 1: A flow invariant synchronous solution , is induced by the automorphism group .

Proof: Assume and for all i’s in cluster i* and for all clusters i* = 1, …, l and for all j’s in cluster j* and for all clusters j* = 1, …, l. It follows that is the same for all i’s in cluster i* and is the same for all j’s in cluster j*. □

Next we assume convergence of the flow invariant solution on the fixed points, for all i’s in cluster i*, and for all j’s in cluster j*. The linearized swing equation, which models the propagation of small disturbances (e.g., affecting the initial condition or affecting the power supplied/demanded at different nodes) [46], is obtained by linearizing (1) and (2) about the stable fixed point , i = 1, … n, , j = 1, …, n, where , , , . The Laplacian matrices with entries and with entries , where δ is the Kronecker delta. Also, represents the total connectivity between generator i and the layer and represents the total connectivity between load j and the layer. Each term and on the right hand side of (6) and (7) represents a small power deviation.

Definition 3: A symmetry for the set of (6) and (7) is defined as a block diagonal permutation matrix of the form, with , , and such that , where the matrix, and s = .

Lemma 2: The set of Eqs. (1) and (2) have the same set of symmetries as Eqs. (6) and (7).

Proof: We need to prove that (a) a matrix that commutes with , commutes also with and (b) viceversa. We first prove (a). The matrix P permutes with one another nodes in layer and permutes with one another nodes in layer . For all v, w = 1, …, n, , where v′(w′) is the node that v(w) is mapped to by P. Now as v is mapped to v′, then v* = v′*; and as w is mapped to w′, then w* = w′*. From this it trivially follows that: (i) if v, , then , (ii) if v, , then and (iii) if and , then . Thus for all v, w = 1, …, n, , where v′(w′) is the node that v(w) is mapped to by P. The proof for (b) is analogous. We start from the observation that where v′(w′) is the node v(w) is mapped to by P, and then by using properties (i), (ii), and (iii) we can prove that for all v and w, . □

We now consider the vectors , , , to rewrite (6) and (7) as, where the diagonal matrix Г = {Г} with diagonal entries Г = γ is the damping matrix for the generators, D is the degree matrix of the generators with respect to the loads and D is the degree matrix of the loads with respect to the generators. D is defined as the n × n diagonal matrix where , if i = j and 0 otherwise. D is defined similarly.

We introduce the vector X = [, ] and rewrite (8) and (9) in the compact form, where , , and .

From knowledge of the group of symmetries , we can compute the irreducible representations (IRRs) of the symmetry group of the network. This defines a transformation T into the so-called IRR coordinate system (see Ref. [36]). The transformation matrix T is orthogonal. Each one of the rows of the matrix T is associated with a specific layer. If a row of the matrix T is associated with layer k, all the i entries of that row are zero for i not in layer S.

Therefore, the matrix T can be cast in the following block diagonal form, where each block T is an n-dimensional square matrix associated with layer k and the symbol ⊕ indicates direct sum of matrices. We will represent these blocks as T for the generator layer and T for the load layer giving a matrix T of the form: Premultiplying (8) and (9) by T, we obtain, and where = T, = T, r = Tp, r = Tp, V = TL(T), V = TL(T), , . Due to the specific diagonal structure of the matrices D, D and Г, they remain unchanged after the transformation.

Remark 1: Both matrices T and T can in turn be written as block-diagonal matrices with the number of blocks equal to the number of clusters in either layer. For our case, T is a matrix with l blocks and T is a matrix with l blocks such that T can be viewed as a block-diagonal matrix with a total of l blocks.

If we define the vector Ω = [, ], (12) and (13) can be recast as, which is in the familiar form,

It is to be noted that the n × n matrix , is block diagonal and is equal to the direct sum , where K is a (generally complex) p × p matrix with p the multiplicity of the u IRR in the permutation representation, U the number of IRRs present and d the dimension of the u IRR, so that ∑ dp = n. The matrix T contains information on which perturbations affecting different clusters get mapped to different IRRs [40]. Due to the block-diagonal structure of the matrix , (14) can be decoupled into a number of lower-dimensional equations, where each equation corresponds to a block of the matrix . There is one representation (labeled u = 1) which we call trivial and the associated block of the matrix corresponds to the dynamics of the quotient network.

Definition 4 (Indicator Matrix): The n × l indicator matrix E has entries E(i, j) = 1 if node i belongs to the cluster j and 0 otherwise,

Definition 5 (Quotient Network): A Quotient Network is a reduced network where redundancies due to symmetries are removed. The dynamics of the quotient network is obtained by pre-multiplying (10) by ((E E)−1E), yielding, where the n × 1 vectors , and the matrices

Hence, the trivial representation is associated with all the clusters. However, it is possible that other IRR representations are only associated with some of the clusters (not all of them.) Each one of these other representations u > 1 describes the deviation/error dynamics of either an isolated cluster or a group of intertwined clusters [36]. A simple interpretation of isolated vs. intertwined clusters is the following. If a cluster is isolated, a perturbation affecting the power of any one of its nodes will not affect the deviation dynamics of other clusters. On the contrary, when a set of two or more clusters are intertwined, a perturbation affecting the power of any of the nodes in a cluster will affect the deviation dynamics of the remaining clusters in the set.

DIAGONALIZATION

The quotient network provides a minimal representation of the full network. Our goal here is to describe the error dynamics between the full network dynamics and the quotient network dynamics. where X(t) is the dynamics of node i which belongs to cluster j of the entire network and is the dynamics of node j of the quotient network, which corresponds to an average of the dynamics of all nodes in cluster j.

To quantify this error, we present an approach to diagonalize the blocks of the matrix with u > 1 and size m > 1. We consider an m-dimensional transverse block from (15), with dynamics, where is an m × m diagonal block of the matrix, is the corresponding m × m block of the matrix and is the corresponding m × m block of the matrix. and are the corresponding m × 1 vectors from the Ω and vectors, respectively. In what follows, we will distinguish between three possible cases: 1) The block represents the error dynamics of intertwined clusters within the generator layer. This means that the matrix will be a diagonal matrix with all non-zero elements in the main diagonal and the matrix can either be homogeneous, i.e, it is a multiple of the identity matrix, or non-homogeneous. This is equivalent to a single layer network with either same or different damping coefficients which has already been studied in [5]. 2) The block represents the error dynamics of intertwined clusters within the load layer. This means that the matrix is a zero matrix and the matrix is the identity matrix. This is a simple case that allows a trivial decomposition. In this paper, we also consider a third, more complex case; 3) The block represents the error dynamics of intertwined clusters in different layers. This means that the matrix is a diagonal matrix with both zero and non-zero elements along the main diagonal, i.e., it is a singular matrix different from the zero matrix. The matrix can be non-homogeneous. To diagonalize such a system, we use the method described in [23]. The advantage of this method is that it is a direct generalization of the modal decoupling for the case that the matrix is invertible, i.e, case 1. Should the system have an invertible mass matrix and be classically damped, i.e, satisfies the commutativity relationship, [8], [23], this method reduces to classical modal analysis. The method is summarized below. We consider the associated quadratic eigenvalue problem (QEP):

The solution of this QEP yields a total of 2m eigenvalues, of which σ = m + r are finite and ϵ = m − r are infinite, where r is the rank of the matrix and m is its size. We randomly assign 2r finite eigenvalues to conjugate pairs: λ and where i = 1, 2, 3, …, r, and we denote the remaining ϵ unpaired eigenvalues, called the “lone eigenvalues” by μ [23]. We now construct the matrices Λ and which are diagonal matrices of dimension r × r and the matrix Ξ which is of dimension ϵ × ϵ. We also construct the matrices V, and W whose columns are the eigenvectors of the matrices Λ, , and Ξ, respectively: v, , i = 1, 2, 3, …, r and w, i = 1, 2, 3, …, ϵ. where here and in what follows, C = [A|B] symbolizes concatenation of two vectors/matrices with the same number of rows. We can now construct the diagonalized matrices: In order to find consistent initial conditions, we also need to define, Multiplying (19) by β2 = 1/λ2 yields: The ϵ infinite eigenvalues from (19) correspond to ϵ zero eigenvalues in (20), β = 0. These ϵ zero eigenvalues and their corresponding eigenvectors must satisfy: where J is an order ϵ Jordan matrix of eigenvalues β = 0, and V is is a n × ϵ matrix of the associated eigenvectors. Since (18) is nondefective, J = 0, which implies . Therefore, V is the matrix of eigenvectors in the null space of . Now, the transformation matrix S can be defined as, where S is a 2m-dimensional real orthogonal matrix, and V is the matrix of eigenvectors in the null space of A2 (analogous to the case of and V). If system (18) is written in the symmetric state space realization, the transformation S will produce the diagonal matrices A0, A1, and A2,

We now derive the two transformation matrices: Using T1 and T2, we can calculate the transformed forcing function,

The diagonalized equation in coordinates can be found using the A0, A1, A2 and g we have already calculated: The initial conditions in the coordinates also need to be calculated from consistent initial conditions in the coordinates, where the σ × m matrices Z, Z are constructed as described in [23]. In what follows, without loss of generality we focus on the forced response and set the initial conditions in the X-coordinates to zero. Hence, we have zero initial conditions in the η-coordinates as well.

Lemma 3: Consider zero initial conditions in X-coordinates. For unit step forcing , only the initial velocity for the second order modes are non-zero; all other initial conditions are zero, that is, if X(0) = 0, , and , where is the vector space for all unit step functions, then,

Proof: We know that if , f (0) = 0 and . Also, we have 0 initial conditions in the X and η coordinates. So, we can calculate our modal initial conditions as:

It is now obvious that χ(0) = 0 ∀ i. Since the non zero term of (29) it follows that when (from (26)). From (25), we can see that the last ϵ columns of T2 are zeros which means that the last ϵ rows of are zeros. Therefore,

After pre-multiplying the expression by V Z, we get the modal initial velocities of the system. This results in the last ϵ elements of to be zero which corresponds to the first order modes, thereby completing the proof. □

Note that we can use the transformation matrices to recover the displacements in η coordinates, We see that (27) can be broken into m independent equations, the solutions of which are called the ‘modes’ of the system. i = 1, …, m. The linear combination of all of the modes, their derivatives, and the forcing function provides the solution for the error as shown in (30).

After diagonalization, the first r modes we obtain are second order modes, (A2) ≠ 0, i = 1, …, r, and the remaining m − r modes are first order modes, (A2) = 0, i = r + 1, …, r + n. In realistic systems with constant damping and inertia, it has been found that all modes are underdamped and propagate through the whole system with , ∀i > 1 [46], [47]. Reference [5] has shown that this assumption holds true even after the assumption of constant damping is removed. We can then characterize the modal responses in terms of well-known first order time responses and underdamped second order time responses to unit step forcing.

Each second order mode can be represented as, where and ζ = (A1)/(2ω). The solution to (32) can be written as, where is the free evolution and is the forced evolution. is given by: where, and . (B1) and (B2) are constants that depend on the initial condition of the modes and can be calculated as follows: We know from Lemma 3 that only are nonzero for second order modes. Therefore we have, We can use this into (34) and obtain, The underdamped forced solution is given by, i = 1, 2, …, r. The full response is the sum of (36) and (35).

On the other hand, each first order mode can be represented as, where, and . The solution is equal to,

LINEAR COMBINATION OF MODES FOR MAXIMUM ERROR QUANTIFICATION

As stated in the previous section, for a unit step forcing, , the linear combination of the modes, their derivatives, and the forcing function can be used to calculate the deviations using (30), Here to ease this calculation, we introduce the concept of a ‘supermode.’

Definition 6: A supermode is a linear combination of modal displacement and velocity. Supermodes can be of two types: First order supermode and second order supermode.

Each supermode can be expressed by the equation, where κ is a real constant. We note that from (30) can be obtained as a linear combination of the , i, j = 1, 2, …, m, from (39). Our definition of supermode in (39) is general, as κ can be seen as a variable parameter that determines each supermode. In what follows we will see how first order supermodes and second order supermodes can be parameterized in a minimal number of parameters.

FIRST ORDER SUPERMODES

First order supermodes, i = r + 1, r + 2, …, m, are obtained by plugging (38) in (39), yielding, which converges to the steady state,

It can be seen from (40) that a first order supermode is completely parameterized by the constant (κ), modal forcing () and time constant (τ).

SECOND-ORDER SUPERMODES

Second order supermodes, i = 1, 2, 3 …, r, are obtained by plugging (36) in (39). We take the first derivative of the solution to get, where . So, we can simplify (39) as follows, where, This converges to a steady state The peak time can be calculated by setting , therefore, all the local extrema of the supermode are at, which can be plugged into (42) to find the associated supermode peak value, . Since there is damping in the system, one may expect the first peak to always be the global maximum. However, it is also possible that this local maximum point is preceded by a local minimum point. So it becomes necessary to check (45) for k = 0, 1, to determine the supermodal peak.

By calculating , an analytical condition can be found to determine the supermodal peak. If is positive (negative), then the first (second) peak is the supermodal peak. If instead , the curvature of the supermode at t = 0 should be assessed, Since ,

Therefore, if the supermode is convex at , the next peak corresponding to k = 1 is the global peak. If the supermode is concave at , the peak corresponding to k = 2 is the global peak. In summary, the global peak is found at, where and are defined in (46) and (48), respectively.

Equation (42) shows that a second order supermode can be completely parameterized by its constant (κ), natural frequency (ω), damped frequency (ζ), modal forcing (g), and initial condition (χ(0)).

LINEAR COMBINATION OF SUPERMODES

Each η can be written as a linear combination of supermodes and of the forcing function, , in the coefficients and which denote how much a certain supermode affects the error dynamics. We note that κ from (39) and , and from (50) are variable constants that can be arbitrarily assigned. These constants can be chosen based on the specific purpose of the analysis. For example, one can pick κ = T2(i, j)/T1(i, j), , and , where i, j = 1, …, m, to obtain from (30).

Since the expressions for steady state values of the supermodes have been derived, we can calculate the steady state error as,

For peak times, we use a slightly modified version of the technique described in [5]. First we take the sum of the supermodal peak values and non-transformed forcing functions for all supermodes. We call this ηL, For first order supermodes, the peak is given by , i = r + 1, …, m, and the peak time can be approximated with the settling time equal to . The ηLs calculated from (52) are then cumulatively added in ascending order of peak time and the peak time corresponding to the largest of these cumulatively added values is then used as the initial guess in the equation: which can be expanded into: where, Solving (53) numerically using our calculated initial time provides the time for the peak error . It is also important to note that the closer the peak times are, the more accurate this initial guess is and the farther apart they are, the more iterations will be required to converge to the actual solution. Once the peak time is known, it can be plug back into (50) to find the peak error. Finally, the maximum error could either be the peak error (if the second order supermode dominates) or the steady state error (if the first order supermode dominates). Thus, we take the largest value between and as the maximum error ,

CASE STUDIES

In this section we provide 2 examples to demonstrate our method: Example 1 shows how our method can be utilized to find the maximum deviation error for a small network and example 2 demonstrates application of the method to the case of the IEEE 145-bus test grid.

EXAMPLE 1

Our first example is the network in Figure 1 with n = 4 nodes in the generator layer and n = 2 nodes in the load layer, with p = [0.1 0.1 0.1 0.2] and p = [−0.2 −0.3].

FIGURE 1.

(a): Full network, (b): Quotient network. Generators are represented by green squares and loads are represented as red circles.

The network dynamics is given by the following equation:

Using the IRR transformation matrix T, (56) is transformed into,

The 2 × 2 diagonal block on the top left of the matrix of (57) corresponds to the quotient network dynamics. The quotient network is depicted in Figure 1(b). The other diagonal blocks represent the error dynamics between the full network and the quotient network, which is what we are interested in. The rightmost 2 blocks are simple so we diagonalize the second 2 × 2 diagonal block from the left to demonstrate our method. The dynamics of this one block is,

The 3 finite eigenvalues for this 2 × 2 system are, λ1 = −0.5193 + 1.5232i, λ2 = −0.5193 − 1.5232i and λ3 = −3.8615. Pairing λ1 and λ2 as conjugates and λ3 as the lone eigenvalue, we can find the normalized eigenvectors,

We can now compute the matrices Λ, , V, , Ξ and W,

We then solve for T1, T2, and T3,

We then calculate the decoupled coefficient matrices A2, A1, A0, transformed power g, and initial condition (0), resulting in the following decoupled equation,

Equation (59) describes the dynamics of two supermodes, a second order supermode (i = 1) and a first order supermode (i = 2). To calcuate the supermodal peak for each, we extract some parameters from the diagonalized system: ζ1 = 0.3227, ω1 = 1.6093, , .

Table 1 provides information on the supermodal peak time and associated peak calculated via our method and compares it to the actual peak time and peak obtained by solving the linear ODE. We also cumulatively add the peaks to calculate the initial guess. Table 2 provides information on the maximum error values for our system calculated using our approach and compares it to the error values obtained by directly solving the linear ODE (rightmost column of the Table.) This is also consistent with the plot of the error dynamics shown in Fig. 2.

Table 1.

Initial guess calculation using a linear combination of supermodes for η1. Supermodal peak time and peak for found by solving the linear ODE are in columns 2 and 3 from the left and the ones calculated using our approach are in columns 4 and 5. The table is arranged in ascending order of peak time and column 6 is the cumulative peak added in that order. The bold peak time represents the initial guess to solve (53).

η index	ODE Peak Time	ODE Peak	Calculated Peak Time	Calculated Peak	Cumulative Peak	Supermode Index
1	1.0359	0.0209	1.0359	0.0209	0.0209	2
	2.4354	0.0105	2.4359	0.0105	0.0314	1
2	1.0359	0.0021	1.0359	0.0020	0.0020	2
	2.1654	0.0377	2.1651	0.0377	0.0397	1

Table 2.

Steady-state error and peak error calculation using an ODE solver in columns 2 and 3 whereas the same calculation using our approach is in columns 4 and 5. The larger of these values is the maximum shown in column 6.

η index	ODE Steady State	ODE Peak	Calculated Steady State	Calculated Peak	Maximum Error
1	0.0283	0.0314	0.0283	0.0314	0.0314
2	0.0300	0.0397	0.0300	0.0397	0.0397

FIGURE 2.

Error vs. time for the network in figure 1. Each curve represents the error dynamics due to the power not respecting the symmetries in the two layers of the network.

Using the IRR transformation matrix we find how this error corresponds to the dynamics of the full network and the quotient network: and where and are the dynamics of the two nodes of the associated quotient network.

EXAMPLE 2

Our second example is an application to the IEEE145 bus network with n = 50 generator nodes, n = 95 load nodes, and 422 transmission lines [50]. The network is shown in Fig. 3a and the corresponding quotient network in Fig. 3b. We assume that the power produced by the generator nodes is equal to and the power absorbed by the load nodes is equal .

FIGURE 3.

(a) Full network representation of IEEE145 test grid. (b) Quotient network representation of IEEE145 test grid. Circles indicate load nodes/clusters and squares indicate generator nodes/clusters. Nodes belonging to the same clusters are colored the same in (a) and the same colors are used for the clusters shown in (b). Nodes colored black belong to a trivial cluster.

Using the IRR transformation, we obtain twentyfour transverse blocks, of which nineteen are scalar, three are 2 × 2, one is 3 × 3 and one is 4 × 4. These blocks represent the deviation dynamics of the individual nodes compared to their clusters. To demonstrate our method, we choose the 3 × 3 blocks, which corresponds to the error dynamics of three intertwined clusters, with two nodes each. The first two clusters are formed of load nodes {121, 122} and {76, 81}, while the third cluster is formed of generator nodes {21, 22}. These three clusters can be seen in Figure 3, where their nodes are represented as red circles, green circles and blue squares on the top right of Fig. 3a, respectively. They map to similarly colored and shaped nodes of the quotient network, all shown in the bottom left of Fig. 3b.

In order to characterize the error dynamics, we increase the power of node 76 to −0.2 to induce an asymmetry in power within that cluster. We also increase the power of node 22 to 0.5053 in order to maintain the network balance. Since nodes in the same clusters do not have the same power, it is expected that the quotient network dynamics will not represent the full network dynamics, which results in the aforementioned error dynamics. As in the previous examples, we are interested in characterizing the maximum overshoot of the error dynamics. Our chosen 3 × 3 block is associated with the following system of equations: Using our method, we diagonalize this system, We also obtain the transformation matrices and initial condition:

Using parameters taken from the diagonalized system and applying the same method as in the previous example, we can calculate the initial guess for each η, which is shown for η1 in Table 3. This initial guess is then used to compute the maximum errors shown in Table 4. The values in the tables are consistent with the error dynamics shown in Fig. 4.

Table 3.

Initial guess calculation using a linear combination of supermodes for η1. Supermodal peak time and peak for found by solving the linear ODE are in columns 1 and 2 from the left and the ones calculated using our approach are in columns 3 and 4. The table is arranged in ascending order of peak time and column 5 is the cumulative peak added in that order. The bold peak time represents the initial guess to solve (53).

ODE Peak Time	ODE Peak	Calculated Peak Time	Calculated Peak	Cumulative Peak	Supermode Index
0.2099	0	0.2099	0	0	3
2.0616	0.0022	2.0616	0.0022	0.0022	2
3.7260	0.0146	3.7255	0.0146	0.0168	1

Table 4.

Calculation of steady-state and peak for all 3 ηs. In columns 2 and 3 we have the steady-state error and peak error calculated using an ODE solver, whereas in columns 4 and 5 we report values calculated using our approach. Column 6 shows the maximum value between the calculated steady-state and peak errors.

η index	ODE Steady State	ODE Peak	Calculated Steady State	Calculated Peak	Maximum Error
1	0.0142	0.0166	0.0142	0.0166	0.0166
2	0.0406	0.0416	0.0406	0.0416	0.0416
3	0.02301	0.2742	0.2301	0.2742	0.2742

FIGURE 4.

Error vs. time for the network in figure 3. Each curve represents the error dynamics in the 3 intertwined clusters presented in example 3.

Like in the previous examples, the ηs can be expressed in terms of displacements of nodes of the full network and the quotient network. , , , where represents the displacement of the i quotient node (the average of the displacements of the individual nodes in that cluster), and X is the displacement of node i of the full network.

CONCLUSION

In this paper, we have specialized the structure-preserving model of the swing equation to a multi-layer network with two layers, one formed of generator nodes and one formed of motor nodes. In the presence of symmetries, a lower dimensional model for the dynamics is provided by the so-called quotient network. However, such representation is often inexact when nodes in the same cluster generate or consume different amounts of power, which makes it important to characterize how much the actual network dynamics deviates from that of the quotient network.

A main difference with large part of the literature is that we have relaxed the commonly used assumption of homogeneous damping coefficients and considered the general case that these can be different from node to node. By using a combination of group representation theory and the solution of the quadratic eigenvalue problem, we have reduced the transient analysis of the linear swing equation in terms of a set of independent first order and second order supermodes. Each supermode is fully defined by a minimal set of parameters. This has led to a characterization of the transient error dynamics for the individual nodes of the multilayer network and the development of a method to compute the maximum value of the transient error. We have presented application of this method to two examples of interest, one of which is the IEEE 145-bus test grid, and demonstrated how it can be used to characterize the deviation between the dynamics of the full network and that of the reduced quotient network.