The Curious Case of Integrator Reach Sets, Part I: Basic Theory

This is the first of a two part paper investigating the geometry of the integrator reach sets, and the applications thereof. In this Part I, assuming box-valued input uncertainties, we establish that this compact convex reach set is semialgebraic, translated zonoid, and not a spectrahedron. We derive the parametric as well as the implicit representation of the boundary of this reach set. We also deduce the closed-form formula for the volume and diameter of this set, and discuss their scaling with state dimension and time. We point out that these results may be utilized in benchmarking the performance of the reach set overapproximation algorithms.


I. INTRODUCTION
Integrators with bounded controls are ubiquitous in systemscontrol.They appear as Brunovsky normal forms for the feedback linearizable nonlinear systems.They also appear frequently as benchmark problems to demonstrate the performance of the reach set computation algorithms.Despite their prominence, specific results on the geometry of the integrator reach sets are not available in the systems-control literature.Broadly speaking, the existing results come in two flavors.On one hand, very generic statements are known, e.g., these reach sets are compact convex sets whenever the set of initial conditions is compact convex, and the controls take values from a compact (not necessarily convex) set [1].On the other hand, several numerical toolboxes [2], [3] are available for tight outer approximation of the reach sets over computationally benign geometric families such as ellipsoids and zonotopes.The lack of concrete geometric results imply the absence of ground truth when comparing the efficacy of different algorithms, and one has to content with graphical or statistical (e.g., Monte Carlo) comparisons.
Building on the preliminary results in [4], this paper undertakes a systematic study of the integrator reach sets.In particular, we answer the following basic questions: Q1. what kind of compact convex sets are these (Section IV)? Q2. how big are these sets (Section V)? Q3. how these results on the geometry of integrator reach sets can be applied in practice (Section VI)?We consider the integrator dynamics having d states and m inputs with relative degree vector r = (r 1 , r 2 , . . ., r m ) ∈ Z m + (vector of positive integers).In other words, we consider a Brunovsky normal form with m integrators where the jth integrator has degree r j for j ∈ [m] := {1, . . ., m}.The dynamics is In (2a), the symbol blkdiag (•) denotes a block diagonal matrix whose arguments constitute its diagonal blocks.In (2b), the notation 0 rj ×1 stands for the r j × 1 column vector of zeros, and e k denotes the kth basis (column) vector in R for k ≤ .The symbol (...|...|...) denotes horizontal concatenation.
Let R (X 0 , t) denote the forward reach set of (1) at time t > 0, starting from a given compact convex set of initial conditions X 0 ⊂ R d , i.e., R (X 0 , t) := measurable u(•)∈U ⊂R m x(t) ∈ R d | (1) and (2) hold, x(t = 0) ∈ X 0 compact convex, U compact .(3) In words, R (X 0 , t) is the set of all states that the controlled dynamics (1)-( 2) can reach at time t > 0, starting from the set X 0 at t = 0, with measurable control u(•) ∈ U compact.Formally, R (X 0 , t) = exp(tA)X 0 t 0 exp ((t − τ )A) BU dτ = exp(tA)X 0 t 0 exp (sA) BU ds, where denotes the Minkowski sum.The set-valued integral [5] in ( 4) is defined for any point-to-set function F (•), as where the summation symbol Σ denotes the Minkowski sum, and • is the floor operator; see e.g., [1].Our objective is to study the geometry of (4) in detail.This paper significantly expands our preliminary works [4], [6]: here we consider multi-input integrators as opposed to the single input case considered in [4].Even for the single input case, while [4,Thm. 1] derived an exact formula for the volume of the reach set, that formula involved limit and nested sums, and in that sense, was not really a closed-form formula [7] -certainly not amenable for numerical computation.In this paper, we derive closed-form formula for the general multi-input case when the input uncertainty is box-valued, i.e., U is hyperrectangle.In the same setting, the present paper addresses previously unexplored directions: the scaling laws for the volume and diameter of integrator reach sets, exact parametric and implicit equations for the boundary, and the classification of these sets.
The paper is structured as follows.After reviewing some preliminary concepts in Sec.II, we consider the integrator reach set resulting from box-valued input set uncertainty in Sec.III.The results on taxonomy and the boundary of the corresponding reach set are provided in Sec.IV.The results on the size of this set are collected in Sec.V.The application of these results for benchmarking the reach set over-approximation algorithms are discussed in Sec.VI.All proofs are deferred to the Appendix.Sec.VII summarizes the paper, and outlines the directions pursued in its sequel Part II.

