JuliaReach: a Toolbox for Set-Based Reachability

JuliaReach: a T oolb ox for Set-Based Reachability Sergiy Bogomolov A ustralian National University Canberra, Australia Marcelo Forets CURE - UdelaR Maldonado, Uruguay Goran Frehse ENST A ParisT ech - U2IS Palaiseau Cedex, France Kostiantyn Potomkin A ustralian National University Canberra, Australia Christian Schilling IST A ustria Klosterneuburg, A ustria ABSTRA CT W e present JuliaReach , a toolbox for set-base d reachability analy- sis of dynamical systems. JuliaRea ch consists of two main pack- ages: 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 conve x sets. The library oers both concr ete and lazy set representations, where the latter stands for the ability to delay set computations until they are needed. The choice of the program- ming 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 com- piled languages. Combining lazy operations in high dimensions and explicit computations in low dimensions, JuliaRea ch can b e applied to solve complex, large-scale problems. CCS CONCEPTS • Security and privacy → Logic and verication ; KEY W ORDS reachability analysis, hybrid systems, lazy computation A CM Reference Format: Sergiy Bogomolov , Marcelo Forets, Goran Frehse, Kostiantyn Potomkin, and Christian Schilling. 2019. JuliaReach: a T oolbox for Set-Based Reachabil- ity . In 22nd A CM International Conference on Hybrid Systems: Computation and Control (HSCC ’19), A pril 16–18, 2019, Montreal, QC, Canada. ACM, New Y ork, N Y, USA, 7 pages. https://doi.org/10.1145/3302504.3311804 1 IN TRODUCTION Set-based reachability analysis is a rigorous approach to reason about the behavior of a dynamical system [ 3 ]. Many well-known approaches to r eachability analysis are based on owpipe construc- tion, i.e., computing a nite co ver of all trajectories starting from a given set of initial states [ 1 , 5 , 10 – 12 , 15 , 22 ]. Such approaches, while being based on dierent theoretical ideas, conceptually shar e underlying machiner y . Thanks to these common grounds, when implementing a owpipe construction approach (in a prototype or in an end-user tool), our vision is that, instead of starting fr om scratch, one can reuse the code from a related implementation. T o HSCC ’19, A pril 16–18, 2019, Montreal, QC, Canada © 2019 This is the author’s version of the work. It is posted here by permission of A CM for your personal use. Not for redistribution. The denitive version was published in ACM ISBN 978-1-4503-6282-5/19/04. . . $15.00 https://doi.org/10.1145/3302504.3311804 make such reuse convenient, it is essential to provide a exible interface for common operations, supporte d by an ecient library with rich functionality . Given these featur es, researchers and their students can then concentrate on the core algorithms and ideally prototype their approach with relatively few lines of code. In prac- tice, howe ver , many of the existing approaches based on owpipe construction are implemented in dierent tools without any shar ed code base. An important element in this context is the choice of program- ming language. In academic prototyping, the choice is usually a consideration between convenience during development and out- put performance in the nal product. On the one hand, compiled languages such as C++ oer high performance, but the compilation overhead is inconvenient for prototyping. On the other hand, in- terpreted languages such as Python oer an interactive session for convenient prototyping, but these languages fall behind in perfor- mance or need to extend the code to work with another lower-lay er program such as Numba or Cython (known as the two-language problem). A compromise between the two w orlds are just-in-time ( JI T) compiled languages such as MA TLAB . Julia , a recent 1 programming language for high-performance scientic computing [ 6 ], r econciles the two advantages of compiled and interpreted languages describe d above, as it comes with an interactive read-evaluate-print loop (REPL) front-end, but is JI T compiled to achieve performance that is comp etitive with com- piled languages such as C . A distinctive feature of Julia is multiple dispatch (i.e., the function to execute is chosen based on each ar- gument type), which allows to write ecient machine code based on a given type, e.g., of the set. As additional features, Julia is plat- form independent, has an ecient interface to C , is supported in Jupyter noteb ooks and well-suited for parallel computing. Julia has a determined and quickly-growing community , especially for scientic tools. 2 W e believe that all this makes Julia an interesting programming language for writing a tool for reachability analysis. In this paper , we present JuliaReach 3 , an open-source tool- box for rapid prototyping of set-based reachability approaches, written in Julia . The JuliaReach toolbox consists of two main packages: The rst package, Rea chability , oers core infrastruc- ture for implementing reachability algorithms for continuous and hybrid systems. As a proof of concept, the package curr ently pro- vides implementations of an algorithm for the analysis of linear time-invariant (LTI) systems base d on decomp osition [ 8 ], and of 1 V ersion 1.0 of Julia was released in A ugust 2018. 2 https://julialang.org 3 http://juliareach.org HSCC ’19, April 16–18, 2019, Montreal, QC, Canada S. Bogomolov et al. the classic interleaving algorithm of continuous and discrete-post operators as used in hybrid system to ols such as Spa ceEx [ 7 , 9 , 12 ]. The second package, LazySets , is a library for state-of-the-art calculus with convex sets. A s the name indicates, a key featur e of LazySets is lazy set representation using support functions. For making the lazy sets concrete , the library oers means to overap- proximate them, e.g., using template dir ections or ε -close approxi- mation. LazySets also supplies concrete set representations, e.g., polyhedra, zonotopes, and ellipsoids. A binary le to use Julia can be downloaded on the ocial web page, and installing the JuliaRea ch packages is just one command. The toolbox puts attention on high eciency and usability , b oth 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 W ork. Several reachability tools are implemented in C++ . For example, Sp aceEx is a C++ package that curr ently integrates several independent reachability algorithms for linear hybrid sys- tems [ 12 ]. Flow ∗ implements T aylor model approximation for non- linear hybrid systems in C++ [11]. There are also other reachability tools that are written in JIT - compiled or interpreted languages. For e xample, CORA implements several algorithms for linear and nonlinear hybrid systems in MA T - LAB [ 1 ], and also oers a standalone set library focused on zono- topes, as well as a librar y for inter val arithmetic. A subset of the algorithms from CORA have been implemented in C++ in the tools CORA/SX and SymRea ch , leading in some instances to a sp eedup of up to a factor of 3 [ 16 ]. A fe w tools are implemented in Python . For instance, Hylaa implements discrete-time reachability algorithms for linear hybrid systems [5]. HyPro is a C++ library for convex set representation with similar aims as the LazySets library [ 23 ]. It provides a common interface for set representations like boxes, convex polyhedra, support func- tions, and zonotop es. Recently , a reachability tool called HyDRA [ 2 ] 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 Se ction 3. In Se ction 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 fr om convex geometry: A set is closed if it contains all its limit p oints. A set S is convex if for any m points v j ∈ S and m non-negative numbers λ j that sum up to 1 we have that Í j λ j v j ∈ S as well. Alternatively , a closed convex set is an intersection of (possibly innitely many) closed half-spaces. Every conve x 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., oating point ( Float64 ) and exact ( Rational ) precision with no additional per- formance penalty: At runtime, Julia uses multiple dispatch on N and JI T -compiles into type-specic code. The library comes with a collection of common set representa- tions, such as balls in dierent norms, ellipsoids, hyperrectangles, polyhedra, p olytopes (i.e., bounded p olyhedra) in constraint and in vertex r epresentation, and zonotopes. One key featur e of LazySets is to apply common binary set op erations 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 con- vex 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. T o realize the lazy paradigm, a set type should provide an implementation of two functions: The rst function has signature σ ( d , S ) in LazySets and returns some furthest point of a set S in a given direction d , i.e., its support vector . W e recall that the set of supp ort vectors corresponds to the optimal points for the support function ρ , i.e., σ S ( d ) : = { x ∈ S | d T x = ρ S ( d )} , ρ S ( d ) : = max x ∈ S d T x . The second function ρ S ( d ) has signature ρ ( d , S ) in LazySets and is exactly the support function. Since the default implementa- tion of the support function is given by the above formula, it suces to dene the function σ ; but if a more ecient 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 eciently without explicitly representing the set result- ing from the operation [ 18 , 19 ], with the exception of intersection. For intersection, we only provide ( over )approximate values for ρ , where we can either use a line search algorithm inspir ed by [ 13 ] or a coarser heuristics given as ρ X ∩ Y ( d ) ≤ min ( ρ X ( d ) , ρ Y ( d )) [18]. LazySets denes 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 repr esentation and their constraint representation. This allows to write generic (or “virtual”) functions for arbitrar y 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 elds that hold the map and the set, respe ctively . 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 · σ ( M T d , S ) , i.e., it asks for the support vector of the wrapped set S . As a second example, consider two p olytopes P 1 and P 2 . The command P 1 ⊕ P 2 instantiates a new MinkowskiSum instance. Sim- ilar to the linear map, the support vector can b e dened recur- sively as σ P 1 ⊕ P 2 ( d ) = σ P 1 ( d ) ⊕ σ P 2 ( d ) . Again, the binary operation MinkowskiSum is dened between any two LazySet types. On the other hand, this operation can be p erformed explicitly using the command minkowski_sum( P 1 , P 2 ) , which calls an external library JuliaReach: a T oolb ox for Set-Based Reachability HSCC ’19, April 16–18, 2019, Montreal, QC, Canada called Polyhedra.jl [ 20 ],that provides a unied interface to well- known implementations of polyhedral computations, such as CDD or LRS . The approach lazy-by-default and optionally-explicit is also avail- able with other set types. For e xample, we can e xecute the concrete linear_map function for Zonotope s, in which case the generators will be transformed. In practice it is convenient to switch from a lazy to a concrete rep- resentation at some point, ideally without an exponential increase in the computational cost when dealing with high-dimensional problems. Below we describe the Approxima tions sub-module, which serves this purpose. 2.2 From lazy to concrete set representation The Approxima tions module implements the transfer from a lazy set representation to a concrete set, generally involving an over- approximation. 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 follow- ing example. A typical set e quation for a discrete approximation model of LTI systems is Y = CH ( e A δ X 0 ⊕ δ BU , X 0 ) , (1) with real matrix A ∈ R n × n , time step δ , initial states X 0 , nondeter- ministic inputs BU , and CH denoting the convex hull. For example, let δ = 0 . 1 , X 0 ⊆ R n be a ball with center ( 1 , . . . , 1 ) and radius 0 . 1 in the innity norm, U ⊆ 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 di- mensions. With LazySets , Eq. (1) for n = 1 , 000 , m = 2 , and random matrix coecients is n = 1 0 0 0 ; m = 2 ; δ = 0 . 1 A = s p r a n d n ( n , n , 0 . 0 1 ) ; B = r a n d n ( n , m ) X 0 = B a l l I n f ( o n e s ( n ) , 0 . 1 ) U = B a l l 2 ( z e r o s ( m ) , 1 . 2 ) Y = C o n v e x H u l l ( S p a r s e M a t r i x E x p ( A * δ ) * X 0 ⊕ δ * B * U , X 0 ) The execution is instantaneous b ecause 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 Approxima tions module can now be used to get informa- tion about Y without having computed any concrete representation; e.g., projection of a set into two dimensions is useful for visualiza- tion or v erication of properties that involve only two variables. For this purpose, the module oers the function overapproximate( S , ε ) where S is a set and ε is an error bound. The function over- approximates a 2D set by adding supporting directions until the error bound is achieved (measured in terms of the Hausdor dis- tance), where the number of dir ections is optimal. The algorithm implements Kamenev’s method [ 17 , 21 ], 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 dene the 2 × n projection matrix and apply it as a linear map (i.e., fr om the left). Then we use the overapproximate method with a specied value of ε : -1.0 -0.5 0.0 0.5 1.0 -1.0 -0.5 0.0 0.5 1.0 Figure 1: Description of the overappro ximation algorithm. W e start with b ox directions ( ± d 1 and ± d 2 ) and then check for each angle if the distance b etween the support vector of the original set and the vertex is less than an error bound ε . 0.4 0.8 1.2 0.8 1.0 1.2 Figure 2: Decomp osition of the 1 , 000 -dimensional set Y , pro- jected onto x 1 and x 50 , with polytop e approximations of dif- ferent precision: ε = Inf (box directions, blue), ε = 0 . 1 (green), and ε = 0 . 001 (red). Runtimes are 55 ms, 195 ms and 2 s, re- spectively . π = s p z e r o s ( 2 , n ) ; π [ 1 , 1 ] = π [ 2 , 5 0 ] = 1 r e s = o v e r a p p r o x i m a t e ( π * Y , ε ) Fig. 2 shows some results for dierent values of ε . For high-dimensional overappr oximation, instead of ε one can alternatively pass a set type (e.g., Hyperrectangle for box approx- imation) or a special iterator typ e for template directions. W e note again that one can call the overapproximate function with an arbi- trary (nested) lazy set, as long as it implements ρ ( d , S ) to evaluate the support function for a given direction. 3 THE REA CHABILI T Y P A CKA GE 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 [ 4 ]. W e assume that the reader is familiar with the basic terminology of hybrid automata. HSCC ’19, April 16–18, 2019, Montreal, QC, Canada S. Bogomolov et al. 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-dened options as key-value pairs that can be passed around to the solver backends. Dynamical systems are passed to Reachability as special wrap- per types. The reachability algorithms themselves are interfaced with post operators, which we will describe later . For hybrid sys- tems, it is common to use clustering algorithms and xe d-point checks. Reachability oers standard convex-hull clustering and checks for a xe d-point (i.e., inclusion) after every jump, either before or after clustering (controlled by the user ). W e support p ost-processing the resulting sets with an output function or a projection to giv en output dimensions. The sets can be visualized with the external Plots.jl to ol using a single command. 3.2 Continuous-Post Operators The interface for continuous-post op erators is rather simple. Based on an :algorithm option, we just call the respective interface func- tion with the continuous mo del, a ReachSet object, and a list of options (the number of steps, the step size, and other algorithm- specic settings). A ReachSet is a wrapper of a set and a correspond- ing 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 ReachSet s, sometimes called a reach tube. For out-of-the-box usage , we provide an ecient implementation of a recent approach for the r eachability analysis of LTI systems [ 2 , 8 ]. The idea is to rst decomp ose the system into low-dimensional blocks and then solve many small reachability problems. While the decomposition comes with an appro ximation 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 discr ete tran- sitions. Given a reach tub e X , a transition with a guard G and assign- ment as д n (·) , and a target invariant I , the discrete-post operator should compute (an ov erapproximation of ) the set as д n (X ∩ G ) ∩ I and pass this set on to the continuous-post operator again. Having both a continuous and a discrete-post op erator , one can implement a classic algorithm that essentially interleaves the two post operators. Reachability oers the function “ solve ” , which implements this algorithm for given post operators. W e implemented two discrete-post operators to instantiate the above algorithm. The rst operator uses (concrete, i.e., non-lazy ) polytopes in constraint representation and relies on the external library Pol yhedra.jl to perform the ab ove-mentioned operations. In particular , the interse ctions are computed explicitly . The se cond discrete-post operator is lazy but oers the options to overappro ximate 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 -0.3 0.3 0.9 -0.25 0.25 0.75 Figure 3: Reachability plot of a two-mode hybrid system. The colors represent location 1 (blue) and location 2 (red), respectively . support vector in one specic 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, w e compare the performance of these tw o operators in a case study . W e will instantiate the lazy discrete-post operator once with both intersections kept lazy and once with both intersections overapproximated. 4 CASE ST UDY W e perform two case studies that investigate some applications. In the rst one, we r e-implement a classic approach for reachabil- ity analysis of LTI systems using the LazySets librar y . W e 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 dierent discrete-post operators outlined in the previous section. 4.1 W orking with the LazySets Library W e wrote a full-edged re-implementation of Girard’s reachability algorithm for LTI systems with zonotopes [ 14 ] in Julia . The goal is to demonstrate prototyping with LazySets and the general appli- cability as a standalone library . The complete Julia sour ce code is given in Algorithm 1, where we underlined those parts that inv olve the LazySets library . W e then extended the example presented in [ 14 , Section 4] to a hybrid system of two modes ℓ i , i = 1 , 2 , with initial states [ 0 . 9 , 1 . 1 ] × [− 0 . 1 , 0 . 1 ] and uncertain inputs from a set u with µ = ∥ u ∥ ∞ = 0 . 001 . The dynamics matrices A i are dened as follows: A 1 =  − 1 − 4 4 − 1  A 2 =  1 4 − 4 − 1  W e add a transition t i from mode ℓ i to ℓ 3 − i with a hyp erplane guard д i : д 1 ≜ x 1 = − 0 . 5 , д 2 ≜ x 2 = − 0 . 3 . LazySets oers order reduction for zonotopes, which we used here with an upp er bound of 10 generators. W e implemente d a custom reachability algorithm JuliaReach: a T oolb ox for Set-Based Reachability HSCC ’19, April 16–18, 2019, Montreal, QC, Canada Algorithm 1: Simple LTI reachability algorithm using zono- topes in Julia . 1 function ReachContinuous(A, X0 , δ , µ , T , max_order) # bloating factors 2 Anorm = norm(A, Inf ) 3 α = (exp( δ ∗ Anorm) − 1 − δ ∗ Anorm)/norm( X0 , Inf ) 4 β = (exp( δ ∗ Anorm) − 1) ∗ µ /Anorm # discretized system 5 n = size(A, 1) 6 ϕ = exp( δ ∗ A) 7 N = oor(Int, T/ δ ) 8 R = V ector{ LazySet }(N) 9 if N == 0 10 return R 11 end # initial reach set in time interval [0, δ ] 12 ϕ p = (I+ ϕ )/2 13 ϕ m = (I − ϕ )/2 14 gens = hcat( ϕ p ∗ X0 .generators, ϕ m ∗ X0 .center , ϕ m ∗ X0 .generators) 15 R[1] = minkowski_sum ( Zonotope ( ϕ p ∗ X0 .center , gens), Zonotope (zeros(n), ( α + β ) ∗ eye(n))) 16 if order(R[1]) > max_order 17 R[1] = reduce_order (R[1], max_order) 18 end # recurrence for [ δ , 2 δ ], . . . , [(N − 1) δ , N δ ] 19 ball β = Zonotop e (zeros(n), β ∗ e ye(n)) 20 for i in 2:N 21 R[i] = minkowski_sum ( linear_map ( ϕ , R[i − 1]), ball β ) 22 if order(R[1]) > max_order 23 R[i] = reduce_order (R[i], max_order) 24 end 25 end 26 return R 27 end for this case study , again only using the LazySets capabilities. For simplicity , a transition is taken as so on as an intersection with the guard is dete cted. The complete source code is given in App en- dix A.W e plot the reachable states for the time interval [ 0 , 4 ] and time step δ = 0 . 001 in Fig. 3. The analysis takes 0.25 seconds. 4.2 W orking with the Reachability Framework In the second case study , we show how the Rea chability frame- work can be use d to quickly evaluate a new approach to reachability analysis. Consider the case wher e a researcher came up with a new strategy for the discrete-post operator and now wants to compare to existing approaches. The researcher w ould just need to imple- ment the dierent approaches as dierent post operators in the Reachability package and then evaluate the performance. Here Û x = − 2 x + 1 . 4 Û y = − y − 0 . 7 Û x 1 = 5 x − 5 x 1 Û x 2 = 5 x 1 − 5 x 2 Û x 3 = 5 x 2 − 5 x 3 . . . Û x m = 5 x m − 1 − 5 x m x ≤ 0 y + 0 . 714286 ∗ x ≥ 0 Û x = − 2 x − 1 . 4 Û y = − y + 0 . 7 Û x 1 = 5 x − 5 x 1 Û x 2 = 5 x 1 − 5 x 2 Û x 3 = 5 x 2 − 5 x 3 . . . Û x m = 5 x m − 1 − 5 x m x ≤ 0 y + 0 . 714286 ∗ x ≤ 0 Û x = − 2 x + 1 . 4 Û y = − y − 0 . 7 Û x 1 = 5 x − 5 x 1 Û x 2 = 5 x 1 − 5 x 2 Û x 3 = 5 x 2 − 5 x 3 . . . Û x m = 5 x m − 1 − 5 x m x ≥ 0 y + 0 . 714286 ∗ x ≥ 0 Û x = − 2 x − 1 . 4 Û y = − y + 0 . 7 Û x 1 = 5 x − 5 x 1 Û x 2 = 5 x 1 − 5 x 2 Û x 3 = 5 x 2 − 5 x 3 . . . Û x m = 5 x m − 1 − 5 x m x ≥ 0 y + 0 . 714286 ∗ x ≤ 0 x ≤ 0 y + 0 . 714286 ∗ x = 0 x = 0 y + 0 . 714286 ∗ x ≤ 0 x ≥ 0 y + 0 . 714286 ∗ x = 0 x = 0 y + 0 . 714286 ∗ x ≥ 0 Figure 4: Hybrid automaton model of a 2-dimensional oscil- lator and an m -dimensional lter . we simulate this situation by using the operators that currently exist in the Reachability package (see Section 3.3). W e compare the p erformance of the dierent discrete-post opera- tors using a model of a lter ed oscillator [ 12 ]. The model represents a parameterized hybrid system consisting of ( i ) a two-dimensional switched oscillator in the variables x and y , and ( i i ) a lter with m state variables x 1 , . . . , x m . W e show the hybrid automaton in Fig. 4. The lter smooths x with x m as an output signal, and the amplitude decreases with increasing dimension of the lter . The initial set is given as x ∈ [ 0 . 2 , 0 . 3 ] , y ∈ [− 0 . 1 , 0 . 1 ] , and x j = 0 for all x j . As a r eference tool, we use Sp aceEx (in version 0.9.8f ) to evaluate the overall performance. T o compare only the performance of the discrete-post operators, we x the number of jumps to four , which corresponds to one loop in the automaton. Reachability oers a simple option for that purpose, but SpaceEx does not. T o have a fair comparison, we do not use our option and instead modify the model in the following way . W e add a new variable b (for bound ) with initial value 1. On the rst 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 ≤ 2 . W e assume the following benchmark settings. The time interval is [ 0 , 20 ] for lters of dimension up to four , and [ 0 , 99 ] for higher dimensions. The time step is δ = 0 . 01 for b oth Reachability and Spa ceEx . W e ran the benchmarks on a computer with a 2.20 GHz CP U and 8 GB RAM. T able 1 outlines the performance of the dif- ferent discrete-post operators and Spa ceEx . For Spa ceEx , we used the support-function algorithm (“ supp ”). Fig. 5 shows some plots of the output variables for the lazy discrete-post operator with over- approximation. (W e use d box approximations in this evaluation.) HSCC ’19, April 16–18, 2019, Montreal, QC, Canada S. Bogomolov et al. T able 1: Runtimes (in seconds) for dierent discrete-post operators in the Reachability package and the SpaceEx tool. Number of lters 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 Spa ceEx 0 . 11 0 . 22 0 . 11 0 . 30 1 . 10 6 . 78 56 . 72 262 . 88 1 , 011 . 17 -0.5 0.0 0.5 -0.5 0.0 0.5 -0.5 0.0 0.5 0.0 0.2 0.4 Figure 5: Output for the ltered oscillator model with dier- ent lter dimension m , using Reachability with the lazy discrete-post operator with overapproximation. Left: m = 4 . Right: m = 16 . As expecte d, the concrete intersection with polytopes is the slowest. (For this model, all computations could be performe d in constraint representation.) The lazy-intersection approaches ( sec- ond 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 W e have presented the JuliaReach toolbox, a new environment for developing reachability algorithms for dynamical systems. Ju- liaReach is written in the mo dern programming language Julia that unies rapid prototyping and high performance. The toolbox consists of a framework for reachability analysis ( Rea chability ) that is built around a standalone library for calculus with conv ex sets ( LazySets ). In the future, we will extend the amount of p ost op erators of- fered 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. A CKNO WLEDGMEN TS This work was partially supported by the Air Force Oce of Sci- entic Research under award number F A2386-17-1-4065, by the ARC project DP140104219 (Robust AI P lanning for Hybrid Sys- tems), by the Austrian Science Fund (F WF) under grants S11402-N23 (RiSE/SHiNE) and Z211-N23 (Wittgenstein A ward), and the Euro- pean Union’s Horizon 2020 resear ch and innovation programme un- der the Marie Skłodowska-Curie grant agreement No 754411. Any opinions, ndings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reect the views of the United States Air Force. REFERENCES [1] Matthias Altho. 2015. An Introduction to CORA 2015. In ARCH@CPSW eek . 120–151. [2] Matthias Altho, Stanley Bak, Xin Chen, Chuchu Fan, Marcelo Forets, Goran Frehse, Niklas Kochdumper , Y angge Li, Sayan Mitra, Rajarshi Ray, Christian Schilling, and Stefan Schupp. 2018. ARCH-COMP18 Category Report: Continuous and Hybrid Systems with Linear Continuous Dynamics. In ARCH18 . 23–52. [3] Rajeev Alur, Costas Courcoubetis, Nicolas Halbwachs, Thomas A. Henzinger, Pei-Hsin Ho, Xavier Nicollin, Alfredo Olivero , Joseph Sifakis, and Sergio Y ovine. 1995. The Algorithmic Analysis of Hybrid Systems. Theor . Comput. Sci. 138, 1 (1995), 3–34. [4] Rajeev Alur , Costas Courcoubetis, Thomas A. Henzinger , and Pei-Hsin Ho. 1992. Hybrid A utomata: An Algorithmic Approach to the Specication and V erication of Hybrid Systems. In Hybrid Systems , V ol. 736. 209–229. [5] Stanley Bak and Parasara Sridhar Duggirala. 2017. Simulation-Equivalent Reach- ability of Large Linear Systems with Inputs. In CA V . 401–420. [6] Je Bezanson, Alan Edelman, Stefan Karpinski, and Viral B. Shah. 2017. Julia: A Fresh Approach to Numerical Computing. SIAM Rev . 59, 1 (2017), 65–98. [7] Sergiy Bogomolov , Alexandre Donzé, Goran Frehse, Radu Grosu, T aylor T . John- son, Hamed Ladan, Andreas Podelski, and Martin W ehrle. 2015. Guided search for hybrid systems based on coarse-grained space abstractions. ST T T (2015). [8] Sergiy Bogomolov , Mar celo Forets, Goran Frehse, Frédéric Viry , Andreas Podelski, and Christian Schilling. 2018. Reach Set Approximation through Decomposition with Low-dimensional Sets and High-dimensional Matrices. In HSCC . 41–50. [9] Sergiy Bogomolov, Goran Frehse, Marius Greitschus, Radu Grosu, Corina S. Pasareanu, Andreas Podelski, and Thomas Strump. 2014. Assume-Guarantee Abstraction Renement Meets Hybrid Systems. In HVC . 116–131. [10] Sergiy Bogomolov , Goran Frehse, Radu Grosu, Hamed Ladan, Andreas Podelski, and Martin W ehrle. 2012. A Box-based Distance between Regions for Guiding the Reachability Analysis of SpaceEx. In CA V . 479–494. [11] Xin Chen, Erika Ábrahám, and Sriram Sankaranarayanan. 2013. Flow*: An Analyzer for Non-linear Hybrid Systems. In CA V . 258–263. [12] Goran Frehse, Colas Le Guernic, Alexandre Donzé, Scott Cotton, Rajarshi Ray , Olivier Lebeltel, Rodolfo Ripado, Antoine Girard, Thao Dang, and Oded Maler . 2011. SpaceEx: Scalable V erication of Hybrid Systems. In CA V . 379–395. [13] Goran Frehse and Rajarshi Ray . 2012. Flowpipe-Guard Intersection for Reacha- bility Computations with Support Functions. In ADHS , V ol. 45. 94–101. [14] Antoine Girar d. 2005. Reachability of Uncertain Linear Systems Using Zonotopes. In HSCC . 291–305. [15] Colas Le Guernic and Antoine Girard. 2009. Reachability Analysis of Hybrid Systems Using Support Functions. In CA V . 540–554. [16] Fabian Immler , Matthias Altho, Xin Chen, Chuchu Fan, Goran Fr ehse, Niklas Kochdumper , Y angge Li, Sayan Mitra, Mahendra Singh T omar, and Majid Zamani. 2018. ARCH-COMP18 Categor y Report: Continuous and Hybrid Systems with Nonlinear Dynamics. In ARCH18 . 53–70. [17] George K. Kamenev . 1996. An algorithm for approximating polyhedra. Computa- tional Mathematics and Mathematical Physics 36, 4 (1996), 533–544. [18] Colas Le Guernic. 2009. Reachability analysis of hybrid systems with linear continuous dynamics . Ph.D. Dissertation. Univ ersité Grenoble 1 - Joseph Fourier . [19] Colas Le Guernic and Antoine Girard. 2010. Reachability analysis of linear systems using support functions. Nonlinear Analysis: Hybrid Systems 4, 2 (2010), 250 – 262. IF AC W orld Congress 2008. [20] Benoît Legat, Robin Deits, Gustavo Goretkin, T wan Koolen, Joey Huchette, Daisuke O yama, Robert Schwarz, Elliot Saba, and Chase Coleman. 2018. Ju- liaPolyhedra/Polyhedra.jl v0.4.4. https://doi.org/10.5281/zenodo.1214290 [21] Alexander V . Lotov and Alexis I. Pospelov . 2008. The modied method of re- ned bounds for polyhe dral approximation of convex polytopes. Computational Mathematics and Mathematical Physics 48, 6 (2008), 933–941. [22] Rajarshi Ray , Amit Gurung, Binayak Das, Ezio Bartocci, Sergiy Bogomolov , and Radu Grosu. 2015. XSp eed: Accelerating reachability analysis on multi-core processors. In HVC . 3–18. [23] Stefan Schupp, Erika Ábrahám, Ibtissem Ben Makhlouf, and Stefan Kowalewski. 2017. HyPro: A C++ Library of State Set Representations for Hybrid Systems Reachability Analysis. In NFM . 288–294. JuliaReach: a T oolb ox for Set-Based Reachability HSCC ’19, April 16–18, 2019, Montreal, QC, Canada A H YBRID REA CHABILI T Y IMPLEMEN T A TION The following Julia function is an emb edding of Algorithm 1 in a simple algorithm for hybrid systems. Apart from LazySets , no other library is r equired. W e underlined the places where LazySets is used. Algorithm 2: Simple hybrid reachability algorithm in Julia . 1 function ReachHybrid(As, T s, init, δ , µ , T , max_order) # start with initial states at time 0 2 queue = [(init[1], init[2], 0.)] 3 res = T uple{ LazySet , Int}[ ] 4 while !isempty(queue) do 5 init, loc, t = pop!(queue) # compute continuous successors 6 R = ReachContinuous(A s[loc], init, δ , µ , T -t, max_order) 7 found_transition = false 8 for i in 1:length(R)-1 9 S = R[i] 10 push!(res, (S, loc)) # check intersection with guards 11 for (guard, tgt_loc) in T s[loc] 12 if ! isdisjoint (S, guard) # nonempty intersection with a guard 13 new_t = t + δ ∗ i 14 push!(queue, (S, tgt_loc, new_t)) 15 found_transition = true 16 end 17 end # stop for first intersection 18 if found_transition 19 br eak 20 end 21 end 22 if !found_transition && length(R) > 0 23 push!(res, (R[end], loc)) 24 end 25 end 26 return res 27 end