Sparse Polynomial Zonotopes: A Novel Set Representation for Reachability Analysis

TL;DR

This paper introduces sparse polynomial zonotopes, a new set representation that efficiently handles non-convex sets and nonlinear system reachability, reducing computational complexity and wrapping effects.

Contribution

The paper presents sparse polynomial zonotopes as a novel set representation that generalizes existing models and improves efficiency in reachability analysis of nonlinear systems.

Findings

Reduces wrapping effect in nonlinear reachability analysis

Decreases computation time compared to zonotopes

Handles non-convex sets effectively

Abstract

We introduce sparse polynomial zonotopes, a new set representation for formal verification of hybrid systems. Sparse polynomial zonotopes can represent non-convex sets and are generalizations of zonotopes, polytopes, and Taylor models. Operations like Minkowski sum, quadratic mapping, and reduction of the representation size can be computed with polynomial complexity w.r.t. the dimension of the system. In particular, for reachability analysis of nonlinear systems, the wrapping effect is substantially reduced using sparse polynomial zonotopes, as demonstrated by numerical examples. In addition, we can significantly reduce the computation time compared to zonotopes when dealing with nonlinear dynamics.

Tables

Table 1. Table 1: Computational complexity with respect to the dimension n 𝑛 n .

Set Operation	SPZ	Quad. Zono. [5]
Multiplication with matrix	$𝒪\left(n^{2}m\right)$	$𝒪\left(n^{2}m\right)$
Mink. add. with zonotope	$𝒪\left(1\right)$	$𝒪\left(n\right)$
Enclosure by zonotope	$𝒪\left(n^2\right)$	$𝒪\left(n^2\right)$
Quadratic map	$𝒪\left(n\log \left(n\right)+n^{3}m\right)$	$𝒪\left(n^{4}m\right)$

Table 2. Table 2: Computation times for the spacecraft rendezvous benchmark. The results for the different tools are taken from [ 36 , Tab. 4] . The computation times are measured on the machines of the participants (see [ 36 , Appendix A] ).

Tool	Comp. Time [s]	Set Rep.	Language
Ariadne [15]	$172$	Taylor models	C++
CORA [6]	$11.8$	Zonotopes	MATLAB
DynIbex [58]	$294$	Zonotopes	C++
Flow* [20]	$18.7$	Taylor models	C++
Isabelle/HOL [34]	$295$	Zonotopes	SML
Our approach	$10.1$	SPZs	MATLAB

Table 3. Table 3: Computation times in seconds for the transcriptional regulator network for different system dimensions.

System Dimension	$𝐧=\mathrm{𝟏𝟐}$	$𝐧=\mathrm{𝟐𝟒}$	$𝐧=\mathrm{𝟑𝟔}$	$𝐧=\mathrm{𝟒𝟖}$
Computation Time [s]	6	20	54	122

Equations

$$\begin{split}\mathcal{PZ}=\bigg{\{}&\sum_{i=1}^{h}\bigg{(}\prod_{k=1}^{p}\alpha_{k}^{E_{(k,i)}}\bigg{)}G_{(\cdot,i)}+\sum_{j=1}^{q}\beta_{j}G_{I(\cdot,j)}~{}\bigg{|}\\ &\alpha_{k},\beta_{j}\in[-1,1]\bigg{\}}.\end{split}$$

$$h=c_{h}n,p=c_{p}n,q=c_{q}n,$$

$$\mathcal{P}\mathcal{Z}={⟨\left[\begin{array}{cccc}4& 2& 1& 2\\ 4& 0& 2& 2\end{array}\right],\left[\begin{array}{c}1\\ 0\end{array}\right],\left[\begin{array}{cccc}0& 1& 0& 3\\ 0& 0& 1& 1\end{array}\right],[12]⟩}_{PZ}$$

$$\begin{split}\mathcal{PZ}=\bigg{\{}&\begin{bmatrix}4\\ 4\end{bmatrix}+\begin{bmatrix}2\\ 0\end{bmatrix}\alpha_{1}+\begin{bmatrix}1\\ 2\end{bmatrix}\alpha_{2}+\begin{bmatrix}2\\ 2\end{bmatrix}\alpha_{1}^{3}\alpha_{2}+\begin{bmatrix}1\\ 0\end{bmatrix}\beta_{1}~{}\bigg{|}\\ &~{}\alpha_{1},\alpha_{2},\beta_{1}\in[-1,1]\bigg{\}}.\end{split}$$

$$\begin{split}\mathcal{PZ}=\bigg{\{}&c+\sum_{j=1}^{p}\alpha_{j}~{}g^{([1],j)}+\sum_{j=1}^{p}\sum_{k=j}^{p}\alpha_{j}\alpha_{k}~{}g^{([2],j,k)}+\dots+\\ &\sum_{j=1}^{p}\sum_{k=1}^{p}\dots\sum_{m=l}^{p}\alpha_{j}\alpha_{k}\dots\alpha_{m}~{}g^{([\mu],j,k,\dots,m)}\\ &+\sum_{i=1}^{q}\beta_{j}~{}G_{I(\cdot,j)}~{}\bigg{|}~{}\alpha_{i},\beta_{j}\in[-1,1]\bigg{\}},\end{split}$$

$$\begin{split}&{\texttt{mergeID}}(\mathcal{PZ}_{1},\mathcal{PZ}_{2})=\big{\{}\langle G_{1},G_{I,1},\overline{E}_{1},\overline{id}\rangle_{PZ},\\ &~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\langle G_{2},G_{I,2},\overline{E}_{2},\overline{id}\rangle_{PZ}\big{\}}\\ &~{}~{}\\ &\text{with}~{}~{}\overline{id}=\begin{bmatrix}id_{1}&id_{2(\cdot,\mathcal{H})}\end{bmatrix},~{}~{}\mathcal{H}=\left\{i~{}|~{}id_{2(i)}\not\in id_{1}\right\},\\ &~{}~{}~{}~{}~{}~{}~{}\overline{E}_{1}=\begin{bmatrix}E_{1}\\ \mathbf{0}^{(|\mathcal{H}|,h_{1})}\end{bmatrix}\in\mathbb{R}^{a\times h_{1}},\\ &~{}~{}~{}~{}~{}~{}~{}\overline{E}_{2(i,\cdot)}=\begin{cases}E_{2(j,\cdot)},~{}\mathrm{if}~{}\exists j~{}\overline{id}_{(i)}=id_{2(j)}\\ \mathbf{0}^{(1,h_{2})},~{}\mathrm{otherwise}\end{cases}i=1\dots a,\end{split}$$

$$\begin{split}&{\texttt{compact}}(\mathcal{PZ})=\langle\overline{G},G_{I},\overline{E},id\rangle_{PZ}\\ &~{}~{}\\ &\text{with}~{}~{}\overline{E}={\texttt{uniqueColumns}}(E)\in\mathbb{N}_{0}^{p\times k},\\ &~{}~{}~{}~{}~{}~{}~{}~{}\mathcal{H}_{j}=\big{\{}i~{}|~{}\overline{E}_{(l,j)}=E_{(l,i)}~{}\forall l\in\{1,\dots,p\}\big{\}},\\ &~{}~{}~{}~{}~{}~{}~{}~{}\overline{G}=\bigg{[}\sum_{i\in\mathcal{H}_{1}}G_{(\cdot,i)}~{}\dots~{}\sum_{i\in\mathcal{H}_{k}}G_{(\cdot,i)}\bigg{]},\end{split}$$

$$\begin{split}&\bigg{\{}\bigg{(}\prod_{k=1}^{p}\alpha_{k}^{e_{(k)}}\bigg{)}G_{(\cdot,1)}+\bigg{(}\prod_{k=1}^{p}\alpha_{k}^{e_{(k)}}\bigg{)}G_{(\cdot,2)}~{}\bigg{|}~{}\alpha_{k}\in[-1,1]\bigg{\}}\\ &=\bigg{\{}\bigg{(}\prod_{k=1}^{p}\alpha_{k}^{e_{(k)}}\bigg{)}\bigg{(}G_{(\cdot,1)}+G_{(\cdot,2)}\bigg{)}~{}\bigg{|}~{}\alpha_{k}\in[-1,1]\bigg{\}}.\end{split}$$

$$\mathcal{Z}=\bigg{\{}c+\sum_{i=1}^{l}\alpha_{i}~{}G_{(\cdot,i)}~{}\bigg{|}~{}\alpha_{i}\in[-1,1]\bigg{\}}.$$

$$\mathcal{Z}={⟨[cG],[],[0^{(n,1)}{I}_l],\text{uniqueID}(l)⟩}_{PZ}.$$

$$\mathcal{P}=\bigg{\{}\sum_{i=1}^{r}\lambda_{i}v_{i}~{}\bigg{|}~{}\lambda_{i}\geq 0,~{}\sum_{i=1}^{r}\lambda_{i}=1\bigg{\}}.$$

$$w_{(i)}\left(x_1,\dots ,x_s\right)=\sum _{j=1}^{m_i}{b}_{i,j}\prod _{k=1}^s{x}_k^{E_{i(k,j)}},$$

$$\mathcal{T}(x)=\left\{\left[\begin{array}{c}{w}_{(1)}(x)\\ ⋮\\ w_{(n)}(x)\end{array}\right]+\left[\begin{array}{c}{y}_{(1)}\\ ⋮\\ y_{(n)}\end{array}\right]|y\in \mathcal{I}\right\},$$

$$\mathcal{T}(x)={⟨\left[\begin{array}{cc}\frac{l_R+u_R}{2}& \hat{G}\end{array}\right],G_I,\left[\begin{array}{cc}{0}^{(s,1)}& \hat{E}\end{array}\right],\text{uniqueID}(s)⟩}_{PZ}$$

$$\begin{split}&\text{with}~{}~{}\widehat{G}=\begin{bmatrix}\big{[}\overline{b}_{1,1}~{}\dots~{}\overline{b}_{1,m_{1}}\big{]}&&\mathbf{0}\\ &\ddots&\\ \mathbf{0}&&\big{[}\overline{b}_{n,1}~{}\dots~{}\overline{b}_{n,m_{n}}\big{]}\end{bmatrix},\\ &~{}~{}~{}~{}~{}~{}~{}~{}\widehat{E}=\begin{bmatrix}\overline{E}_{1}&\dots&\overline{E}_{n}\end{bmatrix},\\ &~{}~{}~{}~{}~{}~{}~{}~{}G_{I}=\begin{bmatrix}\frac{u_{R(1)}-l_{R(1)}}{2}&&\mathbf{0}\\ &\ddots&\\ \mathbf{0}&&\frac{u_{R(n)}-l_{R(n)}}{2}\end{bmatrix}.\end{split}$$