II. PRELIMINARIES
In the following, we summarize some preliminaries which will be useful in the main body and in the Appendix.
1) State transition matrix: For 0 ≤ s < t, the state transition matrix Φ(t, s) associated with (2) is with each diagonal block is upper triangular.Specifically, the jth diagonal block of size r j × r j is written element-wise as where k is the row index, is the column index, and k, ∈ [r j ] for each j ∈ [m].The diagonal entries in (6) are unity.
2) Support function: where •, • denotes the standard Euclidean inner product.Geometrically, h K (y) gives the signed distance of the supporting hyperplane of K with outer normal vector y, measured from the origin.Furthermore, the supporting hyperplane at x bdy ∈ ∂K is y, x bdy = h K (y), and we can write The support function h K (y) is convex in y.For more details on the support function, we refer the readers to [8, Ch.V].
The support function h K (y) uniquely determines the set K. Given matrix-vector pair (Γ, γ) ∈ R d×d × R d , the support function of the affine transform ΓK + γ is Given a function f : R d → R∪{+∞}, its Legendre-Fenchel conjugate is From ( 7)- (10), it follows that h K (y) is the Legendre-Fenchel conjugate of the indicator function Since the indicator function of a convex set is a convex function, the biconjugate This will be useful in Section IV.
To proceed further, we introduce some notations.Since U is compact, let that is, α j and β j are the component-wise minimum and maximum, respectively, of the input vector.Furthermore, let and introduce for j ∈ [m].When t 0 = 0, we simplify the notations as Using (13) and following (7), we deduce Proposition 1 stated next (proof in Appendix A).
Proposition 1. (Support function for compact U) For compact convex X 0 ⊂ R d , and compact U ⊂ R m , the support function of the reach set (4) is where conv(•) denotes the convex hull.
3) Polar dual: The polar dual K • of any non-empty set K ⊂ R d is given by From this definition, it is immediate that K • contains the origin, and is a closed convex set.The bipolar (K • ) • = closure (conv (K ∪ {0})).Thus, if K is compact convex and contains the origin, then we have the involution From ( 7) and (17), notice that K • is the unit support function ball, i.e., In Sec.IV-D, we will mention some properties of the polar dual of the integrator reach set.
4) Vector measure: Let F be a σ-field of the subsets of a set.A countably additive mapping µ : F → R d is termed a vector measure.Here, "countably additive" means that for any sequence Some of the early investigations of vector measures were due to Liapounoff [9] and Halmos [10]; relatively recent references are [11], [12].
5) Zonotope: A zonotope Z ⊂ R d is a finite Minkowski sum of closed line segments or intervals {I i } n i=1 where these intervals are imbedded in the ambient Euclidean space R d .Explicitly, for some positive integer n, we write Thus, a zonotope is the range of an atomic vector measure.Alternatively, a zonotope can be viewed as the affine image of the unit cube.A compact convex polytope is a zonotope if and only if all its two dimensional faces are centrally symmetric [13, p. 182].For instance, the cross polytope {x ∈ R d | x 1 ≤ 1}, is not a zonotope.Standard references on zonotope include [14], [15], [16,Ch. 2.7].

