Interval Reachability of Nonlinear Dynamical Systems with Neural Network
  Controllers
This paper proposes a computationally efficient framework, based on interval
analysis, for rigorous verification of nonlinear continuous-time dynamical
systems with neural network controllers. Given a neural network, we use an
existing verification algorithm to construct inclusion functions for its
input-output behavior. Inspired by mixed monotone theory, we embed the
closed-loop dynamics into a larger system using an inclusion function of the
neural network and a decomposition function of the open-loop system. This
embedding provides a scalable approach for safety analysis of the neural
control loop while preserving the nonlinear structure of the system.
  We show that one can efficiently compute hyper-rectangular
over-approximations of the reachable sets using a single trajectory of the
embedding system. We design an algorithm to leverage this computational
advantage through partitioning strategies, improving our reachable set
estimates while balancing its runtime with tunable parameters. We demonstrate
the performance of this algorithm through two case studies. First, we
demonstrate this method's strength in complex nonlinear environments. Then, we
show that our approach matches the performance of the state-of-the art
verification algorithm for linear discretized systems.
Saber Jafarpour, Akash Harapanahalli, Samuel Coogan
# Interval Reachability of Nonlinear Dynamical Systems with Neural Network Controllers

###### Abstract

This paper proposes a computationally efficient framework, based on interval analysis, for rigorous verification of nonlinear continuous-time dynamical systems with neural network controllers. Given a neural network, we use an existing verification algorithm to construct inclusion functions for its input-output behavior. Inspired by mixed monotone theory, we embed the closed-loop dynamics into a larger system using an inclusion function of the neural network and a decomposition function of the open-loop system. This embedding provides a scalable approach for safety analysis of the neural control loop while preserving the nonlinear structure of the system.

We show that one can efficiently compute hyper-rectangular over-approximations of the reachable sets using a single trajectory of the embedding system. We design an algorithm to leverage this computational advantage through partitioning strategies, improving our reachable set estimates while balancing its runtime with tunable parameters. We demonstrate the performance of this algorithm through two case studies. First, we demonstrate this method's strength in complex nonlinear environments. Then, we show that our approach matches the performance of the state-of-the art verification algorithm for linear discretized systems.

Reachability analysis, Neural feedback loop, Verification of neural networks, Safety verification.

## 1 Introduction

Neural network components are increasingly deployed as controllers in safety-critical applications such as self-driving vehicles and robotic systems. For example, they may be designed based on reinforcement learning algorithms (Zhang et al., 2016) or trained to approximate some dynamic optimization-based controllers (Chen et al., 2018). However, neural networks are known to be vulnerable to small input perturbations (Szegedy et al., 2014); a slight disturbance in their input can lead to a large change in their output. In many applications, these learning-based components are interconnected with nonlinear and time-varying systems. Moreover, they are usually trained with no safety and robustness guarantees. Their learned control policies can therefore suffer from a significant degradation in performance when presented with uncertainties not accounted for during training.
Related works.There is extensive literature on verification of isolated neural networks. Rigorous verification approaches generally fall into three different categories: (i) reachability-based methods which focus on layer-by-layer estimation of reachable sets using interval bound propagation (Mirman et al., 2018; Gowal et al., 2018; Wang et al., 2018), activation function relaxation (Zhang et al., 2018), and symbolic interval analysis (Wang et al., 2018); (ii) optimization-based methods which use linear programming (Wong and Kolter, 2018), semi-definite programming (Fazlyab et al., 2019), or search and optimization (Katz et al., 2017) to estimate the input-output behavior of the neural networks; and (iii) probabilistic methods (Cohen et al., 2019; Li et al., 2019). We refer to the survey (Liu et al., 2021) for a review of neural network verification algorithms. Reachability of nonlinear dynamical systems has been studied using optimization-based methods such as the Hamilton-Jacobi approach (Bansal et al., 2017) and the level set approach (Mitchell and Tomlin, 2000). Several computationally tractable approaches including ellipsoidal (Kurzhanski and Varaiya, 2000) and zonotope methods (Girard, 2005) have been developed for reachability analysis of linear systems. Mixed monotone theory has emerged as a computationally efficient framework for over-approximating the reachable sets of nonlinear systems (Coogan and Arcak, 2015).

Recently, there has been an increased interest in verification methods for closed-loop systems with neural network controllers. However, a direct combination of state-of-the art neural network verification algorithms with the existing reachability analysis toolboxes and can lead to overly-conservative estimates of reachable sets (Dutta et al., 2019, Section 2.1). For linear discrete-time systems, reachable set over-approximation has been studied using semi-definite programming (Hu et al., 2020) and linear programming (Everett et al., 2021, 2021). For nonlinear discrete-time systems, Sidrane et al. (2022) establishes a mixed integer programming framework for reachable set over-approximation using polynomial bounds of the system dynamics. For nonlinear continuous-time systems, Dutta et al. (2019) uses polynomial bounding on the dynamics as well as the neural network to over-approximate reachable sets. Other rigorous approaches for verification of closed-loop neural network controllers include using finite-state abstractions (Sun et al., 2019) and simulation-guided interval analysis (Xiang et al., 2021).

Contributions.In this paper, we use elements of mixed monotone system theory to develop a flexible framework for safety verification of nonlinear continuous-time systems with neural network controllers. First, we employ CROWN--a well-established isolated neural network verification framework (Zhang et al., 2018)--to construct an inclusion function for a pre-trained neural network's output. Then, we use this inclusion function and a decomposition function of the open-loop system to construct suitable embedding systems for the closed-loop system using twice the number of states. Finally, we simulate trajectories of these embedding systems to compute hyper-rectangular over-approximations of the closed-loop reachable sets. Our framework has several advantageous features. It is agnostic to the neural network verifier, only requiring an inclusion function for the neural network. Additionally, our approach is fast and scalable, as only a single simulation of the embedding system is required. This feature, combined with a clever partitioning of the state space, is used to improve the accuracy of our reachable set over-approximations while retaining computational viability. Through several numerical experiments, we study the efficiency of our approach for obstacle avoidance in a simple vehicle model and for perturbation analysis in a linear quadrotor model.
## 2 Notation

For every \(x\in\mathbb{R}^{n}\) and every \(r\in\mathbb{R}_{\geq 0}\) we define the closed ball \(B_{\infty}(x,r)=\{y\in\mathbb{R}^{n}\mid\|y-x\|_{\infty}\leq r\}\). The partial order \(\leq\) on \(\mathbb{R}^{n}\) is defined by \(x\leq y\) if and only if \(x_{i}\leq y_{i}\), for every \(i\in\{1,\ldots,n\}\). For every \(x\leq y\), we can define the interval \([x,y]=\{z\in\mathbb{R}^{n}\mid x\leq z\leq y\}\). The partial order \(\leq\) on \(\mathbb{R}^{n}\) induces the southeast partial order \(\leq_{\mathrm{SE}}\) on \(\mathbb{R}^{2n}\) defined by \([\begin{smallmatrix}\widehat{x}\\ \widehat{y}\end{smallmatrix}]\) if \(x\leq y\) and \(\widehat{y}\leq\widehat{x}\). We define the subsets \(\mathcal{T}_{\geq 0}^{2n},\mathcal{T}_{\leq 0}^{2n},\mathcal{T}^{2n}\subseteq \mathbb{R}^{2n}\) as follows:

\[\mathcal{T}_{\geq 0}^{2n}=\{\lfloor\tfrac{x}{\widehat{x}}\rfloor\in \mathbb{R}^{2n}\mid x\leq\widehat{x}\},\qquad\mathcal{T}_{\leq 0}^{2n}=\{ \lfloor\tfrac{x}{\widehat{x}}\rfloor\in\mathbb{R}^{2n}\mid x\geq\widehat{x}\},\qquad\mathcal{T}^{2n}=\mathcal{T}_{\geq 0}^{2n}\cup\mathcal{T}_{\leq 0}^{2n}.\]

