{'fraction_non_alphanumeric': 0.04944122637645544, 'fraction_numerical': 0.01374992686209116, 'mean_word_length': 4.126374725054989, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimal Control of Large-Scale Networks using Clustering Based Projections

In this paper we present a set of projection-based designs for constructing simplified linear quadratic regulator (LQR) controllers for large-scale network systems. When such systems have tens of thousands of states, the design of conventional LQR controllers becomes numerically challenging, and their implementation requires a large number of communication links. Our proposed algorithms bypass these difficulties by clustering the system states using structural properties of its closed-loop transfer matrix. The assignment of clusters is defined through a structured projection matrix P, which leads to a significantly lower-dimensional LQR design. The reduced-order controller is finally projected back to the original coordinates via an inverse projection. The problem is, therefore, posed as a model matching problem of finding the optimal set of clusters or P that minimizes the H2-norm of the error between the transfer matrix of the full-order network with the full-order LQR and that with the projected LQR. We derive a tractable relaxation for this model matching problem, and design a P that solves the relaxation. The design is shown to be implementable by a convenient, hierarchical two-layer control architecture, requiring far less number of communication links than full-order LQR.


I. INTRODUCTION
A vast majority of practical networked dynamic systems (NDS), ranging from power system networks to wireless networks to social or biological networks, consist of several hundreds to thousands of nodes that are spatially distributed over wide geographical spans. Developing tractable control designs for such large complex networks, and implementing those designs through affordable communication, continue to be a challenge for network designers. Traditionally, control theorists have addressed the problem of controlling largedimensional systems by imposing structure on controllers. The most promising approach, for example, started with the idea of decentralized control [1], followed by techniques such as singular perturbation theory [2], [3], balanced truncation [4], [5], [6], and ν-gap reduction [7] among others. These methods aim to simplify the design of controllers for large systems by exploiting weak coupling between their state variables, and by ignoring states that are 'less important' than others. The tradeoff, however, is that the resulting controllers are often agnostic of the natural coupling between the states, especially the coupling between the closed-loop states, since many of these couplings are forcibly eliminated to facilitate the design itself. N Therefore, extending these methods to facilitate controller designs for networks, especially to NDS whose states may be defined over highly structured topologies such as clustering, is quite difficult. The literature for developing tangible and yet simple low-dimensional controllers that satisfy global stability and dynamic performance requirements of very large NDS is still unfortunately sparse. Ideas on aggregate control [8], glocal control [9] and hierarchical control [10], [11] have recently been proposed to address this challenge. The goal of these designs, however, is to guarantee stability by modular tuning of local controller gains; their degrees of freedom for guaranteeing performance can be limited.
To bridge this gap, in this paper we propose a design method called control inversion. The approach is to cluster the states of an n-dimensional network into r > 0 distinct, non-overlapping groups. We assume n to be a large positive integer, and r ≤ n to be a given design parameter. The grouping is defined by a (r × n) structured projection matrix P whose elements denote the identity of states in the clusters, weighted by certain projection weights. The design thereafter consists of three steps. First, for the full-order network an ndimensional LQR controller is defined for any given choice of Q and R. We refer to this controller as the benchmark LQR. Second, the projection matrix P is used to construct an rdimensional reference model for which an r-dimensional LQR controller is designed. The design matrices for this reducedorder controller, however, are not free; they are constrained by being related to Q and R through P . The important point, however, is that the design dimension reduces to r from n. Finally, this reduced-order controller is projected back to the full-order network by the inverse projection P T . The problem is then to find a projection matrix P that minimizes the H 2norm of the error between the transfer function matrices of the full-order network with the benchmark LQR controller and with the projected LQR controller.
This problem by itself, however, is non-convex even without any structural constraint on P . To bypass this intractability, we introduce a relaxation, which is done in two stages. The first stage relaxes the error minimization to the minimization of its upper bound, while the second stage applies a low-rank approximation. We finally design a P that solves this relaxed problem. Three distinct variants of the design are proposed. In the first case, we optimize over cluster assignment while keeping the projection weights fixed, and establish that this minimization can be posed as an unsupervised clustering problem. We use weighted k-means [12] to solve this minimization. In the second case, we fix the cluster identities, and optimize over the projection weights. Depending on the stability of the open-loop system, we show that this minimization can be posed as finding the dominant eigenvector of the controllability Gramian or as finding the Z-eigenvector of a tensor [13]. In the third case, we propose an iterative method to optimize over both cluster assignment and projection weights. The controllers resulting from all three algorithms are shown to be implementable by a convenient, hierarchical two-layer control architecture, requiring far less number of communication links than full-order LQR as well as sparsity-promoting LQR [14].
Recently, [15] and [16] have used structural projectionbased ideas for model reduction of large networks, but not for control designs. Attention has also been drawn to designing LQR controllers for large systems by finding low-rank solutions of algebraic Riccati equations [17]. However, like most Krylov subspace-based reduction methods such as in [6], the controller in [17] is unstructured, and hence demands as many communication links as the full-order LQR itself. Distributed controllers using model matching [18], sparsitypromoting LQR in [14] and structured LQR in [19], [20] promise to reduce the communication density, but their designs inherit the same dimensionality as the full-order design. Unlike all of these methods, the novelty of our algorithms is in the facilitation of closed-loop control from the perspective of both design and implementation. The recent papers [21], [22] also address both goals, but the dimensionality of their controllers is subject to the sparsity structure of the open-loop network while our design does not necessarily require any such sparsity. Preliminary results on this design have been presented in the recent conference paper [23], but only for a consensus model with specific Q and R matrices.
The remainder of the paper is organized as follows. Section II formulates the problem of clustering-based optimal control. The relaxation for the original problem is derived in Section III, which is then solved by a clustering algorithm based on weighted k-means optimization in Section IV. The design for cluster weights as well as that for the weights and clusters taken together are discussed in Section V. All three algorithms are illustrated via simulations in Section VI. Section VII concludes the paper.
Notation We will use the following notations throughout this paper: |m| absolute value of a scalar m |S| c cardinality of a set S 1 n column vector of size n with all 1 entries with a realization form of g(s) = A B C D . We refer to g(s) as stable if A is Hurwitz, and unstable otherwise. Furthermore, the H 2 and H ∞ norms of a stable transfer matrix g(s) are defined by g(s) H2 = From graph theory, a graph G = (V, E) is defined over a node (vertex) set V = {1, ..., n} and an edge set E ⊂ V × V, which contains two-element subsets of V. If {i, j} ∈ E, we call nodes i and j adjacent, and denote the relation by i ∼ j, or simply ij. The set of nodes adjacent to i ∈ V is noted by N i = {j ∈ V|i ∼ j}. In this paper, G is assumed to be undirected, which implies ij is equivalent to ji, and there are no loops or multiple edges between nodes.

II. PROBLEM FORMULATION
Consider a general LTI system of the form where x(t) ∈ R n , u(t) ∈ R m and y(t) ∈ R p represent the vector of state, control and output variables respectively, and d(t) ∈ R n b is a disturbance entering into the system. We assume (1) to be defined over a network of n s ≤ n interconnected subsystems, with their network topology represented by a connected graph G = (V, E), V = {1, ..., n s }. The dynamics of each subsystem can be written as .., B dns ) and C = diag(C 1 , ..., C ns ). Notice that the dimension of u i (t) can be zero, meaning that the i th subsystem can have no input. In this paper, we consider an LQR design for (1), and assume C = I n for full-state feedback. Given two real-valued matrices Q = Q T 0 and R = R T ≻ 0, the LQR problem is posed as finding a feedback law u(t) = −Kx(t) such that the cost function is minimized. The expression (3), also known as the infinitehorizon continuous-time LQR, can be solved by the following algebraic Riccati equation (ARE) where G = BR −1 B T . The feedback matrix can be found through K = R −1 B T X. For such a solution X to exist, we will adhere to the following assumption throughout this paper.
According to [7], the assumption above guarantees a unique stabilizing solution X = X T ≻ 0. However, finding this solution from (4) in practice is subject to O(n 3 ) computational complexity, which can become unscalable for largescale systems. Moreover, the resulting matrix X is usually an unstructured dense matrix, which demands every subsystem in the network to communicate with every other subsystem for implementing the feedback. These two factors together make both the design and implementation of u = Kx very difficult, especially when G consists of thousands to tens of thousands of nodes. Therefore, we propose a design strategy, which we refer to as control inversion, to repose this LQR problem using a clustering-based projection.

