arXiv:1506.05607v5 [] 22 Aug 2017

1

Unbounded-Time Analysis of Guarded LTI
Systems with Inputs by Abstract
Acceleration
Dario Cattaruzza Alessandro Abate Peter Schrammel Daniel Kroening

Abstract—Reachability analysis of continuous and
discrete time systems is a hard problem that has
seen much progress in the last decades. In many
cases the problem has been reduced to bisimulations
with a number of limitations in the nature of the
dynamics, soundness, or time horizon. In this article
we focus on sound safety verification of UnboundedTime Linear Time-Invariant (LTI) systems with inputs
using reachability analysis. We achieve this by using
Abstract Acceleration, which over-approximates the
reach tube of a system over unbounded time by using
abstraction . The technique is applied to a number of
models and the results show good performance when
compared to state-of-the-art tools.

Reachability in linear programs, however, is a
formidable challenge for automatic analysers: the
problem is undecidable despite the restriction to
linear transformations (i.e., linear dynamics) and
linear guards.
The goal of this article is to push the frontiers
of unbounded-time reachability analysis: we aim at
devising a method that is able to reason soundly
about unbounded trajectories. We present a new
approach for performing abstract acceleration. Abstract acceleration [26], [27], [34] approximates the
effect of an arbitrary number of loop iterations (up to

I. Introduction
Linear loops are an ubiquitous programming pat-

infinity) with a single, non-iterative transfer function
that is applied to the entry state of the loop (i.e., to

tern. Linear loops iterate over continuous variables, the set of initial conditions of the linear dynamics).
which are updated with a linear transformation. Lin- This article extends the work in [34] to systems
ear loops may be guarded, i.e., terminate if a given with non-deterministic inputs elaborating the details
linear condition holds. Inputs from the environment
can be modelled by means of non-deterministic
choices within the loop. These features make linear
loops expressive enough to capture the dynamics of

omitted in [39].
The key contributions of this article are:

many hybrid dynamical models. The usage of such
models in safety-critical embedded systems makes
linear loops a fundamental target for formal methods.
Many high-level requirements for embedded con-

1) We present a new technique to include inputs
(non-determinism) in the abstract acceleration

trol systems can be modelled as safety properties, i.e.
deciding reachability of certain bad states, in which

of general linear loops.
2) We introduce the use of support functions in

the system exhibits unsafe behaviour. Bad states may,
in linear loops, be encompassed by guard assertions.

complex spaces, in order to increase the precision of previous abstract acceleration methods.

2

Since the transformations A and B are linear, and
vector sums preserve convexity, the sets Xn =

II. Preliminaries
A. Linear Loops with Inputs
Simple linear loops are programs expressed in the
form:
while(Gx ≤ h) x := Ax + Bu,
ψ := Gx ≤ h is a linear constraint on the states (with
G ∈ Rr×p and h ∈ Rr ), u ∈ Rq is a non-deterministic

input, and A ∈ R p×p and B ∈ R p×q are linear transformations characterising the dynamics of the system.
In particular, the special instance where ψ = ⊤ (i.e.,

“while true”) represents a time-unbounded loop with
no guards, for which the discovery of a suitable

invariant (when existing) is paramount. As evident
at a semantical level, this syntax can be interpreted
as the dynamics of a discrete-time LTI model with
inputs, under the presence of a guard set which, for
ease of notation, we denote as G = {x | Gx ≤ h}.
In the remaining of this work we will also use the
notation Mi,∗ to represent the rows of a matrix and
M∗, j its columns.
B. Model Semantics
The traces of the model starting from an initial
set X0 ⊆ R p , with inputs restricted to U ⊆ Rq , are
u0
u1
u2
sequences x0 −→ x1 −→ x2 −→ . . ., where x0 ∈ X0
and ∀k ≥ 0, xk+1 = τ(xk , uk ), where

as the union of the reachable sets over n iterations.
S
Moreover, X̂ = n≥0 τn (X0 , U) extends the previous
notion over an unbounded time horizon.
C. Support Functions
1) Support Function Definition: A support function is a convex function on R p which describes
the distance of a supporting hyperplane for a given
geometrical set in R p .
Support functions may be used to describe a set by
defining the distance of its convex hull with respect
to the origin, given a number of directions. More
specifically, the distance from the origin to the hyperplane that is orthogonal to the given direction and
that touches its convex hull at its farthest. Finitely
sampled support functions are template polyhedra in
which the directions are not fixed, which helps avoiding wrapping effects [25]. The larger the number of
directions provided, the more precisely represented
the set will be.
In more detail, given a direction v ∈ R p , the
support function of a non-empty set X ⊆ R p in the

direction of v is defined as
ρX : R p → R,
where x · v =
vectors.

states reached from X by τ in one step:
(2)

We furthermore denote the set of states reached from
X0 via τ in n steps (n-reach set), for n ≥ 0:
0

τ (X0 , U) = X0
n

(4)

(1)

We extend the notation above to convex sets of
states and inputs (X and U), and denote the set of

τ(X, U) = {τ(x, u) | x ∈ X, u ∈ U}

We define the n-reach tube
[
X̂n = τ̂n (X0 , U) =
τk (X0 , U)
k∈[0,n]

where x ∈ R p is a valuation on the state variables,

τ(xk , uk ) = { Axk + Buk | Gxk ≤ h ∧ uk ∈ U}

τn (X0 , U) are also convex.

Pp

i=0

ρX (v) = sup{x · v : x ∈ X} .
xi vi is the dot product of the two

Support functions do not exclusively apply to
convex polyhedra, but in fact to any set X ⊆ R p

represented by a general assertion θ(X). We will
restrict ourselves to the use of convex polyhedra, in
which case the support function definition translates
to solving the linear program

n−1

τ (X0 , U) = τ(τ

(X0 , U) ∩ G, U)

(3)

ρX (v) = max{x · v | Cx ≤ d} .

(5)

3

which is consistent with the real case when θ = 0.
The reason why eiθ cannot be extracted out is be2) Support Functions Properties: Several properties of support functions allow us to reduce operational complexity. The most significant are [23]:
ρkX (v) = ρX (kv) = kρX (v) : k ≥ 0
ρAX (v) = ρX (AT v) : A ∈ R p×p
ρX1 ⊕X2 (v) = ρX1 (v) + ρX2 (v)
ρX (v1 + v2 ) ≤ ρX (v1 ) + ρX (v2 )

ρconv(X1 ∪X2 ) (v) = max{ρX1 (v), ρX2 (v)}
ρX1 ∩X2 (v) ≤ min{ρX1 (v), ρX2 (v)}

As can be seen by their structure, some of these properties reduce complexity to lower-order polynomial
or even to constant time, by turning matrix-matrix
multiplications (O(p3 )) into matrix-vector (O(p2 )), or
into scalar multiplications.
3) Support Functions in Complex Spaces:
The literature does not state, as far as we found
any description of the use of support functions
in complex spaces. Since this is relevant to using
our technique, we extend the definition of support
functions to encompass their operation on complex
spaces.
A support function in a complex vector field
is a transformation:

cause it is a rotation, and therefore follows the same
rules as a matrix multiplication,
 

 cos θ sin θ  
iθ
 v .


ρX (e v) , ρX 
− sin θ cos θ

Since matrices using pseudo-eigenvalues are real, all
other properties remain the same. An important note
is that when using pseudo-eigenvalues, conjugate
eigenvector pairs must be also converted into two
separate real eigenvectors, corresponding to the real
and the imaginary parts of the pair.
III. The Polyhedral Abstract Domain
A. Convex Polyhedra
A polyhedron is a topological element in R p
with flat polygonal (2-dimensional) faces. Each face
corresponds to a hyperplane that creates a halfspace,
and the intersections of these hyperplanes are the
edges of the polyhedron. A polyhedron is said to be
convex if its surface does not intersect itself and a
line segment joining any two points of its surface is
contained in the interior of the polyhedron. Convex
polyhedra are better suited than general polyhedra
as an abstract domain, mainly because they have a
simpler representation and operations over convex

ρX (v) : C p → R = sup{|x · v| | x ∈ X ⊆ C p , v ∈ C p }.

polyhedra are in general easier than for general

The dot product used here is the Euclidean Internal

polyhedra. There are a number of properties of
convex polyhedra that make them ideal for abstract

Product of the vectors, which is commonly defined
in the complex space as:
a·b=

p
X
i=0

ai bi ,

a, b ∈ C p

We are interested in the norm of the complex value,
which is a 1-norm given our definition of dot product:
|a · b| = |re(a · b)| + |im(a · b)|
Returning to our support function properties, we
now have:
ρX (reiθ v) = rρX (eiθ v),

interpretation of continuous spaces, including their
ability to reduce an uncountable set of real points
into a countable set of faces, edges and vertices.
Convex polyhedra retain their convexity across linear
transformations, and are functional across a number
of operations because they have a dual representation [21]. The mechanism to switch between these
two representations is given in section III-B5
1) Vertex Representation: Since every edge in the
polyhedron corresponds to a line between two vertices and every face corresponds to the area enclosed

4

by a set of co-planar edges, a full description of the
polyhedron is obtain by simply listing its vertices.
Since linear operations retain the topological properties of the polyhedron, performing these operations
on the vertices is sufficient to obtain a complete
description of the transformed polyhedron (defined

dron is
V ′ = {v′ = v + t | v ∈ V}
Given an inequality representation X and a translation vector t, the transformed polyhedron corresponds to
X ′ = {x | Cx ≤ d + Ct}

by the transformed vertices). Formally, a polyhedron
is a set V ∈ R p such that v ∈ V is a vertex of the
polyhedron.

2) Linear Transformation: Given a vertex rep-

2) Inequality Representation: The dual of the
Vertex representation is the face representation. Each

resentation V and a linear transformation L, the
transformed polyhedron is

face corresponds to a bounding hyperplane of the

V ′ = LV

polyhedron (with the edges being the intersection
of two hyperplanes and the vertices the intersection

Given an inequality representation X and a linear

of 3 or more), and described mathematically as a
function of the vector normal to the hyperplane. If

transformation L, the transformed polyhedron corresponds to

we examine this description closely, we can see that
it corresponds to the support function of the vector
normal to the hyperplane. Given this description we
formalise the following: A convex polyhedron is a
p

o
n
X ′ = x | C(L+ )T x ≤ ρ′X (L+ CT )
where L+ represents the pseudo-inverse of L. In the
case when the inverse L−1 exists, then

topological region in R described by the set
X = {x ∈ R p | Cx ≤ d, C ∈ Rm×p , d ∈ Rm }
where the rows Ci,∗ for i ∈ [1, m] correspond to the
transposed vectors normal to the faces of the polyhedron and di for i ∈ [1, m] their support functions.
For simplicity of expression, we will extend the use
of the support function operator as follows:
ρ′X : Rm×p → Rm

 ρX (MT )
