Clustering approach to model order reduction of power networks with distributed controllers

This paper considers the network structure preserving model reduction of power networks with distributed controllers. The studied system and controller are modeled as second-order and first-order ordinary differential equations, which are coupled to a closed-loop model for analyzing the dissimilarities of the power units. By transfer functions, we characterize the behavior of each node (generator or load) in the power network and define a novel notion of dissimilarity between two nodes by the ℋ2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal {H}_{2}$\end{document}-norm of the transfer function deviation. Then, the reduction methodology is developed based on separately clustering the generators and loads according to their behavior dissimilarities. The characteristic matrix of the resulting clustering is adopted for the Galerkin projection to derive explicit reduced-order power models and controllers. Finally, we illustrate the proposed method by the IEEE 30-bus system example.


Introduction
Power networks, nowadays, are experiencing the penetration and integration of a wide array of new electronic devices and renewable energy sources (see [6,16] for an overview).In the foreseeable future, power networks will become more and more complex with more variation in the generators, uncertain loads, and denser transmission lines.Centralized power generation is being replaced by a more distributed generation.
The immense size of power grids yields mathematical models with high dimensions complicating further analysis.In most cases, a complete model of the power network is neither practical nor necessary for, e.g., transient analysis, failure detection, distributed controller design, or system simulation.Therefore, we need to construct a reduced-order model that can approximate the behavior of the original complex power system with an acceptable accuracy.More importantly, we desire to preserve the network structure in the reduced-order model such that it can be interpreted as a reduced network with a less complex topology.Specifically, the reduced power system is evolving over a simpler network that consists of fewer generator and load buses and sparser transmission lines.The network structure is necessary for the application of sensor allocation and the management of distributed power generation.The main interest of this paper is to investigate the problem of model order reduction for power networks with distributed controllers with the preservation of network structures.
Literature review Conventional model reduction techniques, including balanced truncation and Krylov subspace methods, have been extended to the dynamic reduction of power systems (see, e.g., [10,23,24,32,42,45,49]).In these papers, the power network systems are modeled in first-order state space representations, and the reduced-order models are constructed within the framework of Petrov-Galerkin projection.Even though these methods provide systematic procedures to produce lower-dimensional models, which benefits the computation and simulation of power systems, the network structure is not retained through the reductions, and the projected states in the reduced-order model often lack a physical interpretation.A structure-preserving approach is proposed in [15], where the network model is firstly reduced by the balanced truncation, which still preserves the semistability such that the reduced-order model can be reformed as a network system.In spite of fewer nodes, the obtained reduced network has a complete interconnection.This restricts the use of the reduced-order models in further applications, especially the implementation of distributed control laws.In a distributed control network, control laws are implemented locally by each generator, and these controllers send and receive data through a communication network, which is only connecting generators and in general, that has a different topology from the power network (see [1,21,31,34,36,47] for more details).In this case, each generator can only receive information from its adjacent generators through the communication network, which means that the control signals have to be generated by using only part of the system states.However, in the reduced-order model resulting from the conventional approaches, each state is composed of a linear combination of all the original states.As a result, the controllers designed based on such models are not to be realized in a distributed fashion.
Another network simplification method is the Kron reduction, which has been widely used in e.g.classic circuit theory, smart grid monitoring, and transient stability assessment.The extension of this approach for model order reduction of electric networks can be found in [11,19,36] and in the references therein.The network structure of a power system is preserved as the Schur complement of the Laplacian matrix of the original network is again a Laplacian matrix, which represents a smaller network.Despite the simplicity of implementation, the Kron reduction modeling only eliminates the load buses in the power network [19].To reduce the complexity of a network of synchronous generators, which is modeled by a second-order swing equation, [20,41] propose a method based on time-scale separation and singular perturbation analysis.This approach identifies the sparsely and densely connected areas of power grids and then aggregates the state variables of the coherent areas.Similar to Kron reduction, the algebraic structure of the Laplacian matrix is maintained through the singular perturbation approximation.However, the method in [41] ignores the fluctuation of power demands from the loads.
Clustering analysis of power networks, derived from the concept of diakoptics, has been intensively explored in the study of coherency recognition, area aggregation, and behavior approximation in the power grid.A large and complex power network could be more easily analyzed and managed when it is decomposed into several smaller components.The coherency-based approach is most popular for the clustering analysis of electric networks.Many related results have been reported in the literature (see [3,16,17,26,33,39,48] and the references therein).The generator coherency is introduced as the tendency of generators when the disturbances are affecting the system.The reduced model of the power network then results from aggregating those coherent generators.However, the coherency identification is a data-driven process, which heavily relies on the accuracy of the sampling data generated from the time-domain voltage angle responses of the generator buses.
Differently, [43] proposes a Petrov-Galerkin approximation framework based on graph clustering, which is applied to physical network systems.However, the methods for cluster selection are still not well-explored.As an extension, [37] suggests to use the almost equitable partition (AEP) as a clustering of the underlying graph, but finding AEPs is generally a very difficult problem and computationally expensive, which limits this method in practical applications.The results in [9,35] can be regarded as the combination of graph clustering and conventional model reduction techniques, balanced truncation, and Krylov subspace methods, respectively.However, it is still a challenge to apply these ideas to second-order systems, which are common model settings in the context of power networks.The approach in [28] offers another feasible solution for the clustering-based model reduction.The notion of cluster reducibility is generalized as the uncontrollability of the system states.Then, the reducible states are merged to reduce the dimension of the system as well as to preserve the network structure.To cope with second-order networks, the results are then extended in [27].Nevertheless, the considered model is a second-order model, which is assumed to be asymptotically stable.This restricts this approach in some applications, such as the coupled swing dynamics in power networks [20].

