Verification of Neural Network Control Systems using Symbolic Zonotopes
  and Polynotopes
Verification and safety assessment of neural network controlled systems
(NNCSs) is an emerging challenge. To provide guarantees, verification tools
must efficiently capture the interplay between the neural network and the
physical system within the control loop. In this paper, a compositional
approach focused on inclusion preserving long term symbolic dependency modeling
is proposed for the analysis of NNCSs. First of all, the matrix structure of
symbolic zonotopes is exploited to efficiently abstract the input/output
mapping of the loop elements through (inclusion preserving) affine symbolic
expressions, thus maintaining linear dependencies between interacting blocks.
Then, two further extensions are studied. Firstly, symbolic polynotopes are
used to abstract the loop elements behaviour by means of polynomial symbolic
expressions and dependencies. Secondly, an original input partitioning
algorithm takes advantage of symbol preservation to assess the sensitivity of
the computed approximation to some input directions. The approach is evaluated
via different numerical examples and benchmarks. A good trade-off between low
conservatism and computational efficiency is obtained.
Carlos Trapiello, Christophe Combastel, Ali Zolghadri
# Verification of Neural Network Control Systems using Symbolic Zonotopes and Polynotopes

###### Abstract

Verification and safety assessment of neural network controlled systems (NNCSs) is an emerging challenge. To provide guarantees, verification tools must efficiently capture the interplay between the neural network and the physical system within the control loop. In this paper, a compositional approach focused on inclusion preserving long term symbolic dependency modeling is proposed for the analysis of NNCSs. First of all, the matrix structure of symbolic zonotopes is exploited to efficiently abstract the input/output mapping of the loop elements through (inclusion preserving) affine symbolic expressions, thus maintaining linear dependencies between interacting blocks. Then, two further extensions are studied. Firstly, symbolic polynotopes are used to abstract the loop elements behaviour by means of polynomial symbolic expressions and dependencies. Secondly, an original input partitioning algorithm takes advantage of symbol preservation to assess the sensitivity of the computed approximation to some input directions. The approach is evaluated via different numerical examples and benchmarks. A good trade-off between low conservatism and computational efficiency is obtained.

 Reachability, neural networks, verification, symbolic zonotopes, polynotopes, nonlinear dynamics.

## I Introduction

The proliferation of data and access to ever-increasing computational power have fueled a renewed interest in deep neural-networks (NNs). These networks have shown a significant ability to address classification/estimation/control tasks that can hardly be formalized and designed from knowledge-based models. However, despite their impressive ability for solving complex problems, it is well known that NNs can be vulnerable to small perturbations or adversarial attacks [1, 2]. This lack of robustness (or fragility) represents a major barrier for their application to safety-critical system where safety assurances are of primary importance. For example, in Guidance, Navigation and Control of flight systems, one must ensure that some output/state trajectories remain inside a flight envelope when some inputs explore a given region. The above issues have fostered a large amount of works that analyze the sensitivity to local disturbances of NNs in isolation (open-loop), as well as the satisfaction of pre-/post- safety conditions [3]. Nevertheless, as reported in [4], reasoning about the safety verification of neural-network control systems (NNCSs), where the NN is used as a feedback controller, still remains a key challenge that requires tractable methods capable of efficiently integrating the heterogeneous components that make up the control loop.

This paper focuses on the reachability analysis of NNCSs which, in turn, allows for formal reasoning about the satisfaction of safety properties (reachability of a target set, or the avoidance of non-secure sets of states). A key challenge in NNCSs reachability analysis is to successfully retain the system-controller interplay by preserving (at each time instant) the dependencies between relevant variables. This, in fact, discourages a direct application of off-the-shelf verification tools which, although able to compute accurate output bounds for elements in isolation, return coarse approximations when iteratively concatenated for the analysis of closed-loop systems since most of (if not all) the I/O dependencies are quickly broken/lost during the computations [5, 6, 7]. Furthermore, effective NNCSs verification tools must be able to assess the system state during (relatively) large time intervals. The above issues motivate the development of computationally efficient analysis methods capable of capturing the interaction between the control loop elements while granting a good scalability both in the system dimensions and in the time horizon length.

Another relevant factor that should be taken into account is the size of the initial state set under study. Mainly, the performance of open- and closed-loop verification techniques that are based on (locally) abstracting the system non-linearities, deteriorates considerably for large initial sets. A common approach to address this issue, particularly in NNCSs verification problems where the number of dimensions is relatively small, is to recast the initial reachability problem into simpler subproblems that analyze a subset of the initial conditions [5, 7, 8, 9, 10, 11, 12, 13]. Nonetheless, the design of efficient and scalable partitioning strategies, specially in closed-loop verification schemes, remains also an open problem.

Related workPreserving dependencies for NNCS verification has spurred on some recent studies. In [14], the authors abstract the I/O mapping of a ReLU NN controller using a polynomial expression (plus an error interval). The polynomial rule is obtained by regression of I/O samples, whereas a sound error term is derived from solving a mixed-integer program (MIP). In a similar fashion, [15] uses Bernstein polynomials to abstract the NN controller. A theoretical and a sampled-based method is proposed to compute the error term based on the Lipschitz constant of the NN. Although both approaches preserve the system-controller interplay, they are computationally expensive, scaling poorly with the number of NN inputs while requiring to be iteratively repeated for each
output. In [9] a NN with differentiable activation functions is transformed into an equivalent hybrid system built upon Taylor models that retain dependencies. However, this approach is not applicable to ReLU functions and the number of states (resp. modes) of the hybrid automaton scales with the number of neurons (resp. layers).

Other approaches preserve system-controller dependencies by formulating the reachable set computation as an optimization problem. The work [16] proposes a semidefinite program for reachability analysis based on the abstraction of the NN non-linearities using quadratic constraints [17]. [11] relies on the tool CROWN [18], and preserves system-controller interaction by solving LP-programs. However, dependencies are broken from one sample to the next. In [10], the closed-loop is firstly abstracted as a conjunction of piecewise-linear functions, and then analyzed using ReLU NNs verification tools like [19], [20].

On the other hand, other works address the reachability problem by chaining different verification tools. In [5], the authors combine a polytopic abstraction of the dynamical system with the tool Sherlock [21], that is used to bound the NN controller outputs. The tool NNV [6], integrates the non-linear dynamics reachability tool CORA [22] with a star sets abstraction of ReLU NN controllers [23]. Besides, [7] combines validated simulation to soundly approximate the dynamical system with common tools for NN output bounding like DeepPoly [24]. In all the above works, dependencies are broken in the switch between the different tools. This latter issue is somehow palliated in [8], where the authors use second order zonotopes (i.e. zonotopes with generators matrix size \(n\times 2n\)) as an interface between system and NN controller analysis tools. Although capable to retain first order dependencies in the system-to-controller (and controller-to-system) set transformations, dependencies in the I/O mapping of the NN controller are broken.

Focusing on partitioning strategies, in [25] the gradient of a ReLU NN (open-loop) is used to decide the next input direction to be bisected, whereas in [26] a uniform grid of the initial set is employed. Other works propose a simulation-based splitting strategy. In [12], the bisection is guided by comparing the interval bound of Monte-Carlo samples with a guaranteed Interval Bound Propagation [27] of the initial subsets. Working in a similar fashion, [13] proposes a simulation-guided framework that unifies standard NN output bounding tools. The decision on the bisection order is based on the distance to the simulation samples enclosure. A closed-loop implementation of the latter algorithm is reported in [11].

ContributionsThis paper takes a new and original direction based on symbolic zonotopes (s-zonotopes) as a generic tool for the closed-loop verification of discrete-time NN controlled systems. The generators (matrix) representation of s-zonotopes enables to efficiently abstract the input-output mappings of the NN controller and non-linear physical system through (inclusion preserving) affine symbolic expressions. The evolution of the closed-loop system can then be bounded in a propagation-based fashion that benefits from the efficient computation of basic operations granted by s-zonotopes, while preserving system-controller linear dependencies. Besides, the computational complexity of the verification tool can be fixed by limiting (reducing) the number of independent symbols. Simulations show the good performance/computational efficiency trade-off granted by this approach.

Furthermore, two extensions are proposed. On the one hand, the use of polynomial symbolic expressions to abstract the input-output mapping of the loop elements is explored. In particular, symbolic polynotope (s-polynotope) structures [28] are used to enclose the NN activation functions graph via the non-convex sets that arise from the polynomial map of interval symbols. Polynomial abstractions enable to reduce the conservatism induced by linear relaxations, at the price of increasing the computation needs.

On the other hand, the symbols preservation throughout the control loop is exploited to develop a smart partitioning strategy of the initial conditions set. The proposed algorithm reasons upon the influence of the input symbols in the output set in order to select which dimension to bisect next, and upon the influence of the (independent) error symbols to assess the quality of each over-approximation.

StructureThe paper is organized as follows. Section II is devoted to some useful preliminaries. Section III introduces the problem statement. Then, Section IV analyzes the closed-loop verification using s-zonotopes. In Section V the use of s-polynotopes is investigated, whereas the input partitioning algorithm is detailed in Section VI. Section VII presents simulation results. Finally, some concluding remarks are provided in Section VIII.

NotationThe following notations are used along this work. \(\mathbb{R}^{n}\), \(\mathbb{R}^{m\times n}\) and \(\mathbb{N}\) denote the \(n\) dimension Euclidean space, the \(m\times n\) dimensional Euclidean space and the set of non-negative integers, respectively. The notation \(v_{i}\) stands for the \(i\)-th element of vector \(v\) and \(M_{[i,\cdot]}\) (\(M_{[\cdot,j]}\)) for the \(i\)-th row (\(j\)-th column) of matrix \(M\). The 1-norm of the (row) vector \(v\) is \(\|v\|_{1}=|v|\textbf{1}\), with \(|.|\) the elementwise absolute value, and \(\mathbf{1}\) a column vector of ones of appropriate size. \(diag(v)\) returns a square diagonal matrix with the elements of vector \(v\) in the main diagonal, whereas \(card(\cdot)\) gives the cardinal.

