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

Shadi Haddad, Abhishek Halder, Senior Member, IEEE Shadi Haddad and Abhishek Halder are with the Department of Applied Mathematics, University of California, Santa Cruz, CA 95064, USA, {shhaddad,ahalder}@ucsc.edu

Abstract

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 over-approximation algorithms.

Keywords: Reach set, integrator, convex geometry, set-valued uncertainty.

I Introduction

Integrators with bounded controls are ubiquitous in systems-control. 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 $\bm{r}=\left(r_{1},r_{2},\ldots,r_{m}\right)^{\top}\in\mathbb{Z}_{+}^{m}$ (vector of positive integers). In other words, we consider a Brunovsky normal form with $m$ integrators where the $j$ th integrator has degree $r_{j}$ for $j\in[m]:=\{1,\ldots,m\}$ . The dynamics is

\displaystyle\dot{\bm{x}}=\bm{A}\bm{x}+\bm{B}\bm{u},\quad\bm{x}\in\mathbb{R}^{d},\quad\bm{u}(\cdot)\in\mathcal{U}\subset\mathbb{R}^{m},

(1)

where $r_{1}+r_{2}+...+r_{m}=d$ , the set $\mathcal{U}$ is compact, and


	$\displaystyle\!\!\bm{A}\!:=\!{\mathrm{blkdiag}}\!\left(\bm{A}_{1},...,\bm{A}_{m}\right),\quad\!\!\bm{B}\!:=\!{\mathrm{blkdiag}}\!\left(\bm{b}_{1},...,\bm{b}_{m}\right),$		(2a)
	$\displaystyle\!\!\bm{A}_{j}:=\left(\bm{0}_{r_{j}\times 1}~{}\bm{\|}~{}\bm{e}_{1}^{r_{j}}\bm{\|}~{}\bm{e}_{2}^{r_{j}}\bm{\|}...\bm{\|}~{}\bm{e}_{r_{j}-1}^{r_{j}}\right),\qquad\!\!\!\bm{b}_{j}:=\bm{e}_{r_{j}}^{r_{j}}.$		(2b)

In (2a), the symbol ${\mathrm{blkdiag}}\left(\cdot\right)$ denotes a block diagonal matrix whose arguments constitute its diagonal blocks. In (2b), the notation $\bm{0}_{r_{j}\times 1}$ stands for the ${r_{j}\times 1}$ column vector of zeros, and $\bm{e}_{k}^{\ell}$ denotes the $k$ th basis (column) vector in $\mathbb{R}^{\ell}$ for $k\leq\ell$ . The symbol $(...|...|...)$ denotes horizontal concatenation.

Notice that an integrator that replaces the ones appearing in the system matrices with arbitrary nonzero reals, is always reducible to the normal form (1)-(2) by renaming the variables. For instance, for any $a,b\in\mathbb{R}\setminus\{0\}$ , the system $\dot{x}_{1}=ax_{2},\dot{x}_{2}=bu_{1},$ is equivalent to $\dot{x}_{1}=\tilde{x}_{2},\dot{\tilde{x}}_{2}=\tilde{u}_{1},$ where $\tilde{x}_{2}:=ax_{2}$ , $\tilde{u}_{1}:=abu_{1}$ .

Let $\mathcal{R}\left(\mathcal{X}_{0},t\right)$ denote the forward reach set of (1) at time $t>0$ , starting from a given compact convex set of initial conditions $\mathcal{X}_{0}\subset\mathbb{R}^{d}$ , i.e.,

	$\displaystyle\mathcal{R}$	$\displaystyle\left(\mathcal{X}_{0},t\right):=\!\!\!\bigcup_{\text{measurable}\>\bm{u}(\cdot)\in\mathcal{U}\subset\mathbb{R}^{m}}\!\!\!\big{\{}\bm{x}(t)\in\mathbb{R}^{d}\mid\text{\eqref{IntegratorDyn} and \eqref{defAB} hold},$
		$\displaystyle\qquad\bm{x}(t=0)\in\mathcal{X}_{0}\;\text{compact convex},\quad\mathcal{U}\;\text{compact}\big{\}}.$		(3)

In words, $\mathcal{R}\left(\mathcal{X}_{0},t\right)$ is the set of all states that the controlled dynamics (1)-(2) can reach at time $t>0$ , starting from the set $\mathcal{X}_{0}$ at $t=0$ , with measurable control $\bm{u}(\cdot)\in\mathcal{U}$ compact. Formally,

	$\displaystyle\!\!\mathcal{R}\left(\mathcal{X}_{0},t\right)\!$	$\displaystyle=\exp(t\bm{A})\mathcal{X}_{0}\dotplus\!\!\int_{0}^{t}\!\!\!\exp\left((t-\tau)\bm{A}\right)\bm{B}\mathcal{U}\>{\rm{d}}\tau$
		$\displaystyle=\exp(t\bm{A})\mathcal{X}_{0}\dotplus\!\!\int_{0}^{t}\!\!\!\exp\left(s\bm{A}\right)\bm{B}\mathcal{U}\>{\rm{d}}s,$		(4)

where $\dotplus$ denotes the Minkowski sum. The set-valued integral [5] in (4) is defined for any point-to-set function $F(\cdot)$ , as

\displaystyle\int_{0}^{t}F(s){\rm{d}}s:=\lim_{\Delta\downarrow 0}\>\sum_{i=0}^{\lfloor t/\Delta\rfloor}\Delta F(i\Delta),

(5)

where the summation symbol $\Sigma$ denotes the Minkowski sum, and $\lfloor\cdot\rfloor$ 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., $\mathcal{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.

II-1 State transition matrix

For $0\leq s<t$ , the state transition matrix $\bm{\Phi}(t,s)$ associated with (2) is

	$\displaystyle\bm{\Phi}(t,s)$	$\displaystyle\equiv\exp(\bm{A}(t-s))$
		$\displaystyle={\mathrm{blkdiag}}\left(\exp(\bm{A}_{1}(t-s)),\ldots,\exp(\bm{A}_{m}(t-s))\right),$

with each diagonal block is upper triangular. Specifically, the $j$ th diagonal block of size $r_{j}\times r_{j}$ is written element-wise as

\displaystyle\exp(\bm{A}_{j}(t-s)):=\begin{cases}\dfrac{(t-s)^{\ell-k}}{(\ell-k)!}&\text{for}\;k\leq\ell,\\ 0&\text{otherwise},\end{cases}

(6)

where $k$ is the row index, $\ell$ is the column index, and $k,\ell\in[r_{j}]$ for each $j\in[m]$ . The diagonal entries in (6) are unity.

II-2 Support function

The support function $h_{\mathcal{K}}(\cdot)$ of a compact convex set $\mathcal{K}\subset\mathbb{R}^{d}$ , is given by

\displaystyle h_{\mathcal{K}}(\bm{y}):=\underset{\bm{x}\in\mathcal{K}}{\sup}~{}\langle\bm{y},\bm{x}\rangle,\quad\bm{y}\in\mathbb{R}^{d},

(7)

where $\langle\cdot,\cdot\rangle$ denotes the standard Euclidean inner product. Geometrically, $h_{\mathcal{K}}(\bm{y})$ gives the signed distance of the supporting hyperplane of $\mathcal{K}$ with outer normal vector $\bm{y}$ , measured from the origin. Furthermore, the supporting hyperplane at $\bm{x}^{\text{bdy}}\in\partial\mathcal{K}$ is $\langle\bm{y},\bm{x}^{\text{bdy}}\rangle=h_{\mathcal{K}}(\bm{y})$ , and we can write

\mathcal{K}=\big{\{}\bm{x}\in\mathbb{R}^{d}\mid\langle\bm{y},\bm{x}\rangle\leq h_{\mathcal{K}}(\bm{y})\;\text{for all}\;\bm{y}\in\mathbb{R}^{d}\big{\}}.

For compact $\mathcal{K}_{1},\mathcal{K}_{2}\subset\mathbb{R}^{d}$ ,

\displaystyle\mathcal{K}_{1}\subseteq\mathcal{K}_{2}\quad\text{if and only if}\quad h_{\mathcal{K}_{1}}(\cdot)\leq h_{\mathcal{K}_{2}}(\cdot).

(8)

The support function $h_{\mathcal{K}}\left(\bm{y}\right)$ is convex in $\bm{y}$ . For more details on the support function, we refer the readers to [8, Ch. V].

The support function $h_{\mathcal{K}}(\bm{y})$ uniquely determines the set $\mathcal{K}$ . Given matrix-vector pair $(\bm{\Gamma},\bm{\gamma})\in\mathbb{R}^{d\times d}\times\mathbb{R}^{d}$ , the support function of the affine transform $\bm{\Gamma}\mathcal{K}+\bm{\gamma}$ is

\displaystyle h_{\bm{\Gamma}\mathcal{K}+\bm{\gamma}}(\bm{y})=\langle\bm{y},\bm{\gamma}\rangle+h_{\mathcal{K}}(\bm{\Gamma}^{\top}\bm{y}).

(9)

Given a function $f:\mathbb{R}^{d}\mapsto\mathbb{R}\cup\{+\infty\}$ , its Legendre-Fenchel conjugate is

\displaystyle f^{*}(\bm{y}):=\underset{\bm{x}\in\operatorname{domain}(f)}{\sup}\>\{\langle\bm{y},\bm{x}\rangle-f(\bm{x})\},\quad\bm{y}\in\mathbb{R}^{d}.

(10)

From (7)-(10), it follows that $h_{\mathcal{K}}(\bm{y})$ is the Legendre-Fenchel conjugate of the indicator function

\mathbf{1}_{\mathcal{K}}(\bm{x}):=\begin{cases}0&\text{if}\quad\bm{x}\in\mathcal{K},\\ +\infty&\text{otherwise.}\end{cases}

Since the indicator function of a convex set is a convex function, the biconjugate $\mathbf{1}_{\mathcal{K}}^{**}(\cdot)=h_{\mathcal{K}}^{*}(\cdot)=\mathbf{1}_{\mathcal{K}}(\cdot)$ . This will be useful in Section IV.

To proceed further, we introduce some notations. Since $\mathcal{U}$ is compact, let

\displaystyle\alpha_{j}:=\underset{\bm{u}\in\mathcal{U}}{\min}\;u_{j},\quad\beta_{j}:=\underset{\bm{u}\in\mathcal{U}}{\max}\;u_{j},\quad j\in[m],

(11)

that is, $\alpha_{j}$ and $\beta_{j}$ are the component-wise minimum and maximum, respectively, of the input vector. Furthermore, let

\displaystyle\mu_{j}:=\frac{\beta_{j}-\alpha_{j}}{2},\quad\nu_{j}:=\frac{\beta_{j}+\alpha_{j}}{2},

(12)

and introduce

\displaystyle\!\!\bm{\xi}(s):=\!\begin{pmatrix}\mu_{1}\bm{\xi}_{1}(s)\\ \vdots\\ \mu_{m}\bm{\xi}_{m}(s)\end{pmatrix},\;\bm{\xi}_{j}(s):=\!\begin{pmatrix}s^{r_{j}-1}/(r_{j}-1)!\\ s^{r_{j}-2}/(r_{j}-2)!\\ \vdots\\ s\\ 1\end{pmatrix},

(13)

for $j\in[m]$ . Also, let

\displaystyle\bm{\zeta}(t_{0},t)=\begin{pmatrix}\mu_{1}\bm{\zeta}_{1}(t_{0},t)\\ \mu_{2}\bm{\zeta}_{2}(t_{0},t)\\ \vdots\\ \mu_{m}\bm{\zeta}_{m}(t_{0},t)\end{pmatrix},\quad\bm{\zeta}_{j}(t_{0},t):=\int_{t_{0}}^{t}\!\bm{\xi}_{j}(s)\>{\rm{d}}s\in\mathbb{R}^{r_{j}},

(14)

for $j\in[m]$ . When $t_{0}=0$ , we simplify the notations as

\displaystyle\!\!\!\!\!\!\bm{\zeta}(t):=\bm{\zeta}(0,t),\;\bm{\zeta}_{j}(t):=\bm{\zeta}_{j}(0,t)\;\text{for all}\;j\in[m].

(15)

Using (13) and following (7), we deduce Proposition 1 stated next (proof in Appendix -A).

Proposition 1.

(Support function for compact $\mathcal{U}$ ) For compact convex $\mathcal{X}_{0}\subset\mathbb{R}^{d}$ , and compact $\mathcal{U}\subset\mathbb{R}^{m}$ , the support function of the reach set (4) is

	$\displaystyle h_{\mathcal{R}\left(\mathcal{X}_{0},t\right)}\left(\bm{y}\right)=\underset{\bm{x}_{0}\in\mathcal{X}_{0}}{\sup}\sum_{j=1}^{m}~{}\langle\bm{y}_{j},\exp\left(t\bm{A}_{j}\right)\bm{x}_{j0}\rangle$
	$\displaystyle+\int^{t}_{0}\underset{\bm{u}\in{\rm{closure}}\left({\rm{conv}}\left(\mathcal{U}\right)\right)}{\sup}\sum_{j=1}^{m}~{}\{\langle\bm{y}_{j},\bm{\xi}_{j}(s)\rangle\>u_{j}\}\>{\rm{d}}{s},$		(16)

where ${\rm{conv}}(\cdot)$ denotes the convex hull.

II-3 Polar dual

The polar dual $\mathcal{K}^{\circ}$ of any non-empty set $\mathcal{K}\subset\mathbb{R}^{d}$ is given by

\displaystyle\mathcal{K}^{\circ}:=\{\bm{y}\in\mathbb{R}^{d}\mid\langle\bm{y},\bm{x}\rangle\leq 1\quad\text{for all}\;\bm{x}\in\mathcal{K}\}.

(17)

From this definition, it is immediate that $\mathcal{K}^{\circ}$ contains the origin, and is a closed convex set. The bipolar $\left(\mathcal{K}^{\circ}\right)^{\circ}={\rm{closure}}\left({\rm{conv}}\left(\mathcal{K}\cup\{\bm{0}\}\right)\right)$ . Thus, if $\mathcal{K}$ is compact convex and contains the origin, then we have the involution $\left(\mathcal{K}^{\circ}\right)^{\circ}=\mathcal{K}$ . From (7) and (17), notice that $\mathcal{K}^{\circ}$ is the unit support function ball, i.e., $\mathcal{K}^{\circ}=\{\bm{y}\in\mathbb{R}^{d}\mid h_{\mathcal{K}}(\bm{y})\leq 1\}$ . In Sec. IV-D, we will mention some properties of the polar dual of the integrator reach set.

II-4 Vector measure

Let $\mathcal{F}$ be a $\sigma$ -field of the subsets of a set. A countably additive mapping $\widetilde{\bm{\mu}}:\mathcal{F}\mapsto\mathbb{R}^{d}$ is termed a vector measure. Here, “countably additive” means that for any sequence $\{\Omega_{i}\}_{i=1}^{\infty}$ of disjoint sets in $\mathcal{F}$ such that their union is in $\mathcal{F}$ , we have $\widetilde{\bm{\mu}}\left(\cup_{i=1}^{\infty}\Omega_{i}\right)=\sum_{i=1}^{\infty}\widetilde{\bm{\mu}}\left(\Omega_{i}\right)<\infty$ . Some of the early investigations of vector measures were due to Liapounoff [9] and Halmos [10]; relatively recent references are [11, 12].

II-5 Zonotope

A zonotope $\mathcal{Z}\subset\mathbb{R}^{d}$ is a finite Minkowski sum of closed line segments or intervals $\{I_{i}\}_{i=1}^{n}$ where these intervals are imbedded in the ambient Euclidean space $\mathbb{R}^{d}$ . Explicitly, for some positive integer $n$ , we write

	$\displaystyle\mathcal{Z}$	$\displaystyle:=I_{1}\dotplus\ldots\dotplus I_{n}$
		$\displaystyle=\bigg{\{}\bm{x}\in\mathbb{R}^{d}\mid\bm{x}=\sum_{i=1}^{n}\bm{x}_{i},\quad\bm{x}_{i}\in I_{i},\quad i\in[n]\bigg{\}}.$

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 $\{\bm{x}\in\mathbb{R}^{d}\mid\|\bm{x}\|_{1}\leq 1\}$ , is not a zonotope. Standard references on zonotope include [14, 15], [16, Ch. 2.7].

The set of zonotopes is closed under affine image and Minkowski sum, but not under intersection. In the systems-control literature, a significant body of work exists on computationally efficient over-approximation of reach sets via zonotopes [17, 18, 19] and its variants such as zonotope bundles [20], constrained zonotopes [21], complex zonotopes [22], and polynomial zonotopes [23, 24].

II-6 Variety and ideal

Let $p_{1},\ldots,p_{n}\in\mathbb{R}[x_{1},\ldots,x_{d}]$ , the vector space of real-valued $d$ -variate polynomials. The (affine) variety $V_{\mathbb{R}[x_{1},\ldots,x_{d}]}(p_{1},\ldots,p_{n})$ is the set of all solutions of the system $p_{1}(x_{1},x_{2},\ldots,x_{d})=\ldots=p_{n}(x_{1},x_{2},\ldots,x_{d})=0$ . Given $p_{1},\ldots,p_{n}\in\mathbb{R}[x_{1},\ldots,x_{d}]$ , the set

I:=\bigg{\{}\sum_{i=1}^{n}\alpha_{i}p_{i}\mid\alpha_{1},\ldots,\alpha_{n}\in\mathbb{R}[x_{1},\ldots,x_{d}]\bigg{\}}

is called the ideal generated by $p_{1},\ldots,p_{n}$ . We write this symbolically as $I=\langle\langle p_{1},\ldots,p_{n}\rangle\rangle$ . Roughly speaking, $\langle\langle p_{1},\ldots,p_{n}\rangle\rangle$ is the set of all polynomial consequences of the given system of $n$ polynomial equations in $d$ indeterminates. We refer the readers to [25, Ch. 1] for detailed exposition of these concepts.

III Box-valued Input Uncertainty

In the remaining of this paper, we characterize the exact reach set (3) when input set $\mathcal{U}\subset\mathbb{R}^{m}$ is box-valued, and remark on the quality of approximation for the same when $\mathcal{U}$ is arbitrary compact.

When $\mathcal{U}\subset\mathbb{R}^{m}$ is box-valued, denote the reach set (3) as $\mathcal{R}^{\Box}$ , i.e., with a box superscript^*^**For the single input ( $m=1$ ) case, we drop the box superscript.. In this case, each of the $m$ single input integrator dynamics with $r_{j}$ dimensional state subvectors for $j\in[m]$ , are decoupled from each other. Then $\mathcal{R}^{\Box}\left(\mathcal{X}_{0},t\right)\subset\mathbb{R}^{d}$ is the Cartesian product of these single input integrator reach sets: $\mathcal{R}_{j}\left(\mathcal{X}_{0},t\right)\subset\mathbb{R}^{r_{j}}$ for $j\in[m]$ , i.e.,

\displaystyle\mathcal{R}^{\Box}=\mathcal{R}_{1}\times\mathcal{R}_{2}\times\ldots\times\mathcal{R}_{m}.

(18)

In what follows, we will sometimes exploit that (18) may also be written as^†^††In general, the Minkowski sum of a given collection of compact convex sets is not equal to their Cartesian product. However, the “factor sets” in (18) belong to disjoint mutually orthogonal $r_{j}$ dimensional subspaces, $j=1,\ldots,m$ , which allows writing this Cartesian product as a Minkowski sum. a Minkowski sum $\mathcal{R}_{1}\dotplus\ldots\dotplus\mathcal{R}_{m}$ . Notice that the decoupled dynamics also allows us to write a Minkowski sum decomposition for the set of initial conditions

\mathcal{X}_{0}=\mathcal{X}_{10}\dotplus\ldots\dotplus\mathcal{X}_{m0},

and accordingly, the initial condition subvectors $\bm{x}_{j0}\in\mathcal{X}_{j0}\subset\mathbb{R}^{r_{j}}$ for $j\in[m]$ . Thus $\bm{x}_{0}=(\bm{x}_{10},\ldots,\bm{x}_{m0})^{\top}$ .

Since the support function of the Minkowski sum is equal to the sum of the support functions, we have

\displaystyle h_{\mathcal{R}^{\Box}\left(\mathcal{X}_{0},t\right)}(\bm{y})=\sum_{j=1}^{m}h_{\mathcal{R}_{j}\left(\mathcal{X}_{j0},t\right)}(\bm{y}_{j}).

(19)

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 $\mathcal{U}$ ) For compact convex $\mathcal{X}_{0}\subset\mathbb{R}^{d}$ , and box-valued input uncertainty set given by

