1

Optimal Control of Large-Scale Networks using
Clustering Based Projections

arXiv:1609.05265v3 [] 4 Oct 2017

Nan Xue, Student Member, IEEE, Aranya Chakrabortty, Senior Member, IEEE

Abstract—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.
Index Terms—Clustering, Large-scale networks, Projection,
LQR, H2 performance.

I. I NTRODUCTION
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. Xue and A. Chakrabortty are with the Department of Electrical and
Computer Engineering, North Carolina State University, Raleigh, NC, 27695
USA, e-mail: nxue@ncsu.edu, achakra2@ncsu.edu
The work is supported partly by the US National Science Foundation (NSF)
under grant ECCS 1054394.

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 H2 norm 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

2

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
1n
column vector of size n with all 1 entries
Ik
identity matrix of size k
Mi,j
the (i, j)th entry of a matrix M
diag(m)
diagonal matrix with vector m on its
principal diagonal
M ⊗N
Kronecker product of M and N
M ◦N
Hadamard product of M and N
tr(M )
trace operation on a matrix M
kM kF
Frobenius pnorm of a matrix M , i.e.
kM kF = tr(M M T )
ker(M )
kernel of a matrix M
σ̄(M ), λ̄(M ) largest singular value, or eigenvalue with
largest real part of a matrix M

σ(M ), λ(M )

smallest singular value, or eigenvalue
with smallest real part of a matrix M
v̄(M )
right eigenvector of λ̄(M )
Given a matrix M = [m1 , ..., mn ] ∈ Rn×n , its vector form is defined by vec(M ) = [mT1 , ..., mTn ]T , with
the inverse operation defined by unvec(vec(M )) = M . A
−1
transfer matrix is defined as g(s) = C(sI
 − A) B + D,
