The Curious Case of Integrator Reach Sets, Part I: Basic Theory
TThe Curious Case of Integrator Reach Sets, Part I: Basic Theory
Shadi Haddad, Abhishek Halder
Abstract —This is the first of a two part paper investigatingthe geometry of the integrator reach sets, and the applicationsthereof. In this Part I, we establish that this compact convex setis semialgebraic, translated zonoid, and not a spectrahedron. Wederive the parametric as well as the implicit representation ofthe boundary of this reach set. We also deduce the closed formformula for the volume and diameter of this set, and discuss theirscaling with state dimension and time. We point out that theseresults may be utilized in benchmarking the performance of thereach set over-approximation algorithms.
Keywords:
Reach set, integrator, convex geometry, set-valueduncertainty. I. I
NTRODUCTION
Integrators with bounded controls are ubiquitous in systems-control. They appear as Brunovsky normal forms for thefeedback linearizable nonlinear systems. They also appearfrequently as benchmark problems to demonstrate the per-formance of the reach set computation algorithms. Despitetheir prominence, specific results on the geometry of theintegrator reach sets are not available in the systems-controlliterature. Broadly speaking, the existing results come in twoflavors. On one hand, very generic statements are known,e.g., these reach sets are compact convex sets whenever theset of initial conditions is compact convex, and the controlstake values from a compact (not necessarily convex) set [1].On the other hand, several numerical toolboxes [2], [3] areavailable for tight outer approximation of the reach sets overcomputationally benign geometric families such as ellipsoidsand zonotopes. The lack of concrete geometric results implythe absence of ground truth when comparing the efficacy ofdifferent algorithms, and one has to content with graphical orstatistical (e.g., Monte Carlo) comparisons.Building on the preliminary results in [4], this paper un-dertakes a systematic study of the integrator reach sets. Inparticular, we answer the following basic questions:
Q1. what kind of compact convex sets are these (Section III)?
Q2. how big are these sets (Section IV)?
Q3. how these results on the geometry of integrator reach setscan be applied in practice (Section V)?We consider the integrator dynamics having d states and m inputs with relative degree vector r = ( r , r , . . . , r m ) (cid:62) ∈ Z m + (vector of positive integers). The dynamics is given by ˙ x = Ax + Bu , x ∈ R d , u ∈ U ⊂ R m , (1)where r + r + ... + r m = d , the set U is compact, and A := blkdiag( A , ..., A m ) , B := blkdiag( b , ..., b m ) , (2a) A j := (cid:16) r j × | e r j | e r j | ... | e r j r j − (cid:17) , b j := e r j r j . (2b) Shadi Haddad and Abhishek Halder are with the Department of Ap-plied Mathematics, University of California, Santa Cruz, CA 95064, USA, { shhaddad,ahalder } @ucsc.edu In (2a), the symbol blkdiag ( · ) denotes a block diagonalmatrix whose arguments constitute its diagonal blocks. In (2b),the notation r j × stands for the r j × column vector ofzeros, and e (cid:96)k denotes the k th basis (column) vector in R (cid:96) for k ≤ (cid:96) .Let R ( X , t ) denote the forward reach set of (1) at time t > , starting from a given compact convex set of initialconditions X ⊂ R d , i.e., R ( X , t ) := (cid:8) x ( t ) ∈ R d | ˙ x = Ax + Bu , x ( t = 0) ∈ X , u ∈ U (cid:9) . (3)In words, R ( X , t ) is the set of all states that the controlleddynamics (1) can reach at time t > , starting from the set X at t = 0 , with control u ( t ) ∈ U compact. Formally, R ( X , t ) = exp( t A ) X (cid:117) (cid:90) t exp (( t − τ ) A ) B U d τ = exp( t A ) X (cid:117) (cid:90) t exp ( s A ) B U d s, (4)where (cid:117) denotes the Minkowski sum. The set-valued Aumannintegral [5] in (4) is defined for any point-to-set function F ( · ) ,as (cid:90) t F ( s )d s := lim ∆ ↓ (cid:98) t/ ∆ (cid:99) (cid:88) i =0 ∆ F ( i ∆) , (5)where the summation symbol Σ denotes the Minkowski sum,and (cid:98)·(cid:99) is the floor operator; see e.g., [1]. Our objective is tostudy the geometric aspects of (4) in detail.This paper significantly expands our preliminary work [4]:here we consider multi-input integrators as opposed to thesingle input case considered in [4]. Even for the single inputcase, while [4, Thm. 1] derived an exact formula for thevolume of the reach set, that formula involved limit and nestedsums, and in that sense, was not really a closed-form formula[6] – certainly not amenable for numerical computation. In thispaper, we derive closed-form formula for the general multi-input case. In addition to these extensions and generalizations,the present paper addresses previously unexplored directions:the scaling laws for the volume and diameter of integratorreach sets, exact parametric and implicit equations for theboundary, and the classification of these sets.The paper is structured as follows. After reviewing somepreliminary concepts in Sec. II, the results on taxonomy andthe boundary of the integrator reach set are provided in Sec.III. The results on the size of the reach set are collected in Sec.IV. The application of these results for benchmarking the reachset over-approximation algorithms are discussed in Sec. V. Allproofs are deferred to the Appendix. Sec. VI summarizes theresults in this paper, and outlines the directions pursued in itssequel Part II. a r X i v : . [ ee ss . S Y ] F e b I. P
RELIMINARIES
In the following, we summarize some preliminaries whichwill be useful in the main body and in the Appendix.
1) State transition matrix:
For ≤ s < t , the statetransition matrix Φ ( t, s ) associated with (2) is Φ ( t, s ) ≡ exp( A ( t − s ))= blkdiag (exp( A ( t − s )) , . . . , exp( A m ( t − s ))) , with each diagonal block being upper triangular, given by exp( A j ( t − s )) := ( t − s ) (cid:96) − k ( (cid:96) − k )! for k ≤ (cid:96), otherwise , (6)for k, (cid:96) = 1 , . . . , r j , for each j = 1 , . . . , m . In particular, thediagonal entries in (6) are unity.
2) Support function:
The support function h K ( · ) of acompact convex set K ⊂ R d , is given by h K ( y ) := sup x ∈K {(cid:104) y , x (cid:105) | y ∈ R d } , (7)where (cid:104)· , ·(cid:105) denotes the standard Euclidean inner product.Geometrically, h K ( y ) gives the signed distance of the support-ing hyperplane of K with outer normal vector y , measuredfrom the origin. Furthermore, the supporting hyperplane at x bdy ∈ ∂ K is (cid:104) y , x bdy (cid:105) = h K ( y ) , and we can write K = (cid:8) x ∈ R d | (cid:104) y , x (cid:105) ≤ h K ( y ) for all y ∈ R d (cid:9) . For more details on the support function, we refer the readersto [7, Ch. V].The support function h K ( y ) uniquely determines the set K .Given matrix-vector pair ( Γ , γ ) ∈ R d × d × R d , the supportfunction of the affine transform Γ K + γ is h Γ K + γ ( y ) = (cid:104) y , γ (cid:105) + h K ( Γ (cid:62) y ) . (8)Given a function f : R d (cid:55)→ R ∪{ + ∞} , its Legendre-Fenchelconjugate is f ∗ ( y ) := sup x ∈ domain( f ) {(cid:104) y , x (cid:105) − f ( x ) | y ∈ R d } . (9)From (7)-(9), it follows that h K ( y ) is the Legendre-Fenchelconjugate of the indicator function K ( x ) := (cid:40) if x ∈ K , + ∞ otherwise.Since the indicator function of a convex set is a convexfunction, the biconjugate ∗∗K ( · ) = h ∗K ( · ) = K ( · ) . This willbe useful in Section III.From (2a), the system matrices are block diagonal. Thus,each of the m single input integrator dynamics with r j dimen-sional state subvectors for j = 1 , . . . , m , are decoupled fromeach other. Hence R ( X , t ) ⊂ R d is the Cartesian product ofthese single input integrator reach sets: R j ( X , t ) ⊂ R r j for j = 1 , . . . , m , i.e., R = R × R × . . . × R m . (10)In what follows, we will sometimes exploit that (10) may alsobe written as a Minkowski sum R (cid:117) . . . (cid:117) R m . Notice that the decoupled dynamics also allows us to write a Mikowskisum decomposition for the set of initial conditions X = X (cid:117) . . . (cid:117) X m , and accordingly, the initial condition subvectors x j ∈ X j ⊂ R r j for j = 1 , . . . , m . Thus x = ( x , . . . , x m ) (cid:62) .Since the support function of the Minkowski sum is equalto the sum of the support functions, we have h R ( X ,t ) ( y ) = m (cid:88) j =1 h R j ( X j ,t ) ( y j ) . (11)To proceed further, we introduce some notations. Since U is compact, for each j = 1 , . . . , m , let α j := min u ∈U u j , β j := max u ∈U u j , (12)that is, α j and β j are the component-wise minimum andmaximum, respectively, of the input vector. Furthermore, let µ j := β j − α j , ν j := β j + α j , (13)and introduce ξ ( s ) := µ ξ ( s ) ... µ m ξ m ( s ) , ξ j ( s ) := s r j − / ( r j − s r j − / ( r j − ... s , (14)for j = 1 , . . . , m . Also, let ζ ( t , t ) = µ ζ ( t , t ) µ ζ ( t , t ) ... µ m ζ m ( t , t ) , ζ j ( t , t ) := (cid:90) tt ξ j ( s ) d s ∈ R r j , (15)for j = 1 , . . . , m . When t = 0 , we simplify the notations as ζ ( t ) := ζ (0 , t ) , ζ j ( t ) := ζ j (0 , t ) for all j = 1 , . . . , m. (16)The following result (proof in Appendix A) will come inhandy in the ensuing development. Theorem 1.
For compact convex X ⊂ R d , and compact U ⊂ R m , the support function of the reach set (4) is h R ( X ,t ) ( y ) = m (cid:88) j =1 (cid:26) sup x j ∈X j (cid:104) y j , exp ( t A ) x j (cid:105) + ν j (cid:104) y j , ζ j ( t ) (cid:105) + µ j (cid:90) t |(cid:104) y j , ξ j ( s ) (cid:105)| d s (cid:27) . (17)
3) Polar dual:
The polar dual K ◦ of any non-empty set K ⊂ R d is given by K ◦ := { y ∈ R d | (cid:104) y , x (cid:105) ≤ for all x ∈ K} . (18)From this definition, it is immediate that K ◦ contains theorigin, and is a closed convex set. The bipolar ( K ◦ ) ◦ =closure (conv ( K ∪ { } )) where conv( · ) denotes the convexhull. Thus, if K is compact convex and contains the origin,then we have the involution ( K ◦ ) ◦ = K . From (7) and(18), notice that K ◦ is the unit support function ball , i.e., K ◦ = { y ∈ R d | h K ( y ) ≤ } . In Sec. III-D, we will mentionsome properties of the polar dual of the integrator reach set. ) Vector measure: Let F be a σ -field of the subsets ofa set. A countably additive mapping (cid:101) µ : F (cid:55)→ R d is termeda vector measure . Here, “countably additive” means that forany sequence { Ω i } ∞ i =1 of disjoint sets in F such that theirunion is in F , we have (cid:101) µ ( ∪ ∞ i =1 Ω i ) = (cid:80) ∞ i =1 (cid:101) µ (Ω i ) < ∞ .Some of the early investigations of vector measures were dueto Liapounoff [8] and Halmos [9]; relatively recent referencesare [10], [11].
5) Zonotope: A zonotope Z ⊂ R d is a finite Minkowskisum of closed line segments or intervals { I j } rj =1 where theseintervals are imbedded in the ambient Euclidean space R d .Explicitly, Z := I (cid:117) . . . (cid:117) I r = (cid:26) x ∈ R d | x = r (cid:88) j =1 x j , x j ∈ I j , j = 1 , . . . , r (cid:27) . Thus, a zonotope is the range of an atomic vector measure.Alternatively, a zonotope can be viewed as the affine imageof the unit cube. A compact convex polytope is a zonotopeif and only if all its two dimensional faces are centrallysymmetric [12, p. 182]. For instance, the cross polytope { x ∈ R d | (cid:107) x (cid:107) ≤ } , is not a zonotope. Standard referenceson zonotope include [13], [14], [15, Ch. 2.7].The set of zonotopes is closed under affine image andMinkowski sum, but not under intersection. In the systems-control literature, a significant body of work exists on com-putationally efficient over-approximation of reach sets viazonotopes [16]–[18] and its variants such as zonotope bundles[19], constrained zonotopes [20], complex zonotopes [21], andpolynomial zonotopes [22], [23].
6) Variety and ideal:
Let p , . . . , p n ∈ R [ x , . . . , x d ] , thevector space of real-valued d -variate polynomials. The (affine) variety V R [ x ,...,x d ] ( p , . . . , p n ) is the set of all solutions of thesystem p ( x , x , . . . , x d ) = . . . = p n ( x , x , . . . , x d ) = 0 .Given p , . . . , p n ∈ R [ x , . . . , x d ] , the set I := (cid:26) n (cid:88) i =1 α i p i | α , . . . , α n ∈ R [ x , . . . , x d ] (cid:27) is called the ideal generated by p , . . . , p n . We write thissymbolically as I = (cid:104)(cid:104) p , . . . , p n (cid:105)(cid:105) . Roughly speaking, (cid:104)(cid:104) p , . . . , p n (cid:105)(cid:105) is the set of all polynomial consequences of thegiven system of n polynomial equations in d indeterminates.We refer the readers to [24, Ch. 1] for detailed exposition ofthese concepts.III. T AXONOMY AND B OUNDARY
For X ⊂ R d compact convex, it is well-known [1, Sec. 2]that R ( X , t ) is a compact convex set for all t > . However,it is not immediate what kind of convex set R is, even forsingleton X ≡ { x } .In this Section, we examine the question “what type ofcompact convex set R ( { x } , t ) is” from several points ofview. In doing so, we also derive the equations for theboundary ∂ R ( { x } , t ) .Notice that for non-singleton X , the taxonomy questionis not well-posed since the classification then will depend on X . Also, setting X ≡ { x } in (4), it is apparentthat R ( { x } , t ) is a translation of the set-valued integral in(4). Thus, classifying R ( { x } , t ) amounts to classifying thesecond summand in (4). A. R ( { x } , t ) is a Zonoid A zonoid is a compact convex set that is defined as therange of an atom free vector measure (see Sec. II-4). Affineimage of a zonoid is a zonoid. Minkowski sum of zonoidsis also a zonoid. We refer the readers to [25]–[27], [28, Sec.I] for more details on the properties of a zonoid. By slightabuse of nomenclature, in this paper we use the term zonoidup to translation, i.e., we refer to the translation of zonoidsas zonoids (instead of using another term such as “zonoidaltranslates”).Let us mention a few examples. Any compact convexsymmetric set in R is a zonoid. In dimensions three or more,all (cid:96) p norm balls for p ≥ are zonoids.An alternative way to think about the zonoid is to viewit as the limiting set (convergence with respect to the two-sided Hausdorff distance, see e.g., [4, Appendix B]) of theMinkowski sum of line segments, i.e., the limit of a sequenceof zonotopes [13], [14], [25]. We will use this viewpoint inSec. IV-A. Our main result in this subsection is the following. Theorem 2.
The reach set (4) with X ≡ { x } is a zonoid. To appreciate Theorem 2 via the limiting viewpoint men-tioned before, one can write R ( { x } , t ) = exp( t A ) x + m (cid:88) j =1 ν j ζ j ( t ) (cid:124) (cid:123)(cid:122) (cid:125) first term (cid:117) m (cid:88) j =1 lim n →∞ n (cid:88) i =0 tn µ j ξ j ( t i ) [ − , (cid:124) (cid:123)(cid:122) (cid:125) second term , (19)where all summation symbols denote Minkowski sums. Thefirst term in (19) denotes a translation. In the secondterm, the outer summation over index j arises by writingthe Cartesian product (10) as the Minkowski sum R (cid:117) . . . (cid:117) R m . Furthermore, uniformly discretizing [0 , t ] into n subintervals [( i − t/n, it/n ) , i = 1 , . . . , n , we write (cid:82) t exp( s A j ) b j [ − µ j , µ j ]d s as the limit of the Minkowski sumover index i . Geometrically, the innermost summands in thesecond term denote non-uniformly rotated and scaled lineintervals in R j . In other words, the second term in (19) isa Minkowski sum of m sets, each of these sets being the limitof a sequence of sets {Z n } comprising of zonotopes Z n := n (cid:88) i =0 tn µ j ξ j ( t i ) [ − , , which are the Minkowski sum of n + 1 line segments.Since lim n →∞ Z n is a zonoid, the second term in (19) isa Minkowski sum of m zonoids, and is therefore a zonoid[25, Thm. 1.5]. The entire right hand side of (19), then, istranslation of a zonoid, and hence a zonoid. emark 1. If X ⊂ R d is not singleton, but instead a zonoid,then R ( X , t ) is still a (translated) zonoid. To see this, noticefrom (4) and (17) that R ( X , t ) = exp( t A ) X (cid:117) R ( { } , t ) , (20) and that exp( t A ) X , being linear image of a zonoid, is azonoid [25, Lemma 1.4]. Thus, (20) being Minkowski sum ofzonoids, is a zonoid too [25, Thm. 1.5], up to translation.B. R ( { x } , t ) is Semialgebraic A set in R d is called basic semialgebraic if it can bewritten as a finite conjunction of polynomial inequalities andequalities, the polynomials being in R [ x , . . . , x d ] . Finiteunion of basic semialgebraic sets is called a semialgebraicset . A semialgebraic set need not be basic semialgebriac; seee.g., [29, Example 2.2].Semialgebraic sets are closed under finitely many unionsand intersections, complement, topological closure, polyno-mial mapping including projection [30], [31], and Cartesianproduct. For details on semialgebraic sets, we refer the readersto [32, Ch. 2]; see [33, Appendix A.4.4] for a short summary.In Proposition 1 below, we derive a parametric representa-tion of x bdy ∈ ∂ R ( { x } , t ) , the boundary of the reach set.Then we use this representation to establish semialgebraicityof R ( { x } , t ) in Theorem 4 that follows. Proposition 1.
For relative degree vector r = ( r , . . . , r m ) (cid:62) ,and fixed x ∈ R d comprising of subvectors x j ∈ R r j where j = 1 , . . . , m , consider the reach set (4) with singleton X ≡ { x } and compact U ⊂ R m . For j = 1 , . . . , m , define µ , . . . , µ m and ν , . . . , ν m as in (12)-(13). Let the indicatorfunction k ≤ (cid:96) := 1 for k ≤ (cid:96) , and := 0 otherwise. Then thecomponents of x bdy = x bdy x bdy ... x bdy m ∈ ∂ R ( { x } , t ) , x bdy j ∈ R r j , j = 1 , . . . , m, admit parameterization x bdy j ( k ) = r j (cid:88) (cid:96) =1 k ≤ (cid:96) t (cid:96) − k ( (cid:96) − k )! x j ( (cid:96) ) + ν j t r j − k +1 ( r j − k + 1)! ± µ j ( r j − k + 1)! (cid:26) ( − r j − t r j − k +1 + 2 r j − (cid:88) q =1 ( − q +1 s r j − k +1 q (cid:27) , (21) where x bdy j ( k ) denotes the k th component of the j th subvector x bdy j for k = 1 , . . . , r j . The parameters ( s , s , . . . , s r j − ) satisfy ≤ s ≤ s ≤ . . . ≤ s r j − ≤ t . The following is an immediate consequence of Proposition 1.
Corollary 3.
The single input integrator reach set R j ( { x } , t ) ⊂ R r j has two bounding surfaces for each j =1 , . . . , m . In other words, there exist p upper j , p lower j : R r j (cid:55)→ R such that R j ( { x } , t ) = { x ∈ R r j | p upper j ( x ) ≤ , p lower j ( x ) ≤ } , Fig. 1:
The “almond-shaped” integrator reach set R ( { x } , t ) ⊂ R with d = 3 , m = 1 , x = (0 . , . , . (cid:62) , U ≡ [ α, β ] = [ − , at t = 2 . . The wireframes correspond to the upper and lower surfaces. with boundary ∂ R j ( { x } , t ) = { x ∈ R r j | p upper j ( x ) = 0 } ∪{ x ∈ R r j | p lower j ( x ) = 0 } . During the proof of Theorem 4 stated below, it will turnout that in fact p upper j , p lower j ∈ R (cid:2) x , . . . , x r j (cid:3) for all j =1 , . . . , m . In words, p upper j , p lower j are real algebraic hypersur-faces for all j = 1 , . . . , m .Let us exemplify the parameterization (21) for the case r =( r , r ) (cid:62) = (2 , (cid:62) . In this case, x bdy (1) x bdy (2) = (cid:32) x (1) + t x (2) + ν ( t / ± µ (cid:0) s − t / (cid:1) x (2) + ν t ± µ (2 s − t ) (cid:33) , (22)and x bdy (1) x bdy (2) x bdy (3) = x (1) + t x (2) + ( t / x (3)+ ν ( t / ± µ (cid:0) t / s / − s / (cid:1) x (2) + t x (3) + ν ( t / ± µ (cid:0) t / s / − s / (cid:1) x (3) + ν t ± µ ( t + 2 s − s ) . (23)In (22), taking plus (resp. minus) signs in each of componentgives the parametric representation of the curve p upper = 0 (resp. p lower = 0 ). These curves are as in [4, Fig. 1(a)], andtheir union defines ∂ R . We note that the parameterization(22) appeared in [34, p. 111].Likewise, in (23), taking plus (resp. minus) signs in each ofcomponent gives the parametric representation of the surface p upper ( x ) = 0 (resp. p lower = 0 ). The resulting set R is thetriple integrator reach set, and is shown in Fig. 1. a) R with x = (0 . , . (cid:62) . (b) R ◦ with x = (0 . , . (cid:62) .(c) R with x = (0 . , . (cid:62) . (d) R ◦ with x = (0 . , . (cid:62) . Fig. 2:
The double integrator reach set R ( { x } , t ) and its polar dual ( R ( { x } , t )) ◦ for different x at t = 2 , U ≡ [ α, β ] = [ − , . Thecurves p upper , p lower defining the reach set boundary (see Corollary 3and the discussion thereafter) are shown too. Now we come to the main result of this subsection.
Theorem 4.
The reach set (4) with X ≡ { x } is semialge-braic. Let us illustrate the bounding curves and surfaces for (22)and (23) respectively, in the implicit form. Eliminating theparameter s from (22) reveals that p upper , p lower are parabolas.Specifically, p upper ( x bdy (1) , x bdy (2)) = 14 (cid:32) x bdy (2) − x (2) − ν tµ + t (cid:33) − x bdy (1) − x (1) − t x (2) − ν t µ − t , (24a) p lower ( x bdy (1) , x bdy (2)) = − (cid:32) − x bdy (2) − x (2) − ν tµ + t (cid:33) − x bdy (1) − x (1) − t x (2) − ν t µ + t . (24b)Similarly, eliminating the parameters s , s from (23) revealsthat p upper , p lower are quartic polynomials, and we have p upper ( x bdy (1) , x bdy (2) , x bdy (3))= 116 (cid:32) x bdy (3) − x (3) − ν tµ − t (cid:33) + 3 (cid:32) x bdy (2) − x (2) − t x (3) − ν t µ − t (cid:33) − (cid:32) x bdy (1) − x (1) − t x (2) − t x (3) − ν t µ − t (cid:33) × (cid:32) x bdy (3) − x (3) − ν tµ − t (cid:33) , (25a) p lower ( x bdy (1) , x bdy (2) , x bdy (3))= − (cid:32) x bdy (3) − x (3) − ν tµ + t (cid:33) − (cid:32) x bdy (2) − x (2) − t x (3) − ν t µ + t (cid:33) + 6 (cid:32) x bdy (1) − x (1) − t x (2) − t x (3) − ν t µ + t (cid:33) × (cid:32) x bdy (3) − x (3) − ν tµ + t (cid:33) . (25b)A natural question is whether one can generalize the implic-itizations (24), (25). This is what we address next. C. Implicitization of ∂ R ( { x } , t ) To derive the implicit equations for the bounding al-gebraic hypersurfaces p upper j , p lower j ∈ R (cid:2) x , . . . , x r j (cid:3) forall j = 1 , . . . , m , we need to eliminate the parameters (cid:0) s , s , . . . , s r j − (cid:1) from (21). For this purpose, it is helpfulto write (21) succinctly as ρ ± j,k = r j − (cid:88) q =1 ( − q +1 s r j − k +1 q , k = 1 , . . . , r j , (26)where ρ ± j,k := ( r j − k + 1)!2 µ j (cid:26) x bdy j ( k ) − r j (cid:88) (cid:96) =1 k ≤ (cid:96) t (cid:96) − k ( (cid:96) − k )! x j ( (cid:96) ) (cid:27) − (cid:26) ± ( − r j − t r j − k +1 + ν j µ j t r j − k +1 (cid:27) . (27)To simplify the rather unpleasant notation ρ ± j,k , we will onlyaddress the m = 1 case. In (26), this allows us to replace r j by d , and to drop the subscript j from the ρ ’s. This doesnot invite any loss of generality in terms of implicitizationsince post derivation, we can replace d by r j to recover therespective p j ’s.With slight abuse of notation, we will also drop the super-script ± from the ρ ’s in (26). Recall that the plus (resp. minus)superscript in the ρ ’s indicates p upper j (resp. p lower j ). From (27),it is clear that in either case, the ρ j,k is affine in x bdy j ( k ) ,which is the k th coordinate of the boundary point for the j thblock. Importantly, for k = 1 , . . . , r j , the quantity ρ j,k doesnot depend on any other component of the boundary pointthan the k th component. Again, the plus-minus superscriptscan be added back post implicitization.Thus, the notationally simplified version of (26) that sufficesfor implicitization, is ρ k = d − (cid:88) q =1 ( − q +1 s d − k +1 q , k = 1 , . . . , d, (28)hich is a system of d homogeneous polynomials in variables ( s , s , . . . , s d − ) . The objective is to derive the implicitizedpolynomial ℘ ( ρ , ρ , . . . , ρ d ) associated with (28).When d = 2 , the parameterization (28) becomes ρ = s ,ρ = s , and we get degree 2 implicitized polynomial ℘ ( ρ , ρ ) = ρ − ρ = 0 . (29)For k = 1 , , substituting for the ρ , ρ in (29) from (27) withappropriate plus-minus signs recovers (24).When d = 3 , the parameterization (28) becomes ρ = s − s ,ρ = s − s ,ρ = s − s , and using elementary algebra, we get degree 4 implicitizedpolynomial ℘ ( ρ , ρ , ρ ) = ρ − ρ ρ + 3 ρ = 0 . (30)As before, for k = 1 , , , substituting for the ρ , ρ , ρ in (30)from (27) with appropriate plus-minus signs recovers (25).However, for d = 4 or higher, it is practically impossible toderive the implicitization via brute force algebra.A principled way to implicitize (28) is due to G. Zaimi [35],and starts with defining λ k := ρ d − k +1 for k = 1 , . . . , d . In-troduce the sequence A k ( s , s , . . . , s d − ) via the generatingfunction (see e.g., [36, Ch. 1]) F ( τ ) = (cid:88) k ≥ A k τ k = (1 − s τ )(1 − s τ ) · · · (1 − s τ )(1 − s τ ) · · · . (31)Taking the logarithmic derivative of (31), and then using thegenerating functions (1 − s q τ ) − = (cid:80) k ≥ ( s q τ ) k for all q =1 , . . . , d − , yields F (cid:48) ( τ ) F ( τ ) = − s (cid:88) k ≥ ( s τ ) k + s (cid:88) k ≥ ( s τ ) k − s (cid:88) k ≥ ( s τ ) k + . . . . (32)Integrating (32) with respect to τ , we obtain F ( τ ) = exp (cid:32) − d (cid:88) k =1 λ k k τ k (cid:33) . (33)Equating (31) and (33) allows us to compute A k as a degree k polynomial of the λ ’s.On the other hand, since the generating function (31) isa rational function with denominator polynomial of degree δ := (cid:98) d − (cid:99) , the following Hankel determinant vanishes * det[ A d − δ + i + j ] δi,j =0 = 0 . (34)Substituting the A k ’s obtained as degree k polynomials of the λ ’s into (34) gives an implicit polynomial in indeterminate ( λ , . . . , λ d ) of degree ( δ + 1)( d − δ ) . Finally, reverting back * This result goes back to Kronecker [37]. See also [38, p. 5, Lemma III]. the λ ’s to the ρ ’s result in the desired implicit polynomial ℘ ( ρ , ρ , . . . , ρ d ) , which is also of degree ( δ + 1)( d − δ ) .For instance, when d = 3 , the relation (34) becomes det (cid:18)(cid:20) A A A A (cid:21)(cid:19) = 0 . (35)In this case, equating (31) and (33) gives A = − λ , A = 12 λ − λ , A = − λ + 12 λ λ − λ . Substituting these back in (35) yields the quartic polyno-mial λ + 3 λ − λ λ = 0 , which under the mapping ( λ , λ , λ ) (cid:55)→ ( ρ , ρ , ρ ) recovers (30), and thus (25).In summary, (34) is the desired implicitization of thebounding hypersurfaces of the single input integrator reachset (up to the change of variables). The Cartesian product ofthese implicit hypersurfaces gives the implicitization in themulti input case. D. Dual of R ( { x } , t ) From convex geometry standpoint, it is natural to ask whatkind of characterization is possible for the polar dual (see Sec.II-3) of the integrator reach set R . We know in general that R ◦ will be compact convex. Depending on the choice of x , U and t , the set R ( { x } , t ) may not contain the origin, and thusthe bipolar ( R ( { x } , t )) ◦◦ = closure (conv ( R ( { x } , t ) ∪ { } )) , that is, we do not have the involution in general.Furthermore, since R ( { x } , t ) is semialgebraic from Sec.III-B, so must be its polar dual ( R ( { x } , t )) ◦ ; see e.g., [33,Ch. 5, Sec. 5.2.2].We also know from Sec. III-A that R ( { x } , t ) is a zonoid.However, the polar of a zonoid is not a zonoid in general [39],[40], and we should not expect ( R ( { x } , t )) ◦ to be one. Fig. 2shows R ( { x } , t ) and ( R ( { x } , t )) ◦ for the double integrator( d = 2 , m = 1 ). E. Summary of Taxonomy
So far we explained that the compact convex set R ( { x } , t ) is semialgebraic, and a translated zonoid. Two well-knownsubclasses of convex semialgebraic sets are the spectrahedra and the spectrahedral shadows . The spectrahedra, a.k.a. linearmatrix inequality (LMI) representable sets are affine slicesof the symmetric positive semidefinite cone. The spectrahe-dral shadows, a.k.a. lifted LMI or semidefinite representablesets are the projections of spectrahedra. The spectrahedralshadows subsume the class of spectrahedra; e.g., the set { ( x , x ) ∈ R | x + x ≤ } is a spectrahedral shadowbut not a spectrahedron. The polar duals of spectrahedra arespectrahedral shadows [33, Ch. 5, Sec. 5.5].We note that the integrator reach set is not a spectrahedron.To see this, we resort to the contrapositive of [41, Thm. 3.1].Specifically, the number of intersections made by a generic linepassing through an interior point of the d -dimensional integra-tor reach set with its real algebraic boundary is not equal tothe degree of the bounding algebraic hypersurfaces, the latter a) Real algebraic curves p upper ,p lower for the double integrator. (b) Real algebraic surfaces p upper ,p lower for the triple integrator. Fig. 3:
The bounding polynomials for the double and triple integratorreach sets at t = 0 . with x = and µ = 1 . R d Compactand convex SemialgebraicZonoids Spectrahedra(lifted LMI representable)(LMI representable)Spectrahedral shadow
Fig. 4:
The summary of the taxonomy for the integrator reach set. we know from Sec. III-C to be ( (cid:98) d − (cid:99) + 1)( d − (cid:98) d − (cid:99) ) . Inother words, the integrator reach set is not rigidly convex, see[41, Sec. 3.1 and 3.2]. Fig. 3 helps visualize this for m = 1 .From Fig. 3a, we observe that a generic line for d = 2 has intersections with the bounding real algebraic curves whereasfrom (24), we know that p upper , p lower are degree polynomials.Likewise, Fig. 3b reveals that a generic line for d = 3 has intersections with the bounding real algebraic surfaces whereasfrom (25), we know that the polynomials p upper , p lower in thiscase, are of degree .Could the integrator reach set be spectrahedral shadow?Some calculations show that sufficient conditions as in [42] donot seem to hold. However, this remains far from conclusive.We summarize our taxonomy results in Fig. 4; the highlightedregion shows where the integrator reach set belongs. To answerwhether this highlighted region can be further narrowed down,seems significantly more challenging.IV. S IZE
We next quantify the “size” of the reach set R ( { x } , t ) by computing two functionals: its d -dimensional volume (Sec.IV-A), and its diameter or maximum width (Sec. IV-B). InSec. IV-C, we discuss how these functionals scale with thestate dimension d . Fig. 5: The integrator reach set R ( { x } , t ) with d = 3 , m = 2 , r = (2 , (cid:62) , x = (1 , , (cid:62) , [ α , β ] = [ − , , [ α , β ] = [ − , at t = 4 . A. Volume
The following result gives the volume formula for theintegrator reach set.
Theorem 5.
Fix x ∈ R d , let X ≡ { x } and compact U ⊂ R m . Consider the integrator dynamics (1)-(2) with d states, m inputs, and relative degree vector r = ( r , r , ..., r m ) (cid:62) . Define µ , . . . , µ m as in (12)-(13). Then the d -dimensional volume ofthe integrator reach set (3) at time t > is vol ( R ( { x } , t )) = 2 d m (cid:89) j =1 (cid:26) µ r j j t r j ( r j +1) / r j − (cid:89) k =1 k !(2 k + 1)! (cid:27) . (36)For a simple illustration of Theorem 5, consider d = 3 , m =2 with r = (2 , (cid:62) . The corresponding reach set R ( { x } , t ) at t = 4 is shown in Fig. 5 for x = (1 , , (cid:62) , U = [ − × × [ − , . Here µ = 5 and µ = 3 .This reach set, being a direct product of the double integra-tor reach set R (cf. Fig. 2) and the single integrator reachset R = { x (3) } (cid:117) [ − µ t, µ t ] , is a cylinder † . In [4], weexplicitly derived that vol ( R ) = µ t , and therefore, thevolume of this cylindrical set must be equal toheight of the cylinder × cross sectional area = 2 µ t × µ t = 43 µ µ t . Indeed, a direct application of the formula (36) recovers theabove expression.
Remark 2.
If the initial set X is not singleton, then comput-ing the volume of the reach set (4) requires us to compute thevolume of a Minkowski sum. Notice that vol (exp( t A ) X ) = | det (exp( t A )) | vol ( X )= exp (trace ( t A )) vol ( X ) † Here, the notation x (3) stands for the third component of vector x . exp m (cid:88) j =1 trace ( t A j ) vol ( X )= vol ( X ) , since from (2b), trace( A j ) = 0 for all j = 1 , . . . , m . There-fore, combining (4), (36) with the classical Brunn-Minkowskiinequality, we have a bound (vol ( R ( X , t ))) /d ≥ (vol ( X )) /d + 2 m (cid:89) j =1 (cid:26) µ r j j t r j ( r j +1) / r j − (cid:89) k =1 k !(2 k + 1)! (cid:27) /d . The above bound holds for any compact X ⊂ R d , notnecessarily convex.B. Diameter We now focus on another measure of the “size” of theintegrator reach set, namely its diameter, or maximal width.By definition, the width [12, p. 42] of the reach set R ( X , t ) , is w R ( X ,t ) ( η ) := h R ( X ,t ) ( η ) + h R ( X ,t ) ( − η ) , (37)where η ∈ S d − (the unit sphere imbedded in R d ), and thesupport function h R ( X ,t ) ( · ) is given by (17). In other words,(37) gives the width of the reach set in the direction η .For singleton X ≡ { x } , combining (17) and (37), we have w R ( { x } ,t ) ( η ) = (cid:90) t (cid:26) |(cid:104) η , ξ ( s ) (cid:105)| + |(cid:104)− η , ξ ( s ) (cid:105)| (cid:27) d s = 2 (cid:90) t |(cid:104) η , ξ ( s ) (cid:105)| d s, (38)where the last equality follows from the fact that ξ ( s ) in (14)is component-wise nonnegative for all ≤ s ≤ t .The diameter of the reach set is its maximal width: diam ( R ( X , t )) := max η ∈ S d − w R ( X ,t ) ( η ) . (39)Notice that (38) is a convex function of η ; see e.g., [43, p.79]. Thus, computing (39) amounts to maximizing a convexfunction over the unit sphere. We next derive a closed formexpression for (39). Theorem 6.
Fix x ∈ R d , let X ≡ { x } and compact U ⊂ R m . Consider the integrator dynamics (1)-(2) with d states, m inputs, and relative degree vector r = ( r , r , ..., r m ) (cid:62) . Define µ , . . . , µ m as in (12)-(13). The diameter of the integratorreach set (3) at time t > is diam ( R ( { x } , t )) = 2 (cid:107) ζ ( t ) (cid:107) = 2 m (cid:88) j =1 µ j (cid:107) ζ j (cid:107) / , (40) wherein ζ ( t ) is defined as in Sec. II-2, and the i th componentof the subvector ζ j ( t ) ∈ R r j is (cid:90) t s ( r j − i ) ( r j − i )! d s = t r j − i +1 ( r j − i + 1)! , i = 1 , , . . . , r j . (41) To illustrate Theorem 6, consider the triple integrator with d = 3 and m = 1 . In this case, U = [ α, β ] , µ := ( β − α ) / ,and we can parameterize the unit vector η ∈ S as η ≡ sin θ cos φ sin θ sin φ cos θ , θ ∈ [0 , π ] , φ ∈ [0 , π ) . Thus (39) reduces to µ max θ ∈ [0 ,π ] φ ∈ [0 , π ) (cid:90) t | s (sin θ cos φ ) / s sin θ sin φ + cos θ | d s. Furthermore, ζ ( t ) = ( t / , t / , t ) (cid:62) , and we obtain η max = sin θ max sin φ max sin θ max cos φ max cos θ max = ± √ t + 9 t + 36 t t , where ± means that either all components are plus or allminus. Thus, the maximizing tuples ( φ max , θ max ) ∈ [0 , π ] × [0 , π ) are given by ( φ max , θ max )= (cid:40)(cid:0) arctan (3 /t ) , arccos (cid:0) / √ t + 9 t + 36 (cid:1)(cid:1) , (cid:0) π + arctan (3 /t ) , arccos (cid:0) − / √ t + 9 t + 36 (cid:1)(cid:1) . (42)Hence, the diameter of the triple integrator reach set at time t is equal to ( µt/ √ t + 9 t + 36 .Fig. 6 shows how the width of the integrator reach setfor d = 3 , m = 1 varies over ( φ, θ ) ∈ [0 , π ] × [0 , π ) ,which parameterize the unit sphere S . The location of themaximizers are given by (42), and are depicted in Fig. 6 viafilled black circle and filled black square.For a visualization of the width and diameter for the doubleintegrator, see [4, Fig. 2]. C. Scaling Laws
We now turn to investigate how the volume and the diameterof the integrator reach set (3) scale with time and the statedimension. For maximal clarity, we focus on the single inputcase.
1) Scaling of the volume:
Fig. 7 plots the volume (36) forthe single input ( m = 1 ) case against time t for varying statespace dimension d . In this case, U = [ α, β ] , and therefore µ := ( β − α ) / . As expected, the volume of the reach setincreases with time for any fixed d .Let us now focus on the scaling of the volume with respectto the state dimension d . For m = 1 , using the knownasymptotic [44] for (cid:81) d − k =1 (2 k + 1)! /k ! , we find the d → ∞ asymptotic for the volume: vol ( R d ( { x } , t )) ∼ (2 µ ) d t d ( d +1) / exp (cid:0) d + (cid:1) c × d − ) d ( d + ) , where c ≈ . . . . is the Glaisher-Kinkelin constant [45,Sec. 2.15].Fig. 7 shows that when t is small, the volume of the largerdimensional reach set stays lower than its smaller dimensionalcounterpart. In particular, given two state space dimensionsig. 6: The width (38) for the single input triple integrator reach set R ( { x } , t ) is shown as a function of ( φ, θ ) ∈ [0 , π ] × [0 , π ) , whichparameterize the unit sphere S . Here U = [ − , and hence µ = 1 .The darker (resp. lighter) hues correspond to the higher (resp. lower)widths. The filled black circle and the filled black square correspondto the maximizers ( φ max , θ max ) given by (42). d, δ with d > δ , and all other parameters kept fixed, thereexists a critical time t cr when the volume of the d dimensionalreach set overtakes that of the δ dimensional reach set.For any d > δ , the critical time t cr satisfies vol ( R d ( { x } , t cr )) (cid:124) (cid:123)(cid:122) (cid:125) d dimensional volume = vol ( R δ ( { x } , t cr )) (cid:124) (cid:123)(cid:122) (cid:125) δ dimensional volume , which together with (36) yields t cr = (2 µ ) − d + δ +1 (cid:32) d − (cid:89) k = δ (2 k + 1)! k ! (cid:33) d − δ )( d + δ +1) . (43)In particular, for δ = d − , we get t cr = (cid:18) µ (2 d − d − (cid:19) /d , d = 2 , , . . . . (44)For instance, when µ = 1 , d = 3 , δ = 2 , we have t cr =(30) / ≈ . . When µ = 1 , d = 4 , δ = 3 , we have t cr = (420) / ≈ . . The dashed vertical lines in Fig. 7show the critical times given by (44).Applying Stirling’s approximation n ! ∼ √ πn ( n/e ) n , weobtain the d → ∞ asymptotic for (44): t cr ∼ e d µ − d − d , where ∼ denotes asymptotic equivalence [46, Ch. 1.4], and e is the Euler number.
2) Scaling of the diameter:
Fig. 8 plots the diameter of(40) for the single input ( m = 1 ) case against time t forvarying state space dimension d . As earlier, U = [ α, β ] , µ :=( β − α ) / . As expected, the diameter of the reach set increaseswith time for any fixed d . Fig. 7: For single input ( m = 1 ), the volume of the integrator reachset R ( { x } , t ) computed from (36) is plotted against time t forstate dimensions d = 2 , , . . . , with U = [ α, β ] = [ − , , µ :=( β − α ) / . The dashed vertical lines show the critical timesgiven by (44). As d → ∞ , the diameter approaches a limiting curve shownby the dotted line in Fig. 8. To derive this limiting curve, noticethat for m = 1 , the formula (40) gives lim d →∞ diam ( R ( { x } , t )) = lim d →∞ µ (cid:118)(cid:117)(cid:117)(cid:116) d (cid:88) j =1 (cid:18) t j j ! (cid:19) . (45)We write the partial sum d (cid:88) j =1 (cid:18) t j j ! (cid:19) = ∞ (cid:88) j =1 (cid:18) t j j ! (cid:19) (cid:124) (cid:123)(cid:122) (cid:125) =: S − ∞ (cid:88) j = d +1 (cid:18) t j j ! (cid:19) (cid:124) (cid:123)(cid:122) (cid:125) =: S , (46)and by ratio test, note that both the sums S , S converge.In particular, S converges to I (2 t ) − , where I ( · ) is thezeroth order modified Bessel function of the first kind. Thisfollows from the very definition of the ν th order modifiedBessel function of the first kind, given by I ν ( z ) := ( z/ ν ∞ (cid:88) j =0 (cid:0) z / (cid:1) j j ! Γ ( ν + j + 1) , ν ∈ R , where Γ( · ) denotes the Gamma function.On the other hand, using the definition of the generalizedhypergeometric function ‡ F ( a ; b , b ; z ) := ∞ (cid:88) n =0 ( a ) n ( b ) n ( b ) n z n n ! , we find that S = t d +1) 1 F (cid:0) d + 2 , d + 2; t (cid:1) (( d + 1)!) . ‡ Here, ( · ) n denotes the Pochhammer symbol [47, p. 256] or rising factorial. ig. 8: For single input ( m = 1 ), the diameter of the integrator reachset R ( { x } , t ) computed from (40) is plotted against time t for statedimensions d = 2 , , . . . , with U = [ α, β ] = [ − , , µ := ( β − α ) / . As d → ∞ , the diameter converges to µ (cid:112) I (2 t ) − ,shown by the dotted line. Therefore, (46) evaluates to S − S = I (2 t ) − − t d +1) 1 F (cid:0) d + 2 , d + 2; t (cid:1) (( d + 1)!) . (47)Combining (45), (46), (47), and using the continuity of thesquare root function on [0 , ∞ ) , we deduce that lim d →∞ diam ( R ( { x } , t )) = 2 µ (cid:113) lim d →∞ ( S − S )= 2 µ (cid:112) I (2 t ) − . (48)That lim d →∞ S exists and equals to zero, follows from (46)and the continuity of the square: lim d →∞ S = lim j →∞ (cid:18) t j j ! (cid:19) = (cid:18) lim j →∞ t j j ! (cid:19) = 0 . To see the last equality, let a j := t j /j ! . By the ratio test, lim sup j →∞ | a j +1 /a j | = lim j →∞ t/j = 0 < , hence { a j } is a Cauchysequence and lim j →∞ a j = 0 .The dotted line in Fig. 8 is the curve (48).V. B ENCHMARKING O VER - APPROXIMATIONS OF I NTEGRATOR R EACH S ETS
In practice, a standard approach for safety and performanceverification is to compute “tight” over-approximation of thereach sets of the underlying controlled dynamical system.Several numerical toolboxes such as [2], [3] are availablewhich over-approximate the reach sets using simple geometricshapes such as zonotopes and ellipsoids. Depending on theinterpretation of the qualifier “tight”, different optimizationproblems ensue, e.g., minimum volume outer-approximation[48]–[55].One potential application of our results in Sec. IV is tohelp quantify the conservatism in different over-approximation . 1
Fig. 9: (Top)
Zonotopic over-approximations of the double integratorreach sets; (bottom) the ratio of the volume of the single inputintegrator reach set R ( t ) and that of its zonotopic over-approximation R approx ( t ) for d = 2 , , , plotted against time t ∈ [0 , . The resultsare computed using the CORA toolbox with µ = 1 , X = { } . algorithms by taking the integrator reach set as a benchmarkcase. For instance, Fig. 9 shows the conservatism in zonotopicover-approximations R approx ( t ) of the single input integratorreach sets R ( { } , t ) ⊆ R approx ( { } , t ) for d = 2 , , with ≤ t ≤ and µ = 1 , computed using the CORAtoolbox [2], [56]. To quantify the conservatism, we used thevolume formula (36) for computing the ratio of the volumes vol( R ) / vol( R approx ) ∈ [0 , . The results shown in Fig. 9were obtained by setting the zonotope order in the CORAtoolbox, which means that the number of zonotopic segmentsused by CORA for over-approximation was ≤ d . As ex-pected, increasing the zonotope order improves the accuracyat the expense of computational speed, but among the differentdimensional volume ratio curves, trends similar to Fig. 9remain. It is possible [28, Thm. 1.1, 1.2] to compute the opti-mal zonotope order as function of the desired approximationaccuracy (i.e., desired Hausdorff distance from the zonoid).For the numerical results shown in Fig. 9, we found thediameters of the over-approximating zonotopes for d = 2 , , ,to be the same as that of the true diameters given by (40) forall times.Fig. 10 depicts the conservatism in ellipsoidal over-approximations R approx ( t ) of the single input integrator reachsets R ( { } , t ) ⊆ R approx ( { } , t ) for d = 2 , , with ≤ t ≤ and µ = 1 , following the algorithms in ellipsoidal toolbox[3]. Specifically, the reach set at time t , is over-approximatedby the intersection of a carefully constructed parameterizedfamily of ellipsoids E (cid:0) q ( t ) , Q (cid:96) i ( t ) ( t ) (cid:1) ig. 10: (Top) Ellipsoidal over-approximations of the double integrator reach sets; (bottom) the ratio of the volume ( left ) and diameter ( right )of the single input integrator reach set R ( t ) and that of its ellipsoidal over-approximation R approx ( t ) for d = 2 , , , plotted against time t ∈ [0 , . Two different ellipsoidal over-approximations are shown: one (in red) based on the S procedure, and the other (in blue) obtainedby scaling the maximum volume inner ellipsoid (MVIE) of the intersection of a parameterized family of ellipsoids. The results are computedfor µ = 1 , X = { } . := { x ∈ R d | ( x − q ( t )) (cid:0) Q (cid:96) i ( t ) ( t ) (cid:1) − ( x − q ( t )) (cid:62) ≤ } , for unit vectors (cid:96) i (0) ∈ R d , i = 1 , . . . , N . The choice of (cid:96) i (0) determines (cid:96) i ( t ) := exp( − A (cid:62) t ) (cid:96) i (0) , which in turnparameterizes the d × d symmetric positive definite shapematrix Q (cid:96) i ( t ) ( t ) ; we refer the readers to [57, Ch. 3.2], [34, Ch.3] where the corresponding evolution equations were derivedusing optimal control. The center vectors q ( t ) ∈ R d , andthe shape matrices Q (cid:96) i ( t ) ( t ) for this parameterized family ofellipsoids are constructed such that ∩ Ni =1 E (cid:0) q ( t ) , Q (cid:96) i ( t ) ( t ) (cid:1) isguaranteed to be a superset of the reach set at time t for anyfinite N , and for N → ∞ , recovers the reach set at that time.For the results shown in Fig. 10, we used N = 20 randomlychosen unit vectors (cid:96) i (0) ∈ R d . Ideally, one would liketo compute the (unique) minimum volume outer ellipsoid(MVOE), a.k.a. the L¨owner-John ellipsoid [58], [59] of theconvex set ∩ i =1 E (cid:0) q ( t ) , Q (cid:96) i ( t ) ( t ) (cid:1) , which is a semi-infiniteprogramming problem [43, Ch. 8.4.1], and has no known exact semidefinite programming (SDP) reformulation. We computedtwo different relaxations of this problem: one based on theS procedure [60, Ch. 3.7.2], and the other by homotheticscaling of the maximum volume inner ellipsoid (MVIE) [58,Thm. III] of the set ∩ i =1 E (cid:0) q ( t ) , Q (cid:96) i ( t ) ( t ) (cid:1) . Both of these lead to solving SDP problems, and both are guaranteed tocontain the L¨owner-John ellipsoid of the intersection of theparameterized family of ellipsoids. These suboptimal (w.r.t.the MVOE criterion) solutions, computed using cvx [61], areshown in Fig. 10.From Fig. 10, we observe that the S procedure resulted inless conservatism compared to the MVIE scaling, in terms ofvolume. While the volume ratio trends in Fig. 10 are similarto that observed in Fig. 9, the approximation quality are lower.In light of the results in Sec. III-A, this is not surprising: theintegrator reach sets being zonoids (i.e., Hausdorff limit ofzonotopes), the zonotopic outer-approximations are expectedto perform better than other over-approximating shape primi-tives.The main point here is that our results in Sec. IV providethe ground truth for the size of the integrator reach set,thereby help benchmarking the performance of reach setapproximation algorithms.VI. E PILOGUE
A. Recap
The present paper initiates a systematic study of integratorreach sets. We showed that this compact convex set is inact semialgebraic (Sec. III-B) as well as a zonoid (range ofan atom free vector measure) up to translation (Sec. III-A).We derived the equation of its boundary in both parametric(Proposition 1) and implicit form (Sec. III-C). We obtained theclosed form formula for the volume (Sec. IV-A) and diameter(Sec. IV-B) of these reach sets. We also derived the scalinglaws (Sec. IV-C) for these quantities. We pointed out that theseresults may be used to benchmark the performance of set over-approximation algorithms (Sec. V).
B. What Next
In the sequel Part II, we will show how the ideas presentedherein enable computing the reach sets for feedback lineariz-able systems. The focus will be in computing the reach setin transformed state coordinates associated with the normalform, and to map that set back to original state coordinatesunder diffeomorphism. This, however, requires nontrivial ex-tension of the basic theory presented here (especially thosein Proposition 1 and Sec III-C) since we will need to handletime-dependent set-valued uncertainty in transformed controlinput even when the original control takes values from a setthat is not time-varying. We will also address how to handlethe state constraints. Several numerical examples will be givento illustrate the results. A
PPENDIX
A. Proof of Theorem 1
Since support function is distributive over sum, we have h R j ( X j ,t ) ( y j ) = sup x j ∈X j (cid:104) y j , exp ( t A j ) x j (cid:105) + h (cid:82) t exp( s A j ) b j [ α j ,β j ]d s ( y j ) . (49)From the definition of support function, we have h b j [ α j ,β j ] ( y j ) = sup u j ∈ [ α j ,β j ] (cid:104) y j , b j u j (cid:105) . (50)The optimizer u opt j in (50) can be written in terms of theHeaviside unit step function H ( · ) as u opt j = α j + ( β j − α j ) H ( (cid:104) y j , b j (cid:105) )= α j + ( β j − α j ) ×
12 (1 + sgn ( (cid:104) y j , b j (cid:105) )) , where sgn( · ) denotes the signum function. Therefore, h b j [ α j ,β j ] ( y j ) = u opt j (cid:104) y j , b j (cid:105) = ν j (cid:104) y j , b j (cid:105) + µ j |(cid:104) y j , b j (cid:105)| . (51)The last equality used (13), and that | x | = x sgn( x ) for anyreal x .Using the property (8), from (51) we get h exp( s A j ) b j [ α j ,β j ] ( y j ) = h b j [ α j ,β j ] (exp (cid:0) s A (cid:62) j (cid:1) y j )= ν j (cid:104) y j , ξ j ( s ) (cid:105) + µ j |(cid:104) y j , ξ j ( s ) (cid:105)| . (52)Applying [4, Proposition 1] to (52), we obtain h (cid:82) t exp( s A j ) b j [ α j ,β j ]d s ( y j ) = (cid:90) t h exp( s A j ) b j [ α j ,β j ] ( y j )d s = (cid:28) y j , ν j ζ j ( t ) (cid:29) + µ j (cid:90) t |(cid:104) y j , ξ j ( s ) (cid:105)| d s, which upon substituting back in (49), and then using (11),completes the proof. (cid:4) B. Proof of Theorem 2
For s ∈ [0 , t ] , let the vector measure (cid:101) µ be definedas d (cid:101) µ ( s ) := ξ ( s )d s where ξ ( s ) is given by (14). Then (cid:82) t |(cid:104) y , ξ ( s ) (cid:105)| d s is exactly in the form of a support functionof a zonoid (see e.g., [25, Sec. 2]). Using the one-to-onecorrespondence between a compact convex set and its supportfunction, the corresponding set is a zonoid.From (8) and (17), notice that R ( { x } , t ) is the translationof a set with support function (cid:82) t |(cid:104) y , ξ ( s ) (cid:105)| d s , i.e., the trans-lation of a zonoid. Thus, we conclude that R ( { x } , t ) is azonoid. (cid:4) C. Proof of Proposition 1
From Sec. II-2, the supporting hyperplane at any x bdy ∈ ∂ R ( { x } , t ) with outward normal y ∈ R d is (cid:104) y , x bdy (cid:105) = h R ( { x } ,t ) ( y ) , and the Legendre-Fenchel conjugate h ∗R ( { x } ,t ) (cid:0) x bdy (cid:1) = 0 . (53)For j = 1 , . . . , m , let y comprise of subvectors y j ∈ R r j . Since the Cartesian product (10) is equivalent to theMinkowski sum R (cid:117) R (cid:117) . . . (cid:117) R m , and the support functionof Minkowski sum is the sum of support functions of thesummand sets [12, p. 48], we have h R ( { x } ,t ) ( y ) = m (cid:88) j =1 h R j ( { x } ,t ) ( y j ) ⇒ h ∗R ( { x } ,t ) (cid:0) x bdy (cid:1) = m (cid:88) j =1 h ∗R j ( { x } ,t ) (cid:16) x bdy j (cid:17) , (54)wherein the last line follows from the property that theLegendre-Fenchel conjugate of a separable sum equals to thesum of the Legendre-Fenchel conjugates [43, p. 95].Therefore, combining (53) and (54), we obtain m (cid:88) j =1 inf y j ∈ R rj (cid:26) (cid:104)− x bdy j + exp ( t A j ) x j + ν j ζ j ( t ) , y j (cid:105) + µ j (cid:90) t |(cid:104) y j , ξ j ( s ) (cid:105)| d s (cid:27) = 0 . (55)For j = 1 , . . . , m , since each objective in (55) involves anintegral of the absolute value of a polynomial in s of degree r j − , that polynomial can have at most r j − roots in theinterval [0 , t ] , i.e., can have at most r j − sign changes in thatinterval. If all r j − roots of the aforesaid polynomial are in [0 , t ] , we denote these roots as s ≤ s ≤ . . . ≤ s r j − , andwrite (cid:90) t |(cid:104) y j , ξ j ( s ) (cid:105)| d s = ± (cid:90) s (cid:104) y j , ξ j ( s ) (cid:105) d s ∓ (cid:90) s s (cid:104) y j , ξ j ( s ) (cid:105) d s ± . . . ± ( − r j − (cid:90) ts rj − (cid:104) y j , ξ j ( s ) (cid:105) d s (cid:104) y j , ± ζ j (0 , s ) ∓ ζ j ( s , s ) ± . . . ± ( − r j − ζ j ( s r j − , t ) (cid:105) . (56)Notice that even if the number of roots in [0 , t ] is strictlyless than § r j − , the expression (56) is generic in the sensethe corresponding summand integrals become zero. Thus,combining (55) and (56), we arrive at m (cid:88) j =1 inf y j ∈ R rj (cid:104)− x bdy j + exp ( t A j ) x j + ν j ζ j ( t ) ± µ j ζ j (0 , s ) ∓ µ j ζ j ( s , s ) ± . . . ± µ j ( − r j − ζ j ( s r j − , t ) , y j (cid:105) = 0 . (57)The left hand side of (57) being the sum of the infimumvalues of linear functions, can achieve zero if and only if eachof those infimum equals to zero, i.e., if and only if x bdy j = exp ( t A j ) x j + ν j ζ j ( t ) ± µ j ζ j (0 , s ) ∓ µ j ζ j ( s , s ) ± . . . ± ( − r j − µ j ζ j ( s r j − , t ) . (58)Using (6), (14) and (15), we simplify (58) to (21), therebycompleting the proof. (cid:4) D. Proof of Corollary 3
From (21), we get two different parametric representationsof x bdy j in terms of ( s , s , . . . , s r j − ) . One parametric repre-sentation results from the choice of positive sign for the ± ap-pearing in (21), and another for the choice of negative sign forthe same. Denoting the implicit representation correspondingto the parametric representation (21) with plus (resp. minus)sign as p upper j ( x ) = 0 (resp. p lower j ( x ) = 0 ), the statementfollows. (cid:4) E. Proof of Theorem 4
We notice that (21) gives polynomial parameterizations ofthe components of x bdy j for all j = 1 , . . . , m . In particular, foreach k ∈ { , . . . , r j } , the right hand side of (21) is a homoge-neous polynomial in r j − parameters ( s , s , . . . , s r j − ) ofdegree r j − k + 1 . By polynomial implicitization [24, p. 134],the corresponding implicit equations p upper j (cid:16) x bdy j (cid:17) = 0 (whenfixing plus sign for ± in (21)) and p lower j (cid:16) x bdy j (cid:17) = 0 (whenfixing minus sign for ± in (21)), must define affine varieties V R [ x ,...,x rj ] ( p upper j ) , V R [ x ,...,x rj ] ( p lower j ) in R [ x , . . . , x d ] .Specifically, denote the right hand sides of (21) as g ± , . . . , g ± r j for all j = 1 , . . . , m , where the superscriptsindicate that either all g ’s are chosen with plus signs, or allwith minus signs. Then write (21) as x bdy j (1) = g ± (cid:0) s , s , . . . , s r j − (cid:1) , x bdy j (2) = g ± (cid:0) s , s , . . . , s r j − (cid:1) , ... x bdy j ( r j ) = g ± r j (cid:0) s , s , . . . , s r j − (cid:1) . § this may happen either because there are repeated roots in [0 , t ] , orbecause some real roots exist outside [0 , t ] , or because some roots are complexconjugates, or a combination of the previous three. Now for each j = 1 , . . . , m , consider the ideal I ± j := (cid:104)(cid:104) x bdy j (1) − g ± , x bdy j (2) − g ± , . . . , x bdy j ( r j ) − g ± r j (cid:105)(cid:105)⊆ R [ s , s , . . . , s r j − , x , x , . . . , x r j ] , and let I ± j,r j − := I ± j ∩ R [ x , ..., x r j ] be the ( r j − thelimination ideal of I ± j . Then for each j = 1 , . . . , m , thevariety V (cid:16) I + j,r j − (cid:17) = V R [ x ,...,x rj ] ( p upper j ) . Likewise, the variety V (cid:16) I − j,r j − (cid:17) = V R [ x ,...,x rj ] ( p lower j ) .Thus, the algebraic boundary (i.e., the Zariski closure of theEuclidean boundary) of R j is ∂ R j = V R [ x ,...,x rj ] (cid:0) p upper j (cid:1) ∪ V R [ x ,...,x rj ] (cid:0) p lower j (cid:1) . Therefore, R j := { x ∈ R r j | p upper j ( x ) ≤ , p lower j ( x ) ≤ } issemialgebraic for all j = 1 , . . . , m .Since the Cartesian product of semialgebraic sets is semi-algebraic, the statement follows from (10). (cid:4) F. Proof of Theorem 5
We organize the proof in three steps.Step 1: From (10), we have vol ( R ( { x } , t )) = vol ( R × R × . . . × R m )= m (cid:89) j =1 vol ( R j ( { x } , t )) . (59)Step 2: Motivated by (59), we focus on deriving the r j -dimensional volume of R j ( { x } , t ) . For this purpose, weproceed as in [4] by uniformly discretizing the interval [0 , t ] into n subintervals (cid:20) ( i − tn , itn (cid:19) , i = 1 , , . . . , n, with ( n + 1) breakpoints { t i } ni =0 , where t i := it/n for i =0 , , . . . , n .From (19), we then have vol ( R j ( { x } , t )) = vol (cid:32) lim n →∞ n (cid:88) i =0 tn exp ( t i A j ) b j [ − µ j , µ j ] (cid:33) = lim n →∞ (cid:18) tn (cid:19) r j vol (cid:32) n (cid:88) i =0 exp ( t i A j ) b j [ − µ j , µ j ] (cid:33) = t r j lim n →∞ n r j vol (cid:32) n (cid:88) i =0 µ j ξ j ( t i )[ − , (cid:33) , (60)where ξ j was defined in (14).At this point, we recognize that the set n (cid:88) i =0 µ j ξ j ( t i )[ − , (61)in (60) is a Minkowski sum of n + 1 intervals, each intervalbeing rotated and scaled in R r j via different linear transfor-mations exp( t i A j ) , i = 0 , , . . . , n . In other words, the set(61) is a zonotope imbedded in R r j .sing the formula for the volume of zonotopes [14, eqn.(57)], [62, exercise 7.19], we can write (60) as vol ( R j ( { x } , t )) = (2 µ j t ) r j lim n →∞ n r j × (cid:88) ≤ i σ ( j ) } , and stands for “the number of”. We will now prove that c ( r j ) = r j − (cid:89) k =1 ( k !) (2 k + 1)! . (68)To this end, we write r j ! · c ( r j ) as an integral over the unitcube [0 , r j : r j ! · c ( r j ) = (cid:90) [0 , rj (cid:89) ≤ a
Verificationof Digital and Hybrid Systems . Springer, 2000, pp. 323–331.[2] M. Althoff, “An introduction to CORA 2015,” in
Proc. of the Workshopon Applied Verification for Continuous and Hybrid Systems , 2015.[3] A. A. Kurzhanskiy and P. Varaiya, “Ellipsoidal toolbox (ET),” in
Proceedings of the 45th IEEE Conference on Decision and Control .IEEE, 2006, pp. 1498–1503.[4] S. Haddad and A. Halder, “The convex geometry of integrator reachsets,” in . IEEE, 2020, pp.4466–4471.[5] R. J. Aumann, “Integrals of set-valued functions,”
Journal of mathemat-ical analysis and applications , vol. 12, no. 1, pp. 1–12, 1965.[6] J. M. Borwein and R. E. Crandall, “Closed forms: what they are andwhy we care,”
Notices of the AMS , vol. 60, no. 1, pp. 50–65, 2013.[7] J.-B. Hiriart-Urruty and C. Lemar´echal,
Convex analysis and minimiza-tion algorithms I: Fundamentals . Springer science & business media,2013, vol. 305.[8] A. A. Liapounoff, “Sur les fonctions-vecteurs completement additives,”
Izvestiya Rossiiskoi Akademii Nauk. Seriya Matematicheskaya , vol. 4,no. 6, pp. 465–478, 1940.[9] P. R. Halmos, “The range of a vector measure,”
Bulletin of the AmericanMathematical Society , vol. 54, no. 4, pp. 416–421, 1948.[10] J. Diestel and B. Faires, “On vector measures,”
Transactions of theAmerican Mathematical Society , vol. 198, pp. 253–271, 1974.[11] N. Dinculeanu,
Vector measures . Elsevier, 2014.[12] R. Schneider,
Convex bodies: the Brunn–Minkowski theory . Cambridgeuniversity press, 2014, no. 151.[13] P. McMullen, “On zonotopes,”
Transactions of the American Mathemat-ical Society , vol. 159, pp. 91–109, 1971.[14] G. Shephard, “Combinatorial properties of associated zonotopes,”
Cana-dian Journal of Mathematics , vol. 26, no. 2, pp. 302–321, 1974.[15] H. S. M. Coxeter,
Regular polytopes . Courier Corporation, 1973.[16] A. Girard, “Reachability of uncertain linear systems using zonotopes,” in
International Workshop on Hybrid Systems: Computation and Control .Springer, 2005, pp. 291–305. [17] A. Girard and C. Le Guernic, “Zonotope/hyperplane intersection forhybrid systems reachability analysis,” in
International Workshop onHybrid Systems: Computation and Control . Springer, 2008, pp. 215–228.[18] M. Althoff, O. Stursberg, and M. Buss, “Computing reachable setsof hybrid systems using a combination of zonotopes and polytopes,”
Nonlinear analysis: hybrid systems , vol. 4, no. 2, pp. 233–249, 2010.[19] M. Althoff and B. H. Krogh, “Zonotope bundles for the efficientcomputation of reachable sets,” in . IEEE, 2011,pp. 6814–6821.[20] J. K. Scott, D. M. Raimondo, G. R. Marseglia, and R. D. Braatz,“Constrained zonotopes: A new tool for set-based estimation and faultdetection,”
Automatica , vol. 69, pp. 126–136, 2016.[21] A. S. Adimoolam and T. Dang, “Using complex zonotopes for stabilityverification,” in . IEEE, 2016,pp. 4269–4274.[22] M. Althoff, “Reachability analysis of nonlinear systems using conserva-tive polynomialization and non-convex sets,” in
Proceedings of the 16thinternational conference on Hybrid systems: computation and control ,2013, pp. 173–182.[23] N. Kochdumper and M. Althoff, “Sparse polynomial zonotopes: Anovel set representation for reachability analysis,”
IEEE Transactionson Automatic Control , 2020.[24] D. Cox, J. Little, and D. OShea,
Ideals, varieties, and algorithms:an introduction to computational algebraic geometry and commutativealgebra . Springer Science & Business Media, 2013.[25] E. D. Bolker, “A class of convex bodies,”
Transactions of the AmericanMathematical Society , vol. 145, pp. 323–345, 1969.[26] R. Schneider and W. Weil, “Zonoids and related topics,” in
Convexityand its Applications . Springer, 1983, pp. 296–317.[27] P. Goodey and W. Wolfgang, “Zonoids and generalisations,” in
Hand-book of convex geometry . Elsevier, 1993, pp. 1297–1326.[28] J. Bourgain, J. Lindenstrauss, and V. Milman, “Approximation ofzonoids by zonotopes,”
Acta mathematica , vol. 162, no. 1, pp. 73–141,1989.[29] C. Vinzant, “The geometry of spectrahedra,” in
Sum of Squares: Theoryand Applications . American Mathematical Society, 2020, pp. 11–36.[30] A. Tarski, “A decision method for elementary algebra and geome-try,” in
Quantifier elimination and cylindrical algebraic decomposition .Springer, 1998, pp. 24–84.[31] A. Seidenberg, “A new decision method for elementary algebra,”
Annalsof Mathematics , pp. 365–374, 1954.[32] J. Bochnak, M. Coste, and M.-F. Roy,
Real algebraic geometry .Springer Science & Business Media, 2013, vol. 36.[33] G. Blekherman, P. A. Parrilo, and R. R. Thomas,
Semidefinite optimiza-tion and convex algebraic geometry . SIAM, 2012.[34] A. B. Kurzhanski and P. Varaiya,
Dynamics and Control of TrajectoryTubes: Theory and Computation . Springer, 2014, vol. 85.[35] G. Zaimi, “A polynomial implicitization,” MathOverflow. [Online].Available: https://mathoverflow.net/q/381335[36] H. S. Wilf, generatingfunctionology . CRC press, 2005.[37] L. Kronecker,
Zur Theorie der Elimination einer Variabeln aus zweialgebraischen Gleichungen . Buchdruckerei der K¨onigl. Akademie derWissenschaften (G. Vogt), 1881.[38] R. Salem,
Algebraic numbers and Fourier analysis . WadsworthPublishing Company, 1983.[39] R. Schneider, “Zonoids whose polars are zonoids,”
Proceedings of theAmerican Mathematical Society , vol. 50, no. 1, pp. 365–368, 1975.[40] Y. Lonke, “On zonoids whose polars are zonoids,”
Israel Journal ofMathematics , vol. 102, no. 1, pp. 1–12, 1997.[41] J. W. Helton and V. Vinnikov, “Linear matrix inequality representationof sets,”
Communications on Pure and Applied Mathematics , vol. 60,no. 5, pp. 654–674, 2007.[42] J. W. Helton and J. Nie, “Semidefinite representation of convex sets,”
Mathematical Programming , vol. 122, no. 1, pp. 21–64, 2010.[43] S. P. Boyd and L. Vandenberghe,
Convex optimization . Cambridgeuniversity press, 2004.[44] OEIS Foundation Inc., “The On-Line Encyclopedia of Integer Se-quences,” http://oeis.org/A107254, 2019.[45] S. R. Finch,
Mathematical constants . Cambridge university press, 2003.[46] N. G. De Bruijn,
Asymptotic methods in analysis . Courier Corporation,1981, vol. 4.[47] M. Abramowitz and I. A. Stegun,
Handbook of mathematical functionswith formulas, graphs, and mathematical tables . US Governmentprinting office, 1970, vol. 55.[48] F. C. Schweppe,
Uncertain dynamic systems . Prentice Hall, 1973.49] F. Chernousko, “Optimal guaranteed estimates of indeterminacies withthe aid of ellipsoids, I,”
Engineering Cybernetics , vol. 18, no. 3, pp.1–9, 1980.[50] ——, “Guaranteed ellipsoidal estimates of uncertainties in controlproblems,”
IFAC Proceedings Volumes , vol. 14, no. 2, pp. 869–874,1981.[51] D. Maksarov and J. Norton, “State bounding with ellipsoidal set de-scription of the uncertainty,”
International Journal of Control , vol. 65,no. 5, pp. 847–866, 1996.[52] C. Durieu, E. Walter, and B. Polyak, “Multi-input multi-output ellip-soidal state bounding,”
Journal of optimization theory and applications ,vol. 111, no. 2, pp. 273–303, 2001.[53] T. Alamo, J. M. Bravo, and E. F. Camacho, “Guaranteed state estimationby zonotopes,”
Automatica , vol. 41, no. 6, pp. 1035–1043, 2005.[54] A. Halder, “On the parameterized computation of minimum volumeouter ellipsoid of Minkowski sum of ellipsoids,” in . IEEE, 2018, pp. 4040–4045.[55] ——, “Smallest ellipsoid containing p-sum of ellipsoids with applicationto reachability analysis,”
IEEE Transactions on Automatic Control , 2020.[56] “CORA: A tool for continuous reachability analysis,” https://tumcps.github.io/CORA/, accessed: 2021-02-14. [57] A. Kurzhanski˘ı and I. V´alyi,
Ellipsoidal calculus for estimation andcontrol . Nelson Thornes, 1997.[58] F. John, “Extremum problems with inequalities as subsidiary condi-tions,”
Studies and Essays: Courant Anniversary Volume, presented toR. Courant on his 60th Birthday , pp. 187–204, 1948.[59] M. Henk, “L¨owner-John ellipsoids,”
Documenta Math , vol. 95, p. 106,2012.[60] S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan,
Linear matrixinequalities in system and control theory . SIAM, 1994.[61] M. Grant and S. Boyd, “CVX: Matlab Software for Disciplined ConvexProgramming, version 2.1,” http://cvxr.com/cvx, Mar 2014.[62] G. M. Ziegler,
Lectures on polytopes . Springer Science & BusinessMedia, 2012, vol. 152.[63] R. A. Horn and C. R. Johnson,
Matrix analysis . Cambridge universitypress, 2012.[64] T. M. Apostol, “An elementary view of Euler’s summation formula,”
TheAmerican Mathematical Monthly , vol. 106, no. 5, pp. 409–418, 1999.[65] E. Hairer and G. Wanner,
Analysis by its History . Springer Science &Business Media, 2006, vol. 2.[66] N. De Bruijn, “On some multiple integrals involving determinants,”