\displaystyle\mathcal{U}:=\left[\alpha_{1},\beta_{1}\right]\times\left[\alpha_{2},\beta_{2}\right]\times\ldots\times\left[\alpha_{m},\beta_{m}\right]\subset\mathbb{R}^{m}\>,

(20)

the support function of the reach set (4) is

	$\displaystyle h_{\mathcal{R}^{\Box}\left(\mathcal{X}_{0},t\right)}\left(\bm{y}\right)=\sum_{j=1}^{m}\bigg{\{}\underset{\bm{x}_{j0}\in\mathcal{X}_{j0}}{\sup}\langle\bm{y}_{j},\exp\left(t\bm{A}\right)\bm{x}_{j0}\rangle$
	$\displaystyle\qquad\qquad\quad+\nu_{j}\langle\bm{y}_{j},\bm{\zeta}_{j}(t)\rangle+\mu_{j}\int_{0}^{t}\!\lvert\langle\bm{y}_{j},\bm{\xi}_{j}(s)\rangle\rvert\>{\rm{d}}s\bigg{\}}.$		(21)

The formula (21) upper bounds (16) resulting from the same initial condition and arbitrary compact $\mathcal{U}\subset\mathbb{R}^{m}$ with $\{\alpha_{j},\beta_{j}\}_{j=1}^{m}$ related to $\mathcal{U}$ via (11). Thus, from (8), the reach set $\mathcal{R}^{\Box}$ with box-valued input uncertainty will over-approximate the reach set $\mathcal{R}$ associated with arbitrary compact $\mathcal{U}$ , at any given $t>0$ , provided $\{\alpha_{j},\beta_{j}\}_{j=1}^{m}$ are defined as (11).

When $\mathcal{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 ${\mathrm{dist}}$ between the convex compact sets $\mathcal{R}^{\Box},\mathcal{R}\subset\mathbb{R}^{d}$ , expressible [13, Thm. 1.8.11] in terms of their support functions $h_{\mathcal{R}^{\Box}}(\cdot),h_{\mathcal{R}}(\cdot)$ as

\displaystyle{\mathrm{dist}}\left(\mathcal{R}^{\Box},\mathcal{R}\right)=\underset{\|\bm{y}\|_{2}=1}{\sup}\quad\big{|}h_{\mathcal{R}^{\Box}}(\bm{y})-h_{\mathcal{R}}(\bm{y})\big{|}.

(22)

Thanks to (8), the absolute value in (22) can be dispensed since $\mathcal{R}\subseteq\mathcal{R}^{\Box}$ with set equality if $\mathcal{U}$ is box, in which case $h_{\mathcal{R}^{\Box}}(\cdot)=h_{\mathcal{R}}(\cdot)$ and ${\mathrm{dist}}=0$ .

It is known that [26, Prop. 6.1] the set $\mathcal{R}\left(\mathcal{X}_{0},t\right)$ resulting from a linear time invariant dynamics such as (1)-(2) remains invariant under the closure of convexification of the input set $\mathcal{U}$ . Therefore, it is possible that $\mathcal{R}=\mathcal{R}^{\Box}$ and ${\mathrm{dist}}=0$ even when the compact set $\mathcal{U}$ is nonconvex. For instance, the reach set $\mathcal{R}\left(\mathcal{X}_{0},t\right)$ resulting from some compact convex $\mathcal{X}_{0}\subset\mathbb{R}^{d}$ and dynamics (1)-(2) with the nonconvex input uncertainty set $\{-1,1\}^{m}$ , is identical to $\mathcal{R}^{\Box}\left(\mathcal{X}_{0},t\right)$ resulting from the same $\mathcal{X}_{0}$ , same dynamics, and the box-valued input uncertainty set (20) with $\alpha_{j}=-1,\beta_{j}=1$ for all $j\in[m]$ .

Likewise, for the same compact convex $\mathcal{X}_{0}\subset\mathbb{R}^{d}$ , the reach set $\mathcal{R}\left(\mathcal{X}_{0},t\right)$ resulting from (1)-(2) with the nonconvex input set $\{\bm{u}\in\mathbb{R}^{m}\mid\|\bm{u}\|_{p}\leq 1\}$ , $0<p<1$ , is the same as that resulting from the cross-polytope $\{\bm{u}\in\mathbb{R}^{m}\mid\|\bm{u}\|_{1}\leq 1\}$ . More generally, for $0<p<\infty$ , suppose $\mathcal{R}_{p}\left(\mathcal{X}_{0},t\right)$ results from the unit $p$ norm ball input uncertainty set $\{\bm{u}\in\mathbb{R}^{m}\mid\|\bm{u}\|_{p}\leq 1\}$ . Let $\bm{M}^{\top}(\tau):=\exp\left(\tau\bm{A}\right)\bm{B}={\mathrm{blkdiag}}\left(\bm{\xi}_{1},\ldots,\bm{\xi}_{m}\right)$ . If $\mathcal{R}^{\Box}\left(\mathcal{X}_{0},t\right)$ results from the same $\mathcal{X}_{0}$ , same dynamics, and input uncertainty set (20) with $\alpha_{j}=-1,\beta_{j}=1$ for all $j\in[m]$ , then using [27, Thm. 1], (22) simplifies to

\displaystyle\!\!\!{\mathrm{dist}}\!\left(\!\mathcal{R}^{\Box},\mathcal{R}_{p}\!\right)\!=\!\!\underset{\|\bm{y}\|_{2}=1}{\sup}\int_{0}^{t}\!\!\!\left(\|\bm{M}(\tau)\bm{y}\|_{1}\!-\!\|\bm{M}(\tau)\bm{y}\|_{q}\right){\rm{d}}\tau

(23)

where $q$ is the Hölder conjugate of $\max\{1,p\}$ , i.e., $\frac{1}{\max\{1,p\}}+\frac{1}{q}=1$ , and $1<q\leq\infty$ . In this case, the positive value (23) quantifies the quality of strict over-approximation $\mathcal{R}_{p}\subset\mathcal{R}^{\Box}$ for $0<p<\infty$ . The objective in (23) being positive homogeneous, admits lossless constraint convexification $\|\bm{y}\|_{2}\leq 1$ , and the corresponding maximal value^‡^‡‡As such, (23) has a difference of convex objective, and by the Weierstrass extreme value theorem, the maximum is achieved. for moderate dimensions $d$ , can be found by direct numerical search.

IV Taxonomy and Boundary

For $\mathcal{X}_{0}\subset\mathbb{R}^{d}$ compact convex, it is well-known [1, Sec. 2] that the reach set $\mathcal{R}$ given by (3) is compact convex for all $t>0$ provided $\mathcal{U}$ is compact. However, it is not immediate what kind of convex set $\mathcal{R}$ is, even for singleton $\mathcal{X}_{0}\equiv\{\bm{x}_{0}\}$ .

In this Section, we examine the question “what type of compact convex set $\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right)$ is” when $\mathcal{U}$ is box-valued uncertainty set of the form (20). In the same setting, we also derive the equations for the boundary $\partial\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right)$ .

Notice that for non-singleton $\mathcal{X}_{0}$ , the taxonomy question is not well-posed since the classification then will depend on $\mathcal{X}_{0}$ . Also, setting $\mathcal{X}_{0}\equiv\{\bm{x}_{0}\}$ in (4), it is apparent that $\mathcal{R}\left(\{\bm{x}_{0}\},t\right)$ is a translation of the set-valued integral in (4). Thus, classifying $\mathcal{R}\left(\{\bm{x}_{0}\},t\right)$ amounts to classifying the second summand in (4).

IV-A $\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right)$ is a Zonoid

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, 29, 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 $\mathbb{R}^{2}$ is a zonoid. In dimensions three or more, all $\ell_{p}$ norm balls for $p\geq 2$ are zonoids.

