# Tight Remainder-Form Decomposition Functions
with Applications to Constrained Reachability and Interval Observer Design

Mohammad Khajenejad, and Sze Zheng Yong M. Khajenejad and S.Z. Yong are with
the School for Engineering of Matter, Transport and Energy, Arizona State
University, Tempe, AZ, USA (e-mail: {mkhajene, szyong}@asu.edu).

###### Abstract

This paper proposes a tractable family of remainder- form mixed-monotone
decomposition functions  that are useful for over-approximating the image set
of nonlinear mappings in reachability and estimation problems. In particular,
our approach applies to a new class  of nonsmooth nonlinear systems that we
call either-sided locally Lipschitz (ELLC) systems, which we show to be a
superset of locally Lipschitz continuous systems,  thus expanding the set of
systems that are formally known to be mixed-monotone.  In addition, we derive
lower and upper bounds for the over-approximation error  and show that the
lower bound is achievable  with our proposed approach. Moreover, we develop a
set inversion algorithm that along with the proposed decomposition functions,
can be used for constrained reachability analysis and interval observer design
for continuous and discrete-time systems with bounded noise.

##  I Introduction

Monotonicity  properties of systems have proven to be very powerful and useful
for analyzing and controlling complex systems [1, 2]. Building upon this idea,
it was shown that certain nonmonotone systems can be lifted to higher
dimensional monotone systems that can be potentially used to deduce critical
information about the original system (see e.g., [3, 4, 5]) by decomposing the
system dynamics into increasing and decreasing components. Systems that are
decomposable in this manner are called _mixed-monotone_ and are significantly
more general than the class of monotone systems.

Furthermore, mixed-monotonicity has also proven to be very beneficial for
system analysis and control. For instance, if mixed-monotonicity holds, it can
be concluded that the original system has global asymptotic stability by
proving the nonexistence of equilibria of the  lifted system, except in a
certain lower-dimensional subspace [6, 7].  Moreover, forward invariant and
attractive sets of the original system can be identified [8] and reachable
sets of the original system can be efficiently  over-approximated  and used
for state estimation and abstraction-based control synthesis [8, 9, 10, 11].
However, the usefulness of these system lifting techniques for analysis and
control is highly dependent on the _tightness_ of the decomposition
approaches, their _computational tractability_ [12, 13] and their
applicability to a broad class of systems; therefore, the capability to
compute or construct _tractable_ and _tight_ mixed-monotone decomposition
functions for a broad range of nonlinear, uncertain and constrained systems is
of great interest and will have a significant impact.

_Literature review_.  Mixed-monotone decomposition functions are generally not
unique, hence several seminal studies have addressed the issue of identifying
and/or constructing appropriate decomposition functions with somewhat
different yet highly related definitions and corresponding sufficient
conditions for mixed-monotonicity [12, 14, 13, 15, 16, 10, 17, 7]. In
particular, recent studies in [12, 13] provided _tight_ decomposition
functions for unconstrained discrete -time and continuous-time dynamical
systems, respectively,  whose _computability_ rely on  the global solvability
of nonlinear optimization programs, which is  only guaranteed in  some
specific cases such as when the vector field is _Jacobian sign-stable_ , or
when all  its _stationary/ critical points_ can be  computed analytically. On
the other hand,  an interesting study in [14] proposed computable and
constructive (but not necessarily tight) decomposition functions for
differentiable vector fields with known bounds for the derivatives.  Building
on these frameworks, we aim to obtain computable/tractable and tighter
decomposition functions for a broad class of nonsmooth systems.

Another relevant body of literature pertains to _interval arithmetic_ [18, 19,
20, 21, 22], which has been successfully applied to problems in numerical
analysis, set estimation, motion planning, etc. Specifically, _inclusion
functions_ and variations thereof (e.g., natural, centered-form and mixed-form
inclusions) that are based on interval arithmetic can be directly related to
decomposition functions and can similarly be used for over-approximation of
the ranges/image sets of functions and for state observer designs [23, 24, 25,
26, 27, 22]. In addition, various properties of inclusion functions have been
studied such as their convergence rates and the _subdivision principle_ [19].
Subsequent studies further introduced refinements of interval over-
approximations and set inversion algorithms that can incorporate new sources
of information about the system such as state constraints,
measurements/observations, manufactured redundant variables, second-order
derivatives, etc., e.g., in [18, 23, 24, 25, 26],  which can be also
beneficial for decomposition function-based constrained reachability and set-
valued estimation problems that we consider in this work.

_Contribution._ In this paper, we introduce a class of mixed-monotone
decomposition functions whose construction is both computationally tractable
and tight for a broad range of discrete-time and continuous-time nonsmooth
nonlinear systems. The proposed mixed-monotone decomposition functions are in
the _remainder form_ (based on the terminology in [19]) and we also show that
they result in _difference of monotone functions_ , which bear some
resemblance with difference of convex functions in DC programming that is also
widely used in optimization and state estimation problems, e.g., [25, 28].

Our work contributes to the literature on mixed-monotone decomposition
functions and more generally, inclusion functions in multiple ways:

  1. 1)

First, our proposed remainder-form decomposition functions apply to a new
broader class of nonsmooth nonlinear systems, which we call _either-sided
locally Lipschitz continuous_ (ELLC) systems, that is proven in the paper to
be a superset of locally Lipschitz continuous (LLC) systems. This new system
class relaxes the (almost everywhere) smoothness and bounded gradient
requirements and allows nonsmooth vector fields with only one-sided bounded
_Clarke Jacobians_.

  2. 2)

Our construction approach for remainder-form decomposition functions is
tractable (i.e., computable in closed-form with a non-iterative algorithm) and
proven to be the tightest for this family of decomposition functions, which we
also show to include [14], thus generalizing and improving on the method in
[14].

  3. 3)

We further show that the over-approximation of the image sets of ELLC systems
converges to the true tightest enclosing interval at least linearly, when the
interval domain width goes to zero and moreover, a _subdivision principle_
applies for improving the enclosure of the range/image set.

  4. 4)

We introduce a novel set inversion algorithm based on mixed-monotone
decomposition functions as an alternative to SIVIA [18] and the refinement
algorithm in [23], which enables the design of algorithms for constrained
reachability and interval observers for systems with known constraints,
modeling redundancy and/or sensor measurements.

Moreover, it is noteworthy that the inclusion function based on mixed-monotone
decomposition functions can be used alongside any existing inclusion
functions, where the “best of them” (by virtue of an intersection property) is
chosen, since we observed that in general, no single inclusion function
consistently outperforms all others. Further more, our proposed inclusion
function can be directly and simply integrated in to existing set inversion,
constrained reachability and interval observer algorithms, e.g., [18, 23, 29],
without much modifications. Finally, to demonstrate the effectiveness of the
proposed algorithms for decomposition function construction and set inversion,
we compare them with existing inclusion/decomposition functions in the
literature.

##  II Background and Problem Formulation

###  II-A Notation

$\mathbb{N}$,  $\mathbb{N}_{a}$, $\mathbb{R}^{n_{z}}$ and $\mathbb{R}^{n\times
m}$ denote the set of positive integers,  the first $a$ positive integers, the
$n_{z}$-dimensional Euclidean space and the space of $n$ by $m$ real matrices,
respectively. Moreover, $\forall
z,\underline{z},\overline{z}\in\mathbb{R}^{n_{z}}$,
$\underline{z}\leq\overline{z}\Leftrightarrow\underline{z}_{i}\leq\overline{z}_{i},\forall
i\in{\color[rgb]{0,0,0}\mathbb{N}_{n_{z}}}$, where ${z}_{i}$ denotes the
$i$-th  element of $z$. Further,
${\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z}\triangleq[\underline{z},\overline{z}]\triangleq{\color[rgb]{0,0,0}\\{}z\in\mathbb{R}^{n_{z}}|\underline{z}\leq
z\leq\overline{z}{\color[rgb]{0,0,0}\\}}$ and
$d(\mathbb{I}\mathcal{Z})\triangleq\|\overline{z}-\underline{z}\|_{\infty}$
are called a closed interval/hyperrectangle in $\mathbb{R}^{n_{z}}$ and the
diameter of ${\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z}$,  respectively, where
$\|z\|_{\infty}\triangleq\max\limits_{i}|z_{i}|$ denotes the
$\ell_{\infty}$-norm of $z\in\mathbb{R}^{n_{z}}$. The set of all intervals in
$\mathbb{R}^{n_{z}}$ is denoted by $\mathbb{IR}^{n_{z}}$.

###  II-B Definitions and Preliminaries

First, we briefly introduce some of the main concepts that we use through out
the paper, as well as some important existing results that will be used  for
comparisons and for deriving our main results. We start by introducing
inclusion  and decomposition functions and some of their  typical instances.

####  II-B1  Inclusion and Decomposition Functions

###### Definition 1 (Inclusion Functions).

[18, Chapter 2.4] Consider a function
$f:\mathcal{Z}\subset\mathbb{R}^{n_{z}}\to\mathbb{R}^{n_{x}}$. The interval
function $T^{f}:\mathbb{IR}^{n_{z}}\to\mathbb{IR}^{n_{x}}$ is an inclusion
function for $f(\cdot)$, if

$\displaystyle\forall{\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z}\in\mathbb{IR}^{n_{z}},f({\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z})\subset
T^{f}({\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z}),$

where $f({\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z})$ is the true image set
(or range) of $f(\cdot)$ for the domain
${\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z}\in\mathbb{IR}^{n_{z}}$.  The
tightest enclosing interval of $f({\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z})$
is denoted by
$T^{f}_{O}(\mathbb{I}\mathcal{Z})\triangleq[\underline{f}^{\text{true}},\overline{f}^{\text{true}}]{\color[rgb]{0,0,0}\triangleq[\min\limits_{z\in\mathbb{I}\mathcal{Z}}f(z),\max\limits_{z\in\mathbb{I}\mathcal{Z}}f(z)]}\supset
f({\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z})$; hence, it is the tightest
inclusion function, i.e.,
$T^{f}_{O}(\mathbb{I}\mathcal{Z})\subseteq{{T}}^{f}({\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z})$.
Further, with a slight abuse of notation, we overload the notation of
$f(\mathbb{I}\mathcal{Z})$ and
$T^{f}({\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z})$ when (1) is a continuous-
time system to represent
$f_{c}(\mathbb{I}\mathcal{Z})\subset\mathbb{R}^{n_{x}}$ and
$T^{f_{c}}({\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z})\in\mathbb{I}\mathbb{R}^{n_{x}}$,
respectively, with $f_{c,i}(\mathbb{I}\mathcal{Z})\triangleq
f_{i}(\mathbb{I}\mathcal{Z}_{c,i})$,
$T^{f_{c}}_{i}(\mathbb{I}\mathcal{Z})\triangleq
T^{f}_{i}(\mathbb{I}\mathcal{Z}_{c,i})$ and
$\mathbb{I}\mathcal{Z}_{c,i}\triangleq\\{{\color[rgb]{0,0,0}[\underline{z},\overline{z}]\in\mathbb{I}\mathcal{Z}}\,|\,[\underline{z}_{j},\overline{z}_{j}]=\mathbb{I}\mathcal{Z}_{j},\forall
j\neq i,\underline{z}_{i}=\overline{z}_{i}=x_{i}\\}$, $\forall
i\in\mathbb{N}_{n_{x}}$.

###### Proposition 1 (Natural ($T^{f}_{N}$) Inclusion Functions).

[18, Theorem 2.2] Consider
${\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z}\triangleq[\underline{z},\overline{z}]\in\mathbb{IR}^{n_{z}}$
and
$f\triangleq[f_{1},\dots,f_{n_{x}}]^{\top}:{\color[rgb]{0,0,0}\mathcal{Z}\subset\mathbb{R}^{n_{z}}}\to\mathbb{R}^{n_{x}}$,
where each $f_{j},\ j\in{\color[rgb]{0,0,0}\mathbb{N}_{n_{x}}}$,  is expressed
as a finite composition of the operators $+,-,\times,/$ and elementary
functions (sine, cosine, exponential, square root,  etc.). A natural inclusion
function
$T^{f}_{N}:\mathbb{IR}^{n_{z}}\to{\color[rgb]{0,0,0}\mathbb{I}}\mathbb{R}^{n_{x}}$
for $f(\cdot)$ is obtained by replacing each real variable
$z_{i},i\in{\color[rgb]{0,0,0}\mathbb{N}_{n_{z}}}$, by its corresponding
_interval variable_
$[z_{i}]\triangleq{\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z}_{i}=[\underline{z}_{i},\overline{z}_{i}]$,
and each operator or function by its interval counterpart by applying
_interval arithmetic_ (cf. [18, Chapter 2] for details).

###### Proposition 2 (Centered ($T^{f}_{C}$) and Mixed Centered ($T^{f}_{M}$)
Inclusion Functions).

[18, Sections 2.4.3–2.4.4] Let
$f\triangleq[f_{1},\dots,f_{n_{x}}]^{\top}:{\color[rgb]{0,0,0}\mathcal{Z}\subset\mathbb{R}^{n_{z}}}\to\mathbb{R}^{n_{x}}$
be differentiable over the  interval
${\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z}\triangleq[\underline{z},\overline{z}]\in\mathbb{IR}^{n_{z}}$.
Then, the interval function

$\displaystyle T^{f}_{C}({\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z})\triangleq
f(m)+{\color[rgb]{0,0,0}\mathbb{I}\mathbf{J}^{f}_{\mathbb{I}\mathcal{Z}}}({\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z}-m),$

is an inclusion function for $f(\cdot)$,  called the _centered inclusion
function_ for $f(\cdot)$ in ${\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z}$, where
$m\triangleq\frac{\underline{z}+\overline{z}}{2}$,
${\color[rgb]{0,0,0}\mathbb{I}\mathbf{J}^{f}_{\mathbb{I}\mathcal{Z}}}$ is an
_interval Jacobian matrix_ with domain $\mathbb{I}\mathcal{Z}$, such that
$J_{f}(z)\in{\color[rgb]{0,0,0}\mathbb{I}\mathbf{J}^{f}_{\mathbb{I}\mathcal{Z}}},\forall
z\in{\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z}$ and $J_{f}(z)$ is the Jacobian
matrix of $f(z)$ at point $z\in{\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z}$.
Moreover,

$\displaystyle
T^{f}_{M}({\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z})\triangleq[T^{f}_{M,i}({\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z})\dots
T^{f}_{M,n_{x}}({\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z})]^{\top},$

with $T^{f}_{M,i}({\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z})\triangleq
f_{i}(m)+\sum_{j=1}^{n_{z}}{\color[rgb]{0,0,0}(\mathbb{I}}\mathbf{J}^{f}_{{\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z}_{1\to
j}})_{i,j}({\color[rgb]{0,0,0}\mathbb{I}\mathcal{Z}_{j}}-m_{j})$ for all
$i\in{\color[rgb]{0,0,0}\mathbb{N}_{n_{x}}}$, is also an inclusion function
for $f(\cdot)$, called the _mixed-centered inclusion function_ , where
${\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z}_{j}\triangleq[\underline{z}_{j}\
\overline{z}_{j}],j\in{\color[rgb]{0,0,0}\mathbb{N}_{n_{z}}}$,
${\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z}_{1\to
j}\triangleq[{\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z}_{1}\dots\color[rgb]{0,0,0}\mathbb{I}\mathcal{Z}_{j}\
m_{j+1}\dots m_{n_{z}}]^{\top}$,
$({\color[rgb]{0,0,0}\mathbb{I}}\mathbf{J}^{f}_{{\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z}_{1\to
j}})_{i,j}$ is the $(i,j)$-th element of the interval  Jacobian matrix
${\color[rgb]{0,0,0}\mathbb{I}}\mathbf{J}^{f}_{{\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z}_{1\to
j}}$ with domain $\mathbb{I}\mathcal{Z}_{1\to j}$, and
$\underline{z}_{j},\overline{z}_{j},m_{j}$ are the $j$-th elements of the
vectors $\underline{z},\overline{z},m$, respectively.

Next, inspired by the work in [19, Section 3], we introduce the notion of
remainder-form (additive) inclusion functions.

###### Definition 2 (Remainder-Form (Additive) Inclusion Functions).

Consider a function
$f:{\color[rgb]{0,0,0}\mathcal{Z}\in\mathbb{R}^{n_{z}}}\to\mathbb{R}^{n_{x}}$.
The interval function $T^{f}_{R}:\mathbb{IR}^{n_{z}}\to\mathbb{IR}^{n_{x}}$ is
an additive (remainder-form) inclusion function for $f(\cdot)$, if there exist
two constituent mappings $g,h:\mathbb{R}^{n_{z}}\to\mathbb{R}^{n_{x}}$, such
that for any
${\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z}\in\mathbb{IR}^{n_{z}}$:

$\displaystyle f({\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z})\subseteq
T^{f}_{R}({\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z})\triangleq
g({\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z})+h({\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z}),$

where addition is based on interval arithmetic  (cf. [18]).

###### Definition 3 (Mixed-Monotonicity  and Decomposition Functions).

[13, Definition 1],[14, Definition 4] Consider the dynamical system with
initial state
$x_{0}\in\mathbb{I}\mathcal{X}_{0}\triangleq[\underline{x}_{0},\overline{x}_{0}]$:

$\displaystyle x_{t}^{+}=\tilde{f}(x_{t},u_{t},w_{t})\triangleq f(z_{t}),$ (1)

where $x_{t}^{+}\triangleq x_{t+1}$ if (1) is a discrete-time system and
$x_{t}^{+}\triangleq\dot{x}_{t}$ if (1) is a continuous-time system,
$\tilde{f}:\mathcal{X}\times\mathcal{U}\times\mathcal{W}\to\mathbb{R}^{n_{x}}$
is the vector field with state $x_{t}\in\mathcal{X}\subset\mathbb{R}^{n_{x}}$,
known input $u_{t}\in\mathcal{U}\subset\mathbb{R}^{n_{u}}$ and disturbance
input
$w_{t}\in\mathcal{W}\subseteq\mathbb{I}\mathcal{W}\triangleq[\underline{w},\overline{w}]\in\mathbb{IR}^{n_{w}}$.
For ease of exposition, we also define
$f:\mathcal{Z}\triangleq\mathcal{X}\times\mathcal{W}\subset\mathbb{R}^{n_{z}}\to\mathbb{R}^{n_{x}}$
as in (1) that is implicitly dependent on $u_{t}$ with the augmented state
$z_{t}\triangleq[x_{t}^{\top}\ w_{t}^{\top}]^{\top}\in\mathcal{Z}$.

Suppose (1) is a discrete-time system. Then, a mapping
$f_{d}:\mathcal{Z}\times\mathcal{Z}\to\mathbb{R}^{n_{x}}$ is a  discrete-time
mixed-monotone decomposition function with respect to $f(\cdot)$, if it
satisfies:

  1. i)

$f(\cdot)$ is embedded on the diagonal of $f_{d}(\cdot,\cdot)$, i.e.,
$f_{d}(z,z)=f(z)$.

  2. ii)

$f_{d}$ is monotone increasing in its first argument, i.e., $\hat{z}\geq
z\implies f_{d}(\hat{z},z^{\prime})\geq f_{d}(z,z^{\prime})$.

  3. iii)

$f_{d}$ is monotone decreasing in its second argument, i.e., $\hat{z}\geq
z\implies f_{d}(z^{\prime},\hat{z})\leq f_{d}(z^{\prime},z).$

Further, if (1) is a continuous-time system, a mapping
$f_{d}:\mathcal{Z}\times\mathcal{Z}\to\mathbb{R}^{n_{x}}$ is a continuous-time
mixed-monotone decomposition function with respect to $f(\cdot)$, if it
satisfies:

  1. i)

$f(\cdot)$ is embedded on the diagonal of $f_{d}(\cdot,\cdot)$, i.e.,
$f_{d}(z,z)=f(z)$.

  2. ii)

$f_{d}$ is monotone increasing in its first argument with respect to “off-
diagonal” arguments, i.e., $\forall i\in\mathbb{N}_{n_{x}}$, $\hat{z}_{j}\geq
z_{j},\forall j\in\mathbb{N}_{n_{z}},\hat{z}_{i}=z_{i}=x_{i}\implies
f_{d,i}(\hat{z},z^{\prime})\geq f_{d,i}(z,z^{\prime})$.

  3. iii)

$f_{d}$ is monotone decreasing in its second argument, i.e., $\hat{z}\geq
z\implies f_{d}(z^{\prime},\hat{z})\leq f_{d}(z^{\prime},z).$

Moreover, systems that admit mixed-monotone decomposition functions are called
mixed-monotone systems.

Moreover, we extend the concept of decomposition functions  to _one-sided
decomposition functions_.

###### Definition 4 (One-Sided Decomposition Functions).

Consider
$f:{\color[rgb]{0,0,0}\mathcal{Z}\subset\mathbb{R}^{n_{z}}}\to\mathbb{R}^{n_{x}}$
and suppose there exist two mixed-monotone mappings
$\overline{f}_{d},\underline{f}_{d}:\mathcal{Z}\times\mathcal{Z}\to\mathbb{R}^{n_{x}}$
such that for any  $\underline{z},z,\overline{z}\in\mathcal{Z}$, the following
statement holds:

$\displaystyle\underline{z}\leq
z\leq\overline{z}\implies\underline{f}_{d}(\underline{z},\overline{z})\leq
f(z)\leq\overline{f}_{d}(\overline{z},\underline{z}).$ (2)

Then, $\overline{f}_{d}$ and $\underline{f}_{d}$ are called upper and lower
decomposition functions for $f$ over $\mathbb{I}\mathcal{Z}$, respectively.

