# Scalable Zonotopic Under-approximation of Backward Reachable Sets for
Uncertain Linear Systems

Liren Yang Necmiye Ozay The authors are with the Dept. of EECS, Univ. of
Michigan, Ann Arbor, MI 48109 {yliren<EMAIL_ADDRESS>
###### Abstract

Zonotopes are widely used for over-approximating forward reachable sets of
uncertain linear systems for verification purposes. In this paper, we use
zonotopes to achieve more scalable algorithms that under-approximate backward
reachable sets of uncertain linear systems for control design. The main
difference is that the backward reachability analysis is a two-player game and
involves Minkowski difference operations, but zonotopes are not closed under
such operations. We under-approximate this Minkowski difference with a
zonotope, which can be obtained by solving a linear optimization problem. We
further develop an efficient zonotope order reduction technique to bound the
complexity of the obtained zonotopic under-approximations. The proposed
approach is evaluated against existing approaches using randomly generated
instances and illustrated with several examples.

##  I Introduction

For autonomous control systems, the control objectives need to be achieved
robustly against system uncertainties. Central to many control synthesis
techniques for uncertain systems is backward reachability analysis. Given an
uncertain control system and a set $X_{0}$ of target states, the backward
reachable set (BRS) consists of the states that can be steered into $X_{0}$ in
finite time, regardless of the system uncertainties. Being able to compute
such sets is important to design controllers with safety or reachability
objectives [4, 19], and is one building block for achieving more complicated
control tasks [6]. Whenever the exact computation is hard, an under-
approximation can be still used to define a conservative control law. A
variety of approaches exists in the literature, including polyhedral
computation [5], interval analysis [17], HJB method [21] and polynomial
optimization [15], just to name a few. For linear dynamics with linear
constraints, polyhedra can be used to represent the BRSs as they are closed
under linear transformation, Minkowski addition and subtraction, and can be
computed leveraging linear optimization tools. However, it is limited to low
dimensional systems (typically, state dimension $\leq 4$) due to an expensive
quantifier elimination step.

One related problem is the forward reachability analysis, where we deal with
uncertain system with no control inputs (e.g., closed-loop systems), and
compute the set of states that can be visited in the future from some initial
state in a given set $X_{0}$. Such forward reachable sets can be computed
offline for verification and online for state prediction [3]. Often times, the
forward reachable sets are overestimated for robustness. For linear systems, a
special polyhedron called zonotope is widely used to represent forward
reachable sets thanks to the favorable complexity of applying linear
transformations (for forward state evolution) and Minkowski additions (to
account for additive uncertainty) to zonotopes (see, e.g., [1, 8]). Algorithms
that compute zonotopic forward reachable sets are more scalable than those
dealing with general polyhedra.

One natural question is: for uncertain linear dynamics, is there a way to
reverse the time so that the efficient zonotopic set computation for forward
reachability analysis can be directly adopted to compute BRSs? Unfortunately,
this is not the case. The main reason is that there lacks a meaningful notion
of two-player game in forward reachability analysis. In the forward case,
there is only one player (i.e., the environment) picking the initial state and
the system uncertainty, whereas in the backward case, there are two players
(i.e., the controller and the environment) picking the control input and the
uncertainty in turn (see Section 4.2 of [20]). Particularly, the existence of
the environment player leads to a Minkowski subtraction step in the sequential
BRS computation, but zonotopes are not closed under Minkowski subtraction [2].
Therefore, combining the idea of time-reversing and efficient computational
tools for forward reachability (e.g., based on zonotopes [18], [10] or
polynomial zonotopes [12]) were explored only for deterministic systems, but
using zonotopes for uncertain systems’ backward reachability, to the best of
our knowledge, is still missing.

In this paper, we use zonotopes to represent and compute BRSs for uncertain
linear systems. The key ingredient is an efficient way to under/over-
approximate the Minkowski difference of two zonotopes by solving convex
optimization problems. While the under approximation allows us to efficiently
compute a subset of the BRS without polyhedral projection, the over-
approximation can be used to quantify how conservative this subset is.
Different from [2], which manipulates a hyperplane-representation, our
approach only deals with the generator-representations of zonotopes, and hence
is more efficient and suitable for sequential computation, but at the cost of
accuracy. The accuracy issue, however, is mitigated by the fact that our
subtrahend zonotope represents the impact of uncertainties and is usually
small compared to the minuend zonotope. Moreover, [2] does not guarantee if
the approximation is an inner one or an outer one. We also leverage the linear
encoding of zonotope-containment problems [23] and derive an alternative
approach for under-approximating the Minkowski difference between zonotopes.
Theoretical analysis and experiments show that our approach scales differently
from this alternative.  In order to upper bound the complexity of each step of
the computation, we further present a way to reduce the order of the obtained
zonotopic BRSs. Zonotope order reduction is extensively studied (e.g., see
[13], [26] and the references therein), but our approach is different: we
search for a lower order zonotope enclosed by the given zonotope, whereas
existing techniques, focusing on forward reachability analysis, all look for
outer approximations. Our approach is evaluated with randomly generated
zonotopes with different dimensions and orders, and its efficacy is
illustrated with several examples.

##  II Preliminaries

Let $G=[g_{1},g_{2},\dots g_{N}]\in\mathbb{R}^{n\times N}$ be a set of
generators, and $c\in\mathbb{R}^{n}$ be a center vector. A zonotope $Z$ with
generator-representation (or G-rep) $(G,c)$ is defined to be the set
$\big{\\{}c+\sum_{i=1}^{N}\theta_{i}g_{i}\mid\theta_{i}\in[-1,1],\ i=1,2,\dots
N\big{\\}}$. With a slight abuse of notation, we will write $Z=(G,c)$. Let
$H\in\mathbb{R}^{L\times n}$ and $h\in\mathbb{R}^{L}$, a polyhedron with
hyperplane-representation (or H-rep) $(H,h)$ is the set
$\\{x\in\mathbb{R}^{n}\mid Hx\leq h\\}$. If polyhedron $X$ is bounded, $X$ is
called a polytope. A set
$V=\\{x_{1},x_{2}\dots,x_{M}\\}\subseteq\mathbb{R}^{n}$ is called a vertex-
representation (or V-rep) of a polytope $X$ if $X$ is the convex hull of $V$,
i.e.,
$X=\text{cvxh}(V):=\big{\\{}\sum_{j=1}^{M}\lambda_{j}x_{j}\mid\sum_{j=1}^{M}\lambda_{j}=1,\lambda_{j}\in[0,1],j=1,2,\dots,M\big{\\}}$,
where cvxh denotes the convex hull. Let $A\in\mathbb{R}^{L\times n}$ and
$X\subseteq\mathbb{R}^{n}$ be a set, $AX$ denotes the set $\\{Ax\mid x\in
X\\}$.

Let $X,Y\subseteq\mathbb{R}^{n}$ be two sets, the Minkowski sum of $X$ and
$Y$, denoted by $X\oplus Y$, is the set $\\{x+y\mid x\in X,y\in Y\\}$.
Whenever $X=\\{x\\}$ is a singleton set, we will write $x+Y$ for $X\oplus Y$.
The Minkowski difference of $X$ and $Y$, denoted by $X\ominus Y$, is defined
to be $\\{z\in\mathbb{R}^{n}\mid z+Y\subseteq X\\}$. For the Minkowski
arithmetics, we assume that the operations are done in order from left to
right, except as specified otherwise by parentheses. The following lemmas will
be useful.

###### Lemma 1.

Let $X,Y,Z\subseteq\mathbb{R}^{n}$.

  * i)