A. Control Inversion
Definition 2.2: Given an integer r, where 0 < r ≤ n, and a non-zero vector w ∈ R n , define r non-empty, distinct, and non-overlapping subsets of the state index set V s = {1, ..., n}, respectively denoted as I = {I 1 , ..., I r }, such that I 1 ∪ ... ∪ I r = V s . A clustering-based projection matrix P ∈ R r×n is defined as where w Ii = [w Ii{1} , ..., w Ii{|Ii|c} ] T is non-zero, and I i {j} denotes the j th element in the set I i . The matrix P has the following three properties: • It is row orthonormal, i.e. P P T = I r ; • Image of P T P lies in the span of w, i.e. P T P w = w; .., r. The construction of P is shown by the following example.
Given P defined over any clustering set I and weight vector w, the control inversion strategy for the LQR problem (3) is then composed of the following three steps.
1) Projection to reduced-order system: Using projection P , we first construct a reduced-order model whereÃ := P AP T ∈ R r×r ,B := P B ∈ R r×m ,B d := P B d ∈ R r×n b andC := P CP T = I r . For this systemx ∈ R r is the state andũ ∈ R m is the control input.
2) Reduced-order LQR design: We similarly project the LQR parameters byQ = P QP T ∈ R r×r , and letR = R such thatG :=BR −1B = P GP T ∈ R r×r . An LQR problem for the reduced-order model (6) is then posed as to minimizẽ with respect toũ(t) = −Kx(t). Here the feedback matrix K =R −1BTX corresponds to the solutionX ∈ R r×r of the reduced-order ARE of (7), which is written as 3) Inverse projection to original coordinates: The solutioñ X from (8) is projected back to the original coordinates through the inverse projection This projected controller can then be implemented in the fullorder model (1) using u = −R −1 B TX x, which implies that the effective feedback gain matrix iŝ

B. Problem Statement
The controllerK is dependent on the projection P through equations (6), (7), and (9). The choice of P is guided in the following way. Consider which is the closed-loop transfer matrix from d to x for (1) with full-order LQR. Similarly, consider which is the closed-loop transfer matrix from d to x for (1) with the projected controller (10). Using (11) and (12), we next state our main problem of interest. Main problem: Given system (1) and an integer r > 0, the problem addressed in this paper is to find a clustering set I and a non-zero vector w such that the corresponding projection matrix P solves the model matching problem 1 However, finding an exact solution for this optimization problem is intractable given that the objective function is an implicit and non-convex function of P , and also because P is defined over a combinatorial structure. Our main contribution, therefore, is finding a tractable relaxation for (13) as a quadratic function of P , and thereafter designing P to solve the relaxed problem. We make the following assumption so that g(s) andĝ(s) both have minimal realization.
Assumption 2.3: The pair (A, B d ) is controllable. Solution Strategy: The outline of our solution strategy is as follows. In Section III, we derive an upper bound relaxation for (13) such that its objective function is quadratic in P . Ideally speaking, one can solve for P from this relaxed problem. The computational complexity for constructing the objective function is, however, O(n 3 ) since it requires the computation of the controllability Gramian of g(s). To bypass this difficulty, a second round of relaxation is applied by exploiting the lowrank (denoted as κ) structure of the controllability Gramian. After these two relaxations, the final objective function, still quadratic in P , can be constructed in O(nκ 2 ) complexity, which is near linear if κ ≪ n. The solution to this optimization is then addressed in two ways -first by finding the clustering set I with a fixed weight vector w (Section IV), and second, by finding w while keeping I fixed (Section V). We also propose to combine these two approaches by an iterative algorithm. The overall design flow and the numerical complexities for each step are previewed in Fig. 1. Detailed explanations of these complexities will be provided in the respective sections to follow.

