Efficient Interaction-Aware Interval Analysis of Neural Network Feedback
  Loops
In this paper, we propose a computationally efficient framework for interval
reachability of systems with neural network controllers. Our approach leverages
inclusion functions for the open-loop system and the neural network controller
to embed the closed-loop system into a larger-dimensional embedding system,
where a single trajectory over-approximates the original system's behavior
under uncertainty. We propose two methods for constructing closed-loop
embedding systems, which account for the interactions between the system and
the controller in different ways. The interconnection-based approach considers
the worst-case evolution of each coordinate separately by substituting the
neural network inclusion function into the open-loop inclusion function. The
interaction-based approach uses novel Jacobian-based inclusion functions to
capture the first-order interactions between the open-loop system and the
controller by leveraging state-of-the-art neural network verifiers. Finally, we
implement our approach in a Python framework called ReachMM to demonstrate its
efficiency and scalability on benchmarks and examples ranging to $200$ state
dimensions.
Saber Jafarpour, Akash Harapanahalli, Samuel Coogan
# Efficient Interaction-Aware Interval Analysis of Neural Network Feedback Loops

###### Abstract

In this paper, we propose a computationally efficient framework for interval reachability of systems with neural network controllers. Our approach leverages inclusion functions for the open-loop system and the neural network controller to embed the closed-loop system into a larger-dimensional embedding system, where a single trajectory over-approximates the original system's behavior under uncertainty. We propose two methods for constructing closed-loop embedding systems, which account for the interactions between the system and the controller in different ways. The interconnection-based approach considers the worst-case evolution of each coordinate separately by substituting the neural network inclusion function into the open-loop inclusion function. The interaction-based approach uses novel Jacobian-based inclusion functions to capture the first-order interactions between the open-loop system and the controller by leveraging state-of-the-art neural network verifiers. Finally, we implement our approach in a Python framework called ReachM to demonstrate its efficiency and scalability on benchmarks and examples ranging to \(200\) state dimensions.

## I Introduction

Problem Description and MotivationsRecent advances in machine learning are bringing learning algorithms into the heart of safety-critical control systems, such as autonomous vehicles, manufacturing sectors, and robotics. For control systems, these learning algorithms often act in the closed-loop setting, as direct feedback controllers [1] or motion planners [2]. Learning-based closed-loop controllers can improve the performance of systems while reducing their computational burden for online implementations as compared with more traditional optimization-based approaches. However, despite their impressive performance, learning models are known to be overly sensitive to input disturbances [3]: a small input perturbation can lead to comparatively large changes in their output. This issue amplifies in feedback settings, as disturbances can accumulate in the closed-loop. As a result, ensuring reliability of learning algorithms is an essential challenge in their integration into safety-critical systems.

Typical learning architectures involve high-dimensional nonlinear function approximators, such as neural networks, necessitating special tools for their input-output analysis. The machine learning and control community have made significant progress in analyzing the safety of neural networks in isolation, including efficient and sound input-output bounds, worst-case adversarial guarantees, and sampling-based stochastic guarantees (cf. [4]). However, most of these existing frameworks for standalone neural networks do not address the unique challenges for closed-loop analysis--namely information propagation, non-stationarity of the bounds, and complex interactions between the system and the controller. In these cases, it is essential to understand and capture the nature of the interactions between the system and the neural network. Recently, several frameworks have emerged for safety verification of learning algorithms in the closed-loop. These frameworks usually incorporate neural network verification algorithms into the closed-loop safety analysis by studying their interactions with the open-loop system. However, most of these existing frameworks are either limited to a specific class of systems or learning algorithms, or they impose significant computational burdens, rendering them unsuitable for runtime safety verification.

Literature ReviewSafety verification of nonlinear dynamical systems without learning-enabled systems is a classical problem and is typically solved using reachability analysis tools, although significant challenges remains even in this setting. Reachability of nonlinear systems has been studied using optimization-based methods such as the Hamilton-Jacobi approach [5] and the level set approach [6]. Several computationally tractable approaches including the ellipsoidal method [7] and the zonotope method [8] have been developed for reachability analysis of dynamical systems. We refer to [9] for a recent review of state-of-the-art reachability analysis techniques. Interval analysis is a classical framework for computing interval functional bounds [10] which has been successfully used for reachability analysis of dynamical systems [11, 12, 13]. A closely-related body of literature is the mixed monotone theory, which extends the classical monotone system theory by separating the cooperative and competitive effect of the dynamics [14, 15]. Recently, mixed monotone theory has been used as an efficient framework for reachability analysis of dynamical systems [16, 17, 18, 19].

Rigorous verification of standalone neural networks has been studied using abstract interpretation techniques such as Reluplex [20] and Neurify [21], using interval bound propagation methods [22, 23], and convex-relaxation approaches such as LipSDP [24] and CROWN [25] and its variants [26]. The simplest method for studying reachability of neural network controlled systems is via an _input-output approach_: using existing techniques for estimating the output of the neural network, the input-output approach performs reachability analysis of the closed-loop system by substituting the full output range of the neural network as the input set
of the open-loop system. Examples of this approach include NNV [27], and simulation-guided interval analysis [28]. This approach is generally computationally efficient and applicable to general forms of neural networks and control systems, but it can lead to overly-conservative estimates of reachable sets [29, Section 2.1]. The reason for this over-conservatism is that this approach essentially treats neural networks controllers as adversaries and, thus, cannot capture the beneficial interactions between the controller and the system. An alternative framework is the _functional approach_, which is based on composing function approximations of the system and the neural network controller. Examples of this approach for linear systems include using linear programming in ReachLP [30], using semi-definite programming in Reach-SDP [31], and using mixed integer programming in [32]. For nonlinear system, the functional approach has been used in ReachNN [33], Sherlock [29], Verisig 2.0 [34], POLAR [35] and JuliaReach [36, 37]. Functional methods are able to capture beneficial interactions between the neural network controller and the system, but they are often computationally burdensome and generally do not scale well to large-scale systems.

ContributionsIn this paper, we introduce a general and computationally efficient framework for studying reachability of continuous-time closed-loop systems with nonlinear dynamics and neural networks controllers. First, we review the notion of inclusion function from classical interval analysis to over-approximate the input-output behavior of functions. Our first minor result states that if an inclusion function is chosen to be minimal, then it provides tighter over-approximations of the output behavior than a Lipschitz bound of the function. For a given function, we study different approaches for constructing inclusion functions. In particular, we propose a novel Jacobian-based cornered inclusion function, which is specially amenable to analysis of closed-loop systems. We further show that this class of inclusion functions plays a critical role in integrating the existing off-the-shelf neural network verification algorithms into our framework.

Second, using the notion of inclusion functions, we develop a general approach for interval analysis of continuous-time dynamical systems. The key idea is to use the inclusion functions to embed the original dynamical system into a higher dimensional embedding system. Then trajectories of the embedding system can be used to provide hyper-rectangular over-approximation for reachable sets of the original system. Our approach unifies and extends several existing approaches for interval analysis of dynamical systems. Noteably, our novel proof technique, based on the classical monotone comparison Lemma, allows for accuracy comparison between different inclusion functions.

As the culminating contribution of this paper, we develop a computationally efficient framework for reachability analysis of nonlinear systems with neural network controllers. The key in our framework is a careful combination of the inclusion function for the open-loop system and the inclusion function for the neural network controller to obtain an inclusion function for the closed-loop system. We propose two methods for combining the open-loop inclusion function and the neural network inclusion function. The first method, called the interconnection-based approach, constructs a closed-loop embedding system using the feedback interconnection of the open-loop embedding system and the inclusion function for the neural network. This method is general, as it can be applied to any open-loop and neural network inclusion function, and is computationally efficient, as it uses the compositional property of inclusion functions. The second method, called the interaction-based approach, constructs a closed-loop embedding system using a Jacobian-based inclusion function for the open-loop system and local affine bounds for the neural network. Compared to the interconnection-based approach, the interaction-based method completely captures the first order interactions between the system and the neural network.

Finally, we implement our approach in a Python toolbox called ReachMM and perform several numerical experiments to show the efficiency and accuracy of our framework. For several linear and nonlinear benchmarks, we show that our method beats the state-of-the-art methods. Moreover, we study scalability of our approach and compare it with the state-of-the-art methods on a platoon of vehicles with up to 200 state variables. Compared to our conference paper [38], this paper provides a more general framework for interval reachability using inclusion functions, introduces the interaction-based approach to capture the first order interactions between the system and the neural network controller, and includes a suite of numerical experiments. We would like to highlight that, in this paper, we focus on continuous-time dynamical systems. However, with minor modifications, our framework can be used to analyze reachability of discrete-time systems.

## II Mathematical Preliminary and Notations

For every two vectors \(v,w\in\mathbb{R}^{n}\) and every subset \(S\subseteq\{1,\ldots,n\}\), we define the vector \(v_{[S:w]}\in\mathbb{R}^{n}\) by \(\big{(}v_{[S:w]}\big{)}_{j}=\begin{cases}v_{j}&j\not\in S\\ w_{j}&j\in S\end{cases}\). Given a matrix \(B\in\mathbb{R}^{n\times m}\), we denote the non-negative part of \(B\) by \([B]^{+}=\max(B,0)\) and the nonpositive part of \(B\) by \([B]^{-}=\min(B,0)\). Given a square matrix \(A\in\mathbb{R}^{n\times n}\) the Metzler and non-Metzler part of \(A\) are denoted by \([A]^{\mathrm{M}}\) and \([A]^{\mathrm{mM}}\), respectively, where \(([A]^{\mathrm{M}})_{ij}=\begin{cases}A_{ij}&A_{ij}\geq 0\text{ or }i=j\\ 0&\text{otherwise,}\end{cases}\) and \([A]^{\mathrm{mM}}=A-[A]^{\mathrm{M}}\). Given a map \(g:\mathbb{R}^{n}\to\mathbb{R}^{m}\), the \(\ell_{\infty}\)-norm Lipschitz bound of \(g\) is the smallest \(L\in\mathbb{R}_{\geq 0}\) such that

\[\|g(x)-g(y)\|_{\infty}\leq L\|x-y\|_{\infty},\qquad\text{ for all }x,y\in\mathbb{R}^{n}.\]

We denote the \(\ell_{\infty}\)-norm Lipschitz constant of \(g\) by \(\mathrm{Lip}_{\infty}(g)\). For a differentiable map \(g:\mathbb{R}^{n}\to\mathbb{R}^{m}\), we denote its Jacobian at point \(x\in\mathbb{R}^{n}\) by \(Dg(x)\). Consider the dynamical system

\[\dot{x}=f(x,w),\qquad\qquad\text{ for all }t\in\mathbb{R}_{\geq 0}, \tag{1}\]

where \(x\in\mathbb{R}^{n}\) is the state of the system and \(w\in\mathbb{R}^{q}\) is the disturbance. Given a piecewise continuous disturbance \(t\mapsto w(t)\) and an initial condition \(x_{0}\in\mathbb{R}^{n}\), the trajectory of (1) starting from \(x_{0}\) is denoted by \(\phi_{f}(t,x_{0},w(t))\). Let
and \(\mathcal{W}\subseteq\mathbb{R}^{q}\), then the reachable set of (1) starting from the initial set \(\mathcal{X}\) and with the disturbance set \(\mathcal{W}\) is defined by:

\[\mathcal{R}_{f}(t,\mathcal{X},\mathcal{W})=\left\{\begin{aligned} &\phi_{f}(t,t_{0},x_{0},w),\,\forall x_{0}\in \mathcal{X},\\ & w:\mathbb{R}\rightarrow\mathcal{W}\text{ piecewise cont.}\end{aligned}\right\} \tag{2}\]

