Contraction-Guided Adaptive Partitioning for Reachability Analysis of
  Neural Network Controlled Systems
In this paper, we present a contraction-guided adaptive partitioning
algorithm for improving interval-valued robust reachable set estimates in a
nonlinear feedback loop with a neural network controller and disturbances.
Based on an estimate of the contraction rate of over-approximated intervals,
the algorithm chooses when and where to partition. Then, by leveraging a
decoupling of the neural network verification step and reachability
partitioning layers, the algorithm can provide accuracy improvements for little
computational cost. This approach is applicable with any sufficiently accurate
open-loop interval-valued reachability estimation technique and any method for
bounding the input-output behavior of a neural network. Using contraction-based
robustness analysis, we provide guarantees of the algorithm's performance with
mixed monotone reachability. Finally, we demonstrate the algorithm's
performance through several numerical simulations and compare it with existing
methods in the literature. In particular, we report a sizable improvement in
the accuracy of reachable set estimation in a fraction of the runtime as
compared to state-of-the-art methods.
Akash Harapanahalli, Saber Jafarpour, Samuel Coogan
Contraction-Guided Adaptive Partitioning for Reachability Analysis of Neural Network Controlled Systems

###### Abstract

In this paper, we present a contraction-guided adaptive partitioning algorithm for improving interval-valued robust reachable set estimates in a nonlinear feedback loop with a neural network controller and disturbances. Based on an estimate of the contraction rate of over-approximated intervals, the algorithm chooses when and where to partition. Then, by leveraging a decoupling of the neural network verification step and reachability partitioning layers, the algorithm can provide accuracy improvements for little computational cost. This approach is applicable with any sufficiently accurate open-loop interval-valued reachability estimation technique and any method for bounding the input-output behavior of a neural network. Using contraction-based robustness analysis, we provide guarantees of the algorithm's performance with mixed monotone reachability. Finally, we demonstrate the algorithm's performance through several numerical simulations and compare it with existing methods in the literature. In particular, we report a sizable improvement in the accuracy of reachable set estimation in a fraction of the runtime as compared to state-of-the-art methods.

## I Introduction

Motivation and Problem StatementNeural networks have become increasingly popular in control systems in recent years due to their relative ease of in-the-loop computation. These learning-based algorithms are known to be vulnerable to input perturbations [1]--small (possibly adversarial) changes in their input can lead to large variations in their output. As such, runtime verification of the safety and performance of neural network controlled systems is essential in safety-critical applications. This task is generally challenging due to the nonlinear and large-scale structure of the neural networks and their interconnection with nonlinear dynamics [2].

A basic ingredient for verifying control systems is the ability to overapproximate the set of reachable states from a given set of initial conditions, possibly in the presence of disturbances. If, for example, this overapproximation avoids obstacles or reaches a goal region, then the system certifiably satisfies the corresponding safety or performance criteria.

Recently, several promising reachability-based methods have been proposed for verifying stand-alone neural networks or feedback systems with neural networks in the control loop, however, these methods either suffer from large computational complexity, large over-approximation error, or lack of generality. For stand-alone neural networks: Interval Bound Propagation (IBP) [3] is fast, but largely over-conservative; CROWN [4] suffers when the input perturbation is large; LipSDP [5] is not scalable to large networks. For neural network closed-loop verification: ReachLP [6] is computationally light, but it can be conservative and only applies to linear systems; ReachSPP [7] can provide tighter bounds, but does not scale well to large networks and only applies to linear systems; POLAR [8], MILP [9], and constrained Zonotope [10] methods can provide very accurate estimates for non-linear systems, but they are too expensive for runtime verification.

Partitioning of the state space is an effective approach for balancing the trade-off between the accuracy of reachable set estimates and the runtime of the verification strategy. In the literature, there are many existing partitioning algorithms. For standalone ReLU networks, [11] splits the input set into partitions based on the stability of each neuron. In [12] and GSG [13], Monte Carlo simulations

Fig. 1: A snapshot of Algorithm 1 is illustrated for \(n=1\). **(top)** There are two initial partitions that both compute the neural network verification step (filled circles in the graph representation). Both partitions are integrated from \(t_{j-1}\) to \(t_{\gamma}=t_{j-1}+\gamma(t_{j}-t_{j-1})\) to estimate the width of the box at \(t_{j}\). Since the estimated blue partition width violates the user-defined maximum width \(\varepsilon\), **(bottom)** the algorithm returns to \(t_{j-1}\) and adds two new partitions to the tree. The maximum neural network verification depth here is \(1\), so the new partitions do not recompute the neural network verification and instead use the same \(E^{c}\) from the initial blue partition.
specific state axes to cut. For closed-loop linear systems, ReachLP-Uniform [6] uses a uniform partitioning of the state space, ReachLP-GSG [6] applies the GSG algorithm to ReachLP, and ReachLipBnB [14] combines a branch and bound algorithm with \(\ell_{2}\)-Lipschitz constants obtained from LipSDP.

