Reachability analysis of linear hybrid systems via block decomposition

TL;DR

This paper introduces a novel approach for reachability analysis of linear hybrid systems that reduces computational complexity by performing most calculations in low-dimensional spaces through block decomposition and strategic operator enhancements.

Contribution

The paper presents an improved method for reachability analysis that optimizes continuous and discrete post operators, enabling efficient analysis of high-dimensional linear hybrid systems.

Findings

Enhanced continuous-post operator reduces computation in high-dimensional space.

Discrete-post operator leverages system structure for low-dimensional computations.

Demonstrated effectiveness on challenging benchmark systems.

Abstract

Reachability analysis aims at identifying states reachable by a system within a given time horizon. This task is known to be computationally expensive for linear hybrid systems. Reachability analysis works by iteratively applying continuous and discrete post operators to compute states reachable according to continuous and discrete dynamics, respectively. In this paper, we enhance both of these operators and make sure that most of the involved computations are performed in low-dimensional state space. In particular, we improve the continuous-post operator by performing computations in high-dimensional state space only for time intervals relevant for the subsequent application of the discrete-post operator. Furthermore, the new discrete-post operator performs low-dimensional computations by leveraging the structure of the guard and assignment of a considered transition. We illustrate the potential of our approach on a number of challenging benchmarks.

Figures

Tables

Table 1. TABLE I : For each benchmark we report the number of dimension and the number of constrained dimensions ( Dim. (constr.) ), the step size of the time discretization ( Step ) and the run time of the different algorithms as described in Section  V-B (e.g., “Deco” refers to the decomposed algorithm presented in this paper) in seconds (where the fastest solution is marked in bold face). For benchmark instances with parenthesized computation time, the safety property could not be proven because the overapproximation was too coarse, in which case the computation terminated as soon as the violation was detected. “TO” and “OOM” indicate a timeout of 5 × 10 4 5 superscript 10 4 5\times 10^{4} seconds and an out-of-memory error, respectively.

Benchmark	Dim. (constr.)	Step	Deco	LazySupp	LazyOptim	SpaceEx LGG	SpaceEx STC
linear_switching	5 (1)	0.0001	$2.50\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟎}}$	$1.27\times {10}^1$	$2.81\times {10}^1$	$2.60\times {10}^1$	$2.30\times {10}^1$
spacecraft_noabort	5 (4)	0.04	$5.30\times {10}^0$	$3.42\times {10}^0$	$2.19\times {10}^2$	$1.18\times {10}^0$	$3.50\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏}}$
spacecraft_120	5 (5)	0.04	$5.30\times {10}^0$	$2.10\times {10}^0$	$4.30\times {10}^1$	$1.91\times {10}^0$	$8.10\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏}}$
platoon_bounded	10 (4)	0.01	$1.30\times {\mathrm{𝟏𝟎}}^{-\mathrm{𝟏}}$	$1.60\times {10}^{-1}$	$5.69\times {10}^0$	$5.55\times {10}^0$	$1.60\times {10}^0$
platoon_unbounded	10 (4)	0.03	$1.08\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟎}}$	$1.16\times {10}^0$	$4.96\times {10}^1$	$3.46\times {10}^1$	$6.50\times {10}^1$
filtered_osc64	67 (3)	0.01	$2.81\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟎}}$	($7.43\times {10}^0$)†	$5.63\times {10}^2$	$2.04\times {10}^1$	$3.25\times {10}^1$
filtered_osc128	131 (3)	0.01	$7.95\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟎}}$	($4.29\times {10}^1$)†	$1.79\times {10}^3$	$1.69\times {10}^2$	$4.67\times {10}^3$
filtered_osc256	259 (3)	0.01	$2.80\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟏}}$	($9.19\times {10}^1$)†	$9.99\times {10}^4$	$8.70\times {10}^3$	OOM
filtered_osc512	515 (3)	0.01	$1.13\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟐}}$	($4.73\times {10}^2$)†	TO	TO	TO
filtered_osc1024	1027 (3)	0.01	$5.09\times {\mathrm{𝟏𝟎}}^{\mathrm{𝟐}}$	($5.11\times {10}^3$)†	TO	TO	TO

Table 2. TABLE II : Evaluation for different time steps on the filtered oscillator model with 64 filters (“filtered_osc64” in Table  I ) and all algorithms (see Section  V-B ) that allow varying the time step. The last row shows the relative change of the fourth row over the first row as base line.

Step	Deco	LazySupp	LazyOptim	SpaceEx LGG
0.01	$2.8\times {10}^0$	$7.4\times {10}^0$	$5.6\times {10}^2$	$2.0\times {10}^1$
0.005	$4.1\times {10}^0$	$1.3\times {10}^1$	$9.5\times {10}^2$	$3.6\times {10}^1$
0.001	$1.7\times {10}^1$	$6.5\times {10}^1$	$4.7\times {10}^3$	$1.6\times {10}^2$
0.0005	$3.2\times {10}^1$	$1.3\times {10}^2$	$8.4\times {10}^3$	$3.2\times {10}^2$
$\times 20.0$	$\times 11.4$	$\times 17.5$	$\times 15.0$	$\times 16.0$

Equations

$$d_{H}^{p}\bigl{(}\mathcal{X},\mathcal{Y}\bigr{)}:=\inf_{\varepsilon\in\mathbb{R}}\left\{\mathcal{Y}\subseteq\mathcal{X}\oplus\varepsilon\mathcal{B}_{p}^{n}\text{ and }\mathcal{X}\subseteq\mathcal{Y}\oplus\varepsilon\mathcal{B}_{p}^{n}\right\}.$$

$${\rho }_{\mathcal{X}}(d):=\underset{x\in \mathcal{X}}{\max }⟨d,x⟩.$$

$$d_H^p(\mathcal{X},\mathcal{Y})=\underset{\Vert d{\Vert }_p\le 1}{\max }{\rho }_{\mathcal{Y}}(d)-{\rho }_{\mathcal{X}}(d).$$

$$\dot{x}(t)=Ax(t)+Bu(t),u(t)\in \mathcal{U}.$$

$${\xi }_{x_0,u}(t)=e^{At}{x}_0+{\int }_0^t{e}^{A(t-s)}Bu(s)ds.$$

$$Pos{t}_C((A,B,\mathcal{U}),{\mathcal{X}}_0):=$$

$$\{{\xi }_{x_0,u}(t)\mid t\ge 0,x_0\in {\mathcal{X}}_0,u(s)\in \mathcal{U}\text{ for all }s\}.$$

$${Post}_D(\ell ,\mathcal{X}):=\bigcup _{{\ell }^{\prime }\in Loc}\{$$

$$(\mathcal{X}\cap Inv(\ell )\cap Grd(\ell ,{\ell }^{\prime }))\cap Inv({\ell }^{\prime }))\}$$

$${\mathcal{R}}_0\cup \bigcup _{(\ell ,\mathcal{X})\in \mathcal{R}}{Post}_D(\ell ,Pos{t}_C(Flow(\ell ),\mathcal{X}))\subseteq \mathcal{R}.$$

$$\mathcal{X}(k):=\Phi \mathcal{X}(k-1)\oplus \mathcal{V}={\Phi }^k\mathcal{X}(0)\oplus {⨁}_{j=0}^{k-1}{\Phi }^j\mathcal{V}.$$

$${\mathcal{X}}_i(k):={⨁}_{j=1}^{\beta }({\Phi }^k{)}_{ij}{\mathcal{X}}_j(0)\oplus {⨁}_{j=0}^{k-1}[({\Phi }^j{)}_{i1}\cdots ({\Phi }^j{)}_{i\beta }]\mathcal{V}.$$

