Forward Invariance in Neural Network Controlled Systems
We present a framework based on interval analysis and monotone systems theory
to certify and search for forward invariant sets in nonlinear systems with
neural network controllers. The framework (i) constructs localized first-order
inclusion functions for the closed-loop system using Jacobian bounds and
existing neural network verification tools; (ii) builds a dynamical embedding
system where its evaluation along a single trajectory directly corresponds with
a nested family of hyper-rectangles provably converging to an attractive set of
the original system; (iii) utilizes linear transformations to build families of
nested paralleletopes with the same properties. The framework is automated in
Python using our interval analysis toolbox $\texttt{npinterval}$, in
conjunction with the symbolic arithmetic toolbox $\texttt{sympy}$, demonstrated
on an $8$-dimensional leader-follower system.
Akash Harapanahalli, Saber Jafarpour, Samuel Coogan
# Forward Invariance in Neural Network Controlled Systems

###### Abstract

We present a framework based on interval analysis and monotone systems theory to certify and search for forward invariant sets in nonlinear systems with neural network controllers. The framework (i) constructs localized first-order inclusion functions for the closed-loop system using Jacobian bounds and existing neural network verification tools; (ii) builds a dynamical embedding system where its evaluation along a single trajectory directly corresponds with a nested family of hyper-rectangles provably converging to an attractive set of the original system; (iii) utilizes linear transformations to build families of nested parallelotopes with the same properties. The framework is automated in Python using our interval analysis toolbox printerval, in conjunction with the symbolic arithmetic toolbox sympy, demonstrated on an \(\delta\)-dimensional leader-follower system.

## I Introduction

Learning enabled components are becoming increasingly prevalent in modern control systems. Their ease of computation and ability to outperform optimization-based approaches make them valuable for in-the-loop usage. However, neural networks are known to be vulnerable to input perturbations--small changes in their input can yield wildly varying results. In safety-critical applications, under uncertainty in the system, it is paramount to verify safe system behavior for an infinite time horizon. Such behaviors are guaranteed through _robust forward invariant sets_, _i.e._, a set for which a system will never leave under any uncertainty.

Forward invariant sets are useful for a variety of tasks. In monitoring, a safety specification extends to infinite time if the system is guaranteed to enter an invariant set. Additionally, for asymptotic behavior of a system, an invariant set can be used as a robustness margin to replace the traditional equilibrium viewpoint in the presence of disturbances. In designing controllers, one can induce infinite-time safe behavior by ensuring the existence of an invariant set containing all initial states of the system and excluding any unsafe regions. There are several classical techniques in the literature for certifying forward invariant sets, such as Lyapunov-based analysis [1], barrier-based methods [2], and set-based approaches [3]. However, a naive application of these methods generally fails when confronted with high-dimensional and highly nonlinear neural network controllers in-the-loop.

_Literature review:_ The problem of verifying the input-output behavior of standalone neural networks has been studied extensively [4]. There is a growing body of literature studying verification of neural networks applied in feedback loops, which presents unique challenges due to the accumulation of error in closed-loop, _i.e._, the wrapping effect. For example, there are functional approaches such as POLAR [5], JuliaReach [6], and ReachMM [7] for nonlinear systems, and linear (resp. semi-definite) programming ReachLP [8] (resp. Reach-SDP [9]) for linear systems. While these methods verify finite time safety, their guarantees do not readily extend to infinite time. In particular, it is not clear how to adapt these tools to search and certify forward invariant sets of neural network controlled systems. There are a handful of papers that directly study forward invariance for neural networks in dynamics: In [10] a set-based approach is used to study forward invariance of a specific class of control-affine systems with feedforward neural network controllers. In [11] an ellipsoidal inner-approximation of a region of attraction of neural network controlled system is obtained using Integral Quadratic Constraints (IQCs). In [12], a Lyapunov-based approach is used to study robust invariance of control systems modeled by neural networks.

_Contributions:_ In this letter, we propose a dynamical system approach for systematically finding nested families of robust forward invariant sets for nonlinear systems controlled by neural networks. Our method uses _localized_ first-order inclusion functions to construct an embedding system, which evaluates the inclusion function separately on the edges of a hyper-rectangle. Our first result is Proposition 1, which certifies (and fully characterizes for some systems) the forward invariance of a hyper-rectangle through the embedding system's evaluation at a single point. Our main result is Theorem 1, which describes how a single trajectory of the dynamical embedding system can be used to construct a nested family of invariant and attracting hyper-rectangles. However, in many applications, hyper-rectangles are not suitable for capturing forward invariant regions, and a simple linear transformation can greatly improve results. In Proposition 2, we carefully construct an accurate localized inclusion function for any linear transformation on the original system, which we use in Theorem 2 to find a nested family of forward invariant paralleletopes. Finally, we implement the framework in Python, demonstrating its applicability to an \(8\)-dimensional
###### Abstract

We consider a framework for bounding the behavior of the uncertain dynamical system (3) using inclusion functions. We first review the notion of inclusion functions in III-A and then use them to construct first-order local bounds for neural network controlled systems in III-B.

###### Contents

* 1 Introduction

[MISSING_PAGE_POST]

* 311
* 311
* 321
* 33 Proofs of Theorem 1.12
* 341
* 35 Proofs of Theorem 1.12
* 36
* 37 Proofs of Theorem 1.12
* 38 Proofs of Theorem 1.12
* 39 Proofs of Theorem 1.12
* 39 Proofs of Theorem 1.12
* 39 Proofs of Theorem 1.12
* 40
* 411
* 421
* 431
* 44