An alternative way to think about the zonoid is to view it as the limiting set (convergence with respect to the two-sided Hausdorff distance, see e.g., [4, Appendix B]) of the Minkowski sum of line segments, i.e., the limit of a sequence of zonotopes [28, 14, 15]. Formally, given a Hausdorff convergent sequence of zonotopes $\!\{\mathcal{Z}_{j}\}$ , the zonoid $\mathcal{Z}_{\infty}$ is

\mathcal{Z}_{\infty}:=\lim_{j\rightarrow\infty}\mathcal{Z}_{j},\;\text{where}\;\mathcal{Z}_{j}:=\sum_{i=1}^{n(j)}\left[\bm{a}_{ij},\bm{b}_{ij}\right],\quad\bm{a}_{ij},\bm{b}_{ij}\in\mathbb{R}^{d},

for some $\bm{a}_{ij}\leq\bm{b}_{ij}$ (element-wise vector inequality), and a suitable mapping $n:\mathbb{Z}_{+}\mapsto\mathbb{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 $\mathcal{R}^{\Box}$ given by (3) with $\mathcal{X}_{0}\equiv\{\bm{x}_{0}\}$ and $\mathcal{U}$ given by (20), is a zonoid.

To appreciate Theorem 2 via the limiting viewpoint mentioned before, let us write

	$\displaystyle\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right)=$	$\displaystyle\underbrace{\exp(t\bm{A})\bm{x}_{0}+\displaystyle\sum_{j=1}^{m}\nu_{j}\displaystyle\bm{\zeta}_{j}(t)\vphantom{\displaystyle\sum_{j=1}^{m}\displaystyle\lim_{n\rightarrow\infty}\displaystyle\sum_{i=0}^{n}\frac{t}{n}\mu_{i}\bm{\xi}_{j}(t_{i})\left[-1,1\right]}}_{\textup{first term}}$
		$\displaystyle\dotplus\underbrace{\displaystyle\sum_{j=1}^{m}\displaystyle\lim_{n\rightarrow\infty}\displaystyle\sum_{i=0}^{n}\frac{t}{n}\mu_{j}\displaystyle\bm{\xi}_{j}(t_{i})\left[-1,1\right]}_{\textup{second term}},$		(24)

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 $\mathcal{R}_{1}\dotplus\ldots\dotplus\mathcal{R}_{m}$ . Furthermore, uniformly discretizing $[0,t]$ into $n$ subintervals $[(i-1)t/n,it/n)$ , $i=1,\ldots,n$ , we write $\int_{0}^{t}\exp(s\bm{A}_{j})\bm{b}_{j}[-\mu_{j},\mu_{j}]{\rm{d}}s$ 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 $\mathbb{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 $\{\mathcal{Z}_{n}\}$ comprising of zonotopes

\mathcal{Z}_{n}:=\displaystyle\sum_{i=0}^{n}\frac{t}{n}\mu_{j}\bm{\xi}_{j}(t_{i})\left[-1,1\right],

which are the Minkowski sum of $n+1$ line segments. Since $\lim_{n\rightarrow\infty}\mathcal{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.

Remark 1.

If $\mathcal{X}_{0}\subset\mathbb{R}^{d}$ is not singleton, but instead a zonoid, then $\mathcal{R}^{\Box}\left(\mathcal{X}_{0},t\right)$ is still a (translated) zonoid. To see this, notice from (4) and (21) that

\displaystyle\mathcal{R}^{\Box}\left(\mathcal{X}_{0},t\right)=\exp(t\bm{A})\mathcal{X}_{0}\dotplus\mathcal{R}^{\Box}\left(\{\bm{0}\},t\right),

(25)

and that $\exp(t\bm{A})\mathcal{X}_{0}$ , being linear image of a zonoid, is a zonoid [28, Lemma 1.4]. Thus, (25) being Minkowski sum of zonoids, is a zonoid too [28, Thm. 1.5], up to translation.

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 $\mathcal{X}_{0}=\{\bm{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) $\mathcal{X}_{0}$ . In this sense, our closed form formulae are as general as one might hope for. For a specific non-singleton $\mathcal{X}_{0}$ , one can use these formulae to first derive the boundary (resp. volume) of $\mathcal{R}^{\Box}\left(\{\bm{0}\},t\right)$ , and then use (25) to get numerical estimates for the boundary (resp. volume) of $\mathcal{R}^{\Box}\left(\mathcal{X}_{0},t\right)$ (cf. Remark 2).

IV-B $\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right)$ is Semialgebraic

A set in $\mathbb{R}^{d}$ is called basic semialgebraic if it can be written as a finite conjunction of polynomial inequalities and equalities, the polynomials being in $\mathbb{R}\left[x_{1},\ldots,x_{d}\right]$ . 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].

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 $\bm{x}^{\text{bdy}}\in\partial\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right)$ , the boundary of the reach set. Then we use this representation to establish semialgebraicity of $\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right)$ in Theorem 4 that follows.

Proposition 2.

For relative degree vector $\bm{r}=(r_{1},\ldots,r_{m})^{\top}$ , and fixed $\bm{x}_{0}\in\mathbb{R}^{d}$ comprising of subvectors $\bm{x}_{j0}\in\mathbb{R}^{r_{j}}$ where $j\in[m]$ , consider the reach set (4) with singleton $\mathcal{X}_{0}\equiv\{\bm{x}_{0}\}$ and $\mathcal{U}$ given by (20). For $j\in[m]$ , define $\mu_{1},\ldots,\mu_{m}$ and $\nu_{1},\ldots,\nu_{m}$ as in (11)-(12). Let the indicator function $\mathbf{1}_{k\leq\ell}:=1$ for $k\leq\ell$ , and $:=0$ otherwise. Then the components of

\bm{x}^{\textup{bdy}}=\begin{pmatrix}\bm{x}^{\textup{bdy}}_{1}\\ \bm{x}^{\textup{bdy}}_{2}\\ \vdots\\ \bm{x}^{\textup{bdy}}_{m}\end{pmatrix}\in\partial\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right),\;\bm{x}^{\textup{bdy}}_{j}\in\mathbb{R}^{r_{j}},\;j\in[m],

admit parametric representation in terms of the parameters $(s_{1},s_{2},\ldots,s_{r_{j}-1})$ satisfying $0\leq s_{1}\leq s_{2}\leq\ldots\leq s_{r_{j}-1}\leq t$ . This parameterization is given by

	$\displaystyle\bm{x}^{\textup{bdy}}_{j}\left(k\right)=\sum_{\ell=1}^{r_{j}}\mathbf{1}_{k\leq\ell}\>\frac{t^{\ell-k}}{(\ell-k)!}\>\bm{x}_{j0}(\ell)+\frac{\nu_{j}\>t^{r_{j}-k+1}}{(r_{j}-k+1)!}$
	$\displaystyle\pm\frac{\mu_{j}}{(r_{j}-k+1)!}\bigg{\{}(-1)^{r_{j}-1}\>t^{r_{j}-k+1}+2\sum_{q=1}^{r_{j}-1}\!\!(-1)^{q+1}\>s_{q}^{r_{j}-k+1}\bigg{\}},$		(26)

where $\bm{x}^{\textup{bdy}}_{j}\left(k\right)$ denotes the $k$ th component of the $j$ th subvector $\bm{x}^{\textup{bdy}}_{j}$ for $k\in[r_{j}]$ .

The following is a consequence of the $\pm$ appearing in (26).

Corollary 3.

The single input integrator reach set $\mathcal{R}_{j}\left(\{\bm{x}_{0}\},t\right)\subset\mathbb{R}^{r_{j}}$ has two bounding surfaces for each $j\in[m]$ . In other words, there exist $p_{j}^{\textup{upper}},p_{j}^{\textup{lower}}:\mathbb{R}^{r_{j}}\mapsto\mathbb{R}$ such that

\mathcal{R}_{j}\left(\{\bm{x}_{0}\},t\right)=\{\bm{x}\in\mathbb{R}^{r_{j}}\mid p_{j}^{\textup{upper}}(\bm{x})\leq 0,\;p_{j}^{\textup{lower}}(\bm{x})\leq 0\},

with boundary $\partial\mathcal{R}_{j}\left(\{\bm{x}_{0}\},t\right)=\{\bm{x}\in\mathbb{R}^{r_{j}}\mid p_{j}^{\textup{upper}}(\bm{x})=0\}\cup\{\bm{x}\in\mathbb{R}^{r_{j}}\mid p_{j}^{\textup{lower}}(\bm{x})=0\}$ .

Refer to caption — Figure 1: The “almond-shaped” integrator reach set $\mathcal{R}(\{\bm{x}_{0}\},t)\subset\mathbb{R}^{3}$ with $d=3$ , $m=1$ , $\bm{x}_{0}=(0.1,0.2,0.3)^{\top}$ , $\mathcal{U}\equiv[\alpha,\beta]=[-1,1]$ at $t=2.1$ . The wireframes correspond to the upper and lower surfaces.

During the proof of Theorem 4 below, it will turn out that in fact $p_{j}^{\text{upper}},p_{j}^{\text{lower}}\in\mathbb{R}\left[x_{1},\ldots,x_{r_{j}}\right]$ for all $j\in[m]$ . In words, $p_{j}^{\text{upper}},p_{j}^{\text{lower}}$ are real algebraic hypersurfaces for all $j\in[m]$ .

Let us exemplify the parameterization (26) for the case $\bm{r}=(r_{1},r_{2})^{\top}=(2,3)^{\top}$ . In this case,

\displaystyle\!\!\begin{pmatrix}\bm{x}^{\text{bdy}}_{1}(1)\\ \vspace*{-0.05in}\\ \bm{x}^{\text{bdy}}_{1}(2)\end{pmatrix}\!\!=\!\!\begin{pmatrix}\bm{x}_{10}(1)+t\bm{x}_{10}(2)+\nu_{1}(t^{2}/2)\pm\mu_{1}\left(s_{1}^{2}-t^{2}/2\right)\\ \vspace*{-0.05in}\\ \bm{x}_{10}(2)+\nu_{1}t\pm\mu_{1}\left(2s_{1}-t\right)\end{pmatrix},

(27)

\displaystyle\!\!\begin{pmatrix}\bm{x}^{\text{bdy}}_{2}(1)\\ \vspace*{-0.05in}\\ \bm{x}^{\text{bdy}}_{2}(2)\\ \vspace*{-0.05in}\\ \bm{x}^{\text{bdy}}_{2}(3)\end{pmatrix}=\begin{pmatrix}\bm{x}_{20}(1)+t\bm{x}_{20}(2)+(t^{2}/2)\bm{x}_{20}(3)\\ +\nu_{2}(t^{3}/6)\pm\mu_{2}\left(t^{3}/6+2s_{1}^{3}/6-2s_{2}^{3}/6\right)\\ \vspace*{-0.05in}\\ \bm{x}_{20}(2)+t\bm{x}_{20}(3)+\nu_{2}(t^{2}/2)\pm\mu_{2}\left(t^{2}/2\right.\\ \left.+2s_{1}^{2}/2-2s_{2}^{2}/2\right)\\ \vspace*{-0.05in}\\ \bm{x}_{20}(3)+\nu_{2}t\pm\mu_{2}\left(t+2s_{1}-2s_{2}\right)\end{pmatrix}.

(28)

In (27), taking plus (resp. minus) signs in each of component gives the parametric representation of the curve $p_{1}^{\text{upper}}=0$ (resp. $p_{1}^{\text{lower}}=0$ ). These curves are as in [4, Fig. 1(a)], and their union defines $\partial\mathcal{R}_{1}$ . We note that the parameterization (27) appeared in [37, p. 111].

Likewise, in (28), taking plus (resp. minus) signs in each of component gives the parametric representation of the surface $p_{2}^{\text{upper}}(\bm{x})=0$ (resp. $p_{2}^{\text{lower}}=0$ ). The resulting set $\mathcal{R}_{2}$ is the triple integrator reach set, and is shown in Fig. 1.

Now we come to the main result of this subsection.

Theorem 4.

The reach set $\mathcal{R}^{\Box}$ given by (3) with $\mathcal{X}_{0}\equiv\{\bm{x}_{0}\}$ and $\mathcal{U}$ as in (20), is semialgebraic.

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_{1}^{\text{upper}},p_{1}^{\text{lower}}$ are parabolas. In particular,

	$\displaystyle p_{1}^{\text{upper}}(\bm{x}^{\text{bdy}}_{1}(1),\bm{x}^{\text{bdy}}_{1}(2))=\frac{1}{4}\left(\frac{\bm{x}^{\text{bdy}}_{1}(2)-\bm{x}_{10}(2)-\nu_{1}t}{\mu_{1}}+t\right)^{\!2}$
	$\displaystyle\quad-\frac{\bm{x}^{\text{bdy}}_{1}(1)-\bm{x}_{10}(1)-t\bm{x}_{10}(2)-\nu_{1}\frac{t^{2}}{2}}{\mu_{1}}-\frac{t^{2}}{2},$		(29)

and the formula for $p_{1}^{\text{lower}}(\bm{x}^{\text{bdy}}_{1}(1),\bm{x}^{\text{bdy}}_{1}(2))$ follows mutatis mutandis.

Similarly, eliminating the parameters $s_{1},s_{2}$ from (28) reveals that $p_{2}^{\text{upper}},p_{2}^{\text{lower}}$ are quartic polynomials. In particular,

	$\displaystyle p_{2}^{\text{upper}}(\bm{x}^{\text{bdy}}_{2}(1),\bm{x}^{\text{bdy}}_{2}(2),\bm{x}^{\text{bdy}}_{2}(3))$
	$\displaystyle=\frac{1}{16}\left(\frac{\bm{x}^{\text{bdy}}_{2}(3)-\bm{x}_{20}(3)-\nu_{2}t}{\mu_{2}}-t\right)^{\!4}$
	$\displaystyle\quad+3\left(\frac{\bm{x}^{\text{bdy}}_{2}(2)-\bm{x}_{20}(2)-t\bm{x}_{20}(3)-\nu_{2}\frac{t^{2}}{2}}{\mu_{2}}-\frac{t^{2}}{2}\right)^{\!2}$
	$\displaystyle\quad-6\left(\frac{\bm{x}^{\text{bdy}}_{2}(1)-\bm{x}_{20}(1)-t\bm{x}_{20}(2)-\frac{t^{2}}{2}\bm{x}_{20}(3)-\nu_{2}\frac{t^{3}}{6}}{\mu_{2}}-\frac{t^{3}}{6}\right)$
	$\displaystyle\qquad\times\left(\frac{\bm{x}^{\text{bdy}}_{2}(3)-\bm{x}_{20}(3)-\nu_{2}t}{\mu_{2}}-t\right),$		(30)

