Spatio-temporal Transfer Function Conditions of Positive Realness for
Translation Invariant Lattice Networks of Interacting Linear Systems∗

arXiv:1802.00185v1 [cs.SY] 1 Feb 2018

Igor G. Vladimirov† ,

Abstract— This paper is concerned with networks of interacting linear systems at sites of a multidimensional lattice. The
systems are governed by linear ODEs with constant coefficients
driven by external inputs, and their internal dynamics and
coupling with the other component systems are translation
invariant. Such systems occur, for example, in finite-difference
models of large-scale flexible structures manufactured from
homogeneous materials. Using the spatio-temporal transfer
function of this translation invariant network, we establish conditions for its positive realness in the sense of energy dissipation.
The latter is formulated in terms of block Toeplitz bilinear
forms of the input and output variables of the composite
system. We also discuss quadratic stability of the network in
isolation from the environment and phonon theoretic dispersion
relations.
Index Terms— Translation invariant networks, spatiotemporal frequency domain, energy balance relations.
MSC Codes — 93C05, 90B10, 37L60, 37K05, 93C80.

I. INTRODUCTION
Complex physical systems can be viewed as a large
number of relatively simple subsystems whose collective
behaviour is a cumulative effect of their interaction rather
than a particular individual structure. Spatially homogeneous
states of matter are modelled as identical building blocks
which interact with each other in a translationally invariant
fashion. A natural example of large-scale composite systems
with translational symmetry is provided by crystalline solids,
where spatially periodic arrangements of constituent particles
result from their interaction and play an important role in
their thermodynamic and mechanical properties (including
the heat transfer and wave propagation) studied in the phonon
theory [10].
Modern engineering exploits translation invariant interconnections in drone swarming, vehicle platooning and artificially fabricated metamaterials [12], such as split ring resonator arrays with unusual electrodynamic characteristics (a
negative refraction index). Nontrivial input-output properties
of such networks of systems (natural or artificial) are not
merely a “sum” of individual internal dynamics of their
constituent blocks and come from a specific structure of
energy flows through the translation invariant interaction.
The energy balance relations, which reflect the energy
conservation and dissipation in isolated and open systems
∗ This work is supported by the Australian Research Council under grant
DP160101121.
† Research School of Engineering, College of Engineering and Computer
Science, Australian National University, ACT 2601, Canberra, Australia.
igor.g.vladimirov@gmail.com, i.r.petersen@gmail.com.

Ian R. Petersen†

(for example, due to electrical resistance and mechanical friction), significantly affect the behaviour of physical systems
and play an increasingly important role in control design
[4], [5], [11]. These equations involve the internal energy
and the work done on the system (which are represented in
the dissipativity theory [14] in terms of storage and supply
rate functions). Work is modelled by using a bilinear form of
the input and output variables, which are interpreted as the
generalised force and velocity respectively. For linear timeinvariant systems with a finite-dimensional internal state,
the properties of being passive, positive real or negative
imaginary (in the case of position variables instead of the
velocity as the output) admit criteria in the form of linear
matrix inequalities for the transfer functions in the frequency
domain or the state-space matrices themselves [7], [15].
The present paper is concerned with similar conditions
for networks of interacting linear systems at sites of a
multidimensional lattice. The composite system is governed
by an infinite set of linear ODEs with constant coefficients
driven by external inputs, and their internal dynamics and
coupling with the other component systems are translation
invariant. These ODEs have block Toeplitz state-space matrices and can be represented in the spatio-temporal frequency
domain by using appropriately modified transfer functions of
several variables. Such systems arise, for example, as finitedifference approximations of PDEs for large-scale flexible
structures made of spatially homogeneous materials. Using the spatio-temporal transfer function of this translation
invariant network, we establish conditions for its positive
realness in the sense of energy dissipation. The latter is
formulated in terms of block Toeplitz bilinear forms of the
input and output variables of the composite system. The
multivariate Laplace and Fourier transform techniques, which
are used for this purpose, are similar to those for distributed
control systems in the classical and quantum settings [9],
[13].
The paper is organised as follows. Section II describes the
class of translation invariant networks under consideration.
Section III represents the network dynamics in terms of the
spatial Fourier transforms of its signals. Section IV discusses
energy balance relations in the case of bilinear supply
rate and quadratic storage functions. Section V establishes
conditions for passivity of the network in terms of its spatiotemporal transfer function, and also discusses quadratic stability bounds for a dissipative network in isolation from
the environment. Section VI considers phonon theoretic
dispersion relations for the isolated network. Section VII

of the lattice forming an exhausting sequence. Since the
matrix-valued functions A, B, C, D (with values in Cn×n ,
Cn×m , Cr×n , Cr×m , respectively) are 2π -periodic with respect
to each of their arguments, then, without loss of generality,
they can be considered on a ν -dimensional torus Tν (where
T can be identified with the interval [−π , π )). These functions are Hermitian in the sense that A(σ ) = A(−σ ) (and
similarly for the other functions B, C, D) for all σ ∈ Rν ,
where (·) is the complexpconjugate. Furthermore,
p
plet kuk2 :=
∑k∈Zν |uk |2 , kxk2 := ∑k∈Zν |xk |2 , kyk2 := ∑k∈Zν |yk |2
denote the norms of the input, state and output of the network
at a fixed but otherwise arbitrary moment of time t > 0 in
the corresponding Hilbert spaces

provides concluding remarks.
II. TRANSLATION INVARIANT NETWORKS
We consider a translation invariant network of coupled
linear systems at sites of the ν -dimensional integer lattice
Zν . For any spatial index k ∈ Zν , the kth system is endowed
with an Rn -valued vector xk of time-varying state variables
(for example, the generalised positions and velocities). Associated with the kth lattice site are vectors uk and yk of
external input and output variables which take values in Rm
and Rr , respectively, and also vary in time (unless indicated
otherwise, vectors are organised as columns). The states and
outputs of these systems are governed by an infinite set of
coupled ODEs
ẋ j =

