JuliaReach: a Toolbox for Set-Based Reachability

TL;DR

JuliaReach is a versatile toolbox for set-based reachability analysis of dynamical systems, combining lazy and explicit set computations in Julia for efficient, large-scale problem solving.

Contribution

It introduces a new Julia toolbox with both reachability algorithms and convex set calculus, enabling efficient analysis of complex systems.

Findings

Achieves runtime performance comparable to compiled languages.

Handles high-dimensional lazy set operations efficiently.

Applicable to large-scale, complex dynamical systems.

Abstract

We present JuliaReach, a toolbox for set-based reachability analysis of dynamical systems. JuliaReach consists of two main packages: Reachability, containing implementations of reachability algorithms for continuous and hybrid systems, and LazySets, a standalone library that implements state-of-the-art algorithms for calculus with convex sets. The library offers both concrete and lazy set representations, where the latter stands for the ability to delay set computations until they are needed. The choice of the programming language Julia and the accompanying documentation of our toolbox allow researchers to easily translate set-based algorithms from mathematics to software in a platform-independent way, while achieving runtime performance that is comparable to statically compiled languages. Combining lazy operations in high dimensions and explicit computations in low dimensions, JuliaReach can be applied to solve complex, large-scale problems.

Tables

Table 1. Table 1. Runtimes (in seconds) for different discrete-post operators in the Reachability package and the SpaceEx tool.

Number of filters $m$	2	4	8	16	32	64	128	196	256
Reachability – concrete	$0.60$	$0.84$	$1.77$	$4.24$	$12.18$	$39.23$	$137.89$	$341.38$	$633.86$
Reachability – lazy	$0.55$	$0.75$	$1.41$	$3.19$	$6.39$	$22.65$	$64.63$	$145.49$	$262.91$
Reachability – lazy with overapproximation	$0.57$	$0.75$	$1.17$	$4.76$	$6.18$	$20.05$	$47.78$	$103.12$	$212.10$
SpaceEx	$0.11$	$0.22$	$0.11$	$0.30$	$1.10$	$6.78$	$56.72$	$262.88$	$1,011.17$

Equations

$$Y=\text{CH}(e^{A\delta }{X}_0\oplus \delta BU,X_0),$$

$$A_1=\left(\begin{array}{cc}-1& -4\\ 4& -1\end{array}\right)A_2=\left(\begin{array}{cc}1& 4\\ -4& -1\end{array}\right)$$

$$\dot{x}$$

$$\dot{y}$$

$${\dot{x}}_1$$

$${\dot{x}}_2$$

$${\dot{x}}_3$$

$$⋮$$

$${\dot{x}}_m$$

$$x$$

$$y+0.714286\ast x$$

$$\dot{x}$$

$$\dot{y}$$

$${\dot{x}}_1$$

$${\dot{x}}_2$$

$${\dot{x}}_3$$

$$⋮$$

$${\dot{x}}_m$$

$$x\le 0$$

$$y+0.714286\ast x\le 0$$

$$\dot{x}$$

$$\dot{y}$$

$${\dot{x}}_1$$

$${\dot{x}}_2$$

$${\dot{x}}_3$$

$$⋮$$

$${\dot{x}}_m$$

$$x$$

$$y+0.714286\ast x$$

$$\dot{x}$$

$$\dot{y}$$

$${\dot{x}}_1$$

$${\dot{x}}_2$$

$${\dot{x}}_3$$

$$⋮$$

$${\dot{x}}_m$$

$$x$$

$$y+0.714286\ast x$$

$$x$$

$$y+0.714286\ast x$$

$$x$$

$$y+0.714286\ast x$$

$$x$$

$$y+0.714286\ast x$$

$$x$$

$$y+0.714286\ast x$$

Full text

reachability analysis, hybrid systems, lazy computation

††journalyear: 2019††copyright: author††conference: 22nd ACM International Conference on Hybrid Systems: Computation and Control; April 16–18, 2019; Montreal, QC, Canada††booktitle: 22nd ACM International Conference on Hybrid Systems: Computation and Control (HSCC ’19), April 16–18, 2019, Montreal, QC, Canada††price: 15.00††doi: 10.1145/3302504.3311804††isbn: 978-1-4503-6282-5/19/04††ccs: Security and privacy Logic and verification

1. Introduction