and again, the formula for $p_{2}^{\text{lower}}(\bm{x}^{\text{bdy}}_{2}(1),\bm{x}^{\text{bdy}}_{2}(2),\bm{x}^{\text{bdy}}_{2}(3))$ 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.

IV-C Implicitization of $\partial\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right)$

To derive the implicit equations for the bounding algebraic hypersurfaces $p_{j}^{\text{upper}},p_{j}^{\text{lower}}\in\mathbb{R}\left[x_{1},\ldots,x_{r_{j}}\right]$ for all $j\in[m]$ , we need to eliminate the parameters $\left(s_{1},s_{2},\ldots,s_{r_{j}-1}\right)$ from (26). For this purpose, it is helpful to write (26) succinctly as

\displaystyle\rho_{j,k}^{\pm}=\displaystyle\sum_{q=1}^{r_{j}-1}(-1)^{q+1}\>s_{q}^{r_{j}-k+1},\qquad k\in[r_{j}],

(31)

where

	$\displaystyle\rho_{j,k}^{\pm}:=$	$\displaystyle\frac{(r_{j}-k+1)!}{2\mu_{j}}\bigg{\{}\bm{x}^{\textup{bdy}}_{j}\left(k\right)-\sum_{\ell=1}^{r_{j}}\mathbf{1}_{k\leq\ell}\>\frac{t^{\ell-k}}{(\ell-k)!}\>\bm{x}_{j0}(\ell)\bigg{\}}$
		$\displaystyle-\frac{1}{2}\bigg{\{}\pm(-1)^{r_{j}-1}\>t^{r_{j}-k+1}+\frac{\nu_{j}}{\mu_{j}}\>t^{r_{j}-k+1}\bigg{\}}.$		(32)

To simplify the rather unpleasant notation $\rho_{j,k}^{\pm}$ , 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 $\rho$ ’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 $\pm$ from the $\rho$ ’s in (31). Recall that the plus (resp. minus) superscript in the $\rho$ ’s indicates $p_{j}^{\text{upper}}$ (resp. $p_{j}^{\text{lower}}$ ). From (32), it is clear that in either case, the $\rho_{j,k}$ is affine in $\bm{x}^{\textup{bdy}}_{j}\left(k\right)$ , which is the $k$ th coordinate of the boundary point for the $j$ th block. Importantly, for $k\in[r_{j}]$ , the quantity $\rho_{j,k}$ does not depend on any other component of the boundary point than the $k$ th component. Again, the plus-minus superscripts can be added back post implicitization.

Thus, the notationally simplified version of (31) that suffices for implicitization, is

\displaystyle\rho_{k}=\displaystyle\sum_{q=1}^{d-1}(-1)^{q+1}\>s_{q}^{d-k+1},\qquad k=1,\ldots,d,

(33)

which is a system of $d$ homogeneous polynomials in variables $\left(s_{1},s_{2},\ldots,s_{d-1}\right)$ . The objective is to derive the implicitized polynomial $\wp(\rho_{1},\rho_{2},\ldots,\rho_{d})$ associated with (33).

When $d=2$ , the parameterization (33) becomes

\displaystyle\rho_{1}=s_{1}^{2},\quad\rho_{2}=s_{1},

and we get degree 2 implicitized polynomial

\displaystyle\wp(\rho_{1},\rho_{2})=\rho_{2}^{2}-\rho_{1}=0.

(34)

For $k=1,2$ , substituting for the $\rho_{1},\rho_{2}$ in (34) from (32) with appropriate plus-minus signs recovers (29).

When $d=3$ , the parameterization (33) becomes

\displaystyle\rho_{1}=s_{1}^{3}-s_{2}^{3},\quad\rho_{2}=s_{1}^{2}-s_{2}^{2},\quad\rho_{3}=s_{1}-s_{2},

elementary algebra gives degree 4 implicitized polynomial

\displaystyle\wp(\rho_{1},\rho_{2},\rho_{3})=\rho_{3}^{4}-4\rho_{3}\rho_{1}+3\rho_{2}^{2}=0.

(35)

As before, for $k=1,2,3$ , substituting for the $\rho_{1},\rho_{2},\rho_{3}$ in (35) from (32) with appropriate plus-minus signs recovers (30). However, for $d=4$ or higher, it is practically impossible to derive the implicitization via brute force algebra.

A principled way to implicitize (33) is due to G. Zaimi [38], and starts with defining $\lambda_{k}:=\rho_{d-k+1}$ for $k=1,\ldots,d$ . Introduce the sequence $A_{k}(s_{1},s_{2},\dots,s_{d-1})$ via the generating function (see e.g., [39, Ch. 1])

\displaystyle F(\tau)=\sum_{k\geq 0}A_{k}\tau^{k}=\frac{(1-s_{1}\tau)(1-s_{3}\tau)\cdots}{(1-s_{2}\tau)(1-s_{4}\tau)\cdots}.

(36)

Taking the logarithmic derivative of (36), and then using the generating functions $(1-s_{q}\tau)^{-1}=\sum_{k\geq 0}\left(s_{q}\tau\right)^{k}$ for all $q=1,\ldots,d-1$ , yields

\displaystyle\!\!\frac{F^{\prime}(\tau)}{F(\tau)}=-s_{1}\!\sum_{k\geq 0}\!\left(s_{1}\tau\right)^{k}\!+s_{2}\!\sum_{k\geq 0}\!\left(s_{2}\tau\right)^{k}\!-s_{3}\!\sum_{k\geq 0}\!\left(s_{3}\tau\right)^{k}\!+\ldots.

(37)

Integrating (37) with respect to $\tau$ , we obtain

\displaystyle F(\tau)=\exp\left(-\sum_{k=1}^{d}\frac{\lambda_{k}}{k}\tau^{k}\right).

(38)

Equating (36) and (38) allows us to compute $A_{k}$ as a degree $k$ polynomial of the $\lambda$ ’s.

On the other hand, since the generating function (36) is a rational function with denominator polynomial of degree $\delta:=\lfloor\frac{d-1}{2}\rfloor$ , the following Hankel determinant vanishes^§^§§This result goes back to Kronecker [40]. See also [41, p. 5, Lemma III].

\displaystyle{\mathrm{det}}[A_{d-2\delta+i+j}]_{i,j=0}^{\delta}=0.

(39)

Substituting the $A_{k}$ ’s obtained as degree $k$ polynomials of the $\lambda$ ’s into (39) gives an implicit polynomial in indeterminate $(\lambda_{1},\ldots,\lambda_{d})$ of degree $(\delta+1)(d-\delta)$ . Finally, reverting back the $\lambda$ ’s to the $\rho$ ’s result in the desired implicit polynomial $\wp(\rho_{1},\rho_{2},\ldots,\rho_{d})$ , which is also of degree $(\delta+1)(d-\delta)$ .

For instance, when $d=3$ , the relation (39) becomes

\displaystyle{\mathrm{det}}\left(\begin{bmatrix}A_{1}&A_{2}\\ A_{2}&A_{3}\end{bmatrix}\right)=0.

(40)

In this case, equating (36) and (38) gives

A_{1}=-\lambda_{1},\;A_{2}=\frac{1}{2}\lambda_{1}^{2}-\frac{1}{2}\lambda_{2},\;A_{3}=-\frac{1}{6}\lambda_{1}^{3}+\frac{1}{2}\lambda_{1}\lambda_{2}-\frac{1}{3}\lambda_{3}.

Substituting these back in (40) yields the quartic polynomial $\lambda_{1}^{4}+3\lambda_{2}^{2}-4\lambda_{3}\lambda_{1}=0$ , which under the mapping $(\lambda_{1},\lambda_{2},\lambda_{3})\mapsto(\rho_{3},\rho_{2},\rho_{1})$ recovers (35), and thus (30).

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.

IV-D Dual of $\mathcal{R}\left(\{\bm{x}_{0}\},t\right)$

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 $\mathcal{R}$ or $\mathcal{R}^{\Box}$ . We know in general that $\mathcal{R}^{\circ}$ will be a closed convex set. Depending on the choice of $\bm{x}_{0},\mathcal{U}$ and $t$ , the set $\mathcal{R}\left(\{\bm{x}_{0}\},t\right)$ may not contain the origin, and thus the bipolar

\left(\mathcal{R}\left(\{\bm{x}_{0}\},t\right)\right)^{\circ\circ}={\rm{closure}}\left({\rm{conv}}\left(\mathcal{R}\left(\{\bm{x}_{0}\},t\right)\cup\{\bm{0}\}\right)\right),

that is, we do not have the involution in general.

Furthermore, since $\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right)$ is semialgebraic from Sec. IV-B, so must be its polar dual $\left(\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right)\right)^{\circ}$ ; see e.g., [36, Ch. 5, Sec. 5.2.2].

We also know from Sec. IV-A that $\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right)$ is a zonoid. However, the polar of a zonoid is not a zonoid in general [42, 43], and we should not expect $\left(\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right)\right)^{\circ}$ to be one. Fig. 2 shows $\mathcal{R}\left(\{\bm{x}_{0}\},t\right)$ and $\left(\mathcal{R}\left(\{\bm{x}_{0}\},t\right)\right)^{\circ}$ for the double integrator ( $d=2$ , $m=1$ ).

IV-E Summary of Taxonomy

So far we explained that the compact convex set $\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right)$ 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})\in\mathbb{R}^{2}\mid x_{1}^{4}+x_{2}^{4}\leq 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 $\mathcal{R}^{\Box}$ 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 $\mathcal{R}^{\Box}$ 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 $(\lfloor\frac{d-1}{2}\rfloor+1)(d-\lfloor\frac{d-1}{2}\rfloor)$ . In other words, the $\mathcal{R}^{\Box}$ is not rigidly convex, see [44, Sec. 3.1 and 3.2]. Fig. 3 helps visualize this for $m=1$ . From Fig. 3(a), 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^{\text{upper}},p^{\text{lower}}$ are degree $2$ polynomials. Likewise, Fig. 3(b) 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^{\text{upper}},p^{\text{lower}}$ in this case, are of degree $4$ .

Could the reach set $\mathcal{R}^{\Box}$ 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 $\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right)$ 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$ .

V-A Volume

The following result gives the volume formula for the integrator reach set $\mathcal{R}^{\Box}$ .

Theorem 5.

Fix $\bm{x}_{0}\in\mathbb{R}^{d}$ , let $\mathcal{X}_{0}\equiv\{\bm{x}_{0}\}$ and $\mathcal{U}$ given by (20). Consider the integrator dynamics (1)-(2) with $d$ states, $m$ inputs, and relative degree vector $\bm{r}=(r_{1},r_{2},...,r_{m})^{\top}$ . Define $\mu_{1},\ldots,\mu_{m}$ as in (11)-(12). Then the $d$ -dimensional volume of the integrator reach set (3) at time $t>0$ is

\displaystyle{\mathrm{vol}}\!\left(\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right)\right)\!=\!2^{d}\displaystyle\prod_{j=1}^{m}\!\bigg{\{}\mu_{j}^{r_{j}}t^{r_{j}(r_{j}+1)/2}\!\prod_{k=1}^{r_{j}-1}\!\dfrac{k!}{(2k+1)!}\bigg{\}}.

(41)

For a simple illustration of Theorem 5, consider $d=3$ , $m=2$ with $\bm{r}=(2,1)^{\top}$ . The corresponding reach set $\mathcal{R}^{\Box}(\{\bm{x}_{0}\},t)$ at $t=4$ is shown in Fig. 5 for $\bm{x}_{0}=(1,1,0)^{\top}$ , $\mathcal{U}=[-5\times 5]\times[-3,3]$ . Here $\mu_{1}=5$ and $\mu_{2}=3$ .

This reach set, being a direct product of the double integrator reach set $\mathcal{R}_{1}$ (cf. Fig. 2) and the single integrator reach set $\mathcal{R}_{2}=\{\bm{x}_{0}(3)\}\dotplus[-\mu_{2}t,\mu_{2}t]$ , is a cylinder^¶^¶¶Here, the notation $\bm{x}_{0}(3)$ stands for the third component of vector $\bm{x}_{0}$ .. In [4], we explicitly derived that ${\mathrm{vol}}\left(\mathcal{R}_{1}\right)=\frac{2}{3}\mu_{1}^{2}t^{3}$ , and therefore, the volume of this cylindrical set must be equal to “height of the cylinder $\times$ cross sectional area”, i.e.,

\displaystyle 2\mu_{2}t\times\frac{2}{3}\mu_{1}^{2}t^{3}=\frac{4}{3}\mu_{1}^{2}\mu_{2}t^{4}.

Indeed, a direct application of the formula (41) recovers the above expression.

Remark 2.

If the initial set $\mathcal{X}_{0}$ is not singleton, then computing the volume of the $\mathcal{R}^{\Box}$ requires us to compute the volume of a Minkowski sum. Notice that

	$\displaystyle{\mathrm{vol}}\left(\exp(t\bm{A})\mathcal{X}_{0}\right)$	$\displaystyle=\|{\mathrm{det}}\left(\exp(t\bm{A})\right)\|{\mathrm{vol}}\left(\mathcal{X}_{0}\right)$
		$\displaystyle=\exp\left({\mathrm{trace}}\left(t\bm{A}\right)\right){\mathrm{vol}}\left(\mathcal{X}_{0}\right)$
		$\displaystyle=\exp\left(\!\sum_{j=1}^{m}{\mathrm{trace}}\left(t\bm{A}_{j}\right)\!\!\right){\mathrm{vol}}\left(\mathcal{X}_{0}\right)$
		$\displaystyle={\mathrm{vol}}\left(\mathcal{X}_{0}\right),$

since from (2b), ${\mathrm{trace}}(\bm{A}_{j})=0$ for all $j=1,\ldots,m$ . Therefore, combining (4), (41) with the classical Brunn-Minkowski inequality, we obtain a lower bound for ${\mathrm{vol}}\left(\mathcal{R}^{\Box}\right)$ as

	$\displaystyle\left({\mathrm{vol}}\left(\mathcal{R}^{\Box}\left(\mathcal{X}_{0},t\right)\right)\right)^{1/d}\geq\left({\mathrm{vol}}\left(\mathcal{X}_{0}\right)\right)^{1/d}$
	$\displaystyle\qquad\qquad+2\left(\displaystyle\prod_{j=1}^{m}\bigg{\{}\mu_{j}^{r_{j}}t^{r_{j}(r_{j}+1)/2}\prod_{k=1}^{r_{j}-1}\dfrac{k!}{(2k+1)!}\bigg{\}}\right)^{\!\!1/d}.$

