1

Subspace Identification of Large-Scale 1D
Homogeneous Networks

arXiv:1702.03539v1 [] 12 Feb 2017

Chengpu Yu, Michel Verhaegen, and Anders Hansson

Abstract—This paper considers the identification of large-scale
1D networks consisting of identical LTI dynamical systems. A
new subspace identification method is developed that only uses
local input-output information and does not rely on knowledge
about the local state interaction. The identification of the local
system matrices (up to a similarity transformation) is done
via a low dimensional subspace retrieval step that enables
the estimation of the Markov parameters of a locally lifted
system. Using the estimated Markov parameters, the state-space
realization of a single subsystem in the network is determined.
The low dimensional subspace retrieval step exploits various key
structural properties that are present in the data equation such
as a low rank property and a two-layer Toeplitz structure in the
data matrices constructed from products of the system matrices.
For the estimation of the system matrices of a single subsystem,
it is formulated as a structured low-rank matrix factorization
problem. The effectiveness of the proposed identification method
is demonstrated by a simulation example.
Index Terms—Large-scale 1D distributed systems, nuclearnorm optimization, two-layer Toeplitz structure.

I. I NTRODUCTION
In recent years, the fundamental system theoretical research
on large-scale interconnected systems has been intensified both
on the topic of distributed control, see e.g. [1]–[4], as well
as on the topic of distributed system identification, see e.g.
[5]–[7]. These fundamental developments are inspired by the
increasing interest from the application side, such as in fluid
mechanics [8], flexible structures [9], large adaptive telescope
mirrors [10], wind turbine farms [11], and so on. In this paper,
we focus on the identification of large-scale 1D homogeneous
networked systems. Such systems may characterize a platoon
of cars or may result by discretizing dynamical systems
described via partial differential equations [3], [12], [13].
When considering the subspace identification of local system dynamics using local data only, there is not only the
missing information about the local system state, but also the
missing state information from the neighboring systems. This
makes the local identification problem more difficult. In [14] a
subspace identification (SID) scheme for large-scale circulant
systems was developed. The specific property of circularity
C. Yu and M. Verhaegen are with the Delft Center for Systems and
Control, Delft University, Delft 2628CD, Netherlands (c.yu-4@tudelft.nl,
m.verhaegen@tudelft.nl)
A. Hansson is with the Division of Automatic Control, Department of Electrical Engineering, Linkoping University, Sweden (anders.g.hansson@liu.se)
Part of the research of M. Verhaegen was done while he was a visiting
professor of the Division of Automatic Control, Department of Electrical
Engineering, Linkoping University, Sweden.
The research leading to these results has received funding from the
European Research Council under the European Union’s Seventh Framework
Programme (FP7/2007-2013) / ERC grant agreement no. 339681.

was exploited to decompose the large model description into
simple subsystems through a global transformation of the
input and output data. In [15] a network of identical systems
interacting in a known interconnection pattern, indicated as
so-called decomposable systems, was considered. The identification of the dynamics of a single subsystem in the network
using transformed input-output data sequences gives rise to
a challenging Bilinear Matrix Inequality (BMI). In [12] the
property that the inverse of the observability Grammian is offdiagonally decaying was used to approximate the neighboring
states, influencing the local system dynamics to be identified,
via an (unknown) linear combination of locally neighboring
inputs and outputs. The selection of these neighboring input
and output quantities requires however an exhaustive search,
which is quite computationally demanding. As a complementary work to [12], a nuclear norm identification solution
was provided in [16] to separate the local dynamics and
global dynamics by exploiting their distinct rank and order
properties; however, it did not consider identification of the
interconnections between the subsystems.
In this paper, we consider the subspace identification of
a local cluster of identical LTI systems connected in a 1D
network. Due to the unknown neighboring-state information,
the considered local identification is a blind identification
problem, where some input sequences are totally unknown. In
this scenario, the existing SID approaches that aim at retrieving
a matrix whose row or column space is of interest, the latter,
e.g., being the extended observability matrix, cannot work. As
a result, we aim to develop a new identification method for
identifying the system matrices of a single subsystem (up to a
similarity transformation) in the network that uses local inputoutput data only.
The subspace identification method presented in this paper
has two major features. First, a novel identification method
is developed that uses similar data equations as many existing
SID variants, but aims for a low dimensional subspace retrieval
that allows the identification of the finite number of Markov
parameters present in this data equation instead. This is a new
alternative to classical SID methods [17] that would fail in
retrieving subspace like the column space of the extended
observability matrix of the system to be identified. This novel
SID approach is applied to the identification of a locally lifted
system with the missing neighboring-state information. By
fully exploiting the two-layer structure of the block Toeplitz
matrices in the data equation of a lifted state-space system,
a low-rank optimization problem is formulated for which the
optimal solution can yield (parts of) the true Markov parameters of the locally lifted system. Second is the retrieval of

2

‫ݑ‬ଵ ሺ݇ሻ

ȭଵ

‫ݕ‬ଵ ሺ݇ሻ

‫ݔ‬௜ିோିଵ ሺ݇ሻ

‫ݑ‬௜ିோ ሺ݇ሻ

‫ݑ‬௜ାோ ሺ݇ሻ

‫ݑ‬௜ ሺ݇ሻ

ȭ௜ାோ

ȭ௜

ȭ௜ିோ

‫ݕ‬௜ିோ ሺ݇ሻ

‫ݔ‬௜ାோାଵ ሺ݇ሻ

‫ݕ‬௜ାோ ሺ݇ሻ

‫ݕ‬௜ ሺ݇ሻ

‫ݑ‬ே ሺ݇ሻ

ȭே

‫ݕ‬ே ሺ݇ሻ

ȭ୧ Ǥ

Fig. 1. Illustration of a cluster of subsystems in a neighborhood of the
subsystem Σi with radius R. The states xi−R−1 (k) and xi+R+1 (k) are
explicitly indicated. They are like all other states unmeasurable.

the system matrices, describing the local LTI dynamics, of an
individual subsystem and its interaction with its neighboring
systems (up to a similarity transformation) from the reliably
estimated Markov parameters. The retrieval of the system
matrices of an individual subsystem is inherently a challenging
structured state-space realization problem [18] for which the
optimal solution can yield the estimates of system matrices
up to a similarity transformation. The second feature also
includes the exploitation of the shifting structure of a timevarying generalized observability matrix.
The rest of the paper is organized as follows. Section II
describes the concerned identification problems and shows the
challenge of dealing with the identification of a small cluster
of subsystems in a large-scale networked system. Section III
shows that the existing subspace identification methods break
down in dealing with the concerned identification problem.
Section IV presents a method for identifying the Markov
parameters of the locally lifted state-space models. Section V
provides a solution to the state-space realization of a single
subsystem. Section VI provides a numerical example. The
conclusions are provided in Section VII.
The following notations are adopted throughout the paper.
The lowercase (uppercase) x (X) is used to denote a vector
(matrix). The superscripts T and −1 are transpose and inverse
operators, respectively. vec (X) represents the vector stacked
by columns of X. tr (X) stands for the trace of X. In denotes
an n × n identity matrix and Om,n is an m × n zero matrix. ⊗
stands for the Kronecker product. kxk2 represents the norm of
x. kXkF and kXk∗ denote the Frobenius norm and nuclear
norm of the matrix X, respectively.
II. P RELIMINARIES

AND

P ROBLEM D EFINITION
{Σi }N
i=1

We consider the LTI systems
connected in a
homogeneous 1D network as shown in Fig. 1. The lifted statespace model of the concerned network is written as:
Σ1 :
Σi :

x1 (k + 1)
y1 (k)
xi (k + 1)

=
=
=

yi (k)
=
i = 2, · · · , N − 1,
ΣN : xN (k + 1) =
yN (k)
=

Ax1 (k) + Ar x2 (k) + Bu1 (k)
Cx1 (k) + e1 (k)
Axi (k) + Al xi−1 (k) + Ar xi+1 (k)
+Bui (k)
Cxi (k) + ei (k)

AxN (k) + Al xN −1 (k) + BuN (k)
CxN (k) + eN (k)
(1)
where xi (k) ∈ Rn×1 , ui (k) ∈ Rm×1 , yi (k) ∈ Rp×1 and
ei (k) ∈ Rp×1 are the state, input, output and measurement

noise of the i-th subsystem, xi−1 (k) and xi+1 (k) are the
neighboring states of the i-th subsystem. For xi (k), the
subscript i denotes the spatial index and k is referred to as the
time index. For a large-scale distributed system, we always
assume that N ≫ n and n > max{p, m}.
By lifting all states xi (k) into the vector x(k) as x(k) =
 T
T
x1 (k) · · · xTN (k) and doing the same for the inputs,
outputs and noises defining resp. the vectors u(k), y(k) and
e(k), the global system in Fig. 1 has the following state space
model:
x(k + 1) = Ax(k) + Bu(k)

(2)

y(k) = Cx(k) + e(k),
where A, B, C are N × N
following forms:

A

 Al
A=






B=


block matrices which have the

A
..
.



B
B
..

.
B



Ar
..

.

..

.
Al






, C = 





,

Ar 
A


C
C
..

.
C



.