## II Symbolic dependencies in set computations

This section provides preliminary concepts which will be used in the following sections. Throughout this article, \(s\) refers to an indexed family of distinct symbolic variables of type unit interval, that is, \(\forall i\in\mathbb{N}\), the symbol \(s_{i}\) (uniquely identified by the integer \(i\)) refers to a scalar real variable the value of which is only known to belong to the unit interval \(\mathcal{D}(s_{i})=[-1,+1]\subset\mathbb{R}\). Also, \(\mathcal{D}(s)=[-1,+1]^{card(s)}\). In other words, the a priori unknown value \(\iota s_{i}\) taken by the symbolic variable \(s_{i}\) satisfies \(\iota s_{i}\in\mathcal{D}(s_{i})\). The generic notation \(\iota\) which reads as "interpretation/valuation of" helps disambiguate between symbols (syntax) and values (semantics) [28]. Note that, in general, several interpretations may coexist. Set-valued interpretations take sets as values. In the following, consistently with the definition domain \(\mathcal{D}(s_{i})\) related to \(s_{i}\), the set-valued interpretation of each symbolic variable \(s_{i}\) will be \(s_{i|_{\iota}}=[-1,+1]\). In addition, the integer-valued vector \(I\) is used to uniquely identify a set of symbols, for example, vector
\(I=[1,\,5,\,3]\) identifies the symbols \(s_{1}\), \(s_{5}\) and \(s_{3}\). For brevity of notation, \(s_{I}\) denotes the column vector \([s_{i}]_{i\in I}\).

**Definition 1** (s-zonotope [28]).: _A symbolic zonotope \(\mathcal{X}_{|s}\) is an affine symbolic function that can be written in the form \(c+Rs_{I}\) where vector \(c\) and matrix \(R\) do not depend on the symbolic variables in \(s_{I}\). Notation: \(\mathcal{X}_{|s}=\langle c,R,I\rangle_{|s}=c+Rs_{I}\)._

**Definition 2** (e-zonotope [28]).: _The e-zonotope \(\mathcal{X}_{|s}\) related to the s-zonotope \(\mathcal{X}_{|s}=\langle c,R,I\rangle_{|s}=c+Rs_{I}\) is the set-valued interpretation of \(\mathcal{X}_{|s}\) as \(\mathcal{X}_{|s}=\langle c,R,I\rangle_{|s}=\{c+R\sigma|\sigma\in\mathcal{D}(s_ {I})\}\)._

A basic example is : given \(i\in\mathbb{N}\), \(\langle 0,1,i\rangle_{|s}=s_{i}\) (symbolic expression corresponding to the \(i\)-th symbol in \(s\)) and \(\langle 0,1,i\rangle_{|s}=\mathcal{D}(s_{i})=[-1,+1]\) (set-valued interpretation of \(s_{i}\)). More generally, s-zonotopes and their interpretation as e-zonotopes make it possible to explicitly perform operations either at symbolic/syntactic level (\(._{|s}\)) or at semantic level (\(._{|s}\)).

**Remark 1**.: _In this work, all symbols being of type unit interval, \(c\) being a real vector and \(R\) a real matrix, \(\mathcal{X}_{|s}\) is a classical zonotope \(\langle c,R\rangle\) with center \(c\) and generator matrix \(R\) (for extensions to other symbol types, see [28]). Note that \(\langle 0,1,i\rangle_{|s}-\langle 0,1,j\rangle_{|s}=s_{i}-s_{j}=0\) for \(i=j\), whereas \(\langle 0,1,i\rangle_{|t}-\langle 0,1,j\rangle_{|t}=[-2,+2]\) for all \((i,j)\), that is, even for \(i=j\). Operating at the symbolic/syntactic level thus permits more accurate set evaluations by preserving trace of symbolic dependencies. This is a key point to prevent from pessimistic outer approximations induced by the so-called dependency problem [29] affecting natural interval arithmetic and other classical set-based operations only considering the semantic level._

From a computational point of view, an s-zonotope is defined by storing the triplet \((c,R,I)\). Due to their affine structure, a key aspect is to efficiently trace the identifier \(i\in I\) of the symbol that multiplies each column of the matrix \(R\). To that end, Matrices with Labeled Columns (MLCs), constitute a data structure featuring columnwise sparsity: It is defined by the pair \((R,I)\) that allows for efficiently recasting standard operations involving s-zonotopes as set-operations on the identifiers vector (\(I\)) and column-wise operations in the projection matrices (\(R\)). For how to translate operations such as sum or linear image onto a computational platform using MLCs the interested reader can refer to [30].

Due to their relevance in further developments, the following operations involving s-zonotopes are briefly recalled.

**Lemma 1** (common symbols [30]).: _Any two s-zonotopes \(\mathcal{X}_{|s}=\langle c_{x},R,I\rangle_{|s}\) and \(\mathcal{Y}_{|s}=\langle c_{y},G,J\rangle_{|s}\), can be rewritten using a common set of symbols \(s_{K}\) as \(\mathcal{X}_{|s}=\langle c_{x},\tilde{R},K\rangle_{|s}\) and \(\mathcal{Y}_{|s}=\langle c_{y},\tilde{G},K\rangle_{|s}\), with_

\[\begin{split}\tilde{R}&=\begin{bmatrix}R_{1},&R_{2},&0 \end{bmatrix},\quad\tilde{G}=\begin{bmatrix}G_{1},&0,&G_{2}\end{bmatrix},\\ K&=\begin{bmatrix}I\cap J;&I\setminus J;&J\setminus I\end{bmatrix}.\end{split} \tag{1}\]

Matrices \((R_{1},G_{1})\) in Lemma 1 may be empty matrices if \(I\cap J\) is empty (similarly for \(R_{2}\) and \(I\setminus J\) or \(G_{2}\) and \(J\setminus I\)).

**Definition 3** (basic operations [30]).: _Given two s-zonotopes \(\mathcal{X}_{|s}\) and \(\mathcal{Y}_{|s}\) with a common set of symbols \(s_{K}\) as in Lemma 1, then their sum and vertical concatenation are the s-zonotopes_

\[\mathcal{X}_{|s}+\mathcal{Y}_{|s} =\left\langle c_{x}+c_{y},[R_{1}+G_{1},\,R_{2},\,G_{2}],K\right\rangle _{|s}, \tag{2}\] \[\left[\mathcal{X}_{|s};\mathcal{Y}_{|s}\right] =\left\langle\begin{bmatrix}c_{x}\\ c_{y}\end{bmatrix},\begin{bmatrix}R_{1}&R_{2}&0\\ G_{1}&0&G_{2}\end{bmatrix},K\right\rangle_{|s}. \tag{3}\]

**Definition 4** (inclusion [28]).: _The s-zonotope \(\mathcal{Y}_{|s}\) is said to include the s-zonotope \(\mathcal{X}_{|s}\), if the set-valued interpretation of \(\mathcal{Y}_{|s}\) includes the set-valued interpretation of \(\mathcal{X}_{|s}\). In other words, the expression \(\mathcal{X}_{|s}\subset\mathcal{Y}_{|s}\) interprets as \(\mathcal{X}_{|s}\subset\mathcal{Y}_{|s}\)._

Definition 4 serves the way for rewriting rules (at symbolic level) that may be either inclusion preserving or inclusion neutral or none of both at set-based evaluation (semantic) level. A more formal treatment of this topic can be found in the definition 27 (rewriting rules and inclusion) in [28], where a definition of inclusion functions is also given in definition 2.

**Definition 5** (reduction [28]).: _The reduction operator \(\downarrow_{q}\) transforms an s-zonotope \(\mathcal{X}_{|s}=\langle c,R,I\rangle_{|s}\) into a new s-zonotope \(\tilde{\mathcal{X}}_{|s}=\downarrow_{q}\mathcal{X}_{|s}=\langle c,G,J\rangle_{|s}\), such that \(\tilde{\mathcal{X}}_{|s}\) includes \(\mathcal{X}_{|s}\) while depending on at most \(q\) symbols, i.e. \(card(J)\leq q\)._

Reduction is thus an inclusion preserving transform. In Definition 5, \(I\cap J\neq\emptyset\) is not mandatory but often useful to prevent from a further propagation of conservative approximations, while controlling the complexity of \(\tilde{\mathcal{X}}_{|s}\) through the maximum number \(q\) of its symbols/generators. In this context, preserving the more significant symbols/dependencies is often beneficial: as in [30], if \(p>q\) a common practice is to replace the \(p-q+1\) less important symbols by a new independent one while guaranteeing the inclusion \(\mathcal{X}_{|s}\subseteq\tilde{\mathcal{X}}_{|s}\). Besides, note that new symbols introduced to characterize independent behaviors must be uniquely identified. Wherever needed, the generation of a vector of \(n\) new unique symbols identifiers is denoted as \(!(n)\). The generation of a pre-specified number of identifiers can be attained by implementing, for example, the Unique Symbols Provider (USP) service introduced in [30].

## III Problem statement

### _System description_

Consider the interconnection of a discrete-time non-linear dynamic model (4) and a neural network. The physical system is modeled as:

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

where \(x(t)\in\mathbb{R}^{n_{x}}\) and \(u(t)\in\mathbb{R}^{n_{u}}\) respectively refer to the state and the control input at time step \(t\in\mathbb{N}\). For all \(t\geq 0\), vector \(w(t)\) accounts for modeling errors and process disturbances and satisfies \(w(t)\in\mathcal{W}=[-1,\,+1]^{n_{u}}\).

