Technische Universit¨at Berlin
Institut f¨ur Mathematik
Data-driven interpolation of dynamical
systems with delay
Philipp Schulze Benjamin Unger
Preprint 17-2015
Preprint-Reihe des Instituts f¨ur Mathematik
Technische Universit¨at Berlin
http://www.math.tu-berlin.de/preprints
Preprint 17-2015 August 2015
Data-driven interpolation of dynamical systems with delay∗
Philipp Schulze†‡ Benjamin Unger†
January 8, 2016
Abstract
We present a data-driven realization for systems with delay, which generalizes the Loewner
framework. The realization is obtained with low computational cost directly from mea-
sured data of the transfer function. The internal delay is estimated by solving a least-
square optimization over some sample data. Our approach is validated by several exam-
ples, which indicate the need for preserving the delay structure in the reduced model.
Keywords: Data-driven model reduction, delay systems, Loewner framework, moment
matching, tangential interpolation, delay recovery, structure-preserving model reduction
AMS(MOS) subject classification: 93A30, 37M99, 34K06, 93A15
1 Introduction
Nowadays, it is common to describe a physical or chemical system by a mathematical sur-
rogate model. Such models are often given by (partial) differential equations and may be
used for analysis, control, and optimization. The demand for high fidelity models results
in large-scale dynamical systems, for which classical numerical methods may be too time
or memory consuming. In such cases, often an analytically justified and numerically stable
approximation of the input-output map is desirable leading to the field of model order re-
duction (MOR) (for an overview see [1, 3] and the references therein). Many of these MOR
methods require access to the internal dynamics of the full system in terms of a state-space
realization. This assumption can be relaxed by employing data-driven MOR techniques that
construct low-dimensional models directly from measurements. Let H(s)∈Cp,m denote the
transfer function of a system, where mand pare the numbers of inputs and outputs, respec-
tively. Since the input-output behavior of a system is characterized by its transfer function,
measurements of Hseem appropriate to construct a low-dimensional model. We assume
measurements H(λi)ri=wiand ℓiH(µi) = vito be given, i. e., we have
right interpolatation data {(λi, ri, wi)|λi∈C, ri∈Cm, wi∈Cp, i = 1, . . . , ρ}and
left interpolation data {(µi, ℓi, vi)|µi∈C, ℓT
i∈Cp, vT
i∈Cm, i = 1, . . . , ρ}.
(1.1)
∗Research supported by the Collaborative Research Center 910 Control of self-organizing nonlinear systems:
Theoretical methods and concepts of application.
†Institut f¨ur Mathematik, TU Berlin, Germany, {pschulze,unger}@math.tu-berlin.de.
‡The author is supported by the Collaborative Research Center 1029 Substantial efficiency increase in gas
turbines through direct use of coupled unsteady combustion and flow dynamics
1
Examples of measurements yielding data in the form (1.1) are scattering parameters for fre-
quency response objects (S-parameters) and admittance parameters for interconnects, which
can be obtained by a vector network analyzer [2].
In this paper we assume that the transfer function is based on a system with (possibly
unknown) delay and study a generalized realization problem with internal delay: Given the
data (1.1), construct matrices Eρ,A1,ρ,A2,ρ,Bρ, and Cρ, such that the transfer function
Hρ(s) = Cρ(sEρ−A1,ρ −e−τsA2,ρ)−1Bρ(1.2)
with delay τ≥0 interpolates the data, i. e.,
wi=H(λi)ri=Hρ(λi)riand vi=ℓiH(µi) = ℓiHρ(µi) for i= 1, . . . , ρ.
The transfer function (1.2) corresponds to a realization Σρ= (Eρ, A1,ρ, A2,ρ, Bρ, Cρ) of the
form
Eρ˙xρ(t) = A1,ρxρ(t) + A2,ρxρ(t−τ) + Bρu(t),(1.3a)
yρ(t) = Cρxρ(t),(1.3b)
which serves as a low-dimensional model.
A generalized realization problem without delay for the data (1.1) is solved in [16], leading
to the Loewner framework. The resulting realization is a generalized state-space representa-
tion of the form −L˙x(t) = −Lσx(t) + V u(t),
y(t) = Wx(t),(1.4)
where Land Lσare the Loewner matrix and shifted Loewner matrix, respectively, and Vand
Ware matrices consisting of the data (for more details see section 2).
The rate of change of realistic models often depends not only on the current time point,
but also on the configuration at previous time instances, which leads to time-delay systems.
Popular examples are nonlinear optics, chemical reactor systems, population dynamics, and
delayed feedback control (cf. [9] and the references within). Finding a realization of a sys-
tem with delay by means of the Loewner framework results in the system (1.4), that does
not feature the delay term and hence cannot reflect the dynamics of the inherently infinite-
dimensional delay system.
In this paper, we propose a generalization of the Loewner realization to systems with
delay. The main contributions are described in the following.
•First, we present a generalization to delay systems based on extensions of the Loewner
matrix and the shifted Loewner matrix. This approach provides a realization as in (1.3)
with the coupling A2,ρ =−Eρ(see section 3).
•Since the coupling A2,ρ =−Eρappears to be rather restrictive, we consider general
conditions for interpolating the transfer function (1.2) in section 4. For this purpose,
the moment matching framework from [4] (see also [10] for an earlier contribution) is
extended to DAE systems with delay.
•Based on the conditions for moment matching, we derive a framework to obtain a
realization (1.3) with the coupling A2,ρ =αEρ+βA1,ρ with scalar parameters αand β
(see section 5). Note that this is a more general form than the realization from section 3,
which is the special case with α=−1 and β= 0.
2
•Since the delay is unknown in general, we propose a methodology to determine a delay
τκ,ρ, which is optimal in the sense that it minimizes the interpolation error for a set of
sampled test data of the transfer function (see section 6). In the same fashion, optimal
parameters ακ,ρ and βκ,ρ may be calculated.
In section 7 we apply the proposed framework to several examples. A comparison to alter-
native techniques from the literature including the original Loewner framework reveals the
necessity of preserving the delay structure in the realization.
2 Notation and Preliminary Results
Recall that the realization (1.3) is given by
Eρ˙xρ(t) = A1,ρxρ(t) + A2,ρxρ(t−τ) + Bρu(t),(2.1a)
yρ(t) = Cρxρ(t),(2.1b)
where xρ(t)∈Rrfor r≤ρ,u(t)∈Rm, and yρ(t)∈Rpdenote, respectively, the state,input,
and output of the model. As common in the delay literature, the right-hand derivative d
dtof
a piecewise smooth function fis denoted by ˙
f[13]. The symbol Instands for the identity
matrix of dimension n×n. The input uis assumed to be sufficiently smooth and the system
(2.1) is equipped with the initial condition (also called history function)
x(t) = ϕ(t) for t∈(−τ, 0],(2.2)
which is assumed to be identically zero, i. e., ϕ≡0. If det(sEρ−A1,ρ −e−τsA2,ρ) is not
vanishing identically, then the Laplace transformation of (2.1) yields the transfer function
Hρ(s) = Cρ(sEρ−A1,ρ −e−τsA2,ρ)−1Bρ.(2.3)
For convenience, we call the set {s∈C|det (sEρ−A1,ρ −e−τsA2,ρ)= 0}the resolvent set of
sEρ−A1,ρ −e−τsA2,ρ. For τ > 0, the system Σρas in (2.1) can be solved on consecutive
time intervals ((K−1)τ, Kτ], such that (2.1a) reduces to the associated differential-algebraic
equation (DAE) [12]
Eρ˙xρ(t) = A1,ρxρ(t) + fρ(t) (2.4)
with inhomogeneity fρ(t) = A2,ρxρ(t−τ) + Bρu(t). Note that we use the terminology DAE
here, since we allow Eρto be singular. This procedure is known as the (Bellman) method of
steps [6, 12]. For τ= 0, (2.1a) reduces to the DAE
Eρ˙xρ(t) = Aρx(t) + Bρu(t) (2.5)
with Aρ:=A1,ρ +A2,ρ. For convenience, we also refer to (2.5) as associated DAE. In both
cases, a unique solution of (2.4) and (2.5) is given if and only if the pencils sEρ−Aρand
sEρ−A1,ρ, respectively are regular [14], i. e., det(sEρ−Aρ) and det(sEρ−A1,ρ) do not vanish
identically. For τ= 0, the regularity of the pencil of (2.5) directly implies the existence of
the transfer function (2.3). In the case τ > 0 the existence of the transfer function (2.3) is
guaranteed by the subsequent lemma.
Lemma 2.1 If τ > 0and the pencil sEρ−A1,ρ is regular, then det(sEρ−A1,ρ−e−τsA2,ρ)≡ 0.
3
Proof. Suppose there exists v= 0 such that (sEρ−A1,ρ −e−τsA2,ρ)v= 0 for all s∈C.
This implies (since (Eρ, A1,ρ) is regular) that
eτsv= (sEρ−A1,ρ)−1A2,ρv
for all snot in the spectrum of (Eρ, A1,ρ). This contradicts v= 0, because (sEρ−A1,ρ)−1is
a rational function in s.
We briefly recall the Loewner realization introduced in [16]. The Loewner matrix L
and shifted Loewner matrix Lσare defined as solutions of the Sylvester equations
LΛ−ML=LW −V R, (2.6)
LσΛ−MLσ=LWΛ −MV R, (2.7)
respectively, where
Λ= diag(λ1, . . . , λρ), M = diag(µ1, . . . , µρ), R =[r1. . . rρ],
W=[w1. . . wρ], L =⎡
⎢
⎣
ℓ1
.
.
.
ℓρ
⎤
⎥
⎦,and V=⎡
⎢
⎣
v1
.
.
.
vρ
⎤
⎥
⎦.
Direct computation shows that if λi=µjfor i, j = 1, . . . , ρ, then the unique solutions of (2.6)
and (2.7), respectively, are given by
[L]i,j =virj−ℓiwj
µi−λj
and [Lσ]i,j =µivirj−λjℓiwj
µi−λj
.(2.8)
Remark 2.2 Both matrices can be assembled efficiently via matrix-matrix operations of size
ρ×ρusing standard tools for scalar, vector, and matrix operations (BLAS). More precisely,
the numerator is given by V R−LW for the Loewner matrix and MV R−LWΛ for the shifted
Loewner matrix. The denominator is implemented as µeT−eλT, where µ, λ ∈Cρare vectors
containing µiand λiand e∈Rρis a vector filled with ones.
Theorem 2.3 ([16, Lemma 5.1]) Let det(˜sL−Lσ)= 0 for all ˜s∈ {λi} ∪ {µi}. Then the
system
−L˙xρ(t) = −Lσxρ(t) + V u(t),
yρ(t) = Wxρ(t)(2.9)
is a minimal realization of an interpolant of the data, i. e., its transfer function
Hρ(s) = W(Lσ−sL)−1V
interpolates the data (1.1).
Let εdenote the machine precision. If det(˜sL−Lσ) = O(ε) for some ˜s∈ {λi} ∪ {µi}, then
one can use the truncated singular value decomposition (SVD) [11] of sL−Lσto truncate
the numerically vanishing singular values as in the next theorem, originally stated in [16].
4
Theorem 2.4 ([16, Theorem 5.1]) Suppose that
rank (˜sL−Lσ) = rank [L Lσ]= rank [L
Lσ]=:rfor all ˜s∈ {λi}∪{µi}.
Then a minimal realization of an interpolant of the data is given by the system
−Y∗LX˙xr(t) = −Y∗LσXxr(t) + Y∗V u(t),
yr(t) = WXxr(t),(2.10)
where Y∈Cρ,r and X∈Cρ,r are the orthogonal factors of the truncated SVD
˜sL−Lσ=Y ΣX∗
for any ˜s∈ {λi}∪{µi}, where Σ∈Cr,r is positive definite and diagonal.
Recently, a generalization of the Loewner framework for a special class of delay systems, where
A1,ρ = 0 in (2.1), was introduced in [18]. There, the interpolant is constructed by means of
Theorem 2.3 leading to the subsequent result, written here in slightly different notation.
Theorem 2.5 ([18, Theorem 3]) Suppose that λieτλi=µjeτµjfor i, j = 1, . . . , ρ. Let
L(τ)and L(τ)
σdenote the Loewner matrix and shifted Loewner matrix, respectively, associated
with the transformed data (λieτλi, ri,e−τλiwi)and (µieτµi, ℓi,e−τµivi)for i= 1, . . . , ρ. If
det (˜sL(τ)−e−τ˜sL(τ)
σ)= 0 for all ˜s∈ {λi}∪{µi}, then the transfer function
Hρ(s) = We−τΛ (e−τsL(τ)
σ−sL(τ))−1e−τM V
of the system
−L(τ)˙xρ(t) = −e−τsL(τ)
σxρ(t−τ)+e−τM V u(t)
yρ(t) = We−τΛxρ(t)
is an interpolant of the original data (1.1).
Again, the condition det(˜sL(τ)−e−τ˜sL(τ)
σ)= 0 can be relaxed by applying the truncated SVD
as in Theorem 2.4.
3 A Loewner Framework for Systems with Delay
In this section, we aim at extending the Loewner matrix and shifted Loewner matrix to the
time delay case. Suppose that µj+ e−τµj=λi+ e−τλifor all i, j = 1, . . . , ρ. Define the
matrices Tand Tσ
[T]i,j :=virj−ℓiwj
µi+ e−τµi−(λj+ e−τλj)and
[Tσ]i,j :=(µi+ e−τµi)virj−(λj+ e−τλj)ℓiwj
µi+ e−τµi−(λj+ e−τλj).
(3.1)
5
Lemma 3.1 Suppose that the denominators in (3.1) are not zero. Then the matrices Tand
Tσare the unique solutions of the Sylvester equations
T(Λ+ e−τΛ)−(M+ e−τM )T=LW −V R, (3.2)
Tσ(Λ+ e−τΛ)−(M+ e−τM )Tσ=LW(Λ+ e−τΛ)−(M+ e−τM )V R, (3.3)
respectively.
Proof. It is well known that under the stated assumptions the Sylvester equations possess
unique solutions (cf. [1]). Multiplication of (3.2) and (3.3) from the left by eT
iand from the
right by ejyields the matrices as in (3.1).
Similar as in [16], we get the following identities.
Corollary 3.2 Suppose that the denominators in (3.1) are not zero. Then we have
Tσ−T(Λ+ e−τΛ) = V R and Tσ−(M+ e−τM )T=LW.
Proof. Multiplying (3.2) by (Λ+ e−τΛ) from the right and subtracting it from (3.3) yields
(Tσ−T(Λ+ e−τΛ)−V R)(Λ+ e−τΛ)
−(M+ e−τM )(Tσ−T(Λ+ e−τΛ)−V R)= 0.(3.4)
Since Λ+ e−τΛ and M+ e−τM have no eigenvalues in common (by assumption), the solution
of the Sylvester equation (3.4) is zero, which yields the first equality. By the same argument,
we get the second identity if we multiply (3.2) by (M+ e−τM ) from the left and subtract it
from (3.3).
Theorem 3.3 Suppose that the denominators in (3.1) are not zero and
det (Tσ−(˜s+ e−τ˜s)T)= 0 for all ˜s∈ {λi}∪{µi}.(3.5)
Then the delay descriptor system given by Σρ= (−T,−Tσ,T, V, W)with associated transfer
function Hρ(s) = W(Tσ−sT−e−τsT)−1Vinterpolates the data, i. e.,
Hρ(λi)ri=wiand ℓiHρ(µi) = vifor i= 1, . . . , ρ.
Moreover, the pencil of the associated DAE of Σρis regular.
Proof. Multiply (3.2) by (s+ e−τs) and subtract it from (3.3) to obtain
(Tσ−T(s+ e−τs))(Λ+ e−τΛ)−(M+ e−τM )(Tσ−T(s+ e−τs))=
LW ((Λ+ e−τΛ)−(s+ e−τs)Iρ)−(M+ e−τM −(s+ e−τs)Iρ)V R. (3.6)
Multiply (3.6) from the right by the i-th unit vector eito obtain
(Tσ−T(s+ e−τs))(λi+ e−τλi)ei−(M+ e−τM )(Tσ−T(s+ e−τs))ei=
LW ((λi+ e−τλi)−(s+ e−τs))ei−(M+ e−τM −(s+ e−τs)Iρ)V ri,
6
which is equivalent to
((λi+ e−τλi)Iρ−(M+ e−τM ))(Tσ−T(s+ e−τs))ei=
LW ((λi+ e−τλi)−(s+ e−τs))ei−(M+ e−τM −(s+ e−τs)Iρ)V ri.
Setting s=λiimplies
((λi+ e−τλi)Iρ−(M+ e−τM ))(Tσ−T(λi+ e−τλi))ei
=((λi+ e−τλi)Iρ−(M+ e−τM ))V ri.(3.7)
By assumption, the matrix (λi+ e−τλi)Iρ−(M+ e−τM ) is nonsingular. In particular, (3.7)
is equivalent to
(Tσ−T(λi+ e−τλi))ei=V ri
and assumption (3.5) implies
ei=(Tσ−T(λi+ e−τλi))−1V ri.
Multiplication by Wfrom the left yields
Hρ(λi)ri=wi=Wei=W(Tσ−T(λi+ e−τλi))−1V ri.
The proof for the left interpolation conditions works analogously by multiplying (3.6) from
the left by eT
iand setting s=µi. The regularity of the pencil of the associated DAE follows
directly from (3.5).
Remark 3.4 Theorem 3.3 is a generalization of Theorem 2.3 in the following sense. If τ= 0,
then system (1.3) reduces to the descriptor system
Eρ˙xρ(t) = Aρxρ(t) + Bρu(t),
yρ(t) = Cρxρ(t)
with Aρ=A1,ρ +A2,ρ and the realization via the matrices Tand Tσyields
[Eρ]i,j =−[T]i,j =−virj−ℓiwj
µi+ 1 −(λj+ 1) =−[L]i,j ,
[Aρ]i,j = [T]i,j −[Tσ]i,j =virj−ℓiwj
µi−λj
−(µi+ 1)virj−(λj+ 1)ℓiwj
µi−λj
=−[Lσ]i,j ,
which are the Loewner matrix and shifted Loewner matrix, respectively.
Remark 3.5 The realization given in Theorem 3.3 is restrictive in the sense that Eρand A2,ρ
coincide apart from a constant factor. Thus, there are delay descriptor systems whose dy-
namical behavior is not captured accurately by this approach. To circumvent this restriction,
we consider a more general framework in the next section.
7
4 Moment Matching for Linear Delay Descriptor Systems
Subsequently, we present a general framework for matching the moments of a linear time-
invariant time-delay descriptor system of the form
E˙x(t) = A1x(t) + A2x(t−τ) + Bu(t),
y(t) = Cx(t)(4.1)
with state-space dimension n, input dimension m, output dimension p, and matrices E,A1,
A2,B, and Cof appropriate size. The transfer function of the system is given by
H(s) = C(sE −A1−e−τsA2)−1B
and we assume that ˜sE −A1−e−τ˜sA2is nonsingular for some ˜s∈C. We generalize the
moment matching approach as introduced in [4] to delay DAE systems of the form (4.1). A
similar method for delay differential equations is derived in [19].
Moment matching is a technique to reduce a large-dimensional system with given state-
space representation by interpolation of the transfer function and its derivatives. In the
next section, we revert this process to derive conditions for a data-driven realization. In the
data-driven approach, we assume that we do not have information about the derivatives of
the transfer function. This is the reason for only considering 0-moments (see Definition 4.1
below) and ignoring moments of higher order. In the case that we have multiple inputs and
multiple outputs (MIMO) the transfer function His a matrix function we interpolate along
interpolation directions. This is known as tangential interpolation [10].
Definition 4.1 (0-moment) Consider λi, µi∈Cin the resolvent set of sE −A1−e−τsA2.
The right 0-moment η(λi, ri)of the system (4.1) at λiin direction riis defined as
η(λi, ri):=C(λiE−A1−A2e−τλi)−1Bri=H(λi)ri.
The left 0-moment η(µi, ℓi)of the system (4.1) at µiin direction ℓiis defined as
η(µi, ℓi):=ℓiC(µiE−A1−A2e−τµi)−1B=ℓiH(µi).
Note that the interpolation data (1.1) consists of the interpolation points λiand µi, the
directions riand ℓi, and the right and left 0-moments. An equivalent notion of moment is
presented in the following proposition.
Proposition 4.2 Consider system (4.1) and complex numbers λi, µi∈Cin the resolvent set
of sE −A1−e−τsA2.
1. The right 0-moment of system (4.1) at λiin direction riis given by
η(λi, ri) = Cπi,
where πiis the unique solution of the linear algebraic system
λiEπi=A1πi+A2πie−τλi+Bri.(4.2)
8
2. The left 0-moment of system (4.1) at µiin direction ℓiis given by
η(µi, ℓi) = ψiB,
where ψiis the unique solution of the linear algebraic system
µiψiE=ψiA1+ψiA2e−τµi+ℓiC. (4.3)
Proof. Since λiis in the resolvent set, there exists a unique solution of (4.2) given by
πi=(λiE−A1−A2e−τλi)−1Briand hence Cπi=η(λi, ri). The second statement follows
analogously.
Remark 4.3 As a direct consequence from Proposition 4.2, the set of right 0-moments cor-
responding to {λ1, . . . , λρ}and {r1, . . . , rρ}is given by CΠ, where Πis the unique solution
of the matrix equation
EΠΛ =A1Π+A2Πe−τΛ +BR, (4.4)
where Λand Rare defined as in section 2. Analogously, the set of left 0-moments correspond-
ing to {µ1, . . . , µρ}and {ℓ1, . . . , ℓρ}is given by ΨB, where Ψis the unique solution of the
matrix equation
MΨE =ΨA1+ e−τM ΨA2+LC. (4.5)
With the help of (4.4) and (4.5), we can derive families of reduced order models achieving
moment matching at the given interpolation points.
Theorem 4.4 Consider the time-delay descriptor system Σρgiven in (1.3) and let right and
left interpolation data of the transfer function of (4.1) be given as in (1.1), where CΠ =W
and ΨB =V(see Proposition 4.2 and Remark 4.3). Furthermore, assume that the matrix
˜sEρ−A1,ρ −e−τ˜sA2,ρ is nonsingular for all ˜s∈ {λ1, . . . , λρ}∪{µ1, . . . , µρ}.
1. System (1.3) matches the right 0-moments CΠ if and only if
CΠ =CρP ,
where Πis the unique solution of (4.4) and Pis the unique solution of the matrix
equation
EρPΛ =A1,ρP+A2,ρPe−τΛ +BρR. (4.6)
2. System (1.3) matches the left 0-moments ΨB if and only if
ΨB =P Bρ,
where Ψis the unique solution of (4.5) and Pis the unique solution of the matrix
equation
MPEρ=PA1,ρ + e−τM PA2,ρ +LCρ.(4.7)
9
Proof. The claim is a straightforward consequence of the notion of moments as introduced
in Proposition 4.2 and Remark 4.3. Following this, the right 0-moments of system (1.3) are
given by CρP, where Pis the unique solution of the matrix equation (4.6). Likewise, the right
0-moments of the original system (4.1) are given by CΠ, where Πis the unique solution of
(4.4). Consequently, the right moments are matched if and only if CΠ =CρP. The matching
of the left 0-moments is proven analogously.
Corollary 4.5
1. The family of reduced order models
Eρ˙xρ(t)=(EρΛ−A2,ρe−τΛ −BρR)xρ(t) + A2,ρxρ(t−τ) + Bρu(t),
yρ(t) = CΠxρ(t)(4.8)
matches the right 0-moments CΠ if ˜sEρ−(EρΛ−A2,ρe−τΛ −BρR)−e−τ˜sA2,ρ is non-
singular for all ˜s∈ {λ1, . . . , λρ}.
2. The family of reduced order models
Eρ˙xρ(t)=(MEρ−e−τM A2,ρ −LCρ)xρ(t) + A2,ρxρ(t−τ) + ΨBu(t),
yρ(t) = Cρxρ(t)(4.9)
matches the left 0-moments ΨB if ˜sEρ−(MEρ−e−τM A2,ρ −LCρ)−e−τ˜sA2,ρ is non-
singular for all ˜s∈ {µ1, . . . , µρ}.
Proof. The assertions follow directly from Theorem 4.4 by setting P=Iρ=P. This
leads to A1,ρ =EρΛ−A2,ρe−τΛ −BρRand Cρ=CΠ for the family matching the right
0-moments and A1,ρ =MEρ−e−τM A2,ρ −LCρand Bρ=ΨB for the family matching the
left 0-moments.
Remark 4.6 Note that the families of reduced order models (4.8) and (4.9) are parameterized
by Eρ,A2,ρ,Bρ, and Cρ. These parameters are only restricted by the generic conditions stated
in Corollary 4.5. This freedom may be exploited in order to tailor the reduced order model to
additional requirements, e. g., preserving structures or properties of the original system (1.3).
To find a family of reduced order models that matches both, the left and right 0-moments,
we compare the coefficient matrices of (4.8) and (4.9). This results in the following theorem.
Theorem 4.7 The family of reduced order models
Eρ˙xρ(t) = (EρΛ−A2,ρe−τΛ −ΨBR)xρ(t) + A2,ρxρ(t−τ) + ΨBu(t)
yρ(t) = CΠxρ(t)(4.10)
matches the 0-moments ΨB and CΠ if the matrix ˜sEρ−(EρΛ−A2,ρe−τΛ −ΨBR)−e−τ˜sA2,ρ
is nonsingular for all ˜s∈ {λi}∪{µi}and the Sylvester equation
EρΛ−MEρ+ e−τM A2,ρ −A2,ρe−τΛ =ΨBR −LCΠ (4.11)
is satisfied.
10
A more general version of Theorem 4.7 may be derived. Instead of coinciding coefficients, it
is sufficient to require that the family (4.8) also matches the left 0-moments of (4.9) or that
the family (4.9) also matches the right 0-moments of (4.8). This leads to an additional degree
of freedom, which we omit for the sake of clarity.
The Sylvester equation (4.11) does not depend on the system matrices E,A1,A2,B
and Cexplicitly, since we have the identities CΠ =Wand ΨB =V. Thus, Theorem 4.7
allows us to construct reduced order models from data only. The obtained realization (4.10)
is parameterized by the matrices Eρand A2,ρ, which are coupled via (4.11).
Corollary 4.8 For the special case A2,ρ =−Eρthe unique solution of the Sylvester equation
(4.11) is given by the matrix Eρ=−T. Moreover, from Lemma 3.2 we have
Aρ=(EρΛ−A2,ρe−τΛ −ΨBR)=Eρ(Λ+ e−τΛ)−V R =−Tσ.
Hence, for the particular case A2,ρ =−Eρ, system (4.10) is equivalent to the realization in
Theorem 3.3.
Remark 4.9 By solving the equations (4.4) and (4.5) for BR and LC, respectively, and
substituting the obtained expressions into the equations (4.10) and (4.11), it is easy to show
that a reduced order model matching the left and right 0-moments is given by the projection
Eρ=ΨEΠ, A1,ρ =ΨA1Π, A2,ρ =ΨA2Π, Bρ=ΨB, and Cρ=CΠ.
Example 4.10 Consider the system
[1 0
0 0][˙x1(t)
˙x2(t)]=[11
2
0 1][x1(t)
x2(t)]+[01
2
0 1][x1(t−π)
x2(t−π)]+[0
1]u(t),
y(t) = [0 1][x1(t)
x2(t)](4.12)
with transfer function H(s) = −(1 + e−πs)−1. We choose ρ= 1, λ= 0, and µ=1
2. Since
(4.12) is a single-input single-output (SISO) system, we set ℓ=r= 1 for the tangential
directions. The projection matrices Ψand Πare given by
Ψ=[0−(1 + e−τ
2)−1]and Π=[1
2
−1
2]
and the reduced order model is given by
0 = 1
2(1+e−τ
2)−1xρ(t) + 1
2(1+e−τ
2)−1xρ(t−τ)−(1+e−τ
2)−1u(t),
yρ(t) = −1
2xρ(t).
Simple calculations show that the left and right 0-moments are matched by the transfer
function of the reduced order model.
Remark 4.11 In [18] the authors consider the special case A1,ρ = 0 (cf. Theorem 2.5 above).
In terms of Corollary 4.5 and Theorem 4.7, this implies that Eρmust satisfy the Sylvester
equation
EρΛeτΛ −eτM MEρ=V ReτΛ −eτM LW
and A2,ρ is given by
A2,ρ = (EρΛ−V R) eτΛ = eτM (Eρ−LW).
Consequently, the result obtained in [18] is a special case of the realization presented in
Theorem 4.7.
11
5 Generalized Delay Loewner Framework
Recall that the reduced models (4.10) allow for choosing either Eρor A2,ρ. Inspired by the
mass- and stiffness-proportional damping, normally referred to as Rayleigh damping [8] in
mechanical systems, we make the ansatz A2,ρ =αEρ+βA1,ρ with some constants αand β.
For the special setting α=−1 and β= 0, we recover the framework presented in section 3.
Lemma 5.1 Suppose that 1 + βe−τ˜s= 0 for all ˜s∈ {λi}∪{µi}. The reduced order model
Eρ˙xρ(t) = A1,ρxρ(t)+(αEρ+βA1,ρ)xρ(t−τ) + Bρu(t),
yρ(t) = Cρxρ(t)(5.1)
matches the 0-moments ΨB and CΠ if the matrix ˜sEρ−A1,ρ −e−τ˜s(αEρ+βA1,ρ)is non-
singular for all ˜s∈ {λi}∪{µi},
Cρ=CΠ, Bρ=ΨB,
A1,ρ =(Eρ(Λ−αe−τΛ)−BρR)(Iρ+βe−τΛ)−1,
and Eρsatisfies the Sylvester equation
Eρ(Λ−αe−τΛ)(Iρ+βe−τΛ)−1−(Iρ+βe−τM )−1(M−αe−τM )Eρ=
BρR(Iρ+βe−τΛ)−1−(Iρ+βe−τM )−1LCρ.(5.2)
Proof. For the right 0-moments, from Corollary 4.5 we get the conditions Cρ=CΠ and
A1,ρ =EρΛ−A2,ρe−τΛ −BρR=EρΛ−(αEρ+βA1,ρ)e−τΛ −BρR,
which is equivalent to
A1,ρ =(Eρ(Λ−αe−τΛ)−BρR)(Iρ+βe−τΛ)−1.(5.3)
Analogously, to match the left 0-moments, the condition Bρ=ΨB must be satisfied and the
matrix A1,ρ is given by
A1,ρ = (Iρ+βe−τM )−1((M−αe−τM )Eρ−LCρ).(5.4)
By comparison of the coefficients (as in the proof of Theorem 4.7), we deduce that (5.3) and
(5.4) must coincide, which yields the Sylvester equation (5.2).
Corollary 5.2 Suppose that the assumptions of Lemma 5.1 are satisfied and the matrices
(Λ−αe−τΛ)(Iρ+βe−τΛ)−1and (Iρ+βe−τM )−1(M−αe−τM )
have no common eigenvalues. Then the unique solution of the Sylvester equation (5.2) is
given by
[Eρ]i,j =
virj
1 + βe−τλj−ℓiwj
1 + βe−τµi
λj−αe−τλj
1 + βe−τλj−µi−αe−τµi
1 + βe−τµi
.(5.5)
12
Proof. Using the identities
Cρ=CΠ =Wand Bρ=ΨB =V
the result follows by multiplying (5.2) from the left by eT
iand from the right by ej.
Lemma 5.3 Under the assumptions of Corollary 5.2, the matrix A1,ρ from Lemma 5.1 sat-
isfies the Sylvester equation
A1,ρ(Λ−αe−τΛ)(Iρ+βe−τΛ)−1−(Iρ+βe−τM )−1(M−αe−τM )A1,ρ =
(Iρ+βe−τM )−1(M−αe−τM )V R(Iρ+βe−τΛ)−1
−(Iρ+βe−τM )−1LW(Λ−αe−τΛ)(Iρ+βe−τΛ)−1
(5.6)
and is given by
[A1,ρ]i,j =
µi−αe−τµi
1 + βe−τµi
virj
1 + βe−τλj−ℓiwj
1 + βe−τµi
λj−αe−τλj
1 + βe−τλj
λj−αe−τλj
1 + βe−τλj−µi−αe−τµi
1 + βe−τµi
.(5.7)
Proof. The proof proceeds analogously to the proof of Lemma 3.2. Multiplying (5.2) with
(Λ−αe−τΛ)(I+βe−τΛ)−1and subtracting it from (5.6) results in a Sylvester equation with a
unique solution that yields (5.4). The result for equation (5.3) follows likewise. Multiplying
(5.6) by eT
ifrom the left and ejfrom the right yields the representation (5.7).
Summarizing the previous discussion, we have shown the following theorem.
Theorem 5.4 Let
det ((˜s−αe−τ˜s)Eρ−(1 + βe−τ˜s)A1,ρ)= 0 for all ˜s∈ {λi}∪{µi}.(5.8)
Then, under the assumptions of Corollary 5.2, the system
Eρ˙xρ(t) = A1,ρxρ(t)+(αEρ+βA1,ρ)xρ(t−τ) + V u(t),
yρ(t) = Wxρ(t)
with Eρand A1,ρ as in (5.5) and (5.7), respectively, is an interpolant of the data, i. e., its
transfer function
Hρ(s) = W((s−αe−τs)Eρ−(1 + βe−τs)A1,ρ)−1V
interpolates the data (1.1). Moreover, the pencil of the associated DAE is regular.
Remark 5.5 Following the idea of [18], we can (formally) derive the same results by only
employing the Loewner realization theorem (Theorem 2.3). Using the proportional ansatz
A2,ρ =αEρ+βA1,ρ with some constants αand β, we can rewrite the transfer function as
H(s) = Cρ(sEρ−A1,ρ −e−τs (αEρ+βA1,ρ))−1Bρ
=Cρ(s−αe−τs
1 + βe−τs Eρ−A1,ρ)−1
Bρ
1
1 + βe−τs
=G(f(s)) 1
1 + βe−τs ,
13
where Gis the transfer function of a generalized state-space system without delay and fis
given by f(s) = s−αe−τs
1+βe−τs . The Loewner matrices are then constructed for the transfer function
G(s) with transformed data (f(λi), ri(1+βe−τλi)−1, wi) and (f(µi),(1+βe−τµi)−1ℓi, vi). Note
that this allows to efficiently implement the system matrices similar as in Remark 2.2.
With Remark 5.5 we can immediately relax assumption (5.8) in Theorem 5.4.
Theorem 5.6 Suppose that the assumptions of Corollary 5.2 are satisfied and that
rank (˜s−αe−τ˜s
1 + βe−τ˜sEρ−A1,ρ)= rank [EρA1,ρ]= rank [Eρ
A1,ρ]=:r
holds for all ˜s∈ {λi}∪{µi}. Then, an interpolant of the data (1.1) is given by the system
−Y∗EρX˙xr(t) = −Y∗A1,ρXxr(t) + Y∗(αEρ+βA1,ρ)Xxr(t−τ) + Y∗V u(t),
yr(t) = WXxr(t),(5.9)
where Y∈Cρ,r and X∈Cρ,r are the orthogonal factors of the truncated SVD
˜s−αe−τ˜s
1 + βe−τ˜sL−Lσ=Y ΣX∗
for any ˜s∈ {λi}∪{µi}, where Σ∈Cr,r is positive definite and diagonal.
Proof. The proof follows directly from Remark 5.5 and Theorem 2.4.
6 Delay Reconstruction
So far, we have assumed that the delay τ≥0 is known, which in general is not the case. In
practical applications, e. g. neuronal processing [20], only a range for the delay time is known
from experiments. Checking the moment matching conditions as stated in Theorem 4.4
implies that the reduced order model matches the moments Vand Wfor any choice of the
delay time. We propose the following strategy to recover a delay term τκ,ρ ≈τ. For some
κ < ρ, split the interpolation data (1.1) as
{(λi, ri, wi)|λi∈C, ri∈Cm, wi∈Cp, i = 1, . . . , κ, κ + 1, . . . , ρ}and
{(µi, ℓi, vi)|µi∈C, ℓT
i∈Cp, vT
i∈Cm, i = 1, . . . , κ, κ + 1, . . . , ρ}(6.1)
and construct the matrices Eκ(τ), A1,κ(τ), A2,κ(τ), Bκ, and Cκonly from the data (λi, ri, wi)
and (µi, ℓi, vi) for i= 1, . . . , κ. Note that the transfer function
Hκ(s, τ) = Cκ(sEκ(τ)−A1,κ(τ)−e−τsA2,κ(τ))−1Bκ
interpolates the data independently from τ. The remaining data (λi, ri, wi) and (µi, ℓi, vi)
for i=κ+ 1, . . . , ρ, in the following called test data, is used to fit the transfer function in a
least-squares sense, i. e., we solve the minimization problem
min
τ∈[τ−,τ+]Jκ,ρ(τ):=1
2
ρ
∑
i=κ+1 (∥Hκ(λi, τ)ri−wi∥2+∥ℓiHκ(µi, τ)−vi∥2)(6.2)
14
over all delay times τin a realistic range [τ−, τ+] and denote the minimizer by τκ,ρ. With
the ‘optimal’ delay τκ,ρ, we can rebuild the realization with the complete data set as in
Theorem 5.6.
Note that the cost functional Jκ,ρ is of nonconvex type and gradient-based optimization
might yield a local instead of the global minimum. Hence, either a good choice of the ini-
tial value for the optimization is required, or global optimization methods like evolutionary
algorithms (see [21] and the references within) can be employed.
In the previous section we have introduced the additional parameters αand β. These
parameters are degrees of freedom, which can be used to tailor the realization to additional
data. Hence, the strategy used to find the delay time may be extended to determine good
values for αand β. More precisely, the cost functional Jκ,ρ is given by
Jκ,ρ(τ, α, β):=1
2
ρ
∑
i=κ+1 (∥Hκ(λi, τ, α, β)ri−wi∥2+∥ℓiHκ(µi, τ, α, β)−vi∥2)
and the minimization in (6.2) can be performed over τ,α, and βsimultaneously. Note that in
contrast to the delay time τ, there are in general no a priori known ranges for the parameters
αand β. The choice of good upper and lower bounds for αand βis not within the scope
of this paper. Instead, we propose to optimize over different intervals, which may be picked
randomly, and then use the result that yields the smallest value of Jκ,ρ.
7 Examples
In all subsequent examples, we denote the transfer function of the original model with H.
The transfer function of the Loewner realization (Theorem 2.3) is denoted by HL. By HT
we denote the transfer function based on transformed data from [18] (Theorem 2.5) and the
generalized delay Loewner realization (Theorem 5.4) is referred to as HD. For details of the
Table 7.1: Transfer functions for the different approaches
Transfer Function Reference
H(s) original model
HL(s) = W(Lσ−sL)−1V[16]
HT(s) = We−τΛ (e−τsL(τ)
σ−sL(τ))−1e−τM V[18]
HD(s) = W(sEρ−A1,ρ −e−τs (αEρ+βA1,ρ))−1VTheorem 5.4
transfer functions, see Table 7.1. Whenever necessary, the redundant parts in the realizations
are truncated by means of the truncated SVD as outlined in Theorem 2.4 and Theorem 5.6.
Example 7.1 Consider again the system from Example 4.10 given by
[1 0
0 0][˙x1(t)
˙x2(t)]=[11
2
0 1][x1(t)
x2(t)]+[01
2
0 1][x1(t−π)
x2(t−π)]+[0
1]u(t),
y(t) = [0 1][x1(t)
x2(t)]
15
together with the interpolation data λ= 0 and µ=1
2. Recall that the transfer function H
is given by H(s) = −(1 + e−πs)−1and note that A2=A1−E, i. e., the considered system
exhibits the structure A2=αE +βA1. Comparison of the transfer functions
HL(s) = −1
sπ(1 −1
e)+2, HD(s) = −1
1+e−πs ,and HT(s) = −1
sπ(1 −1
e) + 2e−πs
of the different realizations shows that the transfer function HDis the only one matching the
full model exactly. This is due to the structure of the system, which can neither be captured
by HLnor by HT.
Example 7.2 For Example 7.1 we employ the methodology of section 6 to recover the delay
time τ. We use the additional interpolation point λ= 10. For the minimization of the cost
functional (6.2), we apply the optimizer fmincon (MATLAB version R2015a with standard
settings) with initial value τ0= 5 to recover the original delay τ=πup to 8 decimal places.
This example fortifies the approach introduced in section 6.
Example 7.3 We test our approach with the delay model from [5] given by the n×nmatrices
E=θIn+T, A1=1
τ(1
ζ+ 1)(T−θIn), A2=1
τ(1
ζ−1)(T−θIn),
where Tis an n×nmatrix with ones on the sub- and superdiagonal, in the (1,1), and in the
(n, n) position and zeros everywhere else. We choose n= 500, τ= 1, ζ= 0.01, and θ= 5. The
input matrix B∈Rnhas ones in the first two components and zeros everywhere else and we
choose C=BT. Note that the matrices A1and A2satisfy the relation A2= (1
ζ+1)/(1
ζ−1)A1
and hence feature the structure A2=αE +βA1. We pick ρ= 10 random interpolation points
λkon the imaginary axis between 10−1ıand 103ıtogether with their complex conjugates µk.
The realizations based on the data from k= 1,2 (i. e., κ= 2 in (6.1)) and exact parameters
τ,α, and βare displayed in Figure 7.1. Clearly, only our framework captures the transfer
function accurately. In contrast, the Loewner realization (HL) and the realization based on
transformed data (HT) yield poor approximations that do not capture one single peak of the
original model.
Next, the test data is used to find estimates for the parameters τ,α, and βfrom the data
by means of minimizing the cost functional (6.2) over the parameter set. The optimization is
performed via the global optimization algorithm particleswarm of the global optimization
toolbox [15]. This algorithm evaluates the cost functional in (6.2) at a collection of sample
points called particles at each steps. After this evaluation, every particle is moved in the
parameter range and reevaluated. This process continues until some stopping criterion is
satisfied. In this case, the algorithm terminated with
τκ,ρ = 0.999367, ακ,ρ = 0.652009,and βκ,ρ = 0.995628,(7.1)
i. e., the delay is captured accurately.
16
(a) Bode plot of H,HL,HD, and HT.
(b) Error plot for HL,HD, and HT.
Figure 7.1: Example 7.3 - Realizations for κ= 2 and exact parameters τ,α, and β.
We illustrate the transfer functions of the realizations based on all interpolation points and
the parameters (7.1) in Figure 7.2. As before, the realization obtained with our methodology
is a very good approximation of the original model. It is worth to emphasize that this real-
ization is purely data-driven, i. e., is constructed from data only. Again, only our realization
approximates the qualitative behavior of the original model reasonably well.
17
Example 7.4 A heated rod with distributed control and homogeneous Dirichlet boundary
conditions, which is cooled by delayed feedback, can be modeled (cf. [7, 17]) via the one-
dimensional partial differential equation
∂v(ξ, t)
∂t =∂2v(ξ, t)
∂ξ2+a1(ξ)v(ξ, t) + a2(ξ)v(ξ, t −1) + u(t) in (0, π)×(0, T],
v(0, t) = v(π, t) = 0 in [0, T].
(7.2)
Discretization of (7.2) via centered finite differences with step size h=π
n+1 yields a system
˙x(t)=(Ln+A1,n)x(t) + A2,nx(t−1) + Bu(t),
y(t) = Cx(t),
where Ln∈Rn,n is the discrete Laplacian and A1,n, A2,n ∈Rn,n are discrete approximations
of the functions a1and a2, respectively. The input matrix B∈Rnis a vector of ones. As
output we use the average temperature of the rod, i. e, C=1
∥B∥BT. For our tests we use
n= 100, κ= 3, and ρ= 8 random interpolation points λkon the imaginary axis between
10−1ıand 102ıtogether with their complex conjugates µk.
We distinguish two cases. First, we set a1≡1≡a2, which yields E=A2, i. e., resulting
in the structure A2=αE +βA1. The computed estimates for τ,α, and βare
τκ,ρ = 1.000196, ακ,ρ = 1.000384, βκ,ρ =−0.002923 (7.3)
and the results are illustrated in Figure 7.3. As before, the generalized Loewner realization
outperforms the other approaches. In the original model from [7, 17], the coefficients a1and
a2are chosen as
a1(ξ) = −2 sin(ξ) and a2(ξ) = 2 sin(ξ).
Note that in this case A2is not a linear combination of Eand A1. Comparing the results
depicted in Figure 7.4, the generalized delay Loewner approach still captures the qualitative
behavior of the full model and is the best approximation in terms of the maximal error within
the considered range of frequencies.
It is worth to note that the performance of the three approaches depends on the choice of
interpolation points. For each method we found at least one selection of interpolation points,
where this method outperforms the others. Qualitatively, only our ansatz captures the main
features of the full model for a small number of interpolation points in all our tests.
19
(a) Bode plot of H,HL,HD, and HT.
(b) Error plot for HL,HD, and HT.
Figure 7.4: Example 7.4 - Realizations for a1(ξ) = −2 sin(ξ), a2(ξ) = 2 sin(ξ) and recovered
τκ,ρ = 1.006960.
8 Conclusions
We have extended the Loewner realization [16] and the moment matching framework [4] to
descriptor systems with internal delay. We used the results obtained for moment matching
21
to construct a realization directly from measurements of the transfer function. The system
matrices can be assembled efficiently by matrix-matrix operations of size ρ×ρ, where 2ρ
is the number of interpolation points. The internal delay is estimated by solving a least-
square optimization over some sample data. Examples show that our approach produces a
low-order model that captures the dynamics of the full model very accurately. Also, the delay
parameter τis recovered almost exactly. Comparing our ansatz with the common Loewner
realization and an extension introduced in [18] reveals the necessity for preserving the delay
structure. Consequently, our approach yields better approximations of the transfer function,
in particular if the state dimension of the realization is small. Open problems are the optimal
choice of the interpolation points λiand µiand estimators for the output error ∥y−yρ∥,
where yand yρdenote the outputs of the original model and the realization, respectively. In
addition, the proportional ansatz A2,ρ =αEρ+βA1,ρ seems rather restrictive compared to
the overall degrees of freedom and further approaches are to be investigated. Moreover, an
extension to multiple delays and more general structures of the transfer function is currently
under investigation.
References
[1] A. C. Antoulas. Approximation of large-scale dynamical systems. Advances in Design
and Control. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA,
2005.
[2] A. C. Antoulas, S. Lefteriu, and A. C. Ionita. A tutorial introduction to the Loewner
framework for model reduction. In P. Benner, A. Cohen, M. Ohlberger, and K. Will-
cox, editors, Model Reduction and Approximation for Complex Systems, ISNM Series.
Birkh¨auser, 2015.
[3] A. C. Antoulas, D. C. Sorensen, and S. Gugercin. A survey of model reduction methods
for large-scale systems. Contemp. Math., 280:193–219, 2001.
[4] A. Astolfi. Model reduction by moment matching for linear and nonlinear systems. IEEE
Trans. Automat. Control, 55(10):2321–2336, 2010.
[5] C. Beattie and S. Gugercin. Interpolatory projection methods for structure-preserving
model reduction. Syst. Control Lett., 58(3):225–232, 2009.
[6] A. Bellen and M. Zennaro. Numerical methods for delay differential equations. Oxford
University Press, 2003.
[7] D. Breda, S. Maset, and R. Vermiglio. Numerical approximation of characteristic values
of partial retarded functional differential equations. Numer. Math., 113(2):181–242, 2009.
[8] R. W. Clough and J. Penzien. Dynamics of Structures. McGraw-Hill, third edition, 2003.
[9] T. Erneux. Applied Delay Differential Equations. Surveys and Tutorials in the Applied
Mathematical Sciences. Springer New York, 2009.
[10] K. Gallivan, A. Vandendorpe, and P. Van Dooren. Model reduction of MIMO systems
via tangential interpolation. SIAM J. Matrix Anal. Appl., 26(2):328–349, 2004.
22
[11] G. H. Golub and C. F. Van Loan. Matrix Computations. Johns Hopkins Studies in the
Mathematical Sciences. Johns Hopkins University Press, Baltimore, Maryland, fourth
edition, 2013.
[12] P. Ha and V. Mehrmann. Analysis and numerical solution of linear delay differential-
algebraic equations. Preprint 42–2014, Institut f¨ur Mathematik, TU Berlin, 2014.
http://www.math.tu-berlin.de/preprints/. Accepted for publication in BIT, 2015.
[13] J. K. Hale and S. M. Verduyn Lunel. Introduction to Functional Differential Equations.
Number 99 in Applied Mathematical Sciences. Springer, New York, 1993.
[14] P. Kunkel and V. Mehrmann. Differential-Algebraic Equations. Analysis and Numerical
Solution. European Mathematical Society, 2006.
[15] MATLAB. Global Optimization Toolbox User’s Guide, Version R2015a. Natick, Mas-
sachusetts, 2015.
[16] A. J. Mayo and A. C. Antoulas. A framework for the solution of the generalized realization
problem. Linear Algebra Appl., 425(2-3):634–662, 2007.
[17] W. Michiels, E. Jarlebring, and K. Meerbergen. Krylov-based model order reduction of
time-delay systems. SIAM J. Matrix Anal. Appl., 32(4):1399–1421, 2011.
[18] I. Pontes Duff, C. Poussot-Vassal, and C. Seren. Realization independent single time-
delay dynamical model interpolation and H2-optimal approximation, 2015. arXiv
Preprint: 1504.06457v1.
[19] G. Scarciotti and A. Astolfi. Model reduction by moment matching for linear time-delay
systems. In Proceedings of the 19th IFAC World Congress, pages 3642 – 3647, 2014.
[20] M. J. Tov´ee. Neuronal processing. How fast is the speed of thought? Curr. Biol.,
4(12):1125–1127, 1994.
[21] X. Yu and M. Gen. Introduction to Evolutionary Algorithms. Decision Engineering.
Springer London, 2010.
23