Contribution
The presented work provides a novel model order reduction scheme for power networks with the associated distributed controllers.The reduction procedure is developed based on our previous results in [12,14], where the reduction of first-order network systems is studied.An extension to second-order models of networks is proposed in [13].This paper explores the application of power networks that contain many generator and load buses, and the proposed method aims to simultaneously simplify the structures of the transmission network and its distributed controller.
By linearizing the nonlinear structure preserving model of the power network [7], we obtain the network model in a second-order representation.Together with the distributed controller, the closed-loop system is derived, which is semistable and contains the Laplacian matrices of the power transmission network and the communication graph of the controllers.The clustering analysis of the networks is based on the closed-loop system.Similar to the idea in [13], that a reduced-order model is generated by clustering the nodes behaving similarly.More specifically, the behaviors of network nodes are characterized by transfer functions, and an H 2 -norm characterization of the dissimilarity of nodal behaviors is adopted such that the hierarchical clustering algorithm in [13] can find suitable clusterings for generator and load nodes, respectively.This notion of dissimilarity is data-independent, which means that we do not need the time-domain sample data to evaluate the differences among the nodes.Compared with [13], the current paper mainly focuses on the application of power networks whose models couple the generator and load dynamics with the distributed controllers.The proposed method need to reduce the power network and communication network simultaneously, instead of only simplifying the network system as in [13].The definition of nodal behaviors cannot be directly applied, since the influence of communication networks of the distributed controllers has to be concerned, and the behaviors of generator and load nodes are supposed to be treated differently.
As the dissimilarity is defined in terms of the H 2 -norm, we use the controllability Gramian of the closed-loop system to evaluate the dissimilarities between all pairs of generator and load buses.A numerical difficulty comes from the computation of the controllability Gramian.Due to the semistability, the corresponding algebraic Lyapunov equation may not have a unique solution [5].This paper extends the idea of [14] to obtain a controllability Gramian by solving two algebraic equations rather than a single Lyapunov equation.
Besides, the hierarchical clustering algorithm is adapted to generate appropriate network clusterings for both generators and loads, so that the nodes with similar behaviors are grouped into the same clusters.Then, the Galerkin method is applied to yield a reduced-order power network with a simplified distributed controller.The characteristic matrices of the resulting clusterings are used as the projection matrices.It is shown that the algebraic structures of Laplacian matrices are retained in the reduced-order model.Consequently, the interconnection topologies of the power transmission network and communication network of the controller are simplified.