Set-based reachability analysis is a rigorous approach to reason about the behavior of a dynamical system (Alur et al., 1995). Many well-known approaches to reachability analysis are based on flowpipe construction, i.e., computing a finite cover of all trajectories starting from a given set of initial states (Guernic and Girard, 2009; Frehse et al., 2011; Chen et al., 2013; Althoff, 2015; Bak and Duggirala, 2017; Bogomolov et al., 2012; Ray et al., 2015). Such approaches, while being based on different theoretical ideas, conceptually share underlying machinery. Thanks to these common grounds, when implementing a flowpipe construction approach (in a prototype or in an end-user tool), our vision is that, instead of starting from scratch, one can reuse the code from a related implementation. To make such reuse convenient, it is essential to provide a flexible interface for common operations, supported by an efficient library with rich functionality. Given these features, researchers and their students can then concentrate on the core algorithms and ideally prototype their approach with relatively few lines of code. In practice, however, many of the existing approaches based on flowpipe construction are implemented in different tools without any shared code base.

An important element in this context is the choice of programming language. In academic prototyping, the choice is usually a consideration between convenience during development and output performance in the final product. On the one hand, compiled languages such as C++ offer high performance, but the compilation overhead is inconvenient for prototyping. On the other hand, interpreted languages such as Python offer an interactive session for convenient prototyping, but these languages fall behind in performance or need to extend the code to work with another lower-layer program such as Numba or Cython (known as the two-language problem). A compromise between the two worlds are just-in-time (JIT) compiled languages such as MATLAB.

Julia, a recent111Version 1.0 of Julia was released in August 2018. programming language for high-performance scientific computing (Bezanson et al., 2017), reconciles the two advantages of compiled and interpreted languages described above, as it comes with an interactive read-evaluate-print loop (REPL) front-end, but is JIT compiled to achieve performance that is competitive with compiled languages such as C. A distinctive feature of Julia is multiple dispatch (i.e., the function to execute is chosen based on each argument type), which allows to write efficient machine code based on a given type, e.g., of the set. As additional features, Julia is platform independent, has an efficient interface to C, is supported in Jupyter notebooks and well-suited for parallel computing. Julia has a determined and quickly-growing community, especially for scientific tools.222https://julialang.org We believe that all this makes Julia an interesting programming language for writing a tool for reachability analysis.

In this paper, we present JuliaReach 333http://juliareach.org, an open-source toolbox for rapid prototyping of set-based reachability approaches, written in Julia. The JuliaReach toolbox consists of two main packages: The first package, Reachability, offers core infrastructure for implementing reachability algorithms for continuous and hybrid systems. As a proof of concept, the package currently provides implementations of an algorithm for the analysis of linear time-invariant (LTI) systems based on decomposition (Bogomolov et al., 2018), and of the classic interleaving algorithm of continuous and discrete-post operators as used in hybrid system tools such as SpaceEx (Frehse et al., 2011; Bogomolov et al., 2014, 2015).

The second package, LazySets, is a library for state-of-the-art calculus with convex sets. As the name indicates, a key feature of LazySets is lazy set representation using support functions. For making the lazy sets concrete, the library offers means to overapproximate them, e.g., using template directions or $\varepsilon$-close approximation. LazySets also supplies concrete set representations, e.g., polyhedra, zonotopes, and ellipsoids.

A binary file to use Julia can be downloaded on the official web page, and installing the JuliaReach packages is just one command. The toolbox puts attention on high efficiency and usability, both from the end-user side and the developer side, i.e., the code is thoroughly tested and documented, including tutorials, and written with extensibility in mind. Code contributions follow a continuous-integration scheme and are peer-reviewed to assure high quality.

Related Work

Several reachability tools are implemented in C++. For example, SpaceEx is a C++ package that currently integrates several independent reachability algorithms for linear hybrid systems (Frehse et al., 2011). Flow ∗implements Taylor model approximation for nonlinear hybrid systems in C++ (Chen et al., 2013).

There are also other reachability tools that are written in JIT-compiled or interpreted languages. For example, CORA implements several algorithms for linear and nonlinear hybrid systems in MATLAB (Althoff, 2015), and also offers a standalone set library focused on zonotopes, as well as a library for interval arithmetic. A subset of the algorithms from CORA have been implemented in C++ in the tools CORA/SX and SymReach, leading in some instances to a speedup of up to a factor of 3 (Immler et al., 2018). A few tools are implemented in Python. For instance, Hylaa implements discrete-time reachability algorithms for linear hybrid systems  (Bak and Duggirala, 2017).