∑ν (A j−k xk + B j−k uk ),

(1)

∑ν (C j−k xk + D j−k uk )

(2)

U := L2 (Zν , Rm ), X := L2 (Zν , Rn ), Y := L2 (Zν , Rr ) (9)

k∈Z

yj =

of appropriately dimensioned square summable real vectorvalued functions f := ( fk )k∈Zν and g := (gk )k∈Zν on the
lattice Zν with the inner product h f , gi := ∑k∈Zν fkT gk . Now,
suppose the state-space matrices in (3) and (4) specify
bounded linear operators A : X → X, B : U → X, C : X → Y,
D : U → Y. This is equivalent to their L2 -induced operator
norms being finite:

k∈Z

for all j ∈
where (˙) is the time derivative (the time
arguments are often omitted for brevity). Their right-hand
sides are organised as convolutions, with the matrices Aℓ ∈
Rn×n , Bℓ ∈ Rn×m , Cℓ ∈ Rr×n , Dℓ ∈ Rr×m depending on the
relative location ℓ ∈ Zν of the lattice sites in accordance
with the translation invariance of the individual dynamics of
the systems and their coupling. By assembling the inputs,
states and outputs into the infinite-dimensional vectors u :=
(uk )k∈Zν , x := (xk )k∈Zν , y := (yk )k∈Zν , the set of ODEs (1)
and (2) can be written as
Zν ,

ẋ = Ax + Bu,

(3)

y = Cx + Du,

(4)

Z t

(Ax(τ ) + Bu(τ ))dτ

e(t−τ )A Bu(τ )dτ ,

kCk∞ = ess sup kC(σ )k < +∞,

(12)

kDk∞ = ess sup kD(σ )k < +∞,

(13)

σ ∈Tν

where the essential supremum is applied to the operator
norms of the matrices in (8). In particular, if the matrices Aℓ ,
Bℓ , Cℓ , Dℓ vanish for all ℓ ∈ Zν with |ℓ| large enough (so that
each of the component systems in (1) and (2) is coupled with
a finite number of the other systems in the network), then
the functions A, B, C, D in (8) are multivariate trigonometric
polynomials, and the conditions (10)–(13) are satisfied in this
case. In general, the fulfillment of (10) and (11) guarantees
that the operator exponential etA and the product e(t−τ )A B
in (6) are bounded block Toeplitz operators acting on the
Hilbert space X and from U to X, respectively, with

(5)

t > 0.

(11)

σ ∈Tν

whose solution is given by
Zt

kBk∞ = ess sup kB(σ )k < +∞,
σ ∈Tν

0

x(t) = etA x(0) +

(10)

σ ∈Tν

where the matrices A := (A j−k ) j,k∈Zν , B := (B j−k ) j,k∈Zν ,
C := (C j−k ) j,k∈Zν , D := (D j−k ) j,k∈Zν are block Toeplitz in
the sense of the additive group structure of the lattice Zν .
The ODE (3) is understood as a Volterra integral equation
(of the second kind)
x(t) =

kAk∞ = ess sup kA(σ )k < +∞,

(6)

ketA k∞ 6 ekAk∞t ,

0

For completeness, we note that (3)–(6) can be endowed with
a rigorous meaning as follows. To this end, the state-space
matrices in (1) and (2) are assumed to be square summable:
 2

Aℓ Bℓ
(7)
∑ Cℓ Dℓ < +∞,
ℓ∈Zν

ke(t−τ )A Bk∞ 6 ekAk∞ (t−τ ) kBk∞ (14)

for any t > τ > 0, in view of the submultiplicativity of the
operator norm and the fact that block Toeplitz matrices form
an algebra. Therefore, if the initial network state is square
summable, that is,
kx(0)k2 < +∞,

where k · k is an arbitrary matrix norm (whose particular
choice is irrelevant in this case). The fulfillment of (7) allows
the Fourier series




T
A(σ ) B(σ )
A Bℓ
:= ∑ e−iℓ σ ℓ
,
(8)
C(σ ) D(σ )
Cℓ Dℓ
ℓ∈Zν

(15)

and the network input is locally absolutely integrable with
respect to time in the sense that
Z T
0

to be defined (in a blockwise fashion) for almost all σ ∈ Rν
as the L2 -limit of appropriate partial sums over finite subsets

ku(t)k2 dt < +∞,

T > 0,

(16)

then these properties are inherited by the subsequent states
of the network. Indeed, a combination of (6) with (14)–(16)
2

where

leads to
Z t


kx(t)k2 6 ekAk∞t kx(0)k2 + kBk∞ ku(τ )k2 dτ < +∞
0
(17)
for all t > 0, and hence,
Z T

kx(t)k2 dt 6T ekAk∞T kx(0)k2
0
Z T

ku(τ )k2 dτ < +∞
+ kBk∞
(18)

U(s, σ ) :=
Here,

F(s, σ ) := C(σ )(sIn − A(σ ))−1 B(σ ) + D(σ )

0

Z T
0

+ kDk∞

0

ku(t)k2 dt < +∞

(19)

for any T > 0. In Section IV, we will replace (16) with a
stronger condition on the network inputs in order to guarantee time-local square integrability for the outputs instead of
(19).
III. NETWORK DYNAMICS IN THE SPATIAL
FREQUENCY DOMAIN
The preservation of the spatial square summability (17)
(provided the input satisfies (16)) allows the network dynamics (3) and (4) to be represented in the spatial frequency
domain as

