Abstract
Reachability analysis of dynamical models is a relevant problem that has seen much progress in the last decades, however with clear limitations pertaining to the nature of the dynamics and the soundness of the results. This article focuses on sound safety verification of unbounded-time (infinite-horizon) linear time-invariant (LTI) models with inputs using reachability analysis. We achieve this using counterexample-guided Abstract Acceleration: this approach over-approximates the reachability tube of the LTI model over an unbounded time horizon by using abstraction, possibly finding concrete counterexamples for refinement based on the given safety specification. The technique is applied to a number of LTI models and the results show robust performance when compared to state-of-the-art tools.
Similar content being viewed by others
1 Introduction
Linear loops are a ubiquitous programming pattern [47]. Linear loops iterate over continuous variables (in the case of physical systems) or discrete variables (in the case of loops in digital programs), which are updated using a linear transformation. Linear loops may be guarded, i.e., halt if a given linear condition holds: at this point the system may either elicit a new behaviour, or simply terminate. 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 many hybrid dynamical models [6, 47]. 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 control systems can be modelled as safety properties, i.e., deciding reachability of certain bad states, for which the model exhibits unsafe behaviour. Sets of bad states may, in linear loops, be described by means of guard assertions, namely (linear) constraints over their continuous variables. Reachability in linear programs, however, is a formidable challenge for automatic analysers: despite the restriction to linear transformations (i.e., linear dynamics) and linear guards, it is undecidable in the general case. The problem has been related to the Skolem Problem, which includes a subset of cases that are known to be decidable (e.g., the orbit problem [46], and for low-order models [54]). The problem is decidable when restricted to finite state spaces, but even then, algorithms are costly (for the infinite-time case in particular).
The goal of this article is to push the frontiers of unbounded-time reachability analysis: we aim to devise a method that performs sound reasoning about unbounded trajectories by abstract acceleration. Abstract acceleration [37, 38, 44] approximates the effect of an arbitrary number of loop iterations (up to infinity) 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). This article extends the work in [44] to models with non-deterministic inputs, elaborating early work in [57] and completing [14, 16].
The key original contributions of this article are as follows:
-
1.
We present a new technique to include time-varying non-deterministic inputs in the abstract acceleration of general linear loops.
-
2.
We extend abstract acceleration to the continuous time case.
-
3.
We introduce a technique to help the analysis of support functions in complex spaces, in order to increase the precision of previous abstract acceleration methods.
-
4.
We develop a counterexample-guided refinement for Abstract Acceleration for safety verification, maximising speed when high precision is not necessary, thus allowing for optimal analysis within a safe region.
-
5.
We provide an implementation of the discussed procedures as a tool, called Axelerator, which is available at http://www.cprover.org/LTI/ We test these novel procedures using a broad set of experiments: our benchmarks (including specifications of the initial states, input ranges and guard sets), are also available.
Finally, in Sect. 10 we provide a thorough review and comparison with related work.
2 Preliminaries
2.1 Linear Loops with Inputs—Syntax
A discrete time LTI model may be described as a simple linear loop. Simple linear loops are functions expressed in the form:
where \(\varvec{x} \in \mathbb {R}^p\) is a valuation on the state variables, \(\varvec{G}\varvec{x} \le \varvec{h}\) is a linear constraint on the states (with \(\varvec{G} \in \mathbb {R}^{r \times p}\) and \(\varvec{h} \in \mathbb {R}^r\)), \(\varvec{u} \in \mathbb {R}^q\) is a non-deterministic input, and \(\varvec{A} \in \mathbb {R}^{p \times p}\) and \(\varvec{B} \in \mathbb {R}^{p \times q}\) are linear transformations characterising the dynamics of the model. 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 =\{\varvec{x} : \varvec{G}\varvec{x}\le \varvec{h}\}\). The purpose of the guard is to restrict the dynamics to a given set, either to ensure safety (think for example of speed limits) and/or to change the behaviour of the model under certain conditions (e.g., once we have reached a certain state we begin a new process).
In particular, the special case where \(G = \top \) (i.e., “while true”) represents a time-unbounded loop with no guards, for which the discovery of a suitable invariant is paramount.
2.2 Model Semantics
The traces of the model starting from an initial set \(X_0\subseteq \mathbb {R}^p\), with inputs restricted to the set \(U \subseteq \mathbb {R}^q\), are sequences \( \varvec{x}_0 \xrightarrow {\varvec{u}_0} \varvec{x}_1 \xrightarrow {\varvec{u}_1} \varvec{x}_2 \xrightarrow {\varvec{u}_2} \ldots \), where \( \varvec{x}_0 \in X_0\) and \(\forall k\ge 0, \varvec{x}_{k+1} = \tau (\varvec{x}_k,\varvec{u}_k) \) and \(\varvec{u}_k \in U\), satisfying:
We extend the point-wise notation above to convex sets of states and inputs (\(X_k\) and U), and denote the set of states reached from set \(X_k\) by \(\tau \) in one step as:
We furthermore denote the set of states reached from \(X_0\) via \(\tau \) in n steps (n-reach set, constrained by G), for \(n\ge 0\):
Since the sets \(X_0\) and U are convex, the transformations \(\varvec{A}\) and \(\varvec{B}\) are linear, and vector sums preserve convexity, the sets \(X_n = \tau ^n(X_0,U)\) are also convex.
We define the n-reach tube
as the union of k-reach sets over n iterations. Moreover, \({\hat{X}} =\bigcup _{n\ge 0} \tau ^n(X_0,U)\) extends the previous notion to an unbounded time horizon (transitive closure).
2.3 Spectral Eigendecomposition
Eigendecomposition [52] denotes the factorisation of a matrix into a canonical form that is characterised by having the non-zero elements (the eigenvalues) only on the main diagonal (note that not all matrices can be factorised this way, as discussed later). Let \(\varvec{A} \in \mathbb {R}^p\) be a diagonalizable square matrix. The eigendecomposition yields the equation \(\varvec{A}=\varvec{S}\varvec{\varLambda }\varvec{S}^{-1}\), where \(\lambda _i=\varvec{\varLambda }_{ii}\) are the eigenvalues of \(\varvec{A}\) and \(\varvec{S}_{*i}\) their corresponding eigenvectors, sharing the known property \(\varvec{A}\varvec{S}_{*i}=\lambda _i\varvec{S}_{*i}\). When a square matrix is not diagonalizable because of repeated eigenvalues, it can be factored into what is known as the Jordan Form, which in addition to the eigenvalues of the matrix, may contain unitary values in the immediate upper diagonal in the case of duplicate eigenvalues: \(\varvec{A}=\varvec{S}\varvec{J}\varvec{S}^{-1}\), where:
In the case of the Jordan Form, the eigenvectors corresponding to repeated eigenvalues are called generalised eigenvectors and have the property that \((\varvec{A}-\lambda _s\varvec{I})^j\varvec{v}_j = 0\), where \(\varvec{v}_j\) is the \(j{\text {th}}\) generalised eigenvector related to eigenvalue \(\lambda _s\).
Finally, the eigenvalues of a real matrix may be complex numbers. This is inconvenient in the subsequent analysis, so rather than using complex arithmetic on these numbers, we choose a different representation over the pseudo-eigenspace. Pseudo-eigendecomposition relies on the observation that complex eigenvalues of a real matrix always come in conjugate pairs. Relaxing the restriction of non-zero values on the main diagonal to include the immediate off-diagonal terms, we leverage the following equivalence:
where \(re(\varvec{v})\) and \(im(\varvec{v})\) are the real and imaginary part of \(\varvec{v}\), respectively. In the case of a non-diagonal Jordan form, the columns are rearranged first (including the upper diagonal ones), and the conversion above is then performed. This representation is also called the Real Jordan Form.
2.4 Support Functions
2.4.1 Definition of Support Functions
A support function is a convex function defined over the vector space \(\mathbb {R}^p\), which describes the distance of a supporting hyperplane to the origin from a given set in \(\mathbb {R}^p\), as illustrated in Fig. 1.
Support function for a polyhedral set in \(\mathbb {R}^2\). The distance between the tangent line and the origin is \(\rho _X\left( \frac{\varvec{v}}{\Vert \varvec{v} \Vert }\right) \)
Support function for a circular set using six directions (\(\varvec{v}_1 \ldots \varvec{v}_6\)). The resulting polyhedron is an over-approximation of the original set (note that directions need not be symmetrical)
Support functions can 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, a support function characterises the distance from the origin to the hyperplane that is orthogonal to the given direction and that touches its convex hull at its farthest. For example, the support function of a sphere centred at the origin given any unit vector \(\varvec{v}\) in \(\mathbb {R}^3\), evaluates to the radius of the sphere.
The intersection of multiple half spaces, each obtained by sampling a support function in a specific direction, can generate a polyhedron (see Fig. 2), as discussed further in the next section. Finitely sampled support functions (i.e., using a limited number of directions) are template polyhedra in which the directions are not fixed, which helps to avoid wrapping effects (wherein sampling in given directions creates an over-approximation of a set that is not aligned with said directions). The larger the number of distinct directions provided, the more precisely represented the set is. In more detail, given a direction \(\varvec{v} \in \mathbb {R}^p\), the support function of a non-empty set \(X \subseteq \mathbb {R}^p\) in the direction of \(\varvec{v}\) is defined as
where \(\varvec{x} \cdot \varvec{v}=\sum _{i=0}^p \varvec{x}_i \varvec{v}_i\) is the dot product of the two vectors. Support functions apply to any non-empty set \(X \subseteq \mathbb {R}^p\), but they are most useful when representing convex sets. We will restrict ourselves to the use of convex polyhedra, in which case the definition of the support function translates to solving the following linear program:
2.4.2 Properties of Support Functions
Several properties of support functions allow us to reduce the complexity of operations. The most significant ones are [34]:
where \(\varvec{v},\varvec{v}_1,\varvec{v}_2 \in \mathbb {R}^p\). As can be seen by their structure, some of these properties reduce the complexity to lower-order polynomial or even to constant time, by turning matrix-matrix multiplications (\({\mathcal {O}}(p^3)\)) into matrix-vector (\({\mathcal {O}}(p^2)\)), or into scalar (\({\mathcal {O}}(p)\)) multiplications.
2.4.3 Support Functions in Complex Spaces
The literature does not state, to the best of our knowledge, any use of support functions in complex spaces. Since we are applying their concept to eigenspaces, which may have complex conjugate eigenvalues, we extend the definition of support functions to encompass the corresponding operations on complex spaces, which are given explicitly.
Theorem 1
A support function in a complex vector field is a transformation:
where \(re(\cdot )\) defines the real part of a complex number. The dot product used here is commonly defined in a complex space as:
where the element \(b_i^*\) is the complex conjugate of \(b_i\).
Proof
Let \(f : \mathbb {C}^p \rightarrow \mathbb {R}^{2p}\) and
Using abstract interpretation [23] we define a Galois connection \(\alpha (\varvec{x})= f(\varvec{x})\) and \(\gamma (\varvec{x}')= f^{-1}(\varvec{x}')\), which is clearly a one-to-one relation. We can therefore establish an equivalence between \(\rho _X(\varvec{v}) = \rho _{X'}(\varvec{v}')\). \(\square \)
As we shall see later, the imaginary part of the dot product is not relevant to the support function, and we therefore disregard it. Using properties of support functions, we now have:
which is consistent with the real case when \(\theta =0\). The reason why \(e^{i\theta }\) cannot be factored out as a constant is because it executes a rotation on the vector, and therefore follows the same rules as a matrix multiplication, namely:
Notice the resemblance of this matrix to a pseudo-eigenvalue representation (see Eq. (6)). Since the vectors we are interested in are conjugate pairs (because they have been created by a transformation into a matrix eigenspace), we can transform our problem into a pseudo-eigenvalue representation. Since this removes the imaginary component from the setup, we can evaluate the support function using standard methods and properties, by transforming the conjugate pairs into separate vectors representing the real and imaginary parts and rotation matrices as in the equation above.
Recently, Adimoolam and Dang [1] have presented a calculus for complex zonotopes that exploits the same idea, namely using complex eigenvalues to represent ellipsoidal and linear boundaries for reach sets. Their projections into real space correspond to our semi-spherical approximations for variable inputs. Apart from the basic difference in domain (zonotopes vs. polyhedra), which in itself changes some aspects of the problem, the main difference is that the authors perform their analysis over the complex space, whereas we ultimately apply the method in the real space by using pseudo-eigenspaces.
2.5 Convex Polyhedra
A polyhedron is a subset of \(\mathbb {R}^p\) with planar faces. Each face is supported by a hyperplane that creates a half-space, and the intersections of these hyperplanes are the edges (and vertices) of the polyhedron. A polyhedron is said to be convex if a line segment joining any two points of its surface is contained within its interior. Convex polyhedra are better suited than general polyhedra to define an abstract domain, mainly because they have a simpler representation and because operations over convex polyhedra are in general easier than over general polyhedra. There are a number of properties of convex polyhedra that make them ideal for abstract interpretation over 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 [32], as detailed next. The algorithm to switch between these two representations is given in Sect. 2.6.5.
2.5.1 Vertex Representation
Since every edge in the polyhedron corresponds to a line between two vertices, and every face corresponds to the area enclosed by a set of co-planar edges, a full description of the polyhedron is obtained simply by 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 by the transformed vertices). Formally, a polyhedron can be described as a set \(V \in \mathbb {R}^p\) such that \(\varvec{v} \in V\) is a vertex of the polyhedron.
2.5.2 Inequality Representation (a.k.a. Face Representation)
The dual of the Vertex representation is the inequality representation, where each inequality represents a face of the polyhedron. Each face corresponds to a bounding hyperplane of the polyhedron (with the edges being the intersection of two hyperplanes and the vertices being the intersection of p or more hyperplanes), and is described mathematically as a function of the vector that is normal to the hyperplane. This representation can be made minimal by eliminating redundant inequalities that do not correspond to any face of the polyhedron. If 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 topological region in \(\mathbb {R}^p\) described by the set
where the rows \(\varvec{C}_{i,*}\) for \(i \in [1,\ldots ,m]\) correspond to the transposed vectors normal to the faces of the polyhedron, and \(\varvec{d}_i\) for \(i \in [1,\ldots ,m]\) to the value of the support function of X in the corresponding direction. For simplicity in the presentation, we will extend the use of the support function operator as follows:
2.6 Operations on Convex Polyhedra
There are a number of operations that we need to be able to perform on convex polyhedra.
2.6.1 Translations
Given a vertex representation V and a translation vector \(\varvec{t}\), the transformed polyhedron is:
Given an inequality representation X and a translation vector \(\varvec{t}\), the transformed polyhedron is:
2.6.2 Linear Transformations
Given a vertex representation V and a linear transformation \(\varvec{L}\), the transformed polyhedron is:
Given an inequality representation X and a linear transformation \(\varvec{L}\), the transformed polyhedron corresponds to
where \(\varvec{L}^+\) represents the pseudo-inverse of \(\varvec{L}\) [55]. In the case when the inverse \(\varvec{L}^{-1}\) exists, then:
From this we can conclude that linear transformations are more efficient when using vertex representation, except when the inverse of the transformation exists and is known a priori. This work makes use of this assumption to avoid alternating between representations.
2.6.3 Set Sums
The addition of two polyhedra is defined such that the resulting set contains the sum of all pairs of points inside the original polyhedra. This operation is commonly known as the Minkowski sum, namely:
Given two vertex representations \(V_1\) and \(V_2\), the resulting polyhedron is
where \(conv(\cdot )\) is the convex hull of the set of vertices contained in the Minkowski sum. Let
be two sets, then
where
Because these sets correspond to systems of inequalities, they can be reduced by removing redundant constraints. Note that if \(\varvec{C}_1=\varvec{C}_2\), then
2.6.4 Set Hadamard Products
Definition 1
Given two vertex representations \(V'\) and \(V''\), we define the set Hadamard product operation using such representations as
where \(\circ \) represents the Hadamard (coefficient-wise) product of the vectors.
Lemma 1
The set \(V=V' \circ V''\) introduced in the previous definition is a convex set that contains all possible combinations of products between elements of sets \(V'\) and \(V''\).
Proof
Given a convex set X with a vertex representation V, by definition we have
which extends to multiple points [13] as
Applying the Hadamard product, we obtain
where \(x_1 \circ x_2 = \sum _{i=1}^{|V'|} k_i\varvec{v}_i' \circ \sum _{j=1}^{|V''|} k_j\varvec{v}_j'' \text { with } \varvec{v}_i' \in V', \varvec{v}_j'' \in V''\). Simplifying, we obtain
where \(\varvec{v}_{ij}=\varvec{v}_i' \circ \varvec{v}_j'' \in V \text { and } \sum _{ij=1}^{|V'||V''|}k_{ij}=\sum _{i=1}^{|V'|}\sum _{j=1}^{|V''|} k_ik_j = 1\). \(\square \)
Note that in the case of the inequality representation, there is no direct result for this product. We therefore enumerate the sets in one of the polyhedra, and use linear solving algorithms to find an over-approximation:
where \(\varvec{t}\) is a template direction for a face in the over-approximation, T is the set of directions selected for the over-approximation, and \(V'\) is the set of vertices of X.
2.6.5 Vertex Enumeration
The vertex enumeration algorithm obtains a list of all vertices of a polyhedron, given a face representation of its bounding hyperplanes. Given the duality of the problem, it is also possible to find the bounding hyperplanes given a vertex description if the chosen algorithm exploits this duality. In this case the description of V is given in the form of a matrix inequality \(\varvec{V}\varvec{x} \le [\begin{array}{cccc}1&1&\cdots&1\end{array}]^\intercal \) with \(\varvec{V}=[\begin{array}{ccc}\varvec{v}_1&\cdots&\varvec{v}_m\end{array}]^\intercal , \varvec{v}_i \in V\). Similarly, \(\varvec{C}\) can be described as a set containing each of its rows. There are two algorithms that efficiently solve the vertex enumeration problem. lrs [4] is a reverse search algorithm, while cdd [32] 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
be the polyhedral cone represented by \(\varvec{C}\). The pair \((\varvec{C},V)\) is said to be a double description pair if
where 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 rays of X if the cone is pointed (i.e., it has a vertex at the origin). This last criterion can be ensured by translating a polyhedral description so that it includes the origin, and then translating the vertices back once they have been discovered (see Sect. 2.6).
We also note that:
The vertex enumeration algorithm starts by finding a base \({\mathcal {C}}_K\) which contains a number of vertices of the polyhedron. This can be done by pivoting over a number of different rows in \(\varvec{C}\) and selecting the feasible visited points, which are known to be vertices of the polyhedron (pivoting p times will ensure at least one vertex is visited if the polyhedron is non-empty). \({\mathcal {C}}_K\) is represented by \(\varvec{C}_K\) which contains the rows used for the pivots. The base \({\mathcal {C}}_K\) is then iteratively expanded to \({\mathcal {C}}_{K+i}\) by exploring the \(i{\text {th}}\) row of \(\varvec{C}\) until \({\mathcal {C}}_K= {\mathcal {C}}\). The corresponding pairs \((\varvec{C}_{K+i},V_{K+i})\) are constructed using the information from \((\varvec{C}_K,V_K)\) as follows.
Let \(\varvec{C}_K \in \mathbb {R}^{n_K \times p}\), \(\varvec{C}_{i,*} \in \mathbb {R}^{1 \times p}\), \(V_K \in \mathbb {R}^p\),
be the spaces outside, inside and on the \(i{\text {th}}\) hyperplane and
the existing vertices lying on each of these spaces. Then [32],
3 Abstract Acceleration—Overview of the Algorithm
Abstract acceleration is a method that seeks to precisely describe the dynamics of a transition system over a number of steps using a concise description between the first and final steps. More precisely, it looks for a direct formula to express the post-image of an unfolded loop from its initial states. Formally, given the dynamics in Eq. (1), an acceleration formula aims to compute the reach tube based on (3) using a function f such that \(f(\cdot )=\tau ^n(\cdot )\). In the case of models without inputs, this equation can be derived from the expression \(\varvec{x}_n=\varvec{A}^n\varvec{x}_0\).
3.1 Overview of the Algorithm
The basic steps required to compute a reach tube using abstract acceleration are given in Fig. 3.
-
1.
The process starts by performing eigendecomposition of the dynamics (based on matrix \(\varvec{A}\)) in order to transform the problem into a simpler one. Since we use unsound arithmetic at this stage, the results are quantities that are marked using a tilde (as in \(\tilde{\varvec{S}}, \tilde{\varvec{J}}\)).
-
2.
The second step involves upper-bounding the rounding errors in order to obtain sound results: bounds on eigenvalues, for example, are well known from the literature and can be obtained as \(|\lambda - {\tilde{\lambda }}| < \Vert \tilde{\varvec{S}}\tilde{\varvec{J}}\tilde{\varvec{S}}^{-1}-\varvec{A} \Vert _2\). In general, a variety of off-the-shelf tools may be used, but since larger problems require numerical algorithms for scalability, all subsequent steps are performed using interval arithmetic in order to maintain soundness: we identify corresponding interval-based quantities with bold symbols (e.g., \({\mathbb {S}}, {\mathbb {J}}\)), as well as subsequent matrix operations (e.g., computing the inverse of \({\mathbb {S}}\)). Thus we obtain \(\Vert {\mathbb {S}}{\mathbb {J}}{\mathbb {S}}^{-1}-\mathbb {A} \Vert _2\) by extending the original unsound matrices by one least-significant bit element-wise. We also note that this equation is symmetric, which is why we use the estimated eigenvectors to calculate the error on the original eigenvalues (see [16] for further details). While bounds on the eigenvectors can be calculated from the eigenvalues, we choose a more complex yet faster over-approximation, which is described in [16].
-
3.
The inverse of the generalised eigenvectors (\(\mathbb {S}^{-1}\)) is calculated soundly (using interval arithmetic).
-
4.
The problem is transformed into canonical form by multiplying both sides of the equation by \(\mathbb {S}^{-1}\) (we use blackboard symbols to indicate interval vectors and matrices, which are employed to ensure sound operations), obtaining
$$\begin{aligned} X_k' {=}\mathbb {J}\left( X_{k-1}' \cap G' \right) +U', \quad \text { where } \quad X_k' {=}\mathbb {S}^{-1}X_k, \quad U'{=}\mathbb {S}^{-1}\varvec{B}U,\quad G'{=}\{\varvec{x} \mid \varvec{G}\mathbb {S}\varvec{x} \le \varvec{h} \}. \end{aligned}$$ -
5.
We calculate the number of iterations n based on the guard set, as explained in Sect. 7. If there are no guards, we set \(n=\infty \). This number need not be exact: if we over-approximate the number of iterations, the resulting reach tube will further over-approximate the desired one.
-
6.
We over-approximate the dynamics subject to general inputs (for parametric inputs or in the absence of inputs this step will be ignored), using the techniques described in Sect. 6.2.
-
7.
We calculate the accelerated dynamics using the techniques described in Sect. 5.1.
-
8.
We transform the input and initial sets into vertex representation, to be used as the source for the reach tube calculation.
-
9.
We employ a sound simplex algorithm [16] to evaluate the convex-set Hadamard product of the abstract dynamics and of the initial set. The two most important elements of the sound simplex are that it uses interval arithmetic to pivot, and that at every stage it verifies the intersection between vertices in order to avoid pivoting on unsound solutions. The latter step is better explained by considering that the simplex algorithm operates by visiting adjacent vertices. The interval arithmetic ensure that the solution at the last visited vertex is sound, but if there is an intersection, the new pivot may choose to start on the intersecting vertex instead (which is not sound), thus, by checking the intersection and extending the interval to encompass both vertices, we retain soundness (see [16] for details).
-
10.
Since we have operated in the eigenspace so far, we transform the reach tube back into the state space via multiplication by \(\mathbb {S}\).
A block diagram describing the different steps used to calculate the abstract reach tube of a model via Abstract Acceleration. The white box is the numerical eigensolver stage. Blue boxes are soundness restoration stages. Red boxes represent linear transformations of the problem. The orange boxes denote abstractions defined in this paper, and the green boxes the reachability computation in the abstract domain. (Color figure online)
4 Abstract Matrices in Abstract Acceleration
We introduce the concept of an abstract matrix.
Definition 2
An abstract matrix \({\mathcal {A}}^n \subseteq \mathbb {R}^{p \times p}\) is an over-approximation of the union of the powers of the matrix \(\varvec{A}^k\) such that \({\mathcal {A}}^n \supseteq \left\{ \varvec{A}^k : k\in [0, \ldots ,n] \right\} \). Its application to the initial set \(X_0\) results in
such that \({\hat{X}}_n^\sharp \supseteq {\hat{X}}_n\) is an over-approximation of the reach tube described in Eq. (4).
Next we explain how to compute such abstract matrices. For simplicity, we first describe this computation for matrices \(\varvec{A}\) with real eigenvalues, whereas the extension to the complex case will be addressed in Sect. 4.1. Similar to [44], we first have to compute the Jordan normal form of \(\varvec{A}\). Let \(\varvec{A}=\varvec{S}\varvec{J}\varvec{S}^{-1}\) where \(\varvec{J}\) is the normal Jordan form of \(\varvec{A}\), and \(\varvec{S}\) is made up by the corresponding generalised eigenvectors. We can then easily compute \(\varvec{A}^n=\varvec{S}\varvec{J}^n\varvec{S}^{-1}\), where given a set of r eigenvalues \(\lambda _s\) with geometric multiplicities \(p_s\) and \(s \in [1,\ldots ,r]\), we have
The abstract matrix \({\mathcal {A}}^n\) is computed as an abstraction over a set of vectors \(\varvec{m}^k \in \mathbb {R}^p, k \in [1,\ldots ,n]\) of distinct entries of \(\varvec{J}^k\), as explained below.
Let \(\varvec{I}_s=[\begin{array}{cccc}1&0&\cdots&0\end{array}] \in \mathbb {R}^{p_s}\). The vector \(\varvec{m}^k\) is obtained by the transformation \(\varphi ^{-1}\) (which is always invertible) as
such that \(\varvec{J}^k = \varphi (\varvec{m}^k)\).
If \(\varvec{J}\) is diagonal [44], then \(\varvec{m}^k\) results in the vector made up of powers of the eigenvalues, \([\begin{array}{ccc}\lambda _1^k&\cdots&\lambda _p^k\end{array}]\). The diagonal entries in the abstract matrix is thus bound by the intervals
We observe that the spectrum of the abstract matrix \(\sigma ({\mathcal {A}}^n)\), which can be derived from its entries in Eq. (12), over-approximates \(\bigcup _{k\in [1,\ldots ,n]} \sigma (\varvec{A}^k)\).
In the case of the \(s\mathrm {th}\) Jordan block \(\varvec{J}_s\) with geometric multiplicity \(p_s > 1\), observe that the first row of \(\varvec{J}_s^n\) contains all (possibly) distinct entries of \(\varvec{J}_s^n\). Hence, the vector section \(\varvec{m}_s\) is the concatenation of the (transposed) first row vectors \(\left( \lambda _s^n , \left( {\begin{array}{c}n\\ 1\end{array}}\right) \lambda _s^{n-1}, \ldots , \left( {\begin{array}{c}n\\ p_s-1\end{array}}\right) \lambda _s^{n-p_s+1}\right) ^\intercal \) of \(\varvec{J}_s^n\).
Since \(\varphi \) transforms the vector \(\varvec{m}\) into the shape of (10) of \(\varvec{J}^n\), it is called a matrix shape [44]. We then define the abstract matrix as
where the constraint \(\varvec{\varPhi }\varvec{m} \le \varvec{f}\) is synthesised from intervals associated to the individual eigenvalues and to their combinations. More precisely, we compute polyhedral relations: for any pair of eigenvalues (or distinct entries) within \(\varvec{J}\), we find an over-approximation of the convex hull containing the points
The reason for evaluating the convex hull over pairs of points is twofold. In the first instance, we note that the set \(\left\{ \varvec{m}^k \mid k \in [1,\ldots ,n] \right\} \) is, in general, not convex. This makes it hard to find its support in arbitrary directions. Ideal directions would be the normals to the gradients of the function, namely \(\nabla m^k\), which would provide the tightest over-approximation at iteration k. However, as will be seen below, when combining negative or complex conjugate eigenvalues, the corresponding hyperplane tangent may intersect the set, and thus it cannot be used to define its convex hull. The second reason for choosing pairwise directions is practical: we need an even distribution of the directions in \(\mathbb {R}^p\), and it is easier to do this in a pairwise manner.Footnote 1
4.1 Abstract Matrices in Complex Spaces
To deal with complex numbers in eigenvalues and eigenvectors, [44] employs the real Jordan form for conjugate eigenvalues \(\lambda = re^{i\theta }\) and \(\lambda ^* = re^{-i\theta }\) (\(\theta \in [0, \pi ]\)), so that
Although this equivalence will be of use once we evaluate the progression of the model, calculating powers under this notation is often more difficult than handling directly the original matrices with complex values.
In the case of real eigenvalues we have abstracted the entries in the power matrix \(\varvec{J}_s^n\) by ranges of eigenvalues \([\min \{\lambda _s^0 \cdots \lambda _s^n\}, \max \{\lambda _s^0 \cdots \lambda _s^n\} ]\), forming a hypercube. In the complex case, where the rotations describe spherical behaviour, we can do something similar by rewriting eigenvalues into the polar form \(\lambda _s = r _s e^{i\theta _s}\) and enclosing the radius in the interval \([0, {\overline{r}}_s]\), where \({\overline{r}}_s=\max \{r_s^k : k \in [0,\ldots ,n]\}\) (in the worst case scenario this is over-approximated by a hyper-box with \(\lambda _s^k \in [-{\overline{r}}_s, {\overline{r}}_s]+[-{\overline{r}}_s, {\overline{r}}_s]i\), but we will introduce tighter bounds in the course of this work).
5 Abstract Acceleration Without Inputs
5.1 Using Support Functions for Abstract Acceleration
As an improvement over [44], the rows in \(\varvec{\varPhi }\) and \(\varvec{f}\) (see (13)) can be obtained by a refined sampling of the support functions of these sets. The choice of directions for these support functions results in an improvement over the logahedral abstractions used in previous work [37, 38, 44] (see Figs. 4,5, 6 and 7). This approach works thanks to the convex properties of the exponential progression. We consider five cases:
-
1.
Positive Real Eigenvalues The exponential curve is cut along the diagonal between the eigenvalues with maximum and minimum range to create a supporting hyperplane. A third point taken from the curve is used to test the direction of the corresponding template vector. An arbitrary number of additional supporting hyperplanes are created by selecting pairs of adjacent points in the curve and creating the corresponding support functions, as illustrated in Fig. 4.
-
2.
Complex Conjugate Eigenvalue Pairs In the case of complex conjugate pairs, the eigenvalue map corresponds to a logarithmic spiral (Fig. 5). In this case, we must first extract the number of iterations (denoted by \({\overline{k}}\)) required for a full cycle. For convergent eigenvalues (\(|{\lambda }|<1\)), only the first \({\overline{k}}\) iterations have an effect on the support functions, while in the divergent case only the last \({\overline{k}}\) iterations are considered (since symmetrically, it corresponds to the reverse spiral case). Support functions are found for adjacent pairs, checking the location of the first point for convergent eigenvalues, and that of the last point for divergent eigenvalues. If a point falls outside of the supporting half-space, we look for an interpolant point that closes the spiral and that is tangent to the origin. This last check is performed as a binary search over the remaining points in the circle (noting that the supporting planes would exclude the considered point) to achieve maximum tightness (Fig. 5).
Care is taken to ensure that subsequent iterations fall within the envelope found on the first/last rotation. This is ensured by extending the support functions outwards by a factor \(f = \max \left( \{1,|\lambda |^{{\hat{n}}}\cos (\theta )^{-1}\}\right) \), where \(\theta \) is the angle of the eigenvalue pair and \({\hat{n}}=n\) for the convergent case or \({\hat{n}}=\frac{1}{n}\) for the divergent case. When this value is too large, we use an interpolation to find better supports. This is achieved by finding a pair such that the first point is obtained from \(\lambda ^k\) and the second from \((\lambda ^{\frac{1}{m}})^{mk+1}\). The relaxation factor then becomes \(\cos \left( \frac{\theta }{m}\right) ^{-1}\).
-
3.
Equal Eigenvalues When two eigenvalues are equivalent, the resulting support functions are those that are orthogonal to the \(x=y\) plane, intersecting the square created by the maximum and minimum values.Footnote 2
-
4.
Jordan Blocks of Non-trivial Size (\(>1\)) In the case of eigenvalues with geometric multiplicities, we find three shapes. When both elements in the pair are convergent (convex sets can be “sharp”), it is important to find the apex of the upper diagonals in order to minimise the over-approximation (Fig. 6). When both elements are divergent, the shape is similar to a positive valued pair since there is no extremum. Finally, when comparing different Jordan blocks, one convergent and one divergent, we evaluate the enclosing hyperbox, thus avoiding the change in convexity at the apex.
-
5.
Negative Eigenvalues and Combinations of Real Eigenvalues with Conjugate Pairs When comparing a positive real eigenvalue to a complex conjugate or a negative one, we must account for the changes of sign in the progression of the latter. We compute envelopes of the progression of the corresponding dynamics, which are obtained via convex over-approximations (cf. Fig. 7). In the case of complex eigenvalues, we use the absolute value in order to determine the envelope. If both eigenvalues have rotating dynamics, we would require full symmetry along the four quadrants, and thus we obtain a hyper-box with vertices at the farthest points from the origin.
Polyhedral faces over \(\mathbb {R}^2\) for pairs of eigenvalues \((\lambda _1^n, \lambda _2^n)\) where \(\lambda _1=2, \lambda _2 =3\), and \(1{\le }n{\le }5\). Bold purple lines represent supports. The dotted grey and dashed red polytopes show logahedral approximations (box and octagon) used in [44]. Note the scales (the sloped dashed lines are parallel to the x = y line, and the dashed red polytope hides two small faces yielding an octagon). (Color figure online)
Polyhedral faces projected onto \(\mathbb {R}^2\) for complex conjugate eigenvalues \((\lambda _1^n, \lambda _2^n)\) where \(\lambda _1=0.8+0.4i, \lambda _2=0.8-0.4i\), and \(1{\le }n{\le }14\). Bold purple lines represent supports. The blue dotted line shows the supporting hyperplane that excludes the point obtained with \(n=1\), which is replaced by a supporting hyperplane tangent to the spiral but touching said origin. (Color figure online)
Polyhedral faces in \(\mathbb {R}^2\) related to a Jordan block \((\lambda _1^n, \left( {\begin{array}{c}n\\ 1\end{array}}\right) \lambda _1^{n-1})\), where \(\lambda _1=0.8\) and \(1{\le }n{\le }15\). Bold purple lines represent supports found in this work. (Color figure online)
Polyhedral faces over \(\mathbb {R}^2\), of different eigenvalue ratios (note that the curves obtained from the progression of the blue and orange dots are convex w.r.t. the \(\lambda _2^n\)-axis, whereas the green and brown are concave). Dotted lines represent convex supports for these layouts. (Color figure online)
An additional drawback of [44] is that calculating the exact Jordan form of any matrix is computationally expensive for high-dimensional matrices. We will instead leverage numerical algorithms that provide an approximation of the Jordan normal form and soundly account for the associated numerical errors. We use properties of eigenvalues to relax \(\varvec{f}\) by finding the maximum error in the calculations, which can be determined by computing the norm \(\delta _{max} = \Vert \hat{\mathbb {S}}\hat{\mathbb {J}}\hat{\mathbb {S}}^{-1}-\mathbb {A}\Vert \), where \(\hat{\mathbb {J}}\) and \(\hat{\mathbb {S}}\) are the eigenvalues and eigenvectors of \(\varvec{A}\) calculated numerically [16]. Recall that the notation above is used to represent interval matrices, and that all operations are performed using interval arithmetic with outward rounding in order to ensure soundness. The constraints in \(\varvec{\varPhi } \varvec{m} \le \varvec{f}\) are then computed by considering the ranges of eigenvalues \(|\lambda _s \pm \delta _{max}|^k\), which are represented in Fig. 4 with blue circles.
The outward relaxation of the support functions (\(\varvec{f}\)) follows a principle similar to that in [33], and reduces the tightness of the over-approximation, while ensuring the soundness of the obtained abstract matrix \({\mathcal {A}}^n\). Graphically, this is equivalent to moving the faces of a polyhedron outward, which practically has a minimal impact due to the small magnitude of \(\delta _{max}\). It is also worth noting that the transformation matrices into and from the eigenspace will also introduce over-approximations due to the intervals, and will exacerbate the over-approximations due to the condition number related to the eigenvalues.
One can still use exact arithmetic with a noticeable improvement over previous work; however, for higher-dimensional models the option of using floating-point arithmetic (taking errors into account and meticulously setting rounding modes) provides a 100-fold improvement, which can render verification practically feasible. For a full description of the numerical techniques see [16].
5.2 Abstract Matrices and Support Functions
Since we are describing operations using abstract matrices and support functions, we briefly review the nature of these operations and the properties that support functions retain within this domain. Let \(X \in \mathbb {R}^p\) be a set and \({\mathcal {A}} \in {\mathcal {R}}^{p \times p}\) an abstract matrix for the same space. From Eq. (13) we have
which leads to
where \( \rho _{\varphi X}(\varvec{v}) = \sup \left\{ \rho _{\varphi }(\varvec{x} \circ \varvec{v}) : \varvec{x} \in X\right\} , \) and \( \rho _{\varphi }(\varvec{v}) = \sup \{\varvec{m} \cdot \varphi ^{-1}(\varvec{v}) : \varvec{\varPhi }\varvec{m} \le \varvec{f}\} \). Here, \(\varvec{x} \circ \varvec{y}\) is the Hadamard product, where \( (\varvec{x} \circ \varvec{y})_{i}=\varvec{x}_i\varvec{y}_i\), and \(\varphi ^{-1}(\cdot )\) is the inverse operation of \(\varphi (\cdot )\). We also define
and, in order to simplify the nomenclature, we write
6 General Abstract Acceleration with Inputs
6.1 Acceleration of Parametric Inputs
Let us now consider the following over-approximation for \(\tau \) on sets:
and add a restriction to constant (also called parametric) inputs, namely \(\varvec{u}_k = \varvec{u}_0, \forall k>0\) and \(\varvec{u}_0 \in U\). Unfolding (3) (ignoring the presence of the guard set G for the time being), we obtain
We further simplify the sum \(\sum _{k=0}^{n-1} \varvec{A}^k\varvec{B}U\), exploiting the following result from linear algebra.
Lemma 2
If \(\varvec{I}-\varvec{A}\) is invertible, then
If furthermore \(\lim \limits _{n\rightarrow \infty } \varvec{A}^n=0\), then \(\lim \limits _{n \rightarrow \infty } \sum \limits _{k=0}^n \varvec{A}^k = (\varvec{I}-\varvec{A})^{-1}.\)
The inverse \((\varvec{I}-\varvec{A})^{-1}\) does not exist for eigenvalues equal to 1, i.e., we need \(1\notin \sigma (\varvec{A})\), where \(\sigma (A)\) is the spectrum (the set of all the eigenvalues) of matrix \(\varvec{A}\). In order to address this problem, we introduce the eigendecomposition of \(\varvec{A} = \varvec{S}\varvec{J}\varvec{S}^{-1}\), (and trivially \(\varvec{I}=\varvec{S}\varvec{I}\varvec{S}^{-1}\)), and by the distributive and transitive properties we obtain
Although \((\varvec{I}-\varvec{J})\) is still not invertible, this representation allows us to accelerate the eigenvalues individually, trivially noticing that \(\sum _{k=0}^{n-1} 1^k = n\) for unitary eigenvalues (thus eliminating the need to calculate said inverse for these eigenvalues). Using the properties above, and translating the problem into the generalised eigenspace to account for unit eigenvalues, we obtain the following representation:
given
where
and where \(gm(\cdot )\) denotes the geometric multiplicity of the given eigenvalue, and \(k=j-i\). With these notions in hand, we next define the abstract acceleration of parametric inputs.
Theorem 2
The abstract acceleration is defined as
where \({\mathcal {B}}^n\supseteq \bigcup _{k \in [1,\ldots ,n]}\varvec{S}(\varvec{D}^k)\varvec{S}^{-1}\varvec{B}\), is an over-approximation of the n-reach tube, namely \({\hat{X}}_n \subseteq {\hat{\tau }}^{\sharp n}(X_0,U)\).
Proof
From Eq. (4) we have
Using Eq. (17), we expand this into
then replace
and finally obtain
The quantities \({\mathcal {A}}^n\) and \({\mathcal {B}}^n\) are calculated using the techniques described in Sect. 5.1, where special consideration is taken to evaluate pairs comprising equal eigenvalues. Figure 8 gives an example of such a pair. Since both functions are monotonic, the set is convex. The techniques applied to positive real eigenvalues (see Sect. 5.1) therefore stands. \(\square \)
Polyhedral faces \((\lambda _1^n, n)\) over \(\mathbb {R}^2\), where \(\lambda _1=2, \lambda _2=3\), and \(1{\le }n{\le }5\). Bold purple lines represent supports found in this work. (Color figure online)
6.2 Acceleration of Time-Varying Inputs
In order to generalise the previous results for use with time-varying inputs, we will over-approximate the term \(\varvec{B}U\) over the eigenspace by a semi-spherical enclosure, namely a set where complex conjugate eigenvalues are enclosed by a spherical support with centre \(\varvec{u}_c'\) and the radius of \(U_b'\), whereas real eigenvalues are enclosed by hyper-rectangles (dashed symbols represent elements in the eigenspace). To this end, we first rewrite
where \(\varvec{u}_c'\) is the centre of the smallest hyperbox (interval hull) containing \(U_J'\), and \(U_d'=\{ \varvec{u} : \varvec{u}+\varvec{u}_c' \in U_J'\}\) as:
We then over-approximate \(U_d'\) via \(U_b'\), by the maximum radius in the directions of the complex eigenvalues (cf. illustration in Fig. 9). Let
be the set of eigenvalues of \(\varvec{A}\), where \(im(\cdot )\) is the imaginary value of a complex number, and conjugate pairs are represented only by one member of the pair. Let us define the function \(f_b:\mathbb {R}^p \rightarrow \mathbb {R}^{p_b}\), where \(p_b\) is the cardinality of \(\varLambda \), such that
\(i \in [1,\ldots ,p]\) and red(\(\cdot \)) is a function that reduces the dimension of a vector by removing the elements where \( \lambda _i \notin \varLambda \) (i.e. the null elements in \(\varvec{v}_b\), such that if for instance \(\varvec{v}_b = [\begin{array}{ccccc}v_1&0&v_3&\ldots&v_p\end{array}]^\intercal \), then \(\text {red} (\varvec{v}_b) = [\begin{array}{cccc}v_1&v_3&\ldots&v_p\end{array}]^\intercal \)). Extending this to matrices we have
where r denotes the number of inequalities describing a set in \(\mathbb {R}^p\). Finally,
Since the description of \(U'_b\) is no longer polyhedral in \(\mathbb {R}^p\), we will also create an over-approximation \(\varvec{J}_b\) of \(\varvec{J}\) in the directions of the complex eigenvectors, in way that is similar to how we generated \(U_b'\) for \(U_d'\). More precisely,
and \(gm(\cdot )\) is the geometric multiplicity of the specific Jordan block.
Definition 3
Given a matrix \(\varvec{A}=\varvec{S}\varvec{J}\varvec{S}^{-1}\) and a vector \(\varvec{x}\), we define the following operations:
Finally, we refer to the accelerated sets
Returning to our original equation for the n-reach set, we obtainFootnote 3
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
Shifting our attention from reach sets to reach tubes, we can now over-approximate the reach tube by abstract acceleration of the summands in (24), as follows.
Theorem 3
The abstract acceleration
where \({\mathcal {A}}^n\supseteq \bigcup _{k \in [1,\ldots ,n]}\varvec{A}^k\), \({\mathcal {B}}^n\supseteq \bigcup _{k \in [1,\ldots ,n]}\sum _{i=0}^{k-1}\varvec{A}^i\varvec{B}\), and \({\mathcal {B}}_b^n\supseteq \bigcup _{k \in [1,\ldots ,n]}F_b^*\left( \sum _{i=0}^{k-1}\varvec{A}^i\varvec{B},\varvec{x}\right) \), denotes an over-approximation of the n-reach tube, namely \({\hat{X}}_n \subseteq {\hat{\tau }}^{\sharp n}(X_0,U)\).
Proof
From Eq. (24) we have that \(X_n \subseteq \varvec{A}^nX_0 \oplus U_{cb}^n = \varvec{A}^nX_0 \oplus U_c^n \oplus U_b^n\). Furthermore, from Eq. (24) we also have \({\hat{\tau }}^{\sharp n}(X_0,U_c) \supseteq \bigcup _{k \in [1,\ldots ,n]} \varvec{A}^kX_0 \oplus U_c^k\). Finally, from the definition of \({\mathcal {B}}_b^n\) we have \({\mathcal {B}}_b^nU_b \supseteq \bigcup _{k \in [1,\ldots ,n]}U_b^k\), hence \({\hat{\tau }}^{\sharp n}(X_0,U) \supseteq {\hat{X}}_n\). \(\square \)
6.3 Combining Abstract Matrices
Calculating the reach set from the set of initial states and that originating from the input set separately, and then adding them together, can result in coarse over-approximations. To minimise this effect, we apply abstract acceleration to the combined input-and-state spaces.
One important property of the abstract matrices \({\mathcal {A}}^n\), \({\mathcal {B}}^n\) and \({\mathcal {B}}_b^n\) is that they are related. In the case of parametric inputs, this correlation is linear and is covered by the acceleration defined in Lemma (). In the case of \({\mathcal {B}}_b^n\) this relationship is not linear (see Eq. 21). However, we can still find a linear over-approximation of the relation between \({\mathcal {B}}_b^n\) and \({\mathcal {A}}^n\) based on the time steps k.
Given two sets \(X \in \mathbb {R}^p\) and \(U \in \mathbb {R}^q\) and a transition equation \(X_{k+1}=\varvec{A}X_k + \varvec{B}U\), which is related to \(\rho _{X_{k+1}}(\varvec{v})=\rho _{\varvec{A}X_k}(\varvec{v}) + \rho _{\varvec{B}U}(\varvec{v})\), we define a set
such that
Accelerating \(X_{k+1}\), we obtain
in the case of parametric inputs. More generally, the diagonal elements of \(\varvec{D}^n\) correspond to the diagonal elements of \(\varvec{A}^n\) and \(\sum _{k=0}^{n-1}{\varvec{A}^k} \varvec{B}\), which means we can construct
where \({\mathcal {A}}^n\) and \({\mathcal {B}}^n\) are the abstract matrices in Eqs. (13) and (19). We can then apply this abstraction to (21) and obtain:
with \(\varvec{J}_b\) defined in (22). This model provides a tighter over-approximation than (25), since the accelerated dynamics of the inputs are now coupled to the acceleration of the dynamical part of the model.
Example 1
In order to illustrate this, let us consider the one-dimensional model \(\varvec{x}_{k+1}=0.5\varvec{x}+1\) with initial state \(\varvec{x}_0=1\). If we calculate \({\mathcal {A}}\) and \({\mathcal {B}}\) separately we get \(\hat{\varvec{x}}=\bigcup _{k=0}^{\infty } \varvec{A}^{k}\varvec{x}_0+\bigcup _{k=0}^{\infty } (1-\varvec{A}^{k})\frac{\varvec{u}}{1-\varvec{A}} = [1, 3]\), however, using \({\mathcal {D}}\) we have \(\hat{\varvec{x}}=\bigcup _{k=0}^{\infty } \varvec{A}^{k}\left( \varvec{x}_0-\frac{\varvec{u}}{1-\varvec{A}}\right) +\frac{\varvec{u}}{1-\varvec{A}} = [1, 2]\). \(\square \)
7 Abstract Acceleration with Guards: Estimation of the Number of Iterations
In the presence of spatial guards G, we are interested in estimating the number of iterations used to calculate the abstract matrices. Since we are dealing with reach sets, we differentiate between sets that are entirely inside the guard, sets that are crossing it, and sets that are entirely outside. The latter reach sets should never be propagated, whereas reach sets crossing guards should be made as tight as possible.
Given a convex polyhedral guard expressed as the assertion \(G=\{\varvec{x} : \varvec{G}\varvec{x} \le \varvec{h}\}\), we define \(G_{i,*}\) as the \(i{\text {th}}\) row of \(\varvec{G}\) and \(h_i\) as the corresponding element of \(\varvec{h}\). We denote the normal vector to the \(i{\text {th}}\) face of the guard as \(\varvec{g}_i=G_i^\intercal \). The distance of the hyperplane defined by the i-th guard to the origin is thus \(\gamma _i = \frac{h_i}{|\varvec{g}_i|}\).
Given a convex set X, we can now describe its position with respect to each face of the guard through the use of its support function alongside the normal vector to the hyperplane (for clarity, we assume that the origin is inside the set X):
Applying this to Eq. (24) we obtain:
From the inequalities above we can determine the number of iterations \(\underline{n_i}\) for which the reach tube remains inside the corresponding hyperplane, and from which iteration \(\overline{n_i}\) the corresponding reach set goes beyond the guard.
Therefore, in order for a reach set to be inside the guard it must 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 \({\underline{n}} = \min \{\ \underline{ n_i }\ \}\) and \({\overline{n}} = \min \{\ \overline{ n_i }\ \}\).
We now discuss why these two cases are important. Looking at the transition in Eq. (1), we can easily derive that if \(\varvec{G}\varvec{x}_k \not \le \varvec{h}\) (i.e., the point lies outside at least one of the faces of the guard set), the post-image of all subsequent iterations of that point must not be included. As such, any over-approximation of the reach set will only add imprecision. Therefore, we will use the bounds \({\underline{n}}\) and \({\overline{n}}\) to create a tighter over-approximation. Let
This double step prevents the set \(\left\{ \varvec{x} : \varvec{x} \in {\hat{X}}_{{\underline{n}}}^{\sharp }, \varvec{x} \notin X_{{\underline{n}}}^{\sharp } \right\} \) from being included in further computations, thus improving the precision of the over-approximation.
Computing the maximum \(\underline{n_i}\) such that Eq. (28) is satisfied is not trivial, because the unknown \(\underline{n_i}\) 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 \(\varvec{g}_i\) into the generalised eigenspace of \(\varvec{A}\). More precisely, let
where \(\varvec{v}_j\) are row vectors of \(\varvec{S}^{-1}\) or \(-\varvec{S}^{-1}\) such that \(k_{ij} \ge 0\), and \( res (\varvec{g}_i)\) is the component of the vector \(\varvec{g}_i\) that lies outside the range of \(\varvec{S}\), namely the subspace spanned by its columns. Notice that since by definition \(\varvec{S}\) always has an inverse, it is full rank and therefore \( res (\varvec{g}_i)=\varvec{0}\) and subsequently not relevant. It is also important to note that \(\varvec{S}\) is the matrix of the generalised eigenvectors of \(\varvec{A}\), and that we are therefore expressing the guard in the generalised eigenspace of \(\varvec{A}\). Thus, we obtain:
7.1 Overestimating the Number of Iterations of a Model Without Inputs
Since rotating dynamics and Jordan shapes will have a complex effect on the behaviour of the model, we seek to transform the Jordan form into a real positive matrix by using the absolute value of the eigenvalues. In such a case, the support function in each direction is monotonically increasing (or decreasing), and it is therefore very easy to find a bound for its progression. We note that the envelope (described by the absolute value) of rotating dynamics will always contain the true dynamics and is therefore a sound over-approximation. We will initially assume that \(\gamma _i\) is positive and then extend to the general case.
Let \(\rho _{X_0}({\varvec{A}^n}^\intercal \varvec{g}_i)=\rho _{X_0'}({\varvec{J}^n}^\intercal \varvec{g}_i')\), so that \(\varvec{g}_i'=\varvec{S}^{-1}\varvec{g}_i\) and
Further, let \(\varLambda _\sigma =\{\lambda _i : i\in [1,\ldots ,p],\ \bigwedge _{j=1}^{i-1} (\lambda _i^*\ne \lambda _j \wedge \lambda _i\ne \lambda _j)\}\), be the set of eigenvalues with distinct values (excluding conjugate pairs and geometric multiplicities). Introduce \(f_\sigma (\varvec{v}):\mathbb {R}^p \rightarrow \mathbb {R}^{p_b}\), where \(p_b\) is the cardinality of \(\varLambda _\sigma \), such that \(f_\sigma (\varvec{v}) = \text {red}(\varvec{v}_\sigma )\), and
and furthermore let \(F_\sigma : \mathbb {R}^{r \times p} \rightarrow \mathbb {R}^{r \times p_b}\) be
Above, red(\(\cdot \)) is a function that reduces the dimension p of a vector to \(p_b\) by removing the elements \( \lambda _i \notin \varLambda _\sigma \). This reduction is not strictly necessary, but it enables a faster implementation. Thus, given \(\varvec{J} = diag\left( \varvec{J}_s, s \in [1,\ldots ,p_b]\right) \), we have
where \({\overline{\sigma }}_s=\sup \limits _{\varvec{x} \ne 0}\frac{\Vert \varvec{J}_s\varvec{x} \Vert _2}{\Vert \varvec{x} \Vert _2}\) is the maximum singular value [48] of the Jordan block \(\varvec{J}_s\). Finally, let
and \(\varvec{v}_\sigma =f_\sigma (\varvec{v})\), where \(\varvec{x}_c'\) is the Cartesian centre of \(X_0'\) and \(X_{c\sigma }'\), an over-approximation of \(X_0'\) centred at \(\varvec{x}_c'\).
Using properties of eigenvalues and of singular values, we obtain \(\rho _{X_0}((\varvec{A}^n)^\intercal \varvec{v}_j) \le {\overline{\sigma }}_j^{\ n}\,\rho _{X_{c\sigma }}((\varvec{v}_\sigma )_j)\), where \(j\in [1,\ldots ,p_b]\), and therefore
where \(k_{ij}\) are the coefficients in Eq. (30).
Since we have assumed to have no inputs, \(\rho _{U_c^n}(\varvec{g}_i) + \rho _{U_b^n}(\varvec{g}_i)=0\), hence we may solve for \({\underline{n}}_i\) as:
In order to separate the divergent parts of the dynamics from the convergent one, let us define
This step will allow us to effectively track which trajectories are likely to hit the guard and when, since it is only the divergent element of the dynamics that can increase the size of the reach tube in a given direction. This condition requires that the set of initial conditions is also inside the guard, which is a reasonable assumption.
Substituting in Eq. (34), we obtain
which allows us to finally formulate the following iteration scheme for under-approximating n.
Proposition 1
An iterative under-approximation of the number of iterations n can be computed by starting with \(\underline{n_i}=1,\) and iterating over
substituting \(n_i=n\) on the right-hand side until we meet the inequality. (If \(n<0\) is obtained at the first iteration, then \({\underline{n}}_i=0\) is the obtained solution.)
Proof
Notice that the sequence \(\underline{n_i}\) is monotonically increasing, before it breaks the inequality. As such, any local minimum represents a sound under-approximation of the number of loop iterations. Note that in the case where \(\gamma _i \le 0\), we must first translate the system coordinates such that \(\gamma _i > 0\). This is simply done by replacing \(\varvec{x}'=\varvec{x}+\varvec{c}\) and operating over the resulting system where \(\gamma _i' = \rho _{\varvec{c}}(\varvec{g}_i)+\gamma _i\).
Mathematically this is achieved as follows: first we get \(\varvec{c}\) by finding the centre of the interval hull (see Eq. (20)) of G.Footnote 4 Next we transform the dynamics into
\(\square \)
7.2 Underestimating the Number of Iterations of a Model Without Inputs
In order to apply a similar technique to (29), we must find an equivalent under-approximation. In the case of Eq. (34), the quantities \({\overline{\sigma }}_j\) ensure that the equation diverges faster than the real dynamics, hence the estimation found is an upper bound to the desired iteration. In this case we want the opposite, hence we look for a model where the dynamics diverge slower. It is easy to show that \({\lambda _b}_j=|\lambda _j|\) represents these slower dynamics,
which reduces to
where \({\underline{k}}_{\ ij}^-=\min \left\{ k_{ij}\ \rho _{X_{c\sigma }}(-(\varvec{v}_\sigma )_j)\ ,\ 0 \right\} \) and \({\underline{k}}_{\ ij}^+=\max \left\{ k_{ij}\ \rho _{X_{c\sigma }}(-(\varvec{v}_\sigma )_j)\ ,\ 0 \right\} \).
An additional consideration must be made regarding the rotational dynamics. In the previous case we did not care about the rotational alignment of the set \(X_n\) with respect to the vector \(\varvec{g}_i\), because any rotation would remain inside the envelope corresponding to the absolute value (\(r^k cos( k \theta ) \le r^k\) ). In the case of an under-approximation, although the magnitude of a complex eigenvalue at a given iteration may be greater than the support of the guard under verification, its angle with respect to the normal to the support vector may cause the corresponding point to remain inside the guard. We must therefore find iterations that are aligned with the normal to the guard, thus ensuring that the chosen point is outside it. 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 rotation to make our equation strictly divergent. Let \({\underline{\theta }}=\min \{\theta _j, j \in [1,\ldots ,p]\}\), where \(\theta _j\) are the angles of the complex conjugate eigenvalues. Let \(n_\theta =\frac{2\pi }{{\underline{\theta }}}\) be the maximum number of iterations needed for any of the dynamics to complete a full turn. Then at any given turn \(|\lambda _j|^{{\overline{n}}_i+n_\theta } \le |\lambda _j|^{{\overline{n}}_i+n}, \text { where } |\lambda _i| \le 1 \text { and } n \in [0, n_\theta ]\). This means that any bound we find on the iterations will be necessarily smaller than the true value. Our problem becomes the solution to:
The problem is simplified by soundly under-approximating the cosines and removing the constants, namely deriving the expressions
The solution to the last equation is
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 \(\theta _j \ge \frac{\pi }{2}\). In such cases, we must explicitly evaluate the set up to \(\frac{2\pi }{\theta _j}+1\) iterations after \({\overline{n}}_i\), in order to ensure that we have covered a full rotation. If the resulting bound does not satisfy the original inequality: \(\rho _{X_0}\left( (\varvec{A}^{\overline{n_i}})^\intercal \varvec{g}_i\right) \ge \gamma _i\), we replace \({\overline{n}}_i=n\) until it does.Footnote 5
Proposition 2
An iterative under-approximation of the number of iterations n can be computed by starting with \(\overline{n_i}'=0\) and iterating over
where k is the result of Eq. (38). We substitute for \({\overline{n}}_i=n\) on the right-hand side as long as the first inequality holds, and then we can find a k such that the second inequality holds.
Since we are explicitly verifying the inequality, there is no further proof required.
7.3 Estimating the Number of Iterations of a Model with Inputs
For an LTI model with inputs, we will use the same paradigm explained in the previous section after transforming the system with inputs into an over-approximating system without inputs.
Let \(X_{c\sigma }', U_{c\sigma }'\) be the corresponding sets of initial states and inputs obtained by applying Eq. (32) to \(X_0'\) and \(U_J'\), and let \(U_{J\sigma }'=(\varvec{I}-\varvec{J}_\sigma )^{-1} U_{c\sigma }'\). The accelerated resulting system may be represented by the equations
Let us now define \((XU)_{\sigma }=\{\varvec{x}-\varvec{u} \mid \varvec{x} \in X_{c\sigma }', \varvec{u} \in U_{J\sigma }' \}\), which allows us to translate the system into
which has the same shape as the equations in the previous section. We may now apply the techniques described above to find the bounds on the iterations.
7.4 Narrowing the Estimation of the Number of Iterations
The estimations above can be conservative, but we may obtain tighter bounds on the number of iterations. In the first instance, note that we have eliminated all negative terms in the sums in Eq. (36). Reinstating these terms can result in loss of monotonicity, but we may still create an iterative approach by fixing the negative value at intermediate stages. Let \({\underline{n}}_i\) be our existing bound for the time horizon before reaching a guard, and \({\underline{k}}_{{\underline{n}}_i}=\sum _{j=1}^p {\underline{k}}_{ij}\left( \frac{ {\overline{\sigma }}_j}{ {\overline{\sigma }}}\right) ^{{\underline{n}}_i}\), \({\overline{k}}_{{\underline{n}}_i}=\sum _{j=1}^p {\overline{k}}_{ij}\left( \frac{ {\overline{\sigma }}_j}{ {\overline{\sigma }}}\right) ^{{\underline{n}}_i}\) the corresponding negative and positive terms of the equation. We can now find upper and lower bounds for \({\underline{n}}_i\) by replacing them in Eq. (39) as:
where \({\underline{n}}_k\) is the bound found in the previous stage. Some steps of this process will provide an unsound result, however, every second step will provide a monotonically increasing sound bound which will be tighter than the one in Eq. (36). Since the elements of the sums are convergent, we have that \(n_i \ge n_k\) implies \({\underline{k}}_{n_i} \ge {\underline{k}}_{n_k} \left( \text {i.e., }|{\underline{k}}_{n_i}| \le |{\underline{k}}_{n_k}| \right) \), thus
which means that \(n_k\) in Eq. (42) is smaller than our n in Eq. (36) (\(n_k \le n \le \underline{n_i} \text { and } \underline{n_i} \ge \underline{n_k})\).
In the case of Eq. (39), the explicit evaluation of the guard at each cycle ensures the behaviour described here.
7.4.1 Maintaining Geometric Multiplicity
A second step in optimising the number of iterations comes from adding granularity to the bounding abstraction by retaining the geometric multiplicity using the matrix \(\varvec{J}_b\) (see Eq. (22)).
Lemma 3
Given a matrix \(\varvec{A}\) with eigenvalues \(\{\lambda _s, s\in [1,\ldots ,r]\}\), where each eigenvalue \(\lambda _s\) has a geometric multiplicity \(p_s\) and corresponding generalised eigenvectors \(\{\varvec{v}_{s,i}, i \in [1,\ldots ,p_s]\}\),
Proof
By definition, given an eigenvector \(\varvec{v}_s\) of \(\varvec{A}\), then \(\varvec{A}\varvec{v}_s=\lambda _s\varvec{v}_s\) [42]. Similarly, a generalised eigenvector \(\varvec{v}_{s,i}\) of \(\varvec{A}\) satisfies the equation \(\left( \varvec{A}-\lambda _s\varvec{I}\right) \varvec{v}_{s,i}=\varvec{v}_{s,i-1}\)Footnote 6. and \(\varvec{v}_{s,1}=\varvec{v}_s\), hence
From here we recursively expand the formula for \(\varvec{A}^{n-j-1}\varvec{v}_{s,i-1}\) and obtain:
\(\square \)
Let \(i'\) denote the position of \(f_b(\lambda _{j})\) within the block \(\varvec{J}_{bs}\) it belongs to, such that its corresponding generalised eigenvector is identified as \(\varvec{v}_{bs,i'}=f_b(\varvec{v}_j)\). Then
In order to manage the product on the right hand side we use slightly different techniques for over- and under-approximations. For \({\underline{n}}_i\) we first find an upper bound \({\underline{n}}_i'\) using Eq. (36) and \(k_{ij}={k_{ij}'}_0+{k_{ij}'}_m\), and then do a second iteration using \(k_{ij}={k_{ij}'}_0+{k_{ij}'}_m\prod \nolimits _{m=0}^{p_s-i'-1} ({\underline{n}}_i'-m)\), which ensures that the true value is below the approximation. In the case of \({\overline{n}}_i\), we also start with \(k_{ij}={k_{ij}'}_0+{k_{ij}'}_m\) and update it during the iterative process.
Example 2
Let us look at the following example, comprising matrices:
with \(\varvec{J}_\sigma \) calculated as in Eq. (31), initial condition \(\varvec{x}_0' = {\left[ \begin{array}{llllll}1&\quad 1&\quad 1&\quad 1&\quad 1&\quad 1\end{array}\right] }\) and guard set \(\varvec{G}\varvec{x} \le 300\) where
The progression of the support function of the reach sets along this vector and the corresponding bounds, as described in the previous section, are shown in Fig. 10.
Progression of the support function of a system for a given guard. The thick blue dots are real values. The dashed green line over-approximates the progression using singular values (Sect. 7.1), the dashed yellow line under-approximates them using eigenvalue norms (Sect. 7.2), whereas the continuous purple lines represent the tighter over-approximation maintaining the geometric multiplicity (Sect. 7.4.1). We can see how the purple line finds a better bound for \({\underline{n}}_i\), while the \({\overline{n}}_i\) bound is conservative for both approaches. (Color figure online)
Changing the eigenvalues to:
we obtain the results in Fig. 11. In this second case we can see that the rotational dynamics force an increase of the initially calculated iteration to account for the effects of the rotation. \(\square \)
Progression of the support function of a rotational system for a given guard. The thick blue dots are real values (negative values are missing due to the log scale). Continuous purple lines represent the over-approximation. The steep vertical line at iteration 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 procedure will either find that this boundary was met at iteration 9 or push it forward to iteration 19. (Color figure online)
7.5 Case Study
We have selected a known benchmark from the literature to illustrate the discussed procedure: the room temperature control problem [28]. The temperature (variable \( temp \)) of a room is controlled via a user-defined input (\( set \)), which can be changed at any discrete time step through a heating (\( heat \)) element, and is affected by ambient temperature (\( amb \)) that is out of the control of the model.
We formalise the description of such a system both via a linear loop and with a dynamical model. Observe that since such a system may be software controlled, Algorithm 1 gives a pseudo-code fragment for the temperature control problem.
We use the read function to represent non-deterministic values between 0 and the maximum value given as its argument. Alternatively, this loop corresponds to the following hybrid dynamical model:
with initial condition \(\left[ \begin{array}{c} temp \\ heat \\ \end{array} \right] _{0} \in \left[ \begin{array}{c} \left[ 5, 40 \right] \\ \left[ 0, 1 \right] \\ \end{array} \right] \), non-deterministic inputs \( \left[ \begin{array}{c} amb \\ set \\ \end{array} \right] _{k} \in \left[ \begin{array}{c} \left[ 5, 40 \right] \\ \left[ 0, 300 \right] \\ \end{array} \right] , \) and guard set \(G = \left\{ \left[ \begin{array}{c} temp \\ heat \\ \end{array} \right] : \left[ \begin{array}{ll} 1 &{}\quad 0\\ 0 &{}\quad 1\\ \end{array} \right] \left[ \begin{array}{c} temp \\ heat \\ \end{array} \right] < \left[ \begin{array}{c} 400\\ 300\\ \end{array} \right] \right\} \).
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 \(\pm 10^{-19}\), moreover the error of the approximate Jordan form computation results in \(\delta _{max}<10^{-17}\). The eigendecomposition of the dynamics is (the values are rounded to three decimal places):
The discussed over-approximations of the reach-sets indicate that the temperature variable intersects the guard set G at iteration \({\underline{n}}=32\). Considering the pseudo-eigenvalue matrix along these iterations, we use Eq. (13) to find that the corresponding complex pair remains within the following boundaries:
The reach tube is calculated by multiplying these abstract matrices with the initial sets of states and inputs, as described in Eq. (25), by the following inequalities:
The negative values represent the lack of restriction in the code on the lower side and correspond to a cooling system (negative heating). The set is displayed in Fig. 12, where for the sake of clarity only 8 directions of the 16 constraints are shown. This results in a rather tight over-approximation which, for comparison’s sake, is not looser than the convex hull of all reach sets obtained by [31] using the given directions. Figure 12 further displays the initial set in black, the collection of reach sets in white, the convex hull of all reach sets in dark blue (as computed by [31]), and finally the abstractly accelerated set in light yellow (dashed lines). The outer lines represent the guards for G.
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. (Color figure online)
8 Application of Abstraction-Refinement to Abstract Acceleration
One of the main limitations of abstract acceleration is that, despite being very fast, it leverages an over-approximation of the actual reach tube for verification. In many cases this over-approximation can be too coarse (i.e., imprecise) for the proof of the safety property of interest. This section deals with methods for refining this over-approximation. Refinements are based on counterexamples, namely vertices of the abstract matrix that lay outside the projection of the safety specification onto the abstract space (calculated using the inverse of the reachability transformations). Our approach can be seen as an instance of the known CounterExample Guided Abstraction Refinement (CEGAR) paradigm [20].
8.1 Finding Counterexample Iterations
Because the objective is to refine the abstract dynamics, we need to find the iterations corresponding to the counterexample (i.e., the ones used to calculate the hyperplanes forming the unsafe vertex). This will allow us to find an interpolant iteration that will reduce the polyhedron in the right direction. Since the abstract dynamics are built over pairs of eigenvalues, it is possible that different eigenvalue pairs provide different results for the counterexample iteration, in which case all of them are used. Let a verification clause explore the solution \(\rho _{\mathcal {A}}(\varvec{v}) = s \le {\overline{s}}\), where \(\varvec{v}\) is the direction we are examining, s its corresponding support function, and \(\rho _{\mathcal {A}}(\varvec{v}) \le {\overline{s}}\) the safety specification. If \(s>{\overline{s}}\) the specification will not be met and we need a refinement. Let \(\varvec{a}_v \in {\mathcal {A}}\) be the vertex at which the maximum is found, i.e., \(\varvec{a}_v \cdot \varvec{v}=s\). The iterations corresponding to this counterexample may be found by analysing the dynamics of each pair of eigenvalues independently, and finding the point closest to the hyperplane whose inequality has been violated. This is done as follows:
-
1.
Conjugate Eigenvalues Since the trajectories along these are circular and centred at the origin, we can find the angle that \(\varvec{a}_v\) forms with the axes of the eigenvalues and use it to calculate the right iteration. Let \(\theta _i\) be the angle of the conjugate eigenvalue pair and \(\theta _{\varvec{a}_v}(i)\) the angle formed by \(\varvec{a}_v\) in the \(i{\text {th}}\) plane (this is equivalent to \(\tan ^{-1}\left( \frac{(\varvec{a}_v)_i}{(\varvec{a}_v)_{i+1}}\right) \)). The corresponding iteration will depend on whether the eigenvalue is convergent or divergent. In the former case, it will be \(\frac{\theta _{\varvec{a}_v}(i)}{\theta _i}\), and in the latter it will be \((n - (n \mod \frac{1}{\theta _i}))+\frac{\theta _{\varvec{a}_v}(i)}{\theta _i}\), where \(\hbox {mod}\) is the modulus operation over the reals.
-
2.
Real Eigenvalues In the case of reals, finding the iteration relies on the direct relation between the given eigenvalue and the target counterexample. Since \((\varvec{a}_v)_i \approx \lambda _i^k \Rightarrow k \approx \log _{\lambda _i}((\varvec{a}_v)_i)\). If the logarithm does not exist, then we presume we cannot further refine using this method.
-
3.
Jordan Blocks with Non-unitary Geometric Multiplicity In the case of larger Jordan blocks we need to examine the nature of the dynamics. Let us look at the equation representing the contribution of a Jordan block to the support:
$$\begin{aligned} \rho _{\lambda _s}(\varvec{v})=\sum _{j=0}^{p_s} \left( {\begin{array}{c}n\\ j\end{array}}\right) \lambda ^{n-j} {\varvec{v}_s}_j. \end{aligned}$$(45)In this case we must use an iterative approximation, as described in Sect. 7.4.1, to find the closest iteration to the unsafe guard. Although this process is more costly than the ones described above, it is also more precise, thus providing a much better refinement. Note that the technique can be applied to the full set of eigenvalues or to any subset of Jordan blocks. This choice is a compromise between precision and speed. We also note that when the refinement process is done in the eigenspace, the new eigenvectors are now the identity set, which makes the problem more tractable.
Since the exclusion of an unsafe vertex from the abstract dynamics does not ensure a sufficiently tight over-approximation, we must perform this step iteratively until either we run out of new refinements or reach a user-defined timeout.
Once the candidate iterations are found, it suffices to add further constraints to the abstract matrix for these iterations as described in Fig. 13. Notice that given the above procedure, it is often faster and more beneficial to begin by performing the refinement over the complex eigenvalues by directly examining the vector of directions \(\varvec{v}\) in the corresponding sub-spaces.
Polyhedral faces from an \(\mathbb {R}^2\) Jordan block subspace \((\lambda _1^n, \left( {\begin{array}{c}n\\ 1\end{array}}\right) \lambda _2^{n-1})\) where \(\lambda _1=0.8, \lambda _2=0.8\) and \(1 \le n \le 15\). The red dot specifies an abstract vertex violating the safety specification (dashed blue line). The closest iteration to the violating vertex is \(n=3\). A new support function (green) based on \(n=3\) eliminates the violating vertex. The new abstract polyhedron meets the safety specification (yellow). (Color figure online)
8.2 Using Bounded Concrete Runs
Owing to the complex interactions in the dynamics, the above procedure may not always find the correct iterations for refinement, or at least not optimal ones. For this reason, a second method is proposed, which in most cases will be more efficient and precise when the dynamics are strictly convergent.
This second approach relies on the direct calculation of the initial \({\overline{k}}\) iterations. Since we operate over the eigenvalues and we limit \({\overline{k}}\) to a conservative bound, this is a relatively inexpensive calculation. The approach leverages the idea that for convergent dynamics, counterexamples are often found in the initial runs. The first step is to directly calculate the trajectory of the counterexample for the first \({\overline{k}}\) iterations, and its corresponding support function in the direction of \(\varvec{v}\). Once again, because this is a single point and a bounded time, this operation is comparatively inexpensive. The second step consists of finding an upper bound for all subsequent iterations, which we can do by using the norms of the eigenvalues and the peaks of each geometrical multiple of a Jordan Block (which relate to these norms). By selecting the larger between these two supports, we ensure soundness over the infinite time horizon. This is equivalent to evaluating the reach tube as
Since the result above is known to be an upper bound for the support in the direction of \(\varvec{v}\), we can directly add it to the inequalities of \({\mathcal {A}}\).
8.3 Case Study
We have taken an industrial benchmark ‘CAFF problem Instance: Building’.Footnote 7 The benchmark consists of a continuous model with 48 state variables and one input. Furthermore, there is an initial state corresponding to a 10-dimensional hyper-rectangle. Time is discretised using a 5 ms sample interval to give us a discrete time model for verification. Notice that the choice of sample time has very little effect on abstract acceleration. It mainly affects the requirement for floating point precision (as very small angles may require higher precision), and may have an effect on counterexample generation which can either decrease precision or increase time (i.e., we may relax the precision of the algorithm to gain speed or tightening at the cost of higher computation times). The provided model requires an analysis on the \(25{\text {th}}\) variable, with a safety specification requiring it to remain below .005. The problem has been verified using SpaceExFootnote 8 for bounded time (20 s) in under three seconds. Axelerator was run on this benchmark using different parameters. We used an Intel 2.6 GHz i7 processor with 8 GB of RAM running on Linux. Although the algorithm lends itself to concurrency, the tool currently supports only single threading. The process itself uses 82 MB on this particular benchmark. The results are summarised in Table 1. Since many tools in this area use unsound arithmetic, we present results for both sound and unsound arithmetic using abstract acceleration to quantify the cost of soundness. It is worth noting that for precisions under 1024 bits the tool returns soundness violation errors when using sound arithmetic.
We note in these results that the performance of Axelerator depends largely on the required level of refinement. State of the art tools can do bounded model checking faster than Axelerator, largely due to implementation optimisations (we expect a better simplex engine would allow us to be more competitive in this regard). This advantage disappears as soon as we require a larger time horizon for verification.
9 Experimental Results
The overall Abstract Acceleration procedure has been implemented in C++ using the eigen-algebra package (v3.2), with multiple precision floating-point arithmetic, and has been tested on a 1.6 GHz Core 2 Duo computer. Unless otherwise specified, we use the sound version of Abstract Acceleration without abstraction-refinement (as per Sect. 8). The tool, called Axelerator, and the corresponding benchmarks used to test it (including specifications of the initial states, input ranges and guard sets), are available at
Since many tools require the specification of a set of directions in which to perform verification, we have chosen to use an octahedral template set. That is the set of vectors that run either along an axis or at a 45 degree angle between two axes (equivalent to a set of inequalities \(\pm x_i \le bound_k\) or \(\pm x_i \pm x_j<bound_k\)). We also make use of interval hulls, which are defined as the smallest hyper-boxes that enclose a set (equivalent to \(\pm x_i \le bound_k\)).
9.1 Comparison with Other Unbounded-Time Approaches
In a first experiment we have benchmarked our implementation against the tools InterProc [43] and Sting [22]. We have tested these tools on different discrete time models, including guarded/unguarded ones, stable/unstable ones, and models with complex/real loops with inputs (details in Table 2). In the first instance, we can see that Axelerator is more successful in finding tight over-approximations that remain very close to the actual reach set. InterProc and Sting are each unable to find nearly 40% of the bounds (i.e., supports), meaning they report an infinite reach-space in the corresponding directions. In these instances, InterProc (owing to limitations related to widening) and Sting (owing to the absence of tight polyhedral, inductive invariants) are unable to infer finite bounds at all. The precision provided by our implementation comes at a reasonable price in terms of computational effort: Axelerator is approximately 10 times slower than InterProc, and in most cases faster than Sting.
Table 3 presents a comparison of our implementation using different levels of precision (long double, 256 bit, and 1024 bit floating-point precision) with the core procedure for abstract acceleration of linear loops without inputs (J) [44] (we use a version without inputs of the parabola, cubic and convoy models). This shows that our implementation gives tighter over-approximations on most benchmarks (column ‘improved’). While on a limited number of instances the current implementation is less precise (the lower right portion of Fig. 4 where the dotted red line crosses into the inside of our support gives a hint why this is happening), the overall increased precision results from mitigating the limitation on chosen directions caused by the use of logahedral abstractions (as done in previous work [44]).
At the same time, our implementation is faster than [44], partly owing to the use of numeric eigendecomposition (as opposed to symbolic), but mostly in view of the cost of increasingly large rational representations in [44] needed to maintain precision when evaluating higher order systems. The fact that many bounds have improved with the new approach, while speed has increased by several orders of magnitude, provides evidence of the advantages of the new approach.
The speed-up is due to the faster Jordan form computation, which takes between 2 and 65 seconds for [44] (using the ATLAS package), whereas our implementation requires at most one second. For the last two benchmarks, the polyhedral computations blow up in [44], whereas our support function approach shows only moderately increasing runtimes. The increase of speed is owing to multiple factors, as detailed in Table 4. The difference in precision of using long double vs. multiple precision floating-point arithmetic is negligible, as all results in the given examples match exactly to 9 decimal places.
9.2 Comparison with Bounded-Time Approaches
In a third experiment, we compare our method with the LGG algorithm [51] used by SpaceEx [31]. Since LGG provides the tightest supports for a given set of directions, this study indicates the quality of our over-approximation. In order to set up a fair comparison, we have provided the implementation of the native algorithm in [51] within our software. We have run both methods on the convoyCar example [44] with inputs, which presents an unguarded, scalable (the dimension can be increased by adding more cars), stable loop with complex dynamics, and we have focused focused on octahedral templates. For convex reach tubes, the approximations computed by Abstract Acceleration are reasonably tight in comparison to those computed by the LGG algorithm (see Table 5). However, by storing finite disjunctions of convex polyhedra, the LGG algorithm is able to generate non-convex reach tubes, which are arguably tighter in case of oscillating or spiralling dynamics. Still, in many applications abstract acceleration can provide a tight over-approximation for the convex hull of (possibly non-convex) reach tubes.
Table 5 provides the quantitative outcomes of this comparison. For simplicity, we present only the projection of the bounds along the variables of interest on a lower dimensional model. As expected, the LGG algorithm performs better in terms of tightness, however unlike our approach, its runtime distinctly increases with the number of iterations. Our implementation of LGG [51] using Convex Polyhedra with octahedral templates becomes slower than the abstractly accelerated version even for small time horizons (our implementation of LGG requires \(\sim \!4\) ms for each iteration on a 6-dimensional problem with octahedral templates). Faster implementations of LGG will change the point at which the speed trade-off happens, but abstract acceleration will always be faster for long enough time horizons.
The evident advantage of abstract acceleration is its speed over finite horizons without much loss in precision, and of course the ability to prove properties for unbounded-time horizons.
9.3 Comparison with Alternative Abstract Acceleration Techniques
Table 6 shows a comparison between our approach and [49]. As can be seen, our approach is not only faster but much more precise. The reasons for this are many-fold. In terms of speed, the fact that [49] uses a dynamical model that is double the dimension of the one presented here and that the algorithmic complexity to manipulate it is \({\mathcal {O}}(n^3)\) (n being the model dimension) will result in slower computations, even when using sparse matrices. In terms of precision, creating an over-approximation around the centre of the input set, as opposed to the origin, makes a considerable difference, which is increased by the fact that circular over-approximations are contained within the interval hulls used in [49], and are therefore up to \(\frac{4}{3}^n\) times smaller in volume. This last consideration is only relevant for rotating dynamics, whereas positive real eigenvalues exhibit the same precision in both approaches. Finally, notice that in the ConvoyCar2 example, where all original eigenvalues are convergent, the interval hull approach [49] creates one divergent eigenvalue, which causes a critical change in the behaviour of the dynamics that leads to unbounded results.
9.4 Scalability
Finally, in terms of scalability, we have an expected \({\mathcal {O}}(p^3)\) complexity with respect to the number of variables p, which derives from the simplex algorithm and the matrix multiplications in Eq. (25). We have parameterised the number of cars in the convoyCar example [44] (also seen in Table 3), and experimented with up to 33 cars (each car after the first requires three state variables, so that for example we obtain \((33-1)\times 3= 96\) variables for the 33-car case), using the same initial states for all cars. We report an average obtained from 10 runs for each configuration. These results demonstrate that our method can scale to industrial-size problems (Table 7).
10 Related Work
There are numerous approaches that tackle the safety verification problem for dynamical models. Time-bounded analysis is unsound in most cases, since it cannot reason about unbounded-time cases (we note that a proof that a given horizon is a completeness threshold [21] would restore such soundness but many tools do not attempt to find such a proof). Unbounded-time solutions are therefore preferred when soundness is required, although they are often either less precise or slower than bounded counterparts.
10.1 Time-Bounded Reachability Analysis
The first approach surrenders exhaustive analysis over the infinite time horizon, and restricts the exploration of model dynamics up to some given finite time bound. Decision procedures for bounded-time reachability problems have made much progress in the past decade, and computational algorithms with related software tools developed. Representatives are STRONG [25], HySon [12], CORA [2], HYLAA [5], Ariadne [7] and SpaceEx [31].
Set-based simulation methods generalise guaranteed integration of ODEs [11, 53] from enclosing intervals to relational domains. They use precise abstractions with low computational cost to over-approximate sets of reachable states up to a given time horizon. Early tools employed polyhedral sets (HyTech [41] and PHAVer [30]), polyhedral flow-pipes [18], ellipsoids [10] and zonotopes [35]. A breakthrough has been achieved by [36, 51], with the representation of convex sets using template polyhedra and support functions. This method is implemented in the tool SpaceEx [31], which can handle dynamical systems with hundreds of variables. It performs computations using floating-point numbers, which is a deliberate choice to boost performance. The downside of this choice is that its implementation is numerically unsound and therefore does not provide genuine formal guarantees, which is at the core of this work. More generally, most tools using eigendecomposition over a large number of variables (more than 10) happen to be numerically unsound owing to the use of unchecked floating-point arithmetic. Some tools (e.g., C2E2 [26]) add a robustness check to the reachability computations in order to restore soundness by over-approximating simulation-based reachability analysis. Another breakthrough in performance has been achieved by HYLAA [5], which is the first tool to solve high-order problems, including several with hundreds of state variables. [7, 17] deal with non-linear dynamics (which are beyond the scope of this work), and the latter does so soundly. Other approaches focus on bounded model checking of hybrid automata, and use specialised constraint solvers (iSAT [27], HySAT [29]), or SMT encodings [19, 39].
10.2 Unbounded Reachability Analysis
The second approach, epitomised by static analysis methods [40], explores unbounded-time horizons. It employs conservative over-approximations to achieve relative completeness (this ensures that we always find a fixpoint) and decidability over infinite time horizons.
Unbounded techniques attempt to infer or compute a loop invariant, i.e., an inductive set that includes all reachable states. If the computed invariant is disjoint from the set of bad states, this proves that the latter are unreachable and hence that the loop is safe. However, analysers frequently struggle to obtain an invariant that is precise enough, with acceptable computational costs. The problem is evidently exacerbated by non-determinism in the loop, which corresponds to the case of open systems (i.e., with inputs). Prominent representatives of this analysis approach include Passel [45], Sting [22] and abstract interpreters such as Astrée [8] and InterProc [43].
Early work in this area has used implementations of abstract interpretation and widening [23], which represent still the underpinnings of many modern tools. The work in [40] uses abstract interpretation with convex polyhedra over piecewise differential inclusions. A more recent approach uses a CEGAR loop to improve the precision of polyhedral abstractions [9] following an inductive sequence of half-space interpolants. Dang and Gawlitza [24] employ optimisation-based (max-strategy iteration) with linear templates for hybrid systems with linear dynamics. Relational abstractions [56] use ad-hoc “loop summarisation” of flow relations, while abstract acceleration focuses on linear relations analysis [37, 38], which is common in program analysis.
10.3 Abstract Acceleration
Abstract acceleration [37, 38, 44] 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 [58] but restricted to subclasses of linear loops, and later to general linear loops but without inputs [44].
The work presented in this article mitigates these limitations by presenting abstract acceleration for general linear loops with inputs [14, 15], developing numeric techniques for scalability and extending the domain to continuous time models.
The work in [50] puts forward improvements on [14], and its outcomes can be compared with the results in [15, 16]. Indeed, [50] proposes an alternative approach to Abstract Acceleration with inputs, which relies on the expansion of the dynamical equations to a model that is four times the original size (dimension). Although this model is mathematically more tractable, it is slower and more imprecise than the one presented in [15]. The reason for this is that while the latter uses a semi-circular over-approximation of the variable inputs, the former uses a rectangular over-approximation of the dynamics, which is in turn expanded by the acceleration. In particular, this rectangular over-approximation of the dynamics may cause a convergent model to become divergent, which leads to unbounded supports. There are cases in which this rectangular over-approximation can give a tighter result, especially if it is implemented using the symmetry properties described in [15], thus a combined approach may be more efficient. Despite the computational overhead, the handling of the guards presented in [50] has potential advantages over our approach: its description is formal, and indeed our procedure is in general not complete, resulting on occasion in an inability to find reasonably precise bounds. However, in a typical case our procedure works well and in the experiments we have run, it has indeed always been able to find the optimal solution.
11 Conclusions and Future Work
We have presented a thorough development of the Abstract Acceleration paradigm to guarded LTI models (namely, conditional linear loops) with inputs. We have extended existing work, which dealt only with autonomous models (i.e., models without inputs). We have decisively shown that the new approach over-competes state-of-the-art tools for unbounded-time reachability analysis in both precision and scalability. In particular, on one hand we claim full soundness (both algorithmic and numerical) of the proposed procedure and of its implementation. On the other hand, the new approach is capable of handling general unbounded-time safety analysis for large-scale open models.
Nested loops are outside of the scope of this contribution, but represent relevant models that ought to be considered. Further work is also needed to extend the approach to non-linear dynamics, which we believe can be explored via hybridisation techniques [3]. We also plan to formalise the framework for general hybrid models, namely dynamical models with multiple discrete locations and location-dependent dynamics under multiple guards.
Notes
A more elaborate search of a set of directions might leverage convex optimisation techniques, which come with a computational overhead: as such, they fall outside the scope of this work.
As an alternative, a space reduction method can be applied to remove one of the eigenspaces, and the original space can be restored after the reduced reach tube has been evaluated.
Note that although we are working in the eigenspace, these sets can be traced back to corresponding sets in the original state space such that \(U_b'=\varvec{S}^{-1}BU_b\), \(U_c'=\varvec{S}^{-1}BU_c\), and \(U_d'=\varvec{S}^{-1}BU_d\). Hence, this inclusion is also valid in the original state space.
If G is open in a given direction we may pick any number in that direction for the corresponding row of c.
This is a tighter value than that in [14], where we over-approximated using \(n_\theta =\frac{(2\pi )^m}{\prod _j \theta _j}\), where m is the number of conjugate pairs.
More generally \(\left( \varvec{A}-\lambda _s\varvec{I}\right) \varvec{v}_{s,i}=k_{s,i}\varvec{v}_{s,i-1}\), since any vector \(k \varvec{v}_{s,i-1}\) is also a generalised eigenvector of \(\varvec{A}\), but we select generalised eigenvectors such that \(k=1\).
References
Adimoolam, A., Dang, T.: Template complex zonotope based stability verification. In: Tarbouriech, S., Girard, A., Hetel, L. (eds.) Control Subject to Computational and Communication Constraints: Current Challenges, pp. 83–96. Springer, Berlin (2018)
Althoff, M.: An introduction to CORA 2015. In: ARCH@CPSWeek, EPiC Series in Computing, vol. 34, pp. 120–151. EasyChair (2015)
Asarin, E., Dang, T., Girard, A.: Hybridization methods for the analysis of nonlinear systems. Acta Informatica 43(7), 451–476 (2007)
Avis, D., Fukuda, K.: A pivoting algorithm for convex hulls and vertex enumeration of arrangements and polyhedra. Discrete Comput. Geom. 8(1), 295–313 (1992)
Bak, S., Duggirala, P.S.: HyLAA: a tool for computing simulation-equivalent reachability for linear systems. In: Alur, R. (ed.) Hybrid Systems: Computation and Control, pp. 173–178. ACM, New York (2017)
Bennett, S.: A brief history of automatic control. IEEE Control Syst. 16(3), 17–25 (1996)
Benvenuti, L., Bresolin, D., Collins, P., Ferrari, A., Geretti, L., Villa, T.: Assume-guarantee verification of nonlinear hybrid systems with Ariadne. Int. J. Robust Nonlinear Control 24(4), 699–724 (2014)
Blanchet, B., Cousot, P., Cousot, R., Feret, J., Mauborgne, L., Miné, A., Monniaux, D., Rival, X.: A static analyzer for large safety-critical software. In: Plan, C. (ed.) Programming Language Design and Implementation, pp. 196–207. ACM, New York (2003)
Bogomolov, S., Frehse, G., Giacobbe, M., Henzinger, T.A.: Counterexample-guided refinement of template polyhedra. In: International Conference on Tools and Algorithms for the Construction and Analysis of Systems, pp. 589–606. Springer, Berlin (2017)
Botchkarev, O., Tripakis, S.: Verification of hybrid systems with linear differential inclusions using ellipsoidal approximations. In: Alur, R. (ed.) Hybrid Systems: Computation and Control, LNCS, pp. 73–88. Springer, Berlin (2000)
Bouissou, O.: Analyse statique par interprétation abstraite de systèmes hybrides. Ph.D. thesis, École Polytechnique (2008)
Bouissou, O., Mimram, S., Chapoutot, A.: Hyson: set-based simulation of hybrid systems. In: Rapid System Prototyping (RSP), pp. 79–85. IEEE (2012)
Carathéodory, C.: Über den Variabilitätsbereich der Koeffizienten von Potenzreihen, die gegebene Werte nicht annehmen. Math. Ann. 64(1), 95–115 (1907)
Cattaruzza, D., Abate, A., Schrammel, P., Kroening, D.: Unbounded-time analysis of guarded LTI systems with inputs by abstract acceleration. In: Static Analysis Symposium, LNCS, vol. 9291, pp. 312–331. Springer, Berlin (2015)
Cattaruzza, D., Abate, A., Schrammel, P., Kroening, D.: Unbounded-time analysis of guarded LTI systems with inputs by abstract acceleration (extended version). Tech. rep., University of Oxford (2015). arXiv:1506.05607
Cattaruzza, D., Abate, A., Schrammel, P., Kroening, D.: Sound numerical computations in abstract acceleration. In: International Workshop on Numerical Software Verification, pp. 38–60. Springer, Berlin (2017)
Chen, X., Ábrahám, E., Sankaranarayanan, S.: Flow*: an analyzer for non-linear hybrid systems. In: Computer Aided Verification, LNCS, pp. 258–263. Springer (2013)
Chutinan, A., Krogh, B.H.: Computing polyhedral approximations to flow pipes for dynamic systems. In: CDC, pp. 2089–2094. IEEE Computer Society (1998)
Cimatti, A., Mover, S., Tonetta, S.: SMT-based verification of hybrid systems. In: AAAI Conference on Artificial Intelligence, p. 2100-2105. AAAI Press (2012)
Clarke, E., Grumberg, O., Jha, S., Lu, Y., Veith, H.: Counterexample-guided abstraction refinement. In: Computer Aided Verification, pp. 154–169. Springer, Berlin (2000)
Clarke, E.M., Kroening, D., Ouaknine, J., Strichman, O.: Computational challenges in bounded model checking. Int. J. Softw. Tools Technol. Transf. 7(2), 174–183 (2005). https://doi.org/10.1007/s10009-004-0182-5
Colón, M.A., Sankaranarayanan, S., Sipma, H.B.: Linear invariant generation using non-linear constraint solving. In: Computer-Aided Verification, pp. 420–432. Springer, Berlin (2003)
Cousot, P., Cousot, R.: Abstract interpretation: A unified lattice model for static analysis of programs by construction or approximation of fixpoints. In: Principles of Programming Languages, pp. 238–252 (1977)
Dang, T., Gawlitza, T.M.: Template-based unbounded time verification of affine hybrid automata. In: Asian Symposium on Programming Languages and Systems, LNCS, pp. 34–49. Springer, Berlin (2011)
Deng, Y., Rajhans, A., Julius, A.A.: STRONG: A trajectory-based verification toolbox for hybrid systems. In: Quantitative Evaluation of Systems, LNCS, vol. 8054, pp. 165–168. Springer, Berlin (2013)
Duggirala, P.S., Mitra, S., Viswanathan, M., Potok, M.: C2E2: A verification tool for Stateflow models. In: International Conference on Tools and Algorithms for the Construction and Analysis of Systems, pp. 68–82. Springer, Berlin (2015)
Eggers, A., Fränzle, M., Herde, C.: SAT Modulo ODE: A direct SAT approach to hybrid systems. In: Automated Technology for Verification and Analysis, LNCS, vol. 5311, pp. 171–185. Springer, Berlin (2008)
Fehnker, A., Ivancic, F.: Benchmarks for hybrid systems verification. In: Hybrid Systems: Computation and Control, LNCS, pp. 326–341. Springer, Berlin (2004)
Fränzle, M., Herde, C.: HySAT: an efficient proof engine for bounded model checking of hybrid systems. Formal Methods Syst. Des. 30(3), 179–198 (2007)
Frehse, G.: PHAVer: Algorithmic verification of hybrid systems past HyTech. In: Hybrid Systems: Computation and Control, LNCS, vol. 3414, pp. 258–273. Springer, Berlin (2005)
Frehse, G., Le Guernic, C., Donzé, A., Ray, R., Lebeltel, O., Ripado, R., Girard, A., Dang, T., Maler, O.: SpaceEx: scalable verification of hybrid systems. In: Computer-Aided Verification, LNCS, vol. 6806, pp. 379–395. Springer, Berlin (2011)
Fukuda, K., Prodon, A.: Double description method revisited. In: Combinatorics and computer science, pp. 91–111. Springer, Berlin (1996)
Gao, S., Avigad, J., Clarke, E.M.: \(\delta \)-complete decision procedures for satisfiability over the reals. In: Automated Reasoning, pp. 286–300. Springer, Berlin (2012)
Ghosh, P.K., Kumar, K.V.: Support function representation of convex bodies, its application in geometric computing, and some related representations. Comput. Vis. Image Underst. 72, 379–403 (1998)
Girard, A.: Reachability of uncertain linear systems using zonotopes. In: Hybrid Systems: Computation and Control, LNCS, vol. 3414, pp. 291–305. Springer, Berlin (2005)
Girard, A., Le Guernic, C., Maler, O.: Efficient computation of reachable sets of linear time-invariant systems with inputs. In: Hybrid Systems: Computation and Control, LNCS, vol. 3927, pp. 257–271. Springer, Berlin (2006)
Gonnord, L., Halbwachs, N.: Combining widening and acceleration in linear relation analysis. In: Static Analysis Symposium, LNCS, pp. 144–160. Springer, Berlin (2006)
Gonnord, L., Schrammel, P.: Abstract acceleration in linear relation analysis. Sci. Comput. Program. 93(Part B), 125–153 (2014)
Gulwani, S., Tiwari, A.: Constraint-based approach for analysis of hybrid systems. In: Computer-Aided Verification, LNCS, vol. 5123, pp. 190–203. Springer, Berlin (2008)
Halbwachs, N., Raymond, P., Proy, Y.E.: Verification of linear hybrid systems by means of convex approximations. In: Static Analysis Symposium, LNCS, vol. 864, pp. 223–237. Springer, Berlin (1994)
Henzinger, T.A., Ho, P.H., Wong-Toi, H.: HyTech: a model checker for hybrid systems. J. Softw. Tools Technol. Transf. 1(1–2), 110–122 (1997)
Horn, R.A., Johnson, C.R.: Matrix Analysis. Cambridge University Press, Cambridge (2012)
Jeannet, B.: Interproc analyzer for recursive programs with numerical variables (2010). http://pop-art.inrialpes.fr/interproc/interprocweb.cgi
Jeannet, B., Schrammel, P., Sankaranarayanan, S.: Abstract acceleration of general linear loops. In: Principles of Programming Languages, pp. 529–540. ACM, New York (2014)
Johnson, T.T., Mitra, S.: Passel: A verification tool for parameterized networks of hybrid automata (2012). https://publish.illinois.edu/passel-tool/
Kannan, R., Lipton, R.J.: Polynomial-time algorithm for the orbit problem. J. ACM (JACM) 33(4), 808–821 (1986)
Knospe, C.: PID control. IEEE Control Syst. 26(1), 30–31 (2006)
Lancaster, P., Tismenetsky, M.: The Theory of Matrices, 2nd edn. Academic Press, Cambridge (1984)
Le Guernic, C.: Toward a sound analysis of guarded LTI loops with inputs by abstract acceleration. In: International Static Analysis Symposium, pp. 192–211. Springer, Berlin (2017)
Le Guernic, C.: Toward a sound analysis of guarded LTI loops with inputs by abstract acceleration (extended version) Working paper or preprint (2017). https://hal.inria.fr/hal-01550767
Le Guernic, C., Girard, A.: Reachability analysis of hybrid systems using support functions. In: CAV, LNCS, vol. 5643, pp. 540–554. Springer, Berlin (2009)
van Loan, C.F., Golub, G.: Matrix Computations. Johns Hopkins Studies in Mathematical Sciences. JHU Press, Baltimore (1996)
Löhner, R.: Einschließung der lösung gewöhnlicher anfangs- und randwertaufgaben und anwendungen. Ph.D. thesis, Universität Karlsruhe (1988)
Ouaknine, J., Worrell, J.: Positivity problems for low-order linear recurrence sequences. In: Proceedings of the twenty-fifth annual ACM-SIAM Symposium on Discrete Algorithms, pp. 366–379. Society for Industrial and Applied Mathematics (2014)
Penrose, R.: A generalized inverse for matrices. In: Mathematical Proceedings of the Cambridge Philosophical Society, vol. 51–3, pp. 406–413. Cambridge University Press, Cambridge (1955)
Sankaranarayanan, S., Tiwari, A.: Relational abstractions for continuous and hybrid systems. In: Computer-Aided Verification, LNCS, vol. 6806, pp. 686–702. Springer, Berlin (2011)
Schrammel, P.: Unbounded-time reachability analysis of hybrid systems by abstract acceleration. In: Embedded Software, pp. 51–54. IEEE (2015)
Schrammel, P., Jeannet, B.: Applying abstract acceleration to (co-)reachability analysis of reactive programs. J. Symbol. Comput. 47(12), 1512–1532 (2012)
Acknowledgements
University of Oxford; We would like to thank Colas Le Guernic for his constructive suggestions and comments on the paper. This work is supported by the Alan Turing Institute, London, UK, by EPSRC Grant EP/J012564/1, ERC Project 280053 (CPROVER), by the H2020 FET OPEN 712689 \(\hbox {SC}^2\) and by Oxford Instruments PLC.
Author information
Authors and Affiliations
Corresponding author
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Cattaruzza, D., Abate, A., Schrammel, P. et al. Unbounded-Time Safety Verification of Guarded LTI Models with Inputs by Abstract Acceleration. J Autom Reasoning 65, 157–203 (2021). https://doi.org/10.1007/s10817-020-09562-z
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s10817-020-09562-z