1,∗

 ρX (MT )
2,∗
ρ′X (MT ) = 
..

.


ρX (MTm,∗ )











B. Operations on Convex Polyhedra
Several operations of interest can be performed on
convex polyhedra
1) Translation: Given a vertex representation V
and a translation vector t, the transformed polyhe-

n
o
X ′ = x | C(L−1 )T x ≤ d)
From this we can conclude that linear transformations are better handled by a vertex representation,
except when the inverse of the transformation exists
and is know a-priori. This works makes use of this
last case to avoid continuous swapping in representations.
3) Set Sums: The addition of two polyhedra is a
slightly more complex matter. The resulting set is
one such that for all possible combinations of points
inside both original polyhedra, the sum is contained
in the result. This operation is commonly known as
the Minkowski sum, namely
A ⊕ B = {a + b | a ∈ A, b ∈ B}
Given two vertex representations V1 and V2 the
resulting polyhedron
V = conv(V1 ⊕ V2 )

5

where conv(·) is the convex hull of the set of vertices
contained in the Minkowski sum.

5) Vertex Enumeration: The vertex enumeration
problem corresponds to the algorithm required to

Let

obtain a list of all vertices of a polyhedron given an
inequality description of its bounding hyperplanes.
X1 = {x | C1 x ≤ d1 }

Given the duality of the problem, it is also possible
to find the bounding hyperplanes given a vertex de-

X2 = {x | C2 x ≤ d2 }

scription if the chosen algorithm exploits this duality.

be two sets, then

In this case the description of V is given in the forms
of a matrix inequality Vx ≤ [ 1 1 · · · 1 ]T with

X = X1 ⊕ X2 = {x | Cx ≤ d},

V = [ v1 · · · vm ]T , vi ∈ V. Similarly, A can be
described as a set containing each of its rows. At

where



 C1 
 d + ρ′X (CT1 )
2
 , d =  1
C = 
C2
d2 + ρ′X1 (CT2 )



 .

Because these sets correspond to systems of inequalities, they may be reduced removing redundant
constraints. Note that if C1 = C2 then
X = X1 ⊕ X2 = {x | C1 x ≤ d1 + d2 },
4) Set Hadamard Product:
Lemma 1. Given two vertex representations V1 and
V2 the resulting polyhedron
V = V1 ◦ V2 = conv({v = v1 ◦ v2 | v1 ∈ V1 , v2 ∈ V2 )
where ◦ represents the Hadamard (coefficient-wise)
product of the vectors, contains all possible combinations of products between elements of each set.
Proof: Given a convex set X, we have:

the time of writing, there are two algorithms that
efficiently solve the vertex enumeration problem. lrs
is a reverse search algorithm, while cdd follows the
double description method. In this work we use the
cdd algorithm for convenience in implementation
(the original cdd was developed for floats, whereas
lrs uses rationals). The techniques presented here can
be applied to either.
Let
C = {x | Ax ≥ 0, A ∈ Rn×p , x ∈ R p }
be the polyhedral cone represented by A. The pair
( A, V) is said to be a double description pair if
C = {λT V | V ∈ R p , λ ∈ R|V|
≥0 }
V is called the generator of X. Each element in V lies
in the cone of X, and its minimal form (smallest m)
has a one-to-one correspondence with the extreme

X ′ = {xi j | xi , x j ∈ X, xi j = txi +(1−t)x j , t ∈ [0, 1]} ⊆ X

rays of X if the cone is pointed (i.e., it has a vertex
at the origin). This last can be ensured by translating

Given xi ∈ X, y j ∈ Y, zi, j = xi ◦ y j ∈ Z

a polyhedral description so that it includes the origin,
and then translating the vertices back once they have

xi j ∈ X ′ , yk ∈ Y, zi,k , z j,k ∈ Z ⇒ zi j,k ∈ Z
xi j ∈ X, ym , yn ∈ Y, zi j,m , zi j,n ∈ Z ⇒ zi j,mn ∈ Z
This equation proves that given v11 , v12 ∈ V1 ,
v21 , v22 ∈ V2 and u, t ∈ [0, 1],

(tv11 + (1 − t)v12 ) (uv21 + (1 − u)v22 ) ∈ V

been discovered (see section III-B).
We will also point out that
n
{x | Ax ≤ b} = x′ | [ −A

 x
p
′
where x ∈ R and x = 
1

b ] x′ ≥ 0


 ∈ R p+1 .

o

The vertex enumeration algorithm starts by finding
a base CK which contains a number of vertices of

6

the polyhedron. This can be done by pivoting over a
number of different rows in A and selecting the feasi-

computing the reachset (3) using a function f such
that f (·) = τn (·). In the case of systems without

ble visited points, which are known to be vertices of
the polyhedron (pivoting p times will ensure at least

inputs, this equation is xn = An x0 . We will use this
property and others derived from it to calculate our

one vertex is visited if the polyhedron is non-empty). abstract matrices.
CK is represented by AK which contains the rows
used for the pivots. The base CK is then iteratively

expanded to CK+i by exploring the ith row of A until
CK = C. The corresponding pairs ( AK+i , VK+i ) are
constructed using the information from ( AK , VK ) as
follows:

(6)

Hi− = {x | Ai,∗ x < 0},

(7)

= {x | Ai,∗ x = 0}

(8)

Hi0

be the spaces outside inside and on the ith hyperplane
and
VK+
VK−

= {v j ∈
= {v j ∈

Hi+ },

Hi− },

VK0 = {v j ∈ Hi0 }

The basic steps required to evaluate a reach tube
using abstract acceleration can be seen in figure 1.
1) The process starts by doing eigendecomposition

Let AK ∈ RnK ×p , Ai,∗ ∈ R1×p , VK ∈ R p ,
Hi+ = {x | Ai,∗ x > 0},

B. Overview of the Algorithm

(9)
(10)
(11)

the existing vertices lying on each of these spaces.
Then

of the dynamics ( A) in order to transform the
problem into a simpler one.
2) A variety of off-the-shelf tools may be used,
but since larger problems require numerical
algorithms for scalability, a second step involves
upper-bounding the error in order to obtain
sound results. In such cases, all subsequent
steps must be performed using interval arithmetic.
3) The inverse of the generalised eigenvectors
must be calculated soundly.
4) The problem gets transformed into canonical
form by multiplying both sides of the equation

VK+i = VK+ ∪ VK− ∪ VK0 ∪ VKi
(12)

VKi = ( Ai,∗ v+ )v− − ( Ai,∗ v− )v+ | v− ∈ V − , v+ ∈ V +
(13)

For the proof see [21].
IV. Abstract Matrices in abstract acceleration
A. Acceleration Techniques
Acceleration of a transition system is a method
that seeks to precisely describe the transition

by S−1 :


′
∩ G′ + U ′ where
Xk′ = J Xk−1

X ′ = S−1 X, U ′ = S−1 BU and G′ = {x | GSx ≤ h}

5) We calculate the number of iterations as explained in section VI. If there are no guards, we
use n = ∞. It is worth noting that this number

need not be exact: if we overapproximate the
number of iterations, the resulting reachtube
will overapproximate the desired one.

relations over a number of steps using a concise
description of the overall transition between the first

6) we overapproximate the dynamics of the variable inputs (for parametric or no inputs this step

and final step. Namely, it looks for a direct formula
to calculate the postimage of a loop from the initial

will be ignored) using the techniques described
in section V-D

states of the loop. Formally, given the dynamics
in equation (1) an acceleration formula aims at

7) we calculate the abstract dynamics using the
techniques described in section V-A

7

A

1. Calculate Ŝ, Ĵ
Eigenspace

2. Restore
Soundness

3. Find
Inverse

S

S−1

J
5. Find Number
of Iterations

G′
X0′
U′

4. Transform
problem into
Eigenspace

J, n

G
X0
BU

X0′ , U ′

7. Get Abstract
6. Semispherical
Dynamics
Jb Approximation
S

X♯

10. Get
Reach Tube

X0′ , Ub′

J
X ′♯

9. Find Eigen
Reach Tube

V0′

8. Find
Vertices

Fig. 1. Block diagram describing the different steps used to
calculate the abstract reach tube of a system.

8) we evaluate the vertices of the combined inputinitial eigenspace to be used as source for the
reachtube calculation
9) we use a sound simplex algorithm to evaluate

made up by the corresponding eigenvectors. We can
then easily compute An = SJ n S−1 , where
 n

 J 1



..

J n = 
(15)
.


n 
Jr
 n  n−p +1 
 n n n−1
s
 λ s 1 λ s

...
p s −1 λ s




..


n n−1
n
.
λs
λs


n
1
J s∈[1,r] =  .

..
..
 .

.
.
 .



λns
(16)

The abstract matrix An is computed as an abstraction
over a set of vectors mk ∈ R p , k ∈ [1 n] of entries of
Jk.

Let I s = [ 1 0 · · · 0 ] ∈ R ps . The vector mk
is obtained by the transformation ϕ−1 :
mk = [ I1 J k1

···

Ir J kr ] ∈ R p

(17)

the convex set product of the abstract dynamics
(used as the tableau) and the initial set (whose

such that J k = ϕ(mk ).

vertices are used as the obejctive functions
alongside a set of template directions for the

If J is diagonal [34], then mk equals the vector
of powers of eigenvalues [λk1 , . . . , λkr ]. An interval

result).

abstraction can thus be simply obtained by computing the intervals [min{λ0s , λns }, max{λ0s , λns }], s ∈

10) since we have calculated our result in the
eigenspace, we transform the reachtube back
into the normal space by multiplying by S .
C. Computation of Abstract Matrices
We define the abstract matrix An as an over-

approximation of the union of the powers of the maS
trix Ak such that An ⊇ k∈[0,n] Ak and its application

to the initial set X0

X̂n♯ = An X0 ⊇ X̂n

(14)

Next we explain how to compute such an abstract
matrix. For simplicity, we first describe this computation for matrices A with real eigenvalues, whereas

[1, r]. We observe that the spectrum of the interval
matrix σ(An ) (defined as intuitively) is an overS
approximation of k∈[0,n] σ( Ak ).

In the case of the sth Jordan block J s with geometric non-trivial multiplicity p s (λi = λi−1 = . . .),
observe that the first row of J ns contains all (possibly)
distinct entries of J ns . Hence, in general, the vector
section ms is the concatenation of the (transposed)
 
 n  n−p +1 T
s
first row vectors λns , n1 λn−1
of
s , . . . , p s −1 λ s
J ns .

Since the transformation ϕ transforms the vector
m into the shape of (16) of J n , it is called a matrix
shape [34]. We then define the abstract matrix as

the extension to the complex case will be addressed
in Section IV-D. Similar to [34], we first have to com-

An = {S ϕ(m) S−1 | Φm ≤ f } ,