[[16], Proposition 3.1, [25], Lemma 4] $X\ominus Y\oplus Z\subseteq X\oplus
Z\ominus Y$, particularly, $X\ominus Y\oplus Y\subseteq X\subseteq X\oplus
Y\ominus Y$.

  * ii)

[From [22]] If $X$, $Y$ and $Z$ are convex, compact and nonempty, then
$X\oplus Z=Y\oplus Z$ implies that $X=Y$.

  * iii)

If $X$, $Y$ are convex, compact and nonempty, then $X\oplus Y\ominus Y=X$.

###### Proof.

To prove iii), note that, by i), $X\oplus Y\ominus Y\oplus Y=X\oplus Y$. Then
applying item ii) yields $X\oplus Y\ominus Y=X$. ∎

###### Lemma 2.

[From [8]] Let $Z=(G,c)\subseteq\mathbb{R}^{n}$ be a zonotope.

  * i)

$Z=\bigoplus_{i=1}^{N}Z_{i}$ where $Z_{i}=\big{(}\\{g_{i}\\},c_{i}\big{)}$
s.t. $\sum_{i=1}^{N}c_{i}=c$.

  * ii)

Let $A\in\mathbb{R}^{L\times n}$, $AZ=\big{(}AG,Ac)$.

  * iii)

Let $Z^{\prime}=(G^{\prime},c^{\prime})$, $Z\oplus
Z^{\prime}=([G,G^{\prime}],c+c^{\prime})$.

##  III Backward Reachable Sets

Consider a discrete-time system in the following form:

$\displaystyle x_{t+1}=Ax_{t}+Bu_{t}+Ew_{t}+K,$ (1)

where $x\in\mathbb{R}^{n_{x}}$ is the state, $u\in U\in\mathbb{R}^{n_{u}}$ is
the control input and $w\in W\in\mathbb{R}^{n_{w}}$ is the disturbance input.
Given a set $X_{0}$ of target states, we want to compute (or to under-
approximate, if exact computation is hard) the $k$-step backward reachable set
$X_{k}$ of set $X_{0}$, defined recursively as

$\displaystyle X_{k+1}$ $\displaystyle=\\{x\in\mathbb{R}^{n_{x}}\mid\exists
u\in U:\forall w\in W:$ $\displaystyle\ \ \ Ax+Bu+Ew+K\in X_{k}\\},\ \
k=0,1,2\dots$ (2)

Set $X_{k}$ contains the states from where it is possible to reach the target
set $X_{0}$ in exactly $k$ steps. A weaker definition of the $k$-step BRS
would require $X_{0}$ to be reached in no more than $k$ steps, whose formal
definition is similar to Eq. (2) except for an extra “$\cup X_{k}$” at the end
of the formula. Here, we adopt the stronger definition in Eq. (2) for
simplicity because the union operation may lead to non-convex sets. There
exists slightly different notions of reachable sets [14], depending on the
order of the quantifiers. We will focus on under-approximating the set defined
by (2) while our approach applies in general.

Suppose that set $U$, $W$, and $X_{0}$ are polytopes, and that the H-rep of
$U$, $X_{0}$ and the V-rep of $W$ is known, one can compute $X_{k}$ as a
polytope in H-rep, i.e.,

$\displaystyle X_{k+1}$
$\displaystyle=\textbf{Proj}_{x}\big{(}\\{x\in\mathbb{R}^{n_{x}},u\in
U\mid\forall w_{j}\in V_{W}:$ $\displaystyle\ \ \ Ax+Bu+Ew_{j}\in
X_{k}\\}\big{)},\ k=0,1,2\dots,$ (3)

where $\textbf{Proj}_{x}(S)=\\{x\mid\exists u:(x,u)\in S\\}$ is the projection
operation (e.g., see [24], Proposition 1). However, polytope projection is
time-consuming, which limits the use of this approach to low dimensional
systems (typically $n_{x}\leq 4$).

In this paper, we consider under-approximating the BRSs of system (1) instead
under the following assumptions.

  * A1.

The target set is a zonotope (denoted by $Z_{0}$ hereafter), whose G-rep is
known.

  * A2.

The disturbance set $W$ is a polytope, whose H-rep $(H,h)$ and V-rep $V$ are
both known.

  * A3.

Matrix $A\in\mathbb{R}^{n_{x}\times n_{x}}$ is invertible. This assumption is
true whenever Eq. (1) is obtained by time-discretizing an underlying
continuous-time linear dynamics.

Finding under-approximation of BRSs is useful in control problems with
reachability objectives and falsification problems against safety requirements
[7].

##  IV Solution Approach

We explore the use of zonotopes in under-approximating the BRS $X_{k}$. This
is based on i) the modest computational complexity of operations on zonotopes
such as Minkowski addition and affine transformation, and ii) the fact that
(2) can be re-written as follows using Minkowski arithmetic [14]:

$\displaystyle X_{k+1}=\\{x\in\mathbb{R}^{n_{x}}\mid Ax\in X_{k}\ominus
EW\oplus-BU-K\\}.$ (4)

In Eq. (4), if $W=\\{0\\}$ and the term “$\ominus EW$” were not there, then
one could show inductively that, under assumption A1-A3, $X_{k+1}$ is a
zonotope whose G-rep can be easily computed from the G-reps of $X_{k}$ and $U$
after Minknowski addition and linear transformation. Whenever $W$ is not a
singleton set, the key step is to efficiently under and over approximate the
Minkowski difference in Eq. (4) with zonotopes in their G-reps. Whereas the
former leads to an inner approximation of $X_{k+1}$, the latter one can be
used to quantify the conservatism of this inner approximation.

###  IV-A Zonotopic Inner/Outer Approximation of $Z\ominus EW$