$$(A,b)\odot (\hat{\mathcal{X}}\cap {\mathcal{G}}^{\ast })$$

$$\hat{\mathcal{X}}\cap \hat{\mathcal{G}}=({\mathcal{X}}_1\times {\mathcal{X}}_2)\cap ({\mathcal{G}}_1\times {\mathcal{G}}_2)=({\mathcal{X}}_1\cap {\mathcal{G}}_1)\times ({\mathcal{X}}_2\cap {\mathcal{G}}_2)$$

$$\hat{\mathcal{X}}\cap \hat{\mathcal{G}}=\varnothing ⟺({\mathcal{X}}_1\cap {\mathcal{G}}_1=\varnothing )\vee ({\mathcal{X}}_2\cap {\mathcal{G}}_2=\varnothing ).$$

$$\hat{\mathcal{X}}\cap \hat{\mathcal{G}}=\varnothing ⟺{\mathcal{X}}_1\cap {\mathcal{G}}_1=\varnothing ,$$

$$\hat{\mathcal{X}}\cap \mathcal{G}\subseteq ({\mathcal{X}}_1\cap {\pi }_1\mathcal{G})\times ({\mathcal{X}}_2\cap {\pi }_2\mathcal{G})$$

$$\begin{array}{rl}\hat{\mathcal{X}}\cap \mathcal{G}=\varnothing & ⟸({\mathcal{X}}_1\cap {\pi }_1\mathcal{G}=\varnothing )\vee ({\mathcal{X}}_2\cap {\pi }_2\mathcal{G}=\varnothing )\\ & ⟸{\mathcal{X}}_1\cap {\pi }_1\mathcal{G}=\varnothing .\end{array}$$

$$d_H^p(\hat{\mathcal{X}}\cap \mathcal{G},\hat{\mathcal{X}}\cap \hat{\mathcal{G}})\le \underset{j}{\max }\min (\Vert {\Delta }_p({\mathcal{X}}_j){\Vert }_p,\Vert {\Delta }_p({\pi }_j\mathcal{G}){\Vert }_p).$$

$$A\hat{\mathcal{X}}\oplus \{b\}\subseteq {\text{bigtimes}}_i{⨁}_j{A}_{ij}{\mathcal{X}}_j\oplus \{b_i\}$$

$$d_H^p(A\mathcal{X},{\text{bigtimes}}_i{⨁}_j{A}_{ij}{\mathcal{X}}_j)$$

$$=$$

$$\le$$

$$\left(\begin{array}{cc}{A}_{11}& 0\\ 0& A_{22}\end{array}\right){\mathcal{X}}_1\times {\mathcal{X}}_2\oplus \{b_1\}\times \{b_2\}$$

$$=(A_{11}{\mathcal{X}}_1\oplus \{b_1\})\times (A_{22}{\mathcal{X}}_2\oplus \{b_2\}).$$

$$\hat{\mathcal{X}}\subseteq \hat{\mathcal{G}}⟺{\bigwedge }_j{\mathcal{X}}_j\subseteq {\mathcal{G}}_j.$$

$$\mathrm{CH}(\hat{\mathcal{X}}\cup \hat{\mathcal{G}})\subseteq {\text{bigtimes}}_j\mathrm{CH}({\mathcal{X}}_j\cup {\mathcal{G}}_j).$$

$$d_H^p(\mathrm{CH}(\hat{\mathcal{X}}\cup \hat{\mathcal{G}}),{\text{bigtimes}}_j\mathrm{CH}({\mathcal{X}}_j\cup {\mathcal{G}}_j))$$

$$\le$$

$$d_H^p(\mathcal{X},{\text{bigtimes}}_j{\pi }_j\mathcal{X})$$

$$=$$

$${\text{bigtimes}}_j{\pi }_j\mathcal{X}\subseteq \{c_{\mathcal{X}}^p\}\oplus \Vert r_{\mathcal{X}}^p{\Vert }_p{\mathcal{B}}_p^n\subseteq \mathcal{X}\oplus \Vert r_{\mathcal{X}}^p{\Vert }_p{\mathcal{B}}_p^n$$

$$d_H^p(\hat{\mathcal{X}}\cap \mathcal{G},\hat{\mathcal{X}}\cap \hat{\mathcal{G}})$$

$$=$$

Full text

I Introduction

A hybrid system [1] is a formalism for modeling cyber-physical systems. Reachability analysis is a rigorous way to reason about the behavior of hybrid systems.

In this paper, we describe a new reachability algorithm for linear hybrid systems, i.e., hybrid systems with dynamics given by linear differential equations and invariants and guards given by linear inequalities. The key feature of our algorithm is that it performs computations in low-dimensional subspaces, which greatly improves scalability. To this end, we integrate (and modify) a recent reachability algorithm for (purely continuous) LTI systems [2], which we call $\mathit{Post}_{C}^{\square}$ in the following, in a new algorithm for linear hybrid systems.

The $\mathit{Post}_{C}^{\square}$ algorithm decomposes the calculation of the reachable states into calculations in subspaces (called blocks). This decomposition has two benefits. The first benefit is that computations in lower dimensions are generally more efficient and thus the algorithm is highly scalable. The second benefit is that the analysis for different subspaces is decoupled; hence one can effectively skip the computations for dimensions that are of no interest (e.g., for a safety property).

Conceptually, reachability analysis for hybrid systems is performed in a “hybrid loop” that interleaves a continuous-post algorithm and a discrete-post algorithm. If we consider $\mathit{Post}_{C}^{\square}$ as a black box, we can plug it into this hybrid loop, which we refer to as $\mathit{Post}_{H}$ (cf. Section III). However, there are two shortcomings of such a naive integration. First, all operations besides $\mathit{Post}_{C}^{\square}$ are still performed in high dimensions, making $\mathit{Post}_{H}$ still suffer from scalability issues. Second, $\mathit{Post}_{H}$ does not utilize the decoupling of $\mathit{Post}_{C}^{\square}$ at all.

We demonstrate that, unlike in $\mathit{Post}_{H}$, it is possible to perform all computations in low dimensions (cf. Section IV). For this purpose, we modify $\mathit{Post}_{C}^{\square}$ as well as $\mathit{Post}_{H}$. In common cases, our algorithm does not introduce an additional approximation error. Furthermore, our algorithm makes proper use of the second benefit of $\mathit{Post}_{C}^{\square}$ by computing the reachable states only in specific dimensions whenever possible.

We implemented the algorithm in JuliaReach, a toolbox for reachability analysis [3, 4], and we evaluate the algorithm on several benchmark problems, including a 1024-dimensional hybrid system (cf. Section V). Our algorithm outperforms the naive $\mathit{Post}_{H}$ and other state-of-the-art algorithms by several orders of magnitude.

To summarize, we show how to modify the decomposed reachability algorithm for LTI systems from [2] in order to efficiently integrate it into a new, decomposed reachability algorithm for linear hybrid systems. The key insights are (1) to exploit the decomposed structure of the algorithm to perform all operations in low dimensions and (2) to only compute the reachable states in specific dimensions whenever possible.

Related work

Reachability analysis for linear hybrid systems

There are several tools available for the analysis of the class of systems we consider in this paper. All participating tools in a recent friendly competition [5] follow the architecture of two post operators, one for continuous time and one for discrete time [6, 7, 8, 3, 9, 10]. We implemented our algorithm in JuliaReach [3] and will consider the mature tool SpaceEx [9] in the evaluation in Section V.