The above bound holds for any compact $\mathcal{X}_{0}\subset\mathbb{R}^{d}$ , not necessarily convex.

V-B Diameter

We now focus on another measure of “size” for the integrator reach set $\mathcal{R}^{\Box}$ , namely its diameter, or maximal width.

By definition, the width [13, p. 42] of $\mathcal{R}^{\Box}\left(\mathcal{X}_{0},t\right)$ , is

\displaystyle w_{\mathcal{R}^{\Box}\left(\mathcal{X}_{0},t\right)}(\bm{\eta}):=h_{\mathcal{R}^{\Box}(\mathcal{X}_{0},t)}\left(\bm{\eta}\right)+h_{\mathcal{R}^{\Box}(\mathcal{X}_{0},t)}\left(-\bm{\eta}\right),

(42)

where $\bm{\eta}\in\mathbb{S}^{d-1}$ (the unit sphere imbedded in $\mathbb{R}^{d}$ ), and the support function $h_{\mathcal{R}^{\Box}(\mathcal{X}_{0},t)}\left(\cdot\right)$ is given by (21). In other words, (42) gives the width of $\mathcal{R}^{\Box}$ in the direction $\bm{\eta}$ .

For singleton $\mathcal{X}_{0}\equiv\{\bm{x}_{0}\}$ , combining (21) and (42), we have

	$\displaystyle w_{\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right)}(\bm{\eta})$	$\displaystyle=\int_{0}^{t}\bigg{\{}\lvert\langle\bm{\eta},\bm{\xi}(s)\rangle\rvert+\lvert\langle-\bm{\eta},\bm{\xi}(s)\rangle\rvert\bigg{\}}\>{\rm{d}}s$
		$\displaystyle=2\int_{0}^{t}\lvert\langle\bm{\eta},\bm{\xi}(s)\rangle\rvert\>{\rm{d}}s,$		(43)

where the last equality follows from the fact that $\bm{\xi}(s)$ in (13) is component-wise nonnegative for all $0\leq s\leq t$ .

The diameter of the reach set $\mathcal{R}^{\Box}$ is its maximal width:

\displaystyle{\rm{diam}}\left(\mathcal{R}^{\Box}\left(\mathcal{X}_{0},t\right)\right):=\underset{\bm{\eta}\in\mathbb{S}^{d-1}}{\max}\>w_{\mathcal{R}^{\Box}\left(\mathcal{X}_{0},t\right)}(\bm{\eta}).

(44)

Notice that (43) is a convex function of $\bm{\eta}$ ; 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).

Theorem 6.

Fix $\bm{x}_{0}\in\mathbb{R}^{d}$ , let $\mathcal{X}_{0}\equiv\{\bm{x}_{0}\}$ and $\mathcal{U}$ given by (20). Consider the integrator dynamics (1)-(2) with $d$ states, $m$ inputs, and relative degree vector $\bm{r}=(r_{1},r_{2},...,r_{m})^{\top}$ . Define $\mu_{1},\ldots,\mu_{m}$ as in (11)-(12). The diameter of the integrator reach set (3) at time $t>0$ is

\displaystyle{\rm{diam}}\!\left(\!\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right)\!\right)\!=\!2\parallel\!\bm{\zeta}(t)\!\parallel_{2}\>=2\!\left(\displaystyle\sum_{j=1}^{m}\mu_{j}^{2}\|\bm{\zeta}_{j}\|^{2}\!\right)^{\!\!\frac{1}{2}}

(45)

wherein $\bm{\zeta}(t)$ is defined as in Sec. II-2, and the $i$ th component of the subvector $\bm{\zeta}_{j}(t)\in\mathbb{R}^{r_{j}}$ is

\displaystyle\!\displaystyle\int_{0}^{t}\!\!\dfrac{s^{(r_{j}-i)}}{(r_{j}-i)!}\>{\rm{d}}s=\dfrac{t^{r_{j}-i+1}}{(r_{j}-i+1)!},\quad i=1,2,\ldots,r_{j}.

(46)

To illustrate Theorem 6, consider the triple integrator with $d=3$ and $m=1$ . In this case, $\mathcal{U}=[\alpha,\beta]$ , $\mu:=(\beta-\alpha)/2$ , and we can parameterize the unit vector $\bm{\eta}\in\mathbb{S}^{2}$ as

\bm{\eta}\equiv\begin{pmatrix}\sin\theta\cos\phi\\ \sin\theta\sin\phi\\ \cos\theta\end{pmatrix}\displaystyle,\quad\theta\in[0,\pi],\quad\phi\in[0,2\pi).

Thus (44) reduces to

2\mu\>\underset{\begin{subarray}{c}\theta\in[0,\pi]\\ \phi\in[0,2\pi)\end{subarray}}{\max}\>\displaystyle\int_{0}^{t}\lvert s^{2}\left(\sin\theta\cos\phi\right)/2+s\sin\theta\sin\phi+\cos\theta\rvert\>{\rm{d}}s.

Furthermore, $\bm{\zeta}(t)=(t^{3}/6,t^{2}/2,t)^{\top}$ , and we obtain

\displaystyle\bm{\eta}^{\max}=\begin{pmatrix}\!\sin\theta^{\max}\sin\phi^{\max}\!\\ \!\sin\theta^{\max}\cos\phi^{\max}\!\\ \!\cos\theta^{\max}\!\end{pmatrix}=\dfrac{\pm 1}{\sqrt{t^{4}+9t^{2}+36}}\begin{pmatrix}t^{2}\\ 3t\\ 6\end{pmatrix},

where $\pm$ means that either all components are plus or all minus. Thus, the maximizing tuples $\left(\phi^{\max},\theta^{\max}\right)\in[0,\pi]\times[0,2\pi)$ are given by

	$\displaystyle\left(\phi^{\max},\theta^{\max}\right)$
	$\displaystyle=\begin{cases}\left(\arctan\left(3/t\right),\arccos\left(6/\sqrt{t^{4}+9t^{2}+36}\right)\right),\\ \left(\pi+\arctan\left(3/t\right),\arccos\left(-6/\sqrt{t^{4}+9t^{2}+36}\right)\right).\end{cases}$		(47)

Hence, the diameter of the triple integrator reach set at time $t$ is equal to $\left(\mu t/3\right)\sqrt{t^{4}+9t^{2}+36}$ .

Fig. 6 shows how the width of the integrator reach set for $d=3$ , $m=1$ varies over $(\phi,\theta)\in[0,\pi]\times[0,2\pi)$ , which parameterize the unit sphere $\mathbb{S}^{2}$ . The location of the maximizers are given by (47), and are depicted in Fig. 6 via filled black circle and filled black square.

For a visualization of the width and diameter for the double integrator, see [4, Fig. 2].

V-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 $\mathcal{R}$ and $\mathcal{R}^{\Box}$ .

V-C1 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, $\mathcal{U}=[\alpha,\beta]$ , and therefore $\mu:=(\beta-\alpha)/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 asymptotic [51] for $\prod_{k=1}^{d-1}(2k+1)!/k!$ , we find the $d\rightarrow\infty$ asymptotic for the volume:

{\mathrm{vol}}\left(\mathcal{R}_{d}\left(\{\bm{x}_{0}\},t\right)\right)\sim(2\mu)^{d}t^{d(d+1)/2}\dfrac{\exp\left(\frac{3}{2}d^{2}+\frac{1}{12}\right)}{c\times 2^{\left(2d^{2}-\frac{1}{12}\right)}\>d^{\left(d^{2}+\frac{1}{12}\right)}},

where $c\approx 1.2824\ldots$ 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^{\prime}$ with $d>d^{\prime}$ , and all other parameters kept fixed, there exists a critical time $t_{\text{cr}}$ when the volume of the $d$ dimensional reach set overtakes that of the $d^{\prime}$ dimensional reach set.

For any $d>d^{\prime}$ , the critical time $t_{\text{cr}}$ satisfies

\displaystyle\underbrace{{\mathrm{vol}}\left(\mathcal{R}_{d}\left(\{\bm{x}_{0}\},t_{\text{cr}}\right)\right)}_{d\;\text{dimensional volume}}=\underbrace{{\mathrm{vol}}\left(\mathcal{R}_{d^{\prime}}\left(\{\bm{x}_{0}\},t_{\text{cr}}\right)\right)}_{d^{\prime}\;\text{dimensional volume}},

which together with (41) yields

\displaystyle t_{\text{cr}}=\left(2\mu\right)^{-\frac{2}{d+d^{\prime}+1}}\left(\prod_{k=d^{\prime}}^{d-1}\dfrac{\left(2k+1\right)!}{k!}\right)^{\frac{2}{(d-d^{\prime})(d+d^{\prime}+1)}}.

(48)

In particular, for $d^{\prime}=d-1$ , we get

\displaystyle t_{\text{cr}}=\left(\frac{1}{2\mu}\frac{(2d-1)!}{(d-1)!}\right)^{1/d},\quad d=2,3,\ldots.

(49)

For instance, when $\mu=1$ , $d=3$ , $d^{\prime}=2$ , we have $t_{\text{cr}}=(30)^{1/3}\approx 3.1072$ . When $\mu=1$ , $d=4$ , $d^{\prime}=3$ , we have $t_{\text{cr}}=(420)^{1/4}\approx 4.5270$ . The dashed vertical lines in Fig. 7 show the critical times given by (49).

Applying Stirling’s approximation $n!\sim\sqrt{2\pi n}(n/e)^{n}$ , we obtain the $d\rightarrow\infty$ asymptotic for (49):

\displaystyle t_{\text{cr}}\sim\frac{4}{e}\>d\>\mu^{-\frac{1}{d}}\>2^{-\frac{3}{2d}},

where $\sim$ denotes asymptotic equivalence [53, Ch. 1.4], and $e$ is the Euler number.

V-C2 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, $\mathcal{U}=[\alpha,\beta]$ , $\mu:=(\beta-\alpha)/2$ . As expected, the diameter of the reach set increases with time for any fixed $d$ .

As $d\rightarrow\infty$ , 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

\displaystyle\displaystyle\lim_{d\rightarrow\infty}{\rm{diam}}\left(\mathcal{R}\left(\{\bm{x}_{0}\},t\right)\right)=\lim_{d\to\infty}2\mu\sqrt{\sum_{j=1}^{d}\left(\dfrac{t^{j}}{j!}\right)^{2}}.

(50)

We write the partial sum

\displaystyle\sum_{j=1}^{d}\left(\dfrac{t^{j}}{j!}\right)^{2}=\underbrace{\sum_{j=1}^{\infty}\left(\dfrac{t^{j}}{j!}\right)^{2}}_{=:S_{1}}-\underbrace{\sum_{j=d+1}^{\infty}\left(\dfrac{t^{j}}{j!}\right)^{2}}_{=:S_{2}},

(51)

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}(\cdot)$ is the zeroth order modified Bessel function of the first kind. This follows from the very definition of the $\nu$ th order modified Bessel function of the first kind, given by

I_{\nu}(z):=\left(z/2\right)^{\nu}\displaystyle\sum_{j=0}^{\infty}\frac{\left(z^{2}/4\right)^{j}}{j!\>\Gamma\left(\nu+j+1\right)},\quad\nu\in\mathbb{R},

where $\Gamma(\cdot)$ denotes the Gamma function.

On the other hand, using the definition of the generalized hypergeometric function^∥^∥∥Here, $(\cdot)_{n}$ denotes the Pochhammer symbol [54, p. 256] or rising factorial.

{}_{1}F_{2}\left(a_{1};b_{1},b_{2};z\right):=\sum_{n=0}^{\infty}\dfrac{(a_{1})_{n}}{(b_{1})_{n}(b_{2})_{n}}\dfrac{z^{n}}{n!},

we find that

S_{2}=\dfrac{t^{2(d+1)}\>_{1}F_{2}\left(1;d+2,d+2;t^{2}\right)}{\left((d+1)!\right)^{2}}.

Therefore, (51) evaluates to

\displaystyle\!\!\!S_{1}\!-\!S_{2}=I_{0}(2t)-1-\dfrac{t^{2(d+1)}\>_{1}F_{2}\left(1;d+2,d+2;t^{2}\right)}{\left((d+1)!\right)^{2}}\!.

(52)

Combining (50), (51), (52), and using the continuity of the square root function on $[0,\infty)$ , we deduce that

	$\displaystyle\!\!\!\!\!\displaystyle\lim_{d\rightarrow\infty}{\rm{diam}}\left(\mathcal{R}\left(\{\bm{x}_{0}\},t\right)\right)$	$\displaystyle=2\mu\sqrt{\lim_{d\to\infty}\left(S_{1}-S_{2}\right)}$
		$\displaystyle=2\mu\sqrt{I_{0}(2t)-1}.$		(53)

That $\lim_{d\rightarrow\infty}S_{2}$ exists and equals to zero, follows from (51) and the continuity of the square:

\lim_{d\rightarrow\infty}S_{2}=\lim_{j\rightarrow\infty}\left(\frac{t^{j}}{j!}\right)^{\!\!2}=\left(\lim_{j\rightarrow\infty}\frac{t^{j}}{j!}\right)^{\!\!2}=0.