∂t X(t, σ ) = A(σ )X(t, σ ) + B(σ )U(t, σ ),
Y (t, σ ) = C(σ )X(t, σ ) + D(σ )U(t, σ ).

xk (t) = Re(est+ik

(20)

T

∑ν e−iℓ σ uℓ(t),

(22)

ℓ∈Z

X(t, σ ) :=

T

∑ν e−iℓ σ xℓ(t),

(23)

ℓ∈Z

Y (t, σ ) :=

T

∑ν e−iℓ σ yℓ(t)

(24)

ℓ∈Z

for almost all σ ∈ Tν at time t > 0. For any given σ ,
the equations (20) and (21) describe an autonomous system
(which is independent of the other systems in this parametric
family with different values of σ ) with a finite-dimensional
internal state X(·, σ ). Their solution can be represented in
terms of the Laplace transform over the time variable as
X(s, σ ) :=

Z +∞
0

Y(s, σ ) :=

Z +∞
0

e

−st

k ∈ Zν , t > 0,

z),

(30)

IV. SUPPLY RATE AND ENERGY BALANCE
A generalised model for the work done by the input u on
the network is provided by a supply rate [14] at time t > 0
in the form

e−st X(t, σ )dt

=(sIn − A(σ ))−1 (X(0, σ ) + B(σ )U(s, σ )),

Tσ

satisfy the ODEs (1) with u = 0. Indeed, substitution of the
latter equality and (30) into the right-hand side of (1) yields

T
T
∑k∈Zν A j−k xk (t) = Re est+i j σ ∑k∈Zν e−i( j−k) σ A j−k z =


T
T
Re est+i j σ A(σ )z = Re sest+i j σ z = ẋ j (t). For purely
imaginary eigenvalues s = iω of the matrix A(σ ) (with
ω ∈ R), the solutions (30) describe persistent oscillations
of frequency ω in the network, which are organised as
(hyper) plane waves of length |2σπ| in Rν . Their wavefronts
are orthogonal to the vector σ and move (along σ ) at
constant phase velocity − |ωσ | , with the direction depending
on the sign of ω . Similarly to the dispersion relations of
the phonon theory, the spectral structure of such oscillations
in the network is represented by the multi-valued map
Tν ∋ σ 7→ ω ∈ R (of the wave vector to the frequency),
which will be discussed in Section VI.

(21)

Here, similarly to (8), the functions U, X, Y on R+ × Tν
with values in Cm , Cn , Cr are the Fourier transforms of the
input, state and output of the network:
U(t, σ ) :=

(28)

where r(·) denotes the spectral radius of a square matrix
(so that ln r(eN ) = max ReS(N) is the largest real part of
the eigenvalues of N, with S(N) denoting its spectrum).
The presence of the upper limit in (29) makes the integral
m
2 ν
in (27) convergent in the Hilbert space
pR L (T , C ) due
2
to the Parseval identity kU(t, ·)k2 :=
Tν |U(t, σ )| dσ =
ν
/2
(2π ) ku(t)k2 for any t > 0. In the case when the input vanishes, so that the network is effectively isolated
from the environment, (20) reduces to a homogeneous linear ODE ∂t X(t, σ ) = A(σ )X(t, σ ). Its solution X(t, σ ) =
et A(σ ) X(0, σ ), considered for different values of σ ∈ Tν ,
admits a direct analogy (formulated in system theoretic
terms) with the phonon theory [10] of collective oscillations
in spatially periodic arrangements of atoms in crystalline
solids. More precisely, for a given σ ∈ Tν , let z ∈ Cn be an
eigenvector of the matrix A(σ ) associated with its eigenvalue
s ∈ C. Then the functions

kx(t)k2 dt

Z T

(27)

is the spatio-temporal transfer function of the network with
values in Cr×m . In view of (10), the relations (25)–(28) are
valid for all s ∈ C satisfying

ln ku(t)k2 
, (29)
Res > max ess sup ln r(eA(σ ) ), lim sup
t
t→+∞
σ ∈Tν

0

ky(t)k2dt 6kCk∞

e−st U(t, σ )dt.

0

for any T > 0. Together with (4), (12), (13), the property (17)
implies that the output norm ky(t)k2 is finite for almost all
t > 0 since so is ku(t)k2 in view of (16). Furthermore, due to
(18), the network output is also locally absolutely integrable
with respect to time:
Z T

Z +∞

S(t) := hu(t), Gy(t)i =

(25)

∑ ν u j (t)TG j−k yk (t)

j,k∈Z

1
=
U(t, σ )∗ G(σ )Y (t, σ )dσ
(2π )ν Tν
1
=
hU(t, ·), GY (t, ·)i,
(2π )ν
Z

Y (t, σ )dt

=C(σ )X(s, σ ) + D(σ )U(s, σ )
=C(σ )(sIn − A(σ ))−1 X(0, σ ) + F(s, σ )U(s, σ ), (26)
3

(31)

network state at time t > 0:

where the Parseval identity is applied to the Fourier transforms U and Y from (22) and (24). The quantity S(t) is
a bilinear function of the current input and output of the
network which is specified by a block Toeplitz matrix G :=
(G j−k ) j,k∈Zν , where Gℓ ∈ Rm×r are given matrices satisfying

∑

ℓ∈Zν

2

kGℓ k < +∞,

1
1
H(t) := hx(t),V x(t)i =
2
2
=

(32)

T

∑ν e−iℓ σ Gℓ.

(33)

It is assumed that the matrix G describes a bounded linear
operator from the output space Y in (9) to the input space
U in (9), so that
(34)