A B
. We rewith a realization form of g(s) =
C D
fer to g(s) as stable if A is Hurwitz, and unstable otherwise. Furthermore, the H2 and H∞ norms of a stable
qR transfer matrix g(s)qare R defined by kg(s)kH2 =
∞
∞
1
∗
∗
−∞ tr[g (t)g(t)]dt =
2π −∞ tr[g (jω)g(jω)]dω and
kg(s)kH∞ = supω σ̄[g(jω)].
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 Ni = {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. P ROBLEM F ORMULATION
Consider a general LTI system of the form
(
ẋ(t) = Ax(t) + Bu(t) + Bd d(t), x(0) = x0
y(t) = Cx(t)

,

(1)

where x(t) ∈ Rn , u(t) ∈ Rm and y(t) ∈ Rp represent the
vector of state, control and output variables respectively, and
d(t) ∈ Rnb is a disturbance entering into the system. We
assume (1) to be defined over a network of ns ≤ n interconnected subsystems, with their network topology represented by
a connected graph G = (V, E), V = {1, ..., ns }. The dynamics
of each subsystem can be written as

P
ẋi (t)=Aii xi (t)+
Aij xj (t)+Bi ui (t)+Bdi di (t)
j∈Ni
(2)
y (t) = C x (t)
i
i i

i = 1, ..., ns , where Aij , Bi , Bdi and Ci are submatrices
with compatible dimensions from A, B = diag(B1 , ..., Bns ),
Bd = diag(Bd1 , ..., Bdns ) and C = diag(C1 , ..., Cns ). Notice
that the dimension of ui (t) can be zero, meaning that the ith
subsystem can have no input.
In this paper, we consider an LQR design for (1), and
assume C = In for full-state feedback. Given two real-valued
matrices Q = QT  0 and R = RT ≻ 0, the LQR problem
is posed as finding a feedback law u(t) = −Kx(t) such that
the cost function
Z ∞
J :=
[xT (t)Qx(t) + uT (t)Ru(t)]dt
(3)
0

is minimized. The expression (3), also known as the infinitehorizon continuous-time LQR, can be solved by the following
algebraic Riccati equation (ARE)
AT X + XA + Q − XGX = 0,

(4)

3

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.
T
1
Assumption 2.1: (Q 2 , A) is observable, and (A, BR− 2 ) is
stabilizable.
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(n3 ) 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 ∈ Rn , define r non-empty, distinct, and
non-overlapping subsets of the state index set Vs = {1, ..., n},
respectively denoted as I = {I1 , ..., Ir }, such that I1 ∪ ... ∪
Ir = Vs . A clustering-based projection matrix P ∈ Rr×n is
defined as
( w
j
j ∈ Ii
,
(5)
Pi,j := kwIi k2
0
otherwise
where wIi = [wIi {1} , ..., wIi {|Ii |c } ]T is non-zero, and Ii {j}
denotes the j th element in the set Ii . The matrix P has the
following three properties:

2) Reduced-order LQR design: We similarly project the
LQR parameters by Q̃ = P QP T ∈ Rr×r , and let R̃ = R
such that G̃ := B̃ R̃−1 B̃ = P GP T ∈ Rr×r . An LQR problem
for the reduced-order model (6) is then posed as to minimize
Z ∞
[x̃T (t)Q̃x̃(t) + ũT (t)R̃ũ(t)]dt
(7)
J˜ :=
0

with respect to ũ(t) = −K̃ x̃(t). Here the feedback matrix
K̃ = R̃−1 B̃ T X̃ corresponds to the solution X̃ ∈ Rr×r of the
reduced-order ARE of (7), which is written as
ÃT X̃ + X̃ Ã + Q̃ − X̃ G̃X̃ = 0.

(8)

3) Inverse projection to original coordinates: The solution
X̃ from (8) is projected back to the original coordinates
through the inverse projection
X̂ = P T X̃P.

(9)

This projected controller can then be implemented in the fullorder model (1) using u = −R−1 B T X̂x, which implies that
the effective feedback gain matrix is
K̂ = R−1 B T X̂.

(10)

B. Problem Statement
The controller K̂ is dependent on the projection P through
equations (6), (7), and (9). The choice of P is guided in the
following way. Consider
g(s) := (sIn − A + BK)−1 Bd ,

(11)

which is the closed-loop transfer matrix from d to x for (1)
with full-order LQR. Similarly, consider
ĝ(s) := (sIn − A + B K̂)−1 Bd ,

(12)

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
The construction of P is shown by the following example.

T and a non-zero vector w such that the corresponding projection
Example 1: Let w = 1 1 1 2 1 1 1 1 1 1 ,
matrix P solves the model matching problem1
I1 = {1, 2}, I2 = {3, 4, 5} and I3 = {6, 7, 8, 9, 10}. Then,
minimize kg(s) − ĝ(s)kH2 .

 1
(13)
P
√
√1
0
0
0
0
0
0
0
0
2
2

√1
√2
√1
However, finding an exact solution for this optimization
0
0
0
0
0
0
P= 0
.
6
6
6
problem is intractable given that the objective function is an
√1
√1
√1
√1
√1
0
0
0
0
0
5
5
5
5
5
implicit and non-convex function of P , and also because P
Given P defined over any clustering set I and weight vector is defined over a combinatorial structure. Our main contribuw, the control inversion strategy for the LQR problem (3) is tion, therefore, is finding a tractable relaxation for (13) as a
quadratic function of P , and thereafter designing P to solve
then composed of the following three steps.
1) Projection to reduced-order system: Using projection P , the relaxed problem. We make the following assumption so
that g(s) and ĝ(s) both have minimal realization.
we first construct a reduced-order model
Assumption 2.3: The pair (A, Bd ) is controllable.
(
˙
x̃(t)
= Ãx̃(t) + B̃ ũ(t) + B̃d d(t), x̃(0) = x̃0
Solution Strategy: The outline of our solution strategy is as
, (6)
follows. In Section III, we derive an upper bound relaxation for
ỹ(t) = C̃ x̃(t)
(13) such that its objective function is quadratic in P . Ideally
where Ã := P AP T ∈ Rr×r , B̃ := P B ∈ Rr×m , B̃d :=
1 The initial condition for the reduced-order model (6) does not need to be
P Bd ∈ Rr×nb and C̃ := P CP T = Ir . For this system x̃ ∈ Rr related to that of the full-order model (1). Our goal is to compare (11) and
is the state and ũ ∈ Rm is the control input.
(12), both of which have zero initial conditions.
•

It is row orthonormal, i.e. P P T = Ir ;
Image of P T P lies in the span of w, i.e. P T P w = w;
Given v 6= 0, P v = 0 only if wITi vIi = 0, i = 1, ..., r.

4

speaking, one can solve for P from this relaxed problem.
The computational complexity for constructing the objective
function is, however, O(n3 ) 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.
E"F&'()2?G
A#@4&"=23"4'0)
>#092?

E"F&'()2???
!"#$%$&'()
*!+,-.

ङ ࢔ࣄ

૛

ङ ࢔ࣄ࢘࢑

this process no coordinator will be able to infer individual
state measurements from other clusters. Finally, the coordinator of Ii computes the control vector ui by taking linear
combinations of the elements of X̃P x. The linear combination
follows from u = −R−1 B T P T X̃P x. The individual elements
of ui are broadcast to the respective input actuators in Ii .
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 xi and a scalar control input ui , the two-layer
control implementation will result
in a much sparser commu
nication topology with n + r2 bidirectionallinks compared to
an LQR controller which would require n2 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.

E"F&'()2??

A()&=(#2?)B"=4'()

?&"=$&'B"
>#092???

*8. C *:D.

E"F&'()2G
/"'01&23"4'0)5
617 89:2*4&$;#"2(<"),#((<.
(=2>#092??2*@)4&$;#"2(<"),#((<.

ङ ࢔

ङ ࢔૝

ङ ࢘࢔૛ ൅ ङ ࢘૜

Figure 1: Step-by-step execution of the proposed designs

C. Benefits of Control Inversion
An important point to note is that the physical meaning of
the state x̃(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
controller X̃ 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(rn2 ),
while that for solving this ARE is O(r3 ). The computational
complexity required to design P through Sections III, IV and
V will be shown to be simpler than the O(n3 ) 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 matrix X̂ =
P T X̃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 T X̂x(t) follows
three steps. First, a coordinator is assigned to each cluster Ii ,
which collects the measurements of all the states belonging
to that cluster. Each coordinator then computes the weighted
averaged state Pi,: x for its cluster, i = 1, 2, .., r. Next, the
coordinators exchange these weighted averages, and each of
them compute the r-dimensional vector X̃P x. Note that in

III. R ELAXATIONS

FOR

M ODEL M ATCHING

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 controller K̂
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 of K̂ is equivalent to two conditions namely, if the reduced-order ARE (8) admits a solution X̃,
and if ĝ(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 soT
1
lution X̃ ≻ 0 if (Ã, G̃ 2 ) is stabilizable and (Q̃ 2 , Ã) is
observable [7]. However, unlike conventional model reduction
techniques that can utilize unstructured projections to preserve
the exact stabilizability and observability properties from
1
T
(A, G 2 ) and (Q 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 solution X̃, and therefore K̂ may
not exist. To bypass this problem, we modify the definitions
of Q̃ and G̃ by using a constant shift. That is, in case the pair
T
1
(P AP T , P G 2 ) is not stabilizable and/or (Q 2 P T , P AP T ) is
not observable, we let
Q̃ = P QP T + αIr

and/or G̃ = P GP T + αIr

(14)

for a small constant shift α > 0. Using (14), Q̃ and G̃ become
1
positive-definite matrices, which yield (Ã, G̃ 2 ) controllable
T
and/or (Q̃ 2 , Ã) observable. The existence of X̃ and K̂ can
thus be guaranteed. The matching error introduced by this shift
will be discussed in the next subsection.

5

C1

‫ݔ‬ଷ ൅ ʹ‫ݔ‬ସ ൅ ‫ݔ‬ହ

‫ݔ‬ଵ ൅ ‫ݔ‬ଶ
ʹ

‫ ଺ݔ‬൅ ‫ ଻ݔ‬൅ ‫ ଼ݔ‬൅ ‫ݔ‬ଽ ൅ ‫ݔ‬ଵ଴

2

‫ݔ‬ଵ଴

10

‫ݔ‬ଽ

9

‫ݔ‬ଵ

1

‫ݔ‬ଷ

‫ݔ‬ସ

4

3

‫଺ݔ‬

6

༠

C1
‫څ‬
ଶ

C3

‫ݔ‬ହ

1

5
‫଻ݔ‬

‫ל‬

C2

C1

ͷ

C3

‫ݔ‬ଶ

‫څ‬

͸

C2

2

༠
ହ

C

2

5

10
9

(a) Step 1 - state averaging P x

‫ל‬
଺

7

10

8

ଶ‫ל‬
଺

6

3

7

7

‫଼ݔ‬

C2

4

5

6

‫ל‬
଺

3 ༠
༠ ହ
༠
༠
ହ
ହ
ହ

1

4
3

‫څ‬
ଶ

8

9

(b) Step 2 - reduced-order control X̃P x =
[⋆ ◦ ⋄]T

8

(c) Step 3 - control inversion u = −X̂x

Figure 2: Cyber-physical architecture for implementing the feedback controller K̂. Here I1 = {1, 2}, I2 = {3, 4, 5} and I3 = {6, 7, 8, 9, 10}.
Solid lines represent physical connections, while dashed lines represent communication links. C1 , C2 , C3 are coordinators for clusters 1, 2
and 3. For the simplicity of illustration, for this example we assume R−1 B T = In .

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:
λ̄(X) ≤ β(A, G, Q) = λ̄(Dt )

λ̄[(Q + KtT Kt )Dt ]
,
λ(F Dt )

(15)

1

where Kt is any matrix stabilizing A + G 2 Kt , and Dt and F
are positive-definite matrices that satisfy
1

1

(A+G 2 Kt )T Dt + Dt (A+G 2 Kt ) ≤ −F.

(16)

From this lemma, given any stabilizing controller Kt 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 ĝ(s).
3) Stability condition: Before stating the main stability
criterion for ĝ(s), we first make the following assumption on
the projection weight w from Definition 2.2.
Assumption 3.2: The weight vector w satisfies wITi vIi 6= 0,
i = 1, ..., r for any Av = λv, Re(λ) ≥ 0.
Note that from Definition 2.2, if wITi vIi = 0, i = 1, ..., r
then P v = 0. This would imply (A − B K̂)v = Av = λv,
Re(λ) ≥ 0, which means that the unstable eigenvalues of A
are retained in the closed-loop. Assumption 3.2 is made to
avoid this situation. With Assumption 3.2, we next state a
sufficient condition for the stability of ĝ(s).
Theorem 3.3: Given Assumption 2.1 and 3.2, the TFM ĝ(s)
is asymptotically stable if
σ(Q) − 2σ̄(X)σ̄(G)σ̄(E) + σ(X)2 σ(G) > 0,

(17)

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 ĝ(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 ĝ(s) can be verified without knowing X.
Lemma 3.4: The stability condition (17) holds if
σ̄(X̃) <

σ(Q)
− β(A, G, Q),
2σ̄(G)β(A, G, Q)

(18)

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 Kt 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 ĝ(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 kg(s)− ĝ(s)kH2 is an intractable problem. To relax
this problem, we find an upper bound for kg(s) − ĝ(s)kH2 .
We denote the error system by ge (s), which can be written as


Ae Be
ge (s) : = g(s) − ĝ(s) =
Ce 0


A − GX
Bd
(19)
=
A − GX̂ −Bd  .
In
In
0

6



I
I
From a similarity transformation of T = n n and T −1 =
In 0


0
In
, (19) yields
In −In




A−GX̂ G(X̂−X) 0
T Ae T −1 T Be
ge (s)=
=
A−GX
Bd 
Ce T −1
0
In
0
0
−1

= −(sIn − A + GX̂)

−1

GE(sIn − A + GX)

Bd ,

(20)

where E = X − X̂. By taking norms on both sides of (20),
we get
kge (s)kH2 ≤ k(sIn − A + GX̂)−1 GkH∞ kEg(s)kH2 ,
which from the bounded real lemma [7] and the definition of
H2 norm reduces to
1

kge (s)kH2 ≤ γkEΦ 2 kF ,

(21)

where γ is any positive real number such that a real-valued
matrix Γ = ΓT  0 exists and satisfies


Γ(A − GX̂) + (A − GX̂)T Γ ΓG
In

GΓ
−γIn
0  ≺ 0,
In
0
−γIn
R ∞ (A−GX)τ
T
T
1
Bd BdT e(A−GX) τ dτ ≻ 0
and Φ := Φ 2 Φ 2 = 0 e
is the solution of the Lyapunov equation
(A − GX)Φ + Φ(A − GX)T + Bd BdT = 0.

(22)

Inequality (21) shows that kge (s)kH2 is linearly bounded by
1
kEΦ 2 kF . Therefore, one way to solve the original model
matching problem (13) will be to find a P that minimizes
1
kEΦ 2 kF . By doing so, σ̄(E) can also be minimized to some
1
1
extent since σ̄(E) ≤ kEkF ≤ σ̄(Φ− 2 )kEΦ 2 kF , 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
1
1
1
1
from K(A, Q 2 , r) = Q 2 AQ 2 · · · Ar−1 Q 2 . 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
1
(21), and derive a new upper bound on kEΦ 2 kF as an explicit
function of P in the following theorem.
1
Theorem 3.5: Denote ξ = k(In − P T P )Φ 2 kF . The norm
1
of the weighted error EΦ 2 satisfies the inequality
1
2

kEΦ kF ≤ f (ξ) = ǫ1 σ̄(Q)ξ 2 + 2ǫ1 ǫ2 ξ + αǫ1 ǫ3 ,

(23)

1

where ǫ1 =

σ̄(Φ− 2 )
1
σ[Φ− 2

1

1

(A−GX)Φ 2 ]

1

, ǫ2 = β̃ σ̄(A)σ̄(Φ 2 )+ σ̄(QΦ 2 ),

ǫ3 = (β̃ 2 + 1)σ̄(Φ), and β̃ = supP β(Ã, 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 kge (s)kH2 ≤ γf (ξ),

due to which we approach the minimization of kge (s)kH2 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
minimize
P

1

1

ξ = kΦ 2 − P T P Φ 2 kF .

(RL-1)

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 kge (s)kH2 in (13)
will be nearly zero as well following (21) and (23).
1
1
The ideal case ξ = 0 will happen when Φ 2 = P T P Φ 2 ,
1
meaning that Φ 2 is invariant to the mapping P T P . This holds
for the trivial case where r = n and P = In . 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


A
−G
H :=
.
(24)
−Q −AT
The eigenvalues of H are symmetric about the imaginary
axis. Suppose
 H is diagonalizable and that the columns of
Y
the matrix
span the stable invariant subspace of H,
Z 2n×n
i.e.
   
Y
Y
H
=
Λ− ,
(25)
Z
Z

−
where Λ− = diag([λ−
1 , ..., λn ]) 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 [6]

Φ = Y (Y −1 Bd BdT Y −T ◦ C)Y T ,
where C ∈ Rn×n is a Cauchy matrix with
i
h
1
−
,
Ci,j = − λ−
+λ
i
j
1

(26)

and subsequently obtain Φ 2 from the Cholesky decomposition.
1
Therefore, to compute Φ and then Φ 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(n3 ) complexity for both computation and memory
[28]. This may defeat the purpose of our design since we

7

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.
Definition 3.6: Define Φκ ∈ Rn×n as
1

T

Φκ := Φκ2 Φκ2 = Y1 (Ω1 Bd BdT ΩT1 ◦ C1,1 )Y1T ,

(27)

−1
where Y1 =Y:,1:κ , Ω1 =Y1:κ,:
and C1,1 =C1:κ,1:κ are respectively
the κ-dimensional partitions of Y , Y −1 and C.
By definition of Φκ , one only needs to compute the first κ
−
eigenvalues λ−
1 , ..., λκ of H, and the Y1 component of the first
κ eigenvectors. Ω1 can be approximated by the pseudo-inverse
of Y1 . 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

minimize
P

1
2

1
2

ξκ = kΦκ − P T P Φκ kF ,

(RL-2)

1

where Φκ2 ∈ Rn×κ from Definition 3.6 can be computed as
1
1
Φκ2 = Y1 (Ω1 Bd BdT ΩT1 ◦ C1,1 ) 2 . 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 Bd has a unitary norm. The minimum
ξ ∗ from (RL-1) and ξκ∗ from (RL-2) satisfy
v
u
n
X
u
1
∗
∗
(28)
− −.
ξ − ξκ ≤ tη 2 nb
2λi
i=κ+1
From Lemma 3.7, the optimality gap between
(RL-1) and (RL-2) will be negligible when the error
q
Pn
η 2 nb i=κ+1 − 2λ1− is kept small. In practice, this
i
situation happens when there exists κ (κ ≪ n) dominant
eigenvalues in the Hamiltonian matrix, i.e., H has the
following spectral gap
0<

|λ−
1|

< ... <

|λ−
κ|

≪

|λ−
κ+1 |

< ... <

|λ−
n |.

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. D ESIGN I: C LUSTER D ESIGN
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
matrix P̄ as
(
1
j ∈ Ii
, i = 1, ..., r.
(29)
P̄i,j := kwIi k2
0
otherwise
From Definition 2.2, P̄ satisfies P = P̄ W , where W =
diag(w). With this notation, we have
1

1

Φκ2 − P T P Φκ2 = W (Ψ − P̄ T P̄ W 2 Ψ),

(30)
1

where Ψ = [ψ1 , ..., ψn ]T denotes the matrix W −1 Φκ2 . Therefore, the entries of the matrix P̄ T P̄ W 2 can be found by
( w2
k
j ∈ Ii & k ∈ Ii
2
T
2
[P̄ P̄ W ]j,k = kwIi k2
0
otherwise
for i = 1, ..., r. Thus the j th row of the matrix P̄ T P̄ W 2 Ψ
can be written as
P
2 T
k∈Ii wk ψk
T
2
T
[P̄ P̄ W Ψ]j,: = ci = P
(31)
2 ,
k∈Ii wk

for j ∈ Ii . It is clear from above that for all the index j that
are assigned to the same cluster Ii , cTi is a weighted average
(or a weighted centroid) of the row vectors ψjT . Moreover, the
matrix P̄ T P̄ 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

ξκ2 =

n
X
j=1

wj2 kψj − ci k22 .

(32)

The optimization problem in (32) is in the same form as
a weighted k-means optimization, which minimizes the Euclidean distance weighted by wj2 between each data point
ψj and its centroid ci . 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
1
clustering for (32) in Algorithm 1. With Φκ2 ∈ Rn×κ 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 get K̂. The

8

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(n2 r) + O(r3 ), which includes
construction of the objective function for (RL-2), execution
of Algorithm 1, computation of reduced-order matrices triple
(Ã, Q̃, G̃), and solving the reduced-order LQR (7), respectively. This complexity can be at most O(n2 r) if κ ≤ r,
which is more tractable compared to the O(n3 ) complexity
of full-order LQR, especially when r and κ are sufficiently
small.
Algorithm 1: H2 closed-loop clustering
1

1

Input : Φκ2 , w and r
Initialization: Assign r random rows from
1

2

Ψ = W −1 Φκ2 = [ψ1 , ..., ψn ]T as the initial centroids
c01 , ..., c0r ;
Find initial clustering sets
I 0 = {j → Ii0 | argmin wj2 kψj − c0i k22 , j = 1, ..., n};
i=1,...,r

3
4
5
6

Update the centroids: c0i =

P

j∈I 0
i

P

8
9

wj2

j∈I 0
i

, i = 1, ..., r ;

k = 1;
while I k−1 6= I k or within maximum iterations do
Update clustering sets I k = {j →
Iik | argmin wj2 kψj − cik−1 k22 , j = 1, ..., n};
i=1,...,r

7

wj2 ψj

Update the centroids: cki =
k =k+1 ;
end
Output: I = I k

P

j∈I k
i

P

wj2 ψj

j∈I k
i

wj2

, i = 1, ..., r ;

V. D ESIGN II: W EIGHT D ESIGN
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
1

minimize
w

1

ξκ = kΦκ2 − P T P Φκ2 kF .

(33)

To solve (33), we consider a binary projection matrix P̂ as
(
1
j ∈ Ii
P̂ (i, j) :=
.
0
otherwise
As can be verified from Definition 2.2, P̂ here also satisfies
P = P̂ Ŵ with Ŵ = diag(ŵ), where this ŵ is defined by
wj
, j ∈ Ii ,
(34)
ŵj =
kwIi k2

such that ŵITi ŵIi = 1, for i = 1, ..., r. Using these notations,
we can rewrite the objective function in (33) as
ξκ2 = tr(Φκ − P T P Φκ ) = tr(Φκ ) − tr(Ŵ P̂ T P̂ Ŵ Φκ )
= tr(Φκ ) − ŵT (P̂ T P̂ ◦ Φκ )ŵ.

Since tr(Φκ ) is a constant number, an equivalent form of (33)
follows as
maximize

ŵT (P̂ T P̂ ◦ Φκ )ŵ

subject to

kŵIi k2 = 1, i = 1, ..., r.

ŵ

(35)

The Hadamard product P̂ T P̂ ◦Φκ preserves the structure from
P̂ T P̂ , or equivalently from the clustering set I in the objective
function. As a result, the optimization problem (35) boils down
to r decoupled optimizations
maximize

ŵITi Φκ[Ii ,Ii ] ŵIi

subject to

kŵIi k2 = 1,

ŵIi

(opt1)

for i = 1, ..., r, where Φκ[Ii ,Ii ] denotes the submatrix of Φ
corresponding to the indices in Ii . This decoupling can be
illustrated by the same example we used before.
Example 2: (Continued) Consider the sets I1 , I2 and I3 ,
and the matrix


1 1 0 0 0 0 0 0 0 0
P̂ = 0 0 1 1 1 0 0 0 0 0 .
0 0 0 0 0 1 1 1 1 1
The objective function in (35) can be block diagonalized as



Φκ[1:2,1:2]
0
0
 ŵ=
0
Φκ[3:5,3:5]
0
ŵ (P̂ P̂ ◦Φκ )ŵ=ŵ 
0
0
Φκ[6:10,6:10]

 T
ŵ1:2 Φκ[1:2,1:2] ŵ1:2
0
0
T
.

0
ŵ3:5 Φκ[3:5,3:5] ŵ3:5
0
T
0
0
ŵ6:10
Φκ[6:10,6:10] ŵ6:10
T

T

T

Since Assumption 3.2 holds trivial for this case, we state
the following theorem.

9

(a) Spatial groups partitioned by the closeness in terms of geometric distances

(b) Node clusters resulting from cluster
design as in Section IV

(c) Imposing projection weights while fixing the clustering sets to the spatial clusters

Figure 3: Illustration of cluster design and weight design.

Theorem 5.1: The global optimum for (opt1) is obtained at
ŵIi = v̄(Φκ[Ii ,Ii ]), for i = 1, .., r.
Proof: Given P̂ T P̂  0 and Φ  0, P̂ T P̂ ◦Φκ 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 ŵ 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 ĝ(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 AT X̂ + X̂A + Q̂ − X̂GX̂ = 0, with
Q̂ := P T Q̃P − (U T U AT X̂ + X̂AU T U ),

(36)

where U is the complement of P . The ARE above implies
−1 T
that our proposed controller
R ∞ K̂T = R BT X̂ is equivalent to
an LQR problem min 0 (x Q̂x + u Ru)dt. Denote the
eigenvalue decomposition of A as



 
 Λ
,
(37)
A V V̄ = V V̄
Λ̄

where Λ̄  0 contains all the unstable eigenvalues of A.
Applying the transformation z = [V V̄ ]−1 x, one can write

 T
V T Q̂V̄
T
T V Q̂V
z.
x Q̂x = z
V̄ T Q̂V V̄ T Q̂V̄

The matrix V̄ T Q̂V̄ 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 by Q̂, we
consider finding w that maximizes tr(V̄ T Q̂V̄ ). Given that
the matrix U T U AT X̂ + X̂AU T U in (36) is bounded but not
sign-definite, we consider maximizing tr(V̄ T Q̂V̄ ) by simply
maximizing its first part as
maximize

fV̄ = tr(V̄ T P T Q̃P V̄ )

subject to

kŵIi k2 = 1, i = 1, ..., r.

ŵ

(opt2)

Combining (opt2) with the main optimization (opt1), one can
formulate the weight design for the unstable case as
maximize

[ŵT (P̂ T P̂ ◦ Φκ )ŵ]2 + ρ · fV̄

subject to

kŵIi k2 = 1, i = 1, ..., r,

ŵ

(opt3)

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 T Q̃P in fV̄
as P T Q̃P ◦ P̂ T P̂ , 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,
maximize
ŵ

subject to

(ŵIi ⊗ ŵIi )T (Φκ[Ii ,Ii ] ⊗ Φκ[Ii ,Ii ]
+ ρQIi ,Ii ⊗ V̄Ii ,: V̄ITi ,: )(ŵIi ⊗ ŵIi )

(opt*)

kŵIi k2 = 1,

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 ŵIi , i = 1, ..., r solved from (opt*)
with ρ > 0, Assumption 3.2 holds for any w satisfying (34).

10

Proof: Note that V̄:,j , j = 1, ..., nv represents the
eigenvector of the j th unstable eigenvalue. The second part
ofPthe objective function in (opt*) can be rewritten as
v
ρ nj=1
(V̄ITi ,j ŵIi )2 (ŵITi QIi ,Ii ŵIi ). Once (opt*) is solved,
each individual scalar (V̄ITi ,j ŵIi )2 has to be maximized away
from 0, meaning V̄ITi ,j ŵIi 6= 0, i = 1, ..., r, j = 1, ..., nv for
AV̄:,j = λV̄:,j , Re(λ) ≥ 0. This satisfies Assumption 3.2.
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 J1 . The maximum J2 of (opt*) satisfies

where Je =

J 1 − J e ≤ J2 ≤ J 1 + J e ,
(38)
Pr
T
max
k=1, k6=i ρkQIi (j),Ik k1 kV̄l,: V̄Ik ,: k1 .

i=1,..,r
j,l=1,...,n

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 = {I1 }. 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 ∈ Rn×n×n×n as
Fi,j,k,l = Φκ[i,j] Φκ[k,l] + ρQi,j Sk,l , i, j, k, l=1, ..., n. (39)

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(n4 ). Although
this computation cost is expensive, the algorithm can be easily
parallelized, and is easier to implement than O(n3 ) full-order
LQR as the memory required is only O(n2 ). Moreover, the
value of n for Algorithm 2 scales down as the number of
clusters increases.
Algorithm 2: Power iteration for projection weight design
1
2
3
4

T

where S denotes the product matrix S = V̄ V̄ . After a few
manipulations, it can be shown that (opt*) is equivalent to the
following problem
maximize

F ⊙ (ŵ ⊗ ŵ ⊗ ŵ ⊗ ŵ)

subject to

kŵk2 = 1,

ŵ

5
6

(opt**)
7

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. Fi,j,k,l = ... = Fl,k,j,i , which fails for (39) as Fi,j,k,l 6=
Fi,k,j,l . However, note that although F is not super symmetric,
F ⊙ (ŵ ⊗ ŵ ⊗ ŵ ⊗ ŵ) 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 ⊙ (ŵ ⊗ ŵ ⊗ ŵ ⊗ ŵ) is identical to

8
9
10
11
12
13

Input : Φκ , Q, V̄ , I, ρ and δ
Partition Φκ[Ii ,Ii ] , QIi ,Ii and V̄Ii based on I;
for i = 1, ..., r do
Find F s corresponding to Ii by Proposition 5.4;
Initialization: Compute the dominant eigenvector of
F s as v̄(F s ), then choose the initial vector v 0 as
v 0 = v̄(unvec(v̄(F s )));
k = 1;
k
k T
) F s (v k ⊗v k )
while (vk−1(v⊗v⊗v
k−1 )T F s (v k−1 ⊗v k−1 ) − 1 > δ or within
maximum iterations do
v k = unvec(F s (v k−1 ⊗ v k−1 ))v k−1 ;
k
v k = kvvk k2 ;
k = k + 1;
end
ŵIi = v k ;
end
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
F ⊙ (ŵ ⊗ ŵ ⊗ ŵ ⊗ ŵ) = (ŵ ⊗ ŵ) F (ŵ ⊗ ŵ),
w iteratively. In this case, one would start with an arbitrarily
chosen w, and minimize ξκ with respect to I using Algorithm
where F s is a super-symmetric tensor specified by
1. Say, the optimal cluster set is given as I ∗ . Thereafter,
1
s
one would fix I to I ∗ , and minimize ξ with respect to w
Fi,j,k,l = (Φκ[i,j] Φκ[k,l] + Φκ[i,k] Φκ[j,l] + Φκ[i,l] Φκ[j,k] )+
3
using Theorem 5.1 or Algorithm 2 depending on whether
1
ρ(Qi,j Sk,l +Qi,k Sj,l +Qi,l Sj,k +Qj,k Si,l +Qj,l Si,k +Qk,l Si,j ) (1) is stable, and so on. The resulting algorithm is shown in
6
Algorithm 3.
2
2
for i, j, k, l = 1, ..., n, and F s ∈ Rn ×n , the matrix unfolding
of F s , can be obtained from
VI. N UMERICAL E XAMPLES
s

s
s
Fn(i−1)+k,n(j−1)+l
= Fi,j,k,l
,

T

s

i, j, k, l = 1, ..., n.

2 This is analogous to an unsymmetric matrix whose quadratic form is inT
variant to the transpose operation, i.e. z T F 2+F z = z T F z = z T F T z ∈ R.

To verify our proposed algorithms, we use a consensus
network model defined over a 500-node graph G500 . The
graph is randomly generated with 0.5 overall probability for
edge attachment, and 6 spatial clusters with a proportion

11

(a) Clusters by coherency clustering

(b) Clusters by H2 open-loop clustering

(c) Clusters by H2 closed-loop clustering with Q1

(d) Clusters by H2 closed-loop clustering with Q2

Figure 4: Clustering of the 500-node network with r = 6 clusters. Nodes assigned to the same cluster are marked by the same color in one
figure. Note that only the node identities are comparable between figures, but not the cluster identities.
Table I: Input specifications for Algorithm 1

Algorithm 3: Iterative algorithm for finding P
Input : A, B, Bd , Q, R and r
1
2
by Definition 3.6;
1 Compute Φκ
0
0
0
2 Choose w = 1n , and compute W = diag(w ) and the

Algorithm
H2 closed-loop clustering
H2 open-loop clustering [16]
Coherency clustering [3]

1

3
4
5
6
7
8

k-means input Ψ0 = (W 0 )−1 Φκ2 ;
k = 1;
while I k−1 6= I k or within maximum iterations do
Solve I k from Algorithm 1 by (Ψk−1 , wk−1 , r) ;
Update wk from Theorem 5.1 or Algorithm 2 by I k ;
k =k+1 ;
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 ≤ Mi,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 364th node,
i.e. Bd equals to the 364th column of I500 . We assume
B = R = I500 and two choices of Q as: a scaled identity
matrix Q1 = 1000 × I500 , and Q2 = [L(G500 )]2 , which is the
square of the unweighted Laplacian matrix of G500 . Both Q1
and Q2 satisfy Assumption 2.1. For this simulation example,
solvability of (8) and stability of ĝ(s) are guaranteed by
Theorem A.2 and Theorem A.5 in Appendix A, respectively.

Data
Φ
Φo a
Φc b

Weight
w
w
1n

Number of Clusters
r
r
r
1

T

Let vcT be the complement of v̄T , and then Φo := Φo2 Φo2 =
R
T T
T
vc [ 0∞ e(vc Avc )τ vcT Bd BdT vc e(vc A vc )τ dτ ]vcT .
b Let the eigenvalues of −L(G) be 0 = λ > λ ≥ ... ≥ λ , Ψ =
n
c
1
2
[v1 , ..., vr ] where vi is the right eigenvector of λi .
a

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, H2 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

12

0.45

0.35
Closed−loop clustering
Open−loop clustering
Coherency clustering

0.4

Closed−loop clustering
Open−loop clustering
Coherency clustering

0.3

0.35
0.25
kg(s)−ĝ(s)k H2
kg(s)k H2

kg(s)−ĝ(s)k H2
kg(s)k H2

0.3
0.25
0.2

0.2
0.15

0.15
0.1
0.1
0.05

0.05
0
0

100

200

300

400

0
0

500

100

r: number of clusters

(a) Design with Q1

300

400

500

(b) Design with Q2
kg(s)−ĝ(s)kH2
kg(s)kH2

nodes), while coherent clusters are decided by both edgeweights and node-weights. The H2 open-loop clusters are
shown in Fig. 4b. As such, they do not follow any definite
pattern with respect to the spatial clusters as they are based
on node aggregation following from the H2 -norm distance of
their output responses. Figures 4c and 4d, on the other hand,
show the clusters identified by our H2 closed-loop algorithm
(Algorithm 1) for Q = Q1 and Q = Q2 , respectively. Both
of these clusters are different from each other for obvious
reasons. They are also different from the spatial clusters,
the coherent clusters as well as the H2 open-loop clusters
as Algorithm 1 is related to the closed-loop controllability
subspace.
We also illustrate the effectiveness of H2 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.
kg(s)−ĝ(s)k
We vary r from 1 to 500, and calculate the ratio kg(s)kH H2
2
resulting from the three clustering algorithms. The results are
shown in Fig. 5. For both Q1 and Q2 the closed-loop clustering
outperforms the other two methods in approaching the H2
performance of g(s). Therefore, even for very small values
of r, the projected controller achieves significantly close H2
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 L1 sparsity-promoting LQR [14]. For
 example,
for this system a standard LQR would require 500
= 124750
2
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 H2 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

with respect to the number of clusters

0.05
Sparsity-promoting LQR
Reduced-order LQR

0.04
kg(s)−ĝ(s)kH2
kg(s)kH2

Figure 5: Performance loss

200

r: number of clusters

0.03
0.02
0.01
0
506

515

528 536 3104

11661

11989

15003

21325

Number of communication links

Figure 6: Projected LQR controller from clustering vs. LQR controller from sparsity-promoting algorithm with Q1 .

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 G500 for implementing the
controller.
B. Weight Design
We next apply Algorithm 2 on G500 to find the optimal
projection weight w while fixing the clusters to those obtained
from coherency. These clusters as previously shown in Fig. 4a
closely resemble the spatial clusters, and their clustering sets
are given by I1 = {1, ..., 167}, I2 = {168, ..., 178}, I3 =
{179, ..., 344}, I4 = {345, ..., 379}, I5 = {380, ..., 463} and
I6 = {464, ..., 500}. After running the algorithm with both
Q1 and Q2 , we plot the two weight vectors compared with
v̄ = v̄(A) in Fig. 7. It can be seen that the weight vectors from
Q1 and Q2 are very different than v̄ or between themselves.
On the other hand, both of these weights at the 364th 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

13

Table II: Results of weight design
Case

Penalty factor ρ

Q1

2
0.011 v̄T Qv̄

Q2

2
0.007 v̄kΦk
T Qv̄

Relative error kg(s) − ĝ(s)kH2 /kg(s)kH2
with w = v̄
with w from Algorithm 2
29.14%
7.35%

kΦk

35.69%

21.79%

1

The convergence threshold for power iteration is chosen as δ = 0.05.
is included in the penalty factor to normalize the two objective
functions in (opt*) to the same scale.
2 kΦk2
v̄ T Qv̄

Table III. The table verifies the O(n3 ) complexity for fullorder LQR vs. the O(n2 r) complexity for reduced-order LQR.
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.
Table III: Scalability Results




w from Q1
w from Q2
v̄



n

Node 364

1000
2000
4000
6000
8000
10000

 






Computation time
Full-order
Reduced-order
16.34 sec
0.79 sec
134.91 sec
2.64 sec
20.36 min
11.61 sec
71.06 min
24.73 sec
Out of memory
48.39 sec
Out of memory
75.90 sec

Perf. loss

kg(s)−ĝ(s)kH2
kg(s)kH2

3.74%
1.06%
0.25%
0.13%
-



VII. C ONCLUSION

ï

Cluster 1
ï 




Cluster 2







Cluster 3




Cluster 4 Cluster 5 Cluster 6








Figure 7: Weight designs from Q1 and Q2

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 Q1 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
1
and κ = 5 for computing Φκ2 , and present the results in

0.25

v̄
w

0.2

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.
A PPENDIX A: S PECIAL C ASE

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 = ns . Each node has a real-valued
node-weight mi > 0, and each edge ij ∈ E has a real-valued
edge-weight aij = aji > 0. A consensus network can then be
modeled in its
simplest form by letting Aij = m−1
i aij and
−1 P
a
in
(2).
The
overall
system
writes as
Aii = − mi
ij
j∈Ni
ẋM (t) = M −1 L(G)xM (t) + BM u(t) + BdM d(t),

where L(G) is the edge-weighted graph Laplacian matrix, and
M = diag([m1 , ..., mn ]) is the matrix of node weights. For
ease of analysis, let us consider a coordinate transformation
1
x = M 2 xM , in which case the state matrices in (1) become
1

kg(s)−ĝ(s)k H2
kg(s)k H2

C ONSENSUS N ETWORK

1

1

1

A = M − 2 L(G)M − 2 , B = M 2 BM , Bd = M 2 BdM . (40)

0.15

Notice that A here is a negative-semidefinite matrix, with only
1
one zero eigenvalue at Av̄(A) = 0, v̄(A) = √ 1 M 2 1n ,

0.1

tr(M)

0.05

0

OF

2

3

4

r: number of clusters

5

Figure 8: Performance loss by weight refinement

6

which we will simply denote as v̄. The state matrix Ã from
the reduced-order system (6) satisfies the following property.
Proposition A.1: Ã is negative-semidefinite (or -definite) if
w ∈ span(v̄) (or w 6∈ span(v̄)).

14

Proof: Denote the complement of the projection matrix
P by U , such that D = [P T U T ]T is unitary, i.e. DDT = In .
Then, Ã = P AP T is a leading principal of the matrix


P AP T P AU T
T
DAD =
,
U AP T U AU T
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 Ã  0.
Moreover, from Definition 2.2, P is defined over w such that
P T P w = w. Then if w ∈ span(v̄), Ã preserves the zero
eigenvalue from A since ÃP v̄ = P Av̄ = 0. Therefore, Ã
is negative-semidefinite if w ∈ span(v̄), or negative-definite
otherwise.
1) Existence Condition: The next theorem shows that the
existence of X̃ for the reduced-order ARE (8) is guaranteed
for any w and I, and matrices Q̃ = P QP T and G̃ = P GP T .
1
T
Theorem A.2: If (Q 2 , A) is detectable and (A, G 2 ) is
stabilizable, then (8) is guaranteed with a unique stabilizing
solution X̃ = X̃ T  0.
Proof: For Q  0 and G  0, similar to Ã in the proof
of Proposition A.1, we have both Q̃  0 and G̃  0. Recall
T
that (Q̃ 2 , Ã) is detectable if and only if for all λ and v that
T
Ãv = λv and λ ≥ 0, Q̃ 2 v 6= 0. If w ∈ span(v̄), Ã  0 has
only one zero eigenvalue, with the corresponding eigenvector
P v̄. Under this situation, v̄ T P T Q̃P v̄ = v̄ T Qv̄ 6= 0 since
T
(Q 2 , A) is detectable. If w 6∈ span(v̄), Ã is negative-definite,
T
which means (Q̃ 2 , Ã) is trivially detectable. The same rational
1
applies to showing stabilizability of (Ã, G̃ 2 ). Therefore, both
stabilizability and detectability are satisfied, and thus (8)
guarantees a unique positive-semidefinite solution.
2) Bound of ARE solution: The following lemma provides
a specific value for β(A, G, Q) in Lemma 3.1.
Lemma A.3: The bound λ̄(X) ≤ β(A, G, Q) holds for
β(A, G, Q) =

σ̄(Q + G)
.
2σ(G − A)

(41)

Proof: The expression for β(A, G, Q) above can be
T
simply obtained by letting Kt = −G 2 , Dt = In and
F = 2G − 2A ≻ 0 in Lemma 3.1.
Note that the value specified by this Lemma only involves
finding two extreme singular values, which can be computed
very efficiently by Lanczos algorithm in O(n) complexity.
3) Stability condition: We next state two sufficient conditions that guarantee the stability of ĝ(s) for consensus
networks.
Theorem A.4: The system ĝ(s) is stable if w = v̄, and I
is an almost equitable partition [25] of the graph G, which
means for k 6= l, the edge weight aij is equal for all i ∈ Ik
and j ∈ Il .
Proof: A similarity transformation of D = [P T U T ]T
and DT on A − GX̂ yields


Ã − G̃X̃
P AU T
. (42)
D(A − GX̂)DT =
U AP T − U GP T X̃ U AU T

From [25], if G admits an almost equitable partition I, the
corresponding matrix P with w = v̄ will satisfy AT P T =
P T N for some N ∈ Rr×r . As a result, P AU T = 0, and thus

ĝ(s) is stable since Ã − G̃X̃ and U AU T are both Hurwitz
given w = v̄.
Theorem A.5: Assume that B is a square invertible matrix.
Then, ĝ(s) is stable if G is similar to αIn for some α > 0.
Proof: Denote the right eigenspace of G by V . If G is
similar to αIn , it follows that V T (GX̂ + X̂G)V = 2αV T X̂V .
This means that GX̂ + X̂G is positive-semidefinite. From the
matrix majorization property [26], we also have
2λ̄(A − GX̂) ≤ z T (A + AT − GX̂ − X̂G)z
hold for any non-zero vector z, with the RHS being nonpositive given A + AT  0 and −GX̂ − X̂G  0. Hence
for λ̄(A − GX̂) to be strictly negative, z T Az = 0 should
contradict z T GX̂z = 0. To prove the contradiction, we assume
z T Az = 0, namely z ∈ span(v̄). Note that z T GX̂z = 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 wITi v̄Ii ,
i = 1, ..., r, which means P v̄ 6= 0 and thus z 6∈ ker(P ).
1
Given that (A, BR− 2 ) is stabilizable, we have v̄ T Gv̄ 6= 0,
i.e. z 6∈ ker(G). To this end, the last condition remaining is
z ∈ span(v̄) 6⇔ P z ∈ ker(X̃) to complete the contradiction.
By assuming a non-zero vector v ∈ ker(X̃), pre- and postmultiplying (8) with v yields
v T ÃT X̃v + v T X̃ Ãv + v T Q̃v − v T X̃ G̃X̃v = 0.
It can be easily verified that ker(X̃) is an Ã-invariant subspace
contained in the null-space of Q̃. Given z ∈ span(v̄), from
the proof of Proposition A.1 we know that P z is Ã-invariant,
i.e., P z is an eigenvector of Ã, when w ∈ span(v̄). As a
T
result, z T P T Q̃P z = z T Qz 6= 0 since (Q 2 , A) is detectable.
This verifies that P z 6∈ ker(Q̃), which proves P z 6∈ ker(X̃).
Therefore, we conclude that z T Az = 0 contradicts z T (GX̂ +
X̂G)z = 0, and thus ĝ(s) is stable.
A PPENDIX B: P ROOFS
A. Proof of Theorem 3.3
Consider a Lyapunov function V (x) = xT Xx > 0, where
X ≻ 0 is the solution of ARE (4). For ĝ(s) to be asymptotically stable, V̇ (x) needs to be negative, or equivalently
(A − GX̂)T X + X(A − GX̂) ≺ 0.
Using the ARE (4), the LMI above reduces to
Q ≻ XGE + EGX − XGX,
which holds if λ(Q) > λ̄(XGE + EGX − XGX). Notice
that the RHS of the eigenvalue inequality follows
λ̄(XGE +EGX −XGX) ≤ λ̄(XGE +EGX)+ λ̄(−XGX),
where we have respectively λ̄(XGE + EGX)
≤
σ̄(XGE + EGX) ≤ 2σ̄(XGE) ≤ 2σ̄(X)σ̄(G)σ̄(E),
and −λ̄(−XGX) = σ(XGX) ≥ σ(X)2 σ(G). Incorporating
these two upper bounds yields the condition in (17).

15

B. Proof of Lemma 3.4
From Lemma 3.1, the ARE solution X is bounded by
σ̄(X) = λ̄(X) ≤ β(A, G, Q). Note that (17) is satisfied if
σ(Q) > 2σ̄(X)σ̄(G)σ̄(E),

(43)

The eigenvalues of B ⊗ A are counted by λi (A)λj (B) with
i, j = 1, ...n, and according to the Weyl’s inequality of
eigenvalues [26], we have the lower bound for σ 2 (L) as
σ 2 (L) ≥ λ(AAT ) + λ(B T B) + 2λ(B ⊗ A)

= σ 2 (A) + σ 2 (B) + 2λ̄(A)λ̄(B) ≥ σ 2 (B).

where the RHS is further bounded by
2σ̄(X)σ̄(G)σ̄(E) ≤ 2β(A, G, Q)σ̄(G)[β(A, G, Q)+σ̄(X̃)]. (44)

Therefore, (43), and then (17) will hold if σ(Q) is greater than
the RHS of (44), which yields (18) in Lemma 3.4.

Combining (52) with (53) then yields the following bound

where ǫ1 =
T

The bound in (23) is derived assuming the worst case from
(14), where Q̃=P QP T +αIr and G̃=P QP T +αIr for α > 0.
We divide the proof into three steps.
1) : We derive an analytical expression for E by recovering
the reduced-order ARE (8) to the full dimension as
P T (ÃT X̃ + X̃ Ã + Q̃ − X̃ G̃X̃)P = 0.

(45)

AT X̂ + X̂A + Q − X̂GX̂ = R,

(47)

with the residue of the approximate ARE denoted by
T

T

T

R := αX̂ + U U A X̂ + X̂AU U + Q − P Q̃P. (48)
By subtracting (4) from (47), we get the Sylvester equation
(A − GX̂)T E + E(A − GX) = −R.

(49)

From (49), we are able to explicitly write E as a function of
A − GX, A − GX̂ and R, and hence obtain an initial bound
1
for kEΦ 2 kF in the next step.
T
1
2) : Pre- and post-multiplying (49) with Φ 2 and Φ 2
respectively, the Sylvester equation takes the form
T

T

1

1

A1 Φ 2 EΦ 2 + Φ 2 EΦ 2 B1 = −Φ 2 RΦ 2 ,
T

(50)
T

where we use the notations A1 = Φ 2 (A − GX̂)T Φ− 2 and
1
1
B1 = Φ− 2 (A−GX)Φ 2 for brevity. It can be easily shown that
λ(A1 ) < 0 and λ(B1 ) < 0, which implies λi (A) + λj (B) 6= 0
for any i, j = 1, ..., n so that (50) is solvable. Therefore, the
1
weighted error EΦ 2 can be expressed as
T

1

T

1

EΦ 2 = Φ− 2 · unvec[L−1 · vec(−Φ 2 RΦ 2 )],

(51)

where L = In ⊗ A + B T ⊗ In is an n2 × n2 matrix. Since
the Frobenius norm is unitary invariant, taking norm on both
1
sides of (51) provides an upper bound on kEΦ 2 kF as
1

T

1

1

]

(54)

> 0 is independent of P . In
T

1

T

1

T

1

T

1

kEΦ 2 kF ≤ σ̄(L−1 )σ̄(Φ− 2 )kΦ 2 RΦ 2 kF .

Note that σ̄(L−1 ) follows σ̄(L−1 ) =
calculated by σ 2 (L) = λ(LLT ) with

1
σ(L) ,

(52)

where σ(L) is

LLT = In ⊗ AAT + B T B ⊗ In + B T ⊗ AT + B ⊗ A.

1

T

1

+ Φ 2 (QU T U + U T U Q − U T U QU T U − αP T P )Φ 2 , (55)

contains the inexplicit functional X̂. We then bypass this term
in the final step.
3) : Taking norm on both sides of (55), and then isolating
the norm of X̂, we can form the bound
1

1

1

+ σ̄(Q)ξ 2 + ασ̄(Φ)[σ̄ 2 (X̂) + 1],

(46)

where U is the complement of P . Thereby substituting ÃP
and P T ÃT in terms of (46), and after a few calculations, (45)
yields the approximated ARE (for details, please see [30])

1

1
σ̄(Φ− 2 )
1
−1
σ[Φ 2 (A−GX)Φ 2

Φ 2 RΦ 2 = αΦ 2 X̂ 2 Φ 2 + Φ 2 U T U AT X̂Φ 2 + Φ 2 X̂AU T U Φ 2

T

ÃP = P A − P AU T U,

T

1

kΦ 2 RΦ 2 kF ≤2[σ̄(A)σ̄(Φ 2 )σ̄(X̂) + σ̄(QΦ 2 )]ξ

Notice that A and Ã are related by

T

T

1

kEΦ 2 kF ≤ ǫ1 kΦ 2 RΦ 2 kF ,

(54), the norm of the weighted residue Φ 2 RΦ 2 , written by

C. Proof of Theorem 3.5

2

(53)

(56)

1
2

with ξ = kU T U Φ kF . Recall that σ̄(X̂) = σ̄(P T X̃P ) =
σ̄(X̃), where X̃ is the solution of the reduced-order ARE
(8). The norm σ̄(X̃) can be further bounded by β(Ã, G̃, Q̃)
through Lemma 3.1. Theorem 3.5, therefore, follows from
(54), (56) and σ̄(X̃) ≤ β(Ã, G̃, Q̃) ≤ supP β(Ã, G̃, Q̃).
D. Proof of Lemma 3.7
To prove the error bound (28), we define a matrix Φ̄ as


Φ̄ = Φ̄1 · · · Φ̄nb ,
(57)
1

where Φ̄i = Y diag(Y −1 bi )C 2 , i = 1, ..., nb and bi is the ith
1
T
column of Bd . The matrix Φ̄ satisfies Φ = Φ̄Φ̄T = Φ 2 Φ 2 .
Besides Y1 and Ω1 defined in Definition 3.6, we further denote
−1
Y2 = Y:,κ+1,n and Ω2 = Yκ+1:n,:
, and partition the Cholesky
1
Decomposition C 2 as
# " 1
#
" 1
2
2
1
0
0
C1:κ,1:κ
C1,1
2
=
.
C =
1
1
1
1
2
2
2
2
C2,1
Cκ+1:n,1:κ
C2,2
Cκ+1:n,κ+1:n
With these notations, Φ̄i in (57)
h can be decomposed
i into
1
2
Φ̄i = Φ̄i,s + Φ̄i,f , where Φ̄i,s = Y1 diag(Ω1 bi )C1,1
0 and
h
i
1
1
2
2
Φ̄i,f = Y2 diag(Ω2 bi )C2,1
, and thus Φκ
Y2 diag(Ω2 bi )C2,2
Pnb
can be rewritten as Φκ = i=1
Φ̄i,s Φ̄Ti,s . Notice that ξ =
1
k(In − P T P )Φ 2 kF = k(In − P T P )Φ̄kF satisfies
1

ξκ∗ ≤ ξ ≤ k(In − P T P )Φκ2 kF + k(In − P T P )Φ̄f kF , (58)


where Φ̄f = Φ̄1,f · · · Φ̄nb ,f . The second norm on the
RHS of (58) is further bounded by k(In − P T P )Φ̄f kF ≤
kΦ̄f kF with
v
v
u nb
u
n
X
u X
1
1
u
1
2
2
2
2
t
2
− −.
kΦ̄f kF ≤ η
(kC2,1 kF + kC2,2 kF ) = tη 2 nb
2λ
i
i=κ+1
i=1

Inserting this along with P = argminP ξκ to the RHS of (58)
yields the error bound (28).

16

E. Proof of Lemma 5.3
Denote ni = |Ii |c , i = 1, ..., r, the objective function fV̄
in (opt2) can be expanded as
fV̄ = tr(V̄ T P T P QP T P V̄ ) = (w ∗ w)T (Q ∗ V̄ V̄ T )(w ∗ w)
where ∗ is the Khatri-Rao product defined by
QI1 ,I1 ⊗ V̄I1 ,: V̄IT1 ,:

..
Q ∗ V̄ V̄ T = 
.
QIr ,I1 ⊗ V̄Ir ,: V̄IT1 ,:



w ∗ w = wIT1 ⊗ wIT1

···
..
.
···

···


QI1 ,Ir ⊗ V̄I1 ,: V̄ITr ,:

..
,
.
QIr ,Ir ⊗ V̄Ir ,: V̄ITr ,:

wITr ⊗ wITr

T

.

Denote the block-diagonal submatrix of Q by Qd , i.e. Qd =
diag(QI1 ,I1 , · · · , QIr ,Ir ), and the off-diagonal by Qo = Q −
Qd . Therefore, we can find the difference between objective
functions of (opt3) and (opt*) as (w∗w)T (ρQo ∗ V̄ V̄ T )(w∗w).
According to the Gershgorin circle theorem, the eigenvalues
of ρQo ∗ V̄ V̄ T are all bounded inside the range of (−kρQo ∗
V̄ V̄ T k1 , kρQo ∗V̄ V̄ T k1 ). Therefore, the theorem follows when
Je = kρQo ∗ V̄ V̄ T k1 .
R EFERENCES
[1] D. D. Siljak, Decentralized Control of Complex Systems, Courier Corporation, 2011.
[2] P. Kokotovic, H. K. Khalil, and J. O’reilly, Singular Perturbation
Methods in Control: Analysis and Design, SIAM, 1999.
[3] J. Chow and P. Kokotovic, “Time Scale Modeling of Sparse Dynamic
Networks,” IEEE Transactions on Automatic Control, vol. 30(8), pp.
714-722, 1985.
[4] B. Moore, “Principal Component Analysis in Linear Systems: Controllability, Observability, and Model Reduction,” IEEE Transactions on
Automatic Control, vol. 26(1), pp. 17-32, 1981.
[5] G. Obinata and B. D. Anderson, Model Reduction for Control System
Design, Springer Science & Business Media, 2012.
[6] A. C. Antoulas, Approximation of Large-Scale Dynamical Systems,
SIAM, 2005.
[7] K. Zhou and J. C. Doyle, Essentials of Robust Control, vol. 180, Prentice
Hall Upper Saddle River, NJ, 1998.
[8] A. M. Boker, T. R. Nudell, and A. Chakrabortty, “On Aggregate Control
of Clustered Consensus Networks,” 2015 American Control Conference
(ACC), pp. 5527-5532, Chicago, IL, USA, 2015.
[9] S. Hara, J. I. Imura, K. Tsumura, T. Ishizaki, and T. Sadamoto, “Glocal
(Global/Local) Control Synthesis for Hierarchical Networked Systems,”
2015 IEEE Conference on Control Applications (CCA), pp. 107-112,
Sydney, NSW, 2015.
[10] T. Sadamoto, T. Ishizaki, and J. I. Imura, “Hierarchical Distributed
Control for Networked Linear Systems,” 53rd IEEE Conference on
Decision and Control, pp. 2447-2452, 2014.
[11] D. Madjidian and L. Mirkin, “Distributed Control with Low-Rank
Coordination,” IEEE Transactions on Control of Network Systems, vol.
1(1), pp. 53-63, March 2014.
[12] S. Lloyd, “Least Squares Quantization in PCM,” IEEE Transactions on
Information Theory, vol. 28(2), pp. 129-137, Mar 1982.
[13] L. Qi, F. Wang, and Y. Wang, “Z-Eigenvalue Methods for a Global
Polynomial Optimization Problem,” Mathematical Programming, vol.
118, pp. 301-316, 2009.
[14] F. Lin, M. Fardad, and M. R. Jovanovic, “Design of Optimal Sparse
Feedback Gains via the Alternating Direction Method of Multipliers,”
IEEE Transactions on Automatic Control, vol. 58(9), pp. 2426-2431,
Sep 2013.
[15] N. Monshizadeh, H. L. Trentelman, and M. K. Camlibel, “ProjectionBased Model Reduction of Multi-Agent Systems using Graph Partitions,” IEEE Transactions on Control of Network Systems, vol. 1(2),
pp. 145-154, Jun 2014.
[16] T. Ishizaki, K. Kashima, A. Girard, J. Imura, L. Chen, and K. Aihara,
“Clustered Model Reduction of Positive Directed Networks,” Automatica, vol. 59, pp. 238-247, 2015.

[17] P. Benner and Z. Bujanovi, “On the Solution of Large-Scale Algebraic
Riccati Equations by Using Low-Dimensional Invariant Subspaces,”
Linear Algebra and its Applications, vol. 488, pp. 430-459, 2016.
[18] M. Rotkowitz and S. Lall, “A Characterization of Convex Problems in
Decentralized Control,” IEEE Transactions on Automatic Control, vol.
51(2), pp. 274-286, Feb. 2006.
[19] X. Wu and M. R. Jovanovic, “Augmented Lagrangian Approach to
Design of Structured Optimal State Feedback Gains,” IEEE Transactions
on Automatic Control, vol. 56, pp. 2923-2929, Dec 2011.
[20] S. Fattahi, G. Fazelnia and J. Lavaei, “Transformation of Optimal
Centralized Controllers into Near-Global Static Distributed Controllers,”
53rd IEEE Conference on Decision and Control, pp. 4915-4922, 2015.
[21] Y. S. Wang, N. Matni, and J. C. Doyle, “Localized LQR Optimal
Control,” 53rd IEEE Conference on Decision and Control, pp. 16611668, 2014.
[22] Y. S. Wang and N. Matni, “Localized LQG Optimal Control for LargeScale Systems,” 2016 American Control Conference (ACC), pp. 19541961, Boston, MA, USA, 2016.
[23] N. Xue and A. Chakrabortty, “H2 -Clustering of Closed-Loop Consensus
Networks under Generalized LQR Designs,” 55th IEEE Conference on
Decision and Control, pp. 5116-5121, Las Vegas, NV, Dec 2016.
[24] A. Helmy, “Architectural Framework for Large-Scale Multicast in Mobile Ad Hoc Networks,” IEEE International Conference on Communications, pp. 2036-2042, vol.4, 2002.
[25] D. M. Cardoso, C. Delorme, and P. Rama, “Laplacian Eigenvectors
and Eigenvalues and Almost Equitable Partitions,” European Journal
of Combinatorics, vol. 28, pp. 665-673, 2007.
[26] J. N. Franklin, Matrix Theory, Courier Corporation, 2012.
[27] R. Davies, P. Shi, and R. Wiltshire, “New Lower Solution Bounds of
the Continuous Algebraic Riccati Matrix Equation,” Linear Algebra and
Its Applications, vol. 427(2-3), pp. 242-255, 2007.
[28] G. H. Golub and C. F. Van Loan, Matrix Computations, JHU Press, vol.
3, 2012.
[29] M. S. Traiantafyllou and F. S. Hover, Maneuvering and Control of
Marine Vehicles, Cambridge, MA, 2002.
[30] J. Sun, “Perturbation Theory for Algebraic Riccati Equations,” SIAM
Journal on Matrix Analysis and Applications, vol. 19, pp. 39-65, 1998.
[31] E. Kofidis and P. A. Regalia, “On the Best Rank-1 Approximation
of Higher-Order Supersymmetric Tensors,” SIAM Journal on Matrix
Analysis and Applications, vol. 23, pp. 863-884, 2002.

Nan Xue (S’15) received his B.E. degree in Electrical Engineering from Xi’an
Jiaotong University, China in 2013. He is currently pursuing his PhD degree
in Electrical Engineering at North Carolina State University, Raleigh, NC. His
research interests include analysis, control and model reduction of large-scale
networked dynamic systems and power systems.

Aranya Chakrabortty (M’08, SM’15) received his PhD degree in Electrical
Engineering from Rensselaer Polytechnic Institute, Troy, NY in 2008. From
2008 to 2009 he was a postdoctoral research associate at the Aeronautics and
Astronautics department of University of Washington, Seattle. He is currently
an Associate Professor in the Electrical and Computer Engineering department
of North Carolina State University, Raleigh, NC, where he is also affiliated
to the FREEDM Systems Center. His research interests are in all branches
of control theory with applications to power systems, especially in wide-area
monitoring and control of large power systems using Synchrophasors. He
received the NSF CAREER award in 2011.

