arXiv:1712.00369v1 [cs.NA] 1 Dec 2017

Reachability Analysis of Large Linear Systems with
Uncertain Inputs in the Krylov Subspace
Matthias Althoff
Abstract. One often wishes for the ability to formally analyze large-scale
systems—typically, however, one can either formally analyze a rather small
system or informally analyze a large-scale system. This work tries to further
close this performance gap for reachability analysis of linear systems. Reachability analysis can capture the whole set of possible solutions of a dynamic
system and is thus used to prove that unsafe states are never reached; this
requires full consideration of arbitrarily varying uncertain inputs, since sensor noise or disturbances usually do not follow any patterns. We use Krylov
methods in this work to compute reachable sets for large-scale linear systems.
While Krylov methods have been used before in reachability analysis, we overcome the previous limitation that inputs must be (piecewise) constant. As a
result, we can compute reachable sets of systems with several thousand state
variables for bounded, but arbitrarily varying inputs, as demonstrated using a
bridge model subject to disturbances.

1. Introduction
Reachability analysis computes the set of possible solutions of dynamic systems
subject to uncertain initial states and inputs. The availability of all possible solutions can be used for many purposes: formal verification of dynamic systems with
discrete and/or continuous dynamics [13,16]; computation of invariance sets [17,41];
computation of the region of attraction [27]; optimization of constrained systems
with uncertainties [46, 55]; set-based observers [21, 45]; and conformance checking [51]. The theory of efficiently computing reachable sets is advancing rapidly and
their usefulness has already been demonstrated for many applications, such as automated driving [4, 28], robotics [49, 57], power systems [26, 61], and analog/mixedsignal circuits [9, 30].
Reachability analysis of linear systems. The reachability analysis in particular of linear continuous systems has been intensely researched. While the main
approach for computing reachable sets of linear systems by using the superposition principle has not significantly changed in recent years [22], much progress
has been made by experimenting with different set representations: ellipsoids [43],
polytopes [20], zonotopes [31], zonotope bundles [6], support functions [32], level
sets [47], and combination of support functions and zonotopes [5]. In particular, when using zonotopes, support functions, or the combination thereof, one can
1

2

MATTHIAS ALTHOFF

efficiently compute systems with more than 100 continuous state variables. Recently, a new technique has been proposed which combines simulation results by
using the superposition principle to represent reachable states via generalized star
sets [15, 25]. Although this technique can compute large systems, it cannot consider uncertain time-varying inputs and requires a formally verified solver for linear
systems. Since one cannot exactly compute the reachable set of linear systems,
except when all eigenvalues are real or imaginary [44], one typically demands tight
over-approximations of reachable sets.
Further use for nonlinear systems. Reachable set computations of linear systems
are often embedded in algorithms for computing reachable sets of nonlinear systems,
either by partitioning the state space into conservatively linearized regions [14] or
by conservatively and continuously linearizing a nonlinear system along its center
trajectory [10, 23]. With conservative linearization, we refer to techniques that
compensate linearization errors by adding uncertainty, e.g., in the form of additive
uncertain inputs.
Order reduction techniques. However, for large-scale linear systems with uncertain time-varying inputs beyond 1000 continuous state variables, even the most
efficient reachability algorithms become too slow for practical use. One of the main
reasons is that the exponential matrix eAt of the system matrix A in the linear system dynamics ẋ = Ax + Bu may be unbearably time-consuming to compute [48].
It should be noted that the dimension is not the only critical parameter; sparsity
of A and the sensitivity of the matrix exponential [48, Sec. 2] are also important
parameters influencing the computation time. For this reason, methods have been
developed to reduce the order of the investigated system, which are generally referred to as order reduction techniques. There exist a vast number of different
techniques surveyed in e.g., [12, 53, 54]. Most order reduction techniques aim at
achieving a similar input/output behavior when the system is initially in a steady
state; they rely on the fact that large-scale systems often have a large number of
state variables but a rather small number of input and output variables. A typical
example is that of infinite-dimensional systems, which are spatially discretized, resulting in high-dimensional differential equations, and controlled by few actuators
and sensed by few sensors. Many applications of reachability analysis, however,
such as formal verification, require that all or many state variables are accurately
approximated as well. For instance, in a large power network, all voltages and
frequencies have to stay within certain bounds. Subsequently, we review previous
work on combining reachability analysis with order reduction techniques.
Order reduction for reachability analysis. To the best knowledge of the author, the first work combining reachability analysis with order reduction techniques
is [37]. There, simulating the solutions of all vertices of an initial set was required
to guarantee an error bound for a set of initial states—this approach is exponential
in the number of state variables when each variable is uncertain within an interval.
At the time the approach in [37] was proposed, tools for reachability analysis had
an exponential complexity as well, so that overall the computation time could be
significantly reduced. However, modern tools such as SpaceEx [29], Flow* [18],
HyLAA [15], XSpeed [36], or CORA [2] have a polynomial complexity, as demonstrated in [3], and thus would most certainly outperform the technique proposed
in [37], even without any order reduction. The same authors later combined reachability analysis with Krylov subspace approximation methods [38]. The advantage

REACHABILITY ANALYSIS WITH UNCERTAIN INPUTS IN THE KRYLOV SUBSPACE

3

of this technique is that it does not scale exponentially in the number of state variables; however, it can only handle linear systems with fixed input, which are also
referred to as affine systems. Recently, the work from [37] was continued in [58],
which considers order reduction in the input/output sense for stable linear systems.
In contrast, this work can also handle unstable systems and reconstruct the whole
set of states, not just the outputs. For nonlinear systems, non-rigorous reduction
techniques have been presented in [19]; however, unlike in this work, the results are
not formal.
Approximate bisimulation. Somewhat related to order reduction techniques is
(approximate) bisimulation [35, 59], which basically shows that two systems have
similar output behavior for the same inputs, while at the same time a relation S
exists so that (x(t), x̃(t)) ∈ S. However, even bisimulation techniques for linear
systems do not scale to the system dimensions handled in this work, since they
require solving linear matrix inequalities [60] to find a bisimulation relation [34].
Contribution. This work extends the state of the art in reachability analysis of
linear systems by proposing the first method with uncertain, time-varying inputs
in the reduced Krylov subspace. This reduction makes it possible to compute
reachable sets of systems with a number of continuous state variables that was
previously infeasible. Please note that input/output order reduction techniques
such as transformation, truncation, and projection (combination of transformation
and truncation) cannot reconstruct the state, in contrast to the reduction presented
in this work. The proposed technique rigorously considers reduction errors for each
consecutive time interval individually. In contrast to [38], our approach a) does
not rely on computing errors by the norm of the system matrix kAk (see (6)),
which would quickly accumulate errors when kAk is greater than 1 (mostly true
in practice), and b) does not require to enlarging the reachable set equally in all
dimensions to account for errors, which can cause large over-approximations. Our
approach is implemented in CORA [2] and will be publicly released with the new
CORA version.
Organization. The paper is organized as follows: Sec. 2 presents preliminaries
from the areas of Krylov subspace approximation, set representation, and reachability analysis. The computation of reachable sets in the Krylov subspace is presented
step by step: The homogeneous solution is described in Sec. 3 and the input solution
in Sec. 4. Then, these are combined in Sec. 5 to demonstrate the overall algorithm.
We close with a detailed numerical example in Sec. 6 and the conclusions in Sec. 7.
2. Preliminaries
Let us first recall some important basics required for this work: Krylov subspace
approximation, representation of continuous sets in high-dimensional spaces, and
computation of reachable sets of linear systems with uncertain inputs.
2.1. Krylov Subspace Approximation. The main obstacle towards reachability analysis of large-scale linear systems is the evaluation of eC v, where C ∈ Rn×n
and v ∈ Rn . To compute eC v more efficiently, we introduce the Krylov subspace
Kξ = span(v, Cv, . . . , C ξ−1 v),
where span(·) returns the linear span of a set of vectors and ξ denotes the dimension
of the subspace. Several possibilities have been developed for approximating eC v in
the Krylov subspace [39, 52, 53]. In this work, we use the simplest approach, which

4

MATTHIAS ALTHOFF

is also one of the most popular: the Arnoldi algorithm, as presented in Alg. 1
(see [52, Sec. 2.1] and [56, Alg. 1]). Please note that w∗ in Alg. 1 denotes the
complex conjugate of w and k.k returns the Euclidean norm. A further reason for
choosing the Arnoldi algorithm is that the tightest a posteriori error bound that
we have found requires this approach [62].
Algorithm 1 Arnoldi iteration
Require: C, v, max order ξ, tolerance tol
Ensure: H, V
1: v (1) = v/kvk
2: for k = 1 . . . ξ do
3:
w = Cv (k)
4:
for j = 1 . . . k do
5:
hj,k = w∗ v (j)
6:
w := w − hj,k v (j)
7:
end for
8:
hk+1,k = kwk
9:
if hk+1,k ≤ tol kCk then
10:
happy-breakdown
11:
end if
12:
v (k+1) = w/hk+1,k
13: end for
14: V = [v (1) , v (2) , . . . , v (ξ) ]
The results of the Arnoldi iteration are an orthogonal basis V = [v (1) , v (2) , . . . , v (ξ) ]
of the Krylov subspace Kξ and the upper ξ × ξ Hessenberg matrix H consisting
of the elements hij from Alg. 1. Please note that we are numbering vectors with
superscripted numbers in parentheses in order to avoid confusion with powers. Using H, V , and e1 = [1, 0, 0, . . . , 0]T , the evaluation of eC v can be approximated
as [52, eq. 3]:
(1)

eC v ≈kvkV eH e1 .

The following rigorous bound for the approximation error exists [52, eq. 29]:
eC v − kvkV eH e1 ≤ 2kvk