similarly to the operator norms in (10)–(13). For example,
the condition (34) holds in the case of a standard supply
rate, when u and y consist of the corresponding force and
velocity variables (with r = m) and G is the identity operator.
Returning to the general case, substitution of (21) into (31)
leads to
1
Re U(t, σ )∗ G(σ )
(2π )ν Tν

× (C(σ )X(t, σ ) + D(σ )U(t, σ )) dσ
Z
1
=
Z(t, σ )∗
2(2π )ν Tν


0
C(σ )∗ G(σ )∗
×
G(σ )C(σ ) G(σ )D(σ ) + D(σ )∗ G(σ )∗




T
x (t)
X(t, σ )
= ∑ e−ik σ k
uk (t)
U(t, σ )
k∈Zν



T

∑ν e−iℓ σ Vℓ

(35)

(36)

N(−σ ) = N(σ )T = N(σ ),

σ ∈ Tν ,

(46)

similarly to (38). In order for the network model to correspond to a real physical system, the dissipation rate in
(44) has to be nonnegative in view of the total energy
conservation, and hence,

takes values in the subspace of complex Hermitian matrices
of order n and satisfies
(38)

S(t) > Ḣ(t).

Under the additional condition
kV k∞ = ess sup kV(σ )k < +∞,

(42)

of order n + m in (44) is obtained by using (35), (42) and
(43), and satisfies

(37)

σ ∈ Tν .

∂t H(t, σ )dσ ,

can be interpreted as the energy dissipation rate. This quantity is part of the supply rate in (31) which is not accounted
for by the rate of change of the internal energy in (42). The
complex Hermitian matrix

−A(σ )∗ V(σ ) − V(σ )A(σ )
N(σ ) :=
G(σ )C(σ ) − B(σ )∗ V(σ )

C(σ )∗ G(σ )∗ − V(σ )B(σ )
(45)
G(σ )D(σ ) + D(σ )∗ G(σ )∗

ℓ∈Z

V(−σ ) = V(σ )T = V(σ ),

Tν

where use is made of (36). If H(t) in (40) describes the
internal energy (usually referred to as the Hamiltonian) of
the network at the current moment of time t > 0, then the
difference
Z
1
Z(t, σ )∗ N(σ )Z(t, σ )dσ
(44)
S(t) − Ḣ(t) =
2(2π )ν Tν

is the Fourier transform of the augmented state-input pair
of the network at time t > 0. Now, let V := (V j−k ) j,k∈Zν be
a symmetric block Toeplitz matrix, which is specified by a
square summable matrix-valued function Zν ∋ ℓ 7→ Vℓ ∈ Rn×n
on the lattice satisfying V−ℓ = VℓT . Then the corresponding
Fourier transform
V(σ ) :=

Z


∂t H(t, σ ) = Re X(t, σ )∗ V(σ )∂t X(t, σ )

= Re X(t, σ )∗ V(σ )(A(σ )X(t, σ ) + B(σ )U(t, σ ))


1
A(σ )∗ V(σ ) + V(σ )A(σ ) V(σ )B(σ )
Z(t, σ ),
= Z(t, σ )∗
0
B(σ )∗ V(σ )
2
(43)

where
Z(t, σ ) :=

(40)

with

Z

× Z(t, σ )dσ ,

H(t, σ )dσ ,

1
Ḣ(t) =
(2π )ν

σ ∈Tν

S(t) =

Tν

1
(41)
H(t, σ ) := X(t, σ )∗ V(σ )X(t, σ ).
2
Here, the Parseval identity is applied to the Fourier transform
X from (23), whose contribution to H(t) at a given σ ∈ Tν
is quantified by the real-valued quantity H(t, σ ) in (41). The
spatial frequency domain representation (20) of the network
state dynamics allows the time derivative of H in (40) to be
computed as

ℓ∈Z

kGk∞ = ess sup kG(σ )k < +∞,

Z

x j (t)TV j−k xk (t)

where

with the Fourier transform
G(σ ) :=

1
(2π )ν

∑

j,k∈Zν

(47)

A sufficient condition for this inequality to hold for arbitrary
admissible inputs u and initial network states x(0) is positive
semi-definiteness of the matrix (45):

(39)

σ ∈Tν

the matrix V describes a self-adjoint operator on the Hilbert
space X in (9). This gives rise to a quadratic form of the

N(σ ) < 0,
4

almost all σ ∈ Tν .

(48)

The condition (54) is stronger than (16) and guarantees
finiteness of the work by such an input on the network over
any bounded time interval. More precisely, (31) implies that
Z T
Z T
2 Z T
|S(t)|dt 6
kGy(t)k22dt
ku(t)k22dt

This condition is also necessary under an additional controllability assumption.
Theorem 1: Suppose the matrix pair (A(σ ), B(σ )) for the
network in (1) and (2), given by (8), is controllable for any
σ ∈ Tν . Then the dissipation rate in (44) is nonnegative
for any square summable initial network states and locally
integrable inputs in (16) if and only if the matrix (45) satisfies
(48).

Proof: In view of (44), the fulfillment of (48) implies
(47) regardless of the controllability of (A, B), thus proving
the sufficiency of (48), mentioned above. In order to show
the necessity, let ζ : Tν → Cn+m be a measurable function
which is Hermitian (that is, ζ (σ ) = ζ (−σ ) for any σ ∈ Tν )
and such that ζ (σ ) is a unit eigenvector associated with the
smallest eigenvalue

λ (σ ) := λmin (N(σ ))

0

6

1
2

(49)

Tν

min(λ (σ ), 0)dσ .

0

ky(t)k22 dt 6
Z T
0

(52)

0

ku(t)k22dt =

Z T

∑

0 k∈Zν

|uk (t)|2 dt < +∞,