III. BOX-VALUED INPUT UNCERTAINTY
In the remaining of this paper, we characterize the exact reach set (3) when input set U ⊂ R m is box-valued, and remark on the quality of approximation for the same when U is arbitrary compact.
When U ⊂ R m is box-valued, denote the reach set (3) as R , i.e., with a box superscript * .In this case, each of the m single input integrator dynamics with r j dimensional state subvectors for j ∈ [m], are decoupled from each other.Then R (X 0 , t) ⊂ R d is the Cartesian product of these single input integrator reach sets: R j (X 0 , t) ⊂ R rj for j ∈ [m], i.e., In what follows, we will sometimes exploit that (18) may also be written as † a Minkowski sum R 1 . . .R m .Notice that the decoupled dynamics also allows us to write a Minkowski sum decomposition for the set of initial conditions and accordingly, the initial condition subvectors x j0 ∈ X j0 ⊂ R rj for j ∈ [m].Thus x 0 = (x 10 , . . ., x m0 ) .Since the support function of the Minkowski sum is equal to the sum of the support functions, we have This leads to the following result (proof in Appendix B) which will come in handy in the ensuing development.
Theorem 1. (Support function for box-valued U) For compact convex X 0 ⊂ R d , and box-valued input uncertainty set given by the support function of the reach set ( 4) is The formula (21) upper bounds (16) resulting from the same initial condition and arbitrary compact U ⊂ R m with {α j , β j } m j=1 related to U via (11).Thus, from (8), the reach set R with box-valued input uncertainty will over-approximate the reach set R associated with arbitrary compact U, at any given t > 0, provided {α j , β j } m j=1 are defined as (11).When U is compact but not box-valued, then we can quantify the quality of the aforesaid over-approximation in terms of the two-sided Hausdorff distance metric dist between the convex compact sets R , R ⊂ R d , expressible [13,Thm. 1.8.11] Thanks to (8), the absolute value in (22) It is known that [26, Prop.6.1] the set R (X 0 , t) resulting from a linear time invariant dynamics such as (1)-(2) remains invariant under the closure of convexification of the input set U. Therefore, it is possible that R = R and dist = 0 even when the compact set U is nonconvex.For instance, the reach set R (X 0 , t) resulting from some compact convex X 0 ⊂ R d and dynamics (1)-( 2) with the nonconvex input uncertainty set {−1, 1} m , is identical to R (X 0 , t) resulting from the same X 0 , same dynamics, and the box-valued input uncertainty set (20) with Likewise, for the same compact convex If R (X 0 , t) results from the same X 0 , same dynamics, and input uncertainty set (20) with where q is the Hölder conjugate of max{1, p}, i.e., 1 max{1,p} + 1 q = 1, and 1 < q ≤ ∞.In this case, the positive value ( 23) quantifies the quality of strict over-approximation R p ⊂ R for 0 < p < ∞.The objective in (23) being positive homogeneous, admits lossless constraint convexification y 2 ≤ 1, and the corresponding maximal value ‡ for moderate dimensions d, can be found by direct numerical search.