kCkξ ekCk
.
ξ!

However, this simple bound is very conservative and may even require choosing
ξ to be larger than the actual system dimension n when a certain threshold for
the approximation error is demanded by the user. For this reason, we use more
sophisticated error computations proposed in [62] for different types of matrices
(skew-Hermetian matrices, positive definite matrices, etc.). In particular, we apply
in our example [62, Corollary 5.3]. For computational reasons, we separate the
error computation of [62] into a part that we refer to as the normalized error ǫnorm
and the norm of the initial value kvk:
(2)

eC v − kvkV eH e1 ≤ kvkǫnorm,

where the computation of ǫnorm is described in detail in Appendix B.

REACHABILITY ANALYSIS WITH UNCERTAIN INPUTS IN THE KRYLOV SUBSPACE

5

2.2. Set Representation. As already summarized in the introduction, several set representations for reachability analysis of linear systems have been proposed. When considering uncertain inputs, zonotopes [31] and support functions
[32] demonstrate the best performance. Recently, [5] showed that combining zonotopes and support functions provides even better benefits: zonotopes can compute
the solution more efficiently, while support functions can represent more general
initial sets. Since most initial sets are multidimensional intervals (special case of
zonotopes), we use zonotopes and neglect their combination with support functions
to focus on the novel aspects of this work.
g

(i)

Definition 2.1 (Zonotope). Given a center c ∈ Rn and so-called generators
∈ Rn , a zonotope is defined as
p
o
n
X
βi g (i) , βi ∈ [−1, 1]
Z := x ∈ Rn x = c +
i=1

We write in short Z = (c, g (1) , . . . , g (p) ) and define the order of a zonotope
as o := np , where p is the number of generators. A zonotope can be seen as the
Minkowski addition of line segments [−1, 1]g (i), which provides an idea of how a
zonotope is constructed; see Fig. 1.

”⊕”

construction
direction

g (1)

”⊕”
g (3)

g (2)

PSfrag replacements
c

c ⊕ g (1)

c ⊕ g (1) ⊕ g (2)

c ⊕ g (1) ⊕ g (2) ⊕ g (3)

Figure 1. Step-by-step construction of a two-dimensional zonotope.
The iterative computation of reachable sets for linear systems requires set-based
addition, which is often referred to as Minkowski addition (X ⊕ Y := {x + y|x ∈
X , y ∈ Y}), and set-based multiplication (X ⊗ Y := {x y|x ∈ X , y ∈ Y}). Note that
the symbol for set-based multiplication is often omitted for simplicity of notation,
and that one or both operands can be singletons. The multiplication with a matrix
M ∈ Ro×n and the Minkowski addition of two zonotopes Z1 = (c, g (1) , . . ., g (p1 ) )
and Z2 = (d, h(1) , . . ., h(p2 ) ) are a direct consequence of the zonotope definition
(see [42]):
(3)

Z1 ⊕ Z2 = (c + d, g (1) , . . . , g (p1 ) , h(1) , . . . , h(p2 ) ),
M ⊗ Z1 = (M c, M g (1) , . . . , M g (p1 ) ).

6

MATTHIAS ALTHOFF

Also, the convex hull of Z1 and Z2 (both having equal order) is required (see [31]):
1
conv(Z1 , Z2 ) ⊆ (c + d, g (1) + h(1) , . . . , g (p1 ) + h(p1 ) ,
2
(4)
c − d, g (1) − h(1) , . . . , g (p1 ) − h(p1 ) ).
For the multiplication of an interval matrix M with a zonotope, the matrix M
is split into a real-valued matrix M ∈ Rn×n and an interval matrix with bound
S ∈ Rn×n , such that M = M ⊕ [−S, S]. After introducing Sj as the j th row of S,
the result is over-approximated as shown in [1, Theorem 3.3] by
MZ1 ⊆(M Z1 ⊕ [−S, S]Z1 )

(i)

hj