pute the Jordan normal form of A. Let A = SJS−1
where J is the normal Jordan form of A, and S is

where the constraint Φm ≤ f is synthesised from
intervals associated to the individual eigenvalues and

(18)

8

to their combinations. More precisely, we compute
polyhedral relations: for any pair of eigenvalues (or

To deal with complex numbers in eigenvalues and
eigenvectors, [34] employs the real Jordan form for
conjugate eigenvalues λ = reiθ and λ∗ = re−iθ (θ ∈

 cos θ
r 
− sin θ


sin θ 
 .
cos θ

Although this equivalence will be of use once we
evaluate the progression of the system, calculating
powers under this notations is often more difficult
than handling directly the original matrices with
complex values.
In Section IV-C, in the case of real
eigenvalues we
in the power

150
100
50

D. Abstract Matrices in Complex Spaces

[0, π]), so that


λ 0 


is replaced by
0 λ∗

n=5
200

λn2

binomials) within J, we find an over-approximation
of the convex hull containing the points
n
o
mk | k ∈ [1, n] ⊆ {m | Φm≤ f }

250

n=4
n=1
n=3
n=2

0
0

5

10

20

15

25

30

35

λn1
Fig. 2. Polyhedral faces from an R2 subspace, where (λn1 , λn2 )
so that λ1 =2, λ2 =3, 1≤n≤5. Bold purple lines represent supports
found by this article. The dotted grey and dashed red polytopes
show logahedral approximations (box and octagon) used in [34].
Note the scales (sloped dashed lines are parallel to the x=y
line, and dashed red polytope hides two small faces yielding an
octagon).

1) Positive Real Eigenvalues
The exponential curve is cut along the diagonal

have abstracted the entries
matrix J ns by ranges of

between the maximum and minimum eigenvalues to create a support function for the

eigenvalues [min{λ0s · · · λns }
max{λ0s · · · λns }].
In the complex case we can do something

corresponding hyperplane. A third point taken
from the curve is used to test the direction of

similar by
form λ s

the corresponding template vector. An arbitrary
number of additional hyperplanes are selected

[min{r0s · · · rns }

rewriting eigenvalues into polar
=
r s eiθs and abstracting by
max{r0s · · · rns }]ei[0

, min(θ s ,2π)]

.

V. General Abstract Acceleration with Inputs
A. Using Support Functions on Abstract Acceleration

by picking pairs of adjacent points in the curve
and creating the corresponding support function
as shown in Figure 2.
2) Complex Conjugate Eigenvalue pairs

As an improvement over [34], the rows in Φ and

In the case of complex conjugate pairs, the
eigenvalue map corresponds to a logarithmic

f (see (18)) are synthesised by discovering support
functions in these sets. The freedom of directions

spiral. In this case, we must first extract the
number of iterations required for a full cycle.

provided by these support functions results in an
improvement over the logahedral abstractions used in

For convergent eigenvalues, only the first n iter-

previous works (see Figures 2 - 5). The mechanism
by which this works follows the convex properties
of the exponential progression. There are four cases
cases to consider 1 :
1 these

explain in detail the procedure alluded in [8]

ations have an effect on the support functions,
while in the divergent case only the last n
iterations are considered. Support functions are
found for adjacent pairs, checking for the location of the origin point (first point for convergent eigenvalues, last for divergent eigenvalues).

9

0.8
2
0.6

n=3

apex
n=3

n=2

0.4

1.5

n=2

λn2

λn2

n=1
0.2
0

n=1

1

−0.2
−0.4
−0.6 −0.4 −0.2

0.5
0

0.2

0.4

0.6

0.8

0

0.2

Fig. 3. Polyhedral faces from an R2 complex conjugate subspace,
where (λn1 , λn2 ) so that λ1 =0.8 + 0.4i, λ2 =0.8 − 0.4i, 1≤n≤14. Bold
purple lines represent supports found by this article. The blue
dotted line shows the support function that excludes the origin
(n=1), which is replaced by the support function projecting from
said origin.

0.4

0.6

0.8

λn1

λn1

Fig. 4. Polyhedral faces from an R2 Jordan block subspace,
where (λn1 , λn2 ) so that λ1 =0.8, λ2 =0.8, 1≤n≤15. Bold purple lines
represent supports found by this article. The blue dotted line
shows the support function that excludes the origin (n=1), which
is replaced by the support function projecting from said origin.

account for both sides of the axis on the latIf the origin falls outside the support function,
we look for an interpolant point that closes the
spiral tangent to the origin. This last check is
performed as a binary search over the remaining
points in the circle (whose supporting planes
would exclude the origin) to achieve maximum
tightness (see Figure 3).
3) Equal Eigenvalues
When two eigenvalues are the same, the resulting support functions are those orthogonal to
the x = y plane, intersecting the square created
by the maximum and minimum values.
4) Jordan Blocks of size > 1
In the case of eigenvalues with geometric mul-

ter. These form mirror images that are merged
during the abstraction. To make matters simple,
we use the magnitude of a complex eigenvalue
and evaluate whether the dynamics are concave
or convex with respect to the mirroring plane.
See Figure 5.
Note that if both eigenvalues are negative and/or
conjugate pairs from a different pair, the mirror
image would be taken on both axes, resulting
in a hyperrectangle. For a tighter bound in the
purely convergent case, we find the convex hull
of a point cloud for a small time horizon and
merge it with the hyperrectangle for the infinite
time horizon thereon.

tiplicities > 1, the shape of the function is
similar to its corresponding unit size block. In

An additional drawback of [34] is that calculating
the exact Jordan form of any matrix is computa-

the convergent case, since the convexity can be

tionally expensive and hard to achieve for large-

sharp, it is important to find the apex of the
upper diagonals in order to minimise the over-

dimensional matrices. We will instead use numerical
algorithms in order to get an approximation of the

approximation. See Figure 4.
5) Negative Eigenvalues and mixed types
When mapping a positive real eigenvalue to a
complex conjugate or negative one, we must

Jordan normal form and account for numerical errors.
In particular, if we examine the nature of (14),
we find out that the numerical operations are not
iterative, therefore the errors do not accumulate with

10

vides a 100-fold plus improvement, which can make
a difference towards rendering verification practically

15

λ2 > |λ1 | > 1

10
λn2

λ2 = |λ1 |

feasible. For a full description on the numerical
processes described here see [9]

|λ1 | > λ2 > 1

5

B. Abstract Matrices in Support Functions
|λ1 | > 1 > λ2

Since we are describing operations using abstract
matrices and support functions, we briefly review the

0
−10

0

10

λn1

Fig. 5. Polyhedral faces from an R2 subspace, with different
convexities (note that the blue and orange plots are convex w.r.t.
the λn2 -axis, whereas the green and brown are concave). Dotted
purple lines represent supports for some of these layouts.

nature of these operations and the properties that the
support functions retain within this domain. Let X ∈