IV. TAXONOMY AND BOUNDARY For
that the reach set R given by ( 3) is compact convex for all t > 0 provided U is compact.However, it is not immediate what kind of convex set R is, even for singleton X 0 ≡ {x 0 }.
In this Section, we examine the question "what type of compact convex set R ({x 0 }, t) is" when U is box-valued uncertainty set of the form (20).In the same setting, we also derive the equations for the boundary ∂R ({x 0 }, t).
Notice that for non-singleton X 0 , the taxonomy question is not well-posed since the classification then will depend on X 0 .Also, setting X 0 ≡ {x 0 } in (4), it is apparent that R ({x 0 }, t) is a translation of the set-valued integral in (4).Thus, classifying R ({x 0 }, t) amounts to classifying the second summand in (4).
A zonoid is a compact convex set that is defined as the range of an atom free vector measure (see Sec. II-4).Affine image of a zonoid is a zonoid.Minkowski sum of zonoids is also a zonoid.We refer the readers to [28]- [30], [31, Sec.I] for more details on the properties of a zonoid.By slight abuse of nomenclature, in this paper we use the term zonoid up to translation, i.e., we refer to the translation of zonoids as zonoids (instead of using another term such as "zonoidal translates").
Let us mention a few examples.Any compact convex symmetric set in R 2 is a zonoid.In dimensions three or more, all p norm balls for p ≥ 2 are zonoids.‡ As such, ( 23) has a difference of convex objective, and by the Weierstrass extreme value theorem, the maximum is achieved.
An alternative way to think about the zonoid is to view it as the limiting set (convergence with respect to the twosided Hausdorff distance, see e.g., [4, Appendix B]) of the Minkowski sum of line segments, i.e., the limit of a sequence of zonotopes [14], [15], [28].Formally, given a Hausdorff convergent sequence of zonotopes {Z j }, the zonoid Z ∞ is for some a ij ≤ b ij (element-wise vector inequality), and a suitable mapping n : Z + → Z + .Our analysis will make use of this viewpoint in Sec.V-A.Our main result in this subsection is the following.
Theorem 2. The reach set R given by (3) with X 0 ≡ {x 0 } and U given by (20), is a zonoid.
To appreciate Theorem 2 via the limiting viewpoint mentioned before, let us write where all summation symbols denote Minkowski sums.The first term in ( 24) denotes a translation.In the second term, the outer summation over index j arises by writing the Cartesian product (18) as the Minkowski sum R 1 . . .R m .Furthermore, uniformly discretizing [0, t] into n subintervals [(i − 1)t/n, it/n), i = 1, . . ., n, we write t 0 exp(sA j )b j [−µ j , µ j ]ds as the limit of the Minkowski sum over index i.Geometrically, the innermost summands in the second term denote non-uniformly rotated and scaled line intervals in R j .In other words, the second term in ( 24) is a Minkowski sum of m sets, each of these sets being the limit of a sequence of sets {Z n } comprising of zonotopes which are the Minkowski sum of n + 1 line segments.Since lim n→∞ Z n is a zonoid, the second term in ( 24) is a Minkowski sum of m zonoids, and is therefore a zonoid [28,Thm. 1.5].The entire right hand side of (24), then, is translation of a zonoid, and hence a zonoid.
In the following, we derive formulae for the boundary (Proposition 2 and Sec.IV-C) and volume (Theorem 5) of the integrator reach set with X 0 = {x 0 } (singleton set).From (25), it is clear that one cannot expect similar closed form formulae for arbitrary compact (or even arbitrary compact convex) X 0 .In this sense, our closed form formulae are as general as one might hope for.For a specific non-singleton X 0 , one can use these formulae to first derive the boundary (resp.volume) of R ({0}, t), and then use (25) to get numerical estimates for the boundary (resp.volume) of R (X 0 , t) (cf.Remark 2).Semialgebraic sets are closed under finitely many unions and intersections, complement, topological closure, polynomial mapping including projection [33], [34], and Cartesian product.For details on semialgebraic sets, we refer the readers to [35,Ch. 2]; see [36,Appendix A.4.4] for a short summary.
In Proposition 2 below, we derive a parametric representation of x bdy ∈ ∂R ({x 0 }, t), the boundary of the reach set.Then we use this representation to establish semialgebraicity of R ({x 0 }, t) in Theorem 4 that follows.
Proposition 2. For relative degree vector r = (r 1 , . . ., r m ) , and fixed x 0 ∈ R d comprising of subvectors x j0 ∈ R rj where j ∈ [m], consider the reach set (4) with singleton X 0 ≡ {x 0 } and U given by (20).For j ∈ [m], define µ 1 , . . ., µ m and ν 1 , . . ., ν m as in ( 11)- (12).Let the indicator function 1 k≤ := 1 for k ≤ , and := 0 otherwise.Then the components of admit parametric representation in terms of the parameters This parameterization is given by where x bdy j (k) denotes the kth component of the jth subvector The following is a consequence of the ± appearing in (26).
Likewise, in (28), taking plus (resp.minus) signs in each of component gives the parametric representation of the surface p upper 2 (x) = 0 (resp.p lower 2 = 0).The resulting set R 2 is the triple integrator reach set, and is shown in Fig. 1.Now we come to the main result of this subsection.
Let us illustrate the bounding curves and surfaces for ( 27) and ( 28) respectively, in the implicit form.Eliminating the parameter s 1 from (27) reveals that p upper 1 , p lower 1 are parabolas.In particular, and the formula for p lower 1 (x bdy 1 (1), x bdy 1 (2)) follows mutatis mutandis.
A natural question is whether one can generalize the implicitizations as in ( 29), (30) to arbitrary state dimensions.This is what we address next.

C. Implicitization of ∂R ({x 0 }, t)
To derive the implicit equations for the bounding algebraic hypersurfaces p upper j , p lower j ∈ R x 1 , . . ., x rj for all j ∈ [m], we need to eliminate the parameters s 1 , s 2 , . . ., s rj −1 from (26).For this purpose, it is helpful to write (26) succinctly as where To simplify the rather unpleasant notation ρ ± j,k , we will only address the m = 1 case.In (31), this allows us to replace r j by d, and to drop the subscript j from the ρ's.This does not invite any loss of generality in terms of implicitization since post derivation, we can replace d by r j to recover the respective p j 's.
With slight abuse of notation, we will also drop the superscript ± from the ρ's in (31).Recall that the plus (resp.minus) superscript in the ρ's indicates p upper j (resp.p lower j ).From (32), it is clear that in either case, the ρ j,k is affine in x bdy j (k), which is the kth coordinate of the boundary point for the jth block.Importantly, for k ∈ [r j ], the quantity ρ j,k does not depend on any other component of the boundary point than the kth component.Again, the plus-minus superscripts can be added back post implicitization.
In summary, (39) is the desired implicitization of the bounding hypersurfaces of the single input integrator reach set (up to the change of variables).The Cartesian product of these implicit hypersurfaces gives the implicitization in the multi input case.