To see the last equality, let $a_{j}:=t^{j}/j!$ . By the ratio test, $\underset{j\rightarrow\infty}{\lim\sup}|a_{j+1}/a_{j}|=\underset{j\rightarrow\infty}{\lim}t/j=0<1$ , hence $\{a_{j}\}$ is a Cauchy sequence and $\underset{j\rightarrow\infty}{\lim}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 [3, 2] 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, 56, 57, 58, 59, 60, 61, 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 $\mathcal{R}_{\text{approx}}(t)$ of the single input integrator reach sets $\mathcal{R}(\{\bm{0}\},t)\subseteq\mathcal{R}_{\text{approx}}(\{\bm{0}\},t)$ for $d=2,3,4$ with $0\leq t\leq 1$ and $\mu=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 ${\mathrm{vol}}(\mathcal{R})/{\mathrm{vol}}(\mathcal{R}_{\text{approx}})\in[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 $\leq 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 over-approximations $\mathcal{R}_{\text{approx}}(t)$ of the single input integrator reach sets $\mathcal{R}(\{\bm{0}\},t)\subseteq\mathcal{R}_{\text{approx}}(\{\bm{0}\},t)$ for $d=2,3,4$ with $0\leq t\leq 1$ and $\mu=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 $\mathcal{E}\left(\bm{q}(t),\bm{Q}_{\bm{\ell}_{i}(t)}(t)\right)$ defined as

\displaystyle\{\bm{x}\in\mathbb{R}^{d}\mid\left(\bm{x}-\bm{q}(t)\right)\left(\bm{Q}_{\bm{\ell}_{i}(t)}(t)\right)^{-1}\left(\bm{x}-\bm{q}(t)\right)^{\top}\leq 1\},

for unit vectors $\bm{\ell}_{i}(0)\in\mathbb{R}^{d}$ , $i=1,\ldots,N$ . The choice of $\bm{\ell}_{i}(0)$ determines $\bm{\ell}_{i}(t):=\exp(-\bm{A}^{\top}t)\bm{\ell}_{i}(0)$ , which in turn parameterizes the $d\times d$ symmetric positive definite shape matrix $\bm{Q}_{\bm{\ell}_{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 $\bm{q}(t)\in\mathbb{R}^{d}$ , and the shape matrices $\bm{Q}_{\bm{\ell}_{i}(t)}(t)$ for this parameterized family of ellipsoids are constructed such that $\cap_{i=1}^{N}\mathcal{E}\left(\bm{q}(t),\bm{Q}_{\bm{\ell}_{i}(t)}(t)\right)$ is guaranteed to be a superset of the reach set at time $t$ for any finite $N$ , and for $N\rightarrow\infty$ , recovers the reach set at that time.

For the results shown in Fig. 10, we used $N=20$ randomly chosen unit vectors $\bm{\ell}_{i}(0)\in\mathbb{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] of the convex set $\cap_{i=1}^{20}\mathcal{E}\left(\bm{q}(t),\bm{Q}_{\bm{\ell}_{i}(t)}(t)\right)$ , 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 $\cap_{i=1}^{20}\mathcal{E}\left(\bm{q}(t),\bm{Q}_{\bm{\ell}_{i}(t)}(t)\right)$ . 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 $\mathcal{R}^{\Box}$ 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

	$\displaystyle h_{\mathcal{R}(\mathcal{X}_{0},t)}\left(\bm{y}\right)$	$\displaystyle=\underset{\bm{x}_{0}\in\mathcal{X}_{0}}{\sup}\langle\bm{y},\exp\left(t\bm{A}\right)\bm{x}_{0}\rangle$
		$\displaystyle\quad\qquad+h_{\int_{0}^{t}\exp(s\bm{A})\bm{B}\mathcal{U}{\rm{d}}s}(\bm{y}).$		(54)

The block diagonal structure of the matrix $\bm{A}$ in (2) implies

\displaystyle\!\underset{\bm{x}_{0}\in\mathcal{X}_{0}}{\sup}\langle\bm{y},\exp\left(t\bm{A}\right)\bm{x}_{0}\rangle\!=\!\underset{\bm{x}_{0}\in\mathcal{X}_{0}}{\sup}\sum_{j=1}^{m}\langle\bm{y}_{j},\exp\left(t\bm{A}_{j}\right)\bm{x}_{j0}\rangle.

(55)

Following the definition of support function and [4, Proposition 1], we then have

		$\displaystyle h_{\int_{0}^{t}\exp(s\bm{A})\bm{B}\>\mathcal{U}\>{\rm{d}}s}(\bm{y})=\int_{0}^{t}\!\!h_{\exp(s\bm{A})\bm{B}\>\mathcal{U}}\>(\bm{y})\;{\rm{d}}{s}$
		$\displaystyle\quad=\int_{0}^{t}\underset{\bm{u}\in\mathcal{U}}{\sup}~{}\langle\bm{y},\exp(s\bm{A})\bm{B}\bm{u}\rangle~{}{\rm{d}}s$
		$\displaystyle\quad=\int^{t}_{0}\underset{\bm{u}\in{\rm{closure}}\left({\rm{conv}}\left(\mathcal{U}\right)\right)}{\sup}~{}\sum_{j=1}^{m}~{}\{\langle\bm{y}_{j},\bm{\xi}_{j}(s)\rangle u_{j}\}\>{\rm{d}}{s}.$		(56)

The last equality in (-A) 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 $\mathcal{U}$ . Substituting (55) and (-A) in (54) yields (16). $\blacksquare$

-B Proof of Theorem 1

Since the uncertainties in (20) along different input co-ordinate axes are mutually independent, the support function of the reach set is of the form (19). Therefore, in this case, (16) takes the form

	$\displaystyle h_{\mathcal{R}(\mathcal{X}_{0},t)}\left(\bm{y}\right)=$	$\displaystyle\sum_{j=1}^{m}\bigg{\{}\underset{\bm{x}_{j0}\in\mathcal{X}_{j0}}{\sup}\langle\bm{y}_{j},\exp\left(t\bm{A}_{j}\right)\bm{x}_{j0}\rangle$
		$\displaystyle+\int^{t}_{0}\underset{u_{j}\in\mathcal{[}\alpha_{j},\beta_{j}]}{\sup}\langle\bm{y}_{j},\bm{\xi}_{j}(s)\rangle\>u_{j}\>{\rm{d}}{s}\bigg{\}}.$		(57)

The optimizer $u_{j}^{\text{opt}}$ of the integrand in the RHS of (57), for $j\in[m]$ , can be written in terms of the Heaviside unit step function $H(\cdot)$ as

	$\displaystyle u_{j}^{\text{opt}}$	$\displaystyle=\alpha_{j}+(\beta_{j}-\alpha_{j})H(\langle\bm{y}_{j},\bm{\xi}_{j}\rangle)$
		$\displaystyle=\alpha_{j}+(\beta_{j}-\alpha_{j})\times\frac{1}{2}\left(1+{\rm{sgn}}\left(\langle\bm{y}_{j},\bm{\xi}_{j}\rangle\right)\right),$

where ${\rm{sgn}}(\cdot)$ denotes the signum function. Therefore,

\displaystyle\underset{u_{j}\in\mathcal{[}\alpha_{j},\beta_{j}]}{\sup}\!\!\!\!\langle\bm{y}_{j},\bm{\xi}_{j}(s)\rangle u_{j}=\nu_{j}\langle\bm{y}_{j},\bm{\xi}_{j}(s)\rangle+\mu_{j}|\langle\bm{y}_{j},\bm{\xi}_{j}(s)\rangle|

(58)

for $0\leq s\leq t$ . Substituting (58) back in (57) and integrating over $s$ completes the proof. $\blacksquare$

-C Proof of Theorem 2

For $s\in[0,t]$ , let the vector measure $\widetilde{\bm{\mu}}$ be defined as ${\rm{d}}\widetilde{\bm{\mu}}(s):=\bm{\xi}(s){\rm{d}}s$ where $\bm{\xi}(s)$ is given by (13). Then $\int_{0}^{t}|\langle\bm{y},\bm{\xi}(s)\rangle|{\rm{d}}s$ 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.

From (9) and (21), $\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right)$ is the translation of a set with support function $\int_{0}^{t}|\langle\bm{y},\bm{\xi}(s)\rangle|{\rm{d}}s$ , i.e., the translation of a zonoid. Thus, $\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right)$ is a zonoid. $\blacksquare$

-D Proof of Proposition 2

From Sec. II-2, the supporting hyperplane at any $\bm{x}^{\textup{bdy}}\in\partial\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right)$ with outward normal $\bm{y}\in\mathbb{R}^{d}$ is $\langle\bm{y},\bm{x}^{\textup{bdy}}\rangle=h_{\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right)}(\bm{y})$ , and the Legendre-Fenchel conjugate

\displaystyle h_{\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right)}^{*}\left(\bm{x}^{\textup{bdy}}\right)=0.

(59)

For $j\in[m]$ , let $\bm{y}$ comprise of subvectors $\bm{y}_{j}\in\mathbb{R}^{r_{j}}$ . Since the Cartesian product (18) is equivalent to the Minkowski sum $\mathcal{R}_{1}\dotplus\ldots\dotplus\mathcal{R}_{m}$ , and the support function of Minkowski sum is the sum of support functions of the summand sets [13, p. 48], we have

	$\displaystyle h_{\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right)}(\bm{y})$	$\displaystyle=\sum_{j=1}^{m}h_{\mathcal{R}_{j}\left(\{\bm{x}_{0}\},t\right)}(\bm{y}_{j})$
	$\displaystyle\Rightarrow h_{\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right)}^{*}\left(\bm{x}^{\textup{bdy}}\right)$	$\displaystyle=\sum_{j=1}^{m}h_{\mathcal{R}_{j}\left(\{\bm{x}_{0}\},t\right)}^{*}\left(\bm{x}^{\textup{bdy}}_{j}\right),$		(60)

wherein the last line follows from the property that the Legendre-Fenchel conjugate of a separable sum equals to the sum of the Legendre-Fenchel conjugates [46, p. 95].

Therefore, combining (59) and (60), we obtain

	$\displaystyle\sum_{j=1}^{m}\>\underset{\bm{y}_{j}\in\mathbb{R}^{r_{j}}}{\inf}\!\bigg{\{}\!\langle-\bm{x}^{\textup{bdy}}_{j}+\exp\left(t\bm{A}_{j}\right)\bm{x}_{j0}+\nu_{j}\bm{\zeta}_{j}(t),\bm{y}_{j}\rangle$
	$\displaystyle\qquad\qquad\qquad\qquad\qquad+\mu_{j}\!\!\int_{0}^{t}\!\!\|\langle\bm{y}_{j},\bm{\xi}_{j}(s)\rangle\|\>{\rm{d}}s\!\bigg{\}}=0.$		(61)

For $j\in[m]$ , since each objective in (61) involves an integral of the absolute value of a polynomial in $s$ of degree $r_{j}-1$ , that polynomial can have at most $r_{j}-1$ roots in the interval $[0,t]$ , i.e., can have at most $r_{j}-1$ sign changes in that interval. If all $r_{j}-1$ roots of the aforesaid polynomial are in $[0,t]$ , we denote these roots as $s_{1}\leq s_{2}\leq\ldots\leq s_{r_{j}-1}$ , and write

	$\displaystyle\int_{0}^{t}\!\!\|\langle\bm{y}_{j},\bm{\xi}_{j}(s)\rangle\|\>{\rm{d}}s=\pm\int_{0}^{s_{1}}\!\!\langle\bm{y}_{j},\bm{\xi}_{j}(s)\rangle\>{\rm{d}}s\mp\int_{s_{1}}^{s_{2}}\!\!\langle\bm{y}_{j},\bm{\xi}_{j}(s)\rangle\>{\rm{d}}s$
	$\displaystyle\qquad\qquad\qquad\qquad\quad\pm\ldots\pm(-1)^{r_{j}-1}\int_{s_{r_{j}-1}}^{t}\!\!\langle\bm{y}_{j},\bm{\xi}_{j}(s)\rangle\>{\rm{d}}s$
	$\displaystyle=\langle\bm{y}_{j},\pm\bm{\zeta}_{j}(0,s_{1})\mp\bm{\zeta}_{j}(s_{1},s_{2})\pm\ldots\pm(-1)^{r_{j}-1}\bm{\zeta}_{j}(s_{r_{j}-1},t)\rangle.$		(62)

Notice that even if the number of roots in $[0,t]$ is strictly less than^**^****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. $r_{j}-1$ , the expression (62) is generic in the sense the corresponding summand integrals become zero. Thus, combining (61) and (62), we arrive at

	$\displaystyle\sum_{j=1}^{m}\>\underset{\bm{y}_{j}\in\mathbb{R}^{r_{j}}}{\inf}\langle-\bm{x}^{\textup{bdy}}_{j}+\exp\left(t\bm{A}_{j}\right)\bm{x}_{j0}+\nu_{j}\bm{\zeta}_{j}(t)\pm\mu_{j}\bm{\zeta}_{j}(0,s_{1})$
	$\displaystyle\mp\mu_{j}\bm{\zeta}_{j}(s_{1},s_{2})\pm\ldots\pm\mu_{j}(-1)^{r_{j}-1}\bm{\zeta}_{j}(s_{r_{j}-1},t),\bm{y}_{j}\rangle=0.$		(63)

The left hand side of (63) being the sum of the infimum values of linear functions, can achieve zero if and only if each of those infimum equals to zero, i.e., if and only if

	$\displaystyle\bm{x}^{\textup{bdy}}_{j}=\exp\left(t\bm{A}_{j}\right)\bm{x}_{j0}+\nu_{j}\bm{\zeta}_{j}(t)\pm\mu_{j}\bm{\zeta}_{j}(0,s_{1})\mp\mu_{j}\bm{\zeta}_{j}(s_{1},s_{2})$
	$\displaystyle\qquad\quad\pm\ldots\pm(-1)^{r_{j}-1}\mu_{j}\bm{\zeta}_{j}(s_{r_{j}-1},t).$		(64)

Using (6), (13) and (14), we simplify (64) to (26), thereby completing the proof. $\blacksquare$

-E Proof of Corollary 3

From (26), we get two different parametric representations of $\bm{x}_{j}^{\text{bdy}}$ in terms of $(s_{1},s_{2},\ldots,s_{r_{j}-1})$ . One parametric representation results from the choice of positive sign for the $\pm$ 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_{j}^{\textup{upper}}(\bm{x})=0$ (resp. $p_{j}^{\textup{lower}}(\bm{x})=0$ ), the result follows. $\blacksquare$

-F Proof of Theorem 4