ContributionsWe introduce a contraction-guided adaptive partitioning algorithm for general nonlinear systems with neural network controllers. The algorithm is (i) _separative_--it decouples the neural network verification calls and reachability partitioning layers, leading to accuracy benefits for little computational cost; (ii) _spatially aware_--it chooses specific sections of the state space to cut based on the contraction rate of the partitions; and (iii) _temporally aware_--it chooses to cut on specific time intervals along the trajectory as needed. Additionally, while our algorithm can be coupled with any interval-based verification method including many of those discussed above, we propose specifically applying it to the recently proposed verification method ReachMM based on mixed monotone system theory [15], for which we call the resulting algorithm ReachMM-CG. In this case, we show that the three features highlighted above are explainable through the lens of contraction theory with provable guarantees. In particular, we prove a connection between the contraction rate of the closed-loop mixed monotone embedding system and the contraction rate of the partitions. Using this observation, we provide upper bounds on the approximation error and develop theoretical guarantees on the algorithm's performance. Moreover, ReachMM-CG is applicable in continuous or discrete time for nonlinear plant models controlled by neural networks with general nonlinear activation functions, which removes several limitations found in most of the above-cited approaches. Finally, we find that ReachMM-CG yields reachable set estimates with a 33% accuracy improvement in a quarter of the time compared to state-of-the-art partitioning algorithms on a benchmark example.

## II Notation and Mathematical Preliminary

We define the partial order \(\leq\) on \(\mathbb{R}^{n}\) as \(x\leq y\iff x_{i}\leq y_{i}\,\forall i\in\{1,\ldots,n\}\). For every \(x\leq y\), we define the interval \([x,y]=\{z:x\leq z\leq y\}\). For \(x,\widehat{x}\in\mathbb{R}^{n}\), let \((x,\widehat{x})\in\mathbb{R}^{2n}\) be their concatenation. The southeast partial order \(\leq_{\text{SE}}\) on \(\mathbb{R}^{2n}\) is induced by \(\leq\) on \(\mathbb{R}^{n}\) as follows: \((x,\widehat{x})\leq_{\text{SE}}(y,\widehat{y})\iff x\leq y\) and \(\widehat{y}\leq\widehat{x}\). Define the following: \(\mathcal{T}_{\geq 0}^{2n}=\{(x,\widehat{x})\in\mathbb{R}^{2n}:x\leq\widehat{x}\}\), \(\mathcal{T}_{\leq 0}^{2n}=\{(x,\widehat{x})\in\mathbb{R}^{2n}:\widehat{x}\leq \widehat{x}\},\ \mathcal{T}^{2n}=\mathcal{T}_{\geq 0}^{2n}\cup\mathcal{T}_{\leq 0}^{2n}\). We define the vector \(x_{[i:\widehat{x}]}\in\mathbb{R}^{n}\) as \(\left(x_{[i:\widehat{x}]}\right)_{j}=\begin{cases}x_{j},&j\neq i\\ \hat{x}_{j}&j=i\end{cases}\). Define the weighted \(\ell_{\infty}\)-norm \(\|x\|_{\infty,\varepsilon}=\|\operatorname{diag}(\varepsilon)^{-1}x\|_{\infty}\). In particular, the weighted maximum width of an interval \([\underline{x},\overline{x}]\) is \(\|\overline{x}-\underline{x}\|_{\infty,\varepsilon}\). Given a matrix \(A\in\mathbb{R}^{n\times n}\), the \(\ell_{\infty}\)-matrix measure of \(A\) is defined by \(\mu_{\infty}(A)=\max_{i\in\{1,\ldots,n\}}\{A_{ii}+\sum_{j\neq i}|A_{ij}|\}\). From [16, Table III], we define the weak pairing \(\llbracket\cdot,\cdot\rrbracket_{\infty}:\mathbb{R}^{n}\times\mathbb{R}^{n} \rightarrow\mathbb{R}\) associated to the norm \(\|\cdot\|_{\infty}\) as follows:

\[\llbracket x,y\rrbracket_{\infty}=\max_{i\in I_{\infty}(y)}y_{i}x_{i},\]

where \(I_{\infty}(x)=\{i\in\{1,\ldots,n\}\ |\ |x_{i}|=\|x\|_{\infty}\}\). Consider the following dynamical system

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

with state vector \(x\in\mathbb{R}^{n}\) and disturbance vector \(w\in\mathcal{W}\subset\mathbb{R}^{q}\). Given a piecewise continuous curve \(w(\cdot)\) where \(w:[t_{0},t]\rightarrow\mathcal{W}\), the trajectory of the system starting from \(x_{0}\) at time \(t_{0}\) is given \(t\rightarrow\phi_{f}(t,t_{0},x_{0},w(\cdot))\). Given an initial set \(\mathcal{X}\), we denote the reachable set of \(f\) at some \(t\geq t_{0}\):

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

The dynamical system (1) is mixed monotone with respect to the decomposition function \(d:\mathcal{T}^{2n}\times\mathcal{T}^{2q}\rightarrow\mathbb{R}^{n}\) if, for every \(i\in\{1,\ldots,n\}\),

1. \(d_{i}(x,x,w,w)=f_{i}(x,w)\), for every \(x\in\mathbb{R}^{n}\), \(w\in\mathbb{R}^{q}\);
2. \(d_{i}(x,\widehat{x},w,\widehat{w})\leq d_{i}(y,\widehat{y},w,\widehat{w})\), for every \(x\leq y\) s.t. \(x_{i}=y_{i}\), and every \(\widehat{y}\leq\widehat{x}\);
3. \(d_{i}(x,\widehat{x},w,\widehat{w})\leq d_{i}(x,\widehat{x},v,\widehat{v})\), for every \(w\leq v\) and every \(\widehat{v}\leq\widehat{w}\).