For every two vectors \(v,w\in\mathbb{R}^{n}\) and every \(i\in\{1,\ldots,n\}\), we define the vector \(v_{[i:w]}\in\mathbb{R}^{n}\) by \(\left(v_{[i:w]}\right)_{j}=\begin{cases}v_{j}&j\neq i\\ w_{j}&j=i.\end{cases}\), for every \(j\in\{1,\ldots,n\}\). 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)\). The Metzler and non-Metzler part of square matrix \(A\in\mathbb{R}^{n\times n}\) are denoted by \(\lceil A\rceil^{\mathrm{Mzl}}\) and \(\lfloor A\rfloor^{\mathrm{Mzl}}\), respectively, where \((\lceil A\rceil^{\mathrm{Mzl}})_{ij}=\begin{cases}A_{ij}&A_{ij}\geq 0\text{ or }i=j\\ 0&\text{otherwise,}\end{cases}\) and \(\lfloor A\rfloor^{\mathrm{Mzl}}=A-\lceil A\rceil^{\mathrm{Mzl}}\). Consider a control system \(\dot{x}=f(x,u)\) on \(\mathbb{R}^{n}\) with a measurable input \(\mathbf{u}:\mathbb{R}_{\geq 0}\to\mathbb{R}^{p}\). We denote the trajectory of the system with the control input \(\mathbf{u}\) starting from \(x_{0}\in\mathbb{R}^{n}\) at time \(t_{0}\) by \(t\mapsto\phi_{f}(t,t_{0},x_{0},u)\). Given a map \(f:\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}^{n}\) and a Lipschitz function \(\mathsf{F}:\mathcal{T}^{2n}\times\mathcal{T}^{2m}\to\mathbb{R}^{n}\), we say that \(\mathsf{F}\) is a _decomposition function for \(f\)_ if

1. \(\mathsf{F}_{i}(x,x,u,u)=f_{i}(x,u)\), for every \(x\in\mathbb{R}^{n}\) and every \(u\in\mathbb{R}^{p}\);
2. \(\mathsf{F}_{i}(x,\widehat{x},u,\widehat{u})\leq\mathsf{F}_{i}(y,\widehat{y}, u,\widehat{u})\), for every \(x\leq y\) such that \(x_{i}=y_{i}\), and every \(\widehat{y}\leq\widehat{x}\);
3. \(\mathsf{F}_{i}(x,\widehat{x},u,\widehat{u})\leq\mathsf{F}_{i}(x,\widehat{x}, v,\widehat{v})\), for every \(u\leq v\) and every \(\widehat{v}\leq\widehat{u}\).

A decomposition function \(\mathsf{F}\) is _tight_, if for every other decomposition function \(d\), we have

\[\left[\begin{smallmatrix}d(x,\widehat{x},u,\widehat{u})\\ d(\widehat{x},x,\widehat{u},u)\end{smallmatrix}\right]\leq_{\mathrm{SE}}\left[ \begin{smallmatrix}\mathsf{F}(x,\widehat{x},u,\widehat{u})\\ \mathsf{F}(\widehat{x},x,\widehat{u},u)\end{smallmatrix}\right],\qquad\text{ for every }x\leq\widehat{x},\ w\leq\widehat{w}.\]

Given a map \(g:\mathbb{R}^{n}\to\mathbb{R}^{m}\), the function \(\left[\frac{\mathsf{G}}{\mathsf{G}}\right]:\mathcal{T}_{\geq 0}^{2n}\to\mathcal{T}_{ \geq 0}^{2m}\) is _an inclusion function for \(g\)_ if, for every \(x\leq\widehat{x}\) and every \(z\in[x,\widehat{x}]\), we have \(\mathsf{G}(x,\widehat{x})\leq g(z)\leq\overline{\mathsf{G}}(x,\widehat{x})\). Note that our definition of inclusion function is closely-related to the notion of _inclusion interval function_ in the interval analysis (Jaulin et al., 2001, Section 2.4.1).

## 3 Problem Statement

We consider a nonlinear plant of the form

\[\dot{x}=f(x,u,w) \tag{1}\]

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

\[\xi^{(i)}(x) =\phi_{i-1}(W^{(i-1)}\xi^{(i-1)}(x)+b^{(i-1)}),\quad i\in\{1, \ldots,k\}\] \[x =\xi^{(0)}\qquad u=W^{(k)}\xi^{(k)}(x)+b^{(k)}:=\mathsf{N}(x), \tag{2}\]
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, and \(\xi^{(i)}(y)\in\mathbb{R}^{n_{i}}\) is the \(i\)-th layer hidden variable and the activation function \(\phi_{i}\) satisfies \(0\leq\frac{\phi_{i}(x)-\phi_{i}(y)}{x-y}\leq 1\). One can show that a large class of activation function including, but not restricted to, ReLU, leaky ReLU, sigmoid, and tanh satisfies this condition (after a possibly re-scaling of their co-domain). We assume that, the neural network (2) is trained to approximate an offline controller with the auxiliary objective of achieving a goal set \(\mathcal{G}\subseteq\mathbb{R}^{n}\) while remaining in a safe set \(\mathcal{S}\subseteq\mathbb{R}^{n}\). Our aim is to consider the closed-loop system given by:

\[\dot{x}=f(x,\mathsf{N}(x),w):=f^{\mathrm{cl}}(x,w). \tag{3}\]

and to verify its safety, i.e., to ensure that the closed-loop system avoids the unsafe domain \(\mathbb{R}^{n}/\mathcal{S}\). We define the reachable set of the closed-loop system (3) by:

\[\mathcal{R}(t,\mathcal{X}_{0},\mathcal{W})=\{\phi_{f^{\mathrm{cl}}}(t,0,x_{0}, \mathbf{w})\mid x_{0}\in\mathcal{X}_{0},\ \ \mathbf{w}:\mathbb{R}_{\geq 0}\to \mathcal{W}\ \ \ \text{is piecewise continuous}\}\]

Thus, our goal is to check if \(\mathcal{R}(t,\mathcal{X}_{0},\mathcal{W})\subseteq\mathcal{S}\) holds for every \(t\in\mathbb{R}_{\geq 0}\). In general, computing the exact reachable sets of the closed-loop system (3) is not computationally tractable. Our approach is based on constructing a computationally efficient over-approximation \(\overline{\mathcal{R}}(t,\mathcal{X}_{0},\mathcal{W})\) of the reachable set of the closed-loop system. Then, avoiding the unsafe set \(\mathbb{R}^{n}/\mathcal{S}\) is guaranteed when \(\overline{\mathcal{R}}(t,\mathcal{X}_{0},\mathcal{W})\subseteq\mathcal{S}\), for every \(t\in\mathbb{R}_{\geq 0}\).

## 4 Input-output reachability of neural networks

In order to estimate the input-output behavior of the neural network controller, we use the verification algorithm called CROWN (Zhang et al., 2018). We first use CROWN to obtain an inclusion function for the input-output map of the neural network. Consider the neural network (2) and let the perturbed input vector \(x\) be in the interval \([\underline{x},\overline{x}]\). For every \(i\in\{1,\ldots,k\}\), we define the pre-activation input to \(i\)th layer as \(z^{(i)}=W^{(i)}\xi^{(i-1)}(x)+b^{(i)}\). When \(x\) is perturbed in the interval \([\underline{x},\overline{x}]\), we assume that \(L^{(i)},U^{(i)}\in\mathbb{R}^{n_{i}}\) are such that \(L^{(i)}\leq z^{(i)}\leq U^{(i)}\). In this case, for the \(j\)th neuron in the \(i\)th layer, there exist \(\alpha^{(i)}_{U,j},\beta^{(i)}_{U,j},\alpha^{(i)}_{L,j},\beta^{(i)}_{L,j}\) such that

\[\alpha^{(i)}_{L,j}(z+\beta^{(i)}_{L,j})\leq\phi_{i}(z)\leq\alpha^{(i)}_{U,j}(z +\beta^{(i)}_{U,j}),\ \ \ \ \text{for every }L^{(i)}_{j}\leq z\leq U^{i}_{j}. \tag{4}\]

For every \(y\in[\underline{y},\overline{y}]\), the output of the neural network is bounded as below:

\[\underline{A}(\underline{x},\overline{x})x+\underline{b}(\underline{x}, \overline{x})\leq\mathsf{N}(x)\leq\overline{A}(\underline{x},\overline{x})x+ \overline{b}(\underline{x},\overline{x}), \tag{5}\]

where matrix-valued functions \(\underline{A},\overline{A}:\mathcal{T}_{\geq 0}^{2n}\to\mathbb{R}^{p\times n}\) and two vector-valued functions \(\underline{b},\overline{b}:\mathcal{T}_{\geq 0}^{2n}\to\mathbb{R}^{p}\) are defined as follows:

\[\overline{A}(\underline{x},\overline{x})=\Lambda^{(0)},\qquad \overline{b}(\underline{x},\overline{x})=\sum\nolimits_{i=1}^{k}\Lambda^{(i)}b ^{(i)}+\mathrm{diag}(\Lambda^{(i)}\Delta^{(i)}),\] \[\underline{A}(\underline{x},\overline{x})=\Omega^{(0)},\qquad \underline{b}(\underline{x},\overline{x})=\sum\nolimits_{i=1}^{k}\Omega^{(i)}b ^{(i)}+\mathrm{diag}(\Omega^{(i)}\Theta^{(i)}) \tag{6}\]

where, for every \(i\in\{1,\ldots,k\}\), \(\Lambda^{(i)},\Delta^{(i)},\Omega^{(i)},\Theta^{(i)}\in\mathbb{R}^{n_{k} \times n_{i}}\) are as defined in (Zhang et al., 2018, Theorem 3.2) for the input perturbation set \(B_{\infty}(\frac{x+\overline{x}}{2},\frac{\overline{x}-\overline{x}}{2})\).
**Theorem 1** (Inclusion functions for Neural Networks): _Consider the \(k\)-layer feed-forward neural network \(u=\mathsf{N}(x)\) given by (2). Then,_

1. _for every input perturbation interval_ \([\underline{x},\overline{x}]\)_, the mapping_ \(\left[\frac{\mathsf{G}_{[\underline{x},\overline{x}]}}{\mathsf{G}_{[ \underline{x},\overline{x}]}}\right]:\mathcal{T}_{\geq 0}^{2n}\to\mathcal{T}_{\geq 0}^{2p}\) _defined by_ \[\frac{\mathsf{G}_{[\underline{x},\overline{x}]}(x,\widehat{x})}{ \overline{\mathsf{G}}_{[\underline{x},\overline{x}]}(x,\widehat{x})} =[\overline{A}(\underline{x},\overline{x})]^{+}x+[\overline{A}( \underline{x},\overline{x})]^{-}\widehat{x}+\underline{b}(\underline{x}, \overline{x}),\] (7) _is an inclusion function for the neural network_ \(\mathsf{N}\) _on_ \([\underline{x},\overline{x}]\)_;_
2. _the mapping_ \(\left[\frac{\mathsf{H}}{\mathsf{H}}\right]:\mathcal{T}_{\geq 0}^{2n}\to\mathcal{T}_{ \geq 0}^{2p}\) _defined by_ \(\underline{\mathsf{H}}(x,\widehat{x})=\mathsf{G}_{[\underline{x},\widehat{x} ]}(x,\widehat{x})\) _and_ \(\overline{\mathsf{H}}(x,\widehat{x})=\overline{\mathsf{G}}_{[\underline{x}, \widehat{x}]}(x,\widehat{x})\) _is an inclusion function for the neural network_ \(\mathsf{N}\) _on_ \(\mathbb{R}^{n}\)_._

**Proof**: Regarding part (i), suppose that \(\eta\leq\widehat{\eta}\in[\underline{x},\overline{x}]\) and \(z\in[\eta,\widehat{\eta}]\). Then,

\[\overline{\mathsf{G}}_{[\underline{x},\overline{y}]}(\eta, \widehat{\eta}) =[\overline{A}(\underline{x},\overline{x})]^{+}\widehat{\eta}+[ \overline{A}(\underline{x},\overline{x})]^{-}\eta+\overline{b}(\overline{x}, \overline{x})\geq[\overline{A}(\underline{x},\overline{x})]^{+}z+[\overline{A }(\underline{x},\overline{x})]^{-}z+\overline{b}(\overline{x},\overline{x})\] \[=[\overline{A}(\underline{x},\overline{x})]z+\overline{b}( \overline{x},\overline{x})\geq\mathsf{N}(z),\]

where the first inequality holds because \(z\in[\eta,\widehat{\eta}]\), the matrix \([\overline{A}(\underline{x},\overline{x})]^{+}\) is non-negative, and the matrix \([\overline{A}(\underline{x},\overline{x})]^{+}\) is non-positive. The second inequality holds by noting the fact that \([\overline{A}(\underline{x},\overline{x})]^{+}+[\overline{A}(\underline{x}, \overline{x})]^{-}=[\overline{A}(\underline{x},\overline{x})]\), for every \(\underline{x}\leq\overline{x}\). Finally, the last inequality holds by (5). Similarly, one can show that \(\underline{\mathsf{G}}_{[\underline{x},\overline{x}]}(\widehat{\eta},\eta) \leq\mathsf{N}(z)\). Regarding part (ii), suppose that \(x\leq\widehat{x}\), for every \(z\in[x,\widehat{x}]\), one can use a similar argument as part (i) to show that

\[\overline{\mathsf{G}}_{[\underline{x},\widehat{x}]}(x,\widehat{x}) =[\overline{A}(x,\widehat{x})]^{+}\widehat{x}+[\overline{A}(x, \widehat{x})]^{-}x+\overline{b}(x,\widehat{x})\geq[\overline{A}(x,\widehat{x} )]^{+}z+[\overline{A}(x,\widehat{x})]^{-}z+\overline{b}(x,\widehat{x})\] \[=[\overline{A}(x,\widehat{x})]z+\overline{b}(x,\widehat{x})\geq \mathsf{N}(z).\]

One can similarly show that \(\underline{\mathsf{G}}_{[\underline{x},\widehat{x}]}(\widehat{x},x)\leq \mathsf{N}(z)\), for every \(z\in[x,\widehat{x}]\). Moreover, suppose that \(\underline{x}=\overline{x}=z\). In this case, both inequalities in (4) are equality with \(\alpha_{L}^{(i)}=\alpha_{U}^{(i)}\) and \(\beta_{L}^{(i)}=\beta_{U}^{(i)}\). Thus, using the equation (6), we get that \(\underline{A}(z,z)=\overline{A}(z,z)\) and \(\underline{b}(z,z)=\overline{b}(z,z)\). This implies that \(\mathsf{N}(z)=\underline{A}(z,z)z+\underline{b}(z,z)=\underline{\mathsf{G}}_{ [z,z]}(z,z)\). Thus, \(\left[\frac{\mathsf{H}}{\mathsf{H}}\right]\) is an inclusion function for \(\mathsf{N}\) on the \(\mathbb{R}^{n}\). \(\blacksquare\)

**Remark 2**: _The following remarks are in order._

1. _For each layer_ \(i\in\{1,\ldots,k\}\)_, the intermediate pre-activation bounds_ \(U^{i},L^{i}\) _can be computed using either the Interval Bound Propagation (IBP)_ _(_Gowal et al._,_ 2018_)_ _or using the CROWN itself_ _(_Zhang et al._,_ 2018_)__. In the former case, the IBP can be considered as a special case of CROWN by choosing_ \(\alpha_{U}^{(i)}=\alpha_{L}^{(i)}=0\) _and_ \(\beta_{L}^{(i)}=L^{(i)}\) _and_ \(\beta_{U}^{(i)}=U^{(i)}\)_, for every_ \(i\in\{1,\ldots,k\}\)_._
2. _The inclusion function_ \(\left[\frac{\mathsf{G}_{[\underline{x},\overline{x}]}}{\mathsf{G}_{[ \underline{x},\overline{x}]}}\right]\) _defined in Theorem_ 1_(i) is linear but it depends on the input perturbation interval_ \([\underline{x},\overline{x}]\)_. On the other hand, the inclusion function_ \(\left[\frac{\mathsf{H}}{\mathsf{H}}\right]\) _defined in Theorem_ 1_(ii) is nonlinear but it is independent of the input perturbation set. For both inclusion functions, it can be shown that the smaller the input perturbation interval, the tighter the relaxation bounds in equation (_4_) are and, thus, the tighter the over- and the under-approximations of the neural network in equation (_5_) are._
## 5 Interval reachability of closed-loop system

In this section, we present a system-level approach for over-approximating the reachable set of the closed-loop system (3) with the neural network controller (2). The key idea is to design a suitable function \(d:\mathcal{T}^{2n}\times\mathcal{T}^{2q}\rightarrow\mathbb{R}^{n}\) and use it to construct the following embedding system associated to the closed-loop system (3):

\[\frac{d}{dt}\begin{bmatrix}x\\ \widehat{x}\end{bmatrix}=\begin{bmatrix}d(x,\widehat{x},w,\widehat{w})\\ d(\widehat{x},x,\widehat{w},w)\end{bmatrix} \tag{8}\]

