Hierarchical H2 Control of Large-Scale Network Dynamic Systems

arXiv:1703.07393v2 [cs.SY] 26 Sep 2017

Nan Xue and Aranya Chakrabortty

Abstract— Standard H2 optimal control of networked dynamic systems tend to become unscalable with network size.
Structural constraints can be imposed on the design to counteract this problem albeit at the risk of making the solution nonconvex. In this paper, we present a special class of structural
constraints such that the H2 design satisfies a quadratic
invariance condition, and therefore can be reformulated as a
convex problem. This special class consists of structured and
weighted projections of the input and output spaces. The choice
of these projections can be optimized to match the closedloop performance of the reformulated controller with that
of the standard H2 controller. The advantage is that unlike
the latter, the reformulated controller results in a hierarchical
implementation which requires significantly lesser number of
communication links, while also admitting model and controller
reduction that helps the design to scale computationally. We
illustrate our design with simulations of a 500-node network.
Index Terms— Optimal Control, Quadratic Invariance, Network Systems, Controller Reduction.

I. I NTRODUCTION
Driven by recent advances in cyber-physical systems,
control synthesis for large-scale network dynamic systems
has become an increasingly prevalent topic in the control and
network communities [1]. Physical networks such as power
systems, wireless sensor networks, or the recently emerging
internet-of-things consist of numerous heterogeneous subsystems that may be defined across complex topologies and
wide geographical spans. The typical numbers of subsystems
in these networks can scale from thousands to millions,
making the design of tractable control mechanisms very
challenging. To tackle the curse of dimension, traditional
model reduction-based techniques such as singular perturbations [2], balanced truncation [3], and controller reduction
[4] have been developed decades ago to facilitate the analysis
and control of large-scale systems, in general. However,
these techniques are not readily applicable to networks as
their control schemes are mostly unstructured, and hence
agnostic of the network topology that contains important
constraints for both design and implementation. One such
pertinent constraint is communication. For example, conventional feedback controllers such as H2 optimal controllers are
defined over unstructured dense transfer matrices, and thus
implementing this control would necessitate an impractically
large number of communication links for even a moderately
sized network.
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.

Starting from the idea of decentralized control [5], several
seminal papers have developed design tools to impose desired implementation structure on controllers. One of them,
the work in [6], states that a decentralized optimal control
admits a convex reformulation if the structural constraint
on the controller is quadratic invariant (QI) under the plant
model. Built around quadratic invariance, papers including
[6], [7] and [8] incorporate sparsity structures in the controller, aiming to reduce the communication density and to
cope with delays. The restriction of these designs is that the
sparsity pattern meeting the QI condition is highly dependent
on the sparsity of the open-loop network, and therefore,
the choice of sparsity can be limited. Moreover, the design
complexities of these controllers are dictated by the order of
the open-loop plant, and thus can become unscalable for very
large-scale networks. A suboptimal design is proposed in [9]
to generalize the choice of sparsity structure by relaxation
algorithms using l1 -weighted norm. The design, however, is
still computationally demanding.
Motivated by these challenges, in this paper we present
a special class of structural constraint for designing H2
controllers. The constraint set is defined by two projections over the input and output spaces, which result in a
hierarchical implementation of the controller. The implementation mechanism works as follows. Selected network
subsystems send their output measurements to a designated
set of coordinators. These coordinators take the average of
these measurements, and exchange the information between
themselves to generate a H2 control law. Each coordinator,
thereafter, broadcasts this control signal to its respective
set of subsystems. The overall execution, therefore, is completely hierarchical, and requires lesser number of communication links. Preliminary work incorporating this hierarchy
has been proposed in our recent paper [10] for a LQG
controller built upon the idea of clustering based projection.
However, the design there is restricted by the choice of
certain design parameters that guarantee closed-loop stability.
In this paper we formalize the design in [10] into a general
H2 framework, where such restrictions do not exist. The
hierarchical structure is proven to meet the QI condition
irrespective of the structure of the plant, as a result of which,
the control problem is reformulated as a standard unconstrained H2 control. The reformulated problem preserves the
hierarchical structure inside its input and output matrices, and
thus allows for model and controller reduction techniques to
reduce the design complexity without breaking the structure
of the controller. The reduction technique adopted in this
paper is Hamiltonian-based approximation for solving the
underlying algebraic Riccati equations (ARE), as proposed

in [11]. The recent paper [12] also addresses similar goals
as ours using receding-horizon control, but the scalability
of their controller is subject to the sparsity structure of the
open-loop network. Related hierarchical designs have been
proposed in [13], [14], however, they do not exploit the
convex reformulation provided by quadratic invariance.
The remainder of this paper is organized as follows.
Section II defines the hierarchical constraint, and formulates
the problem of hierarchical H2 control. The same problem
is parameterized and reformulated into an unconstrained H2
control via quadratic invariance in Section III. Section IV
presents a Hamiltonian-based reduction technique to simplify
the design of the reformulated controller. The optimality gap
between the hierarchical H2 and unconstrained H2 problems
is discussed in Section V, where we also propose a design on
clustering sets to tighten the gap. Validations of our proposed
designs are illustrated via simulations of large-scale networks
in Section VI. Section VII concludes the paper and presents
some future works.
Notation The following notations will be used throughout
this paper: |I|c : cardinality of a set I, diag(m): diagonal
matrix with vector m on its principal diagonal, diag(M, N ):
block-diagonal matrix with matrices M and N on its principal diagonal, tr(M ): trace operation on a matrix
pM , kM kF :
Frobenius norm of a matrix M , i.e. kM kF = tr(M M T ).
A transfer function matrix (TFM) is defined as g(s) =
−1
C(sI
− A)
 B + D, with a realization form of g(s) =

