# Reachability Analysis for Linear Systems with Uncertain Parameters using
Polynomial Zonotopes

Yushen Huang111These authors contributed equally to this work.
<EMAIL_ADDRESS>Ertai Luo22footnotemark: 2<EMAIL_ADDRESS>Stanley Bak<EMAIL_ADDRESS>Yifan Sun<EMAIL_ADDRESS>Stony
Brook University, Stony Brook, NY, USA

###### Abstract

In real world applications, uncertain parameters are the rule rather than the
exception. We present a reachability algorithm for linear systems with
uncertain parameters and inputs using set propagation of polynomial zonotopes.
In contrast to previous methods, our approach is able to tightly capture the
non-convexity of the reachable set. Building up on our main result, we show
how our reachability algorithm can be extended to handle linear time-varying
systems as well as linear systems with time-varying parameters. Moreover, our
approach opens up new possibilities for reachability analysis of linear time-
invariant systems, nonlinear systems, and hybrid systems. We compare our
approach to other state of the art methods, with superior tightness on two
benchmarks including a 9-dimensional vehicle platooning system. Moreover, as
part of the journal extension, we investigate through a polynomial zonotope
with special structure named multi-affine zonotopes and its optimization
problem. We provide the corresponding optimization algorithm and experiment
over the examples obatined from two benchmark systems, showing the efficiency
and scalability comparing to the state of the art method for handling such
type of set representation.

###### keywords:

formal verification , reachability analysis , polynomial zonotopes , linear
systems with uncertain parameters

††journal: Nonlinear Analysis: Hybrid Systems  Figure 1: Different classes of
linear systems explained on the example of the Dubins car model, where
$s_{x}$, $s_{y}$ represent the x- and y-position of the car and $v_{x}$,
$v_{y}$ are the corresponding velocities. The final reachable set is shown in
red and the reachable set for the whole time horizon in gray.

##  1 Introduction

Formal verification is an important aspect of modern cyber-physical systems
since many of them, such as autonomous cars, robots collaborating with humans,
and power system controllers, are safety-critical. Reachability analysis is
one of the most prominent techniques for the formal verification of cyber-
physical systems since it can guarantee robust safety in the presence of
uncertain parameters and inputs. In particular, if the reachable set of the
system does not intersect any unsafe regions, the system is proven to be safe.
The focus of reachability analysis is therefore the computation of a tight
enclosure of the reachable set, because tightness increases the chances to
successfully prove safety and avoid spurious counterexamples. In this paper,
we consider linear systems with uncertain parameters and inputs. The time-
point reachable set for these types of systems is in general non-convex, as it
is visualized for the example of the Dubins car in Fig. 1, where the uncertain
parameter is the turning rate. While previous methods compute convex
enclosures of the reachable set and are consequently limited in their
accuracy, we present a new approach that computes very tight non-convex
enclosures.

There are many ways to perform reachability analysis for linear time-invariant
systems without uncertain parameters (case I. in Fig. 1), which is a well
studied problem. The reachable set is either computed based on simulations
using star sets (Bak and Duggirala, 2017; Duggirala and Viswanathan, 2016), or
based on set-propagation using a variety of set representations such as
support functions (Le Guernic and Girard, 2010; Frehse and et al., 2011),
polytopes (Chutinan and Krogh, 2003), ellipsoids (Kurzhanski and Varaiya,
2000), and zonotopes (Girard, 2005). By using Krylov subspaces (Bak et al.,
2019; Althoff, 2019) or block-decomposition (Bogomolov and et al., 2018), one
can compute the reachable set for very high-dimensional systems with thousands
of states. While algorithm parameters such as time step size often have to be
tuned manually, a recent development in the field is to automate the
corresponding tuning process (Wetzlinger et al., 2020; Frehse et al., 2013).

Fewer approaches exist for linear systems with uncertain parameters (case II.
in Fig. 1). One approach (Ghosh and Duggirala, 2021a) computes approximations
of sets which contain the reachable set with a certain probability. Hereby,
candidate sets are first obtained via reachable set calculation for sampled
parameter values and then verified using statistical verification techniques.
A second method (Lal and Prabhakar, 2015) divides the domain formed by the
parameter space and the considered time horizon into a grid in order to
construct a piecewise affine approximation of the flow function, where the
analytical solution for autonomous linear systems is utilized to obtain the
values for the flow function at the grid points. Since the grid size is
refined until the piecewise affine approximation achieves a certain accuracy,
an enclosure of the reachable set can be obtained by bloating the set defined
by the flow function accordingly. Moreover, for discrete-time systems it is
possible to obtain an enclosure of the reachable set by computing the convex
hull of the reachable sets for all vertices of the parameter set in each time
step (Silvestre, 2022). Other methods use set-propagation techniques to
compute reachable sets (Ghosh and Duggirala, 2019, 2021b; Althoff et al.,
2011a). For discrete-time parametric systems with a special structure, the
reachable set can be represented by a generalized star set with bi-linear
constraints (Ghosh and Duggirala, 2019). The authors extended this approach to
continuous-time systems (Ghosh and Duggirala, 2021b), where they bloat the
reachable set for the linear system corresponding to the average parameter
values. Here, the required bloating is given by an upper bound for the
perturbation of the matrix exponential by the parameter uncertainty, which can
be determined using sensitivity analysis. The method closest to our approach
computes tight enclosures of the reachable set for continuous-time parametric
systems with uncertain inputs using zonotopes (Althoff et al., 2011a). Common
ways to represent the parametric uncertainty are interval matrices (Ghosh and
Duggirala, 2021b; Althoff et al., 2011a), linear matrix equations (Ghosh and
Duggirala, 2021a, 2019), matrix zonotopes (Althoff et al., 2011a; Silvestre,
2022), and matrix polytopes (Lal and Prabhakar, 2015), where linear matrix
equations and matrix zonotopes are different representations of the same set.

All approaches listed above assume that the parameters remain constant over
time. However, there also exist some methods (Serry and Reissig, 2018; Althoff
et al., 2011b) that consider the more general case of time-varying uncertain
parameters (case III. in Fig. 1). The framework of differential inequalities
can be utilized to compute a hyper-rectangular enclosure of the reachable set
efficiently via simulation (Serry and Reissig, 2018). Moreover, ellipsoidal
enclosures can be obtained via optimization with respect to linear matrix
inequality constraints (Zhang and et al., 2020). Another approach (Althoff et
al., 2011b) utilizes the Peano-Baker series, which converges to the transition
matrix for time-varying systems. Once an enclosure of the transition matrix is
obtained, the reachable set can be computed via set-propagation using
zonotopes.

While time-varying linear systems (case IV. in Fig. 1) can be solved with
methods for time-varying parameters by enclosing the time-varying system
matrix with a matrix set in each time step, some approaches specialize on
time-varying systems (Serry and Reissig, 2022; Balandin and Kogan, 2020). In
particular, the first approach (Serry and Reissig, 2022) applies a numerical
approximation of the transition matrix for set-propagation using zonotopes,
where a correction term is added to account for the approximation error.
Moreover, for the special case where the initial set is an ellipsoid, the
reachable set can be represented by an ellipsoid whose shape-defining matrix
is the solution of a linear matrix differential equation (Balandin and Kogan,
2020).

In this paper, we present an approach for determining the reachable sets of
linear systems with uncertain parameters through the use of polynomial
zonotopes. The organization of the paper is as follows: Section 2 outlines the
necessary notations and definitions essential for comprehending the proposed
method. Section 3 examines the fundamental properties of these definitions.
Section 4 describes the algorithm for calculating the reachable set for linear
system with uncertain parameters. Section 5 extends the reachability algorithm
to linear systems with time-varying parameters, linear time-varying systems,
nonlinear systems and explore potential applications in the filed of hybrid
systems. It also provides numerical experiments showing the accuracy and
scalability of our reachability algorithms. Section 6, explores an
optimization algorithm for multi-affine zonotope, which is a special
polynomial zonotope. Leveraging the unique structure of the set, we
experimentally shown its efficiency and scalability comparing to state-of-the-
art method for handling such type of set representation. Finally, Section 7
concludes the paper and suggests directions for future research.

This paper extends the work presented in the conference version Luo et al.
(2023) with significant enhancements, detailed as follows:

  * •

Identification of additional structural properties in polynomial zonotopes
that are beneficial for our algorithm, specifically noticing they are multi-
affine zonotopes.

  * •

Development of improved optimization algorithms that utilized the identified
structural properties for more efficient computation.

  * •

Evaluation of the proposed algorithm with two examples to demonstrate its
performance improvements and practical applicability.

  * •

We provide a detail proof of Lemma 1.

##  2 Preliminaries

###  2.1 Notation

We use lowercase letters to represent vectors $v\in\mathbb{R}^{n}$, uppercase
letters to represent matrices $A\in\mathbb{R}^{m\times n}$, uppercase
calligraphic letter to represent sets $\mathcal{S}\subset\mathbb{R}^{n}$, and
uppercase calligraphic bold letters to represent sets of matrices
$\boldsymbol{\mathcal{A}}\subset\mathbb{R}^{m\times n}$. Given a vector
$v\in\mathbb{R}^{n}$, $v_{(i)}$ represents the $i$-th entry of the vector. For
a matrix $A\in\mathbb{R}^{m\times n}$, $A_{(i,j)}$ represents the element at
$i$-th row and $j$-th column, $A_{(i,~{})}$ denotes the $i$-th row, and
$A_{(~{},j)}$ denotes the $j$-th column. Given two matrices
$A_{1}\in\mathbb{R}^{m\times n_{1}}$ and $A_{2}\in\mathbb{R}^{m\times n_{2}}$,
their horizontal concatenation is represented by
$[A_{1}~{}A_{2}]\in\mathbb{R}^{m\times(n_{1}+n_{2})}$. The identity matrix of
size $n\times n$ is denoted by $I_{n}$, and $\boldsymbol{0}_{m\times n}$ or
$\boldsymbol{1}_{m\times n}$ represents a $m\times n$ matrix consisting of
$0$s or $1$s, respectively. The infinity norm of a matrix set is defined as
$||\boldsymbol{\mathcal{A}}||_{\infty}=\max(||A||_{\infty}~{}|~{}A\in\boldsymbol{\mathcal{A}})$.

###  2.2 Definition

We first present several definitions and operations which are required
throughout the paper. One commonly used set representation for reachability
analysis are zonotopes, since they can represent sets very compactly using
generator vectors:

######  Definition 2.1.

(Zonotope) Given center $c\in\mathbb{R}^{n}$ and a generator matrix
$G\in\mathbb{R}^{n\times h}$, a zonotope is defined as

$\begin{split}\mathcal{Z}=\bigg{\\{}&c+\sum_{i=1}^{h}\alpha_{i}\,G_{(~{},i)}~{}\bigg{|}~{}\alpha_{i}\in[-1,1]\bigg{\\}}.\end{split}$

We use shorthand notation $\mathcal{Z}=\langle c,G\rangle_{Z}$.

Figure 2: Visualization of the polynomial zonotope in Example 1.

While zonotopes are always convex, polynomial zonotopes (Althoff, 2013) can
represent non-convex sets, which allows us to tightly capture the non-convex
shape of the reachable set. In this work, we use the sparse representation of
polynomial zonotopes (Kochdumper and Althoff, 2021)333In contrast to
(Kochdumper and Althoff, 2021, Def. 1), we separate the constant offset $c$
from $G$:

######  Definition 2.2.

(Polynomial Zonotope) Given a constant offset $c\in\mathbb{R}^{n}$, 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}$, a polynomial zonotope
is defined as

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

We use the shorthand notation $\mathcal{PZ}=\langle
c,G,G_{I},E,id\rangle_{PZ}$, where $id\in\mathbb{N}^{p}$ is a list of non-
repeated natural numbers storing a unique identifier for each dependent factor
$\alpha_{k}$.

Let us demonstrate polynomial zonotopes by an example:

######  Example 1.

The polynomial zonotope with sparse representation

$\mathcal{PZ}=\bigg{\langle}\begin{bmatrix}0\\\
0\end{bmatrix},\begin{bmatrix}2&0&1\\\ 1&2&1\end{bmatrix},\begin{bmatrix}1\\\
0.5\end{bmatrix},\begin{bmatrix}1&0&1\\\
0&1&3\end{bmatrix},\begin{bmatrix}1\\\ 2\end{bmatrix}\bigg{\rangle}_{PZ}$

defines the non-convex set:

$\mathcal{PZ}=\bigg{\\{}\begin{bmatrix}0\\\
0\end{bmatrix}+\alpha_{1}\begin{bmatrix}2\\\
1\end{bmatrix}+\alpha_{2}\begin{bmatrix}0\\\
2\end{bmatrix}+\alpha_{1}\alpha_{2}^{3}\begin{bmatrix}1\\\
1\end{bmatrix}+\beta_{1}\begin{bmatrix}1\\\ 0.5\end{bmatrix}\hskip
0.5pt\bigg{|}~{}\alpha_{1},\alpha_{2},\beta_{1}\in[-1,1]\bigg{\\}},$

which is visualized in Fig. 2.

The parameter uncertainty for linear systems can be represented using sets of
matrices. In this work we use the two different matrix set representations
interval matrices and matrix zonotopes.

######  Definition 2.3.

(Interval Matrix) Given a pair of lower and upper matrices
$\underline{A},\overline{A}\in\mathbb{R}^{m\times n}$, an interval matrix is
defined as

$\begin{split}\boldsymbol{\mathcal{A}}=\begin{bmatrix}[\underline{A}_{(1,1)},\overline{A}_{(1,1)}]&\ldots&[\underline{A}_{(1,n)},\overline{A}_{(1,n)}]\\\
\vdots&\ddots&\vdots\\\
[\underline{A}_{(m,1)},\overline{A}_{(m,1)}]&\ldots&[\underline{A}_{(m,n)},\overline{A}_{(m,n)}]\end{bmatrix}.\end{split}$

We use the shorthand notation
$\boldsymbol{\mathcal{A}}=\langle\underline{A},\overline{A}\rangle_{IM}$.

Similar to how zonotopes generalize intervals, matrix zonotopes can be used to
generalize interval matrices.

######  Definition 2.4.

(Matrix Zonotope) Given center matrix $A^{(0)}\in\mathbb{R}^{m\times n}$ and
generator matrices $A^{(1)},\ldots,A^{(w)}\in\mathbb{R}^{m\times n}$, a matrix
zonotope is defined as

$\boldsymbol{\mathcal{A}}=\bigg{\\{}A^{(0)}+\sum_{l=1}^{w}\rho_{l}\,A^{(l)}~{}\bigg{|}\
\rho_{l}\in[-1,1]\bigg{\\}}.$