To identify the system model in (2), the existing SID methods only estimate the triplet (A, B, C) up to a similarity transformation, thereby not preserving the special block-diagonal
and block-tridiagonal structure of these system matrices. In addition, the computational complexity of SID methods, which is
at least O(N 3 ) with N being the number of subsystems in the
network, may easily disqualify their use for the identification
of large-scale networks.
To deal with the high computational complexity for identifying the global system model in (2), we consider the
identification of a cluster of subsystems {Σj }i+R
j=i−R in a
neighborhood of Σi with radius R satisfying R < i < N − R
and R ≪ N , as shown in Fig. 1. The lifted state-space model
of this cluster is represented as
xi (k + 1) = AR xi (k) + B R ui (k) + DR v i (k),
y i (k) = C R xi (k) + ei (k),

(3)

with xi (k), ui (k), y i (k) the subparts of x(k), u(k), y(k) respectively from the block-rows (i−R) to (i+R); AR , B R , C R
are (2R + 1) × (2R + 1) block matrices which have forms
similar to A, B, C in (2), respectively; the (2R + 1) × 2 block
matrix DR and the vector v i (k) are defined as:


Al 0
 0
0 




xi−R−1 (k)
 ..

.
..  , v i (k) =
DR =  .
.
xi+R+1 (k)


 0

0
0 Ar
For the identification of the cluster model in (3), the
following two problems will be addressed.

3

Problem 1: The Markov parameters of the locally lifted
state-space model in (3) will be identified using only the local
input ui (k) and the local output yi (k).
This problem will be addressed in Section IV. The difficulty
of this problem w.r.t. classical subspace identification methods
in [19], [20] or their more recent variants [21], [22] is that
the adjacent state sequences xi−R−1 (k) and xi+R+1 (k) are
missing or the input sequence v i (k) of the considered cluster
model is unknown.
Problem 2: The order the local subsystem Σi , denoted
by n, as well as the local system matrices A, Ar , Al , B, C
are to be identified based on the identified Markov parameters in Problem 1. More specifically, the system matrices are to be identified up to a similarity transformation,
i.e., the estimates Â, Âr , Âl , B̂, Ĉ of A, Ar , Al , B, C satisfy
Â = Q−1 AQ, Âl = Q−1 Al Q, Âr = Q−1 Ar Q, B̂ =
Q−1 B, Ĉ = CQ with Q ∈ Rn×n being a non-singular
ambiguity matrix.
The second problem will be addressed in Section V. It is
noteworthy that, when the system order n is known, the second problem is inherently a structured state-space realization
problem as studied in [18], which turns out to be a challenging
non-convex optimization problem.
The developed method in this paper can be applied to more
general classes of large-scale identification problems than just
a 1D chain of identical LTI systems, such as 1D heterogeneous
networks consisting of clusters of identical LTI dynamical
systems [23].
In addressing the locally lifted identification problem in (3),
we stipulate the following assumption.
Assumption A.1 The global system (A, B, C) and the locally
lifted system (AR , B R , C R ) are assumed to be minimal.
The persistent excitation of the input signal u(k), which
will be used for the identifiability analysis in the sequel, is
defined below.
Definition 1. A time sequence u(k) ∈ RN m is persistently
exciting of order s if there exists an integer h such that the
(block-) Hankel matrix


u(k)
u(k + 1) · · ·
u(k + h − 1)

 u(k + 1)
u(k + 2) · · ·
u(k + h)


,

.
..
..
.
.
.


.
.
.
.
u(k + s − 1) u(k + s) · · · u(k + s + h − 2)
has full row rank for any positive integer k.
III. B REAK

DOWN OF EXISTING SUBSPACE
IDENTIFICATION METHODS

The data equation of the local state-space model (3) is given
as follows:
i
i
i
i
+ Es,h
,
+ TsDR Vs,h
Ys,h
= Os xih + TsBR Us,h
i
In this equation, the (block-) Hankel matrix Ys,h

y i (1) · · ·
yi (h)
 ..
..
i
.
Ys,h =  .
..
.
y i (s) · · · y i (h + s − 1)

(4)

is defined as




with the superscript i being the spatial index of the subsystem
Σi , the subscripts s, h respectively being the number of block
i
,
rows and the number of block columns. Analogous to Ys,h
i
i
i
we define the block-Hankel matrices Us,h , Vs,h , Es,h from the
sequences ui (k), v i (k), ei (k), respectively. The matrix TsBR is
a block Toeplitz matrix defined from the triplet (AR , B R , C R )
as


0


C RBR
0


TsBR = 
,
..
.
.
.
.


.
.
.
s−2
BR · · · C RBR 0
C R AR
and TsDR is defined in a similar way from the triplet
(AR , DR , C R ). The final matrix definitions in (4) are


CR
 C R AR 


Os = 

..


.
s−1
C R AR

and
xih = [xi (k) · · · xi (k + h − 1)] .
Existing subspace identification (SID) methods, see e.g.
[17], break down when retrieving a matrix with a subspace
i
i
like Os or xih from the available data matrices Ys,h
and Us,h
in (4).
To explain this deficiency of existing SID methods, we
consider the simple noise-free case e(k) ≡ 0. Let

−1
i,T
i,T
i
i
Π⊥
Us,h
,
U = I − Us,h Us,h Us,h
i
Then the subspace revealing matrix [17] is Ys,h
Π⊥
U , or
h
iT
i,T
i,T
if we consider the approach in [22].
Ys,h Us,h
When the unknown system input v i (k) in (3) is absent,
i
≡ 0), the column space of
namely v i (k) ≡ 0 (or Vh,s
i
Os can be retrieved from that of Ys,h
Π⊥
U , while the row
i
space of xh can be retrieved from that of compound matrix
iT
h
i,T
i,T
. However, for the case v i (k) 6= 0, the
Ys,h
Us,h
i
presence of the bilinear term TsDR Vs,h
destroys the important
subspace revealing property in existing SID methods.
Next, we shall show that the identification method in [21]
cannot handle the concerned identification problem. One could
consider the unknown input v i (k) in (3) as a ”missing input
variable”. However, for the simple case e(k) ≡ 0 and using
the notation Hv to denote the set block Hankel matrices of
i
the same structure as the matrix Vs,h
, the rank minimization
problem of [21] can be written as
 i 
Ys,h
i

.
min rank Us,h
(5)
V ∈Hv
V

By the following inequality
 i 
 i 
Ys,h
i
 ≥ rank Ys,h
rank  Us,h
,
i
Us,h
V

4

we can see that V = 0 is an optimal solution to (5). It implies
that the whole missing input sequence, not just a few missing
entries, cannot be recovered using the low-rank optimization
of [21].
To resolve this deficiency of SID methods, when considering
the identification of local clusters in a large homogeneous
network, a novel SID approach is developed in this paper.
This method differs from the existing SID methods in two
major ways.
First, the low-rank subspace retrieval step that is characteristic for many SID methods is unable to find an accurate (or
consistent) estimate of the key subspace of Os or xih . In this
paper, the low-rank optimization will be used for the accurate
estimation of (parts of) the Markov parameters in the block
Toeplitz matrix TsBR .
Second, contrary to existing SID methods that are unable
to preserve the structures in the estimated system matrices,
the proposed subspace identification method is capable of
accurately estimating the non-zero block entries (up to a
similarity transformation) of the structured system matrices
in (3), as stated in Problem 2.
IV. I DENTIFYING

THE M ARKOV PARAMETERS OF THE
LOCALLY LIFTED STATE SPACE MODEL (3)

In this section, the Markov parameters of the state-space
model in (3) will be estimated by exploiting the low-rank
property of the unmeasurable-state related terms in (4) and
the specific block Toeplitz structure of TsBR . Subsection IV-A
formulates a low rank optimization problem without considering the specific structure of TsBR . It is shown in Lemma 2 that
this low-rank optimization is unable to recover the true Markov
parameters. Subsection IV-B reformulates another low-rank
optimization problem by incorporating the specific Toeplitz
structure of TsBR . It is shown in Theorem 1 that, under some
mild conditions, the optimal solution to the proposed lowrank optimization problem can yield (parts of) the true Markov
parameters.
A. Formulation of a low-rank optimization problem
Due to the unmeasurable state sequences in a network, the
i
in (4) is unknown. However,
matrix sum Os xih + TsDR Vs,h
this matrix sum has a low-rank property that will be exploited
as a solution in the new subspace identification method.
i
Lemma 1. For the data equation (4), when h > ps or Ys,h
is
i
DR i
a fat matrix, the sum Os xh + Ts Vs,h satisfies the following
rank property

i
rank Os xih + TsDR Vs,h
≤ (2R + 1)n
(6)
+ min{(s − 1)sp, 2(s − 1)n},

where (s − 1)sp and 2(s − 1)n denote the number of non-zero
rows and columns of TsDR , respectively.
Proof: From the structures of Os and TsDR , we can
get that
rank Os xih ≤ rank (Os ) ≤ (2R + 1)n and


i
≤ rank TsDR ≤ min{(s−1)sp, 2(s−1)n}.
rank TsDR Vs,h
Thus, the result in the lemma is straightforward.

From Lemma 1, we can derive a condition to select the
parameter s in the data equation (4) and the cluster radius R,
i
defined above equation (3), such that the sum Os xih +TsDR Vs,h
is of low rank (or rank deficient). This condition reads:
(2R + 1)sp > (2R + 1)n + min{(s − 1)sp, 2(s − 1)n}. (7)
The above condition means that the number of the rows of the
i
is larger than an upper bound of its
matrix Os xih + TsDR Vs,h
rank. In practice, by fixing a value of s satisfying that s > np ,
we can always find a value of R such that the above inequality
holds. Therefore, in the sequel, we assume that the matrix sum
i
Os xih + TsDR Vs,h
is of low-rank (or rank deficient).
Denote the noise-free output as ŷi (k) = y i (k)−ei (k) and its
i
related block Hankel matrix Ŷs,h
. Based on the rank property
discussed above, a low-rank regularized optimization problem
is then proposed as follows:

B

min

i ∈H
Θs R ∈T ,Ŷs,h

h+s−1
X

kŷi (t) − y i (t)k2

t=1

(8)

i
h
i
i
R
+ λ · rank Ŷs,h
− ΘB
s Us,h ,

where T and H denote respectively the set of block Toeplitz
and block Hankel matrices with appropriate block sizes, and
the regularization parameter λ allows to make a trade-off
between the two terms in the cost function. It is remarked
that the block Toeplitz structure of T corresponds to the first
block-Toeplitz layer that is described in Subsection IV-B.
In the absence of measurement noise, it will be shown in
the following lemma that the optimal solution to (8) is nonunique.
Lemma 2. Consider the optimization problem in (8). Suppose
that the following assumptions are satisfied:
1) The global system input u(k) to (2) is persistently
exciting of order N n+s with s being the SID dimension
parameter in (4);
2) The measurement noise is absent, i.e., ŷi (k) = yi (k);
3) The matrix pairs (Al , B) and (Ar , B) have full row
rank.
Then, the optimal solution to the following rank optimization
problem is non-unique:

 i
