# Tractable higher-order under-approximating AE extensions for non-linear
systems

Eric Goubault  Sylvie Putot  LIX, CNRS, Ecole Polytechnique, Institut
Polytechnique de Paris, France (email<EMAIL_ADDRESS>
###### Abstract

We consider the problem of under and over-approximating the image of general
vector-valued functions over bounded sets, and apply the proposed solution to
the estimation of reachable sets of uncertain non-linear discrete-time
dynamical systems. Such a combination of under and over-approximations is very
valuable for the verification of properties of embedded and cyber-physical
controlled systems. Over-approximations prove properties correct, while under-
approximations can be used for falsification. Coupled, they provide a measure
of the conservatism of the analysis. This work introduces a general framework
relying on computations of robust ranges of vector-valued functions. This
framework allows us to extend for under-approximation many precision
refinements that are classically used for over-approximations, such as affine
approximations, Taylor models, quadrature formulae and preconditioning
methods. We end by evaluating the efficiency and precision of our approach,
focusing on the application to the analysis of discrete-time dynamical systems
with inputs and disturbances, on different examples from the literature.

###### keywords:

Uncertain systems, Computer-aided control design

##  1 Introduction

Guaranteed state estimation and reachability analysis are central to many
problems in control, such as robust and optimal control of dynamical systems,
set invariance, safety verification, or control synthesis. This ultimately
relies on computing ranges of functions over a domain, that we have to
approximate since this is an intractable problem.

Much of the existing work focuses on over-approximations of images of
functions, or of reachable sets, generally based on convex set representations
(in particular intervals, ellipsoids, polyhedra). We are interested here in
the much less studied problem of computing under-approximations, that is, sets
of states guaranteed to be reached. Combining over and under approximations is
fundamental for the validation of control systems.When the over-approximation
is not sufficient to prove a property, an under-approximation is helpful to
state the quality of the over-approximation. Additionally, when an under-
approximation of the reachable set intersects the set of error states, it
provides a proof of falsification of the property.

For general controlled systems, the reachability properties will depend on the
initial conditions of the system, but also on the sensitivity of the system to
some control inputs and external disturbances, as reflected by the notions of
minimal and maximal reachability Mitchell (2007). We generalize these notions
here to robust reachability, when both control inputs and adversarial
disturbances are present. The robust image will be the intersection, for all
possible disturbances, of the images of a function or reachable sets of a
system.

#### Contents and contributions

The computation of the reachable set of a dynamical system can be reduced to a
series of images of sets by some vector-valued function. In previous work
Goubault and Putot (2020), we introduced mean-value extensions allowing us to
compute under approximations (also called inner-approximations) of such images
in a very efficient way. In this article, we generalize this approach to
higher-order extensions, and develop quadrature formulas for more precise
under-approximations. We also address many questions that are to be solved for
an accurate and efficient implementation:

  * •

Section 2 recaps the necessary background from previous work. Section 3
generalizes the mean-value extension of Goubault and Putot (2020) for the
under and over-approximation of the robust range of sufficiently smooth real-
valued functions $f:\mathbb{R}^{m}\rightarrow\mathbb{R}$; this generalization
allows us to propose higher-order Taylor extensions for under-approximating
robust ranges of functions. These new extensions, just as the mean-value
extensions, is the basis for under-approximations of elementary vector-valued
functions from $\mathbb{R}^{m}$ to $\mathbb{R}^{n}$, as detailed in Section
2.2, which is instrumental in the reachability analysis of dicrete-time
nonlinear systems proposed later;

  * •

Section 4 proposes a novel approach to subdivisions for mean-value and Taylor
extensions, based on the idea of quadratures (in numerical calculations of
integrals): we show that this improves precision of the computation of under
and over-approximations, while still scaling with the dimension of the system;

  * •

Section 5 applies this work to the approximation of robust reachable sets of
discrete-time dynamical systems; we present results on representative systems
from the literature, demonstrating the tractability and precision of our
approach.

#### Related work

Our approach is related to and partially relies on work on modal intervals and
mean-value extensions, which applications include the computation of under-
approximations of function images Goldsztejn (2012a, b). It is also related to
over-approximations of nonlinear functions and dynamical systems, on which we
rely to compute under-approximations. Many methods for over-approximating
reachable sets for non-linear systems have been developed, among which Taylor
methods Makino and Berz (2003) or polytopes-based methods Guernic and Girard
(2009); Dreossi et al. (2016).

There exist less methods for the harder problem of under-approximating images
of functions or sets of reachable states. Some approaches have been proposed
for linear discrete-time systems Kurzhanski and Varaiya (2000); Girard et al.
(2006); Raković and Fiacchini (2008). Interval-based methods, relying on space
discretization, have been used for under-approximating the image of nonlinear
functions Goldsztejn and Jaulin (2010). They were also used to over and under
approximate solutions of differential systems with uncertain initial
conditions Mézo et al. (2018). Tight approximations for reachable sets of
nonlinear continuous systems can be found via expensive Eulerian methods: the
zero sub-level set of the Lipschitz viscosity solution to a Hamilton-Jacobi
(HJB) partial differential equation gives the (backward) reachable set Chen et
al. (2016). Other approaches, using SoS methods and LMI relaxations have been
proposed for inner approximations, see e.g. Korda et al. (2013). In Xue et al.
(2020), under-approximations for polynomial systems are obtained by solving
semi-definite programs. Taylor models are used on the inverse flow map to
derive under-approximations Chen et al. (2014), but using topological
conditions that are checked with interval constraints solving, which have
difficulties to scale up with dimension. In Xue et al. (2016), the computation
of the under-approximated reachable set is based on a costly analysis of the
boundary of the reachable sets and polytopic approximations. In Kochdumper and
Althoff (2020), some non-convex under-approximations are computed with
polynomial zonotopes, relying on a computation of the outer-approximation of
the reachable set, of an enclosure of the boundary of the reachable set, and a
reduction of the outer-approximation until it is fully included in the region
delimited by the boundary.

#### Notations

For a continuously differentiable vector-valued function
$f:\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$, we note $f_{i}$ its $i$-th
component and $\nabla f=(\nabla_{j}f_{i})_{ij}=(\frac{\partial f_{i}}{\partial
x_{j}})_{1\leq i\leq n,1\leq j\leq m}$ its Jacobian matrix. We note $\langle
x,y\rangle$ the scalar product of vectors $x$ and $y$, and $\lvert x\rvert$
the absolute value extended componentwise.

Intervals are used in many situations to rigorously compute with interval
domains instead of reals, usually leading to over-approximations of function
ranges over boxes. Set valued quantities, whether scalar or vector-valued,
will be noted with bold letters, e.g $x$. We denote
${\mathbb{I}\mathbb{R}}=\\{{\makebox{\boldmath$x$}}=[\underline{x},\overline{x}],\>\underline{x}\in\mathbb{R},\>\overline{x}\in\mathbb{R}\\}$
the set of intervals with real-valued bounds. If $\overline{x}<\underline{x}$,
the interval represents the empty set. For a (possibly vector-valued) interval
${\makebox{\boldmath$x$}}\in{\mathbb{I}\mathbb{R}}^{m}$, we note
$c({{\makebox{\boldmath$x$}}})=(\underline{x}+\overline{x})/2$ its center and
$r({{\makebox{\boldmath$x$}}})=(\overline{x}-\underline{x})/2$ its radius. The
operators over real numbers are lifted in intervals using the same notation.

##  2 Background: AE extensions for computing function images

We recall in this section the notations and results of Goubault and Putot
(2020): we state mean-value over and under-approximating extensions for scalar
and vector-valued functions.

An _over-approximating extension_ , also called _outer-approximating
extension_ , of a function $f:\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$ is a
function
${\makebox{\boldmath$f$}}:\mathcal{P}(\mathbb{R}^{m})\rightarrow\mathcal{P}(\mathbb{R}^{n})$,
such that for all $x$ in $\mathcal{P}(\mathbb{R}^{m})$,
$\mbox{range}(f,{\makebox{\boldmath$x$}})=\\{f(x),x\in{\makebox{\boldmath$x$}}\\}\subseteq{\makebox{\boldmath$f$}}({\makebox{\boldmath$x$}})$.
Dually, under-approximations determine a set of values proved to belong to the
range of the function over some input set. An _under-approximating extension_
, also called _inner-approximating extension_ , of $f$ is a function
${\makebox{\boldmath$f$}}:\mathcal{P}(\mathbb{R}^{m})\rightarrow\mathcal{P}(\mathbb{R}^{n})$,
such that for all $x$ in $\mathcal{P}(\mathbb{R}^{m})$,
${\makebox{\boldmath$f$}}({\makebox{\boldmath$x$}})\subseteq\mbox{range}(f,{\makebox{\boldmath$x$}})$.
Under- and over-approximations can be interpreted as quantified propositions:
$\operatorname{range}(f,{\makebox{\boldmath$x$}})\subseteq{\makebox{\boldmath$z$}}$
can be written $\forall x\in{\makebox{\boldmath$x$}},\,\exists
z\in{\makebox{\boldmath$z$}},\,f(x)=z,$ while
${\makebox{\boldmath$z$}}\,\subseteq\operatorname{range}(f,{\makebox{\boldmath$x$}})$
can be written $\forall z\in\,{\makebox{\boldmath$z$}},\,\exists
x\in{\makebox{\boldmath$x$}},\,f(x)=z.$ Both these propositions are what we
will call _AE propositions_ , for quantified propositions where universal
quantifiers (A) precede existential quantifiers (E).

###  2.1 Mean-value AE extensions for scalar-valued functions

We consider a function $f:\mathbb{R}^{m}\rightarrow\mathbb{R}$. The natural
interval extension consists in replacing real operations by their interval
counterparts in the expression of the function. A generally more accurate
extension relies on a linearization by the mean-value theorem.

####  2.1.1 Mean-value AE extensions

Suppose $f$ is differentiable over the box $x$. The mean-value theorem implies
that

$\forall x^{0}\in{\makebox{\boldmath$x$}},\,\forall
x\in{\makebox{\boldmath$x$}},\exists\xi\in{\makebox{\boldmath$x$}},\,f(x)=f(x^{0})+\langle\nabla
f(\xi),x-x^{0}\rangle.$

If we can bound the range of the gradient of $f$ over $x$, by
${\makebox{\boldmath$\nabla f$}}({\makebox{\boldmath$x$}})$, then we can
derive an interval enclosure, called the mean-value extension. Let us choose
$x^{0}$ to be the center $c({{\makebox{\boldmath$x$}}})$ of $x$ and recall we
note $r({{\makebox{\boldmath$x$}}})=(\overline{x}-\underline{x})/2$ its
radius.

######  Theorem 1 (Thm. 1, Goubault and Putot (2020))

Let $f$ be a continuously differentiable function from $\mathbb{R}^{m}$ to
$\mathbb{R}$ and ${\makebox{\boldmath$x$}}\in{\mathbb{I}\mathbb{R}}^{m}$. Let
${\makebox{\boldmath$f$}}^{0}=[\underline{f^{0}},\overline{f^{0}}]$ include
$f(c({{\makebox{\boldmath$x$}}}))$ and $\nabla$ a vector of intervals
${\makebox{\boldmath$\nabla_{i}$}}=[\underline{\nabla}_{i},\overline{\nabla}_{i}]$
for $i\in\\{1,\ldots,m\\}$ such that
$\left\\{\left\lvert\nabla_{i}f(c({\makebox{\boldmath$x$}}_{1}),\ldots,c({\makebox{\boldmath$x$}}_{i-1}),x_{i},\ldots,x_{m})\right\rvert,x\in{\makebox{\boldmath$x$}}\right\\}\subseteq{\makebox{\boldmath$\nabla$}}_{i}.$
We have the over- and under-approximating extensions

$\displaystyle\mbox{range}(f,{\makebox{\boldmath$x$}})\subseteq[\underline{f^{0}},\overline{f^{0}}]+\langle\overline{\nabla},r({{\makebox{\boldmath$x$}}})\rangle[-1,1]$
(1)
$\displaystyle[\overline{f^{0}}-\langle\underline{\nabla},r({\makebox{\boldmath$x$}})\rangle,\underline{f^{0}}+\langle\underline{\nabla},r({\makebox{\boldmath$x$}})\rangle]\subseteq\mbox{range}(f,{\makebox{\boldmath$x$}})$
(2)

###### Example 2.1

Let us consider the range of $f$ defined by $f(x)=x^{2}-x$ over
${\makebox{\boldmath$x$}}=[2,3]$. We can compute $f(2.5)=3.75$ and
$\lvert\nabla f([2,3])\rvert\subseteq[3,5]$. Then (1) and (2) yield
$3.75+1.5[-1,1]\subseteq\mbox{range}(f,[2,3])\subseteq 3.75+2.5[-1,1],$ from
which we deduce
$[2.25,5.25]\subseteq\mbox{range}(f,[2,3])\subseteq[1.25,6.25]$.

We refer to extensions (1) and (2) as _AE extensions_ , as they can be
interpreted as _AE propositions_. Note that the wider, lesser quality, are the
over-approximations of $f$ and its derivatives, the tighter, less quality, are
the under-approximations. The under-approximation can even become empty if the
width $\overline{f_{0}}-\underline{f_{0}}$ of the approximation of
$f(c({{\makebox{\boldmath$x$}}}))$ exceeds
$2\langle\overline{\nabla},r({{\makebox{\boldmath$x$}}})\rangle$: in this case
the lower bound of the resulting interval is larger than the upper bound,
which by convention we identify with the empty interval. Note also that when
$0\in{\makebox{\boldmath$\nabla_{i}{\makebox{\boldmath$f$}}$}}$, then
$\underline{\nabla_{i}}=0$ and if this is the case for all $i$, the under-
approximation is empty or reduced to a point. A special attention to the
practical evaluation of these extensions over the region $x$ of interest is
thus crucial, this is the object of Section 4.

####  2.1.2 Mean-value AE extensions of the robust range

Mean-value AE extensions can be generalized to compute ranges that are robust
to disturbances, identified as some input components. Let us partition the
indices of the input space in two subsets $I_{\mathcal{A}}$ and
$I_{\mathcal{E}}$, where $I_{\mathcal{A}}$ defines the indices of the inputs
that correspond to disturbances, and $I_{\mathcal{E}}$ the remaining
dimensions. We decompose the input box $x$ accordingly by
${\makebox{\boldmath$x$}}={\makebox{\boldmath$x$}}_{\mathcal{A}}\times{\makebox{\boldmath$x$}}_{\mathcal{E}}$.
We define the robust range of function $f$ on $x$, robust on
${\makebox{\boldmath$x$}}_{\mathcal{E}}$ with respect to disturbances
${\makebox{\boldmath$x$}}_{\mathcal{A}}$, as
$\mbox{range}(f,{\makebox{\boldmath$x$}},I_{\mathcal{A}},I_{\mathcal{E}})=\\{z\,|\,\forall
w\in{\makebox{\boldmath$x$}}_{\mathcal{A}},\,\exists
u\in{\makebox{\boldmath$x$}}_{\mathcal{E}},\,z=f(w,u)\\}$. Intuitively, $u$
will be control components, $w$ disturbances to which the output range should
be robust.

######  Theorem 2 (Thm. 2, Goubault and Putot (2020))

Let $f$ be continuously differentiable function from $\mathbb{R}^{m}$ to
$\mathbb{R}$ and
${\makebox{\boldmath$x$}}={\makebox{\boldmath$x$}}_{\mathcal{A}}\times{\makebox{\boldmath$x$}}_{\mathcal{E}}\in{\mathbb{I}\mathbb{R}}^{m}$.
Let ${\makebox{\boldmath$f$}}^{0}$, ${\makebox{\boldmath$\nabla$}}_{w}$ and
${\makebox{\boldmath$\nabla$}}_{u}$ be vectors of intervals such that
$f(c({\makebox{\boldmath$x$}}))\subseteq{\makebox{\boldmath$f$}}^{0}$,
$\\{\left|\nabla_{w}f(w,c({\makebox{\boldmath$x$}}_{\mathcal{E}}))\right|\,,\,w\in{\makebox{\boldmath$x$}}_{\mathcal{A}}\\}\subseteq{\makebox{\boldmath$\nabla$}}_{w}$
and
$\\{\left|\nabla_{u}f(w,u)\right|\,,w\in{\makebox{\boldmath$x$}}_{\mathcal{A}},\,u\in{\makebox{\boldmath$x$}}_{\mathcal{E}}\\}\subseteq{\makebox{\boldmath$\nabla$}}_{u}$.
We have:

$\mbox{range}(f,{\makebox{\boldmath$x$}},I_{\mathcal{A}},I_{\mathcal{E}})\subseteq[\underline{f^{0}}-\langle\overline{\nabla}_{u},r({{\makebox{\boldmath$x$}}_{\mathcal{E}}})\rangle+\langle\underline{\nabla}_{w},r({{\makebox{\boldmath$x$}}_{\mathcal{A}}})\rangle,\\\
\overline{f^{0}}+\langle\overline{\nabla}_{u},r({{\makebox{\boldmath$x$}}_{\mathcal{E}}})\rangle-\langle\underline{\nabla}_{w},r({{\makebox{\boldmath$x$}}_{\mathcal{A}}})\rangle]$
(3)
$[\overline{f^{0}}-\langle\underline{\nabla}_{u},r({{\makebox{\boldmath$x$}}_{\mathcal{E}}})\rangle+\langle\overline{\nabla}_{w},r({{\makebox{\boldmath$x$}}_{\mathcal{A}}})\rangle,\underline{f^{0}}+\\\
\langle\underline{\nabla}_{u},r({{\makebox{\boldmath$x$}}_{\mathcal{E}}})\rangle-\langle\overline{\nabla}_{w},r({{\makebox{\boldmath$x$}}_{\mathcal{A}}})\rangle]\subseteq\mbox{range}(f,{\makebox{\boldmath$x$}},I_{\mathcal{A}},I_{\mathcal{E}})$
(4)

We refer to Example 2 of Goubault and Putot (2020) for a sample computation.

###  2.2 AE extensions for vector-valued functions

Following Goubault and Putot (2020), we now detail how full n-dimensional
boxes can be included in the image of vector-valued functions
$f:\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$, for $m\geq n$, using AE
extensions of robust ranges. Theorem 3 and Definition 2.2 will be instrumental
in Algorithm 1 for discrete-time reachability of Section 5.

The mean-value extensions of Theorem 1 or the generalization of Theorem 5 give
us under and over-approximations of projections of the image of the function.
The Cartesian product of the over-approximations of each component provides an
over-approximation of a vector-valued function
$f:\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$. This is however not the case for
under-approximation. Suppose for example that we have $\forall
z_{1}\in{\makebox{\boldmath$z$}}_{1},\exists
x_{1}\in{\makebox{\boldmath$x$}}_{1},\,\exists
x_{2}\in{\makebox{\boldmath$x$}}_{2},\>z_{1}=f_{1}(x)$ and $\forall
z_{2}\in{\makebox{\boldmath$z$}}_{2},\exists
x_{1}\in{\makebox{\boldmath$x$}}_{1},\,\exists
x_{2}\in{\makebox{\boldmath$x$}}_{2},\>z_{2}=f_{2}(x)$. We cannot deduce
directly that for all $\forall z_{1}\in{\makebox{\boldmath$z$}}_{1}$ and
$\forall z_{2}\in{\makebox{\boldmath$z$}}_{2}$ there exists $x_{1}$ and
$x_{2}$ such that $z=f(x)$.

Suppose now that we have: $\forall
z_{1}\in{\makebox{\boldmath$z$}}_{1},\forall
x_{1}\in{\makebox{\boldmath$x$}}_{1},\,\exists
x_{2}\in{\makebox{\boldmath$x$}}_{2},\>z_{1}=f_{1}(x)$ and $\forall
z_{2}\in{\makebox{\boldmath$z$}}_{2},\forall
x_{2}\in{\makebox{\boldmath$x$}}_{2},\,\exists
x_{1}\in{\makebox{\boldmath$x$}}_{1},\>z_{2}=f_{2}(x)$ with continuous
selections $x_{2}$ and $x_{1}$. Then there exists functions
$g_{2}(z_{1},x_{1}):{\makebox{\boldmath$z$}}_{1}\times{\makebox{\boldmath$x$}}_{1}\rightarrow{\makebox{\boldmath$x$}}_{2}$
and
$g_{1}(z_{2},x_{2}):{\makebox{\boldmath$z$}}_{2}\times{\makebox{\boldmath$x$}}_{2}\rightarrow{\makebox{\boldmath$x$}}_{1}$
that are continuous in $x_{1}$ (resp. $x_{2}$), and such that
$\forall(z_{1},z_{2})\in{\makebox{\boldmath$z$}}$,
$\forall(x_{1},x_{2})\in{\makebox{\boldmath$x$}}$,
$z_{1}=f_{1}(x_{1},g_{2}(z_{1},x_{1}))$ and
$z_{2}=f_{2}(g_{1}(z_{2},x_{2}),x_{2})$. Using the Brouwer fixed point theorem
on the continuous map
$g~{}:(x_{1},x_{2})\rightarrow(g_{1}(z_{2},x_{2}),g_{2}(z_{1},x_{1}))$ on the
compact set ${\makebox{\boldmath$x$}}_{1}\times{\makebox{\boldmath$x$}}_{2}$,
then $\forall(z_{1},z_{2})\in{\makebox{\boldmath$z$}}$,
$\exists(x^{z}_{1},x^{z}_{2})\in{\makebox{\boldmath$x$}}$ such that
$(z_{1},z_{2})=f(x^{z}_{1},x^{z}_{2})$.

This result can be generalized to functions
$f~{}:\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$ for any $n$, as stated in
Theorem 3.

######  Theorem 3 (Theorem 3 in Goubault and Putot (2020))

Let $f:\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$ be an elementary function and
$\pi:[1\ldots m]\rightarrow[1\ldots n]$ Let us note, for all $i\in[1\ldots
n]$, $J_{E}^{(z_{i})}=\\{j\in[1\ldots m],\;\pi(j)=i\\}$ and
$J_{A}^{(z_{i})}=\\{j\in[1\ldots m]\\}\setminus J_{E}^{(z_{i})}$. Consider the
$n$ AE-extensions $i\in[1\ldots n]$ built from Theorems 2, 4 or 5 and such
that

$\forall z_{i}\in{\makebox{\boldmath$z$}}_{i},\>(\forall
x_{j}\in{\makebox{\boldmath$x$}}_{j})_{j\in J_{A}^{(z_{i})}},\>(\exists
x_{j}\in{\makebox{\boldmath$x$}}_{j})_{j\in J_{E}^{(z_{i})}},\>z_{i}=f_{i}(x)$
(5)

Then
${\makebox{\boldmath$z$}}={\makebox{\boldmath$z$}}_{1}\times{\makebox{\boldmath$z$}}_{2}\times\ldots\times{\makebox{\boldmath$z$}}_{n}$,
if non-empty, is an under-approximation of the image of $f$: $\forall
z\in{\makebox{\boldmath$z$}},\,\exists x\in{\makebox{\boldmath$x$}},\>z=f(x).$

Theorem 3 gives us directly a computation of an under-approximation of
$\mbox{range}(f,{\makebox{\boldmath$x$}})$ for
$f:\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$. It can also be used to compute an
under-approximation of the robust range
$\mbox{range}(f,{\makebox{\boldmath$x$}},I_{\mathcal{A}},I_{\mathcal{E}})$.
For this, we need to choose $\pi:[1\ldots m]\rightarrow([1\ldots n]\setminus
I_{\mathcal{A}})$, which corresponds to the fact that the disturbance part of
the input components will always be quantified universally. We define below
the result of this process, which will be later used in reachability
algorithms for discrete-time dynamical systems.

###### Definition 2.2

Let $f:\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$ and $\pi:[1\ldots
m]\rightarrow([1\ldots n]\setminus I_{\mathcal{A}})$. We define
${\mathcal{I}}(f,{\makebox{\boldmath$x$}},I_{\mathcal{A}},I_{\mathcal{E}},\pi)$
an under-approximation of
$\mbox{range}(f,{\makebox{\boldmath$x$}},I_{\mathcal{A}},I_{\mathcal{E}})$
obtained using Theorem 3 with function $\pi$, in which the under-approximation
of each component is obtained using Theorem 2 (or Corollary 3.2). We define
${\mathcal{O}}(f,{\makebox{\boldmath$x$}},I_{\mathcal{A}},I_{\mathcal{E}},\pi)$
the over-approximation of the robust range obtained using Theorem 2 component-
wise.

##  3 Generalization to new AE extensions

We now introduce new robust AE extensions for a function
$f:\mathbb{R}^{m}\rightarrow\mathbb{R}$, which are no longer necessarily based
on the mean-value theorem.

###### Theorem 4

Suppose we have an approximation function $g$ for $f$, which is an
elementary111Elementary functions are compositions of +, -, $\times$, /, sine,
cosine, log, exp functions in particular. function in the sense of Goldsztejn
(2012a), satisfying $\forall w\in{\bf x}_{\mathcal{A}},\ \forall u\in{\bf
x}_{\mathcal{E}},\ \exists\xi\in{\bf x},\ f(w,u)=g(w,u,\xi)$. Then any under-
approximation (resp. over-approximation) of the robust range of $g$ with
respect to $x_{\mathcal{A}}$ and $\xi$,
${\mathcal{I}}_{g}\subseteq\mbox{range}(g,{\bf x}\times{\bf
x},I_{\mathcal{A}}\cup\\{m+1,\ldots,2m\\},I_{\mathcal{E}})$ is an under-
approximation (resp. over-approximation) of the robust range of $f$ with
respect to $x_{\mathcal{A}}$, i.e.
${\mathcal{I}}_{g}\subseteq\mbox{range}(f,{\bf
x},I_{\mathcal{A}},I_{\mathcal{E}})$.

For instance, for a continuously $(n+1)$-differentiable $f$, the following
$g$, obtained by a Taylor-Lagrange expansion and noting $x=(w,u)$, is an
approximation function for $f$

$g(x,\xi)=f(x^{0})+\sum\limits_{i=1}^{n}\frac{(x-x^{0})^{i}}{i!}D^{i}f(x^{0})\\\
+D^{n+1}f(\xi)\frac{(x-x^{0})^{n+1}}{(n+1)!}$ (6)

where $D^{\alpha}f$ denotes the higher order partial derivative of $f$. For
$n=0$, $g$ is the mean-value approximation.

We focus on the under-approximation. As $g$ is elementary, by Proposition 10.1
of Goldsztejn (2012a), we have a continuous Skolem function
$(w,u)\rightarrow\xi(w,u)$, i.e. a function such that for all
$x=(w,u)\in\mathbb{R}^{m}$, $f(w,u)=g(w,u,\xi(x))$. Consider an under-
approximation ${\mathcal{I}}_{g}$ of the robust range of $g$ with respect to
$\xi$ and $x_{\mathcal{A}}$. It satisfies $\forall z\in{\mathcal{I}}_{g},\
\forall w\in{\bf x}_{\mathcal{A}},\ \forall\xi\in{\bf x},\ \exists u\in{\bf
x}_{\mathcal{E}},\ z=g(w,u,\xi)$. Let $z\in{\mathcal{I}}_{g}$ and $w\in{\bf
x}_{\mathcal{A}}$ be fixed. As $g$ is elementary, we have a corresponding
continuous Skolem function $\xi\rightarrow u(\xi)$, i.e. a function such that
for all $\xi\in{\bf x}$, $z=g(w,u(\xi),\xi)$. For this $z\in{\mathcal{I}}_{g}$
and $w\in{\bf x}_{\mathcal{A}}$, the continuous map $u\rightarrow u(\xi(w,u))$
defined from ${\bf x}_{\mathcal{E}}$ over itself, has a fixed point
$u^{\infty}$, by Brouwer’s theorem. It is such that
$z=f(w,u^{\infty})=g(w,u^{\infty},\xi(w,u^{\infty}))$. Hence $z$ is in the
robust range of $f$ with respect to $x_{\mathcal{A}}$.

###### Example 3.1

Consider function $f(x)=x^{3}+x^{2}+x+1$ on $[-\frac{1}{4},\frac{1}{4}]$. The
exact range is $[0.796875,1.328125]$. Let us approximate $f$ by a quadratic
function, using an order 2 Taylor-Lagrange expansion. We compute
$f^{(1)}(x)=3x^{2}+2x+1$ and $f^{(2)}(x)=6x+2$. By Theorem 4, the range of $f$
over $[-\frac{1}{4},\frac{1}{4}]$ is under (resp. over) approximated by any
under (resp. over) approximation of the robust range with respect to $\xi$ of

$\begin{array}[]{rcl}g(x,\xi)&=&f(x^{0})+(x-x^{0})f^{(1)}(x^{0})+f^{(2)}(\xi)\frac{(x-x^{0})^{2}}{2}\\\
&=&1+x+x^{2}(3\xi+1)\end{array}$

In the general case, it may still be difficult to compute the under-
approximated robust range of $g$. However, Theorem 5 gives a simple way which
is well suited in particular for quadratic Taylor-based approximations.

###### Theorem 5

Let $g$ be an elementary function $g(w,u,\xi)=\alpha(w,u)+\beta(w,u,\xi)$ over
$x=(w,u)\in{\bf x}\subseteq{\mathbb{I}\mathbb{R}}^{m}$ and $\xi\in{\bf x}$.
Let ${\mathcal{I}}_{\alpha}$ be an under-approximation of the robust range of
$\alpha$ with respect to $w$, i.e. $\mbox{range}(\alpha,{\bf
x},I_{\mathcal{A}},I_{\mathcal{E}})$, and ${\mathcal{O}}_{\beta}$ an over-
approximation of the range of $\beta$, i.e. $\mbox{range}(\beta,{\bf
x}\times{\bf x},\emptyset,\\{1,\ldots,2m\\})$.

The robust range of $g$ with respect to $w\in{\bf x}_{\mathcal{A}}$ and
$\xi\in{\bf x}$, i.e. $\mbox{range}(g,{\bf x}\times{\bf
x},I_{\mathcal{A}}\cup\\{m+1,\ldots,2m\\},I_{\mathcal{E}})$, is under-
approximated by
${\mathcal{I}}_{g}=[\underline{\mathcal{I}}_{\alpha}+\overline{\mathcal{O}}_{\beta},\overline{\mathcal{I}}_{\alpha}+\underline{\mathcal{O}}_{\beta}].$

This is the case of Taylor expansions (6), where $\alpha$ is the degree $n$
polynomial, and $\beta$ the degree $n+1$ remainder.

${\mathcal{I}}_{\alpha}$ under-approximating $\mbox{range}(\alpha,{\bf
x},I_{\mathcal{A}},I_{\mathcal{E}})$ means $\forall
a\in{\mathcal{I}}_{\alpha},\ \forall w\in{\bf x}_{\mathcal{A}},\ \exists
u\in{\bf x}_{\mathcal{E}},\ a=\alpha(w,u)$. As $\alpha$ is elementary, we have
a continuous Skolem function $(w,a)\rightarrow u(w,a)$. Moreover, for all $z$
in ${\mathcal{I}}_{g}$, for all $b\in{\mathcal{O}}_{\beta}$, we have
$z-b\in[\underline{\mathcal{I}}_{\alpha}+\overline{\mathcal{O}}_{\beta}-\overline{\mathcal{O}}_{\beta},\overline{\mathcal{I}}_{\alpha}+\underline{\mathcal{O}}_{\beta}-\underline{\mathcal{O}}_{\beta}]={\mathcal{I}}_{\alpha}$.

For given $z\in{\mathcal{I}}_{g}$, $w\in{\bf x}_{\mathcal{A}}$ and $\xi\in{\bf
x}$, consider the continuous function $b\rightarrow\beta(w,u(w,z-b),\xi)$ over
${\mathcal{O}}_{\beta}$. By Brouwer fixed point theorem we have
$b^{\infty}=\beta(w,u(w,z-b^{\infty}),\xi)$. Therefore, for any
$z\in{\mathcal{I}}_{g}$, $\xi\in{\bf x}$, there exist
$a=z-b^{\infty}\in{\mathcal{I}}_{\alpha}$,
$b=b^{\infty}\in{\mathcal{O}}_{\beta}$ and $x=(w,u(w,a))\in{\bf x}$ such that
$z=a+b$. This implies that ${\mathcal{I}}_{g}$ is a robust under-approximation
of $g(x,\xi)=\alpha(x)+\beta(x,\xi)$ with respect to $w$ and $\xi$.

A direct consequence is a simple order 2 under-approximating Taylor method:

###### Corollary 3.2

Consider $f:\mathbb{R}^{m}\rightarrow\mathbb{R}$ a function in $C^{2}$. Let
${\makebox{\boldmath$f$}}^{0}$, ${\makebox{\boldmath$\nabla$}}_{w}^{0}$ and
${\makebox{\boldmath$\nabla$}}_{u}^{0}$ be such that
$f(x^{0})\subseteq{\makebox{\boldmath$f$}}^{0}$,
$\left|\nabla_{w}f(x^{0})\right|\subseteq{\makebox{\boldmath$\nabla$}}_{w}^{0}$
and
$\left|\nabla_{u}f(x^{0})\right|\subseteq{\makebox{\boldmath$\nabla$}}_{u}^{0}$
with $x^{0}=c({\makebox{\boldmath$x$}})$.

Then
$\mbox{range}(f,{\makebox{\boldmath$x$}},I_{\mathcal{A}},I_{\mathcal{E}})$ is
under-approximated by
$[\underline{\mathcal{I}}_{\alpha}+\overline{\mathcal{O}}_{\beta},\overline{\mathcal{I}}_{\alpha}+\underline{\mathcal{O}}_{\beta}]$
where
${\mathcal{I}}_{\alpha}=[\overline{f^{0}}-\langle\underline{\nabla}_{u}^{0},r({{\makebox{\boldmath$x$}}_{\mathcal{E}}})\rangle+\langle\overline{\nabla}_{w}^{0},r({{\makebox{\boldmath$x$}}_{\mathcal{A}}})\rangle,$
$\underline{f^{0}}+\langle\underline{\nabla}_{u}^{0},r({{\makebox{\boldmath$x$}}_{\mathcal{E}}})\rangle-\langle\overline{\nabla}_{w}^{0},r({{\makebox{\boldmath$x$}}_{\mathcal{A}}})\rangle]$
and ${\mathcal{O}}_{\beta}$ is any over-approximation of
$\\{\frac{1}{2}D^{2}f(x)(r({\makebox{\boldmath$x$}}))^{2},x\in{\makebox{\boldmath$x$}}\\}$.

This is a direct application of Theorems 4 and 5 where $g$ is the 2nd order
Taylor approximant, i.e. given by Equation (6) for $n=1$, combined with
Theorem 2 to compute the robust under-approximation of the order 1
approximation. We use a particular case of Theorem 2: the under-approximation
of an order 1 polynomial is almost trivial, we can compute it exactly if the
computation is performed in real numbers with exact evaluation of $f$ and its
gradient at a point. The expression we give here accounts for computation
errors.

###### Example 3.3

We carry on with Example 3.1. The computation of the under approximation of
the range of $1+x$ over $[-\frac{1}{4},\frac{1}{4}]$ yields
$[\frac{3}{4},\frac{5}{4}]$. We also need an over approximation of the range
of $x^{2}(3\xi+1)$ for $x$ and $\xi$ in $[-\frac{1}{4},\frac{1}{4}]$. Standard
interval computation yields
$[0,\frac{1}{16}][\frac{1}{4},\frac{7}{4}]=[0,\frac{7}{64}]$. Overall, we
deduce $[0.859375,1.25]\subseteq\mbox{range}(f,{\makebox{\boldmath$x$}})$. In
comparison, the mean-value AE extension of Theorem 2 would have given us the
less precise under-approximation $[0.875,1.125]$.

##  4 Refinements: preconditioning and quadrature formulae

The n-dimensional inner boxes that we compute with the techniques of Section
2.2 can sometimes be small or empty, even when the projected inner-
approximations on each component are tight. There are different reasons, for
which we propose solutions in this section.

###  4.1 Preconditioning for computing inner skewed boxes

The first difficulty is when the image of the vector-valued function cannot be
precisely approximated by a centered box.

###### Example 4.1

We consider $f(x)=(2x_{1}^{2}-x_{1}x_{2}-1,x_{1}^{2}+x_{2}^{2}-2)^{\intercal}$
with ${\makebox{\boldmath$x$}}=[0.9,1.1]^{2}$. The under-approximated
projections on the two components (respectively $[-0.38,0.38]$ and
$[-0.38,0.38]$) are close to the over-approximated ranges
($[-0.42,0.42]^{2}$), but we only find empty inner boxes.

This problem can be partly solved, as already described in Goubault and Putot
(2020) by computing a skewed box as under-approximation, that is the image of
a box by a linear map, instead of a box. This can be achieved by combining
preconditioning to the mean-value theorem. Let $C\in\mathbb{R}^{n\times n}$ be
a non-singular matrix. If $z$ is an interval vector such that
${\makebox{\boldmath$z$}}\subseteq\mbox{range}(Cf,{\makebox{\boldmath$x$}})$ ,
we can deduce a skewed box to be in the range of $f$, that is
$\\{C^{-1}z|z\in{\makebox{\boldmath$z$}}\\}\subseteq\mbox{range}(f,{\makebox{\boldmath$x$}})$.
A natural choice for $C$ is the inverse of the center of the interval Jacobian
matrix $C=(c({\makebox{\boldmath$\nabla$}}))^{-1}$.

###### Example 4.2

On Example 4.1, using this preconditioning and $pi:(1\rightarrow
1,2\rightarrow 2)$, we obtain for $f_{2}$ as a function of $f_{1}$ the yellow
under-approximating parallelotope of Figure 1(a). We estimate the image
$\mbox{range}(f,{\makebox{\boldmath$x$}})$ by sampling points in the input
domain. This sampling-based estimation is represented as the dark dots-filled
region. The green parallelotope and box are the over-approximations with and
without preconditioning.

(a) Example 4.2: under- (dotted lines) and over-approximation (plain lines)

$x^{0}$$x^{1}$${\makebox{\boldmath$x$}}^{k}={\makebox{\boldmath$x$}}\setminus{\makebox{\boldmath$x$}}^{k-1},\nabla^{k}=[\lvert\nabla
f\rvert]({\makebox{\boldmath$x$}}^{k})$$x_{1}^{-k}$$x_{1}^{-1}$$x_{1}^{0}$$x_{1}^{1}$$x_{1}^{k}$$x_{2}^{-k}$$x_{2}^{0}$$x_{2}^{k}$${\makebox{\boldmath$x$}}^{1},\nabla^{1}$${\makebox{\boldmath$x$}}^{2},\nabla^{2}$
(b) Partitioning the input domain

(c) Example 4.3: approximations for quadrature and order 2 extensions

Figure 1: Illustrations for Sections 2.2 and 4

###  4.2 Quadrature formulae for the mean-value extension

The mean-value interval extension can yield a rough approximation. This is
especially the case when the variation of the gradient is important over the
input range, the extreme case for under-approximation being when this
variation contains zero: the under-approximation is empty or reduced to a
point. Using simple quadrature formulae partially solves this problem.

Let $f:\mathbb{R}^{m}\rightarrow\mathbb{R}$. We partition each dimension
$j=[1\ldots m]$ of the $m$-dimensional input box
${\makebox{\boldmath$x$}}={\makebox{\boldmath$x$}}_{1}\times\ldots\times{\makebox{\boldmath$x$}}_{m}$
in $2k$ sub-intervals and define, for all $j=[1\ldots m]$, $x_{j}^{-k}\leq
x_{j}^{-(k-1)}\leq\ldots\leq x_{j}^{0}\leq\ldots\leq x_{j}^{k}$, with
$x_{j}^{-k}=\underline{{\makebox{\boldmath$x$}}_{j}}$,
$x_{j}^{0}=\mbox{mid}({{\makebox{\boldmath$x$}}_{j}})$,
$x_{j}^{k}=\overline{{\makebox{\boldmath$x$}}_{j}}$. We note
$dx^{i}=x^{i}-x^{i-1}$ the vector-valued deviation.

Let us refine the mean-value AE extensions using such a partition. The first
natural idea is to compute an under-approximation for each sub-box obtained as
product of sub-intervals in each dimension. But in general, the under-
approximating boxes will be non-contiguous, and their convex union is in
general not an under-approximation of
$\mbox{range}(f,{\makebox{\boldmath$x$}})$. Moreover, this approach would not
scale well.

We now propose a scheme that avoids these unions, and remains linear in $k$
with respect to the non-partitioned case. We note
${\makebox{\boldmath$x$}}^{1}=[x_{1}^{-1},x_{1}^{1}]\times[x_{2}^{-1},x_{2}^{1}]\times\ldots\times[x_{m}^{-1},x_{m}^{1}]$,
and for all $i$ between $2$ and $k$,
${\makebox{\boldmath$x$}}^{i}=[x_{1}^{-i},x_{1}^{i}]\times\ldots\times[x_{m}^{-i},x_{m}^{i}]\setminus\mathring{{\makebox{\boldmath$x$}}}^{i-1}$,
where $\setminus$ denotes the set difference and
$\mathring{{\makebox{\boldmath$x$}}}$ the interior of $x$. This partition is
represented in Figure 1(b) for a two-dimensional input space. In practice,
each ”square ring” ${\makebox{\boldmath$x$}}^{i}$ will be decomposed in $2n$
sub-boxes for the Jacobian evaluation.

By the mean-value theorem, $\forall x\in[x^{-1},x^{1}]$,
$\exists\xi^{1}\in[x^{-1},x^{1}],$ $f(x)=f(x^{0})+\langle\nabla
f(\xi^{1}),x-x^{0}\rangle$. Suppose we can compute
${\makebox{\boldmath$f$}}^{0}\supseteq f(x^{0})$ and
${\makebox{\boldmath$\nabla$}}^{i}$ for $i$ in $[1,k]$ such that
$\\{\lvert\nabla
f(x)\rvert,x\in{\makebox{\boldmath$x$}}^{i}\\}\subseteq{\makebox{\boldmath$\nabla$}}^{i}$.
We have
$\mbox{range}(f,{\makebox{\boldmath$x$}}^{1})\subseteq{\makebox{\boldmath$f$}}^{0}+\langle\overline{\nabla}^{1},dx^{1}\rangle[-1,1]$
and
$[\overline{f^{0}}-\langle\underline{\nabla}^{1},dx^{1}\rangle,\underline{f^{0}}+\langle\underline{\nabla^{1}},dx^{1}\rangle]\subseteq\mbox{range}(f,[x^{-1},x^{1}])$.
Let us now take $x\in{\makebox{\boldmath$x$}}^{2}$. We can iterate the mean-
value theorem on the adjacent input subdivision and write that for all
$x\in{\makebox{\boldmath$x$}}^{2}$, there exists
$x^{1}\in{\makebox{\boldmath$x$}}^{1}\cap{\makebox{\boldmath$x$}}^{2}$ (that
is on the border between ${\makebox{\boldmath$x$}}^{1}$ and
${\makebox{\boldmath$x$}}^{2}$), there exists
$\xi^{2}\in{\makebox{\boldmath$x$}}^{2}$ such that
$f(x)=f(x^{1})+\langle\nabla f(\xi^{2}),x-x^{1}\rangle$ and $\lvert
x_{1}-x^{1}_{1}\rvert\leq dx_{1}^{2}$ and $\lvert x_{2}-x^{1}_{2}\rvert\leq
dx_{2}^{2}$. (take for example for $x^{1}$ the intersection of the line from
$x^{0}$ to $x$ with the border between ${\makebox{\boldmath$x$}}^{1}$ and
${\makebox{\boldmath$x$}}^{2}$). We have
$\mbox{range}(f,{\makebox{\boldmath$x$}}^{1}\cup{\makebox{\boldmath$x$}}^{2})\subseteq{\makebox{\boldmath$f$}}^{0}+\langle\overline{\nabla}^{1},dx^{1}\rangle[-1,1]+\langle\overline{\nabla}^{2},dx^{2}\rangle[-1,1].$
There also exists
$(x,x^{1})\in{\makebox{\boldmath$x$}}^{2}\times{\makebox{\boldmath$x$}}^{1}$
such that $\lvert x_{1}-x^{1}_{1}\rvert=dx^{2}_{1}$ and $\lvert
x_{2}-x^{1}_{2}\rvert=dx^{2}_{2}$ (take the corners of the boxes
${\makebox{\boldmath$x$}}^{1}$ and ${\makebox{\boldmath$x$}}^{2}$), so that we
also have
$[\overline{f^{0}}-\langle\underline{\nabla}^{1},dx^{1}\rangle-\langle\underline{\nabla}^{2},dx^{2}\rangle,\underline{f^{0}}+\langle\underline{\nabla}^{1},dx^{1}\rangle+\langle\underline{\nabla}^{2},dx^{2}\rangle]\subseteq\mbox{range}(f,{\makebox{\boldmath$x$}}^{1}\cup{\makebox{\boldmath$x$}}^{2}).$

This generalizes to the $k$ subdivisions:

$\displaystyle\mbox{range}(f,{\makebox{\boldmath$x$}})\subseteq{\makebox{\boldmath$f$}}^{0}+\sum_{i=1}^{k}\langle\overline{\nabla}^{i},dx^{i}\rangle[-1,1]$
(7)
$\displaystyle[\overline{f^{0}}-\sum_{i=1}^{k}\langle\underline{\nabla}^{i},dx\rangle,\underline{f^{0}}+\sum_{i=1}^{k}\langle\underline{\nabla}^{i},dx\rangle]\subseteq\mbox{range}(f,{\makebox{\boldmath$x$}})$
(8)

The same idea applies to the estimation of robust ranges.

Naturally, other schemes can be proposed, relying on the idea that this
technique can be seen as using a quadrature formula for integrating the
Jacobian of a function.

###### Example 4.3

We consider
$f(x)=(2x_{1}^{2}+2x_{2}^{2}-2x_{1}x_{2}-2,x_{1}^{3}-x_{2}^{3}+4x_{1}x_{2}-3)^{\intercal}$
with ${\makebox{\boldmath$x$}}=[0.9,1.1]^{2}$. The results are represented in
Figure 1(c). The sampling-based estimation of the image is the dark dots-
filled region. We choose $pi:(1\rightarrow 1,2\rightarrow 2)$. Using the
preconditioned mean-value extension without partitioning, the over-
approximation is the largest green parallelotope and the under-approximation
for the joint range is empty. The quadrature formula for the mean-value
extension with $k=10$ partitions on one hand, and the order 2 extension of
Corollary 3.2 on the other hand, yield two very similar under-approximating
yellow parallelotopes. They also yield two very similar green over-
approximating parallelotopes. Both approaches actually have comparable
precision on the different examples tested. The light green box is the order 2
over-approximation without preconditioning.

###### Remark 4.4

One could be tempted to use more partitions to improve the quality of the
approximation. However, the computations yield approximations that are
centered at $f(x^{0})$. We can observe that the under-approximating skewed box
is actually already very close to being the largest skewed box entirely
included in the image, given a fixed skewing and a center at $f(x^{0})$. The
order 2 estimation allows for slight decentering, but it would also be
possible to use less basic quadrature formulae for that purpose. Quadrature
can also be combined with order 2 extensions.

Quadrature can also be combined with some classical partition of the inputs,
or with different center choices. However, a possibly disjoint union of
approximations (corresponding to a classical partition of inputs) is
inconvenient, so we would recommend this use only as a property-driven
refinement. Finally, as a sound under-approximation is still obtained by
considering a sub-region of the input set, refinements can be obtained by
detecting and removing sub-regions where Jacobian coefficients are either very
sensitive to the inputs or close to zero.

###  4.3 Bounding the Jacobian matrix

The approach relies on being able to compute over-approximations of $\nabla
f(x)$ over some sub-sets of input box $x$, namely
${\makebox{\boldmath$\nabla$}}^{i}$ for $i$ in $[1,k]$ such that
$\\{\lvert\nabla
f(x)\rvert,x\in{\makebox{\boldmath$x$}}^{i}\\}\subseteq{\makebox{\boldmath$\nabla$}}^{i}$.
Automatic differentiation allows to compute the derivatives, but need to be
combined with set-membership methods to handle uncertainties. We have found
the combination of automatic differentiation with an evaluation in affine
arithmetic to provide a good trade-off between efficiency and precision. All
the experiments presented in this work were performed with an implementation
relying on this combination. Affine forms provide an interesting combination
of parameterization and set-based estimation: a parametric approximate form
for $\nabla f(x)$ valid on all box $x$ is computed, that can be instantiated
on ${\makebox{\boldmath$x$}}^{i}$ to yield tight over-approximations
$\nabla^{i}$, without the need for performing several evaluations of the
differentiation. The full description is out of scope, we give below a flavor
of the use of affine arithmetic for Jacobian estimation on a simple example.

###### Example 4.5

Consider in Example 4.3 the derivative of
$f_{1}(x_{1},x_{2})=2x_{1}^{2}+2x_{2}^{2}-2x_{1}x_{2}-2$ with respect to
$x_{1}$, with $(x_{1},x_{2})\in[0.9,1.1]^{2}$. This derivative is
$\nabla_{1}f_{1}(x_{1},x_{2})=4x_{1}-2x_{2}$. Evaluation with affine
arithmetic first consists in creating a centered form with a fresh noise
symbol for each input: $\hat{x}_{1}=1.+0.1\varepsilon_{1}$ and
$\hat{x}_{2}=1.+0.1\varepsilon_{2}$, with
$(\varepsilon_{1},\varepsilon_{2})\in[-1,1]^{2}$. The gradient evaluated on
these affine forms is
$\hat{\nabla}_{1}\hat{f}_{1}(x)=2+0.4\varepsilon_{1}-0.2\varepsilon_{2}.$ Here
the abstraction is exact. In the general case, affine arithmetic will compute
an approximate affine form and bound the approximation error in a new noise
term. Let us now use this affine form to over-approximate $\nabla_{1}f_{1}(x)$
over all $x$. This amounts to computing the interval bounds when
$\varepsilon_{1}$ and $\varepsilon_{2}$ range in $[-1,1]$, for which we get
${\makebox{\boldmath$\nabla$}}_{1,1}=[1.4,2.6]$. Now we use the same affine
form to over-approximate $\nabla_{1}f_{1}(x)$ over
${\makebox{\boldmath$x$}}^{1}$. This amounts to take $\varepsilon_{1}$ and
$\varepsilon_{2}$ both ranging in $[-\frac{1}{k},\frac{1}{k}]$. Let us
consider $k=10$ subdivisions, we obtain
${\makebox{\boldmath$\nabla$}}^{1}_{1,1}=2+0.4[-0.1,0.1]-0.2[-0.1,0.1]=[1.94,2.06]$.
The process can be iterated to compute
${\makebox{\boldmath$\nabla$}}^{i}_{1,1}$ for $i\in[2\ldots k]$, expressing
the membership to a subset ${\makebox{\boldmath$x$}}^{i}$ of $x$ as
constraints on $\varepsilon_{1}$ and $\varepsilon_{2}$, used to instantiate to
${\makebox{\boldmath$x$}}^{i}$ the estimation of $\nabla_{1}f_{1}(x)$. For
instance, in order to compute ${\makebox{\boldmath$\nabla$}}^{2}_{1,1}$, we
decompose the computation on the 4 rectangles that define
${\makebox{\boldmath$x$}}_{2}$ (see Figure 1(b)):
$(\varepsilon_{1},\varepsilon_{2})\in([\frac{1}{k},\frac{2}{k}]\times[-\frac{2}{k},\frac{2}{k}])\cup([-\frac{2}{k},-\frac{1}{k}]\times[-\frac{2}{k},\frac{2}{k}])\cup([-\frac{1}{k},\frac{1}{k}]\times[-\frac{2}{k},-\frac{1}{k}])\cup([-\frac{1}{k},\frac{1}{k}]\times[\frac{1}{k},\frac{2}{k}])$.

##  5 Application to the reachability of discrete-time systems

We consider discrete-time non-linear dynamical systems with inputs of the form

$\begin{cases}z^{k+1}=f(z^{k},u^{k})\\\
z^{0}\in{\makebox{\boldmath$z$}}^{0}\end{cases}$ (9)

where $f:\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$ is a vector-valued non-
linear function with $m\geq n$, $z\in\mathbb{R}^{n}$ the vector of state
variables, $u\in{{\makebox{\boldmath$u$}}}\subseteq\mathbb{R}^{m-n}$ the input
signal, and ${\makebox{\boldmath$z$}}^{0}$ the initial set.

Given an initial set ${\makebox{\boldmath$z$}}^{0}$, we want to compute the
bounded time reachable set of the dynamical system, i.e, the set of states
visited by the dynamical system up to a fixed time horizon $K\in\mathbb{N}$.
The reachable set can be obtained as the solution of the recursion
${\makebox{\boldmath$z$}}^{k+1}=\\{f(z^{k},u^{k})|z^{k}\in{\makebox{\boldmath$z$}}^{k},u^{k}\in{{\makebox{\boldmath$u$}}}\\}$,
for $k\in[0,K]$. The computation of the reachable set can be seen as a series
of images of sets by vector-valued function $f$. We thus can use the results
of Sections 2 to 4 to compute approximations of these reachable sets.

For conciseness, we consider systems without disturbances and compute maximal
(or classical) reachable sets. The algorithms can be straightforwardly
extended to robust reach set of systems with disturbances, basically replacing
ranges by robust ranges. This allows us to use the lighter notations
${\mathcal{I}}(f,{\makebox{\boldmath$x$}},\pi)$ and
${\mathcal{O}}(f,{\makebox{\boldmath$x$}},\pi)$ to note the under and over-
approximating sets introduced in Definition 2.2.

####  5.0.1 Method 1

the first method consists in iteratively using function image, independently
for under and over-approximation, taking as input the previously computed
approximation of the image. We compute under and over-approximations $I^{k}$
and $O^{k}$ of the reachable set ${\makebox{\boldmath$z$}}^{k}$ by

$\begin{cases}I^{0}={\makebox{\boldmath$z$}}^{0},\;O^{0}=z^{0}\\\
I^{k+1}={\mathcal{I}}(f,I^{k},\pi),\;O^{k+1}={\mathcal{O}}(f,O^{k},\pi)\end{cases}$
(10)

Indeed, at each step $k$, we have
$I^{k+1}\subseteq\mbox{range}(f,I^{k})\subseteq\mbox{range}(f,{\makebox{\boldmath$z$}}^{k})={\makebox{\boldmath$z$}}^{k+1}\subseteq\mbox{range}(f,O^{k})\subseteq
O^{k+1}$.

With this approach, at each step $k$, under and over-approximations $I^{k}$
and $O^{k}$ of the joint range are used as input for the next step. It is thus
particularly important to compute tight under and over-approximations of this
joint range. In particular, using the preconditioning of Section 4.1 will
often be crucial, both for under and over approximation. This yields Algorithm
1.

Algorithm 1 Iterated discrete-time reachability

0: $f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$,
${\makebox{\boldmath$z$}}^{0}\subseteq{\mathbb{I}\mathbb{R}}^{n}$ initial
state, $K\in\mathbb{N}^{+}$, an over-approximating extension $[\nabla f]$ (see
Section 4)

0: $I^{k}$ and $O^{k}$: under and over-approximations of the reachable set
$\mbox{range}(f^{k},{\makebox{\boldmath$z$}}^{0})$ for $k\in[1,K]$

$I^{0}\coloneqq{\makebox{\boldmath$z$}}^{0},O^{0}\coloneqq{\makebox{\boldmath$z$}}^{0}$;
choose $\pi:[1\ldots n]\mapsto[1\ldots n]$

for $k$ from 0 to $K-1$ do

${\makebox{\boldmath$\nabla$}}_{I}^{k}\coloneqq|[\nabla f](I^{k})|$,
${\makebox{\boldmath$\nabla$}}_{O}^{k}\coloneqq|[\nabla f](O^{k})|$

$A_{I}^{k}\coloneqq c({\makebox{\boldmath$\nabla$}}_{I}^{k})$,
$A_{O}^{k}\coloneqq c({\makebox{\boldmath$\nabla$}}_{O}^{k})$ (supposed non-
singular, otherwise taken to identity matrix)

$C_{I}^{k}\coloneqq(A_{I}^{k})^{-1}$, $C_{O}^{k}\coloneqq(A_{O}^{k})^{-1}$

${\makebox{\boldmath$z$}}_{I}^{k+1}\coloneqq{\mathcal{I}}(C_{I}^{k}f,I^{k},\pi)$,
${\makebox{\boldmath$z$}}_{O}^{k+1}\coloneqq{\mathcal{O}}(C_{O}^{k}f,O^{k},\pi)$

if ${\makebox{\boldmath$z$}}_{I}^{k}=\emptyset$ then

return

end

$I^{k+1}\coloneqq A_{I}^{k}{\makebox{\boldmath$z$}}_{I}^{k+1}$,
$O^{k+1}\coloneqq A_{O}^{k}{\makebox{\boldmath$z$}}_{O}^{k+1}$

end for

At each step $k=0\ldots K-1$, ${\makebox{\boldmath$z$}}_{I}^{k+1}$ is an
interval vector such that, if it is non empty,
$I^{k+1}=A_{I}^{k}{\makebox{\boldmath$z$}}_{I}^{k+1}\subseteq\mbox{range}(f,I^{k})\subseteq\mbox{range}(f^{k+1},{\makebox{\boldmath$z$}}^{0})$.
The over-approximation is computed similarly, and is fully decoupled.

####  5.0.2 Method 2

the second method consists in computing the sensitivity to initial state by
approximating the gradient of the iterated function. At each step $k$, we
compute the under and over-approximation of
$\mbox{range}(f^{k},{\makebox{\boldmath$z$}}^{0})$, i.e. the loop body $f$
iterated $k$ times, starting from the initial state
${\makebox{\boldmath$z$}}^{0}$. This yields the schematic Algorithm 2, with
same inputs and hypotheses as in Algorithm 1.

Algorithm 2 Discrete-time reachability computed on $f^{k}$

for $k$ from 0 to $K-1$ do

$I^{k+1}\coloneqq{\mathcal{I}}(f^{k+1},{\makebox{\boldmath$z$}}^{0},\pi)$,
$O^{k+1}\coloneqq{\mathcal{O}}(f^{k+1},{\makebox{\boldmath$z$}}^{O},\pi)$

end for

Here, at each step $k$, the under- and over-approximation are both obtained
from an over-approximation of $f^{k+1}$ evaluated at the center of
${\makebox{\boldmath$z$}}^{0}$ and an over-approximation of the gradient of
$f^{k+1}$ over ${\makebox{\boldmath$z$}}^{0}$: at each step, the gradient can
be obtained by differentiating the gradient from the previous step. Of course,
this can also be combined with preconditioning.

####  5.0.3 Discussion

While relying on the same techniques for range estimation, Algorithm 1 and
Algorithm 2 are different in spirit: Algorithm 2 relies only on the
propagation of over-approximations to deduce under-approximations. In
particular, the under-approximation may be empty at some step $k$, and become
non-empty again at further steps (a similar remark was made in the context of
continuous systems in Goubault and Putot (2017)). In comparison, Algorithm 1
needs at each step an under-approximating box or skew box that is non-empty
for all components. On the other hand, Algorithm 2 is more costly as it
requires a differentiation of the iterated function.

##  6 Implementation and examples

We now present results on small systems from the literature. The approach is
implemented as part of the RINO C++ prototype, available from
https://github.com/cosynus-lix/RINO. The prototype allows to experiment the
function range estimation, but its actual target is discrete and continuous-
time reachability, combining for continuous-time the techniques presented here
with Taylor model methods Goubault and Putot (2019, 2020). The prototype uses
the fadbad++ (http://www.fadbad.com/fadbad.html) automatic differentiation
library and the aaflib (http://aaflib.sourceforge.net/) affine arithmetic
library. The timings are given on a Macbook Pro 2.6GHz Intel Core i7 and 32Gb
of RAM.

##### Test Model

We consider the test model Dreossi et al. (2016):

$\begin{split}x_{1}^{k+1}&=x_{1}^{k}+(0.5(x_{1}^{k})^{2}-0.5(x_{2}^{k})^{2})\Delta\\\
x_{2}^{k+1}&=x_{2}^{k}+2x_{1}^{k}x_{2}^{k}\Delta\end{split}$

with as initial set a box $x_{1}\in[0.05,0.1]$ and $x_{2}\in[0.99,1.00]$, and
$\Delta=0.01$. Figure 2 shows the under and over-approximated reachable sets
(respectively the filled yellow region and green parallelotope) over time up
to 25 steps with Algorithm 1. They are obtained in 0.02 seconds. We can
observe that the under and over-approximations are very close one to another,
confirming the accuracy of the results.

Figure 2:  Skewed box under and over-approximations for 25 steps (Algorithm 1)
for the test model

##### SIR Epidemic Model

We now consider the SIR epidemic model with the parameters of Dreossi et al.
(2016).

$\begin{split}x_{1}^{k+1}&=x_{1}^{k}-\beta x_{1}^{k}x_{2}^{k}\Delta\\\
x_{2}^{k+1}&=x_{2}^{k}+(\beta x_{1}^{k}x_{2}^{k}-\gamma x_{2}^{k})\Delta\\\
x_{3}^{k+1}&=x_{3}^{k}+\gamma x_{2}^{k}\Delta\end{split}$

We compute the reachable set up to 60 steps from the initial box
$(x_{1},x_{2},x_{3})\in[0.79,0.80]\times[0.19,0.20]\times[0,0.1]$. The
parameter values are $\beta=0.34$, $\gamma=0.05$, and $\Delta=0.5$.

The reachable sets computed in 0.05 seconds with Algorithm 1 up to 60 steps
are represented for $(x_{1},x_{2})$ in Figure 3.

Figure 3:  Skewed box under and over-approximations for 60 steps (Algorithm 1)
of the SIR epidemic model

We can note in particular from the zoomed reachable set in the Figure, which
corresponds to last step (60), that the under-approximation (in yellow) is
still of good quality (the purple dots correspond to sample executions).

However, only Algorithm 2 is able to compute non-empty approximations when
taking as initial condition $x_{3}=0$ which is of empty interior, instead of
$x_{3}\in[0,0.1]$. We obtain in 0.05 seconds the very tight approximations of
the projections of the components represented Figure4(a).

(a) SIR epidemic model, 60 steps

(b) Bees model, 1500 steps

Figure 4: Projected under and over-approximations, Algorithm 2

##### Honeybees Site Choice Model

We consider the reachable sets up to 1500 steps of the model studied in
Dreossi et al. (2016):

$\begin{split}x_{1}^{k+1}&=x_{1}^{k}-(\beta_{1}x_{1}^{k}x_{2}^{k}+\beta_{2}x_{1}^{k}x_{3}^{k})\Delta\\\
x_{2}^{k+1}&=x_{2}^{k}+(\beta_{1}x_{1}^{k}x_{2}^{k}-\gamma
x_{2}^{k}+\delta\beta_{1}x_{2}^{k}x_{4}^{k}+\alpha\beta_{1}x_{2}^{k}x_{5}^{k})\Delta\\\
x_{3}^{k+1}&=x_{3}^{k}+(\beta_{2}x_{1}^{k}x_{3}^{k}-\gamma
x_{3}^{k}+\delta\beta_{2}x_{3}^{k}x_{5}^{k}+\alpha\beta_{2}x_{3}^{k}x_{4}^{k})\Delta\\\
x_{4}^{k+1}&=x_{4}^{k}+(\gamma
x_{2}^{k}-\delta\beta_{1}x_{2}^{k}x_{4}^{k}-\alpha\beta_{2}x_{3}^{k}x_{4}^{k})\Delta\\\
x_{5}^{k+1}&=x_{5}^{k}+(\gamma
x_{3}^{k}-\delta\beta_{2}x_{3}^{k}x_{5}^{k}-\alpha\beta_{1}x_{2}^{k}x_{5}^{k})\Delta\\\
\end{split}$

with as initial set the box $x_{1}=500$, $x_{2}\in[390,400]$,
$x_{3}\in[90,100]$, $x_{4}=x_{5}=0$ and the parameter values
$\beta_{1}=\beta_{2}=0.001$, $\gamma=0.3$, $\delta=0.5$, $\alpha=0.7$, and
$\Delta=0.01$.

Algorithm 1 runs very fast, taking only 1.7 seconds for the 1500 steps. But
the under-approximation is very soon empty and the over-approximation tends to
strongly widen after 800 steps. In comparison, Algorithm 2 takes 57 seconds to
complete the 1500 reachability steps. But the projected under-approximations
are very tight, close to the over-approximations, as can be seen on Figure
4(b) which represents under and over-approximations for all components as
functions of steps. These results should also be compared to the much wider
over-approximation of Dreossi et al. (2016) (Figure 7, where time is number of
steps divided by 100), obtained in 81 seconds.

##  7 Conclusion and future work

We focused on new AE under-approximating extensions and their accurate
practical evaluation for non-linear vector-valued functions, and exemplified
their interest for the reachability of discrete-time systems. These techniques
can also be used for the reachability analysis of continuous-time systems,
improving for instance over Goubault and Putot (2019, 2020).

## References

  * Chen et al. (2016) Chen, M., Herbert, S., and Tomlin, C.J. (2016).  Exact and efficient Hamilton-Jacobi-based guaranteed safety analysis via system decomposition.
  * Chen et al. (2014) Chen, X., Sankaranarayanan, S., and Ábrahám, E. (2014).  Under-approximate flowpipes for non-linear continuous systems.  In _FMCAD_. IEEE.
  * Dreossi et al. (2016) Dreossi, T., Dang, T., and Piazza, C. (2016).  Parallelotope bundles for polynomial reachability.  In _Proceedings of the 19th International Conference on Hybrid Systems: Computation and Control_. ACM Press.
  * Girard et al. (2006) Girard, A., Le Guernic, C., and Maler, O. (2006).  Efficient computation of reachable sets of linear time-invariant systems with inputs.  In J.P. Hespanha and A. Tiwari (eds.), _Hybrid Systems: Computation and Control_ , 257–271. Springer Berlin Heidelberg, Berlin, Heidelberg.
  * Goldsztejn (2012a) Goldsztejn, A. (2012a).  Modal intervals revisited, part 1: A generalized interval natural extension.  _Reliable Computing_ , 16, 130–183.
  * Goldsztejn (2012b) Goldsztejn, A. (2012b).  Modal intervals revisited, part 2: A generalized interval mean value extension.  _Reliable Computing_ , 16, 184–209.
  * Goldsztejn and Jaulin (2010) Goldsztejn, A. and Jaulin, L. (2010).  Inner approximation of the range of vector-valued functions.  _Reliable Computing_ , 14.
  * Goubault and Putot (2017) Goubault, E. and Putot, S. (2017).  Forward inner-approximated reachability of non-linear continuous systems.  In _HSCC_. ACM.
  * Goubault and Putot (2019) Goubault, E. and Putot, S. (2019).  Inner and outer reachability for the verification of control systems.  In _HSCC_.
  * Goubault and Putot (2020) Goubault, E. and Putot, S. (2020).  Robust under-approximations and application to reachability of non-linear control systems with disturbances.  _IEEE Control Systems Letters_ , 4(4), 928–933.
  * Guernic and Girard (2009) Guernic, C.L. and Girard, A. (2009).  Reachability analysis of hybrid systems using support functions.  In _Proceedings of 21st International Conference on Computer Aided Verification_ , LNCS.
  * Kochdumper and Althoff (2020) Kochdumper, N. and Althoff, M. (2020).  Computing non-convex inner-approximations of reachable sets for nonlinear continuous systems.  In _CDC_.
  * Korda et al. (2013) Korda, M., Henrion, D., and Jones, C.N. (2013).  Inner approximations of the region of attraction for polynomial dynamical systems.  In _9th IFAC Symposium on Nonlinear Control Systems, NOLCOS 2013_.
  * Kurzhanski and Varaiya (2000) Kurzhanski, A.B. and Varaiya, P. (2000).  Ellipsoidal techniques for reachability analysis.  In N. Lynch and B.H. Krogh (eds.), _HSCC_ , 202–214.
  * Makino and Berz (2003) Makino, K. and Berz, M. (2003).  Taylor models and other validated functional inclusion methods.  _Int. J. Pure Appl. Math_.
  * Mitchell (2007) Mitchell, I.M. (2007).  Comparing forward and backward reachability as tools for safety analysis.  In _HSCC_.
  * Mézo et al. (2018) Mézo, T.L., Jaulin, L., and Zerr, B. (2018).  Bracketing the solutions of an ordinary differential equation with uncertain initial conditions.  _Applied Mathematics and Computation_ , 318.
  * Raković and Fiacchini (2008) Raković, S.V. and Fiacchini, M. (2008).  Approximate reachability analysis for linear discrete time systems using homothety and invariance.  _IFAC Proceedings Volumes_ , 41(2), 15327 – 15332.
  * Xue et al. (2020) Xue, B., Fränzle, M., and Zhan, N. (2020).  Inner-approximating reachable sets for polynomial systems with time-varying uncertainties.  _IEEE Transactions on Automatic Control_ , 65(4), 1468–1483.  10.1109/TAC.2019.2923049.
  * Xue et al. (2016) Xue, B., She, Z., and Easwaran, A. (2016).  Under-approximating backward reachable sets by polytopes.  In _Computer Aided Verification_. Springer.