Let $Z=(G,c)\subseteq\mathbb{R}^{n_{x}}$ be a zonotope, where
$G=[g_{1},g_{2},\dots,g_{N}]$. We formulate two optimization problems, one
computes a zonotopic outer approximation $\overline{\mathfrak{Z}}(Z,EW)$, and
the other computes a zonotopic inner approximation
$\underline{\mathfrak{Z}}(Z,EW)$, of set $EW$ using $Z$ as a “template”. The
obtained outer/inner approximation are also in G-reps. Particularly, their
generators are scaled versions of $Z$’s generators, i.e., in the form of
$\alpha_{i}g_{i}$ for some $\alpha_{i}\in[0,1]$ (see Definition 1). We then
show that the Minkowski difference $Z\ominus\overline{\mathfrak{Z}}(Z,EW)$ and
$Z\ominus\underline{\mathfrak{Z}}(Z,EW)$ can be done element-wise via
generator subtraction. This leads to an efficient way to inner/outer
approximate $Z\ominus EW$ with zonotopes in G-reps. This technique will become
our key ingredient of BRS under-approximation.

###### Definition 1.

Let $Z=(G,c)$ and $Z^{\prime}=(G^{\prime},c^{\prime})$ be zonotopes.
$Z^{\prime}$ is aligned with $Z$ if $G=[g_{1},g_{2},\dots,g_{N}]$ and
$G^{\prime}=[\alpha_{1}g_{1},\alpha_{2}g_{2},\dots,\alpha_{N}g_{N}]$ for some
$\alpha_{i}\in[0,1]$.

####  IV-A1 Outer approximation of $EW$

Consider the following linear programming problem:

$\displaystyle\begin{array}[]{rl}\min_{\theta,\,\alpha,\,c}&\sum_{i=1}^{N}b_{i}\alpha_{i}\\\
\text{s.t.}&\forall w_{j}\in V:c+\sum_{i=1}^{N}\theta_{ij}g_{i}=Ew_{j}\\\ &\ \
\ |\theta_{ij}|\leq\alpha_{i}\leq 1,\ i=0,1,\dots N\\\ \end{array},$ (min-out)

where $b_{i}>0$ are constants and $\theta$ and $\alpha$ are vectors aggregated
from $\theta_{ij}$ and $\alpha_{i}$ respectively. The V-rep $V$ of the
disturbance set $W$, which is available by Assumption A2, is used to formulate
the above problem. Let $N$ be the number of generators in the template
zonotope $Z$, $n_{x}$ be the dimension of the ambient space, and $M$ be the
number of vertices in $V$. In the optimization problem (min-out),

$\displaystyle\begin{array}[]{rl}\\#\text{variables}&=\mathcal{O}(MN+n_{x}),\\\
\\#\text{constraints}&=\mathcal{O}\big{(}M(N+n_{x})\big{)}.\end{array}$ (10)

###### Proposition 1.

Let $(\theta,\overline{\alpha},\overline{c})$ be the minimizer of the
optimization problem (min-out). Define
$\overline{\mathfrak{Z}}(Z,EW)=\big{(}[\overline{\alpha}_{1}g_{1},\overline{\alpha}_{2}g_{2},\dots,\overline{\alpha}_{N}g_{N}],\overline{c}\big{)}$.
We have $EW\subseteq\overline{\mathfrak{Z}}(Z,EW)$.

###### Proof.

By the conditions in (min-out), for any $i$ and $w_{j}\in V$, there exist
$\theta_{ij}\in[-\overline{\alpha}_{i},\overline{\alpha}_{i}]$ s.t.
$Ew_{j}=\overline{c}+\sum_{i=1}^{N}\theta_{ij}g_{i}$. Equivalently, there
exist $\theta_{ij}^{\prime}\in[-1,1]$ s.t.
$Ew_{j}=\overline{c}+\sum_{i=1}^{N}\theta_{ij}^{\prime}\overline{\alpha}_{i}g_{i}$.
Hence
$EV\subseteq\overline{\mathfrak{Z}}(Z,EW)=\big{(}[\overline{\alpha}_{1}g_{1},\overline{\alpha}_{2}g_{2},\dots,\overline{\alpha}_{N}g_{N}],\overline{c}\big{)}$.
It then follows that
$EW=E\text{cvxh}(V)=\text{cvxh}(EV)\subseteq\overline{\mathfrak{Z}}(Z,EW)$
from the convexity of zonotope $\overline{\mathfrak{Z}}(Z,EW)$. ∎

In general, there does not exist a unique minimal (in the set inclusion sense)
zonotopic outer approximation of $EW$ that aligns with the template zonotope
$Z$. We hence minimize a weighted sum of $\alpha_{i}$’s. The weights $b_{i}>0$
can be used for heuristic design to incorporate prior knowledge of disturbance
set $W$. For example, when $W$ is a hyper-rectangle and
$E\in\mathbb{R}^{n_{x}\times n_{w}}$ is full rank, we use
$b_{i}=\|Tg_{i}\|_{1}-\|Tg_{i}\|_{\infty}$, where $T=(E^{\top}E)^{-1}E$ when
$n_{x}\geq n_{w}$ and $T=E^{\top}(EE^{\top})^{-1}$ otherwise. The idea is to
encourage using generators that closely align with vector $Ee_{p}$, where
$e_{p}$ is the $p^{\rm th}$ natural basis of vector space
$\mathbb{R}^{n_{w}}$. A similar criteria was used for zonotope order reduction
in [8].

####  IV-A2 Inner approximation of $EW$

Consider the following optimization problem:

$\displaystyle\begin{array}[]{rl}\max_{\alpha,\,c}&\textstyle\sum_{i=1}^{N}d_{i}\log(\alpha_{i})\\\
\text{s.t.}&Hc+|HG|\alpha\leq h\\\ &0\leq\alpha\leq 1\end{array},$ (max-in)

where $d_{i}\geq 0$ are constants and $|HG|$ is a matrix obtained by taking
element-wise absolute value of matrix $HG$. The H-rep $(H,h)$ of the
disturbance set $W$, which is available by Assumption A2, is used to formulate
the above problem. Suppose that $H$ has $L$ rows. In (max-in), we have

$\displaystyle\begin{array}[]{rl}\\#\text{
variables}&=\mathcal{O}(N+n_{x}),\\\ \\#\text{
constraints}&=\mathcal{O}(N+L).\end{array}$ (16)

###### Proposition 2.