We use shorthand notation $\boldsymbol{\mathcal{A}}=\langle
A^{(0)},A^{(1)},\ldots,A^{(w)},id\rangle_{MZ}$ with $id\in\mathbb{N}^{w}$
storing the unique identifiers for each factor $\rho_{l}$.

The identifiers allow us to preserve dependencies between sets. To achieve
that, we require the operation mergeID (Kochdumper and Althoff, 2021, Prop. 1)
to bring the exponent matrices and id lists into a common format.

######  Definition 2.5.

Given two exponent matrices $E_{1}\in\mathbb{R}^{p_{1}\times h_{1}},$
$E_{2}\in\mathbb{R}^{p_{2}\times h_{2}}$ and two identifier lists
$id_{1}\in\mathbb{R}^{p_{1}}$, $id_{2}\in\mathbb{R}^{p_{2}}$, mergeID returns
two aligned matrices $\overline{E}_{1},~{}\overline{E}_{2}$ and a list of
adjusted identifiers $\overline{id}$:

$\overline{E}_{1},~{}\overline{E}_{2},~{}\overline{id}\leftarrow{\textsc{mergeID}}(id_{1},id_{2},E_{1},E_{2})$

with

$\begin{split}&\mathcal{K}=\big{\\{}i~{}|~{}id_{2(i)}\not\in
id_{1}\big{\\}}:=\\{i_{1},~{}\ldots,~{}i_{k}\\},\\\
&\overline{id}=\big{[}id_{1}^{T}\hskip
5.0ptid_{2(i_{1})}\ldots~{}id_{2(i_{k})}\big{]}^{T},\hskip
5.0pt\overline{E}_{1}=\begin{bmatrix}E_{1}\\\ \boldsymbol{0}_{k\times
h_{1}}\end{bmatrix}\in\mathbb{R}^{(p_{1}+k)\times h_{1}},\\\
&\overline{E}_{2(i,~{})}=\begin{cases}E_{2(j,~{})},&\text{if}\hskip
5.0pt\exists j~{}\,\overline{id}_{(i)}=id_{2(j)}\\\ \boldsymbol{0}_{1\times
h_{2}},&\text{otherwise}\end{cases},~{}i=1,\ldots,p_{1}+k.\end{split}$

Sometimes we also need to destroy dependencies to guarantee the correctness of
our reachability algorithm. For this we introduce the operation
$\textsc{uniqueID}(p)$, which returns a vector of unique and non-existing
identifiers with length $p$. Moreover, we define a convenience operation
$\textsc{fresh}(\mathcal{PZ})$ $=\langle c,G,G_{I},E,$
$\textsc{uniqueID}(p)\rangle_{PZ}$ replacing all identifiers of
$\mathcal{PZ}$. The fresh operator can equally be applied to matrix zonotopes.
For the ease of notation we denote by
$\textsc{fresh}(\mathcal{PZ},\boldsymbol{\mathcal{A}})$ that all dependent
factors of $\mathcal{PZ}$ that are not part of the matrix zonotope
$\boldsymbol{\mathcal{A}}$ are equipped with new unique identifiers. The
operation $\textsc{eval}(\mathcal{PZ},id,val)$ replaces each dependent factor
of $\mathcal{PZ}$ that has the identifier id with a constant value
$val\in[-1,1]$. The operation also applies to matrix zonotopes. As an example,
performing $\textsc{eval}(\mathcal{PZ},2,0.5)$ on the polynomial zonotope from
Example 1 yields the set

$\textsc{eval}(\mathcal{PZ},2,0.5)=\\\ \bigg{\\{}\alpha_{1}\begin{bmatrix}2\\\
1\end{bmatrix}+0.5\begin{bmatrix}0\\\
2\end{bmatrix}+\alpha_{1}(0.5)^{3}\begin{bmatrix}1\\\
1\end{bmatrix}+\beta_{1}\begin{bmatrix}1\\\
0.5\end{bmatrix}~{}\bigg{|}~{}\alpha_{1},\beta_{1}\in[-1,1]\bigg{\\}}.\\\ $

In this paper we use set operations to perform linear maps with either a
matrix or a matrix set, as well as sums that can preserve or drop
dependencies. Given sets
$\mathcal{S}_{1},\mathcal{S}_{2}\subset\mathbb{R}^{n}$, a numerical matrix
$A\in\mathbb{R}^{m\times n}$ and a matrix set
$\boldsymbol{\mathcal{A}}\subset\mathbb{R}^{m\times n}$, these operations are:

$\displaystyle A\,\mathcal{S}_{1}$
$\displaystyle=\\{A\,s~{}|~{}s\in\mathcal{S}_{1}\\},$ (1)
$\displaystyle\boldsymbol{\mathcal{A}}~{}\mathcal{S}_{1}$
$\displaystyle=\\{A\,s~{}|~{}A\in\boldsymbol{\mathcal{A}},s\in\mathcal{S}_{1}\\},$
(2) $\displaystyle\mathcal{S}_{1}\oplus\mathcal{S}_{2}$
$\displaystyle=\\{s_{1}+s_{2}~{}|~{}s_{1}\in\mathcal{S}_{1},s_{2}\in\mathcal{S}_{2}\\}.$
(3)

Given two polynomial zonotopes $\mathcal{PZ}_{1}=\langle
c_{1},G_{1},G_{I1},E_{1},id_{1}\rangle_{PZ}$, $\mathcal{PZ}_{2}=\langle
c_{2},G_{2},G_{I2},E_{2},id_{2}\rangle_{PZ}$ $\subset\mathbb{R}^{n}$ and a
zonotope $\mathcal{Z}=\langle c_{z},G_{z}\rangle_{Z}$ $\subset\mathbb{R}^{n}$,
the Minkowski sum can be computed as (Kochdumper and Althoff, 2021, Prop. 9):

$\displaystyle\mathcal{PZ}_{1}\oplus\mathcal{PZ}_{2}$ $\displaystyle=$
$\displaystyle\hskip
1.0pt\bigg{\langle}c_{1}+c_{2},[G_{1}~{}G_{2}],[G_{I1}~{}G_{I2}],\;\begin{bmatrix}E_{1}&\mathbf{0}\\\
\mathbf{0}&E_{2}\end{bmatrix},\textsc{uniqueID}(id_{1},id_{2})~{}\bigg{\rangle}_{PZ},$
$\displaystyle\mathcal{PZ}_{1}\oplus\mathcal{Z}$ $\displaystyle=$
$\displaystyle\langle
c_{1}+c_{z},G_{1},[G_{I1}~{}G_{z}],E_{1},id_{1}\rangle_{PZ}.$

While Minkowski sum drops dependencies between sets, the exact sum operation
$\boxplus$ explicitly keeps these dependencies (Kochdumper and Althoff, 2021,
Prop. 10):

$\begin{split}&\mathcal{PZ}_{1}~{}\boxplus~{}\mathcal{PZ}_{2}=\big{\langle}c_{1}+c_{2},[G_{1}~{}G_{2}],[G_{I1}~{}G_{I2}],[\overline{E}_{1}~{}\overline{E}_{2}],\overline{id}\big{\rangle}_{PZ},\end{split}$
(4)

Figure 3: Comparison from Example 2 between Minkowski sum
$A~{}\mathcal{PZ}\oplus\mathcal{PZ}$, exact sum
$A~{}\mathcal{PZ}\boxplus\mathcal{PZ}$, and a mixture of both obtained by
destroying the dependencies for the second dependent factor.

where the exponent matrices and identifier are first merged using
$\overline{E}_{1},~{}\overline{E}_{2},~{}\overline{id}\leftarrow\textsc{mergeID}(id_{1},id_{2},E_{1},E_{2})$.
The operation compact (Kochdumper and Althoff, 2021, Prop. 2) is applied after
each exact sum to remove redundancies such as equal columns in the exponent
matrix. Let us demonstrate the difference between Minkowski sum and exact sum:

######  Example 2.

We consider the polynomial zonotope $\mathcal{PZ}=\langle\mathbf{0}_{2\times
1},I_{n},\linebreak[3][~{}],I_{n},[1\hskip 5.0pt2]^{T}\rangle_{PZ}$ together
with the numerical matrix $A=[[1\hskip 5.0pt1]^{T}\hskip
5.0pt\linebreak[3][-1\hskip 5.0pt1]^{T}]$. The comparison of the Minkowski sum
$A~{}\mathcal{PZ}\oplus\mathcal{PZ}$ with the exact sum
$A~{}\mathcal{PZ}\boxplus\mathcal{PZ}$ visualized in Fig. 3 demonstrates that
the Minkowski sum yields an over-approximation of the exact result since the
dependencies between the sets $A~{}\mathcal{PZ}$ and $\mathcal{PZ}$ are
destroyed. Moreover, if we only destroy the dependencies for some factors, as
for example done by the operation
${\textsc{fresh}}(\mathcal{PZ},\boldsymbol{\mathcal{A}})$, we obtain a mixture
between Minkowski sum and exact sum.

We additionally require several other operations that are provided in a
previous paper on polynomial zonotopes Kochdumper and Althoff (2021): The
operator $\textsc{zonotope}(\mathcal{PZ})$ (Kochdumper and Althoff, 2021,
Prop. 5) returns a tight enclosing zonotope of a polynomial zonotope.
Moreover, $~{}\textsc{reduce}(\mathcal{PZ},\rho)$ (Kochdumper and Althoff,
2021, Prop. 16) reduces the representation size and returns a polynomial
zonotope with order smaller or equal to $\rho$ which encloses $\mathcal{PZ}$.

We now introduce a special kind of polynomial zonotope

######  Definition 2.6.

A _multi-affine zonotope_ is a polynomial zonotope in which the exponent
matrix $E$ exclusively composed of $0$ or $1$ elements. The multi-affine
zonotope has similar shorthand notation $\mathcal{MAZ}=\langle
c,G,G_{I},E,id\rangle_{MAZ}$.

##  3 SET OPERATIONS

For our proposed reachability algorithm, we require additional operations on
polynomial zonotopes, which we derive in this section. We begin with the
multiplication of a matrix zonotope with a polynomial zonotope, for which we
consider the case without independent generators first and present the general
case afterwards.

######  Proposition 1.

(Matrix Zonotope Multiplication) Given a matrix zonotope
$\boldsymbol{\mathcal{A}}=\langle
A^{(0)},A^{(1)},\ldots,A^{(w)},id_{mz}\rangle_{MZ}\subset\mathbb{R}^{m\times
n}$ and a polynomial zonotope $\mathcal{PZ}=\langle
c,G,[~{}],E,id_{pz}\rangle_{PZ}\subset\mathbb{R}^{n}$, their multiplication is

$\begin{split}\boldsymbol{\mathcal{A}}~{}\mathcal{PZ}=\hskip
1.0pt&\big{\langle}A^{(0)}c,[A^{(0)}G\hskip
5.0ptA^{(1)}c~{}\ldots~{}A^{(w)}c\hskip 5.0ptA^{(1)}G~{}\ldots~{}A^{(w)}G],\\\
&\hskip 95.0pt[~{}],[\widehat{E}_{1}\hskip 5.0pt\widehat{E}_{2}\hskip
5.0pt\overline{E}_{1}~{}\ldots~{}\overline{E}_{w}],\widehat{id}\big{\rangle}_{PZ}\end{split}$

with

$\overline{E}_{l}=\widehat{E}_{2(~{},l)}\cdot\boldsymbol{1}_{1\times
h}+\widehat{E}_{1},\quad l=1,\ldots,w,$ (5)

where
$\widehat{E}_{1},\widehat{E}_{2},\widehat{id}\leftarrow\textsc{mergeID}(id_{pz},id_{mz},E,I_{w})$
brings the exponent matrices into a common format.

###### Proof.

Inserting the definitions of matrix zonotopes and polynomial zonotopes in Def.
2.2 and Def. 2.4 into the definition of the multiplication with a matrix set
in (2) yields:

$\displaystyle\boldsymbol{\mathcal{A}}~{}\mathcal{PZ}\overset{\eqref{def:matrix
set multiply a set of states}}{=}\big{\\{}As\ \big{|}\
A\in\boldsymbol{\mathcal{A}},s\in\mathcal{PZ}\big{\\}}$
$\displaystyle=\bigg{\\{}\bigg{(}A^{(0)}+\sum_{l=1}^{w}\rho_{l}A^{(l)}\bigg{)}\bigg{(}c+\sum_{i=1}^{h}\bigg{(}\prod_{k=1}^{m}\alpha_{k}^{E_{(k,i)}}\bigg{)}G_{(~{},i)}\bigg{)}~{}\bigg{|}~{}\rho_{l},\alpha_{k}\in[-1,1]\bigg{\\}}$
$\displaystyle=\bigg{\\{}A^{(0)}c+\sum_{i=1}^{h}\bigg{(}\prod_{k=1}^{m}\alpha_{k}^{E_{(k,i)}}\bigg{)}A^{(0)}G_{(~{},i)}+\sum_{l=1}^{w}\rho_{l}A^{(l)}c\hskip
1.0pt+$ $\displaystyle\hskip
20.0pt\sum_{l=1}^{w}\sum_{i=1}^{h}\bigg{(}\prod_{k=1}^{m}\alpha_{k}^{E_{(k,i)}}\bigg{)}\rho_{l}A^{(l)}G_{(~{},i)}~{}\bigg{|}~{}\rho_{l},\alpha_{k},\in[-1,1]\bigg{\\}}$
$\displaystyle=\big{\langle}A^{(0)}c,[A^{(0)}G\hskip
5.0ptA^{(1)}c~{}\ldots~{}A^{(w)}c\hskip
5.0ptA^{(1)}G~{}\ldots~{}A^{(w)}G],[~{}],$ $\displaystyle\hskip
135.0pt[\widehat{E}_{1}\hskip 5.0pt\widehat{E}_{2}\hskip
5.0pt\overline{E}_{1}~{}\ldots~{}\overline{E}_{w}],\widehat{id}\big{\rangle}_{PZ}.$

Here, we exploit that any matrix zonotope has an exponent matrix $I_{w}$
because each generator matrix is only multiplied with a single factor. Since
the operation mergeID brings the exponent matrices to a common representation,
the multiplication of the factors
$\big{(}\prod_{k=1}^{m}\alpha_{k}^{\begin{subarray}{c}E_{(k,i)}\\\
\vspace{-5pt}\end{subarray}}\big{)}~{}\rho_{l}$ can be considered by an
addition of the $l$-th column of exponent matrix $\widehat{E}_{2}$ to $i$-th
column of $\widehat{E}_{1}$ as done in $\eqref{cal:mtimes new exponent
matrices}$. ∎

Based on Prop. 1, the multiplication with powers of matrix sets can be
computed as follows:

######  Corollary 1.

(Higher Order Multiplication) Given a matrix zonotope
$\boldsymbol{\mathcal{A}}$ and a polynomial zonotope $\mathcal{PZ}$, the set

$\big{\\{}A^{k}s~{}\big{|}~{}A\in\boldsymbol{\mathcal{A}},s\in\mathcal{PZ}\big{\\}}=\boldsymbol{\mathcal{A}}^{k}~{}\mathcal{PZ}$

can be computed using $k$ multiplications
$\boldsymbol{\mathcal{A}}\,\big{(}~{}\ldots~{}\big{(}\boldsymbol{\mathcal{A}}~{}(\boldsymbol{\mathcal{A}}~{}\mathcal{PZ})\big{)}\big{)}$.

Any polynomial zonotope can equivalently be represented as a polynomial
zonotope without independent generators by introducing new identifiers for the
independent factors and expanding the exponent matrix (Kochdumper et al.,
2020, Prop. 1). Therefore, also for the case with independent generators the
multiplication with a matrix zonotope can always be performed exactly. In
practice, for computational efficiency and based on the assumption that the
independent generators are relatively small, we compute a tight enclosure of
the multiplication with the independent generators:

######  Proposition 2.

(Matrix Zonotope Multiplication) Given a matrix zonotope
$\boldsymbol{\mathcal{A}}=\langle
A^{(0)},A^{(1)},\ldots,A^{(w)},id_{mz}\rangle_{MZ}\subset\mathbb{R}^{m\times
n}$ and a polynomial zonotope $\mathcal{PZ}$ = $\langle
c,G,G_{I},E,id_{pz}\rangle_{PZ}\subset\mathbb{R}^{n}$, their multiplication
can be tightly enclosed by

$\begin{split}\boldsymbol{\mathcal{A}}~{}\mathcal{PZ}\subseteq\big{\langle}\widehat{c},\widehat{G},[A^{(0)}G_{I}~{}\dots~{}A^{(w)}G_{I}],\widehat{E},\widehat{id}\big{\rangle}_{PZ},\end{split}$

where

$\big{\langle}\widehat{c},\widehat{G},[~{}],\widehat{E},\widehat{id}\big{\rangle}_{PZ}=\boldsymbol{\mathcal{A}}~{}\big{\langle}c,G,[~{}],E,id_{pz}\big{\rangle}_{PZ}$

is computed using Prop. 1.

###### Proof.

The independent generators of $\mathcal{PZ}$ can be viewed as a zero-centered
zonotope, so $\mathcal{PZ}=\mathcal{PZ}_{D}\oplus\mathcal{Z}_{I}$ is given by
the Minkowski sum between the dependent part $\mathcal{PZ}_{D}$ = $\langle
c,G,[~{}],E,\\\ id_{pz}\rangle_{PZ}$ and the independent part
$\mathcal{Z}_{I}=\langle\boldsymbol{0}_{n\times 1},G_{I}\rangle_{Z}$.
Therefore, the multiplication can be over-approximated by:

$\boldsymbol{\mathcal{A}}~{}\mathcal{PZ}\subseteq(\boldsymbol{\mathcal{A}}~{}\mathcal{PZ}_{D})\oplus(\boldsymbol{\mathcal{A}}~{}\mathcal{Z}_{I}),$

where the result for the multiplication between a matrix zonotope with a
zonotope is given by (Althoff et al., 2011b, Prop. 4). ∎

The core operation for reachability analysis is the linear map with the matrix
exponential. In particular, we require the multiplication with the matrix
exponential of two matrix sets $\boldsymbol{\mathcal{A}}$ and
$\boldsymbol{\mathcal{B}}$:

######  Proposition 3.

(Multiplication with Matrix Exponential) Given two matrix zonotopes
$\boldsymbol{\mathcal{A}}$,
$\boldsymbol{\mathcal{B}}\subset\mathbb{R}^{n\times n}$, a polynomial zonotope
$\mathcal{PZ}\subset\mathbb{R}^{n}$, and a Taylor order $\kappa$ satisfying
$\epsilon=\frac{||\boldsymbol{\mathcal{A}}||_{\infty}||\boldsymbol{\mathcal{B}}||_{\infty}}{\kappa+2}<1$,
the set $e^{\boldsymbol{\mathcal{A}}\boldsymbol{\mathcal{B}}}\mathcal{PZ}$ can
be tightly enclosed by

$\begin{split}e^{\boldsymbol{\mathcal{A}}\boldsymbol{\mathcal{B}}}~{}\mathcal{PZ}\subseteq\bigg{(}\operatorname*{\scalerel*{\boxplus}{\sum}}_{i=0}^{\kappa}\frac{\boldsymbol{\mathcal{A}}^{i}\boldsymbol{\mathcal{B}}^{i}}{i!}~{}\mathcal{PZ}\bigg{)}\oplus\Big{(}\boldsymbol{\mathcal{E}}~{}{\textsc{zonotope}}(\mathcal{PZ})\Big{)},\end{split}$

where the interval matrix

$\begin{split}&\boldsymbol{\mathcal{E}}=\langle-\boldsymbol{1}_{n\times
n},\boldsymbol{1}_{n\times
n}\rangle_{IM}\frac{(||\boldsymbol{\mathcal{A}}||_{\infty}||\boldsymbol{\mathcal{B}}||_{\infty})^{\kappa+1}}{(\kappa+1)!}\frac{1}{1-\epsilon}\end{split}$

encloses the remaining terms of the truncated Taylor series. The sets
$\frac{\boldsymbol{\mathcal{A}}^{i}\boldsymbol{\mathcal{B}}^{i}}{i!}~{}\mathcal{PZ}$
are calculated according to Corollary 1, and the multiplication of an interval
matrix with a zonotope using (Althoff et al., 2007, Thm. 4).

###### Proof.

The matrix exponential is defined by its Taylor series:

$\begin{split}&e^{\boldsymbol{\mathcal{A}}\boldsymbol{\mathcal{B}}}~{}\mathcal{PZ}=\bigg{\\{}\bigg{(}\sum_{i=0}^{\infty}\frac{A^{i}B^{i}}{i!}\bigg{)}\,s~{}\bigg{|}~{}A\in\boldsymbol{\mathcal{A}},B\in\boldsymbol{\mathcal{B}},s\in\mathcal{PZ}\bigg{\\}}\overset{\text{\cite[cite]{(\@@bibref{AuthorsPhrase1Year}{Althoff2010a}{\@@citephrase{,
}}{}, Thm.~{}3.2)}}}{=}\\\\[5.0pt]
&\bigg{\\{}\bigg{(}\sum_{i=0}^{\kappa}\frac{A^{i}B^{i}}{i!}\bigg{)}\,s+\underbrace{\bigg{(}\sum_{i=\kappa+1}^{\infty}\frac{A^{i}B^{i}}{i!}\bigg{)}}_{\subseteq\boldsymbol{\mathcal{E}}}s~{}\bigg{|}~{}A\in\boldsymbol{\mathcal{A}},B\in\boldsymbol{\mathcal{B}},s\in\mathcal{PZ}\bigg{\\}}\subseteq\\\
&\bigg{\\{}\bigg{(}\sum_{i=0}^{\kappa}\frac{A^{i}B^{i}}{i!}\bigg{)}\,s~{}\bigg{|}~{}A\in\boldsymbol{\mathcal{A}},B\in\boldsymbol{\mathcal{B}},s\in\mathcal{PZ}\bigg{\\}}\oplus\big{\\{}E\,s~{}\big{|}~{}E\in\boldsymbol{\mathcal{E}},~{}s\in\mathcal{PZ}\big{\\}}\\\\[5.0pt]
&\subseteq\bigg{(}\operatorname*{\scalerel*{\boxplus}{\sum}}_{i=0}^{\kappa}\frac{\boldsymbol{\mathcal{A}}^{i}\boldsymbol{\mathcal{B}}^{i}}{i!}\mathcal{PZ}\bigg{)}\oplus\big{(}\boldsymbol{\mathcal{E}}~{}\textsc{zonotope}(\mathcal{PZ})\big{)}.\end{split}$

Since the remainder interval matrix $\boldsymbol{\mathcal{E}}$ is expected to
be small for a sufficiently large Taylor order, we add the corresponding
vectors to the independent generators by applying a zonotope enclosure. ∎

##  4 Reachability Analysis

We now explain how to compute reachable sets for linear systems with uncertain
parameters (case II. in Fig. 1) using the operations on polynomial zonotopes
derived in the previous section.

###  4.1 Problem Definition

We consider a linear system with uncertain parameters

$\dot{x}(t)=A\,x(t)+B\,u(t),\quad
A\in\boldsymbol{\mathcal{A}},B\in\boldsymbol{\mathcal{B}}~{},$ (6)

where $\boldsymbol{\mathcal{A}}\subset\mathbb{R}^{n\times n}$,
$\boldsymbol{\mathcal{B}}\subset\mathbb{R}^{n\times m}$ are matrix zonotopes.
The initial state $x(0)$ is uncertain within
$\mathcal{X}_{0}\subset\mathbb{R}^{n}$, and the input $u(t)$ is uncertain
within $\mathcal{U}\subset\mathbb{R}^{m}$. The reachable set is defined as
follows:

######  Definition 4.1.

(Reachable Set) Given an initial set $\mathcal{X}_{0}$, an input set
$\mathcal{U}$, and matrix zonotopes $\boldsymbol{\mathcal{A}}$ and
$\boldsymbol{\mathcal{B}}$, the reachable set for the system in (6) at time
$t\geq 0$ is

$\mathcal{R}(t)=\big{\\{}\xi(A,B,t,x(0),u(\cdot))~{}\big{|}~{}A\in\boldsymbol{\mathcal{A}},~{}B\in\boldsymbol{\mathcal{B}},~{}x(0)\in\mathcal{X}_{0},\forall
s\in[0,t]:~{}u(\hskip 1.0pts\hskip 1.0pt)\in\mathcal{U}\big{\\}},$

where $\xi(A,B,t,x(0),u(\cdot))$ denotes the solution of (6) at time t for
numerical matrices $A$ and $B$, initial state $x(0)$, and input signal
$u(\cdot)$.

The goal is to calculate the reachable set for a finite time interval
$[0,t_{end}]$, where time starts at $0$ without loss of generality. Since the
exact reachable set cannot be computed in general, we instead aim to calculate
a tight enclosure of the reachable set. As commonly done in reachability
analysis, we compute the reachable set for consecutive time intervals
$\tau_{k}=[k\Delta t,(k+1)\Delta t]$ with time step $\Delta t$, so that the
reachable set for the whole time horizon is given as
$\mathcal{R}([0,t_{end}])=\bigcup_{k=0}^{(t_{end}/\Delta
t)-1}\mathcal{R}(\tau_{k})$. According to the superposition principle, the
analytical solution to $\eqref{def:probdef}$ for numerical matrices $A$ and
$B$ is given by the addition of the homogeneous solution for the initial state
and the particular solution due to inputs:

$\xi(A,B,t,x(0),u(\cdot))=e^{At}x(0)+\int_{0}^{t}e^{A(t-s)}B\,u(s)\,ds.$ (7)

This concept also carries over to the case with uncertain parameters, so that
we can first compute the homogeneous solution $\mathcal{H}(\tau_{k})$ and the
particular solution $\mathcal{P}(\tau_{k})$ separately, and later combine them
to obtain the reachable set $\mathcal{R}(\tau_{k})$.

###  4.2 Homogeneous Solution

First, we consider the homogeneous solution $\mathcal{H}(\tau_{0})$ in the
first time interval $\tau_{0}=[0,\Delta t]$, which can be calculated as

$\begin{split}\mathcal{H}(\tau_{0})=\big{\\{}e^{At}x(0)~{}\big{|}~{}A\in\boldsymbol{\mathcal{A}},t\in\tau_{0},x(0)\in\mathcal{X}_{0}\big{\\}}=e^{\boldsymbol{\mathcal{A}}\boldsymbol{\mathcal{T}}}\,\mathcal{X}_{0},\end{split}$

where we represent time interval $\tau_{0}=[0,\Delta t]$ by a matrix zonotope

$\begin{split}&\boldsymbol{\mathcal{T}}=[0,\Delta t]\cdot
I_{n}=\big{\langle}T,T,\textsc{uniqueID}(1)\big{\rangle}_{MZ},\quad
T=0.5\,\Delta t\,I_{n}.\end{split}$ (8)

We can compute a tight enclosure of the set
$e^{\boldsymbol{\mathcal{A}}\boldsymbol{\mathcal{T}}}\,\mathcal{X}_{0}$ using
Prop. 3. Given the homogeneous solution for the first time interval, the
solutions for all remaining time intervals can be computed by propagating the
initial solution forward in time via multiplication with the matrix
exponential $e^{\boldsymbol{\mathcal{A}}\Delta t}:$

$\mathcal{H}(\tau_{k})=\underbrace{e^{\boldsymbol{\mathcal{A}}\Delta
t}~{}\ldots~{}e^{\boldsymbol{\mathcal{A}}\Delta
t}}_{k}~{}e^{\boldsymbol{\mathcal{A}}\boldsymbol{\mathcal{T}}}~{}\mathcal{X}_{0}.$

Note that even if $\mathcal{X}_{0}$ is convex, $\mathcal{H}(\tau_{k})$ is
generally non-convex.

###  4.3 Particular Solution

According to (Althoff, 2010, Appendix A), the particular solution at time $t$
for numerical matrices $A\in\boldsymbol{\mathcal{A}}$ and
$B\in\boldsymbol{\mathcal{B}}$ can be enclosed by

$\mathcal{P}_{A,B}(t)=\bigg{\\{}\int_{0}^{t}e^{A(t-s)}~{}B\,u(t)\,ds~{}\bigg{|}~{}u(t)\in\mathcal{U}\bigg{\\}}\subseteq\bigoplus_{i=0}^{\kappa}\bigg{(}\frac{A^{i}t^{i+1}}{(i+1)!}~{}B\,\mathcal{U}\bigg{)}\oplus\big{(}t\,\boldsymbol{\mathcal{E}}_{t}\,\textsc{zonotope}(B\,\mathcal{U})\big{)},$
(9)

where

$\boldsymbol{\mathcal{E}}_{t}=\langle-\boldsymbol{1}_{n\times
n},\boldsymbol{1}_{n\times
n}\rangle_{IM}\frac{(||\boldsymbol{\mathcal{A}}||_{\infty}t)^{\kappa+1}}{(\kappa+1)!}\frac{1}{1-\epsilon},~{}\epsilon=\frac{||A||_{\infty}t}{\kappa+2}<1.$
(10)