The set \(\mathcal{X}\) is _robustly forward invariant_ if \(\mathcal{R}_{f}(t,\mathcal{X},\mathcal{W})\subseteq\mathcal{X}\) for all \(t\geq 0\). The system (1) is continuous-time monotone on \(\mathcal{X}\subseteq\mathbb{R}^{n}\) if, for every \(i\in\{1,\ldots,n\}\), every \(x\leq y\in\mathcal{X}\) with \(x_{i}=y_{i}\), and every \(u\leq v\), we have \(f_{i}(x,u)\leq f_{i}(y,v)\). One can show that if a dynamical system (1) is monotone on \(\mathcal{X}\) and \(\mathcal{X}\) is robustly forward invariant, then its trajectories preserve the standard partial order by time, i.e., for every two trajectories \(t\mapsto x(t)\) and \(t\mapsto y(t)\) of the systems (1), if \(x(0)\leq y(0)\), then \(x(t)\leq y(t)\) for all \(t\in\mathbb{R}_{\geq 0}\)[39]. The _southeast_ partial order \(\leq_{\mathrm{SE}}\) on \(\mathbb{R}^{2n}\) is defined by \(\left[\frac{x}{2}\right]\leq_{\mathrm{SE}}\left[\frac{y}{y}\right]\) if and only if \(x\leq y\) and \(\tilde{y}\leq\tilde{x}\). We also define \(\mathcal{T}_{\geq 0}^{2n}=\{[\left\lceil\frac{x}{2}\right\rceil\in\mathbb{R}^{2n} \mid x\geq\tilde{x}]\text{ and }\mathcal{T}_{\leq 0}^{2n}=\{[\frac{x}{2}]\in \mathbb{R}^{2n}\mid x\geq\tilde{x}\}\text{ and }\mathcal{T}^{2n}=\mathcal{T}_{\geq 0}^{2n} \bigcup\mathcal{T}_{\leq 0}^{2n}\).

## III Problem Statement

We consider a system described by:

\[\dot{x}=f(x,u,w),\qquad\qquad\qquad t\in\mathbb{R}_{\geq 0}, \tag{3}\]

where \(x\in\mathbb{R}^{n}\) is the state of the system, \(u\in\mathbb{R}^{p}\) is the control input, and \(w\in\mathbb{R}^{q}\) is the disturbance. We assume that \(f:\mathbb{R}^{n}\times\mathbb{R}^{p}\times\mathbb{R}^{q}\rightarrow\mathbb{R}^ {n}\) is a parameterized vector field and the state feedback is parameterized by a \(k\)-layer feed-forward neural network controller \(N:\mathbb{R}^{n}\rightarrow\mathbb{R}^{p}\) defined by:

\[\xi^{(0)} =x\] \[\xi^{(i)} =\phi^{(i-1)}(W^{(i-1)}\xi^{(i-1)}+b^{(i-1)}),\quad i\in\{1, \ldots,k\}\] \[u =W^{(k)}\xi^{(k)}(x)+b^{(k)}:=N(x), \tag{4}\]

where \(n_{i}\) is the number of neurons in the \(i\)th layer, \(W^{(i-1)}\in\mathbb{R}^{n_{i}\times n_{i-1}}\) is the weight matrix of the \(i\)th layer, \(b^{(i-1)}\in\mathbb{R}^{n_{i}}\) is the bias vector of the \(i\)th layer, \(\xi^{(i)}(y)\in\mathbb{R}^{n_{i}}\) is the \(i\)th layer hidden variable, and \(\phi^{(i-1)}:\mathbb{R}^{n_{i}}\rightarrow\mathbb{R}^{n_{i}}\) is the \(i\)th layer diagonal activation function satisfying \(0\leq\frac{\phi^{(i-1)}_{\mathrm{}}(x)-\phi^{(i-1)}_{\mathrm{}}(y)}{x-y}\leq 1\) for every \(j\in\{1,\ldots,n_{i}\}\). One can show that a large class of activation functions including ReLU, leaky ReLU, sigmoid, and tanh satisfies this condition (after a possible re-scaling of their co-domain). With this neural network feedback controller, the closed-loop system is given by:

\[\dot{x}=f(x(t),N(x(t)),w):=f^{c}(x(t),w). \tag{5}\]

We assume that \(\mathcal{X}_{0}\subseteq\mathbb{R}^{n}\) is the initial set of states for the closed-loop system (5). This set can represent the uncertainties in the starting state of the system or localization errors in estimating the current state of the system. Moreover, we assume that the disturbance \(w\) belongs to the set \(\mathcal{W}\subseteq\mathbb{R}^{n}\) representing the model uncertainty or exogenous disturbances on the system. In this paper, we focus on the target-avoid problem as defined below:

**Target-Avoid Problem.** Given a target set \(\mathcal{G}\subset\mathbb{R}^{n}\), and an unsafe set \(\mathcal{S}_{\mathrm{unsafe}}\subset\mathbb{R}^{n}\) and a time interval \([0,T]\), check if the closed-loop system (5) reaches the target \(\mathcal{G}\) at time \(T\) while avoiding the unsafe region \(\mathcal{S}_{\mathrm{unsafe}}\). Mathematically,

1. \(\mathcal{R}_{f^{c}}(T,\mathcal{X}_{0},\mathcal{W})\subseteq\mathcal{G}\),
2. \(\mathcal{R}_{f^{c}}(t,\mathcal{X}_{0},\mathcal{W})\bigcap\mathcal{S}_{\mathrm{ unsafe}}=\emptyset\), for all \(t\in[0,T]\).

The target-avoid problem is a classical and prevalent notion of safety in many engineering applications, especially those involving safety-critical systems [9]. Moreover, diverse objectives including multiagent coordination, complex planning specified with formal languages and logics, and classical control-theoretic criteria such as robust invariance and stability can be achieved by concatenating and combining different instantiations of the target-avoid problem.

In general, computing the exact reachable sets of the closed-loop system (5) is not computationally tractable. Our approach to solve the target-avoid safety problem is based on constructing a computationally efficient over-approximation \(\overline{\mathcal{R}}_{f^{c}}(t,\mathcal{X}_{0},\mathcal{W})\) of the reachable set of the closed-loop system. Then, reaching the target is guaranteed by \(\overline{\mathcal{R}}_{f^{c}}(T,\mathcal{X}_{0},\mathcal{W})\subseteq\mathcal{G}\) and avoiding the unsafe region is guaranteed when \(\overline{\mathcal{R}}_{f^{c}}(t,\mathcal{X}_{0},\mathcal{W})\bigcap\mathcal{S}_{ \mathrm{unsafe}}=\emptyset\), for every \(t\in[0,T]\).

## IV Interval Analysis and Inclusion Functions

In this section, we develop a theoretical framework for interval reachability analysis of arbitrary mappings. The key element in our framework is the notion of inclusion function from interval arithmetic [10, 40].

### _Inclusion Functions_

Given a nonlinear input-output map \(g:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}\) and an interval input, the inclusion function of \(g\) provides interval bounds on the output of \(g\).

**Definition 1** (Inclusion function [10]).: Given a map \(g:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}\), the function \(\mathsf{G}=\left[\frac{\mathsf{G}}{\mathsf{G}}\right]:\mathcal{T}_{\geq 0}^{2n}\rightarrow \mathcal{T}_{\geq 0}^{2m}\) is

1. an _inclusion function_ for \(g\) if, for every \(x\leq\widehat{x}\) and every \(z\in[x,\widehat{x}]\), we have \[\underline{\mathsf{G}}(x,\widehat{x})\leq g(z)\leq\overline{\mathsf{G}}(x, \widehat{x}).\] (6)
2. a \([y,\widehat{y}]\)_-localized inclusion function_ if (6) holds for every \([x,\widehat{x}]\subseteq[y,\widehat{y}]\) and every \(z\in[x,\widehat{x}]\).

Moreover, an inclusion function \(\mathsf{G}\) for \(g\) is

1. _monotone_ if, for every \([\frac{x}{2}]\leq_{\mathrm{SE}}\left[\frac{y}{y}\right]\), we have \(\mathsf{G}(x,\widehat{x})\leq_{\mathrm{SE}}\mathsf{G}(y,\widehat{y})\);
2. _thin_ if, for every \(x\in\mathbb{R}^{n}\), \(\underline{\mathsf{G}}(x,x)=\overline{\mathsf{G}}(x,x)=g(x)\);
3. _minimal_ if, for every \(x\leq\widehat{x}\), the interval \([\underline{\mathsf{G}}(x,\widehat{x}),\overline{\mathsf{G}}(x,\widehat{x})]\) is the smallest interval containing \(g([x,\widehat{x}])\).

Given a nonlinear map \(g:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}\), one can show that its minimal inclusion function \(\mathsf{G}^{\min}=\left[\frac{\mathsf{G}^{\min}}{\mathsf{G}^{\min}}\right]\) is given by:

\[\underline{\mathsf{G}}^{\min}_{i}(x,\widehat{x})=\inf_{z\in[x,\widehat{x}]}g_{i}( z),\quad\overline{\mathsf{G}}^{\min}_{i}(x,\widehat{x})=\sup_{z\in[x,\widehat{x}]}g_{i}(z), \tag{7}\]

for every \(i\in\{1,\ldots,n\}\). The next theorem shows that the minimal inclusion function provides sharper over-approximations compared to Lipschitz bounds.
**Theorem 1** (Minimal inclusion functions).: _Consider the function \(g:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}\) and the mapping \(\mathsf{G}^{\min}=\left[\frac{\mathsf{G}^{\min}}{\mathsf{G}^{\min}}\right]\) defined by (7). Then, for every \(x\leq\widehat{x}\), we have_

\[\|\mathsf{G}^{\min}(x,\widehat{x})-\overline{\mathsf{G}}^{\min}(x,\widehat{x} )\|_{\infty}\leq\mathrm{Lip}_{\infty}(g)\|x-\widehat{x}\|_{\infty}.\]

Proof.: Let \(i\in\{1,\ldots,k\}\) be such that \(\|\overline{\mathsf{G}}^{\min}(x,\widehat{x})-\underline{\mathsf{G}}^{\min} (x,\widehat{x})\|_{\infty}=\left|\overline{\mathsf{G}}^{\min}_{i}(x,\widehat{ x})-\underline{\mathsf{G}}^{\min}_{i}(x,\widehat{x})\right|\). Note that since \(g\) is continuous and the box \([x,\widehat{x}]\) is compact, there exist \(\eta^{*},\xi^{*}\in[x,\widehat{x}]\) such that

\[\max_{y\in[x,\widehat{x}]}g_{i}(y)=g_{i}(\eta^{*}),\qquad\min_{y\in[x, \widehat{x}]}g_{i}(y)=g_{i}(\xi^{*}).\]

This implies that \(\|\overline{\mathsf{G}}^{\min}(x,\widehat{x})-\underline{\mathsf{G}}^{\min} (x,\widehat{x})\|_{\infty}=|g_{i}(\eta^{*})-g_{i}(\xi^{*})|\leq\|g(\xi^{*})-g (\eta^{*})\|_{\infty}\leq\mathrm{Lip}_{\infty}(g)\|\xi^{*}-\eta^{*}\|_{\infty} \leq\mathrm{Lip}_{\infty}(g)\|x-\widehat{x}\|_{\infty}\). 

### _Designing Inclusion Functions_

Given a function \(g\), finding a suitable inclusion function for \(g\) is a central problem in interval analysis. We review three approaches for constructing inclusion functions.

Natural inclusion functionAs we showed in the previous section, the minimal inclusion function can be computed using the optimization problem (7). However, for general functions, solving this optimization problem is not tractable. One feasible approach is to find the minimal inclusion function for a handful of elementary functions for which the optimization problem (7) is solvable. Then, for an arbitrary function, a natural inclusion function can be obtained by expressing it as a finite composition of operators and elementary functions [10, Theorem 2][41, Theorem 2]. Several software packages exist for automatically computing natural inclusions functions. In [41], we introduce the software package npinterval1 that implements intervals as a native data type within the Python toolbox numpy[42], yielding flexibility, efficiency through vectorization across arrays, and canonical constructions of natural inclusion functions.

Footnote 1: The code for npinterval is available at [https://github.com/gftastlab/npinterval](https://github.com/gftastlab/npinterval)

Jacobian-based inclusion functionsGiven a continuously differentiable function \(g:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}\), one can use upper and lower bounds on the Jacobian of \(g\) to construct inclusion functions for \(g\). These inclusion functions are often derived from the Taylor expansion of the function \(g\) around a certain point. Most commonly, this point is taken as the midpoint of an interval, leading to, in particular, the centered inclusion function [10, SS2.4.3] and the mixed-centered inclusion function [10, SS2.4.4]. In the next theorem, we develop a Jacobian-based inclusion function obtained by linearization around corners of an interval. As demonstrated in the sequel, such cornered inclusion functions turn out to be particularly amenable to analysis of closed-loop control systems.

**Proposition 1** (Jacobian-based cornered inclusion function).: _Given a continuously differentiable function \(g:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}\) such that \(Dg(z)\in[\underline{J}_{[z,\widehat{x}]},\overline{J}_{[x,\widehat{x}]}]\), for every \(z\in[x,\widehat{x}]\). The function \(\mathsf{G}^{\mathrm{jac}}:\mathcal{T}_{\geq 0}^{2n}\rightarrow\mathcal{T}_{\geq 0}^{2n}\) defined by_

\[\mathsf{G}^{\mathrm{jac}}(x,\widehat{x})=\begin{bmatrix}-[\underline{J}_{[x, \widehat{x}]}]^{-}&-[\underline{J}_{[x,\widehat{x}]}]^{-}\\ -[\overline{J}_{[x,\widehat{x}]}]^{+}&\overline{[J}_{[x,\widehat{x}]}]^{+} \end{bmatrix}\begin{bmatrix}x\\ \widehat{x}\end{bmatrix}+\begin{bmatrix}g(x)\\ g(x)\end{bmatrix}\]

_is an inclusion function for \(g\)._

Proof.: Since \(g\) is continuously differentiable, for every \(z\in[x,\widehat{x}]\), we have \(g(z)=g(x)+\int_{0}^{1}Dg(tz+(1-t)x)(z-x)dt\). Thus, for every \(x\in[x,\widehat{x}]\),

\[g(z) \leq g(x)+\overline{J}_{[x,\widehat{x}]}(z-x)\] \[=g(x)+[\overline{J}_{[x,\widehat{x}]}]^{+}(z-x)+[\overline{J}_{[x, \widehat{x}]}]^{-}(z-x)\] \[\leq g(x)+[\overline{J}_{[x,\widehat{x}]}]^{+}(\widehat{x}-x)= \overline{\mathsf{G}}^{\mathrm{jac}}(x,\widehat{x})\]

where the first inequality holds because \(z-x\geq\mathbb{O}_{n}\), the second equality holds by \(\overline{J}_{[x,\widehat{x}]}=[\overline{J}_{[x,\widehat{x}]}]^{+}+[ \overline{J}_{[x,\widehat{x}]}]^{-}\), and the last inequality holds because \(z-x\leq\widehat{x}-x\) and \(z-x\geq\mathbb{O}_{n}\). Similarly, one can show that, for every \(z\in[x,\widehat{x}]\), \(g(z)\geq g(x)+[\underline{J}_{[x,\widehat{x}]}]^{-}(\widehat{x}-x)=\underline{ \mathsf{G}}^{\mathrm{jac}}(x,\widehat{x})\), completing the proof. 