C. Benefits of Control Inversion
An important point to note is that the physical meaning of the statex(t) of the reduced-order model (6) has no relation to that of the state x(t) of our full-order model (1). This is a key difference of the control inversion design from traditional model-reduction based designs where the reduced-order state vector is typically a direct projection of the full-order state vector. The projection in our design is rather applied on the controllerX instead of x(t). Two natural benefits of this approach are as follow: 1. Tractability of design: The computational complexity for constructing the reduced-order ARE in (8) is O(rn 2 ), while that for solving this ARE is O(r 3 ). The computational complexity required to design P through Sections III, IV and V will be shown to be simpler than the O(n 3 ) complexity of a full-order LQR design. Thus, if r ≪ n, the overall control inversion design becomes numerically more tractable than fullorder LQR.
2. Simplicity in implementation: The projected matrixX = P TX P is a structured r-ranked matrix, which results in a sequential two-layer hierarchical control architecture. The implementation of the feedback u(t) = −R −1 B TX x(t) follows three steps. First, a coordinator is assigned to each cluster I i , which collects the measurements of all the states belonging to that cluster. Each coordinator then computes the weighted averaged state P i,: x for its cluster, i = 1, 2, .., r. Next, the coordinators exchange these weighted averages, and each of them compute the r-dimensional vectorXP x. Note that in this process no coordinator will be able to infer individual state measurements from other clusters. Finally, the coordinator of I i computes the control vector u i by taking linear combinations of the elements ofXP x. The linear combination follows from u = −R −1 B T P TX P x. The individual elements of u i are broadcast to the respective input actuators in I i . Note that since the aggregation is applied to the state x(t), and not on the subsystems, the state vector of any subsystem can be partitioned among different clusters. Thus, in practice, a subsystem may need to transmit its states to more than one coordinator, and also receive control inputs from more than one coordinator. We illustrate the three implementation steps by an example in Fig. 2.
In the worst-case scenario when every subsystem has a scalar state x i and a scalar control input u i , the two-layer control implementation will result in a much sparser communication topology with n + r 2 bidirectional links compared to an LQR controller which would require n 2 number of links. This reduction, combined with standard networking protocols such as multi-casting [24], makes the implementation of our proposed controller convenient and cheap.

III. RELAXATIONS FOR MODEL MATCHING
In this section we describe the theoretical derivation of the relaxation for the model matching problem (13). We start by discussing the well-posedness of the projected controllerK in Section III.A. The final optimization to be solved is then obtained from two stages of relaxations as detailed in Section III.B and III.C respectively. All proofs are presented in the Appendix B.

A. Well-Posedness Conditions
The well-posedness ofK is equivalent to two conditionsnamely, if the reduced-order ARE (8) admits a solutionX, and ifĝ(s) is stable. We discuss these two factors as follows.
1) Existence condition: Similar to the full-order ARE (4), the reduced-order ARE (8) is guaranteed with a unique so- [7]. However, unlike conventional model reduction techniques that can utilize unstructured projections to preserve the exact stabilizability and observability properties from (A, G 1 2 ) and (Q T 2 , A), the structured projection P in (5) does not guarantee that. As a result, for a general system the ARE (8) may not admit a unique solutionX, and thereforeK may not exist. To bypass this problem, we modify the definitions ofQ andG by using a constant shift. That is, in case the pair (P AP T , P G 1 2 ) is not stabilizable and/or (Q T 2 P T , P AP T ) is not observable, we let Q = P QP T + αI r and/orG = P GP T + αI r for a small constant shift α > 0. Using (14),Q andG become positive-definite matrices, which yield (Ã,G 1 2 ) controllable and/or (Q T 2 ,Ã) observable. The existence ofX andK can thus be guaranteed. The matching error introduced by this shift will be discussed in the next subsection.

2) Bound of ARE solution:
We first recall a lemma from [27] about the eigenvalue bound of a general ARE solution.
Lemma 3.1: (Lemma 1.5 in [27]) The solution X from the ARE (4) satisfies the following upper bound: where K t is any matrix stabilizing A + G 1 2 K t , and D t and F are positive-definite matrices that satisfy From this lemma, given any stabilizing controller K t one can always find an upper bound forλ(X) using the function β(A, G, Q). This function will be used next to verify stability ofĝ(s).
3) Stability condition: Before stating the main stability criterion forĝ(s), we first make the following assumption on the projection weight w from Definition 2.2.
Note that from where E = X −X denotes the error between the full-order ARE solution from (4) and the projected solution from (9). Theorem 3.3 provides two options for achieving stability ofĝ(s). In order to meet the condition (17), one can select Q and R such that σ(Q) ≫σ(G) > 0. In practice, this choice of (Q, R) will make the design more robust towards the uncertainties in A and B [29]. As a drawback, a larger Q will also result in a high feedback gain, making the system vulnerable to noise. An alternative to satisfy (17) is to minimizeσ(E). As will be shown in the next subsection,σ(E) is proportional to the objective function of our proposed upper bound relaxation for (13). Therefore, solving the relaxation problem will also assist in enlarging the inequality gap in (17).
The following lemma provides a sufficient condition for (17) by which stability ofĝ(s) can be verified without knowing X.
Lemma 3.4: The stability condition (17) holds if where function β(A, G, Q) is defined in Lemma 3.1. Remark 1: For the case when A is unstable, computing β(A, G, Q) would require the knowledge of K t according to Lemma 3.1. Similarly for w to meet Assumption 3.2, one also needs to know the eigenvectors of all unstable modes of A. Verifying stability ofĝ(s) for an unstable A is, therefore, admittedly more computationally expensive than for a stable A. This computational burden obviously does not exist when A is stable, and may not also exist when A, despite having zero eigenvalues, has specific structural properties. One such example is when A is a weighted Laplacian matrix (i.e., when (1) is a consensus network). We will illustrate this special case of consensus network in Appendix A.