Based on the above definitions, the mixed monotone decomposition function can
be viewed as the special case when the upper and lower decomposition functions
coincide. Moreover, one-sided decomposition functions can be obtained from a
family of (upper and lower) decomposition functions that are generally not
unique via the following result:

###### Corollary 1 (Intersection Property).

Suppose $\overline{f}^{1}_{d},\overline{f}^{2}_{d}$ and
$\underline{f}^{1}_{d},\underline{f}^{2}_{d}$ are pairs of upper and lower
decomposition functions for $f$, respectively. Then,
$\min\\{\overline{f}^{1}_{d},\overline{f}^{2}_{d}\\}$ and
$\max\\{\underline{f}^{1}_{d},\underline{f}^{2}_{d}\\}$ (i.e., their
_intersection_) are also upper and lower decomposition function for
$f(\cdot)$, respectively.

###### Proof.

The results follow from the fact that if two mappings with the same domains
and image spaces are both monotonically increasing or decreasing in their
$i$-th arguments, then their “point-wise” minimum and maximum are also
monotonically increasing or decreasing on their $i$-th arguments. ∎

Further, we can slightly generalize the notion of _embedding system with
respect to $f_{d}$_ in [13, (7)], to the _embedding system with respect to
$\overline{f}_{d},\underline{f}_{d}$_, through the following definition.

###### Definition 5 (Generalized Embedding  Systems).

For an $n$-dimensional system (1)  with any one-sided decomposition functions
$\overline{f}_{d},\underline{f}_{d}$, its _generalized embedding system is
the_ $2n$-dimensional system with initial condition
$\begin{bmatrix}\overline{x}_{0}^{\top}&\underline{x}_{0}^{\top}\end{bmatrix}^{\top}$:

$\displaystyle{\color[rgb]{0,0,0}\begin{bmatrix}{\overline{x}}_{t}^{+}\\\
{\underline{x}}_{t}^{+}\end{bmatrix}=\begin{bmatrix}\overline{f}_{d}(\begin{bmatrix}(\overline{x}_{t})^{\top}\,\overline{w}^{\top}\end{bmatrix}^{\top},\begin{bmatrix}(\underline{x}_{t})^{\top}\,\underline{w}^{\top}\end{bmatrix}^{\top})\\\
\underline{f}_{d}(\begin{bmatrix}(\underline{x}_{t})^{\top}\,\underline{w}^{\top}\end{bmatrix}^{\top},\begin{bmatrix}(\overline{x}_{t})^{\top}\,\overline{w}^{\top}\end{bmatrix}^{\top})\end{bmatrix}.}$
(3)

###### Proposition 3 (State Framer Property).

Let system (1) with initial state
$x_{0}\in\mathbb{I}\mathcal{X}_{0}\triangleq[\underline{x}_{0},\overline{x}_{0}]$
be mixed-monotone with a generalized embedding system (3) with respect to
$\overline{f}_{d},\underline{f}_{d}$. Then, for all $t\geq 0$,
$R^{f}(t,\mathbb{I}\mathcal{X}_{0})\subset\mathbb{I}\mathcal{X}_{t}\triangleq[\underline{x}_{t},\overline{x}_{t}]$,
where
$R^{f}(t,\mathbb{I}\mathcal{X}_{0})\triangleq\\{\phi(t,x_{0},{w}{\color[rgb]{0,0,0}{}_{0:t}})\mid
x_{0}\in\mathbb{I}\mathcal{X}_{0}\text{ and
}w_{t}\in\mathbb{I}\mathcal{W},\forall t\geq 0\\}$ is the reachable set at
time $t$ of (1) when initialized within $\mathbb{I}\mathcal{X}_{0}$ and
$\mathbb{I}\mathcal{X}_{t}\triangleq[\underline{x}_{t},\overline{x}_{t}]$ is
the solution to the generalized embedding system (3). Consequently, the system
state trajectory $x_{t}$ satisfies $\underline{x}_{t}\leq
x_{t}\leq\overline{x}_{t},{\color[rgb]{0,0,0}\forall t\geq 0}$, i.e., is
_framed_ by $\mathbb{I}\mathcal{X}_{t}$.

###### Proof.

The proof for continuous-time systems is similar to the proof in [30,
Proposition 3], while the discrete-time result follows from repeatedly
applying its definition in (2). ∎

Next, the following observation, which follows directly from Definitions 1, 3
and 4, relates the concept of decomposition functions to inclusion functions.

###### Proposition 4 (Decomposition-Based Inclusion Functions).

Given any upper and lower decomposition functions
$\overline{f}_{d},\underline{f}_{d}$ (or any decomposition function
$f_{d}=\overline{f}_{d}=\underline{f}_{d}$) for $f$,

$T^{f_{d}}(\mathbb{I}\mathcal{Z})\triangleq[\underline{f}_{d}(\underline{z},\overline{z}),\overline{f}_{d}(\overline{z},\underline{z})]$

satisfies $f(\mathbb{I}\mathcal{Z})\subset T^{f_{d}}(\mathbb{I}\mathcal{Z})$
with $\mathbb{I}\mathcal{Z}\triangleq[\underline{z},\overline{z}]$ (including
overloading; cf. Definition 1). Consequently,
$T^{f_{d}}(\mathbb{I}\mathcal{Z})$ is an inclusion function (that is based on
decomposition functions).

As noted earlier, (mixed-monotone) decomposition functions defined in
Definition 3 are not unique. Hence, a measure of their tightness is beneficial
for comparing these functions.

###### Definition 6 (Tightness of Decompositions).

[13, Definition 2] A decomposition function $f^{1}_{d}$ for system (1) is
tighter than decomposition function $f^{2}_{d}$, if  for all $z\leq\hat{z}$,

$\displaystyle f^{2}_{d}(z,\hat{z})\leq f^{1}_{d}(z,\hat{z})\ \text{and}\
f^{1}_{d}(\hat{z},{z})\leq f^{2}_{d}(\hat{z},{z}).$ (4)

Then,  $f^{O}_{d}$ is tight, i.e., it is the tightest possible decomposition
function for $f$, if (4) holds with $f^{1}_{d}={\color[rgb]{0,0,0}f^{O}_{d}}$
and any other decomposition function $f^{2}_{d}$.

Furthermore, we define the _metric/measure_ of tightness as the maximum
dimension-wise Hausdorff distance given by:

$\displaystyle
q(f(\mathbb{I}\mathcal{Z}),T^{f_{d}}(\mathbb{I}\mathcal{Z}))\triangleq\max_{i\in\mathbb{N}_{n_{x}}}\tilde{q}(f_{i}(\mathbb{I}\mathcal{Z}),T^{f_{d_{i}}}(\mathbb{I}\mathcal{Z})),$
(5)

where
$\tilde{q}({\color[rgb]{0,0,0}\mathbb{I}}\mathcal{X}_{1},{\color[rgb]{0,0,0}\mathbb{I}}\mathcal{X}_{2})\triangleq\max\\{|\underline{x}_{1}-\underline{x}_{2}|_{\infty},|\overline{x}_{1}-\overline{x}_{2}|_{\infty}\\}$
is the Hausdorff distance between two real intervals
${\color[rgb]{0,0,0}\mathbb{I}}\mathcal{X}_{1}=[\underline{x}_{1},\overline{x}_{1}]$
and
${\color[rgb]{0,0,0}\mathbb{I}}\mathcal{X}_{2}=[\underline{x}_{2},\overline{x}_{2}]$,
both in $\mathbb{IR}$ [19]. Moreover, the above tightest decomposition
function
$T_{O}^{f_{d}}(\mathbb{I}\mathcal{Z})\triangleq[f_{d}^{O}(\underline{z},\overline{z}),f_{d}^{O}(\overline{z},\underline{z})]$
satisfies
$q(f(\mathbb{I}\mathcal{Z}),T_{O}^{f_{d}}(\mathbb{I}\mathcal{Z}))=0$.

Further, by Proposition 4, we can obtain inclusion functions from existing
decomposition functions as shown in Propositions 5 and 6. Note that we
slightly modified the decomposition function in [14, Theorem 2] to make it
also applicable for continuous-time systems and for closed Jacobian intervals.

###### Proposition 5 ($T^{f_{d}}_{L}$ Inclusion Functions).

For any system in the form of (1), suppose that
$f:\mathbb{R}^{n_{z}}\to\mathbb{R}^{n_{x}}$ is differentiable and
$\frac{\partial f_{i}}{\partial z_{j}}(z)\in[a_{ij},b_{ij}],\forall
z\in{\color[rgb]{0,0,0}\mathbb{I}}\mathcal{Z}{\color[rgb]{0,0,0}\triangleq[\underline{z},\overline{z}]}\subseteq{\color[rgb]{0,0,0}\mathbb{I}}\mathbb{R}^{n_{z}}$.
Then, a discrete-time or continuous-time mixed-monotone decomposition function
$f^{L}_{d}=[f^{L}_{d,1}\dots f^{L}_{d,n_{x}}]$ with respect to $f(\cdot)$ and
its corresponding inclusion function
$T^{f_{d}}_{L}(\mathbb{I}\mathcal{Z})=[f^{L}_{d}({\color[rgb]{0,0,0}\underline{z}},{\color[rgb]{0,0,0}\overline{z}}),f^{L}_{d}({\color[rgb]{0,0,0}\overline{z}},{\color[rgb]{0,0,0}\underline{z}})]$
can be constructed as follows:

For all $i\in\mathbb{N}_{n_{x}}$ and $j\in\mathbb{N}_{n_{z}}$,

$\displaystyle
f^{L}_{d,i}({z},\hat{z})=f_{i}({\zeta})+(\alpha_{i}-\beta_{i})({z}-\hat{z}),$
(6)

with $\alpha_{i}=[\alpha_{i1},\dots,\alpha_{in_{z}}]$,
$\beta_{i}=[\beta_{i1},\dots,\beta_{in_{z}}]$ and
$\zeta=[\zeta_{1},\dots,\zeta_{n_{z}}]^{\top}$, where
$\alpha_{ij}=\begin{cases}0,&\text{Cases}\,1,3,4,5,\\\
|a_{ij}|,&\text{Case}\,2,\end{cases}$
$\beta_{ij}=\begin{cases}0,&\text{Cases}\,1,2,4,5,\\\
-|b_{ij}|,&\text{Case}\,3,\end{cases}$ and
$\zeta_{j}=\begin{cases}\underline{z}_{j},\ \text{Cases}\ 1,2,\\\
\overline{z}_{j},\ \text{Cases}\ 3,4,\\\ x_{j},\ \text{Case}\ 5,\end{cases}$
with Cases 1 through 4 for discrete-time systems and Cases 1 through 5 for
continuous-time systems, where the cases as defined as: Case $1:a_{ij}\geq 0$,
Case $2:a_{ij}\leq 0,b_{ji}\geq 0,|a_{ij}|\leq|b_{ij}|$, Case $3:a_{ij}\leq
0,b_{ij}\geq 0,|a_{ij}|\geq|b_{ij}|$, Case $4:b_{ij}\leq 0$ and Case 5: $j=i$.

###### Proposition 6 (Tight Decomposition Functions for Mixed-Monotone
Systems).

[12, Theorem 2],[13, Theorem 1] For any system in the form of (1)  and
$\mathbb{I}\mathcal{Z}\triangleq[\underline{z},\overline{z}]$, a tight
(optimal)  discrete-time or continuous-time mixed-monotone decomposition
function $f^{O}_{d}=[f^{O}_{d,1}\dots f^{O}_{d,n_{x}}]$ and its corresponding
tight inclusion function
$T^{f_{d}}_{O}(\mathbb{I}\mathcal{Z})\triangleq[f^{O}_{d}(\underline{z},\overline{z}),f^{O}_{d}(\overline{z},\underline{z})]$
(i.e.,
$T^{f}_{O}(\mathbb{I}\mathcal{Z})=T^{f_{d}}_{O}(\mathbb{I}\mathcal{Z})$) can
be  constructed as follows:

If (1) is a discrete-time  system, then $\forall i\in\mathbb{N}_{n_{x}}$:

$\displaystyle
f^{O}_{d,i}(z,\hat{z})=\begin{cases}\min\limits_{\zeta\in[z,\hat{z}]}f_{i}(\zeta)\quad\text{if}\
z\leq\hat{z},\\\
\max\limits_{\zeta\in[\hat{z},{z}]}f_{i}(\zeta)\quad\text{if}\
\hat{z}\leq{z}.\end{cases}$ (7)

Moreover, if (1) is a continuous-time system, then $\forall
i\in\mathbb{N}_{n_{x}}$:

$\displaystyle
f^{O}_{d,i}(z,\hat{z})=\begin{cases}\min\limits_{\zeta\in[z,\hat{z}],\zeta_{i}=x_{i}}f_{i}(\zeta)\quad\text{if}\
z\leq\hat{z},\\\\[-2.84544pt]
\max\limits_{\zeta\in[\hat{z},{z}],\zeta_{i}=x_{i}}f_{i}(\zeta)\quad\text{if}\
\hat{z}\leq{z}.\end{cases}$ (8)

###### Remark 1.

On the flip side, the discrete-time mixed-monotone decomposition functions can
also be directly used as inclusion functions; hence, the proposed
decomposition functions can also be relatively easily incorporated into
existing analysis, estimation and planning algorithms that are based on
_interval arithmetic_ , e.g., [18, 19, 20, 21, 22].

Although Proposition 6 provides theoretically tight decomposition functions,
it has some limitations in practice (see [30, Section III] for a detailed
discussion). For instance, exact closed-form solutions to the nonlinear
programs in (7) and (8) may not always be available.  With this in mind, we
define the notion of tractability of decomposition functions as follows.

###### Definition 7 (Tractable Decomposition Functions).

$\overline{f}_{d},\underline{f}_{d}$ are  computationally tractable/computable
one-sided decomposition functions for mapping $f$, if they can be  constructed
in closed-form , i.e., with a finite number of elementary operations and
without differentiation nor an iterative procedure.

Thus, a corollary of Proposition 6 is as follows.

###### Corollary 2 (Tight  and Tractable Decomposition  Functions for JSS
Vector Fields).

Suppose $f(\cdot)$ is continuously differentiable and Jacobian sign-stable
(JSS) [14], i.e.,  $\forall i\in\mathbb{N}_{n_{x}},\forall
j\in\mathbb{N}_{n_{z}},J^{f}_{ij}(z)\triangleq\frac{\partial f_{i}}{\partial
z_{j}}(z)\geq 0,\forall z\in\mathbb{I}\mathcal{Z}$ or
$J^{f}_{ij}(z)\triangleq\frac{\partial f_{i}}{\partial z_{j}}(z)\leq 0,\forall
z\in\mathbb{I}\mathcal{Z}$. Then, the following statements hold :

  1. i)

$\forall i\in\mathbb{N}_{n_{x}},\forall j\in\mathbb{N}_{n_{z}}$,
$f_{i}(\cdot)$ in either monotonically non-decreasing or monotonically non-
increasing in its $j$-th argument $z_{j}$, over the entire domain
$\mathbb{I}\mathcal{Z}$ ;

  2. ii)

The optimization programs in (7) and (8) can be  tractably and exactly solved
by enumerating $f_{i}(\cdot)$ at the vertices of  $\mathbb{I}\mathcal{Z}$
(with fixed $\zeta_{i}=x_{i}$ for continuous-time systems) and choosing the
corresponding optima.

####  II-B2  A Novel Class of Mixed-Monotone Systems

To describe our modeling framework, we formally define a novel class of
nonsmooth systems, which we show to include a wide range of nonlinearities and
prove in Section III to be mixed-monotone.

###### Definition 8 (Either-Sided Locally Lipschitz Systems).

System (1) is _either-sided locally Lipschitz continuous_ (ELLC)  at
$z\in\mathcal{Z}$, if there exists an open neighborhood
$\mathcal{N}_{z}\subset\mathcal{Z}$ of $z$ , and for all
$i\in\mathbb{N}_{n_{x}}$, there exist vectors
$\kappa_{i}\triangleq[\kappa_{i1}\dots\kappa_{in_{z}}]\in\mathbb{R}^{n_{z}}$
with non-zero elements (i.e., $\kappa_{ij}\neq 0,\forall
j\in\mathbb{N}_{n_{z}}$) and constants $\rho_{i}\in\mathbb{R}$ such that
$\forall z^{\prime},z^{\prime\prime}\in\mathcal{N}_{z}$,

$\displaystyle{\color[rgb]{0,0,0}{\color[rgb]{0,0,0}\langle}\kappa_{i}(f_{i}(z^{\prime})-f_{i}(z^{\prime\prime})),z^{\prime}-z^{\prime\prime}{\color[rgb]{0,0,0}\rangle}\leq\rho_{i}\|z^{\prime}-z^{\prime\prime}\|^{2}_{2},\
\forall i\in\mathbb{N}_{n_{x}},}$ (9)

where  $\langle\cdot,\cdot\rangle$ denotes the inner product operator.
Further, if (9) holds $\forall i\in\mathbb{N}_{n_{x}}$, we also call $f$ an
ELLC function on $\mathcal{Z}$.

Note that (9) holds for a very broad range of nonlinear systems. Particularly,
ELLC systems reduce to one-sided locally Lipschitz systems111System (1) with
$n_{x}=n_{z}$ is one-sided locally Lipschitz continuous  if for all
$z\in\mathcal{Z}$, there exist an open neighborhood
$\mathcal{N}_{z}\subset\mathcal{Z}$ of $z$ and a ${\rho}\in\mathbb{R}$ such
that
$(f(z^{\prime})-f(z^{\prime\prime}))^{\top}(z^{\prime}-z^{\prime\prime})\leq{\rho}\|z^{\prime}-z^{\prime\prime}\|^{2}_{2},\forall
z^{\prime},z^{\prime\prime}\in\mathcal{N}_{z}$.,  when $n_{z}=n_{x}$ and
$\kappa_{ii}=\kappa_{jj}\triangleq\kappa>0$ and $\kappa_{ij}=0$, $\forall
i\neq j$, with Lipschitz constant
$\frac{1}{\kappa}\sum_{i=1}^{n_{x}}\rho_{i}$. Moreover, the class of ELLC
systems is a superset of the class of locally Lipschitz continuous (LLC)
systems222System (1) is  LLC if for all $z\in\mathcal{Z}$, there exist an open
neighborhood $\mathcal{N}_{z}$ of $z$ and a ${\rho}\geq 0$ such that
$\|f(z^{\prime})-f(z^{\prime\prime})\|_{2}\leq\rho\|z^{\prime}-z^{\prime\prime}\|_{2}$,
$\forall z^{\prime},z^{\prime\prime}\in\mathcal{N}_{z}$., as shown next.

###### Proposition 7.

An LLC system is also an ELLC system.

###### Proof.

Suppose a system is LLC, then each $f_{i}$ is also LLC with some
$\tilde{\rho}_{i}$. Next, by applying the Cauchy-Schwartz inequality and with
each $f_{i}$ being LLC and any $\kappa_{i}$ with non-zero elements, we have
${\color[rgb]{0,0,0}\langle}\kappa_{i}(f_{i}(z^{\prime})-f_{i}(z)),z^{\prime}-z{\color[rgb]{0,0,0}\rangle}\leq|\kappa_{i}||f_{i}(z^{\prime})-f_{i}(z)|\|z^{\prime}-z\|_{2}\leq\rho_{i}\|z^{\prime}-z\|^{2}_{2}$,
where $\rho_{i}\triangleq|\kappa_{i}|\tilde{\rho}_{i}$. ∎

Next, we  review a notion of generalized gradients used in _nonsmooth_
analysis in systems and control theory [31] when  vector fields of the systems
are not necessarily differentiable.

###### Definition 9 (Clarke Generalized  Directional  Derivatives).

[32, Chapter II] Given a function
$f:\mathcal{Z}\subseteq\mathbb{R}^{n_{z}}\to\mathbb{R}$,

$\displaystyle\begin{array}[]{ll}{f_{C}^{\uparrow}}(z,v)\triangleq\limsup\limits_{t\to
z,\lambda\downarrow 0}\frac{f(t+\lambda
v)-f(t)}{\lambda}=\sup\limits_{\xi\in{\partial}_{C}f(z)}\xi^{\top}v,\\\
{f_{C}^{\downarrow}}(z,v)\triangleq\liminf\limits_{t\to z,\lambda\downarrow
0}\frac{f(t+\lambda
v)-f(t)}{\lambda}=\inf\limits_{\xi\in{\partial}_{C}f(z)}\xi^{\top}v,\end{array}$
(12)

are the (generalized) Clarke upper and lower directional derivatives
/gradients of $f(\cdot)$ at $z\in\mathcal{Z}$ in the direction
$v\in\mathbb{R}^{n_{z}}$, respectively, where the set

$\displaystyle{\partial}_{C}f(z)\triangleq\\{\xi\in\mathbb{R}^{n_{z}}|{f_{C}^{\downarrow}}(z,v)\leq\xi^{\top}v\leq{f_{C}^{\uparrow}}(z,v),\forall
v\in\mathbb{R}^{n_{z}}\\},$

is the Clarke sub-differential (set) of $f(\cdot)$ at $z\in\mathcal{Z}$.

Note that, by definition, ${f_{C}^{\uparrow}}(z,v)$ and
${f_{C}^{\downarrow}}(z,v)$ are the _upper_ and _lower support functions_ of
the set ${\partial}_{C}f(z)$.  Further, as shown in [32,  Appendix I],
$\partial_{C}f(z)$ for an LLC system is nonempty, convex and compact for all
$z\in\mathcal{Z}$, and consequently, at each $z\in\mathcal{Z}$,  the Clarke
directional derivatives are bounded in each direction $v$. However, this does
not hold in general for ELLC systems.  Nonetheless, by the following
proposition , the Clarke directional derivatives  in some specific directions
are bounded from above _or_ from below.

