Model reduction of linear time-varying systems with
applications for moving loads
Maria Cruz Varona† and Boris Lohmann†
Preprint: June 15, 2016

arXiv:1607.02846v1 [math.DS] 11 Jul 2016

Abstract
In this paper we consider different model reduction techniques for systems with moving
loads. Due to the time-dependency of the input and output matrices, the application of timevarying projection matrices for the reduction offers new degrees of freedom, which also come
along with some challenges. This paper deals with both simple methods for the reduction of
particular linear time-varying systems, as well as with a more advanced technique considering
the emerging time derivatives.
Keywords: model reduction; time-varying systems; moving loads

1

Introduction

The detailed modeling of physical and technical phenomena arising in many engineering and
computer science applications may yield models of very large dimension. This is particular
the case in fields such as thermo-fluid dynamics, structural mechanics or integrated circuit
design, where the models are mostly obtained from a spatial discretization of the underlying
partial differential equations. The resulting large systems of ordinary differential equations or
differential-algebraic equations are computationally expensive to simulate and handle. In order
to reduce the computational effort, model reduction techniques that generate reduced-order
models that approximate the dynamic behaviour and preserve the relevant properties of the
original model are required. For the reduction of linear time-invariant (LTI) systems, various
well-established reduction approaches exist (see e.g. [1]). In the past ten years, further model
reduction methods have been developed for linear, parametric and nonlinear systems [4, 3] and
applied in a wide variety of domains.
In this contribution, we investigate model order reduction of linear time-varying (LTV)
systems. Such systems arise in many real-life applications, since dynamical systems often depend
on parameters which vary over time or might alter their behaviour due to ageing, degradation,
environmental changes and time-dependent operating conditions. Another possible application
for LTV systems are moving loads. This particular but still very frequent problem arises, for
example, in working gears, cableways, bridges with moving vehicles or milling processes. Since
the position of the acting force varies over time, systems with sliding components exhibit a timevariant behaviour. The varying load location can be modeled and considered in different ways,
thus yielding diverse alternative representations for systems with moving loads and, according
to this, leading to different approaches to reduce them.
One possibility is to represent moving loads as LTV systems, in which only the input and/or
output matrices are time-dependent. Such systems can be then reduced using balanced truncation model reduction methods developed in [16, 15]. These approaches, however, require a
high computational and storage effort, since two differential Lyapunov equations must be solved.
†
Chair of Automatic Control, Technical University of Munich, Boltzmannstr.
({maria.cruz,lohmann}@tum.de)

1

15, D-85748 Garching

Recently, a practical and efficient procedure of balanced truncation for LTV systems has been
presented in [12]. Note that these aforementioned balanced truncation techniques can be applied
to general LTV systems, where all system matrices are time-dependent. For the reduction of
systems with only time-varying input and output matrices the two-step approach proposed in
[17, 2] can also be pursued. This method consists first on a low-rank approximation of the timedependent input matrix and consequently on applying standard model reduction techniques to
the resulting LTI system with a modified input. The approximation of the input matrix in a
low-dimensional subspace is performed via the solution of a linear least squares minimization
problem.
Systems with moving loads can further be modeled by means of linear switched systems.
Well-known reduction methods such as balanced truncation can then be applied for the reduction
of each LTI subsystem [11].
A last alternative option for describing systems with moving loads is to consider the load
position as a time-dependent parameter of the system model. This results in a linear parametervarying (LPV) system, in which only the input and/or output matrices depend on a time-varying
parameter. In many recent publications, e.g. [8, 11, 9, 2], the parameter is assumed to be timeindependent. Thereby any parametric model order reduction (pMOR) approach [4] can be
applied to the resulting parametric LTI system. In some other recent publications [18, 5, 6, 7]
the time variation of the parameter is taken into account, whereby new time derivative terms
emerge during the time-dependent parametric model reduction process.
This paper deals with different time-varying model reduction techniques for systems with
moving loads. Firstly, LTV systems are considered and the time-dependent projective reduction
framework is briefly explained in section 2. Since moving loads represent particular LTV systems,
we then introduce some straightforward reduction approaches for the resulting special cases in
section 3. In the second part of the paper, we focus on LPV systems and present a timedependent parametric model reduction approach by matrix interpolation [5, 7] in section 4.
Some numerical results for the reduction of systems with moving loads applying the proposed
methods are finally reported and discussed in section 5.

2

Linear Time-Varying Model Order Reduction

In the following we first consider a high-dimensional linear time-varying system of the form
E(t) ẋ(t) = A(t) x(t) + B(t) u(t),
y(t) = C(t) x(t),

(1)