R p be a space and A ∈ R p×p an abstract matrix for
the same space. From the definition we have
[
A=
Sϕ(m)S−1 : Φm ≤ f
which leads to

time. We use properties of eigenvalues to relax f by
finding the maximum error in the calculations that
can be determined by computing the norm δmax =
|A − ŜĴSˆ−1 |, where Ĵ and Ŝ are the numerically



ρAX (v) = ρSϕ(m)S−1 X (v) = ρϕ(m)S−1 X ST v ,

(19)

where
n
o
ρϕX (v) = sup ρϕ (x ◦ v) | x ∈ X

(20)

interval matrices and all operations are performed
using interval arithmetic with outward rounding in

ρϕ (v) = sup{m · ϕ−1 (v) | Φm ≤ f }

(21)

order to ensure soundness. In the following we will
presume exact results and use the regular notation to

Here, x ◦ y is the Hadamard product, where
(x ◦ y)i = xi yi , and ϕ−1 (·) is the reverse operation

calculated eigenvalues and eigenvectors of A. The
notation above is used to represent the matrices as

describe the algorithms. The constraints Φm < f are
then computed by considering the ranges of eigenvalues λ s ±δmax (represented in Fig. 2 as the diameter of
the blue dots). The outward relaxation of the support
functions ( f ), which follows a principle similar to
that introduced in [22], reduces the tightness of the
over-approximation, but ensures the soundness of the
abstract matrix An obtained. It is also worth noting
that the transformation matrices into and from the

eigenspace will also introduce over-approximations
due to the intervals. One can still use exact arithmetic
with a noticeable improvement over previous work;
however, for larger-scale systems the option of using
floating-point arithmetic, while taking into account
errors and meticulously setting rounding modes, pro-

and

of ϕ(·) in order to align the elements on v with the
elements in m. In the case of conjugate pairs this

is equivalent to multiplying the vector section by
h 1 1 i
, and in the case of a Jordan block by an
1
−1
upper triangular matrix of all ones.
We may also define
ρAX (v) = sup {ρ aX (v), ∀a ∈ A}
n
o
= sup Sϕ(m)S−1 x · v, ∀x ∈ X
n
o
= sup ϕ(m)S−1 x · ST v, ∀x ∈ X
n
o
= sup ρϕ (S−1 x ◦ ST v), ∀x ∈ X .

(22)

In order to simplify the nomenclature we write
ρAX (v) = ρX (AT v).

(23)

11

following representation:

C. Acceleration of Parametric Inputs
Let us now consider the following overapproximation for τ on sets:

n−1
X
k=0

♯

τ (X0 , U) = AX0 ⊕ BU

(24)

Unfolding (3) (ignoring the presence of the guard set
G for the time being), we obtain
Xn = An X0 ⊕

X

Ak BU

k∈[0,n−1]

What is left to do is to further simplify the sum
P
k
k∈[0,n−1] A BU. We can exploit the following simple results from linear algebra.
Lemma 2. If I − A is invertible, then
n−1
X
k=0





λ =


k

n
1−λn
1−λ

λ=1
λ,1

⇒ (I − An )(I − A)−1 = SDn S−1
−1k 1 − λni
d(λi , n, k) =
k + 1 (1 − λi )k+1
! n− j−1
k
X
λi
−1k− j n
+
k
−
j
j
−
1
(1
− λi )k− j
j=1



0
j<i






n
i = j ∧ λi = 1




1−λni



i = j ∧ λi , 1

1−λi
Dni, j = 


0
gm(λi ) = 1



n+1




λi = 1


k+1



 d(λi , n, j − i)
λi , 1

(25)

where gm(·) is the geometric multiplicity of the given
Ak = (I − An )(I − A)−1

eigenvalue.

. If furthermore lim An = 0, then
n→∞

lim

n→∞

n
X
k=0

D. Acceleration of Variable Inputs
k

A = (I − A)

−1

The result in the previous section can be only

.
This lemma presents a difficulty in the nature of A.
The inverse (I − A)−1 , does not exist for eigenvalues

directly applied under restricted conditions in the
case of variable inputs. For instance whenever ∀k >

0, uk = uk−1 . In order to generalise it (in particular to
non-constant inputs), we will over-approximate BU

of 1, i.e. we need 1 < σ( A), where σ(A) is the over the eigenspace by a semi-spherical enclosure
spectrum (the set of all the eigenvalues) of matrix A. with centre u′c and radius U ′ . To this end, we first
b
In order to overcome this problem, we introduce the rewrite
eigen-decomposition of A = SJS−1 , (and trivially
U J′ = S−1 BU = {u′c } ⊕ Ud′ ,
I = SIS−1 ), and by the distributive and transitive
properties we obtain
where u′c is the centre of the interval hull of U J′ :
n

(I − A )(I − A)

−1

n

−1 −1

= S(I − J )(I − J) S

.



1
 1
= (ρU′J (vi ) + ρU′J (−vi )) | vi j = 
 0

2

j=i
.
This allows us to accelerate the eigenvalues individj,i
Pn−1 k
ually, using the property k=0 1 = n for eigenvalues
of λ = 1. Using the properties above, and trans- We then over-approximate Ud′ via Ub′ , by the maxiu′ci

lating the problem into the generalised eigenspace
to accounting for unit eigenvalues, we obtain the

mum radius in the directions of the complex eigenvalues as (cf. illustration in Figure 6). Let

12

40

Λ = {λi | i ∈ [1, p], λ∗i , λi−1 }
λi < Λ
λi , λ∗i+1

of a vector by removing the elements where λi < Λ.
Extending this to matrices we have
Fb : R

o×pb

→R

Fb (C) = Cb where (Cb )i,∗ = fb (Ci,∗ )

0
0, 0

λi = λ∗i+1

and red(·) is a function that reduces the dimension

o×p

u′c
λ2




0





|v |
fb (v) = red(vb ) where (vb )i = 

q i



2

 vi + v2i+1

Ud′

20

fb : R p → R pb such that

U′

−20
−40

Ub′
−40

−20

0

20

40

λ1

(26)

Finally
Fig. 6. Relaxation of an input set within a complex subspace, in
order to make it invariant to matrix rotations. Dashed lines and
curves denote translated quantities onto the origin.

Ud′ = {u | C′u u ≤ d′u }

Ud′ ⊆ Ub′ = {u | Fb (C′u ) fb (u) ≤ fb (d′u )}

BU ⊆ Ub ⊕ Uc | Ub = SUb′ and Uc = {Su′c }

(27)

Since the description of Ub′ is no longer polyhedral
in R p , we will also create a semi-spherical overapproximation J b of J in the directions of the complex eigenvectors, in a similar way as we generated
Ub′ for Ud′ . More precisely,


 J b1



..
 where
J b = 
.



n 
J br



 λ j ∈ J b s = |λi | ∈ J s ∩ Λ
∀s ∈ [1 r] 

 gm(J b s ) = gm(J s )

(31)

Shifting our attention from reach sets to tubes,
we can now over-approximate the reach tube by

Theorem 1. The abstract acceleration
(28)

Jordan block.
Definition 1. Given a matrix A = SJS−1 and a
vector x, we define the following operations:


Fb∗ ( A, x) = S fb−1 Fb (J) fb (S−1 x)
(29)
′ 
′
−1
Fb ( A, x) = fb Fb (J) fb (x )
(30)
Finally, we refer to the accelerated sets
o
n
Ubn = Fb∗ ((I − An ), Fb∗ ((I − A)−1 , u)) | u ∈ Ub
n
Ucb
= Ucn ⊕ Ubn

n
Xn ⊆ An X0 ⊕ Ucb

abstract acceleration of the three summands in (31),
as follows.

where gm(·) is the geometric multiplicity of the

Ucn = (I − An )(I − A)−1 Uc

Returning to our original equation for the n-reach
set, we obtain2

τ♯n (X0 , U) =def An X0 ⊕ Bn Uc ⊕ Bnb Ub

(32)

is an over-approximation of the n-reach tube, namely
X̂n ⊆ τ♯n (X0 , U).
Proof: The proof is derived from that in [34]
for An X0 , and extends it as in the developments
presented above.
E. Combining Abstract Matrices
One important property of the abstract matrices
An , Bn and Bnb is that they are correlated. In the
that ∀U b′ , U c′ , U d′ ; ∃U b , U c , U d : U b′ = S−1 BU b so
that
= S−1 BU c and U d′ = S−1 BU d . Hence, this inclusion is
also valid in the original state space.
2 Note

U c′

13

case of parametric inputs, this correlation is linear
and described by the acceleration defined in Lemma

with J b defined by (28). This model provides a
tighter over-approximation than (32) since the accel-

(2). In the case of Bnb this relationship is not linear
(see Eq. 27). However, we can still find a linear over-

erated dynamics of the inputs are now constrained
by the accelerated dynamics of the system.

approximation of the correlation between Bnb and An
based on the time steps k. Given two orthonormal
p

q

spaces X ∈ R ∧ U ∈ R and a transition equation

The most important task remaining is how to

Xk+1 = AXk + BU,

calculate the number of iterations dealing with the
presence of the guard set G.

which is related to

Given a convex polyhedral guard expressed as the
assertion {x | Gx ≤ h}, we define Gi,∗ as the ith row

ρXk+1 (v) = ρ AXk (v) + ρ BU (v),
we define a space
 





 x 

′
X =
|
x
∈
X,
u
∈
U







 Bu


of G and hi as the corresponding element of h. We
denote the normal vector to the ith face of the guard
as gi = GTi . The distance of the guard to the origin
is thus γi = | hgi | .

so that

with

i

Given a convex set X, we may now describe
its position with respect to each face of the guard

 T 


 A v
 = ρX ′ DT v′ ,
ρXk+1 (v) = ρXk′ 
k
v

 A
D = 
0


0
 ,
I

Accelerating Xk+1 , we obtain

through the use of its support function alongside
the normal vector of the hyperplane (for clarity, we

 
v
v =   .
v

assume the origin to be inside set X):

′

ρX ( gi ) ≤ γi ,





ρXn (v) = ρ An X0 (v) + ρ(I− An )(I− A)−1 BU (v) = ρX0′ DnT v′ ,
with
 n
 A
n
D = 
0

0
(I − An )(I − A)−1





in the case of parametric inputs. More generally, the
diagonal elements of Dn correspond to the diagonal
P
k
elements of An and n−1
k=0 A B, which means we can
construct


An 0 
n
(33)
D = 
 | ρXn (v) = ρX0′ (DnT v′ ).
0 Bn

We can then apply this abstraction to (27) and obtain:
ρXn (v) =

′
ρX0′ (DnT
b v )


An
n
Db = 
0

VI. Abstract Acceleration with Guards:
Estimation of the number of Iterations

where



 v 
0 
 , v′ = 

Bnb
fb (v)



Bnb = SFb−1 (I − Jbn )(I − J b )−1 Fb (S−1 )

(34)

inside the hyperplane,

−ρX (−gi ) ≥ γi ,

outside the hyperplane.

Applying this to equation (31) we obtain:
ρXn ( gi ) = ρX0 ( Ani T gi ) + ρUcbn ( gi ) ≤ γi

(35)

ρXn (−gi ) = ρX0 (−A

(36)

ni T

gi ) + ρUcbn (−gi ) ≤ −γi

From the inequalities above we can determine
up to which number of iterations ni the reach tube
remains inside the corresponding hyperplane, and
starting from which iteration ni the corresponding
reach set goes beyond the guard:
In order for a reach set to be inside the guard
it must therefore be inside all of its faces, and we
can ensure it is fully outside of the guard set when
it is fully beyond any of them. Thus, we have n =
min{ ni }, and n = min{ ni }.

We have not however discussed why these two

cases are important. Looking at the transition in
equation (1), we can easily derive that if Gxk  h

14

the postimage of all subsequent iterations is empty.
Therefore, any overapproximation henceforth will

diagonal. In such a case, the progression of the
support function in each direction is monotonically

only add imprecision. We will use the bounds n and
n to create a tighter overapproximation. Let

increasing (or decreasing) and it is therefore very
easy to find a bound for its progression. We note

X̂n♯ = An X0 ⊕ Bn U

that the envelope of rotating dynamics will always
contain the true dynamics and is therefore a sound

(n-reachtube)

Xn♯ = An X0 ⊕ Bn U
(n-reachset)


X̂n♯ | n = τ An−n−1 Xn♯ ⊕ Bn U ∩ G, U
X̂n♯

X̂n♯ | n

∪

overapproximation. We will initially assume that γi
is positive and then extend to the general case.

X̂n♯

(37)
o
n
This double step prevents the set x | x ∈ X̂n♯ , x < Xn♯
to be included in further projections, thus reducing
=

Let ρX0 ( AnT gi ) = ρX0′ (J nT g′i ) such that
g′i = S−1 gi
X0 = {x | CX0 x ≤ dX0 }

the size of the overapproximation.

X0′ = S−1 X0 = {x | SC X0 x ≤ dX0 }

Computing the maximum ni such that (35) is
satisfied is not easy because the unknown ni occurs
in the exponent of the equation. However, since an
intersection with the guard set will always return a
sound over-approximation, we do not need a precise
value. We can over-approximate it by decomposing
gi into the generalised eigenspace of A. Let gi =
Pp
−1
or
j=1 ki j v j +res( gi ), where v j are row vectors of S
−1
−S such that ki j ≥ 0, and res( gi ) is the component

Let
Λσ = {λi : i ∈ [1, p],
p

fσ (v) : R → R

i−1
^
(λ∗i , λ j ∧ λi , λ j )}
j=1

pb

of gi that lies outside the range of S. Notice that
since S has an inverse, it is full rank and therefore

fσ (v) = red(vσ )



0


 r
P
where (vσ )i = 

v2j



j∈[1,p]∧(λ j =λi ∨λ j =λ∗i )

res( gi ) = 0 and subsequently not relevant. It is also

Fσ (C) = Cσ where (Cσ )i,∗ = fσ (Ci,∗ ).

important to note that S is the matrix of generalised
eigenvectors of A and therefore we are expressing
our guard in the generalised eigenspace of A.
 p

p
X
 X


nT
nT

ρX0 ( A gi ) = ρX0  ki j A v j  ≤
ki j ρX0 AnT v j
j=1

j=1

(38)

A. Overestimating the Iterations of a loop without
inputs
We start by looking into the approximation of
the inside bound (i.e. the iterations for which the
reachtube remains fully inside the guard). Since
rotating dynamics and Jordan shapes will have a
complex effect on the behaviour of the equation, we
seek to transform the Jordan form into a real positive

λi < Λσ
λi ∈ Λσ

Fσ : Ro×p → Ro×r

and red(·) is a function that reduces the dimension of a vector by removing the elements where
λi < Λσ . This reduction is not extrictly necessary, but it enables a faster implementation by reducing dimensionality. Correspondingly, given J =
diag ({J s | s ∈ [1, r]})

 σ1

 0
J σ = 
 0

0

0
σ2
0
0

···
···
..
.
···







0 

σr
0
0

(39)

where σ s = ||J s ||2 is the maximum singular value

(hence the induced norm [36]) of the Jordan block
J s.

15

Finally, let
x′c =




 1
vi j = 

 0

1
(ρX ′ (vi ) + ρX0′ (−vi )),
2 0

j=i
j,i

Xσ′ = {x | Fσ (SC X0 ) fσ (x) ≤ fσ (dX0 − SCX0 x′c )}
′
′
) | Xcσ
= { fσ (x′c )} ⊕ Xσ′
X0′ ⊆ fσ−1 (Xcσ

(40)

Proof: This follows from the developments unfolded above. Notice that the sequence ni is monotonically increasing, before it breaks the inequality. As
such any local minimum represents a sound underapproximation of the number of loop iterations. Note
that in the case where γi ≤ 0 we must first translate

the system coordinates such that γi > 0. This is

and vσ = fσ (v).
Using eigenvalue and singular value properties, we
obtain ρX0 ( AnT v j ) ≤ σ j n ρXcσ ((vσ ) j ) | j ∈ [1, r], and
therefore:
p
X
ρX0 ( AnT gi ) ≤
ki j σnj ρXcσ ((vσ ) j ))
(41)
j=1