###### Proposition 8.

Suppose $f:\mathcal{Z}\subseteq\mathbb{R}^{n_{z}}\to\mathbb{R}$ is ELLC on
$\mathcal{Z}$ and let  $e_{j},\forall j\in\mathbb{N}_{n_{z}}$ denote the
standard unit vector /basis in the $j$-th coordinate direction. Then,
$f_{C}^{\uparrow}(z,e_{j})$ is bounded from above or
$f_{C}^{\downarrow}(z,e_{j})$ is bounded from below.

###### Proof.

Setting $z^{\prime}=z+\lambda e_{j}$ in (9), where $\lambda\in\mathbb{R}$ is
sufficiently small, we obtain $\kappa_{ij}(f_{i}(z+\lambda
e_{j})-f_{i}(z))\lambda\leq\rho_{i}\lambda^{2}$. Then, dividing both sides by
$\kappa_{ij}\neq 0$ and taking the $\limsup$ if $\kappa_{ij}>0$, or the
$\liminf$ if $\kappa_{ij}<0$, from both sides when $\lambda\to 0$, imply that
$f_{C}^{\uparrow}(z,e_{j})$ is bounded from above or
$f_{C}^{\downarrow}(z,e_{j})$ is bounded from below, respectively (cf. (12)).
∎

Then, by Proposition 8, we can easily show the following.

###### Corollary 3.

Consider the mapping $f=[f_{1}\dots
f_{n_{x}}]^{\top}:\mathcal{Z}\subset\mathbb{R}^{n_{z}}\to\mathbb{R}^{n_{x}}$,
where $\forall i\in\mathbb{N}_{n_{x}}$, $f_{i}$ is ELLC on $\mathcal{Z}$. We
define upper and lower Clarke Jacobian matrices of $f$ at $z\in\mathcal{Z}$,
$\overline{J}^{f}_{C}(z)=[(\overline{J}^{f}_{C}(z))_{ij}]$ and
$\underline{J}^{f}_{C}(z)=[(\underline{J}^{f}_{C}(z))_{ij}]$, with the upper
and lower partial Clarke derivatives at point $z\in\mathcal{Z}$ for each
$i\in{\color[rgb]{0,0,0}\mathbb{N}_{n_{x}}},j\in{\color[rgb]{0,0,0}\mathbb{N}_{n_{z}}}$,
respectively, defined as:

$\displaystyle(\overline{J}^{f}_{C}(z))_{ij}\triangleq
f^{\uparrow}_{i,C}(z,e_{j}),\ (\underline{J}^{f}_{C}(z))_{ij}\triangleq
f^{\downarrow}_{i,C}(z,e_{j}).$ (13)

Then, $\forall i\in{\color[rgb]{0,0,0}\mathbb{N}_{n_{x}}},\forall
j\in{\color[rgb]{0,0,0}\mathbb{N}_{n_{z}}}$, $\overline{J}^{f}_{C}(z)_{ij}$ is
bounded from above or/and $\underline{J}^{f}_{C}(z)_{ij}$ is bounded  from
below,  i.e., we cannot simultaneously have
$(\overline{J}^{f}_{C}(z))_{ij}=\infty$ and
$(\underline{J}^{f}_{C}(z))_{ij}=-\infty$.

Next, we  provide a slight modification of the results in [32, Proposition
1.12], which plays an important role in our main results later and its proof
goes precisely  along the lines of the proof of [32, Proposition 1.12].

###### Proposition 9.

Let $f:\mathcal{Z}\subseteq\mathbb{R}^{n_{z}}\to\mathbb{R}^{n_{x}}$ be
decomposable into $f=f^{1}+f^{2}$, where
$f^{1},f^{2}:\mathcal{Z}\subseteq\mathbb{R}^{n_{z}}\to\mathbb{R}$ and
$f_{i},f^{1}_{i},f^{2}_{i}$, $\forall
i\in{\color[rgb]{0,0,0}\mathbb{N}_{n_{x}}}$, are ELLC. Then, $\forall
z\in\mathcal{Z},\forall i\in{\color[rgb]{0,0,0}\mathbb{N}_{n_{x}}},\forall
j\in{\color[rgb]{0,0,0}\mathbb{N}_{n_{z}}}$:

$\displaystyle(\overline{J}^{f}_{C}(z))_{ij}$
$\displaystyle\leq(\overline{J}^{f^{1}}_{C}(z))_{ij}+(\overline{J}^{f^{2}}_{C}(z))_{ij},$
$\displaystyle(\underline{J}^{f}_{C}(z))_{i,j}$
$\displaystyle\geq(\underline{J}^{f^{1}}_{C}(z))_{ij}+(\underline{J}^{f^{2}}_{C}(z))_{ij}.$

###### Proof.

The results follow from (13) and the facts that $f^{\uparrow}_{C}(z,v)\leq
f^{1\uparrow}_{C}(z,v)+f^{2\uparrow}_{C}(z,v)$ [32, Chapter II, Proposition
1.12] and $f^{\downarrow}_{C}(z,v)=-f^{\uparrow}_{C}(z,-v)$ [32, Chapter II,
Proposition 1.7], $\forall z\in\mathcal{Z},\forall v\in\mathbb{R}^{n_{z}}$. ∎

It is worth noticing that when $f(\cdot)$ is differentiable, then $\nabla
f(x)\in\partial_{C}f(x)$ and if $f(\cdot)$ is continuously differentiable or
strictly differentiable, then $\partial_{C}f(x)=\\{\nabla f(x)\\}$. Now, we
are ready to explain the notion of _Clarke Jacobian sign-stability_ through
the following definition, which is a generalization of Jacobian sign-stability
in [14] ( cf. Corollary 2).

###### Definition 10 (Clarke Jacobian Sign-Stability).

A mapping $f:\mathcal{Z}\subset\mathbb{R}^{n_{z}}\to\mathbb{R}^{n_{x}}$ is
called Clarke Jacobian sign-stable (CJSS) over $\mathcal{Z}$, if $\forall
z\in\mathcal{Z},\forall i\in{\color[rgb]{0,0,0}\mathbb{N}_{n_{x}}},\forall
j\in{\color[rgb]{0,0,0}\mathbb{N}_{n_{z}}}$,

$\displaystyle(\overline{J}^{f}_{C}(z))_{ij}\leq 0\ \lor\
(\underline{J}^{f}_{C}(z))_{ij}\geq 0.$ (14)

Finally,  we present an extension of Corollary 2, which we will apply later in
our derivations.

###### Proposition 10.

Suppose $f$ is ELLC and CJSS over $\mathbb{I}\mathcal{Z}$. Then, $\forall
i\in\mathbb{N}_{n_{x}},\forall j\in\mathbb{N}_{n_{z}}$, $f_{i}(\cdot)$ in
either monotonically non-decreasing or monotonically non-increasing in its
$j$-th argument $z_{j}$, over the entire domain $\mathbb{I}\mathcal{Z}$ and
consequently, the optima of (7) and (8) are attained at some vertices of
$\mathbb{I}\mathcal{Z}$.

###### Proof.

For $i\in\mathbb{N}_{n_{x}}$, consider any arbitrary $z^{1}$ from the interior
of $\mathbb{I}\mathcal{Z}$ and construct $z^{2}=z^{1}+\lambda
e_{j}\in\mathbb{I}\mathcal{Z}$, for some small enough $\lambda>0$. Then, by an
identifical proof to [32, Chapter II, Theorem 1.3], there exists a
$z_{\theta}\in\mathbb{I}\mathcal{Z}$ on the connecting line between $z^{1}$
and $z^{2}$, such that $\lambda
f^{\downarrow}_{i,C}(z_{\theta},e_{j})=f^{\downarrow}_{i,C}(z_{\theta},\lambda
e_{j})\leq f_{i}(z^{2})-f_{i}(z^{1})\leq f^{\uparrow}_{i,C}(z_{\theta},\lambda
e_{j})=\lambda f^{\uparrow}_{i,C}(z_{\theta},e_{j})$, where the equalities
hold by [32, Chapter II, Proposition 1.5]. Using this and the CJSS assumption
(cf. (14)), $f_{i}(z^{2})\leq f_{i}(z_{2})$ if $f^{\uparrow}_{i,C}(z,\lambda
e_{j})=(\overline{J}^{f}_{C}(z))_{ij}\leq 0$ holds for all
$z\in\mathbb{I}\mathcal{Z}$ (including $z=z_{\theta}$), and similarly,
$f_{i}(z^{2})\geq f_{i}(z_{2})$ if $f^{\downarrow}_{i,C}(z,\lambda
e_{j})=(\underline{J}^{f}_{C}(z))_{ij}\geq 0$ holds for all
$z\in\mathbb{I}\mathcal{Z}$ (including $z=z_{\theta}$). Hence, by moving along
the coordinate directions, one can always increase or decrease each of the
$f_{i}$’s. Using a similar argument, this result also holds for $\lambda<0$.
So, the optimum for each $f_{i}$ is attained at some vertices of the interval
$\mathbb{I}\mathcal{Z}$.  ∎

###  II-C Modeling Framework and Problem Statement

We consider constrained dynamical systems of the form:

$\displaystyle
x_{t}^{+}={\color[rgb]{0,0,0}\tilde{f}(x_{t},u_{t},w_{t})\triangleq
f(z_{t}),}\quad{\color[rgb]{0,0,0}\mu(x_{t},u_{t})\triangleq\nu(x_{t})\in\mathcal{Y}_{t},}$
(15)

where $x_{t}^{+}=\dot{x}_{t}$ if (15) is a continuous-time  system and
$x_{t}^{+}=x_{t+1}$ if (15) is a discrete-time system.
$x_{t}\in\mathcal{X}\subset\mathbb{R}^{n_{x}}$ is the state,
$u_{t}\in\mathcal{U}\subset\mathbb{R}^{n_{x}}$ is the known input and
$w_{t}\in\mathcal{W}\subseteq\mathbb{I}\mathcal{W}\triangleq[\underline{w},\overline{w}]\in\mathbb{IR}^{n_{w}}$
is the augmentation of all exogenous inputs, e.g., bounded disturbance/noise
and internal uncertainties such as uncertain parameters, with known bounds
$\underline{w},\overline{w}$, while
$f:\mathcal{X}\times\mathcal{U}\times\mathcal{W}\to\mathbb{R}^{n_{x}}$ and
$\mu:\mathcal{X}\times\mathcal{U}\to\mathbb{R}^{n_{\color[rgb]{0,0,0}y}}$ are
the nonlinear vector field and the observation/constraint mapping,
respectively,  while $\mathcal{Y}_{t}$ is the known or measured time-varying
constraint/observation set. The mapping $\mu(\cdot)$ and the set
$\mathcal{Y}_{t}\subseteq[\underline{y}_{t},\overline{y}_{t}]$ describe system
constraints that can represent prior or additional knowledge about the system
states, e.g., sensor observations or measurements with bounded noise, known
state constraints, manufactured constraints from modeling redundancy [27] (cf.
Section VI-B for an example), etc. For ease of exposition, we further define
$f:\mathcal{Z}\triangleq\mathcal{X}\times\mathcal{W}\subset\mathbb{R}^{n_{z}}\to\mathbb{R}^{n_{x}}$
and $\nu:\mathcal{X}\subset\mathbb{R}^{n_{x}}\to\mathbb{R}^{n_{y}}$ as in (15)
that is implicitly dependent on $u_{t}$ with the augmented state
$z_{t}\triangleq[x_{t}^{\top}\ w_{t}^{\top}]^{\top}\in\mathcal{Z}$. Further,
we assume the following.

###### Assumption 1.

The mappings $f$ and  $\nu$ are either-sided locally Lipschitz (ELLC; cf.
Definition 8).

###### Assumption 2.

For the mappings $f$ and $\nu$, there exist known bounds on their Clarke
Jacobian matrices,
$\overline{J}^{f}_{C},\underline{J}^{f}_{C}\in\\{\mathbb{R}\cup\pm\infty\\}^{n_{x}\times
n_{z}}$ and
$\overline{J}^{{\color[rgb]{0,0,0}\nu}}_{C},\underline{J}^{{\color[rgb]{0,0,0}\nu}}_{C}\in\\{\mathbb{R}\cup\pm\infty\\}^{\color[rgb]{0,0,0}n_{y}\times
n_{x}}$, that satisfy

$\displaystyle\begin{array}[]{c}((J^{f}_{C}(z))_{ij}\leq({\overline{J}^{f}_{C}})_{ij}<\infty)\vee((J^{f}_{C}(z))_{ij}\geq({\underline{J}^{f}_{C}})_{ij}>-\infty),\\\
((J^{{\color[rgb]{0,0,0}\nu}}_{C}({\color[rgb]{0,0,0}x}))_{pq}\leq({\overline{J}^{{\color[rgb]{0,0,0}\nu}}_{C}})_{pq}<\infty)\vee((J^{{\color[rgb]{0,0,0}\nu}}_{C}({\color[rgb]{0,0,0}x}))_{pq}\geq({\underline{J}^{{\color[rgb]{0,0,0}\nu}}_{C}})_{pq}>-\infty),\end{array}$

$\forall z\in\mathbb{I}\mathcal{Z}$,
$\forall{\color[rgb]{0,0,0}x\in\mathbb{I}\mathcal{X}}$,
$\forall{\color[rgb]{0,0,0}i,q}\in\mathbb{N}_{n_{x}}$, $\forall
j\in\mathbb{N}_{n_{z}}$, $\forall p\in\mathbb{N}_{n_{\color[rgb]{0,0,0}y}}$.

Under the above modeling framework and assumptions, this paper seeks to find
tight and tractable (i.e., computable) remainder-form upper and lower
decomposition functions and their induced inclusion functions (cf. Definitions
1, 2, 4 and  Proposition 4) as well as  to develop set inversion algorithms
based on (mixed-monotone) decomposition functions.

###### Problem 1 (Decomposition Functions).

Suppose Assumptions 1 and 2 hold. Construct and quantify the tightness (via
the metric (5)) of remainder-form decomposition functions, by solving the
following subproblems:

  1. 1.1

Given an ELLC vector field
$f:\mathcal{Z}\subset\mathcal{R}^{n_{z}}\to\mathbb{R}^{n_{x}}$, construct a
tractable family of mixed-monotone remainder-form (i.e., additive)
decomposition functions for $f(\cdot)$.

  2. 1.2

Derive lower and upper bounds for the _tightness_ (quantified via (5)) of the
family of remainder-form decomposition functions obtained in 1.1.

  3. 1.3

Find the tightest decomposition function(s) among the family of remainder-form
decomposition functions obtained in 1.1 and compare them with the
decomposition function in [14] (cf. Proposition 5), and natural, centered and
mixed-centered natural inclusions (cf. Proposition 2).

###### Problem 2 (Set Inversion Algorithm).

Suppose Assumptions 1 and 2 hold. Given a prior/propagated interval
$\mathbb{I}{\color[rgb]{0,0,0}\mathcal{X}}^{p}_{\color[rgb]{0,0,0}t}\in\mathbb{IR}^{n_{\color[rgb]{0,0,0}x}}$,
a constraint/observation function and set,
$\mu({\color[rgb]{0,0,0}x}_{t},u_{t})$ and
$\mathcal{Y}_{t}\subseteq[\underline{y}_{t},\overline{y}_{t}]$ with known
$u_{t}$, develop an algorithm to find an interval superset of all
${\color[rgb]{0,0,0}x}_{t}$ that are compatible with $\mu(\cdot)$,
$[\underline{y}_{t},\overline{y}_{t}]$ and
$\mathbb{I}{\color[rgb]{0,0,0}\mathcal{X}}^{p}_{\color[rgb]{0,0,0}t}$, i.e.,
to find the updated/refined interval
$\mathbb{I}{\color[rgb]{0,0,0}\mathcal{X}}^{u}_{t}$ such that

$\displaystyle\\{{\color[rgb]{0,0,0}x}\in\mathbb{I}{\color[rgb]{0,0,0}\mathcal{X}}^{p}_{\color[rgb]{0,0,0}t}\mid\mu({\color[rgb]{0,0,0}x},u_{t})\in[\underline{y}_{t},\overline{y}_{t}]\\}\subseteq\mathbb{I}{\color[rgb]{0,0,0}\mathcal{X}}^{u}_{t}\subseteq\mathbb{I}{\color[rgb]{0,0,0}\mathcal{X}}^{p}_{\color[rgb]{0,0,0}t}.$
(16)

In the context of constrained reachability analysis and interval observer
design (cf. Section V), the solution of the generalized embedding system (cf.
Definition 5) based on the decomposition functions obtained from solving
Problem 1 provides the unconstrained reachable set (or propagated set),
$\mathbb{I}{\color[rgb]{0,0,0}\mathcal{X}}_{t}^{p}$, while the set inversion
algorithm in Problem 2 finds the constrained reachable set (or updated set)
$\mathbb{I}{\color[rgb]{0,0,0}\mathcal{X}}_{t}^{u}$.

##  III Main Results

To address the aforementioned problems, we first describe  our proposed
construction approach to find a tractable family of mixed-monotone remainder-
form decomposition functions. Then, by characterizing their tightness, we can
determine the tightest decomposition function among the proposed family.
Further, we present a novel set inversion algorithm that serves as an
alternative and improves on existing approaches, e.g., SIVIA in [18, Chapter
3]) and $\mathcal{I}_{G}$ in [23, Algorithm 1]).

###  III-A Remainder-Form Decomposition Functions

To solve Problem 1.1, we provide a constructive procedure for computing the
family of remainder-form decomposition functions in a tractable manner (i.e.,
in closed-form). Intuitively, our approach is based on the idea of decomposing
each ELLC function $f$ into the remainder/additive form, i.e., $f=g+h$, such
that $g$ is a CJSS function (cf. Definition 10; so that Proposition 10
applies) by “shifting” the Clarke directional gradients of $f$ and accounting
for the “error” using $h$. Since there are several “shift” directions, we
obtain a family of decomposition functions. Note that the construction below
is to be independently performed for each dimension of the ELLC function $f$;
hence, without loss of generality, we only consider a scalar ELLC function
$f_{i}$ throughout this subsection.

###### Theorem 1 (A Family of Remainder-Form Decomposition Functions).

Consider an ELLC vector field
$f_{i}:\mathcal{Z}\subset\mathbb{R}^{n_{z}}\to\mathbb{R}$ and suppose that
Assumptions 1 and 2 hold. Then, $f_{i}(\cdot)$ admits a family of mixed-
monotone remainder-form decomposition functions denoted as
$\\{f_{d,i}(z,\hat{z};\mathbf{m},h(\cdot))\\}_{\mathbf{m}\in\mathbf{M}_{i},h(\cdot)\in\mathcal{H}_{\mathbf{M}_{i}}}$
that is parameterized by a supporting vector $\mathbf{m}\in\mathbf{M}_{i}$ and
an ELLC remainder function $h(\cdot)\in\mathcal{H}_{\mathbf{M}_{i}}$, where

$\displaystyle
f_{d,i}(z,\hat{z};\mathbf{m},h(\cdot))=h(\zeta_{\mathbf{m}}(\hat{z},{z}))+f_{i}(\zeta_{\mathbf{m}}(z,\hat{z}))-h(\zeta_{\mathbf{m}}(z,\hat{z})),$
(17)

with
$\zeta_{\mathbf{m}}(z,\hat{z})=[\zeta_{\mathbf{m},1}(z,\hat{z}),\dots,\zeta_{\mathbf{m},n_{z}}(z,\hat{z})]^{\top}$,
where $\forall j\in\mathbb{N}_{n_{z}}$,

$\displaystyle\begin{array}[]{rl}&\zeta_{\mathbf{m},j}(z,\hat{z})=\begin{cases}\hat{z}_{j},\
\text{if}\ \mathbf{m}_{j}\geq\max(({\overline{J}^{f}_{C}})_{ij},{0}),\\\
z_{j},\ \text{if}\
\mathbf{m}_{j}\leq\min(({\underline{J}^{f}_{C}})_{ij},{0}),\end{cases}\end{array}$
(19)

and the set of supporting vectors $\mathbf{M}_{i}$ is defined as:

$\displaystyle\begin{array}[]{rl}\mathbf{M}_{i}\triangleq&\\{\mathbf{m}\in\mathbb{R}^{n_{z}}\
|\ \mathbf{m}_{j}\geq\max(({\overline{J}^{f}_{C}})_{ij},{0})\ \lor\\\ &\ \
\mathbf{m}_{j}\leq\min(({\underline{J}^{f}_{C}})_{ij},{0}),\forall
j\in\mathbb{N}_{n_{z}}\\}\end{array}$ (22)

if (15) is a discrete-time system and with

$\displaystyle\begin{array}[]{rl}&\zeta_{\mathbf{m},j}(z,\hat{z})=\begin{cases}\hat{z}_{j},\
\text{if}\ \mathbf{m}_{j}\geq\max(({\overline{J}^{f}_{C}})_{ij},{0}),\\\
z_{j},\ \text{if}\
\mathbf{m}_{j}\leq\min(({\underline{J}^{f}_{C}})_{ij},{0}),\\\ x_{i},\
\text{if}\ j=i,\end{cases}\end{array}$ (24)
$\displaystyle\begin{array}[]{l}\mathbf{M}_{i}\triangleq\\{\mathbf{m}\in\mathbb{R}^{n_{z}}\
|\ \mathbf{m}_{j}\geq\max(({\overline{J}^{f}_{C}})_{ij},{0})\ \lor\\\ \quad\
\mathbf{m}_{j}\leq\min(({\underline{J}^{f}_{C}})_{ij},{0}),\forall
j\in\mathbb{N}_{n_{z}},j\neq i,\mathbf{m}_{i}=0\\}\end{array}$ (27)