i
R
for R + s < i < N − R − s.
min rank Ys,h
− ΘB
s Us,h
B

Θs R ∈T

(9)

The proof to the above lemma is provided in Appendix A.
The above lemma indicates that the Markov parameters, as the
block entries of TsBR , cannot be recovered even if the global
optimal solution to the low-rank optimization problem in (9)
can be found.
In the next subsection, a solution will be provided by further
constraining the structure of block entries of TsBR , due to the
specific structures of the system matrices (AR , B R , C R ) in
(3).
B. Structure constrained low-rank optimization
The block matrix TsBR in (4) has a two-layer block Toeplitz
structure, which will be utilized in the proposed identification

5

method. The first layer is the block Toeplitz structure of
TsBR with respect to its block entries C R AjR B R . The second
layer is the partial block Toeplitz structure inside the block
entries C R AjR B R , as highlighted in the following example
and lemma.
Example 1. If we take R = 3 and assume that each block in
AR , B R , C R has size 2 × 2, then a visual illustration of the
structures of the matrices {M j = C R AjR B R }3j=1 is given in
Figure 2.

Fig. 2. Illustration of the partial Toeplitz structure of the matrices {M j =
C R AjR B R }3j=1 Top-left for j = 1 with block-bandwidth 1; top-right for
j = 2 with block-bandwidth 2; bottom for j = 3 with block-bandwidth
3. The deep blue color represents zero entries. The parts surrounded by red
curves have block Toeplitz structures.

Lemma 3. Based on the matrices AR , B R , C R defined in (3),
the following hold about the matrix product M j = C R AjR B R :
1) M j is a banded block matrix, with block-bandwidth j.
2) The submatrices of M j , for j < 2R + 1, at the
l−th block row and q−th block column with l, q ∈
{1, · · · , 2R + 1}, are inside the partial block-Toeplitz
region for the index-pair (l, q) satisfying:
i + 1 ≤ l + q ≤ 4R + 3 − i.
Proof: Since B R and C R are block diagonal matrices
with constant diagonal blocks, M j has the same structure
pattern as AjR . The matrix AR can be written as
AR = I ⊗ A + J− ⊗ Al + J+ ⊗ Ar
with J− and J+ block-column shifting matrices such that
post-multiplication of a matrix with J− results in shifting the
block-columns of that matrix to the left and adding a zero
block column in the last column, while post-multiplication of
a matrix with J+ results in the opposite shift and adding a
zero block column in the first column.
The proof is completed by induction. For j = 1 and j = 2
this can be easily checked. Let us assume that AjR has a block
Toeplitz structure for the matrix blocks with the index-pair
(l, q) satisfying that j + 1 ≤ l + q ≤ 4R + 3 − j, and using

as
the above expression for AR , we can express Aj+1
R
Aj+1
= AjR AR = AjR (I ⊗ A + J− ⊗ Al + J+ ⊗ Ar ). (10)
R
By several trivial manipulations of the above matrix product,
j+1
has a partial Toeplitz structure
it can be derived that AR
for the matrix blocks with the index-pair (l, q) satisfying that
(j + 1) + 1 ≤ l + q ≤ 4R + 3 − (j + 1). This completes the
proof.
Let Tf denote the set of two-layer block Toeplitz matrices
that was outlined above for the matrix TsBR . More clearly, the
structure of TsBR is illustrated in the following example.
Example 2. If we take R = 2 and s = 4, the Toeplitz matrix
TsBR can be parameterized in terms of the parameter set
p×m
{ϕj }11
as follows:
j=1 with ϕj ∈ R


0

 M0
0
,
(11)
TsBR = 

 M1 M0
0
M2 M1 M0 0


ϕ1


ϕ1



,
ϕ1
where
M0
=




ϕ1
ϕ1


ϕ3 ϕ4

 ϕ2 ϕ3 ϕ4


, M

ϕ2 ϕ3 ϕ4
=
=
M1
2



ϕ2 ϕ3 ϕ4 
ϕ2 ϕ3


ϕ10 ϕ8 ϕ9
 ϕ6 ϕ7 ϕ8 ϕ9



 ϕ5 ϕ6 ϕ7 ϕ8 ϕ9 .



ϕ5 ϕ6 ϕ7 ϕ8 
ϕ5 ϕ6 ϕ11
Then we specialize the rank-optimization problem in (8) as
B

min

i ∈H
Θs R ∈Tf ,Ŷs,h

h+s−1
X

kŷi (t) − y i (t)k2

t=1

(12)

i
h
i
i
R
+ λ · rank Ŷs,h
U
− ΘB
s
s,h .

The difference between (8) and (12) is in the definition of the
constrained Toeplitz sets indicated by T and Tf , respectively.
Though it is only a minor notational difference, the impact
on the solution is huge. As will be outlined in Theorem 1, it
enables uniqueness of part of the solution.
For the uniqueness property, use will be made of the
following time-varying observability matrix Oj,k , which is a
sub-matrix of the extended observability matrix Oj defined in
(4), consisting of the block rows corresponding to the secondlayer block Toeplitz part of TsBR .
Definition 2. Let Gj be a j × (j + 2) block Toeplitz matrix
of the form


Al A Ar


Al A Ar


(13)
Gj = 
.
..
..
..


.
.
.
Al A Ar

6

Denote C j = Ij ⊗ C. A time-varying observability matrix
Oj,k , for j > 2(k − 1), is defined in terms of the matrix pair
(C j , Gj ) as [24, Chapter 3]:


Cj


C j−2 Gj−2




C j−4 Gj−4 Gj−2
Oj,k = 
.


..


.
C j−2(k−1) Gj−2(k−1) · · · Gj−2

Theorem 1. Suppose that the following assumptions are
satisfied:
1) Assumption A.1 holds and νo is the observability index
of the pair (AR , C R );
2) The cluster radius R of the lifted model in (3) and the
dimension parameter s in (4) satisfy
s > νo ,

R ≥ s − 1;

3) The conditions 1)-3) of Lemma 2 hold;
4) The time-varying observability matrix O2R+1,s−1 , defined in Definition 2, has full column rank.
Then the submatrices of the Markov parameters M j =
C R AjR B R , for j = 0, 1, · · · , s − 2, contained in the matrix
TsBR in (4), for which Lemma 3 has shown that they preserve
the second layer of block-Toeplitz structure, can be computed
in a unique manner from the following low-rank optimization
problem:
min

B
Θs R ∈Tf

h
i
i
i
R
rank Ys,h
− ΘB
s Us,h

for R + s < i < N − R − s.
(14)

The proof of the above theorem can be found in the
Appendix B, where we can see that the (unique) global optimal
solution can yield correct recovery of certain parts of the
Markov parameters M j = C R AjR B R , for j = 0, 1, · · · , s − 2.
In addition, based on the conditions 2) and 4) of the above
theorem, it can be established that the inequality (7) holds,
i
namely the matrix sum Os xih + TsDR Vs,h
is indeed of lowrank (or rank deficient). This in turn indicates that the low-rank
cost function in (12) is a reasonable choice.
Lemma 2 and Theorem 1 demonstrate that the second-layer
Toeplitz structure, as pointed out in Lemma 3, is crucial for
enforcing a (unique) solution to the identification problem 1
of this paper. The considered identification of the Markov
parameters boils down to solving the low-rank regularized
optimization problem in (12).
Since the optimization problem in (12) is non-convex,
it is difficult to obtain an optimal solution under a mild
computational burden. In this paper, the reweighted nuclear
norm optimization method [25] is adopted, which is regarded
as an iterative heuristic for the rank minimization problem
(12).
V. S TATE - SPACE REALIZATION OF

A SINGLE SUBSYSTEM

In this section, we shall study the final realization of
the system matrices {C, A, Al , Ar , B} from the estimated
submatrices of the Markov parameters M j = C R AjR B R , for
j = 0, 1, · · · , s − 2, contained in the matrix TsBR in (4), for