$$w_{(i)}\left({\delta }_1({\alpha }_1),\dots ,{\delta }_s({\alpha }_s)\right):=\sum _{j=1}^{{\overline{m}}_i}{\overline{b}}_{i,j}\prod _{k=1}^s{\alpha }_k^{{\overline{E}}_{i(k,j)}},$$

$$\begin{array}{rl}{\delta }_k({\alpha }_k)=& \frac{l_{D(k)}+u_{D(k)}}{2}+\frac{u_{D(k)}-l_{D(k)}}{2}{\alpha }_k,\\ & {\alpha }_k\in [-1,1],k=1\dots s.\end{array}$$

$$\mathcal{D}=\big{\{}\begin{bmatrix}\delta_{1}(\alpha_{1})&\dots&\delta_{s}(\alpha_{s})\end{bmatrix}^{T}~{}|~{}\alpha_{1},\dots,\alpha_{s}\in[-1,1]\big{\}}.$$

$$\begin{split}&\mathcal{T}(\delta(\alpha))=\left\{\left.\begin{bmatrix}w_{(1)}(\delta(\alpha))\\ \vdots\\ w_{(n)}(\delta(\alpha))\end{bmatrix}+\begin{bmatrix}y_{(1)}\\ \vdots\\ y_{(n)}\end{bmatrix}~{}\right|~{}y\in\mathcal{I}\right\}\overset{\begin{subarray}{c}\scriptscriptstyle\eqref{eq:TaylorModelTransform}\\ \scriptscriptstyle\text{and}\\ \scriptscriptstyle\mathcal{I}=[l_{\scaleto{R\mathstrut}{3pt}},u_{\scaleto{R\mathstrut}{3pt}}]\end{subarray}}{=}\\ &~{}\\ &\Bigg{\{}\sum_{j=1}^{\overline{m}_{1}}\begin{bmatrix}\overline{b}_{1,j}\\ o\end{bmatrix}\prod_{k=1}^{s}\alpha_{k}^{\overline{E}_{1(k,j)}}+\dots+\sum_{j=1}^{\overline{m}_{n}}\begin{bmatrix}o\\ \overline{b}_{n,j}\end{bmatrix}\prod_{k=1}^{s}\alpha_{k}^{\overline{E}_{n(k,j)}}\\ &~{}~{}+\frac{1}{2}\begin{bmatrix}u_{R(1)}-l_{R(1)}\\ o\end{bmatrix}\beta_{1}+\dots+\frac{1}{2}\begin{bmatrix}o\\ u_{R(n)}-l_{R(n)}\end{bmatrix}\beta_{n}\\ &~{}~{}+\frac{l_{R}+u_{R}}{2}~{}\bigg{|}~{}\alpha_{k},\beta_{1},\dots,\beta_{n}\in[-1,1]\Bigg{\}}\\ &~{}\\ &=\left\langle\begin{bmatrix}\frac{l_{R}+u_{R}}{2}&\widehat{G}\end{bmatrix},G_{I},\begin{bmatrix}\mathbf{0}^{(s,1)}&\widehat{E}\end{bmatrix},{\texttt{uniqueID}}(s)\right\rangle_{PZ},\end{split}$$

$$\begin{split}\mathcal{PZ}=\bigg{\langle}&\begin{bmatrix}-0.5&1&0&-1&1\\ -0.5&1&1&1&1\end{bmatrix},[~{}],\\ &\begin{bmatrix}0&1&0&1&2\\ 0&0&1&1&0\end{bmatrix},[1~{}2]\bigg{\rangle}_{PZ}\end{split}$$

$$\begin{split}&{\texttt{zono}}\left(\mathcal{PZ}\right)=\bigg{\langle}\sum_{i\in\mathcal{N}}G_{\left(\cdot,i\right)}+0.5\sum_{i\in\mathcal{H}}G_{\left(\cdot,i\right)},\\ &~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\begin{bmatrix}0.5~{}G_{(\cdot,\mathcal{H})}&G_{(\cdot,\mathcal{K})}&G_{I}\end{bmatrix}\bigg{\rangle}_{Z}\\ &~{}~{}\\ &\text{with}~{}~{}\mathcal{N}=\left\{i~{}\left|~{}E_{(j,i)}=0~{}\forall j\in\{1,\dots,p\}\right.\right\},\\ &~{}~{}~{}~{}~{}~{}~{}~{}\mathcal{H}=\bigg{\{}i~{}\bigg{|}~{}\prod_{j=1}^{p}\left(1-(E_{(j,i)}~{}\mathrm{mod}~{}2)\right)=1\bigg{\}}\setminus\mathcal{N},\\ &~{}~{}~{}~{}~{}~{}~{}~{}\mathcal{K}=\{1,\dots,h\}\setminus(\mathcal{H}\cup\mathcal{N}),\\ &~{}~{}\end{split}$$

$$\begin{array}{rl}\forall i\in \mathcal{H}:& \left(\prod _{k=1}^p[-1,1{]}^{E_{(k,i)}}\right)G_{(\cdot ,i)}=[0,1]G_{(\cdot ,i)}=\\ & =0.5G_{(\cdot ,i)}+[-1,1]0.5G_{(\cdot ,i)}.\end{array}$$

$$\text{poly}(\mathcal{P}\mathcal{Z})=⟨[v_1\dots v_r]{⟩}_P,$$

$$\begin{split}&\overline{\mathcal{PZ}}=\big{\langle}[c_{z}~{}\overline{G}_{(\cdot,\mathcal{K})}],[G_{I}~{}G_{z}],[\mathbf{0}^{(n,1)}~{}\overline{E}_{(\cdot,\mathcal{K})}],id\big{\rangle}_{PZ}\\ &~{}\\ &\text{with}~{}\mathcal{H}=\big{\{}i~{}\big{|}~{}\exists j\in\{1,\dots,p\}:~{}E_{(j,i)}>1\big{\}},\\ &~{}~{}~{}~{}~{}~{}\mathcal{K}=\{1,\dots,h\}\setminus\mathcal{H}\\ &~{}~{}~{}~{}~{}~{}\langle c_{z},G_{z}\rangle_{Z}={\texttt{zono}}\big{(}\langle G_{(\cdot,\mathcal{H})},[~{}],E_{(\cdot,\mathcal{H})},id\rangle_{PZ}\big{)}\end{split}$$

$$\mathcal{P}\mathcal{Z}\subseteq \overline{\mathcal{P}\mathcal{Z}}\stackrel{\text{cite[cite][@@bibrefKochdumper2019a, Alg. 2]}}{\subseteq }⟨[v_1\dots v_r]{⟩}_P.$$

$$s_{\mathcal{S}}(d)=\underset{x\in \mathcal{S}}{\max }{d}^{T}x.$$

$$\text{B}(f(x),\mathcal{I})\supseteq \left[\underset{x\in \mathcal{I}}{\min }f(x),\underset{x\in \mathcal{I}}{\max }f(x)\right]$$

$$\begin{split}&\widehat{s}_{\mathcal{PZ}}(d)=u+\sum_{j=1}^{q}\left|\overline{g}_{I(j)}\right|\\ &~{}\\ &\text{with}~{}~{}\langle\overline{g},\overline{g}_{I},E,id\rangle_{PZ}:=d^{T}\otimes\mathcal{PZ},\\ &~{}~{}~{}~{}~{}~{}~{}~{}[l,u]={\texttt{B}}\big{(}w(\alpha_{1},\dots,\alpha_{p}),[-\mathbf{1}^{(n,1)},\mathbf{1}^{(n,1)}]\big{)},\\ &~{}~{}~{}~{}~{}~{}~{}~{}w(\alpha_{1},\dots,\alpha_{p})=\sum_{i=1}^{h}\left(\prod_{k=1}^{p}\alpha_{k}^{E_{(k,i)}}\right)\overline{g}_{(i)},\end{split}$$

$$\begin{split}d^{T}\otimes\mathcal{PZ}=&\underbrace{\bigg{\{}\sum_{i=1}^{h}\bigg{(}\prod_{k=1}^{p}\alpha_{k}^{E_{(k,i)}}\bigg{)}\overline{g}_{(i)}~{}\bigg{|}~{}\alpha_{k}\in[-1,1]\bigg{\}}}_{\subseteq[l,u]~{}\text{(dependent}~{}\text{part)}}\\ &\oplus\underbrace{\bigg{\{}\sum_{j=1}^{q}\beta_{j}\overline{g}_{I(j)}~{}\bigg{|}~{}\beta_{j}\in[-1,1]\bigg{\}}}_{\overset{\equiv\big{[}-\sum_{j=1}^{q}\left|\overline{g}_{I(j)}\right|,~{}\sum_{j=1}^{q}\left|\overline{g}_{I(j)}\right|\big{]}}{\overset{~{}}{\text{(independent}~{}\text{part)}}}}.\end{split}$$

$$M\otimes \mathcal{P}\mathcal{Z}={⟨MG,M{G}_I,E,id⟩}_{PZ},$$

Full text

{IEEEkeywords}

Reachability analysis, nonlinear dynamics, hybrid systems, sparse polynomial zonotopes.

1 Introduction

\IEEEPARstart

Efficient set representations are highly relevant for controller synthesis and formal verification of hybrid systems, since many underlying algorithms compute with sets; see e.g., [18, 63, 38, 60]. Improvements originating from a new set representation often significantly reduce computation time and improve the accuracy of set-based computations.

1.1 State of the Art

Fig. 1 shows relevant set representations and their relations to each other. Almost all typical set representations are convex, except Taylor models, star sets, level sets, and polynomial zonotopes. Since all convex sets can be represented by support functions, which are closed under Minkowski addition, linear maps, and convex hull operations, they are a good choice for reachability analysis [32, 28, 29, 59, 55]. Ellipsoids and polytopes are special cases of support functions, which are often used for reachability analysis [22, 44, 59, 17, 24] and computations of invariant sets [16, 3, 57, 45]. However, the disadvantage of ellipsoids is that they are not closed under intersection and Minkowski addition; the disadvantage of polytopes is that the Minkowski sum is computationally expensive [62].