**Remark 1**.: (The role of corner points) The proof of Proposition 1 is based on Taylor expansion of \(g\) around the corner point \(x\). The proposition extends straightforwardly to the Taylor expansion of \(g\) around other corner points, i.e., around any point \(y\) such that \(y_{i}\in\{x_{i},\widehat{x}_{i}\}\), resulting in different inclusion functions. Both the Jacobian-based cornered inclusion function (Proposition 1) and the Jacobian-based centered inclusion function [10, SS2.4] are obtained using the Taylor expansion of the original function. In general, they are not comparable (cf. Example 1 and Table I). However, for functions that are monotone in their entries, the Jacobian-based cornered inclusion function will lead to the minimal inclusion function while the Jacobian-based centered inclusion function is often non-minimal. Figure 2 provides a pictorial comparison between the minimal, the centered, and the cornered inclusion functions.

The accuracy of the Jacobian-based cornered inclusion function can be improved at a cost of a slightly more complicated formulation using mixed differentiation in the same way that

Fig. 1: **Left:** npinterval is used to generate interval approximations for a function \(g(x_{1},x_{2})=[(x_{1}+x_{2})^{2},4\sin((x_{1}-x_{2})/4)]^{\top}\) using two different natural inclusion functions. Blue: using the inclusion functions for elementary functions \(x+x^{2}\) and \(x\mapsto\sin(x)\) in [41, Appendix A, Table I], \(\mathrm{G}^{\mathrm{recr}}\) rewriting \(g(x_{1},x_{2})=[x_{2}^{2}+2x_{1}x_{2}+x_{2}^{2},4\sin(x_{1}/4)\cos(x_{2}/4)-4 \cos(x_{1}/4)\sin(x_{2}/4)]^{\top}\) and obtains a different natural inclusion function based on composition of elementary inclusion functions in [41, Appendix A, Table I]. The approximations are generated using the initial set \([-1,1]\times[-1,1]\) and 2000 uniformly sampled outputs are shown in red. **Right:** the same function is analyzed, with the same two natural inclusion functions, but the initial set is partitioned into 1024 uniform sections, and the union of the interval approximations are shown.
the Jacobian-based centered inclusion function is improved with the mixed-centered inclusion function [10, SS2.4.4]. The idea is to use a step-by-step differentiation of the function with respect to each variable.

**Proposition 2** (Mixed Jacobian-based cornered inclusion function).: _Given a continuously differentiable function \(g:\mathbb{R}^{n}\to\mathbb{R}^{m}\) such that \(Dg(z)\in[\underline{\underline{I}}_{[x,\widehat{x}]},\overline{\underline{J}} _{[x,\widehat{x}]}]\), for every \(z\in[x,\widehat{x}]\). The function \(\mathsf{G}^{\mathrm{lac-m}}:\mathcal{T}_{2\geq 0}^{2n}\to\mathcal{T}_{20}^{2n}\) defined by_

\[\mathsf{G}^{\mathrm{jac-m}}(x,\widehat{x})=\begin{bmatrix}-[\underline{M}]^{- }&[\underline{M}]^{-}\\ -[\overline{M}]^{+}&[\overline{M}]^{+}\end{bmatrix}\begin{bmatrix}x\\ \widehat{x}\end{bmatrix}+\begin{bmatrix}g(x)\\ g(x)\end{bmatrix}\]

_is an inclusion function for \(g\), where the \(i\)th columns of \(\underline{M}\) and \(\overline{M}\) are defined by_

\[\underline{M}_{i}=(\underline{J}_{[x,\widehat{x}[h_{i},i]}))_{i},\qquad \overline{M}_{i}=(\overline{J}_{[x,\widehat{x}[h_{i},i]}))_{i},\]

_with \(R_{i}=\{i+1,\ldots,n\}\), for every \(i\in\{1,\ldots,n\}\). Moreover, for every \(x\leq z\leq\widehat{x}\), we have_

\[\underline{\mathsf{G}}^{\mathrm{jac}}(x,\widehat{x})\leq\underline{\mathsf{G} }^{\mathrm{jac-m}}(x,\widehat{x}),\qquad\overline{\mathsf{G}}^{\mathrm{jac- m}}(x,\widehat{x})\leq\overline{\mathsf{G}}^{\mathrm{jac}}(x,\widehat{x}).\]