Let $(\underline{\alpha},\underline{c})$ be the maximizer of optimization
problem (max-in). Define
$\underline{\mathfrak{Z}}(Z,EW)=\big{(}[\underline{\alpha}_{1}g_{1},\underline{\alpha}_{2}g_{2},\dots,\underline{\alpha}_{N}g_{N}],\underline{c}\big{)}$.
We have $\underline{\mathfrak{Z}}(Z,EW)\subseteq EW$.

###### Proof.

We first show that, for $\alpha\geq 0$ and any $c$, $Hc+|HG|\alpha\leq h$ if
and only if

$\displaystyle\forall\theta\in\textstyle\prod_{i=1}^{N}[-\alpha_{i},\alpha_{i}]:H(c+\sum_{i=1}^{N}\theta_{i}g_{i})\leq
h,$ (17)

where $\theta_{i}$ is the $i^{\rm th}$ element of $\theta$. Let $H_{\ell}$ and
$h_{\ell}$ be the $\ell^{\rm th}$ row of $H$ and $h$ respectively. Eq. (17) is
equivalent to

$\displaystyle\forall\ell\in\\{1,2,\dots,L\\}:\begin{array}[]{rl}\max_{\theta}&H_{\ell}(\underline{c}+G\theta)\leq
h_{\ell}\\\ \text{ s.t.
}&\theta\in\textstyle\prod_{i=1}^{N}[-\alpha_{i},\alpha_{i}]\end{array}$ (20)
$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \Updownarrow$
$\displaystyle\forall\ell\in\\{1,2,\dots,L\\}:H_{\ell}c+|H_{\ell}G|\alpha\leq
h_{\ell}.$ (21)

Eq. (21) is equivalent to $Hc+|HG|\alpha\leq h$. Therefore the maximizer
$(\underline{\alpha},\underline{c})$ satisfies Eq. (17), which implies

$\displaystyle\forall\theta^{\prime}\in\textstyle\prod_{i=1}^{N}[-1,1]:H(\underline{c}+\sum_{i=1}^{N}\theta_{i}^{\prime}\underline{\alpha}_{i}g_{i})\leq
h.$ (22)

That is,
$\big{(}[\underline{\alpha}_{1}g_{1},\underline{\alpha}_{2}g_{2},\dots,\underline{\alpha}_{N}g_{N}],\underline{c}\big{)}\subseteq
EW$. ∎

Again, the maximal (in the set inclusion sense) inner approximation does not
exist in general. Here we maximize the volume of a hyper-rectangle in
$\mathbb{R}^{N}$, defined by $d_{i}$ and $\alpha$. In particular, as a
heuristic, we pick $d_{i}=\|g_{i}\|$ for $i=1,2,\dots,N$ throughout the paper.

####  IV-A3 Efficient Minkowski Difference between Aligned Zonotopes

Next, we show that the Minkowski difference amounts to element-wise generator
subtraction when the subtrahend zonotope is aligned with the minuend zonotope.

###### Proposition 3.

Let $Z=(G,c)$ and $Z^{\prime}=(G^{\prime},c^{\prime})$ be zonotopes and
suppose that $Z^{\prime}$ is aligned with $Z$. Then $Z\ominus
Z^{\prime}=\big{(}[(1-\alpha_{1})g_{1},(1-\alpha_{2})g_{2},\dots,(1-\alpha_{N})g_{N}],c-c^{\prime}\big{)}$.

###### Proof.

Let
$\Delta:=\big{(}[(1-\alpha_{1})g_{1},(1-\alpha_{2})g_{2},\dots,(1-\alpha_{N})g_{N}],c-c^{\prime}\big{)}$.
By Lemma 1 iii), $\Delta=\Delta\oplus Z^{\prime}\ominus Z^{\prime}$ as
$\Delta$, $Z^{\prime}$ are convex, compact and nonempty. Also note that
$\Delta\oplus Z^{\prime}=Z$, hence $\Delta=Z\ominus Z^{\prime}$. ∎

We summarize this part by the following proposition.

###### Proposition 4.

Let $Z$ be a zonotope and let $\overline{\mathfrak{Z}}(Z,EW)$,
$\underline{\mathfrak{Z}}(Z,EW)$ be defined by solving (min-out), (max-in)
respectively, then $Z\ominus\overline{\mathfrak{Z}}(Z,EW)\subseteq Z\ominus
EW\subseteq Z\ominus\underline{\mathfrak{Z}}(Z,EW)$. Particularly,
$Z\ominus\overline{\mathfrak{Z}}(Z,EW)$ and
$Z\ominus\underline{\mathfrak{Z}}(Z,EW)$ can be computed efficiently with
generator-wise subtraction.

###### Proof.

It follows from Proposition 1-3 and the fact that both
$\overline{\mathfrak{Z}}(Z,EW)$ and $\underline{\mathfrak{Z}}(Z,EW)$ are
aligned with $Z$. ∎

###  IV-B Approximation of Backward Reachable Sets

We can compute a zonotopic over/under-approximation of the BRS $X_{k}$
recursively as follows:

$\displaystyle\underline{Z}_{0}$ $\displaystyle=\overline{Z}_{0}=Z_{0},$ (23)
$\displaystyle\underline{Z}_{k+1}$
$\displaystyle=A^{-1}\big{(}\underline{Z}_{k}\ominus\overline{\mathfrak{Z}}(\underline{Z}_{k},EW)\oplus-
BU-K\big{)},$ (24) $\displaystyle\overline{Z}_{k+1}$
$\displaystyle=A^{-1}\big{(}\overline{Z}_{k}\ominus\underline{\mathfrak{Z}}(\overline{Z}_{k},EW)\oplus-
BU-K\big{)}.$ (25)

###### Proposition 5.

Let $X_{k}$ be defined by Eq. (2), and $\underline{Z}_{k},\overline{Z}_{k}$ be
defined by Eq. (23)-(25), we have $\underline{Z}_{k}\subseteq
X_{k}\subseteq\overline{Z}_{k}$.

###### Proof.

We prove this by induction. When $k=0$,
$\underline{Z}_{0}=\overline{Z}_{0}=Z_{0}=X_{0}$ by (23). Suppose that
$\underline{Z}_{k}\subseteq X_{k}\subseteq\overline{Z}_{k}$, we have

$\displaystyle\underline{Z}_{k}\ominus\overline{\mathfrak{Z}}(\underline{Z}_{k},EW)$
$\displaystyle\subseteq\underline{Z}_{k}\ominus EW$
$\displaystyle(\text{Proposition }\ref{prop:MinkMinusOverUnder})$
$\displaystyle\subseteq X_{k}\ominus EW.$
$\displaystyle(\underline{Z}_{k}\subseteq X_{k})$ (26)

Combining Eq. (26), (23) and Eq. (4) yields $\underline{Z}_{k+1}\subseteq
X_{k+1}$. Similarly, one can show $X_{k+1}\subseteq\overline{Z}_{k+1}$. ∎