Decomposition

Hybrid systems given as a network of components can be explored efficiently in a symbolic way, e.g., using bounded model checking [11]. We consider a decomposition in the continuous state space. Schupp et al. perform such a decomposition by syntactic independence [12], which corresponds to dynamics with block-diagonal matrices (whereas our decomposition is generally applicable).

For purely continuous systems there exist various decomposition approaches. In this work we build on [2] for LTI systems, which decomposes the system into blocks and exploits the linear dynamics to avoid the wrapping effect. Other approaches for LTI systems are based on Krylov subspace approximations [13], time-scale decomposition [14, 15], similarity transformations [16, 17], projectahedra [18, 19], and sub-polyhedra abstract domains [20]. Approaches for nonlinear systems are based on projections with differential inclusions [21], Hamilton-Jacobi methods [22, 23], and hybridization with iterative refinement [24].

Lazy flowpipe computation

The support-function representation of convex sets can be used to represent a flowpipe (a sequence of sets that covers the behaviors of a system) symbolically [25]. Only sets that are of interest, e.g., those that intersect with a linear constraint, need then be approximated [26]. Using our decomposition approach, we can even avoid the symbolic computation in dimensions that are irrelevant to the intersection. Our approach is independent of the set representation, so it can also be applied in analyses based on, e.g., zonotopes [27, 28, 29, 30]. Given a linear switching system with a hyperplanar state-space partition, one approach computes ellipsoidal approximations of the reach set on the partition borders, without computing the full reach set [31].

Intersection of convex sets

Performing intersections in low dimensions allows for efficient computations that are not possible in high dimensions. For example, checking for emptiness of a polyhedron in constraint representation is a feasibility linear program, which can be solved in weakly polynomial time, but solutions in strongly polynomial time are only known in two dimensions [32].

In the context of hybrid-system reachability, computing intersections is a major challenge because it usually requires a conversion from efficient set representations (like zonotopes, support function, or Taylor models) to polytopes and back, often entailing additional approximation. Below we summarize how other approaches tackle the intersection problem.

A coarse approximation of an intersection of a set $\mathcal{X}$ and another set $\mathcal{Y}$ can be obtained by only detecting a nonempty intersection (which is generally easier to do) and then taking the original set $\mathcal{X}$ as overapproximation [33]. In general, the intersection between a polytope and polyhedral constraints (invariants and guards) can be computed exactly, but such an approach is not scalable [34]. Girard and Le Guernic consider hyperplanar constraints where reachable states are either zonotopes, in which case they work in a two-dimensional projection [28], or general polytopes [35, 36]. The tool SpaceEx approximates the intersection of polytopes and general polyhedra using template directions [9]. Frehse and Ray propose an optimization scheme for the intersection of a compact convex set $\mathcal{X}$, represented by its support function, and a polyhedron $\mathcal{Y}$, and this scheme is exact if $\mathcal{X}$ is a polytope [26]. The problem of performing intersections can also be cast in terms of finding separating hyperplanes [37, 38]. Althoff et al. approximate zonotopes by parallelotopes before considering the intersection [29]. For must semantics, Althoff and Krogh use constant-dynamics approximation and obtain a nonlinear mapping [39]. Under certain conditions, Bak et al. apply a model transformation by replacing guard constraints by time-triggered constraints, for which intersection is easy [40].

II Preliminaries

We introduce some notation. The real numbers are denoted by $\mathbb{R}$. The origin in $\mathbb{R}^{n}$ is written $\mathbf{0}_{n}$. Given two vectors $x,y\in\mathbb{R}^{n}$, their dot product is $\langle x,y\rangle:=\sum_{i=1}^{n}x_{i}\cdot y_{i}$. For $p\geq 1$, the $p$-norm of a matrix $A\in\mathbb{R}^{n\times n}$ is denoted $\|A\|_{p}$. The diameter of a set $\mathcal{X}\subseteq\mathbb{R}^{n}$ is $\Delta_{p}(\mathcal{X}):=\sup_{x,y\in\mathcal{X}}\|x-y\|_{p}$. The $n$-dimensional unit ball of the $p$-norm is $\mathcal{B}_{p}^{n}:=\{x\in\mathbb{R}^{n}\mid 1\geq\|x\|_{p}\}$. An $n$-dimensional half-space is the set $\{\langle a,x\rangle\leq b\mid x\in\mathbb{R}^{n}\}$ parameterized by $a\in\mathbb{R}^{n}$, $b\in\mathbb{R}$. A polyhedron is an intersection of finitely many half-spaces, and a polytope is a bounded polyhedron.

Given sets $\mathcal{X}\subseteq\mathbb{R}^{n}$ and $\mathcal{Y}\subseteq\mathbb{R}^{m}$, scalar $\lambda\in\mathbb{R}$, matrix $A\in\mathbb{R}^{n\times n}$, and vector $b\in\mathbb{R}^{n}$, we use the following operations on sets: scaling $\lambda\mathcal{X}:=\{\lambda x\mid x\in\mathcal{X}\}$, linear map $A\mathcal{X}:=\{Ax\mid x\in\mathcal{X}\}$, Minkowski sum $\mathcal{X}\oplus\mathcal{Y}:=\{x+y\mid x\in\mathcal{X}\text{ and }y\in\mathcal{Y}\}$ (if $n=m$), affine map $(A,b)\odot\mathcal{X}:=A\mathcal{X}\oplus\{b\}$, Cartesian product $\mathcal{X}\times\mathcal{Y}:=\{(x,y)\mid x\in\mathcal{X},y\in\mathcal{Y}\}$, intersection $\mathcal{X}\cap\mathcal{Y}:=\{z\mid z\in\mathcal{X},z\in\mathcal{Y}\}$ (if $n=m$), and convex hull $\mathrm{CH}(\mathcal{X}):=\{\lambda\cdot x+(1-\lambda)\cdot y\mid x,y\in\mathcal{X},0\leq\lambda\leq 1\}$.

Given two sets $\mathcal{X},\mathcal{Y}\subseteq\mathbb{R}^{n}$, the Hausdorff distance is

[TABLE]

Let $\mathcal{C}_{n}\subseteq 2^{\mathbb{R}^{n}}$ be the set of $n$-dimensional compact convex sets. For a nonempty set $\mathcal{X}\in\mathcal{C}_{n}$, the support function $\rho_{\mathcal{X}}:\mathbb{R}^{n}\to\mathbb{R}$ is defined as

[TABLE]

The Hausdorff distance of two sets $\mathcal{X},\mathcal{Y}\in\mathcal{C}_{n}$ with $\mathcal{X}\subseteq\mathcal{Y}$ can also be expressed in terms of the support function as

[TABLE]

We use $\beta$ to denote the number of blocks in a partition. Let $\{\pi_{j}\}_{j=1}^{\beta}$ be a set of (contiguous) projection matrices that partition a vector $x\in\mathbb{R}^{n}$ into $x=[\pi_{1}x,\ldots,\pi_{\beta}x]$. Given a set $\mathcal{X}$ and projection matrices $\{\pi_{j}\}_{j=1}^{\beta}$, we call $\pi_{j}\mathcal{X}$ a block of $\mathcal{X}$ and $\bigtimes_{j}\pi_{j}\mathcal{X}=\pi_{1}\mathcal{X}\times\cdots\times\pi_{\beta}\mathcal{X}$ the Cartesian decomposition with the block structure induced by projections $\pi_{j}$. We write $\widehat{\mathcal{X}}$ to indicate a decomposed set (i.e., a Cartesian product of lower-dimensional sets). For instance, given a nonempty set $\mathcal{X}\in\mathcal{C}_{n}$, its box approximation is the Cartesian decomposition into intervals (i.e., one-dimensional blocks). We can bound the approximation error by the radius of $\mathcal{X}$.