Models of power network and distributed controller
Consider a connected power network consisting of n synchronous generators and m loads.Denote V g = {1, • • • , n} and V l = {n + 1, • • • , n + m} as the index sets of generator and load buses, respectively.Then, the interconnection structure of the power grid can be represented by a connected undirected graph E is the set of unordered pairs (i, j ) representing transmission lines between nodes i and j , which are assumed to be inductive.Notice that The dynamics of generators and loads in a power network are characterized by nonlinear structure-preserving models [7,18] as follows.
(1) Generator bus i ∈ V g : (2) Load bus i ∈ V l : ( All the symbols in the models are described as follows. θ i Voltage phase angle ω i Voltage frequency w.r.t the nominal reference ω * (typically 50 or 60 Hz) -M i > 0 Angular momentum of generator i The inductance of the transmission line connecting nodes i and j -P m i Controllable power generation -P l i Unknown power demand The above structure-preserving power network models are commonly used in stability analysis and controller design of power networks, including the analysis of network synchronization and frequency regulation (see [7,18,47]).For example, [22] studies the full-order description of synchronous generators in a grid setting.However, controlling such model turns out to be rather complicated.Furthermore, it should be remarked that the models in (1) and ( 2) are simplifications for the real power systems based on the assumption that the influence of the windings of the generators is negligible and the power lines are assumed to be lossless.Therefore, only very slow electromechanical transients will be affected by the type of feedback that is proposed for use in this paper.For example, transients as described in [30] would not be applicable.
To obtain a linear model based on (1) and ( 2), we follow [36,44] and assume that the differences of the phase angles are relatively small, i.e., θ i − θ j ≈ 0 for any i, j ∈ V g ∪ V l , which is satisfied in a vicinity of the nominal condition.This assumption then leads to a linearization of the power network model in the following compact second-order form p : M g 0 0 0 where θ g ∈ R n , θ l ∈ R m are the state vectors collecting θ i with i ∈ V g and i ∈ V l , respectively.Matrices M g , D g , and D l are diagonal and positive definite, which are defined as Moreover, P m and P l are collections of P m i with i ∈ V g and P l i with i ∈ V l , respectively.The weighted Laplacian matrix L ∈ R (n+m)×(n+m) of the network topology G is partitioned as The (i, j )-th entry of L is given by which is interpreted as the maximum real power transfer between any two nodes i and j with constant voltage levels.The value of W ij is positive if there is a physical cable directly connecting nodes i and j and 0 otherwise.L has a very special structure as stated in the following lemma [12,40].
Lemma 1 If L is a Laplacian matrix associated to a connected weighted undirected graph G, then 1 .The diagonal entries of L are strictly positive while the off-diagonal components are negative or 0. 2 L T = L ≥ 0 and L has only one 0 eigenvalue. 3 1 T L = 0 and L1 = 0.
It is a crucial task to maintain the frequencies of generators and loads in the power network close to the nominal value.However, with the growth of the network size, centralized controllers become increasingly expensive due to the need for the information of all generators.Recently, distributed averaging PI controllers have been proposed for the frequency control of power networks (see [1,21,36] for the details).The distributed controllers exchange information over a communication network which is assumed to be undirected and connected.It should be emphasized that the topology of the communication network V g is generally different from the topology of the power network G.For each generator i ∈ V g , the controller takes the form Here, N i represents the set of neighboring generators of node i in G c , i.e., the collection of generators that generator i communicates with.The coefficients c ij = c ji > 0 reflect the strengths of the connection between generators i and j .Q i is a positive gain.The controller ( 6) is designed to regulate the frequency deviation to 0 and enforce the controller states ξ i to reach a consensus, i.e., the generated power deviation of all generators becomes equal at steady state.Notice that the controller (6) only need the information from the neighbors of generator i; therefore, it can be implemented in a distributed fashion.We now write the controller in the vector form c : where L c is the Laplacian matrix of the communication network whose (i, j ) entry is given by c ij .The vector ξ ∈ R n is the collection of ξ i , i ∈ V g , and . ω g := θg ∈ R n represent the frequencies of all the generators.By substituting the distributed controller (7) for P m in (3), we obtain the closedloop system (9) Here, d := P l is the power consumption of the m loads, which is uncontrollable and regarded as the stochastic disturbance of the closed-loop system (8).We analyze the stability of cl in the following theorem.

Theorem 1
The closed-loop power system cl is semistable, and the trajectories of its impulse response converge to 0.
Therefore, E is nonsingular, and its inverse reads as which gives By [36], interconnecting the controller c and the power network system p results in a zero frequency deviation and consensus of the controller states.It then implies that, with impulse signals as disturbances of the closed-loop system cl , all the state trajectories converge to constant values, i.e., lim t→∞ e E −1 At exists.Thus, cl is semistable by the definition in [8].Furthermore,using the third property of Laplacian matrices in Lemma 1, it is easy to check that Together with semistability, we conclude that E −1 A only has one 0 eigenvalue at the origin, and all the other eigenvalues have strict negative real parts.Notice that the vectors are the right and left eigenvectors of E −1 A corresponding to the 0 eigenvalue, respectively, i.e., E −1 Av R = 0 and v T L E −1 A = 0. Based on this, we can obtain the following decomposition where U is unitary, and Then, the convergence value of the impulse response is computed as That completes the proof.
Remark 1 Due to the singularity of matrix E −1 A, the closed-loop system cl is not asymptotically stable.However, Theorem 1 indicates that, for any initial condition, the unforced system responses can reach the steady states at 0. Theorem 1 also offers a physical interpretation for the power network.The distributed controller in (7) eliminates the effects of disturbances, namely a sudden impulse change in the power demand, and steers the frequencies of all the nodes (i.e., the generator and load buses) to the nominal value ω * .

Clustering of power networks
To uncover the community structure of generators and loads, a novel clustering method is proposed and applied in this section, where the clusters are constructed based on the differences of node behaviors.

Concepts of clustering
As the physical components of a power network can be interpreted by an undirected weighted graph [38], the techniques of graph analysis help us achieve a better understanding of the essential properties of the power system.Especially, the graph clustering can be applied to identify the communities of generators and loads, which group the units with similar behaviors.In this subsection, we first recapitulate the notion of network clustering (or graph partition) from, e.g., [12,13].
Then, network clustering is to partition V into r (r ≤ k) disjoint clusters which cover all the elements in V.
Mathematically, the clustering of a graph G can be represented by a binary characteristic matrix, which is defined as follows.
The characteristic vector of the cluster C i is defined by binary vector π(C i ) ∈ R k where 1 T k π(C i ) = |C i |, and the j th element of π(C i ) is 1 when j ∈ C i and 0 otherwise.Then, the characteristic matrix of the clustering is a binary matrix denoted by In the following subsections, we propose a new method for cluster selection in power networks.The choices of clustering will determine the accuracy of network approximation.i.e., how close the behaviors of the reduced and original networks are.Hence, the cluster selection algorithm is the most crucial part of the clustering-based model reduction.This paper provides a new way of cluster selection, which involves a particular notion of dissimilarity and an adaption of the hierarchical clustering algorithm.The details of our method are described hereafter.

Characterization and computation of dissimilarity
In the context of clustering, node dissimilarities intuitively describe how different the nodes are from each other.In contrast with the existing results in the literature, such as [16,29,39], the clusters of a power network in this paper are identified without sampling data from a real system.We calculate the dissimilarities of network nodes by the H 2 -norms of transfer function discrepancies, rather than the Euclidean distances of the relative positions.Then, potentially, we can place the nodes with similar behaviors into the same clusters.
More precisely, by considering the state variables of the closed-loop power system cl in (8), the behavior of generator i (i ∈ V g ) is characterized by its responses of the voltage angles θ i , the frequencies ω i , and the controller states ξ i .We include ξ i into behavior characterization of generators since we also approximate the controller c .When the generators are clustered, the communication network of the distributed controllers is automatically simplified.Hence, the effects of the controller states ξ i are also needed to be included.As for the behavior of load i (i ∈ V l ), it can be simply represented by its voltage angle responses θ i with respect to the stochastic disturbances d.
The behaviors of generators and loads can be expressed in the complex frequency domain by the transfer functions from d to the indicative state variables.Denote two new binary vectors where e i is the ith unit base vector.Then, the behavior of generator i is given by where θ i (s), ω i (s), ξ i (s), and d(s) are the Laplace transforms of the states θ i , ω i , and ξ i , and the input d, respectively, in the closed-loop system.Furthermore, the behavior of load i ∈ V l is represented by Thereupon, a pairwise dissimilarity of nodes i and j is defined for generators and loads, respectively.
Moreover, the dissimilarity matrices for the generators and loads are denoted by D g ∈ R n×n and D l ∈ R m×m , which are constructed by collecting D g ij and D l ij .Clearly, D g and D l are symmetric matrices with nonnegative entries and zero diagonal elements.Besides, the boundedness of D g and D l is also guaranteed by Theorem 1.The reasons of the boundedness are explained as follows.
For simplicity, let where E, A, and B are system coefficients of the closed-loop system cl in (8).Note that impulse response of node i in cl is expressed by e At B, which is a bounded smooth function of t and exponentially converges to 0. Thus, the boundedness of D g ij and D l ij can be seen from the the definition of H 2 -norm [2].
and similarly, In both equations, P is the controllability Gramian of the closed-loop system cl .For large-scale networks, ( 23) and ( 24) provide an efficient method to compute matrix D, since we first calculate P and then just apply vector-matrix multiplication to obtain all the entries of D. However, the semistability of cl poses a challenge to implement this idea.For asymptotically stable systems, their controllability Gramians can be uniquely determined by solving the associated continuous-time algebraic Lyapunov equation In the semistable case, A is not Hurwitz, which results in multiple solutions of the above Lyapunov equation [5].To characterize the controllability Gramian from the solutions of ( 25), the following lemma is proven.

Lemma 2
The controllability Gramian P of the semistable power system cl is positive semidefinite, and it is uniquely determined by the combination of ( 25) and where v L is the left eigenvector of A given in (15).
Proof Since v T L A = 0, for any t, we have Thus, (26) holds.It then yields v T L Pv L = 0, which means P is positive semidefinite.Next, we show that solving the algebraic Lyapunov equation (25) together with (26) will yield a unique solution.A contradiction proof is provided.Assume that two distinct symmetric matrices P 1 and P 2 are both the solutions of ( 25) and (26).From (25), we have which leads to Note that lim t→∞ e At = v R v T L .Then, we obtain from (30) that which is equal to 0 due to (26).Thus, it contradicts the fact that P 1 = P 2 .As a result, the common solution of ( 25) and ( 26) is unique.
In the following theorem, we then provide a method to determine the controllability Gramian of the semistable system cl by solving linear matrix equations.
Theorem 2 Let P a be an arbitrary solution that fulfills (25) and ( 26) simultaneously.Then, the controllability Gramian P of the closed-loop power system cl in ( 8) is computed as where J is a constant matrix defined by 3n+m) .(33) Proof Equation ( 27) implies that J PJ T = 0. From (31) in the proof of Lemma 2, we have which leads to (32).
It should be remarked that P in (32) does not depend on the choice of P a , which can be any symmetric real solution of (25).Besides, we can modify the Hessenberg-Schur method (SB04MD (SLICOT)) in [25] to solve the Lyapunov equation with a nonsingular A matrix.The output of the adapted Hessenberg-Schur algorithm will give one solution of ( 25), and we then use the relation in (32) to obtain the real value of P. Hereafter, we apply ( 23) and (24) to determine all the dissimilarities of different pairs of generators or loads in the power grid.