D. Dual of R ({x 0 }, t)
From convex geometry standpoint, it is natural to ask what kind of characterization is possible for the polar dual (see Sec. II-3) of the integrator reach set R or R .We know in general that R • will be a closed convex set.Depending on the choice of x 0 , U and t, the set R ({x 0 }, t) may not contain the origin, and thus the bipolar that is, we do not have the involution in general.Furthermore, since R ({x 0 }, t) is semialgebraic from Sec. IV-B, so must be its polar dual R ({x 0 }, t) • ; see e.g., [36, Ch. 5, Sec.5.2.2].We also know from Sec. IV-A that R ({x 0 }, t) is a zonoid.However, the polar of a zonoid is not a zonoid in general [42], [43], and we should not expect R ({x 0 }, t) • to be one.

E. Summary of Taxonomy
So far we explained that the compact convex set R ({x 0 }, t) is semialgebraic, and a translated zonoid.Two well-known subclasses of convex semialgebraic sets are the spectrahedra and the spectrahedral shadows.The spectrahedra, a.k.a.linear matrix inequality (LMI) representable sets are affine slices of the symmetric positive semidefinite cone.The spectrahedral shadows, a.k.a.lifted LMI or semidefinite  representable sets are the projections of spectrahedra.The spectrahedral shadows subsume the class of spectrahedra; e.g., the set {(x 1 , x 2 ) ∈ R 2 | x 4 1 + x 4 2 ≤ 1} is a spectrahedral shadow but not a spectrahedron.The polar duals of spectrahedra are spectrahedral shadows [36, Ch. 5, Sec.5.5].
We note that R is not a spectrahedron.To see this, we resort to the contrapositive of [44,Thm. 3.1].Specifically, the number of intersections made by a generic line passing through an interior point of the d-dimensional reach set R with its real algebraic boundary is not equal to the degree of the bounding algebraic hypersurfaces, the latter we know from Sec. IV-C to be ).In other words, the R is not rigidly convex, see [44, 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 4 intersections with the bounding real algebraic curves whereas from (29), we know that p upper , p lower are degree 2 polynomials.Likewise, Fig. 3b reveals that a generic line for d = 3 has 6 intersections with the bounding real algebraic surfaces whereas from (30), we know that the polynomials p upper , p lower in this case, are of degree 4.
Could the reach set R be spectrahedral shadow?Some calculations show that sufficient conditions as in [45] do not seem to hold.However, this remains far from conclusive.We summarize our taxonomy results in Fig. 4; the highlighted region shows where the integrator reach set belongs.To answer whether this highlighted region can be further narrowed down, seems significantly more challenging.

V. SIZE
We next quantify the "size" of the reach set R ({x 0 }, t) by computing two functionals: its d-dimensional volume (Sec.V-A), and its diameter or maximum width (Sec.V-B).In Sec.V-C, we discuss how these functionals scale with the state dimension d.

A. Volume
The following result gives the volume formula for the integrator reach set R .Theorem 5. Fix x 0 ∈ R d , let X 0 ≡ {x 0 } and U given by (20).Consider the integrator dynamics (1)-( 2) with d states, m inputs, and relative degree vector r = (r 1 , r 2 , ..., r m ) .Define µ 1 , . . ., µ m as in ( 11)- (12).Then the d-dimensional volume of the integrator reach set (3) For a simple illustration of Theorem 5, consider d = 3, m = 2 with r = (2, 1) .The corresponding reach set R ({x 0 }, t) at t = 4 is shown in Fig. 5 for This reach set, being a direct product of the double integrator reach set R 1 (cf.Fig. 2) and the single integrator reach set R 2 = {x 0 (3)} [−µ 2 t, µ 2 t], is a cylinder ¶ .In [4], we explicitly derived that vol (R 1 ) = 2 3 µ 2 1 t 3 , and therefore, the volume of this cylindrical set must be equal to "height of the cylinder × cross sectional area", i.e., Indeed, a direct application of the formula (41) recovers the above expression.¶ Here, the notation x 0 (3) stands for the third component of vector x 0 .= vol (X 0 ) , since from (2b), trace(A j ) = 0 for all j = 1, . . ., m.Therefore, combining ( 4), ( 41) with the classical Brunn-Minkowski inequality, we obtain a lower bound for vol R as vol R (X 0 , t) The above bound holds for any compact X 0 ⊂ R d , not necessarily convex.

B. Diameter
We now focus on another measure of "size" for the integrator reach set R , namely its diameter, or maximal width.
By definition, the width [13, p. 42] of R (X 0 , t), is where η ∈ S d−1 (the unit sphere imbedded in R d ), and the support function h R (X0,t) (•) is given by (21).In other words, (42) gives the width of R in the direction η.
For singleton X 0 ≡ {x 0 }, combining ( 21) and ( 42), we have where the last equality follows from the fact that ξ(s) in ( 13) is component-wise nonnegative for all 0 ≤ s ≤ t.
The diameter of the reach set R is its maximal width: diam R (X 0 , t) := max Notice that ( 43) is a convex function of η; see e.g., [46, p. 79].Thus, computing (44) amounts to maximizing a convex function over the unit sphere.We next derive a closed form expression for (44).
For a visualization of the width and diameter for the double integrator, see [4, Fig. 2].

C. Scaling Laws
We now turn to investigate how the volume and the diameter of the integrator reach set scale with time and the state dimension.While scaling laws reveal limits of performance of engineered systems (see e.g., [47], [48]), they have not been formally investigated in the context of reach sets even though empirical studies are common [20], [49, Table 1], [50,Fig. 4].
For clarity, we focus on the single input case and hence do not notationally distinguish between R and R .
1) Scaling of the volume: Fig. 7 plots the volume (41) for the single input (m = 1) case against time t for varying state space dimension d.In this case, U = [α, β], and therefore µ := (β − α)/2.As expected, the volume of the reach set increases with time for any fixed d.
Let us now focus on the scaling of the volume with respect to the state dimension d.For m = 1, using the known The darker (resp.lighter) hues correspond to the higher (resp.lower) widths.The filled black circle and the filled black square correspond to the maximizers (φ max , θ max ) given by (47).asymptotic [51] for d−1 k=1 (2k + 1)!/k!, we find the d → ∞ asymptotic for the volume: , where c ≈ 1.2824 . . . is the Glaisher-Kinkelin constant [52, Sec.2.15].Fig. 7 shows that when t is small, the volume of the larger dimensional reach set stays lower than its smaller dimensional counterpart.In particular, given two state space dimensions d, d with d > d , and all other parameters kept fixed, there exists a critical time t cr when the volume of the d dimensional reach set overtakes that of the d dimensional reach set.
2) Scaling of the diameter: Fig. 8 plots the diameter of ( 45) for the single input (m = 1) case against time t for varying state space dimension d.As earlier, U = [α, β], µ := (β −α)/2.As expected, the diameter of the reach set increases with time for any fixed d.
As d → ∞, the diameter approaches a limiting curve shown by the dotted line in Fig. 8.To derive this limiting curve, notice that for m = 1, the formula (45) gives We write the partial sum and by ratio test, note that both the sums S 1 , S 2 converge.In particular, S 1 converges to I 0 (2t) − 1, where I 0 (•) is the zeroth order modified Bessel function of the first kind.This follows from the very definition of the νth order modified Bessel function of the first kind, given by where Γ(•) denotes the Gamma function.
On the other hand, using the definition of the generalized hypergeometric function ||  we find that Therefore, (51) evaluates to Combining ( 50), ( 51), (52), and using the continuity of the square root function on [0, ∞), we deduce that That lim d→∞ S 2 exists and equals to zero, follows from (51) and the continuity of the square: 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 < 1, hence {a j } is a Cauchy sequence and lim j→∞ a j = 0.
The dotted line in Fig. 8 is the curve (53).