One important subclass of polytopes is zonotopes, which can be represented compactly by so-called generators: a zonotope with $l$ generators in $n$ dimensions might have up to ${l\choose n-1}$ halfspaces. More importantly, Minkowski sum and linear maps can be computed cheaply, making them a good choice for reachability analysis [31, 8, 35, 34, 59]. Zonotopes are closely related to affine arithmetic [25] with the zonotope factors being identical to the noise symbols in affine arithmetic. Two relevant extensions to zonotopes are zonotope bundles [10], where the set is represented implicitly by the intersection of several zonotopes, and constrained zonotopes [61], where additional equality constraints on the zonotope factors are considered. Zonotope bundles, as well as constrained zonotopes, are both able to represent any bounded polytope. Both representations make use of lazy computations and thus suffer much less from the curse of dimensionality, as is the case for polytopes [62]. Two other set representations related to zonotopes are complex zonotopes [1] and zonotopes with sub-polyhedric domains [30]. Complex zonotopes are defined by complex valued vectors and are well-suited to verify global exponential stability for systems with complex valued eigenvectors [1]. Zonotopes with sub-polyhedric domains use zones, octagons, and polyhedra instead of intervals to represent the domain for the zonotope factors, which enables the efficient computation of intersections and unions of sets by exploiting lazy computations [30]. A special case of zonotopes are multi-dimensional intervals, which are particularly useful for range bounding of nonlinear functions via interval arithmetic [37], but they are also used for reachability analysis [27, 54]. Since intervals are not closed under linear maps, one often has to split them to reduce the wrapping effect [46].

In general, reachable sets of nonlinear systems are non-convex, so that tight enclosures can only be achieved using non-convex set representations when avoiding the splitting of reachable sets. Taylor models [49], which consist of a polynomial and an interval remainder part, are an example of non-convex set representation. They are typically used for range bounding [50] and reachability analysis [21, 20, 48, 53]. Polynomial zonotopes, another type of non-convex set representation, are introduced in [5] and can equally represent the set defined by a Taylor model, as later shown in this work. Quadratic zonotopes [2] are a special case of polynomial zonotopes. Two other ways to represent non-convex sets are star sets, which are especially useful for simulation-based reachability analysis [26, 14], and level sets of nonlinear functions [51], which are applied to compute reachable sets [52] and controlled invariant regions [43]. While star sets and level sets are very expressive, it is yet unclear how some operations, such as nonlinear mapping, are computed.

1.2 Overview

In this work, we introduce a new non-convex set representation called sparse polynomial zonotopes, which is a non-trivial extension of polynomial zonotopes from [5] and exhibits the following major advantages: a) sparse polynomial zonotopes enable a very compact representation of sets, b) they are closed under all relevant operations, c) many other set representations can be converted to a sparse polynomial zonotope, and most importantly, d) all operations have at most polynomial complexity.

The remainder of this paper is structured as follows: In Sec. 2, the formal definition of sparse polynomial zonotopes is provided and important operations on them are derived. We show how sparse polynomial zonotopes provide substantially better results for reachability analysis in Sec. 3, which is demonstrated in Sec. 4 on four numerical examples.

1.3 Notation

In the remainder of this paper, we will use the following notations: Sets are always denoted by calligraphic letters, matrices by uppercase letters, and vectors by lowercase letters. Given a discrete set $\mathcal{H}\in\{\cdot\}^{n}$, $|\mathcal{H}|=n$ denotes the cardinality of the set. Given a vector $b\in\mathbb{R}^{n}$, $b_{(i)}$ refers to the $i$-th entry. Given a matrix $A\in\mathbb{R}^{n\times m}$, $A_{(i,\cdot)}$ represents the $i$-th matrix row, $A_{(\cdot,j)}$ the $j$-th column, and $A_{(i,j)}$ the $j$-th entry of matrix row $i$. Given a discrete set of positive integer indices $\mathcal{H}=\{h_{1},\dots,h_{|\mathcal{H}|}\}$ with $1\leq h_{i}\leq m~{}\forall i\in\{1,\dots,|\mathcal{H}|\}$, $A_{(\cdot,\mathcal{H})}$ is used for $[A_{(\cdot,h_{1})}~{}\dots~{}A_{(\cdot,h_{|\mathcal{H}|})}]$, where $[C~{}D]$ denotes the concatenation of two matrices $C$ and $D$. The symbols $\mathbf{0}^{(n,m)}\in\mathbb{R}^{n\times m}$ and $\mathbf{1}^{(n,m)}\in\mathbb{R}^{n\times m}$ represent matrices of zeros and ones, and $I_{n}\in\mathbb{R}^{n\times n}$ is the identity matrix. The empty matrix is denoted by $[~{}]$. The left multiplication of a matrix $M\in\mathbb{R}^{m\times n}$ with a set $\mathcal{S}\subset\mathbb{R}^{n}$ is defined as $M\otimes\mathcal{S}=\{Ms~{}|~{}s\in\mathcal{S}\}$, the Minkowski addition of two sets $\mathcal{S}_{1}\subset\mathbb{R}^{n}$ and $\mathcal{S}_{2}\subset\mathbb{R}^{n}$ is defined as $\mathcal{S}_{1}\oplus\mathcal{S}_{2}=\{s_{1}+s_{2}~{}|~{}s_{1}\in\mathcal{S}_{1},s_{2}\in\mathcal{S}_{2}\}$, and the Cartesian product of two sets $\mathcal{S}_{1}\subset\mathbb{R}^{n}$ and $\mathcal{S}_{2}\subset\mathbb{R}^{m}$ is defined as $\mathcal{S}_{1}\times\mathcal{S}_{2}=\{[s_{1}~{}s_{2}]^{T}~{}|~{}s_{1}\in\mathcal{S}_{1},s_{2}\in\mathcal{S}_{2}\}$. We further introduce an $n$-dimensional interval as $\mathcal{I}=[l,u],~{}\forall i~{}l_{(i)}\leq u_{(i)},~{}l,u\in\mathbb{R}^{n}$. The Nabla operator is defined as $\nabla=\sum_{i=1}^{n}e_{i}\frac{\partial}{\partial x_{(i)}}$, with $x\in\mathbb{R}^{n}$ and $e_{i}\in\mathbb{R}^{n}$ being orthogonal unit vectors. For the derivation of computational complexity, we consider all binary operations except concatenations, and initializations are explicitly not considered.

2 Sparse Polynomial Zonotopes

Let us first define sparse polynomial zonotopes (SPZs), followed by derivations of relevant operations on them.

Definition 1:

(Sparse Polynomial Zonotope) Given a generator matrix of dependent generators $G\in\mathbb{R}^{n\times h}$, a generator matrix of independent generators $G_{I}\in\mathbb{R}^{n\times q}$, and an exponent matrix $E\in\mathbb{N}_{0}^{p\times h}$, an SPZ is defined as

[TABLE]

The scalars $\alpha_{k}$ are called * dependent factors* since a change in their value affects multiplication with multiple generators. Consequently, the scalars $\beta_{j}$ are called independent factors because they only affect multiplication with one generator. The number of dependent factors is $p$, the number of independent factors is $q$, and the number of dependent generators is $h$. The order of an SPZ $\rho$ is defined as $\rho=\frac{h+q}{n}$. The independent generators are required for computational reasons: while computations on the dependent generators are exact but computational expensive, computations on the independent generators are often over-approximative but fast. For the derivation of the computational complexity of set operations, we introduce

[TABLE]

with $c_{h},c_{p},c_{q}\in\mathbb{R}_{\geq 0}$. The assumption in (2) is justified by the fact that we limit the order $\rho$ of an SPZ to stay below a desired value $\rho_{d}$. In the remainder of this paper, we call the term $\alpha_{1}^{E_{(1,i)}}\cdot\dotsc\cdot\alpha_{p}^{E_{(p,i)}}\cdot G_{(\cdot,i)}$ a monomial, and $\alpha_{1}^{E_{(1,i)}}\cdot\dotsc\cdot\alpha_{p}^{E_{(p,i)}}$ the variable part of the monomial. In order to keep track of the dependencies between the dependent factors from different SPZs, we introduce an unambiguous integer identifier for each dependent factor $\alpha_{k}$ and store the identifiers for all dependent factors in a row vector $id\in\mathbb{N}_{>0}^{1\times p}$. Using this identifier vector, we introduce the shorthand $\mathcal{PZ}=\langle G,G_{I},E,id\rangle_{PZ}\subset\mathbb{R}^{n}$ for the representation of SPZs. All components of a set $\square_{i}$ have the index $i$, e.g., $p_{i}$, $h_{i}$, and $q_{i}$ belong to $\mathcal{PZ}_{i}$. Since many set operations require the generation of new unique identifiers we introduce the operation ${\texttt{uniqueID}}(m)$, which returns an identifier vector of length $m$ that contains newly generated unique identifiers. The concept of using unique identifiers to keep track of dependencies is also used in [23], so that the operation uniqueID can be implemented as in [23, Tab. 3]. To make SPZs more intuitive, we introduce the following example:

Example 1:

The SPZ

[TABLE]

defines the set

[TABLE]

The construction of this SPZ is visualized in Fig. 2: (a) shows the set spanned by the constant offset vector and the second and third dependent generator, (b) shows the addition of the dependent generator with the mixed term $\alpha_{1}^{3}\alpha_{2}$, (c) shows the addition of the independent generator, and (d) visualizes the final set.

SPZs are a more compact representation of polynomial zonotopes [5], resulting in completely different algorithms for operations on them. In [5, Def. 1], the generators $g^{([i],j,k,\dots,m)}$ for all possible combinations of dependent factors up to a certain polynomial degree $\mu$ are stored:

[TABLE]

with $g^{([\mu],j,k,\dots,m)}\in\mathbb{R}^{n}$, $c\in\mathbb{R}^{n}$, $G_{I}\in\mathbb{R}^{n\times q}$. This results in $h={\mu+p\choose p}$ generators [56, Chapter 3 (3.8)]. For the one-dimensional polynomial zonotope $\mathcal{PZ}=\{\alpha_{1}\cdot\ldots\cdot\alpha_{19}\cdot\alpha_{20}^{10}|\alpha_{i}\in[-1,1]\}$ with $p=20$ dependent factors and with a polynomial degree of $\mu=10$, the number of dependent generators is $h=30045015$. This demonstrates that the number of stored generators can become very large if the polynomial degree and the number of dependent factors are high, which makes computations on the previous set representation very inefficient. Even in comparison with quadratic zonotopes, which correspond to a polynomial order of $\mu=2$, SPZs have lower or equal complexity for all set operations considered in [5] (see Tab. 1). We in turn use a sparse representation, where only the generators for desired factor combinations are stored, which enables the representation of the above polynomial zonotope with only one single generator. Furthermore, our representation does not require limiting the polynomial degree of the polynomial zonotope in advance, which is advantageous for reachability analysis, as shown in Sec. 3.2.

2.1 Preliminary Operations

First, we introduce preliminary set operations that are required for many other operations.

2.1.1 Merging the Set of Identifiers

For all set operations that involve two or more SPZs, the operator mergeID is necessary in order to build a common representation of exponent matrices to fully exploit the dependencies between identical dependent factors:

Proposition 1:

(Merge ID) Given two SPZs, $\mathcal{PZ}_{1}=\langle G_{1},\linebreak[1]G_{I,1},E_{1},id_{1}\rangle_{PZ}$ and $\mathcal{PZ}_{2}=\langle G_{2},\linebreak[1]G_{I,2},$ $E_{2},id_{2}\rangle_{PZ}$, mergeID returns two adjusted SPZs with identical identifier vectors that are equivalent to $\mathcal{PZ}_{1}$ and $\mathcal{PZ}_{2}$, and has a complexity of $\mathcal{O}(n^{2})$:

[TABLE]

where $a=|\mathcal{H}|+p_{1}$.

Proof.

The extension of the exponent matrices with all-zero rows only changes the representation of the set, but not the set itself.

Complexity: The only operation with super-linear complexity with respect to the system dimension $n$ is the construction of the set $\mathcal{H}$ with $\mathcal{O}(p_{1}p_{2})=\mathcal{O}(n^{2})$ using (2). ∎

2.1.2 Transformation to a Compressed Representation

Some set operations result in an SPZ that contains multiple monomials with an identical variable part, which we combine to one single monomial:

Proposition 2:

(Compact) Given an SPZ $\mathcal{PZ}=\langle G,G_{I},\linebreak[0]E,id\rangle_{PZ}$, the operation compact returns a compressed representation of the set $\mathcal{PZ}$ and has a complexity of $\mathcal{O}(n^{2}\log(n))$:

[TABLE]

where the operation uniqueColumns removes identical matrix columns until all columns are unique.

Proof.

For an SPZ where the exponent matrix $E=[e~{}e]$ consists of $2$ identical columns $e\in\mathbb{N}_{0}^{p}$, it holds that

[TABLE]

Summation of the generators for monomials with identical exponents therefore does not change the set, which proves that ${\texttt{compact}}(\mathcal{PZ})\equiv\mathcal{PZ}$. Furthermore, since the number of unique columns $k$ of matrix $E$ is smaller than the number of overall columns $h$, the matrices $\overline{E}$ and $\overline{G}$ are smaller or of equal size compared to the matrices $E$ and $G$, which results in a compressed representation of the set.

Complexity: The operation uniqueColumns in combination with the construction of the sets $\mathcal{H}_{j}$ can be efficiently implemented by sorting the matrix columns followed by an identification of identical neighbors, which can be realized with a worst-case complexity of $\mathcal{O}(ph\log(h))$ [39, Chapter 5]. Since all other operations have at most quadratic complexity, the overall complexity is $\mathcal{O}(n^{2}\log(n))$ using (2). ∎

2.2 Conversion from other Set Representations

This subsection shows how other set representations can be converted to SPZs.

2.2.1 Zonotope and Interval

We first provide the definition of a zonotope:

Definition 2:

(Zonotope) [31, Def. 1] Given a center $c\in\mathbb{R}^{n}$ and a generator matrix $G\in\mathbb{R}^{n\times l}$, a zonotope is defined as

[TABLE]

For a compact notation, we introduce the shorthand $\mathcal{Z}=\langle c,G\rangle_{Z}$. Any zonotope can be converted to an SPZ:

Proposition 3:

(Conversion Zonotope) A zonotope $\mathcal{Z}=\langle c,G\rangle_{Z}$ can be represented by an SPZ:

[TABLE]

The complexity of the conversion is $\mathcal{O}(l)$, with $l$ denoting the number of zonotope generators.

Proof.

If we insert $E=[\mathbf{0}^{(n,1)}~{}I_{l}]$ and $G_{I}=[~{}]$ into (1), we obtain a zonotope (see (3)).

Complexity: The construction of the matrices only involves concatenations, and therefore has complexity $\mathcal{O}(1)$. Generation of $l$ unique identifiers has complexity $\mathcal{O}(l)$, resulting in an overall complexity of $\mathcal{O}(1)+\mathcal{O}(l)=\mathcal{O}(l)$. ∎

Since any interval can be represented as a zonotope [4, Prop. 2.1], their conversion to an SPZ is straightforward.

2.2.2 Polytope

We start with the vertex-representation of a polytope:

Definition 3:

(Polytope) [4, Def. 2.2] Given $r$ polytope vertices $v_{i}\in\mathbb{R}^{n},~{}i\in\{1,\dots,r\}$, a polytope $\mathcal{P}$ is defined as

[TABLE]

We use the shorthand $\mathcal{P}=\langle[v_{1}~{}\dots~{}v_{r}]\rangle_{P}$.

Theorem 1:

(Conversion Polytope) Every bounded polytope can be equivalently represented as an SPZ.

Proof.

As shown in Def. 3, every bounded polytope $\mathcal{P}$ can be described as the convex hull of its vertices. Each vertex $v_{i}\in\mathbb{R}^{n}$ can be equivalently represented as an SPZ without independent generators $v_{i}=\langle v_{i},[~{}],0,{\texttt{uniqueID}}(1)\rangle_{PZ}$. Since SPZs without independent generators are closed under the convex hull operation, as we show later in Prop. 14, computing the convex hull of all vertices results in an SPZ that is equivalent to the bounded polytope $\mathcal{P}$. ∎

An algorithm for the conversion of a polytope in vertex-representation to an SPZ is provided in [40, Alg. 1].

2.2.3 Taylor Model

First, we formally define multi-dimensional Taylor models:

Definition 4:

(Taylor Model) [21, Def. 2.1] Given a vector field $w:\mathbb{R}^{s}\to\mathbb{R}^{n}$, where each sub-function $w_{(i)}:\mathbb{R}^{s}\to\mathbb{R}$ is a polynomial function defined as

[TABLE]

and an interval $\mathcal{I}\subset\mathbb{R}^{n}$, a Taylor model $\mathcal{T}(x)\subset\mathbb{R}^{n}$ is defined as

[TABLE]

where $x=[x_{1}~{}\dots~{}x_{s}]^{T}$, $E_{i}\in\mathbb{N}_{0}^{s\times m_{i}}$ represents an exponent matrix and $b_{i,j}\in\mathbb{R}$ are the polynomial coefficients.

For a concise notation, we introduce the shorthand $\mathcal{T}(x)=\langle w(x),\mathcal{I}\rangle_{T}$. The set defined by any Taylor model can be converted to an SPZ:

Proposition 4:

(Conversion Taylor Model) The set defined by a Taylor model $\mathcal{T}(x)=\langle w(x),\mathcal{I}\rangle_{T}$ on the domain $x\in\mathcal{D}$ with $\mathcal{D}=[l_{D},u_{D}]$ and $\mathcal{I}=[l_{R},u_{R}]$ can be equivalently represented by an SPZ:

[TABLE]

[TABLE]

The coefficients $\overline{b}_{i,j}$ and the matrices $\overline{E}_{i}$ result from the definition

[TABLE]

where $w_{(i)}(\cdot)$ is defined as in (4) and

[TABLE]

The compact operation is applied to remove monomials with an identical variable part. The complexity of the conversion is $\mathcal{O}(n^{4+n}\log(n))$.

Proof.

The auxiliary variables $\delta_{k}(\alpha_{k})$ (see (7)) represent the domain $\mathcal{D}$ with dependent factors $\alpha_{k}\in[-1,1]$:

[TABLE]

Furthermore, the interval $\mathcal{I}=[l_{R},u_{R}]$ can be equivalently represented as a zonotope $\mathcal{I}=\langle 0.5(l_{R}+u_{R}),0.5~{}{\texttt{diag}}(u_{R(1)}-l_{R(1)},\dots,u_{R(n)}-l_{R(n)})\rangle_{Z}$ (see [4, Prop. 2.1]), where the operator diag returns a diagonal matrix. The set defined by the Taylor model $\mathcal{T}(x)$, $x\in\mathcal{D}$ can therefore be equivalently expressed as

[TABLE]

where $\delta(\alpha)=[\delta_{1}(\alpha_{1})~{}\dots~{}\delta_{s}(\alpha_{s})]^{T}$ and $o=\mathbf{0}^{(n-1,1)}$.

Complexity: Let $m=\max(m_{1},\dots,m_{n})$ and $\epsilon=\max(\linebreak[1]\max(E_{1}),\dots,\max(E_{n}))$, where $\max(A)$ returns the maximum entry of a matrix $A$. Since $\delta_{k}(\alpha_{k})=c_{0}+c_{1}\alpha_{k}$, $c_{0},c_{1}\in\mathbb{R}$ (see (7)), it holds that $\delta_{k}(\alpha_{k})^{\epsilon}$ is a polynomial in $\alpha_{k}$ with $\epsilon+1$ polynomial terms. Naive evaluation without intermediate simplification of the function $w_{(i)}\left(\delta_{1}(\alpha_{1}),\dots,\delta_{s}(\alpha_{s})\right)$ with $w_{(i)}(\cdot)$ as defined in (4) therefore results in a multivariate polynomial with $\overline{m}=m(\epsilon+1)^{s}$ terms in the worst-case. From (5) and (6) it can be deduced that the exponent matrix $\widehat{E}$ consequently consists of at most $\widehat{h}=n\overline{m}=nm(\epsilon+1)^{s}$ columns. Since the complexity of compact is $\mathcal{O}(\widehat{p}\widehat{h}\log(\widehat{h}))$ (see Prop. 2) and $\widehat{p}=s$, the subsequent application of the compact operation has complexity $\mathcal{O}(snm(\epsilon+1)^{s}\log(nm(\epsilon+1)^{s}))=\mathcal{O}(snm(\epsilon+1)^{s}(\log(n)+\log(m)+s\log(\epsilon+1)))=\mathcal{O}(n^{4+n}\log(n))$ using (2) and the fact that $s=c_{s}n$, $m=c_{m}n$, and $\epsilon=c_{\epsilon}n$ with $c_{s},c_{m},c_{\epsilon}\in\mathbb{R}_{\geq 0}$ holds. This is also the overall complexity of the conversion, since all other operations have a lower complexity. ∎

2.3 Enclosure by other Set Representations

For computational reasons many algorithms that compute with sets require to enclose sets by simpler set representations. In this subsection we therefore describe how SPZs can be enclosed by other set representations. We consider the SPZ

[TABLE]

as a running example to demonstrate the tightness of the enclosure.

2.3.1 Zonotope

We first show how an SPZ can be enclosed by a zonotope:

Proposition 5:

(Zonotope) Given an SPZ $\mathcal{PZ}=\langle G,G_{I},\linebreak[1]E,id\rangle_{PZ}$, the operation zono returns a zonotope that encloses $\mathcal{PZ}$ with complexity $\mathcal{O}(n^{2})$:

[TABLE]

where $x~{}\mathrm{mod}~{}y,~{}x,y\in\mathbb{R}$ is the modulo operation.