Then, one can use the embedding system (8) to study the propagation of the state and disturbance bounds with time. Suppose that the initial condition set is given by \(\mathcal{X}_{0}\subseteq[\underline{x}_{0},\overline{x}_{0}]\) and the disturbance set is given by \(\mathcal{W}\subseteq[\underline{w},\overline{w}]\). Let \(\mathsf{F}\) be a decomposition function for the open-loop system (1) and \(\left[\frac{\mathsf{G}}{\mathsf{G}}\right]\) and \(\left[\frac{\mathsf{H}}{\mathsf{H}}\right]\) be the inclusion functions of the neural network (2) established in Theorem 1. We introduce "global" function \(d^{\mathrm{G}}\), "hybrid" function \(d^{\mathrm{H}}\), and "local" function \(d^{\mathrm{L}}\), for every \(i\in\{1,\ldots,n\}\), as follows:

\[d^{\mathrm{G}}_{i}(x,\widehat{x},w,\widehat{w}) =\begin{cases}\mathsf{F}_{i}(x,\widehat{x},\underline{\mathsf{H} }(x,\widehat{x}),\overline{\mathsf{H}}(x,\widehat{x}),w,\widehat{w}),&x\leq \widehat{x},w\leq\widehat{w},\\ \mathsf{F}_{i}(\widehat{x},x,\overline{\mathsf{H}}(x,\widehat{x}),\underline {\mathsf{H}}(x,\widehat{x}),\widehat{w},w),&\widehat{x}\leq x,\widehat{w}\leq w,\end{cases}\] \[d^{\mathrm{H}}_{i}(x,\widehat{x},w,\widehat{w}) =\begin{cases}\mathsf{F}_{i}(x,\widehat{x},\underline{\mathsf{G} }_{[x,\widehat{x}]}(x,\widehat{x}_{[i:x]}),\overline{\mathsf{G}}_{[x, \widehat{x}]}(x,\widehat{x}_{[i:x]}),w,\widehat{w}),&x\leq\widehat{x},w\leq \widehat{w},\\ \mathsf{F}_{i}(\widehat{x},x,\overline{\mathsf{G}}_{[x,\widehat{x}]}(x_{[i: \widehat{x}]}),\overline{\mathsf{G}}_{[x,\widehat{x}]}(x_{[i:\widehat{x}]}), \widehat{w},w),&\widehat{x}\leq x,\widehat{w}\leq w,\end{cases}\] \[d^{\mathrm{L}}_{i}(x,\widehat{x},w,\widehat{w}) =\begin{cases}\mathsf{F}_{i}(x,\widehat{x},\underline{\mathsf{H} }(x,\widehat{x}_{[i:x]}),\overline{\mathsf{H}}(x,\widehat{x}_{[i:x]}),w, \widehat{w})&x\leq\widehat{x},w\leq\widehat{w}\\ \mathsf{F}_{i}(\widehat{x},x,\overline{\mathsf{H}}(x_{[i:\widehat{x}]}, \widehat{x}),\underline{\mathsf{H}}(x_{[i:\widehat{x}]},\widehat{x}),\widehat {w},w)&\widehat{x}\leq x,\widehat{w}\leq w.\end{cases}\]

For every \(s\in\{\mathrm{G},\mathrm{L},\mathrm{H}\}\), the trajectory of the embedding system (8) with \(d=d^{s}\) and disturbance \([\begin{smallmatrix}w\\ \widehat{w}\end{smallmatrix}]=[\frac{w}{\overline{w}}]\) starting from \(\left[\frac{\underline{x}_{0}}{\overline{x}_{0}}\right]\) is denoted by \(t\mapsto\left[\frac{\underline{x}^{s}(t)}{\overline{x}^{s}(t)}\right]\).

**Theorem 3** (Interval over-approximation of reachable sets): _Consider the control system (1) with the neural network controller (2). Suppose that the initial condition \(x_{0}\) belongs to \(\mathcal{X}_{0}\subseteq[\underline{x}_{0},\overline{x}_{0}]\) and the disturbance \(w\) belongs to \(\mathcal{W}\subseteq[\underline{w},\overline{w}]\). Let \(\mathsf{F}:\mathcal{T}^{2n}\times\mathcal{T}^{2p}\times\mathcal{T}^{2q} \rightarrow\mathbb{R}^{n}\) be a decomposition function for the open-loop system (1). Then the following statement holds:_

1. _for every_ \(s\in\{\mathrm{G},\mathrm{L},\mathrm{H}\}\) _and every_ \(t\in\mathbb{R}_{\geq 0}\)_, we have_ \[\mathcal{R}(t,\mathcal{X}_{0},\mathcal{W})\subseteq[\underline{x}^{s}(t), \overline{x}^{s}(t)].\]

_If all the activation functions of the neural network \(u=\mathsf{N}(x)\) are ReLU, then_

1. _for every_ \(t\in\mathbb{R}_{\geq 0}\)_, we have_ \[\mathcal{R}(t,\mathcal{X}_{0},\mathcal{W})\subseteq[\underline{x}^{\mathrm{L} }(t),\overline{x}^{\mathrm{L}}(t)]\subseteq[\underline{x}^{\mathrm{H}}(t), \overline{x}^{\mathrm{H}}(t)]\subseteq[\underline{x}^{\mathrm{G}}(t),\overline{ x}^{\mathrm{G}}(t)].\]

* Regarding part (i), we provide the proof for \(s=\mathrm{H}\). The proof of other cases are mutadis mutandis with \(s=\mathrm{H}\). Note that, by (Abate et al., 2021, Theorem 1), the tight decomposition
function \(d^{c}\) for the closed-loop system \(f^{\text{cl}}\) and the tight decomposition function \(d^{o}\) for the open loop system \(f\) can be computed as follows:

\[d^{c}_{i}(x,\widehat{x},w,\widehat{w}) =\begin{cases}\min_{z\in[x,\widehat{x}],\in[w,\widehat{\omega}]}f_ {i}(z,\mathsf{N}(z),\xi),&x\leq\widehat{x},w\leq\widehat{w}\\ \max_{z\in[\widehat{x},x],\in[w,\widehat{\omega}]}f_{i}(z,\mathsf{N}(z),\xi),&\widehat{x}\leq x,\widehat{w}\leq w.\end{cases} \tag{9}\] \[d^{o}_{i}(x,\widehat{x},u,\widehat{u},w,\widehat{w}) =\begin{cases}\min_{z\in[x,\widehat{x}],\in[w,\widehat{\omega}]}f_ {i}(z,\eta,\xi),&x\leq\widehat{x},u\leq\widehat{u},w\leq\widehat{w},\\ \max_{z\in[\widehat{x},x],\in[w,\widehat{\omega}]}f_{i}(z,\eta,\xi),&\widehat {x}\leq x,\widehat{u}\leq u,\widehat{w}\leq w.\end{cases} \tag{10}\]

The trajectory of the embedding system (8) with \(d=d^{c}\) and disturbance \([\frac{w}{\widehat{w}}]=[\frac{w}{\widehat{w}}]\) starting from \(\left[\frac{\overline{x}_{0}}{\overline{x}_{0}}\right]\) is denoted by \(t\mapsto\left[\frac{x^{c}(t)}{\overline{x}^{c}(t)}\right]\). Since \(\mathsf{F}\) is a decomposition function for the open-loop system \(f\), for every \(x\leq\widehat{x}\), \(u\leq\widehat{u}\), \(w\leq\widehat{w}\), and \(i\in\{1,\ldots,n\}\), we get

\[\mathsf{F}_{i}(x,\widehat{x},u,\widehat{u},w,\widehat{w})\leq d^{o}_{i}(x, \widehat{x},u,\widehat{u},w,\widehat{w})=\min_{z\in[x,\widehat{x}],z_{i}=x_{ i}\atop\xi\in[w,\widehat{\omega}],\eta\in[u,\widehat{u}]}f_{i}(z,\eta,\xi). \tag{11}\]

Suppose that \(x\leq\widehat{x}\) and let \(i\in\{1,\ldots,n\}\) and \(y\in[x,\widehat{x}]\) be such that \(y_{i}=x_{i}\). By Theorem 1(i), \(\left[\frac{\underline{\mathsf{G}}_{[x,\widehat{x}]}}{\overline{\mathsf{G}}_{[ x,\widehat{x}]}}\right]\) is an inclusion function for \(\mathsf{N}\) on \([x,\widehat{x}]\). Moreover, \(y\in[x,\widehat{x}_{[i:x]}]\subseteq[x,\widehat{x}]\) and thus

\[\underline{\mathsf{G}}_{[x,\widehat{x}]}(x,\widehat{x}_{[i:x]})\leq\mathsf{N} (y)\leq\overline{\mathsf{G}}_{[x,\widehat{x}]}(x,\widehat{x}_{[i:x]}). \tag{12}\]