Proof.: Since \(g\) is continuously differentiable, for every \(z\in[x,\widehat{x}]\), we have \(g(z_{[R_{1};x]})=g(x)+\int_{0}^{1}D_{x_{1}}g(tz_{[R_{1};x]}+(1-t)x)(z_{1}-x_{ 1})dt\). Thus, for every \(x\in[x,\widehat{x}]\), \(g(z_{[R_{1};x]})\leq g(x)+\overline{J}_{[x,\widehat{x}[h_{1},i]}(z_{1}-x_{1}) \leq g(x)+[\overline{J}_{[x,\widehat{x}[h_{1},i]}]^{+}(\widehat{x}_{1}-x_{1}).\) where the first inequality holds because \(z_{1}-x_{1}\geq\mathbb{0}_{n}\) and the second inequality holds because \(z_{1}-x_{1}\leq\widehat{x}_{1}-x_{1}\) and \(z_{1}-x_{1}\geq\mathbb{0}_{n}\). By successively applying the above procedure on \(g(z_{[R_{1};x]})\), one can show that \(g(z)\leq\underline{\mathsf{G}}^{\mathrm{jac}}(x,\widehat{x})\). Similarly, one can show that \(g(z)\geq\overline{\mathsf{G}}^{\mathrm{jac}}(x,\widehat{x})\) and thus \(\mathsf{G}^{\mathrm{jac}}\) is an inclusion function for \(g\). Finally, one can show that \(\underline{J}_{[x,\widehat{x}]}\leq\underline{M}\) and \(\overline{M}\leq\overline{J}_{[x,\widehat{x}]}\). This implies that \(-[\underline{M}]^{-}x+[\underline{M}]^{-}\widehat{x}\leq-[\underline{J}_{[x, \widehat{x}]}^{-}x+[\overline{J}_{[x,\widehat{x}]}^{-}\widehat{x}.\) As a result \(\underline{\mathsf{G}}^{\mathrm{jac}}(x,\widehat{x})\leq\underline{\mathsf{G} }^{\mathrm{jac-m}}(x,\widehat{x})\). Similarly, one can show that \(\overline{\mathsf{G}}^{\mathrm{jac-m}}(x,\widehat{x})\leq\overline{\mathsf{G }}^{\mathrm{jac}}(x,\widehat{x})\). 

**Example 1** (Cornered inclusion functions).: Consider the function \(f:\mathbb{R}^{2}\to\mathbb{R}^{2}\) defined by \(f(x_{1},x_{2})=\begin{bmatrix}(x_{1}+x_{2})^{2}\\ x_{1}+x_{2}+2x_{1}x_{2}\end{bmatrix}\). We compute the natural inclusion function, the (mixed) Jacobian-based cornered inclusion function, and the (mixed) Jacobian-based centered inclusion function for \(f\) on the interval \([-0.1,0.1]\times[-0.1,0.1]\). The results are shown in Table I, for runtimes averaged over \(10\,000\) runs. We observe that the output interval estimated by the mixed Jacobian-based centered (cornered) is tighter than the interval obtained from the Jacobian-based centered (cornered) inclusion function. This observation is consistent with the comparison between mixed and non-mixed Jacobian-based cornered inclusion functions in Proposition 2. However, there is generally no ranking between other methods, as each can perform better in different scenarios.

Decomposition-based inclusion functionsAnother approach for constructing inclusion functions is by decomposing the function into the sum of increasing and decreasing parts [43, 15]. This separation is often realized using a decomposition function.

**Definition 2** (Decomposition function).: Given a function \(g:\mathbb{R}^{n}\to\mathbb{R}^{m}\), a _decomposition function_ for \(g\) is a map \(d:\mathcal{T}^{2n}\to\mathbb{R}^{m}\) that satisfies

1. \(g(x)=d(x,x)\), for every \(x\in\mathbb{R}^{n}\);
2. \(d(x,y)\leq d(\widehat{x},y)\), for every \(y\in\mathbb{R}^{n}\) and every \(x\leq\widehat{x}\);
3. \(d(x,\widehat{y})\leq d(x,y)\), for every \(x\in\mathbb{R}^{n}\) and every \(y\leq\widehat{y}\).

One can show that every continuous function has at least one decomposition function. Given a map \(g:\mathbb{R}^{n}\to\mathbb{R}^{m}\), we define the tight decomposition function for \(g\) as the mapping \(d^{\mathrm{tight}}:\mathcal{T}^{2n}\to\mathbb{R}^{2m}\) given by

\[d^{\mathrm{tight}}_{i}(x,\widehat{x})=\begin{cases}\min_{z\in[x,\widehat{x}]}g _{i}(z)&x\leq\widehat{x}\\ \max_{z\in[x,\widehat{x}]}g_{i}(z)&\widehat{x}\leq x\end{cases} \tag{8}\]

for every \(i\in\{1,\ldots,n\}\). In general, the decomposition function for function \(g\) is not unique. The next theorem, whose proof is straightforward, shows how to construct inclusion functions from decomposition functions.

**Proposition 3** (Decomposition-based inclusion functions).: _Consider a map \(g:\mathbb{R}^{n}\to\mathbb{R}^{m}\) with a decomposition function \(d:\mathcal{T}^{2n}\to\mathbb{R}^{m}\). Then, the following statements hold:_

1. _the map_ \(\mathsf{G}^{d}:\mathcal{T}_{\geq 0}^{2n}\to\mathbb{R}^{2m}\) _defined by_ \[\mathsf{G}^{d}(x,\widehat{x})=\begin{bmatrix}d(x,\widehat{x})\\ d(\widehat{x},x)\end{bmatrix},\quad\text{ for all }x\leq\widehat{x}\] _is a thin monotone inclusion function for_ \(g\)_._
2. _the mapping_ \(d^{\mathrm{tight}}\) _defined in (_8_) is a decomposition function for_ \(g\) _and_ \[\mathsf{G}^{\min}(x,\widehat{x})=\begin{bmatrix}d^{\mathrm{tight}}(x,\widehat{ x})\\ d^{\mathrm{tight}}(\widehat{x},x)\end{bmatrix},\quad\text{ for all }x\leq\widehat{x},\] _where_ \(\mathsf{G}^{\min}\) _is the minimal inclusion function for_ \(g\)

Fig. 2: Pictorial comparison between the minimal (equation (7)), the Jacobian-based centered ([10, §2.4.3], and the Jacobian-based cornered (Proposition 1) inclusion functions for \(f(x)=x^{2}\) and \(f(x)=x^{3}\) on the interval \([-1,1]\). Note that, for the monotone function \(f(x)=x^{3}\), the intersection of the two Jacobian-based cornered inclusion functions (shown by blue and red dashed lines) leads to the minimal inclusion function.
**Remark 2** (Comparison with the literature).: For a vector field \(f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}\), the decomposition function has already been studied in mixed monotone theory (cf. [44, 16, Definition 2], [45, Definition 4], and [46]). Definition 2 extends this classical notion to arbitrary functions. Theorem 3(ii) shows that the minimal inclusion function, defined by equation (7), is a decomposition-based inclusion function. Moreover, the tight decomposition function \(d^{\rm tight}\) in equation (8) coincides with the tight decomposition function in [47, Definition 2].

One can combine two inclusion functions to obtain another inclusion function with tighter estimates for the output of the original function.

**Proposition 4** (Intersection of inclusion functions).: _Given a map \(g:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}\) with two inclusion functions \(\mathsf{G}_{1}=\left[\frac{\mathsf{G}_{1}}{\mathsf{G}_{1}}\right]:\mathcal{T} _{\geq 0}^{2n}\rightarrow\mathcal{T}_{\geq 0}^{2m}\) and \(\mathsf{G}_{2}=\left[\frac{\mathsf{G}_{2}}{\mathsf{G}_{2}}\right]:\mathcal{T} _{\geq 0}^{2n}\rightarrow\mathcal{T}_{\geq 0}^{2m}\). Define the mapping \(\mathsf{G}_{1}\wedge\mathsf{G}_{2}=\left[\frac{\mathsf{G}_{1}\wedge\mathsf{G }_{2}}{\mathsf{G}_{1}\wedge\mathsf{G}_{2}}\right]:\mathcal{T}_{\geq 0}^{2n} \rightarrow\mathcal{T}_{\geq 0}^{2m}\) by_

\[[\underline{\mathsf{G}}_{1}\wedge\mathsf{G}_{2}](x,\widehat{x}) =\max\{\mathsf{G}_{1}(x,\widehat{x}),\mathsf{G}_{2}(x,\widehat{x})\},\] \[[\overline{\mathsf{G}}_{1}\wedge\overline{\mathsf{G}}_{2}](x, \widehat{x}) =\min\{\overline{\mathsf{G}}_{1}(x,\widehat{x}),\overline{\mathsf{G }}_{2}(x,\widehat{x})\}.\]

_Then \(\mathsf{G}_{1}\wedge\mathsf{G}_{2}\) is an inclusion function for \(g\) and it provides over-approximations of the output of \(g\) which are tighter than both \(\mathsf{G}_{1}\) and \(\mathsf{G}_{2}\)._

Proof.: Note that both \(\mathsf{G}_{1}\) and \(\mathsf{G}_{2}\) are inclusion functions for \(g\). Thus, for every \(x\leq z\leq\widehat{x}\), we have \(\mathsf{G}_{1}(x,\widehat{x})\leq g(z)\) and \(\mathsf{G}_{2}(x,\widehat{x})\leq g(z)\). This implies that \(\max\{\underline{\mathsf{G}}_{1}(x,\widehat{x}),\underline{\mathsf{G}}_{2}(x,\widehat{x})\}=|\underline{\mathsf{G}}_{1}\wedge\underline{\mathsf{G}}_{2} |(x,\widehat{x})\leq g(z)\). Similarly, one can show that \(g(z)\leq[\overline{\mathsf{G}}_{1}\wedge\overline{\mathsf{G}}_{2}](x, \widehat{x})\). 

### _Convergence of inclusion functions_

We now define the convergence rate of an inclusion function which can be used to capture the accuracy of the interval estimations.

**Definition 3** (Convergence rate of inclusion functions).: Consider the mapping \(g:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}\) with an inclusion function \(\mathsf{G}=\left[\frac{\mathsf{G}}{\mathsf{G}}\right]:\mathcal{T}_{\geq 0}^{2n} \rightarrow\mathcal{T}_{\geq 0}^{2m}\). The rate of convergence of \(\mathsf{G}\) is the largest \(\alpha\in\mathbb{R}_{\geq 0}\cup\{\infty\}\) such that

\[\|\overline{\mathsf{G}}(x,\widehat{x})-\underline{\mathsf{G}}(x,\widehat{x})\|_{\infty}-\max_{z,y\in[x,\widehat{x}]}\|g(z)-g(y)\|_{\infty}\] \[\leq\beta\|\widehat{x}-x\|_{\infty}^{\alpha} \tag{9}\]

for some \(\beta\in\mathbb{R}_{\geq 0}\) and for every \(x\leq\widehat{x}\).

Roughly speaking, the rate of convergence of an inclusion function \(\mathsf{G}\) shows how accurately it estimates the range of the original function \(g\). Indeed, for the minimal inclusion function, the interval \([\underline{\mathsf{G}}(x,\widehat{x}),\overline{\mathsf{G}}(x,\widehat{x})]\) is the tightest interval containing \(g([x,\widehat{x}]]\). Thus the left hand side of inequality (9) is always zero and the rate of convergence of the minimal inclusion function is \(\alpha=\infty\). It is shown that, in general, the convergence rate of the natural inclusion functions is \(\alpha=1\)[40, Lemma 4.1] and the convergence rate of the Jacobian-based inclusion functions are \(\alpha=2\)[10, SS2.4.5].

### _Inclusion functions for Neural Networks_

Given a neural network \(N\) of the form (4) and any interval \([y,\widehat{y}]\subseteq\mathbb{R}^{n}\), as described in the sequel, our framework proposes using existing neural network verification algorithms to construct a \([y,\widehat{y}]\)-localized inclusion function of the form \(\mathsf{N}_{[y,\widehat{y}]}=\left[\frac{\mathsf{N}_{[y,\widehat{y}]}}{ \mathsf{N}_{[y,\widehat{y}]}}\right]:\mathcal{T}_{\geq 0}^{2n}\rightarrow\mathcal{T}_{ \geq 0}^{2p}\) for \(N\) satisfying

\[\underline{\mathsf{N}}_{[y,\widehat{y}]}(x,\widehat{x})\leq N(x)\leq\overline {\mathsf{N}}_{[y,\widehat{y}]}(x,\widehat{x}), \tag{10}\]

for any \(x\in[x,\widehat{x}]\subseteq[y,\widehat{y}]\). A large number of the existing neural network verification algorithms can provide bounds of the form (10) for the output of the neural networks, including CROWN [25], LipSDP [24], and IBP [23]. In particular, some neural network verification algorithms can provide _affine_\([y,\widehat{y}]\)-localized inclusion functions for \(N\). Examples of these neural network verification algorithms include CROWN and its variants [25]. Given an interval \([y,\widehat{y}]\), these algorithms provide a tuple \((\underline{C},\overline{C},\underline{d},\overline{d})\) defining affine upper and lower bounds for the output of the neural network

\[\underline{C}(y,\widehat{y})x+\underline{d}(y,\widehat{y})\leq N(x)\leq \overline{C}(y,\widehat{y})x+\overline{d}(y,\widehat{y}), \tag{11}\]

for every \(x\in[y,\widehat{y}]\). Using these linear bounds, we can construct the affine \([y,\widehat{y}]\)-localized inclusion function for \(N\):

\[\underline{\mathsf{N}}_{[y,\widehat{y}]}(x,\widehat{x})=[\underline{C}(y, \widehat{y})]^{+}x+[\underline{C}(y,\widehat{y})]^{-}\widehat{x}+\underline{d}(y,\widehat{y}), \tag{12}\] \[\overline{\mathsf{N}}_{[y,\widehat{y}]}(x,\widehat{x})=[\overline{C} (y,\widehat{y})]^{+}\widehat{x}+[\overline{C}(y,\widehat{y})]^{-}x+\overline{d} (y,\widehat{y})\]

for any \([x,\widehat{x}]\subseteq[y,\widehat{y}]\).

**Remark 3**.:
1. _(Computational complexity):_ the computational complexity of finding the \([y,\widehat{y}]\)-localized inclusion function \(\mathsf{N}_{[y,\widehat{y}]}=\left[\frac{\mathsf{N}_{[y,\widehat{y}]}}{ \mathsf{N}_{[y,\widehat{y}]}}\right]\) depends on the neural network verification algorithm. For instance, for a neural network with \(k\)-layer and \(N\) neurons per layer, the computational complexity of CROWN [25] in finding the bounds of the form (10) is \(\mathcal{O}(k^{2}N^{3})\).
2. _(Inclusion function):_ Since the \([y,\widehat{y}]\)-localized inclusion function \(\mathsf{N}_{[y,\widehat{y}]}=\left[\frac{\mathsf{N}_{[y,\widehat{y}]}}{ \mathsf{N}_{[y,\widehat{y}]}}\right]\) can be obtained for every \([y,\widehat{y}]\), we can define an inclusion function \(\mathsf{N}=\left[\frac{\mathsf{N}}{\mathsf{N}}\right]\) for the neural network \(N\) by \(\mathsf{N}(x,\widehat{x})=\mathsf{N}_{[x,\widehat{x}]}(x,\widehat{x})\), for all \(x\leq\widehat{x}\).

## V Interval Reachability of Dynamical Systems

In this section, we study reachability of the dynamical system (1) using interval analysis. Our key idea is to embed the original dynamical system (1) into an _embedding system_ on \(\mathbb{R}^{2n}\) and use trajectories of the embedding system to study the propagation of the state bounds with time. Our main goal is to use the notion of inclusion function (cf. Definition 1) to design embedding systems for dynamical systems. Consider the dynamical system (1) with an inclusion function \(\mathsf{F}=\left[\frac{\mathsf{F}}{\mathsf{F}}\right]:\mathcal{T}_{\geq 0}^{2n} \times\mathcal{T}_{\geq 0}^{2q}\rightarrow\mathcal{T}_{\geq 0}^{2n}\) for \(f\). We define the associated _embedding system_ on \(\mathbb{R}^{2n}\):

\[\dot{x}_{i} =\mathsf{E}_{i}(x,\widehat{x}_{[i:x]},w,\widehat{w}),\] \[\dot{\widehat{x}}_{i} =\overline{\mathsf{F}}_{i}(x_{[i:\widehat{x}]},\widehat{x},w, \widehat{w}), \tag{13}\]
for every \(i\in\{1,\ldots,n\}\). For instance, using the minimal inclusion function \(\mathsf{F}^{\min}=\left[\frac{\mathsf{F}^{\min}}{\mathsf{F}^{\min}}\right]\) for \(f\), the associated embedding system is given by:

\[\dot{x}_{i} =\min_{\begin{subarray}{c}z\in[x,\overline{x}],\ z_{i}=x_{i},\\ u\in[w,\overline{w}]\end{subarray}}f_{i}(z,u),\] \[\dot{\bar{x}}_{i} =\max_{\begin{subarray}{c}z\in[x,\overline{z}],\ z_{i}=\bar{z} _{i},\\ u\in[w,\overline{w}]\end{subarray}}f_{i}(z,u), \tag{14}\]

for every \(i\in\{1,\ldots,n\}\). Given the disturbance set \(\mathcal{W}=[\underline{w},\overline{w}]\) and the initial set \(\mathcal{X}_{0}=[\underline{x}_{0},\overline{x}_{0}]\), one can use a single trajectory of the embedding system (13) to obtain over-approximations of reachable sets of the system (1). Moreover, it can be shown that among all the inclusion functions used for constructing the embedding system (13), the embedding system obtained from the minimal inclusion function (i.e., the dynamical system (14)) provides the tightest reachable set over-approximations at any time \(t\in\mathbb{R}_{\geq 0}\).

**Proposition 5** (Reachability using embedding systems).: _Consider the dynamical system (1) with the disturbance set \([\underline{w},\overline{w}]\) and the initial set \(\mathcal{X}_{0}=[\underline{x}_{0},\overline{x}_{0}]\). Suppose that \(\mathsf{F}=\left[\frac{\mathsf{F}}{\mathsf{F}}\right]\) is an inclusion function for \(f\), and \(t\mapsto\left[\frac{x}{\overline{x}(t)}\right]\) and \(t\mapsto\left[\frac{x}{\overline{x}^{\min}(t)}\right]\) are the trajectories of embedding systems (13) and (14), starting from \(\left[\frac{x_{0}}{\overline{x}_{0}}\right]\) with fixed disturbance \(\left[\frac{w}{\overline{w}}\right]=[\frac{w}{\overline{w}}]\). Then, for every \(t\in\mathbb{R}_{\geq 0}\), we have_

\[\mathcal{R}_{f}(t,[\underline{x}_{0},\overline{x}_{0}],[\underline{w}, \overline{w}])\subseteq[\underline{x}^{\min}(t),\overline{x}^{\min}(t)] \subseteq[\underline{x}(t),\overline{x}(t)].\]

Proof.: We first show that the dynamical system (14) is a monotone dynamical system with respect to the southeast order \(\leq_{\mathrm{SE}}\). Let \(\left[\frac{x}{\overline{z}}\right]\leq_{\mathrm{SE}}\left[\frac{y}{\overline {y}}\right]\) be such that \(x_{i}=y_{i}\) and let \(\left[\frac{w}{\overline{w}}\right]\leq_{\mathrm{SE}}\left[\frac{v}{\overline {v}}\right]\). Then, we have \([y,\widehat{y}]\subseteq[x,\widehat{x}]\) and \([v,\widehat{v}]\subseteq[w,\widehat{w}]\), and we can compute

\[\mathsf{E}^{\min}_{i}(x,\widehat{x},w,\widehat{w}) =\min_{\begin{subarray}{c}z\in[x,\overline{z}],\ z_{i}=x_{i},\\ u\in[w,\overline{w}]\end{subarray}}f_{i}(z,u)\] \[\leq\min_{\begin{subarray}{c}z\in[y,\widehat{y}],\ z_{i}=y_{i},\\ u\in[v,\overline{v}]\end{subarray}}f_{i}(z,u)=\mathsf{E}^{\min}_{i}(y, \widehat{y},v,\widehat{v}),\]

for every \(i\in\{1,\ldots,n\}\). Similarly, let \(\left[\frac{x}{\overline{z}}\right]\leq_{\mathrm{SE}}\left[\frac{y}{\overline {y}}\right]\) be such that \(\widehat{x}_{i}=\widehat{y}_{i}\) and let \(\left[\frac{w}{\overline{w}}\right]\leq_{\mathrm{SE}}\left[\frac{v}{\overline {v}}\right]\). Then, we have \([y,\widehat{y}]\subseteq[x,\widehat{x}]\) and \([v,\widehat{v}]\subseteq[w,\widehat{w}]\), and similar reasoning implies \(\overline{\mathsf{F}}^{\min}_{i}(x,\widehat{x},w,\widehat{w})\leq\overline{ \mathsf{F}}^{\min}_{i}(y,\widehat{y},v,\widehat{v})\) for every \(i\in\{1,\ldots,n\}\). As a result, the embedding system (14) is continuous-time monotone with respect to \(\leq_{\mathrm{SE}}\). Let \(x_{0}\in[\underline{x}_{0},\overline{x}_{0}]\) and \(u\in[\underline{w},\overline{w}]\) and \(t\mapsto x_{u}(t)\) be the trajectory of the system (1) starting from \(x_{0}\) with disturbance \(u\). Note that \(t\mapsto\left[\frac{x_{u}(t)}{x_{u}(t)}\right]\) is a solution of the embedding system (14) with the initial condition \(\left[\frac{x_{0}}{x_{0}}\right]\) and the disturbance \(\left[\frac{u}{u}\right]\). Moreover, we know that \(\left[\frac{x_{0}}{x_{0}}\right]\leq_{\mathrm{SE}}\left[\frac{x_{0}}{x_{0}}\right]\) and \(\left[\frac{w}{\overline{w}}\right]\leq_{\mathrm{SE}}\left[\frac{u}{u}\right]\). Therefore, by monotonicity of the embedding system (14) with respect to \(\leq_{\mathrm{SE}}\), we get \(\left[\frac{x^{\min}(t)}{x^{\min}(t)}\right]\leq_{\mathrm{SE}}\left[\frac{x_{ u}(t)}{x_{u}(t)}\right]\), for every \(t\in\mathbb{R}_{\geq 0}\). This implies that \(\mathcal{R}_{f}(t,[\underline{x}_{0},\overline{x}_{0}],[\underline{w}, \overline{w}])\subseteq[\underline{x}^{\min}(t),\overline{x}^{\min}(t)]\), for every \(t\in\mathbb{R}_{\geq 0}\). On the other hand, \(\mathsf{F}^{\min}\) is the minimal inclusion function for \(f\) and therefore, for every \(x\leq\widehat{x}\) and every \(w\leq\widehat{w}\), we have \(\left[\frac{\mathsf{F}(x,\widehat{x},w,\widehat{w})}{\mathsf{F}(x,\widehat{x}, w,\widehat{w})}\right]\leq_{\mathrm{SE}}\left[\frac{\mathsf{F}^{\min}(x, \widehat{x},w,\widehat{w})}{\mathsf{F}^{\min}(x,\widehat{x},w,\widehat{w})}\right]\). Thus, for every \(i\in\{1,\ldots,n\}\), every \(x\leq\widehat{x}\), and every \(w\leq\widehat{w}\),

\[\mathsf{E}_{i}(x,\widehat{x}_{[i:x]},w,\widehat{w})\leq\min_{ \begin{subarray}{c}z\in[\underline{x}],\ z_{i}=x_{i},\\ \xi\in[w,\widehat{w}]\end{subarray}}f_{i}(z,\xi),\] \[\mathsf{\overline{F}}_{i}(x_{[i:\widehat{x}]},\widehat{x},w, \widehat{w})\geq\max_{\begin{subarray}{c}z\in[x,\overline{z}],\ z_{i}=\bar{z} _{i},\\ \xi\in[w,\widehat{w}]\end{subarray}}f_{i}(z,\xi).\]

Now, we can use the classical monotone comparison Lemma [48, Theorem 3.8.1] to obtain \(\left[\frac{x}{\overline{x}(t)}\right]\leq_{\mathrm{SE}}\left[\frac{\mathsf{F} ^{\min}_{i}(t)}{x^{\min}(t)}\right]\), for every \(t\in\mathbb{R}_{\geq 0}\). This implies that \([\underline{x}^{\min}(t),\overline{x}^{\min}(t)]\subseteq[\underline{x}(t), \overline{x}(t)]\), for every \(t\in\mathbb{R}_{\geq 0}\), and completes the proof. 

**Remark 4** (Comparison with the literature).:
1. It is straightforward to extend Proposition 5 to the case when the mapping \(\mathsf{F}\) is a \([y,\widehat{y}]\times\mathbb{R}^{q}\)-localized inclusion function for \(f\). In this setting, the results of Proposition 5 hold as long as we have \([\underline{x}(t),\overline{x}(t)]\subseteq[y,\widehat{y}]\).
2. Proposition 5 can be considered as a generalization of [49, Proposition 6] and [18, Proposition 1]. Indeed, in the special case that \(\mathsf{F}\) is a decomposition-based inclusion function constructed from the decomposition function \(d\), one can recover the embedding system \[\dot{x}_{i} =d(x,\widehat{x}_{[i:x]},w,\widehat{w}),\] \[\dot{\bar{x}}_{i} =d(\widehat{x},x_{[i:\widehat{x}]},\widehat{w},w)\] for every \(i\in\{1,\ldots,n\}\). This embedding system is identical to [49, Equation (7)] and [18, Equation (3)]. We highlight that this extension is crucial for our framework as the inclusion functions that we obtain for learning-based systems are neither thin nor decomposition-based.
3. Proposition 5 can alternatively be proved using [12, Theorem 2]. However, our proof of Proposition 5 is different in that it uses monotone system theory [50, 51] and classical monotone comparison Lemma [48, Theorem 3.8.1]. Indeed, compared to [12, Theorem 2], our proof techniques allow to compare accuracy of different over-approximations of reachable sets.

We now study the accuracy of the over-approximations provided in Proposition 5. We first introduce the convergence rate of an embedding system for approximating the reachable set of the dynamical system (1).

**Definition 4** (Convergence rate of embedding systems).: Consider the dynamical system (1) with the embedding system (13). Then the convergence rate of the embedding system (13) is the largest value of \(\alpha\in\mathbb{R
system (13) and the convergence rate of its associated inclusion function \(\mathsf{F}\) (as in Definition 3)? In the special case when \(f\) is a discrete-time monotone vector field, one can choose the inclusion function \(\mathsf{F}(\underline{x},\overline{x})=\left[\begin{smallmatrix}f(\underline{x} )\\ f(\underline{x})\end{smallmatrix}\right]\) for \(f\). In this case, it can be shown that the convergence rate of both the embedding system and the inclusion function are the same and equal to \(\alpha=\infty\). Given a Lipschitz and monotone inclusion function \(\mathsf{F}\) for \(f\), one can show that the convergence rate of the embedding system (13) satisfies \(\alpha\geq 1\)[40, Theorem 4.1]. However, in general, the inclusion function \(\mathsf{F}\) does not reveal more information about the convergence rate of the embedding system (13). Even when the convergent rate of the inclusion function \(\mathsf{F}\) is \(\alpha=\infty\), i.e., the embedding system is given by equations (14), one might get a convergence rate of \(\alpha=1\) for the embedding system.

## VI Interval reachability of learning-based systems

In this section, we develop a computationally efficient framework for studying reachability of neural network controlled systems. The key idea of our approach is to employ the framework of Section V and to capture the interactions between the open-loop system and the neural network controller using a suitable closed-loop embedding system. However, finding a closed-loop embedding system (or, equivalently, computing an inclusion function for the closed-loop vector field) using the approaches in Section IV is often not tractable, as the neural networks can have a large number of parameters. Instead, we propose to use a compositional approach based on an inclusion function for the open-loop dynamics and an inclusion function for the neural network controller. We develop two different approaches for designing the closed-loop inclusion function from the open-loop vector field \(f\) and the localized inclusion function of the neural network.

First, we propose an _interconnection-based_ approach where the key idea is to consider the closed-loop embedding system as the feedback interconnection of the open-loop embedding system and the neural network inclusion function. This method, as elaborated below, can capture some of the stabilizing effects of the neural network controller by approximating the input-output behavior of the neural network separately on the faces of a hyperrectangle. Thus, we consider this approach a semi-input-output approach. The advantages of the interconnection-based approach are its computational speed and that it is amenable to any open-loop inclusion function.

Second, we propose an _interaction-based_ approach where the key idea is to use the Jacobian-based cornered inclusion function (developed in Proposition 1) to more fully capture the interaction between the neural network controller and the system. The advantages of this approach are reduced conservatism compared to the interconnection-based inclusion function while retaining much of the computational efficiencies of interval reachability analysis.

Our methods are applicable to general nonlinear systems with neural network controllers. Nonetheless, the essence of our approaches, and their advantages over a naive input-output approach, are evident even in the simple setting of a linear system with a linear feedback controller, as shown in the following illustrative example.

**Example 2** (Invariant intervals).: Consider the control system \(\dot{x}=Ax+Bu\) where \(x\in\mathbb{R}^{2}\) and \(A=\left[\begin{smallmatrix}-2&1\\ 1&-2\end{smallmatrix}\right]\) and \(B=\left[\begin{smallmatrix}0\\ 1\end{smallmatrix}\right]\). Note that \(A\) is Hurwitz so the origin is the globally asymptotically stable equilibrium point of the open-loop system when \(u\equiv 0\). Consider a state feedback controller \(u=\pi(x)=Kx=k_{1}x_{1}+k_{2}x_{2}\) with \(k_{1}=k_{2}=-3\) designed to, _e.g._, achieve faster convergence to the origin. As a special case of reachability analysis, our goal is to characterize forward invariant intervals of the closed-loop system around the origin. We compare three different approaches for finding invariant intervals.

_Naive input-output approach:_ this approach decouples the system from the controller and only uses knowledge of the output range of the controller. In particular, this approach seeks an interval \(S=[-\xi,\xi]\), \(\xi\in\mathbb{R}^{2}_{\geq 0}\) such that \(S\) is robustly forward invariant for the open loop system under any \(u\in\pi(S)\). However, such a robustly forward invariant set does not exist. By contradiction, suppose \(S=[-\xi,\xi]\) is forward invariant and consider the point \(x=\left[\begin{smallmatrix}0\\ \xi_{2}\end{smallmatrix}\right]\) on the boundary of this set and pick \(u=\pi(-x)=3\xi_{2}\in\pi(S)\), for which \(\dot{x}_{2}=-2\xi_{2}+u=\xi_{2}>0\) and the resulting trajectory immediately leaves \(S\). Despite the stability of the open-loop system and the stabilizing effect of the controller, this approach fails to capture the stability of the closed-loop system.

_Interconnection-based approach:_ this approach also decouples the system from the controller and only assumes knowledge of the output of the controller to characterize forward invariant intervals \(S=[-\xi,\xi]\), for some \(\xi\in\mathbb{R}^{2}_{\geq 0}\). However, it considers the interconnection between the system and the controller on the edges of the interval \(S\). Let \(S^{-}_{1},S^{+}_{1},S^{-}_{2},S^{+}_{2}\) denote the edges of \(S\) defined by \(S^{\pm}_{i}=\{x\in S\mid x_{i}=\pm\xi_{i}\}\) for \(i\in\{1,2\}\). Then a sufficient condition for forward invariance of \(S\) is that \(\dot{x}_{i}\geq 0\) for every \(x\in S^{-}_{i}\), and \(\dot{x}_{i}\leq 0\) for every \(x\in S^{+}_{i}\). Note that, if \(\xi_{1}/\xi_{2}\geq\frac{1}{2}\), then we have \(\dot{x}_{1}\leq 0\) on \(S^{+}_{1}\) and we have \(\dot{x}_{1}\geq 0\) on \(S^{-}_{1}\). On \(S^{+}_{2}\), we have \(\dot{x}_{2}\in[-\xi_{1},\xi_{1}]-2\xi_{2}+u\) and \(u\in-3[-\xi_{1},\xi_{1}]-3\xi_{2}\). This implies that \(\dot{x}_{2}\in[4-\xi_{1},\xi_{1}]-5\xi_{2}\), and, therefore, we need \(\xi_{1}/\xi_{2}\leq\frac{5}{4}\) for \(\dot{x}_{2}\leq 0\). A similar condition is required to ensure \(\dot{x}_{2}\geq 0\) on \(S^{-}_{2}\). This implies that, using this interconnection-based approach, we certify \(S=[-\xi,\xi]\) is forward invariant for any \(\xi\in\mathbb{R}^{2}_{>0}\) satisfying \(\frac{1}{2}\leq\xi_{1}/\xi_{2}\leq\frac{5}{4}\).

_Interaction-based approach:_ this approach uses the knowledge of the controller functional dependency \(\pi(x)=Kx\) and the open-loop dynamics to obtain the closed-loop system \(\dot{x}=(A+BK)x=\left[\begin{smallmatrix}-2&-2\\ -2&-5\end{smallmatrix}\right]x\). One can show that the interval \(S=[-\xi,\xi]\) is forward invariant for this closed-loop system for any \(\xi\in\mathbb{R}^{2}_{\geq 0}\) satisfying \(\frac{1}{2}\leq\xi_{1}/\xi_{2}\leq\frac{5}{2}\).

Now, we develop a mathematical framework for reachability analysis of the learning-based closed-loop system (5). Given an open-loop system of the form (1), we use approaches developed in Section IV-B to construct an inclusion function \(\mathsf{F}^{\rm o}=\left[\begin{smallmatrix}\overline{\mathsf{F}}^{\rm o}\\ \overline{\mathsf{F}}^{\rm o}\end{smallmatrix}\right]:\mathcal{T}^{2n}_{\geq 0} \times\mathcal{T}^{2p}_{\geq 0}\times\mathcal{T}^{2q}_{\geq 0}\to\mathcal{T}^{2n}_{\geq 0}\) for the open-loop vector field \(f\). We also assume that we have access to a neural network inclusion function \(\mathsf{N}\) as developed in Section IV-D.
We construct two classes of inclusion functions for the closed-loop vector field \(f^{c}\) in (5).

### _Interconnection-based approach_

In the first approach, which we refer to as interconnection-based approach, we obtain the closed-loop inclusion function from a feedback interconnection of the open-loop embedding system and the neural network inclusion function.

**Theorem 2** (Closed-loop interconnection-based inclusion function).: _Consider the open-loop system (3) with an inclusion function \(\mathsf{F}^{\rm o}\) and a neural network controller of the form (4) with a \([y,\widehat{y}]\)-localized inclusion function \(\mathsf{N}_{[y,\widehat{y}]}=\left[\frac{\mathbbm{N}_{[y,\widehat{y}]}}{ \mathbbm{N}_{[y,\widehat{y}]}}\right]\). Then, the map \(\mathsf{F}^{\rm con}_{[y,\widehat{y}]}:\mathcal{T}^{2n}_{\geq 0}\times \mathcal{T}^{2q}_{\geq 0}\rightarrow\mathcal{T}^{2n}_{\geq 0}\) defined by_

\[\frac{\mathsf{F}^{\rm con}_{[y,\widehat{y}]}}{\mathsf{F}^{\rm con }_{[y,\widehat{y}]}}(x,\widehat{x},w,\widehat{w}) =\mathsf{\overline{F}}^{\rm o}(x,\widehat{x},\eta,\widehat{\eta},w, \widehat{w}) \tag{15}\]

_is a \([y,\widehat{y}]\times\mathbb{R}^{q}\)-localized inclusion function for the closed-loop vector field \(f^{c}\), where_

\[\xi =\mathbbm{N}_{[y,\widehat{y}]}(x,\widehat{x}),\quad\widehat{\xi} =\overline{\mathsf{N}}_{[y,\widehat{y}]}(x,\widehat{x}),\] \[\eta =\overline{\mathsf{N}}_{[y,\widehat{y}]}(x,\widehat{x}),\quad \widehat{\eta}=\underline{\mathsf{N}}_{[y,\widehat{y}]}(x,\widehat{x}).\]

Proof.: Let \(z\in[x,\widehat{x}]\subseteq[y,\widehat{y}]\) and \(v\in[w,\widehat{w}]\). Since \(\mathsf{N}_{[y,\widehat{y}]}=\left[\frac{\mathbbm{N}_{[y,\widehat{y}]}}{ \mathbbm{N}_{[y,\widehat{y}]}}\right]\) is a \([y,\widehat{y}]\)-localized inclusion function for the neural network controller \(N\), we have \(\underline{\mathsf{N}}_{[y,\widehat{y}]}(x,\widehat{x})\leq N(z)\leq\overline{ \mathsf{N}}_{[y,\widehat{y}]}(x,\widehat{x})\). Therefore, we can compute \(\mathsf{F}^{\rm con}_{[y,\widehat{y}]}(x,\widehat{x},w,\widehat{w})=\mathsf{ \underline{F}}^{\rm o}(x,\widehat{x},\xi,\xi,w,\widehat{w})\leq f(z,N(z),v)=f ^{c}(z,v)\), where the second inequality holds because \(\mathsf{F}^{\rm o}\) is an inclusion function for the open-loop vector field \(f\). Similarly, one can show that \(\overline{\mathsf{F}}^{\rm con}_{[y,\widehat{y}]}(x,\widehat{x},w,\widehat{w} )\geq f^{c}(z,v)\). This implies that \(\mathsf{F}^{\rm con}_{[y,\widehat{y}]}\) is a \([y,\widehat{y}]\)-localized inclusion function for \(f^{c}\). 

### _Interaction-based approach_

The second approach, which we refer to as interaction-based approach, uses a Jacobian-based cornered inclusion function for the open-loop system and an affine neural network bound to capture the first order interaction between the neural network controller and the system.

**Theorem 3** (Closed-loop interaction-based inclusion function).: _Consider the closed-loop system (5) and assume the open-loop vector field \(f\) is continuously differentiable and the neural network controller (4) satisfies affine bounds of the form (11) on the interval \([y,\widehat{y}]\). Assume that, for every \(z\in[x,\widehat{x}]\subseteq[y,\widehat{y}]\), every \(\xi\in[u,\widehat{u}]\), and every \(\eta\in[w,\widehat{w}]\),_

\[D_{\xi}f(z,\xi,\eta) \in[\underline{J}_{[x,\widehat{x}]},\overline{J}_{[x,\widehat{x}]}],\] \[D_{\xi}f(z,\xi,\eta) \in[\underline{J}_{[u,\widehat{u}]},\overline{J}_{[u,\widehat{u}]}],\] \[D_{\eta}f(z,\xi,\eta) \in[\underline{J}_{[u,\widehat{u}]},\overline{J}_{[w,\widehat{u}]}]. \tag{16}\]

_Then, the map \(\mathsf{F}^{\rm act}_{[y,\widehat{y}]}:\mathcal{T}^{2n}_{\geq 0}\times \mathcal{T}^{2q}_{\geq 0}\rightarrow\mathcal{T}^{2n}_{\geq 0}\) defined by_

\[\mathsf{F}^{\rm act}_{[y,\widehat{y}]}(x,\widehat{x},w,\widehat{w})=\left[ \frac{[\underline{H}]^{+}-\underline{J}_{[x,\widehat{x}]}}{[\overline{H}]^{-} }-\overline{J}_{[x,\widehat{x}]},\overline{[\overline{H}]^{+}}\right]\left[ \frac{x}{\widehat{x}}\right]+L\left[\frac{w}{\widehat{w}}\right]+Q, \tag{17}\]

_is a \([y,\widehat{y}]\times\mathbb{R}^{q}\)-localized inclusion function for the closed-loop vector field \(f^{c}\), where_

\[\frac{H}{H} =\underline{J}_{[x,\widehat{x}]}+[\underline{J}_{[u,\widehat{u}]}] ^{+}\underline{C}(y,\widehat{y})+[\underline{J}_{[u,\widehat{u}]}]^{-} \overline{C}(y,\widehat{y})\] \[\overline{H} =\overline{J}_{[x,\widehat{x}]}+[\overline{J}_{[u,\widehat{u}]}] ^{+}\overline{C}(y,\widehat{y})+[\overline{J}_{[u,\widehat{u}]}]^{-} \underline{C}(y,\widehat{y}),\] \[L =\left[\begin{smallmatrix}-[\underline{J}_{u,\widehat{u}]}]^{-} \underline{[J}_{u,\widehat{u}]}^{-}\\ -[\overline{J}_{[u,\widehat{u}]}]^{+}\overline{[J}_{[u,\widehat{u}]}]^{-} \end{smallmatrix}\right],\] \[Q =\left[\begin{smallmatrix}-\underline{J}_{u,\widehat{u}]}u+[ \underline{J}_{u,\widehat{u}]}u+\underline{d}(y,\widehat{y})+[\underline{J}_{u, \widehat{u}]}]^{-}\overline{d}(y,\widehat{y})+f(x,u,w)\\ -\overline{J}_{[u,\widehat{u}]}u+[\overline{J}_{[u,\widehat{u}]}]^{+} \overline{d}(y,\widehat{y})+[\overline{J}_{[u,\widehat{u}]}]^{-}\underline{d} (y,\widehat{y})+f(x,u,w)\end{smallmatrix}\right].\]

Proof.: Let \(z\in[x,\widehat{x}]\) and \(\eta\in[w,\widehat{w}]\). Since \(f\) is continuously differentiable, the fundamental theorem of Calculus gives

\[f^{c}(z,\eta) =f(z,N(z),\eta)=f(x,u,w)+M(z,u,w)(z-x)\] \[+P(z,u,w)(N(z)-u)+R(z,\eta,w)(\eta-w),\]

where \(M(z,u,w)=\int_{0}^{1}D_{x}f(\tau z+(1-\tau)x,u,w)d\tau\), \(P(z,u,w)=\int_{0}^{1}D_{u}f(z,\tau N(z)+(1-\tau)u,w)d\tau\), and \(R(z,\eta,w)=\int_{0}^{1}D_{w}f(z,N(z),\tau\eta+(1-\tau)w)d\tau\). Using the bounds (16), we get \(M(z,u,w)\geq\underline{J}_{[x,\widehat{x}]}\), \(P(z,u,w)\geq\underline{J}_{[u,\widehat{u}]}\), and \(R(z,\eta,w)\geq\underline{J}_{[w,\widehat{w}]}\). Moreover, we have \(N(z)-u\geq 0_{p}\). Therefore, using the affine bounds (11) for the neural network, \(P(z,u,w)(N(z)-u)\geq[P(z,u,w)]^{+}(\underline{C}(y,\widehat{y})z+\underline{d}(y,\widehat{y})-u)+[P(z,u,w)]^{-}(\overline{C}(y,\widehat{y})z+\underline{d}(y,\widehat{y})-u)\). Using the fact that \(z-x\geq 0_{n}\) and \(\eta-w\geq 0_{q}\), we have

\[f^{c}(z,\eta) \geq f(x,u,w)+M(z,u,w)(z-x)\] \[+[P(z,u,w)]^{+}(\underline{C}(y,\widehat{y})z+\overline{d}(y, \widehat{y})-u)\] \[+[P(z,u,w)]^{-}(\overline{C}(y,\widehat{y})z+\overline{d}(y, \widehat{y})-u)\] \[+R(z,\eta,w)(\eta-w)\] \[\geq\underline{H}z-\underline{J}_{[x,\widehat{x}]}x+\underline{J}_{[ w,\widehat{w}]}\eta-\underline{J}_{[w,\widehat{w}]}w\] \[+[\underline{J}_{[w,\widehat{w}]}]^{+}\underline{d}(y,\widehat{y})+[ \underline{J}_{[w,\widehat{w}]}]^{-}\overline{d}(y,\widehat{y})+f(x,u,w),\]

where the second inequality follows from the Jacobian bounds. Now, we note that \(x\leq z\leq\widehat{x}\) and \(w\leq\eta\leq\widehat{w}\) and therefore \(\underline{H}z\geq[\underline{H}]^{+}x+[\underline{H}]^{+}\widehat{x}\) and
expense. Second, as it is clear from the proof of Theorem 3, the Jacobian-based cornered inclusion function is particularly well-suited for integrating the neural network bounds into the closed-loop inclusion function and for studying the interaction between the system and the neural network controller.
2. _(Mixed Jacobian-based inclusion functions):_ Theorem 3 uses the Jacobian-based cornered inclusion function (as in Proposition 1) for the open-loop vector field \(f\). Alternatively, one can use the mixed Jacobian-based cornered inclusion function (as in Proposition 2) for the open-loop vector \(f\) and obtain a different class of inclusion functions for \(f^{c}\). This class of closed-loop inclusion functions can be obtained from (17) by replacing \(\underline{J}_{[x,\widehat{x}]},\overline{J}_{[x,\widehat{x}]},\underline{J}_ {[u,\widehat{u}]},\overline{J}_{[u,\widehat{w}]},\underline{J}_{[w,\widehat{ w}]}\), and \(\overline{J}_{[w,\widehat{w}]}\) with their counterparts \(\underline{M}_{[x,\widehat{x}]},\overline{M}_{[x,\widehat{x}]},\underline{M} _{[u,\widehat{u}]},\overline{M}_{[u,\widehat{w}]},\) and \(\overline{M}_{[w,\widehat{w}]}\) as defined in Proposition 2 and is demonstrated in the numerical experiments below. Our numerical results show that, in most cases, mixed Jacobian-based cornered inclusion functions significantly reduces the conservatism of the over-approximation as compared to their non-mixed counterparts.
3. _(Comparison):_ The interconnection-based approach provides a general and computationally efficient framework for constructing closed-loop inclusion functions. This approach can be applied using any inclusion function for the open-loop system and any inclusion function of the form (10) for the neural network. The computational efficiency of the interconnection-based approach is a direct consequence of its compositional structure. On the other hand, the interaction-based approach fully captures the first-order interaction between the system and the neural networks using a Jacobian-based decomposition. However, this approach requires affine bounds of the form (11) for the neural network and Jacobian bounds of the form (16) for the open-loop system. As a result, computing the closed-loop inclusion functions using the interaction-based approach is generally more computationally expensive.

When the open-loop system (3) is linear, one can significantly simplify the expression for the closed-loop inclusion functions obtained from the interconnection-based and the interaction-based approaches.

**Corollary 1** (Linear open-loop systems).: _Consider the closed-loop system (5) and assume that the open-loop vector field \(f\) is linear with the form \(f(x,u,w)=Ax+Bu+Dw\), where \(A\in\mathbb{R}^{n\times n}\), \(B\in\mathbb{R}^{n\times p}\) and \(D\in\mathbb{R}^{n\times q}\). Then the following statements hold:_

1. _[label=()]_
2. _using the minimal inclusion function for the open-loop vector field_ \(f\) _and affine neural network inclusion functions of the form (_12_), the closed-loop inclusion function_ \(\mathsf{F}^{\mathrm{con}}_{[y,\widehat{y}]}\) _defined in (_15_) is given by:_ \[\mathsf{F}^{\mathrm{con}}_{[y,\widehat{y}]} =\begin{bmatrix}[A]^{+}+[B^{+}\underline{C}+B^{-}\overline{C}]^{+ }&[A]^{+}-[B^{+}\underline{C}+B^{-}\overline{C}]^{-}\\ [A]^{+}+[B^{+}\overline{C}+B^{-}\underline{C}]^{-}&[A]^{+}+[B^{+}\overline{C} +B^{-}\underline{C}]^{+}\end{bmatrix}\begin{bmatrix}\frac{x}{2}\\ \frac{y}{2}\end{bmatrix}\] \[+\begin{bmatrix}D^{+}\ D^{-}\ \big{[}\frac{w}{2}\big{]}+ \begin{bmatrix}B^{+}\underline{d}+B^{-}\overline{d}\\ \frac{y}{2}\end{bmatrix}.\] (18)
3. _the inclusion function_ \(\mathsf{F}^{\mathrm{act}}_{[y,\widehat{y}]}\) _defined in (_17_) is given by_ \[\mathsf{F}^{\mathrm{act}}_{[y,\widehat{y}]} =\begin{bmatrix}[A]^{+}+B^{+}\underline{C}+B^{-}\overline{C}]^{+ }&[A]^{+}+B^{+}\underline{C}+B^{-}\overline{C}]^{-}\\ [A]^{+}+B^{+}\overline{C}+B^{-}\underline{C}]^{-}&[A]^{+}+B^{+}\overline{C}+B^ {-}\underline{C}]^{+}\end{bmatrix}\begin{bmatrix}\frac{x}{2}\\ \frac{y}{2}\end{bmatrix}\] \[+\begin{bmatrix}D^{+}\ D^{-}\ D^{+}\ \big{[}\frac{w}{2}\big{]}+ \begin{bmatrix}B^{+}\underline{d}+B^{-}\overline{d}\\ \frac{y}{2}\end{bmatrix}.\] (19)