B. Relaxation I: Upper Bound Minimization
We present our first stage of relaxation assuming that the sufficient condition for closed-loop stability from Theorem 3.3 holds. As mentioned before, finding a P that exactly minimizes g(s)−ĝ(s) H2 is an intractable problem. To relax this problem, we find an upper bound for g(s) −ĝ(s) H2 . We denote the error system by g e (s), which can be written as From a similarity transformation of T = I n I n I n 0 and T −1 = 0 I n I n −I n , (19) yields where E = X −X. By taking norms on both sides of (20), we get which from the bounded real lemma [7] and the definition of H 2 norm reduces to where γ is any positive real number such that a real-valued matrix Γ = Γ T 0 exists and satisfies Inequality (21) shows that g e (s) H2 is linearly bounded by EΦ 1 2 F . Therefore, one way to solve the original model matching problem (13) will be to find a P that minimizes EΦ 1 2 F . By doing so,σ(E) can also be minimized to some extent sinceσ(E) ≤ E F ≤σ(Φ − 1 2 ) EΦ 1 2 F , which will help in meeting the stability condition (17). This type of bound minimization is common in model and controller reduction, and has been attempted (see [17] and the references therein) under the assumption that P is unstructured, or more specifically P is an r-dimensional Krylov subspace By this assumption, E can be found as an explicit function of P associated with a Householder transformation. In our case, however, P has a structure as in (5), due to which this explicit functional relationship does not hold anymore. We, therefore, apply perturbation theory of ARE to further relax the bound in (21), and derive a new upper bound on EΦ 1 2 F as an explicit function of P in the following theorem. where , andβ = sup P β(Ã,G,Q) are positive scalars that are independent of P , and α is defined in (14). The constant αǫ 1 ǫ 3 in (23) represents the matching error introduced by the constant shift α from (14). This error can be disregarded if the reduced-order ARE (8) admits a solution. From (21) and (23), it then follows that g e (s) H2 ≤ γf (ξ), due to which we approach the minimization of g e (s) H2 by minimizing f (ξ) with respect to P . Since f (ξ) is a monotonic function of ξ, the minimization of f (ξ) is equivalent to minimizing the value of ξ as The optimization (RL-1), therefore, serves as an upper bound relaxation for the original model matching problem (13). Note that, in general, it is impossible to exactly quantify the optimality gap between (13) and (RL-1) since (13) is nonconvex even without posing any combinatorial constraints on P . The optimality gap will be small if the minimum value of ξ is close to zero, in which case the error g e (s) H2 in (13) will be nearly zero as well following (21) and (23). The ideal case ξ = 0 will happen when Φ 2 is invariant to the mapping P T P . This holds for the trivial case where r = n and P = I n . To achieve a sufficiently small minimum for (RL-1) for r < n, we will develop two designs in Sections IV and V.

C. Relaxation II: Low-Rank Approximation
We next discuss the numerical complexity in constructing the optimization problem (RL-1), and how this complexity can be simplified by making appropriate approximations on Φ. In the most general case, Φ required for (RL-1) needs to be computed through the following procedures. First, recall the definition of the Hamiltonian matrix The eigenvalues of H are symmetric about the imaginary axis. Suppose H is diagonalizable and that the columns of the matrix Y Z 2n×n span the stable invariant subspace of H, where consists of all the eigenvalues of H in the left-half plane, i.e. 0 > λ − 1 > ... > λ − n . The stabilizing solution of ARE can thus be found by X = ZY −1 [7]. The first n rows of (25) are expanded as A − GZY −1 = Y Λ − Y −1 , which means Λ − and Y are the eigenvalues and right eigenspace of the closed-loop state matrix A−GX. Then from the Lyapunov equation (22), we can write Φ directly in terms of Y and Λ − as where C ∈ R n×n is a Cauchy matrix with and subsequently obtain Φ 1 2 from the Cholesky decomposition. Therefore, to compute Φ and then Φ 1 2 , one will need to compute the full stable eigenspace Y from H following (25). This computation is as expensive as solving a full-order LQR with O(n 3 ) complexity for both computation and memory [28]. This may defeat the purpose of our design since we want our controller to be numerically much simpler than the full-order LQR. To bypass this difficulty, we next show that Φ can be approximated by a matrix Φ κ that follows from a κ-dimensional (κ < n, not necessarily equal to r) invariant subspace of Y . Ideally κ should be at most r to justify the computational benefit of our design while preserving an acceptable accuracy in the error norm ξ. This matrix Φ κ is constructed as follows.
By definition of Φ κ , one only needs to compute the first κ eigenvalues λ − 1 , ..., λ − κ of H, and the Y 1 component of the first κ eigenvectors. Ω 1 can be approximated by the pseudo-inverse of Y 1 . These κ smallest eigenvalues and eigenvectors can be solved by Krylov subspace-based techniques such as Arnoldi algorithm in O(nκ 2 ) time [28]. Therefore, instead of (RL-1), we consider solving a computationally simpler approximation of (RL-1) as where Φ 1 2 κ ∈ R n×κ from Definition 3.6 can be computed as Φ The optimality gap between optimizations (RL-1) and (RL-2) can be quantified by the following lemma.
Lemma 3.7: Assume Y −1 has a moderate condition number η, and each column of B d has a unitary norm. The minimum ξ * from (RL-1) and ξ * κ from (RL-2) satisfy From Lemma 3.7, the optimality gap between (RL-1) and (RL-2) will be negligible when the error is kept small. In practice, this situation happens when there exists κ (κ ≪ n) dominant eigenvalues in the Hamiltonian matrix, i.e., H has the following spectral gap This gap can exist if the open-loop network (1) exhibits coherent behavior [3]. As a result of this spectral gap, the RHS of (28) can become sufficiently small, in which case the optimal value of (RL-2) will closely resemble that of (RL-1) while the computation of the objective function requiring a much tractable complexity of O(nκ 2 ), κ < n. The idea of utilizing a κ-dimensional subspace for computing Φ κ is similar in spirit to finding an unstructured approximate ARE solution as proposed in [17]. However, it should be noted that unlike [17] where the selection of the κ eigenvectors is undetermined, for our problem the error bound in Lemma 3.7 clearly guides the choice of the κ eigenvectors in terms of tightening the optimality gap between (RL-1) and (RL-2).