Evaluating $\eqref{def: input solution set}$ for matrix sets
$\boldsymbol{\mathcal{A}}$ and $\boldsymbol{\mathcal{B}}$ yields

$\mathcal{P}(t)\subseteq\bigg{(}\operatorname*{\scalerel*{\boxplus}{\sum}}_{i=0}^{\kappa}\frac{\boldsymbol{\mathcal{A}}^{i}t^{i+1}}{(i+1)!}~{}\boldsymbol{\mathcal{B}}\hskip
1.0pt\textsc{fresh}(\mathcal{U})\bigg{)}\oplus\big{(}t\,\boldsymbol{\mathcal{E}}_{t}\,\textsc{zonotope}(\boldsymbol{\mathcal{B}}\,\mathcal{U})\big{)},$
(11)

where we use exact sum to keep the dependencies between the matrix zonotopes
$\boldsymbol{\mathcal{A}}$ and $\boldsymbol{\mathcal{B}}$. At the same time,
we apply the operator $\textsc{fresh}(\mathcal{U})$ to destroy the
dependencies between the input sets $\mathcal{U}$, which ensures that we
maintain the over-approximate nature of $\eqref{def: input solution set}$.
Moreover, based on $\eqref{def: time-point particular solution}$, the
particular solution for the initial time interval is given by

$\mathcal{P}(\tau_{0})\subseteq\bigg{(}\operatorname*{\scalerel*{\boxplus}{\sum}}_{i=0}^{\kappa}\frac{\boldsymbol{\mathcal{A}}^{i}\boldsymbol{\mathcal{T}}^{i+1}}{(i+1)!}~{}\boldsymbol{\mathcal{B}}\hskip
1.0pt\textsc{fresh}(\mathcal{U})\bigg{)}\oplus\big{(}\Delta
t\,\boldsymbol{\mathcal{E}}_{\Delta
t}\,\textsc{zonotope}(\boldsymbol{\mathcal{B}}\,\mathcal{U})\big{)},$ (12)

where the matrix zonotope $\boldsymbol{\mathcal{T}}$ from $\eqref{def: time
matrix zonotope}$ equivalently represents the time interval
$\tau_{0}=[0,\Delta t].$ Given the particular solution for the first time
interval, the particular solutions for the remaining time intervals can
according to (Althoff, 2010, Prop. 3.2) be computed by the following
propagation scheme:

$\mathcal{P}_{A,B}(\tau_{k})=e^{A\Delta
t}\mathcal{P}_{A,B}(\tau_{k-1})\oplus\mathcal{P}_{A,B}(\Delta t),$ (13)

which considers numerical matrices $A$ and $B$. Evaluating $\eqref{def:
particular solution recurrence relation}$ for matrix sets
$\boldsymbol{\mathcal{A}}$ and $\boldsymbol{\mathcal{B}}$ yields

$\begin{split}\mathcal{P}(\tau_{k})\subseteq~{}&e^{\boldsymbol{\mathcal{A}}\Delta
t}~{}\mathcal{P}(\tau_{k-1})\boxplus\textsc{fresh}\big{(}\mathcal{P}(\Delta
t),\boldsymbol{\mathcal{A}},\boldsymbol{\mathcal{B}}\big{)},\end{split}$

where we again use the exact sum to preserve the dependencies between the
matrix sets $\boldsymbol{\mathcal{A}}$ and $\boldsymbol{\mathcal{B}}$, but
destroy dependencies between the particular solutions at different points in
time using fresh to realize the Minkowski sum in $\eqref{def: particular
solution recurrence relation}$.

###  4.4 Reachability Algorithm

Combining the homogeneous solution from Sec. 4.2 and the particular solution
from Sec. 4.3, we obtain the reachability algorithm summarized in Alg. 1:

Algorithm 1 Reachability algorithm (constant parameters) Figure 4: Comparison
between combining the homogeneous solution $\mathcal{H}(\tau_{0})$ and
particular solution $\mathcal{P}(\tau_{0})$ via Minkowski sum (left) and via
exact sum (right).

The algorithm takes advantage of dependency preservation for polynomial
zonotopes and is therefore able to compute very tight non-convex enclosures of
the reachable set. In Line LABEL:algo1:_initizal_homo and
LABEL:algo1:_initizal_input, the homogeneous solution and particular solution
for the first time interval $\tau_{0}=[0,\Delta t]$ are calculated. Since both
of them preserve dependency on time, we can combine them using the exact sum
in Line LABEL:alg1:homo_add_input. As visualized in Fig. 4, this results in a
much tighter set compared to adding the particular solution
$\mathcal{P}(\Delta t)$ at the end of the time interval to every part of the
homogeneous solution, which is the technique applied by other reachability
algorithms Althoff et al. (2011a); Frehse and et al. (2011). Moreover, the
property of time preservation not only holds for the initial step, but also
for all consecutive steps since forward-propagation via multiplication with
$e^{\boldsymbol{\mathcal{A}}\Delta t}$ does not destroy dependency on time.
Therefore, we can also combine the propagated homogeneous and particular
solution for all remaining time intervals with the exact sum in Line
LABEL:alg1:_homo_add_input_in_rounds. Moreover, since dependencies on time are
preserved, we can extract the reachable set $\mathcal{R}(t)$ for a specific
time point $t\in\tau_{k}$ from the reachable set $\mathcal{R}(\tau_{k})$ of
the corresponding time interval $\tau_{k}$ using the eval operation:

$\forall
t\in\tau_{k}:\quad\mathcal{R}(t)=\textsc{eval}\big{(}\mathcal{R}(\tau_{k}),id_{t},\alpha_{t}\big{)},\quad\alpha_{t}=2\,\frac{t-k\,\Delta
t}{\Delta t}-1,$

where $id_{t}$ is the unique id of the dependent factor that represents time
and $\alpha_{t}\in[-1,1]$ maps the time $t$ to the corresponding factor.

Finally, to reduce the over-approximation introduced by the remainder
$\boldsymbol{\mathcal{E}}$ in Prop. 3 one can evaluate additional Taylor terms
with interval arithmetic Jaulin et al. (2006) before truncating the Taylor
series. Moreover, one can apply restructuring as described in (Kochdumper and
Althoff, 2021, Sec. III.A) to move large independent generators to the
dependent generators, which improves accuracy.

##  5 Applications and Reachability Evaluation

There exist several extensions and applications for the reachability algorithm
from the previous section, which we present next. In particular, we first show
how to extend the algorithm to linear systems with time-varying parameters
(Sec. 5.1) and linear time-varying systems (Sec. 5.2). Next, we demonstrate
the advantages our approach has for linear time-invariant systems without
uncertain parameters (Sec. 5.3), before we explain how it can be applied to
hybrid systems (Sec. 5.4) and nonlinear systems (Sec. 5.5).

###  5.1 Time-Varying Parameters

While our basic reachability algorithm presented in Sec. 4 considers systems
where the parameters are constant over time, it can be easily extended to also
handle time-varying parameters. For this, we simply have to replace the state
transition matrix $e^{\boldsymbol{\mathcal{A}}t}$ for the case with constant
parameters by an enclosure of state transition matrix
$\boldsymbol{\mathcal{M}}(t)$ for the case with time varying parameters.
According to (Althoff et al., 2011b, Thm. 1), such an enclosure can be
computed as:

$\begin{split}&\boldsymbol{\mathcal{M}}(t)\subseteq\bigoplus_{i=0}^{\kappa}\overline{\boldsymbol{\mathcal{M}}}_{i}(t)~{}\oplus~{}\boldsymbol{\mathcal{E}}_{t},\quad\overline{\boldsymbol{\mathcal{M}}}_{i}(t)=\frac{t^{i}}{i!}\textsc{conv}\big{(}\textsc{fresh}^{i}(\boldsymbol{\mathcal{A}})\big{)},\end{split}$
(14)

where $\kappa\geq 1$ is the user-defined Taylor order, the interval matrix
$\boldsymbol{\mathcal{E}}_{t}$ is defined as in (10), and operation
$\textsc{conv}(\boldsymbol{\mathcal{A}})$ represents the convex hull of a
matrix set $\boldsymbol{\mathcal{A}}$ which is defined as

$\textsc{conv}(\boldsymbol{\mathcal{A}})=\big{\\{}\lambda\,A_{1}+(1-\lambda)\,A_{2}~{}\big{|}~{}A_{1},A_{2}\in\boldsymbol{\mathcal{A}},~{}\lambda\in[0,1]\big{\\}}.$
(15)

Moreover, the notation

$\begin{split}\textsc{fresh}^{i}(\boldsymbol{\mathcal{A}})=\underbrace{\textsc{fresh}(\boldsymbol{\mathcal{A}})~{}\ldots~{}\textsc{fresh}(\boldsymbol{\mathcal{A}})}_{i}\end{split}$

emphasizes that each matrix zonotope is equipped with different unique
identifiers before the multiplication444Note that the notation
$\textsc{fresh}^{i}(\boldsymbol{\mathcal{A}})$ is different from
$(\textsc{fresh}(\boldsymbol{\mathcal{A}}))^{i}$, where the latter one
represents $\boldsymbol{\mathcal{B}}^{i}$ with
$\boldsymbol{\mathcal{B}}=\textsc{fresh}(\boldsymbol{\mathcal{A}})$. . For our
reachability algorithm, we require the following proposition:

######  Proposition 4.

(Multiplication with Convex Hull) Given a matrix zonotope
$\boldsymbol{\mathcal{A}}\subset\mathbb{R}^{n\times n}$ and a polynomial
zonotope $\mathcal{PZ}\subset\mathbb{R}^{n}$, the multiplication between
${\textsc{conv}}(\boldsymbol{\mathcal{A}})$ and $\mathcal{PZ}$ can be computed
as

$\begin{split}{\textsc{conv}}\big{(}\boldsymbol{\mathcal{A}}\big{)}\,\mathcal{PZ}=\boldsymbol{\mathcal{A}}~{}\boldsymbol{\mathcal{L}}_{1}\,\mathcal{PZ}~{}\boxplus~{}{\textsc{fresh}}(\boldsymbol{\mathcal{A}})~{}\boldsymbol{\mathcal{L}}_{2}~{}\mathcal{PZ}\end{split}$

with

$\boldsymbol{\mathcal{L}}_{1}=\big{\langle}0.5\,I_{n},0.5\,I_{n},{\textsc{uniqueID}}(1)\big{\rangle}_{MZ},~{}~{}\boldsymbol{\mathcal{L}}_{2}=I_{n}-\boldsymbol{\mathcal{L}}_{1}$
(16)

which can be evaluated using Prop. 1.

###### Proof.

According to the definition of the convex hull of a matrix set in (15) we have

$\begin{split}\textsc{conv}(\boldsymbol{\mathcal{A}})&=\big{\\{}\lambda\,A_{1}+(1-\lambda)\,A_{2}~{}\big{|}~{}A_{1},A_{2}\in\boldsymbol{\mathcal{A}},~{}\lambda\in[0,1]\big{\\}}\\\
&=\boldsymbol{\mathcal{A}}~{}\boldsymbol{\mathcal{L}}_{1}~{}\boxplus~{}\textsc{fresh}(\boldsymbol{\mathcal{A}})\,\big{(}I_{n}-\boldsymbol{\mathcal{L}}_{1}\big{)},\end{split}$

where the matrix zonotope $\boldsymbol{\mathcal{L}}_{1}$ defined as in (16)
represents the domain $\lambda\in[0,1]$. While we need to apply the fresh
operation to $\boldsymbol{\mathcal{A}}$ to ensure that we realize two
independent matrices $A_{1},A_{2}\in\boldsymbol{\mathcal{A}}$ as demanded by
the definition in (15), we use the exact sum to keep the dependencies between
$\boldsymbol{\mathcal{L}}_{1}$ and $\boldsymbol{\mathcal{L}}_{2}$. ∎

Combining (14), Prop. 4, and the concept of preserving dependency on time used
in Sec. 4, the homogeneous solution for the first time interval can be
calculated as follows:

$\displaystyle\mathcal{H}(\tau_{0})$
$\displaystyle=\boldsymbol{\mathcal{M}}(\tau_{0})\,\mathcal{X}_{0}$ (17)
$\displaystyle\subseteq\big{(}\boldsymbol{\mathcal{E}}_{\Delta
t}\textsc{zonotope}(\mathcal{X}_{0})\big{)}~{}\oplus\bigg{(}\operatorname*{\scalerel*{\boxplus}{\sum}}_{i=0}^{\kappa}\frac{1}{i!}\textsc{fresh}^{i}(\boldsymbol{\mathcal{A}})\,\boldsymbol{\mathcal{L}}_{1}\,\boldsymbol{\mathcal{T}}^{i}~{}\mathcal{X}_{0}\boxplus\frac{1}{i!}\textsc{fresh}^{i}(\boldsymbol{\mathcal{A}})\,\boldsymbol{\mathcal{L}}_{2}\,\boldsymbol{\mathcal{T}}^{i}~{}\mathcal{X}_{0}\bigg{)},$
(18)

where $\boldsymbol{\mathcal{L}}_{1}$ and $\boldsymbol{\mathcal{L}}_{2}$ are
given by $\eqref{def: lambda matrix zonotope}$, $\boldsymbol{\mathcal{T}}$ is
given by $\eqref{def: time matrix zonotope}$, and
$\boldsymbol{\mathcal{E}}_{\Delta t}$ is given by (10). Based on (Althoff et
al., 2011b, Thm. 2), we can similarly enclose the particular solution at time
$t$ by:

$\begin{split}&\mathcal{P}(t)\subseteq\bigg{(}\operatorname*{\scalerel*{\boxplus}{\sum}}_{i=0}^{\kappa}\frac{t^{i+1}}{(i+1)!}\textsc{fresh}^{i}(\boldsymbol{\mathcal{A}})\,\boldsymbol{\mathcal{L}}_{1}\,\textsc{fresh}(\boldsymbol{\mathcal{B}}\,\mathcal{U})~{}\boxplus\\\
&\frac{t^{i+1}}{(i+1)!}\textsc{fresh}^{i}(\boldsymbol{\mathcal{A}})\,\boldsymbol{\mathcal{L}}_{2}\,\textsc{fresh}(\boldsymbol{\mathcal{B}}\,\mathcal{U})\bigg{)}\oplus\bigg{(}\frac{t}{\kappa+2}\,\boldsymbol{\mathcal{E}}_{t}\,\big{\langle}\widetilde{u},[~{}]\big{\rangle}_{Z}\bigg{)},\end{split}$
(19)

where $\widetilde{u}=[\widetilde{u}_{1}~{}\ldots~{}\widetilde{u}_{n}]^{T}$
with
$\widetilde{u}_{j}=\max\\{|u_{(j)}|~{}|~{}u\in~{}\boldsymbol{\mathcal{B}}\,\mathcal{U}\\}$,
$j=1,\dots,n$. Moreover, based on (LABEL:def:time_varying_input) the
particular solution for the first time interval is then given as