Hierarchical clustering of generator and load buses
The idea of hierarchical clustering has been extensively explored in many fields, including pattern recognition, data compression, and network science (see [4,29]).Hierarchical clustering, in principle, is a greedy algorithm, whose running time grows polynomially with the size of the studied networks.Compared to the K-means algorithm, another popular clustering method, the result of hierarchical clustering is not affected by the initializations of the graph partitions.Due to the above qualities, we choose and adapt the hierarchical algorithm to cluster the generator and load buses for the purpose of model order reduction.
In this paper, the generators and loads are grouped separately on account of different evaluation criteria of dissimilarities (see Section 3.2), yet the same hierarchical clustering procedures can be applied for both.It is because the hierarchical clustering only needs the dissimilarity matrices D g and D l , whose entries already contain the dissimilarities of all pairs of generator and load nodes.Hereby, we only analyze the clustering of generators and the clustering of loads in parallel.Both procedures go in the same way; hence, we only present it for any D ∈ {D g , D l }.
As nodes are merged into small clusters, we must measure how similar two clusters are.In this paper, we consider the average cluster dissimilarity.It characterizes the dissimilarity of two clusters as the average of D ij over all node pairs i and j that belong to distinct clusters [4].More precisely, the dissimilarity of clusters C μ and C ν is defined as follows.
The value of (C μ , C ν ) is computed only using the entries of D that correspond to the nodes in C μ and C ν .Hereupon, the clustering approach links the pairs of nodes that are in close proximity and place them into binary clusters.Then, the newly formed clusters are merged into larger clusters according to the cluster dissimilarity.In Algorithm 1, the pseudocode of hierarchical clustering algorithm is presented.The inputs are the dissimilarity matrix D and the desired number of clusters r.The characteristic matrix of the clustering is considered as the output of the algorithm.

