A sampling-based approach to scalable constraint satisfaction in linear sampled-data systems—Part I: Computation ∗  Shahab Kaynama, Jeremy H. Gillula, Claire J. Tomlin †  (Preprint Submitted for Publication)  Abstract  Sampled-data (SD) systems, which are composed of both discrete- and continuous-time components, are arguably one of the most common classes of cyberphysical systems in practice; most modern controllers are implemented on digital platforms while the plant dynamics that are being controlled evolve continuously in time.   As with all cyberphysical systems, ensuring hard constraint satisfaction is key in the safe operation of SD systems. A powerful analytical tool for guaranteeing such constraint satisfaction is the viability kernel: the set of all initial conditions for which a safety-preserving control law (that is, a control law that satisfies all input and state constraints) exists. In this paper we present a novel sampling-based algorithm that tightly approximates the viability kernel for high- dimensional sampled-data linear time-invariant (LTI) systems. Unlike prior work in this area, our algorithm formally handles both the discrete and continuous characteristics of SD systems. We prove the correctness and convergence of our approximation technique, provide discussions on heuristic methods to optimally bias the sampling process, and demonstrate the results on a twelve-dimensional flight envelope protection problem.  1   Introduction  The mathematical guarantee of satisfaction of hard input and state constraints is an increasingly desirable property that every safety-critical cyberphysical system must implement. The subset of the state space for which this property holds is known as the viability kernel [2], or alternatively, in the infinite horizon case, the maximal controlled-invariant set [3]. Consequently, a tremendous amount of work in the recent literature has been focused on methods for computing the viability kernel. Constrained sampled-data (SD) systems describe a large class of cyberphysical systems. In most practical settings the system evolves continuously in time but the state is measured only at discrete time instants [4]. Consequently, any admissible control policy is piecewise constant; the input can only be applied at the beginning of each sampling interval and is kept constant (under zero-order hold) until the next sampling time. Examples of SD systems include control of blood glucose levels in type-1 diabetes [5] where insulin could only be administered every   30   minutes during the clinical trials. Additionally, restrictions on sampling frequency is not always due to sensory limitations or the particulars of an application.   For instance, in model predictive control (MPC) a higher sampling frequency results in a significantly larger online optimization problem over each prediction horizon—a known limiting factor in embedded control design, e.g. in the automotive industry. Furthermore, without additional rate constraints, a higher sampling frequency would require controllers with much higher bandwidth. Discretizing the dynamics and naively designing control policies in discrete time ignores the behavior of the true system in between two sampling instants. This inter-sample behavior can be crucial in safety-critical systems where the state constraint is associated with “safety” of the system. At the same time, performing continuous-time safety analysis on the system (for example, using the traditional level-set technique [6]), which only requires a mild assumption of  ∗ An earlier version of this paper was partially presented at the 17th Int’l Conference on Hybrid Systems: Computation and Control, April 15-17, 2014, Berlin, Germany [1].  † The authors are with Electrical Engineering & Computer Sciences, University of California at Berkeley, 337 Cory Hall, Berkeley, CA 94720, USA.   { kaynama,   jgillula,   tomlin } @eecs.berkeley.edu  1  arXiv:1405.2363v1 [cs.SY] 9 May 2014
Lebesgue measurability of the control input, cannot provide guarantees about the behavior of the SD system where the input is restricted to draw from the class of piecewise constant signals [7]. All this warrants a technique that can formally handle SD systems. We present a sampling-based approach that yields tight under- and over-approximations of the viability kernel for  high-dimensional SD   linear time-invariant (LTI) dynamics.  1.1   Related Work  The classical numerical schemes for approximating the viability kernel are those based on gridding the state space (and discretizing the dynamics) and appropriately evolving the constraints over this stationary grid. Such schemes include the level-set methods [6] and variants of Saint-Pierre’s viability algorithm [8]. Despite their versatility in handling complex dynamics and sets, the applicability of these schemes has historically been limited to systems of dimensions less than five. Efforts to generalize such grid-based techniques to moderately dimensioned systems include structure decomposition [9, 10] and approximate dynamic programming [11]. Algorithms from within the MPC community have also emerged that enable the computation of the viability kernel for discrete-time LTI systems with polytopic constraints [3]. Due to the fact that these algorithms recursively compute the Minkowski sum, linear transformation, and intersection of polytopes, they can only be applied to low dimensional systems; the number of vertices of the resulting polytope grows rapidly with each subsequent Minkowski sum operation, while the intersection operation at each iteration requires a vertex to facet enumeration of polytopes— an operation that is known to be intractable in high dimensions [12]. In more general contexts (e.g., for continuous-time systems), an ellipsoidal approximation of the region of attraction of the MPC is computed as a (crude) representation of the viability kernel [13]. These, as well as other approximation techniques such as   β -contractive polytopes [14], generally require existence of a stabilizing controller within the constraints—a requirement that may not always be readily satisfied. For LTI systems with convex constraints, [15] (discrete-time) and [16] (discrete- and continuous-time) introduced efficient and scalable algorithms based on support vector representations and piecewise ellipsoidal sets to conser- vatively approximate the viability kernel.   These algorithms follow the flow of the dynamics in the same spirit as Lagrangian techniques for maximal reachability such as [17–19]. Approximating the viability kernel can also be viewed as a search for an appropriate control Lyapunov function subject to additional input and state constraints.   As such, in parallel to the above developments, for polynomial systems with semi-algebraic constraints, various sum-of-squares (SOS) 1   optimization-based techniques have been proposed [20–24] that either directly form a polynomial approximation of the viability kernel, or can be modified to do so. The resulting bilinear SOS program is solved either by alternating search (prone to local optima) or through convex relaxations (introduces additional conservatism). The degree of the SOS multipliers, which directly translates to the presumed degree of the polynomial that is to describe the kernel, is commonly kept low (e.g. quartic), striking a tradeoff between excessive conservatism and computational complexity. A related SOS-based technique is the recent method of occupation measures [25, 26] that, while scalable, can only   over -approximate the desired set (which is insufficient for safety). Though the approximation recovers the true kernel in the limit, the error is not monotonically decreasing with the degree of the multipliers. All of the above algorithms are designed for either discrete- or continuous-time dynamics, and very little effort has been dedicated to   sampled-data   system—precisely the systems that, due to their physically descriptive and realistic nature, stand to benefit the most from safety formalism. To the best of our knowledge, the only other scalable work on computing the viability kernel for SD systems is presented in [7], where a piecewise ellipsoidal algorithm is proposed for LTI dynamics with ellipsoidal constraints. Unfortunately due to the projections and cross-products involved in this algorithm, the quality of the approximation degenerates rapidly with the time horizon.  1.2   Summary of Contributions  We propose a sampling-based approach that directly handles SD systems, albeit under LTI dynamics and convex constraints. The algorithm provably under-approximates the viability kernel in a scalable and efficient manner. Since  1 SOS is a relaxation of the original semidefinite program.  2
the proposed method deals with individual trajectories to infer the evolution of sets of initial conditions, it is not explicitly restricted by a particular shape of the constraints, although its computational complexity does depend on these shapes. Our technique yields a “tight” approximation in the sense that the resulting set touches the boundary of the true viability kernel from the inside with arbitrary precision up to a numerical constant. The algorithm is sampling based, meaning that the points to be included as boundary points of the viability kernel are sampled according to a probability distribution.   We provide a convergence proof that shows that such a sampling-based technique is superior to any alternative deterministic approach. Section 2 formulates the problem we wish to address. Section 3 presents our main under-approximation algorithm, while Section 4 proves its correctness and convergence.   Section 5 discusses the computational complexity of the technique and showcases, via an example, its scalability. To provide a measure of conservatism of our technique, we will also present an over-approximation algorithm in Section 6. This over-approximation algorithm is then utilized in Section 7 to help guide the under-approximation process. This section also describes a few heuristics to bias the random sampling of the state constraint so as to achieve superior convergence and accuracy properties. We study a 12D flight envelope protection problem in Section 8, before providing concluding remarks and future directions in Section 9.  2   Problem Formulation  Consider the LTI system   ̇ x ( t ) =   Ax ( t ) +   Bu ( t )   (1) with state   x ( t )   ∈   R n , control input   u ( t )   ∈ U , where   U   is a compact convex subset of   R m .   A   and   B   are constant matrices of appropriate dimension. The state of the system is measured at every time instant   t k   :=   kδ   for   k   ∈   Z ≥ 0   and fixed sampling interval   δ   ∈   R > 0 . We are concerned with the evolution of the system over   T   := [0 , τ   ]   with arbitrary, finite time horizon   τ   ∈   R > 0 . We denote by   N δ   :=   d τ /δ e   the number of sampling intervals in   T . The input is applied at the beginning of each sampling interval and is kept constant until the next sampling instant. Thus the input signal draws from the set of piecewise constant functions  U   pwc  T   :=   { u   :   T   →   R m   piecewise const. , u ( t k )   ∈ U ∀ k, u ( t ) =   u ( t k )   ∀ t   ∈   [ t k , t k +1 ) } .   (2) For every   x 0   ∈   R n   and   u ( · )   ∈   U   pwc  T   , we denote the (unique) trajectory of (1) by   x u x 0   :   T   →   R n   with initial condition  x u x 0   (0) =   x 0 . When clear from the context, we shall drop the subscript and superscript from the trajectory notation. For a nonempty compact convex state constraint set   K ⊂   R n , deemed   safe , we examine the following construct:  Definition 1   (SD Viability Kernel) .   The finite-horizon SD viability kernel of   K   is the set of all initial states for which there exists a control law such that the trajectories emanating from those states remain in   K   over   T :  Viab sd  T   ( K ) :=   { x 0   ∈ K | ∃ u ( · )   ∈   U   pwc  T   ,   ∀ t   ∈   T , x u x 0   ( t )   ∈ K } .   (3) We seek to find a scalable technique that   under-approximates   the viability kernel. (Any approximation must be conservative, in that at the very least all states for which no admissible input can maintain safety must be excluded.)  2.1   Preliminaries and Notation  We say that the vector field in (1) is bounded on   K   if   ∃ M >   0   such that   ‖ Ax   +   Bu ‖ p   ≤   M   ∀ ( x, u )   ∈ K × U   for some norm   ‖·‖ p . If   K   and   U   are compact, every continuous vector field is bounded on   K . The   ‖·‖ p -distance of a point  x   ∈   R n   from a nonempty set   A ⊂   R n   is defined as   dist p ( x,   A ) := inf a ∈A   ‖ x   −   a ‖ p . The   Minkowski sum   of any two nonempty subsets   A   and   C   is   A ⊕ C   :=   { a   +   c   |   a   ∈ A , c   ∈ C} ; their   Pontryagin difference   (or, the   erosion   of   A   by   C ) is   A   C   :=   { a   |   a   ⊕ C ⊆ A} . We denote by   ∂ C   the boundary of the set   C , by   C c  its complement, and by   vol( C )   its volume.   conv( { v 0 , . . . , v N   } )   denotes the convex hull of a set of points   v 0 , . . . , v N   .  B n p   ( x, a )   denotes the closed   p -norm ball of radius   a   ∈   R ≥ 0   about a point   x   in   R n , and   S n − 1  p   the codimension one boundary   ∂ B n p   (0 ,   1) . A   ray   in   R n   is the set of points   ~ r   =   { r 0   +   sr d   |   s   ∈   R ≥ 0 } , where   r 0   ∈   R n   is the   origin   of the ray, and   r d   ∈   R n   is a unit vector giving the   direction   of the ray. 3
v 0  K  ~ r b r d  (a) Sample a point   r d   ∈ S n − 1 2   and form   ~ r . Find  b   =   ~ r   ∩   ∂ K .  v 0  K  ~ r v 1  b (b) Perform bisection search on   { v 0 , b }   to find  v 1   ∈   Viab sd  T   ( K ) .  v 0  K v 2  v 1  v 7  v 6   v 5  v 4  v 3 (c) Repeat for   N   samples.  K V 7  Viab sd  T   ( K ) (d) Return   conv( { v i } N i =0 ) .  Figure 1: Illustration of Algorithm 1.  3   Methodology  3.1   Overview of the Algorithm  Before we describe our algorithm we first elaborate on some simple subroutines which we make use of, but which we will not formally define due to space constraints. • F IND -I NTERSECTION - ON -B OUNDARY ( C , ~ r )   – Input: A convex compact set   C ⊂   R n , and a ray   ~ r   with origin  r 0   ∈ C . Returns a point   x   along the ray   ~ r   on   ∂ C . We note that since   C   is convex and compact, and since   r 0   ∈ C , there is exactly one such point   x . Runs in time proportional to the number of faces of   C   if   C   is a polytope, and constant time if   C   is an ellipsoid. • S AMPLE -R AY ( x )   – Input: A point   x   ∈   R n . First samples a point   r d   at random from   S n − 1 2   . Returns the ray  ~ r   =   { x   +   sr d } . Runs in time linear in   n . • F EASIBLE ( x,   C )   – Input: A point   x   ∈   R n   and a convex compact set   C ⊂   R n . Returns true if   x   ∈   Viab sd  T   ( C ) , and false otherwise. Further details on the implementation of this subroutine and its time complexity are given in Sections 3.2 and 5. 4
The under-approximation of the viability kernel proceeds as follows. We assume as input a description of   K , 2   as well as some initial point   v 0   ∈   Viab sd  T   ( K ) . 3   We then construct a polytopic approximation of   Viab sd  T   ( K )   by iteratively sampling a direction   r d   and generating a ray   ~ r   centered at   v 0   (Algorithm 1, step 4); finding the point   b   ∈   R n   where  ~ r   intersects the boundary of   K   (Algorithm 1, step 5); and then performing a bisection search along the line segment  { v 0 , b }   until we find a point   v i   such that   v i   ∈   Viab sd  T   ( K )   and   dist p ( v i , ∂   Viab sd  T   ( K ))   <    for some desired accuracy   >   0   in some norm   ‖·‖ p   (Algorithm 2). By repeating for   N   samples and taking the convex hull of the resulting points   { v 0 , . . . , v N   }   we arrive at a polytope   V N   ⊆   Viab sd  T   ( K )   which converges to   Viab sd  T   ( K )   (in a manner we shall formalize in Section 4.2). The conservatism of the algorithm will be formally proven in Section 4.1. Fig. 1 illustrates the under-approximation procedure.  Algorithm 1   Computes a polytopic under-approximation of   Viab sd  T   ( K )   with at most   N   + 1   vertices  1:   procedure   P OLYTOPIC -A PPROX ( K , v 0 , N   )  2:   V 0   ← { v 0 }  3:   for   i   = 1   to   N   do   . N   samples  4:   ~ r   ←   S AMPLE -R AY ( v 0 )  5:   b   ←   F IND -I NTERSECTION - ON -B OUNDARY ( K , ~ r )  6:   v i   ←   B ISECTION -F EASIBILITY ( v 0 , b,   K )  7:   V i   ←   conv( V i − 1   ∪ { v i } )  8:   end for  9:   return   V N  10:   end procedure Algorithm 2   Determines an    -accurate intersection of   ∂   Viab sd  T   ( K )   and the line between   a   and   b .  1:   function   B ISECTION -F EASIBILITY ( a, b,   K )  2:   c   ←   a   + ( b   −   a ) / 2  3:   if   F EASIBLE ( c,   K )   then  4:   if   dist p ( b, c )   <    then  5:   return   c  6:   else  7:   return   B ISECTION -F EASIBILITY ( c, b,   K )  8:   end if  9:   else  10:   return   B ISECTION -F EASIBILITY ( a, c,   K )  11:   end if  12:   end function  3.2   Checking Point Feasibility  The key to our approach is the subroutine F EASIBLE ( x 0 ,   K )   in Step 3 of Algorithm 2, which returns true only if  x 0   ∈   Viab sd  T   ( K ) . An overview of the procedure is as follows: Given an initial condition   x 0 , we verify the existence of a piecewise constant control law that ensures that the trajectory starting from   x 0   belongs to a subset of   K   at every sampling instant. This subset is appropriately chosen at a certain distance from the exterior   K c   such that the inter- sampling portions of the trajectory do not escape   K . The state   x 0   is then labeled as   feasible . We employ forward simulation to determine feasibility.   To do this in a tractable fashion, we use a finite-difference approximation of  2 In theory,   K   can be of any arbitrary (but convex) shape. In practice, this shape directly affects the run time complexity of the subroutines F IND -I NTERSECTION - ON -B OUNDARY   and F EASIBLE   as discussed in the bullet points above.  3 For linear systems if   (0 ,   0)   ∈ K × U   then   0   ∈   Viab sd  T   ( K ) , and thus we can use   v 0   = 0 . Otherwise, we assume we can sample points at random over   K   until we find a point   x   via the subroutine F EASIBLE ( x,   K )   such that   x   ∈   Viab sd  T   ( K ) .  5
0   0.5   1   1.5   2   2.5   3  −1.5  −1  −0.5  0  0.5  1  1.5  Position [m] Velocity [m/s]   Continuous−time  Discrete−time ( δ =0.5)  Discrete−time ( δ =0.2) Figure 2: Simulated trajectory of a mass-spring-damper (unit parameters) under constant force in continuous time vs. its Euler discretization for two sampling times   δ   = 0 . 2 ,   0 . 5 .  x 0  K  x ( t 1 )  K ↓ ( M, δ ) Figure 3: Erosion of   K   ensures that the curvature of the continuous trajectory in between two sampling instants cannot escape safety. the dynamics. Therefore, we also need to take into account the effect of the discretization error (Fig. 2 depicts the significance of this error via a trivial example).  3.2.1   Dealing With Inter-Sampling Behavior  The following lemma ensures that if a control law exists that can keep the trajectory value evaluated at every sampling instant in a certain subset of   K , then the continuous evolution of the system in between sampling instants also maintains safety; that is, the curvature of the trajectory during the sampling intervals does not escape   K   (Fig. 3).  Lemma 1.   Let   (1)   be uniformly bounded on   K   by   M >   0   in some norm   ‖·‖ p 1   . With a sampling interval   δ >   0   define  K ↓ ( M, δ ) :=   { x   ∈ K |   dist p 1   ( x,   K c )   ≥   M δ }  =   K   B n p 1   (0 , M δ ) ,   (4)  For a given initial condition   x 0 , if   ∃ u ( · )   ∈   U   pwc  T   such that   x ( t k )   ∈ K ↓ ( M, δ )   ∀ k   ∈ { 0 , . . . , N δ   } , then   x 0   ∈  Viab sd  T   ( K ) .  The proof is similar to the continuous-time analysis in [16, Proof of Prop. 1] and is omitted here; see [1, Lem. 1] for more detail.  3.2.2   Dealing With Discretization Error  Recall that the solution   x k +1   =   x ( t k +1 )   of the SD system (1) at time   t k +1   for any sampling interval   [ t k , t k   + δ ]   starting from   x k   =   x ( t k )   using a constant input   u k   is  x k +1   =   e Aδ   x k   +  (∫   δ  0  e Aλ dλ   ·   B  )  u k .   (5) 6
x 0 K   ˆ x 1  x 1  e 1  x 2  ˆ x 2  e 2  Figure 4: The mismatch between the true state   x k   at sampling time   t k   and the nominal state   ˆ x k   of the discretized model. The error   e k   propagates in time. If not accounted for in forward simulations, this mismatch could jeopardize safety. Consider the Taylor series expansion   e As   =   ∑ ∞  i =0 ( As ) i /i ! . To approximate the evolution of (1) at every discrete time instant   t k   using a finite-difference equation, we approximate the above infinite sum by a   ζ th order finite sum  A ζ,s   :=   ∑ ζ i =0  ( As ) i  i !   ≈   e As ,   ζ <   ∞   (6) with truncation error  E ζ,s   :=   e As   −   A ζ,s   =   ∑ ∞  i = ζ +1  ( As ) i  i !   .   (7) The trajectory of the resulting finite-difference equation is  ˆ x k +1   =   A ζ,δ   ˆ x k   +  (∫   δ  0  A ζ,λ dλ   ·   B  )  u k .   (8) We refer to   ˆ x   as the   nominal   state of the system. It is important to emphasize that the computation of the matrix exponential has always been a challenging task [27]. An exact computation is generally not possible. Instead, approximate techniques are employed. For instance, Matlab uses the Pad ́ e approximation with scaling and squaring method of [28]. In most practical cases, the forward Euler approximation ( ζ   = 1 ) is used. Truncating the tail   E   of the Taylor series expansion introduces a discretization error that results in mismatch be- tween the true values of the system trajectory at discrete time instants and the values generated by the finite-difference model. That is, the   discretization error   e k   at time   t k   resulting from the truncation error   E   will cause a deviation from the true state   x k   =   x ( t k )   such that  ˆ x k   =   x k   −   e k .   (9) Consequently, to formally guarantee an under-approximation of the SD viability kernel using a time discretized approach—i.e. via simulation of the nominal trajectory—we must take into account the effect of the error   e   and its forward propagation in time (Fig. 4).  Lemma 2.   Suppose that there exist constants   γ k   ≥   0   and a norm   ‖·‖ p 2   such that   ‖ e k ‖ p 2   ≤   γ k   ∀ k . Let  K k  ↓   ( M, δ, γ k ) :=   K ↓ ( M, δ )   B n p 2   (0 , γ k )   6   =   ∅ .   (10)  Then we have that  ˆ x k   ∈ K k  ↓   ( M, δ, γ k )   ⇒   x k   ∈ K ↓ ( M, δ ) .   (11)  Proof.   Regardless of the perturbation caused by the error, since   e k   ∈ B n p 2   (0 , γ k )   we have that   x k   ∈ { ˆ x k }⊕B n p 2   (0 , γ k ) . By enforcing the condition   ˆ x k   ∈ K k  ↓   ( M, δ, γ k )   we guarantee that   x k   ∈   ( K ↓ ( M, δ )   B n p 2   (0 , γ k ) )   ⊕ B n p 2   (0 , γ k )   ⊆ K ↓ ( M, δ )   since for any nonempty sets   X   and   Y   it holds true that   ( X   Y )   ⊕ Y ⊆ X   . For any given initial condition   x 0   we trivially have   e 0   = 0   and thus   K 0  ↓ ( M, δ, γ 0 ) =   K ↓ ( M, δ ) . We now describe a procedure to compute error bounds   γ k   which will be used for   a priori   construction of the sets   K k  ↓   ( M, δ, γ k ) . 7
Using the identity   e As   =   A ζ,s   +   E ζ,s   in (5) yields  x k   =   A ζ,δ   x k − 1   +  (∫   δ  0  A ζ,λ dλ   ·   B  )  u k − 1  ︸   ︷︷   ︸  ˆ x k  +   E ζ,δ   x k − 1   +  (∫   δ  0  E ζ,λ dλ   ·   B  )  u k − 1  ︸   ︷︷   ︸  e k  .   (12) Back-substituting the solutions   x k − 1   into   x k , and   ˆ x k − 1   into   ˆ x k   for every   k   = 1 , . . . , N δ   we can rewrite   e k   =   x k   −   ˆ x k  as  e k   =   ( ( A ζ,δ   +   E ζ,δ   ) k   −   A k ζ,δ  )   x 0   +  k − 1 ∑  i =0  [( ( A ζ,δ   +   E ζ,δ   ) i   −   A i ζ,δ  )  ×  ∫   δ  0  A ζ,λ dλ   + ( A ζ,δ   +   E ζ,δ   ) i  ∫   δ  0  E ζ,λ dλ  ]  Bu k − 1 − i .   (13) We can do so because for any discrete-time LTI system   x k +1   = Φ x k   + Ψ u k , the solution   x k   can be written in terms of the initial condition   x 0   and the past inputs as   x k   = Φ k x 0   +   ∑ k − 1  i =0   Φ i Ψ u k − 1 − i . To compute the upper-bound   γ k   on (13), invoke i) the inequality   ∥ ∥ A k ∥ ∥   ≤ ‖ A ‖ k   which holds for any matrix   A   and positive constant   k , ii) the multiplicative and triangular inequalities, and iii) the binomial expansion for   ( A ζ,δ   +   E ζ,δ   ) k  (which is valid since   A ζ,s E ζ,s   =   E ζ,s A ζ,s   [29]):  ‖ e k ‖ ≤  ∥ ∥ ∥ ∥ ∥  k ∑  l =0  ( k l  )  A l ζ,δ   E k − l ζ,δ   −   A k ζ,δ  ∥ ∥ ∥ ∥ ∥   ‖ x 0 ‖  +  k − 1 ∑  i =0  [∥ ∥ ∥ ∥ ∥  i ∑  l =0  ( i l  )  A l ζ,δ   E i − l ζ,δ   −   A i ζ,δ  ∥ ∥ ∥ ∥ ∥ ∥ ∥ ∥ ∥ ∥ ∫   δ  0  A ζ,λ dλ  ∥ ∥ ∥ ∥ ∥  +   ‖ A ζ,δ   +   E ζ,δ   ‖ i  ∥ ∥ ∥ ∥ ∥ ∫   δ  0  E ζ,λ dλ  ∥ ∥ ∥ ∥ ∥ ]  sup  u ∈U  ‖ Bu ‖   .   (14) For   k   = 1   this inequality is  ‖ e 1 ‖ ≤ ‖ E ζ,δ   ‖ ‖ x 0 ‖   +  ∥ ∥ ∥ ∥ ∥ ∫   δ  0  E ζ,λ dλ  ∥ ∥ ∥ ∥ ∥   sup  u ∈U  ‖ Bu ‖   .   (15) For   k >   1   we find  ‖ e k ‖ ≤  k − 1 ∑  l =0  ( k l  ) ∥ ∥ A l ζ,δ  ∥ ∥   ‖ E ζ,δ   ‖ k − l   ‖ x 0 ‖  +  k − 1 ∑  i =1  [ i − 1 ∑  l =0  ( i l  ) ∥ ∥ A l ζ,δ  ∥ ∥   ‖ E ζ,δ   ‖ i − l  ∥ ∥ ∥ ∥ ∥ ∫   δ  0  A ζ,λ dλ  ∥ ∥ ∥ ∥ ∥  + ( ‖ A ζ,δ   ‖   +   ‖ E ζ,δ   ‖ ) i  ∥ ∥ ∥ ∥ ∥ ∫   δ  0  E ζ,λ dλ  ∥ ∥ ∥ ∥ ∥ ]  sup  u ∈U  ‖ Bu ‖   .   (16) The term   ∫   δ  0   A ζ,λ dλ   in (16) is a definite integral of a finite sum and it can be easily computed:  ∫   δ  0  A ζ,λ dλ   =  ζ ∑  i =0  ∫   δ  0  ( Aλ ) i  i !   dλ   =  ζ ∑  i =0  A i δ i +1  ( i   + 1)!   .   (17) 8
Evaluating the integral   ∫   δ  0   E ζ,λ dλ   is trickier. We can, however, compute an upper-bound on its   ∞ -norm. We will use the following property [30]:  ‖ E ζ,δ   ‖ ∞   ≤   ( ‖ A ‖ ∞   δ ) ζ +1  ( ζ   + 1)!   ·   1 1   −   ε   =:   ψ δ   ,   (18) where the discretization order   ζ   is chosen such that   ε   :=   ‖ A ‖ ∞   δ ζ +2   <   1   (to ensure convergence of the power series  1 +   ε   +   ε 2   +   · · ·   ). Similarly, we can derive  ∥ ∥ ∥ ∥ ∥ ∫   δ  0  E ζ,λ dλ  ∥ ∥ ∥ ∥ ∥ ∞  ≤  ∫   δ  0  ‖ E ζ,λ ‖ ∞   dλ   ≤  ∫   δ  0  ∞ ∑  i = ζ +1  ‖ A ‖ i  ∞   λ i  i !   dλ  =  ∞ ∑  i = ζ +1  ∫   δ  0  ‖ A ‖ i  ∞   λ i  i !   dλ   =  ∞ ∑  i = ζ +1  ‖ A ‖ i  ∞   δ i +1  ( i   + 1)!  ≤ ‖ A ‖ ζ +1  ∞   δ ζ +2  ( ζ   + 2)!   ·   1 1   −   η   ≤   ψ δ   ·   δ ζ   + 2   ,   (19) where   η   =   ε   ·   ( 1   −   1  ζ +3  ) . Note that   ε <   1   ⇒   η <   1 , which automatically ensures the convergence of the power series  1 +   η   +   η 2   +   · · ·   used in derivation of (19). Substituting (17)–(19) into (15)–(16) for the   ∞ -norm we obtain a conservative bound    ̃ γ k   on   ‖ e k ‖ ∞   as   ̃ γ 1   :=   ψ δ   ‖ x 0 ‖ ∞   +  (  ψ δ   ·   δ ζ +2  )  sup  u ∈U  ‖ Bu ‖ ∞   ;   (20)   ̃ γ k   :=  k − 1 ∑  l =0  ( k l  ) ∥ ∥ A l ζ,δ  ∥ ∥ ∞   ψ k − l δ   ‖ x 0 ‖ ∞  +  k − 1 ∑  i =1  [ i − 1 ∑  l =0  ( i l  ) ∥ ∥ A l ζ,δ  ∥ ∥ ∞   ψ i − l δ  ∥ ∥ ∥ ∥  ζ ∑  j =0  A j   δ j +1  ( j   + 1)!  ∥ ∥ ∥ ∥ ∞  +   ( ‖ A ζ,δ   ‖ ∞   +   ψ δ  ) i   ψ δ   ·   δ ζ +2  ]  sup  u ∈U  ‖ Bu ‖ ∞   (21) for   k >   1 . Using the upper-bounds (20)–(21) in Lemma 2 allows us to check for feasibility of a given point   x 0   via only the finite-difference model, while ensuring that safety will not be violated due to discretization. The bound    ̃ γ k   is asymptotically tight in the sense that for any   k ,   lim ζ →∞    ̃ γ k   = 0 . In practice, the chosen order of discretization   ζ   must be large enough so as to ensure convergence of the power series in derivation of (18) as well as non-emptiness of the eroded sets in Lemma 2.  3.2.3   Verifying Feasibility of   x 0   via Forward Simulation  We can now simply use the discretized model  ˆ x k +1   =   A ζ,δ   ˆ x k   +   B ζ,δ   u k   (22) with   B ζ,δ   :=   ∫   δ  0   A ζ,λ dλB   to determine feasibility of a given initial condition   x 0   without worrying about the discretiza- tion error or the inter-sampling behavior of the continuous trajectories of (1) and their potentially negative impact on safety: If there exists a sequence of controls   { u k }   so that the nominal states   ˆ x k   of the closed-loop system belong to the precomputed sets   K k  ↓   ( M, δ,    ̃ γ k )   as described above, then via Lemmas 2 and 1 the trajectories of (1) never exit   K . Let us now construct the prediction equation as in (23) (where we have used the notation   Gx 0   +   H u   to abbreviate the right-hand side of the equality), and formulate the finite horizon feasibility program 9
        ˆ x 0  ˆ x 1  ˆ x 2  . . .  ˆ x N δ          =          I A ζ,δ  A 2  ζ,δ  . . .  A N δ  ζ,δ         ︸   ︷︷   ︸  G  x 0   +          0   0   . . .   0  B ζ,δ   0   . . .   0  A ζ,δ   B ζ,δ   B ζ,δ   . . .   0  . . .   . . .   . . .   . . .  A N δ   − 1  ζ,δ   B ζ,δ   A N δ   − 2  ζ,δ   B ζ,δ   . . .   B ζ,δ         ︸   ︷︷   ︸  H        u 0  u 1  . . .  u N δ   − 1       ︸   ︷︷   ︸  u  (23)  min .  u   0   (24a)  subj .   to   u   ∈ U N δ   (24b)  ˆ x k   ∈ K k  ↓   ( M, δ,    ̃ γ k ) ,   k   = 0 , . . . , N δ   (24c)  [ˆ x 0   · · ·   ˆ x N δ   ] >   =   Gx 0   +   H u .   (24d)  Theorem 1.   If   ∃ u ∗   satisfying   (24) , then   x 0   ∈   Viab sd  T   ( K ) . Proof.   The proof follows directly from Lemmas 1 and 2. More specifically, the prediction equation (24d), for a fixed input sequence   u , generates a forward simulation of the finite-difference model (22) over the desired horizon   [0 , N δ   ] ∩ Z  corresponding to the continuous time horizon   [0 , N δ   δ ] =   T . Constraint (24b) ensures that this input sequence is point- wise admissible (meaning that every member of the sequence belongs to   U ), while constraint (24c) restricts   ˆ x k   so that the trajectory of (1) evaluated at   t k   belongs to   K ↓ ( M, δ )   since via Lemma 2   ˆ x k   ∈ K k  ↓   ( M, δ,    ̃ γ k )   ⇒   x k   =   x ( t k )   ∈ K ↓ ( M, δ ) . Lemma 1 then automatically guarantees that   x ( t )   ∈ K ∀ t   ∈   T   which implies   x 0   ∈   Viab sd  T   ( K ) . The subroutine F EASIBLE   in Algorithm 2 employs Theorem 1 to determine the feasibility of a given sample point  x 0 .   Its computational complexity is proportional to the complexity of (24) which, with polytopic constraints for example, is simply a linear program (LP). 4  4   Conservatism and Convergence  4.1   Algorithm Correctness  Theorem 2.   Given convex sets   K   and   U   and an initial point   v 0   ∈   Viab sd  T   ( K ) ,   V N   =   P OLYTOPIC -A PPROX ( K , v 0 , N   )  is a subset of   Viab sd  T   ( K )   ∀ N   . Proof.   (By induction) First, it is obvious that for   N   = 0 ,   V 0   =   P OLYTOPIC -A PPROX ( K , v 0 ,   0) =   { v 0 } ⊆   Viab sd  T   ( K ) , since we are given that   v 0   ∈   Viab sd  T   ( K ) . Next, assume that   V N   − 1   =   P OLYTOPIC -A PPROX ( K , v 0 , N   −   1) = conv( { v 0 , . . . , v N   − 1 } )   ⊆   Viab sd  T   ( K ) , and that we have   v N   =   B ISECTION -F EASIBILITY ( v 0 , b,   K )   for some point   b   ∈   ∂ K . Let   V N   = conv( V N   − 1   ∪ { v N   } ) . Since B ISECTION -F EASIBILITY   only returns points which are inside   Viab sd  T   ( K ) , we know   v N   ∈   Viab sd  T   ( K ) . Now since   V N   is convex,   ∀ x 0   ∈ V N   ∃ x ′  0   ∈ V N   − 1   and   ∃ θ   ∈   [0 ,   1]   s.t.   x 0   =   θx ′  0   + (1   −   θ ) v N   .   For   x ′  0   we know (by induction hypothesis)   ∃ u x ′  0   ( · )   ∈   U   pwc  T   s.t.   x ′ ( t ) =   e At x ′  0   +   ∫   t  0   e A ( t − r ) Bu x ′  0   ( r ) dr   ∈ K ∀ t   ∈   T . For   v N   we also know  ∃ u v N   ( · )   ∈   U   pwc  T   s.t.   x ′′ ( t ) =   e At v N   +   ∫   t  0   e A ( t − r ) Bu v N   ( r ) dr   ∈ K   ∀ t   ∈   T . Therefore,   ̃ x ( t ) :=   θx ′ ( t ) + (1   −   θ ) x ′′ ( t ) =   e At   ( θx ′  0   + (1   −   θ ) v N   ) +  ∫   t  0  e A ( t − r ) B   ( θu x ′  0   ( r ) + (1   −   θ ) u v N   ( r ) )   dr  =   e At x 0   +  ∫   t  0  e A ( t − r ) Bu x 0   ( r ) dr   ∈ K ∀ t   ∈   T   (25)  4 We note that the vast majority of the calculations for    ̃ γ k   from (20)–(21) can be done once and ahead of time. Online, to be able to construct  K k  ↓   in (24c), two simple operations (a multiplication by   ‖ x 0 ‖   and an addition) are all that is needed to form    ̃ γ k   .  10