Since we have no inputs, ρUcn ( gi ) + ρUbn ( gi ) = 0,
hence we may solve for ni :
ρX0 ( Ani T gi ) ≤

p
X
j=1

n

ki j σ j i ρXcσ ((vσ ) j ) ≤ γi

(42)

simply done by replacing x′ = x + c and operating
over the resulting system where γi′ = ρ c ( gi ) + γi .
Mathematically this is achieved as follows: first
we get c by finding the center of the interval hull
of G (if G is open in a given direction we may pick
any number in that direction for the corresponding
row of c). Next we

 
 xk   A Ac

 = 
1
0
1

where

To separate the divergent element of the dynamics
from the convergent one, let us define


ki j = max ki j ρXcσ ((vσ ) j ) , 0

transform the dynamics into

 


  xk−1   B 
 x
 
+
 uk |  k−1
1
0
1


 


 x   G
G =
|


 1   0
′


Gc   xk−1
 
1
1



 ∈ G′


 
  h 


 ≤ 


1 

σ = max ({σ s | s ∈ [1, p]}) .

This step will allow us to track effectively which

B. Underestimating the Iterations of a loop without

trajectories are likely to hit the guard and when, since
it is only the divergent element of the dynamics that

inputs

can increase the reach tube in a given direction.
Replacing (42), we obtain
p
X

σj
σ

!n

In order to apply a similar techniques to (36) we
must find an equivalent under-approximation. In the
case of equation (42), the σ j esure that the equation

(43)

diverges faster than the real dynamics, hence the
iteration found is an upper bound to the desired

which allows to finally formulate an iteration scheme

iteration. In this case we want the opposite, hence
we look for a model where the dynamics diverge

σ

n

j=1

ki j

≤ γi ,

for approximating n.
Proposition 1. An iterative under-approximation of
the number of iterations n can be computed by
starting with ni = 0 and iterating over
 p
!ni 
X

σ
j
 , (44)
ni ≥ n = logσ (γi ) − logσ  ki j
σ
j=1

substituting ni = n on the right-hand side until we
meet the inequality.

slower. In this case it is easy to demonstrate that
λb j = |λ j | represents these slower dynamics.
ni T

ρX0 (−A

gi ) ≤

p
X
j=1

ki j λb j ni ρXcσ (−(vσ ) j ) ≤ −γi (45)

which reduces to
σ

n

p
X
j=1

k−i j

λb j
σ

!n

+σ

n

p
X
j=1

k+i j ≤ −γi ,

(46)

16

where



k−i j = min ki j ρXcσ (−(vσ ) j ) , 0


k+i j = max ki j ρXcσ (−(vσ ) j ) , 0

cosines and removing the constants:

!
p
 ni X
((n − ni )θ j − αi j )2 

ci j 1 −
max σ
2
j=1
 p

X

⇒ min  ci j ((n − ni )θ j − αi j )2 
j=1
 p

X

2
2

⇒ min  ci j θ j (n − ni ) + ci j αi j θ j (n − ni )
j=1

An additional consideration must also be made
regarding the rotational nature of the dynamics. In
the previous case we did not care about the rotational
alignment of the set Xn with respect to the vector
gi , because any rotation would move the set inside
the guard. In this case, although the magnitude of
the resulting vector is greater than the required one,
the rotation may cause it to be at an angle that
keeps the set inside the guard. We must therefore
account for the rotating dynamics in order to find
the point where the angles align with the guard. In
order to do this, let us first fix the magnitudes of
the powered eigenvalues, in the case of convergent
dynamics we will assume they have converged a full

The solution to this equation is
Pp
j=1 ci j αi j θ j
| n ∈ [ni , ni + nθ ]
n = ni − P p
2 j=1 ci j θ2j

(47)

The second part of the equation is expected to be
a positive value. When this is not the case, the
dominating dynamics will have a rotation θ j ≥ pi2 .
In such cases we must explicitly evaluate the set of
up to 4 iterations after ni . If the resulting
 bound does
T
not satisfy the original inequality: ρX0 Ani gi ≥ γi ,

we replace ni = n until it does 3 .

Proposition 2. An iterative under-approximation of

the number of iterations n can be computed by
rotation in to make our equation strictly divergent.
starting with ni ′ = 0 and iterating over
Let θ = min{θ j | j ∈ [1, p]}, where θ j are the angles

 p
!ni ′ X
p

X
λ
2π
b
j
′
of the complex conjugate eigenvalues. Let nθ = θ
+
ni ≤ n = logσ (γi ) − logσ  k−i j
k+i j 
σ
j=1
j=1
be the maximum number of iterations needed for


′
T
any of the dynamics to complete a full turn. Then
(48)
ni = ni ′ + k | ρX0 A(ni +k) gi ≥ γi
ni +nθ
ni +n
≤ |λ j |
| |λi | ≤ 1, n ∈ 0nθ .
at any given turn |λ j |
where k is the result of equation (47). we substitute
This means that any bound we find on the iterations
for ni = n on the right-hand side until we break
will be necessarily smaller than the true value. Our
the inequality, and then find k such that the second
problem becomes the solution to:
inequality holds.


max σ

−1

ni

p
X
j=1



ci j cos((n − ni )θ j − αi j )

αi j = cos ( gi · v j )

 n

− λb j i


 ki j σ
ci j = 

 k−  λb j ni +nθ

ij σ

Since we are explicitly verifying the inequality,
there is no further proof required.
C. Estimating the Iterations of a loop with inputs
For the case with inputs, we will use the same

|λ j | ≥ 1
|λ j | < 1

The problem is simplified by underapproximating the

paradigm explained in the previous section after performing a mutation that transforms the system with
3 this

is a tighter value than work shown on previous versions
(2π)m
, where
of this paper where we overapproximated using nθ = Q
j θj
m is the number of conjugate pairs.

17

inputs into an over-approximating system without
inputs.
′
′
Let Xcσ
, Ucσ
be the corresponding sets of initial

states and inputs obtained by applying equation (40)
′
′
to X0′ and U J′ , and let U Jσ
= (I − J σ )−1 Ucσ
. The

accelerated resulting system may be represented by
the equations
′
′
⊕ (I − J nσ )U Jσ′
(Xcσ
)n = J nσ Xcσ




nT
′ ) (v) = ρX ′
J nT
ρ(Xcσ
σ v + ρU ′Jσ (v) − ρU ′Jσ J σ v
n
cσ

(49)

′
Xcσ
,u

Proof: Since the elements of the sums are
convergent, we have


ni ≥ nk ⇒ kni ≥ knk i.e. |kni | ≤ |knk |




⇒ logσ kni + kn ≥ logσ kni + kn
k

i

which means that nk in equation (51) is smaller than
or n in equation (44) (nk ≤ n ≤ ni | ni ≥ nk ).

In the case of equation (48), the explicit evaluation

of the guard at each cycle executes the behaviour
described here.

′
U Jσ
}