0

ky(t)k22 dt

(55)

0

(kCx(t)k22 + kDu(t)k22)dt
Z T
0

ku(t)k22 dt,

(56)

V. NETWORK PASSIVITY CONDITIONS IN THE
SPATIAL FREQUENCY DOMAIN
The relations (53) imply that the network is passive in the
sense that any time-locally square integrable input u in (54)
performs a nonnegative work:
Z T

T > 0.

S(t)dt > 0,

(58)

0

The passivity can also be considered irrespective of a specific internal energy (storage) function (provided the initial
network state is zero). For what follows, the network input
u is assumed to be time-space square summable in the sense
that
sZ
rZ
+∞

(53)

|||u||| :=

0

ku(t)k22dt =

+∞

0

∑

k∈Zν

|uk (t)|2 dt < +∞,

(59)
where theR norm ||| · ||| is associated
with
the
inner
product
R
hh f , gii := 0+∞ h f (t), g(t)i = 0+∞ ∑k∈Zν fk (t)∗ gk (t)dt for real
or complex vector-valued functions f (t) := ( fk (t))k∈Zν and
g(t) := (gk (t))k∈Zν of the time Rand space variables. By
the Parseval identity, (59) yields R×Tν |U(iω , σ )|2 dω dσ =

for any T > 0. These relations hold for any network input u
such that ku(t)k2 is a locally square integrable function of
time t > 0 in the sense that
Z T

Z T

where the Cauchy-Bunyakovsky-Schwarz inequality is used
again. The local square integrability of kx(t)k2 as a function
of time t > 0 follows from (57) in view of (15) and (54).
In combination with (56), this ensures
the time-local square
R
integrability of the network output: 0T ky(t)k22 dt < +∞ for
any T > 0. From the last property, (54) and (55), it follows
that the supply rate S(t) is indeed locally integrable with
respect to time.

If the self-adjoint operator V is positive semi-definite
(which is equivalent to that the corresponding Fourier transform V in (37) satisfies V(σ ) < 0 for almost all σ ∈ Tν ),
and the network is initialized at zero state x(0) = 0, then the
integration of both parts of (47) over a time interval [0, T ]
leads to
0

Z T

kx(t)k22 dt + kDk2∞

W (T ) :=

S(t)dt > H(T ) − H(0) = H(T ) > 0

0

ku(t)k22dt

0

If the property (48) does not hold, then the set {σ ∈ Tν :
λ (σ ) < 0} is of positive Lebesgue measure, and the righthand side of (52) is negative, thus contradicting (47). This
proves that, under the controllability assumption on the
pair (A, B), the condition (48) is necessary for the network
dissipativity. 

Z T

Z T

which is obtained by applying the inequality 12 kα + β k2 6
kα k2 + kβ k2 in an arbitrary Hilbert space to the right-hand
side of (4) and using the boundedness (12) and (13) of the
operators C and D. Also, (17) implies that
2 

Z t
kx(t)k22 6 2e2kAk∞t kx(0)k22 + kBk2∞
ku(τ )k2 dτ
Z 0t


2
2
2kAk∞t
kx(0)k2 + tkBk∞ ku(τ )k22 dτ , (57)
6 2e

for any σ ∈ Tν . For such an input, in view of (49)–(51), the
energy dissipation rate (44) takes the form
Z

Z T

6kCk2∞

Note that, due to (46), the function λ in (49) is symmetric.
Therefore, if the above controllability condition is satisfied,
then for any time horizon t > 0, there exists an admissible
network input u on the time interval [0,t] such that the vector
(36) satisfies

ζ (σ ) if λ (σ ) < 0
(51)
Z(t, σ ) =
0
otherwise

1
2(2π )ν

kGk2∞

in view of the Cauchy-Bunyakovsky-Schwarz inequality and
the boundedness (34) of the operator G. Here, the time-local
square integrability of the input allows (19) to be enhanced
as

of the Hermitian matrix N(σ ) in (45), so that, in accordance
with (46),
(50)
N(σ )ζ (σ ) = λ (σ )ζ (σ ).

S(t) − Ḣ(t) =

0

0

T > 0. (54)

5

(2π )ν +1 |||u|||2, where the function U is given by (27). Although (59) is a stronger condition than (54), the input u
produces the same work W (T ) in (58) over a given time
interval [0, T ] as its “truncated” version

u(t) if 0 6 t 6 T
ΠT (u)(t) :=
(60)
0
otherwise

which is a network counterpart of the positive real property. On the other hand, in an extended setting, the block
Toeplitz matrix G can be formed from differential operators
Gℓ (∂t ) with respect to time (whose entries are, for example,
polynomials of ∂t ). In this case, Theorem 2 is modified by
replacing the function G in (61), given by (33), with the
Fourier-Laplace transform

(with ΠT being an orthogonal projection operator for any
T > 0). This property follows from causality of the network
as an input-output operator u 7→ y and from the fact that the
supply rate S(t) in (31) is a bilinear function of u(t) and
y(t).
Theorem 2: Suppose the network, governed by (1) and
(2), satisfies (7) and (10)–(13), and is endowed with the
supply rate (31) subject to (32) and (34). Also, suppose the
matrix A(σ ) in (8) is Hurwitz for any σ ∈ Tν . Then the
network is passive in the sense of (58) for zero initial states
and arbitrary time-locally square integrable inputs u in (54) if
and only if its spatio-temporal transfer function (28) satisfies
∗

∗

E(ω , σ ) := F(iω , σ ) G(σ ) + G(σ )F(iω , σ ) < 0

G(s, σ ) :=