since   K   and   U   are convex and compact. Thus,   u x 0   ( · ) =   θu x ′  0   ( · ) + (1   −   θ ) u v N   ( · )   ∈   U   pwc  T   is safety-preserving and  x 0   ∈   Viab sd  T   ( K ) . Because   x 0   was chosen arbitrarily in   V N   , we conclude that   V N   ⊆   Viab sd  T   ( K ) .  4.2   Algorithm Convergence  One of the striking features of our algorithm is that it is random; additional points on the boundary of the polytopic approximation of   Viab sd  T   ( K )   are iteratively generated based on a random sampling.   This random nature is due to the fact that   ∂   Viab sd  T   ( K )   is unknown   a priori   (and is, in fact, what we are trying to estimate) so it is impossible to know what points to sample to construct a polytope that converges to   Viab sd  T   ( K )   as quickly as possible. In fact, any algorithm which deterministically chose points for which to verify feasibility could be presented with a safe set  K   and system dynamics for which the algorithm would converge arbitrarily poorly. This fact is related to results in the literature of estimating the volume of convex bodies using a   separation oracle 5 . More specifically, part of the literature on algorithms for estimating the volume of convex bodies states that it can be shown that for any algorithm that deterministically queries a separation oracle a polynomial number of times to build a polytopic approximation, the error (the difference in volume between the approximation and the true set) could be exponential in the number of dimensions [31]. Random algorithms, on the other hand, can perform in a provably better manner [32]. To prove our algorithm’s asymptotic convergence, first let   Viab sd  T   ( K )   be the subset of the viability kernel we are actually attempting to approximate, i.e.  Viab sd  T   ( K ) :=   { x 0   | ∃ u ( · )   ∈   U   pwc  T   ,   ∀ k, x ( t k )   ∈ K ↓ ( M, δ ) } .   (26) While the true kernel contains all initial conditions for which a piecewise constant control keeps   x ( t )   ∈ K , the above set only encompasses initial conditions for which   x ( t k )   ∈ K ↓ ( M, δ )   (which is a sufficient, but not necessary, condition to imply   x ( t )   ∈ K ; cf. Lemma 1). Define the volumetric error between these two sets as   cont ( M δ ) := vol(Viab sd  T   ( K ))   −   vol(Viab sd  T   ( K ))   (27) and note that it depends only on the term   M δ , due to the definition of   K ↓ ( M, δ )   in (4). Next, consider the output   V N   of Algorithm 1. Clearly, the accuracy of under-approximation of the viability kernel by the set   V N   is implicitly dependant on the discretization order   ζ   (Section 3.2.2), and the accuracy      of the bisection search (Algorithm 2). To reflect this dependency, we adapt the extended notation   V ζ, N   . 6   Evidently, for fixed values of  ζ   and      as   N   → ∞ , this set only approximates a subset   Viab sd  T   ( K , ζ,  )   of the set   Viab sd  T   ( K ) :  Lim sup  N   →∞  V ζ, N   = Viab sd  T   ( K , ζ,  )   (28)  Lim sup  ζ →∞ , → 0  Viab sd  T   ( K , ζ,  ) = Viab sd  T   ( K )   (29) with   Lim sup   denoting the Kuratowski upper-limit. We are now ready to present our algorithm’s convergence property.  Proposition 1   (Rate of Convergence) .   Let    vol ( N, ζ,  )   be the volumetric error between the viability kernel and the output of our algorithm, minus the error    cont ( M δ )   between the true kernel and   Viab sd  T   ( K ) :   vol ( N, ζ,  ) := vol(Viab sd  T   ( K ))   −   vol( V ζ, N   )   −    cont ( M δ ) .  Then our algorithm converges as  lim  N   →∞  ζ →∞   → 0   vol ( N, ζ,  ) N   2  n − 1   =   c n (Viab sd  T   ( K ) , M δ ) ,   (30)  where   c n (Viab sd  T   ( K ) , M δ )   is a constant dependent on the dimension   n , the shape of the viability kernel (more specifi- cally its Gauss-Kronecker curvature), and the value   M δ .  5 A separation oracle is a function that accepts as input a convex set and a point, and returns whether or not that point is inside the convex set. In our algorithm F EASIBLE ( x,   K )   plays this role.  6 Theorem 2, restated in terms of the extended notation, asserts that   V ζ, N   ⊆   Viab sd  T   ( K )   ∀ ζ, , N   .  11
0   10   20   30   40  0  50  100  150  Dimension (n) Run time [s] Figure 5: Run time of the algorithm for a chain of   n   integrators. The proof requires some background on random algorithms for convex bodies and is provided in Appendix A. Proposition 1 asserts that, for fixed dimension   n , the volumetric error between the outcome of our algorithm and the true viability kernel asymptotically converges, at the exponential rate of   ˆ c n (Viab sd  T   ( K )) /N   2  n − 1   , to a numerical constant due to the sampled-data nature of the system. On the other hand, to keep the accuracy of the approximation the same as   n   →   kn   we would need an increase of   N   →   N   kn − 1  n − 1   . However, the fact that we only store samples on the boundary of the viability kernel to describe that set (as opposed to storing a grid of the entire set   K   and possibly beyond) requires significantly less memory than conventional approaches such as the SD level-set method in [7]. The flexibility in choosing the number of samples strikes a direct tradeoff between accuracy and computational complexity, making our algorithm scalable to high dimensions. The computed approximation is far more accurate (and quite possibly more scalable) than the piecewise ellipsoidal technique also presented in [7]. Again, due to the results in [33], the above convergence rate is optimal (up to a multiplicative constant depending on the probability density function used for sampling); no other algorithm that approximates the kernel by sampling from its boundary will be able to converge at a faster rate.  5   Computational Complexity & Scalability  The run time complexity of our algorithm (for fixed number of sampling intervals   N δ   ) is   O ( N   log( d )Φ( n )) , where  N   is the number of samples/vertices,   d   is the “diameter” of the set   K , and   Φ   is the complexity of the feasibility program (24) as a function of the state dimension   n .   That is, the algorithm runs in time linear in the number of samples   N   , logarithmic in the diameter   d   of   K   due to complexity of the bisection search, and proportional to   Φ   in the complexity of the appropriate feasibility program (24). For instance, with polytopic constraints, the feasibility problem is an LP and thus the algorithm runs in time sub-cubic in   n . This is a direct improvement over existing techniques for approximating the SD viability kernel. Furthermore, since each vertex is processed completely independently of others, our algorithms is highly parallelizable. To demonstrate the scalability of our algorithm, consider the chain of   n   integrators   d n x/dt n   =   u   with constraints  U   = [ − 0 . 15 ,   0 . 15]   and   K   =   { x   | ‖ x ‖ ∞   ≤   0 . 5 } . The state is measured every   δ   = 0 . 05   s and safety is to be maintained over   T   = [0 ,   1] . We use a discretization order of   ζ   = 4 , bisection accuracy of      = 0 . 01   with maximum of three-level bisection depth, and employ YALMIP [34] to implement (24) and MPT [35] for simple operations with polytopes. All of these parameters are kept constant as we increase the dimension   n   and the number of samples   N   = 2 n   to examine the scalability of our algorithm. The results are shown in Fig. 5. The algorithm is implemented in MATLAB R2011b and tested on an Intel Core i7 at   2 . 9   GHz with   16   GB RAM running   64 -bit Windows 7 Pro (without optimizing the code for speed). 12
K  v i  ∂   Viab sd  T   ( K )  r >  d i   x   =   ρ K ( r d i   )  v 0  r d i  (unknown)  y ij   c  Figure 6: An    o -accurate bisection search between   v i   and   c   determines via (33) the supporting hyperplane of the viability kernel.  6   Bounding the Error for Finite Number of Samples: Computing a Tight Over-Approximation  Every convex set can be over-approximated by any finite collection of its   support functions . The support function of a convex compact set   C ⊂   R n   along   `   ∈   R n   is  ρ C   ( ` ) := max  x ∈C   ` > x.   (31) The half-space   { x   |   ` > x   ≤   ρ C   ( ` ) }   contains   C , and the hyper-plane   { x   |   ` > x   =   ρ C   ( ` ) }   is a supporting hyperplane for  C   with normal vector   `   and distance value   ρ C   ( ` ) . It follows that   C ⊆   ⋂  ` ∈L { x   |   ` > x   ≤   ρ C   ( ` ) }   with   L   a finite subset of   R n . Let   r d i   be the direction   r d   along which we have determined the vertex   v i   of the under-approximation set through the   i th iteration of Algorithm 1. It is easy to compute the support function of the set   K   along this direction:   ρ K ( r d i   ) = max x ∈K   r >  d i   x . To find the supporting hyperplane of the true, unknown viability kernel (or rather some appropriate approximation of it) in the direction   r d i   , we move the hyperplane  { x   |   r >  d i   x   =   ρ K ( r d i   ) }   (32) on the interior of   K   until we find at least one feasible point on this plane that belongs to the kernel (or its over- approximation); see Fig. 6. We do so iteratively by first performing an    o -accurate bisection search between the point   v i   and the point of intersection   c   of the ray   ~ r i   passing through   v i   and   v 0   with the hyperplane (32). This gives us points   y ij   indexed by each step   j   of the new bisection search. We then solve the following modified convex feasibility program   for every   y ij  for the given direction   r d i   :  min .  u ,x 0  0   (33a)  subj .   to   u   ∈ U N δ   (33b)  Gx 0   +   H u   ∈ K N δ   +1   (33c)  r >  d i   x 0   =   r >  d i   y ij   .   (33d) Notice that we no longer fix   x 0 ; rather, we implicitly look for a point on the set   { x 0   |   r >  d i   x 0   =   r >  d i   y ij   } ∩ K , when  y ij   varies, that is a feasible point. We also do not erode the constraints as we did before since our goal here is find an  over-approximation   of the true kernel. 13
̂ V N  V N  Figure 7: The viability kernel is sandwiched in between the under-approximation   V ζ, N   and the over-approximation  ̂  V ζ, o  N   . Once a feasible solution to (33) is found for the desired accuracy    o   of the bisection search, we stop the iterations and store two entities: (S1) The last   infeasible   step, i.e. the last value of   y ij   for which (33) was infeasible. Denote this value by   y inf ∗  ij   ; (S2) The feasible solution pair   ( u ∗  i   , x ∗  0 i )   (indexed by   i   to correspond to the direction   r d i   ). We first use entity (S1) to form our over-approximation along   r d i   as the halfspace  { x   |   r >  d i   x   ≤   r >  d i   y inf ∗  ij   } .   (34) Clearly the set   { x   |   r >  d i   x   =   r >  d i   y inf ∗  ij   }   is a supporting hyperplane of   Viab sd  T   ( K )   with an arbitrary (and desirably conservative) error    o . By repeating the above procedure for all   N   directions and forming the set  ̂  V ζ, o  N   :=  N ⋂  i =1  { x   |   r >  d i   x   ≤   r >  d i   y inf ∗  ij   }   (35) we obtain (Fig. 7)  V ζ, N   ⊆   Viab sd  T   ( K )   ⊆   ̂   V ζ, o  N   .   (36) The error   vol(Viab sd  T   ( K ))   −   vol( V ζ, N   )   of our main under-approximation algorithm can then be quantitatively bounded above as   vol( ̂   V ζ, o  N   )   −   vol( V ζ, N   ) .   This upper-bound monotonically decreases as   N   increases, and converges to a numerical constant as   ,  o   →   0   and   N, ζ   → ∞ . The entity (S2) from the feasibility program (33) can help us improve our under-approximation. This is discussed next.  7   Improving the Under-Approximation  We describe two techniques that help improve the quality of our under-approximation for finite number of samples.  7.1   Center of Mass & Gauss-Kronecker Curvature  If we were to naively sample from a uniform distribution on   S n − 1 2   to compute our under-approximating set, two elements would negatively affect the quality of such approximation: (a) the unbalanced distances between the unit ball centered at   v 0   and the boundary points of the true kernel, i.e. the distance between   v 0   and the center of mass of the kernel, and (b) the regions of the boundary of the kernel with high curvature. This is because a uniform distribution 14
Figure 8: Mapping a uniform distribution on   S n − 1 2   ⊆ C   linearly onto   ∂ C   may result in a distribution that is non- uniform.  ̂  V N  V N Figure 9: Injecting additional under-approximation steps guided by the available over-approximation facets improves the quality of the resulting set (compare to Fig. 7). Areas of the kernel with high curvature are now covered at a faster rate. on   S n − 1 2   mapped onto   ∂   Viab sd  T   ( K )   could yield a distribution that is far from uniform depending on severity of (a) or (b); Fig. 8. The former issue can be somewhat mitigated by continually moving the point   v 0   to the center (e.g. in the sense of Chebyshev) of the set   V ζ, N   every time a new vertex is added. The latter issue can be addressed using the information obtained from the over-approximation procedure discussed in the previous section, specifically using the stored en- tity (S2). The idea is that the greater the distance between the   i th under-approximation vertex and over-approximation facet, the higher the Gauss curvature of the boundary of the true kernel in the neighborhood of that unexplored region. To account for this problem, we perform an additional step after the   i th iteration of our combined algorithm, every time a vertex is added to the under-approximation and an over-approximating halfspace is formed along the direction   r d i   : The stored value   x ∗  0 i   in (S2) approximates the   support vector   (the point in which the support function of a convex set touches its boundary) of the viability kernel in the direction   r d i   with    o   accuracy. Therefore, we use   x ∗  0 i  and execute a single instance of our under-approximation procedure this time along not   r d i   , but along the direction  x ∗  0 i   −   v 0   (corresponding to the ray passing through   v 0   and   x ∗  0 i ). By doing so, we allow the over-approximation to “guide” where we should look for the next under-approximation vertex. Moreover, for this new instance of the under-approximation algorithm we can limit the bisection search to a diameter of roughly   2  o   around   x ∗  0 i   (instead of perfoming the search between   v 0   and the point at which the new ray intersects   ∂ K ) since we know that   x ∗  0 i   is already fairly close to the boundary of the true kernel. The resulting vertex is not only in close proximity to the point where the over-approximating halfplane at the   i th iteration has been formed (thus providing a superior confidence that the viability kernel is sandwiched tightly in that area), but also covers the parts of the boundary that could potentially have a high Gauss curvature; Fig. 9. 15
κ 1   κ 2  κ 3  Figure 10: The vMF density function on the unit ball for three different mean directions and concentrations with  κ 1   > κ 2   > κ 3 .  7.2   Biased Random Sampling  The shortcomings of uniform sampling are more pronounced in high dimensions. Thus we additionally seek to   bias  the distribution on the unit ball to mitigate these shortcomings. To this end, we present a few heuristic techniques that are still based on random sampling so as to keep the optimality results of Section 4.2, but could potentially improve the performance of the algorithm. We will make use of the von-Mises Fisher (vMF) distribution [36] whose density function is given by  f vMF ( x ;   μ, κ ) :=   C ( κ ) e κμ >   x   (37) with   concentration   κ   ≥   0 , mean direction   μ   ( ‖ μ ‖   = 1 ), and a normalizing constant   C ( κ ) . The parameter   κ   determines how samples drawn from this distribution are concentrated around the mean direction   μ . For   κ   = 0   the vMF reduces to the uniform density; otherwise, it resembles a normal density (with compact support on the unit ball), centered at   μ  with variance inversely proportional to   κ   (Fig. 10). As   κ   → ∞   the vMF converges to a point distribution. The solutions we propose are by no means exhaustive, and there may be better ways of guiding the sampling process depending on the problem in hand or if we have some   a priori   knowledge of the shape of the viability kernel.  7.2.1   Gradient-Like Methods  The first approach we discuss is related to how some measure of the error varies in consecutive steps of our combined under- and over-approximation algorithm.   We let this change dictate from what vMF distribution should the next sample be drawn. Let   V ζ, i   and   ̂   V ζ, o  i   respectively denote the under-approximation and over-approximation of the viability kernel after adding the   i th sample direction   r d i   , and let the state space be equipped with some metric   d . Define the error function (possibly nonconvex)   err :   R n   →   R ≥ 0   at this   i th iteration as  err( r d i   ) :=   d ( V ζ, i   ,   ̂   V ζ, o  i   ) .   (38) Ideally, the metric   d   is chosen such that   err( r d i   )   is monotonically non-increasing as   i   increases (i.e. as we produce more accurate approximations). Define the normalized quantity  ∇   err( r d i   ) := 1   −   d ( V ζ, i   ,   ̂   V ζ, o  i   ) /d ( V ζ, i − 1 ,   ̂   V ζ, o  i − 1   )  1  π   cos − 1   ( 〈 r d i   , r d i − 1   〉 )   ,   (39) treating the pathological case   0 0   as   0 . The magnitude of   ∇   err( r d i   )   approximates the rate of improvement between two consecutive iterations of our combined algorithm. Thus, we can use this information to draw the next sampling direction   r d i +1   as  r d i +1   ∼   f vMF ( x ;   μ i +1 , κ i +1 )   (40) with mean and concentration parameters updated as  κ i +1   = max   { ν 0 , ω i }   ,   (41)  μ i +1   = sgn( − ν 0   +   ω i )   ·   r d i   ,   (42) 16
where  ω i   :=   ν 1   tanh ( ν 2 ∇   err( r d i   ))   .   (43) Here the positive scalars   ν 0 ,   ν 1 , and   ν 2   are design choices; for example,   ν 2   dictates how quickly the tangent hyperbolic function reaches its maximum   ν 1   where it levels off, while   ν 0   (which is a small scalar) limits how close the resulting vMF distribution can be to the uniform one and, when necessary, sets the sign of the direction   μ i +1   of the distribution to the negative of   r d i   . The intuition behind such an update rule is that when sampling along a given direction   r d i   causes the magnitude of   ∇   err   to become large (thus causing the error to decrease quickly), drawing the next sample from a distribution that is concentrated (proportional to this change) around the same direction may be a good choice to continue to reduce the error:   κ i +1   =   ω i   and   μ i +1   =   r d i   . If, on the other hand, the magnitude of   ∇   err   is small, we are approaching a local minima of the quasilinear function   err( · ) .   Therefore, we would want to sample from a distribution that is  not   concentrated around the current direction   r d i   , and could also be such that it potentially minimizes the likelihood of sampling around   r d i   . In such a case, we would bias the sampling to be slightly concentrated around the opposite direction (and also decrease the likelihood of resampling around   r d i   ) so as to attempt to get us out of this local minima, while not making it to be too far off from a uniform distribution:   κ i +1   =   ν 0   and   μ i +1   =   − r d i   .  Volumetric Measure:   We can use volume as a metric in (38) so that   err( r d i   ) := vol( ̂ V ζ, o  i   )   −   vol( V ζ, i   ) . Of course, computing the volume of a polytope is   # P-hard [37]. However, we can conservatively approximate   err( r d i   )  by calculating analytically the volume of appropriate inscribed and circumscribed ellipsoids: A minimum volume circumscribed ellipsoid (mVCE) containing a polytope that is represented by its vertices can be computed efficiently via a semidefinite program [38]. The same is true for a maximum volume inscribe ellipsoid (MVIE) that is contained in a polytope represented by its facets. (Note that the converse problems are NP-hard.) If we computed the mVCE of  V ζ, i   and shrunk it by a factor of   n   (the dimension), then the resulting ellipsoid would be a subset of   V ζ, i   . Similarly, if we enlarged the MVIE of   ̂   V ζ, o  i   by a factor of   n , then   ̂   V ζ, o  i   would be a subset of the resulting ellipsoid. The difference in the volume of these two ellipsoids would be an upper-bound on   err( r d i   ) . Working with the volume of these extremal ellipsoids does preserve order, in that the relations   vol( V ζ, i +1 )   ≥  vol( V ζ, i   )   and   vol( ̂   V ζ, o  i +1   )   ≤   vol( ̂   V ζ, o  i   )   are also true for their extremal ellipsoids. The reason is that adding constraints (additional vertices or facets) to a convex optimization problem (the SDPs) cannot decrease the value of the optimal objective functions (correlated with ellipsoid volumes). As a result, the error upper-bound, just like the error itself, is monotonically non-increasing as   i   increases. A disadvantage of working with scaled extremal ellipsoids is that the shrinkage/enlargement by a factor of   n   can be too conservative particularly for larger   n .  Hausdorff Distance:   The Hausdorff distance between two compact convex sets   C 1   and   C 2   in   R n   in terms of their support functions is defined as  d H ( C 1 ,   C 2 ) :=   max  ` ∈S n − 1   {| ρ C 1   ( ` )   −   ρ C 2   ( ` ) |}   .   (44) We can use this metric to define the error:   err( r d i   ) =   d H ( V ζ, i   ,   ̂   V ζ, o  i   ) .   Computing this distance between   V ζ, i   = conv( { v j   } i j =0 )   and   ̂   V ζ, o  i   =:   P i x   ≤   p i   (where   P   and   p   respectively are the appropriate matrix and vector corre- sponding to the facet-based outer polytope defined in (35) but with   i   faces) can be cast as a series of LPs: For a fixed direction   `   ∈ S n − 1 , we compute at every iteration   i ρ V ζ, i  ( ` ) =   max  x,   { λ j   ≥ 0 }   ` > x   (45a)  subj .   to   x   =  i ∑  j =0  λ j   v j   ,  i ∑  j =0  λ j   = 1 ,   (45b) as well as  ρ ̂   V ζ, o  i  ( ` ) = max  x   ` > x   (46a)  subj .   to   P i x   ≤   p i .   (46b) 17
The Hausdorff distance can then be approximated over a finite number of directions   L ⊂ S n − 1   as  ˆ d H ( V ζ, i   ,   ̂   V ζ, o  i   ) := max  ` ∈L  {∣ ∣ ∣ ρ ̂   V ζ, o  i  ( ` )   −   ρ V ζ, i  ( ` )  ∣ ∣ ∣ }  ≤   max  ` ∈S n − 1  {∣ ∣ ∣ ρ ̂ V ζ, o  i  ( ` )   −   ρ V ζ, i  ( ` )  ∣ ∣ ∣ }  =   d H ( V ζ, i   ,   ̂   V ζ, o  i   ) .   (47) This approximation is not conservative due to the inequality in (47). But as far as guiding the sampling, it may yield a viable alternative to the volumetric error described above (that could be excessively conservative). Unfortunately, performing these additional LPs at the end of each iteration of our algorithm may undermine its efficiency.  7.2.2   Purely Heuristic Methods  We also propose two purely heuristic methods as alternatives to the gradient-like methods presented above.  Averaged Opposite Direction:   For this approach we associate a given vMF density function to the negative of  each   individual direction vector we have generated so far, and draw our next sampling direction from a convolution of these density functions. That is, in some sense we are drawing at random a sample whose expected value lies in the opposite direction of the samples we have already drawn. The intuition here is that we can approach “true” uniformity at a higher rate by sampling in the direction whose neighborhood we have not yet sampled as densely. We can vary the concentration of the vMF densities as an increasing function of the iteration step   i   (e.g. linearly with the number of vertices), such that at the beginning steps of our algorithm these densities are closer to uniform, and as we progress and generate more vertices/directions their concentrations increase. 7   A design parameter   ν 1   can be a multiplier for the number of vertices. We cap the concentration of the distributions to 100 so as to not lose the randomness properties of the algorithm when   i   grows too large.  Point-to-Plane Distance:   This approach is based on simply calculating the smallest distance from every vertex of  V ζ, i   to facets of   ̂   V ζ, o  i   , and identifying the largest of such distances. This quantity in some sense describes the gaps between the two sets over which we are most uncertain about the boundary of the true kernel. Let    ̄ v j   be the vertex in   V ζ, i   = conv( { v j   } i j =0 )   with the largest point-to-plane distance. To guide the sampling at the next iteration of the algorithm, we calculate the ray that passes through the center   v 0   and the point in the facet of  ̂  V ζ, o  i   that is the closest to    ̄ v j   . To determine our next sampling direction   r d i +1   , we then sample from a vMF distribution whose mean is the unit vector along this ray and whose concentration is dependant on the calculated point-to-plane distance from    ̄ v j   to   ̂   V ζ, o  i   . A design parameter can again be a multiplier   ν 1   for this distance.  7.2.3   Calibration and Comparison on Random Systems  To compare the performance of the four guiding methods described above, we first roughly calibrated the design parameters for each approach by running each guiding technique on 30 randomly generated systems across 2D, 3D, and 4D state and input dimensions 8   for a variety of exponents of 10. For example, we tested the volume extremal ellipsoids and the Hausdorff distance methods for all combinations of   ν 0   = 10 {− 1 , 0 } ,   ν 1   = 10 { 0 , 1 , 2 } , and   ν 2   = 10 { 0 , 1 } /n   (made dimension   n   dependent so that the tangent hyperbolic reaches its maximum   ν 1   more slowly when   n   is larger); and the averaged opposite direction and the point-to-plane distance methods for   ν 1   = 10 {− 2 , − 1 , 0 , 1 } /n   (made dimension dependent so that the vMF concentration increases more slowly when   n   is larger). We found the optimal parameters to be   ν 0   = 0 . 1 ,   ν 1   = 1 ,   ν 2   = 1 /n   for the volume extremal ellipsoids;   ν 0   = 1 ,   ν 1   = 1 ,   ν 2   = 1 /n   for the Hausdorff distance;   ν 1   = 0 . 1 /n   for the averaged opposite direction; and   ν 1   = 1 /n   for the point-to-plane distance. We emphasize that this calibration was meant only to achieve parameters that were roughly the correct order of magnitude; we did no further fine tuning. With the calibrated methods in hand, we then ran all four methods on a fresh set of 60 randomly generated systems across 3D, 4D, and 5D state and input dimensions. We were then able to evaluate the performance of each method by taking advantage of the fact that all the approximations generated are under-approximations.   Thus, the under- approximation with the largest possible volume must be (in a volumetric sense) the closest to the true viability kernel.  7 The density function of the convolution is closer to uniform if (a) the individual densities are closer to uniform, or (b) the existing directions all cancel each other out.  8 Our tests were limited by 5D due to the need to directly compute the volume of polytopes for performance assessment.  18
10   20   30   40   50  −18  −16  −14  −12  −10  −8  −6  −4  −2  0  Number of Samples Percent Volume Difference  Uniform Only  Averaged Opposite Direction  Volume Extremal Ellipsoids  Point to Plane  Hausdorff Figure 11: Percent difference in volume between the result of each guiding method and the best result out of all the guiding methods for a given number of vertices (averaged over 60 randomly generated systems across 3, 4, and 5 state and input dimensions). Using this fact we found for each random system and for a given number of vertices the best under-approximation selected from the results of all four different methods (as well as the basic uniform sampling method). Then, for each system and number of vertices we calculated the percent difference in volume between the best under-approximation and the result of each method. (Note that this percent difference must be negative, since the result from each method must have a volume smaller than or equal to the result of the best method.) Thus a method with a less negative value, i.e. closer to 0, indicates a method that performs better, in a volumetric sense, than a method with a more negative value (which indicates a much smaller and thus less accurate under-approximation in a volumetric sense). Fig. 11 plots the resulting percent volume difference by method as a function of the number of vertices, averaged over the 60 randomly generated systems. On average, the averaged opposite direction (blue) outperforms all methods. More specifically, it improves the resulting under-approximation volume by about 5% over the uniform sampling (red). The Hausdorff distance (yellow) performs consistently better than uniform sampling but since we only compute a crude approximation of the actual Hausdorff distance in (47) (at every step we only use the directions along which we have already generated a vertex so as to maintain monotonicity of the distance approximate while keeping the computation times on par with the other methods), this improvement is expected to be more emphasized for higher number of samples than only 50.   The volume extremal ellipsoids (green) initially performs far better than uniform sampling in 3D and 4D, but then quickly degrades in 5D due to shrinkage/expansion of the ellipsoids by a factor of   n   that is much more exaggerated in higher dimensions. The point-to-plane distance (cyan) simply does not show any improvement over uniform sampling, at least on average in our test setup. We will use the averaged opposite direction as our guiding method of choice for our examples in Section 8. However, we do emphasize again that the performance of these heuristics is problem- dependent. Therefore, while in our limited random tests one method might have outperformed others on average, it may be that another method is even more suitable and performs significantly better for a given problem.  8   Examples  8.1   The Double-integrator  First, consider the simple dynamics    ̈ x   =   u   with   δ   = 0 . 05   s.   The constraints   K   =   { x   | ‖ x ‖ ∞   ≤   0 . 5 }   and   U   = [ − 0 . 15 ,   0 . 15]   are to be respected over   T   = [0 ,   1] . To find a bound   M   on the vector field and construct the eroded set   K ↓ ( M, δ )   we use the following procedure: We first scale the state space by computing for every dimension   d   the quantity   z ∗ ( e d )   −   z ∗ ( − e d )   with   z ∗ ( e d ) := arg max z   e >  d   z   subject to   z   =  ̇ x   and   x   ∈ K ,   u   ∈ U , and where   e d   denotes 19
x   ̇ x  −0.4   −0.2   0   0.2   0.4   0.6  −0.4  −0.2  0  0.2  0.4  0.6 Figure 12: Polytopic under-approximation   V ζ, N   (blue) and over-approximation   ̂   V ζ, o  N/ 2   (lavender) of   Viab sd  T   ( K )   for the double-integrator example with   N   = 20   samples.   The sets   K   and   K ↓ ( M, δ )   are shown in dark and light green, respectively. The SD level-set approximation [7] is also shown (outlined in thick black line). the standard basis vector spanning   d th dimension. Then we divide all such quantities by their minimum value among all dimensions, and transform the dynamics accordingly so that a state increment in all dimensions is equivalent. At this stage, we can calculate   M   as the optimal value of   max x,u   ‖    ̇ x ‖ ∞   subject to   x   ∈ K   and   u   ∈ U . The eroded set  K ↓ ( M, δ )   is now constructed in the scaled state space according to Lemma 1. We used   ζ   = 4 th order discretization and      =    o   = 0 . 01 -accurate bisection search to obtain the under- and over- approximations shown in Fig. 12. The under-approximation is computed using   N   = 20   randomly generated samples via Algorithm 1. The sampling process was guided via the techniques in Section 7. The over-approximation polytope consists of   N/ 2   facets and is computed according to Section 6. The overall computation time was   18   s.  8.2   12D Quadrotor Flight Envelope Protection  We now evaluate our algorithm on the benchmark example described in [39].   Consider the full-order model of a quadrotor based on the nonlinear Newton-Euler rigid body equations of motion. The state vector  x   =   [ x   y   z    ̇ x    ̇ y    ̇ z   φ   θ   ψ    ̇ φ    ̇ θ    ̇ ψ ] >   ∈   R 12   (48) is comprised of translational positions in   [ m ]   with respect to a global origin, their derivatives (linear velocities in   x ,   y ,  z   directions) in   [ m / s ] , the Eulerian angles roll   φ , pitch   θ , and yaw   ψ   in   [ rad ] , and their respective derivatives (angular velocities) in   [ rad / s ] . The control input is the vector   u   = [ u 1   u 2   u 3   u 4 ] >   ∈   R 4   consisting, respectively, of the total thrust in   [ m / s 2 ]   normalized with respect to the mass of the quadrotor ( u 1 ) and the second-order derivatives    ̈ φ ,    ̈ θ ,    ̈ ψ   of the Eulerian angles in   [ rad / s 2 ]   ( u 2   through   u 4 ). The system is under-actuated since there are six degrees of freedom but only four actuators. By linearizing the equations of motion about the hover condition   φ   = 0 ,   θ   = 0 , and   u 1   =   g  (with   g   ≈   9 . 81   being the acceleration of gravity) one would obtain the model    ̇ x   =   Ax   +   Bu   with the state and the input now representing   deviation   from the equilibrium. We follow the example detailed by [40] of an agile quadrotor in which the state is sampled at a frequency of   10   Hz. (See the same reference for values of the system matrices   A   and  B .) For safe operation of the vehicle the Eulerian angles   φ   and   θ   and the speed profile   V   :=   ‖ [  ̇ x  ̇ y  ̇ z] ‖   are bounded as  φ, θ   ∈   [ −   π  4   ,   π  4   ]   and   V   ≤   5 . The angular velocities are constrained as    ̇ φ,    ̇ θ,    ̇ ψ   ∈   [ − 3 ,   3] . We further assume that the vehicle must safely fly within the range of   1   to   7   m above the ground in z direction in an environment that stretches  6   m in each direction in the   x - y   plane. These constraints form the flight envelope   K . The vector   u   is constrained by the hyper-rectangle   U   := [ − g,   2 . 38]   ×   [ − 0 . 5 ,   0 . 5] 3 . 20
x   ̇ x  −5   0   5  −5  0  5 y   ̇ y  −5   0   5  −5  0  5 z   ̇ z  0   2   4   6   8  −4  −2  0  2  4 φ   ̇ φ  −2   0   2  −1  0  1 θ   ̇ θ  −2   −1   0   1   2  −2  −1  0  1 ψ   ̇ ψ  −2   0   2  −2  0  2  ̇ x   ̇ y  −5   0   5  −5  0  5  ̇ φ   ̇ θ  −2   0   2  −2  −1  0  1 Figure 13: Selected 2D projections of the polytopic approximation of   Viab sd  T   ( K )   for the flight envelope example. Under-approximations for   N   = 24 ,   48 ,   96 ,   500 ,   1000 ,   2000 ,   3000   vertices are shown with   N   = 24   in the lightest shade of blue (innermost set), and   N   = 3000   in the darkest shade of blue. An over-approximation (outermost set) with   N/ 2   facets is also shown in lavender. The approximations are tight and touch each other to within a constant accuracy in at least   N/ 2   points—a fact that is obscured by the projections. The quadrotor can travel a distance of roughly half a meter between any two consecutive sampling times despite the relatively high sampling frequency. This fact further warrants the treatment of such a safety-critical system through a SD framework. We wish to compute the set of initial states for which safety can be maintained over   T   = [0 ,   2] . We warm-start our approximation algorithm (in the scaled state space) by first sampling along all axes in order to obtain a full-dimensional object, and then along 72 uniformly-spaced vectors in   x i - x i +3   and   x i +6 - x i +9   subspaces for  i   = 1 ,   2 ,   3 . The remaining samples are generated randomly by guiding the vMF sampler via the averaged opposite direction method. The discretization order and the bisection search accuracies are the same as in the previous example. Fig. 13 shows selected 2D projections of our sampling-based polytopic approximations of   Viab sd  T   ( K )   in the orig- inal unscaled state space. The algorithm requires about   3   s (without optimizing the code for speed) to generate a new vertex of the under-approximation and   5   s to generate a facet of the over-approximation. The under-approximation is tight in the sense that each vertex of the polytope belongs to the boundary of the true viability kernel with some   a priori   known accuracy due to the SD nature of the problem. The over-approximation is also tight in that each facet of the polytope touches the boundary of the true kernel in at least one point. The two approximating sets sandwich the boundary of the viability kernel to within a certain precision in at least half of the number of vertices of the under- approximation, providing an added layer of confidence about the precise location of the kernel. The tightness of the sets are unfortunately unobservable in the projection plots.  9   Conclusions, Extensions, and Future Work  We presented a scalable sampling-based algorithm to generate tight under- and over-approximations of the viability kernel under SD LTI dynamics. We provided correctness and convergence results, discussed a number of heuristics to bias the sampling process for improved performance, and demonstrated the algorithm on a 12D problem of flight envelope protection for an autonomous quadrotor. The extension to discrete-time systems is straightforward: The set   K   is directly used in the feasibility program (24) and thus the vertices of the resulting under-approximation and the facets of the over-approximation touch the boundary of the true viability kernel exactly. The algorithm can also be extended to piecewise affine SD (or discrete-time) system 21
by recasting the optimization problem as a mixed-integer program. Due to space limitations we will present our results on the synthesis of the safety-preserving controllers as well as sufficient conditions under which a generated under- approximation is controlled-invariant (and thus could be used to enforce infinite-horizon safety) in a separate paper. Robustifying our analysis against unknown but bounded disturbances is challenging (due to the minimax nature of the problem), but a sufficient solution is straightforward and can be done in a number of ways. For instance, the disturbance set could be propagated forward in time, according to which the constraint set could be further eroded to take into account the effect of this uncertainty in future time steps. Alternatively, a pre-stabilization technique [41] could be used, under certain additional assumptions, such that a portion of the control is dedicated to managing the growth of the disturbance set propagation while the remaining portion is employed to keep the trajectory in   K . The simulation based nature of our algorithm readily admits incorporation of time delays. We will make our initial results on SD systems with transport delays available as a technical report for the interested reader. Finally, we also mention that a similar sampling-based technique can be formulated to approximate nonconvex maximal reachable tubes under LTI dynamics, leveraging the fact that the underlying reachable sets (time slices of the tube) are convex. Future work for piecewise affine systems includes finding alternative ways to solve the mixed-integer program so as to improve efficiency.   The extension of the algorithm to nonconvex constraints remains another challenging problem.  A   Appendix: Proof of Proposition 1  To prove Proposition 1 we will introduce and make use of a few well-known results from the literature of random algorithms for estimating the boundary of a convex body.  Lemma 3   ([33]) .   Let   C   be a convex body in   R n   with   ∂ C   C 2 , let   f   :   ∂ C →   R +   be a probability density function defined on   ∂ C , let   P f   be the probability measure defined by   f   , and let   E ( f, N   )   be the expected volume of the convex hull of   N  points chosen randomly on   ∂ C   with respect to   P f   . Then  lim  N   →∞   (vol( C )   −   E ( f, N   ))   N   2  n − 1   =   c n ( C ) ,   (49)  where   c n ( C )   is a constant which depends only on the dimension   n , the distribution of   f   , and the shape of   C .  To use this result, we also need the following lemma, which introduces a fictitious source of error between the outcome of our algorithm and the viability kernel:  Lemma 4   ([42]) .   For every compact convex set   C , there exists a compact convex set   C ′   ⊆ C   whose boundary   ∂ C ′   is  C 2 , and   vol( C )   −   vol( C ′ ) =    smooth   for some arbitrarily small positive scalar    smooth .  Define via Lemma 4, a convex body   Viab ′   that is a   C 2   approximation of   Viab sd  T   ( K , ζ,  )   such that  vol(Viab sd  T   ( K , ζ,  ))   −   vol(Viab ′ ) =    smooth .   (50)  Lemma 5.   Any compact convex set   C   with   r 0   ∈ C   is homeomorphic to   B n  2   ( r 0 ,   1)   [43], and in particular we can define the invertible mapping   m   :   S n − 1 2   ( r 0 )   →   ∂ C   as   m ( r d )   ≡   B ISECTION -F EASIBILITY  ( r 0 ,  F IND -I NTERSECTION - ON -B OUNDARY ( C , ~ r ) ,   C ) , where   ~ r   has origin   r 0   ∈ C   and direction   r d .  We are now ready to prove Proposition 1.  Proof of Proposition 1.   Let   f   ( x )   be the probability distribution used by S AMPLE -R AY ( v 0 )   to generate samples on the unit sphere. Then by using the mapping   m   from Lemma 5, we can perform a change of variables to define a new probability density function   g ( m ( r d ))   on   ∂ Viab ′   [44]. Then by Lemma 3 we have  lim  N   →∞  (  vol(Viab ′ )   −   E ( g, N   )  )  N   2  n − 1   =  ̃ c n (Viab ′ )   (51) for some constant    ̃ c n (Viab ′ )   which depends only on the dimension   n , the distribution   g , and the shape of   Viab ′ . 22
Since by Lemma 4   Viab ′   can be made arbitrarily close to   Viab sd  T   ( K , ζ,  )   (i.e. we can take    smooth   →   0 ), the distribution   g   is mapped almost identically on   ∂ Viab sd  T   ( K , ζ,  )   and we can write (51) as  lim  N   →∞  (  vol(Viab sd  T   ( K , ζ,  ))   −   E ( g, N   )  )  N   2  n − 1  =   lim  N   →∞  (  vol(Viab sd  T   ( K , ζ,  ))   −   vol( V ζ, N   )  )  N   2  n − 1  =  ̃ c n (Viab sd  T   ( K , ζ,  ))   (52) since   E ( g, N   ) = vol( V ζ, N   ) . Taking the limit as   ζ   → ∞   and      →   0   on (52) gives  lim  N   →∞  ζ →∞   → 0  (  vol(Viab sd  T   ( K ))   −   vol( V ζ, N   )  )  N   2  n − 1   =  ̃ c n (Viab sd  T   ( K )) .   (53) Now, by (27), we can replace   vol(Viab sd  T   ( K ))   to get  lim  N   →∞  ζ →∞   → 0  (  vol(Viab sd  T   ( K ))   −    cont ( M δ )   −   vol( V ζ, N   )  )  N   2  n − 1   =  ̃ c n (Viab sd  T   ( K )) .   (54) The shape of   Viab sd  T   ( K )   depends only on the shape of   Viab sd  T   ( K )   and the value   M δ . Thus there exists an appro- priate function   c n   such that    ̃ c n (Viab sd  T   ( K )) =   c n (Viab sd  T   ( K ) , M δ ) . Consequently, we arrive at (30).  References  [1] J. H. Gillula, S. Kaynama, and C. J. Tomlin, “Sampling-based approximation of the viability kernel for high- dimensional linear sampled-data systems,” in   Hybrid Systems: Computation and Control , Berlin, Germany, 2014, pp. 173–182. [2] J.-P. Aubin, A. M. Bayen, and P. Saint-Pierre,   Viability Theory: New Directions , 2nd ed.   Springer Verlag, 2011. [3] F. Blanchini and S. Miani,   Set-Theoretic Methods in Control .   Springer, 2008. [4] G. Goodwin, J. Aguero, M. Cea Garridos, M. Salgado, and J. Yuz, “Sampling and sampled-data models: The interface between the continuous world and digital algorithms,”   IEEE Control Systems Magazine , vol. 33, no. 5, pp. 34–53, 2013. [5] L. Magni, D. M. Raimondo, L. Bossi, C. D. Man, G. D. Nicolao, B. Kovatchev, and C. Cobelli, “Model predictive control of type 1 diabetes: An in silico trial,”   Journal of Diabetes Science and Technology , vol. 1, no. 6, pp. 804– 812, 2007. [6] I. M. Mitchell, A. M. Bayen, and C. J. Tomlin, “A time-dependent Hamilton-Jacobi formulation of reachable sets for continuous dynamic games,”   IEEE Transactions on Automatic Control , vol. 50, no. 7, pp. 947–957, July 2005. [7] I. M. Mitchell, S. Kaynama, M. Chen, and M. Oishi, “Safety preserving control synthesis for sampled data systems,”   Nonlinear Analysis: Hybrid Systems , vol. 10, pp. 63–82, 2013. [8] P. Saint-Pierre, “Approximation of the viability kernel,”   Applied Mathematics and Optimization , vol. 29, no. 2, pp. 187–209, Mar 1994. [9] S. Kaynama and M. Oishi, “A modified Riccati transformation for decentralized computation of the viability kernel under LTI dynamics,”   IEEE Transactions on Automatic Control , vol. 58, no. 11, pp. 2878–2892, 2013. 23
[10] I. M. Mitchell, “Scalable calculation of reach sets and tubes for nonlinear systems with terminal integrators: A mixed implicit explicit formulation,” in   Hybrid Systems: Computation and Control , Chicago, IL, 2011, pp. 103–112. [11] P.-A. Coquelin, S. Martin, and R. Munos, “A dynamic programming approach to viability problems,” in   IEEE Symposium on Adaptive Dynamic Programming and Reinforcement Learning , 2007, pp. 178–184. [12] D. D. Bremner, “On the complexity of vertex and facet enumeration for convex polytopes,” Ph.D. dissertation, Montreal, QC, Canada, 1998. [13] W.-H. Chen, D. J. Ballance, and J. O’Reilly, “Optimisation of attraction domains of nonlinear MPC via LMI methods,” in   American Control Conference , Arlington, VA, 2001, pp. 3067–3072. [14] A. Alessio, M. Lazar, A. Bemporad, and W. Heemels, “Squaring the circle: An algorithm for generating polyhe- dral invariant sets from ellipsoidal ones,”   Automatica , vol. 43, no. 12, pp. 2096–2103, 2007. [15] J. Maidens, S. Kaynama, I. M. Mitchell, M. Oishi, and G. A. Dumont, “Lagrangian methods for approximating the viability kernel in high-dimensional systems,”   Automatica , vol. 49, no. 7, pp. 2017–2029, 2013. [16] S. Kaynama, J. Maidens, M. Oishi, I. M. Mitchell, and G. A. Dumont, “Computing the viability kernel using maximal reachable sets,” in   Hybrid Systems: Computation and Control , Beijing, 2012, pp. 55–63. [17] C. Le Guernic and A. Girard, “Reachability analysis of linear systems using support functions,”   Nonlinear Anal- ysis: Hybrid Systems , vol. 4, no. 2, pp. 250–262, 2010. [18] G. Frehse, C. Le Guernic, A. Donze, S. Cotton, R. Ray, O. Lebeltel, R. Ripado, A. Girard, T. Dang, and O. Maler, “SpaceEx: Scalable verification of hybrid systems,” in   23rd International Conference on Computer Aided Veri- fication , G. Gopalakrishnan and S. Qadeer, Eds.   Springer, 2011, pp. 1–16. [19] A. B. Kurzhanski and P. Varaiya, “Ellipsoidal techniques for reachability analysis: internal approximation,”  Systems Control Letters , vol. 41, pp. 201–211, 2000. [20] A. Majumdar and R. Tedrake, “Robust online motion planning with regions of finite time invariance,” in   Al- gorithmic Foundations of Robotics , E. Frazzoli, T. Lozano-Perez, N. Roy, and D. Rus, Eds.   Springer Berlin Heidelberg, 2013, pp. 543–558. [21] M. Tobenkin, I. Manchester, and R. Tedrake, “Invariant funnels around trajectories using sum-of-squares pro- gramming,” in   IFAC World Congress , vol. 18, Milano, Italy, 2011, pp. 9218–9223. [22] R. Tedrake, I. R. Manchester, M. Tobenkin, and J. W. Roberts, “LQR-trees: Feedback motion planning via sums-of-squares verification,”   International Journal of Robotics Research , vol. 29, no. 8, pp. 1038–1052, 2010. [23] W. Tan and A. Packard, “Searching for control lyapunov functions using sums of squares programming,” in  Annual Allerton Conference , 2004, pp. 210–219. [24] S. Prajna and A. Jadbabaie, “Safety verification of hybrid systems using barrier certificates,” in   Hybrid Systems: Computation and Control , 2004, pp. 477–492. [25] D. Henrion and M. Korda, “Convex computation of the region of attraction of polynomial control systems,”  preprint arXiv:1208.1751 , 2012. [26] A. Majumdar, R. Vasudevan, M. M. Tobenkin, and R. Tedrake, “Technical report: Convex optimization of non- linear feedback controllers via occupation measures,”   preprint arXiv:1305.7484 , 2013. [27] C. Moler and C. Van Loan, “Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later,”   SIAM Review , vol. 45, no. 1, pp. 3–49, 2003. 24
[28] N. J. Higham, “The scaling and squaring method for the matrix exponential revisited,”   SIAM Journal of Matrix Analysis and Applications , vol. 26, no. 4, pp. 1179–1193, 2005. [29] C. Standish, “Truncated Taylor series approximation to the state transition matrix of a continuous parameter finite Markov chain,”   Linear Algebra and its Applications , vol. 12, no. 2, pp. 179–183, 1975. [30] M. Liou, “A novel method of evaluating transient response,”   Proceedings of the IEEE , vol. 54, no. 1, pp. 20–23, 1966. [31] I. B ́ ar ́ any and Z. F ̈ uredi, “Computing the volume is difficult,”   Discrete Computational Geometry , vol. 2, no. 1, pp. 319–326, 1987. [32] M. Dyer, “Computing the volume of convex bodies: a case where randomness provably helps,”   Probabilistic Combinatorics and its Applications , vol. 44, pp. 123–170, 1991. [33] C. Sch ̈ utt and E. Werner, “Polytopes with vertices chosen randomly from the boundary of a convex body,” in  Geometric Aspects of Functional Analysis , ser. LNM 1807, V. Milman and G. Schechtman, Eds.   Springer Berlin, 2003, pp. 241–422. [34] J. L ̈ ofberg, “YALMIP: a toolbox for modeling and optimization in MATLAB,” in   Computer Aided Control System Design , 2004, pp. 284–289. [35] M. Kvasnica, P. Grieder, M. Baoti ́ c, and M. Morari, “Multi-Parametric Toolbox (MPT),” in   Hybrid Systems: Computation and Control, LNCS 2993 , R. Alur and G. J. Pappas, Eds.   Berlin, Germany: Springer, 2004, pp. 448–462. [36] K. V. Mardia and P. E. Jupp,   Directional Statistics , 2nd ed.   John Wiley & Sons, Inc, 2000. [37] M. E. Dyer and A. M. Frieze, “On the complexity of computing the volume of a polyhedron,”   SIAM Journal of Computing , vol. 17, no. 5, pp. 967–974, 1988. [38] S. P. Boyd and L. Vandenberghe,   Convex optimization .   Cambridge University Press, 2004. [39] S. Kaynama and C. J. Tomlin, “Benchmark: Flight envelope protection in autonomous quadrotors,” in   Interna- tional Workshop on Applied Verification for Continuous and Hybrid Systems, Part of CPSWeek , Berlin, Germany, April 2014, http://cps-vo.org/group/ARCH. [40] I. Cowling, O. Yakimenko, J. Whidborne, and A. Cooke, “Direct method based control system for an autonomous quadrotor,”   Journal of Intelligent & Robotic Systems , vol. 60, pp. 285–316, 2010. [41] L. Chisci, J. Rossiter, and G. Zappa, “Systems with persistent disturbances: predictive control with restricted constraints,”   Automatica , vol. 37, pp. 1019–1028, 2001. [42] M. Ghomi, “Optimal Smoothing for Convex Polytopes,”   Bulletin of the London Mathematical Society , vol. 36, no. 4, pp. 483–492, July 2004. [43] G. E. Bredon,   Topology and Geometry .   Springer, 1993, vol. 139. [44] J. Pitman,   Probability .   New York: Springer-Verlag, 1993. 25