which Lemma 3 has shown that they preserve the second-layer
of block-Toeplitz structure.
Remark 1. The realization problem considered in this section
is a structured state-space realization problem, which turns out
to be a challenging non-convex optimization problem [18],
[20], [26]. This problem is usually solved by the gradientbased optimization methods, which are sensitive to the selection of the initial condition. Here, we transform the realization
problem into a structured low-rank matrix factorization problem for which the optimal solution enables the estimation of
system matrices (C, Al , A, Ar , B) up to a similarity transformation. This factorization problem is then reformulated into
an equivalent low-rank optimization problem, which is then
numerically solved using the optimization method developed
in [27].
We start the solution by developing expressions of the
submatrices in the second-layer Toeplitz regions in terms of
the system matrices {C, A, Al , Ar , B}. This is done in the
following Lemma.
Lemma 4. Consider the block matrices AR , B R and C R
defined in (3). Let the sequence of non-zero block entries from
left to right of the (j +1)-th block row of the matrix C R AjR B R
be denoted as {Fj,−j , Fj,1−j , · · · , Fj,j−1 , Fj,j }, then these
matrix entries satisfy the following relationship:
j
X

Fj,k z −k = C(Al z −1 + A + Ar z)j B,

(15)

k=−j

where z ∈ C.
Proof: The above result can be derived using the filter
bank theory in [28].
As Fj,k are the Markov parameters inside the second-layer
block Toeplitz part of TsBR , the values of Fj,k for j ∈
{0, 1, · · · , s−2}, k ∈ {−j, · · · , j} are assumed to be available
in this section. Based on these values Fj,k , we address the
problem of estimating the matrices {C, A, Al , Ar , B} up to a
similarity transformation.
Dual to Definition 2, we shall define a time-varying controllability matrix Cj,k , which is a sub-matrix of the extended
controllability matrix determined by (AR , B R ).
Definition 3. Let Γj be a (j + 2) × j block Toeplitz matrix of
the form


Ar

 A Ar




.
..

 Al A


(16)
Γj = 
.
..


.
A
A
r 
l



..

. A 
Al
Denote B j = Ij ⊗ B. A time-varying controllability matrix
Cj,k , for j > 2(k − 1), is defined in terms of the matrix pair
(Γj , B j ):
Cj,k = [ B j

Γj−2 B j−2

···

Γj−2 · · · Γj−2(k−1) B j−2(k−1) ] .

7

In the sequel, we let s be an even integer such that s/2
is an integer as well. The solution to the realization of the
system matrices {C, Al , A, Ar , B} is done in two phases.
In the first phase, the structured time-varying observability
matrix O2R+1,s/2 and the structured time-varying controllability matrix C2R+1,s/2 are to be estimated from the available
matrix values Fi,j . In the second phase, the system matrices {C, A, Al , Ar , B} are derived from these time-varying
observability and controllability matrices. It is remarked that
the subscript s/2 in O2R+1,s/2 (or C2R+1,s/2 ) means that
the maximum moment of the block entries in O2R+1,s/2 (or
C2R+1,s/2 ) is s/2−1. The subscript s/2 is adopted because the
sum of the maximum moments of O2R+1,s/2 and C2R+1,s/2 is
equal to the maximum moment, s− 2, of the available Markov
parameters M j .
A. Determining the time-varying observability and controllability matrices
In this subsection, the determination of O2R+1,s/2 and
C2R+1,s/2 will be formulated as a structured low-rank matrixfactorization problem. More importantly, we show that the
optimal solution to this matrix-factorization problem can
yield the estimates of O2R+1,s/2 and C2R+1,s/2 up to a
block-diagonal ambiguity matrix with identical block-diagonal
entries. This is crucial to identify the system matrices
{C, A, Al , Ar , B} up to a similarity transformation, as stated
in Problem 2.
First, by the definitions of O2R+1,s/2 and C2R+1,s/2 , we
can find that the product of O2R+1,s/2 and C2R+1,s/2 is equal
to a matrix constructed from {Fj,k }jk=−j for j = 0, 1, · · · , s−
2. This is demonstrated in the following example.
Example 3. When R = 1 and s = 4, the product of O3,2 and
C3,2 can be expressed as




C
0
0
 B 0 0 Ar B
 0
C
0
  0 B 0 AB 
O3,2 C3,2 = 
 0
0
C 
0 0 B Al B
CAl CA CAr


F0,0
0
0
F1,1
 0
F1,0 
F0,0
0
.
=
 0
0
F0,0 F1,−1 
F1,−1 F1,0 F1,1 F2,0
The above equation provides a simple example, showing that
the product of O3,2 and C3,2 can be expressed in terms of
{Fj,k }jk=−j for j = 0, 1, 2.
Here, the product of O2R+1,s/2 and C2R+1,s/2 is represented as:
O2R+1,s/2 C2R+1,s/2 = H2R+1,2R+1 ,
(17)
where H2R+1,2R+1 is a (2R + 1) × (2R + 1) block matrix
that is assumed to be known.
Given the matrix H2R+1,2R+1 , the problem of interest is
to determine O2R+1,s/2 and C2R+1,s/2 from equation (17).
According to Definitions 2 and 3, we can see that O2R+1,s/2
and C2R+1,s/2 are structured matrices. These structures are
instrumental to analyzing the properties of the optimal solution

to (17), as shown in Theorem 2. In order to show these
structures, use will be made of the following matrix definition.
For j = 0, 1, · · · , s/2 and z ∈ C, the matrix sequences
{Wj,l }jl=−j and {Ej,l }jl=−j are defined as
j
X

Wj,l z −l = C(Al z −1 + A + Ar z)j

l=−j
j
X

(18)
Ej,l z

−l

= (Al z

−1

j

+ A + Ar z) B.

l=−j

From the defined quantities Wj,l and Ej,l in (18), the matrices
Wj and Ej are defined as


T

T
E0,0
W0,0
T
 E1,−1

 W1,−1 




T




E1,0
W1,0




T




E1,1
W1,1
(19)
Wj = 
 Ej = 
 .
T
 E2,−2

 W2,−2 








..
..




.
.
Wj−1,j−1

T
Ej−1,j−1

It can be seen that Oj,k (or Cj,k ) is a block matrix constructed
from the block components of Wk (or Ek ). This is illustrated
in the next example.
Example 4. The time-varying observability matrix O5,3 is
represented in terms of W3 as follows


W0,0


W0,0




W
0,0




W
0,0


W
O5,3 = 
0,0 
.


 W1,−1 W1,0
W
1,1




W
W
W
1,−1
1,0
1,1



W1,−1 W1,0 W1,1 
W2,−2 W2,−1 W2,0 W2,1 W2,2
Denote by O2R+1,s/2 the set of block matrices having the
same structure as O2R+1,s/2 , as illustrated in Example 4 , and
C2R+1,s/2 the set of block matrices having the same structure
as C2R+1,s/2 . It is remarked that the definitions of both the
sets O2R+1,s/2 and C2R+1,s/2 require knowledge of the system
order n. We then propose the following structured low-rank
matrix factorization problem:
min kH2R+1,2R+1 − OCk2F
O,C

(20)

s.t. O ∈ O2R+1,s/2 , C ∈ C2R+1,s/2 .
It will be shown in the following theorem that the optimal
solution to (20) can yield the estimates of O or C up to a
block-diagonal ambiguity matrix with identical block-diagonal
entries.
Theorem 2. Consider the optimization problem in (20). Suppose that the following assumptions are satisfied:
1) The values of R and s satisfy R ≥ s − 2, and s is a
positive even integer;
2) The system order of a subsystem, n, is known;

8

3) The matrices Oj,s/2 and Cj,s/2 , for any j ≥ min{2s −
3, 2R}, have full column and row rank, respectively;
4) The matrix H2R+1,2R+1 satisfying equation (17) is
known exactly.
Then, any optimal solution pair {Ô, Ĉ} to the optimization
(20) satisfies
Ô = O2R+1,s/2 Q
Ĉ = Q

−1

C2R+1,s/2 ,

(21)

where Q = I2R+1 ⊗ Q with Q ∈ Rn×n being a nonsingular
ambiguity matrix.
The proof of this theorem is given in Appendix C.
B. Rank-constrained form of the optimization problem (20)
In this subsection, the structured low-rank matrix factorization problem in (20) is reformulated into a rank-constrained
optimization problem.
From the matrix quantities defined in (18) and (19), we
can see that O2R+1,s/2 is affine in terms of Ws/2 , while
C2R+1,s/2 is affine in terms of Es/2 . Then, it can be established that product O2R+1,s/2 C2R+1,s/2 is affine in terms of
block entries of Ws/2 Es/2 [29]. This affine operator H(·) is
defined by
O2R+1,s/2 C2R+1,s/2 = H(Ws/2 Es/2 ).

(22)

In the sequel, the affine operator H(·) is assumed to be known.
Instead of O2R+1,s/2 and C2R+1,s/2 , we regard the product
X = Ws/2 Es/2 as the variable to be determined. Then, the
affine function H(X) works on the block entries of X. By
taking into account the low-rank property of the matrix product Ws/2 Es/2 , we propose the following rank-constrained
optimization problem:
min kH2R+1,2R+1 − H(X)k2F
X

s.t. rank(X) = n,

(23)