Let us now define (XU)σ = {x−u | x ∈
∈
which allows us to translate the system into


(50)
ρ((XU)′σ )n (v) = ρ(XU)′σ J nT
σ v
which has the same shape as the equations in the

E. Maintaining Geometric Multiplicity
A second step in optimising the number of itera-

previous section. We may now apply the techniques

tions comes from adding granularity to the bounding
semi-spherical abstraction by retaining the geometric

described above to find the bounds on the iterations.

multiplicity using the matrix J b .

D. Narrowing the estimation of the iterations

Lemma 3. Given a matrix A with eigenvalues
{λ s | s ∈ [1, r]}, where each eigenvalue λ s has a

The estimations above are very conservative, but
we may use further techniques to obtain tighter
bounds on the number of iterations. In the first
instance we note that we have eliminated all negative
terms in the sums in equation (44). Reinstating these

geometric multiplicity p s and corresponding generalised eigenvectors {v s,i | i ∈ [1, p s]},
∀n ≥ 0, An vis = λns v s,i +

terms can cause us to lose monotonicity, but we
may still create an iterative approach by fixing the
negative value at intermediate stages. Let ni be our
existing bound for the time horizon before reaching
 σ ni
 σ ni
P
P
a guard, and kn = pj=1 ki j σj , kni = pj=1 ki j σj
i

the corresponding negative and positive terms of the
equation. We may now find upper and lower bounds
for ni by replacing the equation


ni ≥ nk = logσ (γi ) − logσ kni + kn
k

(51)

where nk is the bound found in the previous stage.
Some stages of this process will provide an unsound
result, but they will also provide an upper bound to
our number of iterations. In fact, every second stage
will provide a monotonically increasing sound bound
which will be tighter than the one in equation (44).

i−1
X
j=1

j
λn−
s

j−1
Y
(n − k)v s,i− j
k=0



j−1
Q


(n
−
k)
i−1


X k=0

n
v
= λ s v s,i +
 (52)
s,i− j 
j


λs
j=1


Proof: By definition, given an eigenvector v s
of A, then Avs = λ s v s [32]. Similarly a generalised eigenvector v s,i of A satisfies the equation
( A − λ s I) v s,i = v s,i−1 and vs,1 = vs hence
Av s,i = λ s vs,i + v s,i−1
An v s,1 = λns v s,1
An v s,i = An−1 (λ s v s,i + v s,i−1 ) = λ s An−1 vs,i + An−1 v s,i−1
= λ2s An−2 v s,i + λ s An−2 v s,i−1 + An−1 vs,i−1
= · · · = λns v s,i +

n−1
X
j=0

λ sj An− j−1 v s,i−1

18

10,000

From here we recursively expand the formula for
An− j−1 v s,i−1 and obtain:

An v s,i = λns vs,i +

n−1
X

j−1
λ sj λn−
v s,i−1
s

+

j=0

n−1 X
n−2
X

λks An−k−2 v s,i−2

ρ( g)

1,000

n( fσ , fσ )

n( fb , fb )
n( fσ , fb )

100

n( fb , fb )

j=0 k=0

= λns vs,i + nλn−1
s v s,i−1 + n

n−2
X

λ sj An− j−2 vs,i−2

10

j=0

= · · · = λns v s,i +

i−1
X

j
λn−
s

j=1

j−1
Y
(n − k)v s,i− j

1

Let i′ denote the position of fb (λ j ) within the
block J bs it belongs to, such that its corresponding
generalised eigenvector is identified as vbs,i′ = fb (v j ).
Then
ρX0′ (J nT g′i )
pb
X


≤
ki j ρX0 J nb T fb (v j )
j=1



k−1
Q


′
(n − m)
pb


i −1
X
X
m=0


′ −k 
v
≤
ki j λb nj ρX0 vbs,i′ +
bs,i

k

λ

b
j
j=1
k=1



k−1
Q


′
(n
−
m)
p


i
−1
b
X
X m=0




′ −k 
v
ρ
≤
ki j λb nj ρX0 vbs,i′ +
bs,i
X

0
k

λ

b
j
j=1
k=1


≤

pb
X

ki′ j0 λb nj

j=1

+

i′
X

ki′ jm λb nj

p sY
−i′ −1
m=0

m=1

(n − m)

(53)

I order to manage the product on the right hand
side we use slightly different techniques for over- and
under-approximations. For ni we first find an upper
bound n′i using equation (44) and ki j = ki′ j 0 + ki′ j m

and then do a second iteration using ki j = ki′ j 0 +
′
p s −i
Q−1 ′
ki′ j m
(ni − m) which ensures the true value is
m=0

under the approximation. In the case of ni , we also
ki′ j 0

start with ki j =
iterative process.

+

ki′ j m

0

2

4

and update it during the

8

6
iteration

k=0

10

Fig. 7. Progression of the support function of a system for a
given guard. Blue dots are real values. The dashed green line overapproximates the progression using singular values (sec VI-A),
the dashed yellow line underapproximates them using eigenvalue
norms (sec VI-B), whereas the continuous purple lines represent
the tighter overapproximation maintaining the gemoetric multiplicity (sec VI-E). We can see how the purple line finds a better bound
for ni , while the ni bound is conservative for both approaches.
Mind the logarithmic scale.

Let us look at the following example:
 3 0 0 0 0 0 


 0 2 1 0 0 0 
 0 0 2 0 0 0 
J =  0 0 0 −2 0 0 


 0 0 0 0 −1 1 




S = 



0

0

1
0
0
0
0
0

0
3
−4
0
0
0




J σ = 

x′0 =



0

0
0
0
1
0
0

0
3
0
0

0
0
2
0

3
0
0
0

1

Gx ≤ 300 | G=

G= 1



−1

0

0
0
1
0
0
0

−1

0
0
0
0
1
1

0
0
0
0
0
1







0
0
0
√
2

1

1

1

1

3

−3

1

1

2









1
2
4

1
4
−3



1



 T
S

The progression of the system along the support
function and corresponding bounds as described in
the previous section are shown in figure 7

19

we assume that part of the system is coded, and
further assume that it is possible to discretise the

1,000,000
100,000

physical environment for simulation. Algorithm 1
shows a pseudo-code fragment for the temperature

ρ( g)

10,000
1,000
n( fb , fb )

control problem. We use the read function to

100

Algorithm 1 Temperature Control Loop
States: temp=temperature, heat=heat output.

n( fb , fb )

10
1

Inputs: set=set-point, amb=ambient temperature.
0

5

10

15

20

iteration
Fig. 8. Progression of the support function of a rotational system
for a given guard. Blue dots are real values (negative values are
missing due to the log scale). Continuous purple lines represent
the overapproximation. The steep vertical line at 19 is due to the
alignment of the rotations with the guard at this point. The point
at iteration 14 appears below the line because of the higher point
at iteration 9. The model will either find that this boundary was
met at iteration 9 or push it forward to 19.

1:

temp=5+read(35);

2:

heat=read(1);
while(temp< 400 && heat< 300)

3:
4:
5:

0

0

0

0

0
0
0
0
1.1e0.5i
0

0
0
0
0
0
1.1e−0.5i









we get the results in figure 8. In this case we can
see that the rotational dynamics force an increase of
the initially calculated iteration to account for the
effects of the rotation.
F. Case Study
We have selected a known benchmark to illustrate
the discussed procedure: the room temperature control problem [17]. The temperature (variable temp)
of a room is controlled to a user-defined set point
(set), which can be changed at any time through a
heating (heat) element, and is affected by ambient
temperature (amb) that is out of the control of the
system.
We formalise the description of such a system both
via a linear loop and via hybrid dynamics. Observe
that since such a system may be software controlled,

amb=5+read(35);

7:

set=read(300);
temp=.97 temp + .02 amb + .1 heat;

8:

heat=heat + .05 set;

6:

9:

Changing the eigenvalues to:
 2e−0.2i
0
0
0

 0
2e0.2i
0
0
√ −0.3i
 0
2e
0
0
√ 0.3i
J =  0
2e
0
0

 0
0
0
0

{

}

represent non-deterministic values between 0 and
the maximum given as argument. Alternatively, this
loop corresponds to the following hybrid dynamical
model:





 0.97 0.1   temp 
 temp 
 


 = 
heat k+1
−0.05 1
heat k



 0.02
0   amb 
 
 ,
+ 
0
0.05
set k

with initial condition

 temp

heat

 

  [5 40] 
 ∈ 
 ,
[0
1]
0

non-deterministic inputs

 
 amb   [5 40]
 ∈ 

set k
[0 300]

and guard set



 temp
G=


 heat

 
  1
 | 
0

0
1



 ,


 

  temp   400 


 < 
 
.


300 
heat

In this model the variables are continuous and
take values over the real line, whereas within the
code they are represented as long double precision floating-point values, with precision of ±10−19 ,

20

moreover the error of the approximate Jordan form
computation results in δmax < 10−17 . Henceforth we

300

focus on the latter description, as in the main text of
this work. The eigen-decomposition of the dynamics

200

A = S JS −1 ⊆ SJS−1 where


 0.798 ± 10−14 0.173 ± 10−15 

S = 
0 ± 10−19
0.577 ± 10−14


 0.985 ± 10−16 0.069 ± 10−17 


J = 
−0.069 ± 10−17 0.985 ± 10−16


 1.253 ± 10−12 −0.376 ± 10−13 
−1
 .

S = 
0 ± 10−18
1.732 ± 10−12

The discussed over-approximations of the reach-sets
indicate that the temperature variable intersects the
guard at iteration n = 32. Considering the pseudoeigenvalue matrix (described in the extended version

heat

is (the values are rounded to three decimal places):

heat = 300

100

0

0

100

200
temp

300

400

Fig. 9. The abstractly accelerated tube (yellow, dashed boundary),
representing an over-approximation of the thermostat reach tube
(dark blue). The set of initial conditions is shown in black,
whereas successive reach sets are shown in white. The guards
and the reach set that crosses them are close to the boundary in
red.

for the case of complex eigenvalues) along these
iterations, we use Equation (18) to find that the corre-

in Figure 9, where for the sake of clarity we display

sponding complex pair remains within the following
boundaries:

only 8 directions of the 16 constraints. This results
in a rather tight over-approximation that is not much

A32


 r
= 
−i


 r
B32 = 
−i

i
r

i
r








 

 














 

 








0.4144
0.0691
0.1082
0.9159
1
0
1
6.145

<
<
<
<
<
<
<
<

r
i
r+i
i−r
r
i
r+i
i−r

<
<
<
<
<
<
<
<

0.985
0.7651
1.247
0.9389
13.41
17.98
29.44
6.514

looser than the convex hull of all reach sets obtained
by [20] using the given directions. In Figure 9, we
can see the initial set in black colour, the collection
of reach sets in white, the convex hull of all reach
sets in dark blue (as computed by [20]), and finally
the abstractly accelerated set in light yellow (dashed
lines). The outer lines represent the guards.

The reach tube is calculated by multiplying these
abstract matrices with the initial sets of states and inputs, as described in Equation (32), by the following
inequalities:




 [5 40] 
 [5 40] 
32 
#
32
X̂32 =A 
 + B 

[0 1]
[0 300]


−24.76 <
temp
< 394.5

"
#



heat
< 253
temp  −30.21 <
=


−40.85
<
temp
+
heat
< 616.6
heat 



 −86.31 < temp − heat < 843.8

G. Calculating the Number of Iterations for Continuous Dynamics

Following the same steps as in Lemma 3 we
develop an equivalent for continuous dynamics.

The negative values represent the lack of restriction

Lemma 4. Given a matrix A with eigenvalues

in the code on the lower side and correspond to system cooling (negative heating). The set is displayed

{λ s | s ∈ [1, r]}, where each eigenvalue λ s has a
geometric multiplicity p s and corresponding gener-

21

alised eigenvectors {v s,i | i ∈ [1, p s]},
∀t ≥ 0, At vis = eλs v s,i +

i−1
X

Table II gives the comparison of our implementation using different levels of precision (long double,

t j eλs vs,i− j

256 bit, and 1024 bit floating-point precision) with
the original abstract acceleration for linear loops

j=1



i−1
X


t j v s,i− j 
= eλs v s,i +

(54)

without inputs (J) [34] (where inputs are fixed to
constants). This shows that our implementation gives

Proof: The proof derives again from the taylor
expansion.

tighter over-approximations on most benchmarks

j=1

∞
X

∞

X tk 

tk
λks v s,i + kλk−1
v s,i =
s v s,i−1
k!
k!
k=0
k=0
∞ k
∞ k
X
X
t k
t
=
λ s v s,i +
kλk−1
s v s,i−1
k!
k!
k=0
k=0

e At v s,i =

Ak

(column ‘improved’). While on a limited number of
instances the current implementation is less precise
(Fig. 2 gives a hint why this is happening), the
overall increased precision is owed to lifting the lim-

(55)

itation on directions caused by the use of logahedral
abstractions.
At the same time, our implementation is faster –

The rest of the proof follows the same expansion as
in 3

even when used with 1024 bit floating-point precision – than the original abstract acceleration (using

Given the similarity of equation (54) with (52) we
may apply exactly the same techniques described in

rationals). The fact that many bounds have improved
with the new approach, while speed has increased by

section VI-E to the continuous case.

several orders of magnitude, provides evidence of the
advantages of the new approach.
The speed-up is due to the faster Jordan form

= eλi t (vs,i + tv s,i−1 )

VII. Experimental Results
The algorithm has been implemented in C++ using the eigen-algebra package (v3.2), with double
precision floating-point arithmetic, and has been
tested on a 1.6 GHz core 2 duo computer.

computation, which takes between 2 and 65 seconds
for [34] (using the ATLAS package), whereas our
implementation requires at most one second. For the
last two benchmarks, the polyhedral computations
blow up in [34], whereas our support function approach shows only moderately increasing runtimes.

A. Comparison with other unbounded-time ap- The increase of speed is owed to multiple factors,
proaches.
as detailed in Table III. The difference of using long
In a first experiment we have benchmarked our double precision floating-point vs. arbitrary precision
implementation against the tools InterProc [33] and
Sting [12]. We have tested these tools on different

arithmetic is negligible, as all results in the given
examples match exactly to 9 decimal places. Note

scenarios, including guarded/unguarded, stable/un-

that, as explained above, soundness can be ensured
by appropriate rounding in the floating-point compu-

stable and complex/real loops with inputs (details
in Table I).4 It is important to note that in many

tations.

instances, InterProc (due to the limitations of widening) and Sting (due to the inexistence of tight

B. Comparison with bounded-time approaches.

polyhedral, inductive invariants) are unable to infer
finite bounds at all.

In a third experiment, we compare our method
with the LGG algorithm [28] used by SpaceEx [20].

4 The tool and
the benchmarks
http://www.cprover.org/LTI/.

In order to set up a fair comparison we have provided the implementation of the native algorithm

are

available

from

22

characteristics
improved
analysis time [sec]
name
type dim inputs bounds IProc Sti IProc
Sti
J+I
parabola i1
¬s,¬c,g 2
1
80
+25 +28 0.007 237
0.049
parabola i2
¬s,¬c,g 2
1
80
+24 +35 0.008 289
0.072
cubic i1
¬s,¬c,g 3
1
120
+44 +50 0.015 704
0.097
cubic i2
¬s,¬c,g 3
1
120
+35 +55 0.018 699
0.124
oscillator i0
s,c,¬g 2
0
56
+24 +24 0.004
0.990 0.021
oscillator i1
s,c,¬g 2
0
56
+24 +24 0.004
1.060 0.024
inv pendulum
s,c,¬g 4
0
16
+8 +8 0.009
0.920 0.012
convoyCar2 i0 s,c,¬g 3
2
12
+9 +9 0.007
0.160 0.043
convoyCar3 i0 s,c,¬g 6
2
24
+15 +15 0.010
0.235 0.513
convoyCar3 i1 s,c,¬g 6
2
24
+15 +15 0.024
0.237 0.901
convoyCar3 i2 s,c,¬g 6
2
24
+15 +15 0.663
0.271 1.416
convoyCar3 i3 s,c,¬g 6
2
24
+15 +15 0.122
0.283 2.103
type: s – stable loop, c – complex eigenvalues, g – loops with guard; dim: system dimension (variables); bounds: nb. of half-planes
defining the polyhedral set;
IProc is [33]; Sti is [12]; J+I is this work;
improved: number of bounds newly detected by J+I over the existing tools (IProc, Sti)
TABLE I
Experimental comparison of unbounded-time analysis tools with inputs

in [28]. We have run both methods on the con- ∼ 4 ms for each iteration on a 6-dimensional problem
voyCar example [34] with inputs, which presents with octagonal abstraction). This can be improved
an unguarded, scalable, stable loop with complex
dynamics, and focused on octahedral abstractions.