The system (4) is controlled by a state-feedback controller \(g(x(t)):\mathbb{R}^{n_{x}}\mapsto\mathbb{R}^{n_{u}}\) parameterized by an \(l\)-layer feed-forward fully connected neural network. The map \(x\mapsto g(x)\) is described by the following recursive equations

\[\begin{split}& x^{(0)}=x,\\ & x^{(k+1)}=\phi^{(k)}(W^{(k)}x^{(k)}+b^{(k)}),\quad k=0,...,l-1, \\ & g(x)=W^{(l)}x^{(l)}+b^{(l)},\end{split} \tag{5}\]
where \(x^{(k)}\in\mathbb{R}^{n_{k}}\) are the outputs (post-activation) of the \(k\)-th layer. The weight matrix \(W^{(k)}\in\mathbb{R}^{n_{k+1}\times n_{k}}\) and bias \(b^{(k)}\in\mathbb{R}^{n_{k+1}}\) define the affine mapping \(z^{(k)}=W^{(k)}x^{(k)}+b^{(k)}\) for the \((k+1)\)-th layer. Besides, the vector-valued function \(\phi^{(k)}:\mathbb{R}^{n_{k+1}}\rightarrow\mathbb{R}^{n_{k+1}}\) is applied element-wise to the pre-activation vector \(z^{(k)}\), that is, \(\phi^{(k)}(z^{(k)})=[\varphi(z_{1}^{(k)}),\cdots,\varphi(z_{n_{k+1}}^{(k)})]^{T}\), where \(\varphi:\mathbb{R}\rightarrow\mathbb{R}\) is the (scalar) activation function. Common activation choices are: ReLU \(\varphi(z)=max(0,z)\); sigmoid \(\varphi(z)=\frac{1}{1+e^{-z}}\); and tanh \(\varphi(z)=tanh(z)\).

The closed-loop system with dynamics (4) and a previously trained neural-network control policy (5), is governed by

\[x(t+1)=f_{g}(x(t),w(t))=f\big{(}x(t),g(x(t)),w(t)\big{)}. \tag{6}\]

Accordingly, given an initial set \(\mathcal{X}_{0}\subset\mathbb{R}^{n_{x}}\), the forward reachable set of (6) at time step \(t\) is denoted as \(\mathcal{X}(t)\). For \(t\geq 1\), this set is defined as:

\[\mathcal{X}(t)=\big{\{}x(t)\,| \exists(x(0),w(0:t-1))\in\mathcal{X}_{0}\times\mathcal{W}\times...\times\mathcal{W}, \tag{7}\] \[\forall\tau\in[0,t-1],x(\tau+1)=f_{g}(x(\tau),w(\tau))\big{\}}.\]

### _Finite-time reach-avoid (RA) verification problem_

Given a goal set \(\mathcal{G}\subset\mathbb{R}^{n_{x}}\) a sequence of avoid sets \(\mathcal{A}(t)\subset\mathbb{R}^{n_{x}}\) and a finite time horizon \(N\in\mathbb{N}^{+}\), it is desired to test whether

\[\mathcal{X}(N)\subseteq\mathcal{G} \tag{8}\] \[\mathcal{X}(t)\cap\mathcal{A}(t)=\emptyset,\quad\forall t=0,...,N-1\]

holds true for the closed loop system (6). In general, the exact evaluation of (8) for a NNCSS is a computationally intractable problem. Thus, the problem is resorted to iteratively compute a tractable over-approximation of the reachable set \(\mathcal{X}(t)\subseteq\mathcal{\bar{X}}(t)\), to test (8) using \(\mathcal{\bar{X}}(t)\) instead. Because of the over-approximation, the proposed verification setting only provides one-sided guarantees, that is, if \(\mathcal{\bar{X}}(t)\) satisfies (8) then it can be guaranteed that (7) will satisfy the RA property, but no sound conclusion about the safety of (7) can be made if the over-approximation \(\mathcal{\bar{X}}(t)\) violates (8). Therefore, the computation of tight over-approximations is of paramount importance, so that a maximum number of truly satisfied specifications can be computationally proven as such.

## IV Closed-loop verification using s-zonotopes

This section presents the methodology for computing a sound over-approximation of the closed-loop system that preserves system-controller linear dependencies. The computation takes advantage of s-zonotopes described in the previous section. The abstraction of the control loop components using affine symbolic expressions is presented below.

### _Initial set_

It is assumed that the initial set can be described by the set-valued interpretation of an s-zonotope \(\mathcal{X}_{|s}(0)=\langle c_{0},R_{0},I_{0}\rangle_{|s}\), where \(c_{0}\in\mathbb{R}^{n_{x}}\) and \(R_{0}\in\mathbb{R}^{n_{x}\times n_{0}}\) and \(I_{0}=!(n_{0})\) is a set of \(n_{0}\) unique identifiers for the interval valued symbols \(s_{l_{0}}\). In other words, it is assumed that \(\mathcal{X}_{0}=\mathcal{X}_{|t}(0)\). Note that, any arbitrary zonotopic set \(\{c+R\xi\,\|\,\|\xi\|_{\infty}\leq 1\}\) can be abstracted as an s-zonotope by characterizing the independent behaviour of the generators through new interval type symbols.

### _NN controller affine abstraction_

For the sake of simplicity of notations, the temporal notation is dropped here. Given a state bounding s-zonotope \(\mathcal{X}_{|s}=\langle c,R,I\rangle_{|s}\) and a NN controller (5), the idea is to abstract the NN behavior through an affine symbolic expression of the form

\[\mathcal{U}_{|s}=\langle C_{u},[G,\,H],[I;\,J]\rangle_{|s}=c_{u}+Gs_{I}+Hs_{J}, \tag{9}\]

such that, it guarantees the local enclosure of the network outputs, i.e. \(g(\mathcal{X}_{|t})\subseteq\mathcal{U}_{|t}\). Note that, expression (9) captures the linear dependencies of the state symbols (identified by \(I\)), plus the addition of new error symbols (identified by \(J\)) that are introduced to guarantee the soundness of the method.

The computation of vector \(c_{u}\), matrices \(G,H\), and the identifiers vector \(J\) is discussed below. The focus is on generating a dependencies-preserving inclusion for an arbitrary layer of NN (5), since a sound enclosure for the whole network follows by induction due its sequential nature. For simplicity, the layer superscript is removed below and the superscript \({}^{+}\) is used to denote the next layer.

**Affine mapping** Given the \(s\)-zonotope \(\mathcal{X}_{|s}=\langle c,R,I\rangle_{|s}\), the affine mapping \(\mathcal{Z}_{|s}=W\mathcal{X}_{|s}+b\) in the layers of (5) yields a (pre-activation) s-zonotope of the form

\[\mathcal{Z}_{|s}=\langle\check{c},\check{R},I\rangle_{|s}, \tag{10}\] \[\check{c}=Wc+b,\quad\check{R}=WR.\]

**Activation functions** Activation functions \(\varphi(\cdot)\) in (5) are applied element-wise to the pre-activation vector. Hence, the projection of \(\mathcal{Z}_{|s}\) onto the \(i\)-th neuron, yields the s-zonotope

\[\mathcal{Z}_{|s}=\langle\check{c}_{i},\check{R}_{[i,:]},I\rangle_{|s}. \tag{11}\]

Notice that, any point belonging to set-valued interpretation \(\mathcal{Z}_{i|t}\) of (11) is confined within an interval \([l_{i},\,u_{i}]\), where, since \(\mathcal{Z}_{i|t}\) is a one dimensional zonotopic set, it follows that \(\mathcal{Z}_{i|t}=[l_{i},\,u_{i}]\) with the lower and upper bounds

\[l_{i}=\check{c}_{i}-\|\check{R}_{[i,:]}\|_{1},\quad u_{i}=\check{c}_{i}+\| \check{R}_{[i,:]}\|_{1}. \tag{12}\]

Therefore, the soundness of the method can be certified by guaranteeing the inclusion (see Definition 4) of the graph of the activation function in the range \([l_{i},\,u_{i}]\). To that end, the activation function \(\varphi(\cdot)\) is abstracted through an affine symbolic function of the form

\[\mathcal{X}_{i|s}^{+}=\alpha_{i}\mathcal{Z}_{i|s}+\beta_{i}+\gamma_{i}s_{j}, \tag{13}\]

where \(s_{j}\) represents a new independent symbol (identified through \(j=!(1)\)) that must be introduced to guarantee the full coverage of the activation function graph on the considered range, that is, in order to satisfy the condition

\[\begin{bmatrix}\mathcal{Z}_{i|t}\\ \varphi(\mathcal{Z}_{i|t})\end{bmatrix}\subseteq\begin{bmatrix}\mathcal{Z}_{i|t} \\ \alpha_{i}\mathcal{Z}_{i|t}+\beta_{i}+\gamma_{i}\mathcal{D}(s_{j})\end{bmatrix}. \tag{14}\]

The \(i\)-th neuron post-activation s-zonotope \(\mathcal{X}_{i|s}^{+}\) in (13) not only guarantees that its set-valued interpretation encloses the neuron output, but it preserves the linear influence of the symbols \(s_{I}\) in the output set. This later point plays a fundamental role since it allows to retain the interplay between the inputs of the neurons at the same layer. Coherently, the
layer post-activation s-zonotope can be computed by vertically concatenating in a recursive fashion the different \(\mathcal{X}_{i|s}^{+}\) after rewriting them using the same set of symbols

\[\mathcal{X}_{|s}^{+}=[\,...\,[\,[\mathcal{X}_{1|s}^{+};\mathcal{X}_{2|s}^{+}]; \mathcal{X}_{3|s}^{+}]\,...\,;\mathcal{X}_{n_{k}|s}^{+}]. \tag{15}\]