where E(t), A(t) ∈ Rn×n , B(t) ∈ Rn×m and C(t) ∈ Rq×n are the time-dependent system
matrices, x(t) ∈ Rn is the state vector and u(t) ∈ Rm , y(t) ∈ Rq represent the inputs and
outputs of the system, respectively. The system matrix E(t) is assumed to be nonsingular for all
t ∈ [0, T ]. Note that it is also possible to consider second-order systems and reformulate them
into the first-order form (1).
2.1

Time-dependent Projective Reduction Framework

In projective model order reduction, we aim to find a reduced-order model by approximating
the state vector x(t) on a subspace of lower dimension r  n. In the time-varying case, the
state vector x(t) might be projected onto a varying subspace spanned by the columns of a timedependent projection matrix V(t) ∈ Rn×r [16, 18]. Therefore, the approximation equations

2

read

x(t) ≈ V(t) xr (t),
ẋ(t) ≈ V̇(t) xr (t) + V(t) ẋr (t),

(2)

whereby the product rule must be considered in this case for the time derivative of the state
vector. Plugging first these both equations into (1), and applying thereon a properly chosen
time-dependent projection matrix W(t) which enforces the Petrov-Galerkin condition leads to
the time-varying reduced-order model


Ar (t)
Br (t)
Er (t)
}|
{
z
}|
{
z
}|
{
z


W(t)T E(t)V(t) ẋr (t) = W(t)T A(t)V(t) −W(t)T E(t)V̇(t) xr (t) + W(t)T B(t) u(t),
(3)
yr (t) = C(t)V(t) xr (t).
| {z }
Cr (t)

It is noteworthy to mention that the system matrix of the reduced-order model (3) not only
comprises the reduced matrix Ar (t), but also includes a further term which depends on the
time derivative V̇(t) of the time-varying projection matrix. This additional term influences the
dynamic behaviour of the reduced-order model and should therefore be taken into account.
The usage of time-dependent projection matrices for the reduction of linear time-varying
systems certainly offers benefits regarding the approximation quality. For their computation,
however, standard reduction methods such as balanced truncation cannot be directly applied,
but must be adapted instead. Furthermore, the time derivative of V(t) should be approximated
numerically (thus increasing the computational effort) and included in the time integration
scheme of the reduced-order model [12].

3

Straightforward Reduction Approaches for particular Linear Time-Varying
Systems

In the previous section we have seen that the application of time-dependent projection matrices
for the reduction of LTV systems comes along with some difficulties and challenges. For the
reduction of particular LTV systems, in which only the input and/or output matrices depend on
time, the usage of time-independent projection matrices V and W might be sufficient. In this
section we discuss some special cases for LTV systems and propose straightforward approaches
to reduce them.
3.1

Case 1: Moving Loads

The first case we want to consider is a high-dimensional LTV system with only time-varying
input matrix, and all other matrices being time-independent:
E ẋ(t) = A x(t) + B(t) u(t),
y(t) = C x(t).

(4)

The time-dependent input matrix describes the position of the moving forces at time t. In the
following we present two straightforward approaches to reduce a system in the form above using
time-independent projection matrices.
Approach 1: Two-step method
The first straightforward reduction method is deducted from the two-step approach presented
in [17, 2]. The method is composed of two steps:
3

1. The time-variability of the input matrix is shifted to the input variables through a lowrank approximation of the input matrix by B(t) ≈ B B̃(t), where B ∈ Rn×m̃ with m̃  n
is a constant matrix and B̃(t) ∈ Rm̃×m . Introducing a new input ũ(t) = B̃(t) u(t), the
original model (4) can be transformed to:
ũ(t)

z }| {
E ẋ(t) = A x(t) + B B̃(t) u(t),

(5)

y(t) = C x(t).
2. The resulting multiple-input multiple-output (MIMO) LTI system (E, A, B, C) can subsequently be reduced by means of balanced truncation, MIMO rational Krylov or MIMOIRKA, for instance. The reduced-order model is then given by
E

A

B

ũ(t)

r
r
r
z }| {
z }| {
z }| { z }| {
T
T
T
W EV ẋr (t) = W AV xr (t) + W B B̃(t) u(t),

(6)

yr (t) = |{z}
CV xr (t),
Cr

where the reduced time-varying input matrix reads Br (t) = Br B̃(t).
For the approximation of the input matrix B(t), other than [17, 2] we simply take the correct
input columns bi (t) with the moving load acting at corresponding nodes i of a coarse finite
element grid and form the low-rank matrix B with them, without performing a least squares
minimization with the basis functions. Note that the two-step approach only provides satisfactory results, if the number of columns m̃ of the low-rank matrix B is sufficiently large [17].
Otherwise the overall approximation error in the output (due to the approximation error in the
input matrix and the model reduction error) can become inadmissibly large. Note also that
this reduction method is limited to systems with a known trajectory of the load before the
simulation.
Approach 2: One-sided reduction with output Krylov subspace
The second straightforward method uses Krylov subspaces for the reduction and exploits the
fact that the only time-varying element in system (4) is the input matrix B(t). Since an input
Krylov subspace would yield a time-varying projection matrix


−1
−1
−1
r−1 A−1 B(t) ,
V(t) := A−1
(7)
s0 B(t) As0 EAs0 B(t) . . . (As0 E)
s0
where As0 = A−s0 E, the idea of this approach is to perform a one-sided reduction with V = W,
where the columns of W form a basis of the output Krylov subspace:


T A−T ET A−T CT . . . (A−T ET )r−1 A−T CT .
W := A−T
(8)
s0 C
s0
s0
s0
s0
Thereby, time-independent projection matrices are obtained for computing the reduced-order
model
Br (t)
Ar
Er
z }| {
z }| {
z }| {
WT EW ẋr (t) = WT AW xr (t) + WT B(t) u(t),
(9)
yr (t) = CW
x
(t).
r
| {z }
Cr

Although only the first r Taylor coefficients (so-called moments) of the transfer function of the
original and the reduced model around the expansion points s0 match due to the application of
a one-sided reduction, we obtain time-independent projection matrices with this approach and
can therefore get rid of the time derivative V̇(t).
4

3.2

Case 2: Moving Sensors

Now we consider a LTV system with only time-varying output matrix
E ẋ(t) = A x(t) + B u(t),

(10)

y(t) = C(t) x(t).

The time-dependent output matrix describes the position of the moving sensors at time t. This
particular LTV system can easily be reduced in the following ways.
Approach 1: Two-step method
1. We shift the time-variability of the output matrix to the output variables through a lowrank approximation by C(t) ≈ C̃(t) C, where C ∈ Rq̃×n with q̃  n is a constant matrix
and C̃(t) ∈ Rq×q̃ . Introducing a new output ỹ(t) = C x(t), the original model (10) can be
transformed to:
E ẋ(t) = A x(t) + B u(t),
(11)

y(t) = C̃(t) C x(t) .
| {z }
ỹ(t)

2. The resulting system (E, A, B, C) can subsequently be reduced by means of any appropriate multiple-input multiple-output LTI reduction technique. The calculated timeindependent projection matrices lead to the reduced-order model
Er

Ar

Br

z }| {
z }| {
z }| {
WT EV ẋr (t) = WT AV xr (t) + WT B u(t),

(12)

yr (t) = C̃(t) |{z}
CV xr (t),
Cr

with the reduced time-varying output matrix Cr (t) = C̃(t) Cr .
The approximation of C(t) is performed by simply taking the output rows with the moving
sensor at the corresponding nodes of a coarse finite element grid. Note that the approximation
of the output matrix yields additional errors in the output [17].
Approach 2: One-sided reduction with input Krylov subspace
Since in this case an output Krylov subspace would lead to a time-varying projection matrix


T A−T ET A−T C(t)T . . . (A−T ET )r−1 A−T C(t)T
W(t) := A−T
(13)
s0 C(t)
s0
s0
s0
s0
due to the time-dependent output matrix C(t), the idea is now to perform a one-sided reduction
with W = V, where the columns of V form a basis of the input Krylov subspace:


−1
−1
−1
r−1 A−1 B .
V := A−1
(14)
s0 B As0 EAs0 B . . . (As0 E)
s0
The reduced model is then given by:
Er

Ar

Br

z }| {
z }| {
z }| {
VT EV ẋr (t) = VT AV xr (t) + VT B u(t),
yr (t) = C(t)V xr (t).
| {z }

(15)

Cr (t)

Due to the application of a one-sided reduction, only r moments are matched. Nevertheless, the
time derivative is avoided, since V and W are time-independent (V̇ = 0).
5

3.3

Case 3: Moving Loads and Sensors

Finally, we consider the combined case with time-varying input and output matrices
E ẋ(t) = A x(t) + B(t) u(t),

(16)

y(t) = C(t) x(t).
If the sensor position coincides with the location of the load, then C(t) = B(t)T .
Approach 1: Two-step method

In this case, the respective two-step techniques explained before have to be combined properly:
1. The time-variability of B(t) is shifted to the input variables and the time-dependency of
C(t) to the output variables, thus obtaining a MIMO LTI system.
2. Time-independent projection matrices are then calculated applying an appropriate model
order reduction method to the resulting system (E, A, B, C). The reduced-order model is
finally given by
E

A

B

ũ(t)

r
r
r
z }| {
z }| {
z }| { z }| {
T
T
T
W EV ẋr (t) = W AV xr (t) + W B B̃(t) u(t),

(17)

yr (t) = C̃(t) |{z}
CV xr (t).
Cr

Approach 2: Reduction with modal truncation
Unfortunately, in this case 3 the application of a one-sided reduction with either an input or an
output Krylov subspace would yield time-varying projection matrices V(t) and W(t) according
to (7) or (13), respectively. A possible alternative to still obtain time-independent projection
matrices and thus get rid of the time derivative V̇ is to use modal truncation as reduction approach. This method only uses the time-independent matrices A and E for computing dominant
eigenvalues (e.g. with smallest magnitude or smallest real part) and eigenvectors, thus yielding
time-independent projection matrices for the reduction.

4

Time-Varying Parametric Model Order Reduction

After having considered linear time-varying systems and presented some straightforward approaches to reduce special cases arising in moving load and sensor problems, in this section we
focus on linear parameter-varying systems of the form
E(p(t)) ẋ(t) = A(p(t)) x(t) + B(p(t)) u(t),
y(t) = C(p(t)) x(t).

(18)

Such systems also exhibit a time-varying dynamic behaviour, since the system matrices explicitly
depend on parameters p(t) which vary over time. Note that moving load and sensor problems
can be represented as LPV systems with only parameter-varying input and/or output matrices,
if the load and sensor location are considered as time-dependent parameters of the system model.
Due to the time-dependency of the parameters, in the next subsection we derive a projectionbased, time-varying parametric model order reduction approach called p(t)MOR, to obtain a
reduced-order model of a LPV system [5, 6]. Based on that, we then adapt the pMOR approach
by matrix interpolation [14] to the parameter-varying case, whereby new time derivative terms
emerge [6, 7]. For the sake of a concise presentation, the time argument t will be omitted in the
state, input and output vectors hereafter.
6

4.1

Projective p(t)MOR

Similarly as explained in subsection 2.1, in the case of projection-based time-dependent parametric model order reduction we aim to approximate the state vector x by x ≈ V(p(t)) xr using a
parameter-varying projection matrix V(p(t)). Plugging the corresponding approximation equations for x and its derivative ẋ in (18), and applying thereon a properly chosen projection matrix
W(p(t)) that imposes the Petrov-Galerkin condition yields the reduced-order model


Er (p(t)) ẋr = Ar (p(t)) − W(p(t))T E(p(t))V̇(p(t)) xr + Br (p(t)) u,
(19)
yr = Cr (p(t)) xr ,
with the time-dependent parametric reduced matrices

Er (p(t)) = W(p(t))T E(p(t))V(p(t)), Ar (p(t)) = W(p(t))T A(p(t))V(p(t)),
Br (p(t)) = W(p(t))T B(p(t)),

(20)

Cr (p(t)) = C(p(t))V(p(t)).

The reduced model comprises an additional term depending on the time derivative V̇(p(t)),
which has to be considered during the extension of the matrix interpolation method to the
parameter-varying case.
4.2

p(t)MOR by Matrix Interpolation

The local pMOR technique of matrix interpolation can be applied to efficiently obtain a parametric reduced-order model from the interpolation of reduced matrices precomputed at different
grid points in the parameter space. Similarly as in the classic method [14], the LPV system (18)
is first evaluated and individually reduced at certain parameter samples pi , i = 1, . . . , k with
respective projection matrices Vi := V(pi ) and Wi := W(pi ). The reduced state vectors xr,i
of the independently calculated reduced models


Er,i ẋr,i = Ar,i − WiT Ei V̇(p(t)) xr,i + Br,i u,
(21)
yr,i = Cr,i xr,i
generally lie in different subspaces and have, therefore, different physical meanings. For this
reason, the direct interpolation of the reduced matrices is not meaningful, and hence the local
reduced models have to be transformed into a common set of coordinates first. This is performed
applying state transformations of the form
xr,i = Ti x̂r,i ,
ẋr,i = Ṫi x̂r,i + Ti x̂˙ r,i ,

(22)

with regular matrices Ti := T(pi ), whereby the product rule is required again for the differentiation of xr,i . These state transformations serve to adjust the different right local bases Vi
to new bases V̂i = Vi Ti . In order to adjust the different left local bases Wi by means of
Ŵi = Wi Mi as well, the reduced models from (21) are subsequently multiplied from the left

7

with regular matrices MTi . The resulting reduced and transformed systems are thus given by
Ânew r,i

z
}|
{
Êr,i
Âr,i
B̂r,i
z }| {
z }| {
z }| {


T
T
T
T
T

MTi Er,i Ti x̂˙ r,i = 
Mi Ar,i Ti −Mi Wi Ei V̇(p(t))Ti − Mi Er,i Ṫi  x̂r,i + Mi Br,i u,

(23)

yr,i = Cr,i Ti x̂r,i .
| {z }
Ĉr,i

One possible way to calculate the transformation matrices Ti and Mi is based on making
the state vectors x̂r,i compatible with respect to a reference subspace spanned by the columns
of the orthogonal matrix R. To this end, the matrices are chosen as Ti := (RT Vi )−1 and
Mi := (RT Wi )−1 , where the columns of R correspond to the r most important directions of
Vall = [V1 . . . Vk ] calculated by a Singular Value Decomposition (SVD) [14].
The resulting system matrix Ânew r,i not only comprises the expected reduced matrix Âr,i ,
but also consists of two further terms that depend on V̇(p(t)) and Ṫi , respectively. The calculation of these time derivatives that are required for the computation of the reduced-order model
will be discussed in the next two sections.
After the transformation of the local models and the computation of the new emerging time
derivatives, a parameter-varying reduced-order model for a new parameter value p(t) is obtained
in the online phase by a weighted interpolation between the reduced matrices from (23) according
to
Xk
Xk
Ẽr (p(t)) =
ωi (p(t))Êr,i , Ãnew r (p(t)) =
ωi (p(t))Ânew r,i ,
i=1
i=1
(24)
Xk
Xk
B̃r (p(t)) =
ωi (p(t))B̂r,i ,
C̃r (p(t)) =
ωi (p(t))Ĉr,i ,
i=1
i=1
Pk
where i=1 ωi (p(t)) = 1. For simplicity, here we use piecewise linear interpolation of the reduced
matrices. Higher order interpolation schemes could also be applied.
4.2.1

Time derivative of V

The time derivative of the projection matrix V(p(t)) can be numerically calculated using a finite
difference approximation. Applying the chain rule first and employing a finite difference method
thereon, the time derivative is given by:
V̇(p(t)) =

∂V
Vt − Vt pt − pt−1
ṗ =
.
∂p
pt − pt
∆t

(25)

pt and pt denote the upper and lower limit of the interval [pt , pt ], in which the parameter vector
pt is located at time instant t. The local bases at these parameter sample points are given by
Vt and Vt , respectively. The partial derivatives ∂V
∂p for each pair of parameter sample points
are calculated in the offline phase of the matrix interpolation approach. In the online phase,
the current time derivative V̇(p(t)) is then computed by multiplying the partial derivative of
the corresponding parameter interval at time instant t with ṗ, which represents the current
velocity of the moving load. Fig. 1 illustrates the aforementioned intervals and the efficient
numerical calculation of the time derivative V̇(p(t)) by a finite difference approximation using
only precomputed local bases.

8

Vt−1

Vt−1
pt−1

Vt

Vt

Vt−1
pt−1
∂V
∂p

pt−1

∂V
∂p

pt

Vt
pt
∂V
∂p

pt

p

Figure 1: Graphical representation of the calculation of the time derivative V̇(p(t)) using the
local bases Vi computed at the parameter sample points pi
4.2.2

Time derivative of T

As explained before, in this paper the transformation matrices Ti are calculated with Ti =
(RT Vi )−1 := K−1 . For the computation of the time derivative Ṫi we make use of the following
definition [10, p. 67]:
Definition 1. Let the matrix K be nonsingular. The time derivative of the inverse matrix is
−1
−1
then given by dKdt = −K−1 dK
dt K .
This leads to:
Ṫi =
4.3

dK−1
= −(RT Vi )−1 RT V̇(p(t))(RT Vi )−1 = −Ti RT V̇(p(t))Ti .
dt

(26)

p(t)MOR by Matrix Interpolation for particular cases

For the reduction of general linear parameter-varying systems the application of time-dependent
parametric projection matrices undoubtedly provides an accurate consideration of the arising
time variability. Their usage, however, involves some difficulties, like the calculation of the
additional derivatives and their incorporation in the numerical simulation of the reduced-order
model. Particular LPV systems with only parameter-varying input and/or output matrices,
arising e.g. in moving load and sensor problems, can efficiently be reduced using the matrix
interpolation approach combined with the usage of parameter-independent projection matrices.
In the following, this technique is briefly explained for some special cases:
Moving Loads
The application of parameter-varying projection matrices V(p(t)) and W(p(t)) for the individual reduction of the local systems within matrix interpolation results in a reduced model, where
all reduced matrices vary with the time-dependent parameter, although the original LPV system
only contains variations in the input matrix. In order to get rid of the emerging derivatives and
only have to interpolate the input matrix in the online phase of matrix interpolation, one-sided
reductions with an output Krylov subspace W = span(W) should be employed.
Moving Sensors
In a similar manner, for the case of a LPV system with only parameter-varying output matrix
C(p(t)) one-sided projections with a single input Krylov subspace V = span(V) computed with
the input matrix should be performed for the reduction of the sampled models during matrix
interpolation. In this way, we obtain parameter-independent projection matrices V = W and
only have to interpolate the output matrix, thus reducing the computational effort in the online
phase.
9

Moving Loads and Sensors
For the combined moving load and sensor example the application of one-sided projections with
either input or output Krylov subspaces is not helpful, since both the input and output matrices
are parameter-varying in this case. Therefore, parameter-independent projection matrices can
only be calculated using modal truncation. By doing so, the reduced-order model only contains
parameter variations in the input and output reduced matrices like in the original system.

5

Numerical Examples

In this section, we present some numerical results for systems with moving loads.
5.1

Timoshenko beam

The presented reduction approaches are first applied to the finite element model of a simply
supported Timoshenko beam of lenght L subjected to a moving load.

F(t)
z

y
x

Figure 2: A simply supported Timoshenko beam subjected to a moving force F (t)
Since the moving force F (t) is applied in the negative z-direction and we are only interested in
the vertical displacement of the beam, the model described in [13] is adapted from a 3D to a 1D
finite element model. Furthermore, both the moving load and/or sensor case are incorporated
into the model, yielding time-dependent input and/or output matrices. The resulting singleinput single-output second-order system is reformulated into a LTV first-order model of the
form
b(t)
E
A
ẋ
x
z }| {
z }| { z}|{
z
}|
 