Proof.: Regarding part (i), it is easy to see that the minimal inclusion function of the linear open-loop vector field \(f\) is

\[\mathsf{F}^{\mathrm{o}}(x,\widehat{x},u,\widehat{u},w,\widehat{w}) =\begin{bmatrix}A^{+}\ A^{-}\ A^{+}\ \big{[}\frac{x}{2}\big{]}+\begin{bmatrix}B^{+}\ B^{-}\ B^{-}\end{bmatrix} \begin{bmatrix}\frac{w}{2}\\ \frac{w}{2}\end{bmatrix}\] \[+\begin{bmatrix}D^{+}\ D^{-}\ D^{+}\ \big{[}\frac{w}{2}\big{]}\,. \end{bmatrix}.\]

One can then replace the above minimal inclusion function and the affine bounds (12) into equation (15) and use Theorem 2 to obtain the result. Regarding part (ii), we use the bounds \(\underline{J}_{[x,\widehat{x}]}=\overline{J}_{[x,\widehat{x}]}=A\), \(\underline{J}_{[u,\widehat{u}]}=\overline{J}_{[u,\widehat{u}]}=B\), and \(\underline{J}_{[w,\widehat{w}]}=\overline{J}_{[w,\widehat{w}]}=D\) in Theorem 3. Thus, using the notation of Theorem 3, we compute