We notice that (26) gives polynomial parameterizations of the components of $\bm{x}_{j}^{\text{bdy}}$ for all $j\in[m]$ . In particular, for each $k\in[r_{j}]$ , the right hand side of (26) is a homogeneous polynomial in $r_{j}-1$ parameters $(s_{1},s_{2},\ldots,s_{r_{j}-1})$ of degree $r_{j}-k+1$ . By polynomial implicitization [25, p. 134], the corresponding implicit equations $p_{j}^{\text{upper}}\left(\bm{x}_{j}^{\text{bdy}}\right)=0$ (when fixing plus sign for $\pm$ in (26)) and $p_{j}^{\text{lower}}\left(\bm{x}_{j}^{\text{bdy}}\right)=0$ (when fixing minus sign for $\pm$ in (26)), must define affine varieties $V_{\mathbb{R}[x_{1},...,x_{r_{j}}]}(p_{j}^{\text{upper}}),V_{\mathbb{R}[x_{1},...,x_{r_{j}}]}(p_{j}^{\text{lower}})$ in $\mathbb{R}\left[x_{1},\ldots,x_{d}\right]$ .

Specifically, denote the right hand sides of (26) as $g_{1}^{\pm},\ldots,g_{r_{j}}^{\pm}$ for all $j\in[m]$ , where the superscripts indicate that either all $g$ ’s are chosen with plus signs, or all with minus signs. Then write (26) as

	$\displaystyle\bm{x}^{\text{bdy}}_{j}(1)$	$\displaystyle=g_{1}^{\pm}\left(s_{1},s_{2},\ldots,s_{r_{j}-1}\right),$
		$\displaystyle\;\;\vdots$
	$\displaystyle\bm{x}^{\text{bdy}}_{j}(r_{j})$	$\displaystyle=g_{r_{j}}^{\pm}\left(s_{1},s_{2},\ldots,s_{r_{j}-1}\right).$

Now for each $j\in[m]$ , consider the ideal

	$\displaystyle I_{j}^{\pm}:=$	$\displaystyle\langle\langle\bm{x}^{\text{bdy}}_{j}(1)-g_{1}^{\pm},\bm{x}^{\text{bdy}}_{j}(2)-g_{2}^{\pm},\ldots,\bm{x}^{\text{bdy}}_{j}(r_{j})-g_{r_{j}}^{\pm}\rangle\rangle$
		$\displaystyle\subseteq\mathbb{R}[s_{1},s_{2},\ldots,s_{r_{j}-1},x_{1},x_{2},\ldots,x_{r_{j}}],$

and let $I_{j,r_{j}-1}^{\pm}:=I_{j}^{\pm}\cap\mathbb{R}[x_{1},...,x_{r_{j}}]$ be the $(r_{j}-1)$ th elimination ideal of $I_{j}^{\pm}$ . Then for each $j\in[m]$ , the variety

V\left(I_{j,r_{j}-1}^{+}\right)=V_{\mathbb{R}[x_{1},...,x_{r_{j}}]}(p_{j}^{\text{upper}}).

Likewise, the variety $V\left(I_{j,r_{j}-1}^{-}\right)=V_{\mathbb{R}[x_{1},...,x_{r_{j}}]}(p_{j}^{\text{lower}})$ .

Thus, the algebraic boundary (i.e., the Zariski closure of the Euclidean boundary) of $\mathcal{R}_{j}$ is

\partial\mathcal{R}_{j}=V_{\mathbb{R}[x_{1},...,x_{r_{j}}]}\left(p_{j}^{\text{upper}}\right)\cup V_{\mathbb{R}[x_{1},...,x_{r_{j}}]}\left(p_{j}^{\text{lower}}\right).

Therefore, $\mathcal{R}_{j}:=\{\bm{x}\in\mathbb{R}^{r_{j}}\mid p_{j}^{\textup{upper}}(\bm{x})\leq 0,\;p_{j}^{\textup{lower}}(\bm{x})\leq 0\}$ is semialgebraic for all $j\in[m]$ .

Since the Cartesian product of semialgebraic sets is semialgebraic, the statement follows from (18). $\blacksquare$

-G Proof of Theorem 5

We organize the proof in three steps.
Step 1: From (18), we have

	$\displaystyle{\mathrm{vol}}\left(\mathcal{R}^{\Box}\left(\{\bm{x}_{0}\},t\right)\right)$	$\displaystyle={\mathrm{vol}}\left(\mathcal{R}_{1}\times\mathcal{R}_{2}\times\ldots\times\mathcal{R}_{m}\right)$
		$\displaystyle=\displaystyle\prod_{j=1}^{m}{\mathrm{vol}}\left(\mathcal{R}_{j}\left(\{\bm{x}_{0}\},t\right)\right).$		(65)

Step 2: Motivated by (65), we focus on deriving the $r_{j}$ -dimensional volume of $\mathcal{R}_{j}\left(\{\bm{x}_{0}\},t\right)$ . For this purpose, we proceed as in [4] by uniformly discretizing the interval $[0,t]$ into $n$ subintervals $[(i-1)t/n,it/n)$ , $i=1,\ldots,n$ , with $(n+1)$ breakpoints $\{t_{i}\}_{i=0}^{n}$ , where $t_{i}:=it/n$ for $i=0,1,\ldots,n$ .

From (24), we then have

	$\displaystyle{\mathrm{vol}}\left(\mathcal{R}_{j}\left(\{\bm{x}_{0}\},t\right)\right)={\mathrm{vol}}\left(\displaystyle\lim_{n\rightarrow\infty}\displaystyle\sum_{i=0}^{n}\frac{t}{n}\exp\left(t_{i}\bm{A}_{j}\right)\bm{b}_{j}[-\mu_{j},\mu_{j}]\right)$
	$\displaystyle=\!\displaystyle\lim_{n\rightarrow\infty}\!\left(\dfrac{t}{n}\right)^{\!\!r_{j}}\!\!{\mathrm{vol}}\left(\displaystyle\sum_{i=0}^{n}\exp\left(t_{i}\bm{A}_{j}\right)\bm{b}_{j}[-\mu_{j},\mu_{j}]\right)$
	$\displaystyle=t^{r_{j}}\displaystyle\lim_{n\rightarrow\infty}\dfrac{1}{n^{r_{j}}}{\mathrm{vol}}\left(\displaystyle\sum_{i=0}^{n}\mu_{j}\bm{\xi}_{j}(t_{i})[-1,1]\right),$		(66)

where $\bm{\xi}_{j}$ was defined in (13). We recognize that the set $\sum_{i=0}^{n}\mu_{j}\bm{\xi}_{j}(t_{i})[-1,1]$ in (66) is a Minkowski sum of $n+1$ intervals, i.e., is a zonotope imbedded in $\mathbb{R}^{r_{j}}$ , wherein each interval is rotated and scaled in $\mathbb{R}^{r_{j}}$ via different linear transformations $\exp(t_{i}\bm{A}_{j})$ , $i=0,1,\ldots,n$ .

Using the formula for the volume of zonotopes [15, eqn. (57)], [69, exercise 7.19], we can write (66) as

	$\displaystyle\mathrm{vol}$	$\displaystyle\left(\mathcal{R}_{j}\left(\{\bm{x}_{0}\},t\right)\right)=(2\mu_{j}t)^{r_{j}}\displaystyle\lim_{n\rightarrow\infty}\dfrac{1}{n^{r_{j}}}\>\times$
		$\displaystyle\displaystyle\sum_{0\leq i_{1}<i_{2}<\ldots<i_{r_{j}}\leq n}\!\!\!\!{\mathrm{det}}\left(\bm{\xi}_{j}(t_{i_{1}})\|\bm{\xi}_{j}(t_{i_{2}})\|\ldots\|\bm{\xi}_{j}(t_{i_{r_{j}}})\right).$		(67)

To compute the summand determinants in (67), let

\Delta_{j}\left(i_{1},i_{2},\ldots,i_{r_{j}}\right):={\mathrm{det}}\left(\bm{\xi}_{j}(t_{i_{1}})|\bm{\xi}_{j}(t_{i_{2}})|\ldots|\bm{\xi}_{j}(t_{i_{r_{j}}})\right),

where $0\leq i_{1}<i_{2}<\ldots<i_{r_{j}}\leq n$ . In the matrix list notation, let us use the vertical bars $|\cdot|$ to denote the absolute value of determinant. From (13), $\Delta_{j}\left(i_{1},i_{2},\ldots,i_{r_{j}}\right)$ equals

	$\displaystyle\begin{vmatrix}\dfrac{\left(i_{1}t/n\right)^{r_{j}-1}}{(r_{j}-1)!}&\dfrac{\left(i_{2}t/n\right)^{r_{j}-1}}{(r_{j}-1)!}&\ldots&\dfrac{\left(i_{r_{j}}t/n\right)^{r_{j}-1}}{(r_{j}-1)!}\\ &&&\\ \dfrac{\left(i_{1}t/n\right)^{r_{j}-2}}{(r_{j}-2)!}&\dfrac{\left(i_{2}t/n\right)^{r_{j}-2}}{(r_{j}-2)!}&\ldots&\dfrac{\left(i_{r_{j}}t/n\right)^{r_{j}-2}}{(r_{j}-2)!}\\ &&&\\ \vdots&\vdots&\vdots&\vdots\\ &&&\\ i_{1}t/n&i_{2}t/n&\ldots&i_{r_{j}}t/n\\ &&&\\ 1&1&\ldots&1\end{vmatrix}$
	$\displaystyle=\dfrac{(t/n)^{1+2+\ldots+(r_{j}-1)}}{1!\times 2!\times\ldots\times(r_{j}-1)!}\begin{vmatrix}i_{1}^{r_{j}-1}&i_{2}^{r_{j}-1}&\ldots&i_{r_{j}}^{r_{j}-1}\\ &&&\\ i_{1}^{r_{j}-2}&i_{2}^{r_{j}-2}&\ldots&i_{r_{j}}^{r_{j}-2}\\ &&&\\ \vdots&\vdots&\vdots&\vdots\\ &&&\\ i_{1}&i_{2}&\ldots&i_{r_{j}}\\ &&&\\ 1&1&\ldots&1\end{vmatrix},$
	$\displaystyle=\dfrac{(t/n)^{r_{j}(r_{j}-1)/2}}{\displaystyle\prod_{k=1}^{r_{j}-1}k!}\begin{vmatrix}1&1&\ldots&1\\ &&&\\ i_{1}&i_{2}&\ldots&i_{r_{j}}\\ &&&\\ \vdots&\vdots&\vdots&\vdots\\ &&&\\ i_{1}^{r_{j}-2}&i_{2}^{r_{j}-2}&\ldots&i_{r_{j}}^{r_{j}-2}\\ &&&\\ i_{1}^{r_{j}-1}&i_{2}^{r_{j}-1}&\ldots&i_{r_{j}}^{r_{j}-1}\\ \end{vmatrix},$		(68)

where we used the properties of elementary row operations.

Notice that the determinant appearing in the last step of (68) is the Vandermonde determinant (see e.g., [70, p. 37])

\displaystyle\displaystyle\prod_{1\leq a<b\leq r_{j}}\left(i_{b}-i_{a}\right).

(69)

Combining (67), (68) and (69), we obtain

	$\displaystyle\mathrm{vol}$	$\displaystyle\left(\mathcal{R}_{j}\left(\{\bm{x}_{0}\},t\right)\right)=\dfrac{(2\mu_{j})^{r_{j}}t^{r_{j}(r_{j}+1)/2}}{\displaystyle\prod_{k=1}^{r_{j}-1}k!}\displaystyle\lim_{n\rightarrow\infty}\dfrac{1}{n^{r_{j}(r_{j}+1)/2}}\>\times$
		$\displaystyle\displaystyle\sum_{0\leq i_{1}<i_{2}<\ldots<i_{r_{j}}\leq n}\quad\displaystyle\prod_{1\leq a<b\leq r_{j}}\left(i_{b}-i_{a}\right).$		(70)

Step 3: Our next task is to simplify (70). Observe that the sum

\displaystyle\displaystyle\sum_{0\leq i_{1}<i_{2}<\ldots<i_{r_{j}}\leq n}\quad\displaystyle\prod_{1\leq a<b\leq r_{j}}\left(i_{b}-i_{a}\right),

(71)

returns a polynomial in $n$ of degree $r_{j}(r_{j}+1)/2$ , and hence the limit in (70) is always well-defined. Specifically, the limit extracts the leading coefficient of this polynomial.

Let us denote the leading coefficient of the sum (71) as $c(r_{j})$ . By the Euler-Maclaurin formula [71], [72, Chap. II.10]:

\displaystyle\!\!\!\!c(r_{j})=\!\!\!\displaystyle\int\limits_{0\leq y_{1}<y_{2}<\ldots<y_{r_{j}}\leq 1}\;\prod\limits_{1\leq\alpha<\beta\leq r_{j}}\!\!(y_{a}-y_{b})\cdot\prod\limits_{a=1}^{r_{j}}{\rm{d}}y_{a}.

(72)

One way to unpack (72) is to write it as a sum over the symmetric permutation group $\mathfrak{S}_{r_{j}}$ of the finite set $[r_{j}]$ , i.e.,

\!\!c(r_{j})=\displaystyle\sum_{\sigma\in\mathfrak{S}_{r_{j}}}\text{sgn}(\sigma)\displaystyle\frac{1}{\prod_{k=1}^{r_{j}}(\sigma_{1}+\sigma_{2}+...+\sigma_{k})},