HyPro is a C++ library for convex set representation with similar aims as the LazySets library (Schupp et al., 2017). It provides a common interface for set representations like boxes, convex polyhedra, support functions, and zonotopes. Recently, a reachability tool called HyDRA (Althoff et al., 2018) has been implemented on top of HyPro. At the time of writing, HyDRA is not publicly available.

In the remainder of the paper, we describe the LazySets library in Section 2 and the Reachability package in Section 3. In Section 4, we demonstrate the viability of JuliaReach in two case studies.

2. The LazySets Library

This section outlines the LazySets library. In the library we consider representations of closed convex sets in the usual sense from convex geometry: A set is closed if it contains all its limit points. A set $S$ is convex if for any $m$ points $v_{j}\in S$ and $m$ non-negative numbers $\lambda_{j}$ that sum up to 1 we have that $\sum_{j}\lambda_{j}v_{j}\in S$ as well. Alternatively, a closed convex set is an intersection of (possibly infinitely many) closed half-spaces.

Every convex set type in the library inherits from the parametric abstract type LazySet{N}, where N is a parameter for the numeric type. This way one can easily switch between, e.g., floating point (Float64) and exact (Rational) precision with no additional performance penalty: At runtime, Julia uses multiple dispatch on N and JIT-compiles into type-specific code.

The library comes with a collection of common set representations, such as balls in different norms, ellipsoids, hyperrectangles, polyhedra, polytopes (i.e., bounded polyhedra) in constraint and in vertex representation, and zonotopes. One key feature of LazySets is to apply common binary set operations in a lazy fashion, i.e., to not construct the result of an operation explicitly. For that purpose, we use special wrapper types to represent operations such as convex hull, Minkowski sum, linear and exponential map, intersection, and Cartesian product. Since these lazy set operations themselves are subtypes of LazySet, they can also be nested.

New subtypes of LazySet can be easily added. To realize the lazy paradigm, a set type should provide an implementation of two functions: The first function has signature $\sigma$($d$, $S$) in LazySets and returns some furthest point of a set $S$ in a given direction $d$, i.e., its support vector. We recall that the set of support vectors corresponds to the optimal points for the support function $\rho$, i.e., $\sigma_{S}(d):=\{x\in S\mid d^{\mathrm{T}}x=\rho_{S}(d)\}$, $\rho_{S}(d):=\max_{x\in S}d^{\mathrm{T}}x$.

The second function $\rho_{S}(d)$ has signature $\rho$($d$, $S$) in LazySets and is exactly the support function. Since the default implementation of the support function is given by the above formula, it suffices to define the function $\sigma$; but if a more efficient implementation is available, it will be used automatically by the dispatch machinery in Julia. For the above lazy operations, the support vector can be evaluated efficiently without explicitly representing the set resulting from the operation (Le Guernic and Girard, 2010; Le Guernic, 2009), with the exception of intersection. For intersection, we only provide (over)approximate values for $\rho$, where we can either use a line search algorithm inspired by (Frehse and Ray, 2012) or a coarser heuristics given as $\rho_{\mathcal{X}\cap\mathcal{Y}}(d)\leq\min(\rho_{\mathcal{X}}(d),\rho_{\mathcal{Y}}(d))$ (Le Guernic, 2009).

LazySets defines several abstract supertypes (sometimes called interfaces in other programming languages) to bundle common functionality. For instance, all polytopic set types should subtype the AbstractPolytope type and implement methods to obtain their vertex representation and their constraint representation. This allows to write generic (or “virtual”) functions for arbitrary sets that belong to the AbstractPolytope family. If later a new set type is added to this family, all such generic functions are available.

2.1. The lazy paradigm

For illustration, consider the linear map of a set $S$ by a matrix $M$. In JuliaReach, we can either write LinearMap($M$, $S$) or use the short hand $M$ * $S$ for convenience. The LinearMap type has two fields that hold the map and the set, respectively. If $S$ is itself a LinearMap instance, the constructor multiplies the two matrices immediately. Otherwise, it creates a LinearMap instance that wraps $M$ and $S$ instead of computing the map explicitly. Given a direction $d$, the support vector computation is based on the formula $M\cdot\sigma(M^{T}d,S)$, i.e., it asks for the support vector of the wrapped set $S$.