Eq. (24), (25) only involve Minkowski addition, linear transformation of
zonotopes and Minkowski difference where the subtrahend zonotope is aligned
with the minuend zonotope. The above three operations can be done efficiently
with G-rep manipulations. The time for computing $\underline{Z}_{k}$ grows
modestly with $k$ because the number of $\underline{Z}_{k}$’s generators,
denoted by $N_{k}$, increases linearly with $k$. In fact,
$N_{k+1}=N_{k}+N_{U}$ where $N_{U}$ is the (constant) number of generators of
the zonotopic control input set $U$. In what follows, we introduce an order
reduction technique to upper bound the time complexity of computing
$\underline{Z}_{k}$.

####  IV-B1 Zonotope Order Reduction

The order of an $n$-dimensional zonotope with $N$ generators is defined to be
$N/n$. Zonotope order reduction problem concerns approximating a given
zonotope with another one with lower order. Most of the existing techniques
focus on finding outer approximations because zonotopes are typically used to
overestimate forward reachable sets. Whereas in this paper, we find inner
approximations using the following fact.

###### Proposition 6.

Let $Z=\big{(}G=[g_{1},g_{2},\dots,g_{N}],c\big{)}$ be a zonotope. Define
$G_{1}$ to be the matrix after removing arbitrary two columns $g_{i}$, $g_{j}$
from $G$ and appending $g_{i}+g_{j}$, and define $G_{2}$ to be the matrix
after removing columns $g_{i}$, $g_{j}$ from $G$ and appending $g_{i}-g_{j}$.
Then $Z_{1}=(G_{1},c)\subseteq Z$ and $Z_{2}=(G_{2},c)\subseteq Z$.

###### Proof.

Let $\l(g_{k}):=\\{\theta g_{k}\mid\theta\in[-1,1]\\}$, then
$Z=c+\bigoplus_{k=1}^{N}l(g_{k})$, $Z_{1}=c+\bigoplus_{k\neq
i,j}^{N}l(g_{k})\oplus l(g_{i}+g_{j})$. Since $l(g_{i}+g_{j})=\\{\theta
g_{1}+\theta
g_{2}\mid\theta\in[-1,1]\\}\subseteq\\{\theta_{1}g_{1}+\theta_{2}g_{2}\mid\theta_{1},\theta_{2}\in[-1,1]\\}=l(g_{i})\oplus
l(g_{j})$, $Z_{1}\subseteq Z$. Similarly, $Z_{2}\subseteq Z$. ∎

Note that, in Proposition 6, the number of generators of $Z_{1}$ (or $Z_{2}$)
is fewer than that of $Z$ by one. Our zonotope order reduction procedure will
keep replacing two generators $g_{i},g_{j}$ by their combination (either
$g_{i}+g_{j}$ or $g_{i}-g_{j}$) until the order of the resulting zonotope is
small enough. Particularly, we use the following heuristic to select
$g_{i},g_{j}$:

$\displaystyle(i,j)=\text{arg\,}\text{min}_{1\leq i<j\leq
N}\|g_{i}\|_{2}\|g_{j}-\hat{g}_{i}g_{j}^{\top}\hat{g}_{i}\|_{2},$ (27)

where $\hat{g}_{i}=\tfrac{g_{i}}{\|g_{i}\|_{2}}$. Then we will add
$g_{i}+g_{j}$ if
$\|(g_{i}+g_{j})^{\top}G^{\prime}\|_{2}\geq\|(g_{i}-g_{j})^{\top}G^{\prime}\|_{2}$,
and add $g_{i}-g_{j}$ otherwise, where $G^{\prime}$ is the transpose of the
right inverse of the generator matrix after removing columns $g_{i}$, $g_{j}$.
The idea is to combine two generators that are either closely aligned or small
in 2-norm, and the combined generator should be larger and more perpendicular
to the remaining generators.

####  IV-B2 Deriving Reachability Control Law using $\underline{Z}_{k}$

Once zonotopic inner approximations $\underline{Z}_{k}$ of the BRSs are
computed, checking if a state $x$ belongs to $\underline{Z}_{k}$ amounts to
solving a linear program. Moreover, for any state $x\in\underline{Z}_{k+1}$,
we can find a control input $u\in U(x,\underline{Z}_{k})$ that brings $x$ to
$\underline{Z}_{k}$ in one step, where $U(x,\underline{Z}_{k})$ is defined to
be

$\displaystyle\\{u\in U\mid\forall w\in W:Ax+Bu+Ew+K\in\underline{Z}_{k}\\}$
$\displaystyle=$
$\displaystyle\textbf{Proj}_{u}\left\\{(u,\theta)\,\Bigg{|}\,\begin{array}[]{l}Ax+Bu+K=\\\
c^{(k)}+\textstyle\sum_{i=1}^{N_{k}}\theta_{i}g_{i}^{(k)},\\\ u\in
U,\theta_{i}\in[-1,1]\end{array}\right\\},$ (31)

where $\big{(}[g_{1}^{(k)},g_{2}^{(k)},\dots,g_{N_{k}}^{(k)}],c^{(k)}\big{)}$
is the G-rep of
$\underline{Z}_{k}\ominus\overline{\mathfrak{Z}}(\underline{Z}_{k},EW)$, which
can be saved during the computation (see Eq. (24)). We do not need to
explicitly perform the projection step in Eq. (31) as it is sufficient to find
one $u\in U(x,\underline{Z}_{k})$ by solving a linear program. For any initial
state $x_{0}\in\underline{Z}_{k}$, iteratively applying $u_{t}\in
U(x_{t},\underline{Z}_{k-t-1})$ yields a feedback control strategy, which
generates a sequence $u_{0},u_{1},\dots u_{k-1}$ and drives the initial state
$x_{0}$ into the target set $\underline{Z}_{0}=Z_{0}$ in precisely $k$-steps,
regardless of the disturbance inputs.

##  V Evaluation & Discussion

###  V-A Comparisons

We compare our approach for under-approximating $Z\ominus EW$ with two other
methods: one by Althoff [2] and one based on the work by Sadraddini and
Tedrake [23]. Whenever the disturbance set $W$ is a zonotope in its G-rep,
$Z\ominus EW$ can be estimated by [2], but the result is not guaranteed to be
an under-approximation. This approach outperforms the exact computation but is
still expensive due to an H-rep manipulation. Alternatively, using the linear
encoding of zonotope-containment problems [23], we derive the linear program
below that under-approximates $Z\ominus EW$:

$\displaystyle\begin{array}[]{rl}\max_{\Gamma,\gamma,\alpha,c}&\textstyle\sum_{i=1}^{N}\alpha_{i}\\\
\text{s.t.}&[G_{Z}\text{diag}(\alpha),EG_{W}]=G_{Z}\Gamma\\\
&c_{Z}-(c+Ec_{W})=G_{Z}\gamma\\\ &\|[\Gamma,\gamma]\|_{\infty}\leq 1,\
0\leq\alpha\leq 1\end{array},$ (36)

where $(G_{W},c_{W})$ and $(G_{Z},c_{Z})$ are the G-reps of $W$ and $Z$
respectively. Similar to our approach, the solution of (36) also gives a
zonotopic under-approximation $(G_{Z}\text{diag}(\alpha),c)$ of $Z\ominus EW$
that aligns with the template $Z$. The linear program (36) scales differently
from (min-out), which dominates the time of computing BRSs. Let $N_{W}$ and
$N$ be the number of generators of $W$ and $Z$ respectively. For (36),

$\displaystyle\begin{array}[]{rl}\\#\text{variables}&=\mathcal{O}\big{(}N(N+N_{W})+n_{x}\big{)},\\\
\\#\text{constraints}&=\mathcal{O}(N+n_{x}).\end{array}$ (39)

The size of (36) is independent of the number of $W$’s vertices and grows with
$N_{W}$, the number of generators of $W$. Thus (36) is more advantageous than
(min-out) whenever $W$ is a high dimensional zonotope with a small order.  On
the other hand, the number of variables in (min-out) is linear in $N$, whereas
that in (36) is quadratic in $N$.

We randomly generate about 2000 test cases, each case consists of a zonotope
$Z\subseteq\mathbb{R}^{n_{x}}$, a hyper-rectangular
$W\subseteq\mathbb{R}^{n_{x}}$ and a square matrix
$E\in\mathbb{R}^{{n_{x}}\times{n_{x}}}$. The Minkowski difference $Z\ominus
EW$ is estimated using the three different methods. Fig. 1 shows the
computation time w.r.t. the dimension and the order of zonotope $Z$. Each dot
represents the time for a specific case, and the surface is plotted with
averaged values. All the experiments are run on a 1.80 GHz laptop with 16 GB
RAM. The computation time of Althoff’s approach grows fast w.r.t. the order
and the dimension of $Z$ (in fact, we could not finish running any one of the
higher-order cases after hours). Our approach scales better with the order of
$Z$, but still grows relatively fast with the dimension $n_{x}$ because the
number of $W$’s vertices grows exponentially with $n_{x}$ since we choose $W$
to be hyper-rectangles in this example. Somewhat surprisingly, the computation
time of Sadraddini’s approach grows very slowly w.r.t. the order and the
dimension of $Z$. This is consistent with the big-O analysis: in the largest
test case, $n_{x}=10$ and $N=100$, but $W$ has about $10^{3}$ vertices
$(M=1000)$. Hence (min-out) has approximately ten times more variables than
(36).

Figure 1: Upper: computation time for estimating $Z\ominus EW$. Lower: volume
ratio distribution.

Another metric is the size of the obtained estimation. The volumes of the
obtained zonotopic estimations are comparable. Define
$r_{1}=\left(\tfrac{V_{\rm Althoff}}{V_{\rm min-out}}\right)^{1/n_{x}}$ and
$r_{1}=\left(\tfrac{V_{\rm Sadraddini}}{V_{\rm min-out}}\right)^{1/n_{x}}$,
the statistics of $r_{1}$, $r_{2}$ are given in the table below.

| mean | std. | min | max |  confidence of $[0.95,1.05]$
---|---|---|---|---|---
$r_{1}$ | $1.0017$ | $0.0577$ | $0.9900$ | $1.3856$ | $98.83\%$
$r_{2}$ | $0.9678$ | $0.1891$ | $0.8372$ | $1.7498$ | $95.10\%$

###  V-B Order Reduction

We evaluate our order reduction technique with 29000 randomly generated
zonotopes with different dimensions and orders. The approach introduced in
Section IV-B1 is used to reduce the order of each testing zonotope by one. As
shown in Fig. 2 (upper), the computation time grows modestly with the
zonotope’s dimension and order.

Figure 2: Upper: averaged computation time for reducing a zonotope’s order by
one. Lower: volume ratio between the reduced-order zonotope and the one before
reduction.

The quality of the reduced-order zonotope is measured by the ratio between its
volume and that of the original zonotope before reduction, defined in Fig. 2
(lower). We are able to run this evaluation for lower-dimensional cases
because computing the exact volume of a zonotope is difficult for high-
dimensional case due to the combinatorial complexity [9]. In Fig. 2 (lower),
the volume ratio increases with the the original zonotope’s order because
higher order means more freedom in selecting the generators to combine. In the
presented cases, the ratio is close to one if the original zonotope’s order is
greater than three.

##  VI Case Studies

Figure 3: Backward reachable set computation for lateral dynamics. Left:
computation time. Right: set volume.

Figure 4: Backward reachable set computation for longitudinal dynamics. Left:
computation time. Right: set volume.

###  VI-A Aircraft Position Control.

With an aircraft position control system, we illustrate the overall BRS
computation approach that combines Minkowski difference and order reduction to
implement the iterations in Eq. (23)-(25). The linearized 6D lateral dynamics
and the 6D longitudinal dynamics of the aircraft are in the form of Eq. (1),
whose $A$, $B$ matrices are given in Eq. (40). For both systems, $E_{\rm
lat}=E_{\rm long}=I$. The states of the lateral and longitudinal dynamics are
$x_{\rm lat}=[v,p,r,\phi,\psi,y]^{\top}$ and $x_{\rm
long}=[u,w,q,\theta,x,h]^{\top}$ respectively, and control inputs are $u_{\rm
lat}=[\delta_{a},\delta_{r}]^{\top}$ and $u_{\rm
long}=[\delta_{e},\delta_{t}]^{\top}$ respectively (see TABLE I and Fig. 5).
We assume that the disturbance sets are hyper-boxes and their G-rep are
$W_{\rm lat}=(\text{diag}([0.037,0.00166,0.0078,0.00124,0.00107,$
$0.07229]),0)$ and $W_{\rm long}=(\text{diag}([0.3025,0.4025,0.01213,$
$0.006750,1.373,1.331]),0)$.