\[\underline{H} =A+B^{+}\underline{C}(y,\widehat{y})+B^{-}\overline{C}(y,\widehat {y}),\] \[\overline{H} =A+B^{+}\overline{C}(y,\widehat{y})+B^{-}\underline{C}(y,\widehat {y}),\] \[L =\begin{bmatrix}-D^{-}\ D^{-}\end{bmatrix},\qquad Q=\begin{bmatrix} Ax+Dw+B^{+}\underline{d}(y,\widehat{y})+B^{-}\overline{d}(y,\widehat{y})\\ Ax+Dw+B^{+}\overline{d}(y,\widehat{y})+B^{-}\underline{d}(y,\widehat{y})\end{bmatrix}.\]

The result then follows by replacing the above terms into equation (17) and using Theorem 3. 

**Remark 6** (The role of interactions).: The closed-form expressions (18) and (19) demonstrate the difference between the interconnection-based and the interaction-based approach. In both cases, one can interpret the term \(A\) as the effect of open-loop dynamics and the terms \(B^{+}\underline{C}+B^{-}\overline{C}\) and \(B^{+}\overline{C}+B^{-}\underline{C}\) as the effect of the neural network controller. The interconnection-based approach considers the interconnection of the cooperative and competitive effect of the open-loop system and the neural network controller, separately. This leads to the terms \([A]^{\pm}+[B^{+}\underline{C}+B^{-}\overline{C}]^{\pm}\) and \([A]^{\pm}+[B^{+}\overline{C}+B^{-}\underline{C}]^{\pm}\). On the other hand, the interaction-based approach first considers the interaction between the open-loop system and the neural network controller via the terms \(A+B^{+}\underline{C}+B^{-}\overline{C}\) and \(A+B^{+}\overline{C}+B^{-}\underline{C}\) and then the cooperative and competitive effect of these interactions are separated.

