Piecewise Robust Barrier Tubes for Nonlinear Hybrid Systems with Uncertainty

Piecewise Robust Barrier T ub es for Nonlinear Hybrid Systems with Uncertain t y Hui K ong 1 , Ezio Barto cci 3 , Y u Jiang 4 , Thomas A. Henzinger 2 1 Max-Planc k-Institute for Soft ware Systems, Kaiserslautern, Germany 2 IST Austria, Klosterneuburg, Austria 3 TU Wien, Vienna, Austria 4 T singhua Univ ersity , Beijing, China Abstract. Piecewise Barrier T ub es (PBT) is a new technique for flo w- pip e o verappro ximation for nonlinear systems with p olynomial dynam- ics, which leverages a combination of barrier certificates. PBT has ad- v an tages ov er traditional time-step based metho ds in dealing with those nonlinear dynamical systems in which there is a large difference in sp eed b et ween tra jectories, pro ducing an ov erapproximation that is time inde- p enden t. How ever, the existing approach for PBT is not efficient due to the application of interv al metho ds for enclosure-b ox computation, and it can only deal with contin uous dynamical systems without uncertaint y . In this paper, we extend the approach with the abilit y to handle b oth con tinuous and h ybrid dynamical systems with uncertain ty that can re- side in parameters and/or noise. W e also improv e the efficiency of the metho d significan tly , by a voiding the use of in terv al-based methods for the enclosure-b ox computation without lo osing soundness. W e hav e de- v elop ed a C++ prototype implementing the prop osed approach and we ev aluate it on sev eral b enc hmarks. The exp erimen ts show that our ap- proac h is more efficient and precise than other metho ds in the literature. 1 In tro duction Hybrid systems (HS) [21] are a suitable mathematical framework to mo del dy- namical systems with b oth discrete and contin uous dynamics. This formalism has b een successfully adopted to design cyb er-ph ysical systems (CPS) whose b e- ha vior is characterized b y an embedded soft ware monitoring and/or con trolling a physical substratum. F ormal verification of HS has indeed a practical impact in engineering by assuring imp ortan t safety-critical requiremen ts at design-time. Despite the great effort to adv ance the state-of-the-art, reachabilit y analy- sis of HS remains one of the most c hallenging verification tasks. Although the problem of reac hability analysis is in general undecidable [21] for HS, in the last decade sev eral efficient and scalable semidecidable approaches hav e b een prop osed to analyse HS with linear dynamics [14,15,34,39,19]. This researc h was supp orted in part b y the Austrian Science F und (FWF) under gran ts S11402-N23, S11405-N23 (RiSE/SHiNE), ADynNet (P28182), and Z211-N23 (Wittgenstein A ward) and the Deutsc he F orsc hungsgemeinsc haft pro ject 389792660- TRR 248. 2 Hui Kong, Ezio Barto cci, Y u Jiang, Thomas A. Henzinger HS with nonlinear ordinary differential equations (ODEs) remains still v ery c hallenging to solve b ecause these ODEs do not hav e a closed form solution in general. One common strategy to tackle this problem is to compute an o ver- appro ximation (also called flowpipe) that contains all the p ossible tra jectories originating from an initial set of states within a b ounded-time horizon [6,1,11,9,10]. If the o verappro ximation do es not intersect with the unsafe set of states, then the system is safe. How ever, if the ov erapproximation is to o coarse, it ma y in- tersect the unsafe set of states only due to the approximation errors and then the verdict ab out safety may b e inconclusiv e. Th us, one of the main problem to address is how to efficiently c ompute tight over-appr oximations of the r e achable set of states for nonline ar c ontinuous and hybrid systems . T o o vercome this problem, in a recent pap er [22], w e hav e in tro duced the notion of Piecewise Barrier T ub es (PBT), a new flowpipe ov erappro ximation for nonlinear systems with p olynomial dynamics. The main idea of this approach is that for each segmen t of a flowpipe, it constructs a coarse b o x that is big enough to contain the segment and then it computes in the b o x a set of barrier functions [32,24] whic h work together to form a tub e surrounding the flowpipe. PBT has adv an tages o ver traditional time-step based metho ds in dealing with those nonlinear dynamical systems in whic h there is a large difference in speed b et w een tra jectories, pro ducing a tight o ver-appro ximation that is time inde- p enden t. Ho wev er, the approach in [22] cannot handle uncertaint y and h ybrid systems. In addition, the use of interv al metho d for enclosure-box computation reduces its efficiency . In this paper, we extend the approach with the ability to handle b oth con- tin uous and h ybrid dynamical systems with uncertaint y which can reside in pa- rameters and/or noise. W e improv e the efficiency of the metho d significan tly , by a voiding the use of interv al metho d for enclosure-b ox computation without lo os- ing soundness. W e hav e developed a C++ prototype implementing the prop osed approac h and we ev aluate it on sev eral b enc hmarks. The exp erimen ts sho w that our approac h is more efficient and precise than other metho ds proposed in the literature. The other existing techniques used to compute a b ounded flo wpip e are mainly based on in terv al metho d [30] or T aylor model [4]. Interv al metho d is quite ef- ficien t even for high dimensional systems [30], but it suffers from the wr apping effe ct that arises due to an uncon trollable growth of the interv al enclosure that accum ulates ov erapproximation errors. The use of T aylor mo del is more pre- cise b ecause it uses a v ector of p olynomials plus a vector of small interv als to sym b olically represent the flowpipe. Ho wev er, c hecking the intersection with the unsafe region requires generally the use of in terv al metho d that brings back the wrapping effect. In particular, the wrapping effect can explo de easily when the flowpipe segment ov er a time in terv al is stretched drastically due to a large difference in sp eed b et ween individual tra jectories. Only recen tly , to ols such as CLR T [9,10], Flo w* [6], MathSA T SMT solv er [8,7], HySA T/iSA T [12], dReach [25], C2E2 [11] and CORA [1], hav e made some pro- gresses in verifying nonlinear contin uous and hybrid mo dels. Some of these to ols [7,25,12] are based on decision pro cedures that o vercome the theoretical limits in nonlinear the ories o ver the reals. The main idea is to encode the reacha- bilit y problem for nonlinear systems as first-order logic formulas ov er the real n umbers. A satisfiability mo dulo theories (SMT) solv er implementing such pro- cedures can return either a verdict of unsatisfiabilit y when the unsafe region is Title Suppressed Due to Excessive Length 3 not reached or an inconclusive verdict [25,12] such as δ -sat if the problem is satisfiable given a certain precision δ (the same problem may result unsatisfiable b y increasing the precision). How ev er, in the case of unsatisfiabilit y these to ols generally do not provide a reachable set represen tation that explains the v erdict. Other techniques for reachabilit y analysis of nonlinear systems include inv ariant generation [29,40,23,36,37], abstraction and h ybridization [35,26,2,31,16,5]. The paper is organized as follows. Section 2 presents the necessary prelimi- naries. Section 3 sho ws how to compute robust barrier certificates using linear programming, while in Section 4 we present our approach to address the reach- abilit y analysis problem of nonlinear contin uous and hybrid systems with uncer- tain ty . Section 5 provides our exp erimen tal results and we conclude in Section 6. 2 Preliminaries In this section, we recall some concepts used throughout the pap er. W e first clarify some notation conv entions. If not sp ecified otherwise, we use b oldface lo wer case letters to denote v ectors, w e use R for the real num b ers field and N for the set of natural num b ers, and w e consider m ultiv ariate p olynomials in R [ x ] , where the comp onen ts of x act as indeterminates. In addition, for all the polynomials B ( c , x ) , w e denote b y c the vector comp osed of all the c i and denote by x the vector comp osed of all the remaining v ariables x i that o ccur in the p olynomial. W e use R ≥ 0 and R > 0 to denote the domain of nonnegative real n umber and positive real n umber resp ectively . With an abuse of notation, w e sometimes use B ( x ) = 0 for the semialgebraic set it defines. ∂ S denotes the b oundary of compact set S . Next, we present the notation of the Lie deriv ative, which is widely used in the discipline of differential geometry . Let f : R n → R n b e a contin uous v ector field suc h that ˙ x i = f i ( x ) where ˙ x i is the time deriv ativ e of x i ( t ) . Definition 1 (Lie deriv ative). F or a given p olynomial p ∈ R [ x ] over x = ( x 1 , . . . , x n ) and a c ontinuous system ˙ x = f , wher e f = ( f 1 , . . . , f n ) , t he Lie derivative of p ∈ R [ x ] along f is define d as L f p = P n i =1 ∂ p ∂ x i · f i . Essen tially , the Lie deriv ative of p is the time deriv ative of p , i.e., reflects the c hange of p ov er time. In this pap er, w e fo cus on semialgebraic systems with uncertain ty , which is describ ed by the following ODE. ˙ x = f ( x ( t ) , u ( t )) (1) where f is a vector of p olynomial functions, x ( t ) is a solution of the system, u ( t ) is the v ector of uncertain parameters and/or p erturbation and u ( t ) is Lipschitz con tinuous. Note that we do not make a distinction b et ween uncertain parame- ters and p erturbation since we deal with them uniformly . F ormally , semialgebraic system with uncertain ty is defined as follows. Definition 2 (Semialgebraic system with uncertain ty). A semialgebr aic system with unc ertainty is a 5 -tuple M def = h X, f , X 0 , I , U i , wher e 1. X ⊆ R n is the state sp ac e of the system M , 4 Hui Kong, Ezio Barto cci, Y u Jiang, Thomas A. Henzinger 2. f ∈ R [ x , u ] n is lo c al ly Lipschitz c ontinuous ve ctor function defining the ve c- tor flow as in ODE (1) , 3. X 0 ⊆ X is the initial set, which is semialgebr aic [41], 4. I is the invariant or domain of the system, 5. U is a domain for the unc ertain p ar ameters and p erturb ation, i.e., u ( t ) ∈ U The lo cal Lipschitz contin uity guarantees the existence and uniqueness of the differen tial equation ˙ x = f lo cally . A tra jectory of a semialgebraic system with uncertain ty is defined as follows. Definition 3 (T ra jectory). Given a semialgebr aic system with unc ertainty M , a tr aje ctory originating fr om a p oint x 0 ∈ X 0 to time T > 0 is a c ontinu- ous and differ entiable function ζ ( x 0 , t ) : [0 , T ) → R n such that 1) ζ ( x 0 , 0) = x 0 , and 2) ∃ u ( · ) : ∀ τ ∈ [0 , T ) : d ζ dt   t = τ = f ( ζ ( x 0 , τ ) , u ( τ )) , wher e u ( · ) : [0 , T ) → U . T is assume d to b e within the maximal interval of existenc e of the solution fr om x 0 . F or ease of readability , w e also use ζ ( t ) for ζ ( x 0 , t ) if it is clear from the con text. Definition 4 (Safety). Given an unsafe set X us ⊆ X , a semialgebr aic system with unc ertainty M is said to b e safe if no tr aje ctory ζ ( x 0 , t ) of M satisfies that ∃ τ ∈ R ≥ 0 : ζ ( x 0 , τ ) ∈ X us , wher e x 0 ∈ X 0 . 3 Robust Barrier Certificate b y Linear Programming A barrier certificate for a contin uous dynamics system is a real-v alued function B ( x ) such that 1) the initial set and the unsafe set are lo cated on differen t sides of the hyper-surface H = { x ∈ R n | B ( x ) = 0 } resp ectiv ely , and 2) no tra- jectory originating from the same side of H as the initial set can cross through H to reach the other side. Therefore, the existence of suc h a function B ( x ) can guarantee the safet y of the system. The ab o ve condition can b e formalized using an infinite sequence of higher order Lie deriv atives [42]. Unfortunately , this formalization cannot b e applied directly to barrier certificate computation. Therefore, a couple of sufficient conditions for the ab o ve condition hav e b een prop osed [24,28,17]. Most recently , based on the sufficient condition in [32], a new approach was prop osed to o verappro ximate the flowpipe of nonlinear contin- uous dynamical systems using combination of barrier certificates [22]. How ev er, the approach is limited to contin uous dynamical systems without uncertaint y . T o tackle this problem, we extend the approach to deal with contin uous and h ybrid systems with uncertaint y . Similarly , we adopt the same barrier certificate condition as [22], but we introduce the uncertain ty in the barrier certificate con- dition. Note that in order to distinguish it from barrier certificate for dynamical system without uncertaint y , we call a barrier certificate satisfying the follo wing condition r obust b arrier c ertific ate . Theorem 1. Given an unc ertain semialgebr aic system M = h X , f , X 0 , I , U i , let X us b e the unsafe set, the system is guar ante e d to b e safe if ther e exists a r e al-value d function B ( x ) such that ∀ x ∈ X 0 : B ( x ) > 0 (2) ∀ ( x , u ) ∈ I × U : L f B > 0 (3) ∀ x ∈ X us : B ( x ) < 0 (4) Title Suppressed Due to Excessive Length 5 The most common approac h to barrier certificate computation is b y SOS programming [32,24]. The idea of this kind of approach is to first relax the origi- nal constraints like (2)–(4) into a set of p ositive semidefinite (PSD) p olynomials b y applying Putinar representation [33], whic h is further relaxed by requiring ev ery PSD p olynomial has a sum-of-squares decomp osition, which can b e solv ed b y SOS programming in p olynomial time. How ever, constructing automatically a set of consistent templates for the barrier certificate as well as the auxiliary p olynomials is not trivial. In addition, SOS programming metho d can yield fake solution sometimes due to n umerical error. An alternative to SOS programming based approaches is to use linear pro- gramming based approaches. This class of approaches relies on an LP-relaxation to the original constraint. In [38,3], to compute Ly apunov function, an LP- relaxation was obtained by applying Handelman represen tation to the original constrain t. Recen tly , this kind of LP-relaxation w as adopted in [22] to com- pute piecewise barrier tub es. In [43], an extended version of Handelman rep- resen tation, called Krivine representation [27], was employ ed for barrier certifi- cate computation. Compared to Handelman representation, which can only deal with conv ex p olytopes, Krivine representation can deal with more general com- pact semialgebraic sets. How ever, Krivine representation requires normalizing the p olynomials inv olved, which is exp ensiv e. In this pap er, we adopt the same represen tation as in [22], i.e., Handelman represen tation as our LP-relaxation scheme for Theorem 1. W e assume that the initial set X 0 , the unsafe set X us , the in v ariant I , the parameter and/or p ertur- bance space are all conv ex and compact p olyhedra, i.e., X 0 = { x ∈ R n | p 1 ( x ) ≥ 0 , · · · , p m 1 ( x ) ≥ 0 } , I = { x ∈ R n | q 1 ( x ) ≥ 0 , · · · , q m 2 ( x ) ≥ 0 } , U = { u ∈ R l | w 1 ( u ) ≥ 0 , · · · , w m 3 ( u ) ≥ 0 } and X us = { x ∈ R n | r 1 ( x ) ≥ 0 , · · · , r m 4 ( x ) ≥ 0 } where p i ( x ) , q i ( x ) , r k ( x ) and w i ( u ) , are all linear p olynomials. Then, Theorem 1 can b e relaxed as follows. Theorem 2. Given a semialgebr aic system with unc ertainty M = h X , f , X 0 , I , U i , let X 0 , X us , I and U b e define d as ab ove, the system is guar ante e d to b e safe if ther e exists a r e al-value d p olynomial function B ( x ) such that B ( x ) ≡ X | α |≤ M 1 λ α m 1 Y i =1 p α i i +  1 (5) L f B ≡ X | β |≤ M 2 λ β m 2 Y i =1 q β i i m 3 Y j =1 w β m 2 + j j +  2 (6) − B ( x ) ≡ X | γ |≤ M 3 λ γ m 4 Y i =1 r γ i i +  3 (7) wher e α = ( α k ) , β = ( β k ) , γ = ( γ k ) , λ α , λ β , λ γ ∈ R ≥ 0 ,  i ∈ R > 0 and M i ∈ N , i = 1 , · · · , 3 . R emark 1. Theorem 2 implies that the system M can b e prov ed to b e safe as long as we can find a real-v alued polynomial function B ( x ) such that B ( x ) , − B ( x ) and L f B can be written as a nonnegativ e com bination of the prod- ucts of the p o wers of the p olynomials defining X 0 , X us and I × U resp ectiv ely . 6 Hui Kong, Ezio Barto cci, Y u Jiang, Thomas A. Henzinger This theorem provides us with a solution to solv e barrier certificate by linear programming. Given a p olynomial template B ( c , x ) for B ( x ) , where c is the co efficien ts of the monomials to b e decided in B ( c , x ) , we substitute B ( c , x ) for B ( x ) o ccurring in the conditions (5)–(7) to obtain three p olynomial identities in R [ x ] with linear p olynomials in R [ c , λ ] as their coefficients, where λ is a vector comp osed of all the λ α , λ β , λ γ o ccurring in (5)–(7). Since (5)–(7) are identities, then all the co efficien ts of the corresp onding monomials on b oth sides of the iden tities must be identical. By collecting the corresp onding co efficien ts of the monomials on b oth sides of the identities and let them equal resp ectively , w e obtain a system S of linear equations and inequalities on c , λ . Now, finding a robust barrier certificate is conv erted to finding a feasible solution for S , which can b e solv ed by linear programming efficiently . Since the degree of B ( c , x ) is k ey to the expressiv e p ow er of B ( c , x ) , in our implemen tation, w e attempt to solv e a barrier certificate from a group of templates with different degrees. Due to the page limit, we do not elab orate on our algorithm for barrier cer- tificate computation, but we demonstrate how it w orks in the follo wing example. Example 1. Given a 2D system defined by ˙ x = 2 x + 3 y + u 1 , ˙ y = − 4 x + 2 y + u 2 , let X 0 = { ( x, y ) ∈ R 2 | p 1 = x + 100 ≥ 0 , p 2 = − 90 − x ≥ 0 , p 3 = y + 45 ≥ 0 , p 4 = − 40 − y ≥ 0 } , I = { ( x, y ) ∈ R 2 | q 1 = x + 110 ≥ 0 , q 2 = − 80 − x ≥ 0 , q 3 = y + 45 ≥ 0 , q 4 = − 20 − y ≥ 0 } , U = { ( u 1 , u 2 ) ∈ R 2 | w 1 = u 1 + 50 . 0 ≥ 0 , w 2 = 50 . 0 − u 1 ≥ 0 , w 3 = u 2 + 50 . 0 ≥ 0 , w 4 = 50 . 0 − u 2 ≥ 0 } and X us = { ( x, y ) ∈ R 2 | r 1 = x + 98 ≥ 0 , r 2 = − 90 − x ≥ 0 , r 3 = y + 24 ≥ 0 , r 4 = − 20 − y ≥ 0 } . Assume B ( c , x ) = c 0 + c 1 x + c 2 y , M i =  i = 1 for i = 1 , · · · , 3 , then w e obtain the follo wing p olynomial identities according to Theorem 2 c 1 + c 2 x + c 3 y − 4 X i =1 λ 1 i p i −  1 ≡ 0 c 2 (2 x + 3 y + u 1 ) + c 3 ( − 4 x + 2 y + u 2 ) − 4 X j =1 λ 2 j q j − 4 X j =1 λ 3 j w j −  2 ≡ 0 − ( c 1 + c 2 x + c 3 y ) − 4 X k =1 λ 4 k r k −  3 ≡ 0 where λ ij ≥ 0 for i, j = 1 , · · · , 4 . If we collect the co efficien ts of x, y, u 1 , u 2 in the ab o v e p olynomials and let them b e 0 , we obtain a system S of linear p olynomial equations and inequalities o ver c i , λ ij . By solving S using linear programming, w e obtain a feasible solution with c 1 = − 1263 . 5 , c 2 = − 11 . 5 , c 3 = − 5 . 85 . 4 Piecewise Robust Barrier T ub es The idea of piecewise robust barrier tubes (PRBT s) is to use robust barrier tub es (RBT s) to piecewise ov erapproximate the flowpipe segmen ts of nonlinear h ybrid systems with uncertain ty , where eac h RBT is essentially a cluster of robust barrier certificates which are situated around the flowpipe to form a tight tub e enclosing the flowpipe. The basic idea of PRBT computation is shown in Algorithm 1. Title Suppressed Due to Excessive Length 7 Algorithm 1: PRBT computation input : f : d ynamics of the system; X 0 : Initial set; U : set of uncertaint y; N : n umber of robust barrier tub es ( RBT ) in PRBT ; ( θ min , d min ) : parameters for simulation output: PRBT : piecewise robust barrier tub e 1 PRBT ← empty queue; 2 while Length ( PRBT ) < N do 3 [ F ound , θ, d ] ← [ f alse, θ 0 , d 0 ] ; 4 while θ > θ min do 5 E ← construct a coarse enclsoure-b o x for X 0 b y ( θ , d ) -simulation; 6 [ F ound , RBT , X 0 0 ] ← compute RBT inside E and obtain a set X 0 0 ⊇ ( RBT ∩ ∂ E ) ; 7 if not F ound then 8 ( θ , d ) ← 1 / 2 ∗ ( θ , d ) ; // to shrink E 9 con tinue; 10 else 11 PRBT ← Push ( PRBT , RBT ) ; // add RBT to the queue of PRBT 12 X 0 ← X 0 0 ; // update X 0 for computing next RBT 13 break; 14 if not F ound then break ; 15 return PRBT ; 4.1 Construction of the Enclosure-b o x A key step in PRBT computation is the construction of enclosure-box for a giv en compact initial set. Note that here an enclosur e-b ox is a hyp err e ctangle that entir ely c ontains a flowpip e se gment . In principle, the smaller the enclosure- b o x, the easier it is to compute a barrier tub e. Ho wev er, to make full use of the p o w er of nonlinear o verappro ximation, it is desirable to ha ve as big enclosure-b ox as p ossible so that fewer barrier tub es are needed to cov er a flowpipe. In [22], in terv al metho d was adopted to build an enclosure-b o x. How ever, the main problem with interv al metho d is that the enclosure-box thus compute d is usually very small which will result in a big num b er of barrier tub es for a fixed length of flowpipe. On the one hand, this will lead to an increasing burden on barrier tube computation. On the other hand, the capability of barrier tube in o verappro ximating complex flo wpip e can not b e fully released. F or these reasons, w e choose to use a purely simulation-based approach without losing soundness. A k ey concept inv olved in our sim ulation-based enclosure-b ox construction is twisting of tr aje ctory , whic h is a measure of maximal b ending of tra jectories in a b o x. F or the con venience of presen tation, w e present the formal definition of t wisting of tra jectory as follo ws. Definition 5 (T wisting of tra jectory). L et M b e a c ontinuous system and ζ ( t ) b e a tr aje ctory of M . Then, ζ ( t ) is said to have a twisting of θ on the time interval I = [ T 1 , T 2 ] , written as ξ I ( ζ ) , if it satisfies that ξ I ( ζ ) = θ , wher e ξ I ( ζ ) def = sup t 1 ,t 2 ∈ I arccos  h ˙ ζ ( t 1 ) , ˙ ζ ( t 2 ) i k ˙ ζ ( t 1 ) kk ˙ ζ ( t 2 ) k  . Then, w e hav e Algorithm 2 to compute enclosure-b o x. 8 Hui Kong, Ezio Barto cci, Y u Jiang, Thomas A. Henzinger Algorithm 2: Construct enclosure-b o x input : f ( x , u ) : system dynamics; X 0 : initial set; U : uncertain parameters; θ : t wisting of sim ulation; θ min : minimal theta for simulation; d : maxim um distance of simulation; output: E : an enclosure-b o x containing X 0 ; P : plane where flowpipe exits ; G : range of intersection of F low f ( X 0 ) with plane P by simulation 1 u c ← center p oin t of U ; 2 f c ( x ) ← the center dynamic f ( x , u c ) ; 3 S 0 ← sample a set of p oin ts from X 0 ; 4 select a p oin t x 0 ∈ S 0 ; 5 succ ← f alse ; 6 while θ ≥ θ min do 7 x e ← end p oin t of ( θ , d ) -simulation of f c ( x ) for x 0 ; 8 foreac h x i e : plane in the i ’th dimension of x e do 9 do simulation for all the p oin ts in S 0 , up date G and E ; 10 if al l the simulations hit x i e then 11 P ← x i e ; 12 succ ← tr ue ; 13 if succ then 14 bloat E s.t F low f ( X 0 ) exits from E only through the facet in P ; 15 return [ E , P , G ] ; 16 else 17 [ θ , d ] ← 1 / 2 ∗ [ θ , d ] ; R emark 2. In this pap er, w e assume that b oth X 0 and U are defined by hy- p errectangles. The basic idea of enclosure-b o x construction is that, given a con- tin uous dynamical system with uncertaint y , w e first remov e the uncertain ty by taking the center p oin t u c of U for the dynamics (line 1–2). Then, w e sample a set S 0 of p oin ts from X 0 for sim ulation (line 3). Prior to doing simulation for S 0 , w e first select a p oin t x 0 (usually the center p oin t of X 0 ) to do ( θ , d ) -simulation to obtain the end p oin t x e of the simulation (line 7). A ( θ , d ) -simulation is a sim- ulation that stops either when the twisting of the simulation r e aches θ or when the Euclide an distanc e b etwe en x 0 and x e r e aches d . The motiv ation to get the end p oin t x e is that, there are n planes of the form x i = x i e (the i’th element of x e ) intersecting at x e , so we wan t to c heck if one of the n planes, say P , was hit by all the simulations that start from S 0 , and if yes, it is very likely that P cut through the entire flowpipe. Then, w e take P as one of the facets of the desired enclosure-b o x E . In addition, during the simulations, we simultaneously k eep updating 1) the b oundary where th e simulations can reach and use that range as our candidate enclosure-b o x E , and 2) the b oundary range G where the sim ulations in tersect with the plane P : x i = x i e . If we end up finding such a plane P , we will push the other facets of E outw ards to make the flo wpip e exit only from this sp ecific facet of E . Of course, this ob jective cannot b e guaranteed only by simulation and pushing, we need to further chec k if the flowpipe do es not in tersect the other facets of E , which can b e done according to Theorem 3. Theorem 3. Given an unc ertain semialgebr aic system M = h X , f , X 0 , I , U i , assume E ⊂ I is an enclosur e-b ox of X 0 and F i is a fac et of E . The flowpip e of Title Suppressed Due to Excessive Length 9 M fr om X 0 do es not interse ct F i , i.e, ( F l ow f ( X 0 ) ∩ F i ) ∩ E = ∅ if ther e exists a b arrier c ertific ate B i ( x ) for F i inside E . R emark 3. Theorem 3 can b e easily prov ed by the definition of barrier certificate, whic h is ignored here. In order to make sure that the flowpipe ev ades a facet F i of E , according to Theorem 1, we only need to find a barrier certificate for F i . In the case of no barrier certificate b eing found, further bloating to the facet of E will b e p erformed. If bloating facet still end up with failure, we keep shrinking E b y setting ( θ , d ) to ( θ / 2 , d/ 2) un til barrier certificates are found for all the facet of E or θ gets less than some threshold θ min . 4.2 Computation of Robust Barrier T ube An ideal application scenario of barrier certificate is when w e can prov e the safet y prop erty using a single barrier certificate. Unfortunately , this is usually not true b ecause the flo wpip e can b e very complicated so that no p olynomial function of a sp ecified degree satisfies the constraint. In the previous subsection, w e introduce how to obtain for an initial set X 0 an enclosure-b o x E in which the system dynamics is simple enough so that a robust barrier certificate B ( x ) can b e easily computed. Therefore, we can compute a set of robust barrier certificates, whic h we call Robust Barrier T ub e (RBT), to create a tight ov erappro ximation for the flowpipe provided that there is a set of auxiliary sets serving as unsafe sets. F ormally , w e define RBT as follows. Definition 6 (Robust Barrier T ub e (RBT)). Given a semialgebr aic sys- tem M = h X , f , X 0 , I , U i , let E b e an enclosur e-b ox of X 0 and X AS = { X i AS : X i AS ⊆ E } b e a set of auxiliary sets (AS), an RBT is a set of r e al-value d func- tions Φ = { B i ( x ) , i = 1 , · · · , m } such that for al l B i ( x ) ∈ Φ : i) ∀ x ∈ X 0 : B i ( x ) > 0 , ii) ∀ ( x , u ) ∈ E × U : L f B i > 0 , and iii) ∀ x ∈ X i AS : B i ( x ) < 0 . The precision of RBT dep ends closely on the set X AS of ASs. Therefore, to deriv e a goo d barrier tub e, we need to first construct a set of high quality ASs. The factors that could affect the quality of the set X AS of ASs include 1) the n umber of ASs, and 2) the p osition, size and shap e of AS. Roughly sp eaking, the more ASs we hav e, if positioned prop erly , the more precise the RBT would b e. Regarding the p osition, size and shap e of AS, a desirable AS should 1) b e as close to the flo wpip e as p ossible, 2) spread widely around the flowpipe, and 3) be shap ed lik e a shell for the flo wpip e. In tuitively , a high quality set of ASs could b e shap ed lik e a ring around a h uman finger so that the barrier tube is tightly confined in the narrow space b et w een the ring and the finger. With the key factors aforemen tioned in mind, we developed Algorithm 3 for RBT computation. R emark 4. In principle, the more barrier certificates w e use, the b etter o verap- pro ximation w e ma y ac hieve. Ho wev er, using more barrier certificates also means more computation time. Therefore, we hav e to make a trade-off b et ween preci- sion and efficiency . In Algorithm 3, we choose to use RBT consisting of 2( n − 1) barrier certificates for n dimensional dynamical systems, which means w e need to construct 2( n − 1) ASs. W e use the same scheme as in [22] to construct ASs. Recall that we get a coarse region G where the flowpipe intersects with one of the facets of E during the construction of the enclosure-b o x E . Since G is an 10 Hui Kong, Ezio Barto cci, Y u Jiang, Thomas A. Henzinger Algorithm 3: Compute robust barrier tub e input : f : system dynamics; X 0 : Initial set; E : enclosure-b o x of X 0 ; U : set of uncertain ty; P : plane where flowpipe exits from enclosure-b o x E ; G : b o x approx. of ( P ∩ F low f ( X 0 ) ) by simulation;  : difference b et ween AS ’s (auxiliary set) output: RBT : barrier tub e; X 0 0 : b o x ov er-approx. of ( RBT ∩ E ) 1 foreac h G ij : a facet of G do 2 AS ← − CreateAS ( G , P , G ij ) ; 3 while true do 4 [ f ound, B ij ] ← − ComputeRBC ( f , X 0 , E , AS , U ) ; 5 if f ound then AS ’ ← − Expand ( AS ) ; 6 else AS ’ ← − Contract ( AS ) ; 7 if Diff ( AS 0 , AS ) ≤  then 8 break; 9 AS ← − AS ’; 10 if f ound then 11 RBT ← − Push ( RBT , B ij ) ; 12 break; 13 else 14 return F AIL 15 return SUCCEED; n − 1 dimensional b ox, the RBT must con tain G . Therefore, we choose to con- struct 2( n − 1) ASs whic h are able to form a tigh t hollow hyper-rectangle around G . The idea is that for each facet G ij of G , we construct an n − 1 dimensional h yp er-rectangle b et ween G ij and E ij as an AS (line 2), where E ij is the n − 1 di- mensional face of E that corresp onds to G . Then, we use Algorithm 3 to compute an RBT(line 4). In the while lo op 3, we try to find the best barrier certificate b y adjusting the width of AS (line 5 and 6) iteratively until the difference in width b et ween tw o consecutive AS s is less than the sp ecified threshold  . T o b e in tuitive, we provide Figure 1 to demonstrate the pro cess. 4.3 PRBT for Con tinuous Dynamics The idea of computing PRBT is straightforw ard. Given an initial set X 0 , we first construct a coarse enclosure-b ox E containing X 0 and then we further compute an RBT inside E to get a muc h more precise ov erapproximation for the flowpipe. Mean while, we obtain a hyper-rectangle R formed b y ASs with a hollow X 0 0 in the middle. Since the in tersection of the RBT and the facet of E is contained en tirely in the hollow X 0 0 of R , w e use X 0 0 as a new initial set and rep eat the en tire pro cess to compute a PRBT step b y step. Since our approac h is time indep enden t, the length of a PRBT cannot b e measured by the length of time horizon. Hence, in our implementation, we try to compute a sp ecified num b er of RBT s. Title Suppressed Due to Excessive Length 11 ! : enclosure (a) More simulation ! : enclosure (b) bloat ! : enclosure (c) : enclosure (d) : enclosure bloated ! " (e) : enclosure bloated ! " " " #$%&'( )*+,'- . (f ) ! " # " : enclosure bloated (g) Fig. 1: (a) → (g): demonstration of RBT computation 4.4 PRBT for Hybrid Dynamics T o extend our approach with the ability to deal with hybrid systems, w e need to handle tw o problems i) compute the in tersection of RBT and guard set, and ii) compute the image of the intersection after discrete jump. In general, these t wo issues can b e very hard dep ending on what kind of guard sets and transitions are defined for the h ybrid systems. In this pap er, we make some assumptions on the hybrid systems under con- sideration. Let a discrete transition τ b e defined as follo ws. τ ll 0 = h Guard ll 0 , Trans ll 0 i (8) where l and l 0 are the lo cations of the dynamics b efore and after a discrete tran- sition resp ectiv ely , Guard ll 0 = { x ∈ R n | x i ∼ b i , ∼∈ {≤ , ≥}} and Trans ll 0 : x 0 = A x , where A is an n -dimensional matrix. Based on this assump tion, the problem of computing the in tersection of RBT and guard set is reduced to computing the in tersection of RBT with a plane of x i = b i , which can be handled using a similar strategy to computing the in tersection of RBT with the facet of enclosure-b o x. Hence, we hav e Algorithm 4 to deal with discrete transition of a hybrid system. R emark 5. The strategy to deal with discrete transitions of h ybrid systems is that every time we obtain an enclosure-b o x E , we first detect whether E intersects with some guard set Guard ll 0 . If no, we pro ceed with the normal pro cess of PRBT computation. Otherwise, we switch to the pro cedure of Algorithm 4 in which the input X l 0 is the last state set whose enclosure-b o x intersects with Guard ll 0 . Since the flowpipe ma y not cross through the guard plane en tirely , we use the while lo op in line 2 to compute an ov erapproximation for the in tersection. The basic idea of the while lo op is that, given a state set X l 0 , we first construct an enclosure- b o x E b y simulation(line 4), if E intersects with Guard ll 0 , we shrink E by cutting off the part of E that lies in the guard set (line 7). As a result of this op eration, the flo wpip e could exit from E not only through the guard plane but also through other facets of E . F or eac h of those facets, we compute an o verappro ximation for its intersection with the flowpipe using simulation and barrier certificate 12 Hui Kong, Ezio Barto cci, Y u Jiang, Thomas A. Henzinger Algorithm 4: handle discrete transition of hybrid system input : X l 0 : intermediate initial set at lo cation l ; Guard ll 0 : guard set of transition τ ll 0 ; T rans ll 0 : image mapping of transition τ ll 0 output: X l 0 0 : image of transition τ ll 0 1 InitQ ← Push ( InitQ , X l 0 ) ; 2 while InitQ not empty do 3 X l 0 ← Pop ( InitQ ) ; 4 E ← construct enclosure-b o x for X l 0 ; 5 if E ∩ Guard ll 0 == ∅ then 6 con tinue; 7 E ← E ∩ Guard ll 0 ; 8 X Φ ∩ E ← do simulation and barrier certificate computation to find an o verappro ximation for the region where the flowpipe Φ intersects with the guard plane x i = b i ; 9 Q Φ ∩ E ← Push ( Q Φ ∩ E , X Φ ∩ E ) ; 10 foreac h E ij : facet of E except guar d plane do 11 X ij 0 ← do simulation and barrier certificate computation to an o verappro ximation for the region where the barrier tub e intersects with E ij ; InitQ ← Push ( InitQ , X ij 0 ) ; 12 X Φ ∩ E ← b o x ov erapprox. Q Φ ∩ E ; 13 X l 0 0 ← T rans ll 0 X Φ ∩ E ; computation(line 8 and 11). In addition, since thos e intersections X ij 0 that do not lie in the guard plane could still reach the guard plane later, we therefore push them in to a queue for further exploration. 5 Implemen tation and Exp eriments W e ha ve dev elop ed PRBT, a soft ware prototype written in C++ that imple- men ts the concepts and the algorithms presented in this paper. PRBT computes piecewise robust barrier tub es for nonlinear con tinuous and hybrid systems with p olynomial dynamics. W e compare our approach in efficiency and precision with the state-of-the-art to ols Flo w* and CORA using several b enchmarks of non- linear contin uous and hybrid systems. Note that since C2E2 do es not supp ort uncertain ty , so w e cannot compare with it. The exp erimen ts were carried out on a desktop computer with a 3 . 6 GHz Intel 8 Cor e i7-7700 CPU and 32 GB memory . 5.1 Nonlinear Contin uous Systems W e consider six nonlinear benchmark systems with p olynomial dynamics for whic h their mo dels and settings are provided in T able 1. The exp erimental results are rep orted in T able 2. Since our approach is time indep enden t, which is different from Flow* and CORA, to make the comparison fair enough, we choose to compute a slightly longer flowpipe than the other t wo to ols. Note that there are t wo columns for time for Flo w*. The reason wh y we Title Suppressed Due to Excessive Length 13 T able 1: Contin uous dynamical mo del definitions Mo del Dynamics Uncertain ty X 0 Con troller 2D ˙ x = d 1 xy + y 3 + 2 d 1 ∈ [0 . 95 , 1 . 05] x ∈ [29 . 9 , 30 . 1] ˙ y = d 2 x 2 + 2 x − 3 y d 2 ∈ [0 . 95 , 1 . 05] y ∈ [ − 38 , − 36] V an der Pol ˙ x = y + d 1 d 1 ∈ [ − 0 . 01 , 0 . 01] x ∈ [1 , 1 . 5] Oscillator ˙ y = y − x − x 2 y + d 2 d 2 ∈ [ − 0 . 01 , 0 . 01] y ∈ [2 . 40 , 2 . 45] Lotk a-V olterra ˙ x = x (1 . 5 − y ) + d 1 d 1 ∈ [ − 0 . 01 , 0 . 01] x ∈ [4 . 6 , 5 . 5] ˙ y = − y (3 − x ) − d 2 d 2 ∈ [ − 0 . 01 , 0 . 01] y ∈ [1 . 6 , 1 . 7] Buc kling ˙ x = y + d 1 d 1 ∈ [ − 0 . 01 , 0 . 01] x ∈ [ − 0 . 5 , − 0 . 4] Column ˙ y = 2 x − x 3 − 0 . 2 y + 0 . 1 + d 2 d 1 ∈ [ − 0 . 01 , 0 . 01] y ∈ [ − 0 . 5 , − 0 . 4] Jet Engine ˙ x = − y − 1 . 5 x 2 − 0 . 5 x 3 − 0 . 5 d 1 ∈ [ − 0 . 005 , 0 . 005] x ∈ [1 . 19 , 1 . 21] + d 1 d 2 ∈ [ − 0 . 005 , 0 . 005] y ∈ [0 . 8 , 1 . 0] ˙ y = 3 x − y + d 2 ˙ x = 10( y − x ) + d 1 d 1 ∈ [ − 0 . 001 , 0 . 001] x ∈ [1 . 79 , 1 . 81] Con troller 3D ˙ y = x 3 + d 2 d 2 ∈ [ − 0 . 001 , 0 . 001] y ∈ [1 . 0 , 1 . 1] ˙ z = xy − 2 . 667 z z ∈ [0 . 5 , 0 . 6] ha ve an extra time column for Flo w* is that it can b e v ery fast and precise to compute the T aylor mo del for a given system. How ev er, T a ylor mo dels cannot b e in general applied directly to solve the safety verification problem. Checking their intersection with the unsafe set requires their transformation in to simpler geometric form (e.g. b o x), which has an exp onen tial complexit y in the n umber of the dimensions and it needs to b e considered in the ov erall running time. T able 2: Exp erimen tal results for b enc hmark systems. #v ar: num b er of v ariables; T: computing time for flo wpip e; TT: computing time in- cluding b o x transformation; N: n umber of flowpipe segmen ts; D: can- didate degrees for template p olynomial (for PRBT only); TH: time horizon for flowpipe (for Flo w* and CORA only). F/E: failed to ter- minate under 30min or exception happ ened. PRBT Flo w* CORA Mo del #v ar T N D TH T TT N T N Con troller 2D 2 35.73 12 3 0.012 48.8 417.64 240 F 1200 V an der Pol 2 221.62 17 4 6.74 23.88 1111.05 135 E – Lotk a-V olterra 2 30.10 9 4 3.20 6.06 405.32 320 40.30 160 Buc kling Column 2 74.02 35 3 14.00 F – – 734.81 1400 Jet Engine 2 240.98 18 4 9.50 F – – 1.69 190 Con troller 3D 3 774.58 20 3 0.55 F – – F – R emark 6. T able 2 shows us how brutal the reality of reac hability analysis of nonlinear systems is and this gets ev en worse in the presence of uncertain t y and large initial set. As can b e seen in T able 2, b oth Flo w* and CORA failed to give a solution for half of the benchmarks either due to timeout or due to exception. This phenomenon may b e alleviated if smaller initial sets are provided or uncertaint y is remov ed. In terms of computing time T , PRBT do es not alwa ys 14 Hui Kong, Ezio Barto cci, Y u Jiang, Thomas A. Henzinger (a) PRBT (b) Flow* (c) CORA (d) PRBT Fig. 2: Lotk a-V olterra: (a), (b), (c); controller 3D: (d) (a) PRBT (b) Flo w* (c) PRBT (d) Flow* Fig. 3: Controller 2d: (a), (b); V an der P ol Oscillator: (c), (d) (a) PRBT (b) CORA (c) PRBT (d) CORA Fig. 4: Buckling Column: (a), (b); Jet Engine: (c), (d) outp erform the other tw o to ols. A ctually , Flow* or CORA can b e muc h faster in some cases. Ho wev er, when the b ox transformation time for T aylor model w as tak en into account, the total computing time T T of Flow* increased significantly . One p oin t to note here is that, PRBT, in general, pro duces a muc h smaller n umber N of flowpipe segmen ts than the other tw o, whic h means that the time used to chec k the intersection of flo wpip e with the unsafe set can b e reduced considerably . In addition, as shown in Figure 2–4, PRBT is more precise than the other t wo on av erage. 5.2 Nonlinear Hybrid System W e use the tunnel dio de oscillator (TDO) circuit (with different setting) intro- duced in [20] to illustrate the application of our approach to h ybrid system. The tw o state space v ariables are the v oltage x 1 = V C across capacitor and the curren t x 2 = I L through the inductor. The system dynamics is describ ed as follo ws, ˙ x 1 = 1 C ( − h ( x 1 ) + x 2 ) ˙ x 2 = 1 L ( − x 1 − x 2 G + V in ) where h ( x 1 ) describes the tunnel dio de c haracteristic and V in = 0 . 3 V , G = 5 m ohm − 1 , L = 0 . 5 µH and C = 2 pF . Title Suppressed Due to Excessive Length 15 . ) ( 2 x f x = . 3 . 0 1 ³ x 3 . 0 1 £ x ) ( 1 x f x = 1 l 2 l 5 . 0 3 . 0 1 £ £ x 3 . 0 0691 . 0 1 £ £ x ) ( 3 x f x = 3 l 0691 . 0 0 1 £ £ x 0691 . 0 1 £ x 0691 . 0 1 ³ x . Fig. 5: Hybridized mo del of TDO (a) PRBT −0.2 0 0.2 0.4 0.6 0.8 1 1.2 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 x2 x1 (b) Flow* (c) CORA Fig. 6: Flowpipe of hybridized TDO F or this mo del, we wan t to define an initial state region X 0 whic h can guar- an tee the oscillating behaviour for the system. Due to the highly nonlinear b e- ha viour of the system, a common strategy to deal with this mo del is to use a h ybridized mo del to approximate the dynamics system and then apply formal v erification to the hybrid mo del[13,18]. In our exp erimen t, we use three cubic equations to appro ximate the curve of h ( x 1 ) . h ( x 1 ) =    0 . 000847012 + 35 . 2297 x 1 − 395 . 261 x 2 1 + 1372 . 29 x 3 1 , 0 ≤ x 1 ≤ 0 . 0691 1 . 242 − 0 . 033 x 1 − 47 . 4311 x 2 1 + 116 . 48 x 3 1 , 0 . 0691 ≤ x 1 ≤ 0 . 3 − 16 . 544 + 139 . 64 x 1 − 389 . 245 x 2 1 + 359 . 948 x 3 1 , 0 . 3 ≤ x 1 ≤ 0 . 50 F rom the piecewise function h ( x 1 ) , we can derive a 3 -mo de hybrid system whic h is sho wn in Figure 5. The system switches b et w een the lo cations as the v alue of x 1 c hanges. Let the initial set b e X 0 = { ( x 1 , x 2 ) ∈ R 2 | 0 . 40 ≤ x 1 ≤ 0 . 48 , 0 . 38 ≤ x 2 ≤ 0 . 39 } on lo cation l 3 , w e compute an ov erapproximation for the flo wpip e using PRBT, Flow* and CORA resp ectiv ely . As illustrated in Figure 6, b oth PRBT and CORA found an inv ariant with roughly the same precision, which indicates the mo del oscillates for the initial set, while Flow* ran into an error. 6 Conclusion W e prop ose a nov el metho d to compute efficiently a tight o verappro ximation (flo wpip e) of the reachable set for nonlinear contin uous and hybrid systems with uncertain ty . Our approach av oid the use of interv al metho d for enclosure-b o x computation, which can pro duce bigger enclosure-b o xes and fewer flo wpip e seg- men ts. Experiments on sev eral benchmark systems sho w that our metho d is more efficien t and precise than other p opular approaches prop osed in literature. 16 Hui Kong, Ezio Barto cci, Y u Jiang, Thomas A. Henzinger References 1. Matthias Althoff and Dmitry Grebenyuk. Implementation of interv al arithmetic in CORA 2016. In Pr o c. of ARCH , volume 43 of EPiC Series in Computing , pages 91–105. EasyChair, 2017. 2. Eugene Asarin, Thao Dang, and Antoine Girard. Hybridization methods for the analysis of nonlinear systems. A cta Informatic a , 43(7):451–476, 2007. 3. Mohamed Amin Ben Sassi, Sriram Sank aranaray anan, Xin Chen, and Erik a Ábrahám. Linear relaxations of p olynomial p ositivit y for p olynomial ly apunov function synthesis. IMA Journal of Mathematic al Contr ol and Information , 33(3):723–756, 2015. 4. Martin Berz and Ky oko Makino. V erified integration of odes and flo ws using differen tial algebraic metho ds on high-order ta ylor mo dels. R eliable Computing , 4(4):361–369, 1998. 5. Sergiy Bogomolov, Christian Schilling, Ezio Barto cci, Grégory Batt, Hui Kong, and Radu Grosu. Abstraction-based parameter synthesis for multiaffine systems. In Pr o c. of HVC: the 11th International Haifa V erific ation Confer enc e on Har dwar e and Software: V erific ation and T esting , volume 9434 of LNCS , pages 19–35, 2015. 6. Xin Chen, Erik a Ábrahám, and Sriram Sank aranaray anan. Flo w*: An analyzer for non-linear hybrid systems. In Pr o c. of CA V 2013: the 25th International Con- fer enc e on Computer Aide d V erific ation , v olume 8044 of LCNS , pages 258–263. Springer, 2013. 7. Alessandro Cimatti, Alb erto Griggio, Ahmed Irfan, Marco Ro veri, and Rob erto Sebastiani. Exp erimen ting on solving nonlinear integer arithmetic with incremen tal linearization. In Pro c. of SA T 2018: the 21st International Conferenc e on The ory and Applic ations of Satisfiability T esting , volume 10929 of LNCS , pages 383–398. Springer, 2018. 8. Alessandro Cimatti, Alb erto Griggio, Ahmed Irfan, Marco Ro veri, and Rob erto Sebastiani. Incremen tal linearization for satisfiability and verification mo dulo nonlinear arithmetic and transcendental functions. ACM T rans. Comput. L o g. , 19(3):19:1–19:52, 2018. 9. Jacek Cyrank a, Md. Ariful Islam, Greg Byrne, Paul L. Jones, Scott A. Smolk a, and Radu Grosu. Lagrangian reachabililt y . In Pro c. of CA V 2017: the 29th In- ternational Confer enc e on Computer Aide d V erific ation , volume 10426 of LNCS , pages 379–400. Springer, 2017. 10. Jacek Cyrank a, Md. Ariful Islam, Scott A. Smolk a, Sicun Gao, and Radu Grosu. Tigh t con tinuous-time reach tubes for lagrangian reac hability . In Pr o c. of CDC 2018: 57th IEEE Confer enc e on De cision and Contr ol , page to app ear. IEEE, 2018. 11. P arasara Sridhar Duggirala, Sa yan Mitra, Mahesh Viswanathan, and Matthew P otok. C2E2: A verification to ol for stateflow mo dels. In Pr o c. of T ACAS 2015: the 21st International Conferenc e on T ools and Algorithms for the Construction and Analysis of Systems , v olume 9035 of LNCS , pages 68–82. Springer, 2015. 12. Martin F ränzle, Christian Herde, Tino T eige, Stefan Ratschan, and T obias Sc hu- b ert. Efficient solving of large non-linear arithmetic constraint systems with com- plex b oolean structure. JSA T , 1(3-4):209–236, 2007. 13. G F rehse, BH Krogh, and RA Ruten bar. V erification of hybrid systems using iterativ e refinemen t. Pr o c. SR C TECHCON 2005, Portland, USA, Oct. 24-26 , 2005. 14. Goran F rehse, Colas Le Guernic, Alexandre Donzé, Scott Cotton, Ra jarshi Ray , Olivier Lebeltel, Ro dolfo Ripado, Antoine Girard, Thao Dang, and Oded Maler. SpaceEx: Scalable v erification of hybrid systems. In Pr o c. of CA V 2011: the 23r d International Confer enc e on Computer A ide d V erific ation , volume 6806 of LNCS , pages 379–395. Springer, Springer, 2011. Title Suppressed Due to Excessive Length 17 15. An toine Girard and Colas Le Guernic. Efficient reac hability analysis for linear systems using supp ort functions. Pr o c. of IF AC W orld Congr ess , 41(2):8966–8971, 2008. 16. Radu Grosu, Grégory Batt, Flavio H. F enton, James Glimm, Colas Le Gu ernic, Scott A. Smolk a, and Ezio Barto cci. F rom cardiac cells to genetic regulatory net works. In Pr o c. of CA V 2011: the 23r d International Confer enc e Computer Aide d V erific ation , volume 6806 of LNCS , pages 396–411. Springer, 2011. 17. S. Gulwani and A. Tiwari. Constraint-based approach for analysis of hybrid sys- tems. In Pro c. of CA V 2008: the 20th International Confer enc e on Computer Aide d V erific ation , volume 5123 of LNCS , p ages 190–203. Springer, 2008. 18. Smriti Gupta, Bruce H Krogh, and Rob A Rutenbar. T ow ards formal verification of analog and mixed-signal designs. TECHCON, 2003. 19. Amit Gurung, Ra jarshi Ra y , Ezio Bartocci, Sergiy Bogomolov, and Radu Grosu. P arallel reachabilit y analysis of h ybrid systems in xspeed. International Journal on Software T o ols for T e chnolo gy T r ansfer , to app ear:1–23, 2018. 20. W alter Hartong, Lars Hedrich, and Erich Barke. Mo del chec king algorithms for analog verification. In Pro c e e dings of the 39th annual Design Automation Confer- enc e , pages 542–547. ACM, 2002. 21. T.A. Henzinger. The theory of h ybrid automata. In Pro c. IEEE Symp. L o gic in Computer Science , pages 278–292, 1996. 22. Hui K ong, Ezio Barto cci, and Thomas A. Henzinger. Reac hable set ov er- appro ximation for nonlinear systems using piecewise barrier tubes. In Pr o c. of CA V 2018: the 30th International Confer enc e on Computer A ide d V erific ation , v olume 10981 of LNCS , pages 449–467. Springer, 2018. 23. Hui Kong, Sergiy Bogomolov, Christian Schilling, Y u Jiang, and Thomas A. Hen- zinger. Safety verification of nonlinear hybrid systems based on inv ariant clusters. In Pr o c. of HSCC 2017: the 20th International Confer enc e on Hybrid Systems: Computation and Control , pages 163–172. ACM, 2017. 24. Hui Kong, F ei He, Xiaoyu Song, William NN Hung, and Ming Gu. Exp onen tial- condition-based barrier certificate generation for safety verification of hybrid sys- tems. In Pro c. of CA V 2013: the 25th International Confer enc e on Computer Aide d V erific ation , volume 8044 of LNCS , p ages 242–257. Springer, 2013. 25. So onho Kong, Sicun Gao, W ei Chen, and Edmund M. Clarke. dReac h: δ - reac hability analysis for hybrid systems. In Pro c. of T ACAS’15: the 21st Interna- tional Confer enc e on T o ols and Algorithms for the Construction and Analysis of Systems , volume 9035 of LNCS , pages 200–205. Springer, 2015. 26. T omas Krilavicius. Hybrid T e chniques for Hybrid Systems . PhD thesis, Universit y of T wen te, Enschede, Netherlands, 2006. 27. Jean B Lasserre. P olynomial programming: Lp-relaxations also conv erge. SIAM Journal on Optimization , 15(2):383–393, 2005. 28. J. Liu, N. Zhan, and H. Zhao. Computing semi-algebraic inv ariants for p olynomial dynamical systems. In Pr o c. of EMSOFT 2011: the 11th International Confer enc e on Embe dde d Softwar e , pages 97–106. ACM, 2011. 29. Nadir Matringe, Arnaldo Vieira Moura, and Rachid Rebiha. Generating in v ariants for non-linear hybrid systems by linear algebraic metho ds. In Pr o c. of SAS 2010: the 17th International Symp osium on Static Analysis , volume 6337 of LNCS , pages 373–389. Springer, 2010. 30. Nedialk o S Nedialko v. Interv al tools for ODEs and DAEs. In Pr o c. of SCAN 2006: the 12th GAMM - IMACS International Symp osium on Scientific Comput- ing, Computer Arithmetic and V alidate d Numerics , pages 4–4. IEEE, 2006. 31. P avithra Prabhak ar and Miriam García Soto. Hybridization for stabilit y analysis of switched linear systems. In Pr o c. of HSCC 2016: of the 19th Intern. Confer enc e on Hybrid Systems: Computation and Contr ol , pages 71–80. ACM, 2016. 32. S. Pra jna and A. Jadbabaie. Safety verification of hybrid systems using barrier cer- tificates. Pr o c. of HSCC 2004: the 7th International W orkshop on Hybrid Systems: Computation and Control , 2993:271–274, 2004. 18 Hui Kong, Ezio Barto cci, Y u Jiang, Thomas A. Henzinger 33. Mihai Putinar. P ositive polynomials on compact semi-algebraic sets. Indiana University Mathematics Journal , 42(3):969–984, 1993. 34. Ra jarshi Ray , Amit Gurung, Binay ak Das, Ezio Barto cci, Sergiy Bogomolov, and Radu Grosu. Xsp eed: A ccelerating reachabilit y analysis on m ulti-core pro cessors. In Pr o c. of HV C 2015: the 11th International Haifa V erific ation Confer enc e on Har dwar e and Softwar e: V erific ation and T esting , v olume 9434 of LNCS , pages 3–18. Springer, 2015. 35. Nima Ro ohi, Pa vithra Prabhak ar, and Mahesh Viswanathan. Hybridization based CEGAR for h ybrid automata with affine dynamics. In Pr o c. of T ACAS 2016: the 22nd International Confer enc e on T o ols and Algorithms for the Construction and Analysis of Systems , volume 9636 of LNCS , pages 752–769. Springer, 2016. 36. S. Sank aranaray anan. Automatic in v ariant generation for hybrid systems using ideal fixed p oin ts. In Pr o c. of HSCC 2010: the 13th ACM International Confer enc e on Hybrid Systems: Computation and Contr ol , pages 221–230. ACM, 2010. 37. S. Sank aranara yanan, H. Sipma, and Z. Manna. Constructing inv ariants for hybrid systems. In Pr o c. of HSCC 2004: the 7th Intern. W orkshop on Hybrid Systems: Computation and Control , volume 2993 of LNCS , pages 539–554. Springer, 2004. 38. Sriram Sank aranaray anan, Xin Chen, et al. Lyapuno v function synthesis using handelman representations. IF AC Pr o c e e dings V olumes , 46(23):576–581, 2013. 39. Stefan Sch upp, Erik a Ábrahám, Ibtissem Ben Makhlouf, and Stefan K ow alewski. Hypro: A C++ library of state set representations for hybrid systems reachabilit y analysis. In Pr o c. of NFM 2017: the 9th International Symp osium on NASA F ormal Metho ds , volume 10227 of LNCS , pages 288–294, 2017. 40. Andrew Sogok on, Khalil Ghorbal, Paul B Jackson, and André Platzer. A metho d for in v ariant generation for p olynomial contin uous systems. In Pr o c. of VMCAI 2016: the 17th International Confer enc e on V erific ation, Model Che cking, and Ab- str act Interpr etation , v olume 9583 of LNCS , pages 268–288. Springer, 2016. 41. Gilb ert Stengle. A Nullstellensatz and a Positivstellensatz in semialgebraic geom- etry . Mathematische Annalen , 207(2):87–97, 1974. 42. A. T aly and A. Tiwari. Deductive verification of contin uous dynamical systems. In FSTTCS , volume 4, pages 383–394, 2009. 43. Zhengfeng Y ang, Chao Huang, Xin Chen, W ang Lin, and Zhiming Liu. A linear programming relaxation based approach for generating barrier certificates of hybrid systems. In International Symp osium on F ormal Metho ds , pages 721–738. Springer, 2016.