$\begin{split}\mathcal{P}(\tau_{0})\subseteq&~{}\big{(}\Delta
t/(\kappa+2)\,\boldsymbol{\mathcal{E}}_{\Delta
t}\,\langle\widetilde{u},[~{}]\rangle_{Z}\big{)}~{}\oplus\\\
&~{}\bigg{(}\operatorname*{\scalerel*{\boxplus}{\sum}}_{i=0}^{\kappa}\frac{1}{(i+1)!}\textsc{fresh}^{i}(\boldsymbol{\mathcal{A}})\,\boldsymbol{\mathcal{T}}^{i+1}\,\boldsymbol{\mathcal{L}}_{1}\,\textsc{fresh}(\boldsymbol{\mathcal{B}}\,\mathcal{U})\\\
&\hskip
20.0pt\boxplus\frac{1}{(i+1)!}\textsc{fresh}^{i}(\boldsymbol{\mathcal{A}})\,\boldsymbol{\mathcal{T}}^{i+1}\,\boldsymbol{\mathcal{L}}_{2}\,\textsc{fresh}(\boldsymbol{\mathcal{B}}\,\mathcal{U})\bigg{)}.\end{split}$
(20)

Finally, a tight enclosure of the reachable set can be computed using Alg. 2,
which is a modification of (Althoff et al., 2011b, Alg. 1). Note that since we
have to destroy dependencies on $\boldsymbol{\mathcal{A}}$ and
$\boldsymbol{\mathcal{B}}$ to ensure correctness, we can use a simpler
propagation scheme compared to Alg. 1.

Algorithm 2 Reachability algorithm (time-varying parameters)

###  5.2 Linear Time-Varying Systems

Based on the reachability algorithm for linear systems with time-varying
parameters presented in the previous subsection, we can also compute tight
enclosures of reachable sets for time-varying linear systems. In particular,
given a time-varying linear system

$\dot{x}(t)=A(t)\,x(t)+B(t)\,u(t)$

with $A:~{}\mathbb{R}_{\geq 0}\to\mathbb{R}^{n\times n}$ and
$B:~{}\mathbb{R}_{\geq 0}\to\mathbb{R}^{n\times m}$, we can for each
reachability time step $\tau_{k}$ abstract it by the linear parametric system

$\forall
t\in\tau_{k}:~{}~{}\dot{x}(t)\in\boldsymbol{\mathcal{A}}\,x(t)+\boldsymbol{\mathcal{B}}\,u(t),$

where

$\boldsymbol{\mathcal{A}}=\big{\\{}A(t)~{}\big{|}~{}t\in\tau_{k}\big{\\}},~{}~{}\boldsymbol{\mathcal{B}}=\big{\\{}B(t)~{}\big{|}~{}t\in\tau_{k}\big{\\}}.$
(21)

To tightly enclose the matrix sets $\boldsymbol{\mathcal{A}}$ and
$\boldsymbol{\mathcal{B}}$ in (21) by matrix zonotopes, range bounding using
affine arithmetic de Figueiredo and Stolfi (2004) can be applied. Finally, we
can simply use Alg. 2 to compute the reachable set, where we update the
matrices $\boldsymbol{\mathcal{A}}$ and $\boldsymbol{\mathcal{B}}$ according
to (21) in each time step. We need to use the algorithm for time-varying
parameters instead of Alg. 1 since $A(t)$ and $B(t)$ change over time.

Figure 5: Reachable set for the linear system from Example 3 (red), where the
initial set is depicted in white with a black border.

###  5.3 Linear Time-Invariant Systems

Even though reachability analysis for linear time-invariant systems
$\dot{x}(t)=A\,x(t)+B\,u(t)$ is a well studied topic, our approach provides
unique advantages for this problem. As summarized in a recent survey Forets
and Schilling (2022), existing methods compute the homogeneous time interval
solution either via optimization Frehse and et al. (2011) or by calculating
the convex hull of the time-point solutions Girard (2005); Le Guernic and
Girard (2010), where an additional correction term is used to account for the
trajectory curvature. All of these methods are only accurate if the time step
size $\Delta t$ is small. With our approach, on the other hand, we can compute
the homogeneous time interval solution very elegantly by simply evaluating
$\mathcal{H}(\tau_{0})=e^{A\boldsymbol{\mathcal{T}}}\,\mathcal{X}_{0}$ with
$\boldsymbol{\mathcal{T}}$ from (8) using Prop. 3. Moreover, the
representation size required to compute the set with an arbitrary user-defined
precision remains moderate even for extremely large time intervals, as we
demonstrate by the following lemma:

######  Lemma 1.

Fix an arbitrary error bound $\psi>0$. Then the number of Taylor terms
$\kappa$ required to ensure that the over-approximation error for the
enclosure of the set $\mathcal{H}([0,\Delta
t])=e^{A\boldsymbol{\mathcal{T}}}\,\mathcal{X}_{0}$ (calculated using Prop. 3)
stays below $\psi$ only grows linearly with respect to the time step size
$\Delta t$.

###### Proof.

According to Prop. 3, the enclosure is calculated as follows:

$e^{A\boldsymbol{\mathcal{T}}}\,\mathcal{X}_{0}\subseteq\underbrace{\bigg{(}\operatorname*{\scalerel*{\boxplus}{\sum}}_{i=0}^{\kappa}\frac{A^{i}\boldsymbol{\mathcal{T}}^{i}}{i!}~{}\mathcal{X}_{0}\bigg{)}}_{\text{exact}}\oplus\underbrace{\Big{(}\boldsymbol{\mathcal{E}}~{}\textsc{zonotope}(\mathcal{X}_{0})\Big{)}}_{\begin{subarray}{c}\text{over-
approximation}\\\ \vspace{-4pt}\\\ \text{error}\end{subarray}}.$ (22)

Since the first term is evaluated exactly using the exact sum, the over-
approximation error is given by the second term. Note that
$||\boldsymbol{\mathcal{T}}||_{\infty}=\Delta t$ holds since
$\boldsymbol{\mathcal{T}}=[0,\Delta t]\cdot I_{n}$ according to (8).
Therefore, the interval matrix $\boldsymbol{\mathcal{E}}$ is according to
Prop. 3 defined as

$\boldsymbol{\mathcal{E}}=\langle-\boldsymbol{1}_{n\times
n},\boldsymbol{1}_{n\times n}\rangle_{IM}\frac{||A||_{\infty}^{\kappa+1}\Delta
t^{\kappa+1}}{(\kappa+1)!}\frac{1}{1-\epsilon}$ (23)

with the convergence requirement for the remainder term, such that

$\epsilon=\frac{\|A\|_{\infty}\Delta t}{\kappa+2}<1.$

For any $\psi>0$, define $\zeta=\min(\frac{\psi}{2},0.5)$. By selecting
$\kappa$ to satisfy

$\kappa>\frac{\|A\|_{\infty}\Delta t}{\zeta\frac{1}{e}}-1,$ (24)

we obtain

$(\kappa+1)!\overset{(a)}{>}\left(\frac{1}{e}\right)^{\kappa+1}(\kappa+1)^{\kappa+1}\overset{\eqref{eq:condition}}{>}\frac{\|A\|_{\infty}^{\kappa+1}\Delta
t^{\kappa+1}}{\zeta^{\kappa+1}},$

where (a) follows from Stirling’s approximation,

$\frac{1}{e}\cdot(\kappa+1)<\frac{1}{e}(2\pi(\kappa+1))^{\frac{1}{2(\kappa+1)}}\cdot
e^{\frac{1}{12(\kappa+1)^{2}+(\kappa+1)}}\cdot(\kappa+1)<\sqrt[\kappa+1]{(\kappa+1)!}.$

This leads to that

$\frac{\|A\|_{\infty}^{\kappa+1}\Delta
t^{\kappa+1}}{(\kappa+1)!}<\zeta^{\kappa+1}<\zeta.$

Additionally, considering the definition of $\epsilon$, we find

$\epsilon=\frac{\|A\|_{\infty}\Delta t}{\kappa+2}<\frac{\|A\|_{\infty}\Delta
t}{\kappa+1}\overset{\eqref{eq:condition}}{<}\frac{\zeta}{e}<\frac{1}{2}<1,$

which implies that $\frac{1}{1-\epsilon}<2$. As a result, it leads to the
desired result

$\frac{||A||_{\infty}^{\kappa+1}\Delta
t^{\kappa+1}}{(\kappa+1)!}\frac{1}{1-\epsilon}<\psi.$

∎

Since the matrix zonotope $\boldsymbol{\mathcal{T}}$ has only a single
generator according to (8), it holds that the enclosure for
$\mathcal{H}(\tau_{0})=e^{A\boldsymbol{\mathcal{T}}}\,\mathcal{X}_{0}$
calculated using Prop. 3 has $(\kappa+1)h$ dependent generators if
$\mathcal{X}_{0}$ has $h$ dependent generators. According to Lemma 1, the
number of generators required to enclose the homogeneous solution
$\mathcal{H}([0,\Delta t])$ with a certain precision therefore only grows
linearly with respect to the time step. Consequently, we can make the time
step size $\Delta t$ very large (for example, the entire time bound) since the
representation size grows slowly.

######  Example 3.

We consider the linear system

$\begin{bmatrix}\dot{x}_{1}\\\ \dot{x}_{2}\\\
\dot{x}_{3}\end{bmatrix}=\begin{bmatrix}0&0&-0.9\\\ 0&-0.5&0\\\
0.9&0&0\end{bmatrix}\begin{bmatrix}x_{1}\\\ x_{2}\\\ x_{3}\end{bmatrix}$

together with the initial set $\mathcal{X}_{0}=[1,2]\times[3,4]\times[1,2]$.
As shown in Fig. 5, with our approach we can represent the reachable set for
the entire time horizon of $t_{end}=5$s by a single polynomial zonotope.

###  5.4 Hybrid Systems

Hybrid systems divide the state space into several regions with different
dynamic equations, where the regions are separated by so-called guard sets.

Figure 6: The left plot visualizes the intersection (red) between the
reachable set of the linear system from Example 3 (gray) and the guard set in
(25) (green), where the initial set is depicted in white with a black border.
The resulting guard intersection enclosure computed with different set
representations visualized on the plane defined by the guard set is shown on
the right.

As a simple example for a hybrid system we consider the linear system from
Example 3 with the guard set

$\mathcal{G}=\big{\\{}x\in\mathbb{R}^{3}~{}\big{|}~{}x_{(2)}-0.2\,x_{(3)}=1.1\big{\\}},$
(25)

which is visualized in Fig. 6. The common procedure for reachability analysis
of hybrid systems is to use a reachability algorithm for continuous dynamics
to compute the reachable set until it has fully crossed the guard set and is
therefore located outside of the current region. The initial set for the next
region is then determined by calculating the intersection between the
reachable set and the guard set, before the whole procedure is repeated until
the time horizon is reached. Handling those guard intersections efficiently is
one of the main challenges for hybrid system reachability analysis. In
particular, most reachability algorithms for continuous systems require a very
small time step size $\Delta t$ to provide accurate results. Consequently,
often the reachable sets for many consecutive time intervals will intersect
the guard set. If the intersection between the guard set and each of these
sets would be treated as a separate initial set for the next region, one would
have to solve not only one but many reachability problems, which is
computationally demanding. Moreover, since usually multiple consecutive guard
intersections occur during the overall time horizon, the number of parallel
sets would grow exponentially in this case Duggirala and Bak (2019). To avoid
these issues, existing approaches for hybrid system reachability unite the
sets for different time intervals by enclosing the corresponding guard
intersections by convex set representations such as template polyhedra Girard
and Le Guernic (2008), zonotopes Althoff and Krogh (2012), or zonotope bundles
Althoff and Krogh (2011). However, as shown in Fig. 6, even for very simple
linear dynamics the intersection between the reachable set and the guard set
often yields a highly non-convex set, so that convex enclosures are often very
conservative.

Figure 7: Reachable set enclosure for the Lotka-Volterra system in Example 4
at time $t=0.25$s (left) and time $t=0.5$s (right).

As we discussed in Sec. 5.3, with our reachability approach we can represent
the reachable set for the entire time horizon with a single set. Consequently,
with our approach one can completely avoid the unification step when computing
guard intersections. Moreover, by converting the polynomial zonotope that
represents the reachable set to a constrained polynomial zonotope Kochdumper
and Althoff (2020), the intersection between the reachable set and the guard
set can be computed exactly as a single set, even for nonlinear guard set such
as polynomial level sets (Kochdumper, 2022, Prop. 3.2.24). Due to the
similarity of polynomial zonotopes and constrained polynomial zonotopes, the
reachability algorithm presented in this paper can easily be extended to
compute with constrained polynomial zonotopes, so that we can simply use the
resulting constrained polynomial zonotope as the initial set for the next
region. Finally, checking if a polynomial zonotope intersects a guard set can
use the computationally efficient polynomial zonotope refinement approach Bak
and et. al (2022).

###  5.5 Nonlinear Systems

The approach presented in this paper can also be used to construct a new
reachability algorithm for nonlinear systems. This algorithm linearizes the
system in each time step and captures the linearization error with uncertain
parameters in order construct a guaranteed enclosure of the reachable set. In
particular, given the nonlinear system

$\dot{x}(t)=f\big{(}x(t),u(t)\big{)}$

with $f:~{}\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}^{n}$, we can for
each reachability time step $\tau_{k}$ enclose it by a linear parametric
system using a first-order Taylor series expansion of the function $f$:

$\forall t\in\tau_{k}:~{}~{}\dot{x}(t)\in
f(x_{l},u_{l})+\boldsymbol{\mathcal{A}}\big{(}x(t)-x_{l}\big{)}+\boldsymbol{\mathcal{B}}\big{(}u(t)-u_{l}\big{)},~{}$

where

$\begin{split}\boldsymbol{\mathcal{A}}=&\bigg{\\{}\frac{\partial
f(x,u)}{\partial
x}~{}\bigg{|}~{}x\in\mathcal{R}(\tau_{k}),u\in\mathcal{U}\bigg{\\}},\\\
\boldsymbol{\mathcal{B}}=&\bigg{\\{}\frac{\partial f(x,u)}{\partial
u}~{}\bigg{|}~{}x\in\mathcal{R}(\tau_{k}),u\in\mathcal{U}\bigg{\\}}.\end{split}$
(26)