### _Interval Reachability of Closed-loop Systems_

Our culminating conclusion is that, using inclusion functions constructed by the interconnection-based and the interaction-based approaches, we obtain computationally efficient over-approximations of reachable sets of the closed-loop system.

**Theorem 4** (Reachability of closed-loop system).: _Let \(\mathrm{c}\in\{\mathrm{con},\mathrm{act}\}\), the disturbance set \(\mathcal{W}=[\underline{w},\overline{w}]\), and the initial set \(\mathcal{X}_{0}=[\underline{x}_{0},\overline{x}_{0}]\). For the closed-loop dynamical system (5) with \(t\mapsto\left[\frac{\underline{x}^{\mathrm{c}}(t)}{\mathbb{E}^{\mathrm{c}}(t)}\right]\) being the trajectory of the embedding system_

\[\dot{x}_{i} =\left(\mathsf{E}^{\mathrm{c}}_{[x,\widehat{x}]}(x,\widehat{x} _{[i:x]},\underline{w},\overline{w})\right)_{i},\] \[\dot{\hat{x}}_{i} =\left(\overline{\mathsf{F}}^{\mathrm{c}}_{[x,\widehat{x}]}(x_{[ i:\widehat{x}]},\widehat{x},\overline{w},\underline{w})\right)_{i}, \tag{20}\]
_for every \(i\in\{1,\ldots,n\}\), starting from \(\left[\frac{x}{\mathbb{F}_{0}}\right]\), the following statements hold:_

1. _for all_ \(t\in\mathbb{R}_{\geq 0}\)_,_ \[\mathcal{R}_{f^{c}}(t,[\underline{x}_{0},\overline{x}_{0}],[\underline{w}, \overline{w}])\subseteq[\underline{x}^{\mathrm{c}}(t),\overline{x}^{\mathrm{c }}(t)].\]
2. _if the open-loop system (_3_) is linear, then for all_ \(t\in\mathbb{R}_{\geq 0}\)_,_ \[[\underline{x}^{\mathrm{act}}(t),\overline{x}^{\mathrm{act}}(t)]\subseteq[ \underline{x}^{\mathrm{con}}(t),\overline{x}^{\mathrm{con}}(t)].\]

Proof.: Regarding part (i), using Theorem 2 and Theorem 3, we show that \(\mathsf{F}^{\mathrm{con}}_{[y,\widehat{y}]}\) and \(\mathsf{F}^{\mathrm{act}}_{[y,\widehat{y}]}\) are \([y,\widehat{y}]\times\mathbb{R}^{q}\)-localized inclusion functions for the closed-loop vector field \(f^{c}\). Thus, one can use the localized version of Proposition 5 (see Remark 4(i)) with \(y=x\) and \(\widehat{y}=\widehat{x}\) to obtain the result. Regarding part (ii), for linear open-loop systems, using Theorem 1, the embedding system associated to the interconnection-based inclusion function function has the following form:

\[\frac{d}{dt}\left[\frac{x}{\widehat{x}}\right]=\left[\begin{smallmatrix}[A^{ \mathrm{M}}+B^{\mathrm{t}}\underline{C}+B^{-}\overline{C}]^{\mathrm{M}}&[A^{ \mathrm{M}}]^{\mathrm{M}}+[B^{\mathrm{t}}\underline{C}+B^{-}\overline{C}]^{ \mathrm{M}}\\ [A^{\mathrm{M}}+[B^{\mathrm{t}}\underline{C}+\overline{C}]^{\mathrm{M}}]&[A^{ \mathrm{M}}]^{\mathrm{M}}+[B^{\mathrm{t}}\underline{C}+B^{-}\overline{C}]^{ \mathrm{M}}\end{smallmatrix}\right]\left[\frac{x}{\widehat{x}}\right]\] \[+\left[\begin{smallmatrix}D^{+}D^{-}\\ D^{-}D^{+}\end{smallmatrix}\right]\left[\frac{w}{\widehat{w}}\right]+\left[ \begin{smallmatrix}B^{\mathrm{t}}\underline{d}+B^{-}\overline{d}\\ B^{-}\underline{d}+B^{+}\overline{d}\end{smallmatrix}\right]=\mathsf{F}^{ \mathrm{con}}(x,\widehat{x},w,\widehat{w}), \tag{21}\]

and the embedding system associated to the interaction-based inclusion function has the following form:

\[\frac{d}{dt}\left[\frac{x}{\widehat{x}}\right]=\left[\begin{smallmatrix}[A^{ \mathrm{B}}+B^{\mathrm{t}}\underline{C}+B^{-}\overline{C}]^{\mathrm{M}}&[A^{ \mathrm{B}}+\underline{C}+B^{-}\overline{C}]^{\mathrm{M}}\\ [A^{\mathrm{B}}+\underline{C}+B^{-}\overline{C}]^{\mathrm{M}}&[A^{\mathrm{B}}+ B^{-}\overline{C}]^{\mathrm{M}}\end{smallmatrix}\right]\left[\frac{x}{\widehat{x}}\right]\] \[+\left[\begin{smallmatrix}D^{+}D^{-}\\ D^{-}D^{+}\end{smallmatrix}\right]\left[\frac{w}{\widehat{w}}\right]+\left[ \begin{smallmatrix}B^{\mathrm{t}}\underline{d}+B^{-}\overline{d}\\ B^{-}\underline{d}+B^{+}\overline{d}\end{smallmatrix}\right]=\mathsf{F}^{ \mathrm{act}}(x,\widehat{x},w,\widehat{w}). \tag{22}\]

It is easy to check that both embedding systems (21) and (22) are continuous-time monotone. Moreover, we have

\[[A+B^{+}\underline{C}+B^{-}\overline{C}]^{\mathrm{M}}\leq[A]^{ \mathrm{M}}+[B^{+}\underline{C}+B^{-}\overline{C}]^{\mathrm{M}},\] \[[A+B^{+}\underline{C}+B^{-}\overline{C}]^{\mathrm{M}}\geq[A]^{ \mathrm{M}}+[B^{+}\underline{C}+B^{-}\overline{C}]^{\mathrm{nM}}.\]

Therefore, for every \(x\leq\widehat{x}\), we have

\[[A+B^{+}\underline{C}+B^{-}\overline{C}]^{\mathrm{M}}x+[A+B^{+} \underline{C}+B^{-}\overline{C}]^{\mathrm{nM}}\widehat{x}\] \[\geq([A]^{\mathrm{M}}+[B^{+}\underline{C}+B^{-}\overline{C}]^{ \mathrm{M}})x+([A]^{\mathrm{nM}}+[B^{+}\underline{C}+B^{-}\overline{C}]^{ \mathrm{nM}})\widehat{x}.\]

This implies that, \(\mathsf{F}^{\mathrm{con}}(x,\widehat{x},w,\widehat{w})\leq_{\mathrm{SE}} \mathsf{F}^{\mathrm{act}}(x,\widehat{x},w,\widehat{w})\), for every \(x\leq\widehat{x}\) and every \(w\leq\widehat{w}\). Now, we can use the classical monotone comparison Lemma [48, Theorem 3.8.1] to obtain \(\left[\frac{x^{\mathrm{con}}(t)}{\overline{x}^{\mathrm{con}}(t)}\right]\leq_{ \mathrm{SE}}\left[\frac{x^{\mathrm{an}}(t)}{\overline{x}^{\mathrm{an}}(t)}\right]\) for every \(t\in\mathbb{R}_{\geq 0}\). As a result, we get \([\underline{x}^{\mathrm{act}}(t),\overline{x}^{\mathrm{act}}(t)]\subseteq[ \underline{x}^{\mathrm{con}}(t),\overline{x}^{\mathrm{con}}(t)]\), for every \(t\in\mathbb{R}_{\geq 0}\). 

**Remark 7**.:
1. _(Comparison with the literature):_ For linear open-loop systems, a straightforward computation shows that Theorem 4 with \(c=\mathrm{act}\) coincides with [30, Lemma IV.3] with \(p=\infty\) and \(q=1\).
2. _(Computational complexity):_ from a computational perspective, the framework presented in Theorem 4 consists of two main ingredients: (i) the neural network verification algorithm for computing the inclusion function for the neural network (either in the form (10) or in the form (12)), and (ii) integration of a single trajectory of the embedding system (13). Part (i) includes querying the neural network verification once per integration step and its runtime depends on the computational complexity of the associated algorithm. The runtime of the part (ii) depends on the integration method and the form of the open-loop decomposition function \(\mathsf{F}\).

## VII Numerical Experiments

We demonstrate the effectiveness of our reachability analysis, implemented as an open-source Python toolbox called ReachMM2, using numerical experiments3 for (i) a nonlinear bicycle model, (ii) the linear double integrator, (iii) several existing benchmarks in the literature, and (iv) a platoon of double integrator vehicles moving in a plane. We first briefly mention several algorithmic techniques that are implemented in the ReachMM toolbox and are used in several examples below to broaden the applicability and to improve the accuracy of our reachable set over-approximations.