With a valid decomposition function \(d\), one can construct an embedding system associated to (1) as follows:

\[\frac{d}{dt}\left[\frac{x}{\overline{x}}\right]=\begin{bmatrix}d(\underline{x},\overline{x},\underline{w},\overline{w})\\ d(\overline{x},\underline{x},\overline{w},\underline{w})\end{bmatrix}:=E( \underline{x},\overline{x},\underline{w},\overline{w}). \tag{3}\]

## III Problem Statement

Consider a nonlinear continuous-time dynamical system of the form

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

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

\[\xi^{(i)}=\sigma^{(i-1)}\left(W^{(i-1)}\xi^{(i-1)}+b^{(i-1)}\right), \,i=1,\ldots k \tag{5}\] \[\xi^{(0)}=x,\quad N(x)=W^{(k)}\xi^{(k)}+b^{(k)}\]

where \(m_{i}\) is the number of neurons in the \(i\)-th layer, \(W^{(i)}\in\mathbb{R}^{m_{i}\times m_{i-1}}\) is the weight matrix on the \(i\)-th layer, \(b^{(i)}\in\mathbb{R}^{m_{i}}\); is the bias vector on the \(i\)-th layer, \(\xi^{(i)}\in\mathbb{R}^{m_{i}}\) is the \(i\)-th layer hidden variable and \(\sigma_{i}\) is the activation function for the \(i\)-th layer. In practice, \(N(x(t))\) cannot be evaluated at each instance of time \(t\), and, instead, the control must be implemented via, _e.g._, a zero-order hold strategy between sampling instances. We assume that there exists an increasing sequence of control time instances \(\{t_{0},t_{1},t_{2},\ldots\}\) in which the control input is updated. Thus, the closed-loop system with the neural-network feedback controller is given by:

\[\dot{x}(t)=f(x(t),N(x(t_{j})),w(t)):=f^{c}(j,x(t),w(t)),\]
where \(t\in\mathbb{R}_{\geq 0}\) and \(j\in\mathbb{Z}_{\geq 0}\) is such that \(t\in[t_{j},t_{j+1}]\). In our analysis when \(j\) is clear from context or does not affect the result, we drop \(j\) as an argument of \(f^{c}\).

The goal of this paper is to verify the behavior of the closed-loop system (6). To verify the safety of a system under uncertainty, one needs to verify the entire reachable set. However, in general, computing the reachable set exactly is not computationally tractable--instead, approaches typically compute an over-approximation \(\overline{\mathcal{R}}_{f^{c}}(t,t_{0},\mathcal{X},\mathcal{W})\supseteq \mathcal{R}_{f^{c}}(t,t_{0},\mathcal{X},\mathcal{W})\). Therefore, the main challenge addressed in this paper is to develop an approach for providing tight over-approximations of reachable sets while remaining computationally tractable for runtime computation.

## IV Interval Reachability of Neural Network Controlled Systems

### _General Framework_

We assume we have access to an off-the-shelf dynamical system reachability tool that supports _interval analysis_.

**Assumption 1** (Open-loop interval reachability).: _Given a dynamical system of the form (4), any intervals \(\mathcal{X}_{0}=[\underline{x}_{0},\overline{x}_{0}]\subseteq\mathbb{R}^{n}\), \(\mathcal{U}=[\underline{u},\overline{u}]\subseteq\mathbb{R}^{p}\), and \(\mathcal{W}=[\underline{w},\overline{w}]\subseteq\mathbb{R}^{q}\), and some initial time \(t_{0}\), there exists a reachability algorithm that returns a valid interval approximation \([\underline{x}(t),\overline{x}(t)]\) satisfying_

\[\mathcal{R}_{f}(t,t_{0},\mathcal{X}_{0},\mathcal{U},\mathcal{W})\subseteq[ \underline{x}(t),\overline{x}(t)],\quad\forall t\geq t_{0}. \tag{7}\]

Reachability analysis of dynamical systems is a classical and well-studied research field with several off-the-shelf toolboxes with this capability, including Flow*[17], Hamilton-Jacobi approach [18], level set toolbox [19], CORA [20], and mixed monotonicity [21].

We also assume access to off-the-shelf neural network verification algorithms that can provide _interval inclusion functions_ as follows:

**Assumption 2** (Neural network verification).: _Given a neural network \(N\) of the form (5) and any interval \([\underline{y},\overline{y}]\subseteq\mathbb{R}^{n}\), there exists a neural network verification algorithm that returns a valid inclusion function \(\left[\frac{\underline{N}_{[\underline{y},\overline{y}]}}{\underline{N}_{[ \underline{y},\overline{y}]}}\right]:\mathcal{T}_{\geq 0}^{2n}\to\mathcal{T}_{\geq 0}^{2p}\) satisfying_

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

_for any \(x\in[\underline{x},\overline{x}]\subseteq[\underline{y},\overline{y}]\)._

A large number of the existing neural network verification algorithms can provide bounds of the form (8) for the output of the neural networks, including CROWN (and its subsequent variants) [4], LipSDP [22], and IBP [3].