by the use of zonotopes, or by careful selection
of the directions along the eigenvectors, but this

For convex reach sets, the approximations computed
by abstract acceleration are quite tight in comparison

comes at a cost on precision. Even when finding
combinations that outperform our approach, this will

to those computed by the LGG algorithm. However,
storing finite disjunctions of convex polyhedra, the

only allow the time horizon of the LGG approach
to be slightly extended before matching the analysis

LGG algorithm is able to generate non-convex reach
tubes, which are arguably more proper in case of

time from abstract acceleration, and the reachable
states will still remain unknown beyond the extended

oscillating or spiralling dynamics. Still, in many
applications abstract acceleration can provide a tight

time horizon.

over-approximation of the convex hull of those nonconvex reach sets.
Table IV gives the results of this comparison. For
simplicity, we present only the projection of the

The evident advantage of abstract acceleration is
its speed over finite horizons without much precision
loss, and of course the ability to prove properties for
unbounded-time horizons.

bounds along the variables of interest. As expected,

C. Scalability

the LGG algorithm performs better in terms of
tightness, but its runtime increases with the number

Finally, in terms of scalability, we have an expected O(n3 ) complexity worst-case bound (from

of iterations. Our implementation of LGG using
Convex Polyhedra with octagonal templates is slower
than the abstractly accelerated version even for small
time horizons (our implementation of LGG requires

the matrix multiplications in equation (32)). We
have parameterised the number of cars in the convoyCar example [34] (also seen in Table II), and
experimented with up to 33 cars (each car after

23

characteristics
improved
analysis time (sec)
name
type
dim
bounds
tighter
looser
J (jcf)
mpfr+(jcf)
mpfr
ld
parabola i1
¬s,¬c,g
3
80
+4(5%)
0(0%)
2.51
( 2.49) 0.16 (0.06) 0.097 0.007
parabola i2
¬s,¬c,g
3
80
+4(5%)
0(0%)
2.51
( 2.49) 0.26 (0.06) 0.101 0.008
cubic i1
¬s,¬c,g
4
120
0(0%)
0(0%)
2.47
( 2.39) 0.27 (0.20) 0.110 0.013
cubic i2
¬s,¬c,g
4
120
0(0%)
0(0%)
2.49
( 2.39) 0.32 (0.20) 0.124 0.014
oscillator i0
s,c,¬g
2
56
0(0%)
-1(2%)
2.53
( 2.52) 0.12 (0.06) 0.063 0.007
oscillator i1
s,c,¬g
2
56
0(0%)
-1(2%)
2.53
( 2.52) 0.12 (0.06) 0.078 0.008
inv pendulum
s,c,¬g
4
12
+8(50%)
0(0%)
65.78
(65.24) 0.24 (0.13) 0.103 0.012
convoyCar2 i0
s,c,¬g
5
12
+9(45%)
0(0%)
5.46
( 4.69) 3.58 (0.22) 0.258 0.005
convoyCar3 i0
s,c,¬g
8
24
+10(31%)
-2(6%)
24.62
(11.98) 3.11 (1.01) 0.552 0.051
convoyCar3 i1
s,c,¬g
8
24
+10(31%)
-2(6%)
23.92
(11.98) 4.94 (1.01) 0.890 0.121
convoyCar3 i2
s,c,¬g
8
24
+10(31%)
-2(6%)
1717.00
(11.98) 6.81 (1.01) 1.190 0.234
convoyCar3 i3
s,c,¬g
8
24
+10(31%)
-2(6%)
1569.00
(11.98) 8.67 (1.01) 1.520 0.377
type: s – stable loop, c – complex eigenvalues, g – loops with guard; dim: system dimension (including fixed inputs); bounds: nb.
of half-planes defining the polyhedral set; improved: number of bounds (and percentage) that were tighter (better) or looser (worse)
than [34];
J is [34]; mpfr+ is this article using 1024bit mantissas (e < 10−152 ); mpfr uses a 256bit mantissa (e < 10−44 ); ld uses a 64bit
mantissa (e < 10−11 ); here e is the accumulated error of the dynamical system; jcf: time taken to compute Jordan form
TABLE II
Experimental comparison with previous work

the first requires 3 variables, so that for example
(33 − 1) × 3 = 96 variables), and have adjusted

the time bounded analysis is in most cases unsound
since it cannot reason about the unbounded time case

the initial states/inputs sets. We report an average
of 10 runs for each configuration. These results

(we not that a proof of the existence of a fix-point for
the given horizon would restore such soundness by

demonstrate that our method scales to industrial-size

many tools do not attempt to find such proof which

problems.

is left to the user). Unbounded-time solutions are
therefore preferred when such soundness is required,

# of variables
runtime (s)

3
6
12
24
48 96
0.004 0.031 0.062 0.477 5.4 56

although they are often either less precise or slower
than their bounded counterparts.

VIII. Related Work
There are several approaches that solve the safety
problem for the linear and other cases such as hybrid
systems. They are broadly divided into two categories due to the inherent nature of these. Namely
Optimisation
Eigen vs. ATLAS 5
Support functions vs. generators
long double vs. multiple precision arithmetic
interval vs. regular arithmetic
Total

Speed-up
2–10
2–40
5–200
.2–.5
4–80000

TABLE III
Performance improvements by feature

A. Time-Bounded Reachability Analysis
The first approach is to surrender exhaustive analysis over the infinite time horizon, and to restrict
the exploration to system dynamics up to some
given finite time bound. Bounded-time reachability
is decidable, and decision procedures for the resulting satisfiability problem have made much progress
in the past decade. The precision related to the
bounded analysis is offset by the price of uncertainty:
behaviours beyond the given time bound are not
considered, and may thus violate a safety requirement. Representatives are STRONG [15], HySon [7],

24

name
run time
car acceleration
car speed
car position

this article
100 iterations
unbounded
166 ms
166 ms
[-0.820 1.31]
[-1.262 1.31]
[-1.013 5.11]
[-4.515 6.15]
[43.7 83.4]
[40.86 91.9]

LGG
100 iterations 200 iterations 300 iterations
50 ms
140 ms
195 ms
[-0.815 1.31] [-0.968 1.31] [-0.968 1.31]
[-1.013 4.97] [-3.651 4.97] [-3.677 4.97]
[44.5 83.4]
[44.5 88.87] [44.5 88.87]

TABLE IV
Comparison on convoyCar2 benchmark, between this work and the LGG algorithm [28]

CORA [1], HYLAA [3] and SpaceEx [20].

B. Unbounded Reachability Analysis
The second approach, epitomised in static anal-

Set-based simulation methods generalise guaran- ysis methods [30], explores unbounded-time horiteed integration [6], [37]from enclosing intervals to zons. It employs conservative over-approximations to
relational domains. They use precise abstractions
with low computational cost to over-approximate
sets of reachable states up to a given time horizon.

achieve completeness and decidability over infinite
time horizons.
Unbounded techniques attempt to infer a loop

Early tools used polyhedral sets (HyTech [31] and
PHAVer [19]), polyhedral flow-pipes [10], ellip-

invariant, i.e., an inductive set of states that includes
all reachable states. If the computed invariant is

soids [5] and zonotopes [24]. A breakthrough was

disjoint from the set of bad states, this proves that the

been achieved by [25], [28], with the representation of convex sets using template polyhedra and

latter are unreachable and hence that the loop is safe.
However, analysers frequently struggle to obtain an

support functions. This method is implemented in
the tool SpaceEx [20], which can handle dynamical

invariant that is precise enough with acceptable computational cost. The problem is evidently exacerbated

systems with hundreds of variables. Although it
may use exact arithmetic to maintain soundness, it

by non-determinism in the loop, which corresponds
to the case of open systems. Prominent representa-

performs computations using floating-point numbers:
this is a deliberate choice to boost performance,

tives of this analysis approach include Passel [35],
Sting [12], and abstract interpreters such as Astrée

which, although quite reasonable, its implementation
is numerically unsound and therefore does not pro-

[4] and InterProc [33]. Early work in this area has
used implementations of abstract interpretation and

vide genuine formal guarantees. In fact, most tools
using eigendecomposition over a large number of

widening [13], which are still the foundations of
most modern tools. The work in [30] uses abstract

variables (more than 10) are numerically unsound
due to the use of unchecked floating-point arithmetic.

interpretation with convex polyhedra over piecewiseconstant differential inclusions. Dang and Gawl-

Another breakthrough in performance was done by

itza [14] employ optimisation-based (max-strategy

HYLAA [3] which was the first tool to solve all
high order problems of hundreds and thousands

iteration) with linear templates for hybrid systems
with linear dynamics. Relational abstractions [38]

dimensions. Other approaches use specialised con- use ad-hoc “loop summarisation” of flow relations,
straint solvers (HySAT [18], iSAT [16]), or SMT while abstract acceleration focuses on linear relations
encodings [11], [29] for bounded model checking
of hybrid automata.

analysis [26], [27], which is common in program
analysis.

25

C. Abstract Acceleration

References

Abstract acceleration [26], [27], [34] captures the
effect of an arbitrary number of loop iterations with a
single, non-iterative transfer function that is applied
to the entry state of the loop (i.e., to the set of initial
conditions of the linear dynamics). Abstract acceleration has been extended from its original version
to encompass inputs over reactive systems [40] but
restricted to subclasses of linear loops, and later to
general linear loops but without inputs [34].
The work presented in this article lifts these limitations by presenting abstract acceleration for general
linear loops with inputs [8], developing numeric

[1] M. Althoff. An introduction to cora 2015. In ARCH@
CPSWeek, pages 120–151, 2015.
[2] E. Asarin, T. Dang, and A. Girard. Hybridization methods
for the analysis of nonlinear systems. Acta Informatica,
43(7):451–476, 2007.
[3] S. Bak and P. S. Duggirala. Hylaa: A tool for computing
simulation-equivalent reachability for linear systems. In
Proceedings of the 20th International Conference on Hybrid
Systems: Computation and Control, HSCC 2017, Pittsburgh,
PA, USA, April 18-20, 2017, pages 173–178, 2017.
[4] B. Blanchet, P. Cousot, R. Cousot, J. Feret, L. Mauborgne,
A. Miné, D. Monniaux, and X. Rival. A static analyzer
for large safety-critical software. In PLDI, pages 196–207.
ACM, 2003.

techniques for scalability and extending the domain
to continuous time systems.

[5] O. Botchkarev and S. Tripakis. Verification of hybrid
systems with linear differential inclusions using ellipsoidal
approximations. In HSCC, LNCS, pages 73–88. Springer,
2000.

IX. Conclusions and Future Work

[6] O. Bouissou. Analyse statique par interprétation abstraite
de systèmes hybrides. PhD thesis, École Polytechnique,
2008.

We have presented an extension of the Abstract
Acceleration paradigm to guarded LTI systems (linear loops) with inputs, overcoming the limitations of
existing work dealing with closed systems. We have
decisively shown the new approach to over-compete
state-of-the-art tools for unbounded-time reachability
analysis in both precision and scalability. The new
approach is capable of handling general unboundedtime safety analysis for large scale open systems
with reasonable precision and fast computation times.
Conditionals inside loops and nested loops are out
of the scope of this paper.
Work to be done is extending the approach to nonlinear dynamics, which we believe can be explored
via hybridisation techniques [2], and to formalise the
framework for general hybrid models with multiple
guards and location-dependent dynamics, with the
aim to accelerate transitions across guards rather than
integrate individual accelerations on either side of the
guards.
a) Acknowledgments.: We would like to thank
Colas Le Guernic for his constructive suggestions
and comments on the paper.