IV. DESIGN I: CLUSTER DESIGN
In this section we present an algorithm to design P by solving (RL-2). Note that P has two degrees of freedom -I and w. For the design in this section, we keep w fixed, and minimize (RL-2) over I. Although inherently this is an NPhard problem, fortunately the specific structure of the objective function ξ κ , together with the structure imposed on P in (5), enables (RL-2) to be solved by efficient numerical algorithms such as weighted k-means [12]. We show these results as follows.
To establish the equivalency of (RL-2) to the weighted kmeans optimization, it is useful to borrow a nominal projection matrixP as From Definition 2.2,P satisfies P =P W , where W = diag(w). With this notation, we have where Ψ = [ψ 1 , ..., ψ n ] T denotes the matrix W −1 Φ 1 2 κ . Therefore, the entries of the matrixP TP W 2 can be found by .., r. Thus the j th row of the matrixP TP W 2 Ψ can be written as for j ∈ I i . It is clear from above that for all the index j that are assigned to the same cluster I i , c T i is a weighted average (or a weighted centroid) of the row vectors ψ T j . Moreover, the matrixP TP W 2 Ψ will have identical rows for those whose indices are inside the same cluster. Therefore, (RL-2) can be posed as an unsupervised clustering problem minimize I1,...,Ir The optimization problem in (32) is in the same form as a weighted k-means optimization, which minimizes the Euclidean distance weighted by w 2 j between each data point ψ j and its centroid c i . Thus, data points which are close to each other in the weighted distance are assigned to the same cluster. A standard method for solving this problem is Lloyd's algorithm [12], using which we present the weighted k-means clustering for (32) in Algorithm 1. With Φ 1 2 κ ∈ R n×κ as the input, the running time of Lloyd's algorithm is O(nκrk), where k is the total number of iterations. As in any heuristic algorithm, Algorithm 1 does not guarantee convergence to a global minimum. Hence, if one is not satisfied by the k-means solution one can apply any state-of-art clustering algorithms for solving (32). Using the resulting clustering set I and the fixed weight w, we can then construct the projection P and execute the control inversion design to getK. The clustering weight w can be selected as any vector that satisfies Assumption 3.2. When A is Hurwitz, a simple choice of w can be the vector of all ones.
We conclude this section by summarizing the total numerical complexity for our design based on Algorithm 1. The chain of approximations involved in this design till the output of Algorithm 1 follows the sequence of equations: (21), (23), (RL-1), (RL-2) and (32). The total complexity amounts to O(nκ 2 ) + O(nκrk) + O(n 2 r) + O(r 3 ), which includes construction of the objective function for (RL-2), execution of Algorithm 1, computation of reduced-order matrices triple (Ã,Q,G), and solving the reduced-order LQR (7), respectively. This complexity can be at most O(n 2 r) if κ ≤ r, which is more tractable compared to the O(n 3 ) complexity of full-order LQR, especially when r and κ are sufficiently small. We next state a variant of our proposed controller where we solve (RL-2) by fixing the clustering set I and varying the projection weights in w instead. This type of a controller may be needed when a subsystem prefers to have all of its states assigned to the same cluster. Multiple subsystems in the network may also exhibit spatial clustering based on their geographical proximities, in which case they may prefer to have a dedicated coordinator for themselves. This scenario commonly arises in power system networks. Utility companies always prefer to send the state information of their generators to only their own local control centers instead of sharing that information with any other company. Therefore, in this section we develop a new set of algorithms where we fix I to represent the identities of the desired clusters, and minimize ξ κ in (RL-2) over w. A visual interpretation of this approach and its comparison to the design in Section IV are shown in Fig. 3. Note that the optimal values of w so obtained denote the relative importance of the network nodes in the closedloop system with the projected controller. In Fig. 3c we show this relative importance by shrinking or expanding the size of the nodes. We describe the design for two cases depending on whether (1) is stable or not.