Proof.

We over-approximate all monomial variable parts in (1) with additional independent factors, which yields a zonotope (3). Since monomials with exclusively even exponents ($i\in\mathcal{H}$) are strictly positive, we can enclose them tighter using

[TABLE]

For all other monomials ($i\in\mathcal{K}$), evaluation of the monomial variable part directly results in the interval $[-1,1]$. A dependent factor affects all monomials that contain the dependent factor. Since the over-approximation of the monomial variable parts with new independent factors destroys this dependence between different monomials (e.g., $\{\alpha_{1}\alpha_{2}^{2}+\alpha_{1}^{3}\alpha_{2}|\alpha_{1},\alpha_{2}\in[-1,1]\}\subseteq\{\beta_{1}+\beta_{2}|\beta_{1},\beta_{2}\in[-1,1]\}$), the resulting zonotope encloses $\mathcal{PZ}$, because removing dependence results in an over-approximation [37].

Complexity: The calculation of the set $\mathcal{H}$ has complexity $\mathcal{O}(ph)$, and the construction of the zonotope is $\mathcal{O}(nh)$ in the worst-case where all exponents are exclusively even, resulting in an overall complexity of $\mathcal{O}(ph)+\mathcal{O}(nh)$, which is $\mathcal{O}(n^{2})$ using (2). ∎

The enclosing zonotope for the SPZ in (8) calculated according to Prop. 5 is visualized in Fig. 3 (left).

2.3.2 Polytope

Next, we show how to enclose an SPZ with a polytope:

Proposition 6:

(Polytope) Given an SPZ $\mathcal{PZ}=\langle G,G_{I},E,id\rangle_{PZ}$, the operation poly returns a polytope that encloses $\mathcal{PZ}$ with complexity $\mathcal{O}(2^{n^{2}})$:

[TABLE]

where the vertex-representation of the polytope $\langle[v_{1}~{}\dots~{}v_{r}]\rangle_{P}$ is computed by applying [40, Alg. 2] to the SPZ

[TABLE]

Proof.

Alg. 2 in [40] computes an enclosing polytope for an SPZ for which the exponent matrix has only zeros or ones as entries. We therefore first split $\mathcal{PZ}$ into one part $\langle G_{(\cdot,\mathcal{K})},G_{I},E_{(\cdot,\mathcal{K})},id\rangle_{PZ}$ with only zeros or ones in the exponent matrix, and one remainder part $\langle G_{(\cdot,\mathcal{H})},[~{}],E_{(\cdot,\mathcal{H})},id\rangle_{PZ}$. In order to remove exponents that are greater than one the remainder part is enclosed by a zonotope (see (9)) using Prop. 5. Combination of the two parts finally yields the SPZ $\overline{\mathcal{PZ}}$ which satisfies

[TABLE]

Complexity: The calculation of the sets $\mathcal{H}$ and $\mathcal{K}$ in (9) has complexity $\mathcal{O}(ph)$, and the computation of an enclosing zonotope has according to Prop. 5 complexity $\mathcal{O}(n^{2})$. According to [40, Prop. 7] the complexity of [40, Alg. 2] is $\mathcal{O}((2^{p})^{\lfloor n/2\rfloor+1}+4^{p}(p+n))$. Using (2), the overall complexity is therefore $\mathcal{O}(2^{n^{2}})$. ∎

The enclosing polytope for the SPZ in (8) calculated according to Prop. 9 is visualized in Fig. 3 (middle).

2.3.3 Support Function, Interval, and Template Polyhedra

Let us first derive the support function of an SPZ.

Definition 5:

(Support Function) [32, Def. 1] Given a set $\mathcal{S}\subset\mathbb{R}^{n}$ and a direction $d\in\mathbb{R}^{n}$, the support function $s_{\mathcal{S}}:\mathbb{R}^{n}\to\mathbb{R}$ of $\mathcal{S}$ is defined as

[TABLE]

If $\mathcal{S}$ is convex, its support function is an exact representation; otherwise, an over-approximation is returned. Since SPZs are non-convex in general, one can only over-approximate them by support functions. To compute support functions for SPZs, we need to introduce the range bounding operation. Given a function $f:\mathbb{R}^{m}\to\mathbb{R}$ and an interval $\mathcal{I}\subset\mathbb{R}^{m}$, the range bounding operation

[TABLE]

returns an over-approximation of the exact bounds.

Proposition 7:

(Support Function) An SPZ $\mathcal{PZ}=\langle G,G_{I},E,id\rangle_{PZ}$ is over-approximated by the support function

[TABLE]

where the vector $d\in\mathbb{R}^{n}$ specifies the direction. The calculation of $\widehat{s}_{\mathcal{PZ}}(d)$ has complexity $\mathcal{O}(n^{2})+\mathcal{O}({\texttt{B}})$, where $\mathcal{O}({\texttt{B}})$ denotes the computational complexity of the range bounding operation.

Proof.

We first project the SPZ onto the direction $d$, and then divide the one-dimensional projected SPZ into one part with dependent generators and one with independent generators:

[TABLE]

The bounds for the independent part calculated by the sum of absolute values are exact [32, Sec. 2]. However, the lower bound $l$ and the upper bound $u$ of the dependent part are over-approximative since the range bounding operation B returns an over-approximation, so that $\widehat{s}_{\mathcal{PZ}}(d)$ over-approximates $\mathcal{PZ}$.

Complexity: The calculation of the projection onto $d$ has a complexity of $\mathcal{O}(nh)+\mathcal{O}(nq)$ (see Prop. 8), which results in an overall complexity of $\mathcal{O}(n^{2})+\mathcal{O}({\texttt{B}})$ using (2) since all other operations have linear complexity. ∎

Note that the tightness of $\widehat{s}_{\mathcal{PZ}}(d)$ solely depends on the tightness of the bounds of the function $w(\cdot)$ obtained by one of the range bounding techniques, e.g., interval arithmetic [37], Bernstein polynomials [19], and verified global optimization [50]. A comparison of different techniques can be found in [9].

A template polyhedron enclosing an SPZ can easily be constructed by evaluating the support function $\widehat{s}_{\mathcal{PZ}}(d)$ for a discrete set of directions $\mathcal{D}=\{d_{1},\dots,d_{r}\},~{}d_{i}\in\mathbb{R}^{n},~{}i=1\dots r$. The over-approximation with an interval represents a special case where $\mathcal{D}=\{I_{n(\cdot,1)},\linebreak[1]\dots,\linebreak[1]I_{n(\cdot,n)},\linebreak[1]-I_{n(\cdot,1)},\linebreak[1]\dots,\linebreak[1]-I_{n(\cdot,n)}\}$. The enclosing interval for the SPZ in (8) calculated by using Bernstein polynomials for range bounding is visualized in Fig. 3 (right).

2.4 Basic Set Operations

This subsection derives basic operations on SPZs.

2.4.1 Multiplication with a Matrix

The left-multiplication with a matrix is obtained as:

Proposition 8:

(Multiplication) Given an SPZ $\mathcal{PZ}=\langle G,\linebreak[1]G_{I},E,id\rangle_{PZ}\subset\mathbb{R}^{n}$ and a matrix $M\in\mathbb{R}^{m\times n}$, the left-multiplication is computed as

[TABLE]

which has complexity $\mathcal{O}(mn^{2})$.

Proof.

The result follows directly from inserting the definition of SPZs in (1) into the definition of the operator $\otimes$ (see Sec. 1.3).

Complexity: The complexity results from the complexity of matrix multiplications and is therefore $\mathcal{O}(mnh)+\mathcal{O}(mnq)=\mathcal{O}(mn^{2})$ using (2). ∎

2.4.2 Minkowski Addition

Even though every zonotope can be represented as an SPZ, we provide a separate definition for the Minkowski addition of an SPZ and a zonotope for computational reasons.

Proposition 9:

(Addition) Given two SPZs, $\mathcal{PZ}_{1}=\langle G_{1},\linebreak[1]G_{I,1},E_{1},id_{1}\rangle_{PZ}$ and $\mathcal{PZ}_{2}=\langle G_{2},G_{I,2},$ $E_{2},id_{2}\rangle_{PZ}$, as well as a zonotope $\mathcal{Z}=\langle c_{z},G_{z}\rangle_{Z}$, their Minkowski sum is defined as

[TABLE]

[TABLE]

where (10) has complexity $\mathcal{O}(n)$ and (11) has complexity $\mathcal{O}(1)$.

Proof.

The result is obtained by inserting the definition of zonotopes (3) and SPZs (1) into the definition of the Minkowski sum (see Sec. 1.3). For two SPZs, we generate new identifiers for all factors since the Minkowski sum per definition removes all dependencies between the two added sets.

Complexity: The construction of the resulting SPZs only involves concatenations and therefore has complexity $\mathcal{O}(1)$. For two SPZs, $p_{1}+p_{2}$ unique identifiers have to be generated, resulting in the complexity $\mathcal{O}(p_{1}+p_{2})$, which is $\mathcal{O}(n)$ using (2). ∎

2.4.3 Exact Addition

By computation of the Minkowski sum of two SPZs, $\mathcal{PZ}_{1}$ and $\mathcal{PZ}_{2}$, as defined in (10), possible dependencies between $\mathcal{PZ}_{1}$ and $\mathcal{PZ}_{2}$ due to common dependent factors get lost. We therefore introduce the exact addition $\mathcal{PZ}_{1}\boxplus\mathcal{PZ}_{2}$ of two SPZs, which explicitly considers the dependencies between $\mathcal{PZ}_{1}$ and $\mathcal{PZ}_{2}$. To bring the exponent matrices to a common representation, we apply mergeID prior to the computation.

Proposition 10:

(Exact Addition) Given two SPZs, $\mathcal{PZ}_{1}=\langle G_{1},\linebreak[1]G_{I,1},E_{1},id\rangle_{PZ}$ and $\mathcal{PZ}_{2}=\langle G_{2},G_{I,2},$ $E_{2},id\rangle_{PZ}$ with a common identifier vector $id$, their exact addition is defined as

[TABLE]

which has complexity $\mathcal{O}(n^{2}\log(n))$. The compact operation is applied to remove monomials with an identical variable part.

Proof.

The result is identical to the one for Minkowski addition of two SPZs in (10), with the difference that the common identifier vector $id$ is used instead of newly generated unique identifiers.

Complexity: Merging the identifier vectors has complexity $\mathcal{O}(n^{2})$ (see Prop. 1). The construction of the resulting SPZ only involves concatenations and therefore has complexity $\mathcal{O}(1)$. Subsequent application of the compact operation has complexity $\mathcal{O}(p_{1}(h_{1}+h_{2})\log(h_{1}+h_{2}))$ (see Prop. 2), which is $\mathcal{O}(n^{2}\log(n))$ using (2). The overall complexity is therefore $\mathcal{O}(n^{2})+\mathcal{O}(1)+\mathcal{O}(n^{2}\log(n))=\mathcal{O}(n^{2}\log(n))$. ∎