In particular, if r = m and the extended operator G acts on
the network output as Gy = ẏ, then (64) yields G(s, σ ) =
sIm . In accordance with the structure of the supply rate
(31), this describes the setting when u and y consist of
the corresponding force and position (rather than velocity)
variables. In this case, (61) takes the form

Since the spatio-temporal transfer function F in (28) satisfies
F(−iω , −σ ) = F(iω , σ ) for all ω ∈ R, σ ∈ Rν , and the complex conjugation of a Hermitian matrix preserves positive
semi-definiteness, then (65) is equivalent to

(61)

1
ω ∈ R+ , σ ∈ Tν . (66)
(F(iω , σ ) − F(iω , σ )∗ ) 4 0,
i
Similarly to (63), the condition (66) is a network version of
the negative-imaginary property [7], [15]. In the case when
the matrix A(σ ) in (8) is Hurwitz for all σ ∈ Tν (as assumed
in Theorem 2), there exists a positive definite Hermitian
Cn×n -valued function V in (37) on the torus Tν satisfying
(38), (39) and such that
A(σ )∗ V(σ ) + V(σ )A(σ ) ≺ 0,

S(t)dt = hhu, Gyii
0
Z  Z +∞

1
=
U(iω , σ )∗ G(σ )Y(iω , σ )dω dσ
ν
+1
(2π )
Tν
−∞
Z
1
=
U(iω , σ )∗ G(σ )F(iω , σ )U(iω , σ )dω dσ
(2π )ν +1 R×Tν
Z
1
=
U(iω , σ )∗ E(ω , σ )U(iω , σ )dω dσ , (62)
2(2π )ν +1 R×Tν

σ ∈ Tν .

(67)

In isolation from the environment (when u = 0), the quadratic
stability of the network can be formulated by enhancing the
positive definiteness as

µ := ess inf
λmin (V(σ )) > 0.
ν
σ ∈T

(68)

Then the corresponding block Toeplitz matrix V describes a
positive definite self-adjoint operator on L2 (Zν , Rn ) whose
inverse V −1 is also such an operator and its norm is related
to the quantity µ in (68) by

where the function F in (28) is well-defined on the set
(iR) × Tν (with iR the imaginary axis) due to the matrix
A(σ ) being Hurwitz for all σ ∈ Tν . Here, in view of (29), the
Parseval identity is used together with the Fourier transforms
(26), (27) and the Hermitian matrix E(ω , σ ) ∈ Cm×m in
(61). In view of (62), the positive semi-definiteness of E
almost everywhere in R × Tν implies that W (T ) > 0 for all
admissible inputs vanishing outside the time interval [0, T ].
The necessity of the matrix inequality in (61) for the network
passivity can be obtained by letting T → +∞ and considering
all possible time-space square summable inputs satisfying
(59). 

kV −1 k∞ =

1
.
µ

(69)

Furthermore, the Hamiltonian H, associated with V by (40),
admits the following bounds in terms of the standard L2 norm of the network state:
µ
kV k∞
kx(t)k22 6 H(t) 6
kx(t)k22 .
(70)
2
2
Similarly to finite-dimensional settings, the dissipativity (47)
(or a nonstrict version of the inequality (67)) implies that the
Hamiltonian of the isolated network does not increase with
time (that is, Ḣ 6 0), which, in combination with (70) leads
to
2
kV k∞
2
kx(0)k22 , t > 0 (71)
kx(t)k22 6 H(t) 6 H(0) 6
µ
µ
µ

Note that in the case of r = m and the standard supply rate
mentioned above, G in (33) is the identity matrix, and (61)
reduces to

ω ∈ R, σ ∈ Tν ,

ω ∈ R, σ ∈ Tν . (65)

iω (F(iω , σ ) − F(iω , σ )∗ ) < 0,

Z +∞

F(iω , σ ) + F(iω , σ )∗ < 0,

(64)

ℓ∈Z

for almost all (ω , σ ) ∈ R × Tν . Here, the function G is given
by (33).

Proof: Assuming that the network is initialized at zero,
consider the work W (T ) up until a given time horizon T > 0.
Since an arbitrary admissible input u (in the sense of (54))
can be replaced with ΠT (u) in (60) without affecting the
work, we will assume, without loss of generality, that u
vanishes beyond the time interval [0, T ], and hence, so also
does the supply rate S in (31). For any such input, substitution
of (31) into (58) leads to
W (T ) =

T

∑ν e−iℓ σ Gℓ(s).

(63)
6

Since M ≻ 0, then s ∈ C is an eigenvalue of the matrix A(σ ) in (78) if and only if s2 is an eigenvalue of
−M −1/2 K(σ )M −1/2 . Therefore, if the stiffness operator K
is positive semi-definite (so that 0 delivers a minimum to
the potential energy), then A(σ ) has a purely imaginary
spectrum:
 √
S(A(σ )) = ± i λ : λ ∈ S(K(σ )M −1 )
(80)

(similar bounds, arising from quadratic Hamiltonians, are
also used for quantum harmonic oscillators [6, Eq. (22)]).
In view of (69), the factor kVµk∞ = kV k∞ kV −1 k∞ > 1 on
the right-hand side of (71) is the condition number of the
operator V (which quantifies how far V is from scalar
operators).
VI. PHONON THEORETIC DISPERSION
RELATIONS
We will now consider a class of translation invariant
networks, which, in isolation from the environment, manifest
phonon-like dynamics mentioned in Section III. Suppose
the state vectors xk of the component systems consist of
the conjugate position and momentum variables which are
assembled into n2 -dimensional vectors qk and pk , respectively
(with n being even):
 
q
k ∈ Zν .
(72)
xk = k ,
pk