if (15) is a continuous-time system, while the set of remainder functions
$\mathcal{H}_{\mathbf{M}_{i}}$ is the family of all ELLC remainder functions
whose Clarke sub-differential set over $\mathcal{Z}$ (cf. Definition 9) is a
subset of ${\mathbf{M}_{i}}$ given by:

$\displaystyle\mathcal{H}_{\mathbf{M}_{i}}\triangleq\\{{h}:\mathcal{Z}\to\mathbb{R}\
|\
[\underline{J}^{{h}}_{C}(z),\overline{J}^{{h}}_{C}(z)]\subseteq\mathbf{M}_{i},\forall
z\in\mathcal{Z}\\}.$ (28)

The resulting family of decomposition-based inclusion functions is given by:

$\displaystyle
T^{f_{d,i}}_{\underline{\mathbf{m}},\overline{\mathbf{m}}}(\mathbb{I}\mathcal{Z})\triangleq[{f}_{d,i}(\underline{z},\overline{z};\underline{\mathbf{m}},\underline{h}(\cdot)),{f}_{d,i}(\overline{z},\underline{z};\overline{\mathbf{m}},\overline{h}(\cdot))],$
(29)

for all $\underline{\mathbf{m}},\overline{\mathbf{m}}\in\mathbf{M}_{i}$, and
the corresponding
$\underline{h}(\cdot),\overline{h}(\cdot)\in\mathcal{H}_{\mathbf{M}_{i}}$.

It may be worth noting that the small difference in the definitions of
$\mathbf{M}$ for discrete-time versus continuous-time systems in Theorem 1
originates from the subtle difference between the definitions of decomposition
functions for discrete-time and continuous-time systems (cf. Definition 3).
To prove the above theorem, we first prove the following lemma.

###### Lemma 1.

Any remainder function $h(\cdot)\in\mathcal{H}_{\mathbf{M}_{i}}$ (cf. (28)) is
CJSS and the function $g_{i}(\cdot)\triangleq f_{i}(\cdot)-h(\cdot)$ is also
CJSS. Moreover, the pair $(g_{i}(\cdot),-h(\cdot))$ is aligned, i.e., $\forall
j\in\mathbb{N}_{n_{z}}$,

  1. (i)

$(\overline{J}^{g_{i}}_{C}(z))_{j}\leq 0$ if and only if
$(\underline{J}^{h}_{C}(z))_{j}\geq 0$, or

  2. (ii)

$(\underline{J}^{g_{i}}_{C}(z))_{j}\geq 0$ if and only if
$(\overline{J}^{h}_{C}(z))_{j}\leq 0$.

###### Proof.

Since $\mathbf{m}\in\mathbf{M}_{i}$, by construction of $\mathbf{M}_{i}$ (cf.
(22), (27)), $h(\cdot)\in\mathcal{H}_{\mathbf{M}_{i}}$ is a CJSS function.
Next, by applying Proposition 9 to $g_{i}(\cdot)\triangleq
f_{i}(\cdot)-h(\cdot)$, we know that $\forall z\in\mathcal{Z},\forall
j\in\mathbb{N}_{n_{z}}$:

$\displaystyle(\underline{J}^{g_{i}}_{C}(z))_{j}$
$\displaystyle\leq(\underline{J}^{f}_{C}(z))_{ij}-(\underline{J}^{h}_{C}(z))_{j},$
(30) $\displaystyle(\overline{J}^{g_{i}}_{C}(z))_{j}$
$\displaystyle\geq(\overline{J}^{f}_{C}(z))_{ij}-(\overline{J}^{h}_{C}(z))_{j}.$
(31)

Then, since $h(\cdot)\in\mathcal{H}_{\mathbf{M}_{i}}$, according to (19)–(28),
we consider the following two cases (the case with $\mathbf{m}_{i}=0$ for
continuous-time systems is trivial):

  1. (i)

$\forall
z\in\mathcal{Z},(\underline{J}^{h}_{C}(z))_{j}\geq\max(({\overline{J}^{f}_{C}})_{ij},{0})\geq
0$. Then, from (30), we obtain $\forall
z\in\mathcal{Z},(\underline{J}^{g_{i}}_{C}(z))_{j}\leq
0\Rightarrow{(\underline{J}^{g_{i}}_{C}})_{j}\leq 0$.

  2. (ii)

$\forall
z\in\mathcal{Z},(\overline{J}^{h}_{C}(z))_{j}\leq\min(({\underline{J}^{f}_{C}})_{ij},{0})\leq
0$. Then, from (31), we have $\forall
z\in\mathcal{Z},(\overline{J}^{g_{i}}_{C}(z))_{j}\geq
0\Rightarrow({\overline{J}^{g_{i}}_{C}})_{j}\geq 0$.

The reverse can be similarly deduced. Finally, since
$({\overline{J}^{g_{i}}_{C}})_{j}\geq 0$ or
${(\underline{J}^{g_{i}}_{C}})_{j}\leq 0$ holds, $g_{i}(\cdot)$ is CJSS by
Definition 10. ∎

###### Remark 2.

Since the pair $(g_{i},-h)$ is aligned and $f_{i}=g_{i}+h$, the proposed
remainder-form decomposition function can also be viewed as the decomposition
of $f_{i}$ into a _difference of monotone functions_ , which is similar in
spirit to difference of convex functions in DC programming, e.g., [25, 28].

###### Proof of Theorem 1.

Armed by Lemma 1, we now prove that (17) is mixed-monotone (cf. Definition 3).
Having defined $g_{i}(\cdot)\triangleq f_{i}(\cdot)-h(\cdot)$, (17) can be
rewritten as

$\displaystyle
f_{d,i}(z,\hat{z};\mathbf{m},h)=h(\zeta_{\mathbf{m}}(\hat{z},{z}))+g_{i}(\zeta_{\mathbf{m}}({z},\hat{z})).$
(32)

First, it directly follows from (19) that $\zeta_{\mathbf{m}}(z,{z})=z$ and
$f_{d,i}(z,{z};\mathbf{m},h(\cdot))=f_{i}(z)$. Hence, it remains to show that
$f_{d,i}(z,\hat{z};\mathbf{m},h(\cdot))$ is non-decreasing in $z$ and non-
increasing in $\hat{z}$. To do this, consider
$z,\tilde{z},\hat{z}\in\mathcal{Z}$, where $\tilde{z}\geq z$. Let
$j_{0}\in\mathbb{N}_{n_{z}}$, and suppose case (i) in the proof of Lemma 1
holds for dimension $j_{0}$. Then, by the first case in (19),

$\displaystyle\begin{array}[]{c}\zeta_{\mathbf{m},j_{0}}(\hat{z},\tilde{z})=\tilde{z}_{j_{0}}\geq{z}_{j_{0}}=\zeta_{\mathbf{m},j_{0}}(\hat{z},z),\\\
\zeta_{\mathbf{m},j_{0}}(\tilde{z},\hat{z})=\hat{z}_{j_{0}}=\zeta_{\mathbf{m},j_{0}}(z,\hat{z}).\end{array}$
(35)

Next, we define ${z}^{1}\in\mathcal{Z}$ as follows:
${z}^{1}_{j_{0}}=\tilde{z}_{j_{0}}$ and ${z}^{1}_{j}=z_{j},\forall j\neq
j_{0}$. Thus, $z^{1}\geq z$, and by (35),
$\zeta_{\mathbf{m}}(z^{1},\hat{z})=\zeta_{\mathbf{m}}(z,\hat{z})$,
$\zeta_{\mathbf{m},j}(\hat{z},z^{1})=\zeta_{\mathbf{m},j}(\hat{z},z),\forall
j\neq j_{0}$ and
$\zeta_{\mathbf{m},j_{0}}(\hat{z},z^{1})\geq\zeta_{\mathbf{m},j_{0}}(\hat{z},z)$.
Moreover, by case (i) in the proof of Lemma 1 and Proposition 10, $h(\cdot)$
is non-decreasing in the dimension $j_{0}$ and thus,
$h(\zeta_{\mathbf{m}}(\hat{z},{z}^{1}))\geq
h(\zeta_{\mathbf{m}}(\hat{z},{z}))$ and
$g_{i}(\zeta_{\mathbf{m}}({z}^{1},\hat{z}))=g_{i}(\zeta_{\mathbf{m}}({z},\hat{z}))$.
Then, it follows from (32) that
$f_{d,i}(z^{1},\hat{z};\mathbf{m},h(\cdot))\geq
f_{d,i}(z,\hat{z};\mathbf{m},h(\cdot))$. Repeating this procedure sequentially
for all dimensions $j$ for which case (i) in Lemma 1 holds (where $\tau$ is
the size of this set), we obtain:

$\displaystyle\begin{array}[]{c}f_{d,i}(z^{\tau},\hat{z};\mathbf{m},h(\cdot))\geq
f_{d,i}(z^{\tau-1},\hat{z};\mathbf{m},h(\cdot))\geq\dots\\\ \hskip
35.56593pt\geq f_{d,i}(z^{1},\hat{z};\mathbf{m},h(\cdot))\geq
f_{d,i}(z,\hat{z};\mathbf{m},h(\cdot)).\end{array}$ (38)

Next, we consider the rest of the dimensions $j^{\prime}$ that satisfy case
(ii) in Lemma 1. It follows from the second case in (19) that for such a
dimension $j_{0}^{\prime}\in\mathbb{N}_{n_{z}}$,

$\displaystyle\begin{array}[]{c}\zeta_{\mathbf{m},{j^{\prime}_{0}}}(\hat{z},\tilde{z})=\hat{z}_{j^{\prime}_{0}}=\zeta_{\mathbf{m},{j^{\prime}_{0}}}(\hat{z},z^{\tau}),\\\
\zeta_{\mathbf{m},{j^{\prime}_{0}}}(\tilde{z},\hat{z})=\tilde{z}_{j^{\prime}_{0}}\geq{z}^{\tau}_{j^{\prime}_{0}}=\zeta_{\mathbf{m},{j^{\prime}_{0}}}(z^{\tau},\hat{z}).\end{array}$
(41)

Repeating a similar procedure as for case (i), we define
${z}^{\tau+1}\in\mathcal{Z}$ as follows:
${z}^{\tau+1}_{j^{\prime}_{0}}=\tilde{z}_{j^{\prime}_{0}}$ and
${z}^{\tau+1}_{j}=z^{\tau}_{j},\forall j\neq j^{\prime}_{0}$. Thus,
$z^{\tau+1}\geq z^{\tau}$, and by (35),
$\zeta_{\mathbf{m}}(\hat{z},z^{\tau+1})=\zeta_{\mathbf{m}}(\hat{z},z^{\tau})$,
$\zeta_{\mathbf{m},j}(z^{\tau+1},\hat{z})=\zeta_{\mathbf{m},j}(z^{\tau},\hat{z}),\forall
j\neq j_{0}^{\prime}$ and
$\zeta_{\mathbf{m},j_{0}^{\prime}}(z^{\tau+1},\hat{z})\geq\zeta_{\mathbf{m},j_{0}^{\prime}}(z^{\tau},\hat{z})$.
Moreover, by case (ii) in the proof of Lemma 1 and Proposition 10,
$g_{i}(\cdot)$ is non-decreasing in the dimension $j_{0}^{\prime}$ and thus,
$g_{i}(\zeta_{\mathbf{m}}({z}^{\tau+1},\hat{z}))\geq
g_{i}(\zeta_{\mathbf{m}}({z}^{\tau},\hat{z}))$ and
$h(\zeta_{\mathbf{m}}(\hat{z},{z}^{\tau+1}))=h(\zeta_{\mathbf{m}}(\hat{z},{z}^{\tau}))$.
Then, it follows from (32) that
$f_{d,i}(z^{\tau+1},\hat{z};\mathbf{m},h(\cdot))\geq
f_{d,i}(z^{\tau},\hat{z};\mathbf{m},h(\cdot))$.

Repeating this procedure sequentially for all dimensions $j$ for which case
(ii) in Lemma 1 holds, we obtain:

$\displaystyle\begin{array}[]{c}f_{d,i}(z^{\tau},\hat{z};\mathbf{m},h(\cdot))\leq
f_{d,i}(z^{\tau+1},\hat{z};\mathbf{m},h(\cdot))\leq\dots\\\ \hskip
21.33955pt\leq f_{d,i}(\tilde{z},\hat{z};\mathbf{m},h(\cdot)),\end{array}$
(44)

where the last term is $f_{d,i}(\tilde{z},\hat{z};\mathbf{m},h(\cdot))$ since
there exist only two possible cases (i) or (ii) for each dimension. Combining
(38) and (44) yields $f_{d,i}(\tilde{z},\hat{z};\mathbf{m},h(\cdot))\geq
f_{d,i}({z},\hat{z};\mathbf{m},h(\cdot))$, which means that $f_{d,i}$ is non-
decreasing in its first argument. A n almost identical argument shows that
$f_{d,i}$ is non-increasing in its second argument. Thus, $f_{d,i}$ is mixed-
monotone. ∎

Theorem 1 mathematically introduces a family of decomposition functions (cf.
Definition 3), but the results are not yet tractable (cf. Definition 7), since
to build such a family, we have to search over $\mathbf{M}_{i}$ (cf. (22) and
(27)), which is an unbounded and infinite set, as well as over
$\mathcal{H}_{\mathbf{M}_{i}}$ (cf. (28)), which is an infinite-dimensional
space of functions. To overcome this problem, we propose tractable upper and
lower decomposition functions that only need to search over a _finite_ set of
supporting vector $\mathbf{M}_{i}^{c}\subset\mathbf{M}_{i}$ with the choice of
linear remainder functions
$h(\zeta)=\langle\mathbf{m},\zeta\rangle=\mathbf{m}^{\top}\zeta$, and prove
that these tractable decomposition functions are the tightest among the family
of decomposition functions in Theorem 1.

###### Theorem 2 (Tight and Tractable Remainder-Form Upper and Lower
Decomposition Functions).

Consider an ELLC vector field
$f_{i}:\mathcal{Z}\subset\mathbb{R}^{n_{z}}\to\mathbb{R}$ and let Assumptions
1 and 2 hold. Then, the tightest tractable (mixed-monotone) remainder-form
upper and lower decomposition functions with $z\geq\hat{z}$ are

$\displaystyle\begin{array}[]{l}\overline{f}_{d,i}(z,\hat{z})=\min\limits_{\mathbf{m}\in\mathbf{M}^{c}}f_{i}(\zeta_{\mathbf{m}}^{+})+\mathbf{m}^{\top}(\zeta_{\mathbf{m}}^{-}-\zeta_{\mathbf{m}}^{+}),\\\
\underline{f}_{d,i}(\hat{z},z)=\max\limits_{\mathbf{m}\in\mathbf{M}^{c}}f_{i}(\zeta_{\mathbf{m}}^{-})+\mathbf{m}^{\top}(\zeta_{\mathbf{m}}^{+}-\zeta_{\mathbf{m}}^{-}),\end{array}$
(47)

where $\zeta_{\mathbf{m}}^{+}\triangleq\zeta_{\mathbf{m}}(z,\hat{z})$ and
$\zeta_{\mathbf{m}}^{-}\triangleq\zeta_{\mathbf{m}}(\hat{z},z)$ with
$\zeta_{\mathbf{m}}(\cdot,\cdot)$ defined in (19) and the _finite_ set of
supporting vectors $\mathbf{M}_{i}^{c}$ defined as:

$\displaystyle\begin{array}[]{rl}\mathbf{M}_{i}^{c}\triangleq&\\{\mathbf{m}\in\mathbb{R}^{n_{z}}\
|\ \mathbf{m}_{j}=\max(({\overline{J}^{f}_{C}})_{ij},{0})\ \lor\\\ &\
\mathbf{m}_{j}=\min(({\underline{J}^{f}_{C}})_{ij},{0}),\forall
j\in\mathbb{N}_{n_{z}}\\}\end{array}$ (50)

if (15) is a discrete-time system and

$\displaystyle\begin{array}[]{l}\mathbf{M}_{i}^{c}\triangleq\\{\mathbf{m}\in\mathbb{R}^{n_{z}}\
|\ \mathbf{m}_{j}=\max(({\overline{J}^{f}_{C}})_{ij},{0})\ \lor\\\ \quad\
\mathbf{m}_{j}=\min(({\underline{J}^{f}_{C}})_{ij},{0}),\forall
j\in\mathbb{N}_{n_{z}},j\neq i,\mathbf{m}_{i}=0\\}\end{array}$ (53)

if (15) is a continuous-time system.

Moreover, we call the resulting inclusion function for an interval domain
$\mathbb{I}\mathcal{Z}=[\underline{z},\overline{z}]$ the remainder-form
inclusion function
$T^{f_{d,i}}_{R}\triangleq[\underline{f}_{d,i}(\underline{z},\overline{z}),\overline{f}_{d,i}(\overline{z},\underline{z})]$
(cf. Definition 2).

We will prove the above theorem in two steps with the help of the following
lemmas, where the two steps show that we can restrict our search for the
tightest upper and lower decomposition functions to a _finite_ set of
supporting vectors and the set of _linear_ remainder functions, respectively,
without introducing any conservatism.

###### Lemma 2 (Finite Set of Supporting Vectors).

Suppose the assumptions in Theorem 2 hold. Then, $\forall
z,\hat{z}\in\mathcal{Z}$, $\forall h(\cdot)\in\mathcal{H}_{\mathbf{M}_{i}}$
and for both $\operatorname*{opt}\in\\{\min,\max\\}$,

$\displaystyle\operatorname*{opt}\limits_{\begin{subarray}{c}\mathbf{m}\in\mathbf{M}_{i},\\\
h(\cdot)\in\mathcal{H}_{\mathbf{M}_{i}}\end{subarray}}f_{d,i}(z,\hat{z};\mathbf{m},h(\cdot))=\operatorname*{opt}\limits_{\begin{subarray}{c}\mathbf{m}\in\mathbf{M}_{i}^{c},\\\
h(\cdot)\in\mathcal{H}_{\mathbf{M}_{i}}\end{subarray}}f_{d,i}(z,\hat{z};\mathbf{m},h(\cdot)),$

where $f_{d,i}(z,\hat{z};\mathbf{m},h(\cdot))$ is defined in (17),
$\mathbf{M}_{i}$ in (22) or (27), $\mathbf{M}_{i}^{c}$ in (50) or (53), and
$\mathcal{H}_{\mathbf{M}_{i}}$ in (28).

###### Proof.

We consider $\mathbf{m}\in\mathbf{M}_{i}$ and construct
$\tilde{\mathbf{m}}\in\mathbf{M}_{i}^{c}$ as follows for all
$j\in\mathbb{N}_{n_{z}}$:
$\tilde{\mathbf{m}}_{j}=\max(({\overline{J}^{f}_{C}})_{ij},0)$ if
$\mathbf{m}_{j}\geq\max(({\overline{J}^{f}_{C}})_{ij},0)$ and
$\tilde{\mathbf{m}}_{j}=\min(({\underline{J}^{f}_{C}})_{ij},0)$ if
$\mathbf{m}_{j}\leq\min(({\underline{J}^{f}_{C}})_{ij},0)$. Then, it can be
easily verified from (19) that $\forall
z^{1},z^{2}\in\mathcal{Z},{\zeta}_{\mathbf{m}}(z^{1},z^{2})={\zeta}_{\tilde{\mathbf{m}}}(z^{1},z^{2})$,
and hence by (17),
$f_{d,i}(z^{1},z^{2};\mathbf{m},h(\cdot))=f_{d,i}(z^{1},z^{2};\tilde{\mathbf{m}},h(\cdot)){\color[rgb]{0,0,0},\forall
h(\cdot)\in\mathcal{H}_{\mathbf{M}_{i}}}$. Hence, for any
$\mathbf{m}\in\mathbf{M}_{i}$, there exists
$\tilde{\mathbf{m}}\in\mathbf{M}_{i}^{c}$ that admits an equivalent
decomposition function and correspondingly, the optimization over
$\mathbf{M}_{i}$ and $\mathbf{M}_{i}^{c}$ are equivalent. ∎

###### Lemma 3 (Linear Remainder Functions).

Suppose the assumptions in Theorem 2 hold. Then, $\forall
z,\hat{z}\in\mathcal{Z},z\geq\hat{z}$,

$\displaystyle\begin{array}[]{l}\min\limits_{\begin{subarray}{c}\mathbf{m}\in\mathbf{M}_{i}^{c},\\\
h(\cdot)\in\mathcal{H}_{\mathbf{M}_{i}}\end{subarray}}f_{d,i}(z,\hat{z};\mathbf{m},h(\cdot))=\min\limits_{\mathbf{m}\in\mathbf{M}^{c}}f_{i}(\zeta_{\mathbf{m}}^{+})+\mathbf{m}^{\top}(\zeta_{\mathbf{m}}^{-}-\zeta_{\mathbf{m}}^{+}),\\\
\max\limits_{\begin{subarray}{c}\mathbf{m}\in\mathbf{M}_{i}^{c},\\\
h(\cdot)\in\mathcal{H}_{\mathbf{M}_{i}}\end{subarray}}f_{d,i}(\hat{z},z;\mathbf{m},h(\cdot))=\max\limits_{\mathbf{m}\in\mathbf{M}^{c}}f_{i}(\zeta_{\mathbf{m}}^{-})+\mathbf{m}^{\top}(\zeta_{\mathbf{m}}^{+}-\zeta_{\mathbf{m}}^{-}),\vspace{-0.2cm}\end{array}$
(56)