As demonstrated later in Sec. 3, using exact addition instead of the Minkowski sum for reachability analysis results in a significantly tighter enclosure of the reachable set.

2.4.4 Cartesian Product

Even though every zonotope can be represented as an SPZ, we provide a separate definition for the Cartesian product of an SPZ and a zonotope for computational reasons.

Proposition 11:

(Cartesian Product) Given two SPZs, $\mathcal{PZ}_{1}=\langle G_{1},\linebreak[1]G_{I,1},E_{1},id_{1}\rangle_{PZ}\subset\mathbb{R}^{n}$ and $\mathcal{PZ}_{2}=\langle G_{2},G_{I,2},$ $E_{2},id_{2}\rangle_{PZ}\subset\mathbb{R}^{m}$, as well as a zonotope $\mathcal{Z}=\langle c_{z},G_{z}\rangle_{Z}\subset\mathbb{R}^{m}$, their Cartesian product is defined as

[TABLE]

[TABLE]

where (12) has complexity $\mathcal{O}(n)$ and (13) has complexity $\mathcal{O}(1)$.

Proof.

The result is obtained by inserting the definition of zonotopes (3) and SPZs (1) into the definition of the Cartesian product (see Sec. 1.3).

Complexity: The construction of the resulting SPZs only involves concatenations and therefore has complexity $\mathcal{O}(1)$. For two SPZs, $p_{1}+p_{2}$ unique identifiers have to be generated, resulting in the complexity $\mathcal{O}(p_{1}+p_{2})$, which is $\mathcal{O}(n)$ using (2). ∎

2.4.5 Quadratic Map

For reachability analysis based on the conservative polynomialization approach [5], a polynomial abstraction of the nonlinear dynamic function is calculated, requiring quadratic and higher-order maps. Here, we derive the equations for the quadratic map.

Definition 6:

(Quadratic Map) [5, Theorem 1] Given a set $\mathcal{S}\subset\mathbb{R}^{n}$ and a discrete set of matrices $Q_{i}\in\mathbb{R}^{n\times n},i=1\dots m$, the quadratic map of $\mathcal{S}$ is defined as

[TABLE]

For SPZs, we first consider the special case without independent generators, and later present the general case.

Proposition 12:

Given the SPZ $\widehat{PZ}\linebreak[1]=\langle\widehat{G},[~{}],\widehat{E},\widehat{id}\rangle_{PZ}$ and a discrete set of matrices $Q_{i}\in\mathbb{R}^{n\times n},i=1\dots m$, the result of the quadratic map is

[TABLE]

where

[TABLE]

The compact operation is applied to remove monomials with an identical variable part. The overall complexity is $\mathcal{O}(n^{3}\log(n))+\mathcal{O}(n^{3}m)$.

Proof.

The equations are obtained by substitution of $s$ in Def. 6 with the definition of an SPZ from (1), which yields

[TABLE]

Note that only the generator matrix, but not the exponent matrix, is different for each dimension $x_{(i)}$.

Complexity: The construction of the matrices $\overline{E}_{j}$ has complexity $\mathcal{O}(\widehat{h}^{2}\widehat{p})$, and the construction of the matrices $\overline{G}_{j}$ has complexity $\mathcal{O}(n^{2}\widehat{h}m)+\mathcal{O}(n\widehat{h}^{2}m)$ if the results for $Q_{i}\widehat{G}$ are stored and reused. Since the matrices $\overline{E}$ and $\overline{G}$ both consist of $\overline{h}=\widehat{h}^{2}$ columns, and because the complexity of the compact operation is $\mathcal{O}(\widehat{p}\overline{h}\log(\overline{h}))$ (see Prop. 2), the complexity of the subsequent application of compact is $\mathcal{O}(\widehat{p}\widehat{h}^{2}\log(\widehat{h}^{2}))$. The resulting overall complexity is therefore $\mathcal{O}(n^{2}\widehat{h}m)+\mathcal{O}(n\widehat{h}^{2}m)+\mathcal{O}(\widehat{p}\widehat{h}^{2}\log(\widehat{h}^{2}))$, which is $\mathcal{O}(n^{3}\log(n))+\mathcal{O}(n^{3}m)$ using (2). ∎

We now extend Prop. 12 to the general case including independent generators, for which we compute an over-approximation for computational reasons.

Proposition 13:

(Quadratic Map) Given an SPZ $\mathcal{PZ}=\langle G,G_{I},E,id\rangle_{PZ}\subset\mathbb{R}^{n}$ and a discrete set of matrices $Q_{i}\in\mathbb{R}^{n\times n},i=1\dots m$,

[TABLE]

where the zonotope enclosure is calculated by applying Prop. 5. The matrices $\overline{G}$ and $\overline{E}$ are computed according to (14) with the extended generator and exponent matrices $\widehat{G}$ and $\widehat{E}$, as well as the extended identifier vector $\widehat{id}$ defined as

[TABLE]

so that $\widehat{p}=p+q$. The compact operation is applied to remove monomials with an identical variable part. The complexity of the calculations is $\mathcal{O}(n^{3}\log(n))+\mathcal{O}(n^{3}m)$.

Proof.

With the extended generator and exponent matrices $\widehat{G}$ and $\widehat{E}$ and the extended identifier vector $\widehat{id}$ from (16), $\mathcal{PZ}$ can be represented equivalently as an SPZ without independent generators

[TABLE]

which enables the computation of the quadratic map according to Prop. 12. For computational reasons, the resulting matrices $\overline{G}$ and $\overline{E}$ from Prop. 12 are divided into a dependent part that contains the dependent factors $\alpha_{1},\dots,\alpha_{p}$ only ($i\in\mathcal{H}$), and a remainder part that contains all remaining monomials ($i\in\mathcal{K}$):

[TABLE]

Since the part containing the remaining monomials is enclosed by a zonotope, it holds that the resulting SPZ encloses the quadratic map.

Complexity: With the extended matrices $\widehat{G}$ and $\widehat{E}$ from (16), the calculation of the matrices $\overline{E}_{j}$ and $\overline{G}_{j}$ has complexity $\mathcal{O}((h+q)^{2}(p+q))$ and $\mathcal{O}(n^{2}(h+q)m)+\mathcal{O}(n(h+q)^{2}m)$, respectively, which is $\mathcal{O}(n^{3}m)$ using (2). As for the case without independent generators, the complexity from the subsequent application of the compact operation is $\mathcal{O}(ph^{2}\log(h^{2}))=\mathcal{O}(n^{3}\log(n))$ using (2). The resulting overall complexity is therefore $\mathcal{O}(n^{3}\log(n))+\mathcal{O}(n^{3}m)$, since all other operations have lower complexity. ∎

For the extension to cubic and higher-order maps of sets, which is omitted at this point due to space limitations, the results from Prop. 13 can be reused. We demonstrate the tightness of the quadratic map enclosure with an example:

Example 2:

We consider the SPZ

[TABLE]

and the matrices

[TABLE]

A comparison between the exact quadratic map and the over-approximation computed by Prop. 13 is shown in Fig. 4.

2.4.6 Convex Hull

Another important set operation is the convex hull of two sets:

Definition 7:

(Convex Hull) Given two sets $\mathcal{S}_{1}\subset\mathbb{R}^{n}$ and $\mathcal{S}_{2}\subset\mathbb{R}^{n}$, the convex hull of them is defined as

[TABLE]

We first consider the special case without independent generators and later present the general case.

Proposition 14:

Given $\mathcal{PZ}_{1}=\langle G_{1},[~{}],E_{1},\linebreak[1]id_{1}\rangle_{PZ}$ and $\mathcal{PZ}_{2}=\langle G_{2},[~{}],$ $E_{2},id_{2}\rangle_{PZ}$, the convex hull is computed as

[TABLE]

which has complexity $\mathcal{O}(n^{2})$.

Proof.

We first generate new unique identifiers for all dependent factors to remove possible dependencies between the two SPZs. For SPZs, the definition of the convex hull from Def. 7 can be equivalently written as

[TABLE]

Evaluation of the exact additions and the Minkowski sum in (18) for the SPZs according to Prop. 10 and Prop. 9 results in the equations in (17), where we substituted the parameter $\lambda$ with an additional dependent factor $\alpha_{p_{1}+p_{2}+1}=\lambda$. Since $\lambda\in[-1,1]$ and $\alpha_{p_{1}+p_{2}+1}\in[-1,1]$, this substitution does not change the set.

Complexity: The construction of the matrix $\overline{G}$ has complexity $\mathcal{O}(2n(h_{1}+h_{2}))$. Generation of $p_{1}+p_{2}+1$ unique identifiers for $\overline{id}$ has complexity $\mathcal{O}(p_{1}+p_{2}+1)$. The overall complexity is therefore $\mathcal{O}(2n(h_{1}+h_{2}))+\mathcal{O}(p_{1}+p_{2}+1)$, which is $\mathcal{O}(n^{2})$ using (2). ∎

We now extend Prop. 14 to the general case including independent generators, for which we compute an over-approximation for computational reasons.

Proposition 15:

(Convex Hull) Given the two SPZs $\mathcal{PZ}_{1}=\langle G_{1},G_{I,1},E_{1},id_{1}\rangle_{PZ}$ and $\mathcal{PZ}_{2}=\langle G_{2},G_{I,2},$ $E_{2},id_{2}\rangle_{PZ}$,

[TABLE]

where ${\texttt{conv}}(\overline{\mathcal{PZ}}_{1},\overline{\mathcal{PZ}}_{2})$ is calculated according to Prop. 14, and an over-approximation of the convex hull of two zonotopes is computed according to [4, Eq. (2.2)]:

[TABLE]

where $\mathcal{K}_{1}=\{1,\dots,q_{1}\}$ and $\mathcal{K}_{2}=\{1,\dots,q_{2}\}$. The overall complexity is $\mathcal{O}(n^{2})$.

Proof.

Each SPZ $\mathcal{PZ}=\langle G,G_{I},E,id\rangle_{PZ}$ can be represented equivalently as the Minkowski sum of an SPZ and a zonotope $\mathcal{PZ}=\langle G,[~{}],E,id\rangle_{PZ}\oplus\langle\mathbf{0}^{(n,1)},G_{I}\rangle_{Z}$ (see (11)). Given sets $S_{1}$, $S_{2}$, $S_{3}$, $S_{4}\subset\mathbb{R}^{n}$, it holds that ${\texttt{conv}}(S_{1}\oplus S_{2},S_{3}\oplus S_{4})\subseteq{\texttt{conv}}(S_{1},S_{3})\oplus{\texttt{conv}}(S_{2},S_{4})$, which can be derived from the definition of the convex hull from Def. 7:

[TABLE]

[TABLE]

An enclosure of the convex hull of two SPZs can therefore be obtained by computing the convex hull of two SPZs without independent generators according to Prop. 14, and the computation of the convex hull of two zonotopes according to [4, Eq. (2.2)].

Complexity: As shown in Prop. 14, the complexity for computing the convex hull of two SPZs without independent generators is $\mathcal{O}(n^{2})$. The complexity for the computation of the convex hull of two zonotopes as defined in (19) is $\mathcal{O}(2n\min(q_{1},q_{2}))$, resulting from the matrix additions and subtractions that are required for the construction of the matrix $\widehat{G}_{1}$ or $\widehat{G}_{2}$. The overall complexity is therefore $\mathcal{O}(n^{2})+\mathcal{O}(2n\min(q_{1},q_{2}))$, which is $O(n^{2})$ using (2). ∎

We demonstrate the tightness of the convex hull enclosure with an example:

Example 3:

We consider the SPZs

[TABLE]

and

[TABLE]

A comparision between the exact convex hull and the over-approximation computed by Prop. 15 is shown in Fig. 5.

2.5 Auxiliary Operations

This subsection derives useful auxiliary operations on SPZs.

2.5.1 Order Reduction

Many set operations, such as Minkowski addition or quadratic maps, increase the number of generators and consequently also the order $\rho$ of the SPZ. Thus, for computational reasons, it is necessary to repeatedly reduce the zonotope order during reachability analysis. We propose a reduction operation for SPZs that is based on the order reduction of zonotopes (see e.g., [42]).

Proposition 16:

(Reduce) Given an SPZ $\mathcal{PZ}=\langle G,G_{I},\linebreak[1]E,id\rangle_{PZ}$ and a desired zonotope order $\rho_{d}\geq 1+1/n$, the operation reduce returns an SPZ with an order smaller than or equal to $\rho_{d}$ that encloses $\mathcal{PZ}$:

[TABLE]

where the zonotope enclosure is calculated by applying Prop. 5. For reduction, the

[TABLE]

smallest generators are selected:

[TABLE]

where $\widehat{G}=[G~{}G_{I}]$ and $d\in\mathbb{N}_{>0}^{h+q}$ is a vector of indices for which each next entry represents a longer generator than the previous entry. The complexity is $\mathcal{O}(n^{2})+\mathcal{O}({\texttt{reduce}})$, where $\mathcal{O}({\texttt{reduce}})$ denotes the complexity of the zonotope reduction, which depends on the selected method.

Proof.

The definition of $a$ (see (20)) ensures that $|\widehat{\mathcal{K}}|+|\widehat{\mathcal{H}}|+n+1\leq\rho_{d}n$ for values $\rho_{d}\geq 1+1/n$. Furthermore, ${\texttt{reduce}}(\mathcal{PZ},\rho_{d})\supseteq\mathcal{PZ}$ since the operators zono and reduce are both over-approximative, and therefore ${\texttt{reduce}}({\texttt{zono}}(\cdot))$ is over-approximative, too.

Complexity: Sorting the generators has a complexity of $\mathcal{O}(n(h+q))+\mathcal{O}((h+q)\log(h+q))$, which is $\mathcal{O}(n^{2})$ using (2). In the worst-case where all dependent generators get reduced, the enclosure with a zonotope has complexity $\mathcal{O}(ph)+\mathcal{O}(nh)$ (see. Prop. 5), which is $\mathcal{O}(n^{2})$ using (2). The overall complexity is therefore $\mathcal{O}(n^{2})\linebreak[1]+\mathcal{O}({\texttt{reduce}})$. ∎

After reduction, we remove possibly generated all-zero rows in the exponent matrix.

2.5.2 Restructure

Due to the repeated order reduction and Minkowski addition during reachability analysis, the volume spanned by independent generators grows relative to the volume spanned by dependent generators. As explained later in Sec. 3, this has a negative effect on the tightness of the reachable sets. We therefore define the operation restructure, which introduces new dependent generators that over-approximate the independent ones:

Proposition 17:

(Restructure) Given an SPZ $\mathcal{PZ}=\langle G,\linebreak[1]G_{I},E,id\rangle_{PZ}$, restructure returns an SPZ that encloses $\mathcal{PZ}$ and removes all independent generators:

[TABLE]

The overall complexity is $\mathcal{O}(n)+\mathcal{O}({\texttt{reduce}})$, where $\mathcal{O}({\texttt{reduce}})$ is the complexity of the zonotope reduction.

Proof.

The result of the restructure operation encloses the original set since reduce is over-approximative, and the redefinition of independent generators as new dependent generators just changes the set representation, but not the set in (21) itself:

[TABLE]

Complexity: The generation of $n$ new unique identifiers has complexity $\mathcal{O}(n)$. Since the construction of the matrices only involves concatenations, the overall complexity equals $\mathcal{O}(n)+\mathcal{O}({\texttt{reduce}})$, where $\mathcal{O}({\texttt{reduce}})$ is the complexity of the zonotope reduction. ∎

We demonstrate the effectiveness of Prop. 17 by numerical examples in Sec. 4. To save computation time, we define an upper bound $p_{d}$ of factors for the SPZ after restructuring so that $p+n\leq p_{d}$ holds, where independent factors are removed first to maintain as much dependence as possible. If required, dependent factors are removed by enclosing the corresponding generators with a zonotope (see Prop. 5).

3 Reachability Analysis

In this section, we demonstrate how SPZs can be used to improve reachability analysis for nonlinear systems. We consider nonlinear systems of the form

[TABLE]

where $x$ is the vector of system states, $u$ is the input vector, and $f:~{}\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}^{n}$ is a Lipschitz continuous function. The reachable set of the system is defined as follows:

Definition 8:

(Reachable Set) Let $\xi(t,x_{0},u(\cdot))$ denote the solution to (22) for an initial state $x(0)=x_{0}$ and the input trajectory $u(\cdot)$. The reachable set for an initial set $\mathcal{X}_{0}\subset\mathbb{R}^{n}$ and a set of possible input values $\mathcal{U}\subset\mathbb{R}^{m}$ is

[TABLE]

The superscript $e$ on $\mathcal{R}_{\mathcal{X}_{0}}^{e}(t)$ denotes the exact reachable set, which cannot be computed for general nonlinear systems. Therefore, we compute a tight over-approximation $\mathcal{R}(t)\supseteq\mathcal{R}_{\mathcal{X}_{0}}^{e}(t)$. Furthermore, we calculate the reachable set for consecutive time intervals $\tau_{s}=[t_{s},t_{s+1}]$ with $t_{s+1}=t_{s}+\Delta t$ so that the reachable set for a time horizon $t_{f}$ is given as $\mathcal{R}([0,t_{f}])=\bigcup_{s=0}^{t_{f}/\Delta t-1}\mathcal{R}(\tau_{s})$, where $t_{f}$ is a multiple of $\Delta t$.

3.1 Reachability Algorithm

Our algorithm in based on the conservative polynomialization approach in [5], which is an extension to the conservative linearization approach in [12]. The principle of conservative linearization is quite simple: In each time interval $\tau_{s}$ the nonlinear function $f(\cdot)$ is linearized and the set of linearization errors is treated as an additional uncertain input to the system, which then allows to calculate the reachable set with a reachability algorithm for linear systems. The conservative polynomialization approach improves this method by abstracting the nonlinear function $f(\cdot)$ by a Taylor expansion of order $\kappa$:

[TABLE]

where $z=[x^{T}~{}u^{T}]^{T}$, the set $\mathcal{L}_{(i)}(t)$ is the Lagrange remainder, defined in [5, Eq. (2)], and the vector $z^{*}\in\mathbb{R}^{n+m}$ is the expansion point for the Taylor series.

In order to fully exploit the advantages of SPZs, Alg. 1 is slightly modified from [5]. We only specify the algorithm for the Taylor order $\kappa=2$ for simplicity, since the extension to higher orders is straightforward. The operation taylor returns the matrices

[TABLE]

storing the coefficients of the Taylor series expansion for the nonlinear function $f(\cdot)$ at the expansion point $z^{*}$, and the operation enlarge enlarges a set by a given scalar factor $\lambda\in\mathbb{R}_{>1}$. The definitions of the operations ${\texttt{post}}^{\Delta}$ [5, Sec. 4.1], varInputs [5, Sec. 4.2], and lagrangeRemainder [5, Sec. 4.1] are identical to the ones in [5]. The post operator changed since we precompute some of the sets in our algorithm, and since we use the exact addition as defined in Prop. 10 instead of the Minkowski addition to add $\mathcal{F}_{1}$ to $\mathcal{F}_{2}$:

[TABLE]

where $\Gamma(\Delta t)$ is defined as in [5, Sec. 3.2], and $\mathcal{R}^{p,\Delta}(\cdot)$ is defined as in [5, Eq. (9)]. We heuristically trigger the restructure process (see Sec. 2.5.2) when

[TABLE]

where $\mathcal{PZ}=\langle G,G_{I},E,id\rangle_{PZ}$, interval computes an interval enclosure, and volume calculates the volume of a multi-dimensional interval.

The while-loop in lines 3-24 of Alg. 1 iterates over all time intervals $\tau_{s}$ until the time horizon $t_{f}$ is reached. For each time interval we first abstract the nonlinear equation $f(\cdot)$ by a Taylor expansion in Line 4. In lines 5-16 we then compute the set of linearization errors $\Psi(\tau_{s})$ on the time interval reachable set $\mathcal{R}(\tau_{s})$. The problem we are facing here is that we need $\mathcal{R}(\tau_{s})$ to calculate $\Psi(\tau_{s})$, but on the other hand we also need $\Psi(\tau_{s})$ to calculate $\mathcal{R}(\tau_{s})$. To resolve this mutual dependence we first compute $\mathcal{R}(\tau_{s})$ using the set of linearization errors from the previous time step ($\Psi(\tau_{s+1})=\Psi(\tau_{s})$, see Line 23) as an initial guess in Line 13. Next, we use $\mathcal{R}(\tau_{s})$ to calculate $\Psi(\tau_{s})$ in Line 15. In the next iteration of the repeat-until loop we then enlarge the set of linearization errors in Line 9 and calculate $\mathcal{R}(\tau_{s})$ and $\Psi(\tau_{s})$ using the enlarged set $\overline{\Psi}(\tau_{s})$. We repeat this process until the enlarged set of linearization errors $\overline{\Psi}(\tau_{s})$ contains the set of actual linearization erros $\Psi(\tau_{s})$ (see Line 16), which guarantees that $\mathcal{R}(\tau_{s})$ and $\Psi(\tau_{s})$ are both over-approximative. For the benchmarks in Sec. 4 this procedure results in 1 to 3 iterations of the repeat-until loop. After we computed the set of linearization errors, we finally calculate the reachable set $\mathcal{R}(t_{s+1})$ for the next point in time in Line 17, which we then use as the new initial set for the next time interval $\tau_{s+1}$.