The expansion points $x_{l}$ and $u_{l}$ can be heuristically chosen as
$x_{l}=x_{c}+0.5\,\Delta tf(x_{c},u_{c})$ and $u_{l}=u_{c}$, where $x_{c}$ is
the center of $\mathcal{R}(t_{k})$ and $u_{c}$ is the center of $\mathcal{U}$.
To tightly enclose the matrix sets $\boldsymbol{\mathcal{A}}$ and
$\boldsymbol{\mathcal{B}}$ in (26) by matrix zonotopes, we can apply range
bounding using affine arithmetic de Figueiredo and Stolfi (2004). One problem
we are facing here is that we require the time interval reachable set
$\mathcal{R}(\tau_{k})$ to compute $\boldsymbol{\mathcal{A}}$ and
$\boldsymbol{\mathcal{B}}$, but also require $\boldsymbol{\mathcal{A}}$ and
$\boldsymbol{\mathcal{B}}$ to compute $\mathcal{R}(\tau_{k})$. To resolve this
mutual dependence, we apply the following strategy Althoff (2013): We compute
$\boldsymbol{\mathcal{A}}$ and $\boldsymbol{\mathcal{B}}$ using an estimate
$\widehat{\mathcal{R}}(\tau_{k})$ for $\mathcal{R}(\tau_{k})$, and then
compute $\mathcal{R}(\tau_{k})$ using the resulting $\boldsymbol{\mathcal{A}}$
and $\boldsymbol{\mathcal{B}}$. If the resulting set satisfies
$\mathcal{R}(\tau_{k})\subseteq\widehat{\mathcal{R}}(\tau_{k})$, we have a
guaranteed enclosure of the reachable set. If not, we update our estimate by
bloating the resulting set
$\widehat{\mathcal{R}}(\tau_{k})=x_{c}+\eta\,(\mathcal{R}(\tau_{k})-x_{c})$,
where $\eta>1$ is a bloating factor and $x_{c}$ is the center of
$\mathcal{R}(\tau_{k})$. This process is repeated until
$\mathcal{R}(\tau_{k})\subseteq\widehat{\mathcal{R}}(\tau_{k})$ is satisfied
at some point. As an initial estimate $\widehat{\mathcal{R}}(\tau_{k})$ we can
use the final reachable set from the previous time step $\mathcal{R}(k\Delta
t)$. Since containment checks as well as affine arithmetic cannot be directly
applied to polynomial zonotopes, we use an interval or oriented hyper-
rectangle enclosure of $\mathcal{R}(\tau_{k})$ to compute
$\boldsymbol{\mathcal{A}}$ and $\boldsymbol{\mathcal{B}}$. Moreover, we have
to use the algorithm for time-varying parameters in Alg. 2 to compute the
reachable set since the abstraction error represented by the matrix sets
$\boldsymbol{\mathcal{A}}$ and $\boldsymbol{\mathcal{B}}$ does not remain
constant over time. Let us demonstrate reachability analysis for nonlinear
systems by an example:

Figure 8: Comparison of reachable set enclosures for the Dubins car computed
with our approach and with the zonotope method Althoff et al. (2011a), where
the time interval reachable set is shown on the top and the final reachable
set is shown on the bottom.

######  Example 4.

We consider the Lotka-Volterra system

$\begin{split}\dot{x}_{1}&=3\,x_{1}-3\,x_{1}\,x_{2}\\\
\dot{x}_{2}&=x_{1}\,x_{2}-x_{2}\end{split}$

together with the initial set $\mathcal{X}_{0}=[1.1,1.5]\times 1$ and the time
horizon $t_{end}=0.5$\mathrm{s}$$. The results visualized in Fig. 7 show that
our approach is able to compute a non-convex enclosure of the reachable set.

In this section we compare our reachability algorithm to other state of the
art methods using two benchmark systems. We implemented our approach in
MATLAB, and all computations are carried out on a Apple M2 Pro processor with
16GB memory.

Figure 9: Comparison of enclosures for the final reachable set for the Dubins
car with time-varying parameters computed with our approach and with the
zonotope method Althoff et al. (2011b), where the exact reachable set for the
constant parameter case is shown for comparison. Figure 10: Comparison of the
reachable set enclosure of the platoon benchmark computed with our approach
and the CORA toolbox Althoff (2015), where random simulations are shown in
black and the required safe distance of $42$m is visualized in red.

###  5.6 Dubins Car

First, we consider the example of a Dubins car shown in Fig. 1 to demonstrate
how well our reachabilty algorithm captures the non-convexity of the reachable
set. According to (Bak and Tran, 2022, Sec. 2.2), the dynamics of the Dubins
car can be described as a four-dimensional linear parametric system. We
consider a time horizon of $t_{end}=2$\mathrm{s}$$ together with the initial
state
$\mathcal{X}_{0}=[0~{}0~{}0~{}10$\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$]$. Fig.
8 compares the reachable set enclosure for case with constant parameters
computed with our algorithm and the zonotope method Althoff et al. (2011a),
which is implemented in the reachablity toolbox CORA Althoff (2015). Our
polynomial zonotope approach tightly captures the non-convexity of the
reachable set, while the zonotope method yields a large convex over-
approximation. Moreover, for our approach we can compute the reachable set
using a single time step together with $\kappa=60$ Taylor terms and a desired
zonotope order of $\rho_{d}=60$, which takes $0.25$\mathrm{s}$$. For the
zonotope approach, we obtained the tightest results for a time step size of
$\Delta t=0.05$\mathrm{s}$$, $\kappa=6$ Taylor terms, and desired zonotope
order $\rho_{d}=100$, which takes $0.18$\mathrm{s}$$ to compute.

As a second experiment, we consider the scenario where the parameters are
time-varying. Since reachability analysis is more challenging in this case, we
reduce the time horizon to $t_{end}=1$\mathrm{s}$$. We compare the approach
presented in Sec. 5.1 to the zonotope method Althoff et al. (2011b), which
again is implemented in the reachability toolbox CORA Althoff (2015). For the
zonotope method we obtained the tightest result with the settings $\Delta
t=0.05$\mathrm{s}$$, $\kappa=4$, and $\rho_{d}=30$, for which the calculation
of the reachable set takes $0.069$\mathrm{s}$$. For our approach, we use the
same time step size but increase the other settings to $\kappa=6$, and
$\rho_{d}=80$. While the calculation of the reachable set takes with
$4$\mathrm{s}$$ longer than for the zonotope method, our approach yields an
enclosure that is in most regions more accurate as shown in Fig. 9. Moreover,
while the zonotope method only computes a very rough convex enclosure of the
reachable set, our approach still captures the non-convexity of the reachable
set to some extent.

###  5.7 Vehicle Platoon

To demonstrate the scalability of our approach, we now consider the
PLAA01-BND42 instance of the $9$-dimensional platoon benchmark from the 2021
ARCH competition (Althoff and et al., 2021, Sec. 3.7), which describes a
platoon consisting of four vehicles. The verification task for the benchmark
is to show that even in case of a communication loss the vehicles still keep a
safe distance of 42m from each other. This can be formulated as a reachability
problem, where the corresponding parametric system encloses both the system
with and without communication loss. The system has one uncertain input, which
is the acceleration of the leading vehicle. For our experiments we slightly
modify the benchmark by considering constant instead of time-varying
parameters. We compare our approach with the state of the art reachability
toolbox CORA Althoff (2015), which uses the zonotope method Althoff et al.
(2011a) and requires a computation time of $15$\mathrm{s}$$. For our approach,
calculation of the reachable set using Alg. 1 with $\Delta
t=0.036$\mathrm{s}$$, $\kappa=6$, and $\rho_{d}=50$ takes $9$\mathrm{s}$$, and
is therefore faster than the CORA toolbox. On top of that, as shown in Fig. 10
the reachable set enclosure computed with our approach is significantly
tighter than the one computed with CORA. In summary, our approach is therefore
both faster and more accurate than the existing state of the art reachability
tool.

##  6 Scalable Optimization for Multi-Affine Zonotopes

###  6.1 Reachability sets are multi-affine zonotopes

Let us reconsider Algorithm 2 from section 5.1. Let us make a (mild)
assumption that the matrix sets $\mathbfcal{A}$ and $\mathbfcal{B}$, initial
set $\mathcal{X}_{0}$, input set $\mathcal{U}$ are indenpendent (true if they
do not share a common uncertain factor). Then, the resulting polynomial
zonotope has a curiously useful property; the exponents over all the factors,
minus the time factor, are at most 1. This means that if the time factor is
replaced with numerical values (as they are, after all, known), then _the
resulting reachable set is a multiaffine zonotope (Def 2.6)._ Intuitively,
this is because we call fresh operator on the power of $\mathbfcal{A}$, and
because the particular solution is added through Minkowski sum for each time
interval.

Since multi-affine zonotopes are a type of polynomial zonotopes, it is of
course possible to apply existing set operations for polynomial zonotopes;
however, usually more efficient methods exist when restricted to this special
class. Theoretically, as shown in Huang et al. (2023), optimization and
intersection checking of polynomial zonotopes is NP-complete. This is also
true of plotting polynomial zonotopes of dimension 2. This is because both
tasks rely on a splitting algorithm Huang et al. (2023); Bak and Tran (2022)
which models polynomial zonotopes as a union of zonotopes, created by
sequentially splitting the original nonconvex set. Given the importance of
this splitting method, we now investigate how its efficiency can be improved
over multiaffine zonotopes.

Figure 11: Time point reachable set of the Dubins car with time-varying
parameters plotted through the default polynomial zonotope approach with 60
splits costing around 2 hours.

###  6.2 Splitting methods for polynomial zonotopes

The core idea behind splitting methods for polynomial zonotopes is the
observation that any nonconvex set can be represented (or overapproximated) as
a union of zonotopes. This is a useful representation, since tasks such as
intersection checking and 2-D plotting are much easier for zonotopes (can be
done through quick index search) than polynomial zonotopes (requires nonconvex
optimization). However, a drawback is that, even with careful implementation,
the number of zonotopes may increase significantly, and as mentioned in (Huang
et al., 2023), it computationally demanding to achieve an accurate
overrepresentation; in some implementations, it may even lead to an unbounded
number of sets by not monotonically decreasing the overrepresentation set size
in each round.

As we see in Fig.11, applying the built-in plotting algorithm for polynomial
zonotope relying on splitting to plot the same time instance set as shown in
Fig. 9 takes around 2 hours, and may still result in a bad over-approximation.
One idea to improve the plotting process is to model the convex hull of the
set as an intersection of hyperspaces by finding the supporting hyperplanes at
various points on the boundary of the set. These supporting hyperplanes can be
derived by solving multi-affine optimization problems.

###  6.3 Improved splitting methods for multiaffine zonotopes

In this section we propose a scalable optimization algorithm for minimizing
linear functions over multiaffine zonotopes. Let us start by defining the
multi-affine optimization problem

######  Definition 6.1.

Consider the polynomial defined as:

$\displaystyle
p(\alpha_{1},\alpha_{2},\cdots,\alpha_{n})=\sum_{I\subseteq\\{1,2,\cdots,n\\}}g_{I}\prod_{i\in
I}\alpha_{i},$ (27)

where $a_{I}\in\mathbb{R}$. The multi-affine polynomial optimization problem
is:

$\displaystyle\min_{\begin{subarray}{c}\alpha_{1},...,\alpha_{n}\\\
\alpha_{i}\in[-1,1]\end{subarray}}p(\alpha_{1},\alpha_{2},\cdots,\alpha_{n}).$
(28)

Given a set $C\subseteq\\{1,2,\ldots,n\\}$ and a multi-affine polynomial $p$,
as defined in Equation (27), we formalize the restriction of $p$ to the set
$C$, denoted $p|_{C}$, with the following expression:

$p|_{C}=\sum_{I\subseteq C}g_{I}\prod_{i\in I}\alpha_{i}$ (29)

Note that we can use multi-affine optimization to find a tight over-
approximation of multi-affine zonotope. Suppose we have a multi-affine
zonotope $\mathcal{MAZ}$. Let us start with a vector $\textbf{d}_{1}$. To find
the support, we minimize the following:

$\textbf{x}_{1}=\arg\min_{\textbf{x}\in\mathcal{MAZ}}\textbf{d}_{1}^{T}\textbf{x}$
(30)

Therefore, the half-space $\mathcal{H}_{1}$ outlined below include the multi-
affine zonotope, and the associated supporting hyperplane, denoted as
$\partial\mathcal{H}_{1}$, will be

$\mathcal{H}_{1}=\\{\textbf{x}~{}|~{}\textbf{d}_{1}^{T}(\textbf{x}_{1}-\textbf{x})\leq
0\\},\quad\partial\mathcal{H}_{1}=\\{\textbf{x}~{}|~{}\textbf{d}_{1}^{T}(\textbf{x}_{1}-\textbf{x})=0\\}$

Subsequently we can calculate $\mathcal{H}_{2}$. This iterative process will
continue until we calculate $\mathcal{H}_{m}$. The specific procedure of how
to choose $\textbf{d}_{1},\cdots,\textbf{d}_{m}$ is so called the Kamenev
method and being discussed in Lotov et al. (2004). Hence this can be
formalized to some incremental approach to find the directions to perform
optimization, until the consecutive optimal solution values are close to each
other under some small tolerance. As a result,
$\mathcal{H}_{1}\bigcap\mathcal{H}_{2}\cdots\bigcap\mathcal{H}_{m}$ will be an
over-approximation of the multi-affine zonotope. This type of convex over-
approximation that can be combined with other methods, like polynomial
zonotopes splitting. This combination can make the single polynomial zonotope
split method more effective, as we show in Section 6.4. By merging these
approaches, we can get a tighter approximation than using the method based on
polynomial zonotope splits alone. A visualization of the above procedure is in
figure 12.

Figure 12: In the diagram, we initiate with a predefined series of directions
identified as $d_{1}$. Subsequently, we identify a sequence of supports
$\partial\mathcal{H}_{1},\partial\mathcal{H}_{2},\ldots,\partial\mathcal{H}_{m}$
and the respective points of intersection
$\textbf{x}_{1},\textbf{x}_{2},\ldots,\textbf{x}_{m}$ with the multi-affine
zonotope, label as green dots. The red dots are used to mark the points where
adjacent supports intersect, specifically at
$\partial\mathcal{H}_{1}\bigcap\partial\mathcal{H}_{2},\partial\mathcal{H}_{2}\bigcap\partial\mathcal{H}3,\ldots,\partial\mathcal{H}_{m-1}\bigcap\partial\mathcal{H}_{m}$.

We now present our method for solving (28). Since all variables in a multi-
affine optimization problem have an exponent of one, the partial derivative
along each variable cannot change sign. This means that the optimal value must
occur on one of the corners of the $n$-dimensional box of the domain, and it
is sufficient to consider this finite set when optimizing.