## 1 Introduction

### 1.1 Introduction

In this section, we introduce a framework for bounding the behavior of the uncertain dynamical system (3) using inclusion functions. We first review the notion of inclusion functions in III-A and then use them to construct first-order local bounds for neural network controlled systems in III-B.

### Inclusion Functions

Interval analysis aims to provide interval bounds on the output of a function given an interval of possible inputs [13]. Given a function \(f:\mathbb{R}^{n}\to\mathbb{R}^{m}\), the function \(\mathsf{F}=\left(\frac{\mathsf{F}}{\mathsf{F}}\right):\mathcal{T}_{\geq 0}^{2n} \to\mathcal{T}_{\geq 0}^{2m}\) is called an _inclusion function_ for \(f\) if

\[\mathsf{E}(\underline{x},\overline{x})\leq f(x)\leq\overline{\mathsf{F}}( \underline{x},\overline{x}),\quad\text{for every $x\in[\underline{x},\overline{x}]$}, \tag{4}\]

for every interval \([\underline{x},\overline{x}]\subset\mathbb{R}^{n}\), and is an \(\mathcal{S}\)-_localized inclusion function_ if the bounds (4) are valid for every interval \([\underline{x},\overline{x}]\subseteq\mathcal{S}\), in which case we instead write \(\mathsf{F}_{\mathcal{S}}\). Additionally, an inclusion function is

1. _monotone_ if \(\mathsf{F}(\underline{x},\overline{x})\geq_{\mathrm{SE}}\mathsf{F}(\underline{ y},\overline{y})\), for any \([\underline{x},\overline{x}]\subseteq[\underline{y},\overline{y}]\);
2. _thin_ if for any \(x\), we have \(\mathsf{E}(x,x)=f(x)=\overline{\mathsf{F}}(x,x)\);
3. _minimal_ if \(\mathsf{F}\) returns the tightest possible interval, _i.e._ for each \(i\in\{1,\ldots,m\}\), \[\mathsf{E}_{i}^{\min}(\underline{x},\overline{x})=\inf_{x\in[\underline{x}, \overline{x}]}f_{i}(x),\quad\overline{\mathsf{F}}_{i}^{\min}(\underline{x}, \overline{x})=\sup_{x\in[\underline{x},\overline{x}]}f_{i}(x).\]

Given the one-to-one correspondence between \(\mathcal{T}_{\geq 0}^{2n}\) and intervals of \(\mathbb{R}^{n}\), an inclusion function is often interpreted as a mapping from intervals to intervals--given an inclusion function \(\mathsf{F}=\left(\frac{\mathsf{F}}{\mathsf{F}}\right):\mathcal{T}_{\geq 0}^{2n} \to\mathcal{T}_{\geq 0}^{2m}\), we use the notation \([\mathsf{F}]=[\underline{\mathsf{F}},\overline{\mathsf{F}}]:\mathbb{R}^{n} \to\mathbb{R}^{m}\) to denote the equivalent interval-valued function with interval argument.

In this paper, we focus on two main methods to construct inclusion functions: (1) Given a composite function \(f=f_{1}\circ f_{2}\circ\cdots\circ f_{N}\), and inclusion functions \(\mathsf{F}_{i}\) for \(f_{i}\) for every \(i\in\{1,\ldots,N\}\), the _natural inclusion function_

\[\mathsf{F}^{\mathrm{nat}}(\underline{x},\overline{x}):=(\mathsf{F}_{1}\circ \mathsf{F}_{2}\cdots\circ\mathsf{F}_{N})(\underline{x},\overline{x}) \tag{5}\]

provides a simple but possibly conservative method to build inclusion functions using the inclusion functions of simpler mappings; (2) Given a differentiable function \(f\), with an inclusion function \(\mathsf{J}_{x}\) for its Jacobian matrix, the _jacobian-based inclusion function_

\[[\mathsf{F}^{\mathrm{jac}}(\underline{x},\overline{x})]:=[\mathsf{J}_{x}( \underline{x},\overline{x})]([\underline{x},\overline{x}]-\dot{x})+f(\dot{x}), \tag{6}\]
for any \(\hat{x}\in[\underline{x},\overline{x}]\), can provide better estimates by bounding the first order taylor expansion of \(f\). Both of these inclusion functions are monotone (resp. thin) assuming the inclusion functions used to build them are also monotone (resp. thin).

In our previous work [14], we introduce the open source package npinterval2, which extends the popular package numpy to include an interval data-type. The package automatically builds natural inclusion functions using minimal inclusion functions for standard functions and operations [14, Table 1], and can be used to construct jacobian-based inclusion functions when coupled with a symbolic toolbox like sympy.