Footnote 2: [https://github.com/gftactslab/ReachMM_TAC2023](https://github.com/gftactslab/ReachMM_TAC2023)

Footnote 3: All the experiments are performed using an AMD Ryzen 5 5600X CPU and 32 GB of RAM. Unless otherwise specified, runtimes are averaged over 100 runs, with mean and standard deviations reported.

_Partitioning:_ The computational efficiency of our method makes it well-suited for partitioning to reduce compounding over-approximation error caused by the wrapping effect [10, Section 2.2.4]. In Section VII-B, we apply the partitioning methods proposed in [38, 52] to improve accuracy.

_Redundant Variable Refinement:_ We accommodate a technique introduced in [13] to refine interval over-approximations of reachable sets using redundant state variables of the form \(y=Ax+b\) for some \(A\in\mathbb{R}^{m\times n}\) and \(b\in\mathbb{R}^{m}\). By augmenting the dynamical system (5) with the additional dynamics \(\dot{y}=A\dot{x}\), we create a new dynamical system with the new states \(z=\left[\begin{smallmatrix}x\\ y\end{smallmatrix}\right]\in\mathbb{R}^{n+m}\) and the affine constraints \(Mz=b\) where \(M=\left[-A&I_{m}\right]\in\mathbb{R}^{m\times(n+m)}\). Following [13], we use the interval bounds on \(z\) to improve the accuracy of the interval bounds on \(x\). We apply this technique in Section VII-C for the TORA benchmark.

_Zero-order hold control:_ In some applications, we wish to explicitly account for the practical restriction that the control \(u(t)=N(x(t))\) cannot be continuously updated and, instead, the control must be implemented via, _e.g._, a zero-order hold strategy between sampling instances. It is straightforward to extend the interconnection-based approach to be able to capture the zero-order hold policy [52]. Adopting the zero-order hold policy in the interaction-based approach requires over-bounding the error and we omit this calculation for the sake of brevity. We refer to our code for details of this adaptation. We use this framework in Section VII-C.

_Implementation:_ Our current implementation for all examples is in standard Python. For the interconnection-based approach, the inclusion function for the open-loop system is the natural inclusion functions computed using our package npinterval [41]. For neural networks, the affine inclusion functions are obtained from CROWN [25] and computed using autoLiRPA [53]. Further performance improvements of ReachMM are likely possible by, _e.g._, implementing as compiled code, and are subject of ongoing and future work.

### _Nonlinear bicycle model_

In the first experiment, we compare the interconnection-based and the interaction-based approach. Consider the non
linear dynamics of a bicycle adopted from [54]:

\[\dot{p_{x}} =v\cos(\phi+\beta(u_{2})) \dot{\phi} =\frac{v}{\ell_{r}}\sin(\beta(u_{2}))\] \[\dot{p_{y}} =v\sin(\phi+\beta(u_{2})) \dot{v} =u_{1} \tag{23}\]

where \([p_{x},p_{y}]^{\top}\in\mathbb{R}^{2}\) is the displacement of the center of mass in the \(x-y\) plane, \(\phi\in[-\pi,\pi)\) is the heading angle in the plane, and \(v\in\mathbb{R}_{\geq 0}\) is the speed of the center of mass. Control input \(u_{1}\) is the applied force, input \(u_{2}\) is the angle of the front wheel, and \(\beta(u_{2})=\arctan\left(\frac{\ell_{f}}{\ell_{f}+\ell_{r}}\tan(u_{2})\right)\) is the slip slide angle where the parameter \(\ell_{f}\) (\(\ell_{r}\)) is the distance between the center of mass and the front (rear) wheel. In this example, for the sake of simplicity, we set \(\ell_{f}=\ell_{r}=1\). Let \(x=[p_{x},p_{y},\phi,v]^{\top}\) and \(u=[u_{1},u_{2}]^{\top}\). We apply the neural network controller (\(4\times 100\times 100\times 2\), ReLU activations) defined in [38], which was trained to mimic an MPC that stabilizes the vehicle to the origin while avoiding a circular obstacle centered at \((4,4)\) with a radius of \(2\). The dynamics are simulated using Euler integration with a step size of \(0.125\). In Figure 3, we compare the accuracy and runtime of different reachability approaches for the bicycle model.

DiscussionFigure 3 shows that the naive input-output approach is the fastest approach with low accuracy of over-approximation. Using the interconnection-based approach with \(\mathsf{F}^{\rm con}\) improves the accuracy with a slight increase in the runtime. In comparison, the interaction-based approach with the Jacobian-based cornered inclusion function \(\mathsf{F}^{\rm act}\) is notably slower due to the computation of Jacobian bounds (16) in this approach. Over short time horizons, the interaction-based approach is more accurate, however, its accuracy deteriorates for longer horizons, which can be attributed to the decrease in accuracy of the Jacobian bounds (16). As has been shown in Proposition 4, the accuracy of the intersection \(\mathsf{F}^{\rm con}\wedge\mathsf{F}^{\rm act}\) is better than both \(\mathsf{F}^{\rm con}\) and \(\mathsf{F}^{\rm act}\). Finally, using the mixed Jacobian-based cornered inclusion function (cf. Remark 5(ii)) significantly improves the accuracy of the interaction-based approach with little effect on its runtime.

### _Double Integrator Model_

In the second experiment, we focus on reachability of linear open-loop systems with neural network controllers. We study the accuracy and efficiency of the interaction-based approach (19) for the double integrator benchmark system

\[x(t+1)=\begin{bmatrix}1&1\\ 0&1\end{bmatrix}x(t)+\begin{bmatrix}0.5\\ 1\end{bmatrix}u(t). \tag{24}\]

For this example, we use the discrete-time version of our framework as derived in [52, Section VII.B]. We apply the neural network controller (\(2\times 10\times 5\times 1\), ReLU activations) defined in [30] and consider the interaction-based approach with the closed-loop inclusion function \(\mathsf{F}^{\rm act}\) defined in (19). We use contraction-guided adaptive partitioning proposed in [52] to improve the accuracy. Additionally, we compare our proposed ReachMM to state-of-the-art algorithms for linear discrete-time systems: ReachLP [30] with uniform partitioning (ReachLP-Unif) and greedy simulation guided partitioning (ReachLP-GSG), and ReachLipBnB [55] (branch-and-bound using LipSDP [24]). Each algorithm is run with two different sets of hyper-parameters, aiming to compare their performances across various regimes. The setup for ReachMM is \((\varepsilon,\,D_{p},\,D_{\text{N}})\); ReachLP-Unif is \(\#\) initial partitions, ReachLP-GSG is \(\#\) of total propogator calls, ReachLipBnB is \(\varepsilon\). The results are shown in Figure 4 and the runtime comparisons are shown in Table II.

DiscussionFigure 4 and Table II show that, when the open-loop system is linear, our interaction-based approach (Corollary 1(ii)) combined with a suitable partitioning scheme beats state-of-the-art approaches in both accuracy and runtime.

Fig. 4: The over-approximated reachable sets of the closed-loop double integrator model (24) are compared for three different algorithms on two different runtime regimes, for the initial set \([2.5,3]\times[-0.25,0.25]\) and final time \(T=5\). The experiment setup and performances are reported in Table II. 200 true trajectories are shown in red. The horizontal axis is \(x_{1}\) and the vertical axis is \(x_{2}\).

Fig. 3: Accuracy and runtime comparison between different reachability approaches for the bicycle model (23) from the initial set \([7.95,8.05]\times[6.95,7.05]\times[-2\pi/3-0.01,-2\pi/3+0.01]\times[1.99,2.01]\): (a) a naive input-output approach combining the natural inclusion function for the open-loop dynamics and the affine inclusion function for the neural network, (b) the interconnection-based approach with the closed-loop inclusion function \(\mathsf{F}^{\rm con}\) defined in (15) constructed from the natural inclusion function for the open loop dynamics and affine inclusion function for the neural network, (c) the interaction-based approach with Jacobian-based cornered inclusion function \(\mathsf{F}^{\rm act}\) defined in (17), (d) the intersection of the interconnection-based inclusion and interaction-based approach \(\mathsf{F}^{\rm con}\wedge\mathsf{F}^{\rm act}\), (e) the interaction-based approach with _mixed states_ Jacobian-based cornered inclusion function defined in (17) and Remark 5(ii), and (f) the interaction-based approach with _mixed states and control_ Jacobian-based cornered inclusion function defined (17) and Remark 5(ii). The blue boxes are hyper-rectangular over-approximation of reachable sets and 100 simulated trajectories of the system are shown in red.
### _ARCH-COMP Benchmarks_

In the third experiment, we analyze three of the benchmarks from ARCH-COMP22 [56].

Adaptive Cruise Control (ACC)We consider the Adaptive Cruise Control benchmark from [56, Equation (1)]. The neural network controller (\(5\times 20\times 20\times 20\times 20\times 20\times 1\), ReLU activations) with the input \((x_{\rm lead}-x_{\rm ego},v_{\rm lead}-v_{\rm ego},a_{\rm lead}-a_{\rm ego},v_{ \rm set},T_{\rm gap})^{\top}\in\mathbb{R}^{5}\) is applied with a zero-order holding of \(0.1\) seconds [56]. Our goal is to verify that from the initial set \(\mathcal{X}_{0}=[90,110]\times[32,32.2]\times[0,0]\times[10,11]\times[30,30.2 ]\times[0,0]\), the collision specification

\[x_{\rm lead}-x_{\rm ego}\geq D_{\rm default}+T_{\rm gap}x_{\rm ego}\]

is never violated in the next \(5\) seconds, where \(D_{\rm default}=10\) and \(T_{\rm gap}=1.4s\). To verify the specification, we define \(D_{\rm rel}=x_{\rm lead}-x_{\rm ego}\) and \(D_{\rm safe}=D_{\rm default}+T_{\rm gap}x_{\rm ego}\) and ensure that \(D_{\rm rel}-D_{\rm safe}\geq 0\) for the given time horizon.

We consider the interconnection-based approach with the closed-loop inclusion function \(\mathsf{F}^{\rm con}\) defined in (15) constructed from the natural inclusion function for open-loop system and an affine inclusion function for the neural network. We use this interconnection-based inclusion function and the Euler integration with a step-size of \(0.01\) to compute upper and lower bounds on the states of the system starting from the initial set \(\mathcal{X}_{0}\). These bounds are then used to obtain upper and lower bounds on \(D_{\rm rel}\) and \(D_{\rm safe}\). The results are shown in Figure 5.
the state-of-the-art methods POLAR [35] and JuliaReach [37]. On the other benchmarks in [56], our method is unable to verify the desired specification. In the next section, we show that our interconnection-based and interaction-based methods enjoy better scalability compared to POLAR and JuliaReach.

### _Vehicle Platooning_

In the fourth experiment, we investigate the scalability of our method. We consider a platoon of \(N\) vehicles \(\mathcal{V}=\{\mathcal{V}_{j}\}_{j=1}^{N}\), each with the two-dimensional double integrator dynamics

\[\dot{p}_{x}^{j}=v_{x}^{j}, \dot{v}_{x}^{j}=\sigma(u_{x}^{j})+w_{x}^{j},\] \[\dot{p}_{y}^{j}=v_{y}^{j}, \dot{v}_{y}^{j}=\sigma(u_{y}^{j})+w_{y}^{j}, \tag{26}\]

where \(p^{j}=(p_{x}^{j},p_{y}^{j})\in\mathbb{R}^{2}\) is the displacement of the center of mass of \(\mathcal{V}_{j}\), \(v^{j}=(v_{x}^{j},v_{y}^{j})\in\mathbb{R}^{2}\) is the velocity of the center of mass of \(\mathcal{V}_{j}\), \((u_{x}^{j},u_{y}^{j})\in\mathbb{R}^{2}\) are desired acceleration inputs limited by the softmax operator \(\sigma(u)=u_{\text{lim}}\tanh(u/u_{\text{lim}})\) with \(u_{\text{lim}}=5\), and \(w_{x}^{j},w_{y}^{j}\sim\mathcal{U}([-0.001,0.001])\) are uniformly distributed disturbances. We consider a leader-follower structure for the system, where the first vehicle \(\mathcal{V}_{1}\) chooses its control \(u=(u_{x}^{1},u_{y}^{1})\) as the output of a neural network (\(4\times 100\times 100\times 2\), ReLU activations), and the rest of the platoon \(\{\mathcal{V}_{j}\}_{j=2}^{N}\) applies a PD tracking control input

\[u_{\text{d}}^{j}=k_{p}\left(p_{\text{d}}^{j-1}-p_{\text{d}}^{j}-r\frac{v_{ \text{d}}^{j-1}}{\|v^{j-1}\|_{2}}\right)+k_{v}(v_{\text{d}}^{j-1}-v_{\text{d}} ^{j}), \tag{27}\]

for each \(\mathbf{d}\in\{x,y\}\) with \(k_{p}=k_{v}=5\) and \(r=0.5\). The neural network was trained by gathering data from an offline MPC control policy for the leader only (\(N=1\)). The offline policy minimized a quadratic cost aimed at stabilizing to the origin while avoiding a circular obstacle centered at \((4,4)\) with radius \(2.25\), implemented as a hard constraint with \(33\%\) padding and a slack variable.

We consider the interconnection-based approach with closed-loop inclusion function \(\mathsf{F}^{\mathrm{con}}\) defined in (15) constructed from the natural inclusion function for open-loop system. We also consider the interaction-based approach with the closed-loop mixed Jacobian-based cornered inclusion function \(\mathsf{F}^{\mathrm{act}}\) as defined in (17) constructed for four corners \((x,u),(u,\widehat{u}),(\widehat{x},u),(\widehat{x},\widehat{u})\) (cf. Remark 5). We perform reachability analysis for platoons of \(N\) vehicles with \(N\in\{1,4,9,20,50\}\). For \(j\)th vehicle in the platoon, we compute the reachable sets for the time frame \(t\in[0,1.5]\) using the closed-loop inclusion functions \(\mathsf{F}^{\mathrm{con}}\) and \(\mathsf{F}^{\mathrm{act}}\) starting from the initial set \(([7.225+0.5(j-1)\cos(\pi/3),7.275+0.5(j-1)\cos(\pi/3)]\times[5.725+0.5(j-1) \sin(\pi/3),5.775+0.5(j-1)\sin(\pi/3)]\times[-0.5,-0.5]\times[-5,-5])\). All integration is performed using Euler discretization with a step size of \(0.0125\). Figure 8 visualizes the result for \(N=9\). Runtimes of our approach as well as POLAR [35] and JuliaReach [37] for platoons of varying size are reported in IV. Note that the implementations in POLAR and JuliaReach omit the disturbances \(w_{\text{d}}^{j}\) and the softmax \(\sigma\) in equation (26).

DiscussionThis experiment demonstrates one of the key features of our interval analysis framework--its scalability to large-scale systems. In general, Figure 8 shows how \(\mathsf{F}^{\mathrm{act}}\) outperforms \(\mathsf{F}^{\mathrm{con}}\) in accuracy of the reachable set over-approximation. However, Table IV shows how \(\mathsf{F}^{\mathrm{con}}\) outperforms \(\mathsf{F}^{\mathrm{act}}\) in terms of scalability with respect to state dimensions. Moreover, both \(\mathsf{F}^{\mathrm{con}}\) and \(\mathsf{F}^{\mathrm{act}}\) demonstrate better runtime scalability than POLAR and JuliaReach.

Fig. 8: **Left:** The over-approximation of the reachable sets computed using \(\mathsf{F}^{\mathrm{con}}\) for each individual unit in the platooning example with \(N=9\) are shown on \(p_{x}\)-\(p_{y}\) planes, with the circular obstacle in red. **Right:** The over-approximation of the reachable sets computed using \(\mathsf{F}^{\mathrm{act}}\) for each individual unit in the platooning example with \(N=9\) are shown on \(p_{x}\)-\(p_{y}\) planes, with the circular obstacle in red.

Fig. 7: The reachable set for the TORA benchmark is shown on the \(x_{1}\)-\(x_{2}\) plane. The goal set is pictured in green, and 100 true trajectories of the system are pictured in red. **Left:** Every tenth over-approximation of the reachable set is pictured in blue. **Right:** The plot is zoomed into the goal set, and all over-approximations of the reachable set are shown in blue.
## VIII Conclusions

We present a framework based on interval analysis for safety verification of neural network controlled systems. The main idea is to embed the closed-loop system into a higher dimensional system whose over-approximation of reachable sets can be easily computed using its extreme trajectory. Using an inclusion function of the open-loop system and interval bounds for neural networks, we proposed an interconnection-based approach and an interaction-based approach for constructing the closed-loop embedding system. We show how these approaches can capture the interactions between the nonlinear plant and the neural network controllers and we provide several numerical experiments comparing them with the state-of-the-art algorithms in the literature.