where n is the system order of a subsystem in (1) and X has
the same size as Ws/2 Es/2 .
Next, we will show that the rank-constrained optimization
problem in (23) is an equivalent formulation of (20). To show
this result, the following lemma is required.
Lemma 5. Suppose that Oj,s/2 and Cj,s/2 , for j > s − 2,
have full column and row rank, respectively. Then, Ws/2 and
Es/2 have full column and row rank, respectively, and both
ranks are n.
Proof: Since Oj,s/2 has full column rank, the block
iT
h
T
T
T
, which is constructed
vector W0,0
, W1,−1
, · · · , Ws/2,−s/2
by stacking the non-zero block entries in the first block column
of Oj,s/2 , has full column rank as well. Since this block vector
is a sub-part of Ws/2 , it can be derived that Ws/2 has full
column rank. Similarly, it can be proven that the block vector
Es/2 has full row rank as well.
The next theorem shows the equivalence between the optimization problems in (20) and (23).

Theorem 3. Consider the optimization problems in (20) and
(23). Suppose that all the conditions of Theorem 2 hold. Then,
the optimization problems (20) and (23) are equivalent in
the
can yield the estimate
 sense that their optimal solutions

Ŵs/2 , Ês/2 of Ws/2 , Es/2 up to an ambiguity matrix,
i.e., there exist a nonsingular matrix Q ∈ Rn×n such that
Ŵs/2 = Ws/2 Q, Ês/2 = Q−1 Es/2 .
Proof: First, we will show that the optimal
 solution
to (20) can yield the estimate of Ws/2 , Es/2 up to an
ambiguity matrix.
From the matrix quantities in (18) and (19), we can see
that structured matrices O2R+1,s/2 and C2R+1,s/2 are linearly
parameterized by the block components of Ws/2 and Es/2 ,
respectively. By the result of Theorem 2, it can be derived
that the optimal
 solution to (20) can yield the estimate of
Ws/2 , Es/2 up to an ambiguity matrix.
Second, we will show that the optimal
solution to (23) can

yield the estimate of Ws/2 , Es/2 up to an ambiguity matrix.
By equation (22), the optimization problem (20) can be
reformulated as
min

W,E

kH2R+1,2R+1 − H(WE)k2F ,

(24)

where the variables W and E have the same sizes of Ws/2 and
Es/2 , respectively. It is noted that equation (24) is a parameterized form of the structured low-rank matrix factorization
problem in (20).
By Theorem 2 and Lemma
 5, the optimal solution
(Ŵ, Ê) to (24) satisfies rank ŴÊ = n and H(ŴÊ) =

H2R+1,2R+1 . It can be derived that X = ŴÊ is an optimal
solution to (23) and the criterion of (23) becomes zero.
Let X̂ be an optimal solution to (23). It should satisfy
H(X̂) = H2R+1,2R+1 and rank X̂ = n. Let the SVD

decomposition of X̂ be given by X̂ = Û ΣV̂ T , where Û and
V̂ are constructed by n orthogonal columns, and Σ ∈ Rn×n
is a diagonal matrix with positive diagonal entries. Then,
(Ŵ, Ê) = (Û , ΣV̂ T ) is an optimal solution
 to (24). Therefore,
(Û , ΣV̂ T ) is an estimate of Ws/2 , Es/2 up to an ambiguity
matrix. The proof is then completed.
The rank-constrained optimization problem in (23) is nonconvex and NP-hard. Following our previous work [27], the
rank-constrained optimization problem (23) is recast into
a difference-of-convex optimization problem which is then
solved using the sequential convex programming method.
Let X̂ be an optimal solution to (23). The SVD decomposition of X̂ is given by
 T 

 Σ1

V1
,
(25)
X̂ = U1 U2
Σ2
V2T

where U1 and V1 consists of n orthogonal columns, and the
diagonal matrix Σ1 ∈ Rn×n has larger diagonal entries than
Σ2 . The estimates of Ws/2 and Es/2 can then be obtained as
follows:
Ŵs/2 = U1 ,
Ês/2 = Σ1 V1T .

9

Since O2R+1,s/2 and C2R+1,s/2 are respectively affine in
terms of Ws/2 and Es/2 , the variables O and C in (20) can
be estimated from Ŵs/2 and Ês/2 , respectively.
C. Determining the system matrices {A, Al , Ar , B, C}
In view of the theoretical result in Theorem 2, we assume
that the obtained estimates Ô2R+1,s/2 and Ĉ2R+1,s/2 satisfy
Ô2R+1,s/2 = O2R+1,s/2 Q,
Ĉ2R+1,s/2 = Q−1 C2R+1,s/2 ,

(26)

where Q = I2R+1 ⊗ Q with Q ∈ Rn×n being nonsingular.
Based on these estimates, we will address the identification
of the system matrices {A, Al , Ar , B, C} up to a similarity
transformation.
First, the shifting structure of the time-varying observability
matrix O2R+1,s/2 will be explored. Denote


C j−2k1 Gj−2k1 · · · Gj−2
 C j−2(k +1) Gj−2(k1 +1) · · · Gj−2 
1


Oj,k1 :k2 = 

..


.
C j−2k2 Gj−2k2 · · · Gj−2
where 0 ≤ k1 < k2 ≤ s/2 − 2 and 2k2 ≤ j ≤ 2R + 1. The
matrix Oj,k1 :k2 above is constructed by the block rows of Oj,k
with block-row indices from k1 to k2 . Then, the structureshifting property of O2R+1,s/2 can be represented as
O2R−1,0:s/2−2 G2R−1 = O2R+1,1:s/2−1 ,

(27)

where O2R−1,0:s/2−2 and O2R+1,1:s/2−1 are sub-matrices of
O2R+1,s/2 , and G2R−1 is a block Toeplitz matrix defined in
Definition 2. This is illustrated in the following example.
Example 5. When R = 2 and s = 6, the structure-shifting
property in (27) can be explicitly written as

Lemma 6. Let Ô2R+1,s/2 satisfy equation (26). Assume
that Ô2R−1,0:s/2−2 has full column rank. Then, the optimal
solution Ĝ to the optimization problem in (28) satisfies

(29)
Ĝ = I2R−1 ⊗ Q−1 G2R−1 (I2R+1 ⊗ Q) ,

where Q ∈ Rn×n is a nonsingular matrix.

The above lemma can be derived straightforwardly based
on equation (26) and the optimization formulation (28).
The condition that the time-varying observability matrix
Ô2R−1,0:s/2−2 has full column rank is similar to condition 3)
of Theorem 2. Lemma 6 implies that the matrices Al , A, Ar
can be determined up to a similarity transformation, i.e.,
Âl = Q−1 Al Q, Â = Q−1 AQ, Âr = Q−1 Ar Q.
In addition, according to equation (26), the estimates Ĉ and B̂
can be extracted respectively from Ô2R+1,s/2 and Ĉ2R+1,s/2 ,
satisfying that
Ĉ = CQ, B̂ = Q−1 B.
To ease the reference, the proposed local network identification algorithm is summarized in Algorithm 1.
Algorithm 1: Local identification for 1D distributed systems
Construct a spatially stacked state-space model (3) and its
Step 1
temperally stacked equation (4) based on local observations;
B
Step 2
Estimate Ts R from the optimization problem (12);
Step 4
Estimate O2R+1,s/2 and C2R+1,s/2 using
the method described in Subsection V-B;
Step 5
Extract the estimates of C and B from the estimates of
O2R+1,s/2 and C2R+1,s/2 , respectively;
Step 6
Estimate Al , A, Ar by solving (28).

VI. N UMERICAL

SIMULATION

In this section, numerical simulations are provided to
demonstrate
the effectiveness of the proposed identification



C
method
–
Algorithm
1. In the simulation, the distributed
A
A Ar
C
 l


system
is
constructed
by
connecting 40 identical subsystems
Al A Ar

C 
Al A Ar
in
a
line,
and
the
identification
for the 20-th subsystem
CAl CA CAr

 is performed. The system matrices (A, Al , Ar , B, C) with
CAl
CA
CAr
3×3
, B ∈ R3×2 and C ∈ R2×3 are randomly
CAl
CA
CAr

 A, Al , Ar ∈ R
,
=

generated
such
that
Assumption A.1 is satisfied and the 1D
CAl
CA
CAr
CA2l C(AAl + Al A)
∗
C(AAr + Ar A) CA2r
networked system is stable.
To construct an augmented state-space model in (3), we set
where ∗ is used to represent the term C(Al Ar + A2 + Ar Al ). R = 5 and s = 8. The system input and the measurement noise
Based on equation (27), we formulate the following struc- are generated as white noise sequences and the data length is
tured least-squares optimization problem to identify the ma- set to 800.
To evaluate the identification performance against the noise
trices Al , A, Ar based on the estimate Ô2R+1,s/2 :
effect, the criterion signal-to-noise ratio (SNR) is adopted,
which is defined as
min kÔ2R−1,0:s/2−2 G − Ô2R+1,1:s/2−1 k2F
G


(28)
var(yi (k) − ei (k))
s.t. G ∈ G2R−1 ,
.
SNR = 10 log
var(ei (k))
where G2R−1 denotes a set of matrices having the same
In the sequel, we shall carry out numerical simulations with
structure as G2R−1 , as shown in Definition 2.
The optimal solution to (28) has properties shown in the the SNR ranging from 0 dB to 95 dB.
We use the criterion impulse-response fitting to evaluate
following lemma.
the performance of the proposed identification method. The

10

normalized fitting error of the impulse-response sequence
CAi B is defined by
T P10
1 X i=0 kCAi B − Ĉj Âij B̂j kF
,
P10
i
T j=1
i=0 kCA BkF