VI. BENCHMARKING OVER-APPROXIMATIONS OF INTEGRATOR REACH SETS
In practice, a standard approach for safety and performance verification is to compute "tight" over-approximation of the reach sets of the underlying controlled dynamical system.Several numerical toolboxes such as [2], [3] are available which over-approximate the reach sets using simple geometric shapes such as zonotopes and ellipsoids.Depending on the interpretation of the qualifier "tight", different optimization problems ensue, e.g., minimum volume outer-approximation [55]- [62].
One potential application of our results in Sec.V is to help quantify the conservatism in different over-approximation algorithms by taking the integrator reach set as a benchmark case.For instance, Fig. 9 shows the conservatism in zonotopic over-approximations R approx (t) of the single input integrator reach sets R({0}, t) ⊆ R approx ({0}, t) for d = 2, 3, 4 with 0 ≤ t ≤ 1 and µ = 1, computed using the CORA toolbox [2], [63].To quantify the conservatism, we used the volume formula (41) for computing the ratio of the volumes vol(R)/vol(R approx ) ∈ [0, 1].The results shown in Fig. 9 were obtained by setting the zonotope order 50 in the CORA toolbox, which means that the number of zonotopic segments used by CORA for over-approximation was ≤ 50d.As expected, increasing the zonotope order improves the accuracy at the expense of computational speed, but among the different dimensional volume ratio curves, trends similar to Fig. 9 remain.It is possible [31,Thm. 1.1,1.2] to compute the optimal zonotope order as function of the desired approximation accuracy (i.e., desired Hausdorff distance from the zonoid).
For the numerical results shown in Fig. 9, we found the diameters of the over-approximating zonotopes for d = 2, 3, 4, to be the same as that of the true diameters given by ( 45) for all times.Fig. 10 depicts the conservatism in ellipsoidal overapproximations R approx (t) of the single input integrator reach sets R({0}, t) ⊆ R approx ({0}, t) for d = 2, 3, 4 with 0 ≤ t ≤ 1 and µ = 1, following the algorithms in ellipsoidal toolbox [3].Specifically, the reach set at time t, is over-approximated by the intersection of a carefully constructed parameterized family of ellipsoids E q(t), Q i(t) (t) defined as for unit vectors i (0) ∈ R d , i = 1, . . ., N .The choice of i (0) determines i (t) := exp(−A t) i (0), which in turn parameterizes the d × d symmetric positive definite shape matrix Q i(t) (t); we refer the readers to [64,Ch. 3.2], [37,Ch. 3] where the corresponding evolution equations were derived using optimal control.The center vectors q(t) ∈ R d , and the shape matrices Q i(t) (t) for this parameterized family of ellipsoids are constructed such that ∩ N i=1 E q(t), Q i(t) (t) is guaranteed to be a superset of the reach set at time t for any finite N , and for N → ∞, recovers the reach set at that time.
For the results shown in Fig. 10, we used N = 20 randomly chosen unit vectors i (0) ∈ R d .Ideally, one would like to compute the (unique) minimum volume outer ellipsoid (MVOE), a.k.a. the Löwner-John ellipsoid [65], [66]   i=1 E q(t), Q i(t) (t) .Both of these lead to solving SDP problems, and both are guaranteed to contain the Löwner-John ellipsoid of the intersection of the parameterized family of ellipsoids.These suboptimal (w.r.t. the MVOE criterion) solutions, computed using cvx [68], are shown in Fig. 10.Fig. 10 shows that the S procedure entail less conservatism compared to the MVIE scaling, in terms of volume.While the volume ratio trends in Fig. 10 are similar to that observed in Fig. 9, the approximation quality is lower.In light of the results in Sec.IV-A, this is not surprising: the integrator reach sets being zonoids (i.e., Hausdorff limit of zonotopes), the zonotopic outer-approximations are expected to perform better than other over-approximating shape primitives.
The main point here is that our results in Sec.V provide the ground truth for the size of the integrator reach set, thereby help benchmarking the performance of reach set approximation algorithms.