Therefore, for every \(i\in\{1,\ldots,n\}\), we get

\[d^{c}_{i}(x,\widehat{x},w,\widehat{w}) =\min_{z\in[x,\widehat{x}],\in[w,\widehat{\omega}],z_{i}=x_{i} \atop z_{i}=x_{i}}f_{i}(z,\mathsf{N}(z),\xi)\geq\min_{z\in[x,\widehat{x}], \in[w,\widehat{\omega}],z_{i}=x_{i}\atop u\in[\underline{\mathsf{G}}_{[x, \widehat{x}]}(x,\widehat{x}_{[i:x]}),\overline{\mathsf{G}}_{[x,\widehat{x}]}( x,\widehat{x}_{[i:x]})]}f_{i}(z,u,\xi)\] \[\geq\mathsf{F}_{i}(x,\widehat{x},\underline{\mathsf{G}}_{[x, \widehat{x}]}(x,\widehat{x}_{[i:x]}),\overline{\mathsf{G}}_{[x,\widehat{x}]}( x,\widehat{x}_{[i:x]}),w,\widehat{w})=d^{\mathrm{H}}_{i}(x,\widehat{x},w, \widehat{w}). \tag{13}\]

where the first equality holds by definition of \(d^{c}\), the second inequality holds by equation (12), the third inequality holds by equation (11), and the fourth inequality holds by definition of \(d^{\mathrm{H}}\). Similarly, one can show that \(d^{c}_{i}(\widehat{x},x,\widehat{w},w)\leq d^{\mathrm{H}}_{i}(\widehat{x},x, \widehat{w},w)\), for every \(i\in\{1,\ldots,n\}\). This implies that \(\left[\begin{smallmatrix}d^{\mathrm{H}}(x,\widehat{x},w,\widehat{w})\\ d^{\mathrm{H}}(\widehat{x},x,\widehat{w},w)\end{smallmatrix}\right]\leq_{ \mathrm{SE}}\left[\begin{smallmatrix}d^{c}(x,\widehat{x},w,\widehat{w})\\ d^{c}(\widehat{x},x,\widehat{w},w)\end{smallmatrix}\right]\), for every \(x\leq\widehat{x}\) and every \(w\leq\widehat{w}\). Note that, by (Abate et al., 2021, Theorem 1), the vector field \(\left[\begin{smallmatrix}d^{c}(x,\widehat{x},w,\widehat{w})\\ d^{c}(\widehat{x},x,\widehat{w},w)\end{smallmatrix}\right]\) is monotone with respect to the southeast order \(\leq_{\mathrm{SE}}\) on \(\mathbb{R}^{2n}\). Now, we can use (Michel et al., 2008, Theorem 3.8.1), to deduce that \(\left[\frac{x^{\mathrm{H}}(t)}{\overline{x}^{\mathrm{H}}(t)}\right]\leq\left[ \frac{x^{c}(t)}{\overline{x}^{c}(t)}\right]\), for every \(t\in\mathbb{R}_{\geq 0}\). This implies that \([\underline{x}^{c}(t),\overline{x}^{c}(t)]\subseteq[\underline{x}^{\mathrm{H} }(t),\overline{x}^{\mathrm{H}}(t)]\). On the other hand, by (Abate et al., 2021, Theorem 2), we know that \(\mathcal{R}(t,\mathcal{X}_{0},\mathcal{W})\subseteq[\underline{x}^{c}(t), \overline{x}^{c}(t)]\), for every \(t\in\mathbb{R}_{\geq 0}\). This lead to \(\mathcal{R}(t,\mathcal{X}_{0},\mathcal{W})\subseteq[\underline{x}^{\mathrm{H} }(t),\overline{x}^{\mathrm{H}}(t)]\), for every \(t\in\mathbb{R}_{\geq 0}\).

Regarding part (ii), suppose that all the activation functions are ReLU and let \(L^{(i)}_{[x,\widehat{x}]}\) and \(U^{(i)}_{[x,\widehat{x}]}\) are the intermediate bounds associated to the input perturbation \([x,\widehat{x}]\). Then, for every \(\left[\frac{x}{\widehat{x}}\right]\leq_{\mathrm{SE}}\left[\frac{y}{\widehat{y} }\right]\in\mathcal{T}_{\geq 0}^{2n}\), the intermediate bounds in CROWN satisfy:

\[L^{(i)}_{[x,\widehat{x}]}\leq L^{(i)}_{[y,\widehat{y}]},\qquad U^{(i)}_{[y, \widehat{y}]}\leq U^{(i)}_{[x,\widehat{x}]},\qquad\text{ for every }i\in\{1,\ldots,k\}.\]
Since all the activation functions are ReLU, we can use (Zhang et al., 2018, Table 1 and Table 2) to show that

\[\alpha^{(i)}_{L,[x,\widehat{x}]} \leq\alpha^{(i)}_{L,[y,\widehat{y}]},\qquad\alpha^{(i)}_{U,[y, \widehat{y}]}\leq\alpha^{(i)}_{U,[x,\widehat{x}]},\] \[\beta^{(i)}_{L,[x,\widehat{x}]} \leq\beta^{(i)}_{L,[y,\widehat{y}]},\qquad\beta^{(i)}_{U,[y, \widehat{y}]}\leq\beta^{(i)}_{U,[x,\widehat{x}]},\qquad\text{ for every }i\in\{1,\ldots,k\}\]

Therefore, using the equations (4) and (6), and the formula in (Zhang et al., 2018, Theorem 3.2) for \(\Lambda^{(i)},\Delta^{(i)},\Omega^{(i)},\Theta^{(i)}\), we get

\[\underline{A}(x,\widehat{x}) \leq\underline{A}(y,\widehat{y}),\qquad\overline{A}(y,\widehat{y })\leq\overline{A}(x,\widehat{x}),\] \[\underline{b}(x,\widehat{x}) \leq\underline{b}(y,\widehat{y}),\qquad\overline{b}(y,\widehat{y })\leq\overline{b}(x,\widehat{x}),\]

This implies that

\[\underline{\mathsf{G}}_{[x,\widehat{x}]}(w,\widehat{w}) =[\underline{A}(x,\widehat{x})]^{+}w+[\underline{A}(x,\widehat{ x})]^{-}\widehat{w}\leq[\underline{A}(y,\widehat{y})]^{+}w+[\underline{A}(y, \widehat{y})]^{-}\widehat{w}\] \[\leq[\underline{A}(y,\widehat{y})]^{+}v+[\underline{A}(y, \widehat{y})]^{-}\widehat{v}=\underline{\mathsf{G}}_{[y,\widehat{y}]}(v, \widehat{v}).\]

Similarly, we can show that \(\overline{\mathsf{G}}_{[y,\widehat{y}]}(v,\widehat{v})\leq\overline{\mathsf{ G}}_{[x,\widehat{x}]}(w,\widehat{w})\). As a result, for every \([\frac{x}{x}]\leq_{\mathrm{SE}}\left[\frac{y}{y}\right]\in\mathcal{T}_{\geq 0}^{2n}\) and \([\frac{w}{\widehat{w}}]\leq_{\mathrm{SE}}\left[\frac{v}{v}\right]\in\mathcal{T }_{\geq 0}^{2n}\),

\[\left[\frac{\underline{\mathsf{G}}_{[x,\widehat{x}]}(w,\widehat{w})}{ \overline{\mathsf{G}}_{[x,\widehat{x}]}(w,\widehat{w})}\right]\leq_{\mathrm{ SE}}\left[\frac{\underline{\mathsf{G}}_{[y,\widehat{y}]}(v,\widehat{v})}{ \overline{\mathsf{G}}_{[y,\widehat{y}]}(v,\widehat{v})}\right] \tag{14}\]

Note that, using property (14) and definition of decomposition function, one can easily check that \(d^{\mathrm{G}},d^{\mathrm{L}},d^{\mathrm{H}}\) are decomposition functions for \(f^{\mathrm{cl}}\) and thus, they are monotone with respect to the southeast order \(\leq_{\mathrm{SE}}\) on \(\mathbb{R}^{2n}\). On the other hand, we have \(\widehat{x}_{[i:x]}\leq\widehat{x}\). Therefore, using (14),