Proposition 1

Let $\mathcal{X}\in\mathcal{C}_{n}$ be nonempty, $p=\infty$, $r_{\mathcal{X}}^{p}$ be the radius of the box approximation of $\mathcal{X}$, and let $\pi_{j}$ be appropriate projection matrices. Then $d_{H}^{p}(\mathcal{X},\bigtimes_{j}\pi_{j}\mathcal{X})\leq\|r_{\mathcal{X}}^{p}\|_{p}$.

II-A LTI systems

An $n$-dimensional LTI system $(A,B,\mathcal{U})$, with matrices $A\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{n\times m}$, and input domain $\mathcal{U}\in\mathcal{C}_{m}$, is a system of ODEs of the form

[TABLE]

We denote the set of all $n$-dimensional LTI systems by $\mathcal{L}_{n}$. From now on, a vector $x\in\mathbb{R}^{n}$ is also called a (continuous) state. Given an initial state $x_{0}\in\mathbb{R}^{n}$ and an input signal $u$ such that $u(t)\in\mathcal{U}$ for all $t$, a trajectory of (1) is the unique solution $\mathit{\xi}_{x_{0},u}:\mathbb{R}_{\geq 0}\rightarrow\mathbb{R}^{n}$ with

[TABLE]

Given an LTI system $(A,B,\mathcal{U})$ and a set $\mathcal{X}_{0}\in\mathcal{C}_{n}$ of initial states, the continuous-post operator, $\mathit{Post_{C}}$, computes the set of reachable states for all input signals $u$ over $\mathcal{U}$:

[TABLE]

II-B Linear hybrid systems

We briefly introduce the syntax of linear hybrid systems used in this work and refer to the literature for the semantics [41, 1]. An $n$-dimensional linear hybrid system is a tuple $\mathcal{H}=(\mathit{Var},Loc,\mathit{Flow},\mathit{Inv},\mathit{Grd},\mathit{Asgn})$ with variables $\mathit{Var}=\{x_{1},\dots,x_{n}\}$, a finite set of locations $Loc$, two functions $\mathit{Flow}:Loc\to\mathcal{L}_{n}$ and $\mathit{Inv}:Loc\to\mathcal{C}_{n}$ that respectively assign continuous dynamics and an invariant to each location, and two functions $\mathit{Grd}:Loc\times Loc\to\mathcal{C}_{n}$ and $\mathit{Asgn}:Loc\times Loc\to\mathbb{R}^{n\times n}\times\mathbb{R}^{n}$ that respectively assign a guard and an assignment in the form of a deterministic affine map to each pair of locations. If $\mathit{Grd}(\ell,\ell^{\prime})\neq\emptyset$, we call $(\ell,\ell^{\prime})$ a (discrete) transition.

Let $\mathcal{H}=(\mathit{Var},Loc,\mathit{Flow},\mathit{Inv},\mathit{Grd},\mathit{Asgn})$ be a linear hybrid system. A (symbolic) state of $\mathcal{H}$ is a pair $(\ell,\mathcal{X})\in\mathit{Loc\times 2^{\mathbb{R}^{n}}}$. The discrete-post operator, $\mathit{Post}_{D}$, maps a symbolic state to a set of symbolic states by means of discrete transitions:

[TABLE]

The reach set of $\mathcal{H}$ from a set of initial symbolic states $\mathcal{R}_{0}$ of $\mathcal{H}$ is the smallest set $\mathcal{R}$ of symbolic states such that

[TABLE]

We note that the usual definition of the reach set of a hybrid system is stricter wrt. invariants. However, practical analysis tools use the above definition for efficiency reasons.

Safety properties can be given as a set of symbolic error states that should be avoided and be encoded as the guard of a transition to a new error location. In our implementation we can also check inclusion in the safe states (the complement) and do not require an encoding with additional transitions.

III Reachability analysis of linear hybrid systems

Our reachability algorithm for linear hybrid systems is centered around the algorithm from [2], called $\mathit{Post}_{C}^{\square}$ for convenience, which implements $\mathit{Post_{C}}$ for LTI systems in a compositional way. In this section, we first recall the $\mathit{Post}_{C}^{\square}$ algorithm. Two important properties of $\mathit{Post}_{C}^{\square}$ are that (1) the output is a sequence of decomposed sets and that (2) this sequence is computed efficiently in low dimensions.

After explaining the algorithm $\mathit{Post}_{C}^{\square}$, we incorporate it in a standard reachability algorithm for linear hybrid systems. However, this standard reachability algorithm will not make use of the above-mentioned properties. This motivates our new algorithm (presented in Section IV), which is a modification of both this standard reachability algorithm and $\mathit{Post}_{C}^{\square}$ to make optimal use of these properties.

III-A Decomposed reachability analysis of LTI systems

The decomposition-based approach [2] follows a flowpipe-construction scheme using time discretization, which we shortly recall here. Given an LTI system $(A,B,\mathcal{U})$ and a set of initial (continuous) states $\mathcal{X}_{0}$, by fixing a time step $\delta$ we first compute a set $\mathcal{X}(0)$ that overapproximates the reach set up to time $\delta$, a matrix $\Phi=e^{A\delta}$ that captures the dynamics of duration $\delta$, and a set $\mathcal{V}$ that overapproximates the effect of the inputs up to time $\delta$. We obtain an overapproximation of the reach set in time interval $[k\delta,(k+1)\delta]$, for step $k>0$, with

[TABLE]

Algorithm $\mathit{Post}_{C}^{\square}$ decomposes this scheme. Fixing some block structure, let $\widehat{\mathcal{X}}(0):=\mathcal{X}_{1}(0)\times\dots\times\mathcal{X}_{\beta}(0)$ be the corresponding Cartesian decomposition of $\mathcal{X}(0)$. We compute a sequence $\widehat{\mathcal{X}}(k):=\mathcal{X}_{1}(k)\times\dots\times\mathcal{X}_{\beta}(k)$ such that $\mathcal{X}(k)\subseteq\widehat{\mathcal{X}}(k)$ for every $k$. Each low-dimensional set $\mathcal{X}_{i}(k)$ is computed as

[TABLE]

The above sequences $\mathcal{X}(k)$ resp. $\widehat{\mathcal{X}}(k)$ are called flowpipes.

Example

We illustrate the algorithms with a running example throughout the paper. For illustration purposes, the example is two-dimensional (and hence we decompose into one-dimensional blocks, but we emphasize that the approach also generalizes to higher-dimensional decomposition) and we consider a hybrid system with only a single location and one transition (a self-loop). Figure 1 depicts the flowpipe construction for the example.

III-B Reachability analysis of linear hybrid systems

We now discuss a standard reachability algorithm for hybrid systems. Essentially, this algorithm interleaves the operators $\mathit{Post_{C}}$ and $\mathit{Post}_{D}$ following (3) until it finds a fixpoint. Here we use $\mathit{Post}_{C}^{\square}$ as the continuous-post operator.

We first compute a flowpipe $\widehat{\mathcal{X}}=\widehat{\mathcal{X}}(0),\dots,\widehat{\mathcal{X}}(N)$ using $\mathit{Post}_{C}^{\square}$ as described above. Then we use $\mathit{Post}_{D}$ to take a discrete transition. According to (2), we want to compute $((A,b)\odot(\widehat{\mathcal{X}}\cap\mathcal{I}_{1}\cap\mathcal{G}))\cap\mathcal{I}_{2}$, where $(A,b)$ is a deterministic affine map $\mathcal{I}_{1}$ and $\mathcal{I}_{2}$ are the source and target invariant, and $\mathcal{G}$ is the guard. Frehse and Ray showed that for such maps the term can be simplified to