VII. EPILOGUE
Recap: The present paper initiates a systematic study of integrator reach set.When the input uncertainty set is hyperrectangle, we showed that the corresponding compact convex reach set R is in fact semialgebraic (Sec.IV-B) as well as a zonoid (range of an atom free vector measure) up to translation (Sec.IV-A).We derived the equation of its boundary in both parametric (Proposition 2) and implicit form (Sec. IV-C).We obtained the closed form formula for the volume (Sec.V-A) and diameter (Sec.V-B) of these reach sets.We also derived the scaling laws (Sec.V-C) for these quantities.We pointed out that these results may be used to benchmark the performance of set over-approximation algorithms (Sec.VI).What Next: In the sequel Part II, we will show how the ideas presented herein enable computing the reach sets for feedback linearizable systems.The focus will be in computing the reach set in transformed state coordinates associated with the normal form, and to map that set back to original state coordinates under diffeomorphism.This, however, requires nontrivial extension of the basic theory presented here (especially those in Proposition 2 and Sec IV-C) since we will need to handle time-dependent set-valued uncertainty in transformed control input even when the original control takes values from a set that is not time-varying.

A. Proof of Proposition 1
Since support function is distributive over sum, we have The block diagonal structure of the matrix A in (2) implies The last equality in (56) follows from (13), and from the fact [26, Prop.6.1] that the reach set remains invariant under the closure of convexification of the input set U. Substituting ( 55) and ( 56) in ( 54) yields (16).