Combining Assumptions 1 and 2 leads naturally to an algorithm for over-approximating solutions to (6). In particular, starting with \(j=0\) and an initial interval of states \(\mathcal{X}_{j}=[\underline{x}(t_{j}),\overline{x}(t_{j})]\), first obtain \(\underline{N}_{[\underline{y},\overline{y}]}\) and \(\overline{N}_{[\underline{y},\overline{y}]}\) for some \([\underline{y},\overline{y}]\supseteq\mathcal{X}_{j}\) from Assumption 2. Next, set \(\mathcal{U}=[\underline{N}_{[\underline{y},\overline{y}]}(\underline{x}(t_{j}),\overline{x}(t_{j})),\overline{N}_{[\underline{y},\overline{y}]}(\underline{ x}(t_{j}),\overline{x}(t_{j}))]\). Finally, compute the reachable set of the closed-loop system on the interval \(t\in(t_{j},t_{j+1}]\) as the interval reachable set obtained from the algorithm of Assumption 1, set \(\mathcal{X}_{j+1}\) as the reachable set at time \(t_{j+1}\), increment \(j\gets j+1\), and iterate. This iteration serves as the backbone of our proposed paritioning-based algorithm below.

### _Mixed Monotone Framework_

One specific framework that satisfies both Assumptions 1 and 2 is developed in [15], where interval reachability of neural network controlled systems is studied through the lens of mixed monotone dynamical systems theory. Suppose that we have access to a decomposition function \(d\) for the open-loop system (4) with the open-loop embedding system

\[\frac{d}{dt}\begin{bmatrix}\underline{x}\\ \overline{x}\end{bmatrix}=\begin{bmatrix}d(\underline{x},\overline{x}, \underline{u},\overline{u},\underline{w},\overline{w})\\ d(\overline{x},\overline{x},\overline{u},\underline{u},\overline{w}, \overline{w})\end{bmatrix}:=E(\underline{x},\overline{x},\underline{u}, \overline{u},\underline{w},\overline{w}). \tag{9}\]

Given \(\mathcal{X}_{0}=[\underline{x}_{0},\overline{x}_{0}]\subseteq\mathbb{R}^{n}\), \(\mathcal{U}=[\underline{u},\overline{u}]\subseteq\mathbb{R}^{p}\), and \(\mathcal{W}=[\underline{w},\overline{w}]\subseteq\mathbb{R}^{q}\), the trajectory of (9) starting from \((\underline{x}_{0},\overline{x}_{0})\) provides the inclusion (7) from Assumption 1.

Following the treatment in [15], one can use CROWN [4] to obtain the desired bounds from Assumption 2. Given an interval \([\underline{y},\overline{y}]\), the algorithm provides an efficient procedure for finding a tuple \((\underline{C},\overline{C},\underline{d},\overline{d})\) defining linear upper and lower bounds for the output of the neural network

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

for every \(x\in[\underline{y},\overline{y}]\). Using these linear bounds, we can construct the inclusion function for any \([\underline{x},\overline{x}]\subseteq[\underline{y},\overline{y}]\):

\[\begin{split}&\underline{N}_{[\underline{y},\overline{y}]}( \underline{x},\overline{x})=[\underline{C}(\underline{y},\overline{y})]^{+} \underline{x}+[\underline{C}(\underline{y},\overline{y})]^{-}\overline{x}+ \underline{d}(\underline{y},\overline{y}),\\ &\overline{N}_{[\underline{y},\overline{y}]}(\underline{x},\overline{x })=[\overline{C}(\underline{y},\overline{y})]^{+}\overline{x}+[\overline{C}( \underline{y},\overline{y})]^{-}\underline{x}+\overline{d}(\underline{y}, \overline{y}).\end{split} \tag{11}\]

By combining these two tools, one can construct a new closed-loop embedding system to over-approximate the reachable sets of the closed-loop system (6). From [15], with the following definitions,

\[\begin{split}&\underline{\eta}_{j}=\underline{N}_{[\underline{y}, \overline{y}]}(\underline{x}(t_{j}),\overline{x}(t_{j})_{[i:\underline{x}(t_{j})]}), \\ &\overline{\eta}_{j}=\overline{N}_{[\underline{y},\overline{y}]}( \underline{x}(t_{j}),\overline{x}(t_{j})_{[i:\underline{x}(t_{j})]}),\\ &\underline{\nu}_{j}=\underline{N}_{[\underline{y},\overline{y}]}( \overline{x}(t_{j}),\underline{x}(t_{j})_{[i:\overline{x}(t_{j})]}),\\ &\overline{\nu}_{j}=\overline{N}_{[\underline{y},\overline{y}]}( \overline{x}(t_{j}),\underline{x}(t_{j})_{[i:\overline{x}(t_{j})]}),\end{split} \tag{12}\]

we use the following "hybrid" function

\[\begin{split}\big{(}d^{c}_{[\underline{y},\overline{y}]}(j, \underline{x},\overline{x},\underline{w},\overline{w})\big{)}_{i}=\begin{cases}d_{i}( \underline{x},\overline{x},\overline{\eta}_{j},\overline{\eta}_{j},\underline{w}, \overline{w}),&\underline{x}\leq\overline{x}\\ d_{i}(\overline{x},\overline{x},\overline{\nu}_{j},\underline{\nu}_{j}, \overline{w},\overline{w}),&\overline{x}\leq\underline{x}\end{cases}\end{split} \tag{13}\]