A. Case I: Stable Open-Loop
We consider the same optimization objective as in (RL-2) but now minimize it with respect to w as To solve (33), we consider a binary projection matrixP aŝ As can be verified from Definition 2.2,P here also satisfies P =PŴ withŴ = diag(ŵ), where thisŵ is defined bŷ such thatŵ T Iiŵ Ii = 1, for i = 1, ..., r. Using these notations, we can rewrite the objective function in (33) as Since tr(Φ κ ) is a constant number, an equivalent form of (33) follows as maximizê wŵ T (P TP • Φ κ )ŵ subject to ŵ Ii 2 = 1, i = 1, ..., r.
The Hadamard productP TP •Φ κ preserves the structure from P TP , or equivalently from the clustering set I in the objective function. As a result, the optimization problem (35) boils down to r decoupled optimizations  Since Assumption 3.2 holds trivial for this case, we state the following theorem. Proof: GivenP TP 0 and Φ 0,P TP •Φ κ is positivesemidefinite according to the Schur product theorem. Also, the objective function in (opt1) is a standard Rayleigh quotient for symmetric eigenvalue problem. Therefore, the maximum of objective function in (opt1) is obtained at the largest singular value of Φ κ [Ii,Ii] . Since Φ κ [Ii,Ii] 0, its largest singular value is the same as its largest eigenvalueλ(Φ κ[Ii,Ii] ), and hence the optimum is obtained at its dominant eigenvector v(Φ κ[Ii,Ii] ).
Given that Φ κ[Ii,Ii] is symmetric, the eigenvector of its largest eigenvalue can be computed efficiently by Krylov subspace-based techniques, e.g. Lanczos algorithm [28] in this case. Solving (opt1) for i = 1, ..., r, therefore, requires a worst case complexity of O(n) in total. Onceŵ is solved from Theorem 5.1, one can then choose any w that satisfies (34), and this w would serve as global optimum for (33).

B. Case II: Unstable Open-Loop
For an unstable open-loop system (1), w solved from Theorem 5.1 is not guaranteed to conform to Assumption 3.2, as a result of whichĝ(s) can become unstable. To avoid such a hazardous situation, we add an extra penalty on (opt1) to restrict w under Assumption 3.2. This penalty term is formulated as follows.
Recall the reduced-order ARE (8). By pre-and postmultiplying it with P T and P , and after a few calculations, we get A TX +XA +Q −XGX = 0, witĥ where U is the complement of P . The ARE above implies that our proposed controllerK = R −1 B TX is equivalent to an LQR problem min ∞ 0 (x TQ x + u T Ru)dt. Denote the eigenvalue decomposition of A as whereΛ 0 contains all the unstable eigenvalues of A. Applying the transformation z = [VV ] −1 x, one can write The matrixV TQV corresponds to the LQR weight for the unstable dynamics of the state trajectories. Therefore, to assure that the unstable modes of A are indeed penalized byQ, we consider finding w that maximizes tr(V TQV ). Given that the matrix U T U A TX +XAU T U in (36) is bounded but not sign-definite, we consider maximizing tr(V TQV ) by simply maximizing its first part as subject to ŵ Ii 2 = 1, i = 1, ..., r. (opt2) Combining (opt2) with the main optimization (opt1), one can formulate the weight design for the unstable case as where ρ > 0 is a penalty factor. In (opt3), the objective function from (opt1) is squared to match with the order of fV . This optimization problem is in the form of a fourthorder sum of squares (SOS) over r > 1 sphere constraints, for which finding even a local optimal is very difficult. One way to bypass this can be to approximate matrix P TQ P in fV as P TQ P •P TP , meaning only to retain the block-diagonal component of Q only. In this way, (opt3) is reduced to r SOS sub-problems with one sphere constraint for each as, for i = 1, ..., r. While this approximation will follow naturally if Q is block-diagonal, the upshot is that the closed-loop performance of the projected system may suffer if Q has dominant off-block-diagonals. In practical networks, however, it is quite common to simply minimize the energy of a node itself, or the energy within a cluster, which implies that Q is very commonly a diagonal or block diagonal matrix. In fact, Q would indeed be preferred as block-diagonal for the scenario described in this section since network operators will always try to discourage closed-loop coupling of their own cluster nodes with other clusters. In those cases, (opt3) and (opt*) become equivalent problems, yielding the same solution. However, irrespective of whether Q is block-diagonal or not, the following theorem shows that the solution of (opt*) will satisfy Assumption 3.2. Theorem 5.2: Givenŵ Ii , i = 1, ..., r solved from (opt*) with ρ > 0, Assumption 3.2 holds for any w satisfying (34).
The following lemma shows the performance trade-off due to the approximation in (opt*) when Q is not block-diagonal. Lemma 5.3: Let the maximum of (opt3) be J 1 . The maximum J 2 of (opt*) satisfies where J e = max i=1,..,r j,l=1,...,n r k=1, k =i ρ Q Ii(j),I k 1 V l,:V T I k ,: 1 .
Proof: The proof follows directly from the Gershgorin circle theorem [28], and is shown in the Appendix.
We next present the solution for (opt*). Since (opt*) is a set of r decoupled problems, we illustrate the solution for just one cluster I = {I 1 }. This will also allow us to drop the subscripts in all the variables used in (opt*), making the notations easier to follow. We define a fourth-order tensor F ∈ R n×n×n×n as where S denotes the product matrix S =VV T . After a few manipulations, it can be shown that (opt*) is equivalent to the following problem where ⊙ denotes the element-wise product. It has been studied in [13] that such a polynomial optimization is equivalent to finding the largest Z-eigenvalue of F , if F is super-symmetric. From the definition in [13], a super-symmetric tensor is one whose entries are invariant to any permutation to the index, i.e. F i,j,k,l = ... = F l,k,j,i , which fails for (39) as F i,j,k,l = F i,k,j,l . However, note that although F is not super symmetric, F ⊙ (ŵ ⊗ŵ ⊗ŵ ⊗ŵ) is a one-dimensional polynomial which is invariant to any index permutations. 2 Following this logic, we rewrite the objective function in (opt**) as follows. Proposition 5.4: Given the fourth-order tensor F specified by (39), the polynomial F ⊙ (ŵ ⊗ŵ ⊗ŵ ⊗ŵ) is identical to where F s is a super-symmetric tensor specified by for i, j, k, l = 1, ..., n, and F s ∈ R n 2 ×n 2 , the matrix unfolding of F s , can be obtained from F s n(i−1)+k,n(j−1)+l = F s i,j,k,l , i, j, k, l = 1, ..., n.
2 This is analogous to an unsymmetric matrix whose quadratic form is in- The proof is omitted as the equations above can be easily verified by matching the coefficients of the polynomials on both sides. In summary, the optimization problem (opt**) can be approached by substituting F with a super-symmetric tensor F s . One can, thereafter, solve (opt**) using techniques developed for Z-eigenvalue problems. We solve (opt**) using the tensor power iteration method [31] in Algorithm 2. The convergence properties of this algorithm can be found in [31]. Due to the super symmetry of F s , the worst case (only one cluster) complexity for each iteration of Algorithm 2 is O(n 4 ). Although this computation cost is expensive, the algorithm can be easily parallelized, and is easier to implement than O(n 3 ) full-order LQR as the memory required is only O(n 2 ). Moreover, the value of n for Algorithm 2 scales down as the number of clusters increases. Initialization: Compute the dominant eigenvector of F s asv(F s ), then choose the initial vector v 0 as v 0 =v(unvec(v(F s ))); 10 end 11ŵ Ii = v k ; 12 end 13 Construct w and then P by (34); Output: P