$\displaystyle A_{\rm lat}=$
$\displaystyle\left[\begin{smallmatrix}1.004&0.1408&0.3095&-0.3112&0&0\\\
0.03015&1.177&0.6016&-0.6029&0&0\\\ -0.02448&-0.1877&0.3803&0.5642&0&0\\\
-0.01057&-0.09588&-0.3343&1.277&0&0\\\
0.0003943&0.0095901&-0.005341&-0.007447&1&0\\\
-0.2579&-23.32&-51.03&61.35&-37.86&1\\\ \end{smallmatrix}\right],$
$\displaystyle A_{\rm long}=$
$\displaystyle\left[\begin{smallmatrix}0.9911&-0.04858&-0.01709&-0.4883&0&0\\\
0.0005870&0.9968&0.5168&-0.0001398&0&0\\\
0.0002070&-0.001123&0.9936&-5.092\times 10^{-5}&0&0\\\
1.907&-1.032&0.01832&1&0&0\\\ -0.04601&0.001125&0.0002638&0.01130&1&0\\\
-5.095\times 10^{-5}&-0.1874&-0.01185&4.004&0&1\\\ \end{smallmatrix}\right],$
$\displaystyle B_{\rm lat}=$
$\displaystyle\left[\begin{smallmatrix}-0.1189&0.007812\\\ -0.1217&0.2643\\\
0.01773&-0.2219\\\ -0.02882&-0.09982\\\ -0.0005607&0.002437\\\
0.1120&-0.5785\\\ \end{smallmatrix}\right],$ $\displaystyle B_{\rm long}=$
$\displaystyle\left[\begin{smallmatrix}1.504&7.349\times 10^{-5}\\\
-0.04645&-3.421\times 10^{-6}\\\ -0.009812&-1.488\times 10^{-6}\\\
-9.080\times 10^{-5}&-1.371\times 10^{-8}\\\ -0.03479&-1.700\times 10^{-6}\\\
0.004171&2.913\times 10^{-7}\\\ \end{smallmatrix}\right].$ (40)

Figure 5: Illustration of the states and control inputs. TABLE I: Variables in the aircraft model variable | physical meaning | range | unit
---|---|---|---
$v$ | velocity | $[-1,1]$ | m/s
$p$ | roll angular rate | $[-1,1]$ | rad/s
$r$ | yaw angular rate | $[-1,1]$ | rad/s
$\phi$ | roll angle | $[-\pi/5,\pi/5]$ | rad
$\psi$ | yaw angle | $[-\pi/5,\pi/5]$ | rad
$y$ | lateral deviation | $[-2,2]$ | m
$u$ | velocity | $[40,60]$ | m/s
$w$ | velocity | $[0,10]$ | rad/s
$q$ | pitch angular rate | $[-0.1,0.1]$ | rad/s
$\theta$ | pitch angle | $[-\pi,\pi]$ | rad
$x$ | horizontal displacement | $[0,800]$ | rad
$h$ | altitude | $[260,390]$ | m
$\delta_{a}$ | aileron deflection | $[-\pi,\pi]$ | m
$\delta_{r}$ | rudder deflection | $[-\pi,\pi]$ | m
$\delta_{e}$ | elevator deflection | $[-0.262,0.524]$ | m
$\delta_{t}$ | throttle control | $[0,10^{4}]$ | m

For both the lateral and longitudinal dynamics, we can efficiently compute
their $k$-step BRSs using the proposed approach for reasonably large horizons
$k$, whereas the computation gets stuck at $k=3$ using the exact Minkowski
difference provided by MPT3 [11], or the approximation function implemented in
CORA. Fig. 4, 4 show the results for the lateral dynamics and the longitudinal
dynamics, respectively. In each figure, the left (right, resp.) plot shows the
cpu time for computing $\underline{Z}_{k}$ (the size of $\underline{Z}_{k}$,
resp.) versus $k$, the number of backward expansion steps. The red curves are
for our approach and the blue ones for the approach using Sadraddini’s
zonotope containment encoding.

####  VI-A1 Cpu time plots

The solid (dotted, resp.) lines correspond to the computation time with
(without, resp.) zonotope order reduction. Using the order reduction technique
(actived at $k=39$), our approach and Sadraddini’s approach give comparable
results. Without order reduction, the computation time of Sadraddini’s
approach (dotted blue) grows faster w.r.t. $k$ than ours (dotted red). This is
consistent with the big-O analysis because in our approach, the time-dominant
Minkwoski difference step amounts to solving a linear program whose number of
variables is linear in $k$ (proportional to the $\underline{Z}_{k}$’s order),
whereas the number of variables is quadratic in $k$ in Sadraddini’s
formulation. Although our approach scales well even without order reduction,
order reduction is still important in efficiently storing the zonotopic BRSs
and deriving the control law.

####  VI-A2 Volume plots

The solid lines correspond to the results with order reduction “in the loop”
(i.e., in the $k^{\rm th}$ step, $\underline{Z}_{k}$ is reduced to a certain
order before $\underline{Z}_{k+1}$ is computed). Whereas the dotted lines
correspond to the results with order reduction after all $\underline{Z}_{k}$’s
are computed (Ideally, we would like to compute $\underline{Z}_{k}$’s volume
without order reduction at all, but this is impossible with the off-the-shelf
volume computation tools in CORA because the complexity is combinatorial in
$\underline{Z}_{k}$’s order). The two approaches give comparable results with
order reduction. Moreover, since the dotted red line and the solid red line
are close to each other, this indicates that the “wrapping effect” due to the
order reduction in-the-loop is relatively small.

###  VI-B Double Integrator with Uncontrollable Subspace

With a 10D system, we show the effectiveness of the reachability controller
derived from the zonotopic BRSs as described in Section IV-B2. The system
consists of a double-integrator dynamics in the 3D space and a 4D
uncontrollable subspace (the uncontrollable part affects the controllable
part). The continuous-time dynamics is

$\displaystyle\dot{x}_{1}$ $\displaystyle=x_{2}+x_{7}+x_{10}+w_{1},\ \
\dot{x}_{2}=u_{1}+w_{2},$ $\displaystyle\dot{x}_{3}$
$\displaystyle=x_{4}-x_{8}+w_{3},\ \ \dot{x}_{4}=u_{2}+w_{4},$
$\displaystyle\dot{x}_{5}$ $\displaystyle=x_{6}+x_{9}+w_{5},\ \
\dot{x}_{6}=u_{3}+w_{6},$ (41) $\displaystyle\dot{x}_{7}$
$\displaystyle=-0.01x_{7}+x_{8}+w_{7},\ \ \dot{x}_{8}=-x_{8}-0.01x_{7}+w_{8},$
$\displaystyle\dot{x}_{9}$ $\displaystyle=-10^{-4}x_{7}+2x_{10}+w_{9},$
$\displaystyle\dot{x}_{10}$ $\displaystyle=-2x_{9}-10^{-4}x_{10}+w_{10}.$