A B
. Rsp and Rp denote respectively the set of realC D
rational strictly proper TFMs and the set of real-rational
proper TFMs. RH∞ denotes the set of all stable TFMs from
Rp , and RH2 denotes the set of all stable TFMs from Rsp .
The H∞ norm of g(s) ∈ RH∞ is defined by kg(s)kH∞ =
supω σ̄[g(jω)], where σ̄ denotes the largest singular value.
The
RH2 is defined by kg(s)kH2 =
qR H2 norm of g(s) ∈ q
R∞
∞
1
∗
∗
2π −∞ tr[g (jω)g(jω)]dω.
−∞ tr[g (t)g(t)]dt =
Furthermore, a Lyapunov equation is a matrix equation in
the form of
AΦ + ΦAT + BB T = 0,
and an algebraic Riccati equation (ARE) is written as
AT X + XA + C T C + XBD−1 B T X = 0.
We use Φ = LYAP(A, B) and X = ARE(A, B, C, D) to
respectively denote solutions from these two equations.
II. P ROBLEM F ORMULATION
We motivate our design from the standard formulation of
H2 optimal control. Consider a transfer function matrix G(s)
in a realization form




B2
A B1
G11 G12
=  C1
G(s) =
0
D12  , (1)
G21 G22
0
C2 D21

where G22 is the plant model with A ∈ Rn×n , B2 ∈ Rn×nu
and C2 ∈ Rny ×n . The plant G22 is defined over ns ≤ n

interconnected subsystems with each modeled by
ẋi (t) = Aii xi (t) +

ns
X

Aij xj (t) + B2,ii ui (t)

i6=j

yi (t) = C2,ii xi (t),

,

(2)

i = 1, ..., ns

where Aii , B2,ii and C2,ii are submatrices with compatible
dimensions from A, B2 = diag(B2,11 , ..., B2,ns ns ) and
C2 = diag(C2,11 , ..., C2,ns ns ). Note that for ny , nu ≤ n,
B2,ii and C2,ii can be null matrices, that is, a subsystem can
have no input or output. Consider a standard H2 optimal
control for G(s), we begin with the following assumptions.
A1. (A, B2 ) is stabilizable and (C2 , A) is detectable
T
T
D12 ≻ 0
A2. D21 D21
≻ 0 and D12
A3. (A, B1 ) and (C1 , A) have no uncontrollable or unobservable modes on the imaginary axis
T
T
=0
A4. D12
C1 = 0 and B1 D21
Assumptions (A1-A3) are standard assumptions to formulate
H2 optimal control [15]. Assumption (A4) is made to simplify the notations for presenting the technical results. Our
proposed design is still valid if this assumption is removed.
n ×n
Let K(s) ∈ Rp u y denote the controller of interest. The
lower linear fractional transformation (LFT) of G(s) and
K(s) is defined by
f (G, K) := G11 + G12 K(I − G22 K)−1 G21 .

(3)

With these notations, the constrained H2 optimal control
problem considered in this paper can be posed as
minimize

kf (G, K)kH2

subject to K stabilizes G,
K∈S
n ×n

(4)

where S ⊆ Rp u y represents the subset of all strictly
proper transfer function matrices that satisfy certain strucn ×n
tural constraints. Note that when S = Rp u y , (4) becomes
the unconstrained H2 design. In practice, the unconstrained
case usually yields a controller defined over a dense transfer
function matrix. As a result, implementation of this controller for any large-scale network dynamic systems where
n, nu and ny may scale up to millions, would necessitate
an impractically large number of communication links. To
bypass this challenge, in this paper we propose a hierarchical
structure for the constraint S, that can significantly simplify
the implementation of K(s). Another advantageous property
of this structure is that it allows (4) to be reformulated as a
convex optimization problem as will be unfolded in Section
III. The basic mechanism behind the hierarchical formulation
is the partitioning of subsystems into a set of non-overlapping
groups or clusters. This grouping strategy can be arbitrary,
and is not necessarily dictated by any system property. The
constraint S for this hierarchy can then be formulated as
follows.
Definition 2.1: Given the index set of subsystems as V =
{1, ..., ns }, and an integer r, where 0 < r ≤ ns , define
r non-empty, distinct, and non-overlapping subsets of V
respectively denoted as I = {I1 , ..., Ir }, such that I1 ∪

