Invariant, Viability and Discriminating Kernel Under-Approximation via Zonotope Scaling

Invariant, Viability and Discriminating K ernel Under- Appro ximation via Zonotope Scaling Ian M. Mitchell Department of Computer Science The University of British Columbia V ancouver , British Columbia , Canada mitchell@cs.ubc.ca Jacob Budzis Engineering Physics Program The University of British Columbia V ancouver , British Columbia , Canada jrbudzis@gmail.com Andriy Bolyachevets Department of Computer Science The University of British Columbia V ancouver , British Columbia , Canada andriy .bolyachevets@alumni.ubc.ca ABSTRA CT Scalable safety verication of continuous state dynamic systems has been demonstrated through both r eachability and viability analyses using parametric set representations; howev er , these two analyses are not interchangable in practice for such parametric repr esenta- tions. In this paper we consider viability analysis for discrete time ane dynamic systems with adversarial inputs. Given a set of state and input constraints, and treating the inputs in best-case and/or worst-case fashion, we construct invariant, viable and discriminat- ing sets, which must therefore under-approximate the invariant, viable and discriminating kernels respectively . The sets are con- structed by scaling zonotop es represented in center-generator form. The scale factors are found through ecient convex optimizations. The results ar e demonstrated on two toy examples and a six dimen- sional nonlinear longitudinal model of a quadrotor . KEY W ORDS reachability , viability , controlled invariance, r obust safety analysis, set-theoretic methods, convex optimization, zonotopes A CM Reference format: Ian M. Mitchell, Jacob Budzis, and Andriy Bolyachevets. 2016. Invariant, Viability and Discriminating Kernel Under- Approximation via Zonotope Scaling. In Proceedings of preprint, arXiv .org, 2019, 10 pages. DOI: 10.1145/nnnnnnn.nnnnnnn 1 IN TRODUCTION Reachability analysis is a rigorous alternative to sampling-base d verication of dynamic systems, and at least for linear (or ane) sys- tems there have been recent demonstrations of techniques capable of handling thousands of continuous state space dimensions [ 6 , 7 ]. Reachable sets and tub es—or more typically over-approximations of them to ensure soundness—are an eective tool for demonstrating safety: If the for ward / backward reach set or tub e does not intersect the unsafe / initial set respectively then the system is safe. Any input or parameter uncertainty is typically treated in a worst-case fashion to make the reach set or tube larger and hence the system less likely to be judged safe. In this paper we adopt the alternativ e appr oach to proving safety advocated in viability theor y [ 5 ] but more commonly framed as Permission to make digital or hard copies of part or all of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for prot or commercial advantage and that copies bear this notice and the full citation on the rst page. Copyrights for third-party components of this work must be honor ed. For all other uses, contact the owner /author(s). preprint, arXiv .org © 2019 Copyright held by the owner/author(s). 978-x-xx xx-xxx x-x/YY/MM. . . $15.00 DOI: 10.1145/nnnnnnn.nnnnnnn various versions of invariant sets: Find the set of states fr om which trajectories of the system are guarantee d to satisfy a specie d safety constraint. A critical feature of viable or controlled invariant sets is that a control input may b e chosen in a b est-case fashion to keep the trajectories safe. Robust versions of these sets (aka discriminating sets) also allow an adv ersarial disturbance input which is tr eated in a worst-case fashion to drive the trajectories out of the safety constraints. In every case we must under-appr oximate the results to ensure soundness of the safety analysis. While viability analysis can be reduced to reachability analysis and vice versa in theor y , the parametric representations which can handle high dimensional systems do not support such reductions; for example, the reductions r equire set complements but the para- metric representations are usually r estricted to convex sets whose complements are non-convex. The ne ed to develop scalable algorithms for viability / invari- ance in addition to those for reachability is therefore practical: The former require under-approximation while the latter over- approximation, and some analyses are more naturally amenable to parametric representations in one formulation or the other , but rarely both. Although we do not have space to explore it in this paper , an example of an application which naturally ts into the via- bility framework is testing at run-time whether an exogenous input signal—such as might arise from a human-in-the-loop contr ol—will maintain safety; for example, see [17]. The focus of this paper is therefore development of mor e scalable formulations for (robust) invariance / viability based analysis of ane continuous state dynamic systems. Scalability is achieved by framing the calculations as convex optimizations to nd ecient parametric representations in the form of zonotopes. The specic contributions of this paper are to: • Show how a nite horizon invariance kernel of an ane system with disturbance input can b e underapproximated with a zonotope via a linear program. • Extend the formulation to allow control inputs and thereby underapproximate nite horizon viability and discriminat- ing kernels with zonotopes via a convex program. • Demonstrate the use of the discriminating kernel to com- pute a larger robust controlled invariant set for a moderate dimensional nonlinear quadr otor model than was achieved using an ellipsoidal representation in [17]. 2 RELA TED W ORK Reachability and viability have been applied to a broad variety of dierent dynamic systems resulting in a huge range of dierent algorithms. W e focus here on parametric approaches for linear or preprint, 2019, arXiv .org Ian M. Mitchell, Jacob Budzis, and Andriy Bolyachevets ane dynamics. Ellipsoid and support vector parametric represen- tations of viability constructs were explored in [ 15 , 16 ]. Zonotopic representations for reachability were introduced in [ 9 ] and have since been extensively explored; for e xample [ 3 , 4 , 10 ]. Our work was inspired, how ever , by the papers [ 18 , 20 ] and in particular [ 19 ], which utilizes a convex optimization to select zonotope generator weights to construct a control scheme that will drive a set of initial states into the smallest possible set of nal states. Also similar to this work is [ 14 ] in which the authors se ek a linear fee dback control input which will maximize the size of the backward reachable set and arrive at a bilinear matrix inequality based on containment of one zonotope within another . Signicantly , unlike most other work these papers treat the input in a best-case fashion. The dierence with the work b elow lies in the reachability vs viability formula- tion and the fact that our approach constructs a set-valued viable feedback control from a convex optimization. 3 PRELIMINARIES For a matrix M , let | M | denote the elementwise absolute value and ( M ) i , j denote the element in row i and column j ; therefore, ( M s ) i , j denotes the element in row i and column j of the matrix power M s . Also dene 1 d and 0 d to be the vectors in R d of all ones and all zeros respectively . 3.1 System Dynamics Consider a discrete time dynamic system for times t = 0 , 1 , 2 , . . . : x ( t + 1 ) = A x ( t ) + B u ( t ) + C v ( t ) + w (1) where • The state is x ∈ R d x . • The control input is u ∈ U ⊂ R d u . The control input constraint U is an interval hull (or hyperrectangle) dene d by the elementwise inequalities U = { u ∈ R d u | u ≤ u ≤ u } . (2) • The disturbance input is v ∈ V ⊂ R d v . The disturbance input constraint V is a zonotop e (see section 3.3). • The drift w ∈ R d x is constant. W e note that w , 0 can alternatively be treated by osetting the center of the zono- tope V , but we will carr y w through separately so as to support a drift term for disturbance-free systems. Given an initial condition x ( 0 ) and input signals u ( ·) and v (·) , the solution of (1) is given by x ( t ) = A t x ( 0 ) + t − 1 Õ s = 0 A t − 1 − s ( B u ( s ) + C v ( s ) + w ) . (3) 3.2 Sets of Interest for Safety V erication W e will formulate our safety verication problem as keeping the system state within a constraint set X . For reasons of notational simplicity , we will assume in the rest of the paper that X is an interval hull or box constraint of the form X = { x ∈ R d x | x ≤ x ≤ x } . (4) It is possible to relax this assumption; see se ction 7 for further comments. W e seek to approximate various subsets of the constraint set. The most general is a nite horizon discriminating or robust controlled invariant set Disc ( [ 0 , T ] , X ) , ( x ( 0 ) ∈ X      ∃ u (·) , ∀ v (·) , ∀ t ∈ [ 0 , T ] , x ( t ) ∈ X ) , (5) where u ( t ) ∈ U and v ( t ) ∈ V for all t = 0 , . . . , T . Note that the control input u ( t ) tries to keep the system state within the constraint set X , while the disturbance input v ( t ) is treated in a worst-case or adversarial fashion and tries to drive the system state outside the constraint set. W e will also consider two special cases of the discriminating set. For systems which lack a control input, we se ek a nite horizon invariant set Inv ( [ 0 , T ] , X ) , ( x ( 0 ) ∈ X      ∀ v (·) , ∀ t ∈ [ 0 , T ] , x ( t ) ∈ X ) , (6) while for systems which lack a disturbance input w e seek a nite horizon viable or controlled invariant set Viab ( [ 0 , T ] , X ) , ( x ( 0 ) ∈ X      ∃ u (·) , ∀ t ∈ [ 0 , T ] , x ( t ) ∈ X ) . (7) Our algorithm for computing these previous sets will often make use of the forward reach set of some specied set S . R ( t , S ) , n x ( t ) ∈ R d x | ∃ w ( ·) , x ( 0 ) ∈ S o , (8) where the choice of input signal w (·) = u ( ·) or w (·) = v ( ·) should be clear from context. Unlike the discriminating, invariant and viable sets, the reach set is dene d at a single time rather than over an interval, and the set S is an initial condition rather than a constraint. 3.3 Set Representation: Zonotopes W e will use zonotopes as our parametric representation of invariant, viable and/or discriminating sets. A zonotope S ⊆ R d is a polytope dened by a center c ( S ) ∈ R d and a nite number of generators д i ( S ) ∈ R d for i = 1 , . . . , p ( S ) : S =        c ( S ) + p ( S ) Õ i = 1 λ i д i ( S )       − 1 ≤ λ i ≤ + 1        . (9) In most cases it is more convenient to work with the center-generator tuple (or G-repr esentation for a zonotope S rather than (9) S = h c ( S ) | д 1 ( S ) , д 2 ( S ) , . . . , д p ( S ) ( S ) i , = h c ( S ) | G ( S ) i , where the generator matrix G ( S ) is formed by horizontal concate- nation of the generator vectors G ( S ) =  д 1 ( S ) д 2 ( S ) · · · д p ( S ) ( S )  ∈ R d × p ( S ) . When it is necessar y to refer to individual elements of a generator matrix or vector , we will use the notation д j , i ( S ) to specify the ele- ment in the j t h row and i t h column of matrix G ( S ) , or equivalently the j t h element of vector д i ( S ) . Invariant, Viability and Discriminating Kernel Under- Approximation via Zonotope Scaling preprint, 2019, arXiv .org With the generator matrix notation, we can write (9) more com- pactly as S = n c ( S ) + G ( S ) λ    − 1 p ( S ) ≤ λ ≤ + 1 p ( S ) o . (10) For a vector x ∈ S , dene λ ( x , S ) such that x = c ( S ) + G ( S ) λ ( x , S ) and − 1 p ( S ) ≤ λ ( x , S ) ≤ + 1 p ( S ) (11) and note that by (10) λ ( x , S ) exists but is not necessarily unique. The coordinate bounds for a zonotope (or equivalently the inter- val hull containing that zonotope) are easily computed; for example , see [4, 10]. In particular , for x ∈ S , c j ( S ) − p ( S ) Õ i = 0   д j , i ( S )   ≤ x j ≤ c j ( S ) + p ( S ) Õ i = 0   д j , i ( S )   , (12) for all j = 1 , . . . , d , which we can write in compact form as c ( S ) − | G ( S ) | 1 p ( S ) ≤ x ≤ c ( S ) + | G ( S ) | 1 p ( S ) . (13) Lemma 3.1. Consider an interval hull B = { b ∈ R d | b ≤ b ≤ b } . For zonotope S ⊂ R d , the containment S ⊆ B is equivalent to the constraints c j ( S ) − p ( S ) Õ i = 0   д j , i ( S )   ≥ b j , c j ( S ) + p ( S ) Õ i = 0   д j , i ( S )   , ≤ b j (14) for all j = 1 , . . . , d . More compactly , c ( S ) − | G ( S ) | 1 p ( S ) ≥ b , c ( S ) + | G ( S ) | 1 p ( S ) ≥ b . (15) Proof. A straightforward conse quence of the bounds (12) and (13).  The class of zonotopes is close d under linear transformation [ 9 ], and the eect of a linear transform on a zonotope is easily imple- mented using the center-generator tuple representation M S = h M c ( S ) | MG ( S ) i , where S is a zonotope in dimension d and M ∈ R m × d is the matrix representing the linear transformation. Proposition 3.2. For systems without control inputs but with disturbance inputs constrained by the zonotope V = h c ( V ) | G ( V ) i , the center , generator count and generator matrix of the exact reach set R ( t , S ) are c ( R ( t , S )) = A t c ( S ) + t − 1 Õ s = 0 A t − 1 − s ( C c ( V ) + w ) , p ( R ( t , S )) = p ( S ) + t p ( V ) , G ( R ( t , S )) =  A t G ( S ) A t − 1 CG ( V ) · · · CG ( V )  (16) Proof. See [1, 10]  3.4 Computational Environment The examples in the paper were implemented in Ma tlab (R2018a) using CVX [ 12 , 13 ] ( version 2.1) with the SDPT3 solver ( version 4.0) for the convex optimizations and CORA [ 2 ] (2018 release) for zono- tope visualization. Timings were taken on a Lenovo ThinkPad X1 Y oga 1st Signature Edition laptop with an Intel Core i7-6500U (dual cor e) at 2.5 GHz with 8 GB memory running Windows 10 Pro (version 1803). 4 COMP U TING INV ARIAN T SETS As dened in (6), invariant sets are computed for systems without control input, so in this section we focus on the case where U = ∅ . In order to minimize notational complexity , we rst sketch the algorithm for systems with no disturbance input before proceeding to the more general case. 4.1 With No Inputs In this section we will derive conditions under which a parameter- ized zonotope I is invariant with respect to the state constraints (4). Dene I = h α | γ 1 д 1 ( I ) , γ 2 д 2 ( I ) , . . . , γ p ( I ) д p ( I ) ( I ) i , = h α | G ( I ) Γ i , (17) where Γ =          γ 1 0 · · · 0 0 γ 2 · · · 0 . . . . . . . . . . . . 0 0 · · · γ p ( I )          ∈ R p ( I ) × p ( I ) is the diagonal matrix with vector γ along its diagonal. In this parameterization, G ( I ) is a spe cied set of generators, but the center vector α ∈ R d x and generator scaling vector γ ∈ R p ( I ) with γ ≥ 0 are free parameters. Each element of γ is associated with a generator of I and can be thought of intuitively as the “width” of I in that generator direction. Proposition 4.1. Assuming a system with no disturbance or con- trol input, the reach set for an initial state space zonotope I can b e represented by a center , generator count and generator matrix given by c ( R ( t , I )) = A t α + t − 1 Õ s = 0 A t − 1 − s w p ( R ( t , I )) = p ( I ) G ( R ( t , I )) = A t G ( I ) Γ , (18) where the rst equation can b e written out elementwise as c j ( R ( t , I )) = d x Õ k = 1 ( A t ) j , k α k + t − 1 Õ s = 0 ( A t − 1 − s ) j , k w k ! (19) and the last equation can be written out for the i t h generator elemen- twise as д j , i ( R ( t , I )) = d x Õ k = 1 ( A t ) j , k д k , i ( I ) γ i . (20) Proof. Straightforward substitution of V = ∅ and (17) into (16).  preprint, 2019, arXiv .org Ian M. Mitchell, Jacob Budzis, and Andriy Bolyachevets Knowing how the reach set R ( t , I ) evolves allo ws us to ensure that trajectories which start within I stay within the constraint set X . Proposition 4.2. Assuming a system with no disturbance or con- trol input, I ⊆ Inv ( [ 0 , T ] , X ) if d x Õ k = 1 ( A t ) j , k α k + t − 1 Õ s = 0 ( A t − 1 − s ) j , k w k ! − p ( I ) Õ i = 0       d x Õ k = 1 ( A t ) j , k д k , i ( I )       γ i ≥ x j d x Õ k = 1 ( A t ) j , k α k + t − 1 Õ s = 0 ( A t − 1 − s ) j , k w k ! + p ( I ) Õ i = 0       d x Õ k = 1 ( A t ) j , k д k , i ( I )       γ i ≤ x j (21) for all j = 1 , . . . , d x and t = 0 , . . . , T . More compactly , A t α + t − 1 Õ s = 0 A t − 1 − s w −   A t G ( I )   γ ≥ x , A t α + t − 1 Õ s = 0 A t − 1 − s w +   A t G ( I )   γ ≤ x . (22) for all t = 0 , . . . , T . Proof. T o show that I ⊆ Inv ( [ 0 , T ] , X ) , by (6) we ne ed to show for x ( 0 ) ∈ I that x ( t ) ∈ X . If x ( 0 ) ∈ I then x ( t ) ∈ R ( t , I ) by (8). Plugging (19) and (20) into (14) yields (21), and so by (4) and Lemma 3.1, the constraint (21) implies R ( t , I ) ⊆ X for all t = 0 , . . . , T ; (23) consequently , x ( t ) ∈ X . Note that w e can move γ i outside the abso- lute value because γ i ≥ 0 . Rearranging (21) and taking advantage of the fact Γ 1 p ( I ) = γ yields (22).  Remark 4.3. The constraints (22) are linear in α and γ . Ideally , w e would then seek the set I of maximum volume which satises (22). Although an analytic formula for zonotope volume exists [ 11 ], it is combinatorially complex in the number of gener- ators and hence we settle for the simpler heuristic of maximizing the sum of the elements of γ (and thereby the sum of the “widths” of the generators). Our algorithm can therefore b e written as an optimization problem max α , γ 1 T p ( I ) γ such that γ ≥ 0 and (22) holds ∀ t ∈ 0 , . . . , T (24) Remark 4.4. The optimization (24) is a linear program with d x + p ( I ) decision variables, p ( I ) non-negativity constraints, and 2 d x ( T + 1 ) general constraints from (22). 4.2 With Disturbance Inputs For systems with uncertainty in the form of a disturbance input v ∈ V , ∅ , we must include the ee ct of V on R ( t , I ) when we construct the constraints for our optimization. -1 -0.5 0 0.5 1 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 2 -1 -0.5 0 0.5 1 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 2 -1 -0.5 0 0.5 1 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 2 -1 -0.5 0 0.5 1 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 2 Figure 1: Computed invariant set I (red thick line) for the rotation example with varying numb ers of generators in I . Also shown are the constraint set X (blue thick line) and R ( t , I ) (thin green lines) for t = 1 , . . . T . T op left: T wo gener- ators (the coordinate axes). T op right: Four generators (co- ordinate axes plus diagonals). Bottom left: Nine generators (equally spaced around top half circle). Bottom right: Six- teen random generators (only sev en have scaling γ i > 0 . 01 ). Proposition 4.5. Assuming a system with no control input, I ⊆ Inv ( [ 0 , T ] , X ) if ©       « A t α + t − 1 Õ s = 0 A t − 1 − s ( C c ( V ) + w ) −   A t G ( I )   γ − t − 1 Õ s = 0    A t − 1 − s CG ( V )   1 p ( V )  ª ® ® ® ® ® ® ¬ ≥ x , ©       « A t α + t − 1 Õ s = 0 A t − 1 − s ( C c ( V ) + w ) +   A t G ( I )   γ + t − 1 Õ s = 0    A t − 1 − s CG ( V )   1 p ( V )  ª ® ® ® ® ® ® ¬ ≤ x . (25) for all t = 0 , . . . , T . Proof. Plugging (16) into (15) demonstrates that (25) implies R ( t , I ) ⊆ X when V , ∅ . The remainder of the proof follows that of Proposition 4.2.  Note that the generators arising from the disturbance input appear in the constraints but are not scaled. The constraints (25) are still linear in α and γ , so we can compute an invariant set using the linear program optimization (24) with (25) substituted for (22). Invariant, Viability and Discriminating Kernel Under- Approximation via Zonotope Scaling preprint, 2019, arXiv .org -1 -0.5 0 0.5 1 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 2 Figure 2: Computed invariant set for the rotation example with disturbance input using eight equally spaced genera- tors for I . 4.3 Example: Rotational Dynamics T o demonstrate our algorithm and the eects of the choice of G ( I ) , we consider a system with rotational dynamics. Starting from the continuous time system Û x =  0 − 1 + 1 0  x we use the matrix e xponential with time step 0 . 2 to construct the discrete time system x ( t + 1 ) =  + 0 . 9801 − 0 . 1987 + 0 . 1987 + 0 . 9801  x ( t ) . (26) W e use the optimization (24) with x = − 1 2 , x = + 1 2 and T = 32 to compute an invariant set after slightly more than one full rotation. The true invariance kernel in this case is the largest circle contained in X . Figure 1 shows the computed invariant sets with dierent numbers of generators in I . Run times for the optimization were below 3 seconds for all of these examples. Figure 2 illustrates the results of the calculation if we introduce a disturbance with C =  1 0 0 1  V =   0 0       0 . 05 0 0 0 . 05   into the rotational dynamics (26). 5 COMP U TING VIABLE SETS Accommodating the disturbance input for invariant sets was no- tationally complicated but conceptually straightfor ward: The con- straints on the reach sets of the initial set had to take into account the eect of an input which could take on any value in V at any state at any time. Achieving viability typically requires that the con- trol input be chosen based on the curr ent state; consequently , we need to tackle the control inputs in a dierent manner . T o simplify the notation we work thr ough the algorithm without disturbance inputs in this section, and consequently compute viable sets. W e also omit the drift term for now . 5.1 A ugmenting Generators with Control Dene the function J (I , β , Φ ) =   α β       γ 1 д 1 ( I ) ϕ 1  ,  γ 2 д 2 ( I ) ϕ 2  , . . . ,  γ p ( I ) д p ( I ) ( I ) ϕ p ( I )  ,  , =   α β       G ( I ) Γ Φ   , (27) where β ∈ R d u , ϕ i ∈ R d u and Φ =  ϕ 1 ϕ 2 · · · ϕ p ( I )  ∈ R d u × p ( I ) . Before proceeding, we note: • The value of J (I , β , Φ ) is a zonotope in R d x + d u . • The set h β | Φ i is a zonotope in R d u , but J (I , β , Φ ) , I × h β | Φ i . • Unlike the diagonal matrices Γ (and Ψ encountered below ), the matrix Φ is dense: All entries may b e nonzero. • In the remainder of this w ork we will often refer to a col- lection of vectors and matrices { β ( s ) , Φ ( s )} t s = 0 . The fact that the colle ction is parameterized by time does not imply any direct temporal dependence between its elements. Proposition 5.1. Given { β ( s ) , Φ ( s )} t − 1 s = 0 and assuming a system with no disturbance input or drift, the reach set for an initial state space zonotope I can be represented by a center , generator count and generator matrix given by c ( R ( t , I )) = A t α + t − 1 Õ s = 0 A t − 1 − s B β ( s ) p ( R ( t , I )) = p ( I ) G ( R ( t , I )) = A t G ( I ) Γ + t − 1 Õ s = 0 A t − 1 − s B Φ ( s ) (28) Proof. For a system with V = ∅ and w = 0 , the dynamics (1) can be written as the linear transform x ( s + 1 ) =  A B   x ( s ) u ( s )  . And therefore  x ( 0 ) u ( 0 )  ∈ J (I , β ( 0 ) , Φ ( 0 )) ⇔ x ( 1 ) ∈  A B  J (I , β ( 0 ) , Φ ( 0 )) , ⇔ R ( 1 , I ) =  A B  J (I , β ( 0 ) , Φ ( 0 )) . Since a linear transformation of a zonotope is a zonotope, we can apply the linear transformation to (27) to determine c ( R ( 1 , I )) = A α + B β ( 0 ) , p ( R ( 1 , I )) = p ( I ) , G ( R ( 1 , I )) = AG ( I ) Γ + B Φ ( 0 ) , . The result (28) can then be derived by induction on s = 1 , . . . , T − 1 .  With this evolution formula for the initial zonotope I , we can deduce the necessar y constraints for viability . preprint, 2019, arXiv .org Ian M. Mitchell, Jacob Budzis, and Andriy Bolyachevets Proposition 5.2. Given { β ( t ) , Φ ( t )} T − 1 t = 0 and assuming a system with no disturbance input or drift, I ⊆ Viab ( [ 0 , T ] , X ) if ©       « A t α + t − 1 Õ s = 0 A t − 1 − s B β ( s ) −      A t G ( I ) Γ + t − 1 Õ s = 0 A t − 1 − s B Φ ( s )      1 p ( I ) ª ® ® ® ® ® ® ¬ ≥ x , ©       « A t α + t − 1 Õ s = 0 A t − 1 − s B β ( s ) +      A t G ( I ) Γ + t − 1 Õ s = 0 A t − 1 − s B Φ ( s )      1 p ( I ) ª ® ® ® ® ® ® ¬ ≤ x (29) and β ( t ) − | Φ ( t ) | 1 p ( I ) ≥ u , β ( t ) + | Φ ( t ) | 1 p ( I ) ≤ u (30) for t = 0 , . . . , T . Proof. Observe that x ( 0 ) ∈ I implies x ( t ) ∈ R ( t , I ) by (8). Dene the (not necessarily unique) input u ( t ) = β ( t ) + Φ ( t ) λ ( x ( t ) , R ( t , I ) ) . (31) By (11), u ( t ) ∈ h β ( t ) | Φ ( t )i . By (2) and Lemma 3.1, the con- straint (30) implies h β ( t ) | Φ ( t )i ⊆ U . Therefore, u ( t ) ∈ U for all t = 0 , . . . , T . By (4), (28) and Lemma 3.1, the constraint (29) implies R ( t , I ) ⊆ X and consequently x ( t ) ∈ X . W e have therefore proved for any x ( 0 ) ∈ I that there exists feasible input signal u ( ·) such that the trajector y x (·) starting from x ( 0 ) generated by u ( ·) satises x ( t ) ∈ X for t = 0 , . . . , T . By (7), I ⊆ Viab ( [ 0 , T ] , X ) .  By Proposition 5.2, nding a viable set I reduces to nding α , γ and { β ( t ) , Φ ( t )} T − 1 t = 0 to satisfy (29) and (30). W e can simply substitute (29) and (30) for (22) in the optimization (24) and add { β ( t ) , Φ ( t )} T − 1 t = 0 to the decision variables. A feasible solution to the resulting optimization problem will dene a set of viable states through α and γ , and a (time dependent) set of viable controls through { β ( t ) , Φ ( t )} T − 1 t = 0 ; however , the set of viable controls may be very small: By (31) the range of u ( t ) is directly proportional to | Φ ( t ) | , but reducing | Φ ( t ) | always makes it easier to satisfy (30) and the objective from (24) provides no direct incentive to increase | Φ ( t ) | (although a nonzero value may be necessar y to satisfy (29)). In order to achie ve a larger set of viable contr ols we would like to maximize | Φ ( t ) | , but maximizing an absolute value is a non-convex objective and we are not willing to destroy the convexity of our optimization to directly incorporate such a term. W e also cannot constrain Φ ( t ) ≥ 0 elementwise, since the sign of elements of ϕ i ( t ) relative to the sign of elements of the corresponding д i ( I ) may be critical to achieving viability (more discussion in section 5.3). Another approach is needed to seek broader control authority . 5.2 Control as Scaled Disturbance Intuitively , we w ould like to characterize the range of control input which could be applied while maintaining the viability of I . That control input authority might cause the reach set to be larger than it would be otherwise, but such growth may be accommodated in regions where the state constraints are not tight. Furthermore, we already have a method of mathematically characterizing the eect of such a priori indeterminant input: Tr eat it as a disturbance. With that in mind, we introduce a parameterized zonotope F to capture this control input authority which is independent of the authority in the equations above. Dene F = h 0 d u | ψ 1 д 1 ( F ) , ψ 2 д 2 ( F ) , . . . , ψ p ( F ) д p ( F ) ( F ) i , = h 0 d u | G ( F ) Ψ i , (32) where Ψ ∈ R p ( F ) × p ( F ) is the diagonal matrix with vector ψ along its diagonal. Like I , we will x the generators G ( F ) ∈ R d u × p ( F ) but leave the generator scaling vector ψ ∈ R p ( F ) with ψ ≥ 0 as a free parameter . W e then allow for time-dep endent F ( t ) and compute the result of input u ( t ) ∈ F ( t ) on the reach set according to Proposition 3.2. Proposition 5.3. Given { β ( s ) , Φ ( s ) , ψ ( s ) } t − 1 s = 0 and assuming a sys- tem with no disturbance input or drift, the reach set for an initial state space zonotope I can be represented by a center , generator count and generator matrix given by c ( R ( t , I )) = A t α + t − 1 Õ s = 0 A t − 1 − s B β ( s ) p ( R ( t , I )) = p ( I ) + t − 1 Õ s = 0 p ( F ( s ) ) G ( R ( t , I )) =  F I F 0 F 1 · · · F t − 1  (33) where F I = A t G ( I ) Γ + t − 1 Õ s = 0 A t − 1 − s B Φ ( s ) F s = A t − 1 − s BG ( F ( s ) ) Ψ ( s ) for s = 0 , . . . , t − 1 (34) Proof. Combine (16) and (28) by superposition.  W e note in passing that the equation for c ( R ( t , I )) in (33) demon- strates why it is sucient to x the center of F ( t ) at the origin in (32): The parameter β ( t ) already provides a mechanism to shift the center of the input set at time t . Proposition 5.4. Given { β ( t ) , Φ ( t ) , ψ ( t ) } T − 1 t = 0 and assuming a sys- tem with no disturbance input or drift, I ⊆ Viab ( [ 0 , T ] , X ) if ©             « A t α + t − 1 Õ s = 0 A t − 1 − s B β ( s ) −      A t G ( I ) Γ + t − 1 Õ s = 0 A t − 1 − s B Φ ( s )      1 p ( I ) − t − 1 Õ s = 0    A t − 1 − s BG ( F ( s ) )   ψ ( s )  ª ® ® ® ® ® ® ® ® ® ® ® ® ¬ ≥ x , ©             « A t α + t − 1 Õ s = 0 A t − 1 − s B β ( s ) +      A t G ( I ) Γ + t − 1 Õ s = 0 A t − 1 − s B Φ ( s )      1 p ( I ) + t − 1 Õ s = 0    A t − 1 − s BG ( F ( s ) )   ψ ( s )  ª ® ® ® ® ® ® ® ® ® ® ® ® ¬ ≤ x (35) Invariant, Viability and Discriminating Kernel Under- Approximation via Zonotope Scaling preprint, 2019, arXiv .org and β ( t ) − | Φ ( t ) | 1 p ( I ) − | G ( F ( t ) ) | ψ ( t ) ≥ u , β ( t ) + | Φ ( t ) | 1 p ( I ) + | G ( F ( t ) ) | ψ ( t ) ≤ u (36) for t = 0 , . . . , T . W e have shifted Ψ ( t ) outside the absolute value in (35) and (36) b e- cause ψ ( t ) ≥ 0 , and then taken advantage of the fact that Ψ ( t ) 1 p ( F ( t ) ) = ψ ( t ) . Proof. Let ˜ U ( t ) = h β ( t ) |  Φ ( t ) G ( F ( t ) ) Ψ ( t )  i , ˜ Φ ( t ) = h Φ ( t ) 0 d u × Í t − 1 s = 0 p ( F ( s ) ) i , and then choose any ρ ( t ) ∈ R p ( F ( t ) ) such that − 1 p ( F ( t ) ) ≤ ρ ( t ) ≤ + 1 p ( F ( t ) ) . (37) Dene the input u ( t ) = β ( t ) + ˜ Φ ( t ) λ ( x ( t ) , R ( t , I ) ) + G ( F ( t ) ) Ψ ( t ) ρ ( t ) . (38) W e augment ˜ Φ ( t ) with zero columns / generators to account for the extra generators (33) in R ( t , I ) arising from F ( s ) for s < t . Note that these extra generators have no dir ect eect on the choice of input for step t in (38) because they are zero vectors, although we do need to account for their continuing eect on the choice of λ ( x ( t ) , R ( t , I ) ) through these extra columns in ˜ Φ ( t ) . By (11), the fact that the extra generators in ˜ Φ ( t ) compared to Φ ( t ) are all zero vectors, and (37), u ( t ) ∈ ˜ U ( t ) . By (2) and Lemma 3.1, the constraint (36) implies ˜ U ( t ) ⊆ U . Therefore, u ( t ) ∈ U for all t = 0 , . . . , T . By (4), (33) and Lemma 3.1, the constraint (35) implies R ( t , I ) ⊆ X and conse quently x ( t ) ∈ X . The remainder is the same as the proof of Proposition 5.2.  Remark 5.5. The constraints (35) and (36) are linear in α and { β ( t ) , ψ ( t ) } T − 1 t = 0 , and are convex in γ and { Φ ( t )} T − 1 t = 0 . Our viability optimization problem is then written as max α , γ , { β ( t ) , Φ ( t ) , ψ ( t ) } T − 1 t = 0 1 T p ( I ) γ + η T − 1 Õ t = 0 1 T p ( F ( t ) ) ψ ( t ) such that γ ≥ 0 ψ ( t ) ≥ 0 and (35), (36) hold ∀ t ∈ 0 , . . . , T (39) where η > 0 is a weighting parameter used to trade o the relative importance of large γ (to encourage a larger set of viable states) against large ψ ( t ) (to encourage a larger set of viable controls). Remark 5.6. The optimization (39) is a convex program with d x + p ( I ) + T d u ( p ( I ) + 1 ) + Í t p ( F ( t ) ) decision variables, p ( I ) + Í t p ( F ( t ) ) non-negativity constraints, and 2 ( d x + d u )( T + 1 ) general convex constraints from (35) and (36). 5.3 Example: Double Integrator T o demonstrate the viability algorithm we consider the traditional double integrator . Starting from the continuous time system Û x =  0 + 1 0 0  x +  0 1  u -1 -0.5 0 0.5 1 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 2 -1 -0.5 0 0.5 1 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 2 Figure 3: Computed viable set I (red thick line) for the double integrator example with (left) and without (right) F . Eight generators e qually spaced in the north-west quadrant are used (only ve have scaling γ i > 0 . 01 ). Also shown are the constraint set X (blue thick line), R ( I , t ) (thin green lines) for t = 1 , . . . T , and a collection of sample traje ctories (star shows the initial conditions for each). 0 5 10 15 20 25 30 -1 0 1 x 1 0 5 10 15 20 25 30 -1 0 1 x 2 0 5 10 15 20 25 30 time step -1 0 1 u 1 +1.0 0 5 10 15 20 25 30 -1 0 1 x 1 0 5 10 15 20 25 30 -1 0 1 x 2 0 5 10 15 20 25 30 time step -1 0 1 u 1 +1.0 Figure 4: Sample viable trajectory (shown as the thick black trajectory in gure 3) with (left) and without (right) F . Con- trol at each time is chosen as close to + 1 as possible subje ct to (38) on the left and (31) on the right. The dotted cur ve in the control subplot shows the scaling ψ 1 ( t ) (which is always zero on the right). Any dierence between the actual control u 1 ( t ) and scaling ψ 1 ( t ) arises from Φ ( t ) . we use the matrix exponential with time step 0 . 1 and Ma tlab ’s integral() adaptive quadrature r outine to construct the discrete time system x ( t + 1 ) =  + 1 . 0000 + 0 . 1000 0 + 1 . 0000  x ( t ) +  + 0 . 0050 + 0 . 1000  u ( t ) . W e use the optimization (39) with d x = 2 , d u = 1 , x = − 1 d x , x = + 1 d x , u = − 1 , u = + 1 , and T = 30 to compute a viable set, and in the process demonstrate the importance of including characterizations of control authority from both sections 5.1 and 5.2. Figure 3 shows the computed viable set with and without the additional control authority enabled by the technique from section 5.2 (with F ( t ) equal to the 2 × 2 identity matrix for all t ). Run time for the optimization was below 5 seconds for both cases. No viable set is found if we omit Φ ( t ) from section 5.1, so that case is not shown. Some intuition for the failure of this latter case to nd any viable set can b e found by examining ϕ 1 ( t ) for the other cases. This input component corresponds to the state generator д 1 ( I ) =  0 1  T , which is the generator whose scaling to a large extent determines the height of the viable set. For the rst half of preprint, 2019, arXiv .org Ian M. Mitchell, Jacob Budzis, and Andriy Bolyachevets -1 -0.5 0 0.5 1 x 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 u 1 -1 -0.5 0 0.5 1 x 1 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 u 1 Figure 5: Projection of the set of viable inputs at t = 20 with (left) and without (right) F . the time horizon, ϕ 1 ( t ) < 0 ; in other words, states that hav e a large positive velocity x 2 ( t ) will be forced (through λ ( x ( t ) , R ( t , I ) ) in (31) or (38)) to choose a correspondingly large negativ e control input u ( t ) . Without the coupling of state and input signal through Φ ( t ) , viability is infeasible. Figure 4 shows a single trajector y’s states and control input componentwise over time and displays just this behavior: Although the desired input is + 1 , the early input signal u ( t ) is forced to be negative because x 2 ( t ) is quite positive. At later times ϕ 1 ( t ) does be come slightly positive, but it is at these times toward the end of the viability horizon that the b enets of section 5.2 are visible. Although the sets of viable states in gure 3 are the same with or without F , the set of viable controls is notably larger if we include a non-empty F ; for example, gure 5 shows that the set of viable controls at t = 20 is much larger on the left, while on the left of gure 4 u 1 ( t ) is much more able to approach the target value of + 1 in the latter half of the horizon b ecause of the growth of ψ 1 ( t ) . 6 COMP U TING DISCRIMINA TING SETS In this se ction we incorporate the control, disturbance and drift terms by combining the derivations in Sections 4.2 and 5.2 to han- dle the full generality of the dynamics (1). The longer equations referenced in this section are collected in gure 6. Proposition 6.1. Given { β ( s ) , Φ ( s ) , ψ ( s ) } t − 1 s = 0 , the reach set for an initial state space zonotope I can be represented by a center , generator count and generator matrix given by (40). Proof. Combine (16) and (33) by superposition.  Proposition 6.2. Given { β ( t ) , Φ ( t ) , ψ ( t ) } T − 1 t = 0 , if (36) and (41) hold then I ⊆ Disc ( [ 0 , T ] , X ) . Proof. Redene ˜ Φ ( t ) = h Φ ( t ) 0 d u × Í t − 1 s = 0 p ( F ( s ) ) 0 d u × t p ( V ) i , to account for the extra columns / generators in G ( R ( t , I )) aris- ing from the disturbance inputs, and then combine the proofs of Propositions 4.5 and 5.4 by superposition.  6.1 Example: Nonlinear Quadrotor T o demonstrate the utility of the full algorithm, we compute a discriminating set for a partially nonlinear six dimensional lon- gitudinal mo del of a quadrotor taken from [ 8 ]. The state space dimensions are: • horizontal position x 1 [m] (positive rightward), • vertical position x 2 [m] (positive upward), • horizontal velocity x 3 [m/s], • vertical velocity x 4 [m/s], • roll x 5 [rad] (positive clockwise), • roll velocity x 6 [rad/s]. The control input dimensions are: • total thrust u 1 , • desired roll angle u 2 . The nonlinear continuous time plant dynamics model is: Û x 1 = x 3 , (43a) Û x 2 = x 4 , (43b) Û x 3 = u 1 K sin x 5 , (43c) Û x 4 = − д + u 1 K cos x 5 , (43d) Û x 5 = x 6 , (43e) Û x 6 = − d 0 x 5 − d 1 x 6 + n 0 u 2 , (43f ) W e adopt the constraint set used in [ 17 ], except that we broaden the range of x 5 : x 1 ∈ [− 1 . 7 , + 1 . 7 ] , x 2 ∈ [ + 0 . 3 , + 2 . 0 ] , x 3 ∈ [− 0 . 8 , + 0 . 8 ] , x 4 ∈ [− 1 . 0 , + 1 . 0 ] , x 5 ∈ [− π 12 , + π 12 ] , x 6 ∈ [− π 2 , + π 2 ] . (44) W e also broaden the range of allow ed controls: u 1 ∈ [− 1 . 5 , + 1 . 5 ] + ¯ u 1 , u 2 ∈ [− π 12 , + π 12 ] + ¯ x 5 , (45) where ¯ u 1 = д / K and ¯ x 5 = 0 . W e then linearize (43) about ¯ u 1 and ¯ x 5 to arrive at the linear model (42). Note that unlike [ 17 ] we do not hybridize the dynamics in order to linearize ab out multiple operat- ing points; the entire constraint set is handled with a single linear model. For the range of x 5 in (44) and u 1 in (45), the linearization errors are error in Û x 3 ∈ [− 0 . 2760 , + 0 . 2760 ] , error in Û x 4 ∈ [ 0 . 0000 , + 0 . 3668 ] . For conservativeness, we choose V to be this error rectangle dilated by 10% . Note that the disturbance ae cts only x 3 and x 4 because the dynamics of the remaining dimensions are exactly linear . The parameter values use d in simulation are taken from [ 17 ]: K = 0 . 89 / 1 . 4 , d 0 = 70 , d 1 = 17 and n 0 = 55 . W e use Ma tlab ’s expm() and integral() routines to construct the discrete time version of (42) with time step 0 . 05 , and then com- pute for a horizon T = 40 . T o choose generators, we note that for ¯ x 5 = 0 the pairs of dimensions {( x 1 , x 3 ) , ( x 2 , x 4 ) , ( x 3 , x 5 ) , ( x 5 , x 6 )} look like double integrators. For these pairs we create ve gener- ators in the north-west quadrant, while for the remaining pairs of states we create only two generators along the diagonal and anti-diagonal. T o these we add the six coordinate axes for a total of p ( I ) = 48 generators. W e also tried running the optimization with an additional 52 randomly oriented generators to see whether coupling between more than two dimensions w ould improve the Invariant, Viability and Discriminating Kernel Under- Approximation via Zonotope Scaling preprint, 2019, arXiv .org c ( R ( t , I )) = A t α + t − 1 Õ s = 0 A t − 1 − s ( B β ( s ) + C c ( V ) + w ) , p ( R ( t , I )) = p ( I ) + t − 1 Õ s = 0 p ( F ( s ) ) + t p ( V ) G ( R ( t , I )) =  F I F 0 F 1 · · · F t − 1 A t − 1 CG ( V ) A t − 2 CG ( V ) · · · CG ( V )                       where F I and { F s } t − 1 s = 0 are given in (34) . (40) ©       « A t α + t − 1 Õ s = 0 A t − 1 − s ( B β ( s ) + C c ( V ) + w ) −      A t G ( I ) Γ + t − 1 Õ s = 0 A t − 1 − s B Φ ( s )      1 p ( I ) − t − 1 Õ s = 0    A t − 1 − s BG ( F ( s ) )   ψ ( s )  − t − 1 Õ s = 0    A t − 1 − s CG ( V )   1 p ( V )  ª ® ® ® ® ® ® ¬ ≥ x , ©       « A t α + t − 1 Õ s = 0 A t − 1 − s ( B β ( s ) + C c ( V ) + w ) +      A t G ( I ) Γ + t − 1 Õ s = 0 A t − 1 − s B Φ ( s )      1 p ( I ) + t − 1 Õ s = 0    A t − 1 − s BG ( F ( s ) )   ψ ( s )  + t − 1 Õ s = 0    A t − 1 − s CG ( V )   1 p ( V )  ª ® ® ® ® ® ® ¬ ≤ x                                        for t = 0 , . . . , T . (41) Figure 6: Some long equations for Propositions 6.1 and 6.2.             Û x 1 Û x 2 Û x 3 Û x 4 Û x 5 Û x 6             = A z }| {             0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 + K ¯ u 1 cos ¯ x 5 0 0 0 0 0 − K ¯ u 1 sin ¯ x 5 0 0 0 0 0 0 1 0 0 0 0 − d 0 − d 1             x z }| {             x 1 x 2 x 3 x 4 x 5 x 6             + B z }| {             0 0 0 0 K sin ¯ x 5 0 K cos ¯ x 5 0 0 0 0 n 0             u z }| {  u 1 u 2  + C z }| {             0 0 0 0 1 0 0 1 0 0 0 0             v z }| {  v 1 v 2  + w z }| {             0 0 − K ¯ u 1 sin ¯ x 5 + K ¯ u 1 cos ¯ x 5 − д 0 0             (42) Figure 7: Linearize d quadrotor dynamics. discriminating set; however , none of the randomly oriente d gen- erators achieved a scaling γ i ≥ 0 . 01 and so we discar ded them in subsequent runs. W e run the optimization problem (39) with constraint (41) substi- tuted for (35). Run time for the optimization is just over 5 minutes. Only 11 of the 48 generators achieved a scaling factor γ i ≥ 0 . 01 . Figure 8 shows projections of the resulting discriminating set. Note that this set is discriminating for both the linear and nonlinear models, since we conservatively capture all linearization error in the disturbance input bounds. Although we hav e insucient space to make a detaile d compari- son with the ellipsoid-based results from [ 17 ], we can observe that in roughly the same computational time (albeit on a slightly faster laptop) our zonotope-based algorithm is able to nd a much larger discriminating set over twice the time horizon with double the time resolution. Furthermore, the zonotope representation is much less conservative with its treatment of the disturbance input, and hence we do not need to hybridize or to restrict the range of x 5 and u 1 so severely . W e hyp othesize that most of the improv ement in accuracy arises from the fact that zonotopes can exactly r epresent the rect- angular form of typical state and input constraints, while ellipsoids are forced to adopt dramatic under- or over-appr oximations. -2 0 2 x 1 0 0.5 1 1.5 2 x 2 -2 0 2 x 1 -1 -0.5 0 0.5 1 x 3 0 1 2 x 2 -1 -0.5 0 0.5 1 x 4 -0.2 0 0.2 x 5 -2 -1 0 1 2 x 6 Figure 8: Projections of the computed discriminating set I (red thick line) for the quadrotor model. preprint, 2019, arXiv .org Ian M. Mitchell, Jacob Budzis, and Andriy Bolyachevets 7 DISCUSSION The formulations developed above leav e a number of parameters to be chosen by the user and make some assumptions which could be relaxed. W e briey discuss these issues in this section. The key parameter open to the user is the choice of generators G ( I ) and G ( F ( t ) ) . Unless there is some reason to b elieve that direct coupling of the inputs is benecial, the latter will typically be chosen as an identity matrix. Choosing the former is trickier; howev er , the experience in section 6.1 indicates that examination of the sparsity pattern of A may allow one to choose these vectors more eciently than simply trying to cover the unit hypersphere in R d x . One factor that is perhaps not so obvious when choosing gener- ators is that this choice impacts the quality of the resulting sets not just directly through the generators but also indirectly through the heuristic obje ctive function in (39). If G ( I ) is orthonormal then max 1 T p ( I ) γ is not an unreasonable heuristic to make I large, but as the generators lose perpendicularity (inevitable as the number of generators grows) the quality of this heuristic decreases. It should be noted that this heuristic appears to encourage sparse γ , which could be a signicant benet for downstream uses of I . Finally , the restriction to interval hulls of the control input set U and constraint set X was driven by the need to include constraints in the optimization which conrmed that (projections of ) zonotopes were contained within those sets. That restriction can be relaxed to any class of sets for which a reasonable number of constraints can conrm zonotope containment; for example, intersections of slabs, or even conv ex polygons with a modest number of faces. It may even be possible to allow full zonotopes using [ 14 , Lemma 3], albeit at the cost of swapping a small number of constraints (such as (36)) for a full linear matrix inequality . 8 CONCLUSIONS AND F U T URE W ORK W e have derived conve x optimizations whose solutions represent invariant, viable or discriminating sets for discrete time, continuous state ane dynamical systems, and demonstrated the results on a simple rotation, a double integrator and a six dimensional nonlinear longitudinal model of a quadrotor respectively . The optimizations can be solv ed with modest computing p ower in seconds to minutes. A key shortcoming of the curr ent formulation is its restriction to discrete time. Unfortunately , the typical approaches to soundly mapping continuous time reachability into discrete time reachabil- ity (for example , see [ 3 ]) cannot be used in our appr oach, so we are exploring alternatives. W e are hopeful that some of the computational eciency tech- niques demonstrated in [ 6 , 7 ] could be applied to improve the scalability of the optimization problems derived here. Although there are pruning heuristics which could be applied to reduce the number of generators, the size of the optimization inevitably grows with ner time discretization and/or longer horizons. W e have con- sidered approaches to replace one large optimization with many smaller ones, although there are tradeos in such schemes. Finally , we are exploring the use of the viable and discriminating sets for classifying and ltering exogenous control inputs at run- time. A CKNO WLEDGMEN TS This work was supp orted by National Science and Engineering Research Council of Canada (NSERC) Undergraduate Student Re- search A wards and Discovery Grant #298211. REFERENCES [1] Matthias Altho. 2010. Reachability Analysis and its Application to the Safety Assessment of Autonomous Cars . P h.D. Dissertation. [2] Matthias Altho. 2015. An introduction to CORA 2015. In In Proc. of the W orkshop on A pplied V erication for Continuous and Hybrid Systems . 120–151. [3] Matthias Altho and Goran Frehse. 2016. Combining zonotopes and support functions for ecient reachability analysis of linear systems. In IEEE Conference on Decision and Control (CDC) . 7439–7446. [4] Matthias Altho and Bruce H. Krogh. 2011. Zonotope Bundles for the Ecient Computation of Reachable Sets. In IEEE Conference on De cision and Control (CDC) . 6814–6821. [5] Jean-Pierre A ubin, Alexandre M. Bayen, and Patrick Saint-Pierre. 2011. Viability Theory: New Directions . Springer . https://doi.org/10.1007/978- 3- 642- 16684- 6 [6] Stanley Bak and Parasara Duggirala. 2017. HyLAA: A T ool for Computing Simulation-Equivalent Reachability for Linear Systems. In Hybrid Systems: Com- putation and Control (HSCC) . [7] Sergiy Bogomolov , Marcelo Forets, Goran Frehse, Frédéric Viry, Andreas Podelski, and Christian Schilling. 2018. Reach Set Appro ximation Through Decomposition with Low-dimensional Sets and High-dimensional Matrices. In Hybrid Systems: Computation and Control (HSCC) . 41–50. https://doi.org/10.1145/3178126.3178128 [8] Patrick Bouard. 2012. On-board Model Predictive Control of a Quadrotor Heli- copter: Design, Implementation, and Experiments . T echnical Report UCB/EECS- 2012-241. Department of Electrical Engineering and Computer Science, Univer- sity of California at Berkeley . http://w ww .eecs.berkeley .edu/Pubs/T echRpts/ 2012/EECS- 2012- 241.html [9] Antoine Girard. 2005. Reachability of Uncertain Linear Systems using Zonotopes. In Hybrid Systems: Computation and Control (HSCC) , Manfred Morari and Lothar Thiele (Eds.). Number 3414 in Lecture Notes in Computer Science. Springer V erlag, 291–305. https://doi.org/10.1007/978- 3- 540- 31954- 2_19 [10] Antoine Girard, Colas Le Guernic, and Oded Maler . 2006. Ecient Computation of Reachable Sets of Linear Time-Invariant Systems with Inputs. In Hybrid Systems: Computation and Control (HSCC) , João Hespanha and Ashish Tiwari (Eds.). Numb er 3927 in Lecture Notes in Computer Science. Springer V erlag, 257–271. [11] Eugene Gover and Nishan Krikorian. 2010. Determinants and the volumes of parallelotopes and zonotopes. Linear Algebra A ppl. 433, 1 (2010), 28 – 40. https://doi.org/10.1016/j.laa.2010.01.031 [12] Michael Grant and Stephen Boyd. 2008. Graph implementations for nonsmooth convex programs. In Recent Advances in Learning and Control , V . Blondel, S. Boyd, and H. Kimura (Eds.). Springer- V erlag Limited, 95–110. http://stanford.edu/ ~boyd/graph_dcp.html. [13] Michael Grant and Stephen Boyd. 2017. CVX: Matlab Software for Disciplined Convex Programming, version 2.1. http://cvxr .com/cv x. (Dec. 2017). [14] Dongkun Han, Albert Rizaldi, Ahmed El-Guindy , and Matthias Altho. 2016. On enlarging backward reachable sets via Zonotopic set membership. In IEEE Int. Symp. on Intelligent Control (ISIC) . 1–8. [15] John N. Maidens, Shahab Kaynama, Ian M. Mitchell, Me eko M. K. Oishi, and Guy A. Dumont. 2013. Lagrangian Methods for Approximating the Viability Kernel in High-Dimensional Systems. A utomatica 49, 7 (July 2013), 2017–2029. https://doi.org/10.1016/j.automatica.2013.03.020 [16] Ian M. Mitchell and Shahab Kaynama. 2015. An Improved Algorithm for Robust Safety Analysis of Sampled Data Systems. In Hybrid Systems: Computation and Control (HSCC) . 21–30. https://doi.org/10.1145/2728606.2728619 [17] Ian M. Mitchell, Jerey Y eh, Forrest J. Laine, and Claire J. T omlin. 2016. Ensuring Safety for Sampled Data Systems: An Ecient Algorithm for Filtering Potentially Unsafe Input Signals. In IEEE Conference on Decision and Control (CDC) . Las V egas, NV , 7431–7438. https://doi.org/10.1109/CDC.2016.7799417 [18] Bastian Schürmann and Matthias Altho. 2017. Convex interpolation control with formal guarantees for disturbed and constrained nonlinear systems. In Hybrid Systems: Computation and Control (HSCC) . 121–130. https://doi.org/10. 1145/3049797.3049800 [19] Bastian Schürmann and Matthias Altho. 2017. Guaranteeing constraints of disturbed nonlinear systems using set-base d optimal control in generator space. In International Federation of Automatic Control (IFA C) W orld Congress . [20] Bastian Schürmann and Matthias Altho. 2017. Optimal control of sets of solu- tions to formally guarantee constraints of disturbed linear systems. In A merican Control Conference (A CC) . 2522–2529. https://doi.org/10.23919/A CC.2017.7963332