**Proposition 1** (NN s-zonotope).: _Given the s-zonotope \(\mathcal{X}_{|s}^{(0)}=\langle c^{(0)},R^{(0)},I^{(0)}\rangle_{|s}\), and let \(\alpha^{(k)},\beta^{(k)},\gamma^{(k)}\in\mathbb{R}^{n_{k}+1}\) be some parameter vectors that guarantee the inclusion of the \(n_{k+1}\) activation functions in the \(k\)-th layer, then the enclosure of the NN output set \(g(\mathcal{X}_{|_{i}}^{(0)})\subseteq\mathcal{U}_{|_{i}}=\langle c_{u},[G,\,H ],[I;\,J]\rangle_{\iota}\) is guaranteed for the s-zonotope in (9) with parameters_

\[c^{(k+1)}=diag(\alpha^{(k)})(W^{(k)}c^{(k)}+b^{(k)})+\beta^{(k)}, \tag{16a}\] \[\tilde{H}^{(k+1)}=\big{[}diag(\alpha^{(k)})W^{(k)}\tilde{H}^{(k)},\quad diag(\gamma^{(k)})\big{]}\,,\] (16b) \[\tilde{G}^{(k+1)}=diag(\alpha^{(k)})W^{(k)}\tilde{G}^{(k)},\quad k =1,...,l-1,\] (16c) \[c_{u}=W^{(l)}c^{(l)}+b^{(l)},\] (16d) \[H=W^{(l)}\tilde{H}^{(l)},\] (16e) \[K=W^{(l)}\tilde{K}^{(l)},\] (16f) \[J=[!(n_{1});\,...;!(n_{l})], \tag{16g}\]

_where \(\tilde{G}^{(1)}=diag(\beta^{(0)})W^{(0)}R^{(0)}\) and \(\tilde{H}^{(1)}=diag(\gamma^{(0)})\)._

Proof.: Expressions (16a)-(16f) result from the recursive application of Lemma 1 and the vertical concatenation of the post-activation s-zonotopes (13) for the \(n_{k+1}\) neurons of the \(k\)-th layer. Besides, (16g) reflects the symbols identifier update of the noise terms introduced at the neurons of each layer.

Regarding the output inclusion, starting with an initial set \(\mathcal{X}_{|_{i}}^{(0)}\), by induction, given the pre-activation s-zonotope \(\mathcal{X}_{|_{i}}^{(k)}\), the operations at the \(k\)-th layer are: affine mapping; linear abstraction (inclusion preserving for appropriate triplet \((\alpha_{i}^{(k)},\beta_{i}^{(},\gamma_{i}^{(k)}))\)); and vertical concatenation. Thus, the composition of inclusion functions being an inclusion function, the proof follows. 

For each neuron, the triplet of parameters \((\alpha,\beta,\gamma)\) must be appropriately designed to satisfy (14), while minimizing the conservatism induced by using an affine relaxation. In this regard, a relevant heuristic consists in minimizing the magnitude of the error symbol introduced to guarantee the activation function graph enclosure, i.e. to minimize \(|\gamma|\). Due to the independent behaviour of the error symbol, this can be reformulated as minimizing the area of the enclosing parallelogram [31].

**Lemma 2**.: _Given the bounds \([l,\,u]\) in (12) with \(l<u\), the triplet of parameters \((\alpha^{*},\beta^{*},\gamma^{*})\) that minimizes \(|\gamma|\) while guaranteeing the satisfaction of (14) are:_

* _ReLU function_ \(\varphi(x)=max(0,x)\)__ \[\alpha^{*}=\frac{\varphi(u)-\varphi(l)}{u-l},\quad\beta^{*}=\gamma^{*}=\frac{ \varphi(l)-\alpha^{*}\cdot l}{2}.\] (17)
* _S-shaped functions_
* _Sigmoid_ \(\varphi(x)=\frac{1}{1+e^{-x}}\) _with_ \(\varphi^{\prime}(x)=\varphi(x)(1-\varphi(x))\)__
* _tanh_ \(\varphi(x)=tanh(x)\) _with_ \(\varphi^{\prime}(x)=1-\varphi(x)^{2}\)__ \[\alpha^{*}=\min(\varphi^{\prime}(l),\varphi^{\prime}(u)),\] \[\beta^{*}=\frac{\varphi(u)+\varphi(l)-\alpha^{*}\cdot(u+l)}{2},\] (18) \[\gamma^{*}=\frac{\varphi(u)-\varphi(l)-\alpha^{*}\cdot(u-l)}{2}.\]

**Remark 2**.: _The proposed NN abstraction method shares a similar structure with the zonotope abstraction based on affine arithmetic presented in [31] for the (open-loop) NN output bounding. However, here, the explicitly computed affine symbolic expression (9) will further play a key role in closed-loop verification, and an efficient computation of the projection matrices exploiting the generators (matrix) structure of s-zonotopes, is also used._

### _Dynamical system affine abstraction_

Similar to the NN controller dynamics (5), the function (4) that describes the state evolution at time \((t+1)\) can be abstracted by means of an (inclusion preserving) affine mapping. The resulting s-zonotope will depend on the symbols that define the state at time \(t\), plus some extra symbols that account for: I) NN controller non-linearities; II) abstract system non-linearities; III) the uncertainty sources.

For the computation of a state bounding s-zonotope, it is assumed that the function \(f(\cdot)\) in (4) results from the composition of elementary functions and operators for which an affine symbolic expression (s-zonotope) can be computed. Note that this is not much restrictive since (linear) operations such as linear image, sum or vertical concatenation are closed (i.e. they return s-zonotopes) under affine mappings. Besides, any univariate locally continuous differentiable function can be abstracted through an affine mapping.

**Lemma 3**.: _Let \(h:[l,\,u]\rightarrow\mathbb{R}\) be a class \(\mathcal{C}^{1}\) function on a given interval \([\,l,\,u]\subset\mathbb{R}\). Then, the function \(\tilde{h}(x,\epsilon)=\alpha x+\beta+\gamma\epsilon\) satisfies that \(\forall x\in[\,l,\,u],\,\,\exists\epsilon\in[-1,\,+1],h(x)=\tilde{h}(x,\epsilon)\) for the triplet of parameters:_

\[\alpha=\frac{h(u)-h(l)}{u-l},\quad\beta=\frac{h(\underline{x})+h( \bar{x})-\alpha(\underline{x}+\bar{x})}{2},\] \[\gamma=\frac{h(\bar{x})-h(\underline{x})+\alpha(\underline{x}- \bar{x})}{2},\]

_where, defining \(\xi(x)=h(x)-\alpha x\), then_

\[\bar{x}=\operatorname*{arg\,max}_{x\in\{\delta_{1},...,\delta_{n},u\}}\xi(x), \quad\,\underline{x}=\operatorname*{arg\,min}_{x\in\{\delta_{1},...,\delta_{n}, u\}}\xi(x),\]

_with \(\delta_{1},...,\delta_{n}\) the stationary-points of \(\xi(\cdot)\) in \([l,\,u]\)._

Proof.: See Appendix A. 

Lemma 3 provides a method to propagate (inclusion preserving) s-zonotopes through univariate non-linearities. Besides, the approach in Lemma 3 returns an optimal, in the sense of minimizing the magnitude \(|\gamma|\) of the error symbol, set of parameters for convex/concave differentiable functions [28]. On the other hand, the interaction between multiple variables can be handled through the sum operation (2) or by over-approximating the product of two s-zonotopes.
**Lemma 4**.: _Given two 1-D s-zonotopes \(\mathcal{X}_{|s}=\langle c_{x},r^{T},K\rangle_{|s}\) and \(\mathcal{Y}\rangle_{|s}=\langle c_{y},g^{T},K\rangle_{|s}\) with a common set of symbols \(s_{K}\) (with \(n=card(s_{K})\)), then the product \(\mathcal{X}_{|s}\times\mathcal{Y}_{|s}\) is included by the s-zonotope \(\mathcal{L}_{|s}=\langle c_{l},[l^{T},\,m],[K;\,j]\rangle_{|s}\) with_

\[\begin{split}& c_{l}=c_{x}c_{y}+\frac{1}{2}\sum_{i=1}^{n}r_{i}g_{i}, \qquad l=c_{x}g+c_{y}r,\\ & m=\frac{1}{2}\sum_{i=1}^{n}|r_{i}g_{i}|+\sum_{i=1}^{n}\sum_{l> i}^{n}|r_{i}g_{l}+r_{l}g_{i}|,\end{split} \tag{19}\]

_and \(j=!(1)\)._

**Example 1**.: _Consider the system \(x^{+}=\sin(x)-u+0.1w\), with an initial set \(\mathcal{X}_{0}=[0,\,1]\) described by \(\mathcal{X}_{|s,0}=0.5+0.5s_{1}\). This system is controlled by a NN with 1 layer of 2 neurons that, for \(\mathcal{X}_{0}\), is abstracted as \(\mathcal{U}_{|s}=0.1+0.2s_{1}-0.1s_{2}+0s_{3}\). The non-linear function \(h(x)=\sin(x)\) is abstracted, for \(\mathcal{X}_{0}\), as \(\hat{\mathcal{X}}_{|s}=0.45+0.42s_{1}+0.03s_{4}\). Besides, the independent behaviour of the disturbances is captured by \(\mathcal{W}|_{s}=s_{5}\). Accordingly, a dependency preserving over-approximation of the successor state is given by \(\mathcal{X}_{|s}^{+}=0.35+0.22s_{1}+0.1s_{2}+0.03s_{4}+0.1s_{5}\). The a priori knowledge on the number of error symbols introduced at each abstraction (e.g. \(\mathcal{U}_{|s}\) introduces up to two error symbols, one per neuron) allows to directly store the generators matrix of each s-zonotope by taking into account the common set of symbols, thus providing an efficient computation of required operations as shown below:_