[TABLE]

where the set $\mathcal{G}^{*}$ can be statically precomputed [26], which is usually easy because the sets $\mathcal{I}_{1}$, $\mathcal{G}$, and $\mathcal{I}_{2}$ are given as polyhedra in constraint representation (also called H-representation, as opposed to the (vertex) V-representation). Hence we ignore invariants in the rest of the presentation for simplicity.

Example

We continue with the flowpipe from Figure 2. The guard $\mathcal{G}$ is a half-space that is constrained in dimension $x_{1}$ and unconstrained in dimension $x_{2}$. Only the last two sets in the flowpipe intersect with the guard. The assignment here is a translation in dimension $x_{1}$. The resulting intersection, before and after the translation, is depicted in Figure 3.

Finally, the algorithm checks for a fixpoint, i.e., for inclusion of the symbolic states we computed with $\mathit{Post}_{D}$ in previously-seen symbolic states.

Example

The green set in Figure 3 shows that the two sets we obtained from $\mathit{Post}_{D}$ are contained in $\widehat{\mathcal{X}}(0)$ computed before. (Recall that in this example we only consider a single location; hence the inclusion of the sets implies inclusion of the symbolic states.)

The steps outlined above describe one iteration of the standard reachability algorithm. Each symbolic state for which the fixpoint check was negative spawns a new flowpipe. Since this can lead to a combinatorial explosion, one typically applies a technique called clustering (cf. [9]), where symbolic states are merged after the application of $\mathit{Post}_{D}$. Here we consider clustering with a convex hull.

Example

Assume that the fixpoint check above was negative for both sets that we tested. In Figure 4 we show the convex hull of the sets in purple.

Up to now, we have seen a standard incorporation of an algorithm for the continuous-post operator $\mathit{Post_{C}}$ (for which we used $\mathit{Post}_{C}^{\square}$) into a reachability algorithm for hybrid systems. Observe that $\mathit{Post_{C}}$ was used as a black box. Consequently, we could not make use of the properties of the specific algorithm $\mathit{Post}_{C}^{\square}$. In particular, apart from $\mathit{Post}_{C}^{\square}$, we performed all computations in high dimensions. In the next section, we describe a new algorithm that instead performs all computations in low dimensions.

IV Decomposed reachability analysis

We now present a new, decomposed reachability algorithm for linear hybrid systems. The algorithm uses a modified version of $\mathit{Post}_{C}^{\square}$ for computing flowpipes and has two major performance improvements over the algorithm seen before.

Recall that $\mathit{Post}_{C}^{\square}$ computes flowpipes consisting of decomposed sets. The first improvement is to exploit the decomposed structure to perform all other operations (intersection, affine map, inclusion check, and convex hull) in low dimensions.

The second improvement is to compute flowpipes in a sparse way. Roughly speaking, we are only interested in those dimensions of a flowpipe that are relevant to determine an intersection with a guard. We only need to compute the other dimensions of the flowpipe if we detect such an intersection.

We defer proofs in this section to Appendix References.

The algorithm starts as before: given an initial (symbolic) state, we compute $\mathcal{X}(0)$ (discretization) and decompose the set to obtain $\widehat{\mathcal{X}}(0)$. Next we want to compute a flowpipe, and this is where we deviate from the previous algorithm.

IV-A Computing a sparse flowpipe

We modify $\mathit{Post}_{C}^{\square}$ in order to control the dimensions of the flowpipe. Recall that the black-box version of $\mathit{Post}_{C}^{\square}$ computed the flowpipe $\widehat{\mathcal{X}}(k)=\mathcal{X}_{1}(k)\times\dots\times\mathcal{X}_{\beta}(k)$ for $k=1,\dots,N$, i.e., in all dimensions. This is usually not necessary for detecting an intersection with a guard. We will discuss this formally below, but want to establish some intuition first. Recall the running example from before. The guard was only constrained in dimension $x_{1}$. This means that the bounds of the sets $\widehat{\mathcal{X}}(k)=\mathcal{X}_{1}(k)\times\mathcal{X}_{2}(k)$ in dimension $x_{2}$ are irrelevant. Consequently, we do not need to compute the sets $\mathcal{X}_{2}(k)$ at all (at least for the moment). We only compute those dimensions of a flowpipe that are necessary to determine intersection with the guards. Identifying these dimensions and projecting the guards accordingly has to be performed only once per transition and is often just a syntactic procedure.

In particular, Algorithm 1 describes the new implementation of $\mathit{Post}_{C}^{\square}$ to compute a sparse flowpipe. The algorithm starts the same way as the original algorithm in [2] with a decomposed set in line 1 and iteratively computes a sparse flowpipe only for the dimensions of interest (line 1). However, if we detect an intersection with a guard constraint (line 1), we compute the full-dimensional flowpipe for the corresponding time frame. The computation of the inputs $\mathcal{\hat{V}_{\text{tmp}}}$ remains the same as in the original algorithm.

Example

As discussed, we only compute the flowpipe $\mathcal{X}_{1}(1),\dots,\mathcal{X}_{1}(4)$ for the first block (in dimension $x_{1}$), i.e., a sequence of intervals, which is depicted in Figure 5(a). Projecting the guard to $x_{1}$, we obtain a ray $\mathcal{G}_{1}$. As expected, we observe an intersection with the guard for the same time steps as before (namely steps $k=3$ and $k=4$).

IV-B Decomposing an intersection

Computing the intersection of two $n$-dimensional sets $\mathcal{X}$ and $\mathcal{G}$ in low dimensions is generally beneficial for performance; yet it is particularly interesting if one of the sets is a polytope that is not represented by its linear constraints. Common cases are the V-representation or zonotopes represented by their generators, which are used in several approaches [27, 28, 29, 30]. To compute the (exact) intersection of such a polytope $\mathcal{X}$ with a polyhedron $\mathcal{G}$ in H-representation, $\mathcal{X}$ needs to be converted to H-representation first. A polytope with $m$ vertices can have $O\!\left(\binom{m-n/2}{m-n}\right)$ linear constraints [42]. (For two polytopes in V-representation in general position there is a polynomial-time intersection algorithm [43], but this assumption is not practical.) A zonotope with $m$ generators can have $O\!\left(m\binom{m}{n-1}\right)$ linear constraints [29]. If $\mathcal{G}$ is a polytope in H-representation, checking disjointness of $\mathcal{X}$ and $\mathcal{G}$ can also be solved more efficiently in low dimensions, e.g., for $m$ linear constraints in $O(m)$ for $n\leq 3$ [44].

Next we discuss how to apply low-dimensional reasoning to the intersection $\widehat{\mathcal{X}}\cap\mathcal{G}$ of a decomposed set $\widehat{\mathcal{X}}$ and a polyhedron $\mathcal{G}$, or respectively detect emptiness of the intersection $\widehat{\mathcal{X}}\cap\mathcal{G}$ (which can often be achieved more efficiently). The key idea is to exploit that $\widehat{\mathcal{X}}$ is decomposed. For ease of discussion, we consider the case of two blocks (i.e., $\widehat{\mathcal{X}}=\mathcal{X}_{1}\times\mathcal{X}_{2}$). Below we discuss the two cases that $\mathcal{G}$ is decomposed or not.