⊆(M c1 , M g (1) , . . . , M g (p1 ) , h(1) , . . . , h(n) ),
(
Pp1
|g|(k) ), for i = j
Sj (|c| + k=1
=
.
0, for i 6= j

Another important operation is the enclosure of a zonotope by an axis-aligned box
(see [1, Proposition 2.2]):
box(Z) :=(c, ĥ(1) , . . . , ĥ(n) ),
(P
p1
( k=1
|g|(k) )j , for i = j
(i)
ĥj =
0, for i 6= j

(5)

.

We also require some new operations on zonotopes in the Krylov subspace,
which are introduced later in Sec. 3.
2.3. Reachability Analysis. Reachable set computations are typically performed iteratively for short time intervals
τk := [tk , tk+1 ].
In this work, constant-size time intervals tk := k δ are used to focus on the main
innovations, where k ∈ N is the time step and δ ∈ R+ is referred to as the time
increment. An extension to variable time increments is described in [29].
We recapitulate the reachability analysis of a linear differential inclusion
ẋ ∈ Ax(t) ⊕ U,

(6)
n

n×n

where x ∈ R , A ∈ R
, and U ⊂ Rn is a set of uncertain inputs. Please note
that this also includes the form ẋ = Ax(t) + Bu(t), u(t) ∈ Ũ often used in control
theory, since one could easily choose U = B Ũ. For linear systems, the main task
is to compute the reachable set of the first time interval [0, δ]. Most of this paper
will deal with computing the reachable set for the initial time interval since the
propagation for later time intervals is rather simple, as shown later in Alg. 2. We
take advantage of the superposition principle of linear systems by computing the
following reachable sets separately and later joining them together: the reachable
set of the homogeneous solution Rh (t), the reachable set of the particular solution
Rp (t) due to the uncertain input U, and the reachable set Rǫ correcting the initial
assumption that trajectories are straight lines within [0, δ]. According to [1], the
reachable set for [0, δ] is computed as shown in Fig. 2:
(1) Starting from X0 , compute the set of all homogeneous solutions Rh (δ).
(2) Obtain the convex hull of X0 and Rh (δ) to approximate the reachable set
for the time interval [0, δ].

REACHABILITY ANALYSIS WITH UNCERTAIN INPUTS IN THE KRYLOV SUBSPACE

7

S
(3) Compute R([0, δ]) := t∈[0,δ] R(t) by considering uncertain inputs by
adding Rp (δ) and accounting for the curvature of trajectories by adding
Rǫ .
R([0, δ])
Rh (δ)
PSfrag replacements
convex hull of
X0 , Rh (δ)
X0
enlargement
➀

➁

➂

Figure 2. Steps for computing the reachable set of a linear system
for the first time interval.

3. Homogeneous Solution in the Krylov Subspace
In this section, the basic idea of computing reachable sets as presented in
Sec. 2.3 is extended so that reachable sets can be computed in the Krylov subspace.
We first present new techniques for the homogeneous solution of points in time and
time intervals. Subsequently, we consider for the first time how reachable sets can
be computed for arbitrarily changing input trajectories within the Krylov subspace.
3.1. Solution for a Point in Time. The well-known homogeneous solution
of a linear time-invariant system with initial state x0 is
xh (t) = eAt x0 .
We can bound the exact solution using the lemma below. For that lemma and
n
subsequent derivations we introduce [−1, 1] as an n-dimensional vector whose
entries are intervals [−1, 1]. Analogously, we write [−1, 1]n×m to represent an
n × m matrix whose entries are intervals [−1, 1].
Lemma 3.1 (Single state homogeneous solution). After obtaining V , H from
Alg. 1 with inputs C = A, v = x0 , we can bound the homogeneous solution xh (t) =
eAt x0 by
n

(7)

xh (t) ∈ kx0 kV eH t e1 ⊕ [−1, 1] kx0 kǫnorm t
= x̂h (t) ⊕ Ered (t, x0 )

where
(8)

x̂h (t) := kx0 kV eH t e1
n

Ered (t, x0 ) = [−1, 1] kx0 kǫnorm t.

and ǫnorm is computed as described in Appendix B.

8

MATTHIAS ALTHOFF

Proof. The lemma directly follows from enlarging the approximate solution
in (1) by the error kx0 kǫnorm from (2). Since we have for a vector a ∈ Rn and a
scalar b ∈ R that kak < b ⇒ a ∈ [−1, 1]n b, one obtains from (2) that
n

eA t x0 − kx0 kV eH t e1 ∈ [−1, 1] kx0 kǫnorm t,
which proves the lemma. Please note that the error can simply be multiplied by t
when the exponential matrix has the form eA t (see [62, Corollary 5.3]).

For reachability analysis, one has to compute the homogeneous solution for a
set of initial states. Replacing the single initial state in Lemma 3.1 by a set of
initial states X0 is not trivial since the matrices V and H depend on each initial
state x0 ∈ X0 . To indicate this dependency, we write V (x0 ) and H(x0 ) from now
on to stress that those matrices have been obtained from the state x0 . Since X0
is represented as a zonotope in this work, the homogeneous reachable set can be
computed by the following theorem:
Theorem 3.1 (Homogeneous solution for a point in time). The reachable
set of the homogeneous solution Rh (t) := eA t X0 for the initial zonotope X0 =
(c, g (1) , . . . , g (p) ) can be over-approximated by the zonotope
Rh (t) ⊆ (ĉ, ĝ (1) , . . . , ĝ (p) ) ⊕ Rh,err ,
ĉ = kckV (c)eH(c) t e1 ,
(i)

ĝ (i) = kg (i) kV (g (i) )eH(g ) t e1 ,
p


X
n
kg (i) k ǫnorm t.
Rh,err = [−1, 1] kck +
i=1

Proof. Inserting the definition of a zonotope (Def. 2.1) into Rh (t) = eA t X0
and using (3), we obtain
Rh (t) = e

At

c⊕

p
M

[−1, 1]eA t g (i) .

i=1

Using Lemma 3.1 yields
n

Rh (t) ⊆kckV (c)eH(c) t e1 ⊕ [−1, 1] kckǫnorm t
p 
M
(i)
[−1, 1]kg (i)kV (g (i) )eH(g ) t e1
⊕
i=1

n

⊕ [−1, 1] kg (i) kǫnorm t

=kckV (c)eH(c) t e1
⊕



p
M
(i)
[−1, 1]kg (i)kV (g (i) )eH(g ) t e1
i=1

p


X
kg (i) k ǫnorm t .
⊕ [−1, 1] kck +
n

i=1

|

{z

=Rh,err

}

REACHABILITY ANALYSIS WITH UNCERTAIN INPUTS IN THE KRYLOV SUBSPACE

9

This results in the zonotope of the theorem. Since Rh,err is a multidimensional
interval, and thus also a zonotope, the addition of the two zonotopes Rh,err and
(ĉ, ĝ (1) , . . . , ĝ (p) ) results in the zonotope Rh (t).

We propose two different representations of Rh,err . The first one uses the
equivalence
p
n


M
X
[−1, 1]ei ǫ̂norm t,
kg (i) k ǫnorm t =
[−1, 1]n kck +

(9)

i=1

i=1

|

{z

=:ǫ̂norm t

}

where ei is a unit vector with the ith element being 1 and all others 0. Let us also
introduce 0n and 1n as an n-dimensional vector of zeros and ones, respectively.
Using (9), we can write Rh,err = (0n , ĥ(1) , . . . , ĥ(n) ) with
(
ǫ̂norm t, for i = j
(i)
ĥj =
,
0, for i 6= j
which makes it possible to obtain (see Theorem 3.1)
Rh (t) ⊆(ĉ, ĝ (1) , . . . , ĝ (p) ) ⊕ Rh,err
=(ĉ, ĝ (1) , . . . , ĝ (p) , ĥ(1) , . . . , ĥ(n) ).
The advantage of this technique is that we can iteratively apply Lemma 3.1 when
propagating the homogeneous solution for the next time steps (as later described
in Sec. 5) since the result is in zonotope representation. However, this approach
adds new generators, which slows down the computation for the next time step.
For this reason, we introduce the possibility to add the error Rh,err as an
interval vector. This makes it possible to apply Theorem 3.1 without adding new
generators. Therefore, we have to enclose eA t [b, b] by an n-dimensional interval
within the Krylov subspace to further propagate the error as detailed in Sec. 5.
Lemma 3.2 (Interval vector multiplication). Using [b, b] = bc ⊕ [−b∆, b∆ ] we
can provide the bound
eA t [b, b] ⊆kbc kV (bc )eH(bc ) t e1 ⊕ [−µ, µ],
where
µ = kb∆ kV̄ (|b∆ |)eH̄(|b∆ |) t e1 + 1n (kb∆ kǭnorm t + kbc kǫnorm t)
and V̄ , H̄ are obtained from Alg. 1 and ǭnorm t from Appendix B using C = |A|.
Please note that the absolute values are computed elementwise, i.e., |A|i,j = |Ai,j |.
Proof. We first show that eA t [b, b] can be over-approximated by
eA t [b, b] =eA t (bc ⊕ [−b∆ , b∆ ])
⊆eA t bc ⊕ eA t [−b∆ , b∆ ]
⊆eA t bc ⊕ [−|eA t ||b∆ |, |eA t ||b∆ |]
⊆eA t bc ⊕ [−e|A|t |b∆ |, e|A|t |b∆ |],

10

MATTHIAS ALTHOFF

where the over-approximation achieved in the last line directly follows from the
Taylor series of eA t . Using the above result and Lemma 3.1, we obtain
eA t bc ⊕ [−e|A|t |b∆ |, e|A|t |b∆ |]

(10)

n

=kbc kV (bc )eH(bc ) t e1 ⊕ [−1, 1] kbc kǫnorm t ⊕ [−γ, γ],

where
γ = kb∆ kV̄ (|b∆ |)eH̄(|b∆ |) t e1 + 1n kb∆ kǭnorm t.
The result in (10) can be simplified to
kbc kV (bc )eH(bc ) t e1 ⊕ [−µ, µ],
with µ as in the lemma.



3.2. Solution for a Time Interval. In the previous subsection, we overapproximated the homogeneous solution for points in time. This subsection overapproximates xh (t) for a time interval [0, δ]. Since xh (t) ∈ x̂h (t) ⊕ Ered (t, x0 ) as
shown in (7), where Ered (t, x0 ) is monotonically increasing according to (8), we have
that ∀t ∈ [0, δ] : Ered (t, x0 ) ⊆ Ered (δ, x0 ). It remains to approximate x̂h (t) within
time intervals:
t
∀t ∈ [0, δ] : x̂h (t) ≈ x0 + (x̂h (δ) − x0 ).
δ
To consider the error of this approximation, we use a finite Taylor expansion of the
exponential matrix of η th order with error matrix E (see [1, eq. (3.3)]):
eH(x0 ) t =

(11)

η
X
1
(H(x0 ) t)i + E(t, x0 ),
i!
i=0

where E is enclosed by an interval matrix:

n×n

E(t, x0 ) ∈ E(t, x0 ) = [−1, 1]

(12)

η+1

Φ(t, x0 ) =

(kH(x0 )k∞ t)
(η + 1)!

1
,
1−ǫ

ǫ=

Φ(t, x0 ),
kH(x0 )k∞ t !
< 1.
η+2

!

Please note that < means that the above inequality has to be enforced by choosing
η large enough. We are proposing an enclosure of the error xh (t) − x̂h (t), which is
multiplicative with the initial state, since large initial states result in larger errors
for a given time horizon.
Lemma 3.3 (Correction matrix F ). The interval matrix
F (x0 ) =

η
X

−i

−1

[(i i−1 − i i−1 )δ i , 0]

i=2

H i (x0 )
⊕ E(δ, x0 ),
i!

with E(δ, x0 ) and η according to (12) ensures the enclosure of the exact solution:
(13)

∀t ∈ [0, δ] :

t
x̂h (t) ∈ x0 + (x̂h (δ) − x0 )
δ
⊕ kx0 kV (x0 )F (x0 )e1 .

Proof. Let us start by rearranging (13):
t
x̂h (t) − x0 − (x̂h (δ) − x0 ) ∈ kx0 kV (x0 )F (x0 )e1 .
δ

REACHABILITY ANALYSIS WITH UNCERTAIN INPUTS IN THE KRYLOV SUBSPACE 11

After replacing x̂h (t) by kx0 kV (x0 )eH(x0 ) t e1 from (8) and using kx0 kV (x0 )e1 = x0
(this follows from the fact that the first column of V (x0 ) is x0 /kx0 k according to
Alg. 1), the above inclusion becomes
kx0 kV (x0 )eH(x0 ) t e1 − kx0 kV (x0 )e1

t
kx0 kV (x0 )eH(x0 ) δ e1 − kx0 kV (x0 )e1
−
δ
∈kx0 kV (x0 )F (x0 )e1 .
Next, we use from linear algebra that A(B + C)D = ABD + ACD so that we can
simplify the above equation to


t  H(x0 ) δ
H(x0 ) t
kx0 kV (x0 ) e
−I −
e
− I e1
δ
∈kx0 kV (x0 )F (x0 )e1 .
By comparing the inner part  of kx0 kV (x0 )e1 , one obtains
t
∀t ∈ [0, δ] : eH(x0 ) t − I − (eH(x0 ) δ − I) ∈ F (x0 ).
δ
Substituting eH(x0 ) t by (11) and canceling linear terms yields
η
X

1
t
(ti − t δ i−1 ) H i (x0 ) + E(t, x0 ) − E(δ, x0 ) ∈ F (x0 ).
i!
δ
i=2

The interval of ti −t δ i−1 for t ∈ [0, δ] is obtained exactly by computing the minimum
d i
and maximum of ti −t δ i−1 for which only one extreme value exists: dt
(t −t δ i−1 ) =
1
0 ⇒ tmin = i− i−1 δ. This means that the maximum values are to be found at the
borders of t ∈ [0, δ], which are both 0 for t = 0 and t = δ. Thus,
−i

−1

[(i i−1 − i i−1 )δ i , 0] = {ti − t δ i−1 |t ∈ [0, δ]}.
It remains to bound E(t, x0 ) − δt E(δ, x0 ) for t ∈ [0, δ]. The error can be bounded
η+1

n×n

0 )k∞ t)
1
by E(t, x0 ) ∈ [−1, 1]
φ(t, x0 ) using φ(t, x0 ) = (kH(x(η+1)!
1−ǫ from (12). As
φ(t, x0 ) is strictly increasing with time, it follows that φ(t, x0 ) ∈ [0, 1]φ(δ, x0 ) for
t ∈ [0, δ]. From this follows that φ(t, x0 ) − δt φ(δ, x0 ) ∈ [0, 1]φ(δ, x0 ) − [0, 1]φ(δ, x0) =
[−1, 1]φ(δ, x0 ). Thus, ∀t ∈ [0, δ] : E(t, x0 ) − δt E(δ, x0 ) ∈ E(δ, x0 ), which completes
the proof.


When the initial state x0 is substituted by the set of initial states X0 , Lemma 3.3
can be generalized as described in the subsequent theorem.
Theorem 3.2 (Solution for time intervals). The reachable set for t ∈ [0, δ] is
over-approximated by
[

Rh (t) =conv X0 , Rh (δ) ⊕ N ,
Rh ([0, δ]) :=
t∈[0,δ]

where
N = kckV (c)F (c)e1 ⊕

p
M
i=1

kg (i) kV (g (i) )F (g (i) )e1 .

12

MATTHIAS ALTHOFF

Proof. Let us start by formulating the reachable set for t ∈ [0, δ] using xh (t) ∈
x̂h (t)⊕ Ered (t, x0 ) from (7) and x̂h (t) from (13); we also use ∀t ∈ [0, δ] : Ered (t, x0 ) ⊆
Ered (δ, x0 ):
o
n
t
Rh ([0, δ]) ⊆ x0 + (x̂(δ) − x0 ) x0 ∈ X0
[δ
⊕
kx0 kV (x0 )F (x0 )e1
(14)
x0 ∈X0
[
Ered (δ, x0 ).
⊕
x0 ∈X0

Since
o
n
t
∀t ∈ [0, δ] : x0 + (x̂(δ) − x0 ) x0 ∈ X0
δ
[
⊕
Ered (δ, x0 )
x0 ∈X0

⊆ conv(X0 , Rh (δ))

we can simplify (14) to
Rh ([0, δ]) ⊆conv(X0 , Rh (δ))
[
⊕
kx0 kV (x0 )F (x0 )e1
x0 ∈X0

⊆conv(X0 , Rh (δ)) ⊕ N̂ ,

S

where N̂ := x0 ∈X0 kx0 kV (x0 )F (x0 )e1 . It remains to over-approximate N̂ when
the initial set is a zonotope:
[
kx0 kV (x0 )F (x0 )e1
x0 ∈X0

⊆ kckV (c)F (c)e1 ⊕
|
{z
}
=:N (0)

p
M
i=1

[−1, 1]kg (i)kV (g (i) )F (g (i) )e1 =: N ,
|
{z
}
=kg(i) kV (g(i) )F (g(i) )e1 =:N (i)

and the interval vectors N (i) are added using standard interval arithmetic [40].
Thus, N is an interval vector, which is converted to a zonotope and added to the
convex hull. The interval [−1, 1] in the above computation can be removed since
all entries of F (g (i) ) have the form [−a, a], a ∈ R+ .

Next, we derive the set of solutions due to uncertain inputs.
4. Input Solution in the Krylov Subspace
In this section, we obtain for the first time the set of input solutions for
uncertain, time-varying inputs in the Krylov subspace. We first consider overapproximations for solutions of a single constant input. Next, we generalize this
result to uncertain but constant inputs. Finally, we derive an over-approximation
for arbitrarily-varying and uncertain inputs. The first lemma provides the overapproximation for a constant, known input.

REACHABILITY ANALYSIS WITH UNCERTAIN INPUTS IN THE KRYLOV SUBSPACE 13

Lemma 4.1 (Krylov error of constant input solution). After obtaining Ṽ (u),
H̃(u) from Alg. 1 and ǫ̃norm δ from Appendix B with inputs


 n
A
u
0
C = 1×n
, v=
,
0
0
1
we can bound the particular solution (aka input solution) for constant inputs
Z δ
Z δ
A(δ−t)
(15)
xp,const (δ) =
e
dt u =
eAt dt u
0

0

by

n

xp,const (δ) ∈ P Ṽ (u)eH̃(u) δ e1 ⊕ [−1, 1] ǫ̃norm δ
= x̃p (δ) ⊕ Ẽred (δ),
where
P = [I, 0n ]
x̃p (δ) = P Ṽ (u)eH̃(u) δ e1

(16)

n

Ẽred (δ) = [−1, 1] ǫ̃norm δ.
Proof. Using [56, Sec. 5] we first rewrite the solution of (15) as


 
 n

ẋp,const
0
A u xp,const
,
x̃(0)
=
x̃
=
.
=
0
x̌
1
0 0
x̌˙
| {z } | {z } | {z }
x̃˙

=:Ã(u)

x̃

After inserting the projection matrix P = [I, 0n ] (see (16)), we can write
xp,const (t) = P eÃ(u)t x̃0 .
This makes it possible to reformulate the Krylov error of the input solution:
k P eÃ(u) t x̃0 −
| {z }
xp,const (t)

≤kP k ke

Ã(u) t

P kx̃0 kṼ (u)eH̃(u) t e1
{z
}
|

k

Krylov approximation of xp,const (t)

x̃0 − kx̃0 kṼ (u)eH̃(u) t e1 k,

(kP k = kx̃0 k = 1)

=keÃ(u) t x̃0 − Ṽ (u)eH̃(u) t e1 k ≤ ǫ̃norm t
and ǫ̃norm t is obtained as in Appendix B.



In order to generalize the previous results to arbitrarily-varying inputs, we
require the following corollary on the input solution for general time bounds:
Corollary 4.1 (Input solution for general time bounds). The partial input
solution for
Z te
eA(δ−t) dt u,
xp,const ([t0 , te ]) :=
t0

where 0 ≤ t0 ≤ te ≤ δ, can be over-approximated as

xp,const ([t0 , te ]) ∈P Ṽ (u)(eH̃(u) (δ−t0 ) − eH̃(u) (δ−te ) )e1
n

⊕ [−1, 1] ǫ̃norm (te − t0 ).

14

MATTHIAS ALTHOFF

Proof. Let us rewrite xp,const ([t0 , te ]) as
xp,const ([t0 , te ]) :=

Z

te

e

A(δ−t)

Z

dt u =

Z

δ−t0

eAt dt −

Z


eAt dt u.

0

0

eAt dt u

δ−te
δ−te

t0

=

δ−t0

After applying Lemma 4.1, we obtain
xp,const ([t0 , te ])
=P Ṽ (u)eH̃(u) (δ−t0 ) e1 − P Ṽ (u)eH̃(u) (δ−te ) e1
n

⊕ [−1, 1] ǫ̃norm (δ − t0 − (δ − te ))
=P Ṽ (u)(eH̃(u) (δ−t0 ) − eH̃(u) (δ−te ) )e1
⊕ [−1, 1]n ǫ̃norm (te − t0 ).

Next, we over-approximate the solution for constant inputs which are uncertain
within the set U. In many cases, it is desired to have the solution for constant inputs,
e.g., if a control system with zero-order hold is considered.
Theorem 4.1 (Reachable set for constant inputs). The reachable set of the
input solution
R̂p (δ) :=

Z

δ

e

A(δ−t)

dt U =

0

(1)

Z

δ

eAt dt U,
0

(q)

where U = (cu , gu , . . . , gu ) is a zonotope, can be over-approximated by the zonotope
R̂p (δ) ⊆ (c̃u , g̃u(1) , . . . , g̃u(q) ) ⊕ R̂p,err ,
H̃(cu ) δ
c̃(i)
e1 ,
u = P Ṽ (cu )e
(i)

g̃u(i) = P Ṽ (gu(i) )eH̃(gu ) δ e1 ,
q


X
n
ǫ̃norm (gu(i) ) δ.
R̂p,err = [−1, 1] ǫ̃norm (cu ) +
i=1

Rδ
0

Proof. After inserting the definition of a zonotope (Def. 2.1) into R̂p (δ) :=
e
dt U and using (3), we obtain
A(δ−t)

R̂p (δ) =

Z

0

δ

e

A(δ−t)

Z
q
M
[−1, 1]
dt cu ⊕
i=1

0

δ

eA(δ−t) dt gu(i) .

REACHABILITY ANALYSIS WITH UNCERTAIN INPUTS IN THE KRYLOV SUBSPACE 15

Using Lemma 4.1, one receives
n

R̂p (δ) ⊆P Ṽ (cu )eH̃(cu ) δ e1 ⊕ [−1, 1] ǫ̃norm (cu ) δ
q 
M
(i)
[−1, 1]P Ṽ (gu(i) )eH̃(gu ) δ e1
⊕
i=1

⊕ [−1, 1]n ǫ̃norm (gu(i) ) δ

=P Ṽ (cu )eH̃(cu ) δ e1
⊕



q
M
(i)
[−1, 1]P Ṽ (gu(i) )eH̃(gu ) δ e1
i=1

q


X
ǫ̃norm (gu(i) ) δ .
⊕ [−1, 1]n ǫ̃norm (cu ) +
i=1

|

{z

=R̂p,err

}

This results in the zonotope of the theorem. Since R̂p,err is a multidimensional
interval, and thus also a zonotope, the addition of the two zonotopes R̂p,err and
(1)
(q)

(c̃u , g̃u , . . . , g̃u ) results in the zonotope R̂p (δ).
The above derivations only hold for constant inputs. We generalize the previous
results to arbitrary input trajectories in the following theorem.
Theorem 4.2 (Reachable set for varying inputs). The reachable set due to
uncertain inputs
Z δ
n
Rp (δ) = xp (δ) xp (δ) =
eA(δ−t) u(t) dt,
0
o
∀t ∈ [0, δ] : u(t) ∈ U

can be over-approximated by

η n
M
o
H̃ j (u)
Rp (δ) ⊆P
Ṽ (u)
u ∈ U δj
j!
j=1
n
o
n×n
⊕[−1, 1]
Φ̃(δ, u) u ∈ U e1
n
o
⊕[−1, 1]n ǫ̃norm (u) u ∈ U δ.

The proof can be found in Appendix A.
Next, we consider a particular solution of Theorem 4.2 when U is a zonotope.
Corollary 4.2 (Varying inputs with zonotopic bounds). The reachable set
due to uncertain inputs
Z δ
n
Rp (δ) = xp (δ) xp (δ) =
eA(δ−t) u(t) dt,
0
o
∀t ∈ [0, δ] : u(t) ∈ U

16

MATTHIAS ALTHOFF
(1)

(q)

when U = (cu , gu , . . . , gu ) is a zonotope, can be over-approximated by
Rp (δ) =(c̃u , g̃u(1,1) , . . . , g̃u(q,1) , g̃u(1,2) , . . . , g̃u(q,2) , . . . , g̃u(q,η) )
⊕ R̂p,err ,

X
η 
H̃ j (cu )  j
δ e1 ,
c̃u = P
Ṽ (cu )
j!
j=1
(i)

H̃ j (gu ) j
δ e1 ,
j!

q


X
n×n
Φ̃(δ, gu(i) ) e1
= P [−1, 1]
Φ̃(δ, cu ) +

g̃u(i,j) = P Ṽ (gu(i) )
Rp,err

i=1

q


X
n
ǫ̃norm (gu(i) ) δ.
⊕ [−1, 1] ǫ̃norm (cu ) +
i=1

Proof. After inserting the definition of a zonotope (Def. 2.1) into the result
of Theorem 4.2, we obtain
M
η 
H̃ j (cu )
Ṽ (cu )
Rp (δ) ⊆P
j!
j=1
⊕

q
M

[−1, 1]Ṽ (gu(i) )

i=1

⊕ [−1, 1]

n×n

(i)
H̃ j (gu )  j
δ
j!

q


M
[−1, 1]Φ̃(δ, gu(i) ) e1
Φ̃(δ, cu ) ⊕
i=1

|



{z

=[−1,1]n×n Φ̃(δ,cu )+

Pq

(i)
i=1 Φ̃(δ,gu )

}



q


M
[−1, 1]ǫ̃norm(gu(i) ) δ ,
⊕ [−1, 1]n ǫ̃norm (cu ) ⊕
i=1

|

n

=[−1,1]



{z

ǫ̃norm (cu )+

Pq

(i)
i=1 ǫ̃norm (gu )



δ

}

which can be simplified to
X
η 


H̃ j (cu )  j
δ e1
Rp (δ) ⊆P
Ṽ (cu )
j!
j=1


q
η M
j (i) 
M
(i) H̃ (gu )
j
[−1, 1] P Ṽ (gu )
⊕
δ e1 ⊕ Rp,err .
j!
j=1 i=1
This results in the zonotope of the theorem. Since Rp,err is a multidimensional
interval, and thus also a zonotope, the addition of the two zonotopes Rp,err and
(1,1)
(q,1) (1,2)
(q,2)
(q,η)
(c̃u , g̃u , . . . , g̃u , g̃u , . . . , g̃u , . . . , g̃u ) results in the zonotope Rp (δ). 

REACHABILITY ANALYSIS WITH UNCERTAIN INPUTS IN THE KRYLOV SUBSPACE 17

5. Propagation
So far, we have only described the computation of a reachable set for the first
time interval. This section presents how the initial results are propagated for further
consecutive time intervals τi .
5.1. Overall Algorithm. To grasp the implemented propagation scheme more
easily, we start with one of the simplest propagation procedures (see [31]):
(17)

Rh (τk+1 ) = eAδ Rh (τk ),

(18)

Rp (τk+1 ) = eAδ Rp (τk ) ⊕ Rp (τ0 ).

In order to keep the number of generators of zonotopes that are multiplied by the
matrix exponential small, the wrapping-free approach from [33] is used in this work
as a basis; we later introduce modifications to make the best use of the proposed
computation in the Krylov subspace. The wrapping-free approach modifies the
above procedure by introducing the auxiliary reachable set Rb , which is enclosed
by an axis-aligned box: Rp (τ0 ) = box(Rb (τ0 )) (see line 3 of Alg. 2 and (5)). This has
the effect that the representation of Rp (τk+1 ) in (18) does not grow in complexity
due to the Minkowski addition, but stays a simple axis-aligned box. The modified
computation of Rp (τk+1 ) according to [33] is:
(19)

Rb (τk+1 ) = eAδ Rb (τk ),

(20)

Rp (τk+1 ) = Rp (τk ) ⊕ box(Rb (τk+1 )).

Another modification, which is not proposed in [33], is that we change (17): instead of computing the homogeneous solution for the first time interval only and
then propagating it, we apply Theorem 3.2 in each time step. The reason is that
for large systems, one saves much computation time when the matrix exponential
multiplication is performed with Rh (tk ) instead of Rh (τk ) due to a fewer number
generators, while the computation in Theorem 3.2 is negligible when using zonotopes. Since the sets in Alg. 2 are indexed without explicitly indicating time, we
distinguish sets representing points in time by an asterisk (see e.g., R∗h,1 ) from the
ones representing time intervals.
To sum up, in the Krylov subspace, the computation of (17) and (19) is carried
out by Theorem 3.1 (see line 6 and line 7 in Alg. 2). The propagation of the
particulate solution (20) is realized in line 8 of Alg. 2 and the homogeneous solution
for a time interval is obtained in line 9. The aggregation of the homogeneous and
the particulate solution are performed in line 4 and line 10 of Alg. 2.
When only constant inputs are considered, corollary 4.2(A, B, U, δ) in line 2
is replaced by theorem 4.1(A, B, U, δ). Also, when the center uc of the set of
uncertain inputs U is large compared to the deviation from the center U∆ := U ⊕
(−uc ), one should move the solution of uc inside the convex hull computation
performed in theorem 3.2(A, X0 , δ) in line 9 (see e.g., [7, Alg. 1]). This, however,
is independent from the extension to compute in the Krylov subspace and thus not
discussed in this work.
5.2. Computational Complexity. Let us first consider the computational
complexity when applying Alg. 2 without utilizing the Krylov subspace: For each
time step, we have to map zonotopes using the matrix exponential, compute the
box enclosure as well as the convex hull of zonotopes, and add two zonotopes. We
introduce the number of generators of R∗h,1 in Alg. 2 as ph and the number of

18

MATTHIAS ALTHOFF

Algorithm 2 Compute R([0, tf ])
Require: State matrix A, input matrix B, initial set X0 , input set U, time step δ,
time horizon tf
Ensure: R([0, tf ])
1: R∗
h,1 = theorem 3.1(X0 , A, δ)
2: R∗
b,1 = corollary 4.2(U, A, B, δ)
3: Rp,1 = box(Rb,1 )
4: R1 = Rh,1 ⊕ Rp,1
5: for k = 1 . . . ⌈tf /δ − 1⌉ do
6:
R∗h,k+1 = theorem 3.1(R∗h,k , A, δ)
7:
R∗b,k+1 = theorem 3.1(R∗b,k , A, δ)
8:
Rp,k+1 = Rp,k ⊕ box(R∗b,k+1 )
9:
Rh,k+1 = theorem 3.2(R∗h,k , A, δ)
10:
Rk+1 = Rh,k+1 ⊕ Rp,k+1
11: end for
Stf /δ
12: R([0, tf ]) = k=1 Rk
generators of R∗b,1 as pb . Thus, the orders of those zonotopes are oh = pnh and ob =
pb
n . The number of required binary operations for each of the previously mentioned
high-level operands are listed in Tab. 1. We are interested in the complexity with
respect to the dimension n, since the number of time steps is typically fixed or
given by reaching a fixed-point. As can be seen from Tab. 1, the complexity with
respect to n is cubic for linear maps (when M is quadratic as for eAδ ), linear for
addition, quadratic for over-approximating the convex hull, and quadratic for the
box enclosure. Thus, the overall complexity is dominated by the linear map, which
has complexity O((oh + ob ) n3 ). Note that according to Alg. 2, the order of the
involved zonotopes does not grow compared to other propagation algorithms (see
e.g., [1, 31]).
The only difference when computing in the Krylov subspace is that the complexity of computing eAδ Z changes. As discussed above, this is also the operation
Table 1. Required operations.
operands
(1)

(p )

(1)

(pp )

(1)

(p )

Zh = (ch , gh , . . . , gh h ), Ẑh = (ĉh , ĝh , . . . , ĝh h ) ,
Zp = (cp , gp , . . . , gp

) ⊂ Rn , M ∈ Rm×n

operation

nr. of binary operations

M Zh , see (3)

2mn(oh n + 1)

Zh ⊕ Zp , see (3)

n

conv(Zh , Ẑh ), see (4) oh n2
box(Zh ), see (5)

oh n2

REACHABILITY ANALYSIS WITH UNCERTAIN INPUTS IN THE KRYLOV SUBSPACE 19

which determines the overall computational complexity. We are discussing the case
when the Krylov order is chosen such that the exponential matrix is computed up
to machine precision so that the results have the same accuracy and the comparison is fair. The Arnoldi iteration requires ξ matrix-vector multiplications (ξ is
the dimension of the Krylov subspace; see Alg. 1). However, since the involved
matrices are sparse, we can add a sparsity constant s << 1 so that we obtain
O(sξn2 ) operations. Further, we require ξ 2 /2 inner products (O(ξ 2 n)) and ξ 2 /2
other operations with O(n) resulting in O(ξ 2 n), so that the overall complexity is
O(ξ 2 n) + O(sξn2 ) = O(sξn2 ) with respect to n.
The Arnoldi iteration has to be computed for (oh + ob )n generators, so that
the computational complexity amounts to O(sξ(oh + ob )n3 ) concerning the Arnoldi
iteration. The required computation of each exponential matrix eHδ in the Krylov
subspace only depends on ξ and not on n. Finally, we require matrix/matrix
computations between V and the exponential matrix eHδ with complexity O(nξ 2 ).
Thus, the overall complexity using the Krylov technique is O(sξ(oh + ob )n3 ) +
O(nξ 2 ) = O(sξ(oh +ob )n3 ). As a result, both the classical approach and the Krylov
approach are cubic in the number of continuous state variables. However, the
Krylov method can be much faster (as demonstrated in the next section) depending
on the sparsity of the system matrix A, since typically sξ << 1 (compare O((oh +
ob ) n3 ) for the classical approach with O(sξ(oh + ob )n3 ) for the Krylov approach).
Please note that the sparsity typically increases the larger a system is (see e.g., [63,
Fig. 3]. Also, it should be noted that although A is often very sparse, all entries of
eAδ are typically non-zero.
6. Numerical Experiments
To demonstrate the usefulness of the presented approach, we consider the formal analysis of a bridge as a representative of a safety-critical structure. To the best
knowledge of the author, no work on formally bounding values of a safety-critical
structure exists, whose safety can determine life or death; the most prominent
example is the collapsed Tacoma Narrows Bridge [11].
6.1. Bridge Model. The model of the Roosevelt Lake Bridge (Arizona) is
taken from a student work [50], which investigated its structural dynamics. A
picture1 of the bridge and the corresponding finite element model are shown in
Fig. 3. The bridge is mostly made out of steel, except for the roadway deck. The
dam structure in front of the bridge (see Fig. 3(a)) creates a special aerodynamic
situation, which can excite certain natural frequencies of the bridge [24]. Using the
proposed reachability analysis, we can investigate the entire range of possible time
responses for all kinds of frequencies of the exciting forces at once. Thus, we do
not miss a single possible frequency that could cause a problem.
The bridge model consists of 445 nodes connecting bars and beams. The dynamics of the bridge can be formulated as
M q̈(t) + Dq̇(t) + Kq(t) = g(t),
where M is the mass matrix, D is the damping matrix, K is the stiffness matrix,
q is the vector of generalized coordinates, and g(t) summarizes the external loads.
1The picture is taken from commons.wikimedia.org

20

MATTHIAS ALTHOFF

The above model can be directly transformed to the state space format used in this
work by introducing x1...m = q, xm+1...2 m = q̇ so that
ẋ(t) = A x(t) + B u(t),
where




 m×m 
0m×m
I
0
A=
.
,
B=
M −1
−M −1 K −M −1 D
Since the bridge has 2520 degrees of freedom, the dimension of the linear system is
n = 5040. The exact model is provided as an example in CORA [2].

PSfrag replacements
reference
trajectory
vehicle I
(a) Picture of the bridge with the dam accelerating winds in the background.

60
PSfrag replacements
reference 40
trajectory
vehicle I20
ego vehicle
previously verified part 0
ego vehicle
new intended part 0

-150

-100

-50

0

50

100

150

(b) Finite-element model of the bridge (distances are in meters).

Figure 3. Picture and finite-element model of the Roosevelt Lake
Bridge (Arizona).
6.2. Excitation of the Bridge. We consider two excitations of the bridge:
(1) The bridge is excited by lateral forces simulating the winds acting on the bridge
and the street deck and (2) vertical excitation is generated by the street deck caused
by moving traffic on the bridge. It should be noted that we consider all kinds of
frequencies of these excitations so that no possible solution is missed. Due to the
lengths of the input vector u (5040 entries), we provide u within the uploaded
example in CORA.
6.3. Results. The results below are obtained using a standard laptop with
an Intel i7-3520M CPU running at 2.90GHz and 8GB of RAM.

REACHABILITY ANALYSIS WITH UNCERTAIN INPUTS IN THE KRYLOV SUBSPACE 21

Numerical simulation. Since we consider the full model and do not reduce it by
removing higher frequencies, the numerical simulation of the bridge model is already
computationally demanding. For a time interval [0, 10−4 ] s, the average time step
of the ode45 solver of MATLAB is 5.06 · 10−8 s and the average computation time
for one simulation is 122 s. Please note that eigenfrequencies of steel constructions
can be very high; in our case, the highest frequency is 1.64 · 106 .
Reachability analysis. For the reachability analysis we have selected a time
step δ = 1e − 7 that is a little larger than the average step size for the numerical
simulation. To be as precise as possible, we have chosen the order reduction so that ǫ
is below the precision of floating point numbers in MATLAB, which is 2.2204·10−16.
Given the proposed accuracy and time step, the reduced order ξ for the system
matrix A is 120 and for Ã we obtain 124, which is a significant reduction from the
original dimension 5040.
The presented technique requires 101 s for the initial time interval and 54 s for
each subsequent one. The original approach (see e.g., [8, 33]) requires 1912 s for
the initial time interval and 489 s per time step. It should also be mentioned that
our approach can be easily parallelized since the Arnoldi algorithm (Alg. 1), which
consumes almost all of the computation time, can be executed independently for
each generator, as shown in Theorem 3.1.
In Fig. 4, selected projections of the 5040-dimensional state space are shown,
which is entirely computed, unlike in methods that only consider the output of a
system.
Discussion. As demonstrated, the proposed approach is roughly 10 times faster
while providing the same accuracy up to the floating point precision of MATLAB.
Further, MATLAB is not a suitable language for performing the Arnoldi algorithm
since it cannot be vectorized (i.e., removing for-loops with vector-based operations).
Thus, we expect more significant advantages with languages that are compiled, such
as C++. It should also be noted that the reachable set computation takes as long as
443 simulations, while the problem with just 10 uncertain initial states has already
1024 vertices.

7. Conclusions
We have presented the first work for computing the over-approximative reachable set of linear systems in the Krylov subspace for arbitrarily-varying, bounded
inputs. Unlike most other work on applying reduction techniques for reachability
analysis, which compute an output abstraction, we perform a state abstraction and
can fully reconstruct the complete reachable sets. When using output abstractions,
only the reachable outputs can be over-approximated—if those output-abstraction
techniques would be used to reconstruct the whole reachable set, no reduction would
be achieved. Due to the strict consideration of error bounds, our approach can be
used for formal verification and other formal techniques, such as computation of
invariance sets, computation of the region of attraction, optimization of constrained
systems with uncertainties, set-based observers, and conformance checking. In our
numerical example we have seen a 10 times speed up although the approximation
was accurate up to floating point precision.
The only other work computing state abstraction, to the best knowledge of
the author, is [38]. Besides being able to handle arbitrarily varying inputs, the

22

MATTHIAS ALTHOFF

Figure 4. Reachable set for selected projections. Black lines show
random simulations, the gray area shows the reachable set.

approach presented here extends the previous work by not requiring affine set representations, which uniformly bloat the result equally in all directions and often
cause unnecessarily large over-approximations. Also, we do not use accumulating
norm bounds. The usefulness of our approach is demonstrated for a realistic model
of an existing bridge.

REACHABILITY ANALYSIS WITH UNCERTAIN INPUTS IN THE KRYLOV SUBSPACE 23

Acknowledgment
This work was partly supported by the German Research Foundation (DFG)
under grant number AL 1185/5-1.

Appendix A. Proof of Theorem 4.2
Proof. We first compute an over-approximative reachable set resulting from
inputs when assuming constant inputs for time intervals [tk−1 , tk [, where 0 = t0 <
t1 < . . . < tl−1 < tl = δ, of equal duration γ = tk − tk−1 , so that

u(0)




u(t1 )
u(t) = .
..





u(tl−1 )

for t ∈ [0, t1 [
for t ∈ [t1 , t2 [
for t ∈ [tl−1 , δ].

In a second step, we let γ → 0 to obtain an over-approximation for arbitrarily
varying inputs. The solution of piecewise constant inputs is obtained from

xp,pw (δ) =

Z

t1

eA(δ−t) dt u(0) +

Z

t2

eA(δ−t) dt u(t1 )

t1

0

+ ...+

Z

δ

eA(δ−t) dt u(tl−1 ).

tl−1

When the set of inputs is uncertain within U, we obtain

(21)

R̂p,l (δ) =

Z

t1

eA(δ−t) dt U ⊕ . . . ⊕

0

Z

δ

eA(δ−t) dt U.

tl−1

In order to obtain not only an approximation, but an over-approximation, the
solution for a time interval [tk−1 , tk [ is further abstracted. From Corrollary 4.1 we
have that
Z

te

eA(δ−t) dt U
R̂p ([t0 , te ]) :=
t0
n
o
H̃(u) (δ−t0 )
⊆P Ṽ (u)(e
− eH̃(u) (δ−te ) ) u ∈ U e1
n
o
n
⊕ [−1, 1] ǫ̃norm (u) u ∈ U (te − t0 ).

24

MATTHIAS ALTHOFF

The above over-approximation is rewritten using a finite Taylor series (see (11)):
n

H̃ 1 (u)
R̂p ([t0 , te ]) ⊆ P Ṽ (u) I +
(δ − t0 )1 + . . .
1!

H̃ 1 (u)
H̃ η (u)
(δ − t0 )η + Ẽ(δ − t0 , u) − I +
(δ − te )1
+
η!
1!

H̃ η (u)
(δ − te )η + Ẽ(δ − te , u)
u ∈ U,
+ ...+
η!
o
n
o
n
Ẽ(t, u) ∈ Ẽ(t, u) e1 ⊕ [−1, 1] ǫ̃norm (u) u ∈ U (te − t0 )

n
H̃ 1 (u)
[(δ − t0 )1 − (δ − te )1 ] + . . .
=P Ṽ (u)
1!
H̃ η (u)
+ Ṽ (u)
[(δ − t0 )η − (δ − te )η ]
η!

+ Ẽ(δ − t0 , u) − Ẽ(δ − te , u) u ∈ U,
n
o
o
Ẽ(t, u) ∈ Ẽ(t, u) e1 ⊕ [−1, 1]n ǫ̃norm (u) u ∈ U (te − t0 ).

Finally, the uncertainty of the input is moved inwards, which results in a further
over-approximation:

R̂p ([t0 , te ]) ⊆ P

(22)

o
n
H̃ 1 (u)
u ∈ U [(δ − t0 )1 − (δ − te )1 ] ⊕ . . .
Ṽ (u)
1!
n
o
H̃ η (u)
⊕ Ṽ (u)
u ∈ U [(δ − t0 )η − (δ − te )η ]
η!
n
o
⊕ Ẽ(δ − t0 , u) − Ẽ(δ − te , u) u ∈ U, Ẽ(t, u) ∈ Ẽ(t, u)

n
o
e1 ⊕ [−1, 1]n ǫ̃norm (u) u ∈ U (te − t0 ).

After introducing

one can rewrite (22) as

o
n
H̃ j (u)
u∈U ,
D(j) := Ṽ (u)
j!
Z

(23)

te

eA(δ−t) dt U
R̂p ([t0 , te ]) =
t0

⊆P D(1) [(δ − t0 )1 − (δ − te )1 ] ⊕ . . .

⊕ D(η) [(δ − t0 )η − (δ − te )η ]
n
o
⊕ Ẽ(δ − t0 , u) − Ẽ(δ − te , u) u ∈ U, Ẽ(t, u) ∈ Ẽ(t, u)

n
o
n
e1 ⊕ [−1, 1] ǫ̃norm (u) u ∈ U (te − t0 ).

REACHABILITY ANALYSIS WITH UNCERTAIN INPUTS IN THE KRYLOV SUBSPACE 25

By repeatedly inserting (23) into (21), rearranging time intervals (δ − tk ) = tl−k ,
n×n
and using Ẽ(t, u) ∈ [−1, 1]
Φ̃(t, u) from (12), we obtain

R̂p,l (δ) =P
⊕
⊕
⊕
⊕
(24)

n

(D(1) (tl − tl−1 )
⊕ . . . ⊕D(η) (tηl − tηl−1 )
(D(1) (tl−1 − tl−2 ) ⊕ . . . ⊕D(η) (tηl−1 − tηl−2 )
...
⊕ . . . ⊕ D(η) (tη1 − tη0 )
(D(1) (t1 − t0 )
{z
}
|
|
{z
}
L
L
n×n

[−1, 1]

⊕
⊕
⊕

u∈U
⊕
⊕
⊕
⊕


Φ̃(tl , u) − Φ̃(tl−1 , u)


n×n
[−1, 1]
(Φ̃ tl−1 , u) − Φ̃(tl−2 , u)
...

n×n
[−1, 1]
Φ̃(t1 , u) − Φ̃(t0 , u)
|
{z
}
P

o
e1

n
o
[−1, 1]n ǫ̃norm (u) u ∈ U (tl − tl−1 )
n
o
n
[−1, 1] ǫ̃norm (u) u ∈ U (tl−1 − tl−2 )
...
n
o
n
[−1, 1] ǫ̃norm (u) u ∈ U (t1 − t0 ) .
}
| {z
P

The summation symbols indicate that the terms written in one column are summed
up. Since the expressions tjk − tjk−1 are positive scalars, the following statement can
be used for the summation: For any two positive scalars a, b ∈ R+ and the convex
set S, one can state that {a s + b s|s ∈ S} = aS + bS = (a + b)S. From this follows
that
(25)

D(j) (tj1 − tj0 ) ⊕ D(j) (tj2 − tj1 ) ⊕ . . . ⊕ D(j) (tjl − tjl−1 )
= D(j) (tjl − tj0 ) = D(j) tjl = D(j) δ j

and similarly for the Taylor remainder terms


[−1, 1] Φ̃(t1 , u) − Φ̃(t0 , u) ⊕ [−1, 1] Φ̃(t2 , u) − Φ̃(t1 , u)

⊕ . . . ⊕ [−1, 1] Φ̃(tl , u) − Φ̃(tl−1 , u)

= [−1, 1] Φ̃(tl , u) − Φ̃(t0 , u) = [−1, 1]Φ̃(δ, u)

and the Krylov error terms

n
o
ǫ̃norm (u) u ∈ U (t1 − t0 ) ⊕ . . .
n
o
n
[−1, 1] ǫ̃norm (u) u ∈ U (tl − tl−1 )
n
o
n
= [−1, 1] ǫ̃norm (u) u ∈ U δ.
n

[−1, 1]

(26)

26

MATTHIAS ALTHOFF

Inserting (25) - (26) into (24) yields

R̂p,l (δ) =P D(1) δ ⊕ . . . ⊕ D(η) δ η
n
o
n×n
⊕[−1, 1]
Φ̃(δ, u) u ∈ U e1
n
o
⊕[−1, 1]n ǫ̃norm (u) u ∈ U δ.

Since the above result is independent of the number of intermediate time steps l,
we can choose l → ∞, meaning that
lim R̂p,l (δ) = Rp (δ).

l→∞


Appendix B. Error Bound for the Approximation of Matrix
Exponentials
We summarize the complete procedure to obtain the error bound from [62,
Corollary 5.3] for stable (i.e., positive definite) matrices C, which satisfies
keC v − kvkV eH e1 k ≤ kvkǫnorm.
To obtain ǫnorm , we require from [62, Eq. (5.1)]
1
a = min{λi ( (C + C ∗ ))},
i
2
1
b = max{λi ( (C + C ∗ ))},
i
2
1
c = max{|λi ( (C − C ∗ ))|},
i
2
where λi () returns the ith eigenvalue and C ∗ is the conjugate transpose of C. These
values are used to determine the next auxiliary variable m, which is obtained by
solving [62, Eq. (4.3)]:
(27)

ν :=

E(1 − m) − mK(1 − m)
E(m) − (1 − m)K(m)
=
,
c
0.5(b − a)

using the Jacobi elliptic integrals (see [62, Eq. (2.6)])
Z π2
1
K(m) =
(1 − m sin2 θ)− 2 dθ,
0

E(m) =

Z

π
2

1

(1 − m sin2 θ) 2 dθ.

0

As a byproduct of solving (27), we obtain ν, which is another auxiliary variable.
Next, an optimization is performed according to Alg. 3 to find the best compromise
for two opposing error sources, which solves [62, Eq. (5.9)].
Finally, ǫnorm is obtained according to [62, Corollary 5.3] as
ǫnorm = 22.16kCk

q ξ−1 −(a−1/(2ν)(1/q−q))
e
.
1−q

REACHABILITY ANALYSIS WITH UNCERTAIN INPUTS IN THE KRYLOV SUBSPACE 27

Algorithm 3 Compute optimal value for q
Require: reduced dimension ξ, ν, m
Ensure: q
1: q = 0,
q = 0, q = 1, res → ∞
2: C̃ = 1/(2ν)
3: while |res| > eps do
p
4:
res = (ξ − 1)q + (2 − ξ)q 2 − C̃(1 − q) (1 − q 2 )2 + 4mq 2
5:
if res < 0 then
6:
q = q, q = 0.5(q + q)
7:
else
q = q, q = 0.5(q + q)
8:
9:
end if
10: end while

References
[1] M. Althoff. Reachability Analysis and its Application to the Safety Assessment of
Autonomous Cars. Dissertation, Technische Universität München, 2010. http://nbnresolving.de/urn/resolver.pl?urn:nbn:de:bvb:91-diss-20100715-963752-1-4.
[2] M. Althoff. An introduction to CORA 2015. In Proc. of the Workshop on Applied Verification
for Continuous and Hybrid Systems, pages 120–151, 2015.
[3] M. Althoff, S. Bak, D. Cattaruzza, X. Chen, G. Frehse, R. Ray, and S. Schupp. ARCHCOMP17 category report: Continuous and hybrid systems with linear continuous dynamics.
In Proc. of the 4th International Workshop on Applied Verification for Continuous and
Hybrid Systems, pages 143–159, 2017.
[4] M. Althoff and J. M. Dolan. Online verification of automated road vehicles using reachability
analysis. IEEE Transactions on Robotics, 30(4):903–918, 2014.
[5] M. Althoff and G. Frehse. Combining zonotopes and support functions for efficient reachability analysis of linear systems. In Proc. of the 55th IEEE Conference on Decision and
Control, pages 7439–7446, 2016.
[6] M. Althoff and B. H. Krogh. Zonotope bundles for the efficient computation of reachable sets.
In Proc. of the 50th IEEE Conference on Decision and Control, pages 6814–6821, 2011.
[7] M. Althoff and B. H. Krogh. Reachability analysis of nonlinear differential-algebraic systems.
IEEE Transactions on Automatic Control, 59(2):371–383, 2014.
[8] M. Althoff, C. Le Guernic, and B. H. Krogh. Reachable set computation for uncertain timevarying linear systems. In Hybrid Systems: Computation and Control, pages 93–102, 2011.
[9] M. Althoff, A. Rajhans, B. H. Krogh, S. Yaldiz, X. Li, and L. Pileggi. Formal verification
of phase-locked loops using reachability analysis and continuization. Communications of the
ACM, 56(10):97–104, 2013.
[10] M. Althoff, O. Stursberg, and M. Buss. Reachability analysis of nonlinear systems with
uncertain parameters using conservative linearization. In Proc. of the 47th IEEE Conference
on Decision and Control, pages 4042–4048, 2008.
[11] O. H. Amman, T. von Kármán, and G. B. Woodruff. The failure of the Tacoma Narrows
Bridge. Technical report, Federal Works Agency, 1941.
[12] A. C. Antoulas. Approximation of Large-Scale Dynamical Systems. Society for Industrial and
Applied Mathematics, 2005.
[13] E. Asarin, T. Dang, G. Frehse, A. Girard, C. Le Guernic, and O. Maler. Recent progress in
continuous and hybrid reachability analysis. In Proc. of the IEEE Conference on Computer
Aided Control Systems Design, pages 1582–1587, 2006.
[14] E. Asarin, T. Dang, and A. Girard. Hybridization methods for the analysis of nonlinear
systems. Acta Informatica, 43:451–476, 2007.
[15] S. Bak, S. Bogomolov, and M. Althoff. Time-triggered conversion of guards for reachability analysis of hybrid automata. In Proc. of the 15th International Conference on Formal
Modelling and Analysis of Timed Systems, 2017.

28

MATTHIAS ALTHOFF

[16] A. Banerjee, K. K. Venkatasubramanian, T. Mukherjee, and S. K. S. Gupta. Ensuring safety,
security, and sustainability of mission-critical cyber-physical systems. Proceedings of the
IEEE, 100(1):283–299, 2012.
[17] F. Blanchini. Set invariance in control. Automatica, 35(11):1747 – 1767, 1999.
[18] X. Chen, E. Ábrahám, and S. Sankaranarayanan. Flow*: An analyzer for non-linear hybrid
systems. In Proc. of Computer-Aided Verification, LNCS 8044, pages 258–263. Springer,
2013.
[19] Y. Chou, X. Chen, and S. Sankaranarayanan. A study of model-order reduction techniques
for verification. In Proc. of Numerical Software Verification, pages 98–113, 2017.
[20] A. Chutinan and B. H. Krogh. Computational techniques for hybrid system verification. IEEE
Transactions on Automatic Control, 48(1):64–75, 2003.
[21] C. Combastel. A state bounding observer for uncertain non-linear continuous-time systems
based on zonotopes. In Proc. of the 44th IEEE Conference on Decision and Control, and the
European Control Conference, pages 7228–7234, 2005.
[22] T. Dang. Vérification et synthèse des systèmes hybrides. PhD thesis, Institut National Polytechnique de Grenoble, 2000.
[23] T. Dang, O. Maler, and R. Testylier. Accurate hybridization of nonlinear systems. In Hybrid
Systems: Computation and Control, pages 11–19, 2010.
[24] A. G. Davenport and J. P. C. King. The influence of topography on the dynamic wind loading
of long span bridges. Journal of Wind Engineering and Industrial Aerodynamics, 36(2):1373–
1382, 1990.
[25] P. S. Duggirala and M. Viswanathan. Parsimonious, simulation based verification of linear
systems. In Proc. of International Conference on Computer Aided Verification, pages 477–
494, 2016.
[26] A. El-Guindy, D. Han, and M. Althoff. Formal analysis of drum-boiler units to maximize the load-following capabilities of power plants. IEEE Transactions on Power Systems,
31(6):4691–4702, 2016.
[27] A. El-Guindy, D. Han, and M. Althoff. Estimating the region of attraction via forward reachable sets. In Proc. of the American Control Conference, pages 1263–1270, 2017.
[28] P. Falcone, M. Ali, and J. Sjöberg. Predictive threat assessment via reachability analysis and
set invariance theory. IEEE Transactions on Intelligent Transportation Systems, 12(4):1352–
1361, 2011.
[29] G. Frehse, C. Le Guernic, A. Donzé, S. Cotton, R. Ray, O. Lebeltel, R. Ripado, A. Girard,
T. Dang, and O. Maler. SpaceEx: Scalable verification of hybrid systems. In Proc. of the
23rd International Conference on Computer Aided Verification, LNCS 6806, pages 379–395.
Springer, 2011.
[30] G. Frehse, B. H. Krogh, and R. A. Rutenbar. Verifying analog oscillator circuits using forward/backward abstraction refinement. In Proc. of Design, Automation and Test in Europe,
2006.
[31] A. Girard. Reachability of uncertain linear systems using zonotopes. In Hybrid Systems:
Computation and Control, LNCS 3414, pages 291–305. Springer, 2005.
[32] A. Girard and C. Le Guernic. Efficient reachability analysis for linear systems using support
functions. In Proc. of the 17th IFAC World Congress, pages 8966–8971, 2008.
[33] A. Girard, C. Le Guernic, and O. Maler. Efficient computation of reachable sets of linear
time-invariant systems with inputs. In Hybrid Systems: Computation and Control, LNCS
3927, pages 257–271. Springer, 2006.
[34] A. Girard and G. J. Pappas. Approximate bisimulation relations for constrained linear systems. Automatica, 43:1307–1317, 2007.
[35] A. Girard and G. J. Pappas. Approximate bisimulation: A bridge between computer science
and control theory. European Journal of Control, 17(5-6):568–578, 2011.
[36] A. Gurung, A. Deka, E. Bartocci, S. Bogomolov, R. Grosu, and R. Ray. Parallel reachability
analysis for hybrid systems. In Proc. of the ACM/IEEE International Conference on Formal
Methods and Models for System Design, pages 12–22, 2016.
[37] Z. Han and B. Krogh. Reachability analysis of hybrid control systems using reduced-order
models. In Proc. of the American Control Conference, pages 1183–1189, 2004.
[38] Z. Han and B. H. Krogh. Reachability analysis of large-scale affine systems using lowdimensional polytopes. In Hybrid Systems: Computation and Control, LNCS 3927, pages
287–301. Springer, 2006.

REACHABILITY ANALYSIS WITH UNCERTAIN INPUTS IN THE KRYLOV SUBSPACE 29

[39] Marlis Hochbruck and Christian Lubich. On Krylov subspace approximations to the matrix
exponential operator. SIAM Journal on Numerical Analysis, 34(5):1911–1925, 1997.
[40] L. Jaulin, M. Kieffer, and O. Didrit. Applied Interval Analysis. Springer, 2006.
[41] I. Kolmanovsky and E. G. Gilbert. Theory and computation of disturbance invariant sets for
discrete-time linear systems. Mathematical Problems in Engineering, 4:317–367, 1998.
[42] W. Kühn. Rigorously computed orbits of dynamical systems without the wrapping effect.
Computing, 61:47–67, 1998.
[43] A. A. Kurzhanskiy and P. Varaiya. Ellipsoidal techniques for reachability analysis of discretetime linear systems. IEEE Transactions on Automatic Control, 52(1):26–38, 2007.
[44] G. Lafferriere, G. J. Pappas, and S. Yovine. A new class of decidable hybrid systems. In
Hybrid Systems: Computation and Control, LNCS 1569, pages 137–151. Springer, 1999.
[45] V. T. H. Le, C. Stoica, T. Alamo, E. F. Camacho, and D. Dumur. Zonotopes: From Guaranteed State-estimation to Control. Wiley, 2013.
[46] D. Limon, I. Alvarado, T. Alamo, and E. F.Camacho. Robust tube-based MPC for tracking of
constrained linear systems with additive disturbances. Journal of Process Control, 20(3):248–
260, 2010.
[47] I. M. Mitchell, A. M. Bayen, and C. J. Tomlin. A time-dependent Hamilton–Jacobi formulation of reachable sets for continuous dynamic games. IEEE Transactions on Automatic
Control, 50:947–957, 2005.
[48] C. Moler and C. Van Loan. Nineteen dubious ways to compute the exponential of a matrix,
twenty-five years later. SIAM Review, 45(1):3–49, 2003.
[49] A. Pereira and M. Althoff. Overapproximative arm occupancy prediction for human-robot
interaction built from archetypal movements. In International Conference on Robotics and
Automation, pages 1394–1401, 2016.
[50] Michele Pisaroni. Dynamic response of a 3D structure: Roosevelt lake bridge. Technical
University of Munich, June 2012.
[51] H. Roehm, J. Oehlerking, M. Woehrle, and M. Althoff. Reachset conformance testing of
hybrid automata. In Proc. of Hybrid Systems: Computation and Control, pages 277–286,
2016.
[52] Y. Saad. Analysis of some Krylov subspace approximations to the matrix exponential operator. SIAM Journal on Numerical Analysis, 29(1):209–228, 1992.
[53] Y. Saad. Iterative Methods for Sparse Linear Systems. Society for Industrial and Applied
Mathematics, 2003.
[54] W. H. A. Schilders, H. A. van der Vorst, and J. Rommes. Model Order Reduction: Theory,
Research Aspects and Applications. Springer, 2008.
[55] B. Schürmann and M. Althoff. Convex interpolation control with formal guarantees for disturbed and constrained nonlinear systems. In Proc. of Hybrid Systems: Computation and
Control, pages 121–130, 2017.
[56] R. B. Sidje. Expokit: A software package for computing matrix exponentials. ACM Transactions on Mathematical Software, 24(1):130–156, 1998.
[57] H. Täubig, U. Frese, C. Hertzberg, C. Lüth, S. Mohr, E. Vorobev, and D. Walter. Guaranteeing functional safety: design for provability and computer-aided verification. Autonomous
Robots, 32(3):303–331, 2012.
[58] H.-D. Tran, L. V. Nguyen, W. Xiang, and T. T. Johnson. Order-reduction abstractions
for safety verification of high-dimensional linear systems. Discrete Event Dynamic Systems,
27(2):443–461, 2017.
[59] A. J. van der Schaft. Equivalence of dynamical systems by bisimulation. IEEE Transactions
on Automatic Control, 49(12):2160–2172, 2004.
[60] J. G. VanAntwerp and R. D. Braatz. A tutorial on linear and bilinear matrix inequalities.
Journal of Process Control, 10:363–385, 2000.
[61] H. N. Villegas Pico and D. C. Aliprantis. Voltage ride-through capability verification of wind
turbines with fully-rated converters using reachability analysis. IEEE Transactions on Energy
Conversion, 29(2):392–405, 2014.
[62] H. Wang and Q. Ye. Error bounds for the Krylov subspace methods for computations of
matrix exponentials. arXiv:1603.07358v1, 2016.
[63] S. Williams, L. Oliker, R. Vuduc, J. Shalf, K. Yelick, and J. Demmel. Optimization of sparse
matrix-vector multiplication on emerging multicore platforms. Parallel Computing, 35:178–
194, 2009.

30

MATTHIAS ALTHOFF

Department of Computer Science, Technische Universität München, 85748 Garching, Germany
E-mail address: althoff@tum.de