\[\begin{split}&\qquad\qquad\qquad c\quad\quad s_{1}\quad\quad s_{2} \quad\quad s_{3}\quad\quad s_{4}\quad\quad s_{5}\\ -&\left[\begin{array}{c|c|ccccc}0.1&0.2&-0.1&0&0&0 \\ 0.45&0.42&0&0&0.03&0\\ \end{array}\right]&\leftrightarrow&\mathcal{U}_{|s}\\ +&\left[\begin{array}{c|c|ccccc}0&0&0&0&0&0.1\\ \end{array}\right]&\leftrightarrow&0.1\cdot\mathcal{W}_{|s}\\ =&\left[\begin{array}{c|ccccc}0.35&0.22&0.1&0&0.03&0.1\\ \end{array}\right]&\leftrightarrow&\mathcal{X}_{|s}^{+}\end{split}\]

```
1:\(\mathcal{X}_{0}\), NN param \((\boldsymbol{W},\boldsymbol{b})\), \(\mathcal{G},\mathcal{A}(i)\), \(N\), \(q\)
2:isRAok, \(t_{err}\), \(\mathcal{X}_{|s}(j)\)\(j=0,...,min(t_{err},N)\).
3:Initialize: Generate \(\mathcal{X}_{|s}(0)\); set \(t_{err}\leftarrow\infty\)
4:for i = 0 to \(N-1\)do
5:if\(\mathcal{X}_{|i}(i)\cap\mathcal{A}(i)\)then
6:\(t_{err}\gets i\)
7:break all
8:else
9:\(\mathcal{U}_{|s}(i)\gets controller(\mathcal{X}_{|s}(i),\boldsymbol{W}, \boldsymbol{b})\)
10:\(\bar{\mathcal{X}}_{|s}(i+1)\gets system(f(\cdot),\mathcal{X}_{|s}(i), \mathcal{U}_{|s}(i),\mathcal{W}_{|s}(i))\)
11:\(\mathcal{X}_{|s}(i+1)\leftarrow_{\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! \!\!
their relevant information (generators, symbol identifiers and order of the monomials) into matrices.

**Definition 6** (s-polynotope [28]).: _A symbolic polynotope \(\mathcal{P}_{|s}\) is a polynomial function that can be written in the form \(c+Rs_{I}^{E}\) where vector \(c\) and matrices \(R\) and \(E\) do not depend on the symbolic variables in \(s_{I}\). Notation: \(\mathcal{P}_{|s}=\langle c,R,I,E\rangle_{|s}=c+Rs_{I}^{E}\)._

Definition 6 uses the exponential matrix notation (as in Definitions 23-25 in [28]), where the usually sparse matrix \(E\) accounts for exponents of the symbols involved in each monomial. As an example, \(s_{I}=[s_{1},s_{2}]^{T}\) and \(E=[1\,0\,3;0\,2\,4]\), yields \(s_{I}^{E}=[s_{1},s_{2}^{2},s_{1}^{3}s_{2}^{4}]^{T}\). Similar to s-zonotopes, a distinction is made between an s-polynotope as defined in Definition 6 and its set-valued interpretation defined as the (possibly non-convex) set \(\mathcal{P}_{|u}=\{c+R\sigma^{E}\,|\,\sigma\in\mathcal{D}(s_{I})\}\).

Symbolic polytopes obviously extend s-zonotopes (obtained from \(E=\mathcal{I}\), i.e. with identity as exponent matrix) and are closed under the extension of basic operations already defined for s-zonotopes like linear image, sum or concatenation. The reader is referred to [28] for further details on how to define and operate on s-polynotopes.

### _NN controller polynomial abstraction_

The abstraction of the I/O map of a NN controller of the form (5) using s-polynotopes is presented below. In particular, given a state bounding s-polynotope \(\mathcal{X}_{|s}=\langle c,R,I,E\rangle_{|s}\) (note that any initial s-zonotope in Section IV-A can be directly transformed into an equivalent s-polynotope) the idea is to compute a polynomial symbolic map of the form

\[\mathcal{U}_{|s}=\langle c_{u},G,Q,E_{u}\rangle=c_{u}+Gs_{Q}^{E_{u}}, \tag{20}\]

such that, the enclosure of the network outputs is guaranteed, i.e. \(g(\mathcal{X}_{|s})\subseteq\mathcal{U}_{|t}\). The vector of identifiers in (20) has the structure \(Q=[I,\,J]\), thus involving the symbols in the state bounding s-polynotope (identified by \(I\)) as well as error symbols (identified by \(J\)). Notice that the exponent matrix \(E_{u}\) may also capture cross terms involving symbols with identifiers in both \(I\) and \(J\).

Similar to section IV-B, the computation of (20) can be obtained from a forward propagation of \(\mathcal{X}_{|s}\) through (in this case) a polynomial relaxation of the activation functions. The pre-activation s-polynotope \(\mathcal{Z}_{|s}=W\mathcal{X}_{|s}+b\) and its projection onto the \(i\)-th neuron are given by

\[\begin{split}\mathcal{Z}_{|s}&=\langle\bar{c},\tilde {R},I,E\rangle_{|s},\qquad\mathcal{Z}_{|s}=\langle\bar{c}_{i},\tilde{R}_{[i, ]},I,E\rangle_{|s},\\ \dot{c}&=Wc+b,\qquad\qquad\qquad\tilde{R}=WR.\end{split} \tag{21}\]

Bounding the set-valued interpretation of \(\mathcal{Z}_{|s}\) within the interval \([l_{i},\,u_{i}]\), then the polynomial structure of s-polynotopes enables to obtain a sound over-approximation of the NN output by locally covering the activation function graph through an \(n\)-order polynomial expression of the form

\[\mathcal{X}_{i|s}^{+}=\sum_{m=1}^{n}\alpha_{i,m}(\mathcal{Z}_{i|s})^{m}+\beta _{i}+\gamma_{i}s_{j}, \tag{22}\]

where \(s_{j}\) represents a new independent symbol (identified through \(j=!(1)\)) introduced to guarantee the enclosure of the activation function graph in the range \([l_{i},\,u_{i}]\). Therefore, since \(\mathcal{X}_{i|s}^{+}\) results in an s-polynotope arising from the polynomial mapping of s-polynotopes, the layer post-activation s-polynotope \(\mathcal{X}_{s}^{+}\) is computed by vertically concatenating the neuron post-activation s-zonotopes.

Polynomial over linear abstractions of the activation function not only allow to reduce the conservatism introduced by the error symbols, but also enable the to compute input-output symbolic relationships that better fit the activation behaviour.

**Example 2**.: _Suppose that the projection onto a ReLU neuron is given by the s-polynotope \(\mathcal{Z}_{|s}=0.5-0.5s_{1}+s_{1}s_{2}\), whose set-valued interpretation is bounded/included in the interval \([l,\,u]=[-1,\,2]\). Then, the ReLU function can be locally abstracted over this range using an \(n=2\)-order polynomial of the form \(\mathcal{X}_{|s}^{+}=\alpha_{2}(\mathcal{Z}_{|s})^{2}+\alpha_{1}\mathcal{Z}_{| s}+\beta+\gamma s_{3}\), with \((\alpha_{2},\,\,\alpha_{1})=(0.25,\,0.5),\,\,\beta=\gamma=0.125\). This, in turn, generates the post-activation s-polynotope_

\[\begin{split}\mathcal{X}_{|s}^{+}&=0.25(0.5-0.5s_{ 1}+s_{1}s_{2})^{2}+0.5(0.5-0.5s_{1}+s_{1}s_{2})\\ &\quad+0.125+0.125s_{3}\\ &=0.4375-0.375s_{1}+0.0625s_{1}^{2}+0.125s_{3}+0.75s_{1}s_{2}\\ &\quad-0.25s_{1}^{2}s_{2}+0.25s_{1}^{2}s_{2}^{2}.\end{split}\]

_Figure 1 depicts the non-convex local enclosure generated by the set-valued interpretation of \([\mathcal{Z}_{|s};\mathcal{X}_{|s}^{+}]\)._

Example 2 evidences the complexity/accuracy trade-off inherent to using \(n\)-order polynomial abstractions: a high \(n\) grants an accurate representation of the activation functions; whereas, on the other hand, it increases the computational complexity due to the increased number of monomials. A reduction strategy is thus used to manage the representation complexity. It can consist in either truncating the maximum degree of the polynomial approximation (20), or limiting the maximum number of monomials involved. To address this latter issue, an approach consists in (independently) assessing the monomials relevance based on the \(2\)-norm of the generators (matrix columns) [32], and use natural interval extension [33] to bound the list of selected monomials through a reduced number of independent symbols.

The ability of the final s-polynotope (20) to generate a sound over-approximation of the network outputs is guaranteed by selecting (for each neuron) a triplet \((\alpha,\beta,\gamma)\), where

Fig. 1: 2-nd order polynomial enclosure of ReLu function (Example 2)
\(\alpha=[\alpha_{n},...,\alpha_{1}]^{T}\) is an \(n\)-dimensional vector, that ensures the (local) coverage of the activation function. In the case of the commonly used ReLU activation functions, as shown in Example 2 their convex nature allows describe them more accurately than with the sole affine dependencies by using \(2\)-nd order polynomial expressions. The reduction in the magnitude related to the error symbol is especially significant in those situations where an affine approximation yields a rough description of ReLU function, i.e. for \(|l|\approx u\) (if \(l<0<u\)).

**Proposition 2**.: _Given the interval \([l,\,u]\) with \(l<0<u\) and a ReLU activation function \(\varphi(x)=max(x,0)\). The set of parameters_

\(\alpha_{2}=\frac{1}{2u},\ \alpha_{1}=1-\alpha_{2}u,\ \beta=\gamma=\frac{\alpha_{2}u^{2}}{8}\) _if_ \(|l|\leq u\leq 2|l|\)__

\(\alpha_{2}=\frac{\gamma}{2},\ \alpha_{1}=-\alpha_{2}l,\ \beta=\gamma=\frac{ \alpha_{2}l^{2}}{8}\) _if_ \(u<|l|\leq 2u\)__

_guarantees that \(\eta(x,\epsilon)=\alpha_{2}x^{2}+\alpha_{1}x+\beta+\gamma\epsilon\) satisfies that \(\forall x\in[l,\,u],\exists\epsilon[-1,\,1],\ \varphi(x)=\eta(x,\epsilon)\) with \(|\gamma|\leq\frac{3}{8}|\gamma_{aff}^{*}|\), where \(\gamma_{aff}^{*}\) is the error symbol introduced by the affine abstraction in Lemma 2._

Proof.: See Appendix B. 

### _Finite-time RA using s-polynotopes_

The main aspects in addressing the RA verification problem using s-polynotopes are discussed below. In general, the same steps presented in Algorithm 1 can be used while adapting the NN controller and dynamical system abstraction to an s-polynotope formulation. In this case, the computation of the s-polynotope \(\mathcal{U}_{|s}\) in Step 7 of Algorithm 1 has already been presented in section V-B. Regarding the abstraction of the nonlinear function \(f(\cdot)\) in Step 8, since s-polynotopes constitute an extension of s-zonotopes, \(f(\cdot)\) can always be abstracted using (at least) the affine dependency preserving method in Lemma 3. Note that, s-polynotopes also enable the description of (multivariate) polynomial equations without the need to over-approximate them, at least in all intermediate symbolic compositions and up to some tunable computation load.

It must be taken into account that several operations in a s-polynotope formulation of Algorithm 1 such as bounding the projection of an s-polynotope onto a neuron, intersection/inclusion of an s-polynotope with an avoid/reach set or the reduction operator (in Step 9), in turn require the computation of interval bounds from (multivariate) interval polynomial expressions. If computationally affordable, the range bounds computed by a (simple and fast) interval extension may be refined either by iteratively bisecting the variables domain, or resorting to numerically reliable optimization-based methods.

## VI Input partitioning strategy

In general, the conservatism induced by abstraction-based verification tools strongly depends on the size of the initial set. It is thus extremely useful to assess the regions of the initial/input space for which meaningful (i.e. not too coarse) over-approximations of the closed-loop system evolution can be obtained. To that end, this section presents an algorithm to split the initial set of a NNCSs verification problem in a smartly guided way. More precisely, the proposed splitting strategy relies on and benefits from the dependency modeling and tracing used in Section IV. In particular, the algorithm assesses the sensitiveness of the initial/input directions on the satisfaction of a safety property by an s-zonotopic over-approximation through the analysis of the relative influence of the initial symbols. The principle of the algorithm is introduced in Section VI-A. Then, some relevant notions are detailed in section VI-B, whereas the algorithm pseudo-code for a RA problem implementation is reported in Section VI-C. Finally, some further discussion on different settings is presented.

### _Splitting principle_

The main idea of the proposed algorithm is to keep a linear increase in the number of subsets by only splitting at each iteration the sole initial/input set symbol that has greater influence on the satisfaction of the safety property \(\mathcal{S}\) to be verified. To that end, notice that given any initial s-zonotope \(\mathcal{X}_{|s}(0)=\langle c_{0},R_{0},I\rangle_{|s}\) as defined in Section IV-A, then the successive computation of forward reachable sets returns over-approximating s-zonotopes structured as

\[\mathcal{X}_{|s}(t)=c_{f}+R_{f}s_{I}+G_{f}s_{J}, \tag{23}\]

where the matrix \(R_{f}\) (resp. \(G_{f}\)) reflects the impact of the initial (resp. error) symbols identified by \(I\) (resp. \(J\)) on the computed over-approximation at time \(t\). Typically, testing \(\mathcal{S}(\mathcal{X}_{|s}(t))\) boils down to a metric/size evaluation on the set-valued interpretation of \(\mathcal{X}_{|s}(t)\) (e.g. to check threshold preassigning). Hence, due to the linearity of (23), the influence of each input symbol \(s_{i}\)\((i\in I)\) can be assessed using a metric that gauges the generator (column of \(R_{f}\)) size that is related to \(s_{i}\), whereas the accuracy of an s-zonotope approximation to evaluate \(\mathcal{S}\) can be determined by measuring the zonotope \(\langle 0,G_{f}\rangle\) spanned by the error symbols.

Therefore, at each iteration of the algorithm, an input s-zonotope \(\mathcal{X}_{|s}(0)\), such that the corresponding output/final s-zonotope does not satisfy \(\mathcal{S}\), is split into two new input s-zonotopes that are later evaluated on the satisfaction of \(\mathcal{S}\). The algorithm may run until the satisfaction of the safety property, or until the accuracy of the method (gauged through \(\langle 0,G_{f}\rangle\)) is below a certain threshold.

### _Relevant notions_

Some relevant notions for the s-zonotope based partitioning algorithm are discussed below.

#### Vi-B1 Accuracy assessment

considering the evaluation of a safety property for a s-zonotope of the form (23), the accuracy of the over approximation can be assessed by gauging the zonotope \(\langle 0,G_{f}\rangle\) spanned by a (set-valued) interpretation of the error symbols. In particular, further implementations make use of the zonotope \(F\)-radius [34], that is,the Frobenius norm of the generators matrix \(\|G_{f}\|_{F}\) to reasoning upon the quality of the affine approximation.

#### Vi-B2 Input symbols relative influence

the sensitivity of an input symbol \(s_{i}\)\((i\in I)\) is computed based on the F-radius ratios of the I/O zonotopes spanned by \(\iota s_{i}\). That is, through the ratio \(\|R_{0}^{[i]}\|_{2}/\|R_{0}^{[i]}\|_{2}\) where \(R_{0}^{[i]}\) (and \(R_{f}^{[i]}\)) denote the columns
of \(R_{0}\) (and \(R_{f}\)) that multiply the symbol \(s_{i}\). This relation is used to quantify how a variation on \(s_{i}\) at the input s-zonotope \(\mathcal{X}_{|s}(0)\) affects the output s-zonotope.

#### Vi-B3 Symbol bisection

bisecting a unit interval symbol \(s_{i}\ (i\in I)\) is done by rewriting it as \(s_{i}\to 0.5+0.5s_{j}\) and \(s_{1}\to-0.5+0.5s_{k}\), where \(j=!(1)\) and \(k=!(1)\), thus generating two new s-zonotopes.

#### Vi-B4 Polyhedral RA sets

checking the empty intersection and/or the inclusion of a state bounding set with/within a polyhedron in half-space representation of the type \(\{h_{i}^{T}x\leq r_{i},\ i=1,...,m\}\) can be done by evaluating the infimum/supremum of the projections of the bounding set onto the directions \(h_{i}\in\mathbb{R}^{n_{x}}\)[35]. For a state bounding s-zonotope of the form (23), the supremum of the dot product with \(h\) is computed as

\[\sup_{x\in\mathcal{X}_{|s}(t)}h^{T}x=h^{T}c_{f}+\|h^{T}R_{f}\|_{1}+\|h^{T}G_{f} \|_{1}. \tag{24}\]

### _Algorithm implementation for finite-time RA_

Algorithm 2 reflects the pseudo-code of the proposed partitioning strategy to check the satisfaction of a RA problem over a time horizon \(N\). Algorithm 2 uses square brackets to label the different s-zonotopes that arise after splitting. As an example, \(\mathcal{X}_{|s}[0]\) reads as the initial s-zonotope (i.e. \(\mathcal{X}_{|s}[0]=\mathcal{X}_{|s}^{0}\)), which is split onto a second \(\mathcal{X}_{|s}[1]\) and a third \(\mathcal{X}_{|s}[2]\) s-zonotopes (with \(\mathcal{X}_{|s}[1]\cup\mathcal{X}_{|s}[2]=\mathcal{X}_{|s}[0]\)). Besides, \(L\) denotes a set of integer labels/indices (for the above example \(L=\{0,1,2\}\)), and \(\mathcal{X}_{|s}[L]\) is a shorthand for the set of s-zonotopes \(\{\mathcal{X}_{|s}[l]\ |\ l\in L\}\).

At each iteration, the routine reach runs a slightly modified version of Algorithm 1, that, in this case, returns the last time instant (and the corresponding s-zonotope) for which the RA problem is not satisfied. These times-to-last-error are managed by vector \(T\). The algorithm iteratively selects the label of the initial s-zonotope that yields the largest time-to-last-error (Step 14). The use of this backward management of the information that prioritizes to split until the RA constraints are satisfied at time \(t\), then at time \(t-1\), etc., will be further discussed in the next paragraph. Once the \(l\)-th (with \(l\in L\)) s-zonotope has been selected, sym-select returns the initial symbol identifier \(i\in I\) that has greater relative influence over the violated property. The symbol \(s_{i}\) of the \(l\)-th set is split by the routine sym-split that returns two new initial subsets (Step 6). The times-to-last-error for the new s-zonotopes are computed and the set \(L\) and vector \(T\) updated (Steps 7-11). In particular, Algorithm 2 runs either until the RA problem is satisfied for the whole set \(L\), or until a maximum number \(n_{max}\) of splits is reached.

```
0: same as Algorithm 1, \(n_{max}\)
0: isRAok, set of s-zonotopes \(\mathcal{X}_{|s}[L]\)
1:Initialize:\(l=n=0\); \(L=\{l\}\); \(\mathcal{X}_{|s}[0]=\mathcal{X}_{0}\)
2:\((t,\mathcal{X}_{|s}(t)[0])\leftarrow\texttt{reach}(\mathcal{X}_{|s}[0],N)\)
3:\(T\gets append(t)\)
4:while\((\max(T)>0)\ \vee\ (n/2==n_{max})\)do
5:\(i\leftarrow\texttt{sym-select}(\mathcal{X}_{|s}[l],\mathcal{X}_{|s}(T(l))[l])\)
6:\((\mathcal{X}_{|s}[n+1],\mathcal{X}_{|s}[n+2])\leftarrow\texttt{sym-split}( \mathcal{X}_{|s}[l],i)\)
7:for\(j=1\) to \(2\)do
8:\((t,\mathcal{X}_{|s}(t)[n+j])\leftarrow\texttt{reach}(\mathcal{X}_{|s}[n+j], \max(T))\)
9:\(L\gets L\cup\{n+j\}\)
10:\(T\gets append(t)\)
11:endfor
12:\(L\gets L\setminus\{l\}\)
13:\(T(l)\gets 0\)
14:\(l\leftarrow\arg\max(T(L))\)
15:\(n\gets n+2\)
16:endwhile
17: isRAok = (\(\max(T)==0\))
```

**Algorithm 2** Input partitioning for finite-time RA

Algorithm 2 manages the information in a backward fashion, that is, it selects a s-zonotope with the higher time-to-last-error. This usually gives better results than working in a forward fashion (that is, selecting the s-zonotope with lower time-to-first-error) since it avoids to get stuck by exhaustively splitting up to the satisfaction of a constraint at time \(t\), which, then, may have a small impact on the constraint satisfaction at \(t+1\). On this subject, the algorithm can be straightforwardly adapted to handle the forward case by directly using Algorithm 1 (instead of reach), using \(T(l)\gets N+1\) (instead of \(T(l)\gets 0\)) in Step 13 and selecting the minimum (instead of the maximum) of vector \(T\). Besides, the reduction operation used in Algorithm 1 must not truncate the initial symbols even if their relevance decreases with time, so that the input-output mapping of the symbols identified by \(I\) is preserved.

### _Other possible settings and applications_

Other choices for the proposed input partitioning strategy are as follows:

* The strategy in Algorithm 2 can be adapted to handle open-loop verification problems (e.g. elements like the NN in isolation). In this case, the reach routine will only compute the output s-zonotope for the isolated element for a number of forward steps \(N=1\).
* The maximum number of splits stopping criterion in Algorithm 2 can be modified/complemented with a tolerance on the accuracy assessment (see Section VI-B). In other words, if the accuracy tolerance is fulfilled and a property is still violated, then the algorithm should stop to prevent from further splitting and the safety property is considered as unsatisfied up to the accuracy tolerance.
* Another interesting application is to modify the s-zonotope split decision rule (Step 14 of Algorithm 2) to focus the split in those regions for which the accuracy of using an affine abstraction is low (i.e. high \(\|G_{f}\|_{F}\)). This tends to return a set of initial s-zonotopes such that each locally provides an accurate (affine) abstraction of the system behavior.

## VII Simulations

### _Benchmark description_

The numerical simulations consist in the discrete-time version of some of the verification problems proposed in the ARCH-COMP 2021 [4]. Five dynamical systems are assessed,
namely: single pendulum (**S**), TORA (**T**), unicycle car (**C**), adaptive cruise control (**ACC**) and double pendulum (**D**). The above systems have been discretized using the forward Euler method with sampling period \(\Delta T\), and they are controlled by a NN controller with control period \(\Delta h\). The NN controllers are the ones provided in [4] to control the continuous-time version of the models. To address this issue, the dynamical models have been analyzed under sampling times \(\Delta T(\leq\Delta h)\) chosen sufficiently small for the discretization to have negligible impact in the model responses. Under this context, the same safety constraints and initial conditions than the ones proposed in [4] have been re-used to setup the reported simulations. Note that, as discussed in Remark 3, the use of symbolic approaches supports the variation of \(\Delta T\) (for some \(\Delta h\)) without inducing conservatism due to a loss of dependencies between repeated control inputs. A detailed description of the systems dynamics can be found in [4], whereas the main parameters that define each safety verification problem are shown in Table I.

### _Benchmark results using \(s\)-zonotopes_

All the results reported below were obtained on a standard laptop with Intel Core i7-8550U@1.8GHz\(\times\)4 processor and 16GB RAM running Windows 10. Table II shows the set of initial states, the safety constraints (with their time horizon), as well as the time required by an s-zonotope implementation to verify each problem. The reduction order is \(q=200\) in all the experiments. Some particularities are discussed below:

* _Single pendulum_ (**S**): in a discrete-time setting, the constraint \(x_{1}\in[0,\,1]\) is guaranteed to be satisfied for problem **S1** (with \(\Delta T=0.05\)s) for the time interval \(t\in[0.55,\,1]\) (that is, for samples \(\{11,...,20\}\)), whereas in **S2** (with \(\Delta T=0.001\)s) the constraint satisfaction is guaranteed for the time interval \(t\in[0.516,\,1]\).
* _TORA_ (**T**): in **T1**, the closed-loop system is not stable for the discrete-time model obtained for \(\Delta T=1\)s. In this case, an unambiguous constraint violation is achieved at \(t=3\) in \(0.036\)s. On the other hand, the closed-loop model obtained in **T2** and **T3** is stable, and the s-zonotope method verifies the satisfaction of the safety constraint in both problems without resorting to split the input set.
* _Unicycle car_ (**C**): the model under study considers the addition of an unknown-but-bounded disturbance \(w\in 10^{-4}[-1,\,+1]\) affecting the fourth state. The safety properties are verified for both **C1** and **C2**. In particular, Figure 2 shows the envelope computed for **C2** in the time interval \(t\in[0,\,10]\) and how the outer-approximation lies within the goal set at \(t=10\).
* _Adaptative cruise control_ (**ACC**): both problems **ACC1** and **ACC2** are verified for the given time horizon.
* _Double pendulum_ (**D**): the set of constraints in problems **D1-2** are violated by the closed-loop system. An unambiguous constraint violation is achieved for **D1** at \(t=0.25\) and for **D2** at \(t=0.278\). On the other hand, the problem **D3** cannot be verified from a simple affine abstraction: the accumulated error indeed increases in the reachability analysis of **D3**, not allowing to guarantee the constraints satisfaction or their unambiguous violation, and thus motivating further extensions.

The results presented above show how, despite their low computational complexity, s-zonotopes yield a high performance in NNCSs verification, being able to verify almost all the benchmark problems without splitting the input set. It is also remarkable the scalability of this approach. As an example, for problem **T3** with \(\Delta T=0.001\)s, \(\Delta h=1\)s and time horizon \(t\in[0,20]\), the proposed tool only requires of \(1.515\)s to compute and assess \(N=20\mathrm{s}/\Delta T=20.000\) forward iterations.

### _Use of s-polynotopes_

The capability of s-polynotopes to capture the non-convex map of NNs is illustrated below. To that end, the set of randomly generated neural networks used in [17] are analyzed. All the NNs consists of 2 inputs, 2 outputs and they differ on the number of hidden layers and neurons per layer. The

Fig. 2: Problem **C2**: framed zonotopes represent the computed bounds at each \(\Delta h=.2\)s; blurred lines represent the bounds update at \(\Delta T=.001\)s.
first four NNs present \(l=\{1,2,3,4\}\) hidden layers, each having \(n_{k}=100\) ReLU neurons per layer. The examined NN input set is \(\mathcal{X}_{0}=[0.9,\,1.1]\times[0.9,\,1.1]\). Figure 3 shows the set-valued interpretation of the output bounding s-polynotopes obtained by abstracting the activity functions of active neurons with second order polynomials i.e. with \(n=2\) in (22). The computation times are \(\{0.178,0.240,2.021,3.329\}s\) for the NNs with \(l=1\) to \(l=4\) hidden layers, respectively. Similarly, another set of NNs with \(l=\{7,8,9,10\}\) hidden layers and \(n_{k}=10\) ReLU neurons per layer is evaluated for the same input set. Figure 4 represents the interpretation of the resulting s-polynotopes that are computed in \(\{0.2389,0.108,0.786,0.155\}\)s, respectively. Those examples (Figure 3 and Figure 4) taken from [17] illustrate the ability of s-polynotopes composition to accurately generate inclusion preserving polynomial I/O mappings of NNs. As a byproduct, an efficient implicit description of possibly non convex output sets is obtained.

### _Partitioning strategy_

Firstly, in order to show the performance of the partitioning algorithm, it will be applied to the open-loop robotic arm example used in [12, 13]. Particularly, the non-linear dynamics of a 2 DOF robot arm are modeled by a \((2,5,2)\) NN with \(tanh\) activations. The considered set of joint angles are extended to \((\theta_{1},\theta_{2})\in[\frac{\pi}{3},\frac{4\pi}{3}]^{2}\). An implementation of Algorithm 2 adapted to analyze the NN in isolation is executed in order to iteratively minimize the \(F\)-radius1 of the zonotope spanned by the error symbols (\(\|G_{f}\|_{F}\)) for a fixed number of \(n_{max}=400\) splits. The computation time of the algorithm is \(0.097\)s. Figure 4(b) shows the resulting pattern of \(401\) input subsets, whereas Figure 4(a) represents the corresponding s-zonotope interpretations obtained in the output space altogether with an exhaustive evaluation of the NN (blue dots). This latter figure shows how Algorithm 2 achieves an accurate description of the non-convex output set by focusing the splitting effort in those regions of the input space for which an affine abstraction granted by s-zonotopes is not accurate enough.

Footnote 1: The \(F\)-radius of a zonotope is the Frobenius norm of its generator matrix (see Definition 3 in [34]).

Furthermore, considering the initial set \((\theta_{1},\theta_{2})\in[\frac{\pi}{3},\frac{2\pi}{3}]^{2}\), Algorithm 2 is set to split up to the satisfaction of the safety constraint \(y_{1}\leq d\) (where \(y_{1}\) denotes the first output). Table III reflects the number of splits and the time required by Algorithm 2 to satisfy the above safety constraint for different values of \(d\). Besides, Table III also shows, for a fixed number of splits, the number of existing possible combinations of set selections and symbols bisections, as well as how many among them are able to satisfy the property. As an example, for \(d=1.2\), Algorithm 2 requires \(8\) splits. For the same

Fig. 4: NNs with 10 neurons per layer and \(l=\{7,8,9,10\}\) hidden layers. Set-valued interpretation of the over-approximating s-polynotope (red set); exhaustive evaluation of the NNs (blue dots).

Fig. 3: NNs with 100 neurons per layer and \(l=\{1,2,3,4\}\) hidden layers. Set-valued interpretation of the over-approximating s-polynotope (red set); exhaustive evaluation of the NNs (blue dots).
problem, there are no possible combinations of less than 7 splits for which the property can be proven; there exist 18 out of \(5.491\cdot 10^{4}\) possibilities that satisfy it with 7 splits (\(0.032\%\)); and 336 out of \(3.66\cdot 10^{5}\) possibilities that satisfy it with 8 splits (\(0.0918\%\)).

Regarding the closed-loop examination of Algorithm 2, this is applied to assess the satisfaction of problem **D3**, which cannot be satisfied by a simple s-zonotope abstraction. To that end, Algorithm 2 is set to split up to the satisfaction of the safety constraints in Table II. The algorithm requires a total of 19 splits (i.e. 20 subsets) computed in \(5.12\)s. Figure 6 shows the time evolution of the interval enclosure of the resulting 20 reachable sets (light blue background), altogether with 50 random simulations of the closed-loop system (blue dots).

## VIII Conclusions

A compositional approach focused on inclusion preserving long term symbolic dependency modeling is introduced in this work for the analysis of NNCSs, where such long term is to be understood both in time iterations (regarding the controlled system dynamics) and in layer iterations (regarding the sole NNs). This results in a generic method that has been developped in several ways. Firstly, the matrix structure of s-zonotopes enables to compute (fast and simple) affine symbolic mappings to abstract the I/O mapping of the control loop components. Two further extensions are also proposed: the use of s-polynotopes to compute inclusion preserving polynomial mappings capable of accurately describing the non-convex map of NNs, and an input partitioning algorithm that benefits from the ability granted by s-zonotopes to preserve linear dependencies between the loop elements. Simulations show the comparative efficiency of the proposals and support the prevalence of dependency preserving methods for closed-loop analysis over the use of accurate, but dependency breaking, output bounding verification tools. Future works should address the integration with the analysis of continuous-time dynamical systems, as well as the study of optimized affine/polynomial abstractions for achieving better performance in verifying specific safety properties.

## Acknowledgments

The authors would like to thank Prof. Vicenc Puig (UPC, Barcelona, Spain) who initially prompted this collaboration and his support during the research stay of the first author at the University of Bordeaux. This work was in part supported by the Margarita Salas grant from the Spanish Ministry of Universities funded by the European Union NextGenerationEU.

Fig. 5: Robot arm example: input partitioning for \((\theta_{1},\theta_{2})\in[\frac{\pi}{2},\frac{4\pi}{3}]^{2}\); yellow-green colors characterize the corresponding I/O set pairs; exhaustive evaluation of the NN (blue dots).

Fig. 6: Input partitioning problem **D3**: interval enclosure of the resulting 20 reachable sets (light blue); random simulations (blue dots).
## Appendix A Proof of the Lemma 3

Define \(\alpha=\frac{h(u)-h(l)}{u-l}\) and form \(\xi(x)=h(x)-\alpha x\) which is continuous in \([l,\,u]\), differentiable in \((l,\,u)\) and satisfies \(\xi(l)=\xi(u)=0\). Then, the maximum \(\xi(\bar{x})\) (resp. minimum \(\xi(\underline{x})\)) of \(\xi(x)\) on \([l,\,u]\) must be for \(\bar{x}\) (resp. \(\underline{x}\)) in the boundary points \(\{u,l\}\) or in its identity points denoted as \(\{\delta_{1},...,\delta_{n}\}\), that is, the solutions of \(\xi^{\prime}(x)=0\to h^{\prime}(x)=\alpha\).

Given \(\xi(\bar{x})\) and \(\xi(\underline{x})\), then for all \(x\in[l,\,u]\)\(\xi(\underline{x})\leq\xi(\underline{x})\leq\xi(\underline{x})\implies\underline{y}(x)\leq h (x)\leq\bar{y}(x)\), with \(\underline{y}(x)=\alpha(x-\underline{x})+h(\underline{x})\) and \(\bar{y}(x)=\alpha(x-\bar{x})+h(\bar{x})\). Thus, it follows that \(\forall x\in[l,\,u],\exists\epsilon\in[-1,+1]\) such that

\[h(x)=\frac{\bar{y}(x)+\underline{y}(x)}{2}+\frac{\bar{y}(x)-\underline{y}(x) }{2}\epsilon\]

or equivalently \(h(x)=\tilde{h}(x,\epsilon)=\alpha x+\beta+\gamma\epsilon\) with

\[\beta =\frac{\bar{y}(x)+\underline{y}(x)}{2}-\alpha x=\frac{h( \underline{x})+h(\bar{x})-\alpha(\underline{x}+\bar{x})}{2},\] \[\gamma =\frac{\bar{y}(x)-\underline{y}(x)}{2}=\frac{h(\bar{x})-h( \underline{x})+\alpha(\underline{x}-\bar{x})}{2}.\]

## Appendix B Proof of the Proposition 2

_Inclusion preservation_: the sign criterion \(\gamma\geq 0\) is chosen below. Given \(\bar{y}(x)=\alpha_{2}x^{2}+\alpha_{1}x+\beta+\gamma\) and \(\underline{y}(x)=\alpha_{2}x^{2}+\alpha_{1}x+\beta-\gamma\), then parameters \((\alpha_{2},\alpha_{1},\beta,\gamma)\) ensure local coverage of \(\varphi(x)\) if \(\underline{y}(x)\leq\varphi(x)\leq\bar{y}(x),\forall x\in[l,\,u]\).

Consider firstly the scenario \(|l|\leq u\leq 2|l|\). In this case, \(\alpha_{2}=\frac{1}{2u}>0\) and thus \(\underline{y}(x),\bar{y}(x)\) are strictly convex.

On the one hand, \(\beta=\gamma\) and \(\alpha_{1}=-\alpha_{2}u\) impose that, \(y(0)=0=\varphi(0)\) and \(\underline{y}(u)=u=\varphi(u)\), whereas \(\underline{y}(l)=\frac{1}{2}(\frac{l}{u}+l)\leq 0=\varphi(l)\) for \(u\geq|l|\) and \(l<0\). Therefore, from \(\underline{y}(l)\leq\varphi(l)\), \(\underline{y}(0)=\varphi(0)\), \(\underline{y}(u)=\varphi(u)\) and the convexity of \(\underline{y}(x)\) wrt \(x\), it follows that \(\underline{y}(x)\leq\varphi(x),\forall x\in[l,\,u]\).

On the other hand, from \(\alpha_{1}=1-\alpha_{2}u\) and \(\beta=\gamma=\frac{\alpha_{2}u^{2}}{8}\), then \(\bar{y}(x)\) is tangent to the positive region of \(\varphi(x)\) in \(\hat{x}=\frac{u}{2}\) (that is, \(\bar{y}(\hat{x})=\varphi(\hat{x})=\hat{x}\) and \(\bar{y}^{\prime}(\hat{x})=\varphi^{\prime}(\hat{x})=1\)), and thus, since \(\bar{y}(x)\) is convex, it follows that \(\bar{y}(x)\geq\varphi(x),\forall x\geq 0\). Additionally, for \(\alpha_{2}=\frac{1}{2u}\) the (global) minimum of \(\bar{y}(x)\) is \(\bar{y}(x^{*})=0\) for \(x^{*}=\frac{-u}{2}\) (that is, \(\bar{y}^{\prime}(x^{*})=0\)), and thus \(\bar{y}(x)\geq\varphi(x),\forall x\leq 0\).

A similar reasoning can be used to show the inclusion preservation for the scenario \(u<|l|\leq 2u\).

_Conservatism reduction:_ For the scenario \(|l|\leq u\leq 2|l|\), the parameter \(\gamma\) has the value \(\gamma=\frac{\alpha_{2}u^{2}}{8}=\frac{u}{16}\). On the other hand, for a ReLU function \(\varphi(x)=max(0,x)\) the triplet for an affine abstraction in Lemma 2 yields \(\gamma^{*}_{aff}=\frac{u|l|}{2(u+|l|)}\). Therefore, for \(u\leq 2|l|\) the following inequality is obtained

\[\gamma^{*}_{aff}=\frac{u|l|}{2(u+|l|)}\geq\frac{u|l|}{2(2|l|+|l|)}=\frac{u}{6} >\frac{u}{16}=\gamma\]

and thus \(\gamma\leq\frac{3}{8}\gamma^{*}_{aff}\sim|\gamma|\leq\frac{3}{8}|\gamma^{*}_{aff}|\) (since \(\gamma,\gamma^{*}_{aff}>0\)). A similar reasoning can be used to prove the case \(u<|l|\leq 2u\).