for any σ ∈ Tν . Such a network is organised as an infinitedimensional spring-mass system whose vibrations are linear superpositions
of “basis” phonons given by (30) with
√
s = ±i λ for
different
values of σ ∈ Tν . The matrix
√
−1/2
−1
K(σ )M −1/2 M −1/2 in (80), which is
K(σ )M = MM
−1/2
isospectral to M
K(σ )M −1/2 < 0, corresponds to the
stiffness-to-mass ratio for one-mode harmonic oscillators.
The dynamic properties of phonons with large wavelengths
2π
|σ | (manifesting themselves in highly correlated motions of
distant subsystems in the network) depend on the asymptotic
behaviour of the matrix K(σ ) in (79) in a small neighbourhood of σ = 0.
Theorem 3: Suppose the isolated translation invariant network (75) with the Hamiltonian (73) has a positive semidefinite stiffness operator K satisfying

Accordingly, the network is assumed to have the following
Hamiltonian:


1
H := ∑ kp j k2M−1 + qTj ∑ K j−k qk .
(73)
2 j∈Zν
k∈Zν
Here, M is a real positive definite symmetric mass matrix of
order 2n , √and use is√made of a weighted Euclidean norm
kvkL := vT Lv = | Lv| of a vector v specified by such
a matrix L. Also, K := (K j−k ) j,k∈Zν is a symmetric block
n n
Toeplitz stiffness operator whose blocks Kℓ ∈ R 2 × 2 specify
the potential energy part of the Hamiltonian. In view of (72)
and (73), the blocks of the corresponding block Toeplitz
matrix V in (40) are given by


Kℓ
0
(74)
Vℓ =
,
ℓ ∈ Zν ,
0 δ0ℓ M −1

∑

∑

ℓ∈Zν

J :=




0 1
⊗ In/2
−1 0

(76)

is the symplectic structure matrix associated with the partitioning (72), with ⊗ the Kronecker product of matrices. A
comparison of (75) with (1) allows the blocks of the matrix
A in (3) to be expressed in terms of (74) and (76) as


0
δ0ℓ M −1
Aℓ = JVℓ =
,
ℓ ∈ Zν .
(77)
−Kℓ
0

K(σ ) =

K(σ ) :=

T

∑ν e−iℓ σ Kℓ .

1
2

ν

∑
j,k=1

= |σ |

Upon substituting (77) into (8), the function A takes the form


0
M −1
,
(78)
A(σ ) =
0
−K(σ )
where

(81)

|ℓ| kKℓ k < +∞.

(82)

Then the phase velocities of the phonons in the network are
uniformly bounded:
o
nω
: σ ∈ Tν \ {0}, ω ∈ R+ , iω ∈ S(A(σ )) < +∞.
sup
|σ |
(83)

Proof: The condition (82) ensures twice continuous
differentiability of the function K in (79). Since K < 0 is
equivalent to K(σ ) < 0 for all σ ∈ Tν , then (81) implies that
K(0) = 0. Hence, K has zero Frechet derivative at the origin:
K′ (0) = 0 (that is, ∑ℓ∈Zν ℓ ⊗ Kℓ = 0). Indeed, otherwise, the
leading term of the linearly truncated Taylor series expansion
K(σ ) = ∑νj=1 σ j ∂σ j K(σ ) σ =0 + o(|σ |) would fail to be a
positive semi-definite matrix for arbitrary σ := (σk )16k6ν .
Therefore, the quadratically truncated Taylor series expansion of K in a small neighbourhood of the origin reduces
to

where δ jk is the Kronecker delta. In application to the
isolated network, the Hamiltonian equations of motion take
the form
 

 

∂p j H
q̇ j
M −1 p j
= J ∂x j H =
ẋ j =
(75)
=
ṗ j
−∂q j H
− ∑k∈Zν K j−k qk
for all j ∈ Zν , where

Kℓ = 0,

ℓ∈Zν
2

2

σ j σk Γ jk + o(|σ |2 )
1
2

ν

∑
j,k=1

!

θ j θk Γ jk + o(1) ,

as σ → 0. (84)

Here, we have associated with σ 6= 0 the unit direction vector
θ := (θk )16k6ν by
σ
θ :=
,
(85)
|σ |

(79)

ℓ∈Z

7

n

n

and the matrices Γ jk ∈ R 2 × 2 specify the second Frechet
derivative of the function K as
Γ jk := ∂σ j ∂σk K(σ )

σ =0

=−

∑

ℓ j ℓk Kℓ

Here, g := ( f (t, jh, kh)) j,k∈Z is a discretization of the forcing
term, and L is a self-adjoint Toeplitz operator which approximates the Laplacian ∆ as
1
(L℧) jk := 2 (q j+1,k + q j−1,k + q j,k+1 + q j,k−1 − 4q jk ). (92)
h
A comparison of (90) and (91) with (75) shows that the
scalar ρ plays the role of the mass matrix M, and β2 L2 is
the stiffness operator K. The latter implies that the Fourier
T
transforms (79) and L(σ ) := ∑ℓ∈Z2 e−iℓ σ Lℓ are related by
K = β2 L2 , and hence, K inherits the property (81) from the
discretized Laplacian L in (92), since L(0) = ∑ℓ∈Z2 Lℓ = 0.

(86)

ℓ∈Zν

in view of (79). The leading term ∑νj,k=1 σ j σk Γ jk in (84)
is a real positive semi-definite symmetric matrix of order n2
for any σ ∈ Rν . Now, since the spectrum (as a set-valued
function) of a matrix depends continuously on it [1], [2],
then, in view of (80),
q
sup r(A(σ )) = maxν λmax (K(σ )M −1 ) < +∞,
(87)
σ ∈Tν

σ ∈T