As a second example, consider two polytopes $P_{1}$ and $P_{2}$. The command $P_{1}\oplus P_{2}$ instantiates a new MinkowskiSum instance. Similar to the linear map, the support vector can be defined recursively as $\sigma_{P_{1}\oplus P_{2}}(d)=\sigma_{P_{1}}(d)\oplus\sigma_{P_{2}}(d)$. Again, the binary operation MinkowskiSum is defined between any two LazySet types. On the other hand, this operation can be performed explicitly using the command minkowski_sum($P_{1}$, $P_{2}$), which calls an external library called Polyhedra.jl (Legat et al., 2018),that provides a unified interface to well-known implementations of polyhedral computations, such as CDD or LRS.

The approach lazy-by-default and optionally-explicit is also available with other set types. For example, we can execute the concrete linear_map function for Zonotopes, in which case the generators will be transformed.

In practice it is convenient to switch from a lazy to a concrete representation at some point, ideally without an exponential increase in the computational cost when dealing with high-dimensional problems. Below we describe the Approximations sub-module, which serves this purpose.

2.2. From lazy to concrete set representation

The Approximations module implements the transfer from a lazy set representation to a concrete set, generally involving an overapproximation. As a particular feature, the module can be used to combine lazy high-dimensional sets with explicit low-dimensional approximations of projections. For illustration, consider the following example. A typical set equation for a discrete approximation model of LTI systems is

[TABLE]

with real matrix $A\in\mathbb{R}^{n\times n}$, time step $\delta$, initial states $X_{0}$, nondeterministic inputs $BU$, and CH denoting the convex hull. For example, let $\delta=0.1$, $X_{0}\subseteq\mathbb{R}^{n}$ be a ball with center $(1,\dots,1)$ and radius $0.1$ in the infinity norm, $U\subseteq\mathbb{R}^{m}$ be a ball centered in the origin with radius $1.2$ in the 2-norm, and $B$ be a linear map of appropriate dimensions. With LazySets, Eq. (1) for $n=1{,}000$, $m=2$, and random matrix coefficients is

⬇

n = 1000; m = 2; $\delta$ = 0.1

A = sprandn(n, n, 0.01); B = randn(n, m)

X0 = BallInf(ones(n), 0.1)

U = Ball2(zeros(m), 1.2)

Y = ConvexHull(SparseMatrixExp(A*$\delta$)*X0$\,\oplus\,\delta$BU, X0)

The execution is instantaneous because we just created a nested lazy set. Note the use of SparseMatrixExp($M$), a wrapper around the matrix exponential $e^{M}$; the evaluation of the support vector for the lazy matrix exponential relies on computing the matrix action, a technique taken from the numerical ODE/PDE domain.

The Approximations module can now be used to get information about Y without having computed any concrete representation; e.g., projection of a set into two dimensions is useful for visualization or verification of properties that involve only two variables. For this purpose, the module offers the function overapproximate($S$, $\varepsilon$) where $S$ is a set and $\varepsilon$ is an error bound. The function overapproximates a 2D set by adding supporting directions until the error bound is achieved (measured in terms of the Hausdorff distance), where the number of directions is optimal. The algorithm implements Kamenev’s method (Kamenev, 1996; Lotov and Pospelov, 2008), which we sketch in Fig. 1.

Suppose now that we are interested in observing the projection of Y onto the variables 1 and 50. First we define the $2\times n$ projection matrix and apply it as a linear map (i.e., from the left). Then we use the overapproximate method with a specified value of $\varepsilon$:

⬇

$\pi$ = spzeros(2, n); $\pi$[1, 1] = $\pi$[2, 50] = 1

res = overapproximate($\pi$*Y, $\varepsilon$)

Fig. 2 shows some results for different values of $\varepsilon$.

For high-dimensional overapproximation, instead of $\varepsilon$ one can alternatively pass a set type (e.g., Hyperrectangle for box approximation) or a special iterator type for template directions. We note again that one can call the overapproximate function with an arbitrary (nested) lazy set, as long as it implements $\rho$($d$, $S$) to evaluate the support function for a given direction.

3. The Reachability Package

The purpose of the Reachability package is to provide the basic infrastructure for implementing a custom reachability algorithm for continuous and hybrid systems, modeled as a hybrid automaton (Alur et al., 1992). We assume that the reader is familiar with the basic terminology of hybrid automata.

3.1. Infrastructure

Reachability provides the basic infrastructure that is needed or at least convenient to have for writing reachability algorithms.