$\displaystyle\min_{\begin{subarray}{c}\alpha_{1},...,\alpha_{n}\\\
\alpha_{i}\in[-1,1]\end{subarray}}p(\alpha_{1},\alpha_{2},\cdots,\alpha_{n})=\min_{\begin{subarray}{c}\alpha_{1},...,\alpha_{n}\\\
\alpha_{i}\in\\{-1\\}\cup\\{1\\}\end{subarray}}p(\alpha_{1},\alpha_{2},\cdots,\alpha_{n}).$

The above problem is, in general, NP-complete due to the following theorem:

######  Theorem 1.

(Huang et al., 2023) Optimization over bilinear polynomial optimization is NP-
complete where bilinear polynomial is defined as

$\displaystyle\min_{\begin{subarray}{c}\alpha_{1},...,\alpha_{n}\\\
\alpha_{i}\in\\{-1\\}\cup\\{1\\}\end{subarray}}\sum_{i=1}^{n}\sum_{j=i+1}^{n}g_{i,j}\alpha_{i}\alpha_{j}.$
(31)

We can design a brute force algorithm to find the optimal value of the multi-
affine polynomial.

Algorithm 3 Brute Force Minimization of Multi-Affine Polynomial

1:function Minimize($p$, $n$)

2: $m\leftarrow\infty$

3: for $i\leftarrow 0$ to $2^{n}-1$ do

4: for $j\leftarrow 1$ to $n$ do

5: $\alpha_{j}=2\left(\left\lfloor\frac{i}{2^{j}}\right\rfloor\mod 2\right)-1$

6: end for

7: $m\leftarrow\min(m,p(\alpha_{1},\alpha_{2},\ldots,\alpha_{n}))$

8: end for

9: return $m$

10:end function

The algorithm, as defined, has exponential running time, aligning with the NP-
Completeness of the problem at hand. In the worst case, it may check all
$2^{n}$ possible values for $\\{\alpha_{i}\\}_{i}^{n}$. However, it is
possible that some problems can be solved more efficiently. Specifically, the
work in (Del Pia and Di Gregorio, 2023) highlights that the complexity of the
multi-affine polynomial optimization problem is closely linked to the
associated dependency graph structure.

######  Definition 6.2.

Given a multiaffine polynomial $p$ as defined in (27), the associated factor
dependency graph $G=(V,E)$ is constructed as follows:

  * •

The vertex set $V=\\{1,...,n\\}$ represents the indices of the variables
$\alpha_{1},\alpha_{2},\ldots,\alpha_{n}$ in the polynomial $p$.

  * •

The edge set $E\subset V\times V$ contains an edge $e=\\{i,j\\}\in E$ if there
exists a monomial in the polynomial $p$ (27) that includes both coefficients
$a_{i}$ and $a_{j}$.

This graph represents the pairwise dependencies between variables in the
polynomial optimization problem, where an edge signifies that the connected
variables jointly affect the value of the polynomial.

Suppose we have a multi-affine polynomial $p$ and its associated dependency
graph $G=(V,E)$. If $G$ is disconnected, e.g. consisting of disconnected
components
$G_{1}=(V_{1},E_{1}),G_{2}=(V_{2},E_{2}),\ldots,G_{m}=(V_{m},E_{m})$, (e.g.
where there is no edge $\\{i,j\\}\in E$ where $i\in V_{k}$, $j\in V_{l}$,
$k\neq l$), then the optimization problem can be decomposed over each
connected component:

$\min_{\alpha_{i}\in[-1,1],i\in V}p=\min_{\alpha_{i}\in[-1,1],i\in
V_{1}}p_{1}+\min_{\alpha_{i}\in[-1,1],i\in
V_{2}}p_{2}+\cdots+\min_{\alpha_{i}\in[-1,1],i\in V_{m}}p_{m}$

Here, $p_{i}$ only contains the monomials from $p$ that involve coefficients
represented by nodes in $V_{i}$. Then, the runtime reduces from $2^{n}$ to
$\sum_{i}2^{|V_{i}|}$, which is often significant in practice. Let’s
illustrate this with the following example:

######  Example 5.

Consider a multi-affine polynomial $p$ and its associated dependency graph
$G=(V,E)$ defined as follows:

$\displaystyle p$
$\displaystyle=g_{1}\alpha_{1}\alpha_{2}\alpha_{3}+g_{2}\alpha_{3}\alpha_{4}+g_{3}\alpha_{5}\alpha_{6}+g_{4}\alpha_{7},$
$\displaystyle V$
$\displaystyle=\\{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4},\alpha_{5},\alpha_{6},\alpha_{7}\\},$
$\displaystyle E$
$\displaystyle=\\{\\{\alpha_{1},\alpha_{2}\\},\\{\alpha_{1},\alpha_{3}\\},\\{\alpha_{2},\alpha_{3}\\},\\{\alpha_{3},\alpha_{4}\\},\\{\alpha_{5},\alpha_{6}\\}\\}.$

The dependency graph $G$ can be divided into the following connected
components:

  * •

$G_{1}=(V_{1},E_{1})$ with
$V_{1}=\\{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4}\\}$ and
$E_{1}=\\{\\{\alpha_{1},\alpha_{2}\\},\\{\alpha_{1},\alpha_{3}\\},\\{\alpha_{2},\alpha_{3}\\},\\{\alpha_{3},\alpha_{4}\\}\\}$,

  * •

$G_{2}=(V_{2},E_{2})$ with $V_{2}=\\{\alpha_{5},\alpha_{6}\\}$ and
$E_{2}=\\{\\{\alpha_{5},\alpha_{6}\\}\\}$,

  * •

$G_{3}=(V_{3},E_{3})$ consisting of the isolated vertex
$V_{3}=\\{\alpha_{7}\\}$ with no edges.

The polynomial $p=p_{1}+p_{2}+p_{3}$ can be decomposed into sub-polynomials
$p_{1}$, $p_{2}$, and $p_{3}$ corresponding to these components:

$\displaystyle p_{1}$
$\displaystyle=g_{1}\alpha_{1}\alpha_{2}\alpha_{3}+a_{2}\alpha_{3}\alpha_{4},$
(32) $\displaystyle p_{2}$ $\displaystyle=g_{3}\alpha_{5}\alpha_{6},$ (33)
$\displaystyle p_{3}$ $\displaystyle=g_{4}\alpha_{7}.$ (34)

Thus, the original optimization problem can be decomposed into:

$\displaystyle\min_{\alpha_{i}\in[-1,1],i\in
V}p=\min_{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4}\in[-1,1]}p_{1}+\min_{\alpha_{5},\alpha_{6}\in[-1,1]}p_{2}+\min_{\alpha_{7}\in[-1,1]}p_{3}.$

Without decomposition, Alg. 3 would enumerate $2^{7}=128$ possibilities;
however, through decomposition, we only need to enumerate $2^{4}+2^{2}+2=22$
possibilities, which is a significant reduction.

Now suppose an associated graph is _almost_ disconnected; for example, there
is only one vertex, that, if removed, would disconnect the graph. Then we
propose to remove this vertex so that the dependency graph can be partitioned
into multiple connected components. Let us illustrate it by considering the
following example:

######  Example 6.

Consider the multi-affine polynomial optimization problem defined by the
polynomial:

$p(\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4},\alpha_{5},\alpha_{6},\alpha_{7},\alpha_{8})=(g_{1}\alpha_{1}\alpha_{2}\alpha_{3}+a_{2}\alpha_{3}\alpha_{4}+g_{3}\alpha_{5}\alpha_{6},+g_{4}\alpha_{7})\alpha_{8}.$

where $g_{1},g_{2},g_{3},$ and $g_{4}$ are real coefficients. The associated
dependency graph $G=(V,E)$ for this polynomial is given by:

  * •

$V=\\{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4},\alpha_{5},\alpha_{6},\alpha_{7},\alpha_{8}\\}$

  * •

$E=\\{\\{\alpha_{1},\alpha_{2}\\},\\{\alpha_{1},\alpha_{3}\\},\\{\alpha_{2},\alpha_{3}\\},\\{\alpha_{3},\alpha_{4}\\},\\{\alpha_{5},\alpha_{6}\\}\\{\alpha_{1},\alpha_{8}\\},$
$\\{\alpha_{2},\alpha_{8}\\},\\{\alpha_{3},\alpha_{8}\\},\\{\alpha_{4},\alpha_{8}\\},\\{\alpha_{5},\alpha_{8}\\},\\{\alpha_{6},\alpha_{8}\\},\\{\alpha_{7},\alpha_{8}\\}\\}$

In this dependency graph, vertex $\alpha_{8}$ is a key node connecting
multiple edges. The removal of $\alpha_{8}$ leads to partitioning the graph
into separate connected components. The remaining graph, after removing
$\alpha_{8}$, can be separated into three distinct connected components

$V_{1}=\\{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4}\\},\quad
V_{2}=\\{\alpha_{5},\alpha_{6}\\}\quad V_{3}=\\{\alpha_{7}\\}.$

More specifically, if $\alpha_{8}=1$

$\displaystyle\min_{\alpha_{i}\in[-1,1],i\in
V}p_{+}=\min_{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4}\in[-1,1]}p_{1}+\min_{\alpha_{5},\alpha_{6}\in[-1,1]}p_{2}+\min_{\alpha_{7}\in[-1,1]}p_{3}.$

where

$p_{1}=g_{1}\alpha_{1}\alpha_{2}\alpha_{3}+g_{2}\alpha_{3}\alpha_{4},\qquad
p_{2}=g_{3}\alpha_{5}\alpha_{6},\qquad p_{3}=g_{4}\alpha_{7}.$

Similarly, if $\alpha_{8}=-1$, we simply flip signs where needed, and minimize

$\displaystyle\min_{\alpha_{i}\in[-1,1],i\in
V}p_{-}=\min_{\alpha_{1},\alpha_{2},\alpha_{3},\alpha_{4}\in[-1,1]}(-p_{1})+\min_{\alpha_{5},\alpha_{6}\in[-1,1]}(-p_{2})+\min_{\alpha_{7}\in[-1,1]}(-p_{3}).$

Then the minimum of $p$ is $\min\\{p_{+},p_{-}\\}$.

(a) Factor dependency graph for Example 5

(b) Factor dependency graph for Example 6

Figure 13: The visualization shows the dependency graph. By removing the node
$\alpha_{8}$ from the graph on the right, we see that the graph is split into
three separate parts, as shown on the left.

Given the above example, we design the Algorithm 4 and 5 to solve the multi-
affine optimization problem efficiently.

Algorithm 4 Minimize a Multi-affine Polynomial Based on Splitting

1:function SplitMin($p$, $g$, $t$)

2: Input: A multi-affine polynomial $p$, associated dependency graph
$g=(v,e)$, threshold $t$.

3: Output: Optimal value of polynomial $p$.

4: if $|v|\leq t$ then

5: return $\textsc{Minimize}(p)$

6: end if

7: $m\leftarrow 0$

8: $gc\leftarrow\textsc{ConnectedComponents}(g)$

9: for all $c\in gc$ do

10: $m\leftarrow m+\textsc{MinOneComponent}(p,g,c,t)$

11: end for

12: return $m$

13:end function

Algorithm 5 Function for Minimizing Over a Single Component

1:function MinOneComponent($p$, $g$, $c$, $t$)

2: Input: A multi-affine polynomial $p$, associated dependency graph
$g=(v,e)$, threshold $t$ and a connected component $c$ from graph $g$.

3: Output: Optimal value of polynomial $p|_{c}$ defined in (29)

4: if $|c|\leq t$ then

5: return $\textsc{Minimize}(p|_{c})$

6: else

7: $ct\leftarrow\textsc{FindMinVertexCuts}(g,c)$555The algorithm returns the
minimum vertex cut West et al. (2001) of the connected component
$c$.$\triangleright$ Cut on the connected component c from g

8: Initialize a stack $s$ with $p|_{c}$ at the top

9: for all $nd\in ct$ do $\triangleright$ Each vertex return by the cut

10: New stack $s^{\prime}$ is empty

11: while $s$ is not empty do

12: $p^{\prime}\leftarrow s.pop()$

13: for $i\in\\{-1,1\\}$ do

14: $p^{\prime}\leftarrow\textsc{Substitute}(p^{\prime},nd,i)$

15: $s^{\prime}.push(p^{\prime})$

16: end for

17: end while

18: $s\leftarrow s^{\prime}$

19: end for

20: $m\leftarrow\infty$

21: while $s$ is not empty do

22: $p^{\prime}\leftarrow s.pop()$

23: $g^{\prime}\leftarrow\textsc{RemoveVertices}(g,ct)$

24: $m\leftarrow\min(m,\textsc{SplitMin}(p^{\prime},g^{\prime},t))$

25: end while

26: return $m$

27: end if

28:end function

###  6.4 Evaluation

In this section, we perform two experiments to demonstrate how the multi-
affine zonotope specialized optimization algorithm 4 and 5 can help relax the
difficulties existing in several crucial operations, which are based on the
splitting stragety Bak and Tran (2022), of polynomial zonotopes, including 2-d
plotting, support function computation and intersection checking. We compare
our method with the newest CORA v2024.1.3, which is a state-of-the-art
polynomial zonotopes tool box. Our optimization method is implemented in
Python 3.9 and experiments run on an Apple M2 Pro processor with 16GB memory.
666In contrast, the reachability evaluation in section 5 compares against CORA
toolbox which is the one at the time of the HSCC conference paper getting
reviewed, for the purpose of consistency.

For the first experiment, we notice that in figure 11 the 2-d reachable set
projection plotted by the CORA built-in splitting based approach for
polynomial zonotope does not provide a tight over-approximation of the
reachable set after running around 2 hours. Instead, as we notice that time-
point reachable set is actually a multi-affine zonotope, we can apply our
multi-affine optimization approach to provide a tighter over-approximation
which is given in figure 14. The grey region corresponds to CORA polynomial
zonotope plotting approach with 60 splits, which takes time more than 2 hours.
The blue contour corresponds to getting the convex enclosure taking the idea
from Sec. 6.3 and apply through our multi-affine optimization approach. After
that we intersect our convex contour with the non-convex over-approximated
polygon obtained through the CORA plotting and splitting process with 20
splits. The computation time of getting the convex contour, including the
Kamenev approach to find the direction to optimize with tolerance of
$10^{-7}$, and the graph based multi-affine optimization approach itself,
combined with the computation time of the CORA built-in function with 20
splits takes $1063s$ in total. Through this experiment, we demonstrate from a
visualized perspective that the splitting based algorithm would cost huge
computation time when the split number is high for the purpose of getting a
tighter over-approximation enclosure. However, the reality is that though huge
computation budget is spent, the gain is limited.