C. Proof of Theorem 2
For s ∈ [0, t], let the vector measure µ be defined as d µ(s) := ξ(s)ds where ξ(s) is given by (13).Then t 0 | y, ξ(s) |ds is exactly in the form of a support function of a zonoid (see e.g., [28,Sec. 2]).Using the one-to-one correspondence between a compact convex set and its support function, the corresponding set is a zonoid.

E. Proof of Corollary 3
From (26), we get two different parametric representations of x bdy j in terms of (s 1 , s 2 , . . ., s rj −1 ).One parametric representation results from the choice of positive sign for the ± appearing in (26), and another for the choice of negative sign for the same.Denoting the implicit representation corresponding to the parametric representation (26) with + (resp.−) sign as p upper j (x) = 0 (resp.p lower j (x) = 0), the result follows.

F. Proof of Theorem 4
We notice that (26) gives polynomial parameterizations of the components of x bdy j for all j ∈ [m].In particular, for each k ∈ [r j ], the right hand side of ( 26) is a homogeneous polynomial in r j − 1 parameters (s 1 , s 2 , . . ., s rj −1 ) of degree r j − k + 1.By polynomial implicitization [25, p. 134], the corresponding implicit equations p upper Specifically, denote the right hand sides of (26) as g ± 1 , . . ., g ± rj for all j ∈ [m], where the superscripts indicate that either all g's are chosen with plus signs, or all with minus signs.Then write (26) as x bdy j (1) = g ± 1 s 1 , s 2 , . . ., s rj −1 , ** this may happen either because there are repeated roots in [0, t], or because some real roots exist outside [0, t], or because some roots are complex conjugates, or a combination of the previous three.
called basic semialgebraic if it can be written as a finite conjunction of polynomial inequalities and equalities, the polynomials being in R [x 1 , . . ., x d ].Finite union of basic semialgebraic sets is called a semialgebraic set.A semialgebraic set need not be basic semialgebriac; see e.g., [32, Example 2.2].

Corollary 3 .
The single input integrator reach set R j ({x 0 }, t) ⊂ R rj has two bounding surfaces for each j ∈ [m].In other words, there exist p upper j , p lower j
(a) Real algebraic curves p upper , p lower for the double integrator.(b) Real algebraic surfaces p upper , p lower for the triple integrator.

Remark 2 .
If the initial set X 0 is not singleton, then computing the volume of the R requires us to compute the volume of a Minkowski sum.Notice that vol (exp(tA)X 0 ) = |det (exp(tA)) |vol (X 0 ) = exp (trace (tA)) vol (X 0 )

Fig. 10 :
Fig. 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 Rapprox(t) for d = 2, 3, 4, plotted against time t ∈ [0, 1].Two different ellipsoidal over-approximations are shown: one (in red) based on the S procedure, and the other (in blue) obtained by scaling the maximum volume inner ellipsoid (MVIE) of the intersection of a parameterized family of ellipsoids.The results are computed for µ = 1, X0 = {0}.
of the convex set ∩ 20 i=1 E q(t), Q i(t) (t) , which is a semi-infinite programming problem [46, Ch. 8.4.1], and has no known exact semidefinite programming (SDP) reformulation.We computed two different relaxations of this problem: one based on the S procedure [67, Ch. 3.7.2],and the other by homothetic scaling of the maximum volume inner ellipsoid (MVIE) [65, Thm.III] of the set ∩