Algorithm 1 Hierarchical clustering
1: Compute the dissimilarity matrix D by ( 23) or (24) 2: Place each node into its own singleton cluster, that is Set m to be an arbitrary large number 6: Compute (C i , C j ) by (35) 8: end for 12: Merge cluster μ and ν into a single cluster 13: In principle, hierarchical clustering does not require preliminary knowledge about the number and the size of clusters.In practice, it generates a dendrogram that organizes a hierarchy of clusters in a tree structure.Any cut of the hierarchical tree offers a potentially valid clustering of the studied network.
It is also worth mentioning that he formation of clusters in Algorithm 1 does not intend to only cluster the adjacent nodes.Even if two nodes are not neighbors, their dissimilarity still can be computed and thus their difference of behaviors is measured, and two nonadjacent vertices can be clustered when they have relatively similar behaviors, as the obtained aggregated network will achieve a small approximation (see some numerical examples in e.g.[12,13]).However, our method does not exclude the possibility of modifying the algorithm to only aggregate adjacent nodes.To this end, we need to update the computation of dissimilarity matrix, i.e., we add a penalty factor if two nodes are not neighbors.

Reduced model of power network
In Section 3.2, the dissimilarities are calculated in the closed-loop systems, and the generator and load buses are grouped separately by the hierarchical clustering in Section 3.3.Then, the controllers are automatically clustered as each of them is associating with a generator.
Suppose that we acquire r clusters of generator buses and q clusters of load buses.Correspondingly, the characteristic matrices are denoted by g ∈ R n×r and l ∈ R m×q , respectively.Let Then, a simplified weighted power network Ĝ can be obtained by aggregating all the nodes with the same clusters in the original complex network.Mathematically, the reduced Laplacian matrix of Ĝ is generated by the following projection.
where L is the weighted Laplacian matrix representing the original network.Similarly, the inertia and damping coefficients in the reduced network are given by Mg = T g M g g , Dg = T g D g g ∈ R r×r , and Dl = T l D l l ∈ R q×q .(38) This clustering-based projection also has a physical interpretation, as shown in Fig. 1, which shows the aggregation of generator and load buses, respectively.Notice that the nodes in the same cluster can be aggregated even if they do not have a direct connection with each other.
The transmission lines contained in a cluster are neglected, while those connecting two clusters are aggregated.Moreover, the maximum real power transfer on the new branch between the clusters C μ and C ν is given by where W ij is the maximum real power transfer between nodes i and j .In the reduced network, a node represents a cluster of nodes in the original network, and its angular momentum and damping coefficient are the sum of these parameters of generators and loads in the original power network.Specifically, the inertia and damping coefficients of the cluster C i is computed as  (37) and (38), respectively.θg ∈ R r and θl ∈ R q are the voltage phase angles of the aggregated generators and loads in the reduced-order power network, which are used to approximate the states in the original system: θ g ≈ g θg , ω g ≈ g ωg = g θg , and θ l ≈ l θl . ( Based on (41), the reduced-order nonlinear power system can be also constructed, which follows the same form as (1) or (2).Here, the reduced-order nonlinear model is omitted.
As the generators are clustered, the size of the communication network G c for the distributed controller in ( 7) is also simplified.As a result, the dimension of the controller is reduced simultaneously.Analogously, the representation of the lowerdimensional controller is written as where Q := T g Q g and ξ ∈ R r are states of the reduced-order controller.Lc := T g L c g represents the simplified communication network.
Combining the reduced versions of the power network ˆ p and the distributed controller ˆ c , we derive the reduced-order closed-loop system as follows.
where xT := θT (45) Here, d := P l is the uncontrollable power demand in (8).The Galerkin projection matrix for the closed-loop system cl is which satisfies Ê = T E , Â = T A , and B = T B. (47) Next, the properties of the reduced-order models are discussed in the following theorem.

Theorem 3
The reduced power network p in (41) and the distributed controller c in (43) preserve the network structures.The closed-loop system cl in (44) is also semistable, and the trajectories of its impulse response converge to 0.
Proof In [13], we have proven that the Galerkin projection based on the characteristic matrix of network clustering which can preserve the algebraic structure of a Laplacian matrix, which means L in (37) and Lc in (43) are the reduced Laplacian matrices representing the simplified power network and communication links.Besides, (40) shows that Mg , Dg , and Dl are diagonal positive definite.Therefore, the reduced models in (41) and (43) have the same structures as ( 8) and (7), respectively.
Following the same reasoning line of Theorem 1, the closed-loop system cl is semistable.It has only one pole at the origin, and all the other poles are located in the open-left half plane.The vectors are the left and right eigenvectors of Ê−1 Â corresponding to the only zero eigenvalue, which can be verified by vT L Ê−1 Â = 0, and Ê−1 Â vR = 0. Note that the characteristic matrix of the graph clustering has the property g 1 r = 1 n .Hence, vR 2 = 1, and vT L vR = 1.We further obtain (49) That completes the proof.

Case study
We illustrate the proposed method on the IEEE 30-bus test system [46], which contains 6 generator buses and 41 transmission lines.The graph representation of this power system is depicted in Fig. 2a.In order to fit the network data into the index setting of the states in our model (3), we use the bus numberings that are different from the original system date.Now, the generator buses are with the indices from 1 to 6. Assume that the distributed controllers of the generators are identical and connected based on the communication network in Fig. 2b.The control parameters in (6) are given by Q i = 0.1 for all generators and c ij = 1.Thereby, the closedloop power system in (44) is established, which include 42 states, namely the voltage phase angles of generators and loads, θ g ∈ R 6 and θ l ∈ R 24 , the generator voltage frequencies ω g ∈ R 6 , and the controller states ξ ∈ R 6 .The positive semidefinite controllability Gramian is then obtained by Theorem 2. Consequently, the dissimilarity matrices of generator and load behaviors are computed using ( 23) and (24).Due to space reasons, we only give the results for generators as follows.
We then group the generator and load nodes by Algorithm 1, which divides the node sets in the original network into several subsets that contain nodes with small dissimilarities to each other.The clustering results of the generator and load nodes are straightforwardly interpreted by the dendrograms in Fig. 3.The leaves are the bottom vertical lines representing the buses, and the clustering of two nodes is indicated by merging two leaves into a single branch.The dissimilarity can be read from the horizontal position of each fusion.In this example, we cluster node sets of the generators and loads as in Table 1.The clustering results are reasonable owing to the hierarchical structures of the dendrograms.
Using the characteristic matrices of the resulting clusterings for Galerkin projection, we obtain the reduced-order models of the power networks and distributed controller in forms of ( 41) and (43).The simplified networks are depicted in Fig. 4.  Therefore, the network structures are preserved through the reduction process, which means that we can use the power network with much smaller size to approximate the behavior of the original one.
Next, the quality of the approximation is evaluated.We consider the outputs of the closed-loop system cl as θ g , θ l , ω g , and ξ , respectively.Correspondingly, they are approximated by g θg , l θl , g ωg , and g ξ , which are assumed to be the outputs of reduced-order model and controller.Based on this, we compute the approximation errors of those variables in terms of H 2 -norms (see Table 2).The errors are not significant if we think of the values of the node dissimilarities in (50) and the dimension of the reduced closed-loop model which is considerably lower compared to the original system.
Next, the performance of the reduced-order model is further demonstrated in the time domain.The states of both the original and simplified networks are initialized at zero.From 0 to 20 s, the power demands of the 24 loads, P l , are assumed by a random vector with all the entries in the range [0, 3].After 20 s, the demands are set to be a new random constant vector, which is generated from [0, 1].The unit of P l is 100 MVA.In Fig. 5, the state trajectories of both systems are compared.Evolving over time, each state trajectory of the reduced network shows a similar directional tendency to the trajectories representing a cluster of nodes in the original network.Furthermore, the most of the response curves of the full-order system are approximated by those of the reduced model with acceptable errors.Besides, from  Fig. 5c, we can see how the full-order and reduced-order distributed controllers regulate the generator frequency deviations back to 0 in both systems.Moreover, Fig. 5d illustrates that the states of both controllers are synchronized to the same trajectories.
We implement this numerical test by MatLab 2016a in the environment of a 64bit operating system with Intel Core i5-3470 CPU at 3.20 GHz, RAM 8.00 GB.We observe that the evaluation of dissimilarities takes 0.0564 s, while the clustering algorithm only uses 0.0033 s to find three clusters for both generator and load buses.Therefore, the total computation time is mainly consumed by the calculation of the controllability Gramian for dissimilarity evaluation.
Fig. 5 The comparisons of the state of the full-order and reduced-order power systems, where the solid and dashed lines are representing the trajectories of the original and simplified networks

Conclusion
We have considered the reduced-order modeling of power networks and distributed controllers, which are expressed as semistable second-order and first-order differential algebraic equations, respectively.By exploiting the controllability Gramian of the semistable closed-loop system, we have proposed a novel notion of node ity and apply a hierarchical clustering approach to divide the generator and load buses into several subsets.Towards the preservation of network structures in the reducedorder models, the characteristic matrices of the resulting clusters are adopted for the Galerkin projections of both power system and its controller.The explicit reducedorder models are established in the same forms of the original models, which inherit a network interpretation for the interconnections of the power units.A numerical example, at last, has shown the performance of the reduced-order model and controller.This paper provides an idea to explore clusterings of controlled power networks based on dynamical models.As a further extension, other representations of interconnections, generators, loads, and controllers may be investigated, e.g., including resistors on the transmission lines, or more complicated nonlinearities in the generators, or different types of controllers.

Fig.
Fig. Illustration of aggregation of generator buses (a) and load buses (b).The nodes in the same cluster are included in the dashed boxes

Fig. 2 a
Fig. 2 a The topology of IEEE 30-bus test system.The generators and load buses are represented by circles and squares, respectively.b The communication network that links the distributed controllers on 6 generators

Fig. 3
Fig.3 Dendrograms showing the clusterings of 6 generator buses (a) and 24 load buses.The horizontal axis are labeled by bus numberings, and the dissimilarity data are read from vertical axis

Fig. 4
Fig. 4 The reduced topology of the power transmission lines and communication links, which are represented by solid and dashed edges, respectively

Table 1
Clustering results of generator and load buses

Table 2
The approximation errors evaluated by H 2 -norms