C. Optimizing (RL-2) with respect to both I and w
The designs proposed in Section IV and this section can be combined to optimize (RL-2) as a function of both I and w iteratively. In this case, one would start with an arbitrarily chosen w, and minimize ξ κ with respect to I using Algorithm 1. Say, the optimal cluster set is given as I * . Thereafter, one would fix I to I * , and minimize ξ with respect to w using Theorem 5.1 or Algorithm 2 depending on whether (1) is stable, and so on. The resulting algorithm is shown in Algorithm 3.

VI. NUMERICAL EXAMPLES
To verify our proposed algorithms, we use a consensus network model defined over a 500-node graph G 500 . The graph is randomly generated with 0.5 overall probability for edge attachment, and 6 spatial clusters with a proportion  Solve I k from Algorithm 1 by (Ψ k−1 , w k−1 , r) ; 6 Update w k from Theorem 5.1 or Algorithm 2 by I k ; 7 k = k + 1 ; 8 end Output: P of 100 for the number of edges within clusters versus the number of edges across clusters. We also apply a random weight 1 ≤ M i,i ≤ 2 on each node. The resulting state matrix A follows the expression (40) presented in Appendix A. The disturbance is assumed to enter from the 364 th node, i.e. B d equals to the 364 th column of I 500 . We assume B = R = I 500 and two choices of Q as: a scaled identity matrix Q 1 = 1000 × I 500 , and Q 2 = [L(G 500 )] 2 , which is the square of the unweighted Laplacian matrix of G 500 . Both Q 1 and Q 2 satisfy Assumption 2.1. For this simulation example, solvability of (8) and stability ofĝ(s) are guaranteed by Theorem A.2 and Theorem A.5 in Appendix A, respectively.  [16] Φo a w r Coherency clustering [3] Φc b 1n r a Let v T c be the complement ofv T , and then Φo : where v i is the right eigenvector of λ i .