for any \([\underline{y},\overline{y}]\supseteq[\underline{x},\overline{x}]\), to create the closed-loop embedding system:

\[\begin{split}&\frac{d}{dt}\begin{bmatrix}\underline{x}\\ \overline{x}\end{bmatrix}=\begin{bmatrix}d^{c}_{[\underline{y},\overline{y}]}(j, \underline{x},\overline{x},\underline{w},\overline{w})\\ d^{c}_{[\underline{y},\overline{y}]}(j,\overline{x},\underline{w}, \overline{w})\end{bmatrix}:=E^{c}_{[\underline{y},\overline{y}]}(j,\underline{x}, \overline{x},\underline{w},\overline{w}).\end{split} \tag{14}\]

In our analysis when \(j\) is clear from context or does not affect the result, we drop \(j\) as an argument of \(d^{c}\) and \(E^{c}\). For an input set \(\mathcal{X}_{0}\subseteq[\underline{x}_{0},\overline{x}_{0}]\) and a disturbance set \(\mathcal{W}\subseteq[\underline{w},\overline{w}]\), with \(\left[\frac{x}{x}(t)\right]=\phi_{E^{c}_{[\underline{y},\overline{y}]}}\left(t _{0},\left[\frac{\underline{x}_
following inclusion holds [15, Theorem 1]:

\[\mathcal{R}_{f^{c}}(t,t_{0},\mathcal{X},\mathcal{W})\subseteq[\underline{x}(t), \overline{x}(t)]. \tag{15}\]

Thus, running one single trajectory of \(E^{c}\) can over-approximate the reachable set of the closed-loop system (6).

### _Partitioning_

While interval analysis techniques can often be computationally inexpensive, as the size of the uncertain sets grow, the reachable set estimates become overly conservative. In a closed-loop, these effects can compound significantly as the error continues to accumulate. This phenomenon is known in the literature as the _wrapping effect_. To mitigate over-conservatism in interval-based techniques, one can split uncertain regions into smaller partitions. A valid over-approximation of the reachable set can be found by taking the union of the over-approximations of the reachable sets for each partition. Due to the locality of the smaller partitions, this technique can drastically reduce the conservatism of interval-based techniques.

## V Contraction-Guided Adaptive Partitioning

In this section, we introduce a contraction-guided adaptive partitioning algorithm to improve the results of interval reachability methods satisfying Assumptions 1 and 2.

### _Algorithm Description_

At a particular control step \(t_{j}\), define a partition

\[P=((\underline{x}_{0},\overline{x}_{0}),\mathrm{N},\mathcal{S}), \tag{16}\]

where \((\underline{x}_{0},\overline{x}_{0})\in\mathcal{T}_{\geq 0}^{2n}\) is the initial condition of the partition, \(\mathrm{N}\in\{\mathrm{True},\mathrm{False}\}\) is whether or not the partition should compute the neural network verification (8), and \(\mathcal{S}\) is a set of tuples representing its subpartitions. Algorithm 1 starts by initializing a tree at time \(t_{0}\) with a single root node \(\mathcal{P}=((\underline{X}_{0},\overline{X}_{0}),\mathrm{True},\emptyset)\) representing the initial set.

step( ) _Procedure:_ Along a particular control interval \([t_{j-1},t_{j}]\), every partition undergoes the following procedure. If specified, the embedding system \(E^{c}\) is computed using (11, 13, 14). Initially, \(E^{c}\) is integrated to a fraction \(\gamma\in(0,1]\) of the control interval, where the contraction rate is computed by comparing its width to its initial condition. Then, the width of the final box at \(t_{j}\) is estimated. If it violates hyper-parameter \(\varepsilon\), then the algorithm returns to time \(t_{j-1}\) and divides the initial condition into \(2^{n}\) sub-partitions. This process is repeated for each of the sub-partitions. Finally, if the estimated width does not violate \(\varepsilon\), \(E^{c}\) is fully integrated to \(t_{j}\). This procedure is visualized in Figure I, and formalized as step() in Algorithm 1, which repeats for every control step until some final time \(T\).

```
Input: initial set: \(\mathcal{X}_{0}=[\underline{X}_{0},\overline{X}_{0}]\subset\mathbb{R}^{n}\), control time instances: \(\{t_{0},t_{1},\dots,\}\), final time: \(T\) Parameter: desired width \(\varepsilon\in\mathbb{R}_{>0}^{n}\), check contraction factor \(\gamma\in(0,1]\), max partition depth \(D_{p}\in\mathbb{Z}_{\geq 0}\), max neural network verification depth \(D_{\mathrm{N}}\in\mathbb{Z}_{\geq 0}\) Output: Over-approximated reachable set trajectory for \(t\in[t_{0},T]\); \(\overline{\mathcal{R}}(t)\)
1\(\overline{\mathcal{R}}(t_{0})=\mathcal{X}_{0}\):
2\(\mathcal{P}=((\underline{X}_{0},\overline{X}_{0}),\mathrm{True},\emptyset)\);
3\(t_{m}=\) smallest control time instance \(\geq T\);
4for\(j=1,\dots,m\)do
5\((\overline{\mathcal{R}}(\cdot)|_{[t_{j-1},t_{j}]},\mathcal{P}_{+})\leftarrow \texttt{step}\left(\mathcal{P},\emptyset,j\right)\);
6\(\mathcal{P}\leftarrow\mathcal{P}_{+}\);
7
8 end for
9return\(\overline{\mathcal{R}}(t)\)\(\forall t\in[t_{0},T]\);
10
11Procedure step\(((\underline{x}_{0},\overline{x}_{0}),\mathrm{N},\mathcal{S}),E^{c},j)\)
12ifN is Truethen // NN Verification
13\(E^{c}\gets E^{c}\) as (14, 13) using (11) on \((\underline{x}_{0},\overline{x}_{0})\);
14if\(\mathcal{S}=\emptyset\)then // No Subpartitions
15\(d=\mathrm{get\_depth}()\) ; // Tree Depth
16if\(d<D_{p}\)then // Can Partition
17\(t_{\gamma}\gets t_{j-1}+\gamma(t_{j}-t_{j-1})\);
18\((\underline{x},\overline{x})(t_{\gamma})\leftarrow\phi_{E^{c}}(t_{\gamma},t_ {j-1},(\underline{x}_{0},\overline{x}_{0}))\);
19\(C=\frac{\|\overline{x}(t_{\gamma})-\underline{x}(t_{\gamma})\|_{\infty, \varepsilon}}{\|\underline{x}_{0}-\underline{x}_{0}\|_{\infty,\varepsilon}}\); // Contr Rate
20if\(C^{1/\sqrt{\|}}\overline{x}_{0}\max_{\emptyset}\|_{\infty,\varepsilon}>1\)then
21\(\{\underline{x}_{0}^{k},\overline{x}_{0}^{k}\}_{k=1}^{2}\leftarrow\) uni_div\((\underline{x}_{0},\overline{x}_{0})\);
22\(\mathcal{S}\leftarrow\{((\underline{x}_{0},\overline{x}_{0}^{k}),d<D_{\mathrm{N}}, \emptyset)\}_{k=1}^{2n}\);
23\(\mathrm{N}\leftarrow\mathrm{N}\land((d+1)>D_{\mathrm{N}})\);
24return step\(((\underline{x}_{0},\overline{x}_{0}),\mathrm{N},\mathcal{S}),E^{c},j)\)
25\(((\underline{x},\overline{x})(\cdot)|_{[t_{j-1},t_{j}]}\leftarrow\phi_{E^{c}}( \cdot,t_{j-1},(\underline{x}_{0},\overline{x}_{0}))\);
26return\(([\underline{x},\overline{x}](\cdot)|_{[t_{j-1},t_{j}]},((\underline{x}, \overline{x})(t_{f}),\mathrm{N},\emptyset))\)
27else // Iterate Subpartitions
28\(\{(\overline{\mathcal{R}}^{k}(\cdot),P_{+}^{k})\leftarrow\texttt{step}\left( P^{k},E^{c},j\right)\}_{P^{k}\in\mathcal{S}}\);
29\((\underline{x}_{0},\overline{x}_{0})_{+}\leftarrow\) tightest \([\underline{x},\overline{x}]\supseteq\bigcup_{k}\overline{\mathcal{R}}^{k}(t_{f})\);
30return\((\bigcup_{k}\overline{\mathcal{R}}^{k}(\cdot),((\underline{x}_{0}, \overline{x}_{0})_{+},\mathrm{N},\{P_{+}^{k}\}_{k=1}^{2n}))\);
```

**Algorithm 1**Contraction-Guided Adaptive Partitioning (ReachMM-CG)

Hyper-parametersThere are several user-defined parameters that are important for the algorithm's performance. The choice of the check-contraction factor \(\gamma\in(0,1]\) is purely for computational benefit, and in particular, \(\gamma=1\) checks the true width of the box at \(t_{j}\) instead of an estimate. The maximum partition depth \(D_{p}\in\mathbb{Z}_{\geq 0}\) specifies how deep new partitions can be added in the tree to integrate \(E^{c}\). The maximum neural network verification depth \(D_{\mathrm{N}}\in\mathbb{Z}_{\geq 0}\) specifies how deep neural network verification is performed in the partition tree. \(D_{p}\) and \(D_{\mathrm{N}}\) are directly related to the first two terms of the right-hand side of Therorem 1 in the next section--every partitioning layer will improve the first term, while only partitions computing the neural network verification will improve the second term. Finally, the choice of desired width \(\varepsilon\in\mathbb{R}_{\geq 0}^{n}\) specifies how wide partitions can grow before being cut: if chosen too small, the algorithm will partition early in the start (\(\varepsilon_{i}=0\) for any \(i\) implies a uniform initial partitioning); if chosen too large, the algorithm will not partition at all (\(\varepsilon_{i}=\infty\) implies no condition on the \(i\)-th component).
### _Discussion_

Algorithm 1 has three key features.

SeparationIn practice, the neural network verification step usually contributes the most computational expense to the algorithm. Introducing the maximum depths \(D_{p}\) and \(D_{\text{N}}\) allows us to separate the calls to the neural network verification algorithm and the total number of partitions. This allows the algorithm to improve its accuracy without any additional calls to the neural network verifier, leaving it to the user to control the trade-off between computational complexity and accuracy.

Spatial awarenessBased on the dynamics of the system, the over-approximation error can be highly sensitive to the spatial location of each partition. By estimating the contraction rate of each partition separately, the algorithm can localize partitioning exactly _where_ necessary.

Temporal awarenessOne of the main drawbacks of interval reachability frameworks is that as time evolves, the over-approximation error tends to compound exponentially. As such, by partitioning along trajectories, the algorithm can save computations to improve reachable set estimates exactly _when_ necessary.

Generality of Algorithm 1The algorithm, as written, specifically uses the Mixed Monotone framework (developed in Section V.B) for the interval reachability of the dynamical system, as well as CROWN (specified by (10)) for the neural network verification. We will exclusively use this setting for the rest of the paper, and refer it as "ReachMM-CG". However, it is easy to see that one can replace Lines 11, 16, and 23 with any other interval reachability setting satisfying Assumptions 1 and 2.

## VI Contraction-Based Guarantees on Reachable Set Over-Approximation

In this section, we use contraction theory to provide rigorous bounds on the accuracy of the mixed-monotone reachable set over-approximations. For this section, we remove the assumption that the control is applied in the piecewise constant fashion and instead use a continuously applied neural network. This is done for the simplicity of analysis.

**Theorem 1** (Accuracy guarantees for mixed-monotone reachability).: _Consider the closed-loop system (6) with the neural network controller \(u=N(x)\) given by (5) and let \(t\mapsto x(t)\) be a trajectory of the closed-loop system (6). Suppose that \(d\) is the decomposition function for the open-loop system (4), there is a neural network verification algorithm which provides the bounds of the form (8). Let and be a trajectory of the closed-loop embedding system (14) and \(t\mapsto\left[\frac{y}{y}\right]\) is a curve such that, for every \(t\in\mathbb{R}_{\geq 0}\). Then,_

(17)

_where_

\[c_{x} =\sup_{\underline{z},\overline{y}\in\Omega_{t}}\mu_{\infty}(D_{[ \frac{z}{2}]}E^{\text{c}}_{[\underline{y}(t),\overline{y}(t)]}(\underline{z}, \overline{z},\underline{w},\overline{w})),\] \[\ell^{o}_{u} =\sup_{\begin{subarray}{c}x\in\Omega_{t},w\in[\underline{w}, \overline{w}],\\ \underline{u},\overline{w}\in[\underline{N}(x,x),\overline{N}(x,x)]\end{subarray} }\|D_{[\frac{z}{u}]}E^{\text{c}}(x,x,\underline{u},\overline{w},w,w)\|_{\infty},\] \[\ell^{c}_{w} =\sup_{x\in\Omega_{t},\underline{w},\overline{z}\in[\underline{ w},\overline{w}],\\ \|D_{[\frac{z}{u}]}E^{\text{c}}_{[\underline{y}(t),\overline{y}(t)]}(x,x, \underline{z},\overline{z})\|_{\infty},\]

_and the set \(\Omega_{t}\) is defined by \(\Omega_{t}=\bigcup_{\tau\in[0,t]}[\underline{y}(\tau),\overline{y}(\tau)] \subseteq\mathbb{R}^{n}\)._

Proof.: We first make the following conventions:

\[r(t) =\begin{bmatrix}\underline{x}(t)\\ \overline{x}(t)\end{bmatrix}-\begin{bmatrix}x(t)\\ x(t)\end{bmatrix},\] \[\text{1} =E^{\text{c}}_{[\underline{y},\overline{y}]}(\underline{x}, \overline{x},\underline{w},\overline{w})-E^{\text{c}}_{[\underline{y}, \overline{y}]}(x,x,\underline{w},\overline{w}),\] \[\text{2} =E^{\text{c}}_{[\underline{y},\overline{y}]}(x,x,\underline{w}, \overline{w})-E^{\text{c}}_{[\underline{y},\overline{y}]}(x,x,w,w),\] \[\text{3} =E^{\text{c}}_{[\underline{y},\overline{y}]}(x,x,w,w)-\begin{bmatrix} f^{\text{c}}(x,w)\\ f^{\text{c}}(x,w)\end{bmatrix}.\]

Therefore, we can compute

\[\frac{d}{dt}\|r(t)\|_{\infty}^{2} =\left[E^{\text{c}}_{[\underline{y},\overline{y}]}(\underline{x},\overline{x},\underline{w},\overline{w})-\begin{bmatrix}f^{\text{c}}(x,w)\\ f^{\text{c}}(x,w)\end{bmatrix},r(t)\right]_{\infty}\] \[\leq\left[\text{1},r(t)\right]_{\infty}+\left[\text{2},r(t) \right]_{\infty}+\left[\text{3},r(t)\right]_{\infty},\]

where the first equality is by the curve-norm derivative property of the weak pairing [16, Theorem 25(ii)] and the second inequality is by the subadditivity of weak pairing [16, Definition 15, property (i)]. By [16, Theorem 18], the first term on the RHS of the above equation can be estimated. By [16, Definition 15, property (iv)], the second term in the RHS can be estimated

Fig. 2: A partition tree structure of Algorithm 1 for a run on the double integrator system (20) for \(\varepsilon=0.1\), \(D_{p}=10\), \(D_{\text{N}}=2\). The blue color represents the algorithm’s _separation_—only nodes filled with blue compute the neural network verification step, and all the integrations are performed only on the leaf nodes. The imbalanced structure is a consequence of the algorithm’s _spatial awareness_. As a consequence of the algorithm’s _temporal awareness_, the structure of the partition tree deepens as time increases.
as

\[\left[\vbox{\hbox{\includegraphics[height=14.226378pt]{figs/
a step size of \(0.01\), and CROWN is computed using auto_LiRPA [25]. In Table I, the performance of ReachMM-CG with different selections of hyper-parameters (all with check contraction factor \(\gamma=0.1\)) is compared to ReachMM, the non-adaptive purely initial set partitioning strategy from [15]. The runtimes are averaged over 100 runs, with mean and standard deviation reported. The volumes for the tightest over-bounding box for the approximated reachable set at the final time step is reported. The adaptive strategies see an improved run-time over the non-adaptive counterparts, and in some cases see a tighter approximation.

### _Linear Discrete-Time Double Integrator_

Consider a zero-order hold discretization of the classical double integrator with a step-size of \(1\) (from [6]):

\[x_{t+1}=\underbrace{\begin{bmatrix}1&1\\ 0&1\end{bmatrix}}_{A}x_{t}+\underbrace{\begin{bmatrix}0.5\\ 1\end{bmatrix}}_{B}u_{t} \tag{20}\]

Here, we consider a fixed actuation step size of \(1\)--the same as the integration step. In this special case (discrete-time LTI systems), [15, Corollary 4] shows that the following is a valid closed-loop embedding system that provides tighter bounds than (14):

\[\begin{bmatrix}\underline{x}_{t+1}\\ \overline{x}_{t+1}\end{bmatrix}=\begin{bmatrix}\underline{M}^{+}\\ \overline{M}^{-}\end{bmatrix}\begin{bmatrix}\underline{x}_{t}\\ \overline{x}_{t}\end{bmatrix}+\begin{bmatrix}B^{+}&B^{-}\\ B^{-}&B^{+}\end{bmatrix}\begin{bmatrix}\underline{d}_{t}\\ \overline{d}_{t}\end{bmatrix}\]

with \(\underline{C},\overline{C},\underline{d},\overline{d}\) taken from CROWN as (10) on \((\underline{x}_{t},\overline{x}_{t})\), and \(\underline{M}=A+B^{+}\underline{C}+B^{-}\overline{C}\) and \(\overline{M}=A+B^{+}\overline{C}+B^{-}\underline{C}\).

We use the neural network from [6] (\(2\times 10\times 5\times 1\), ReLU). The performance of Algorithm 1 for various choices of \(\varepsilon\), \(D_{p}\), and \(D_{\text{N}}\) is shown in Figure 4, for the initial set \([2.5,3]\times[-0.25,0.25]\) and a final time of \(T=5\).

Comparing with the literatureAdditionally, we compare the proposed ReachMM-CG to state-of-the-art partitioning algorithms for linear discrete-time systems: ReachMM [15] (uniform initial partitioning), ReachLP-Uniform [6] (ReachLP with uniform initial partitioning), ReachLP-GSG (ReachLP with greedy sim-guided partitioning) [6], and ReachLipBnB [14] (branch-and-bound using LipSDP [5]).

Each algorithm is run with two different sets of hyper-parameters, aiming to compare their performances across various regimes. The setup for ReachMM-CG is \((\varepsilon,\,D_{p},\,D_{\text{N}})\); ReachMM is \((D_{p},\,D_{\text{N}})\); ReachLP-Uniform is \(\#\) initial partitions, ReachLP-GSG is \(\#\) of total propogator calls, ReachLipBnB is \(\varepsilon\). Their runtimes are averaged over 100 runs, with mean and standard deviation reported. The true areas of the reachable sets are computed using Python packages (Shapely, polytope). The performance is outlined in Table II, and notable reachable sets are displayed in Figure 5. ReachMM-CG outperforms SOTA across the board: for both setups, ReachMM-CG is significantly faster than the other methods while returning a tighter reachable set.

Fig. 4: The over-approximated reachable sets of the closed-loop double integrator model (20) are shown in blue for the initial set \([2.5,3]\times[-0.25,0.25]\) and final time \(T=5\). They are computed using Algorithm 1 with the specified parameters in the title of each plot. The average runtime and standard deviation across 100 runs is reported, as well as the true area of the final reachable set. 200 true trajectories are shown in red. The horizontal axis is \(x_{1}\) and the vertical axis is \(x_{2}\).

Fig. 3: The over-approximated reachable sets of the closed-loop nonlinear vehicle model (18) in the \((p_{x},p_{y})\) coordinates are shown in blue for the 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]\) over the time interval \([0,1.25]\). They are computed using Algorithm 1 with \(\varepsilon=[0.2,0.2,\infty,\infty]^{\top}\), \(D_{p}=2\), and \(D_{\text{N}}=1\). The average runtime across 100 runs with standard deviation is reported, as well as the volume of the over-bounding box at the final time \(T=1.25\). 200 true trajectories of the system are shown in the time-varying yellow line.
## VIII Conclusions

In this paper, we propose an adaptive partitioning approach for interval reachability analysis of neural network controlled systems. The algorithm uses _separation_ of the neural network verifier and the dynamical system reachability tool, is _spatially aware_ in choosing the right locations of the state space to partition, and is _temporally aware_ to partition along trajectories rather than merely on the initial set. Using contraction theory for mixed monotone reachability analysis, we provide formal guarantees for the algorithm's performance. Finally, we run simulations to test the algorithm and show significant improvement in both runtime and accuracy as compared to existing partitioning approaches.