Intersection of two decomposed sets

We first consider the case that $\mathcal{G}$ is also decomposed and agrees with $\widehat{\mathcal{X}}$ on the block structure, i.e., $\mathcal{G}=\widehat{\mathcal{G}}=\mathcal{G}_{1}\times\mathcal{G}_{2}$ and $\mathcal{X}_{1},\mathcal{G}_{1}\subseteq\mathbb{R}^{n_{1}}$ for some $n_{1}$. Clearly

[TABLE]

because the Cartesian product and intersection distribute; thus

[TABLE]

Now consider the second disjunct in (5) and assume that $\mathcal{G}_{2}$ is universal. We get $\mathcal{X}_{2}\cap\mathcal{G}_{2}=\emptyset\iff\mathcal{X}_{2}=\emptyset.$ In our context, $\widehat{\mathcal{X}}$ (and hence $\mathcal{X}_{2}$) is nonempty by construction. Hence (5) simplifies to

[TABLE]

so we never need to compute $\mathcal{X}_{2}$ to determine whether the intersection is empty. In practice, the set $\mathcal{G}^{*}$ from (4) takes the role of $\widehat{\mathcal{G}}$ and is often only constrained in some dimensions (and hence decomposed and universal in all other dimensions). We illustrate this algorithm in Figure 6 b).

Intersection of a decomposed and a non-decomposed set

If $\mathcal{G}$ is not decomposed in the same block structure as $\mathcal{X}$, we can still decompose it, at the cost of an approximation error. Let $\pi_{1}$ and $\pi_{2}$ be suitable projection matrices. Then

[TABLE]

and hence

[TABLE]

From (6) we obtain 1) a sufficient test for emptiness of $\widehat{\mathcal{X}}\cap\mathcal{G}$ in terms of only $\mathcal{X}_{1}$ and 2) a more precise sufficient test in terms of $\mathcal{X}_{1}$ and $\mathcal{X}_{2}$ in low dimensions. If both tests fail, we can either fall back to the (exact) test in high dimensions or conservatively assume that the intersection is nonempty.

The precision of the above scheme highly depends on the structure of $\widehat{\mathcal{X}}$ and $\mathcal{G}$. If several (but not all) blocks of $\mathcal{G}$ are constrained, instead of decomposing $\mathcal{G}$ into the low-dimensional block structure, one can alternatively compute the intersection for medium-dimensional sets to avoid an approximation error; we apply this strategy in the evaluation (Section V). If $\mathcal{G}$ is compact, the following proposition shows that the approximation error is bounded by the maximal entry in the diameters of $\widehat{\mathcal{X}}$ and $\mathcal{G}$, and this bound is tight. This strategy is illustrated in Figure 6 c).

Proposition 2

Let $\widehat{\mathcal{X}}=\bigtimes_{j}\mathcal{X}_{j}\in\mathcal{C}_{n}$, $\mathcal{G}\in\mathcal{C}_{n}$, $\widehat{\mathcal{X}}\cap\mathcal{G}\neq\emptyset$, $\widehat{\mathcal{G}}:=\bigtimes_{j}\pi_{j}\mathcal{G}$ for appropriate projection matrices $\pi_{j}$ corresponding to $\mathcal{X}_{j}$, and $p=\infty$. Then

[TABLE]

Example

Consider Figure 5(b). Previously we have already identified the intersection with the flowpipe for time steps $k=3$ and $k=4$. The resulting sets are $\widehat{\mathcal{X}}(k)\cap\mathcal{G}=\mathcal{X}_{1}(k)\cap\mathcal{G}_{1}\times\mathcal{X}_{2}(k)$, where $\mathcal{G}_{1}$ was the projection of $\mathcal{G}$ to $x_{1}$. We emphasize that we compute the intersections in low dimensions, that we need not compute $\mathcal{X}_{2}(1)$ and $\mathcal{X}_{2}(2)$ at all, and that in this example all computations are exact (i.e., we obtain the same sets as in Figure 3).

Now consider the case that $\mathcal{G}$ is not decomposed in the same structure as $\widehat{\mathcal{X}}$, e.g., the hyperplane $x_{1}=x_{2}$. One option is to decompose $\mathcal{G}$ to the blocks of $\widehat{\mathcal{X}}$, i.e., $\mathcal{G}_{1}:=\pi_{1}\mathcal{G}$, $\mathcal{G}_{2}:=\pi_{2}\mathcal{G}$ and compute the block-wise intersection $(\mathcal{X}_{1}(k)\cap\mathcal{G}_{1})\times(\mathcal{X}_{2}(k)\cap\mathcal{G}_{2})$. Here $\mathcal{G}_{1}$ and $\mathcal{G}_{2}$ are universal, so we obtain a coarse approximation of the intersection (namely $\widehat{\mathcal{X}}(k)$ itself). Alternatively, computing the intersection $\widehat{\mathcal{X}}(k)\cap\mathcal{G}$ is exact but computationally more demanding. If we assume that the system has higher dimension, e.g., $10$, then computing the intersection in two dimensions (i.e., with a 2D projection $(\mathcal{X}_{1}(k)\times\mathcal{X}_{2}(k))\cap\pi\mathcal{G}$) is still exact and yet cheaper than computing the intersection in full dimensions (i.e., $\widehat{\mathcal{X}}(k)\cap\mathcal{G}$).

We exemplify the possible difference between the intersection algorithms in Figure 7. There we compare the three algorithms visualized in Figure 6. The high-dimensional intersection is the most precise as expected because it does not suffer from a projection error, however, the result is a high-dimensional set, instead of a decomposed high-dimensional set. In addition, this algorithm intersects the sets in three dimensions (including $x_{3}$) while the other algorithms ignore the unconstrained dimensions. The medium-dimensional intersection corresponds to the box approximation of the true intersection. We observe that if we cannot decompose the polyhedron to the same block structure as the decomposed set without a projection error (e.g. projecting $x_{1}+x_{2}\leq 4$ onto $x_{1}$ and $x_{2}$ individually results in one-dimensional universes), then the low-dimensional intersection can produce a large approximation error. Thanks to the second linear constraint $x_{1}\leq 1.5$ and the fact that we treat each half-space separately, we are still able to compute a nontrivial intersection.

IV-C Decomposing an affine map

The next step after computing the intersection with the guard is the application of the assignment. We consider an affine map $A\widehat{\mathcal{X}}\oplus\{b\}$ with $A\in\mathbb{R}^{n\times n}$ and $b\in\mathbb{R}^{n}$. Affine-map decomposition has already been presented as part of the operator $\mathit{Post}_{C}^{\square}$ [2]:

[TABLE]

where $A_{ij}$ is the block in the $i$-th block row and the $j$-th block column. We recall an error estimation.

Proposition 3

[2, Prop. 3]** Let $\mathcal{X}=\bigtimes_{j=1}^{\beta}\mathcal{X}_{j}\in\mathcal{C}_{n}$ be nonempty, $A\in\mathbb{R}^{n\times n}$, $q_{j}:=\arg\max_{i}\|A_{ij}\|_{p}$ (the index of the block with the largest matrix norm in the $j$-th block column) so that $\alpha_{j}:=\max_{i\neq q_{j}}\|A_{ij}\|_{p}$ is the second largest matrix norm in the $j$-th block column. Let $\alpha_{\text{max}}:=\max_{j}\alpha_{j}$ and $\Delta_{\text{sum}}:=\sum_{j}\Delta_{\infty}(\mathcal{X}_{j})$. Then

[TABLE]