\[\left[\frac{\underline{\mathsf{G}}_{[x,\widehat{x}]}(x,\widehat{x}_{[i:x]})}{ \overline{\mathsf{G}}_{[x,\widehat{x}]}(x,\widehat{x}_{[i:x]})}\right]\leq_{ \mathrm{SE}}\left[\frac{\underline{\mathsf{G}}_{[x,\widehat{x}]}(x,\widehat{x} )}{\overline{\mathsf{G}}_{[x,\widehat{x}]}(x,\widehat{x})}\right],\qquad \text{ for every }x\leq\widehat{x} \tag{15}\]

This implies that, for every \(x\leq\widehat{x}\) and every \(w\leq\widehat{w}\),

\[d^{\mathrm{H}}(x,\widehat{x},w,\widehat{w}) =\mathsf{F}_{i}(x,\widehat{x},\underline{\mathsf{G}}_{[x, \widehat{x}]}(x,\widehat{x}_{[i:x]}),\overline{\mathsf{G}}_{[x,\widehat{x}]}( x,\widehat{x}_{[i:x]}),w,\widehat{w})\] \[\leq\mathsf{F}_{i}(x,\widehat{x},\underline{\mathsf{G}}_{[x, \widehat{x}]}(x,\widehat{x}),\overline{\mathsf{G}}_{[x,\widehat{x}]}(x, \widehat{x}),w,\widehat{w})=d^{\mathrm{G}}(x,\widehat{x},w,\widehat{w})\]

where the inequality holds using equation 15 and the fact that \(\mathsf{F}\) is a decomposition function for \(f^{\mathrm{cl}}\). Similarly, we can show that \(d^{\mathrm{H}}(\widehat{x},x,\widehat{w},w)\leq d^{\mathrm{G}}(\widehat{x},x,\widehat{w},w)\), for every \(x\leq\widehat{x}\) and every \(w\leq\widehat{w}\). This implies that \(\left[\frac{d^{\mathrm{G}}(x,\widehat{x},w,\widehat{w})}{d^{\mathrm{G}}(x,x,\widehat{w},w)}\right]\leq_{\mathrm{SE}}\left[\frac{d^{\mathrm{H}}(x,\widehat{x },w,\widehat{w})}{d^{\mathrm{H}}(\widehat{x},x,\widehat{w},w)}\right]\), for every \(x\leq\widehat{x}\) and every \(w\leq\widehat{w}\). Note that \(d^{\mathrm{G}}\) is a decomposition function for \(f^{\mathrm{cl}}\) and thus the vector field \(\left[\frac{d^{\mathrm{G}}(x,\widehat{x},w,\widehat{w})}{d^{\mathrm{G}}(x, \widehat{x},\widehat{w},w)}\right]\) is monotone with respect to the southeast order \(\leq_{\mathrm{SE}}\) on \(\mathbb{R}^{2n}\). Now, we can use (Michel et al., 2008, Theorem 3.8.1), to deduce that \(\left[\frac{x^{\mathrm{G}}(t)}{x^{\mathrm{G}}(t)}\right]\leq\left[\frac{x^{ \mathrm{H}}(t)}{x^{\mathrm{H}}(t)}\right]\), for every \(t\in\mathbb{R}_{\geq 0}\). This implies that \([\underline{x}^{\mathrm{H}}(t),\overline{x}^{\mathrm{H}}(t)]\subseteq[ \underline{x}^{\mathrm{G}}(t),\overline{x}^{\mathrm{G}}(t)]\). The proofs for the other inclusions are similar and we remove it for the sake of brevity. 

For the special case when the open-loop system (1) is linear in state and affine in control, one can find closed-from expressions for \(d^{\mathrm{G}}\), \(d^{\mathrm{H}}\), and \(d^{\mathrm{L}}\).
**Corollary 4** (Linear systems): _Consider the control system (1) with \(f(x,u,w)=Ax+Bu+Cw\) with \(A\in\mathbb{R}^{n\times n}\), \(B\in\mathbb{R}^{n\times p}\), and \(C\in\mathbb{R}^{n\times q}\) and with the neural network controller (2). Then we recover all the results of Theorem (3) with the following "global", "hybrid", and "local" functions:_

\[d^{\rm G}(x,\widehat{x},w,\widehat{w}) =\left(\lceil A\rceil^{\rm Mzl}+\left[R(x,\widehat{x})\right]^{+} \right)x+\left(\lfloor A\rfloor^{\rm Mzl}+\left[S(x,\widehat{x})\right]^{-} \right)\widehat{x}+C^{+}w+C^{-}\widehat{w}\] \[d^{\rm H}(x,\widehat{x},w,\widehat{w}) =\lceil A+R(x,\widehat{x})\rceil^{\rm Mzl}x+\lfloor A+S(x, \widehat{x})\rceil^{\rm Mzl}\widehat{x}+C^{+}w+C^{-}\widehat{w}\] \[d^{\rm L}_{i}(x,\widehat{x},w,\widehat{w}) =\lceil A+R(x,\widehat{x}_{[i:x]})\rceil^{\rm Mzl}x+\lfloor A+S(x,\widehat{x}_{[i:x]})\rfloor^{\rm Mzl}\widehat{x}+C^{+}w+C^{-}\widehat{w}, \quad\forall i\in\{1,\ldots,n\},\]

_where \(R(x,\widehat{x})=B^{+}[\underline{A}(x,\widehat{x})]+B^{-}[\overline{A}(x, \widehat{x})]\) and \(S(x,\widehat{x})=B^{+}[\overline{A}(x,\widehat{x})]+B^{-}[\underline{A}(x, \widehat{x})]\)._

**Proof** Note that a decomposition function for the open-loop system is given by \(\mathsf{F}(x,\widehat{x},u,\widehat{u},w,\widehat{w})=\lceil A\rceil^{\rm Mzl }x+\lfloor A\rfloor^{\rm Mzl}\widehat{x}+B^{+}u+B^{-}\widehat{u}+C^{+}w+C^{-} \widehat{w}\)(Coogan, 2020, Example 3). Using simple algebraic manipulations, one can show that for every \(i\in\{1,\ldots,n\}\) and every \(M\in\mathbb{R}^{n\times n}\),

\[\left(M^{+}x+M^{-}\widehat{x}_{[i:x]}\right)_{i}=\left(\lceil M\rceil^{\rm Mzl }x+\lfloor M\rfloor^{\rm Mzl}\widehat{x}\right)_{i}. \tag{16}\]

Thus, using Theorem 3 and the identity (16), for every \(i\in\{1,\ldots,n\}\), we have

\[d^{\rm H}_{i}(x,\widehat{x},w,\widehat{w}) =\left(\lceil A\rceil^{\rm Mzl}x+\lfloor A\rfloor^{\rm Mzl} \widehat{x}+B^{+}[\underline{A}]^{+}(x,\widehat{x})x+B^{+}[\underline{A}]^{- }(x,\widehat{x})\widehat{x}_{[i:x]}\right)_{i}\] \[\quad\quad+\left(B^{-}[\overline{A}]^{+}(x,\widehat{x})x+B^{-}[ \overline{A}]^{-}(x,\widehat{x})\widehat{x}_{[i:x]}+C^{+}w+C^{-}\widehat{w} \right)_{i}\] \[=\left(\lceil A\rceil^{\rm Mzl}x+\lfloor A\rfloor^{\rm Mzl} \widehat{x}+\lceil B^{+}\underline{A}(x,\widehat{x})\rceil^{\rm Mzl}x+ \lfloor B^{+}\underline{A}(x,\widehat{x})\rceil^{\rm Mzl}\widehat{x}\right)_ {i}\] \[\quad\quad+\left(\lfloor B^{-}\overline{A}(x,\widehat{x})\rfloor^ {\rm Mzl}x+\lceil B^{-}\overline{A}(x,\widehat{x})\rceil^{\rm Mzl}\widehat{x}+ C^{+}w+C^{-}\widehat{w}\right)_{i}\] \[=\left(\lceil A+B^{+}\underline{A}(x,\widehat{x})+B^{-}\overline {A}(x,\widehat{x})\rceil^{\rm Mzl}+\lfloor A+B^{+}[\overline{A}(x,\widehat{ x})]+B^{-}[\underline{A}(x,\widehat{x})]\rfloor^{\rm Mzl}\right)_{i}\] \[\quad\quad+\left(C^{+}w+C^{-}\widehat{w}\right)_{i}\] \[=\left(\lceil A+R(x,\widehat{x})\rceil^{\rm Mzl}x+\lfloor A+S(x, \widehat{x})\rfloor^{\rm Mzl}\widehat{x}+C^{+}w+C^{-}\widehat{w}\right)_{i}.\]