A. Cluster Design
We start by finding the closed-loop clustering set I with respect to a fixed weight w =v =v(A) with r = 6. For comparison, we also apply two other popular clustering algorithms, namely, H 2 open-loop clustering proposed in [16], and coherency-based clustering proposed in [3]. Both of these clustering algorithms can be transformed into Algorithm 1, with their equivalent inputs as summarized in Table I. Note that these two algorithms capture only the open-loop characteristics of the network, and hence do not depend on the choice of Q and R. Fig. 4a shows that the clusters identified by coherency based clustering closely resemble the spatial clusters of the open-loop network except for a few discrepancies. For example, two distant groups of nodes are assigned to the same cluster shown in yellow. These discrepancies arise from the fact that the spatial clusters are only based on the edge-weights (that model geographical distance between two We also illustrate the effectiveness of H 2 closed-loop clustering with respect to the number of clusters r. As evident from the design, the error between the transfer matrices in (11) and (12) will be minimal when r = n, and will degrade with decreasing r while improving tractability of the design. We vary r from 1 to 500, and calculate the ratio resulting from the three clustering algorithms. The results are shown in Fig. 5. For both Q 1 and Q 2 the closed-loop clustering outperforms the other two methods in approaching the H 2 performance of g(s). Therefore, even for very small values of r, the projected controller achieves significantly close H 2 performance as the full-order LQR controller. In terms of implementation, the projected controller needs far less number of communication links than a full-order standard LQR as well as a full-order L 1 sparsity-promoting LQR [14]. For example, for this system a standard LQR would require 500 2 = 124750 links. Meanwhile as shown in Fig. 6, a sparsity promoting LQR requires from 3104 to as many as 21325 links to retain a performance loss under 5%. By choosing r ≤ 9, the similar performance loss can be maintained by our design using at most 536 links.
It is also noted that the H 2 closed-loop clusters do not need to strictly follow the spatial geometric clustering of the network. For example, in both Figures (4c) and (4d), a cluster can be one single node as shown by the red, or can be scattered over the network such as yellow. In practice, this means that to implement the proposed control law, nodes from different geographical locations may need to be part of the same cluster for the closed-loop model, i.e., nodes that belong to two different spatial clusters in open-loop may need to collaborate and send their states to a common coordinator. The assignment, therefore, encourages system-wide participation from nodes at various corners of G 500 for implementing the controller.

B. Weight Design
We next apply Algorithm 2 on G 500 to find the optimal projection weight w while fixing the clusters to those obtained from coherency. These clusters as previously shown in  Fig. 7. It can be seen that the weight vectors from Q 1 and Q 2 are very different thanv or between themselves. On the other hand, both of these weights at the 364 th node, i.e., the node where the disturbance enters, show a sudden jump in magnitude from the rest of the nodes. To verify the closed-loop performance, we construct the P matrices using these two vectors, and summarize the error ratios with some design parameters in Table II. As expected, by applying the weight design, the closed-loop errors as shown in Table II are significantly reduced from w =v. Despite the fact that these  two errors are much larger than what we get from closedloop clustering (which yields an error of 0.68%), the weight design still grants us with significant improvement over the hard constraint on I.
Finally, we compare the closed-loop performance of the iterative Algorithm 3 (where both w and I are free) with Algorithm 1 in Fig. 8. The comparison is shown for Q 1 and only r ≤ 6 as the error ratio already becomes under 1% after r = 6. For this example, it is worth mentioning that Algorithm 3 turns out to be surprisingly efficient as it converges right after the first iteration. In this sense, the iterative process reduces to a single weight design after the clustering. Fig. 8 verifies that Algorithm 3 achieves better matching between g(s) and g(s) than Algorithm 1.

C. Scalability Results
To verify the scalability of Algorithm 1, we increase the size of the network, and compare the computation time with that of solving a full-order LQR controller. We let r = 5 and κ = 5 for computing Φ  When n = 8000, the full-order LQR is already beyond the capability of our computation facility. The reduced-order LQR design, however, requires remarkably less computation time, while still providing a close performance match to the fullorder LQR controller.

VII. CONCLUSION
In this paper we developed a set of projection-based algorithms that improve the dynamic response of large-scale network systems with reduced-order LQR controllers. The advantage of these reduced-order controllers is that they are structured, and significantly easier to design and implement compared to regular full-order LQR controllers. Our future work will be to address the robustness of this approach to communication delays, to exploit additional input-output properties such as passivity to further improve performance, and to inspect the influence of network heterogeneity on clustering.

APPENDIX A: SPECIAL CASE OF CONSENSUS NETWORK
The well-posedness conditions from Section III.A can be relaxed if system (1) is a consensus network. Consensus is commonly used to model the dynamic behavior of many practical networks such as social networks, power networks, and wireless networks. For the same network graph G = (V, E) as in (1), we suppose n = n s . Each node has a real-valued node-weight m i > 0, and each edge ij ∈ E has a real-valued edge-weight a ij = a ji > 0. A consensus network can then be modeled in its simplest form by letting A ij = m −1 i a ij and A ii = − m −1 i j∈Ni a ij in (2). The overall system writes aṡ which we will simply denote asv. The state matrixÃ from the reduced-order system (6) satisfies the following property.
Proof: Denote the complement of the projection matrix P by U , such that D = [P T U T ] T is unitary, i.e. DD T = I n . Then,Ã = P AP T is a leading principal of the matrix which is similar to A. Given the leading principal of a symmetric positive-definite (or-semidefinite) matrix still being positive-definite (or-semidefinite), it follows thatÃ 0. Moreover, from Definition 2.2, P is defined over w such that P T P w = w. Then if w ∈ span(v),Ã preserves the zero eigenvalue from A sinceÃPv = P Av = 0. Therefore,Ã is negative-semidefinite if w ∈ span(v), or negative-definite otherwise.
1) Existence Condition: The next theorem shows that the existence ofX for the reduced-order ARE (8)  3) Stability condition: We next state two sufficient conditions that guarantee the stability ofĝ(s) for consensus networks.
Theorem A.4: The systemĝ(s) is stable if w =v, and I is an almost equitable partition [25] of the graph G, which means for k = l, the edge weight a ij is equal for all i ∈ I k and j ∈ I l .
Proof From [25], if G admits an almost equitable partition I, the corresponding matrix P with w =v will satisfy A T P T = P T N for some N ∈ R r×r . As a result, P AU T = 0, and thusĝ (s) is stable sinceÃ −GX and U AU T are both Hurwitz given w =v.
Theorem A.5: Assume that B is a square invertible matrix. Then,ĝ(s) is stable if G is similar to αI n for some α > 0.
Proof: Denote the right eigenspace of G by V . If G is similar to αI n , it follows that V T (GX +XG)V = 2αV TX V . This means that GX +XG is positive-semidefinite. From the matrix majorization property [26], we also have 2λ(A − GX) ≤ z T (A + A T − GX −XG)z hold for any non-zero vector z, with the RHS being nonpositive given A + A T 0 and −GX −XG 0. Hence forλ(A − GX) to be strictly negative, z T Az = 0 should contradict z T GXz = 0. To prove the contradiction, we assume z T Az = 0, namely z ∈ span(v). Note that z T GXz = 0 holds if and only if either one of the following three conditions is satisfied: (1) z ∈ ker(P ), (2) z ∈ ker(G), or (3) P z ∈ ker(X). First of all, Assumption 3.2 requires w T Iiv Ii , i = 1, ..., r, which means Pv = 0 and thus z ∈ ker(P ). Given that (A, BR − 1 2 ) is stabilizable, we havev T Gv = 0, i.e. z ∈ ker(G). To this end, the last condition remaining is z ∈ span(v) ⇔ P z ∈ ker(X) to complete the contradiction. By assuming a non-zero vector v ∈ ker(X), pre-and postmultiplying (8) with It can be easily verified that ker(X) is anÃ-invariant subspace contained in the null-space ofQ. Given z ∈ span(v), from the proof of Proposition A.1 we know that P z isÃ-invariant, i.e., P z is an eigenvector ofÃ, when w ∈ span(v). As a result, z T P TQ P z = z T Qz = 0 since (Q T 2 , A) is detectable. This verifies that P z ∈ ker(Q), which proves P z ∈ ker(X). Therefore, we conclude that z T Az = 0 contradicts z T (GX + XG)z = 0, and thusĝ(s) is stable.