Normalized fitting error

10-1

10-2

10-3

0.02

Normalized fitting error

where T is the number of randomly generated networked
systems and {Ĉj , Âj , B̂j } are the estimates of {C, A, B} of
the j-th generated networked system, respectively. Similarly,
we can also define the normalized fitting error for the impulseresponse sequences CAil B and CAir B.
Fig. 3 shows the identification performance of the proposed
method against the SNR. The normalized fitting errors are
calculated by averaging the results of 100 randomly generated
networked systems, and the regularization parameter λ in (8)
is set to λ = 10−3 . It can be observed from Fig. 3 that, using
the proposed identification method, the normalized fitting error
decreases along with the increase of the SNR. When the SNR
is larger than 50 dB, the normalized fitting errors can be as
small as 10−4 .

0.025

0.015

0.01

0.005

0
10-8

10

20

30

40

50

60

70

80

90

VII. C ONCLUSION
The local identification of 1-D large-scale distributed systems has been studied. Compared with the classical system
identification problems, the challenging point of the local
system identification is that there are two unknown system
inputs which are the states of its neighboring subsystems. By
exploiting both the spatial and temporal structures of the distributed system, especially the two-layer Toeplitz structure of
the Markov-parameter matrix, a low-rank optimization problem has been provided for identifying the Markov parameters
of a local cluster of identical subsystems, where the associated
optimal solution can yield (parts of) the true Markov parameters. Moreover, the system realization of a structured statespace model is formulated as a structured low-rank matrix

100

A PPENDIX A
P ROOF OF L EMMA 2

100

Fig. 4 shows the impulse-response fidelity of individual
subsystems against the regularization parameter λ, where
SNR is set to SNR=40 dB. The normalized fitting errors are
computed by averaging the results of 100 randomly generated
networked systems. We can observe from Fig. 4 that good
identification performance can be achieved by choosing the
regularization parameter λ around 10−2 .

10-2

factorization problem, showing that the system matrices can
be determined up to a similarity transformation by enforcing
the structure of the generalized observability/controllability
matrix.
Although we only consider the identification of 1D homogeneous networked systems, it can be extended to 2D homogeneous networks using the same identification framework,
namely exploiting both the spatial and temporal structures
of the concerned distributed system. In our future work,
the identification of large heterogeneous networks will be
investigated.

SNR (dB)

Fig. 3. Normalized fitting errors under different noise levels with λ = 10−3 .

10-4

Fig. 4. Normalized fitting errors against the regularization parameter λ with
SNR=40 dB.

10-4
0

10-6

The proof of the lemma is divided
 into three steps.
xih
i
 has full row rank under
Step I: we shall show that  Vs,h
i
Us,h
the given assumptions in the lemma.
It can be derived from (2) that






 
I
···
 A
··· 
 
 = .
· · ·   ..
As−1
x(k + s − 1) · · ·

x(k)
x(k + 1)
..
.



×


B
..
.
As−2 B
x(k)
u(k)
..
.


..

.
···
···
···

B




.
··· 
u(k + s − 2) · · ·






(30)

Let
Ξ = [h O(2R+1)n,(i−R−1)n I(2R+1)n O(2R+1)n,(N−i−R)n ], i
O
In
On,(2R+1)n
On,n
On,(N −i−R−1)n
Ω = On,(i−R−2)n
On,n
On,(2R+1)n
In
On,(N −i−R−1)n
n,(i−R−2)n
and
Π = [ O(2R+1)m,(i−R−1)m I(2R+1)m O(2R+1)m,(N−i−R)m ]

11

be selection matrices such that

Ξ

Ω

 ΩA
ΩB

 i  
..
..
..

.
xh
.
.

i
s−1
s−2
 Vs,h  =  ΩA
ΩA
B
·
·
·
ΩB

i

Us,h
0
Π


..
..

.
.


0
Π
0


x(k)
···

u(k)
··· 




.
..
×
.
·
·
·


 u(k + s − 2) · · · 
u(k + s − 1) · · ·

and









0 





 (31)

Π

By Lemma 10.4 in [20], under the assumption that u(k) is
persistently exciting of order nN + s, it can be established
that


x(k)
···

u(k)
··· 




..

.
··· 
u(k + s − 1) · · ·

has full row rank. In addition, by explicitly unfolding the
expressions of ΩAj B, it can be verified that, when (Al , B)
and (Ar , B) have full row rank and R ≥ s − 1, the coefficient
matrix on the right-hand side of (31) has full row rank. For
the sake of brevity, we only show this for the case with
s = 3. The associated coefficient matrix in (31) is shown
in (32) with ∗ denoting unimportant entries for analyzing the
concerned rank property. Under the condition that (Al , B) and
(Ar , B) have full row rank, it is easy to see that the submatrix stacked by the first five and last two block rows of
(32) has full row rank. It can be further verified that the whole
coefficient matrix in (32) has full row rank. This
 is a conse
 2


 Al B 0
quence of Al Al B B = Al B
0 0 I



 Ar B 0
 2
with
and Ar Ar B B = Ar B
0 0 I




Al B and Ar B having full row rank.
Now that the coefficient matrix in (31)has full row rank,
xih
i
 has full row rank
Sylvesters’ inequality shows that  Vs,h
i
Us,h
as well.
R
R
Step II: Let TB
∈ T be the true value of ΘB
s
s . Next, we
BR
shall show that Ts is an optimal solution
(9). Let ∆s ∈ T .
 i to 
xh
i
, we can derive
Under the full row rank property of  Vs,h
i
Us,h
that

 i
i
R
rank Ys,h
− TB
s Us,h

 i 


xh
Os TsDR
=rank
i
Vs,h


=rank Os TsDR

 i 
 i
R
− ∆s Us,h
rank Ys,h
− TB
s
 i 

xh


i

=rank  Os TsDR ∆s  Vs,h
i
Us,h


DR
∆s .
=rank Os Ts

By the inequality



rank Os TsDR ≤ rank Os

TsDR

∆s

we can obtain that

 i
i
R
rank Ys,h
− TB
s Us,h
 i 
 i
R
− ∆s Us,h
− TB
≤ rank Ys,h
s



(33)

,

(34)

R
is an optimal solution to (9).
for any ∆s ∈ T . Thus, TB
s
Step III: In this step, we shall show that the optimal solution
to (9) is non-unique. Let

∆s = TsDR Gs ,

(35)

R
where Gs = Is ⊗GR with GR ∈ R2n×(2R+1)n . Denote TB
=
s
BR
BR
Ts − ∆s . It is easy to verify that Ts ∈ T .
Since ∆s = TsDR Gs , it can be established that