This completes the proof for \(s={\rm H}\). The proofs of the other cases follow similarly.

**Remark 5**: _The following remarks are in order._

1. _(Tight decomposition function) One can use the tight decomposition function of the closed loop system (_3_) in the embedding system (_8_) to obtain hyper-rectangular over-approximations for reachable set of the system. However, finding the tight decomposition function requires solving a computationally complicated optimization problem and generally difficult. In this context, Theorem_ 3 _uses the decomposition function of the open-loop system (_1_) and the inclusion function of the neural network to establish three under-approximation the tight decomposition function of the closed-loop system (_8_)._
2. _(Comparison with the literature) For linear systems, Corollary_ 4 _can be used to show that the forward Euler integration of the embedding system (_8_) with_ \(d=d^{\rm H}\) _and with a small enough time-step will lead to the identical over-approximation sets as_ _(Everett et al._,_ 2021b_, Lemma IV.3_) applied to the forward Euler discretization of the linear system._
_
3. _(Generality of the approach) Theorem_ 3 _proposes a general embedding-based framework for verification of the closed-loop system (_3_) which is based on combining mixed monotone reachability of the open-loop system with interval analysis of the neural network. Using this perspective, our approach can be applied to arbitrary neural network verification algorithms as long as one can construct an inclusion function for the neural network. This is in contrast with most of the existing approaches for neural network closed-loop reachability analysis, which are heavily dependent on the neural network verification algorithm (see for instance_ (Everett et al., 2021)_) and_ _(Sidrane et al., 2022)__)._
4. _(Computational complexity)_ _From a computational perspective, the framework presented in Theorem_ 3 _consists of two main ingredients: (i) evaluating CROWN to compute the inclusion function of the neural network as in Theorem_ 1_, and (ii) integrating the embedding dynamical system (_8_). For a neural network with_ \(k\)_-layer and_ \(N\) _neuron per layer, the complexity of CROWN is_ \(\mathcal{O}(k^{2}N^{3})\)__(Zhang et al., 2018)_. Moreover, the functions_ \(d^{\mathrm{G}}\) _and_ \(d^{\mathrm{H}}\) _call CROWN once per integration step, while the function_ \(d^{\mathrm{L}}\) _calls CROWN_ \(n\) _times per integration step. The run time of the integration process depends on the form of the open-loop decomposition function_ \(\mathsf{F}\)_._

## 6 Efficient reachability analysis via partitioning

In this section, we develop a suitable partitioning of the state space and combine it with Theorem (3) to obtain a computationally efficient algorithm for generating reachable set over-approximations of the closed-loop system. Interval methods are known to suffer from large over-approximation error due to the wrapping effect (Jaulin et al., 2001, Section 2.2.4). However, their ease of computation allows for combining them efficiently with partitioning strategies. Our partitioning strategy consists of two main components: (i) a uniform division of the state space of the embedding system (8) to compute the neural network inclusion functions using CROWN, and (ii) a uniform subdivision to implement the integration on the embedding system (8). We first pick an \(s\in\{\mathrm{G},\mathrm{H},\mathrm{L}\}\) and start with the initial perturbation set \(\mathcal{X}_{0}\). In the first step, we find the smallest hyper-rectangle containing \(\mathcal{X}_{0}\) by computing \((\overline{x}_{0})_{j}=\max_{x\in\mathcal{X}_{0}}x_{j}\) and \((\underline{x}_{0})_{j}=\min_{x\in\mathcal{X}_{0}}x_{j}\), for every \(i\in\{1,\ldots,n\}\). We then divide the hyper-rectangle \([\underline{x}_{0},\overline{x}_{0}]\) into \(D_{a}\) partitions and obtain the set \(\{[\underline{x}^{1},\overline{x}^{1}],\ldots,[\underline{x}^{D_{a}}, \overline{x}^{D_{a}}]\}\). For every \(k\in\{1,\ldots,D_{a}\}\), we compute the trajectory of the embedding system (8) with \(d=d^{\mathrm{s}}\) and the initial condition \(\left\lfloor\frac{\underline{x}^{k}}{\underline{x}^{k}}\right\rfloor\) at time \(\Delta t\) as in Theorem 3. For every \(k\in\{1,\ldots,D_{a}\}\), we divide each states of the hyper-rectangle \([\underline{x}^{1},\overline{x}^{1}]\) into \(D_{s}\) partitions, obtaining the subpartitions \(\{[\underline{x}^{k,1},\overline{x}^{k,1}],\ldots,[\underline{x}^{k,D_{s}}, \overline{x}^{k,D_{s}}]\}\). For every \(k\in\{1,\ldots,D_{a}\}\) and \(l\in\{1,\ldots,D_{s}\}\), we compute the trajectory of the embedding system (8) with \(d=d^{\mathrm{Bs}}\), where, for every \(i\in\{1,\ldots,n\}\),

\[d^{\mathrm{Bs}}_{i}(x,\widehat{x},w,\widehat{w})=\mathsf{F}_{i}(x,\widehat{x},\underline{\mathsf{G}}_{[\underline{x}^{k},\overline{x}^{k}]}(x,\widehat{x} _{[i:x]}),\overline{\mathsf{G}}_{[\underline{x}^{k},\overline{x}^{k}]}( \widehat{x}_{[i:x]},x),w,\widehat{w}) \tag{17}\]

and the initial condition \(\left\lfloor\frac{\underline{x}^{k,l}}{\underline{x}^{k,l}}\right\rfloor\) at time \(\Delta t\) as in Theorem 3. We then set \(\mathcal{X}_{1}=\bigcup_{k=1}^{D_{a}}\bigcup_{l=1}^{D_{s}}[\underline{x}^{l,k} (\Delta t),\overline{x}^{l,k}(\Delta t)]\) and repeat this procedure on the initial set \(\mathcal{X}_{1}\). We keep repeating this algorithm until we get to the final time \(T\). Note that our partitioning approach is different from (Everett et al., 2021) and (Xiang et al., 2021) in that we re-partitioning the state-space at every time step. A summary of the above procedure is presented in Algorithm 1.
```
1:\(s\in\{\mathrm{G},\mathrm{H},\mathrm{L}\}\), the initial set \(\mathcal{X}_{0}\), the final time \(T\), the actuation step \(\Delta t\), divisions \(D_{a}\), subdivision \(D_{s}\) Output: the over-approximation of the reachable set \(\overline{\mathcal{R}}(T,0,\mathcal{X}_{0},\mathcal{W})\)
2:\(j\gets 0\)
3:while\(j<\lfloor\frac{T}{\Delta t}\rfloor\)do
4:\(\underline{x}_{i}\leftarrow\min_{x\in\mathcal{X}_{j}}x_{i}\), for every \(i\in\{1,\ldots,n\}\)
5:\(\overline{x}_{i}\leftarrow\max_{x\in\mathcal{X}_{j}}x_{i}\), for every \(i\in\{1,\ldots,n\}\)
6:\(\{[\underline{x}^{1},\overline{x}^{1}],\ldots,[\underline{x}^{D_{a}}, \overline{x}^{D_{a}}]\}\leftarrow\) uniform_partition\(([\underline{x},\overline{x}],D_{a})\)
7:for\(k=\{1,\ldots,D_{a}\}\)do
8: Compute \(\underline{A}(\underline{x}^{k},\overline{x}^{k})\), \(\overline{A}(\underline{x}^{k},\overline{x}^{i})\), \(\underline{b}(\underline{x}^{k},\overline{x}^{k})\), and \(\overline{b}(\underline{x}^{k},\overline{x}^{k})\) using CROWN and (6).
9:\(\{[\underline{x}^{k,1},\overline{x}^{k,1}],\ldots,[\underline{x}^{k,D_{s}}, \overline{x}^{k,D_{s}}]\}\leftarrow\) uniform_partition\(([\underline{x}^{k},\overline{x}^{k}],D_{s})\)
10:for\(l=\{1,\ldots,D_{s}\}\)do
11: Compute \(\left[\frac{\underline{x}^{k,l}(\Delta t)}{\overline{x}^{k,l}(\Delta t)}\right]\) for system (8) with \(d=d^{\mathrm{Bs}}\) and initial condition \(\left[\frac{\underline{x}^{k,l}}{\overline{x}^{k,l}}\right]\).
12:endfor
13:endfor
14:\(\mathcal{X}_{j+1}=\bigcup_{k=1}^{D_{a}}\bigcup_{l=1}^{D_{s}}[\underline{x}^{l,k}(\Delta t),\overline{x}^{l,k}(\Delta t)]\)
15:endwhile
16:return\(\overline{\mathcal{R}}(T,0,\mathcal{X}_{0},\mathcal{W})\leftarrow\mathcal{X} _{\lfloor\frac{T}{\Delta t}\rfloor}\)
```