An Options type stores user-defined options as key-value pairs that can be passed around to the solver backends.

Dynamical systems are passed to Reachability as special wrapper types. The reachability algorithms themselves are interfaced with post operators, which we will describe later. For hybrid systems, it is common to use clustering algorithms and fixed-point checks. Reachability offers standard convex-hull clustering and checks for a fixed-point (i.e., inclusion) after every jump, either before or after clustering (controlled by the user).

We support post-processing the resulting sets with an output function or a projection to given output dimensions. The sets can be visualized with the external Plots.jl tool using a single command.

3.2. Continuous-Post Operators

The interface for continuous-post operators is rather simple. Based on an :algorithm option, we just call the respective interface function with the continuous model, a ReachSet object, and a list of options (the number of steps, the step size, and other algorithm-specific settings). A ReachSet is a wrapper of a set and a corresponding uncertain time interval. Initially we would call the algorithm with the initial states and the time point 0 as in (X0, [0, 0]), but in the context of a hybrid system, the uncertain time interval will grow with each jump. The expected return type is just a sequence of ReachSets, sometimes called a reach tube.

For out-of-the-box usage, we provide an efficient implementation of a recent approach for the reachability analysis of LTI systems (Bogomolov et al., 2018; Althoff et al., 2018). The idea is to first decompose the system into low-dimensional blocks and then solve many small reachability problems. While the decomposition comes with an approximation error, the approach is suited for large-scale sparse systems.

3.3. Discrete-Post Operators

For analyzing hybrid systems, in addition to a post operator for continuous systems, one needs a post operator for the discrete transitions. Given a reach tube $\mathcal{X}$, a transition with a guard $\mathcal{G}$ and assignment $asgn(\cdot)$, and a target invariant $\mathcal{I}$, the discrete-post operator should compute (an overapproximation of) the set $asgn(\mathcal{X}\cap\mathcal{G})\cap\mathcal{I}$ and pass this set on to the continuous-post operator again.

Having both a continuous and a discrete-post operator, one can implement a classic algorithm that essentially interleaves the two post operators. Reachability offers the function “solve”, which implements this algorithm for given post operators.

We implemented two discrete-post operators to instantiate the above algorithm. The first operator uses (concrete, i.e., non-lazy) polytopes in constraint representation and relies on the external library Polyhedra.jl to perform the above-mentioned operations. In particular, the intersections are computed explicitly.

The second discrete-post operator is lazy but offers the options to overapproximate any of the intersections for obtaining a concrete set. Since evaluating a nested lazy intersection is a computationally complex task, this operator uses the line-search algorithm (see Section 2) only on the highest nesting level, and uses the coarser heuristics on lower nesting levels. If one is only interested in the support vector in one specific direction, this purely lazy approach scales very well. For more directions, however, there is no clear recipe which of the intersections should be performed lazily and which one should be overapproximated.

In the next section, we compare the performance of these two operators in a case study. We will instantiate the lazy discrete-post operator once with both intersections kept lazy and once with both intersections overapproximated.

4. Case Study

We perform two case studies that investigate some applications. In the first one, we re-implement a classic approach for reachability analysis of LTI systems using the LazySets library. We then add a custom extension for hybrid systems, demonstrating the rapid-prototyping capabilities of Julia. In the second case study, we consider Reachability and compare the performance of different discrete-post operators outlined in the previous section.

4.1. Working with the LazySets Library

We wrote a full-fledged re-implementation of Girard’s reachability algorithm for LTI systems with zonotopes (Girard, 2005) in Julia. The goal is to demonstrate prototyping with LazySets and the general applicability as a standalone library. The complete Julia source code is given in Algorithm 1, where we underlined those parts that involve the LazySets library.

We then extended the example presented in (Girard, 2005, Section 4) to a hybrid system of two modes $\ell_{i}$, $i=1,2$, with initial states $[0.9,1.1]\times[-0.1,0.1]$ and uncertain inputs from a set $u$ with $\mu=\|u\|_{\infty}=0.001$. The dynamics matrices $A_{i}$ are defined as follows:

[TABLE]