where $\zeta_{\mathbf{m}}^{+}\triangleq\zeta_{\mathbf{m}}(z,\hat{z})$ and
$\zeta_{\mathbf{m}}^{-}\triangleq\zeta_{\mathbf{m}}(\hat{z},z)$ with
$f_{d,i}(z,\hat{z};\mathbf{m},h(\cdot))$ is defined in (17),
$\zeta_{\mathbf{m}}(\cdot,\cdot)$ defined in (19), $\mathbf{M}_{i}^{c}$ in
(50) or (53), and $\mathcal{H}_{\mathbf{M}_{i}}$ in (28).

###### Proof.

Consider any $h(\cdot)\in\mathcal{H}_{\mathbf{M}_{i}}$,
$\mathbf{m}\in\mathbf{M}_{i}$ and from (17),

$\displaystyle f^{+}_{d,i}$ $\displaystyle\triangleq
f_{d,i}(z,\hat{z};\mathbf{m},h(\cdot))=f_{i}(\zeta_{\mathbf{m}}^{+})+\Delta
h_{\mathbf{m}},$ (57) $\displaystyle f^{-}_{d,i}$ $\displaystyle\triangleq
f_{d,i}(\hat{z},z;\mathbf{m},h(\cdot))=f_{i}(\zeta_{\mathbf{m}}^{-})-\Delta
h_{\mathbf{m}},$ (58)

where $\Delta h_{\mathbf{m}}\triangleq
h(\zeta_{\mathbf{m}}^{-})-h(\zeta_{\mathbf{m}}^{+})$. Then, applying the
Clarke Mean Value Theorem [32, Chapter II, Theorem 1.3] to $\Delta
h_{\mathbf{m}}$, there exists
$\xi\in[{\underline{J}^{h}_{C}},{\overline{J}^{h}_{C}}]\subset\mathbf{M}_{i}$
such that

$\displaystyle\Delta
h_{\mathbf{m}}=\langle\xi,(\zeta_{\mathbf{m}}^{+}-\zeta_{\mathbf{m}}^{-})\rangle.$
(59)

Since $\xi\in\mathbf{M}_{i}$, by (22) and (27), we know that
$\xi_{j}\leq\min(({\underline{J}^{f}_{C}})_{ij},{0})$ or
$\xi_{j}\geq\max(({\overline{J}^{f}_{C}})_{ij},{0})$, $\forall
j\in\mathbb{N}_{n_{z}}$ (excluding $j=i$ for continuous-time systems where
$\xi_{j}=0$).

Then, similar to the proof of Lemma 1, $\forall j\in\mathbb{N}_{n_{x}}$, we
can consider two cases corresponding to the two cases in (19):

  1. (i)

$\xi_{j}\geq\max(({\overline{J}^{f}_{C}})_{ij},{0})\geq 0$: From (19),
$\zeta_{\mathbf{m},j}^{-}=\zeta_{\mathbf{m},j}(\hat{z},{z})={z}_{j}$,
$\zeta_{\mathbf{m},j}^{+}=\zeta_{\mathbf{m},j}(z,\hat{z})=\hat{z}_{j}$ and
${z}_{j}\geq\hat{z}_{j}$, thus, we have
$\zeta_{\mathbf{m},j}^{-}-\zeta_{\mathbf{m},j}^{+}\geq 0$ and
$\xi_{j}(\zeta_{\mathbf{m},j}^{-}-\zeta_{\mathbf{m},j}^{+})\geq\max(({\overline{J}^{f}_{C}})_{ij},{0})(\zeta_{\mathbf{m},j}^{-}-\zeta_{\mathbf{m},j}^{+})$.
Then, the minimum of $f^{+}_{d,i}$ in (57) and the maximum of $f^{-}_{d,i}$ in
(58) are attained in (56) when
$\xi_{j}=\max(({\overline{J}^{f}_{C}})_{ij},{0})\in(\mathbf{M}^{c}_{i})_{j}$.

  2. (ii)

$\xi_{j}\leq\min(({\underline{J}^{f}_{C}})_{ij},{0})\leq 0$: From (19),
$\zeta_{\mathbf{m},j}^{-}=\zeta_{\mathbf{m},j}(\hat{z},{z})=\hat{z}_{j}$,
$\zeta_{\mathbf{m},j}^{+}=\zeta_{\mathbf{m},j}(z,\hat{z})={z}_{j}$ and
${z}_{j}\geq\hat{z}_{j}$, thus, we have
$\zeta_{\mathbf{m},j}^{-}-\zeta_{\mathbf{m},j}^{+}\leq 0$ and
$\xi_{j}(\zeta_{\mathbf{m},j}^{-}-\zeta_{\mathbf{m},j}^{+})\geq\min(({\underline{J}^{f}_{C}})_{ij},{0})(\zeta_{\mathbf{m},j}^{-}-\zeta_{\mathbf{m},j}^{+})$.
Then, the minimum of $f^{+}_{d,i}$ in (57) and the maximum of $f^{-}_{d,i}$ in
(58) are attained in (56) when
$\xi_{j}=\min(({\underline{J}^{f}_{C}})_{ij},{0})\in(\mathbf{M}^{c}_{i})_{j}$.

Finally, we can restrict our search to the class of linear remainder functions
$h(\zeta)=\langle\mathbf{m},\zeta\rangle=\mathbf{m}^{\top}\zeta$ with
$\mathbf{m}\in\mathbf{M}^{c}_{i}$, since it can achieve the optima in (56). ∎

###### Proof of Theorem 2.

First, by repeatedly applying Corollary 1 on all the decomposition functions
in the family (17) and the fact that the upper and lower decomposition
functions can be optimized independently, it can be seen that the tightest
upper and lower decomposition functions with
$z,\hat{z}\in\mathcal{Z},z\geq\hat{z}$ are

$\displaystyle\begin{array}[]{l}\overline{f}_{d,i}(z,\hat{z})=\min\limits_{\begin{subarray}{c}\mathbf{m}\in\mathbf{M}_{i},h(\cdot)\in\mathcal{H}_{\mathbf{M}_{i}}\end{subarray}}f_{d,i}(z,\hat{z};\mathbf{m},h(\cdot)),\\\
\underline{f}_{d,i}(\hat{z},z)=\max\limits_{\begin{subarray}{c}\mathbf{m}\in\mathbf{M}_{i},h(\cdot)\in\mathcal{H}_{\mathbf{M}_{i}}\end{subarray}}f_{d,i}(\hat{z},z;\mathbf{m},h(\cdot)),\end{array}$

Then, by Lemmas 2 and 3, we obtain the tractable and tight upper and lower
decomposition functions in (47). ∎

Theorem 2 guarantees that in order to obtain the tightest possible
decomposition function in the form of (17), it is sufficient to only search
over a _finite_ set of supporting vectors
$\mathbf{M}_{i}^{\color[rgb]{0,0,0}c}$ and the class of _linear_ remainder
functions with gradients from this same set
$\mathbf{M}_{i}^{\color[rgb]{0,0,0}c}$, i.e.,
$h(\zeta)={\color[rgb]{0,0,0}\langle}\mathbf{m},\zeta{\color[rgb]{0,0,0}\rangle}=\mathbf{m}^{\top}\zeta,\forall\mathbf{m}\in\mathbf{M}_{i}^{c}$,
where the search space is the finite and countable set $\mathbf{M}^{c}_{i}$.
Hence, the optimal search for the tightest decomposition functions is
computable/tractable according to Definition 7.

Further, the result in Theorem 2 can be applied to each $f_{i}$,
$i\in\mathbb{N}_{n_{x}}$ of the function $f$ to obtain the tightest remainder-
form decomposition functions from the family of _remainder-form CJSS
decomposition functions_ in (17). This is summarized in Algorithm 1, which
takes an interval domain $\mathbb{I}\mathcal{Z}=[\underline{z},\overline{z}]$,
the function $f$ and its Clarke  Jacobians, $\overline{J}_{C}^{f}$,
$\underline{J}_{C}^{f}$, as inputs, and  outputs the remainder-form inclusion
function
$T^{f_{d}}_{R}\triangleq[\underline{f}_{d}(\underline{z},\overline{z}),\overline{f}_{d}(\overline{z},\underline{z})]$
(cf. Definition 2).

Algorithm 1 Remainder-Form Decomposition Functions

1:function
$T^{f_{d}}_{R}$($f(\cdot),\overline{J}_{C}^{f},\underline{J}_{C}^{f},\overline{z},\underline{z}$)

2: Initialize:
$\overline{f}_{d}\leftarrow\infty,\underline{f}_{d}\leftarrow-\infty$;

3: for $i=1$ to $n_{x}$ do

4: if (15) is a discrete-time system then

5: $\mathbf{M}_{i}^{c}\triangleq\\{\mathbf{m}\in\mathbf{M}_{i}|\
\mathbf{m}_{j}=\min(({\underline{J}_{C}^{f}})_{ij},0)\ \vee$

6: $\mathbf{m}_{j}=\max(({\overline{J}_{C}^{f}})_{ij},{0}),\forall
j\in\mathbb{N}_{n_{z}}\\}$.

7: end if

8: if (15) is a continuous-time system then

9: $\mathbf{M}_{i}^{c}\triangleq\\{\mathbf{m}\in\mathbf{M}_{i}|\
\mathbf{m}_{j}=\min(({\underline{J}_{C}^{f}})_{ij},0)\ \vee\mathbf{m}_{j}=$

10: $\max(({\overline{J}_{C}^{f}})_{ij},{0}),\forall j\in\mathbb{N}_{n_{z}}$,
$j\neq i$ and $\mathbf{m}_{i}=0\\}$.

11: end if

12: for $\mathbf{m}\in\mathbf{M}^{c}_{i}$ do

13: for $j=1$ to $n_{z}$ do

14: if (15) is a continuous-time system $\land(j=i)$ then

15: $\zeta^{+}_{\mathbf{m},j}\leftarrow x_{i}$;
$\zeta^{-}_{\mathbf{m},j}\leftarrow x_{i}$;

16: else if $\mathbf{m}_{j}=\min(({\underline{J}_{C}^{f}})_{ij},{0})$ then

17: $\zeta^{+}_{\mathbf{m},j}\leftarrow\overline{z}_{j}$;
$\zeta^{-}_{\mathbf{m},j}\leftarrow\underline{z}_{j}$;

18: else

19: $\zeta^{+}_{\mathbf{m},j}\leftarrow\underline{z}_{j}$;
$\zeta^{-}_{\mathbf{m},j}\leftarrow\overline{z}_{j}$;

20: end if

21:
$\overline{f}_{d,i}\leftarrow\min(\overline{f}_{d,i},f_{i}(\zeta^{+}_{\mathbf{m}})+\mathbf{m}^{\top}(\zeta^{-}_{\mathbf{m}}-\zeta^{+}_{\mathbf{m}}))$;

22:
$\underline{f}_{d,i}\leftarrow\min(\underline{f}_{d,i},f_{i}(\zeta^{-}_{\mathbf{m}})+\mathbf{m}^{\top}(\zeta^{+}_{\mathbf{m}}-\zeta^{-}_{\mathbf{m}}))$;

23: end for

24: end for

25: end for

26: return $\overline{f}_{d},\underline{f}_{d}$;

27:end function

###  III-B Error Bounds

Next, we more formally characterize the tightness of our proposed family of
remainder-form decomposition functions in (17), where we use the
metric/measure of tightness in (5) that is based on the Hausdorff distance
function. In particular, we derive lower and upper bounds on the _over-
approximation error_ of the image set/range of function $f(\cdot)$, where the
lower bound is achievable by our tight and tractable decomposition functions
in Theorem 2.

###### Theorem 3 (Error Bounds).

Suppose that all the assumptions in Theorem 1 are satisfied for each
$f_{i}:\mathcal{Z}\subset\mathbb{R}^{n_{z}}\to\mathbb{R},i\in\mathbb{N}_{n_{x}}$.
Let
$T^{f}_{O}(\mathbb{I}\mathcal{Z})\triangleq[\underline{f}^{\text{true}},\overline{f}^{\text{true}}]{\color[rgb]{0,0,0}\triangleq[\min\limits_{z\in\mathbb{I}\mathcal{Z}}f(z),\max\limits_{z\in\mathbb{I}\mathcal{Z}}f(z)]}=T^{f_{d}}_{O}(\mathbb{I}\mathcal{Z})$
be the tightest enclosing interval of the true image set/range of $f(\cdot)$
over
$\mathbb{I}\mathcal{Z}=[\underline{z},\overline{z}]\in\mathbb{IR}^{n_{z}}$
(cf. Definition 1),
$T_{\underline{\mathbf{m}},\overline{\mathbf{m}}}^{f_{d{\color[rgb]{0,0,0},i}}}(\mathbb{I}\mathcal{Z})\supseteq
T^{f_{\color[rgb]{0,0,0}i}}_{O}(\mathbb{I}\mathcal{Z})$,
${\color[rgb]{0,0,0}\forall i\in\mathbb{N}_{n_{x}},}$ given in (29), be
inclusion functions using the family of decomposition functions in (17),  and
$T^{f_{d}}_{\mathbb{M}}(\mathbb{I}\mathcal{Z})\supseteq
T^{f}_{O}(\mathbb{I}\mathcal{Z})$ be any inclusion function such that
$T^{f_{d,i}}_{\mathbb{M}}(\mathbb{I}\mathcal{Z})=T_{\underline{\mathbf{m}},\overline{\mathbf{m}}}^{f_{d,i}}(\mathbb{I}\mathcal{Z})$,
$\forall i\in\mathbb{N}_{n_{x}}$, for some
$\underline{\mathbf{m}},\overline{\mathbf{m}}\in\mathbf{M}_{i}$. Then, the
following inequalities hold:

  1. (i)

$\underline{q}_{f_{d}}(\mathbb{I}\mathcal{Z})\leq
q(T_{{\color[rgb]{0,0,0}\mathbb{M}}}^{f_{d}}(\mathbb{I}\mathcal{Z}),T^{f}_{O}(\mathbb{I}\mathcal{Z}))$,
and

  2. (ii)

$\underline{q}_{f_{d}}(\mathbb{I}\mathcal{Z})\leq\overline{q}_{f_{d}}(\mathbb{I}\mathcal{Z})\leq\hat{\overline{q}}_{f_{d}}(\mathbb{I}\mathcal{Z})$,

with the tightness metric
$q([\underline{v},\overline{v}],[\underline{w},\overline{w}])=\max\limits_{i\in\mathbb{N}_{n_{x}}}\max\\{|\overline{v}_{i}-\overline{w}_{i}|,|\underline{v}_{i}-\underline{w}_{i}|\\}$
defined in (5), and with the bounds
$\underline{q}_{f_{d}}(\mathbb{I}\mathcal{Z})\triangleq\max_{i\in\mathbb{N}_{n_{x}}}\underline{q}_{f_{d,i}}(\mathbb{I}\mathcal{Z})$,
$\overline{q}_{f_{d}}(\mathbb{I}\mathcal{Z})\triangleq\max_{i\in\mathbb{N}_{n_{x}}}\overline{q}_{f_{d,i}}(\mathbb{I}\mathcal{Z})$
and
$\hat{\overline{q}}_{f_{d}}(\mathbb{I}\mathcal{Z})\triangleq\max_{i\in\mathbb{N}_{n_{x}}}\hat{\overline{q}}_{f_{d,i}}(\mathbb{I}\mathcal{Z})$,
where for each $i\in\mathbb{N}_{n_{x}}$,

$\displaystyle\begin{array}[]{rl}\underline{q}_{f_{d,i}}(\mathbb{I}\mathcal{Z})&\triangleq\max\\{\min\limits_{\mathbf{m}\in\mathbf{M}^{c}_{i}}\Delta^{1}_{i,\mathbf{m}}-\overline{f}_{i}^{\text{true}},\underline{f}_{i}^{\text{true}}-\max\limits_{\mathbf{m}\in\mathbf{M}^{c}_{i}}\Delta^{2}_{i,\mathbf{m}}\\},\\\
\hat{\overline{q}}_{f_{d,i}}(\mathbb{I}\mathcal{Z})&\triangleq\min\limits_{\mathbf{m}\in\mathbf{M}_{i}^{c}}\Delta^{3}_{i,\mathbf{m}},\\\
{\overline{q}}_{f_{d,i}}(\mathbb{I}\mathcal{Z})&\triangleq\min\\{\min\limits_{\mathbf{m}\in\mathbf{M}_{i}^{c}}\Delta^{3}_{i,\mathbf{m}},\min\limits_{\mathbf{m}\in\mathbf{M}_{i}^{c}}\Delta^{3}_{i,\mathbf{m}}+\Delta^{4}_{i,\mathbf{m}}\\},\end{array}$

with $\mathbf{M}_{i}$ and $\mathbf{M}_{i}^{c}$ in (22), (50) or  (27), (53),
respectively,

$\displaystyle\begin{array}[]{ll}\Delta^{1}_{i,\mathbf{m}}\triangleq
f_{i}(\zeta_{i,\mathbf{m}}^{+})+\Delta^{3}_{i,\mathbf{m}},&\Delta^{2}_{i,\mathbf{m}}\triangleq
f_{i}(\zeta_{i,\mathbf{m}}^{-})-\Delta^{3}_{i,\mathbf{m}},\\\
\Delta^{3}_{i,\mathbf{m}}\triangleq
h_{i}(\zeta_{i,\mathbf{m}}^{-})-h_{i}(\zeta_{i,\mathbf{m}}^{+}),&\Delta^{4}_{i,\mathbf{m}}\triangleq
f_{i}(\zeta_{i,\mathbf{m}}^{+})-f_{i}(\zeta_{i,\mathbf{m}}^{-}),\end{array}$
(62)

as well as
$\zeta_{i,\mathbf{m}}^{+}\triangleq\zeta_{i,\mathbf{m}}(\overline{z},\underline{z})$,
$\zeta_{i,\mathbf{m}}^{-}\triangleq\zeta_{i,\mathbf{m}}(\underline{z},\overline{z})$
and $\zeta_{i,\mathbf{m}}(.,.)$ defined in Theorem 1. Further, without loss of
tightness (cf. Lemma 3), we can replace $\Delta^{3}_{i,\mathbf{m}}$ in (62)
with

$\displaystyle\Delta^{3}_{i,\mathbf{m}}\triangleq\mathbf{m}^{\top}(\zeta_{i,\mathbf{m}}^{-}-\zeta_{i,\mathbf{m}}^{+}),$
(63)

and the lower bound $\underline{q}_{f_{d}}(\mathbb{I}\mathcal{Z})$ is attained
by the upper and lower decomposition functions (i.e., $T_{R}^{f_{d}}$) in
Theorem 2.

###### Proof.

First, from (17), (29) and (62), we find

$\displaystyle\tilde{q}(T_{\underline{\mathbf{m}},\overline{\mathbf{m}}}^{f_{d,i}}(\mathbb{I}\mathcal{Z}),T^{f_{i}}_{O}(\mathbb{I}\mathcal{Z}))=\max\\{\Delta^{1}_{i,\mathbf{m}}-\overline{f}_{i}^{\text{true}},\underline{f}_{i}^{\text{true}}-\Delta^{2}_{i,\mathbf{m}}\\}$
$\displaystyle\geq\max\\{\min_{\mathbf{m}\in\mathbf{M}^{c}_{i}}\Delta^{1}_{i,\mathbf{m}}-\overline{f}_{i}^{\text{true}},\underline{f}_{i}^{\text{true}}-\max_{\mathbf{m}\in\mathbf{M}^{c}_{i}}\Delta^{2}_{i,\mathbf{m}}\\},$
(64)

where the inequality in (i) follows from independently searching over
${\color[rgb]{0,0,0}\mathbf{M}^{c}_{i}\subset\mathbf{M}_{i}}$ to minimize each
argument of the maximization,  as well as the fact that by Lemmas 2 and 3, we
can apply (63) and search only over
$\mathbf{M}^{c}_{i}\subset{\color[rgb]{0,0,0}\mathbf{M}_{i}}$ without any
conservatism, and it can be verified that by construction, the lower bound
$\underline{q}_{f}(\mathbb{I}\mathcal{Z})$ is attained by the decomposition
functions in Theorem 2.

To obtain (ii), we apply [19, Theorem 4-(b)], which proved that for any
remainder-form inclusion functions with remainder function $r_{i}(\cdot)$
satisfying
$r_{i}(\mathbb{I}\mathcal{Z})\subset[\underline{r}_{i},\overline{r}_{i}]$,
$\tilde{q}(W^{R}_{f_{d,i}}(\mathbb{I}\mathcal{Z}),V_{f_{i}}(\mathbb{I}\mathcal{Z}))\leq\overline{r}_{i}-\underline{r}_{i}$
holds. For the first inequality  in (ii), only $h_{i}(\cdot)$ is considered as
the remainder function, while in the second inequality, both $h_{i}(\cdot)$
and $g_{i}(\cdot)\triangleq f_{i}(\cdot)-h_{i}(\cdot)$ are considered as
remainder functions, separately, with the minimum chosen as the bound.
Moreover, since $h_{i}(\cdot)$ is CJSS, by (19), $\overline{h}_{i}$ and
$\underline{h}_{i}$ are attained at the corner points given by
$\zeta^{-}_{i,\mathbf{m}}$ and $\zeta^{+}_{i,\mathbf{m}}$, respectively, with
$\Delta^{3}_{i,\mathbf{m}}\triangleq\overline{h}_{i}-\underline{h}_{i}$.
Further, since $g_{i}(\cdot)$ is aligned with $-h_{i}(\cdot)$ by Lemma 1,
$\overline{g}_{i}$ and $\underline{g}_{i}$ are attained at the corner points
given by $\zeta^{+}_{i,\mathbf{m}}$ and $\zeta^{-}_{i,\mathbf{m}}$,
respectively, with
$\Delta^{3}_{i,\mathbf{m}}+\Delta^{4}_{i,\mathbf{m}}\triangleq\overline{g}_{i}-\underline{g}_{i}$.
∎