VII. CONCLUSION
We have considered a class of translation invariant networks of interacting linear systems at sites of a multidimensional lattice. Energy-based input-output properties of
such a network (including passivity, positive realness and the
negative imaginary property) have been considered in terms
of its transfer function representation in the spatio-temporal
frequency domain. We have also discussed quadratic stability
of the network in the isolated regime, and dispersion relations
for phonons in the isolated network. The results of the paper
are applicable to quadratic regulation problems for large
flexible structures with a sparse collocation of sensors and
actuators.

where λmax (·) is the largest eigenvalue of a matrix with a real
spectrum. Hence, the left-hand side of (83) can be unbounded
only because of the asymptotic behaviour of r(A(σ )) as σ →
0 (which corresponds to long-wave phonons). However, it
follows from (80) and (84)–(86) that
p
λmax (K(σ )M −1 )
r(A(σ ))
= lim sup
lim sup
|σ |
|σ |
σ →0
σ →0
v
u
 ν

u1
= max t λmax ∑ θ j θk Γ jk M −1 < +∞. (88)
2
θ ∈Rν : |θ |=1
j,k=1

A combination of (87) and (88) leads to (83). 

R EFERENCES

The speed of propagation of disturbances is described in
physics literature in terms of their group velocity ∂σ ω (see,
for example, [3]). Note that the differentiability is ensured
for those phonons whose frequency ω corresponds to a nondegenerate eigenvalue (of multiplicity one) of the matrix
K(σ )M −1 in (80). We will now provide an example which
illustrates the condition (81). Consider the movement of a
thin isotropic plate of density ρ > 0 and bending stiffness
β > 0, governed by the Kirchhoff-Love equation [8]:

[1] N.J.Higham, Functions of Matrices, SIAM, 2008.
[2] R.A.Horn, and C.R.Johnson, Matrix Analysis, Cambridge University
Press, New York, 2007.
[3] L.D.Landau, and E.M.Lifshitz, Statistical Physics, 3rd Ed.,
Butterworth-Heinemann, Oxford, 1980.
[4] R.Ortega, A.J.van der Schaft, I.Mareels, and B.Maschke, Putting
energy back in control, IEEE Contr. Sys., vol. 21, no. 2, pp. 18–33,
2001.
[5] R.Ortega, A.van der Schaft, B.Maschke, and G.Escobar, Interconnection and damping assignment passivity-based control of port-controlled
Hamiltonian systems, Automatica, vol. 38, no. 4, 2002, pp. 585–596.
[6] I.R.Petersen, A direct coupling coherent quantum observer, 2014
IEEE Conference on Control Applications (CCA), 8–10 October 2014
(arXiv:1408.0399v2 [quant-ph], 8 August 2014).
[7] I.R.Petersen, and A.Lanzon, Feedback control of negative-imaginary
systems, IEEE Contr. Sys., vol. 30, no. 5, 2010, pp. 54–72.
[8] J.N.Reddy, Theory and Analysis of Elastic Plates and Shells, 2nd Ed.,
CRC Press, Taylor and Francis, London, 2007.
[9] A.Kh.Sichani, I.G.Vladimirov, and I.R.Petersen, Covariance dynamics
and entanglement in translation invariant linear quantum stochastic
networks IEEE 54th Annual Conference on Decision and Control
(CDC), December 15-18, 2015, Osaka, Japan, pp. 7107–7112.
[10] G.P.Srivastava, The Physics of Phonons, Taylor & Francis, 1990.
[11] A. van der Schaft, and D. Jeltsema, Port-Hamiltonian Systems Theory:
An Introductory Overview, Foundations and Trends R in Systems and
Control, vol. 1, no. 2–3, 2014, pp. 173–378.
[12] V.Veselago, L.Braginsky, V.Shklover, and C.Hafner, Negative refractive index metamaterials, J. Comp. Theor. Nanosci., vol. 3, 2006, pp.
1–30.
[13] I.G.Vladimirov, and I.R.Petersen, Physical realizability and mean
square performance of translation invariant networks of interacting
linear quantum stochastic systems, 21st Int. Symp. Math. Theor.
Networks Syst., Groningen, The Netherlands, July 7–11, 2014, pp.
1881–1888.
[14] J.C.Willems, Dissipative dynamical systems. Part I: general theory,
Part II: linear systems with quadratic supply rates, Arch. Rational
Mech. Anal., vol. 45, no. 5, 1972, pp. 321–351, 352–393.
[15] J.Xiong, I.R.Petersen, and A.Lanzon, A negative imaginary lemma and
the stability of interconnections of linear negative imaginary systems,
IEEE Trans. Automat. Contr., vol. 55, no. 10, 2010, pp. 2342–2347.

β
ρ∂t2 w = − ∆2 w + f .
(89)
2
Here, w(t, ξ ) is the local deviation of the plate from the
equilibrium (unbent) position, which is a real-valued
 function

ξ1
of time t and the two-dimensional vector ξ := ξ ∈ R2
2

of Cartesian coordinates on the plane, and ∆ := ∂ξ2 + ∂ξ2
1
2
is the Laplacian over the spatial variables. Also, f (t, ξ )
is the density of the external force acting transversally
on the plate. A finite-difference scheme for the numerical
solution of the PDE (89) with spatial step size h > 0 can
be viewed as a translation invariant network on the twodimensional integer lattice Z2 . By denoting ℧ := (℧ jk ) j,k∈Z
the spatial discretization of the solution of (89) at time t > 0
(so that ℧ jk (t) approximates w(t, jh, kh)), the corresponding
functions
(90)
q := ℧,
p := ρ∂t ℧
on the set R+ × Z2 are evolved according to the Hamilton
equations
β
1
q̇ = p,
(91)
ṗ = − L2 q + g.
ρ
2
8