We discretize the above dynamics with a sampling period $\Delta t=0.5$s, and
define the disturbance set $W$ so that $w_{\\{1,3,5\\}}\in[-0.12,0.12]$,
$w_{\\{2,4,6\\}}\in[-0.2,0.2]$, $w_{\\{7,8,9,10\\}}\in[-0.1,0.1]$, and the
control set $U=[-0.5,0.5]^{3}$. Starting from a randomly picked initial
condition in $\underline{Z}_{50}$, our goal is to reach a final state for
which $x_{i}\in[9.5,10.5]$ for $i\in\\{1,3,5\\}$ and $x_{i}\in[-0.5,0.5]$ for
the remaining $i$’s. We defined a controller as described in Section IV-B2.
Fig. 6 shows a closed-loop trajectory under random disturbances. The small
target set is reached despite the oscillating uncontrollable dynamics.

Figure 6: A closed-loop trajectory for the double-integrator dynamics. The red
box is the target set.

##  VII Conclusion

In this paper, we develop an approach that under-approximates the backward
reachable sets for uncertain linear systems using zonotopes. The main
technical ingredients are i) under-approximating the Minkowski difference
between two zonotopes and ii) an order reduction technique tailored to
enclosed zonotopes. These developments were evaluated with randomly generated
instances and two case studies. Experiments show that our method is more
scalable than the off-the-shelf tools (MPT3, CORA) and scales differently from
the approach based on Sadraddini’s zonotope-inclusion technique. In our
method, the dominant Minkowski subtraction step requires solving a linear
program whose size is linear in the zonotope’s order, while that dependency is
quadratic in Sadraddini’s approach. We will investigate extending our approach
to nonlinear systems in the future.

Acknowledgments: The authors would like to thank Yuhao Zhang from the
University of Wisconsin-Madison and Sara Shoouri and Jiahong Xu from the
University of Michigan for sharing the aircraft model.

## References

  * [1] M. Althoff.  An introduction to CORA 2015.  In 1st and 2nd Intl. Workshop on ARCH, pages 120–151, 2015.
  * [2] M. Althoff.  On computing the minkowski difference of zonotopes.  arXiv preprint arXiv:1512.02794, 2015.
  * [3] M. Althoff, G. Frehse, and A. Girard.  Set propagation techniques for reachability analysis.  Annual Review of Control, Robotics, and Autonomous Systems, 4:369–395, 2021.
  * [4] D. Bertsekas.  Infinite time reachability of state-space regions by using feedback control.  IEEE TAC, 17(5):604–613, 1972.
  * [5] F. Blanchini and S. Miani.  Set-theoretic methods in control.  Springer, 2008.
  * [6] M. Chen, Q. Tam, S. C. Livingston, and M. Pavone.  Signal temporal logic meets reachability: Connections and applications.  In Intl. WAFR, pages 581–601. Springer, 2018.
  * [7] G. Chou, Y. E. Sahin, L. Yang, K. J. Rutledge, P. Nilsson, and N. Ozay.  Using control synthesis to generate corner cases: A case study on autonomous driving.  IEEE TCAD, 37(11):2906–2917, 2018.
  * [8] A. Girard.  Reachability of uncertain linear systems using zonotopes.  In Intl. Workshop on HSCC, pages 291–305. Springer, 2005.
  * [9] E. Gover and N. Krikorian.  Determinants and the volumes of parallelotopes and zonotopes.  Linear Algebra Its Appl., 433(1):28–40, 2010.
  * [10] D. Han, A. Rizaldi, A. El-Guindy, and M. Althoff.  On enlarging backward reachable sets via zonotopic set membership.  In ISIC, pages 1–8. IEEE, 2016.
  * [11] M. Herceg, M. Kvasnica, C. Jones, and M. Morari.  Multi-Parametric Toolbox 3.0.  In Proc. of ECC, pages 502–510, Zürich, Switzerland, July 17–19 2013.  http://control.ee.ethz.ch/~mpt.
  * [12] N. Kochdumper and M. Althoff.  Computing non-convex inner-approximations of reachable sets for nonlinear continuous systems.  In 59th CDC, pages 2130–2137. IEEE, 2020.
  * [13] A.-K. Kopetzki, B. Schürmann, and M. Althoff.  Methods for order reduction of zonotopes.  In 56th CDC, pages 5626–5633. IEEE, 2017.
  * [14] A. A. Kurzhanskiy and P. Varaiya.  Reach set computation and control synthesis for discrete-time dynamical systems with disturbances.  Automatica, 47(7):1414–1426, 2011.
  * [15] J. B. Lasserre.  Tractable approximations of sets defined with quantifiers.  Mathematical Programming, 151(2):507–527, 2015.
  * [16] Y. Li.  Robustly complete temporal logic control synthesis for nonlinear systems.  PhD thesis, University of Waterloo, 2019.
  * [17] Y. Li and J. Liu.  Invariance control synthesis for switched nonlinear systems: An interval analysis approach.  IEEE TAC, 63(7):2206–2211, 2017.
  * [18] L. Liebenwein, W. Schwarting, C.-I. Vasile, J. DeCastro, J. Alonso-Mora, S. Karaman, and D. Rus.  Compositional and contract-based verification for autonomous driving on road networks.  In Robotics Research, pages 163–181. Springer, 2020.
  * [19] J. Lygeros, C. Tomlin, and S. Sastry.  Controllers for reachability specifications for hybrid systems.  Automatica, 35(3):349–370, 1999.
  * [20] I. M. Mitchell.  Comparing forward and backward reachability as tools for safety analysis.  In Intl. Workshop on HSCC, pages 428–443. Springer, 2007.
  * [21] I. M. Mitchell, A. M. Bayen, and C. J. Tomlin.  A time-dependent hamilton-jacobi formulation of reachable sets for continuous dynamic games.  IEEE TAC, 50(7):947–957, 2005.
  * [22] L. Montejano.  Some results about minkowski addition and difference.  Mathematika, 43:265–273, 1996.
  * [23] S. Sadraddini and R. Tedrake.  Linear encodings for polytope containment problems.  In 58th CDC, pages 4367–4372. IEEE, 2019.
  * [24] S. W. Smith, P. Nilsson, and N. Ozay.  Interdependence quantification for compositional control synthesis with an application in vehicle safety systems.  In 55th CDC, pages 5700–5707. IEEE, 2016.
  * [25] L. Yang and N. Ozay.  Efficient safety control synthesis with imperfect state information.  In 59th CDC, pages 874–880. IEEE, 2020.
  * [26] X. Yang and J. K. Scott.  A comparison of zonotope order reduction techniques.  Automatica, 95:378–384, 2018.