The above result holds for both discrete-time and continuous-time systems
(with overloading described in Definition 1). Further, note that lower bound
$\underline{q}_{f,d}(\mathbb{I}\mathcal{Z})$ is attainable by $T^{f_{d}}_{R}$
but since it is a function of the unknown $\overline{f}^{true}$ and
$\underline{f}^{true}$, it cannot be computed. Thus, its upper bounds
${\overline{q}}_{f_{d}}(\mathbb{I}\mathcal{Z})$ and
$\hat{\overline{q}}_{f_{d}}(\mathbb{I}\mathcal{Z})$ in (ii) that are
independent of $\overline{f}^{true}$ and $\underline{f}^{true}$ are more
useful, e.g., as worst case function over-approximation error bounds in
reachability and robust control problems.

###  III-C Convergence Rate and Subdivision Principle

In this subsection, we study the convergence rate of our proposed
$T^{f_{d}}_{R}$, i.e., the rate at which its approximation error goes to zero,
when the domain interval diameter $d(\mathbb{I}\mathcal{Z})$ shrinks. We show
that when using $T^{f_{d}}_{R}$, the error converges at least linearly, which
is also the convergence rate of natural inclusions $T^{f}_{N}$ [20, Chapter
6]. Further, we show that the _subdivision principle_ introduced in [21] can
be applied to improve the convergence. We first introduce the notion of
convergence rate, inspired by [21], and then, we present the convergence rate
and the subdivision principle for our proposed $T^{f_{d}}_{R}$.

###### Definition 11 (Convergence Rate).

An inclusion function $T^{f}:\mathbb{IR}^{n_{z}}\to\mathbb{IR}^{n_{x}}$ for an
ELLC vector field
$f:\mathcal{Z}\subset\mathbb{R}_{n_{z}}\to\mathbb{R}^{n_{x}}$ has a
convergence rate $\alpha>0$, if

$\displaystyle
q(T^{f}(\mathbb{I}\mathcal{Z}),T^{f}_{O}(\mathbb{I}\mathcal{Z}))\leq\beta\,d(\mathbb{I}\mathcal{Z})^{\alpha},$
(65)

for some $\beta>0$, where $T{\color[rgb]{0,0,0}{}^{f}}(\mathbb{I}\mathcal{Z})$
is the interval over-approximation of the range of $f(\cdot)$ over
$\mathbb{I}\mathcal{Z}$ (i.e., an inclusion function),
$T^{f}_{O}(\mathbb{I}\mathcal{Z})$ is the tightest inclusion function (cf.
Definition 1), $q(\cdot,\cdot)$ is defined in (5) and
$d(\mathbb{I}\mathcal{Z})\triangleq\|\overline{z}-\underline{z}\|_{\infty}$.

###### Theorem 4 (Convergence Rate and Subdivision Principle for
$T^{f_{d}}_{R}$).

The $T^{f_{d}}_{R}$ inclusion function for any ELLC vector field
$f:\mathcal{Z}\subset\mathbb{R}_{n_{z}}\to\mathbb{R}^{n_{x}}$ satisfies:

$\displaystyle{\color[rgb]{0,0,0}\underline{q}_{f_{d}}}(\mathbb{I}\mathcal{Z})=q(T_{R}^{f_{d}}(\mathbb{I}\mathcal{Z}),T^{f}_{O}(\mathbb{I}\mathcal{Z}))\leq\beta^{f}_{R}d(\mathbb{I}\mathcal{Z}),$
(66)

with convergence rate $\alpha=1$ for some $\beta^{f}_{R}>0$. Moreover,
applying the subdivision principle, we have

$\displaystyle
q(T_{R}^{f_{d}}(\mathbb{I}\mathcal{Z};k),T^{f}_{O}(\mathbb{I}\mathcal{Z};k))\leq\frac{\gamma^{f}_{R}}{k},$
(67)

where $\mathbb{I}\mathcal{Z}$ is subdivided into $k^{n_{z}}$ interval vectors
$\mathbb{I}\mathcal{Z}^{l},l\in\mathbb{N}_{k^{n_{z}}}$ (i.e., with $k$
divisions in each dimension such that
$d(\mathbb{I}\mathcal{Z}^{l}_{j})=\frac{d(\mathbb{I}\mathcal{Z}_{j})}{k}$ for
$j\in\mathbb{N}_{n_{z}}$, $l\in\mathbb{N}_{k^{n_{z}}}$),
$T_{R}^{f_{d}}(\mathbb{I}\mathcal{Z};k)\triangleq\bigcup_{l=1}^{k^{n_{z}}}T_{R}^{f_{d}}(\mathbb{I}\mathcal{Z}^{l})$
and
$T^{f}_{O}(\mathbb{I}\mathcal{Z};k)\triangleq\bigcup_{l=1}^{k^{n_{z}}}T^{f}_{O}(\mathbb{I}\mathcal{Z}^{l})$.

###### Proof.

For a linear remainder function $h(z)=\mathbf{m}^{\top}z$ with
$\mathbf{m}\triangleq\begin{bmatrix}{\mathbf{m}}_{1}\ \ldots\
{\mathbf{m}}_{n_{x}}\end{bmatrix}$ (cf. Lemma 3),
$\hat{\overline{q}}_{f_{d}}(\mathbb{I}\mathcal{Z})$ in Theorem 3 can be upper
bounded by triangle inequality by (66) with
$\beta^{f}_{R}=\max_{i\in\mathbb{N}_{n_{x}}}\|\mathbf{m}_{i}\|_{\infty}$.
Further, the proof of (67) follows the same lines as the proof of [33, Theorem
4.1]. ∎

###  III-D Set Inversion Algorithm

The remainder-form decomposition functions returned by Algorithm 1 can be used
with the generalized embedding system in Definition 5 to over-approximate
(unconstrained) reachable sets of a dynamic system governed by the vector
field $f(\cdot)$, which corresponds to the propagated/predicted sets in state
observers/estimators. However, when additional state constraint information is
available (e.g., sensor observations/measurements in state estimation
problems, known safety constraints from system design and manufactured
constraints from modeling redundancy [23, 24, 25, 26, 27]), an additional set
inversion (also known as update or refinement) step will allow us to take the
advantage of the constraints to shrink the propagated sets, i.e., to obtain a
tighter subset of the propagated set that is compatible/consistent with the
given constraints. Further details about the application of the
reachable/propagated set and set inversion algorithms will be described in
Section V.

Formally, given the constraint/observation function
$\mu({\color[rgb]{0,0,0}x},u)$ in (15) with known $u$, a
constraint/observation set $\mathcal{Y}$ with maximal and minimal values
$\overline{y},\underline{y}$ (satisfying
$\mu({\color[rgb]{0,0,0}x},u)\subset\mathcal{Y}\subseteq[\underline{y},\overline{y}]$)
and a prior (or propagated/predicted) reachable interval
$\mathbb{I}\mathcal{{\color[rgb]{0,0,0}X}}^{p}=[\underline{{\color[rgb]{0,0,0}x}}^{p},\overline{{\color[rgb]{0,0,0}x}}^{p}]$,
we wish to find an updated/refined interval
$\mathbb{I}\mathcal{{\color[rgb]{0,0,0}X}}^{u}\subseteq\mathbb{I}\mathcal{{\color[rgb]{0,0,0}X}}^{p}$,
such that (16) holds, i.e., to solve Problem 2. Finding
$\mathbb{I}\mathcal{{\color[rgb]{0,0,0}X}}^{u}$ in (16) is called the _set
inversion_ problem [18]. To our best knowledge, existing set inversion
algorithms/operators either compute subpavings (i.e., unions of intervals)
instead of an interval using (conservative) natural inclusions (SIVIA [18,
Chapter 3]) or only applies if relatively restrictive monotonicity assumptions
hold ($\mathcal{I}_{G}$ [23, Algorithm 1]).

In this section, leveraging our proposed _discrete-time_ decomposition-based
inclusion functions for an ELLC function
${\nu}({\color[rgb]{0,0,0}x})\triangleq\mu({\color[rgb]{0,0,0}x},u)$ with
known $u$, i.e., $T^{{\nu}_{d}}_{R}$ , we develop a novel set inversion
algorithm that solves Problem 2, which is summarized in Algorithm 2. The main
idea behind Algorithm 2 is based on the observation that a candidate interval
$\mathbb{I}\Xi\triangleq[\underline{\xi},\overline{\xi}]\subseteq{\mathbb{I}\mathcal{{\color[rgb]{0,0,0}X}}}^{p}$
that satisfies
$T_{R}^{\nu_{d}}(\mathbb{I}\Xi)\cap[\underline{y},\overline{y}]=\emptyset$
(i.e., if $\underline{\nu}_{d}(\underline{\xi},\overline{\xi})>\overline{y}$
_or_ if $\overline{\nu}_{d}(\overline{\xi},\underline{\xi})<\underline{y}$) is
_incompatible/inconsistent_ with the set
$\\{{\color[rgb]{0,0,0}x}\in\mathbb{I}\mathcal{{\color[rgb]{0,0,0}X}}^{p}\,|\,\underline{y}\leq\mu({\color[rgb]{0,0,0}x},u)\leq\overline{y}\\}$
and can be eliminated/ruled out from
$\mathbb{I}\mathcal{{\color[rgb]{0,0,0}X}}^{p}$ and  thus, shrinking
$\mathbb{I}\mathcal{{\color[rgb]{0,0,0}X}}^{u}$.

Using this idea, starting from the prior/propagated interval and using
bisection for each dimension, Algorithm 2 shrinks the compatible interval from
below and/or above if $\underline{\nu}_{d}$ is strictly greater than
$\overline{y}$ or if $\overline{\nu}_{\color[rgb]{0,0,0}d}$ is strictly
smaller than $\underline{y}$ (cf. Lines 8 and 18). Repeating this procedure
along with  bisections with a threshold $\epsilon$, the candidate intervals
that are determined to be _inconsistent_ with the constraint/observation set
are ruled out. Note that the ordering of the dimensions in the ‘for’ loop on
Line 3 may have an impact on the tightness of the returned interval
$\mathbb{I}\mathcal{{\color[rgb]{0,0,0}X}}^{u}=[\underline{{\color[rgb]{0,0,0}x}}^{u},\overline{{\color[rgb]{0,0,0}x}}_{u}]$
and  so, it may be desirable to tailor the order to the problem at hand, to
randomize the order or to repeat the algorithm with the previous
$\mathbb{I}\mathcal{{\color[rgb]{0,0,0}X}}^{u}$ as
$\mathbb{I}\mathcal{{\color[rgb]{0,0,0}X}}^{p}$ multiple times.

The following result shows that Algorithm 2 returns
$\mathbb{I}\mathcal{{\color[rgb]{0,0,0}X}}^{u}=[\underline{{\color[rgb]{0,0,0}x}}^{u},\overline{{\color[rgb]{0,0,0}x}}_{u}]$
that satisfies (16), i.e., solves Problem 2.

Algorithm 2 Set Inversion based on Discrete-Time $T^{\nu_{d}}_{R}$

1:function Set-
Inv($\nu(\cdot),\overline{J}_{C}^{\nu},\underline{J}_{C}^{\nu},\overline{{\color[rgb]{0,0,0}x}}^{p},\underline{{\color[rgb]{0,0,0}x}}^{p},\overline{y},\underline{y},\epsilon$)

2: Initialize:
$\overline{{\color[rgb]{0,0,0}x}}^{u}\leftarrow\overline{{\color[rgb]{0,0,0}x}}^{p},\underline{{\color[rgb]{0,0,0}x}}^{u}\leftarrow\underline{{\color[rgb]{0,0,0}x}}^{p}$;

3: for ${\color[rgb]{0,0,0}i}=1$ to $n_{{\color[rgb]{0,0,0}x}}$ do

4:
$\underline{\zeta}\leftarrow\underline{{\color[rgb]{0,0,0}x}}^{u}_{{\color[rgb]{0,0,0}i}}$;
$\overline{\zeta}\leftarrow\overline{{\color[rgb]{0,0,0}x}}^{u}_{{\color[rgb]{0,0,0}i}}$;

5: while $\overline{\zeta}-\underline{\zeta}>\epsilon$ do

6: $\zeta_{m}\leftarrow\frac{1}{2}(\overline{\zeta}+\underline{\zeta})$;
$\overline{\xi}\leftarrow\overline{{\color[rgb]{0,0,0}x}}^{u}$;
$\underline{\xi}\leftarrow\underline{{\color[rgb]{0,0,0}x}}^{u}$;
$\underline{\xi}_{{\color[rgb]{0,0,0}i}}\leftarrow\zeta_{m}$;

7: $(\overline{\nu}_{d},\underline{\nu}_{d})\leftarrow
T^{\nu_{d}}_{R}(\nu(\cdot),\overline{J}_{C}^{\nu},\underline{J}_{C}^{\nu},\overline{\xi},\underline{\xi})$;
(Algorithm 1)

8: if $(\overline{\nu}_{d}<\underline{y})\ \vee\
(\underline{\nu}_{d}>\overline{y})$ then

9: $\overline{\zeta}\leftarrow\zeta_{m}$;
$\overline{{\color[rgb]{0,0,0}x}}^{u}_{{\color[rgb]{0,0,0}i}}\leftarrow\overline{\zeta}$;

10: else

11: $\underline{\zeta}\leftarrow\zeta_{m}$;

12: end if

13: end while

14:
$\underline{\zeta}\leftarrow\underline{{\color[rgb]{0,0,0}x}}^{u}_{{\color[rgb]{0,0,0}i}}$;
$\overline{\zeta}\leftarrow\overline{{\color[rgb]{0,0,0}x}}^{u}_{{\color[rgb]{0,0,0}i}}$;

15: while $\overline{\zeta}-\underline{\zeta}>\epsilon$ do

16: $\zeta_{m}\leftarrow\frac{1}{2}(\overline{\zeta}+\underline{\zeta})$;
$\overline{\xi}\leftarrow\overline{{\color[rgb]{0,0,0}x}}^{u}$;
$\underline{\xi}\leftarrow\underline{{\color[rgb]{0,0,0}x}}^{u}$;
$\overline{\xi}_{{\color[rgb]{0,0,0}i}}\leftarrow\zeta_{m}$;

17: $(\overline{\nu}_{d},\underline{\nu}_{d})\leftarrow
T^{\nu_{d}}_{R}(\nu(\cdot),\overline{J}_{C}^{\nu},\underline{J}_{C}^{\nu},\overline{\xi},\underline{\xi})$;
(Algorithm 1)

18: if $(\overline{\nu}_{d}<\underline{y})\ \vee\
(\underline{\nu}_{d}>\overline{y})$ then

19: $\underline{\zeta}\leftarrow\zeta_{m}$;
$\underline{{\color[rgb]{0,0,0}x}}^{u}_{{\color[rgb]{0,0,0}i}}\leftarrow\underline{\zeta}$;

20: else

21: $\overline{\zeta}\leftarrow\zeta_{m}$;

22: end if

23: end while

24: end for

25: return
$\overline{{\color[rgb]{0,0,0}x}}^{u},\underline{{\color[rgb]{0,0,0}x}}^{u}$;

26:end function

###### Proposition 11.

Suppose Assumptions 1 and 2 hold and consider an ELLC constraint/observation
function
$\nu:\mathcal{{\color[rgb]{0,0,0}X}}\subset\mathbb{R}^{n_{{\color[rgb]{0,0,0}x}}}\to\mathbb{R}^{n_{{\color[rgb]{0,0,0}\mu}}}$,
where $\nu({\color[rgb]{0,0,0}x})\triangleq\mu({\color[rgb]{0,0,0}x},u)$ with
known $u$, a constraint/observation set
$\mathcal{Y}_{t}\subseteq[\underline{y}_{t},\overline{y}_{t}]$ and a
prior/propagated interval
$\mathbb{I}\mathcal{{\color[rgb]{0,0,0}X}}^{p}\triangleq[\underline{{\color[rgb]{0,0,0}x}}_{p},\overline{{\color[rgb]{0,0,0}x}}_{p}]\in\mathbb{IR}^{n_{{\color[rgb]{0,0,0}x}}}$.
Then, the updated/refined interval
$\mathbb{I}\mathcal{{\color[rgb]{0,0,0}X}}^{u}\triangleq[\underline{{\color[rgb]{0,0,0}x}}^{u},\overline{{\color[rgb]{0,0,0}x}}^{u}]$
returned by Algorithm 2 satisfies (16).

###### Proof.

Obviously,
$\mathbb{I}\mathcal{{\color[rgb]{0,0,0}X}}^{u}\subseteq{\mathbb{I}\mathcal{{\color[rgb]{0,0,0}X}}^{p}}$
(i.e.,
$\underline{{\color[rgb]{0,0,0}x}}^{p}\leq\underline{{\color[rgb]{0,0,0}x}}^{u}$
and
$\overline{{\color[rgb]{0,0,0}x}}^{p}\geq\overline{{\color[rgb]{0,0,0}x}}^{u}$)
by initialization and construction (cf. Lines 2, 9 and 19). Further, we show
that
$\mathbb{I}\mathcal{{\color[rgb]{0,0,0}X}}^{u}\supseteq\mathbb{I}\mathcal{{\color[rgb]{0,0,0}X}}^{*}\triangleq\\{{\color[rgb]{0,0,0}x}\in\mathbb{I}\mathcal{{\color[rgb]{0,0,0}X}}^{p}\,|\,\underline{y}\leq\nu({\color[rgb]{0,0,0}x})\leq\overline{y}\\}$.
To use contradiction, suppose that it does not hold. Then,
$\exists\zeta\in\mathbb{I}\mathcal{{\color[rgb]{0,0,0}X}}^{*}$ such that
$\zeta\notin\mathbb{I}\mathcal{{\color[rgb]{0,0,0}X}}^{u}$, i.e.,
$\exists{\color[rgb]{0,0,0}i}\in\mathbb{N}_{n_{{\color[rgb]{0,0,0}x}}}$ such
that
$\zeta_{{\color[rgb]{0,0,0}i}}>\overline{{\color[rgb]{0,0,0}x}}^{u}_{{\color[rgb]{0,0,0}i}}$
or
$\zeta_{{\color[rgb]{0,0,0}i}}<\underline{{\color[rgb]{0,0,0}x}}^{u}_{{\color[rgb]{0,0,0}i}}$.
Without loss of generality, suppose the first case holds, i.e.,
$\zeta_{{\color[rgb]{0,0,0}i}}>\overline{{\color[rgb]{0,0,0}x}}^{u}_{{\color[rgb]{0,0,0}i}}$
(the proof for
$\zeta_{{\color[rgb]{0,0,0}i}}<\underline{{\color[rgb]{0,0,0}x}}^{u}_{{\color[rgb]{0,0,0}i}}$
is similar). Then,
$\zeta\in[\underline{{\color[rgb]{0,0,0}x}}^{m},\overline{{\color[rgb]{0,0,0}x}}^{p}]$,
where
$\underline{{\color[rgb]{0,0,0}x}}^{m}_{{\color[rgb]{0,0,0}i}}>\overline{{\color[rgb]{0,0,0}x}}^{u}_{{\color[rgb]{0,0,0}i}}$
and
$\underline{{\color[rgb]{0,0,0}x}}^{m}_{i^{\prime}}=\underline{{\color[rgb]{0,0,0}x}}^{p}_{i^{\prime}},\forall
i^{\prime}\neq i$. Hence,

$\displaystyle\underline{\nu}_{d}^{R}(\underline{{\color[rgb]{0,0,0}x}}^{m},\overline{{\color[rgb]{0,0,0}x}}^{p})\leq{\nu}(\zeta)\leq\overline{\nu}_{d}^{R}(\overline{{\color[rgb]{0,0,0}x}}^{p},\underline{{\color[rgb]{0,0,0}x}}^{m}),$
(68)

where $\overline{\nu}_{d}^{R}(\cdot,\cdot)$ and
$\underline{\nu}_{d}^{R}(\cdot,\cdot)$ are the proposed upper and lower
remainder-form decomposition functions in Algorithm 1. On the other hand, note
that
$\mathcal{{\color[rgb]{0,0,0}X}}^{u}\cap[\underline{{\color[rgb]{0,0,0}x}}^{m},\overline{{\color[rgb]{0,0,0}x}}^{p}]=\emptyset$,
hence the interval
$[\underline{{\color[rgb]{0,0,0}x}}^{m},\overline{{\color[rgb]{0,0,0}x}}^{p}]$
has been “ruled out” by Algorithm 2. In other words, one of the “or”
conditions in line 8 of Algorithm 2 must hold for this interval, i.e.,
$\overline{\nu}_{d}^{R}(\overline{{\color[rgb]{0,0,0}x}}^{p},\underline{{\color[rgb]{0,0,0}x}}^{m})<\underline{y}\
\lor\
\underline{\nu}_{d}^{R}(\underline{{\color[rgb]{0,0,0}x}}^{m},\overline{{\color[rgb]{0,0,0}x}}^{p})<\overline{y}$.
Combining this and (68), we obtain
$\nu(\zeta)<\underline{y}\lor{\mu}(\zeta)>\overline{y}$, which contradicts
with $\zeta\in\mathbb{I}{\mathcal{{\color[rgb]{0,0,0}X}}}^{*}$ (i.e.,
$\underline{y}\leq\nu(\zeta)\leq\overline{y}$). ∎