[7] O. Bouissou, S. Mimram, and A. Chapoutot. Hyson: Setbased simulation of hybrid systems. In Rapid System Prototyping (RSP), 2012 23rd IEEE International Symposium on,
pages 79–85. IEEE, 2012.
[8] D. Cattaruzza, A. Abate, P. Schrammel, and D. Kroening.
Unbounded-time analysis of guarded LTI systems with inputs by abstract acceleration. In SAS, volume 9291 of LNCS,
pages 312–331. Springer, 2015.
[9] D. Cattaruzza, A. Abate, P. Schrammel, and D. Kroening.
Sound numerical computations in abstract acceleration. In
International Workshop on Numerical Software Verification,
pages 38–60. Springer, 2017.
[10] A. Chutinan and B. H. Krogh. Computing polyhedral
approximations to flow pipes for dynamic systems. In CDC,
pages 2089–2094. IEEE Computer Society, 1998.
[11] A. Cimatti, S. Mover, and S. Tonetta. SMT-based verification of hybrid systems. In AAAI Conference on Artificial
Intelligence. AAAI Press, 2012.
[12] M. A. Colón, S. Sankaranarayanan, and H. B. Sipma. Linear
invariant generation using non-linear constraint solving. In
CAV, pages 420–432. Springer, 2003.
[13] P. Cousot and R. Cousot. Abstract interpretation: A unified
lattice model for static analysis of programs by construction
or approximation of fixpoints. In POPL, pages 238–252,
1977.
[14] T. Dang and T. M. Gawlitza. Template-based unbounded
time verification of affine hybrid automata. In APLAS,
LNCS, pages 34–49. Springer, 2011.
[15] Y. Deng, A. Rajhans, and A. A. Julius. STRONG: A
trajectory-based verification toolbox for hybrid systems. In

26

Quantitative Evaluation of Systems, volume 8054 of LNCS,
pages 165–168. Springer, 2013.
[16] A. Eggers, M. Fränzle, and C. Herde. SAT Modulo ODE:
A direct SAT approach to hybrid systems. In ATVA, volume
5311 of LNCS, pages 171–185. Springer, 2008.
[17] A. Fehnker and F. Ivancic. Benchmarks for hybrid systems
verification. In HSCC, pages 326–341. Springer, 2004.
[18] M. Fränzle and C. Herde. HySAT: An efficient proof engine
for bounded model checking of hybrid systems. Formal
Methods in System Design, 30(3):179–198, 2007.
[19] G. Frehse. PHAVer: Algorithmic verification of hybrid
systems past HyTech. In HSCC, volume 3414 of LNCS,
pages 258–273. Springer, 2005.
[20] G. Frehse, C. L. Guernic, A. Donzé, R.
R. Ripado, A. Girard, T. Dang, and O.
Scalable verification of hybrid systems.
6806 of LNCS, pages 379–395. Springer,

Ray, O. Lebeltel,
Maler. SpaceEx:
In CAV, volume
2011.

[21] K. Fukuda and A. Prodon. Double description method
revisited. In Combinatorics and computer science, pages
91–111. Springer, 1996.
[22] S. Gao, J. Avigad, and E. M. Clarke. δ-complete decision
procedures for satisfiability over the reals. In Automated
Reasoning, pages 286–300. Springer, 2012.
[23] P. K. Ghosh and K. V. Kumar. Support function representation of convex bodies, its application in geometric
computing, and some related representations. Computer
Vision and Image Understanding, 72:379–403, 1998.
[24] A. Girard. Reachability of uncertain linear systems using
zonotopes. In HSCC, volume 3414 of LNCS, pages 291–
305. Springer, 2005.
[25] A. Girard, C. L. Guernic, and O. Maler. Efficient computation of reachable sets of linear time-invariant systems with
inputs. In HSCC, volume 3927 of LNCS, pages 257–271.
Springer, 2006.
[26] L. Gonnord and N. Halbwachs. Combining widening and
acceleration in linear relation analysis. In SAS, LNCS, pages
144–160. Springer, 2006.
[27] L. Gonnord and P. Schrammel. Abstract acceleration in
linear relation analysis. Science of Computer Programming,
93(Part B):125–153, 2014.
[28] C. L. Guernic and A. Girard. Reachability analysis of hybrid
systems using support functions. In CAV, volume 5643 of
LNCS, pages 540–554. Springer, 2009.
[29] S. Gulwani and A. Tiwari. Constraint-based approach for
analysis of hybrid systems. In CAV, volume 5123 of LNCS,
pages 190–203. Springer, 2008.
[30] N. Halbwachs, P. Raymond, and Y.-E. Proy. Verification of
linear hybrid systems by means of convex approximations.
In SAS, volume 864 of LNCS, pages 223–237. Springer,
1994.
[31] T. A. Henzinger, P.-H. Ho, and H. Wong-Toi. HyTech: A
model checker for hybrid systems. Journal on Software
Tools for Technology Transfer, 1(1-2):110–122, 1997.

[32] R. A. Horn and C. R. Johnson. Matrix analysis. Cambridge
university press, 2012.
[33] B. Jeannet.
Interproc analyzer for recursive
programs
with
numerical
variables,
2010.
http://pop-art.inrialpes.fr/interproc/interprocweb.cgi.
[34] B. Jeannet, P. Schrammel, and S. Sankaranarayanan. Abstract acceleration of general linear loops. In POPL, pages
529–540. ACM, 2014.
[35] T. T. Johnson and S. Mitra.
Passel: A verification
tool for parameterized networks of hybrid automata, 2012.
https://publish.illinois.edu/passel-tool/.
[36] P. Lancaster and M. Tismenetsky. The Theory of Matrices.
Academic Press, 2nd edition, 1984.
[37] R. Löhner. Einschließung der Lösung gewöhnlicher Anfangsund Randwertaufgaben und Anwendungen. PhD thesis,
Universität Karlsruhe, 1988.
[38] S. Sankaranarayanan and A. Tiwari. Relational abstractions
for continuous and hybrid systems. In CAV, volume 6806
of LNCS, pages 686–702. Springer, 2011.
[39] P. Schrammel. Unbounded-time reachability analysis of
hybrid systems by abstract acceleration. In Embedded
Software, pages 51–54. IEEE, 2015.
[40] P. Schrammel and B. Jeannet. Applying abstract acceleration
to (co-)Reachability analysis of reactive programs. Journal
of Symbolic Computation, 47(12):1512–1532, 2012.