... ∪ Ir = V. We denote the collection of subsystems in Ii ,
i = 1, ..., r as a cluster.
We next define two structured projection matrices Pu and Py
for clusters I1 , ..., Ir as follows.
Definition 2.2: Given a non-zero vector wu ∈ Rnu and
the input index set Vu ={1, ..., nu }, define r non-overlapping
and non-empty subsets denoted as I u = {I1u , ..., Iru }, such
that Iiu includes indices of all the inputs in cluster Ii for
i = 1, ..., r. The input projection matrix Pu ∈ Rr×nu is
defined by
( wu,[j]
, j ∈ Iiu
u
,
(5)
Pu,[ij] := kwu,[Ii ] k2
0, otherwise
where the vector wu,[Iiu ] is non-zero. The output projection
matrix Py is defined in the same way by vector wy ∈
Rny , output index set Vy ={1, ..., ny } and subsets I y =
{I1y , ..., Iry }.
Based on the projections Pu and Py , the constrained
subspace S of all hierarchical controllers can be formally
stated as follows.
n ×n
Definition 2.3: The subspace S ∈ Rp u y admits a hierarchical structure defined over projection matrices Pu ∈
Rr×nu and Py ∈ Rr×ny if there exists a lower-dimensional
transfer matrix S̃ ∈ Rr×r
such that
p
S = PuT S̃Py .
(6)
The controllers in S contains the projection structures of
Py and Pu , which lead to a sequential two-layer hierarchical
control architecture. We illustrate this architecture along with
the implementation of controllers in S by the following
example.
Example 2.4: Consider a networked system shown in the
bottom layer of Fig. 1, where four subsystems are indexed by
V = {1, 2, 3, 4}, and are connected by a simple line graph.
All four subsystems are assumed to have scalar outputs,
and only the first three are equipped with inputs, i.e. Vy =
{1, 2, 3, 4} and Vu = {1, 2, 3}. The system is partitioned
into two clusters I1 = {1, 2} and I2 = {3, 4}. Let both wu
and wy to be vectors of all ones. Then given I1u = {1, 2},
I2u = {3}, I1y = {1, 2} and I2y = {3, 4}, the projection
matrices can be constructed as
#
"

 1
√1
√1
√
√1
0
0
0
2
2
2
2
.
, Py =
Pu =
√1
√1
0
0
0
0 1
2
2


s s2
, the hierarchical subspace S follows
Denoting S̃ = 1
s3 s4
as
 s1 s1 s2 
2

 s1

S = PuT S̃Py =  s23
2
s3
2

2
s1
2
s3
2
s3
2

√

2
s2 
√
2
s4  .
√

2
s4
√
2

Therefore, to implement the controller in S, we designate
two coordinators, one for each cluster. These are denoted as
O1 and O2 in Fig. 1. The implementation of S can then be
shown by the following three steps with illustrations in Fig.
1.

•

Step 1 - Output averaging (ȳ = Py y): Each coordinator receives output measurements from its designated cluster, and computes the averaged output ȳi =
P
wy,[j] yj
for i = 1, ..., r.
j∈I y kw
y k2
i

y,[I ]
i

Step 2 - Lower-dimensional control (ū = S̃ ȳ): Next,
the coordinators communicate with each other, and
exchange the averaged outputs ȳi , i = 1, ..., r. Each
coordinator then computes a lower-dimensional H2
control input ūi = S̃i,: ȳ, i = 1, ..., r in a distributed
manner.
T
• Step 3 - Control inversion (u = Pu ū): In the final
step, each coordinator broadcasts the control signal uIiu
to subsystems in its cluster by a simple scaling of ūi ,
i = 1, ..., r.
By applying such a hierarchical structure, the total number of communication links required for implementing the
r
control would be
 at most ns + 2 . This can be much
ns
cheaper than 2 required for the all-to-all communication
in unstructured controllers. This structure also fits many
common cyber-physical networks such as power systems
where the subsystems may represent synchronous generators
and loads, the clusters may represent utility companies, and
coordinators may represent the control centers of each respective company. The proposed hierarchical implementation
can also preserve data privacy between coordinators as only
an averaged output ȳ is shared between them. Therefore,
no coordinator can infer the output measurements from
subsystems assigned to other coordinators.
Compared to [10], where the LQG design is facilitated by
only a single output projection matrix, the hierarchy in S
is defined over both input and output projections. We next
present the main results of this paper, starting with the convex
reformulation for the constrained design (4).
•

III. C ONVEX R EFORMULATION
I NVARIANCE

BY

Q UADRATIC

In this section, we present the convex reformulation for
the hierarchical H2 problem (4) facilitated by the quadratic
invariant property of S. The following well-known results of
Youla parameterization will be used to unfold the solution.
Theorem 3.1: [15] Let F and L be such that A+LC2 and
A + B2 F are Hurwitz. Then all controllers that internally
stabilize G(s) can be parameterized by
K(s) = f (Knom , Q),
with any Q ∈ RH∞ and Knom defined as


A + B2 F + LC2 −L B2
F
0
I .
Knom = 
−C2
I
0

(7)

(8)

Furthermore, the closed-loop transfer function equals to
f (G, K) = T11 + T12 QT21 ,

(9)

Ȫଵ

1

‫ݕ‬തଵ ൌ

‫ݕ‬ଵ

Ȫଶ

‫ݕ‬ଵ ൅ ‫ݕ‬ଶ

2

ʹ

‫ݕ‬ଶ

3

‫ݕ‬തଶ ൌ

‫ݕ‬ଷ

Ȫଵ

‫ݕ‬ଷ ൅ ‫ݕ‬ସ

4

‫ݕ‬തଶ

‫ݑ‬തଵ ൌ ‫ݏ‬ଵ ‫ݕ‬തଵ ൅ ‫ݏ‬ଶ ‫ݕ‬തଶ

ʹ

‫ݕ‬ସ

2

1

(a) Step 1 - Output averaging

‫ݕ‬തଵ

Ȫଵ

Ȫଶ

‫ݑ‬ത ଶ ൌ ‫ݏ‬ଷ ‫ݕ‬തଵ ൅ ‫ݏ‬ସ ‫ݕ‬തଶ

3

‫ݑ‬ଵ ൌ ‫ݑ‬ଶ ൌ

4

‫ݑ‬ଵ

1

(b) Step 2 - Lower-dimensional control

2

Ȫଶ

‫ݑ‬തଵ

‫ݑ‬ଶ

ʹ

‫ݑ‬ଷ ൌ ‫ݑ‬ത ଶ

3

‫ݑ‬ଷ

4

(c) Step 3 - Control inversion

Fig. 1: Hierarchical architecture for implementing any controller in S.

X T  0 and Y = Y T  0 are unique solutions of AREs

where T is given by

T
T = 11
T21




Â
T12
=  Ĉ1
T22
Ĉ2

A + B2 F

0

 C1 + D12 F
0

−B2 F
A + LC2
−D12 F
C2

B̂1
0
D21


B̂2
D12 
0

B1
B1 + LD21
0
D21

X = ARE(A, B2 PuT , C1 , R1 ),
T

Y = ARE(A

B2
0 
 . (10)
D12 
0

Given J, T11 , T12 and T21 from Theorem 3.1, one can
rewrite the original H2 problem (4) into a model matching
problem with respect to the Youla parameter Q as
minimize
subject to

kT11 + T12 QT21 kH2
Q ∈ RH∞

.

(11)

K(s) = f (Knom , Q) ∈ S
Note that K(s) = f (Knom , Q) ∈ S is not an affine
constraint in Q, in general. However, it has been shown in
[6] that K ∈ S is equivalent to the affine constraint Q ∈ S in
Youla domain if the subspace S is quadratic invariant under
the plant model G22 . We introduce the notion of quadratic
invariance according to [6] as follows.
Definition 3.2: The subspace S is called quadratic invariant under G22 if KG22 K ∈ S for all K ∈ S.
From this definition, it can be easily verified that for any
controller K = PuT K̃Py ∈ S, quadratic invariance holds for
our hierarchical constraint S under G22 given that
KG22 K = PuT (K̃Py G22 PuT K̃)Py ∈ S
for K̃Py G22 PuT K̃ ∈ Rr×r
sp . As K ∈ S is equivalent to Q ∈
S, (11) becomes
minimize
subject to

kT11 + T12 QT21 kH2
Q ∈ RH∞ , Q ∈ S

.

(12)

To this end, we can find the optimal solution for (4) as
follows.
Theorem 3.3: If (A, B2 PuT ) is stabilizable and (Py C2 , A)
is detectable, the hierarchical H2 problem (4) admits an
optimal solution Kopt = PuT K̃opt Py ,


A + B2 PuT F2 + L2 Py C2 −L2
,
(13)
K̃opt =
F2
0
T
T
where R1 = Pu D12
D12 PuT , R2 = Py D21 D21
PyT , F2 =
−1
T
T T −1
−R1 Pu B2 X and L2 = −Y C2 Py R2 . Matrices X =

1

(14)

, C2T PyT , B1T , R2 ).

(15)

Proof: From (6), the constraint Q ∈ S implies that
there exists Q̃ ∈ Rr×r
such that Q = PuT Q̃Py . By replacing
p
T
Q with Pu Q̃Py , (12) becomes
minimize

kT11 + T12 PuT Q̃Py T21 kH2

subject to

Q̃ ∈ RH∞

.

Now consider an intermediate TFM T̄ defined as


T11
T12 PuT
=
T̄ =
Py T21
T22

A + B2 F
−B2 F
B1
B2 PuT

0
A + LC2 B1 + LD21
0

 C1 + D12 F −D12 F
0
D12 PuT
Py D21
0
0
Py C2

(16)




.


Given that (A, B2 PuT ) is stabilizable and (Py C2 , A) is detectable, there exists matrices F̄ and L̄ such that A+B2 PuT F̄
and A + L̄Py C2 are both Hurwitz. Therefore, by choosing F
and L as PuT F̄ and L̄Py respectively in T̄ , the TFM within
the H2 norm in (16) can be rewritten as the LFT form
T11 + T12 PuT Q̃Py T21 = f (Ḡ, K̄),
where Ḡ is specified by

A
Ḡ =  C1
Py C2

B1
0
Py D21


B2 PuT
D12 PuT 
0

(17)

(18)

and K̄ is any stabilizing controller for Ḡ parameterized
by Theorem 3.1. By (17), the problem (16) becomes an
unconstrained H2 design for Ḡ, which yields the optimal
controller K̃opt as in (13). The optimal hierarchical controller, therefore, follows as Kopt = PuT K̃opt Py .
Note that the hierarchical structure of Kopt is defined
by pre- and post-projections PuT and Py over a lowerdimensional controller K̃opt . This allows applications of
conventional unstructured model reduction techniques for
designing K̃opt , without breaking the hierarchical structure
of Kopt . We present one such reduction technique in the next
section to simplify the design of Kopt .
1 Removing Assumption A4 will add extra coupled terms to A, R and
1
R2 . Interested readers can refer to [15] for the expressions. This parameter
change, however, does not affect design procedures presented in this paper.

IV. H AMILTONIAN -BASED A PPROXIMATION
Constructing the optimal controller Kopt requires two
matrices X and Y from AREs (14) and (15). Given that X
and Y are defined by the same ARE operator, we will present
the results in this section based on X only. We define the
Hamiltonian matrix for (14) as


A
−M
H:=
,
(19)
−C1T C1 −AT

where M = B2 PuT R1−1 Pu B2T . The eigenvalues of H are
symmetric about the imaginary axis [15]. Supposethat H is
Z
diagonalizable and that the columns of the matrix 1
Z2 2n×n
span the stable invariant subspace of H, i.e.
   
Z
Z
(20)
H 1 = 1 Λ− ,
Z2
Z2

−
where Λ− = diag([λ−
1 , ..., λn ]) consists of all the eigenvalues of H in the left-half plane. The matrix X can then
be found by X = Z2 Z1−1 [15]. Computing the full stable
eigenspaces above, however, is subject to O(n3 ) complexity,
and thus can become unscalable for a large network. In order
to make the hierarchical design applicable, we next consider
the problem of finding an approximate and computationally
simpler solution for X, noted as X̄ (and Ȳ for Y ).
We introduce a few notations based on the stable invariant
κ×κ
subspace of H defined in (20). Let Λ−
be a principal
κ ∈ R
−
submatrix of Λ , and denote the corresponding subspace
partitions by Z1κ ∈ Rn×κ and Z2κ ∈ Rn×κ , i.e.
 −





Λκ
−
Λ =
, Z1 = Z1κ Z1κ̄ , Z2 = Z2κ Z2κ̄ .
Λκ̄−

Using the lower-dimensional Hamiltonian subspaces Z1κ and
Z2κ , one can construct an approximate solution for (14) as
T
T
X̄ = Z2κ (Z2κ
Z1κ )−1 Z2κ
.

(21)

This approximation form has been previously proposed in
[11], although the choice of the κ eigenvalues contained in
Λ−
κ was not determined. In contrast, motivated by conventional frequency-weighted model reduction, we decide the
partitioning of Λ− from the minimization problem
minimize
X̄

k(X − X̄)Wx kH2

,

(22)

where Wx = (sI − A + M X)−1 B1 is a frequency weight
representing the closed-loop output response. Due to the nonconvex nature of the model reduction, it is intractable to find
an exact minimum for (22). Therefore, we approach (22) by
an upper bound minimization derived as follows.
T
T
Theorem 4.1: Denote Eκ =Z2κ̄ −Z2κ (Z2κ
Z1κ )−1 Z2κ
Z1κ̄ .
The objective function in (22) is upper bounded by
(23)
k(X − X̄)Wx kH2 ≤ ǫkEκ kF ,
qP
n
n×n
where ǫ :=
is a Cauchy-like
i=κ+1 Cii , and C ∈ R
matrix defined by
1
Cij := − −
[Z1−1 B1 B1T Z1−T ]ij .
(24)
λi + λ−
j

Proof: See Appendix.
From Theorem 4.1, the norm of interest k(X − X̄)Wx kH2
is linearly bounded by a scalar ǫ. The minimization in (22)
can then be approached by minimizing ǫ with respect to the
choice of Λ−
κ . It is worth emphasizing that the value of Cii
is inversely proportional to the magnitude of the closed-loop
eigenvalue λ−
i for i = κ + 1, ..., n. Therefore, one heuristic
way for minimizing ǫ can be to order the eigenvalues of Λ−
−
−
as 0 < |λ−
1 | ≤ ... ≤ |λn |, and let Λκ contain the first κ
eigenvalues with the smallest magnitudes. In that case, λ−
κ+1
till λ−
n will have larger magnitudes, and thus both Cii and ǫ
will be smaller positive numbers. Hence, computation of X̄
will require only the eigenspace of these first κ eigenvalues,
which can be completed in O(nκ2 ) time using Krylovsubspace based techniques, i.e. Arnoldi algorithm [17]. This
complexity is more tractable than O(n3 ) of finding an exact
solution X. The matching between X and X̄ becomes better
−
as the gap between λ−
κ and λκ+1 increases. If
−
−
−
0 < |λ−
1 | ≤ ... ≤ |λκ | ≪ |λκ+1 | ≤ ... ≤ |λn |,

(25)

which means that closed-loop network is naturally clustered
into κ coherent groups [1], ǫ becomes negligible, indicating
a close matching between X and X̄. The computational
complexity for X̄ can be further reduced if κ is a sufficiently
small number.
The approximate solution X̄, however, does not come
without a drawback. In general, X̄ can only stabilize A −
M X̄ over the subspace Z1κ , and hence a small κ will
not necessarily stabilize an unstable plant. Therefore, we
briefly comment on the stabilizability condition of X̄ to
conclude this section. First, some useful properties of X̄ are
summarized as follows.
Lemma 4.2: [11] The approximate solution X̄ given by
(21) satisfies X̄  X, and its residue can be written as
R(X̄) := AT X̄ + X̄A + C1T C1 − X̄M X̄ = C̄1T C̄1 , (26)
T
T
]C1T .
Z1κ )−1 Z1κ
where C̄1T = [I − Z2κ (Z2κ
A simple stability condition then follows immediately from
Lemma 4.2.
Theorem 4.3: The approximate solution X̄ is stabilizing
if C1T C1 − C̄1T C̄1  0, and (C1T C1 − C̄1T C̄1 , A) has no
unobservable modes on the imaginary axis.
Proof: By moving C̄1T C̄1 to the left hand side of
(26), (26) becomes an ARE in the form of AT X̄ + X̄A +
C1T C1 − C̄1T C̄1 − X̄M X̄ = 0. The stabilizability for X̄
then follows from standard conditions in [15] based on
(C1T C1 − C̄1T C̄1 , A).
Theorem 4.3 implies that one sufficient condition for
stabilizing an unstable system is to choose a positive-definite
matrix C1T C1 with a large norm. The flip side, however,
is that this would lead to a larger control effort, thereby
amplifying the measurement noise in the H2 design. An
alternative is to gradually increase the value of κ, and check
for closed-loop stability. Theorem 4.3 can thus be used as
a stability test, requiring computation of only Yκ and Zκ
without any need for checking the eigenvalues of A − M X̄.

V. D ESIGN

OF

H IERARCHICAL C ONSTRAINT S

The results presented so far produces a stabilizing controller Kopt , which minimizes the values of kf (G, K)kH2
for any given choice of S in (4). In this section, we show
that it is also possible to design S such that this norm is
close to the unconstrained H2 problem, that is the problem
in (4) without the constraint K ∈ S.
A. Quantification of the Optimality Gap
The unconstrained H2 problem can be written in model
matching form as
minimize
subject to

kJ1 (Q)kH2 = kT11 + T12 QT21 kH2

. (27)

Q ∈ RH∞

We start with the following two lemmas.
Lemma 5.1: The optimal solution of (27) is given by
Q∗ = −WL W̄R = −W̄L WR ,
where WL , W̄L , WR and W̄R are defined as



Â + B̂2 F̂
Â + B̂2 F̂ I
, W̄L =
WL =
0
F̂
B̂2 F̂



Â + L̂Ĉ2
Â + L̂Ĉ2 L̂
, W̄R =
WR =
I
0
L̂Ĉ2

(28)

F̂
F̂
L̂
L̂



,
,

T −1
T
L̂= − Ŷ Ĉ2T (D21 D21
) , Ŷ =ARE(ÂT , Ĉ2T , B̂1T , D21 D21
),

and matrices Â, B̂1 , B̂2 , Ĉ1 and Ĉ2 are defined in (10).
Proof: The solution for the model matching problem
(27) follows from spectral factorization. See [15].
Lemma 5.2: The constrained H2 problem (12) is equivalent to the unconstrained problem

subject to

kJ2 (Q)kH2 = kT12 PuT Pu QPyT Py T21
Q ∈ RH∞

+ T11 kH2

. (29)

Proof: See Appendix.
With Lemma 5.1 and 5.2, we can quantify the optimality
gap between the constrained and unconstrained H2 problems
as follows.
Theorem 5.3: Denote the optimal values of (27) and (29)
by J1∗ and J2∗ respectively. The optimal value J2∗ is upper
bounded by
2
2
J2∗
≤ J1∗
+ 2ξJ1∗ + ξ 2 ,

ξ = ǫ1 ξu + ǫ2 ξy + ǫ2 ξy , (30)

where ǫ1 = kT12 kH∞ kT21 kH∞ kW̄R kH∞ and ǫ2 =
kT12 kH∞ kT21 kH∞ kW̄L kH∞ are bounded scalars, and ξu
and ξy are defined by
1

ξu = k(I−PuT Pu )F̂ Φu2 kF ,

minimize
Pu

ξu

&

1

ξy = k(I−PyT Py )L̂T Φy2 kF ,

Py

ξy

.

(31)

The projections Pu and Py are defined over two variables, i.e. the clustering sets (I u , I y ) and clustering weights
(wu , wy ). Under the scope of this paper, we consider solving
only the clustering sets (I u , I y ) from (31), as these sets
define the combinatorial structure of S. To gain full degree of
freedom in designing (I u , I y ), we make the following choice
of (wu , wy ) to satisfy the stabilizability and detectability
conditions required for Theorem 3.3.
Lemma 5.4: (A, B2 PuT ) is stabilizable and (Py C2 , A) is
detectable if wu = Vl vu and wy = Vr vy for any vu and vy
that satisfy
VrT C2T C2 Vr vy 6= 0,

(32)

where 0 is a vector with all zero entries. Matrices Vl and
Vr denote the left and right eigenvector matrices such that
AT Vl = Vl Λ+ , AVr = Vr Λ+ , and Λ+  0 is a diagonal matrix that contains all unstable eigenvalues (including
imaginary eigenvalues) of A.
The proof is omitted as this lemma follows directly from
PBH test. The implication of Lemma 5.4 is that when the
open-loop system (1) is not stable, vectors wu and wy can
be chosen to meet linear inequalities in (32), which can be
done by linear programming. When the open-loop system
(1) is stable, the stabilizability and detectability conditions
become trivial. In this case, one can choose wu and wy to
be any non-zero vectors, e.g. a vector of all ones.
Once (wu , wy ) are fixed, it has been shown in our recent
paper [1] that the minimization (31) with respect to (I u , I y )
is equivalent to an unsupervised clustering optimization. A
simple yet efficient heuristic algorithm to solve this problem
is weighted k-means [16]. Therefore, one can design the clus1
1
tering sets by simply providing data matrices (F̂ Φu2 , L̂T Φy2 ),
weight vectors (wu , wy ), and number of clusters r as inputs
to the k-means algorithm, i.e.
1

1

I u = kmeans(F̂ Φu2 , wu , r) & I y = kmeans(L̂T Φy2 , wy , r).
(33)
VI. N UMERICAL E XAMPLE
In this section we verify the performance of our proposed
designs using a first-order consensus network. The plant G22
is described by ns = 500 integrators that are interconnected
by the network topology shown in Fig. 2, where each
node represents a single integrator. The dynamics at each
integrator is modeled by

Φu = LYAP(Â + B̂2 F̂ , I), Φy = LYAP((Â + L̂Ĉ2 )T , I).
The equality of (30) is attained at ξu = ξy = 0.
Proof: See Appendix.
Smaller values of ξu and ξy will result in a smaller value
of ξ, and therefore, in a tighter gap between J1∗ and J2∗ .

minimize

B. Design of Pu and Py

VlT B2 B2T Vl vu 6= 0,


T
T
F̂ = −(D12
D12 )−1 B̂2T X̂, X̂ = ARE(Â, B̂2 , Ĉ1 , D12
D12 ),

minimize

Also note that ξu and ξy are monotonic functions of Pu and
Py , respectively. Therefore, as a relaxation we design Pu and
Py from

ẋi = ui +

ns
X

aij (xj − xi ),

i = 1, ..., ns ,

j=1,j6=i

where ui is a scalar control input, and aij ≥ 0 represents
the connection strength between the ith and j th oscillators.

V4

1.3
1.2

1

1.1
0.5

1
0.9

Normalized H2 norm

Computation time /s

V1

1.4

1.5

0
0

1

2

3

4

5

6

7

κ

Fig. 3: Performance results for Hamiltonian-based reduction.
V3
10

V2

The scalability of our design is further tested upon consensus networks whose dimensions range from 100 to 10000
with aij generated randomly. The resulting computation

Computation time /s

We assume the state xi to be measurable, i.e., yi = xi . The
design parameters C1 and B1 are chosen as identity matrices
scaled by 10, and D12 and D21 are chosen as identity
matrices. For the hierarchical constraint S, the weight vectors
(wu , wy ) are chosen as vectors of all ones, which allow
Pu and Py to satisfy the stabilizability and detectability
conditions in Theorem 3.3. The clustering sets I u and I y
are selected as four groups of integrators indexed by V1 =
{1, ..., 125}, V2 = {126, ..., 250}, V3 = {251, ..., 375} and
V4 = {376, ..., 500}, where inside the groups integrators are
densely connected while the groups themselves are sparsely
connected, as shown in Fig. 2. This results in four naturally
coherent clusters. Given this constraint S, the optimal hierarchical controller Kopt requires 3.6418 seconds to solve, and
yields an H2 norm kf (G, Kopt )kH2 = 2.3191.
We first verify the Hamiltonian-based approximation with
respect to choices of κ from 1 to 6, and plot the corresponding computation costs and H2 norms yielded by the
approximated controllers in Fig. 3. For comparison, H2
norms shown on the right axis are normalized/divided by
2.3191 from the optimal case. As κ increases, the normalized
H2 norm decreases, and becomes close to 1 when κ ≥
4, indicating a close performance matching between the
approximated controller and Kopt . On the left axis, the
computation cost for the approximated controller scales up
as κ increases. When κ ≤ 4, the computation costs are all
under 0.3 seconds, which are much lower compared with
3.6418 seconds for Kopt . Therefore, at the intersection point
κ = 4, one can design an approximated controller in a much
lower cost, yet still get the similar performance as Kopt .
This performance is mainly facilitated by the 4 smallest
dominant eigenvalues resulting from the natural clustering
of the open-loop network shown in Fig. 2. The 4 dominant
eigenvalues also explain the sudden increase in computation
time at κ = 5, where the convergence of Arnoldi algorithm
become slower in searching the one extra non-dominant
eigenvalue. It is worth mentioning that if the open-loop
system is not naturally clustered, the H2 norm can increase.
The computation cost, however, will not change significantly.

Approximate
Optimal

10

2

10

0

10

-2

10 2

10 3

10 4

Dimension: n

Fig. 4: Scalability results for computation time

costs for solving the Hamiltonian-based approximation versus that of Kopt are plotted in Fig. 4. As n scales up,
the approximate approach becomes remarkably cheaper than
the optimal case. For example, at n = 4000, solving the
optimal case is already beyond our computing capacity, while
the approximate approach can be completed within only
17 seconds. This verifies the tractability of our design in
handling large-scale networks.
We finally illustrate the design of S for the same 500node network. The clustering sets (I u , I y ) in this example
are determined from (33) with respect to r = 1, ..., 6, and
then used to design the controller Kopt . We show in Fig. 5 the
ratio between J2∗ and J1∗ , which as defined in Section V are
the resulting optimal values of constrained and unconstrained
H2 problems respectively. The ratio J2∗ /J1∗ as shown in
Fig. 5 decreases monotonically, and becomes close to 1 at
r = 4. The clusters generated by r = 4 actually resembles
the four coherent groups in Fig. 2. This verifies our proposed
design in finding a hierarchical structure that tightens the gap
between constrained and unconstrained H2 problems.
VII. C ONCLUSION
In this paper we presented a hierarchical H2 optimal
controller that provides benefits for both design and implementation for large-scale network dynamic systems. The
1.8

H2 norm: J2∗ /J1∗

Fig. 2: Network topology for a 500-node network.

4

1.6
1.4
1.2
1
0

1

2

3

4

5

Number of clusters: r

Fig. 5: Optimality results for design of S

6

advantage of the hierarchical structure is that it satisfies
quadratic invariance property for any generic network model,
and allows a convex reformulation for the original problem.
The reformulated design can be simplified by conventional
controller reduction techniques without breaking the hierarchical structure of the controller. In addition, we proposed a
preliminary design of the clustering sets to tighten the gap
between unconstrained and constrained H2 problems. Our
future work will be to extend this design to the clustering
weights.
A PPENDIX
Proof of Theorem 4.1
By applying the definition of H2 norm, the objective
function in (22) can be rewritten as
1

k(X − X̄)Wx kH2 = k(X − X̄)Φ 2 kF ,
1

T

where Φ = Φ 2 Φ 2 = LYAP(A − M X, B1 ). Note that from
(20) we can find that A − M X = Z1 Λ− Z1−1 , which means
Λ− contains all the eigenvalues of closed-loop state matrix
A − M X on the eigenspace Z1 . Hence, the controllability
Gramian Φ follows directly from [18] as Φ = Z1 CZ1T . Also
given the expansion of the error matrix (X − X̄) as


T
T
Z1κ )−1 Z2κ
Z1κ̄ Z1−1 ,
(X − X̄) = 0 Zκ̄ − Z2κ (Z2κ
1

the norm k(X − X̄)Φ 2 kF follows
q
1
k(X − X̄)Φ 2 kF = tr[(X − X̄)Φ(X − X̄)]
q
1
2
= tr(Eκ Cκ+1:n,κ+1:n EκT ) = kEκ Cκ+1:n,κ+1:n
kF
q
1
2
kF = kEκ kF tr(Cκ+1:n,κ+1:n ).
≤ kEκ kF kCκ+1:n,κ+1:n
Expanding the trace above then yields (23).
Proof of Lemma 5.2

Given that Pu and Py are orthonormal matrices, it can
be verified from PBH test that (A, B2 PuT Pu ) is stabilizable
and (PyT Py C2 , A) is detectable. Following the same rationale
in the proof of Theorem 3.3, J2 can be rewritten as J2 =
f (Ğ, K̆), where Ğ follows


A
B1
B2 PuT Pu
Ğ = 
C1
0
D12 PuT Pu 
T
T
Py Py C2 Py Py D21
0
and K̆ is any stabilizing controller for Ğ parameterized by
Theorem 3.1. An important property is that the inversion of
T
T
the matrix Pu D12
D12 PuT (or Py D21 D21
PyT ) follows
T
T
(Pu D12
D12 PuT )−1 = Pu (PuT Pu D12
D12 PuT Pu )−1 PuT .

Given this equality, we can find that
T
ARE(A, B2 PuT Pu , C1 , PuT Pu D12
D12 PuT Pu ) = X,
T
ARE(AT , C2T PyT Py , B1T , PyT Py D21 D21
PyT Py ) = Y,

from which it can be verified that the state space solution of
(29) is same as Kopt from Theorem 3.3. Therefore, (29) is
equivalent to (12).

Proof of Theorem 5.3
Denote the complement of matrices Pu and Py by P̄u and
P̄y respectively, the function J1 (Q) can be written as
J1 (Q) = T11 + T12 (PuT Pu + P̄uT P̄u )Q(PyT Py + P̄yT P̄y )T21 .
By letting Q = Q∗ , we can write
1

1

2
+ kJe k2H2 + 2kJe2 J12 (Q∗ )k2H2 ,
kJ2 (Q∗ )k2H2 = J1∗

where Je =T12 PuT Pu Q∗ P̄yT P̄y T21 +T12 P̄uT P̄u Q∗ PyT Py T21 +
T12 P̄uT P̄u Q∗ P̄yT P̄y T21 . Note that Je = 0 when ξu = ξy = 0,
or equivalently P̄uT P̄u Q∗ = Q∗ P̄yT P̄y = 0. As a result, the
equation above yields the equality J1∗ = kJ2 (Q∗ )k2H2 =
J2∗ . The bound in (30) simply follows the inequality
1

1

2
2
+ kJe k2H2 + 2kJe2 J12 (Q∗ )k2H2
J2∗
≤ kJ2 (Q∗ )k2H2 = J1∗

2
2
+ ξ 2 + 2ξJ1∗ .
≤ J1∗
+ kJe k2H2 + 2J1∗ kJe kH2 ≤ J1∗

R EFERENCES
[1] N. Xue and A. Chakrabortty, “Optimal Control of LargeScale Networks using Clustering Based Projections,” submitted to Transactions on Automatic Control, available online:
http://arxiv.org/abs/1609.05265.
[2] P. Kokotovic, H. K. Khalil, and J. Oreilly, Singular Perturbation
Methods in Control: Analysis and Design, SIAM, 1999.
[3] 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.
[4] G. Obinata and B. D. Anderson, Model Reduction for Control System
Design, Springer Science & Business Media, 2012.
[5] D. D. Siljak, Decentralized Control of Complex Systems, Courier
Corporation, 2011.
[6] 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.
[7] A. Lamperski and L. Lessard, “Optimal Decentralized State-Feedback
Control with Sparsity and Delays,” Automatica, vol. 58, pp. 143-151,
2015.
[8] S. Sabau and N. C. Martins, “Youla-Like Parametrizations Subject to
QI Subspace Constraints,” IEEE Transactions on Automatic Control,
vol. 59(6), pp. 1411-1422, June 2014.
[9] 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,
2013.
[10] N. Xue and A. Chakrabortty, “LQG Control of Large Networks: A
Clustering-Based Approach,” in proceedings of the American Control
Conference (ACC), Seattle, WA, 2017.
[11] 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.
[12] Y. S. Wang and N. Matni, “Localized LQG Optimal Control for LargeScale Systems,” 2016 American Control Conference (ACC), pp. 19541961, Boston, MA, USA, 2016.
[13] D. H. Nguyen and S. Hara, “Hierarchical Decentralized Stabilization
for Networked Dynamical Systems by LQR Selective Pole Shift,”
IFAC Proceedings Volumes, vol. 47(3), pp. 5778-5783, 2014.
[14] T. Sadamoto, T. Ishizaki, and J. I. Imura, “Hierarchical Distributed
Control for Networked Linear Systems,” 53rd IEEE Conference on
Decision and Control (CDC), pp. 2447-2452, Los Angeles, CA, 2014.
[15] K. Zhou and J. C. Doyle, Essentials of Robust Control, vol. 180,
Prentice Hall Upper Saddle River, NJ, 1998.
[16] S. Lloyd, “Least Squares Quantization in PCM,” IEEE Transactions
on Information Theory, vol. 28(2), pp. 129-137, Mar 1982.
[17] G. H. Golub and C. F. Van Loan, Matrix Computations, JHU Press,
vol. 3, 2012.
[18] A. C. Antoulas, Approximation of Large-Scale Dynamical Systems,
SIAM, 2005.