It is noteworthy that our set inversion algorithm can also be used with any
applicable inclusion functions (such as
$T^{f}_{N},T^{f}_{C},T^{f}_{M},T^{f_{d}}_{L},T^{f_{d}}_{O}$) or the best of
them (i.e., by independently computing the reachable sets of all inclusion
functions and intersecting them; cf. Corollary 1) in place of $T^{f_{d}}_{R}$
in Lines 7 and 17. On the other hand, the proposed $T^{f_{d}}_{R}$ (as well as
$T^{f_{d}}_{L},T^{f_{d}}_{O}$) can also be directly used in place of or in
combination with natural inclusions within SIVIA [18, Chapter 3] to obtain
subpavings  (i.e., a union of intervals).

##  IV Comparison with Existing Inclusion Functions

###  IV-A Comparison with the $T^{f_{d}}_{L}$ Inclusion Function

In this subsection, we compare the performance of the proposed $T^{f_{d}}_{R}$
with $T^{f_{d}}_{L}$ (cf. Proposition 5) through the following Theorem 5. We
show that the decomposition function $f^{L}_{d}$ introduced in [23] and
recapped in Proposition 5, belongs to the family of the remainder-form
decomposition functions in (17) and hence, $T^{f_{d}}_{L}$ cannot be tighter
than $T^{f_{d}}_{R}$, which is the tightest decomposition function that
belongs to (17).

###### Theorem 5 ($T^{f_{d}}_{L}$ vs $T^{f_{d}}_{R}$).

Suppose all the assumptions in Theorem 1 hold. Then, the following statements
are true.

  1. (i)

$f^{L}_{d}(\cdot,\cdot)$ belongs to the family of decomposition functions in
(17),  i.e., for each $i\in\mathbb{N}_{n_{x}}$, a specific pair of
$\mathbf{m}_{i}^{L}\in\mathbf{M}_{i}$ and
$h_{i}^{L}(\cdot)\in\mathcal{H}_{\mathbf{M}_{i}}$ corresponds to the
decomposition function $f^{L}_{d}(\cdot,\cdot)$ in [14, Theorem 2] (cf.
Proposition 5).

  2. (ii)

The optimal remainder-form decomposition function $T^{f}_{R}$ is always
tighter than (at least as good as) the inclusion function $T^{f}_{L}$, induced
by the decomposition function $f^{L}_{d}$.

###### Proof.

To prove (i), consider a specific decomposition function from the family of
remainder functions in (17) that is constructed, for each
$i\in\mathbb{N}_{n_{x}}$, with a supporting vector $\mathbf{m}_{i}^{L}$ and a
linear remainder function
${h}^{L}_{i}(\cdot)={\color[rgb]{0,0,0}\langle}{\mathbf{m}}^{L}_{i},\cdot{\color[rgb]{0,0,0}\rangle}$
as follows:

$\displaystyle(\mathbf{m}^{L}_{i})_{j}$ (69)
$\displaystyle=\begin{cases}\min(({\underline{J}^{f}_{C}})_{ij},0),&\text{if}\,|\min(({\underline{J}^{f}_{C}})_{ij},0)|\leq|\max(({\overline{J}^{f}_{C}})_{ij},0)|,\\\\[-0.71114pt]
\max(({\overline{J}^{f}_{C}})_{ij},0),&\text{if}\,|\max(({\overline{J}^{f}_{C}})_{ij},0)|<|\min(({\underline{J}^{f}_{C}})_{ij},0)|,\end{cases}$

for all $j\in\mathbb{N}_{n_{z}}$. Clearly,
$\mathbf{m}^{L}_{i}\in\mathbf{M}^{c}_{i}\subset\mathbf{M}_{i}$ by its
definition. Furthermore, it is easy to observe that $\mathbf{m}^{L}_{i}$ can
be rewritten as

$\displaystyle(\mathbf{m}^{L}_{i})_{j}=\begin{cases}0,&\text{if}\ (a_{ij}\geq
0)\lor(b_{ij}\leq 0)\lor(j=i),\\\ a_{ij},&\text{if}\
(a_{ij}<0)\land(b_{ij}>0)\land(|a_{ij}|\leq|b_{ij}|),\\\ b_{ij}&\text{if}\
(a_{ij}<0)\land(b_{ij}>0)\land(|b_{ij}|\leq|a_{ij}|),\end{cases}$ (70)

where $a_{ij}\triangleq({\underline{J}^{f}_{C}}_{ij})$,
$b_{ij}\triangleq({\overline{J}^{f}_{C}})_{ij}$. Recall that $(a_{ij}\geq 0)$,
$(a_{ij}<0)\land(b_{ij}>0)\land(|a_{ij}|\leq|b_{ij}|)$,
$(a_{ij}<0)\land(b_{ij}>0)\land(|b_{ij}|\leq|a_{ij}|)$, $(b_{ij}\leq 0)$ and
$j=i$ (continuous-time systems only) correspond to Cases $1$–$5$ in
Proposition 5, respectively. Then, by (69) and (70), we find that
$\zeta_{\mathbf{m}^{L}_{i},j}(z,\hat{z})$ in (19) coincides with $\zeta_{j}$
in Proposition 5. Moreover, by (70),

$\displaystyle\begin{array}[]{l}{h}^{L}_{i}(\zeta_{\mathbf{m}^{L}_{i}}(\hat{z},{z}))=\langle\mathbf{m}^{L}_{i},\zeta_{\mathbf{m}^{L}_{i}}(\hat{z},{z})\rangle=\sum_{j=1}^{\mathbb{N}_{n_{z}}}\phi^{i}_{j},\\\
{h}^{L}_{i}(\zeta_{\mathbf{m}^{L}_{i}}({z},\hat{z}))=\langle\mathbf{m}^{L}_{i},\zeta_{\mathbf{m}^{L}_{i}}({z},\hat{z})\rangle=\sum_{j=1}^{\mathbb{N}_{n_{z}}}\psi^{i}_{j},\end{array}$
(73)

where  $\phi^{i}_{j}\triangleq\begin{cases}0,&\text{Cases}\,1,4,5,\\\
a_{ij}\hat{z}_{j},&\text{Case}\ 2,\\\ b_{ij}z_{j},&\text{Case}\ 3,\end{cases}$
and $\psi^{i}_{j}\triangleq\begin{cases}0,&\text{Cases}\,1,4,5,\\\
a_{ij}{z}_{j},&\text{Case}\ 2,\\\ b_{ij}\hat{z}_{j},&\text{Case}\
3.\end{cases}$ Consequently,

$\displaystyle{h}^{L}_{i}(\zeta_{\mathbf{m}^{L}_{i}}(\hat{z},{z}))-{h}^{L}_{i}(\zeta_{\mathbf{m}^{L}_{i}}({z},\hat{z}))=\textstyle\sum_{j=1}^{\mathbb{N}_{n_{z}}}\phi^{i}_{j}-\psi^{i}_{j}=\sum_{j=1}^{\mathbb{N}_{n_{z}}}\theta^{i}_{j},$

where $\theta^{i}_{j}=\begin{cases}0,&\text{Cases}\,1,4,5,\\\
a_{ij}(\hat{z}_{j}-z_{j}),&\text{Case}\ 2,\\\
b_{ij}(z_{j}-\hat{z}_{j}),&\text{Case}\ 3.\end{cases}$ Then, defining two
indicator functions $\alpha^{i},\beta^{i}\in\mathbb{R}^{n_{z}}$, where for all
$j\in\mathbb{N}_{n_{z}}$,
$\alpha^{i}_{j}\triangleq\begin{cases}0,&\text{Cases}\,1,3,4,5,\\\
|a_{ij}|,&\text{Case}\ 2,\end{cases}$ and
$\beta^{i}_{j}\triangleq\begin{cases}0,&\text{Cases}\,1,2,4,5,\\\
-|b_{ij}|,&\text{Case}\ 3,\end{cases}$ $\theta^{i}_{j}$ can be rewritten as
$\theta^{i}_{j}=(\alpha^{i}_{j}-\beta^{i}_{j})(z_{j}-\hat{z}_{j})$, and hence

$\displaystyle{h}^{L}_{i}(\zeta_{\mathbf{m}^{L}_{i}}(\hat{z},{z}))-{h}^{L}_{i}(\zeta_{\mathbf{m}^{L}_{i}}({z},\hat{z}))=\textstyle\sum_{j=1}^{\mathbb{N}_{n_{z}}}\theta^{i}_{j}=\langle\alpha^{i}-\beta^{i},z-\hat{z}\rangle.$

Finally, since $\zeta_{\mathbf{m}^{L}_{i}}({z},\hat{z})$ coincides with
$\zeta$ in Proposition 5, by (17),
$f^{i}_{d}(z,\hat{z};\mathbf{m}^{L}_{i},h^{L}_{i}(\cdot))=h^{L}_{i}(\zeta_{\mathbf{m}^{L}_{i}}(\hat{z},{z}))+f_{i}(\zeta_{\mathbf{m}^{L}_{i}}({z},\hat{z}))-h_{i}(\zeta_{\mathbf{m}^{L}_{i}}({z},\hat{z})))=f_{i}(\zeta)+\langle\alpha^{i}-\beta^{i},z-\hat{z}\rangle=f^{L}_{d,i}(z,\hat{z})$,
where $f^{L}_{d,i}(z,\hat{z})$ is the decomposition function introduced in
Proposition 5 and is defined in (6).

(ii) The result directly follows from (i) and Theorem 2. ∎

From the above, we know that $T^{f_{d}}_{L}$, which only considers a specific
$\mathbf{m}_{i}\in\mathbf{M}_{i}$ for all $i\in\mathbb{N}_{n_{x}}$, cannot be
tighter than $T^{f_{d}}_{L}$, which considers all
$\mathbf{m}_{i}\in\mathbf{M}_{i}$. Nonetheless, since $T^{f_{d}}_{L}$ requires
less computation, it can still be useful for system s with large dimensions,
and can also be tighter than $T^{f}_{N},T^{f}_{C}$ and $T^{f}_{M}$ (see
Examples 1 and 2 below). Further, this suggests that when computational
resources are limited, it is also possible to consider a strict subset of
$\mathbf{M}_{i}$ on top of the one in $T^{f_{d}}_{L}$ to obtain a tighter
decomposition function than $T^{f_{d}}_{L}$.  In addition, the above theorem
indirectly proves that $T^{f_{d}}_{L}$ in Proposition 5 also applies to ELLC
systems.

###  IV-B Comparison with $T^{f}_{N}$, $T^{f}_{C}$ and $T^{f}_{M}$ Inclusion
Functions

In this subsection, we compare the performance of natural inclusions and some
of their modifications, i.e., $T^{f}_{N}$, $T^{f}_{C}$, $T^{f}_{M}$ with the
(discrete-time) $T^{f_{d}}_{R}$, via computing the over-approximation for the
range of some example functions. It is worth mentioning that we were not able
to derive any theoretical results that show the superiority of one over the
others. In fact, our simulation results showed that depending on function and
its corresponding considered domain, one of them can be tighter than the
others in some cases and the opposite holds for other cases. However, in some
cases, reflected in the following examples, the $T^{f_{d}}_{R}$ typically
returns tighter intervals.

####  IV-B1 Composition of Non-Elementary Functions

In cases where the considered vector field is not a composition of “elementary
functions” (e.g., simple monomials, $\sin(\cdot)$, $\cos(\cdot)$, monotone
functions, etc .), $T^{f}_{N}$, $T^{f}_{C}$ and $T^{f}_{M}$ are known to be
hard to compute and conservative over-approximations for bounding the
constituent functions are often needed, which lead to poor inclusion
functions, i.e., large errors. In these cases, it is most likely that
$T^{f_{d}}_{R}$ returns better bounds. The following example describes one
such function.

###### Example 1.

Consider $f(x)=x\arctan{(x^{2}-2x+5)}$, which is composed of non-elementary
functions, and an interval domain $\mathbb{I}\mathcal{X}=[1,3]$. In this case,
$T^{f}_{N}$, $T^{f}_{C}$, $T^{f}_{M}$, $T^{f_{d}}_{L}$ and $T^{f_{d}}_{R}$
return $[-4.7124,4.7124]$, $[1.3258,4.3393]$, $[1.3187,4.2475]$,
$[1.2835,2.9461]$ and $[1.1760,2,7468]$, respectively, where the final
interval (corresponding to $T^{f_{d}}_{R}$) is a subset of all others.

####  IV-B2 “Almost” Sign-Stable Functions

In cases where $f(\cdot)$ can be decomposed into a CJSS constituent and a
relatively small additive perturbation, $T^{f}_{R}$ will most likely return
tighter bounds  than the bounds returned by $T^{f}_{N}$, $T^{f}_{C}$ and
$T^{f}_{M}$. For instance, consider the following example.

###### Example 2.

Consider $f(x)=x^{3}-0.1x$, which is a monotone increasing (and hence CJSS)
function on its interval domain $\mathbb{I}\mathcal{X}=[-1,3]$, except on the
short interval $[-\sqrt{\frac{0.1}{3}},\sqrt{\frac{0.1}{3}}]$. For this
example, $T^{f}_{N}$, $T^{f}_{C}$, $T^{f}_{M}$, $T^{f_{d}}_{L}$ and
$T^{f_{d}}_{R}$ return
$[-8.9000,27.0100]$,$[-49.9000,54.7000]$,$[-49.9000,54.7000]$,
$[-1.0300,26.0100]$ and $[-1.0300,26.0100]$, respectively, where
$T^{f_{d}}_{L}$ and $T^{f_{d}}_{R}$ are much tighter than $T^{f}_{N}$,
$T^{f}_{C}$ and $T^{f}_{M}$.

####  IV-B3 Vector Fields with Several Additive Terms

It is also well-known in the literature that natural, centered and mixed-
centered inclusions perform worse for the functions with many additive terms,
compared to the ones with fewer additive terms [18, 20]. This is not
necessarily true for the performance of $T^{f_{d}}_{R}$. The following example
illustrates this fact, where a function with several additive terms is
considered.

###### Example 3.

Consider
$f(x)=x_{1}x_{2}x_{3}+x_{1}^{2}x_{2}+x_{2}^{2}x_{3}+x_{3}^{2}x_{1}+x_{1}^{2}x_{3}+x_{3}^{2}x_{2}+x_{2}^{2}x_{1}+x_{1}^{3}+x_{2}^{3}+x_{3}^{3}$
with an interval domain
$\mathbb{I}\mathcal{X}=[-2,2]\times[-2,2]\times[-2,2]$. Then, $T^{f}_{N}$,
$T^{f}_{C}$, $T^{f}_{M}$, $T^{f_{d}}_{L}$ and $T^{f_{d}}_{R}$ return
$[-80,80]$, $[-76.45,76.45]$, $[-73.62,73.62]$, $[-176,176]$ and
$[-54.4,54.4]$, respectively, where the final interval from $T^{f_{d}}_{R}$ is
the tightest among all.

####  IV-B4 Existence of Closed-Form Decomposition Functions

Finally, it is notable that our proposed $T^{f_{d}}_{R}$ approach (and
$T^{f_{d}}_{L}$) enables us to find closed-form inclusion functions for a wide
class of vector fields. This, can be analytically beneficial, e.g., in
convergence analysis for reachable sets or stability analysis in interval
observer designs [34, 35]. This is in contrast to natural, centered and mixed-
centered inclusions (and also $T^{f_{d}}_{O}$ in general), where a closed-form
inclusion function for general classes of functions, is often not available.

##  V Applications

###  V-A Application to Constrained Reachability Analysis

Consider the following constrained bounded-error system:

$\displaystyle
x_{t}^{+}=f(x_{t},u_{t},w_{t}),\quad\mu(x_{t},u_{t})\in\mathcal{Y}_{t},$ (74)

where $x_{t}^{+}\triangleq x_{t+1}$ if (74) is a discrete-time  system (with
sampling time $\delta t$) and $x_{t}^{+}\triangleq\dot{x}_{t}$ if (74) is a
continuous-time system, $x_{t}\in\mathbb{R}^{n_{x}}$ with
$x_{0}\in[\underline{x}_{0},\overline{x}_{0}]$ and
$u_{t}\in\mathbb{R}^{n_{u}}$ are state and known input signals,
$w_{t}\in[\underline{w},\overline{w}]\in\mathbb{I}\mathbb{R}^{n_{w}}$ is a
bounded process disturbance signals,
$\mathcal{Y}_{t}\subseteq[\underline{y}_{t},\overline{y}_{t}]\in\mathbb{IR}^{n_{y}}$
is the time-varying, _uncertain_ state interval constraint and
$f:\mathbb{R}^{n_{x}+n_{u}+n_{w}}\to\mathbb{R}^{n_{x}},\mu:\mathbb{R}^{n_{x}+n_{u}}\to\mathbb{R}^{n_{{\color[rgb]{0,0,0}\mu}}}$
are known vector fields. The following proposition shows how to apply
Algorithms 1–2, i.e., the mixed-monotone remainder-form decomposition function
construction and the set inversion algorithms, to compute  over-approximations
of the reachable sets /framers of the states for the system in (74).

###### Proposition 12.

Consider the system (74) with initial state
$x_{0}\in\mathbb{I}\mathcal{X}_{0}\triangleq[\underline{x}_{0},\overline{x}_{0}]$
and let  $f(z_{t})\triangleq\tilde{f}(x_{t},u_{t},w_{t})$ and
$\nu(x_{t})\triangleq\mu(x_{t},u_{t})$ with
$z_{t}\triangleq[x_{{\color[rgb]{0,0,0}t}}^{\top}\ w_{t}^{\top}]^{\top}$ and
known $u_{t}$. Suppose that Assumptions 1 and 2 hold and $\epsilon$ is a
chosen small positive threshold. Then, for all $t\geq
0,CR^{f}(t,\mathbb{I}\mathcal{X}_{0})\subset\mathbb{I}\mathcal{X}^{u}_{t}\triangleq[\underline{x}^{u}_{t},\overline{x}^{u}_{t}]$,
where
$CR^{f}(t,\mathbb{I}\mathcal{X}_{0})\triangleq\\{\phi(t,x_{0},{w}_{{\color[rgb]{0,0,0}t}})\mid
x_{0}\in\mathbb{I}\mathcal{X}_{0},\mu(x_{t},u_{t})\in{\color[rgb]{0,0,0}\mathcal{Y}_{t}\subseteq[\underline{y}_{t},\overline{y}_{t}]}\text{
and }w_{t}\in\mathbb{I}\mathcal{W},\forall t\geq 0\\}$ is the constrained
reachable set at time $t$ of (74) when initialized within
$\mathbb{I}\mathcal{X}_{0}$ and
$\mathbb{I}\mathcal{X}^{u}_{t}\triangleq[\underline{x}^{u}_{t},\overline{x}^{u}_{t}]$
is the solution to the following constrained embedding system:

$\displaystyle\begin{bmatrix}{\overline{x}}_{t}^{p+}\\\
{\underline{x}}_{t}^{p+}\end{bmatrix}$
$\displaystyle=\begin{bmatrix}\overline{{\color[rgb]{0,0,0}f}}_{d}(\begin{bmatrix}(\overline{x}_{t}^{u})^{\top}\,\overline{w}^{\top}\end{bmatrix}^{\top},\begin{bmatrix}(\underline{x}_{t}^{u})^{\top}\,\underline{w}^{\top}\end{bmatrix}^{\top})\\\
\underline{{\color[rgb]{0,0,0}f}}_{d}(\begin{bmatrix}(\underline{x}_{t}^{u})^{\top}\,\underline{w}^{\top}\end{bmatrix}^{\top},\begin{bmatrix}(\overline{x}_{t}^{u})^{\top}\,\overline{w}^{\top}\end{bmatrix}^{\top})\end{bmatrix},\
\begin{bmatrix}\overline{x}^{p}_{0}\\\
\underline{x}^{p}_{0}\end{bmatrix}=\begin{bmatrix}\overline{x}_{0}\\\
\underline{x}_{0}\end{bmatrix},$
$\displaystyle(\overline{x}^{u}_{t},\underline{x}^{u}_{t})$
$\displaystyle=\text{SET-
INV}({\nu}(\cdot),\overline{J}^{{\nu}}_{C},\underline{J}^{{\nu}}_{C},\overline{x}^{p}_{t},\underline{x}^{p}_{t},\overline{y}_{t},\underline{y}_{t},\epsilon),$

where
$(\overline{{\color[rgb]{0,0,0}f}}_{d}(\cdot,\cdot),\underline{{\color[rgb]{0,0,0}f}}_{d}(\cdot,\cdot))=T^{{\color[rgb]{0,0,0}f}_{d}}_{R}({\color[rgb]{0,0,0}f}(\cdot),\overline{J}^{{{\color[rgb]{0,0,0}f}}}_{C},\underline{J}^{{\color[rgb]{0,0,0}f}}_{C},\cdot,\cdot)$
is the discrete-time or continuous-time decomposition-based inclusion function
in Algorithm 1 and the SET-INV function in Algorithm 2 is based on the
discrete-time
$T^{\nu_{d}}_{R}(\nu(\cdot),\overline{J}^{{\nu}}_{C},\underline{J}^{{\nu}}_{C},\cdot,\cdot)$
from Algorithm 1. Consequently, the constrained system state trajectory
$x_{t}$ satisfies $\underline{x}^{u}_{t}\leq x_{t}\leq\overline{x}^{u}_{t}$ at
all times $t$.

###### Proof.

The results directly follow from applying  Propositions 3–4, Theorems 1–2 and
Lemmas 3–11. ∎

###  V-B Application to Interval Observer Design

Now, we consider the interval observer design problem for the following
bounded-error system:

$\displaystyle x^{+}=f(x_{t},u_{t},w_{t}),\quad
y_{t}={\mu}(x_{t},u_{t})+Vv_{t},$ (75)