In particular, if only one block per block column of matrix $A$ is nonzero, the approximation is exact [2]. For example, consider a two-block scenario and a block-diagonal matrix $A$, i.e., $A_{12}=A_{21}=0$. Then

[TABLE]

In practice, affine maps with such a structure are unrealistic for the $\mathit{Post}_{C}^{\square}$ operator but typical for assignments. Prominent cases include resets, translations, and scalings, for which $A$ is even diagonal and hence block diagonal for any block structure.

Example

Recall that, after computing the intersections, we ended up with the same sets as in Figure 3. In our example, the assignment is a translation in dimension $x_{1}$. Hence, as mentioned above, the application of the decomposed assignment is also exact. In particular, the translation only affects $\mathcal{X}_{1}(k)$ and we obtain the same result as in Figure 3.

IV-D Inclusion check for decomposed sets

Our algorithm is now fully able to take transitions. Observe that all sets ever occurring in scheme (3) using the algorithm are decomposed. The following proposition gives an exact low-dimensional fixpoint check under this condition.

Proposition 4

Let $\widehat{\mathcal{X}}=\bigtimes_{j}\mathcal{X}_{j}\in\mathcal{C}_{n},\widehat{\mathcal{G}}=\bigtimes_{j}\mathcal{G}_{j}\in\mathcal{C}_{n}$ be nonempty sets with identical block structure. Then

[TABLE]

IV-E Decomposing a convex hull

As the last part of the algorithm, we decompose the computation of the convex hull. We exploit that all sets in the same flowpipe share the same block structure.

Proposition 5

Let $\widehat{\mathcal{X}}=\bigtimes_{j}\mathcal{X}_{j}\in\mathcal{C}_{n},\widehat{\mathcal{G}}=\bigtimes_{j}\mathcal{G}_{j}\in\mathcal{C}_{n}$ be nonempty sets with identical block structure. Then

[TABLE]

For the decomposition operations proposed before (intersection, affine map, and inclusion), there are common cases where the approximations were exact. In these cases it is always beneficial to perform the decomposed operations instead of the high-dimensional counterparts. The decomposition of the convex hull, however, always incurs an approximation error, which we can bound by the radius of the box approximation and by the block-wise difference in bounds.

Proposition 6

Let $\widehat{\mathcal{X}}=\bigtimes_{j}\mathcal{X}_{j}\in\mathcal{C}_{n},\widehat{\mathcal{G}}=\bigtimes_{j}\mathcal{G}_{j}\in\mathcal{C}_{n}$ be nonempty sets with identical block structure and let $r^{\infty}$ be the radius of the box approximation of $\mathrm{CH}(\widehat{\mathcal{X}}\cup\widehat{\mathcal{G}})$. Then

[TABLE]

Example

Figure 4 shows the decomposed convex hull of the sets $\widehat{\mathcal{X}}(3)\cap\mathcal{G}$ and $\widehat{\mathcal{X}}(4)\cap\mathcal{G}$ after applying the translation. Since each block is one-dimensional in our example, we obtain the box approximation.

IV-F Discussion

The combination of the decomposed set operations outlined above results in a reachability algorithm that is sound. This immediately follows from the fact that all operations compute an overapproximation.

Proposition 7** (Soundness)**

The analysis using Algorithm 1 and the aforementioned procedures to apply the decomposed set operations (intersection, affine map, inclusion check, and convex hull) results in an overapproximation of the reach set.

When using the above algorithm, the question about the choice of the decomposition arises. Unfortunately, there is generally no automatic way to find the “best” decomposition.

Regarding LTI systems, different decompositions generally result in different flowpipes [2]. Unless the LTI system consists of fully decoupled sub-components (which is usually not the case), the only “best” decomposition is the degenerate case of a single big block.

Regarding hybrid systems given a fixed decomposition of the LTI systems, there is indeed a unique “best” decomposition: a single block for all constrained dimensions. Finding this “best” decomposition is trivial. However, this block may be too big in terms of computational cost.

We can generally say that a finer partition in the decomposition (i.e., using more and smaller blocks) only ever degrades precision. However, while lower-dimensional set operations have lower computational cost, a decrease in precision may actually affect the rest of the algorithm in a negative way. For instance, we may not be able to prove that a transition cannot be taken. Thus statements about computational cost are generally not possible in the case of hybrid systems.

In the next section we shall see that in practical use cases the decomposition is always beneficial and that using the block size of the constrained dimensions is a good heuristic.

Our reachability algorithm benefits greatly from a low number of constrained dimensions, especially in high-dimensional models. High-dimensional models are scarcely available, but those models that we found supported the hypothesis that this number is typically low. In particular, the ARCH-COMP competition [9] represents the state of the art in reachability analysis and considers a scalable model of a powertrain from [39] where only one dimension is constrained. In the next section we consider another scalable model with just three constrained dimensions.

V Evaluation

We implemented the ideas presented in Section IV in JuliaReach [4, 3]. The code is publicly available [4]. We performed the experiments presented in this section on a Mac notebook with an Intel i5 CPU@3.1 GHz and 16 GB RAM.

V-A Benchmark descriptions

We evaluate our implementation on a number of benchmarks taken from the HyPro model library [45], from the ARCH-COMP 2019 competition [5], and a scalable model from [9]. To demonstrate the qualitative performance of our approach, we verify a safety property for each benchmark, which requires precise approximations in each step of the algorithm. We briefly describe the benchmarks below.

Linear switching system. This five-dimensional model taken from [45] is a linear hybrid system with five locations of different controlled, randomly-generated continuous dynamics stabilized by an LQR controller. The discrete structure has a ring topology with transitions determined heuristically from simulations. The safety property for this system is $x_{1}>-1.2$.

Spacecraft rendezvous. This model with five dimensions represents a spacecraft steering toward a passive target in space [46]. We use a linearized version of this model with three locations. We consider two scenarios, one where the spacecraft successfully approaches the target and another one where a mission abort occurs at $t=120$ min. For the safety properties we refer to [5].

Platooning. This ten-dimensional model with two locations represents a platoon of three vehicles with communication loss at deterministic times [47]. The safety property enforces a minimum distance $d$ between the vehicles: $\bigwedge_{x\in\{x_{1},x_{4},x_{7}\}}x\geq-d$. We consider both a time-bounded setting with $d=42$ and a time-unbounded setting with $d=50$; note that in the latter setting a fixpoint must be found.

Filtered oscillator. This scalable model consists of a two-dimensional switched oscillator (dimensions $x$ and $y$) and a parametric number of filters (here: 64–1024) which smooth $x$ [9]. We fixed the maximum number of transitions to five by adding a new variable. The safety property is $y<0.5$.

V-B Tool descriptions

We compare our implementation in JuliaReach to two other algorithms available in the same tool. All three algorithms use the same decomposition-based continuous-post operator $\mathit{Post}_{C}^{\square}$ [3], so the main difference between these algorithms is the intersection operation with discrete jumps, which allows for a direct comparison of the approach presented in this paper. The existing implementations can be seen as instances of the algorithm in Section III. Furthermore, we compare the implementation to two algorithms available in SpaceEx [9], which is an efficient and mature verification tool for linear hybrid systems. We summarize the different approaches below.

Deco. This algorithm implements the decomposed approach presented in this paper. To compute the (low-dimensional) intersections, we use a polyhedra library [48].

LazySupp. This algorithm uses a (lazy) support-function-based approximation of the intersection operation using the simple heuristic $\rho_{X\cap Y}(\ell)\leq\min({\rho_{X}(\ell),\rho_{Y}(\ell)})$. This heuristic is fast but not precise enough to verify the safety property of the filtered oscillator model.