{ z}|{
 
0
F 0
ż
0
F
z
(t) =
(t) +
F (t),
0 M z̈
−K −D ż
b̂(t)
(27)
 


z
y(t) = ĉ(t)T 0T
(t),
|
{z
} ż
c(t)T

where the arbitrary nonsingular matrix F ∈ R2N ×2N is chosen in our case to F = K for the aim
of stability preservation using a one-sided reduction (cf. [13, 6]). The dimension of the original
model is then given by n = 2 · 2N with N finite elements.
Moving load case
We first consider the reduction of a beam of length L = 1 m subjected to a point force moving
from the tip to the supporting with a constant velocity v and an amplitude of F (t) = 20 N.
For the numerical simulation we use an implicit Euler scheme with a step size of dt = 0.001 s.
In Fig. 3 the simulation results for the different proposed reduction methods are presented.
We first apply the standard matrix interpolation (MatrInt) approach using k = 76 equidistantly
10

−2
−4
−6
0.2

·10−2

0.4
0.6
time [s]

0.8

FOM
ROM MatrInt V(p(t)) 1-RK
ROM MatrInt V̇, Ṫ 1-RK
ROM MatrInt W 1-RK
ROM approx B + IRKA

−4

0

z-deflection of beam tip [m]

−2
−4
0.4
0.6
time [s]

0.2

·10−2

0

0.2

−2

v = 10 ms

2

0

0

−6

1

v = 5 ms

·10−2

2

0

0

z-deflection of beam tip [m]

v = 1 ms

z-deflection of beam tip [m]

z-deflection of beam tip [m]

·10−2

0.8

1

0.4
0.6
time [s]

0.8

1

0.8

1

v = 20 ms

2

0

−2
0

0.2

0.4
0.6
time [s]

Figure 3: Simulation results for the Timoshenko beam with moving load for different reduction
methods and velocities. Original dimension n = 2 · 2 · 451 = 1804, reduced dimension r = 10.
Krylov-based reductions performed with expansion points s0 = 0
distributed local models with corresponding current input, which are individually reduced applying one-sided projections with input Krylov subspaces (V(p(t))) for r = 10. The consideration
of the theorically emerging derivatives V̇ and Ṫ according to (23) in the matrix interpolation
scheme only yields better results than the standard MatrInt method for large velocities of the
moving load. In any case, the application of a single time-independent output Krylov subspace
(W) during MatrInt and the two-step method (m̃ = 76) combined with MIMO-IRKA yields the
best results (see Table 1).
Moving load and sensor case
Now we consider a larger beam of length L = 50 m with both moving load and sensor. The
observation of the z-deflection of the beam coincides at any time with the position of the moving
load, meaning that c(p(t))T = b(p(t)). First we apply the matrix interpolation approach and
use modal truncation for the individual reduction of the k = 201 sampled models constructed
with the input and output vectors corresponding to each parameter sample point. Since modal
truncation only considers the matrices A and E for the reduction and these matrices do not vary
over time, we only have to compute one single pair of time-independent projection matrices V
11

Table 1: Absolute L2 output error norms ky − yr kL2
MatrInt V(p(t))

MatrInt V̇, Ṫ

MatrInt W

approx B + IRKA

v = 1m
s

4.8e−4

1.5e−2

3.0e−5

1.0e−4

v = 5m
s

3.0e−3

6.8e−2

1.8e−5

2.0e−4

v = 10 m
s

2.6e−3

1.0e−3

1.7e−5

5.8e−5

v = 20 m
s

2.2e−2

5.4e−3

1.2e−5

3.9e−5

and W in the offline phase. During the online phase, only the parameter-varying input and
output vectors have to be interpolated in order to obtain a reduced-order model for each current
position of the load/sensor.
Next, we further apply the aforementioned two-step method for the reduction. To this end,
the time-varying input and output vectors are first approximated by low-rank matrices B and
C on a coarse finite element grid. To ensure a proper comparability with MatrInt, we choose
the same m̃ = 201 nodes where local models were constructed before. The herewith obtained
approximated output y(t) and approximation errors are depiced in Fig. 4. One can see that the
number of chosen columns m̃ is sufficiently large, since the approximation error is adequately
small. After that, we both apply two-sided MIMO rational Krylov (2-RK) and MIMO-IRKA
for the reduction of the resulting LTI system. Fig. 4 shows the simulated output for the different explained reduction methods as well as the corresponding absolute and relative L2 errors.
Although all results show a similar behaviour, the matrix interpolation approach combined with
modal truncation together with the two-step method by IRKA lead to the smallest errors. Simulations were also conducted with the extended p(t)MOR approach by matrix interpolation
considering the time derivatives like in 23. Unfortunately, these additional terms make the
pencils (Ânew r,i , Êr,i ) often unstable, yielding unstable interpolated systems and results.
5.2

Beam with moving heat source

We now apply the presented techniques on a second example [12], which describes the heat
transfer along a beam of length L with a moving heat source. The temperature is observed
at the same position as the heat source, thus c(t)T = b(t). In our case, we consider a system
dimension of n = 2500, apply an input heating source of u(t) = 50 ◦ C and use an implicit
Euler scheme with a step size of dt = 1 s for the time integration. Fig. 5 shows the simulation
results, and the absolute and relative errors for the different employed reduction methods. One
interesting observation is that in this case the application of the extended MatrInt approach
with the consideration of the time derivatives yields a slightly better approximation than the
classic MatrInt combined with modal truncation (k = 84). In general, this fact could also
be observed for the previous and some other numerical experiments with higher velocities, as
long as the overall interpolated systems were stable. This slightly better approximation can be
explained through the more accurate consideration of the arising time variability using timedependent projection matrices, as opposed to modal truncation which does not consider the
moving interactions. The approximation of the time-dependent input and output vectors by
low-rank matrices using m̃ = 84 nodes, and the subsequent application of balanced truncation
(TBR) or MIMO-IRKA for the reduction shows a similar behaviour. Although the extended
MatrInt shows in this case the best results, it is difficult to clearly identify a superior method,
since all presented approaches are suitable for the reduction of systems with moving loads.

12

z-deflection at current sensor position [m]

v = 11 ms

·10−2

0

−11

FOM
ROM MatrInt V, W MODAL
FOM approx B, C
ROM approx B, C + 2-RK
ROM approx B, C + IRKA

−2
−11.2
−6
−11.4
−10
−11.6
−14
−11.8
−18

0

1

2
3
time [s]

4

3.95

5

v = 11 ms

·10−3

4.05

4
relative error

2

0

4
time [s]
v = 11 ms

·10−2

4
absolute error

v = 11 ms

·10−2

0

1

2
3
time [s]

4

2

0

5

0

1

2
3
time [s]

4

5

Figure 4: Simulation results for the Timoshenko beam with moving load and sensor for different
reduction methods. Original dimension n = 2 · 2 · 1001 = 4004, reduced dimension r = 80.
Krylov-based reductions performed with expansion points s0 = 0

13

temperature at current sensor position [◦ C]

·10−3

5

·10−3

4
4.9
FOM
ROM MatrInt V, W MODAL
ROM MatrInt V̇, Ṫ 1-RK
FOM approx B, C
ROM approx B, C + TBR
ROM approx B, C + IRKA

4.8

2

4.7
0

0

20

40
60
time [s]

80

100

40

·10−4

1.5

45

50
time [s]

20

40
60
time [s]

55

60

·10−2

10

1

relative error

absolute error

8

0.5

6
4
2

0

20

40
60
time [s]

80

100

0

0

80

100

Figure 5: Simulation results for the 1D beam with moving heat source for different reduction
methods. Original dimension n = 2500, reduced dimension r = 40. Krylov-based reductions
performed with expansion points s0 = 0

14

6

Conclusions

In this paper, we have presented several time-varying model reduction techniques for systems
with moving loads. Such particular, but still frequent problems lead to high-dimensional systems
with time-varying input and/or output matrices. For their reduction, time-dependent projection
matrices can be applied, thus offering an accurate consideration of the time variation, but leading
also to an additional derivative in the reduced model which have to be taken into account. Since
moving load problems represent particular LTV systems, we have presented straightforward
reduction approaches for some special cases, where time-independent projection matrices are
calculated and therefore the emerging time derivative is avoided. Systems with moving loads
can also be modeled as special LPV systems, where the input and/or output matrices depend
on a time-varying parameter describing the position of the load. In this context we have derived
a projection-based, time-varying parametric model reduction approach and extended the matrix
interpolation scheme to the parameter-varying case. With the appropriate combination of this
method with the application of parameter-independent projection matrices, special LPV systems
can be efficiently reduced avoiding the time derivatives. The proposed methods have been tested
on two different beam models for both the moving load and/or sensor cases. All techniques have
provided similar satisfactory results, showing that all methods are suitable for the reduction of
systems with moving loads. In particular, the presented straightforward approaches using timeindependent projection matrices are very simple, but may be absolutely sufficient for certain
problems. They provide a basis for comparison with more complex techniques that consider the
time variability using time-dependent projection matrices. These advanced techniques should be
investigated more deeply in the future, especially concerning general LTV systems, the increased
computational effort due to the time-dependent projection matrices and derivatives, fast-varying
load variations and stability preservation in the reduced-order model.
Acknowledgements
The authors would like to thank N. Lang, T. Stykel and J. Saak for kindly proving us the 1D
heat transfer model used for the numerical example in section 5.2. Furthermore, we thank the
former and current members of our model reduction lab for the fruitful discussions.

References
[1] A. C. Antoulas. Approximation of Large-Scale Dynamical Systems. SIAM, Philadelphia,
PA, 2005.
[2] M. Baumann, A. Vasilyev, T. Stykel, and P. Eberhard. Comparison of two model order
reduction methods for elastic multibody systems with moving loads. In Proceedings of the
Institution of Mechanical Engineers, Part K: Journal of Multi-body Dynamics, 2016.
[3] U. Baur, P. Benner, and L. Feng. Model order reduction for linear and nonlinear systems: a system-theoretic perspective. Archives of Computational Methods in Engineering,
21(4):331–358, 2014.
[4] P. Benner, S. Gugercin, and K. Willcox. A survey of projection-based model reduction
methods for parametric dynamical systems. SIAM Review, 57(4):483–531, 2015.
[5] M. Cruz Varona, M. Geuss, and B. Lohmann. p(t)MOR: Time-varying parametric model
order reduction and applications for moving loads. In Proceedings of the 8th Vienna Conference on Mathematical Modelling (MATHMOD), volume 48, pages 677–678. Elsevier, 2015.
15

[6] M. Cruz Varona, M. Geuss, and B. Lohmann. Zeitvariante parametrische Modellordnungsreduktion am Beispiel von Systemen mit wandernder Last. In Methoden und Anwendungen der Regelungstechnik – Erlangen-Münchener Workshops 2013 und 2014, pages
57–70. B. Lohmann und G. Roppenecker (Hrsg.), Shaker-Verlag, 2015.
[7] M. Cruz Varona and B. Lohmann. Time-varying parametric model order reduction by matrix interpolation. In Proceedings of the MoRePaS 2015–Model Reduction of Parametrized
Systems III, 2015.
[8] M. Fischer and P. Eberhard. Application of parametric model reduction with matrix interpolation for simulation of moving loads in elastic multibody systems. Advances in Computational Mathematics, pages 1–24, 2014.
[9] M. Fischer and P. Eberhard. Interpolation-based parametric model order reduction for
material removal in elastic multibody systems. In Proceedings ECCOMAS Thematic Conference on Multibody Dynamics, 2015.
[10] G. H. Golub and C. F. Van Loan. Matrix Computations. Johns Hopkins University Press,
Baltimore, 4th edition, 2013.
[11] N. Lang, J. Saak, and P. Benner. Model order reduction for systems with moving loads.
at-Automatisierungstechnik, 62(7):512–522, June 2014.
[12] N. Lang, J. Saak, and T. Stykel. Balanced truncation model reduction for linear timevarying systems. MCMDS, to appear.
[13] H. Panzer, J. Hubele, R. Eid, and B. Lohmann. Generating a parametric finite element model of a 3D cantilever Timoshenko beam using MATLAB. TRAC, Lehrstuhl für
Regelungstechnik, Technische Universität München, 09 2009.
[14] H. Panzer, J. Mohring, R. Eid, and B. Lohmann. Parametric model order reduction by
matrix interpolation. at–Automatisierungstechnik, 58(8):475–484, 8 2010.
[15] H. Sandberg and A. Rantzer. Balanced truncation of linear time-varying systems. IEEE
Transactions on Automatic Control, 49(2):217–229, 2004.
[16] S. Shokoohi, L. M. Silverman, and P. M. van Dooren. Linear time-variable systems: balancing and model reduction. IEEE Transactions on Automatic Control, 28(8):810–822,
1983.
[17] T. Stykel and A. Vasilyev. A two-step model reduction approach for mechanical systems
with moving loads. Journal of Computational and Applied Mathematics, 297:85–97, 2016.
[18] T. Tamarozzi, G. Heirman, and W. Desmet. An on-line time dependent parametric model
order reduction scheme with focus on dynamic stress recovery. Computer Methods in Applied
Mechanics and Engineering, 268(0):336 – 358, 2014.

16