where its state equation is similar to the system (74), but instead of a known
state constraint, an observation/measurement signal
$y_{t}\in\mathbb{R}^{n_{\color[rgb]{0,0,0}{y}}}$ is known/measured with a
known observation function
${\mu}(\cdot):\mathbb{R}^{n_{x}+n_{u}}\to\mathbb{R}^{n_{{\color[rgb]{0,0,0}y}}}$,
known matrix
$V\in\mathbb{R}^{n_{{\color[rgb]{0,0,0}y}}\times{\color[rgb]{0,0,0}n}_{v}}$
and measurement noise signal
$v_{t}\in[\underline{v},\overline{v}]\in\mathbb{IR}^{n_{v}}$. It can be easily
verified that the observation equation can be equivalently written as
${\nu}(x_{t})\triangleq{\mu}(x_{t},u_{t})\in\mathcal{Y}_{t}\triangleq[y_{t}-\overline{s},y_{t}-\underline{s}]$,
where $\overline{s}=V^{+}\overline{v}-V^{-}\underline{v}$,
$\underline{s}=V^{+}\underline{v}-V^{-}\overline{v}$,
$V^{+}\triangleq\max(V,\mathbf{0}_{n_{v}})$ and $V^{-}\triangleq V^{+}-V$,
with $\mathbf{0}_{n_{v}}$ being a zero vector in $\mathbb{R}^{n_{v}}$ [36].
This transforms the system (75) into the form of (74) and hence, Proposition
12 directly applies, where $[\underline{x}^{p}_{t},\overline{x}^{p}_{t}]$ is
the state interval from the prediction/propagation step and
$[\underline{x}^{u}_{t},\overline{x}^{u}_{t}]$ is the updated state interval
after a measurement update step. Further, if (75) is a sampled-data system,
i.e., the system dynamics is continuous and the observations are sampled in
discrete-time, our proposed approach still directly applies  with minor
modifications.

##  VI Simulations

In this section, we compare the performance s of $T^{f}_{N}$ (natural
inclusions; cf. Proposition 1), $T^{f}_{C},T^{f}_{M}$ (centered and mixed-
centered inclusions; cf. Proposition 2), $T^{f_{d}}_{L}$ (decomposition
functions proposed in [14]; cf. Proposition 5), $T^{f_{d}}_{R}$ (our proposed
remainder-form decomposition function in Algorithm 1), $T^{f}_{S_{1}}$ (the
first proposed bounding approach in [23,  Theorem 1], if applicable),
$T^{f}_{S_{2}}$ (the second proposed approach in [23,  Theorem 2], if
applicable) and $T^{f_{d}}_{O}$ (the tight decomposition functions proposed in
[12, Theorem 2] for discrete-time and in [13, Theorem 1] for continuous-time
systems, when computable; cf. Proposition 6) in computing the reachable sets
of several unconstrained and constrained dynamical systems in the form of
(15). Further, by the intersection property in Corollary 1, we can also
intersect the reachable sets of all applicable inclusion functions (except
$T^{f_{d}}_{O}$) for comparison. Note that we only compare methods in the
literature that use _first-order_ (generalized) gradient information. The
consideration of higher-order information is a subject of our future work.

###  VI-A Van Der Pol System

Consider the following discretized Van der Pol system [27]:

$\displaystyle\begin{array}[]{ll}x_{1,t+1}=x_{1,t}+\delta t\,x_{2,t},\\\
x_{2,t+1}=x_{2,t}+\delta t\,((1-x^{2}_{1,t})x_{2,t}-x_{1,t}),\end{array}$ (78)

with $\delta t=0.1$ s and $x_{0}\in[1.15,1.4]\times[2.05,2.3]$. Starting from
the initial intervals, the results for computing reachable intervals, using
several applicable inclusion functions, are depicted in  Figure 1.
Unfortunately, $T^{f}_{S_{1}}$ and $T^{f}_{S_{2}}$ are not valid (applicable)
here, due to the lack of required monotonicity assumptions (cf. [23,
conditions (6) and (16)]), while as expected, $T^{f_{d}}_{O}$ (that is
computable here, since the corresponding optimization problems can be
analytically solved by computing and comparing all stationary points of the
vector fields and the boundaries) returns the tightest bounds. Moreover, our
proposed $T^{f_{d}}_{R}$ notably returns tighter bounds than all applicable
methods, except $T^{f_{d}}_{O}$ and the best of $T^{f}_{N}$–$T^{f_{d}}_{R}$,
and seems to offer significant improvement over natural, centered and mixed-
centered inclusions as well as $T^{f_{d}}_{L}$. Further, “the best of
$T^{f}_{N}$–$T^{f_{d}}_{R}$” approach is able to get very close to the
tightest set obtained by $T^{f_{d}}_{O}$.

Figure 1: Upper and lower bounds on ${x}_{1}$ and $x_{2}$ in System (78)
(Example A: Van der Pol system), when applying $T^{f}_{N}(--)$, $T^{f}_{C}$
($\circ$), $T^{f}_{M}$ ($\diamond$), $T^{f_{d}}_{L}$ ($\square$),
$T^{f_{d}}_{R}$ ($\ast$), the best of $T^{f}_{N}$–$T^{f_{d}}_{R}$ ($\cdot$-)
and $T^{f_{d}}_{O}$ (+), as well as the midpoint trajectory (–).

###  VI-B Example 3 in [23] with Modeling Redundancy

(a) System B without observations (79)–(81).

(b) System B with observations (79)–(82).

Figure 2: Upper and lower bounds/framers for ${x}_{1}$ and $x_{2}$ in Example
B, when applying $T^{f}_{N}(--)$, $T^{f}_{S_{1}}$ ($\triangle$), $T^{f}_{C}$
($\circ$), $T^{f}_{M}$ ($\diamond$), $T^{f_{d}}_{L}$ ($\square$),
$T^{f_{d}}_{R}$ ($\ast$), $T^{f_{d}}_{R}$ with manufactured redundancies
($\times\cdot$) and $T^{f_{d}}_{O}$ (+), as well as the midpoint trajectory
(–).

Now consider the following discrete-time dynamical system from [24, 23] with
bounded noise:

$\displaystyle x_{t+1}=\begin{bmatrix}0&-0.5\\\
1&1+0.3{\color[rgb]{0,0,0}\lambda}_{t}\end{bmatrix}x_{t}+0.02\begin{bmatrix}-6\\\
1\end{bmatrix}w_{t},$ (79)

where $w_{t}$, ${\color[rgb]{0,0,0}\lambda}_{t}\in[-0.001,0.001]$ and
$x_{0}\in[-0.55,-0.445]\times[0.145,0.248]$. The approximated reachable sets
are depicted in  Figure 2(a). Here, $T_{S_{1}}$ is applicable and returns the
exact same bounds as $T^{f}_{N}$ (natural inclusion), but $T^{f}_{S_{2}}$ is
not applicable due to the lack of monotonicity (cf. [23,  (16)]). The tight
$T^{f_{d}}_{O}$ is again computable and, as expected, returns the tightest
possible intervals. Again, our proposed approach, $T^{f_{d}}_{R}$ outperforms
all applicable ones, except for the tight $T^{f_{\color[rgb]{0,0,0}d}}_{O}$.
To further improve our results and inspired by the literature on modeling
redundancy, e.g., [26, 23], we consider a manufactured _redundant_ state
$z_{t}$ chosen as:

$\displaystyle z_{t}$ $\displaystyle=x_{1,t}+6x_{2,t},$ (80)

which implies that:

$\displaystyle z_{t+1}$ $\displaystyle=z_{t}+5x_{1,t}+(1.8v_{t}-0.5)x_{2,t}.$
(81)

Augmenting (81) with the original system (79), we consider the over-
approximation of the reachable sets for the augmented system, subject to (80)
as a constraint, using our proposed $T^{f_{d}}_{R}$ for the propagation step,
as well as Algorithm 2 for set inversion (refinement/update). The results in
Figure 2(a) show an improvement when considering the manufactured redundant
variable  when compared to using $T^{f_{d}}_{R}$ without any redundancies. It
can also be shown that different manufactured states result in different
amounts of improvement. Thus, as future work, we will also consider the
development of a _principled_ approach for choosing manufactured/redundant
variables.

Next, we consider (79)–(81) with an observation equation:

$\displaystyle y_{t}=1.6x_{1,t}+0.3x_{2,t}+v_{t},$ (82)

where $y_{t}$ is a _known_ measurement/observation signal at time step $t$
with measurement noise signal $v_{t}\in[-.05,0.05]$.  Figure 2(b) shows the
approximated upper and lower bounds for the states of the system (82), using
the same inclusion functions that we applied to the system (79), as well as
applying Algorithm 2 for the set inversion/refinement procedure. As expected
and as can be seen by comparing Figures 2(a) and 2(b), the additional
measurement information enables us to tighten the resulting intervals for all
inclusion functions.

###  VI-C Example 2.11 in [18]

Next, we consider the following discrete-time dynamical system in [18, Example
2.11]:

$\displaystyle\begin{array}[]{ll}x_{1,t+1}=x_{1,t}^{2}+x_{1,t}e^{x_{2,t}}-x_{2,t}^{2},\\\
x_{2,t+1}=x_{1,t}^{2}-x_{1,t}e^{x_{2,t}}+x_{2,t}^{2},\end{array}$ (85)

where $x_{0}\in[0.12,0.121]\times[0.182,0.185]$. Here, $T^{f}_{S_{1}}$ and
$T^{f}_{S_{2}}$ approaches are not applicable, due to lack of the required
monotonicity conditions, while the tight $T^{f_{d}}_{O}$ is not computable
since the stationary/critical points are not  analytically solvable. On the
other hand, our proposed remainder-form inclusion function, $T^{f_{d}}_{R}$,
can be tractably computed.  Figure 3 depicts the resulting reachable sets,
where we observe that $T^{f_{d}}_{R}$ outperforms all other applicable and
computable approaches, and  is comparable to the “best of all” approach (via
intersection).

Figure 3: Upper and lower bounds on ${x}_{1}$ and $x_{2}$ in System (85)
(Example C), when applying $T^{f}_{N}(--)$, $T^{f}_{C}$ ($\circ$), $T^{f}_{M}$
($\diamond$), $T^{f_{d}}_{L}$ ($\square$), $T^{f_{d}}_{R}$ ($\ast$) and the
best of $T^{f}_{N}$–$T^{f_{d}}_{R}$ ($\cdot$-), as well as the midpoint
trajectory (–).

###  VI-D Continuous-Time System in [13]

As the next example, we consider the following continuous-time dynamical
system from [13]:

$\displaystyle\begin{bmatrix}\dot{x}_{1,t}\\\ \dot{x}_{2,t}\\\
\dot{x}_{3,t}\end{bmatrix}=\begin{bmatrix}w_{1{\color[rgb]{0,0,0},t}}x_{2,t}^{2}-x_{2,t}+w_{2,t}\\\
x_{3,t}+2\\\ x_{1,t}-x_{2,t}-w_{1,t}^{3}\end{bmatrix},$ (86)

with ${x}_{0}\in[-\frac{1}{2},\frac{1}{2}]^{3}$ and
$w_{t}\in[-\frac{1}{4},0]\times[0,\frac{1}{4}]$. Figure 4 depicts the
approximations of the reachable sets when applying $T^{f}_{N}$, $T^{f}_{C}$,
$T^{f}_{M}$, $T^{f_{d}}_{L}$, $T^{f_{d}}_{R}$, the best of
$T^{f}_{N}$–$T^{f_{d}}_{R}$ and $T^{f_{d}}_{O}$. As expected, the tight
$T^{f_{d}}_{O}$, which is implemented using the corresponding embedding
functions given in [13, Section VI], returns the tightest sets/framers.
Further, among the other inclusion functions, $T^{f_{d}}_{R}$ has the the best
performance and it  can be slightly improved  by using the best of
$T^{f}_{N}$–$T^{f_{d}}_{R}$.

Figure 4: Upper and lower bounds for ${x}_{1}$, $x_{2}$ and $x_{3}$ in System
(86) (Example D), when applying $T^{f}_{N}(--)$, $T^{f}_{C}$ ($\circ$),
$T^{f}_{M}$ ($\diamond$), $T^{f_{d}}_{L}$ ($\square$), $T^{f_{d}}_{R}$
($\ast$), the best of $T^{f}_{N}$–$T^{f_{d}}_{R}$ ($\cdot$-) and
$T^{f_{d}}_{O}$ (+), as well as the midpoint trajectory (–).

###  VI-E Unicycle System

Next, we are interested in computing the reachable sets for a well-known
continuous-time system, namely the unicycle-like mobile robot, e.g., in [37,
38], with two driving wheels, mounted on the left and right sides of the
robot, with their common axis passing through the center of the robot. The
dynamics of such a system can be described as in [38]:

$\displaystyle\begin{array}[]{c}\dot{s}_{x,t}=\phi_{\omega,t}\cos\theta_{t}+w_{x,t},\
\dot{s}_{y,t}=\phi_{\omega,t}\sin\theta_{t}+w_{y,t},\\\
\dot{\theta}_{t}=\phi_{\theta,t}+w_{\theta,t},\end{array}$ (89)

with state $x_{t}\triangleq[s_{x,t}\ s_{y,t}\ \theta_{t}]^{\top}$, where
$s_{x,t}$ and $s_{y,t}$ are the $X$ and $Y$ coordinates of the main axis mid-
point between the two driving wheels and $\theta_{t}$ is the angle between the
robot forward axis and the $X$-direction, while $\phi_{\omega,t}$ and
$\phi_{\theta,t}$ are the displacement and angular velocities of the robot (as
known inputs), respectively, and $w_{t}=[w_{x,t}\ w_{y,t}\
w_{\theta,t}]^{\top}$ is the process noise vector. Setting
$\phi_{\omega,t}=0.3$, $\phi_{\theta,t}=0.15$,
$w_{x,t}=0.2(0.5\rho_{x_{1,t}}-0.3)$, $w_{y,t}=0.2(0.3\rho_{x_{2,t}}-0.2)$ and
$w_{\theta,t}=0.2(0.6\rho_{x_{3,t}}-0.4)$, with $\rho_{x_{l,t}}\in[0,1]$
$(l=1,2,3)$ and initial state $x_{0}=[0.1\ 0.2\ 1]^{\top}$, Figure 5(a) shows
the over-approximations of the reachable sets for system (89), using the
various methods discussed in this paper. As can be observed, $T^{f_{d}}_{O}$,
which is computable for this example, returns the tightest intervals.
Moreover, in some intervals, natural inclusions and their modifications,
return tighter bounds than $T^{f_{d}}_{R}$. However, by taking their
intersection and computing the best of $T^{f}_{N}$–$T^{f_{d}}_{R}$, we obtain
further improvements.

Then, to illustrate the effectiveness of our proposed set inversion algorithm
for continuous-time systems, we consider observations/measurements similar to
[37], as follows. In the $X$–$Y$ plane, two known points, denoted as
$(s_{x_{i}},s_{y_{i}})$ $(i=1,2)$, are chosen as markers/landmarks. Then, the
distance from the robot’s planar Cartesian coordinates $(s_{x,t},s_{y,t})$ to
each marker $(s_{x_{i}},s_{y_{i}})$ can be expressed as
$d_{i,t}=\sqrt{(s_{x_{i}}-s_{x,t})^{2}+(s_{y_{i}}-s_{y,t})^{2}}$. Furthermore,
the azimuth $\phi_{i,t}$ at time $t$ can be related to the current system
state variables $s_{x,t},s_{y,t}$ and $\theta_{t}$ as
$\phi_{i,t}=\theta_{t}-\arctan(\frac{s_{y_{i}}-s_{y,t}}{s_{x_{i}}-s_{x,t}})$.
Treating both the distance $d_{i,t}$ and the azimuth $\phi_{i,t}$ as the
measurements, as well as considering measurement disturbances /noise $v_{t}$,
the  nonsmooth, nonlinear observation/measurement equation can be written as:

$\displaystyle y_{t}=[d_{1,t}\ \phi_{1,t}\ d_{2,t}\ \phi_{2,t}]^{\top}+v_{t},$
(90)

with $v_{1,t}=0.02\rho_{{y_{1}},t}-0.01$, $v_{2,t}=0.03\rho_{{y_{2}},t}-0.01$,
$v_{3,t}=0.03\rho_{{y_{3}},t}-0.02$, $v_{4,t}=0.05\rho_{{y_{4}},t}-0.03$ and
$\rho_{{y_{k}},t}\in[0,1]$ $(k=1,2,3,4)$. Now, applying all the methods along
with the set inversion approach in Algorithm 2 to the constrained system
(89)–(90), one can observe that by taking advantage of the observations, the
reachable set approximations have been significantly improved in Figure 5(b)
(with observations), when compared to Figure 5(a) (without observations).

(a) Unicycle system without observation (89).

(b) Unicycle system with observations (89)–(90).

Figure 5: Upper and lower bounds on ${x}_{1}$ and $x_{2}$ in the unicycle
system (Example E), when applying $T^{f}_{N}(--)$, $T^{f}_{C}$ ($\circ$),
$T^{f}_{M}$ ($\diamond$), $T^{f_{d}}_{L}$ ($\square$), $T^{f_{d}}_{R}$
($\ast$), the best of $T^{f}_{N}$–$T^{f_{d}}_{R}$ ($\cdot$-) and
$T^{f_{d}}_{O}$ (+), as well as the midpoint trajectory (–).

(a) GTM system without observations (94)

(b) GTM system with observations (94)–(95).

Figure 6: Upper and lower bounds on ${x}_{1}$ and $x_{2}$ in the GTM system
(Example F), when applying $T^{f}_{N}(--)$, $T^{f}_{C}$ ($\circ$), $T^{f}_{M}$
($\diamond$), $T^{f}_{L}$ ($\square$), $T^{f}_{R}$ ($\ast$) and the best of
$T^{f}_{N}$–$T^{f}_{R}$ ($\cdot$-), as well as the midpoint trajectory (–).

###  VI-F Generic Transport Longitudinal Model

Finally, we consider NASA’s Generic Transport Model (GTM) [39], a remote-
controlled 5.5% scale commercial aircraft [40], with the following main
parameters: wing area $S=5.902\ \text{ft}^{2}$, mean aerodynamic chord
$\overline{c}=0.9153$ ft, mass $m=1.542$ slugs, pitch axis moment of inertia
$I_{yy}=4.254\ \text{slugs/ft}^{2}$, air density $\rho=0.002375\
\text{slugs/ft}^{3}$ and gravitational acceleration $g=32.17\
\text{ft/s}^{2}$. The longitudinal dynamics of the GTM can be described as the
following continuous-time dynamical system:

$\displaystyle\begin{array}[]{ll}\dot{V}_{t}&=\frac{-D_{t}-mg\sin(\theta_{t}-\alpha_{t})+T_{x,t}\cos\alpha_{t}+T_{z,t}\sin\alpha_{t}}{m},\\\
\dot{\alpha}_{t}&=q_{t}+\frac{-L_{t}+mg\cos(\theta_{t}-\alpha_{t})-T_{x,t}\sin\alpha_{t}+T_{z,t}\cos\alpha_{t}}{mV_{t}},\\\
\dot{q}_{t}&=\frac{M_{t}+T_{m,t}}{I_{yy}},\
\dot{\theta}_{t}=q_{t},\end{array}$ (94)

where $V_{t}$, $\alpha_{t}$, $q_{t}$ and $\theta_{t}$ are air speed (ft/s),
angle of attack (rad), pitch rate (rad/s) and pitch angle (rad), respectively.
Moreover, $T_{x,y}$ (lbs), $T_{z,t}$ (lbs), $T_{m,t}$ (lbs–ft), $D_{t}$ (lbs),
$L_{t}$ (lbs) and $M_{t}$ (lbs) denote the projection of the total engine
thrust along the body’s  $X$-axis, the projection of the total engine thrust
along the body’s  $Z$-axis, the pitching moment due to both engines, the drag
force, the lift force and the aerodynamic pitching moment, respectively, with
their nominal values given in [41]. Defining $x_{t}\triangleq[V_{t}\
\alpha_{t}\ q_{t}\ \theta_{t}]^{\top}$ with
${x}_{0}\in[147,158]\times[0.04,0.05]\times[0.1,0.2]\times[0.04,0.05]$, Figure
6(a) depicts the reachable set approximations for $x_{1,t}$ and $x_{2,t}$ of
the system (94). For this system , $T^{f_{d}}_{O}$ is not computable, since
the stationary/critical points of the vector fields cannot be obtained
analytically, and as shown in Figure 6(a), $T^{f_{d}}_{R}$ obtains a tighter
over-approximation than $T^{f}_{N},T^{f}_{C}$, $T^{f}_{M}$ and
$T^{f_{d}}_{L}$, with the best of $T^{f}_{N}$–$T^{f_{d}}_{R}$ showing further
improvement.

Next, we consider an additional set of measurements, in the form of a linear
observation equation:

$\displaystyle y_{t}=x_{1,t}+x_{2,t}-x_{3,t}+v_{t},$ (95)

with $v_{t}\in[-0.01,0.01]$. Then applying $T^{f}_{N}$–$T^{f_{d}}_{R}$ along
with the set inversion (update) approach in Algorithm 2 to the constrained
system (94)–(95), we observe considerably tighter approximations for all
approaches in Figure 6(b) (with observations) when compared to the
approximations of the reachable sets in Figure 6(a) (without observations).

##  VII Conclusion And Future Work

A tractable family of remainder- form mixed-monotone decomposition functions
was proposed in this paper for a relatively large class of nonsmooth systems
called either-sided locally Lipschitz systems that is proven to include
locally Lipschitz continuous systems. We also characterized the lower and
upper bounds for the over-approximation errors when using the proposed
remainder-form decomposition functions to over-approximate the true
range/image set of a nonsmooth, nonlinear mapping, where the lower bound is
achieved by our proposed tight and tractable decomposition function, which was