We add a transition $t_{i}$ from mode $\ell_{i}$ to $\ell_{3-i}$ with a hyperplane guard $g_{i}$: $g_{1}\triangleq x_{1}=-0.5$, $g_{2}\triangleq x_{2}=-0.3$. LazySets offers order reduction for zonotopes, which we used here with an upper bound of 10 generators. We implemented a custom reachability algorithm for this case study, again only using the LazySets capabilities. For simplicity, a transition is taken as soon as an intersection with the guard is detected. The complete source code is given in Appendix A.We plot the reachable states for the time interval $[0,4]$ and time step $\delta=0.001$ in Fig. 3. The analysis takes 0.25 seconds.

4.2. Working with the Reachability Framework

In the second case study, we show how the Reachability framework can be used to quickly evaluate a new approach to reachability analysis. Consider the case where a researcher came up with a new strategy for the discrete-post operator and now wants to compare to existing approaches. The researcher would just need to implement the different approaches as different post operators in the Reachability package and then evaluate the performance. Here we simulate this situation by using the operators that currently exist in the Reachability package (see Section 3.3).

We compare the performance of the different discrete-post operators using a model of a filtered oscillator (Frehse et al., 2011). The model represents a parameterized hybrid system consisting of ($i$) a two-dimensional switched oscillator in the variables $x$ and $y$, and ($ii$) a filter with $m$ state variables $x_{1},\dots,x_{m}$. We show the hybrid automaton in Fig. 4. The filter smooths $x$ with $x_{m}$ as an output signal, and the amplitude decreases with increasing dimension of the filter. The initial set is given as $x\in[0.2,0.3],y\in[-0.1,0.1]$, and $x_{j}=0$ for all $x_{j}$.

As a reference tool, we use SpaceEx (in version 0.9.8f) to evaluate the overall performance. To compare only the performance of the discrete-post operators, we fix the number of jumps to four, which corresponds to one loop in the automaton. Reachability offers a simple option for that purpose, but SpaceEx does not. To have a fair comparison, we do not use our option and instead modify the model in the following way. We add a new variable $b$ (for bound) with initial value 1. On the first transition (at the bottom in Fig. 4) we add an assignment $b:=2*b$. Finally, we augment the invariant of the target location (lower right) by the constraint $b\leq 2$.

We assume the following benchmark settings. The time interval is $[0,20]$ for filters of dimension up to four, and $[0,99]$ for higher dimensions. The time step is $\delta=0.01$ for both Reachability and SpaceEx. We ran the benchmarks on a computer with a 2.20 GHz CPU and 8 GB RAM. Table 1 outlines the performance of the different discrete-post operators and SpaceEx. For SpaceEx, we used the support-function algorithm (“supp”). Fig. 5 shows some plots of the output variables for the lazy discrete-post operator with overapproximation. (We used box approximations in this evaluation.)

As expected, the concrete intersection with polytopes is the slowest. (For this model, all computations could be performed in constraint representation.) The lazy-intersection approaches (second and third row in the table) scale much better, both roughly in the same order. Thanks to the fast continous-post algorithm, all operators outperform SpaceEx for the biggest model instances.

5. Conclusion

We have presented the JuliaReach toolbox, a new environment for developing reachability algorithms for dynamical systems. JuliaReach is written in the modern programming language Julia that unifies rapid prototyping and high performance. The toolbox consists of a framework for reachability analysis (Reachability) that is built around a standalone library for calculus with convex sets (LazySets).

In the future, we will extend the amount of post operators offered in Reachability. In particular, we are working on a new discrete-post operator that exploits the low-dimensional output of the continuous-post operator in the hybrid loop.

Acknowledgements.

This work was partially supported by the Sponsor Air Force Office of Scientific Research http://www.wpafb.af.mil/afrl/afosr/ under award number Grant #FA2386-17-1-4065, by the ARC project Grant #DP140104219 (Sponsor Robust AI Planning for Hybrid Systems https://cs.anu.edu.au/research/research-projects/robust-ai-planning-hybrid-systems), by the Sponsor Austrian Science Fund (FWF) https://www.fwf.ac.at/ under grants Grant #S11402-N23 (RiSE/SHiNE) and Grant #Z211-N23 (Wittgenstein Award), and the Sponsor European Union’s Horizon 2020 research and innovation programme https://ec.europa.eu/programmes/horizon2020/ under the Marie Skłodowska-Curie grant agreement No Grant #754411. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the United States Air Force.

Appendix A Hybrid reachability implementation

The following Julia function is an embedding of Algorithm 1 in a simple algorithm for hybrid systems. Apart from LazySets, no other library is required. We underlined the places where LazySets is used.