**Algorithm 1** Over-approximation of reachable sets of (3)

## 7 Numerical Simulations

In this section, we show the efficiency of our reachability analysis using numerical experiments on a nonlinear vehicle model and a linear quadrotor model1.

Footnote 1: All the code is available at

[https://github.com/gtfactslab/L4DC2023_NNControllerReachability](https://github.com/gtfactslab/L4DC2023_NNControllerReachability)

### Nonlinear Vehicle model

We consider the dynamics of a vehicle adopted from (Polack et al., 2017) satisfying the following nonlinear ordinary differential equation:

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

where \([p_{x},p_{y}]^{\top}\in\mathbb{R}^{2}\) is the displacement of the center of mass, \(\phi\in[-\pi,\pi)\) is the heading angle in the plane, \(v\in\mathbb{R}^{+}\) is the speed of the center of mass. Control input \(u_{1}\) is the applied force subject to disturbance \(w\), input \(u_{2}\) is the angle of the front wheels, and \(\beta(u_{2})=\arctan\left(\frac{\ell_{f}}{\ell_{f}+\ell_{r}}\tan(u_{2})\right)\) is the slip slide angle. We set \(x=[p_{x},p_{y},\phi,v]^{\top}\) and \(u=[u_{1},u_{2}]^{\top}\). We use the following tight decomposition function for the open-loop system:

\[\mathsf{F}(x,\widehat{x},u,\widehat{u},w,\widehat{w})=\begin{bmatrix}d^{b_{1}b_ {2}}\left([v,d^{\cos}(\phi+\beta(u_{2}),\widehat{\phi}+\beta(\widehat{u}_{2})) ]^{\top},[\widehat{v},d^{\cos}(\widehat{\phi}+\beta(\widehat{u}_{2}),\phi+ \beta(u_{2}))]^{\top}\right)\\ d^{b_{1}b_{2}}\left([v,d^{\sin}(\phi+\beta(u_{2}),\widehat{\phi}+\beta( \widehat{u}_{2}))]^{\top},[\widehat{v},d^{\sin}(\widehat{\phi}+\beta(\widehat{u }_{2}),\phi+\beta(u_{2}))]^{\top}\right)\\ d^{b_{1}b_{2}}\left([v,d^{\sin}(\beta(u_{2}),\beta(\widehat{u}_{2}))]^{\top},[ \widehat{v},d^{\sin}(\beta(\widehat{u}_{2}),\beta(u_{2}))]^{\top}\right)\\ u_{1}+w\end{bmatrix},\]
where \(d^{b_{1}b_{2}}\), \(d^{\text{cos}}\), and \(d^{\text{sin}}\) are defined in (Cao et al., 2022).

We designed an offline nonlinear model predictive controller in Python using Casadi (Andersson et al., 2019) to steer the vehicle to the origin while avoiding obstacles. We use a fixed horizon of 20 with an actuation step of 0.25 seconds, and a quadratic cost function with \(Q=\operatorname{diag}(1,1,0,0)\), \(Q_{hor}=\operatorname{diag}(100,100,0,1)\), and other regularizing terms tuned empirically. Additionally, we add circular obstacles with 25% padding as hard constraints with slack variables; in Figures 1 and 2, we consider one centered at \((4,4)\) with a radius of 2.4. We simulated 65000 real trajectories (5s, 20 control actions) with initial conditions uniformly sampled from a specified region, and aggregated the data into a set of \(1.3M\) training pairs \((x,u)\in\mathbb{R}^{4}\times\mathbb{R}^{2}\). A neural network \(u=\mathsf{N}(x)\) with \(2\) hidden layers with \(100\) neurons per layer and ReLU activation was trained in Pytorch to approximate the model predictive controller under a scaled Mean Squared Error loss. We use Algorithm 1 with \(D_{a}=16\) and \(D_{s}=1\) to provide over-approximations for the reachable sets of the vehicle model (18), comparing functions \(d^{\text{G}}\), \(d^{\text{L}}\), and \(d^{\text{H}}\) from Theorem 3. Algorithm 1 line 7 is computed using auto_LiRPA (Xu et al., 2020). The results are shown in Figure 1.

### Linear 6D quadrotor model

We use the linear 6D quadrotor model adopted from (Ivanov et al., 2019) with the dynamics \(\dot{x}=Ax+Bu+c\) where \(A,B,c,u\) are as define in (Everett et al., 2021) and \(x=[p_{x},p_{y},p_{z},v_{x},x_{y},v_{z}]^{\top}\) is the state of the system with \(p_{x},p_{y}\), and \(p_{y}\) being the linear displacement in the \(x,y\), and \(z\) directions, respectively and \(v_{x},v_{y}\), and \(v_{z}\) the velocity in the \(x,y\), and \(z\) directions respectively. The neural network controller \(u=\mathsf{N}(x)\) consists of \(2\) layers with \(32\) neurons in each layer and is identical to (Everett et al., 2021, Section H). One can compute the following tight decomposition function for the open-loop system \(\mathsf{F}(x,\widehat{x},u,\widehat{u},w,\widehat{w})=Ax+B^{+}u+B^{-}\widehat {u}+c\). We compare the

Figure 1: Performance of the three different functions \(d^{\text{G}}\), \(d^{\text{H}}\), and \(d^{\text{L}}\) in Theorem (3) for over-approximation of the reachable set of the system (18) with neural network \(u=\mathsf{N}(x)\) trained to approximate an offline model predictive controller. The \(p_{x}-p_{y}\) plot of the motion of the vehicle is shown starting from an initial set \([7.9,8.1]^{2}\times[-\frac{2\pi}{3}-0.01,-\frac{2\pi}{3}+0.01]\times[1.99,2.01]\). In terms of accuracy, the size of the over-approximations obtained using the functions \(d^{\text{H}}\) and \(d^{\text{L}}\) are close to each other, but are much smaller than the size of the over-approximations obtained using the function \(d^{\text{G}}\). On the other hand, finding the over-approximations using the functions \(d^{\text{G}}\) and \(d^{\text{H}}\) take about the same amount of time and are much faster than finding the over-approximations using the functions \(d^{\text{L}}\). The runtimes are averaged over \(10\) instances and mean and standard deviation are reported.
over-approximation of the reachable sets using different integration techniques for Algorithm 1. The results are shown in Figure 3. Note that, our Algorithm 1 with the forward Euler method yields identical results as the linear programming approach in (Everett et al., 2021b).

Figure 3: The time evolution of the reachable sets of the 6D quadrotor with the neural network controller \(u=\mathsf{N}(x)\) in \(p_{x},p_{y},p_{z}\) coordinate. The blue hyper-rectangles are the over-approximations of the reachable sets of the system computed using the function \(d^{\mathrm{H}}\) in Theorem 3. **Left plot**: The ODEs are integrated using a Runge-Kutta method. **Middle plot**: The ODEs are integrated using forward Euler method with time step \(0.01\). **Right plot**: The ODEs are integrated using forward Euler method with time step \(0.1\) to match the actuation time step. Note that the reachable set estimates of this plot are identical to those in (Everett et al., 2021b, Figure 10(a)) where the over-approximations are computed using a linear program. We observed that, on the same computer, this linear program takes more than twice as long, \(0.113\pm 0.004\) seconds. The runtimes are averaged over \(10\) instances and mean and standard deviation are reported.

Figure 2: Performance of the Algorithm (1) with \(s=\mathrm{H}\) for different partitions \(D_{s}\) and sub-partitions \(D_{a}\) for over-approximation of the reachable set of the system (18) with neural network \(u=\mathsf{N}(x)\) trained to approximate an offline model predictive controller. All four figures show the \(p_{x}-p_{y}\) plot of the motion of the vehicle. The runtimes are averaged over \(10\) instances and the mean and standard deviation are reported.
## 8 Conclusions

We presented a fast and scalable method, based on mixed monotone theory, to over-approximate the reachable sets of nonlinear systems coupled with neural network controllers. We introduce three methods to intertwine neural network inclusion functions from CROWN (Zhang et al., 2018) and open-loop decomposition functions with varying empirical results. For future research, we plan to study the role of the neural network verification algorithms in the tightness of our approximations.