where $\text{sgn}(\sigma):=(-1)^{\nu}$ , $\nu:=\{\#(i,j)\mid i<j,\sigma(i)>\sigma(j)\}$ , and $\#$ stands for “the number of”. We will now prove that

\displaystyle c(r_{j})=\displaystyle\prod_{k=1}^{r_{j}-1}\dfrac{\left(k!\right)^{2}}{(2k+1)!}.

(73)

To this end, we write $r_{j}!\cdot c(r_{j})$ as an integral over $[0,1]^{r_{j}}$ :

\displaystyle r_{j}!\cdot c(r_{j})=\int_{[0,1]^{r_{j}}}\prod_{1\leq a<b\leq r_{j}}|y_{a}-y_{b}|\>{\rm{d}}y_{1}...{\rm{d}}y_{r_{j}}.

(74)

In 1955, de Bruijn [73, see toward the end of Sec. 9] used certain Pfaffians to evaluate

		$\displaystyle\int_{[0,1]^{r_{j}}}\prod_{1\leq a<b\leq r_{j}}\|y_{a}-y_{b}\|\>{\rm{d}}y_{1}...{\rm{d}}y_{r_{j}}$
	$\displaystyle=$	$\displaystyle\frac{r_{j}!\cdot\{1!\times 2!\times...\times(r_{j}-1)!\}^{2}}{1!\times 3!\times...\times(2r_{j}-1)!},\quad r_{j}=2,3,\ldots,$

which upon substitution in (74), indeed yields (73).

Combining (70) and (73), we arrive at

	$\displaystyle\!\!\!\mathrm{vol}\left(\mathcal{R}_{j}\left(\{\bm{x}_{0}\},t\right)\right)$	$\displaystyle=\dfrac{(2\mu_{j})^{r_{j}}t^{r_{j}(r_{j}+1)/2}}{\displaystyle\prod_{k=1}^{r_{j}-1}k!}c(r_{j})$
		$\displaystyle=(2\mu_{j})^{r_{j}}t^{r_{j}(r_{j}+1)/2}\displaystyle\prod_{k=1}^{r_{j}-1}\dfrac{k!}{(2k+1)!}.$		(75)

Finally, substituting (75) in (65), and recalling that $r_{1}+r_{2}+\ldots+r_{m}=d$ , the expression (41) follows. $\blacksquare$

-H Proof of Theorem 6

From (13), the subvector $\bm{\xi}_{j}(s)$ , where $j\in[m]$ , is component-wise nonnegative for all $s\in[0,t]$ .

Therefore, by triangle inequality, we have

\displaystyle\!\int_{0}^{t}\!\lvert\langle\bm{\eta},\bm{\xi}(s)\rangle\rvert\>{\rm{d}}s\leq\!\int_{0}^{t}\!\sum_{j=1}^{m}\langle\lvert\bm{\eta}_{j}\rvert,\mu_{j}\bm{\xi}_{j}(s)\rangle=\langle|\bm{\eta}|,\bm{\zeta}(t)\rangle,

(76)

where $|\bm{\eta}_{j}|$ denotes the $j$ th subvector with component-wise absolute values. Let us call $|\bm{\eta}|$ as the “absolute unit vector”.

The upper bound in (76) is convex in $\bm{\eta}$ , and is maximized by an absolute unit vector collinear with $\bm{\zeta}(t)$ given by

\displaystyle\bm{\eta}=\pm\dfrac{\bm{\zeta}(t)}{\parallel\bm{\zeta}(t)\parallel_{2}},

(77)

i.e., the unit vectors associated with $\bm{\zeta}(t)$ up to plus-minus sign permutations among its components.

Out of the $2^{d}$ unit vectors given by (77), the “all plus” and “all minus” unit vectors achieve equality in (76), and hence must be the maximizers of (43). The inequality (76) remains strict for the remaining $2^{d}-2$ unit vectors in (77), thus are suboptimal for (43). Therefore, the maximizers in (44) are

\bm{\eta}^{\max}=\bm{\zeta}(t)/\parallel\bm{\zeta}(t)\parallel_{2},\quad-\bm{\zeta}(t)/\parallel\bm{\zeta}(t)\parallel_{2},

which upon substitution in (43), results in (45). $\blacksquare$

References

[1] P. Varaiya, “Reach set computation using optimal control,” in Verification of Digital and Hybrid Systems. Springer, 2000, pp. 323–331.
[2] M. Althoff, “An introduction to CORA 2015,” in Proc. of the Workshop on 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 reach sets,” in 2020 American Control Conference (ACC). IEEE, 2020, pp. 4466–4471.
[5] R. J. Aumann, “Integrals of set-valued functions,” Journal of Mathematical Analysis and Applications, vol. 12, no. 1, pp. 1–12, 1965.
[6] S. Haddad and A. Halder, “Boundary and taxonomy of integrator reach sets,” in 2022 American Control Conference (ACC). IEEE, 2022, pp. 4133–4138.
[7] J. M. Borwein and R. E. Crandall, “Closed forms: what they are and why we care,” Notices of the AMS, vol. 60, no. 1, pp. 50–65, 2013.
[8] J.-B. Hiriart-Urruty and C. Lemaréchal, Convex analysis and minimization algorithms I: Fundamentals. Springer Science & Business media, 2013, vol. 305.
[9] A. A. Liapounoff, “Sur les fonctions-vecteurs completement additives,” Izvestiya Rossiiskoi Akademii Nauk. Seriya Matematicheskaya, vol. 4, no. 6, pp. 465–478, 1940.
[10] P. R. Halmos, “The range of a vector measure,” Bulletin of the American Mathematical Society, vol. 54, no. 4, pp. 416–421, 1948.
[11] J. Diestel and B. Faires, “On vector measures,” Transactions of the American Mathematical Society, vol. 198, pp. 253–271, 1974.
[12] N. Dinculeanu, Vector measures. Elsevier, 2014.
[13] R. Schneider, Convex bodies: the Brunn–Minkowski theory. Cambridge University Press, 2014, no. 151.
[14] P. McMullen, “On zonotopes,” Transactions of the American Mathematical Society, vol. 159, pp. 91–109, 1971.
[15] G. Shephard, “Combinatorial properties of associated zonotopes,” Canadian Journal of Mathematics, vol. 26, no. 2, pp. 302–321, 1974.
[16] H. S. M. Coxeter, Regular polytopes. Courier Corporation, 1973.
[17] A. Girard, “Reachability of uncertain linear systems using zonotopes,” in International Workshop on Hybrid Systems: Computation and Control. Springer, 2005, pp. 291–305.
[18] A. Girard and C. Le Guernic, “Zonotope/hyperplane intersection for hybrid systems reachability analysis,” in International Workshop on Hybrid Systems: Computation and Control. Springer, 2008, pp. 215–228.
[19] M. Althoff, O. Stursberg, and M. Buss, “Computing reachable sets of hybrid systems using a combination of zonotopes and polytopes,” Nonlinear Analysis: Hybrid Systems, vol. 4, no. 2, pp. 233–249, 2010.
[20] M. Althoff and B. H. Krogh, “Zonotope bundles for the efficient computation of reachable sets,” in 2011 50th IEEE Conference on Decision and Control and European Control Conference. IEEE, 2011, pp. 6814–6821.
[21] J. K. Scott, D. M. Raimondo, G. R. Marseglia, and R. D. Braatz, “Constrained zonotopes: A new tool for set-based estimation and fault detection,” Automatica, vol. 69, pp. 126–136, 2016.
[22] A. S. Adimoolam and T. Dang, “Using complex zonotopes for stability verification,” in 2016 American Control Conference (ACC). IEEE, 2016, pp. 4269–4274.
[23] M. Althoff, “Reachability analysis of nonlinear systems using conservative polynomialization and non-convex sets,” in Proceedings of the 16th International Conference on Hybrid Systems: Computation and Control, 2013, pp. 173–182.
[24] N. Kochdumper and M. Althoff, “Sparse polynomial zonotopes: A novel set representation for reachability analysis,” IEEE Transactions on Automatic Control, 2020.
[25] D. Cox, J. Little, and D. OShea, Ideals, varieties, and algorithms: an introduction to computational algebraic geometry and commutative algebra. Springer Science & Business Media, 2013.
[26] J. Yong and X. Y. Zhou, Stochastic controls: Hamiltonian systems and HJB equations. Springer Science & Business Media, 1999, vol. 43.
[27] S. Haddad and A. Halder, “Certifying the intersection of reach sets of integrator agents with set-valued input uncertainties,” IEEE Control Systems Letters, 2022.
[28] E. D. Bolker, “A class of convex bodies,” Transactions of the American Mathematical Society, vol. 145, pp. 323–345, 1969.
[29] R. Schneider and W. Weil, “Zonoids and related topics,” in Convexity and its Applications. Springer, 1983, pp. 296–317.
[30] P. Goodey and W. Wolfgang, “Zonoids and generalisations,” in Handbook of Convex Geometry. Elsevier, 1993, pp. 1297–1326.
[31] J. Bourgain, J. Lindenstrauss, and V. Milman, “Approximation of zonoids by zonotopes,” Acta Mathematica, vol. 162, no. 1, pp. 73–141, 1989.
[32] C. Vinzant, “The geometry of spectrahedra,” in Sum of Squares: Theory and Applications. American Mathematical Society, 2020, pp. 11–36.
[33] A. Tarski, “A decision method for elementary algebra and geometry,” in Quantifier elimination and cylindrical algebraic decomposition. Springer, 1998, pp. 24–84.
[34] A. Seidenberg, “A new decision method for elementary algebra,” Annals of Mathematics, pp. 365–374, 1954.
[35] J. Bochnak, M. Coste, and M.-F. Roy, Real algebraic geometry. Springer Science & Business Media, 2013, vol. 36.
[36] G. Blekherman, P. A. Parrilo, and R. R. Thomas, Semidefinite optimization and convex algebraic geometry. SIAM, 2012.
[37] A. B. Kurzhanski and P. Varaiya, Dynamics and Control of Trajectory Tubes: Theory and Computation. Springer, 2014, vol. 85.
[38] G. Zaimi, “A polynomial implicitization,” MathOverflow. [Online]. Available: https://mathoverflow.net/q/381335
[39] H. S. Wilf, generatingfunctionology. CRC Press, 2005.
[40] L. Kronecker, Zur Theorie der Elimination einer Variabeln aus zwei algebraischen Gleichungen. Buchdruckerei der Königl. Akademie der Wissenschaften (G. Vogt), 1881.
[41] R. Salem, Algebraic numbers and Fourier analysis. Wadsworth Publishing Company, 1983.
[42] R. Schneider, “Zonoids whose polars are zonoids,” Proceedings of the American Mathematical Society, vol. 50, no. 1, pp. 365–368, 1975.
[43] Y. Lonke, “On zonoids whose polars are zonoids,” Israel Journal of Mathematics, vol. 102, no. 1, pp. 1–12, 1997.
[44] J. W. Helton and V. Vinnikov, “Linear matrix inequality representation of sets,” Communications on Pure and Applied Mathematics, vol. 60, no. 5, pp. 654–674, 2007.
[45] J. W. Helton and J. Nie, “Semidefinite representation of convex sets,” Mathematical Programming, vol. 122, no. 1, pp. 21–64, 2010.
[46] S. P. Boyd and L. Vandenberghe, Convex optimization. Cambridge University Press, 2004.
[47] L.-L. Xie and P. R. Kumar, “A network information theory for wireless communication: Scaling laws and optimal operation,” IEEE transactions on information theory, vol. 50, no. 5, pp. 748–767, 2004.
[48] F. Xue, P. R. Kumar et al., “Scaling laws for ad hoc wireless networks: an information theoretic approach,” Foundations and Trends® in Networking, vol. 1, no. 2, pp. 145–270, 2006.
[49] C. Fan, J. Kapinski, X. Jin, and S. Mitra, “Locally optimal reach set over-approximation for nonlinear systems,” in 2016 International Conference on Embedded Software (EMSOFT). IEEE, 2016, pp. 1–10.
[50] Y. Meng, D. Sun, Z. Qiu, M. T. B. Waez, and C. Fan, “Learning density distribution of reachable states for autonomous systems,” in Conference on Robot Learning. PMLR, 2022, pp. 124–136.
[51] OEIS Foundation Inc., “The On-Line Encyclopedia of Integer Sequences,” http://oeis.org/A107254, 2019.
[52] S. R. Finch, Mathematical constants. Cambridge University Press, 2003.
[53] N. G. De Bruijn, Asymptotic methods in analysis. Courier Corporation, 1981, vol. 4.
[54] M. Abramowitz and I. A. Stegun, Handbook of mathematical functions with formulas, graphs, and mathematical tables. US Government Printing Office, 1970, vol. 55.
[55] F. C. Schweppe, Uncertain dynamic systems. Prentice Hall, 1973.
[56] F. Chernousko, “Optimal guaranteed estimates of indeterminacies with the aid of ellipsoids, I,” Engineering Cybernetics, vol. 18, no. 3, pp. 1–9, 1980.
[57] ——, “Guaranteed ellipsoidal estimates of uncertainties in control problems,” IFAC Proceedings Volumes, vol. 14, no. 2, pp. 869–874, 1981.
[58] D. Maksarov and J. Norton, “State bounding with ellipsoidal set description of the uncertainty,” International Journal of Control, vol. 65, no. 5, pp. 847–866, 1996.
[59] C. Durieu, E. Walter, and B. Polyak, “Multi-input multi-output ellipsoidal state bounding,” Journal of Optimization Theory and Applications, vol. 111, no. 2, pp. 273–303, 2001.
[60] T. Alamo, J. M. Bravo, and E. F. Camacho, “Guaranteed state estimation by zonotopes,” Automatica, vol. 41, no. 6, pp. 1035–1043, 2005.
[61] A. Halder, “On the parameterized computation of minimum volume outer ellipsoid of Minkowski sum of ellipsoids,” in 2018 IEEE Conference on Decision and Control (CDC). IEEE, 2018, pp. 4040–4045.
[62] ——, “Smallest ellipsoid containing p-sum of ellipsoids with application to reachability analysis,” IEEE Transactions on Automatic Control, 2020.
[63] “CORA: A tool for continuous reachability analysis,” https://tumcps.github.io/CORA/, accessed: 2021-02-14.
[64] A. Kurzhanskiĭ and I. Vályi, Ellipsoidal calculus for estimation and control. Nelson Thornes, 1997.
[65] F. John, “Extremum problems with inequalities as subsidiary conditions,” Studies and Essays: Courant Anniversary Volume, presented to R. Courant on his 60th Birthday, pp. 187–204, 1948.
[66] M. Henk, “Löwner-John ellipsoids,” Documenta Math, vol. 95, p. 106, 2012.
[67] S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan, Linear matrix inequalities in system and control theory. SIAM, 1994.
[68] M. Grant and S. Boyd, “CVX: Matlab Software for Disciplined Convex Programming, version 2.1,” http://cvxr.com/cvx, Mar 2014.
[69] G. M. Ziegler, Lectures on polytopes. Springer Science & Business Media, 2012, vol. 152.
[70] R. A. Horn and C. R. Johnson, Matrix analysis. Cambridge University Press, 2012.
[71] T. M. Apostol, “An elementary view of Euler’s summation formula,” The American Mathematical Monthly, vol. 106, no. 5, pp. 409–418, 1999.
[72] E. Hairer and G. Wanner, Analysis by its History. Springer Science & Business Media, 2006, vol. 2.
[73] N. De Bruijn, “On some multiple integrals involving determinants,” J. Indian Math. Soc, vol. 19, pp. 133–151, 1955.