i
i
R
rank [Yk,s,h
− TB
s Uk,s,h


= rank Os TsDR ∆s




 I 0 0
DR
Os Ts
(36)
= rank
0 I Gs


= rank Os TsDR

i
i
R
= rank [Ys,h
− TB
s Us,h .
R
It can then be observed from the above equation that both TB
s
BR
and Ts are optimal solutions to (9); therefore, the optimal
solution to (9) is non-unique.

A PPENDIX B
P ROOF OF T HEOREM 1
Following the proof of Lemma 2, the following matrix in
(31) has full row rank:
 i 
xh
i
.
 Vs,h
i
Us,h
Let Tf denote the set of the two-layer block Toeplitz
R
matrices, defined in Subsection IV-B. Let TB
∈ Tf be the
s
R
true value of the estimate of ΘB
in
(14).
Since
the matrix set
s
Tf belongs to the matrix set T , the minimal value of the cost
function in (9) is smaller than or equal to that in (14), i.e.

 i
i
R
min rank Ys,h
− ΘB
s Us,h
B

Θs R ∈T

≤

 i

i
R
min rank Ys,h
− ΘB
s Us,h .

B
Θs R ∈Tf

R
∈ Tf is an
As shown in Lemma 2, the true value TB
s
optimal solution to (9). Then, it can be observed from the
R
above inequality that TB
∈ Tf is also an optimal solution to
s

12















0
0
0
0
0
A2l
0

0
0
0
Al
0
∗
0

0
In
0
A
0
∗
0

| In
| 0
| 0
| Ar
| 0
| ∗
| 0
|
|

∗
∗
∗
∗
∗
∗
∗

In
0
0
0
Al
0
∗

| 0
| 0
| In
| 0
| A
| 0
| ∗
|
|

0
0
0
0
Ar
0
∗

0
0
0
0
0
0
A2r

0
0
Al B
0

B
0
AB
0

|
|
|
0 |
∗
0 |
∗
Ar B |
∗
0 |
∗
| I(2R+1)m
|
0

|
|
|
| 0
| 0
| 0
| Al B
|
|

0
B
0
AB

0
0
0
Ar B

0
0

|
|
|
|
|
0|
∗
0|
∗
|
0
| I(2R+1)m

B
0

|
|
|
|
|
|0
|0
|
|










0 

0 



0
B

(32)

(14). As a consequence, the corresponding minimal value of
the cost function of (14) is
rank



TsDR

Os



.

BR
R
R
Let TB
∈ Tf be such that ∆s = TB
is a perturbation
s
s −Ts
we seek of the true value such that it retains the minimal value
of the cost function in (14). Using the full row rank condition
of the matrix in (31), we obtain that


i
i
R
rank Ys,h
− TB
s Us,h


i
i
+ ∆s Us,h
= rank Os xih + TsDR Vs,h


= rank Os TsDR ∆s .

(37)

∗
 ..
 .

 ∗

 0




 0



0 




0 
.
∗ 

.. 
. 

(38)

∗

(b) Permutation of ∆s : Let Ls denote the matrix that
permutes all zero block rows of TsDR to top such that
Ls TsDR =
2



TsDR

0
T2



,

Ls Os =



Os1
∗



following special structures



0
∆0,1
..
.
···

0
..

.

∆1,s−3
∗

..

.
∆0,s−2





,


0 

,
(40)

Now, we seek that part of ∆s ∈ Tf that retains the minimal
value of the cost function in (14) only when it is zero. For that
reason, we permute the block rows of ∆s compatible with the
zero block-row pattern of the matrix TsDR .
(a) Structure of TsDR : The structure of TsDR is determined
by C R AjR DR for j = 0, 1, · · · , s − 2. By the definitions of
AR , C R and DR in (3), we can find that C R AjR DR is a (2R+
1) × 2 block matrix. When R ≥ (s − 1), the block rows of
C R AjR D R indexed from (j + 2) to (2R − j) are zero. For
instance, C R AjR D R has the block structure as shown below
with ∗ representing non-zero block entries:


∆s and Os , yielding the
 1 
∆
Ls ∆s =
∗

0

∆
0,0


∆
1,0

=
..

.

 ∆s−2,0

(39)

where T is a sub-matrix of
by removing the corresponding zero block rows, as illustrated in (38). It is noteworthy that
the zero block rows of TsDR correspond to the second-layer
block Toeplitz part of TsBR , as shown in Lemma 3.
Then, the permutation matrix Ls is applied to the matrices

such that Os1 = O2R+1,s and ∆j,l , for j = 0, 1, · · · , s − 2 and
l = 0, · · · , s − 2 − j, is a [2(R − j − l) − 1] × [2R + 1] block
Toeplitz matrix in the form


0 δj,−j δj,1−j · · ·
δj,j
0


..
..
..
..
..
..
∆j,l = 
,
.
.
.
.
.
.
0

δj,−j

δj,1−j

···

δj,j

0

with {δj,k }jk=−j being block entries of appropriate sizes. Note
that ∆j,l for l ≥ 1 is a sub-matrix of ∆j,0 , which is constructed
by stacking the block rows indexed from (l+1) to (2R−1−l).
Now, we are ready to show for what values of ∆j,l the
following rank constraint holds




rank Os TsDR ∆s = rank Os TsDR .
(41)
It can be derived from the above equality that



 Q1
DR
= Ls ∆s ,
Ls Os Ts
Q2
or
  1 


∆
O2R+1,s 0
Q1
=
,
Q2
∗
T2
∗

(42)

(43)

where Q1 and Q2 are coefficient matrices of appropriate sizes,
and ∗ denotes non-zero block entries caused by the block-row
permutation. So we focus on the top part of (43) and rewrite
it as


(44)
∆1 = O2R+1,s Q1,0 Q1,1 · · · Q1,s−1 ,

s−1
where {Q1,j }j=0
are sub-matrices of Q1 satisfying that Q1 =
[Q1,0 Q1,1 · · · Q1,s−1 ].
By considering the last block column of equation (44),
under the assumption that O2R+1,s−1 has full column rank, we
can obtain that Q1,s−1 = 0. Similarly, by considering the second to the last block column, it can be derived that Q1,s−2 = 0

13

and ∆0,s−2 = 0, further ∆0,j = 0 for j = 0, 1, · · · , s − 3.
Next, for the third to the last block column, since ∆0,s−3 = 0
and O2R+1,s−1 has full column rank, we can obtain that
Q1,s−3 = 0 and ∆1,j = 0 for j = 0, · · · , s − 3. By repeating
this procedure, it can be established that equation (44) holds
only when ∆1 = 0 and Q1,j = 0 for j = 0, 1, · · · , s − 1. In
other words, the matrix ∆s with a nontrivial sub-matrix
∆1

DR
cannot be linearly represented by Os Ts
.
Since the block rows of ∆1 correspond to the secondlayer Toeplitz part of ∆s , we can see that the rank
R
Os TsD
∆s
is strictly larger than that of
of


DR
Os Ts
as long as the second-layer Toeplitz part of
∆s is non-zero. Thus, the proof is completed.

To illustrate this, the first block-column equation of (47) can
be equivalently written as follows:











Ĉ
ˆl
CA
0
0
0
0
0
0




C
  CAl
 
  0
  0
=
  0
 
  0

0
0

0
CA
C
0
0
0
CAl
0

0
CAr
0
C
0
0
CA
CAl

0
0
0
0
C
0
CAr
CA

0
0
0
0
0
C
0
CAr












Q1,1
Q2,1 
Q3,1 
.
Q4,1
Q5,1



It is obvious that, since O4,2 (or O2R,s/2 ) has full column
rank, it has

A PPENDIX C
P ROOF OF T HEOREM 2
Q2,1 = Q3,1 = Q4,1 = Q5,1 = 0,
ˆ l = CAl Q1,1 .
Ĉ = CQ1,1 , CA

Assume that the matrix pair {Ô2R+1,s/2 , Ĉ2R+1,s/2 } is one
of the optimal solution pairs to (20), it can be established that
Ô2R+1,s/2 and O2R+1,s/2 has the same structure and the same
column space, which can be algebraically represented as
Ô2R+1,s/2 = O2R+1,s/2 Q,

(45)

where Q ∈ R(2R+1)n×(2R+1)n is a (2R + 1, 2R + 1) block
matrix which is nonsingular. In the sequel, we denote Qi,j the
(i, j)-th block entry of Q. In order to prove the theorem, it is
sufficient to prove that the estimate Ŵs/2 of Ws/2 satisfies
that
Ŵs/2 = Ws/2 Q,
(46)
where Q ∈ Rn×n is a nonsingular matrix and Ws/2 is defined
in (19). In the sequel, we denote by Ŵi,j , for j = −i, · · · , i,
the estimate of Wi,j which is a block component of Ws/2 .
To illustrate the structure of equation (45), it is expanded
into equation (47) by setting R = 2 and s = 4.
The following proof will be divided into three steps.
Step 1. We will show that, under the assumption that
O2R,s/2 has full column rank, the following relations can be
derived from equation (45):
i = 0 : Ŵi,−i = Wi,−i Q1,1 ,
Wi,−i Q1,j = 0 for j = 2, · · · , s − 1;
i = 1 : Ŵi,−i = Wi,−i Q1,1 ,
Ŵi,2−i = Wi,−i Q1,3 + Wi,1−i Q1,2 + Wi,2−i Q1,1 ,
Wi,−i Q1,j+2 + Wi,1−i Q1,j+1 + Wi,2−i Q1,j = 0
for j = 2, · · · , s − 3;
···

Qi,2 = 0 for i = 3, · · · , 2R + 1;

CQ1,2 = 0;

Ŵi,1−i = Wi,−i Q1,2 + Wi,1−i Q1,1 for i = 1, · · · , s/2 − 1.
(51)
To illustrate this, by using the relations in (50), the second
block-column equation of (47) can be equivalently written as












ˆ
CA
0
CQ1,1
0
0
0
CAl Q1,1
0



 CA
l
  C
 
 
 
 
=
 
 
 


CA

CAr



C
C
C
C
CAl

CA
CAl

CAr
CA

CAr











Q1,2
Q2,2
Q3,2
Q4,2
Q5,2

We can observe that the bottom subpart of the vector on the
left-hand side equals the product of the first block column of
O4,2 and Q1,1 . By the assumption that O4,2 (or O2R,s/2 ) has
full column rank, we can derive that
Q3,2 = Q4,2 = Q5,2 = 0; CQ1,2 = 0;
ˆ = CAl Q1,2 + CAQ1,1 .
CA

i = s/2 − 1 : Ŵi,−i = Wi,−i Q1,1 ,
···
Ŵi,s−2−i = Wi,−i Q1,s−1 + · · · + Wi,s−2−i Q1,1 .
(48)

Let’s consider the equation corresponding to the first blockcolumn of (45). Since O2R,s/2 has full column rank, it can be
derived that
Ŵi,−i = Wi,−i Q1,1 for i = 0, 1, · · · , s/2 − 1.

Next, by substituting the relations in (49) into the equation
corresponding to the second block-column equation of (45)
and under the assumption that O2R,s/2 has full column rank,
we can obtain that

(52)

Ŵi,1−i = Wi,−i Q1,2 + Wi,1−i Q1,1 ,

Qi,1 = 0 for i = 2, 3, · · · , 2R + 1;

(50)

(49)

Q2,2 = Q1,1 ;

(53)

By iteratively considering from the first to the (s − 1)-th
block-column equation of (45) using the same strategy, we
can derive the relations in equation (48).
Step 2. Under the assumption that O2R,s/2 has full column
rank, the following relations can be derived from equation (45)





.


14





Ĉ









 ˆ
 CAl



Ĉ
Ĉ
Ĉ
Ĉ
ˆ
CA
ˆl
CA

ˆr
CA
ˆ
CA
ˆl
CA

ˆr
CA
ˆ
CA

ˆr
CA

 C

 
 
 
 
 
=
 
  CA
l
 





C
C
C
C
CA
CAl

CAr
CA
CAl

as well:

CAr
CA

CAr











Q1,1
Q2,1
Q3,1
Q4,1
Q5,1

Q1,2
Q2,2
Q3,2
Q4,2
Q5,2

Q1,3
Q2,3
Q3,3
Q4,3
Q5,3

Q1,4
Q2,4
Q3,4
Q4,4
Q5,4

Q1,5
Q2,5
Q3,5
Q4,5
Q5,5







(47)

equations (55)-(56) can be explicitly written as

i = 0 : Ŵi,i = Wi,i Q2R+1,2R+1 ,
Wi,i Q2R+1,j = 0 for j = 2R + 3 − s, · · · , 2R;
i = 1 : Ŵi,i = Wi,i Q2R+1,2R+1 ,
Ŵi,i−1 = Wi,i Q2R+1,2R + Wi,i−1 Q2R+1,2R+1 ,
Ŵi,i−2 = Wi,i Q2R+1,2R−1 + Wi,i−1 Q2R+1,2R
+ Wi,i−2 Q2R+1,2R+1 ,
Wi,i Q2R+1,j + Wi,i−1 Q2R+1,j+1
+ Wi,i−2 Q2R+1,j+2 = 0
for j = 2R + 3 − s, · · · , 2R − 2;
···
i = s/2 − 1 : Ŵi,i = Wi,i Q2R+1,2R+1 ,
···
Ŵi,i−s+2 = Wi,i Q2R+1,2R+3−s + · · ·
+ Wi,i−s+2 Q2R+1,2R+1 .
(54)

The relations in (54) are derived using the same strategy in
Step 1 by iteratively considering from the (2R + 1)-th to the
(2R + 3 − s)-th block-column equation of (45) in a reverse
order.
Step 3. We will show Ŵs/2 = Ws/2 Q1,1 using the relations
in (48) and (54).
Denote by Ô2s−3,s/2 (:, s−1) the (s−1)-th block column of
the estimate Ô2s−3,s/2 , which is a sub-matrix of Ô2R+1,s/2 .
The relations in (48) can be compactly written as:


Q1,s−1


..


.


 Q1,1 
,
Ô2s−3,s/2 (:, s − 1) = O2s−3,s/2 
(55)


0




..


.
0

while those in (54) can be compactly written as:

0

..

.


0
Ô2s−3,s/2 (:, s − 1) = O2s−3,s/2 
 Q2R+1,2R+1


..

.

Q2R+1,2R+3−s







 (56)





To illustrate this using the specific example shown in (47), the














0
0
Ĉ
0
0
ˆr
CA
ˆ
CA
ˆl
CA








Q1,3


 Q1,2

 = O5,2 
 Q1,1

 0



0






,














0
0
Ĉ
0
0
ˆr
CA
ˆ
CA
ˆl
CA






0


 0

 = O5,2 
 Q5,5

 Q

5,4


Q5,3






.

(57)

In view of equations (55) and (56), we can obtain that




Q1,s−1
0




..
..




.
.




 Q1,2 


0







O2s−3,s/2  Q1,1  = O2s−3,s/2  Q2R+1,2R+1 





0
Q2R+1,2R








..
..




.
.
0

Q2R+1,2R+3−s

By assumption 3) of the theorem, the matrix O2s−3,s/2 has
full column rank. As a result, it can be established from the
above equation that
Q1,1 = Q2R+1,2R+1 ;

Q1,2 = · · · = Q1,s−1 = 0;

Q2R+1,2R+3−s = · · · = Q2R+1,2R = 0.

(58)

Since the block vector Ô2s−3,s/2 (:, s − 1) contains all the
block components of Ŵs/2 , we can obtain that
Ŵs/2 = Ws/2 Q1,1 .
Therefore, the proof is completed.
R EFERENCES
[1] B. Bamieh, F. Paganini, and M. Dahleh, “Distributed control of spatially
invariant systems,” Automatic Control, IEEE Transactions on, vol. 47,
pp. 1091–1107, Jul 2002.
[2] N. Motee and A. Jadbabaie, “Optimal control of spatially distributed
systems,” Automatic Control, IEEE Transactions on, vol. 53, pp. 1616–
1629, Aug 2008.
[3] J. Rice and M. Verhaegen, “Distributed control in multiple dimensions:
A structure preserving computational technique,” Automatic Control,
IEEE Transactions on, vol. 56, pp. 516–530, March 2011.
[4] A. Mahajan, N. C. Martins, M. Rotkowitz, and S. Yüksel, “Information
structures in optimal decentralized control.,” in CDC, pp. 1291–1306,
2012.
[5] M. Wozny and W. Carpenter, “A characteristics approach to parameter
identification in distributed systems,” Automatic Control, IEEE Transactions on, vol. 14, pp. 422–422, Aug 1969.
[6] D. Gorinevsky and M. Heaven, “Performance-optimized applied identification of separable distributed-parameter processes,” Automatic Control,
IEEE Transactions on, vol. 46, pp. 1584–1589, Oct 2001.

15

[7] J. Rice and M. Verhaegen, “Efficient system identification of heterogeneous distributed systems via a structure exploiting extended kalman
filter,” Automatic Control, IEEE Transactions on, vol. 56, pp. 1713–
1718, July 2011.
[8] T. Bewley and S. Liu, “Optimal and robust control and estimation of
linear paths to transition,” Journal of Fluid Mechanics, vol. 365, pp. 305–
349, Jun 25 1998.
[9] F. Wu and S. Yildizoglu, “Distributed parameter-dependent modeling
and control of flexible structures,” Journal of Dynamic Systems Measurement and Control-Transactions of The ASME, vol. 127, pp. 230–239,
Jun 2005.
[10] D. G. MacMynowski, R. Heimsten, and T. Andersen, “Distributed Force
Control of Deformable Mirrors,” European Journal of Control, vol. 17,
pp. 249–260, May-Jun 2011.
[11] R. R. Pena, R. D. Fernandez, and R. J. Mantz, “Passivity control via
Power Shaping of a wind turbine in a dispersed network,” International
Journal of Hydrogen Energy, vol. 39, pp. 8846–8851, May 2014.
[12] A. Haber and M. Verhaegen, “Subspace identification of large-scale
interconnected systems,” Automatic Control, IEEE Transactions on,
vol. 59, pp. 2754–2759, Oct 2014.
[13] J. Kulkarni, R. DAndrea, and B. Brandl, “Application of distributed
control techniques to the adaptive secondary mirror of cornells large atacama telescope,” in SPIE Astronomical Telescopes and Instrumentation
Conf., pp. 750–756, 2002.
[14] P. Massioni and M. Verhaegen, “Subspace identification of circulant
systems,” Automatica, vol. 44, no. 11, pp. 2825 – 2833, 2008.
[15] P. Massioni and M. Verhaegen, “Subspace identification of distributed,
decomposable systems,” in Decision and Control, 2009 held jointly
with the 2009 28th Chinese Control Conference. CDC/CCC 2009.
Proceedings of the 48th IEEE Conference on, pp. 3364–3369, Dec 2009.
[16] N. Matni and A. Rantzer, “Low-rank and low-order decompositions for
local system identification,” arXiv:1403.7175, 2014.
[17] M. Verhaegen, “Subspace techniques in system identification,” Encyclopedia of Systems and Control, pp. 1386–1396, 2015.
[18] G. Mercere, O. Prot, and J. Ramos, “Identification of parameterized
gray-box state-space systems: From a black-box linear time-invariant
representation to a structured one,” Automatic Control, IEEE Transactions on, vol. 59, pp. 2873–2885, Nov 2014.
[19] L. Ljung, System Identification: Theory for the user (Second Edition).
Prentice Hall PTR, 1999.
[20] M. Verhaegen and V. Verdult, Filtering and System Identification: A
Least Squares Approach. Cambridge University Press, 2007.
[21] Z. Liu, A. Hansson, and L. Vandenberghe, “Nuclear norm system
identification with missing inputs and outputs,” Systems and Control
Letters, vol. 62, no. 8, pp. 605 – 612, 2013.
[22] M. Verhaegen and A. Hansson, “Nuclear norm subspace identification
(n2sid) for short data batches,” CoRR abs/1401.4273, 2014.
[23] C. Yu and M. Verhaegen, “Subspace identification of distributed clusters
of homogeneous systems,” IEEE Transactions on Automatic Control,
vol. 62, no. 1, pp. 463–468, 2017.
[24] P. Dewilde and A.-J. Van der Veen, Time-varying systems and computations. Springer Science & Business Media, 2013.
[25] K. Mohan and M. Fazel, “Reweighted nuclear norm minimization with
application to system identification,” in American Control Conference
(ACC), 2010, pp. 2953–2959, IEEE, 2010.
[26] L.-L. Xie and L. Ljung, “Estimate physical parameters by black box
modeling,” in Proceedings of the 21st Chinese Control Conference,
pp. 673–677, 2002.
[27] C. Yu, M. Verhaegen, S. Kovalsky, and R. Basri, “Identification of
structured lti mimo state-space models,” in Decision and Control (CDC),
2015 IEEE 54th Annual Conference on, pp. 2737–2742, IEEE, 2015.
[28] G. Strang, “Fast transforms: Banded matrices with banded inverses,”
Proceedings of the National Academy of Sciences of the United States
of America, vol. 107, no. 28, pp. 12413 – 12416, 2009.
[29] S. Choudhary and U. Mitra, “Identifiability scaling laws in bilinear
inverse problems,” arXiv preprint arXiv:1402.2637, 2014.