Footnote 2: The code for npinterval can be found at [https://github.com/qftfactslab/npinterval](https://github.com/qftfactslab/npinterval).

### Localized Closed-Loop Inclusion Functions

One of the biggest challenges in neural network controlled system verification is correctly capturing the interactions between the system and the controller. For invariance analysis, it is paramount to capture the stabilizing nature of the controller, which can easily be lost with naive overbounding of the input-output interactions between the system and controller. We make the following assumption throughout.

**Assumption 1** (Local affine bounds of neural network).: For the neural network (2), there exists an algorithm that, for any interval \([\underline{\xi},\overline{\xi}]\), produces a local affine bound \((C_{[\underline{\xi},\overline{\xi}]},\underline{d}_{[\underline{\xi}, \overline{\xi}]},\overline{d}_{[\underline{\xi},\overline{\xi}]})\) such that for every \(x\in[\underline{\xi},\overline{\xi}]\),

\[C_{[\underline{\xi},\overline{\xi}]}x+\underline{d}_{[\underline{\xi}, \overline{\xi}]}\leq N(x)\leq C_{[\underline{\xi},\overline{\xi}]}x+ \overline{d}_{[\underline{\xi},\overline{\xi}]}.\]

Many off-the-shelf neural network verification frameworks can produce the linear estimates required in Assumption 1, and in particular, we focus on CROWN [15]. For ReLU and otherwise piecewise linear networks, one can setup a mixed integer linear program similar to [16], which is tractable for small-sized networks. The bounds from Assumption 1 can be used to construct a \([\underline{\xi},\overline{\xi}]\)-localized inclusion function for \(N(x)\):

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

where \((C^{+})_{i,j}=\max(C_{i,j},0)\), and \(C^{-}=C-C^{+}\).

Using the localized first-order bounds of the neural network, we propose a general framework for constructing closed-loop inclusion functions for neural network-controlled systems that capture the first-order stabilizing effects of the controller. First, assuming \(f\) is differentiable, with inclusion functions \(\mathsf{J}_{x},\mathsf{J}_{u},\mathsf{J}_{w}\) for its Jacobians, one can construct a closed-loop Jacobian-based inclusion function \(\mathsf{F}^{\mathsf{c}}\). With the Jacobian bounds evaluated on the input \((\underline{z},\overline{z},\underline{\mathbb{N}}_{[\underline{z}, \overline{z}]}(\underline{z},\overline{z}),\overline{\mathbb{N}}_{[\underline{z },\overline{z}]}(\underline{z},\overline{z}),\underline{w},\overline{w})\), define

\[[\mathsf{F}^{\mathsf{c}}_{[\underline{z},\overline{z}]}( \underline{x},\overline{x},\underline{w},\overline{w})]=([\mathsf{J}_{x}]+[ \mathsf{J}_{u}]C_{[\underline{x},\overline{\pi}]})[\underline{x},\overline{x}] \tag{8}\] \[+[\mathsf{J}_{u}][\underline{d}_{[\underline{x},\overline{x}]}, \overline{d}_{[\underline{x},\overline{x}]}]+[\mathsf{R}_{[\underline{z}, \overline{z}]}(\underline{w},\overline{w})],\]

where \(\hat{x}\in[\underline{x},\overline{x}]\subseteq[\underline{z},\overline{z}]\), \(\hat{u}=N(\hat{x})\), \(\hat{w}\in[\underline{w},\overline{w}]\), and \([\mathsf{R}_{[\underline{z},\overline{u}]}(\underline{w},\overline{w})]:=-[ \mathsf{J}_{x}]\hat{x}-[\mathsf{J}_{u}]\hat{u}+[\mathsf{J}_{w}]([\underline{w },\overline{w}]-\hat{w})+f(\hat{x},\hat{u},\hat{w})\). In Proposition 2, we provide a more general result which proves that (8) is a \([\underline{z},\overline{z}]\)-localized inclusion function for \(f^{\mathsf{c}}\) (take \(T=I\)).

In the case that \(f\) is not differentiable, or finding an inclusion function for its Jacobian is difficult, as long as a (monotone) inclusion function \(\mathsf{F}\) for the open-loop system \(f\) is known, a (monotone) closed-loop inclusion function for \(f^{\mathsf{c}}\) from (3) can be constructed using the natural inclusion approach in (5) with \(\mathsf{N}\) from (7)

\[\mathsf{F}^{\mathsf{c}}(\underline{x},\overline{x},\underline{w},\overline{w})= \mathsf{F}(\underline{x},\overline{\mathbb{N}}_{[\underline{x},\overline{x}]}( \underline{x},\overline{x}),\overline{\mathbb{N}}_{[\underline{x},\overline{ \pi}]}(\underline{x},\overline{x}),\underline{w},\overline{w}). \tag{9}\]

## 4 A Dynamical Approach to Set Invariance

Using the closed-loop inclusion functions developed in Section 3, we embed the uncertain dynamical system (3) into a larger certain system that enables computationally tractable approaches to verifying and computing families of invariant sets. Consider the closed-loop system (3) with an \(\mathcal{S}\)-localized inclusion function \(\mathsf{F}^{\mathsf{c}}_{S}:\mathcal{T}_{\geq 0}^{2n}\times\mathcal{T}_{\geq 0}^{2q} \rightarrow\mathcal{T}_{\geq 0}^{2n}\) for \(f^{\mathsf{c}}\) constructed via, _e.g._, (8) or (9), with the disturbance set \(\mathcal{W}\subseteq[\underline{w},\overline{w}]\). Then \(\mathsf{F}^{\mathsf{c}}_{S}\)_induces_ an _embedding system_ for (3) with state \((\frac{\overline{x}}{\overline{x}})\in\mathcal{T}_{\geq 0}^{2n}\) and dynamics defined by

\[\begin{split}\dot{\underline{x}}_{i}&=\Big{(} \mathsf{E}_{S}(\underline{x},\overline{x},\underline{w},\overline{w})\Big{)}_{i} :=\Big{(}\mathsf{E}_{S}^{\mathsf{c}}(\underline{x},\overline{x}_{i\underline{ x}},\underline{w},\overline{w})\Big{)}_{i},\\ \dot{\overline{x}}_{i}&=\Big{(}\mathsf{E}_{S}(\underline{x },\overline{x},\underline{w},\overline{w})\Big{)}_{i}:=\Big{(}\mathsf{F}_{S}^{ \mathsf{c}}(\underline{x}_{i\underline{x},\overline{x}},\underline{w}, \overline{w})\Big{)}_{i},\end{split} \tag{10}\]

where \(\mathsf{E}_{S}:\mathcal{T}_{\geq 0}^{2n}\times\mathcal{T}_{\geq 0}^{2q} \rightarrow\mathbb{R}^{2n}\). One of the key features of the embedding system, which evolves on \(\mathcal{T}_{\geq 0}^{2n}\), is that the inclusion function is evaluated separately on each face of the hyper-rectangle \([\underline{x},\overline{x}]\), represented by \([\underline{x},\overline{x}_{i\underline{x}}]\) and \([\underline{x}_{i\underline{x}},\overline{x}]\) for each \(i\in\{1,\ldots,n\}\). One favorable consequence of this separation is that the neural network verification steps from (8) and (9) are computed separately on each face, _i.e._, the \((C_{[\underline{x},\overline{x}]},\underline{d}_{[\underline{x},\overline{x}]}, \overline{d}_{[\underline{x},\overline{x}]})\) tuple from Assumption 1 is localized separately to each face. As seen in Proposition 1, this meshes nicely with Nagumo's Theorem [3], which allows us to guarantee forward invariance by checking the boundary of the hyper-rectangle through one evaluation of the embedding system (10).

**Proposition 1** (Forward invariance in hyper-rectangles).: _Consider the neural network controlled system (3) with the disturbance set \(\mathcal{W}=[\underline{w},\overline{w}]\). Let \(\mathsf{F}^{\mathsf{c}}_{S}\) be a \(\mathcal{S}\)-localized inclusion function for \(f^{\mathsf{c}}\), e.g. (8) or (9), and \(\mathsf{E}_{\mathcal{S}}\) be the embedding system (10) induced by \(\mathsf{F}^{\mathsf{c}}_{S}\). Consider some \([\underline{x}^{*},\overline{x}^{*}]\subseteq\mathcal{S}\). If_

\[\mathsf{E}_{S}(\underline{x}^{*},\overline{x}^{*},\underline{w},\overline{w}) \geq_{\mathrm{SE}}0, \tag{11}\]

_then \([\underline{x}^{*},\overline{x}^{*}]\) is a \([\underline{w},\overline{w}]\)-robustly forward invariant set. Moreover, if \(\mathsf{F}^{\mathsf{c}}_{S}\) is the minimal inclusion function, the condition (11) is also necessary for \([\underline{x}^{*},\overline{x}^{*}]\) to be a \([\underline{w},\overline{w}]\)-robustly forward invariant set._

Proof.: For brevity, since \([\underline{x}^{*},\overline{x}^{*}]\subseteq\mathcal{S}\), we drop \(\mathcal{S}\) from the notation. Consider the set \([\underline{x}^{*},\overline{x}^{*}]\), and suppose that \(\mathsf{E}(\underline{x}^
on the \(i\)-th upper face of the hyperrectangle \([\underline{x}^{\star}_{i;\overline{x}^{\star}},\overline{x}^{\star}]\). Since this holds for every \(i\in\{1,\ldots,n\}\), by Nagumo's theorem [3, Theorem 3.1], the closed set \([\underline{x}^{\star},\overline{x}^{\star}]\) is forward invariant since for every point \(x\) along its boundary \(\bigcup_{i}([\underline{x}^{\star},\overline{x}^{\star}_{i;\underline{x}^{ \star}}]\cup[\underline{x}^{\star}_{i;\overline{x}^{\star}},\overline{x}^{ \star}])\), the vector field \(f^{\mathrm{c}}(x,w)\) points into the set, for every \(w\in[\underline{w},\overline{w}]\). Now, suppose that \(\mathsf{E}\) is the embedding system induced by the minimal inclusion function \(\mathsf{F}^{\min}\) for \(f^{\mathrm{c}}\), and \([\underline{x}^{\star},\overline{x}^{\star}]\) is a hyper-rectangle such that \(\mathsf{E}(\underline{x}^{\star},\overline{x}^{\star},\underline{w},\overline{ w})\not\geq_{\mathrm{SE}}0\). Then there exists \(i\in\{1,\ldots,n\}\) such that either

\[\mathsf{E}^{\min}_{i}(\underline{x}^{\star},\overline{x}^{\star} _{i;\underline{x}^{\star}},\underline{w},\overline{w}) =\inf_{x\in[\underline{x}^{\star},\overline{x}^{\star}_{i; \underline{x}^{\star}}],w\in[\underline{w},\overline{w}]}f(x,w)<0,\text{ or}\] \[\mathsf{F}^{\min}_{i}(\underline{x}^{\star}_{i;\overline{x}^{ \star}},\overline{x}^{\star},\underline{w},\overline{w}) =\sup_{x\in[\underline{x}^{\star}_{i;\overline{x}^{\star}}, \overline{x}^{\star}],w\in[\underline{w},\overline{w}]}f(x,w)>0.\]

If the first case holds, then there exists \(x^{\prime}\in[\underline{x}^{\star},\overline{x}^{\star}_{i;\underline{x}}],w \in[\underline{w},\overline{w}]\) such that \(f^{\mathrm{c}}(x^{\prime},w)<0\) along the \(i\)-th lower face of the hyper-rectangle. If the second case holds, then there exists \(x^{\prime}\in[\underline{x}^{\star}_{i;\overline{x}^{\star}},\overline{x}^{ \star}],w\in[\underline{w},\overline{w}]\) such that \(f^{\mathrm{c}}(x^{\prime},w)>0\) along the \(i\)-th upper face of the hyper-rectangle. Thus, by Nagumo's theorem, the set \([\underline{x}^{\star},\overline{x}^{\star}]\) is not \([\underline{w},\overline{w}]\)-robustly forward invariant, as there exists a point along its boundary such that the vector field \(f^{\mathrm{c}}\) points outside the set. 

**Remark 1** (Linear systems).: For the special case of the linear system \(\dot{x}=Ax+Bu+Dw\) controlled by a neural network \(u=N(x)\) with piecewise linear activations, ReLU or Leaky ReLU, one can compute the minimal inclusion function using a Mixed Integer Linear Program (MILP) similar to the formulation from [16] as

\[\mathsf{E}^{\min}_{i} =\min_{x\in[\underline{x},\overline{x}],w\in[\underline{w}, \overline{w}]}(Ax+BN(x)+Dw)_{i},\] \[\mathsf{F}^{\min}_{i} =\max_{x\in[\underline{x},\overline{x}],w\in[\underline{w}, \overline{w}]}(Ax+BN(x)+Dw)_{i},\]

for \(i\in\{1,\ldots,n\}\). When using the minimal inclusion function, Proposition 1 can be used to additionally guarantee that a hyper-rectangle is _not_ forward invariant.

The next Theorem shows how monotonicity of the dynamical embedding system can be used to define a family of nested robustly forward invariant sets using the condition from Proposition 1.

**Theorem 1** (A nested family of invariant sets).: _Consider the neural network controlled system (3) with the disturbance set \(\mathcal{W}=[\underline{w},\overline{w}]\). Let \(\mathsf{F}^{\mathrm{c}}_{\mathcal{S}}\) be a \(\mathcal{S}\)-localized monotone inclusion function for \(f^{\mathrm{c}}\), e.g. (8) or (9), and \(\mathsf{E}_{\mathrm{S}}\) be the embedding system induced by \(\mathsf{F}^{\mathrm{c}}_{\mathcal{S}}\). Consider \([\underline{x}_{0},\overline{x}_{0}]\subseteq\mathcal{S}\) and let \(t\mapsto\left(\frac{\underline{x}^{(t)}}{\overline{x}(t)}\right)\) be the trajectory of (10) with \(\left(\frac{\underline{x}^{(0)}}{\overline{x}_{0}}\right)=\left(\frac{\underline {x}_{0}}{\overline{x}_{0}}\right)\). If_

\[\mathsf{E}_{\mathcal{S}}(\underline{x}_{0},\overline{x}_{0},\underline{w}, \overline{w})\geq_{\mathrm{SE}}0,\]

_then the following statements hold:_

1. \([\underline{x}(t),\overline{x}(t)]\) _is a_ \([\underline{w},\overline{w}]\)_-robustly forward invariant set for the system (_3_) for every_ \(t\geq 0\)_, and_ \[[\underline{x}(\tau),\overline{x}(\tau)]\subseteq[\underline{x}(t),\overline{ x}(t)],\ \ \text{for every}\ t\leq\tau;\]
2. \(\lim_{t\to\infty}\left(\frac{\underline{x}^{(t)}}{\overline{x}(t)}\right)= \left(\frac{\underline{x}^{\star}}{\overline{x}^{\star}}\right)\)_, where_ \(\left(\frac{\underline{x}^{\star}}{\overline{x}^{\star}}\right)\in\mathcal{T}_{ \geq 0}^{2n}\) _is an equilibrium point of the embedding system (_10_) and_ \([\underline{x}^{\star},\overline{x}^{\star}]\) _is a_ \([\underline{w},\overline{w}]\)_-attracting set for the system (_3_) with region of attraction_ \([\underline{x}(t^{\prime}),\overline{x}(t^{\prime})]\)_._

Proof.: Consider \((\frac{x}{\overline{x}})\leq_{\mathrm{SE}}(\frac{\overline{x}}{\overline{x}})\). This implies that \(\left(\frac{\underline{x}}{\overline{x}_{i;\underline{x}}}\right)\leq_{ \mathrm{SE}}\left(\frac{\underline{x}}{\overline{x}_{i;\underline{x}}}\right)\) and \(\left(\frac{\underline{x}}{\overline{x}}\right)\approx_{\mathrm{SE}}\left( \frac{\underline{x}_{i;\overline{x}}}{\overline{x}^{\star}}\right)\). Since \(\mathsf{F}^{\mathrm{c}}_{\mathcal{S}}\) is monotone,

\[\left(\mathsf{E}^{\mathrm{c}}_{\mathcal{S}}(\underline{x}, \overline{x}_{i;\underline{x}},\underline{w},\overline{w})\right)_{i} \leq\left(\mathsf{E}^{\mathrm{c}}_{\mathcal{S}}(\underline{z},\overline{z}_{i; \underline{x}},\underline{w},\overline{w})\right)_{i},\] \[\left(\overline{\mathsf{F}}^{\mathrm{c}}_{\mathcal{S}}(\underline{z}, \overline{x},\overline{w},\overline{w})\right)_{i} \geq\left(\overline{\mathsf{F}}^{\mathrm{c}}_{\mathcal{S}}(\underline{z}, \overline{z},\overline{w},\overline{w})\right)_{i},\]

for every \(i\in\{1,\ldots,n\}\). This implies that the embedding system \(\mathsf{E}^{\mathrm{c}}_{\mathcal{S}}\) is monotone w.r.t. the southeast partial order \(\leq_{\mathrm{SE}}\)[17, 18]. Now using [19, Proposition 2.1], the set \(\mathcal{P}_{+}=\{(\frac{\overline{x}}{\overline{x}}):(\frac{\overline{x}}{ \overline{x}})\geq_{\mathrm{SE}}0\}\) is a forward invariant set for \(\mathsf{E}^{\mathrm{c}}_{\mathcal{S}}\). Since \(\left(\frac{\overline{x}}{\overline{x}_{0}}\right)\in\mathcal{P}_{+}\), forward invariance implies \(\left(\frac{\underline{x}^{(t)}}{\overline{x}(t)}\right)\in\mathcal{P}_{+}\) for every \(t\geq 0\). Therefore, using Proposition 1, \([\underline{x}(t),\overline{x}(t)]\) is forward invariant for the closed-loop system (3) for every \(t\geq 0\). Additionally, the curve \(t\mapsto\left(\frac{\underline{x}^{(t)}}{\overline{x}(t)}\right)\) is nondecreasing with respect to the partial order \(\leq_{\mathrm{SE}}\)[19, Proposition 2.1]. This means that for every \(t\leq\tau\), \(\left(\frac{\underline{x}^{(t)}}{\overline{x}(t)}\right)\leq_{\mathrm{SE}}\left( \frac{\underline{x}^{(\tau)}}{\overline{x}(\tau)}\right)\), which implies that \([\underline{x}(\tau),\overline{x}(\tau)]\subseteq[\underline{x}(t),\overline{x}(t)]\). Regarding part (ii), since \(t\mapsto\left(\frac{\underline{x}^{(t)}}{\overline{x}(t)}\right)\) is nondecreasing w.r.t. \(\leq_{\mathrm{SE}}\), for every \(i\in\{1,\ldots,n\}\), the curves \(t\mapsto\underline{x}_{i}(t)\) (resp. \(t\mapsto\overline{x}_{i}(t)\)) are nondecreasing (resp. nonincreasing) w.r.t. \(\leq\) and bounded on \(\mathbb{R}\). This implies that there exists \(\
where \(N^{\prime}(y):=N(T^{-1}y)\) is the neural network from (2), with an extra initial layer with weight matrix \(T^{-1}\) and linear activation \(\phi(x)=x\). There is a one-to-one correspondence between the transformed system (12) and the original system (3), in the sense that every trajectory \(t\mapsto y(t)\) of (12) uniquely corresponds with the trajectory \(t\mapsto\Phi^{-1}(y(t))\) of (3). Given an interval \([\underline{y},\overline{y}]\), the set \(\Phi^{-1}([\underline{y},\overline{y}])=\{T^{-1}y:y\in[\underline{y}, \overline{y}]\}\) defines a paralleletope in standard coordinates.

Using the definitions from (8), we construct a localized closed-loop Jacobian-based inclusion function \(\mathsf{G}^{\mathsf{c}}_{\Phi([\underline{x},\overline{y}])}\) as

\[[\mathsf{G}^{\mathsf{c}}_{\Phi([\underline{x},\overline{y}])}( \underline{y},\overline{y},\underline{w},\overline{w})] =T([\mathsf{J}_{x}]+[\mathsf{J}_{u}](C^{\prime}_{[\underline{y}, \overline{y}]})T)^{-1}[\underline{y},\overline{y}]\] \[\quad+T[\mathsf{J}_{u}][\underline{d}^{\prime}_{[\underline{y}, \overline{y}]},\overline{d}^{\prime}_{[\underline{y},\overline{y}]}]+T[ \mathsf{R}_{[\underline{x},\overline{y}]}(\underline{w},\overline{w})], \tag{13}\]

where \([\mathsf{J}_{\mathbf{s}}]\) are inclusion functions on the Jacobians of the original system evaluated on \((\underline{z},\overline{z},\overline{\mathsf{N}}_{[\underline{z},\overline{ z}]}(\underline{z},\overline{z}),\overline{\mathsf{N}}_{[\underline{z},\overline{z}]}( \underline{z},\underline{z}),\underline{w},\overline{w})\) for \(\mathbf{s}\in\{x,u,w\}\), \((C^{\prime},\underline{d}^{\prime},\overline{d}^{\prime})\) evaluated over \([\underline{y},\overline{y}]\) on the transformed neural network \(N^{\prime}(y)=N(T^{-1}y)\), \(\hat{x}\in[\underline{z},\overline{z}]\), \(\hat{\bar{u}}=N(\hat{x})\), and \(\hat{\bar{w}}\in[\underline{w},\overline{w}]\).

**Proposition 2**.: _Consider the neural network controlled system (3) with (monotone) inclusion function (8) and let \(T\in\mathbb{R}^{n\times n}\) be an invertible matrix transforming the system into (12). Then (13) is a \(\Phi([\underline{z},\overline{z}])\)-localized (monotone) inclusion function for \(g^{\mathsf{c}}\)._

Proof.: Using the mean value theorem with the inclusion functions \(\mathsf{J}_{\mathbf{s}}\) for \(\mathbf{s}\in\{x,u,w\}\) valid on the input \((\hat{\underline{z}},\overline{z},\underline{\mathsf{N}}_{[\underline{z}, \overline{z}]}(\underline{z},\overline{z}),\overline{\mathsf{N}}_{[\underline{ z},\overline{z}]}(\underline{z},\underline{z}),\underline{w},\overline{w})\),

\[g^{\mathsf{c}}(y,w) \in T[\mathsf{J}_{x}](T^{-1}y-\hat{\bar{x}})+T[\mathsf{J}_{u}](N( T^{-1}y)-\hat{\bar{u}})\] \[\quad+T[\mathsf{J}_{u}](w-\hat{\bar{w}})+Tf(\hat{\underline{z}}, \hat{\bar{u}},\hat{\bar{w}}).\]

Using the tuple \((C^{\prime}_{[\underline{y},\overline{y}]},\underline{d}^{\prime}_{[ \underline{y},\overline{y}]},\overline{d}^{\prime}_{[\underline{y},\overline{ y}]})\) from Assumption 1,

\[g^{\mathsf{c}}(x,w) \in T([\mathsf{J}_{x}]+[\mathsf{J}_{u}](C^{\prime}_{[\underline{y}, \overline{y}]}T))T^{-1}[\underline{y},\overline{y}]\] \[\quad+T[\mathsf{J}_{u}][\underline{d}^{\prime}_{[\underline{y}, \overline{y}]},\overline{d}^{\prime}_{[\underline{y},\overline{y}]}]+T[ \mathsf{R}_{[\underline{x},\overline{z}]}(\underline{w},\overline{w})],\]

implying that (13) is an inclusion function for \(g^{\mathsf{c}}\). Regarding monotonicity, if the Jacobian inclusion functions are monotone, since addition and matrix multiplication are also monotone, (8) is a composition of monotone inclusion functions and is therefore monotone. 

It is important to note that the inclusion functions for the Jacobian matrices in (13) are obtained in the original coordinates. In principle, one could instead symbolically write \(g\) as a new system and directly apply the closed-loop Jacobian-based inclusion function from (8). In practice, however, these transformed dynamics often have complicated expressions that lead to excessive conservatism when using natural inclusion functions and are not suitable for characterizing invariant sets. In the next Theorem, we link forward invariant hyperrectangles in transformed coordinates to forward invariant paralleletopes in standard coordinates.

**Theorem 2** (Forward invariance in paralleletopes).: _Consider the neural network controlled system (3) and let \(T\in\mathbb{R}^{n\times n}\) be an invertible matrix transforming the system into (12). Let \(\mathsf{G}^{\mathsf{c}}_{\mathcal{S}}\) be a \(\mathcal{S}\)-localized inclusion function for \(g^{\mathsf{c}}\), e.g. (13), and let \(\mathsf{E}_{T,\mathcal{S}}\) be the embedding system (10) induced by \(\mathsf{G}^{\mathsf{c}}_{\mathcal{S}}\) and \([\underline{y}_{0},\overline{y}_{0}]\subseteq\mathcal{S}\). If_

\[\mathsf{E}_{T,\mathcal{S}}(\underline{y}_{0},\overline{y}_{0},\underline{w}, \overline{w})\geq_{\mathrm{SE}}0,\]

_then the paralleletope \(\Phi^{-1}([\underline{y}_{0},\overline{y}_{0}])\) is a \([\underline{w},\overline{w}]\)-robustly forward invariant set for the neural network controlled system (3)._

**Remark 3**.:
1. _(Evaluation on faces):_ When using (13) as \(\mathsf{G}^{\mathsf{c}}_{\mathcal{S}}\), the neural network verification step from Assumption 1 to find \((C^{\prime},\underline{d}^{\prime},\overline{d}^{\prime})\) is evaluated separately on each face of the \([\underline{y},\overline{y}]\) hyperrectangle, corresponding to each face of the paralleletope \(\Phi^{-1}([\underline{y},\overline{y}])\).
2. _(Nested invariant paralleletopes):_ Using the dynamical approach from Theorem 1, one can obtain a nested family of invariant hyper-rectangles for the transformed system 12. Using Theorem 2, this nested family corresponds to a nested family of paralleletope invariant sets for the original system (3).
3. _(Choice of transformation):_ In practice, there are principled ways of choosing good transformations \(T\). One approach is to take the transformation associated with the Jordan decomposition of the linearization of the system at an equilibrium point, as is done in the example below.

## 6 Numerical Experiments

Consider two vehicles \(\mathrm{L}\) and \(\mathrm{F}\) each with dynamics

\[\begin{split}\dot{p}^{j}_{x}=v^{j}_{x},&\dot{v}^{j} _{x}=\sigma(u^{j}_{x})+w^{j}_{x},\\ \dot{p}^{j}_{y}=v^{j}_{y},&\dot{v}^{j}_{y}=\sigma(u^{j }_{y})+w^{j}_{y},\end{split} \tag{14}\]

for \(j\in\{\mathrm{L},\mathrm{F}\}\), where \(p^{j}=(p^{j}_{x},p^{j}_{y})\in\mathbb{R}^{2}\) is the displacement of the center of mass of \(j\) in the plane, \(v^{j}=(v^{j}_{x},v^{j}_{y})\in\mathbb{R}^{2}\) is the velocity of the center of mass of \(j\), \((u^{j}_{x},u^{j}_{y})\in\mathbb{R}^{2}\) are desired acceleration inputs limited by the nonlinear softmax operator \(\sigma(u)=u_{\text{lim}}\tanh(u/u_{\text{lim}})\) with \(u_{\text{lim}}=20\), and \(w^{j}_{x},w^{j}_{y}\sim\mathcal{U}([-0.01,0.01])\) are uniformly distributed bounded disturbances on \(j\). Denote the combined state of the system \(x=(p^{\mathrm{L}},v^{\mathrm{L}},p^{\mathrm{F}},v^{\mathrm{F}})\in\mathbb{R}^{8}\). We consider a leader-follower structure for the system, where the leader vehicle \(\mathrm{L}\) chooses its control \(u=(u^{\mathrm{L}}_{x},u^{\mathrm{L}}_{y})\) as the output of a neural network \((4\times 10\times 100\times 2\), ReLU activations), and the follower vehicle \(\mathrm{F}\) follows the leader with a \(\mathrm{D}\) controller

\[u^{\mathrm{F}}_{\mathbf{d}}=k_{p}(p^{\mathrm{L}}_{\mathbf{d}}-p^{\mathrm{F}}_{ \mathbf{d}})+k_{v}(v^{\mathrm{L}}_{\mathbf{d}}-v^{\mathrm{F}}_{\mathbf{d}}), \tag{15}\]

for each \(\mathbf{d}\in\{x,y\}\) with \(k_{p}=6\) and \(k_{v}=7\). The neural network was trained by gathering \(5.7\)M data points from an offline MPC control policy for the leader only, with control limits implemented as hard constraints rather than \(\sigma\). The offline policy minimized a quadratic cost aiming to stabilize to the origin while avoiding four circular obstacles centered at \((\pm 4,\pm 4)\) with radius \(2.25\) each, implemented as hard constraints with \(33\%\) padding and a slack variable.

First, a trajectory of the undisturbed system is run until it reaches the equilibrium point \(x^{\star}\approx[0.01,0,0,0,0.01,0,0,0]^{T}\). Then, using sympy, the Jacobian matrices
using auto_LiRPA [21] along the interval \([\underline{z},\overline{z}]:=x^{\star}+[-0.06,0.06]^{4}\times[-0.25,0.25]^{2} \times[-0.325,0.325]^{2}\) yielding \((C_{[\underline{z},\overline{z}]},\underline{d}_{[\underline{z},\overline{z} ]},\overline{d}_{[\underline{z},\overline{z}]})\). Then, the Jordan decomposition \(T^{-1}JT=(Df_{x}(x^{\star},N(x^{\star}),0)+Df_{x}(x^{\star},N(x^{\star}),0)C_{[ \underline{z},\overline{z}]})\) yields the transformation \(T\), and the matrix \(J\approx\mathrm{diag}(-6,-6,-4.12,-4.26,-0.93,-0.95,-1,-1)\) is filled with negative real eigenvalues, signifying that the equilibrium \(x^{\star}\) is locally stable. The \(T\)-transformed system (12) is analyzed with the embedding system induced by (13). The interval \([\underline{y}_{0},\overline{y}_{0}]=Tx^{\star}+[-0.1,0.1]^{6}\times[-0.15,0. 15]\times[-0.2,0.2]\) yields \(\Phi^{-1}([\underline{y}_{0},\overline{y}_{0}])\subseteq[\underline{z}, \overline{z}]\), and \(E_{T}(\underline{y}_{0},\overline{y}_{0},\underline{w},\overline{w})\geq \mathrm{se}\ 0\). Thus, using Theorem 2, the paralleletope \(\Phi^{-1}([\underline{y}_{0},\overline{y}_{0}])\) is an invariant set of (14). The embedding system in \(y\) coordinates is simulated forwards using Euler integration with a step-size of \(0.1\) for \(90\) time steps, and at each step the localization \([\underline{z},\overline{z}]=T^{-1}[\underline{y}(t),\overline{y}(t)]\) is refined. Starting from \(\left(\frac{y_{0}}{\overline{y}_{0}}\right)\), the forward-time embedding system converges to a point \(\left(\frac{y^{\star}}{\overline{y}}\right)\), where \(E_{T}(\underline{y}^{\star},\overline{y}^{\star},\underline{w},\overline{w})=0\). The transformed embedding system is also simulated backwards using Euler integration with a step-size of \(0.05\) while the condition \(\left(\frac{y^{\star}(t)}{\overline{y}(t)}\right)\geq_{\mathrm{SE}}0\) is satisfied (call the final time \(t^{\prime}\)). Using Theorems 1 and 2, the collection \(\{\Phi^{-1}([\underline{y}(t),\overline{y}(t)])\}_{t\geq t^{\prime}}\) consists of nested invariant paralleletopes converging to the attractive set \(\Phi^{-1}([\underline{y}^{\star},\overline{y}^{\star}])\) with region of attraction \(\Phi^{-1}([\underline{y}(t^{\prime}),\overline{y}(t^{\prime})])\). The initial paralleletope takes \(0.38\) seconds to verify, and the entire nested family takes \(40.28\) seconds to compute \(95\) paralleletopes. These results are visualized in Figure 1.

## VII Conclusions

Using interval analysis and neural network verification tools, we propose a framework for certifying hyper-rectangle and paralleletope invariant sets in neural network controlled systems. The key component of our approach is the dynamical embedding system, whose trajectories can be used to construct a nested family of invariant sets. This work opens up an avenue for future work in designing safe learning-enabled controllers.