At the end of each time step, we apply the operation reduce in Line 18 to reduce the zonotope order to the desired order $\rho_{d}$. Since the zonotope order is defined as $\rho=\frac{h+q}{n}$ (see Sec. 2), this ensures that the SPZs that represent the reachable set contain at most $\rho_{d}n$ generators, where $n$ is the system dimension. While order reduction limits the growth of the size of the exponent matrix, the growth of the integer entries in the exponent matrix is not contolled explicitly. However, since the generators that belong to large exponents are usually small they are automatically removed by the operation reduce.

Since Alg. 1 contains several tuning parameters we shortly discuss their influence on the performance. The parameter with the largest effect is the time step size $\Delta t$; a smaller time step size improves the accuracy, but also prolongs the computation time. In addition, the maximum volume ratio $\mu_{d}$ can have a large impact. A smaller ratio potentially increases the accuracy, but a very small ratio can have a negative effect since every restructure operation results in an over-approximation. All other parameters have a smaller influence on the performance. Good default values are $\lambda=0.1$ for the enlargement factor, $\rho_{d}=50$ for the maximum zonotope order, and $p_{d}=50$ for the maximum number of dependent factors. Furthermore, the method in [31, Sec. 3.4] (Girard’s method) is applied for zonotope reduction, and we use principal-component-analysis-based order reduction in combination with Girard’s method [42, Sec. III.A] for the reduction during the restructure operation. For the future we aim to develop strategies to tune the parameters automatically.

We proceed with a discussion of the main advantages resulting from SPZs.

3.2 Advantages of using Sparse Polynomial Zonotopes

As mentioned earlier, one of the main advantages of SPZs is that they reduce the dependency problem in Alg. 1. We demonstrate this with a short example:

Example 4:

We consider the one-dimensional system $\dot{x}=f(x)=-x+x^{2}$, the initial set $\mathcal{R}(0)=\{\alpha_{1}|\alpha_{1}\in[-1,1]\}$, and the time step size $\Delta t=1$. Computation of the Taylor expansion at $z^{*}=0$ in Line 4 of Alg. 1 results in the parameter values $w=f(z^{*})=0$, $A=\frac{\partial f}{\partial x}|_{z^{*}}=-1$, and $D=\frac{\partial^{2}f}{\partial x^{2}}|_{z^{*}}=2$. The quadratic map in Line 6 evaluates to $\frac{1}{2}{\texttt{sq}}(D,\mathcal{R}(0))=\{\alpha_{1}^{2}|\alpha_{1}\in[-1,1]\}$ for SPZs. On the other hand, if we use zonotopes, then the quadratic map has to be over-approximated with $\frac{1}{2}{\texttt{sq}}(D,\mathcal{R}(0))=\{0.5+0.5\alpha_{2}|\alpha_{2}\in[-1,1]\}$. Furthermore, the exact addition as defined in Prop. 10 is not possible for zonotopes. Therefore, the sets $\mathcal{F}_{1}$ and $\mathcal{F}_{2}$ in (23) have to be added using the Minkowski sum, which results in an additional over-approximation due to the loss of dependency. With zonotopes, we obtain for (23)

[TABLE]

With SPZs, however, we obtain the exact set

[TABLE]

Using zonotopes for reachability analysis therefore leads to a significant over-approximation error in each time step. A similar problem occurs with the polynomial zonotope representation from [5], since this requires limiting the maximum polynomial degree in advance.

3.3 Hybrid Systems

In reachability analysis for hybrid systems, the main difficulty is the calculation of the intersection between the reachable set and the guard sets. For SPZs, three different strategies exist:

We calculate the intersection with a zonotope over-approximation of the SPZ. By doing so, it is possible to directly apply the well-developed techniques for the computation of guard intersections with zonotopes, like e.g., the ones from [10, Sec. VI] or [33]. 2. 2.

We calculate the intersections with the guard sets by using the guard-mapping approach in [11] or the time-scaling approach in [13]. Both methods require basic set operations only and can therefore by applied to SPZs. 3. 3.

We apply the method in [41] which tightly encloses the intersection of the reachable set with guard sets represented by nonlinear level sets with an SPZ.

Which method performs best depends on the system.

4 Numerical Examples

In this section we demonstrate the improvements to reachability analysis due to using SPZs on four benchmark systems. The computations for our approach are carried out in MATLAB on a 2.9GHz quad-core i7 processor with 32GB memory. Our implementation of SPZs will be made publicly available with the next release of the CORA toolbox [6].

4.1 Van-der-Pol Oscillator

The system considered first is the Van-der-Pol oscillator taken from the 2019 ARCH competition [36, Sec. 3.1]:

[TABLE]

For this system, we compare the results for the computation of the reachable set with Alg. 1 using zonotopes, the quadratic zonotopes from [5], and our SPZ representation. We consider the initial set $x_{1}\in[1.23,1.57]$ and $x_{2}\in[2.34,2.46]$, and use the parameter values $\Delta t=0.005$ s, $\rho_{d}=50$, $\lambda=0.1$, $\mu_{d}=0.01$, and $p_{d}=100$. For a fair comparison, we use the same parameter values for every set representation.

The resulting reachable sets are shown in Fig. 6. It is clearly visible that the stability of the limit cycle can only be verified with SPZs when sets are not split. The computation time is $9.33$ seconds for zonotopes, $13.38$ seconds for quadratic zonotopes, and $16.52$ seconds for SPZs.

An impression on how tight the reachable set can be over-approximated with SPZs is provided in Fig. 7, where the reachable set after $t=3.15$ seconds computed with a time step size of $\Delta t=0.0001$ s and a maximum volume ratio of $\mu_{d}=0.001$ is compared to the exact reachable set of the system. The figure also demonstrates how well the SPZ approximates the shape of the exact reachable set.

4.2 Drivetrain

For the second numerical example, we examine a drivetrain, which is also a benchmark from the ARCH 2019 competition [7, Sec. 3.3]. We consider the case with two rotating masses, resulting in a system dimension of $n=11$. The model is a hybrid system with linear dynamics. However, we apply the novel approach from [13] for calculating the intersections with guard sets, which is based on time-triggered conversion of guards and results in a significant nonlinearity due to the time-scaling process. The initial set is given by $\mathcal{R}(0)=0.5(\mathcal{X}_{0}-{\texttt{center}}(\mathcal{X}_{0}))+{\texttt{center}}(\mathcal{X}_{0})$, where $\mathcal{X}_{0}$ is defined as in [7, Sec. 3.3], and we consider the same extreme acceleration maneuver as in [7, Sec. 3.3]. As a specification, we require that the engine torque after $1.5$ seconds is at least $59Nm$, which can be formally specified as $T_{m}\geq 59Nm~{}\forall t\geq 1.5s$.

The results for the drivetrain model are shown in Fig. 8. We explicitly considered the possibility of splitting the reachable sets along the largest generator vector so that the specification could be verified with all set representations. However, splitting sets prolongs the computation time: with quadratic zonotopes, the verification took $93$ seconds, and $221$ seconds with zonotopes. Only with SPZs was it possible to verify the specification without splitting, resulting in a computation time of $15$ seconds, which is six times faster than with quadratic zonotopes and more than $14$ times faster than with zonotopes. Compared to other non-zonotopic set representations, the speed-up is even larger.

4.3 Spacecraft Rendezvous

As a third numerical example we consider the docking-maneuver of a spacecraft taken from the ARCH 2019 competition [36, Sec. 3.4]. The model is a hybrid system with $n=4$ states and nonlinear dynamics. The three discrete modes are approaching, rendezvous attempt, and aborting. We consider the same initial set and the same specifications as in [36, Sec. 3.4]. The specifications are that in mode rendezvous attempt the spacecraft is located inside the line-of-sight cone and the absolute velocity stays below 3.3 m/min. Furthermore, in mode aborting the spacecraft should not collide with the space station.

We apply Alg. 1 with the parameter values $\Delta t=0.2$ min (mode approaching and abortion), $\Delta t=0.05$ min (mode rendezvous attempt), $\rho_{d}=10$, $\lambda=0.1$, $\mu_{d}=1$, and $p_{d}=10$. To calculate the intersection between the reachable set and the guard sets we use the method in [41]. The resulting reachable satisfies all specifications. To compare the performance of SPZs with other reachability tools we consider the results from the ARCH 19 competition [36]. The comparison in Tab. 2 shows that using SPZs resulted in the smallest computation time.

4.4 Transcriptional Regulator Network

To demonstrate the scalability of our approach, we consider the benchmark in [47, Sec. VIII.D] describing a transcriptional regulator network with $N$ genes. For a network with $N$ genes the system has $n=2N$ dimensions. We consider the case without artificial guard set so that the benchmark represents a continuous nonlinear system with uncertain inputs. Furthermore, we consider the same initial set, time horizon, and set of uncertain inputs as in [47, Sec. VIII.D].

We compute the reachable set with SPZs using Alg. 1 with the parameter values $\Delta t=0.1$ min, $\rho_{d}=10$, $\lambda=0.1$, $\mu_{d}=1$, and $p_{d}=50$. The reachable set is visualized in Fig. 9, and the computation times for different system dimensions are listed in Tab. 3. Even for a system dimension of $n=48$ the computation of reachable set with SPZs takes only $122$ seconds, which demonstrates how well our approach scales with the system dimension.

5 Conclusions

We have introduced sparse polynomial zonotopes, a new non-convex set representation. The sparsity results in several advantages compared to previous representations of polynomial zonotopes: sparse polynomial zonotopes enable a compact representation of sets, they are closed under all relevant set operations, and all operations have at most polynomial complexity. The fact that sparse polynomial zonotopes are a generalization of several other set representations like Taylor models, polytopes, and zonotopes further substantiates the relevance of this new representation.

The main application for sparse polynomial zonotopes is reachability analysis for nonlinear systems. Numerical examples demonstrate that using sparse polynomial zonotopes results in much tighter over-approximations of reachable sets compared to using zonotopes or quadratic zonotopes. Due to the improved accuracy, splitting can be avoided, resulting in a significant reduction of the computation time, since splitting of sets results in an exponential number of sets to be propagated with respect to the system dimension.