In this second experiment, we are going to make the same argument about
splitting based approaches, which have the imbalance of the accuracy
improvement with the high computation cost, through numerical evaluation of
the supports functions. To demonstrate the inadequacies of the splitting
approach not only exist in the time point reachable set from the Dubin’s car
dynamic with time-varying parameters, but also is a more general problem. We
experiment over a different 5-D benchmark system with time-varying parameters
from section 5.1 of (Althoff et al., 2011a). We perform reachability analysis
through Algorithm. 2 by using a high zonotope order of 500. Then we get the
third time point reachable set, and reduce it to a polynomial zonotope (which
is multi-affine) with $495$ dependent generators, $5$ independent generators,
$96$ dependent factors and project to the first two dimension as our
experiment target. The projection would not prevent us from showing the
scalability of our result since one always need to project the set
representation onto some direction to perform optimization process (by eq.
30).

We run the incremental Kamenev approach to find 6 directions to optimize with
tolerance of $10^{-2}$. Then test our optimization approach for each direction
and comparing against the CORA splitting based support function computation
approach for $10$, $20$, $30$ splits. The results are given in Table 1. Since
multi-affine optimization always returns the optimal solution, we set the
optimal values we get as a comparison standard $f^{*}$ for the CORA method.
Let the optimal value given by CORA support function as $f$. We define the
relative error as $\frac{|f-f^{*}|}{|f^{*}|}$. In our study, we found that
when we increase the number of splits in the support function, which makes the
computation time longer, the improvement in error reduction is not as
significant as expected. This result shows a clear imbalance between the high
computational cost from more splits and the small gains in accuracy. On the
other hand, directly apply our multi-affine optimization approach, guarantees
to find the optimal solution while having very few computation time.

Figure 14: Time point reachable set of the Dubins car dynamics with time
varying parameters. The grey region corresponds to CORA polynomial zonotope
plotting approach with 60 splits, takes time more than 2 hours. The blue line
corresponds to by using our multi-affine optimization approach with polynomial
approach with 20 splits. Table 1: Relative error and Running Time for
different methods.

Method | Case 1 (Error/ Time) | Case 2 (Error/ Time) | Case 3 (Error/ Time)
---|---|---|---
10 Splits |  $0.41\%$ / 0s |  $0.6\%$ / 1s |  $0.49\%$ / 1s
20 Splits |  $0.36\%$ / 9s |  $0.54\%$ / 54s |  $0.49\%$ / 344s
30 Splits |  $0.34\%$ / 227s |  $0.52\%$ / 5818s |  $0.49\%$ / 3087s
Ours |  0 / 16s |  0 / 16s |  0 / 15s

Method | Case 4 (Error / Time) | Case 5 (Error/ Time) | Case 6 (Error/ Time)
---|---|---|---
10 Splits |  $0.98\%$ / 0s |  $0.74\%$ / 1s |  $0.79\%$ / 0s
20 Splits |  $0.88\%$ / 43s |  $0.72\%$ / 358s |  $0.64\%$ / 330s
30 Splits |  $0.83\%$ / 1156s | Exceeds limit |  $0.59\%$ / 521s
Ours |  0 / 16s |  0 / 15s |  0 / 15s

##  7 Conclusion

In this work, we used polynomial zonotopes—a set representation originally
developed to analyze nonlinear systems—in order to examine four different
flavors of linear systems. Since polynomial zonotopes preserve dependencies
among both symbolic variables and time, the method computes tight non-convex
enclosures of reachable sets. Comparisons on two benchmark systems illustrate
that our reachability algorithm yields much tighter enclosures than other
state-of-the-art methods. For linear time invariant systems, our method can
represent the reachable set for the entire time horizon with a single
polynomial zonotope. We extend the conference paper in to this journal version
by introducing an efficient and scalable algorithm to optimize over multi-
affine zonotopes, which is a special case of polynomial zonotopes. Through
experiments we demonstrate the efficiency and scalability of our algorithm and
reveal the problem of existing state-of-the-art method. However, it still
faces some problems. For instance, if the way we split the data (vertex cut)
is uneven, the algorithm takes longer to run. Creating an algorithm that
always splits the data evenly (balanced cut) is also difficult. For future
work, we want to explore the possibility of making a random algorithm that can
split the data evenly, or we might try a different method that doesn’t require
checking every possible option for each factor.

## References

  * Althoff (2010) M. Althoff.  _Reachability Analysis and its Application to the Safety Assessment of Autonomous Cars_.  Dissertation, Technische Universität München, 2010.
  * Althoff (2013) M. Althoff.  Reachability analysis of nonlinear systems using conservative polynomialization and non-convex sets.  In _Proc. of the International Conference of Hybrid Systems: Computation and Control_ , pages 173–182, 2013.
  * Althoff (2015) M. Althoff.  An introduction to CORA 2015.  In _Proc. of the Workshop on Applied Verification for Continuous and Hybrid Systems_ , pages 120–151, 2015.
  * Althoff (2019) M. Althoff.  Reachability analysis of large linear systems with uncertain inputs in the Krylov subspace.  _Transactions on Automatic Control_ , 65(2):477–492, 2019.
  * Althoff and et al. (2021) M. Althoff and et al.  Arch-comp21 category report: Continuous and hybrid systems with linear continuous dynamics.  In _Proc. of the International Workshop on Applied Verification of Continuous and Hybrid Systems_ , pages 1–31, 2021.
  * Althoff and Krogh (2011) M. Althoff and B. H. Krogh.  Zonotope bundles for the efficient computation of reachable sets.  In _Proc. of the International Conference on Decision and Control_ , pages 6814–6821, 2011.
  * Althoff and Krogh (2012) M. Althoff and B. H. Krogh.  Avoiding geometric intersection operations in reachability analysis of hybrid systems.  In _Proc. of the International Conference on Hybrid Systems: Computation and Control_ , pages 45–54, 2012.
  * Althoff et al. (2007) M. Althoff, O. Stursberg, and M. Buss.  Reachability analysis of linear systems with uncertain parameters and inputs.  In _Proc. of the International Conference on Decision and Control_ , pages 726–732, 2007.
  * Althoff et al. (2011a) M. Althoff, B. H. Krogh, and O. Stursberg.  _Modeling, Design, and Simulation of Systems with Uncertainties_ , chapter Analyzing Reachability of Linear Dynamic Systems with Parametric Uncertainties, pages 69–94.  Springer, 2011a.
  * Althoff et al. (2011b) M. Althoff, C. Le Guernic, and B. H. Krogh.  Reachable set computation for uncertain time-varying linear systems.  In _Proc. of International Conference on Hybrid Systems: Computation and Control_ , pages 93–102, 2011b.
  * Bak and Duggirala (2017) S. Bak and P. S. Duggirala.  Simulation-equivalent reachability of large linear systems with inputs.  In _Proc. of International Conference on Computer Aided Verification_ , pages 401–420, 2017.
  * Bak and et. al (2022) S. Bak and et. al.  Reachability of koopman linearized systems using random fourier feature observables and polynomial zonotope refinement.  In _Proc. of the International Conference on Computer Aided Verification_ , pages 490–510, 2022.
  * Bak and Tran (2022) S. Bak and H.-D. Tran.  Neural network compression of acas xu early prototype is unsafe: Closed-loop verification through quantized state backreachability.  In _Proc. of the NASA Formal Methods Symposium_ , pages 280–298, 2022.
  * Bak et al. (2019) S. Bak, H.-D. Tran, and T. T. Johnson.  Numerical verification of affine systems with up to a billion dimensions.  In _Proc. of the International Conference on Hybrid Systems: Computation and Control_ , pages 23–32, 2019.
  * Balandin and Kogan (2020) D. V. Balandin and M. M. Kogan.  Control and estimation in linear time-varying systems based on ellipsoidal reachability sets.  _Automation and Remote Control_ , 81(8):1367–1384, 2020.
  * Bogomolov and et al. (2018) S. Bogomolov and et al.  Reach set approximation through decomposition with low-dimensional sets and high-dimensional matrices.  In _Proc. of the International Conference on Hybrid Systems: Computation and Control_ , pages 41–50, 2018.
  * Chutinan and Krogh (2003) A. Chutinan and B. H. Krogh.  Computational techniques for hybrid system verification.  _Transactions on Automatic Control_ , 48(1):64–75, 2003.
  * de Figueiredo and Stolfi (2004) L. H. de Figueiredo and J. Stolfi.  Affine arithmetic: Concepts and applications.  _Numerical Algorithms_ , 37(1-4):147–158, 2004.
  * Del Pia and Di Gregorio (2023) Alberto Del Pia and Silvia Di Gregorio.  On the complexity of binary polynomial optimization over acyclic hypergraphs.  _Algorithmica_ , 85(8):2189–2213, 2023.
  * Duggirala and Bak (2019) P. S. Duggirala and S. Bak.  Aggregation strategies in reachable set computation of hybrid systems.  _Transactions on Embedded Computing Systems_ , 18(5s), 2019.
  * Duggirala and Viswanathan (2016) P. S. Duggirala and M. Viswanathan.  Parsimonious, simulation based verification of linear systems.  In _Proc. of International Conference on Computer Aided Verification_ , pages 477–494, 2016.
  * Forets and Schilling (2022) M. Forets and C. Schilling.  Conservative time discretization: A comparative study.  In _Proc. of the International Conference on Integrated Formal Methods_ , pages 149–167, 2022.
  * Frehse and et al. (2011) G. Frehse and et al.  SpaceEx: Scalable verification of hybrid systems.  In _Proc. of the International Conference on Computer Aided Verification_ , pages 379–395, 2011.
  * Frehse et al. (2013) G. Frehse, R. Kateja, and C. Le Guernic.  Flowpipe approximation and clustering in space-time.  In _Proc. of the International Conference on Hybrid Systems: Computation and Control_ , pages 203–212, 2013.
  * Ghosh and Duggirala (2019) B. Ghosh and P. S. Duggirala.  Robust reachable set: Accounting for uncertainties in linear dynamical systems.  _Transactions on Embedded Computing Systems_ , 18(5s), 2019.
  * Ghosh and Duggirala (2021a) B. Ghosh and P. S. Duggirala.  Reachability of linear uncertain systems: Sampling based approaches.  _arXiv preprint arXiv:2109.07638_ , 2021a.
  * Ghosh and Duggirala (2021b) B. Ghosh and P. S. Duggirala.  Robustness of safety for linear dynamical systems: Symbolic and numerical approaches.  _arXiv preprint arXiv:2109.07632_ , 2021b.
  * Girard (2005) A. Girard.  Reachability of uncertain linear systems using zonotopes.  In _Proc. of the International Conference on Hybrid Systems: Computation and Control_ , pages 291–305, 2005.
  * Girard and Le Guernic (2008) A. Girard and C. Le Guernic.  Zonotope/hyperplane intersection for hybrid systems reachability analysis.  In _Proc. of the International Conference on Hybrid Systems: Computation and Control_ , pages 215–228, 2008.
  * Huang et al. (2023) Yushen Huang, Ertai Luo, Stanley Bak, and Yifan Sun.  On the difficulty of intersection checking with polynomial zonotopes.  In _International Symposium on Automated Technology for Verification and Analysis_ , pages 51–71. Springer, 2023.
  * Jaulin et al. (2006) L. Jaulin, M. Kieffer, and O. Didrit.  _Applied Interval Analysis_.  Springer, 2006.  ISBN 1852332190.
  * Kochdumper (2022) N. Kochdumper.  _Extensions of Polynomial Zonotopes and their Application to Verification of Cyber-Physical Systems_.  Dissertation, Technische Universität München, 2022.
  * Kochdumper and Althoff (2020) N. Kochdumper and M. Althoff.  Constrained polynomial zonotopes.  _arXiv preprint arXiv:2005.08849_ , 2020.
  * Kochdumper and Althoff (2021) N. Kochdumper and M. Althoff.  Sparse polynomial zonotopes: A novel set representation for reachability analysis.  _Transactions on Automatic Control_ , 66(9):4043–4058, 2021.
  * Kochdumper et al. (2020) N. Kochdumper, B. Schürmann, and M. Althoff.  Utilizing dependencies to obtain subsets of reachable sets.  In _Proc. of the International Conference on Hybrid Systems: Computation and Control_ , 2020.
  * Kurzhanski and Varaiya (2000) A. B. Kurzhanski and P. Varaiya.  Ellipsoidal techniques for reachability analysis.  In _Proc. of the International Conference on Hybrid Systems: Computation and Control_ , pages 202–214, 2000.
  * Lal and Prabhakar (2015) R. Lal and P. Prabhakar.  Bounded error flowpipe computation of parameterized linear systems.  In _Proc. of the International Conference on Embedded Software_ , pages 237–246, 2015.
  * Le Guernic and Girard (2010) C. Le Guernic and A. Girard.  Reachability analysis of linear systems using support functions.  _Nonlinear Analysis: Hybrid Systems_ , 4(2):250–262, 2010.
  * Lotov et al. (2004) Alexander V Lotov, Vladimir A Bushenkov, and Georgy K Kamenev.  _Interactive decision maps: Approximation and visualization of Pareto frontier_ , volume 89.  Springer Science & Business Media, 2004.
  * Luo et al. (2023) Ertai Luo, Niklas Kochdumper, and Stanley Bak.  Reachability analysis for linear systems with uncertain parameters using polynomial zonotopes.  In _Proceedings of the 26th ACM International Conference on Hybrid Systems: Computation and Control_ , pages 1–12, 2023.
  * Serry and Reissig (2018) M. Serry and G. Reissig.  Hyper-rectangular over-approximations of reachable sets for linear uncertain systems.  In _Proc. of the International Conference on Decision and Control_ , pages 6275–6282, 2018.
  * Serry and Reissig (2022) M. Serry and G. Reissig.  Overapproximating reachable tubes of linear time-varying systems.  _Transactions on Automatic Control_ , 67(1):443–450, 2022.
  * Silvestre (2022) Daniel Silvestre.  Accurate guaranteed state estimation for uncertain lpvs using constrained convex generators.  In _Proc. of the International Conference on Decision and Control_ , pages 4957–4962, 2022.
  * West et al. (2001) Douglas Brent West et al.  _Introduction to graph theory_ , volume 2.  Prentice hall Upper Saddle River, 2001.
  * Wetzlinger et al. (2020) M. Wetzlinger, N. Kochdumper, and M. Althoff.  Adaptive parameter tuning for reachability analysis of linear systems.  In _Proc. of the International Conference on Decision and Control_ , pages 5145–5152, 2020.
  * Zhang and et al. (2020) Z. Zhang and et al.  Reachable set estimation for uncertain linear systems with known inputs.  _International Journal of Control, Automation and Systems_ , 18(10):2445–2455, 2020.