LazyOptim. In contrast to the coarse intersection in LazySupp, this algorithm uses a more precise implementation of the intersection operation based on line search [26].

SpaceEx LGG. This is an efficient implementation of the algorithm by Le Guernic and Girard [35]. The algorithm uses a support-function representation and can hence check low-dimensional conditions efficiently, but operations such as intersections work in high dimensions.

SpaceEx STC. The STC algorithm is an extension of the LGG algorithm with automatic time-step adaptation [49].

All tools use a support-function representation of sets. We use template polyhedra with box constraints to overapproximate these sets, which roughly corresponds to a decomposition into one-dimensional blocks. This approximation is fast.

For the algorithms implemented in JuliaReach, we use one-dimensional block structures for all models. For SpaceEx we use the options given in the benchmark sources. For algorithms that use a fixed time step, we fix this parameter to the same value for each benchmark. The SpaceEx STC algorithm does not have such a parameter, so we instead fix the parameter “flowpipe tolerance,” which controls the approximation error, to the following values: Linear switching system: $0.01$, Spacecraft rendezvous: $0.2$, Platooning: $1$, Filtered oscillator: $0.05$.

V-C Experimental results

The results are presented in Table I. We generally observe a performance boost of the Deco algorithm for models with more than five dimensions, and only a minor overhead for “small” models. This demonstrates the general scalability improvement by performing operations, especially the intersection, in low dimensions. Since all models have a small number of constrained dimensions in their guards and invariants, choosing an appropriate block structure results in very low-dimensional sets for computing the intersections, for which concrete polyhedral computations are very efficient and most precise. We note that such an intersection computation does not scale with the dimension, and so other algorithms must resort to approximation techniques.

Moreover, we found that our approach scales more favorably compared to the high-dimensional approaches when decreasing the time step $\delta$ (cf. Section III-A). We demonstrate this observation for the filtered oscillator model in Table II and explain it as follows. Recall that we only compute those sets in the flowpipe in high dimensions for which we have detected an intersection in low dimensions. With a smaller time step, the total number of sets increases and hence the savings due to our approach become more dominant. To give an example, for the filtered oscillator with time step $0.0005$, out of the $9{,}661$ sets in total we only computed $1{,}400$ sets in high dimensions. Consequently, making the time step $20$ times smaller only makes the run time $11.4$ times slower for our algorithm as opposed to at least $15$ times slower for other algorithms.

Furthermore, while in general our algorithm may induce an additional approximation error with the block structure, in the evaluation it was always precise enough to prove the safety properties. In Table I we report the total amount of constrained dimensions in guards, invariants, and safety properties for each benchmark instance. These are the dimensions that determine our low-dimensional flowpipes. For some of the models, the constrained dimensions differ between the invariants/guards and the safety property, but a one-dimensional block structure was still sufficient. Note that we need to compute intersections only with invariants and guards; for safety properties we just need to check inclusion.

In the linear switching model, only one dimension is constrained, so our algorithm, together with one-dimensional block structure, is appropriate, especially because a small time step is required to prove the property. In the platooning model, the safety property constrains three dimensions but the invariants and guards constrain just one dimension.

The invariants and guards of the models spacecraft and filtered oscillator constrain two dimensions, so the natural choice for the decomposition is to keep these dimensions in the same block. However, we chose to decompose into one-dimensional blocks for better run-time comparison in Table I. In the implementation, we follow the strategy to perform the intersection in two dimensions and then project back to one dimension (cf. Section IV-B). If instead we decide to employ a two-dimensional block structure, we can further gain precision by using different template polyhedra, e.g., with octagon constraints. In Figure 8 we visualize this alternative for the filtered oscillator. As expected, the one-dimensional analysis is less precise as it corresponds to a box approximation. However, we note that a one-dimensional block structure also inherently reduces the precision of $\mathit{Post}_{C}^{\square}$, so the additional approximation error does not only stem from the handling of discrete transitions by to our approach. In addition, one can observe that the approximation using a one-dimensional block structure with low-dimensional intersection is coarser compared with the medium-dimensional strategy and we cannot prove the safety property anymore.

V-D Scaling the number of constrained dimensions

We now investigate the effect of increasing the number of constrained dimensions in two experiments.

In the first experiment we focus on the intersection operation. We fix a hypercube $\mathcal{X}$ of dimension $1{,}024$, of radius $4$ and centered in the origin, and a half-space $\mathcal{G}$ $x_{1}+\dots+x_{k}\leq 2$ for parameter $k$ that controls the number of constrained dimensions. The half-space properly cuts the hypercube for any $k$. In Figure 9 we show the run times for different intersection algorithms. We note that these algorithms compute different approximations of the true intersection: the “DecoLow” and “DecoMedium” algorithms implement the respective ideas from Figure 6; the “LazySupp” and “LazyOptim” algorithms use the ideas as in the reachability algorithms of the same name (described in Section V-B); the “Concrete” algorithm uses a polyhedra library (note that the H-representation of a hypercube of dimension $n$ has only $2n$ constraints). We can see that the high-dimensional intersection algorithms are not affected by the number of constrained dimensions $k$ whereas the decomposition-based algorithms scale gracefully with $k$.

In the second experiment we consider the full reachability setting. For that we modify the filtered oscillator benchmark with $128$ filters (filtered_osc128 in Table I) by adding small nonzero entries to $k$ previously unconstrained dimensions in the invariants and guards. We consider all the reachability algorithms that are precise enough to verify the safety property (which is still satisfied in our modified benchmark instances), i.e., all the algorithms from Table I except for LazySupp. Note that for SpaceEx algorithms we extrapolate the results we report for this benchmark instance in Table I to a range of values of parameter $k$ as the reachability algorithms implemented in SpaceEx do not depend on the number of constrained dimensions. The run times are plotted in Figure 10 as a function of $k$. Again the high-dimensional algorithms are not affected by the number $k$ whereas the decomposition-based algorithm scales well with $k$ (note the log scale). For $k\approx 120$ (close to $n=131$) we see a cross-over with the SpaceEx LGG run time, which is expected due to the overhead of the decomposition.

VI Conclusion

We have presented a schema that integrates a reachability algorithm based on decomposition for LTI systems in the analysis loop for linear hybrid systems. The key insight is that intersections with polyhedral constraints can be efficiently detected and computed (approximately or often even exactly) in low dimensions. This enables the systematic focus on appropriate subspaces and the potential for bypassing large amounts of flowpipe computations. Moreover, working with sets in low dimensions allows to use precise polyhedral computations that are infeasible in high dimensions.

Future work

An essential step in our algorithm is the fast computation of a low-dimensional flowpipe for the detection of intersections. In the presented algorithm, we recompute the flowpipe for the relevant time frames in high dimensions using the same decomposed algorithm with the same time step. However, this is not necessary. We could achieve higher precision by using different algorithmic parameters or even a different, possibly non-decomposed, algorithm (e.g., one that features arbitrary precision [25]). This is particularly promising for LTI systems because one can avoid recomputing the homogeneous (state-based) part of the flowpipe [26].

In the benchmarks considered in the experimental evaluation, it was not necessary to change the block structure when switching between locations. In general, different locations may constrain different dimensions, so tracking the “right” dimensions may be necessary to maintain precision. While it is easy to merge different blocks, subsequent computations would become more expensive. Hence one may also want to split blocks again for optimal performance. Since this comes with a loss in precision, heuristics for rearranging the block structure, possibly in a refinement loop, are needed.