# Safety verification of Neural-Network-based controllers:
a set invariance approach

Louis Jouret, Adnane Saoud and Sorin Olaru This work is supported by the ANR
PIA funding: ANR-20-IDEES-0002.Louis Jouret is with Swiss Federal School of
Technology in Lausanne-EPFL, Switzerland<EMAIL_ADDRESS>Saoud is
with the College of Computing, University Mohammed VI Polytechnic, Benguerir,
Morocco<EMAIL_ADDRESS>Olaru is with CentraleSupelec, University
Paris-Saclay, Gif-sur-Yvette, France<EMAIL_ADDRESS>Object Identifier (DOI): 10.1109/LCSYS.2023.3342088

# Safety verification of Neural-Network-based controllers: a set invariance
approach

Louis Jouret, Adnane Saoud and Sorin Olaru This work is supported by the ANR
PIA funding: ANR-20-IDEES-0002.Louis Jouret is with Swiss Federal School of
Technology in Lausanne-EPFL, Switzerland<EMAIL_ADDRESS>Saoud is
with the College of Computing, University Mohammed VI Polytechnic, Benguerir,
Morocco<EMAIL_ADDRESS>Olaru is with CentraleSupelec, University
Paris-Saclay, Gif-sur-Yvette, France<EMAIL_ADDRESS>Object Identifier (DOI): 10.1109/LCSYS.2023.3342088

###### Abstract

This paper presents a novel approach to ensure the safety of continuous-time
linear dynamical systems controlled by a neural network (NN) based state-
feedback. Our method capitalizes on the use of continuous piece-wise affine
(PWA) activation functions (e.g. ReLU) which render the NN a PWA continuous
function. By computing the affine regions of the latter and applying Nagumo’s
theorem, a subset of boundary points can effectively verify the invariance of
a potentially non-convex set. Consequently, an algorithm that partitions the
state space in affine regions is proposed. The scalability of our approach is
thoroughly analyzed, and extensive tests are conducted to validate its
effectiveness.

©2023 IEEE. Personal use of this material is permitted. Permission from IEEE
must be obtained for all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes,
creating new collective works, for resale or redistribution to servers or
lists, or reuse of any copyrighted component of this work in other works.

##  I Introduction

Machine learning, particularly NNs, has had a transformative impact on various
scientific fields, including control systems. Two main approaches have
emerged: the first involves approximating a complex control law such as Model
Predictive Control (MPC) using a NN [1], making the controller more memory
efficient and enabling faster computations. The second approach entails
synthesizing a NN controller through reinforcement learning, which has gained
popularity due to its capability to learn intricate control strategies from
data generated by the system [2]. However, when applying a NN to safety-
critical systems [3], there is often skepticism and valid concerns regarding
their black box nature and the inherent difficulty in interpreting their
behavior. In control systems, safety refers to the property of a system
remaining in a set of safe states for all future time instances. NNs,
characterized by their numerous neurons and nonlinear activation functions,
present computational challenges in explicitly representing the input-output
relationship. This complexity hampers the ability to interpret the actions
taken by a NN-controller and verify its safe operation. Extensive research has
been conducted to address this issue using different approaches. One approach
involves estimating the reachable set of the NN-controlled system [4, 5, 6].
However, existing reachability-based approaches are limited to discrete-time
systems and finite time safety, while in the present work we are dealing with
continuous-time systems and infinite-time safety properties. Another approach
focuses on finding a barrier certificate [7, 8] for the closed-loop system.
The main challenge lies in computing the barrier certificate. Recent
advancements involve training a separate NN that acts as barrier certificates
to the closed-loop system [7, 8]. This approach leverages the universal
approximation capabilities of the NN, enabling them to estimate a barrier
certificate if one exists. Nevertheless, a notable limitation of this method
is its lack of completeness. In other words, if the process fails to identify
a barrier certificate, it remains uncertain whether the failure stems from an
inability to discover the correct certificate or the absence of a valid
certificate altogether. This paper presents a novel approach to ensure safety
of linear dynamical systems controlled by a NN. The method revolves around
Nagumo’s condition, which states that a set is invariant if the vector field
at every point on the boundary points back inside the set [9]. By utilizing
continuous PWA activation functions within the NN, like $ReLU$, the output of
the latter can be explicitly expressed as a continuous composite PWA function.
Nagumo’s condition is then used for the invariance of linear systems
controlled by PWA controllers within a polytopic set. Consequently, the
verification of a small subset of boundary points is sufficient to prove the
set’s invariance. Although the number of regions within which the PWA
controller is affine increases non-polynomially with the number of neurons
[10], the advantage of the present approach relies on the fact that the
calculations are limited to the regions connected to the set’s boundaries. To
this end, an algorithm that leverages the automatic differentiation [11]
offered by modern deep-learning libraries like PyTorch [12] is described and
analysed.

Section II introduces the class of systems and NN considered in the paper. In
Section III, the condition to guarantee the safety of the considered system is
presented. Section IV proposes an algorithmic procedure to apply the presented
method. Finally, Section V focuses on the scalability of the proposed
approach.

Notation: $\mathbb{R}$ and $\mathbb{R}^{+}_{0}$ are the set of reals and of
non-negative reals, respectively. Given a set
$\mathcal{O}\subset\mathbb{R}^{n}$, its border set is denoted
$\partial\mathcal{O}$. For a matrix $M$, we denote $M_{i,j}$ the element in
the $i^{th}$ row and $j^{th}$ column and $M_{i,*}$ the $i^{th}$ row vector.
The matrix $diag(v)$ designates a diagonal matrix having the scalars
$v_{i}\in\mathbb{R}$ on its diagonal. For $v,w\in\mathbb{R}^{n}$, we write
$v\preccurlyeq w$ the element-wise inequality of two vectors such that
$v_{i}\leq w_{i},\forall i$. To represent the multiplication of multiple
matrices, denoted as $M_{1},M_{2},\ldots,M_{n}$, we introduce the notation
${\overset{\curvearrowleft}{\prod}^{n}_{k=1}}M_{k}$. In this notation, the
left arrow signifies that the matrices are multiplied from right to left,
starting with $M_{n}$ and ending with $M_{1}$, resulting in the final matrix
product. Considering a NN with $L$ hidden layers and $N$ neurons per hidden
layer, we denote $n_{in}$, $n_{out}$ and $n_{in}\times N^{(L)}\times n_{out}$
the network’s input dimension, output dimension and architecture respectively.
The number of neurons on layer $l$ will be written as $n^{(l)}$. We write
$\binom{n}{k}$ the binomial coefficient $\frac{n!}{k!(n-k)!}$.

##  II Preliminaries and Problem Formulation

###  II-A Linear Dynamical Systems

In this paper, we consider a linear system $\Sigma$ described by

$\displaystyle\dot{x}=Ax+Bu$ (1)

where $x\in\mathbb{R}^{n}$ is the state and $u\in\mathbb{R}^{m}$ is the
control input.

###  II-B Feed-forward Neural Networks

A feed-forward neural network (NN) consists of interconnected layers, with
each layer containing individual neurons [13]. The connection weights between
layer $l-1$ and layer $l$ are denoted as $W^{(l)}$, and layer $l$’s biases are
represented as $b^{(l)}$. The output $z^{(l)}$ of layer $l$, of a NN with an
input $x$, can be expressed as follows:

$z^{(l)}=\sigma^{(l)}(W^{(l)}z^{(l-1)}+b^{(l)}),\quad
l\in\\{1,2,...,L\\},\quad z^{(0)}=x,$

with $z^{(l-1)}\in\mathbb{R}^{n^{(l-1)}}$,
$W^{(l)}\in\mathbb{R}^{n^{(l)}\times n^{(l-1)}}$ and
$b^{(l)}\in\mathbb{R}^{n^{(l)}}$ and the componentwise activation function
$\sigma^{(l)}:\mathbb{R}^{n^{(l)}}\rightarrow\mathbb{R}^{n^{(l)}}$ introduces
a nonlinear behavior to the neuron’s output. The overall output of the NN is
$\mathcal{N}(x;\theta)=z^{(L)}$, where $\theta$ are the parameters of the NN,
i.e. the weights and biases.

###  II-C Convex Polytopes

This section recalls key concepts related to polytopes [14], highlighting
their main geometric properties for this study.

######  Definition II.1.

Let $C\in\mathbb{R}^{m\times n}$ and $d\in\mathbb{R}^{m}$. A closed convex
polyhedral set is defined in its $\mathcal{H}$-representation as follows:
$\mathcal{R}=\\{x\in\mathbb{R}^{n}\mid Cx\leq d\\}$.

When the inequality is strict (i.e., $Cx<d$), the polyhedron is referred to as
open. A polytope is a bounded polyhedron.

######  Definition II.2.

Let $\mathcal{V}=\\{v_{1},\ldots,v_{n}\\}$ be a finite set of points in
$\mathbb{R}^{d}$. The convex hull of $\mathcal{V}$ defines a polytope as
follows:

$\displaystyle\mathcal{R}$
$\displaystyle=\left\\{\sum_{k=1}^{n}\lambda_{i}v_{i}\mid\lambda_{i}\geq
0,\sum_{k=1}^{n}\lambda_{i}=1,v_{i}\in\mathcal{V}\right\\}$

This is also known as the $\mathcal{V}$-representation of a polytope and
$\mathcal{V}(\mathcal{R})$ is the set of vertices of the polytope
$\mathcal{R}$.

######  Definition II.3.

Let $\mathcal{R}\subset\mathbb{R}^{n}$ be a convex polytope defined by
$\mathcal{R}=\\{x\in\mathbb{R}^{n}\mid Cx\leq d\\}$ where
$C\in\mathbb{R}^{m\times n}$ and $d\in\mathbb{R}^{m}$. Consider the hyperplane
$\mathcal{H}_{i}=\\{x\in\mathbb{R}^{n}\mid
C_{i,*}x=d_{i}\\},i\in\\{1,\ldots,m\\}$. The face $\mathcal{F}_{i}$ of
$\mathcal{R}$ is defined as
$\mathcal{F}_{i}(\mathcal{R})=\mathcal{H}_{i}\cap\mathcal{R}$ and
$\mathcal{F}(\mathcal{R})$ will denote the set of all the faces of
$\mathcal{R}$.

###  II-D Neural Network controlled systems

The linear system $\Sigma$ in (1) is considered to be controlled by a NN,
represented by the state-feedback function $u(x)=\mathcal{N}(x;\theta)$ [3].
The NN receives the system’s states and produces the corresponding control
input for the system. This class of systems is referred to as NN-controlled
systems and has been extensively explored recently [4, 5, 8]. In this case,
the closed-loop system is given by

$\displaystyle\dot{x}(t)=f(x(t))=Ax(t)+B\mathcal{N}(x(t);\theta)$ (2)

###  II-E Problem formulation

In this paper, we consider the following problem:

######  Problem II.1.

Consider the closed-loop NN-controlled system in (2). Let $\mathcal{S}$ be a
polytopic set of admissible states in $\mathbb{R}^{n}$, and let
$\mathcal{O}_{1},\mathcal{O}_{2},\ldots,\mathcal{O}_{m}\subset\mathcal{S}$ be
$m$ open polytopic unsafe sets. Consider the safe set
$\mathcal{X}:=\mathcal{S}\setminus\bigcup_{i=1}^{m}\mathcal{O}_{i}$. The
objective is to verify that for any trajectory
$x:\mathbb{R}_{0}^{+}\rightarrow\mathbb{R}^{n}$ satisfying
$x(0)\in\mathcal{X}$, we have that $x(t)\in\mathcal{X}$, for all
$t\in\mathbb{R}_{0}^{+}$.

Intuitively, the objective is to prove the positive invariance of the
potentially non-convex set $\mathcal{X}$ for the NN-controlled system in (2).

##  III Main Results

###  III-A Piece-wise affine Neural Networks

Consider the activation function $\sigma:\mathbb{R}\rightarrow\mathbb{R}$ to
be PWA:

$\displaystyle\sigma(x)$ $\displaystyle=\left\\{\begin{matrix}c_{1}\cdot
x+d_{1}&\text{if }x\leq m_{1}\\\ c_{2}\cdot x+d_{2}&\text{if }m_{1}<x\leq
m_{2}\\\ \vdots\\\ c_{n}\cdot x+d_{n}&\text{if }m_{n-1}<x\\\
\end{matrix}\right.$ (3)

where $c_{i},d_{i},m_{i}\in\mathbb{R}$, $i\in\\{1,\ldots,n\\}$. Commonly used
functions that satisfy this condition are ReLU, leaky-ReLU, PReLU or the
Binary Step function. Note that the approach presented is specifically
tailored for NNs with PWA activation functions. The exploration to accommodate
nonlinear activation functions will be a subject of future research. The
activation function $\sigma^{(l)}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ of
a layer $l$ with $n$ neurons can then be written as:

$\displaystyle\sigma^{(l)}(x)=[\sigma(x_{1}),\ldots,\sigma(x_{n})]^{T}=C_{\Phi(x)}^{(l)}x+d_{\Phi(x)}^{(l)},\forall
x\in\mathbb{R}^{n}$

where $x=(x_{1},x_{2},\ldots,x_{n})^{T}\in\mathbb{R}^{n}$,
$C_{\Phi(x)}^{(l)}\in\mathbb{R}^{n\times n}$ is a square diagonal matrix and
$d_{\Phi(x)}^{(l)}\in\mathbb{R}^{n}$.

Note that the elements of $C_{\Phi(x)}^{(l)}$ and $d_{\Phi(x)}^{(l)}$ are the
first order coefficients of $\sigma$ for $x_{1},x_{2},...,x_{n}$ respectively,
i.e. $\forall i\in\\{1,\ldots,n\\}$,

$\displaystyle\exists j\text{ s.t. }m_{j-1}<x_{i}\leq m_{j}\text{
and}\begin{cases}C_{{\Phi(x)}^{i,i}}^{(l)}=c_{j}\\\
d_{{\Phi(x)}^{i}}^{(l)}=d_{j}\end{cases}$

The subscript $\Phi(x)$ will throughout this paper represent the index of the
affine region of the NN when given an input $x$. Assume that
$x\in[\underline{m},\overline{m}]$, i.e,
$\underline{m}_{i}<x_{i}\leq\overline{m}_{i}$, $i=1,\ldots,n$, where
$\underline{m},\overline{m}\in\mathbb{R}^{n}$. We refer to
$\underline{m},\overline{m}$ as delimiters of the vector $x$. These vectors
essentially define the boundary points for each element in $x$ and when an
individual element crosses these boundaries, $C_{\Phi(x)}^{(l)}$ and
$d_{\Phi(x)}^{(l)}$ are altered accordingly.
Given an input $x$, the output of layer $l$ can be written as:

$\displaystyle z^{(l)}$ $\displaystyle=\sigma^{(l)}(W^{(l)}z^{(l-1)}+b^{(l)})$
$\displaystyle=(C_{\Phi(x)}^{(l)}W^{(l)})z^{(l-1)}+(C_{\Phi(x)}^{(l)}b^{(l)}+d_{\Phi(x)}^{(l)})$
(4)

Thus, we can write the output of layer $l$ w.r.t. the input $x$:

$\displaystyle z^{(l)}$
$\displaystyle=\mathcal{E}^{(l)}_{\Phi(x)}x+\mathcal{G}^{(l)}_{\Phi(x)}\text{
, where}$ (5)
$\displaystyle\mathcal{E}^{(l)}_{\Phi(x)}={\overset{\curvearrowleft}{\prod}^{l}_{k=1}}(C_{\Phi(x)}^{(k)}W^{(k)})$
$\displaystyle\mathcal{G}^{(l)}_{\Phi(x)}=\sum_{k=1}^{l}[{\overset{\curvearrowleft}{\prod}^{l}_{j=k+1}}C_{\Phi(x)}^{(j)}W^{(j)}](C_{\Phi(x)}^{(k)}b^{(k)}+d_{\Phi(x)}^{(k)})$

We call $\mathcal{E}^{(l)}_{\Phi(x)}$ and $\mathcal{G}^{(l)}_{\Phi(x)}$ the
active parameters of $x$.

######  Definition III.1.

Let $\mathcal{N}$ be a feed-forward NN, utilizing solely PWA activation
functions, defined by (3). Let $\mathcal{X}$ be the set of all possible inputs
to $\mathcal{N}$. A linear region $\mathcal{R}$ of the layer $(l)$ of
$\mathcal{N}(x;\theta)$ is defined as the set:

$\displaystyle\mathcal{R}=\big{\\{}x\in\mathcal{X}\mid\underline{$m$}\preccurlyeq
W^{(l)}z^{(l-1)}+b^{(l)}\preccurlyeq\overline{m},$ $\displaystyle
z^{(l-1)}=\mathcal{E}^{(l-1)}_{\Phi(x)}x+\mathcal{G}^{(l-1)}_{\Phi(x)}\big{\\}}$

where $\underline{$m$}$, $\overline{m}$ are the delimiters of
$W^{(l)}z^{(l-1)}+b^{(l)}$ in $\sigma^{(l)}$.

It is proven that each linear region is a convex polytope [15]. Intuitively, a
linear region is the convex set of states that share the same active
parameters. Moreover, if the PWA activation functions are all continuous, the
output will be continuous as well. The closed-loop system in (2) becomes:

$\displaystyle\dot{x}=Ax+B\mathcal{N}(x;\theta)=(A+B\mathcal{E}^{(L)}_{\Phi(x)})x+B\mathcal{G}^{(L)}_{\Phi(x)}$
(6)

Note that $\Phi(x)=j$ if $x\in\mathcal{R}_{j}$ where $R_{j}$ is a polytope
corresponding to a linear region and is defined as $\mathcal{R}_{j}=\\{z\mid
C_{j}z\leq d_{j}\\}$. Hence the dynamic of the closed-loop system is
continuous PWA and will change when the state crosses into another linear
region.

###  III-B Invariance condition

For the PWA dynamical system to remain inside a safe set according to Nagumo
theorem also [9, Theorem 4.7], the vector field has to point tangentially or
inward at every point along the set’s boundary. Considering our problem
statement, we can take advantage of the linear constraints of the polytopic
linear regions and the linearity of the closed-loop system within these
regions to simplify this condition. Indeed, we show that it is sufficient to
examine solely the vertices of the set we want to verify and the points on its
border where a transition between linear regions occurs.

######  Theorem III.1.

Consider a polytopic set $\mathcal{S}=\\{x\in\mathbb{R}^{n}\mid
C^{\mathcal{S}}x\leq d^{\mathcal{S}}\\}$ and let
$\mathcal{F}(\mathcal{S})=\\{\mathcal{F}_{1},\ldots,\mathcal{F}_{N}\\}$ the
faces of $\mathcal{S}$. Consider the linear switched system $\Sigma$ in (6)
defined by $\dot{x}=f(x)=A_{\Phi(x)}x+b_{\Phi(x)}$ where $\Phi(x)=j$ if
$x\in\mathcal{R}_{j}=\\{z\mid C_{j}z\leq d_{j}\\}$, $1\leq j\leq M$.
Furthermore, assume $\mathcal{S}\subseteq\bigcup_{j=1}^{M}\mathcal{R}_{j}$.
Consider the set

$\displaystyle\mathcal{D}_{ij}=\mathcal{F}_{i}\cap\mathcal{R}_{j},i\in\\{1,\ldots,N\\},j\in\\{1,\ldots,M\\}$
(7)

Then, for any trajectory $x:\mathbb{R}_{0}^{+}\rightarrow\mathbb{R}^{n}$
satisfying $x(0)\in\mathcal{S}$, we have $x(t)\in\mathcal{S},\forall
t\in\mathbb{R}_{0}^{+}$ for the system $\Sigma$ if and only if the following
condition is satisfied:

$\displaystyle C^{\mathcal{S}}_{i,*}\cdot(A_{j}v+b_{j})\leq 0,$
$\displaystyle\forall i\in\\{1,\ldots,N\\},\forall
j\in\\{1,\ldots,M\\},\forall v\in\mathcal{V}(\mathcal{D}_{ij})$ (8)

###### Proof.

Let us rewrite $\mathcal{S}=\\{x\in\mathbb{R}^{n}\mid
g^{\mathcal{S}}_{i}(x)\leq 0,i\in\\{1,\ldots,N\\}\\}$ where
$g^{\mathcal{S}}_{i}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ defined as
$g^{\mathcal{S}}_{i}(x)=C^{\mathcal{S}}_{i,*}x-d_{i}$. Consider a point
$x\in\partial\mathcal{S}$. Then there exists $i\in\\{1,\ldots,N\\}$, such that
$x\in\mathcal{F}_{i}(\mathcal{S})$. Using the fact that
$\mathcal{S}\subseteq\bigcup_{j=1}^{M}\mathcal{R}_{j}$ there exists
$j\in\\{1,\ldots,M\\}$ such that $x\in\mathcal{D}_{ij}$. Hence, one gets

$\displaystyle x=\sum_{k=1}^{l}\lambda_{k}v_{k}\text{ where
}\lambda_{k}\in[0,1],\sum_{k=0}^{l}\lambda_{k}=1,v_{k}\in\mathcal{V}(\mathcal{D}_{ij})$

where $l$ is the number of vertices of the polytope $\mathcal{D}_{ij}$. Thus,

$\displaystyle\nabla g^{\mathcal{S}}(x)^{T}f(x)$
$\displaystyle=(C^{\mathcal{S}}_{i,*})^{T}\cdot(A_{j}x+b_{j})=(C^{\mathcal{S}}_{i,*})^{T}\cdot(A_{j}\sum_{k=0}^{l}\lambda_{k}v_{k}+b_{j})$
$\displaystyle=\sum_{k=0}^{l}\underbrace{\lambda_{k}[(C^{\mathcal{S}}_{i,*})^{T}\cdot(A_{j}v_{k}+b_{j})}_{\leq
0}]\leq 0$

implying that $f(x)\in T_{S}(x)$111$T_{\mathcal{S}}(x)$ denotes the tangent
cone to the set $\mathcal{S}$ at the point $x$ (see [9, Definition 4.6]).. It
follows from [9, Theorem 4.7] that $\mathcal{S}$ is invariant for the
considered system. Now assume that $\mathcal{S}$ is an invariant for the
system $\Sigma$. We have from [9, Theorem 4.7] that $f(x)\in T_{S}(x)$ for all
$x\in\partial S$. Now consider $v_{ij}\in\mathcal{D}_{ij}$,
$i\in\\{1,\ldots,N\\}$, $j\in\\{1,\ldots,M\\}$. Since $v_{ij}\in\partial S$,
it follows that $f(v_{ij})\in T_{S}(v_{ij})$, which in turn implies (8). ∎

Figure 1: Every linear region has a different color. A finite set of vertices
must be verified to assess the invariance of
$\mathcal{S}\setminus\mathcal{O}$.

For an open convex polytopic unsafe set $\mathcal{O}=\\{x\in\mathbb{R}^{n}\mid
C^{\mathcal{O}}x<d^{\mathcal{O}}\\}$, the inequality of equation (8) changes,
such that for any trajectory $x:\mathbb{R}_{0}^{+}\rightarrow\mathbb{R}^{n}$
satisfying $x(0)\notin\mathcal{O}$, we have $x(t)\notin\mathcal{O},\forall
t\in\mathbb{R}_{0}^{+}$ for the system $\Sigma$ if and only if the following
condition is satisfied:

$\displaystyle C^{\mathcal{O}}_{i,*}\cdot(A_{j}v+b_{j})\geq 0,$
$\displaystyle\forall i\in\\{1,\ldots,N\\},\forall
j\in\\{1,\ldots,M\\},\forall v\in\mathcal{V}(\mathcal{D}_{ij})$ (9)

We call (8) and (9) the invariance condition. An intuitive 2-dimensional
example is depicted in Figure 1.

##  IV Algorithmic Implementation

###  IV-A Segmentation into Linear Regions

The division of the state space into linear regions comes from the breakpoints
in the PWA activation function of each neuron. The hyperplane separating two
regions can be written w.r.t the output of the previous layer $z^{(l-1)}$:

$\\{z^{(l-1)}\in\mathbb{R}^{n^{(l-1)}}\mid
W^{(l)}_{n,*}z^{(l-1)}+b^{(l)}_{n}=m_{k}\\}$ (10)

where $W^{(l)}_{n,*}$ are the weights linked to neuron $n$ in layer $l$,
$b^{(l)}_{n}$ its bias and $m_{k}$ the $k^{th}$ breakpoint of its activation
function. In a multi-layer NN, every layer of the network will cut every
region $\mathcal{R}_{j}$ coming from the previous layer with (10). As the
linear regions should be expressed w.r.t. the state space, we can use (5) to
rewrite (10) as:

$\displaystyle\mathcal{H}^{(l)}_{n,k}\\{\mathcal{R}_{j}\\}=$
$\displaystyle\\{x\in\mathbb{R}^{n_{in}}\mid
W^{(l)}_{n,*}[\mathcal{E}_{j}^{(l-1)}x+\mathcal{G}_{j}^{(l-1)}]+b_{n}^{(l)}=m_{k}\\}$
(11)

This hyperplane represents the border between two linear regions of layer
$(l)$ within region $\mathcal{R}_{j}$. Note that the hyperplane depends on the
active parameters of the linear region it cuts. As a result, the hyperplane is
only valid inside its associated linear region. Equation (11) also shows that
the computation of a region of a given layer is dependent on a region of the
previous layer. Consequently, we compute all the linear regions of a layer
before moving to the next layer.

The active parameters can be computed using an iterative formula of (3).
However, modern deep-learning libraries (e.g. PyTorch [12]) build a
computation-graph using automatic differentiation [11] making the computation
of the gradient of any element in the network w.r.t. another possible and
efficient. As the output of a layer inside a linear region has by definition a
constant derivative, the Jacobian matrix only needs to be evaluated for a
single point inside the region. We take the Chebyshev center $x_{c}$, i.e. the
center of the largest Euclidean ball that is entirely contained within the
polytope, because it is by definition inside the region and its computation
scales well for higher dimensions. We can express the active parameters in a
more simple way than (5):

$\displaystyle\mathcal{E}^{(l)}_{\Phi(x_{c})}$
$\displaystyle=\left.\frac{\partial z^{(l)}}{\partial
x}\right|_{x_{c}};\mathcal{G}^{(l)}_{\Phi(x_{c})}$
$\displaystyle=\left.z^{(l)}\right|_{x_{c}}-\left.\frac{\partial
z^{(l)}}{\partial x}\right|_{x_{c}}\cdot x_{c}$

Figure 2 depicts the segmentation algorithm. It takes as input the NN
$\mathcal{N}$, a closed polytopic set $\mathcal{S}\subset\mathbb{R}^{n_{in}}$
and
$\mathcal{O}=\\{\mathcal{O}_{1},\ldots,\mathcal{O}_{m}\\}\subset\mathcal{S}$ a
list of open polytopic sets and returns the list of linear regions of the
state-space. Note that the segmentation algorithm will provide a unique
decomposition of the set $\mathcal{S}\setminus\mathcal{O}$ under a given NN.
Figure 3 shows how the segmentation can be seen as a tree-graph. If we
consider $\mathcal{R}_{4}$ of no interest, we can ignore the computation of
all its subregions and hence save considerable time. Thus, in our approach, a
substantial portion of the computation time can be saved since we do not need
to process linear regions that do not intersect the border of the set to be
analysed. An example would be the purple region in Figure 1.

Figure 2: Segmentation algorithm: $\mathcal{L}_{\mathcal{R}}^{(l)}$ is the
list of all the regions of layer $l$. The algorithm cuts every region of
$\mathcal{L}_{\mathcal{R}}^{(l)}$ with their associated hyperplanes of layer
$l+1$. Figure 3: Segmentation of the green polytope in $\mathbb{R}^{2}$ by a
NN having for every layer $l$, $\sigma^{(l)}(x)=max(0,x)$. When a region
$\mathcal{R}_{j}$ of layer $(l)$ is being cut by the forthcoming layer
$(l+1)$, we first compute the active parameters $\mathcal{E}_{j}^{(l)}$ and
$\mathcal{G}_{j}^{(l)}$ of the region. These are then used to compute all the
hyperplanes $\mathcal{H}^{(l+1)}\\{\mathcal{R}_{j}\\}$ in layer $(l+1)$.
Finally, these hyperplanes cut (only) the region $\mathcal{R}_{j}$ in
subregions.

###  IV-B Safety Verification

Once we have divided the state-space into linear regions, we use Algorithm 1
to verify the safety of the closed-loop system. It assesses whether the vector
field guides the system back inside the safe set for each vertex of the
invariance condition.

Algorithm 1 Verifying the safety of a NN-controller

0: neural network $\mathcal{N}$, closed polytopic set
$\mathcal{S}\subset\mathbb{R}^{n_{in}}$,
$\mathcal{O}=\\{\mathcal{O}_{1},\ldots,\mathcal{O}_{n}\\}\subset\mathcal{S}$ a
list of open polytopic sets and the closed-loop dynamical system $\Sigma$ in
(2).

0: a boolean safety Initialisation : $\texttt{safety}\leftarrow\texttt{TRUE}$

1:
$\mathcal{L}_{\mathcal{R}}\leftarrow\texttt{segmentation\\_algorithm}(\mathcal{N},\mathcal{S},\mathcal{O})$

2: for each polytope $\mathcal{R}$ in $\mathcal{L}_{\mathcal{R}}$ do

3: for each face $\mathcal{F}_{i}$ in $\mathcal{F}(\mathcal{S})$ do

4: $\mathcal{P}=\mathcal{R}\cap\mathcal{F}_{i}$

5: if $\mathcal{P}\neq\emptyset$ then

6: for each vertex $v$ in $\mathcal{V}(\mathcal{P})$ do

7: if $C^{\mathcal{S}}_{i,*}f(v)>0$ then

8: $\texttt{safety}\leftarrow\texttt{FALSE}$

9: for each polytope $\mathcal{O}_{k}$ in $\mathcal{O}$ do

10: for each face $\mathcal{F}_{i}$ in $\mathcal{F}(\partial\mathcal{O}_{k})$
do

11: $\mathcal{P}=\mathcal{R}\cap\mathcal{F}_{i}$

12: if $\mathcal{P}\neq\emptyset$ then

13: for each vertex $v$ in $\mathcal{V}(\mathcal{P})$ do

14: if $C^{\mathcal{O}}_{i,*}f(v)<0$ then

15: $\texttt{safety}\leftarrow\texttt{FALSE}$

16: return safety

##  V Numerical Examples

###  V-A Mobile robot

We consider a 2-dimensional environment in which an agent navigates. The agent
behaves as an integrator and is represented by: $\dot{x}_{2\times
1}=u_{2\times 1}$. The safety of the controlled system will be verified for
$\mathcal{S}_{1}=\\{(x,y)\mid-5\leq x\leq 5,-5\leq y\leq 5\\}$ and two
polytopic unsafe sets $\mathcal{O}_{1}$ and $\mathcal{O}_{2}$. Figure 4
depicts the environment and the linear regions identified by the segmentation
algorithm illustrated in Figure 2. The NN is trained using Reinforcement
Learning, more specifically DDPG [2], to reach a target position while
avoiding obstacles. The architecture of the NN controller is $2\times
128^{(2)}\times 2$. The activation function utilized in the hidden layers is
the $\textit{leaky-relu}(x)=\max(0.01x,x)$ and on the last layer we use a
linearized version of $tanh(x)$.

(a) 2D environment with the agent in blue and the unsafe polytopic sets
$\mathcal{O}_{1}$ and $\mathcal{O}_{2}$ in red

(b) Segmentation of the state space into 1002 regions by a $2\times
128^{(2)}\times 2$ NN controller

Figure 4: Visualization of the Mobile Robot

###  V-B Spring-mass-damper

The system comprises $n$ interconnected wagons, with each wagon having two
variables: position and velocity. Figure 5 depicts the system. The system’s
representation is as follows:

$\displaystyle\dot{x}=\begin{bmatrix}[0]_{n\times n}&I_{n}\\\
-M^{-1}K&M^{-1}C\\\ \end{bmatrix}x+\begin{bmatrix}[0]_{n\times n}\\\ M^{-1}\\\
\end{bmatrix}[F_{u_{1}},\ldots,F_{u_{n}}]^{T}$

where

$\displaystyle K=\left[\begin{smallmatrix}(k_{1}+k_{2})&-k_{2}&0&0&\ldots&0\\\
-k_{2}&(k_{2}+k_{3})&-k_{3}&0&\ldots&0\\\
0&-k_{3}&(k_{3}+k_{4})&-k_{4}&\ldots&0\\\
\vdots&\ddots&\ddots&\ddots&\ddots&\vdots\\\
0&\ldots&0&-k_{(n-1)}&(k_{(n-1)}+k_{n})&-k_{n}\\\
0&\ldots&\ldots&0&-k_{n}&k_{n}\\\ \end{smallmatrix}\right],$

$M=diag[m_{1},m_{2},\ldots,m_{n}]$,
$x=[z_{1},\ldots,z_{n},\dot{z}_{1},\ldots,\dot{z}_{n}]$, where $z_{i}$ is the
position of wagon $i$, $i=1,\ldots,n$ and $C$ has the same structure as $K$.
The system is well-suited for the investigation of the scalability of our
approach w.r.t. the dimensionality of the systems. Algorithm 1 is applied to
the set $\mathcal{S}_{2}=\\{x\in\mathbb{R}^{2n}\mid 0\leq x_{i}\leq
1,|x_{n+i}|\leq|x_{i}|,i=1,2,\ldots,n\\}$. We train the NN by approximating an
MPC controller. During training, the initial state is chosen from
$\\{x\in\mathbb{R}^{2n}\mid|x_{i}|<1,i=1,\ldots,n\\}$ and we impose $-1\leq
u_{i}\leq 1,i=1,\ldots,n$.

Figure 5: Spring-mass-damper with 3 wagons, i.e. 6 states

###  V-C Scalability

The investigation of the time complexity relies on estimating the count of
linear regions in a multi-layer, multi-breakpoint NN, which unfortunately
remains an open research question [10, 16]. Nevertheless, we aim to provide a
clear idea of how computational time scales.

Considering the segmentation algorithm, the worst-case time complexity for
partitioning a region by a layer is directly related to the maximum number of
subregions that a set of $k$ hyperplanes can partition a $d$-dimensional
space. We can define $K(n^{(l)};d)$ as the time complexity associated with
segmenting a region $\mathcal{R}$ with dimensionality $d$ by layer $l$:

$\displaystyle
K(n^{(l)};d)=O\left(1+n^{(l)}+\binom{n^{(l)}}{1}+\ldots+\binom{n^{(l)}}{d}\right)$

Let $\\#\mathcal{R}(l)$ represent the count of regions in layer $l$ that
intersect the border of the verification set. Consequently, the time
complexity of the segmentation algorithm becomes:

$\displaystyle\centering
O\left(K(n^{(1)};d)+K(n^{(2)};d)\\#\mathcal{R}(1)+\ldots+K(n^{(L)};d)\\#\mathcal{R}(L-1)\right)\@add@centering$

The time complexity of Algorithm 1 is contingent upon the count of linear
regions identified and their number of vertices. The time complexity of the
vertex enumeration problem amounts to $O(f^{2}dv)$, where $f$ denotes the
number of faces of the region and $v$ signifies the number of vertices.
Importantly, the number of vertices exhibits exponential growth with the
dimensionality $d$. The time complexity of verifying all the vertices becomes:
$O\left(\\#\mathcal{R}(L)\cdot 2^{d}\right)$.

Due to the absence of a formal estimation concerning the count of linear
regions, we conduct empirical tests. We investigate the effects of varying the
depth and width of the NN and the dimension of the system. The data presented
in Table V-C, V-C, and V-C includes essential metrics: $\\#N$, $\\#\theta$,
$\\#\mathcal{R}$, and $t_{v}$ representing the number of hidden neurons, the
number of network parameters, the number of linear regions and the
verification time of Algorithm 1, respectively. Table V-C and V-C provide
results from experiments conducted on the Mobile Robot, while Table V-C from
the Spring-mass-damper system. [htb]  Architecture $\\#N$ $\\#\theta$
$\\#\mathcal{R}$ $t_{v}$ [s] $2\times 16^{(2)}\times 2$ 32 354 125 28 $2\times
32^{(2)}\times 2$ 64 1218 248 92 $2\times 64^{(2)}\times 2$ 128 4482 477 307
$2\times 128^{(2)}\times 2$ 256 17154 973 1263 $2\times 256^{(2)}\times 2$ 512
67074 1835 4541 Scalability w.r.t. the NN’s width [htb]  Architecture $\\#N$
$\\#\theta$ $\\#\mathcal{R}$ $t_{v}$ [s] $2\times 32^{(1)}\times 2$ 32 162 122
31 $2\times 32^{(2)}\times 2$ 64 1218 248 92 $2\times 32^{(4)}\times 2$ 128
3330 412 263 $2\times 32^{(6)}\times 2$ 192 5442 597 566 $2\times
32^{(8)}\times 2$ 256 7554 862 1185 Scalability w.r.t. the NN’s depth [htb]
Architecture $\\#N$ $\\#\theta$ $\\#\mathcal{R}$ $t_{v}$ [s] $2\times
8^{(2)}\times 1$ 16 105 37 6 $4\times 8^{(2)}\times 2$ 16 130 253 62 $6\times
8^{(2)}\times 3$ 16 155 1937 384 $8\times 8^{(2)}\times 4$ 16 180 8120 1629
Scalability w.r.t. the system’s dimension When comparing, it is important to
clarify that our approach is currently tailored to linear systems, while
existing barrier-based approaches [7, 8] can extend to nonlinear systems.
However, their numerical examples are limited to systems of dimension 3. Our
approach exhibits favorable scalability when it comes to network size and we
demonstrate that although our method is limited to linear systems, we manage
to treat systems with dimensionality up to 8 in a reasonable time on an
ordinary machine222All experiments were conducted on a 2016 MacBook Pro with
16 GB RAM, a 3.3 GHz Intel Core i7 dual-core processor. All the processing was
done on the CPU. The codebase (Python) is available at
https://github.com/LouisJouret/Neural-Control-Invariance-Checker.. Our
approach provides finite-time safety assessment, a capability lacking in
existing barrier-based methods [7, 8]. Furthermore, if the computational time
of these methods becomes excessive, they may require a trial-and-error
procedure consisting of successive attempts with different NN architectures
until potentially a successful NN-barrier-function is found. Conversely, our
approach scales in a more predictable manner with the number of neurons and
the dimensionality of the system.

##  VI Conclusion and Future Work

This paper presents a new method for evaluating the safety of linear dynamical
systems controlled by a NN within a potentially non-convex region of the state
space. Our algorithm is complete and exhibits excellent scalability properties
regarding the number of neurons and the system’s dimensionality. Our future
research will concentrate on extending this approach to nonlinear systems.
This will entail approximating the nonlinear system using a NN with PWA
activation functions, effectively converting the closed-loop system into a
linear representation.

## References

  * [1] E.T. Maddalena, C.G. da S. Moraes, G. Waltrich and C.N. Jones  “A Neural Network Architecture to Learn Explicit MPC Controllers from Data” 21st IFAC World Congress  In _IFAC-PapersOnLine_ 53.2, 2020, pp. 11362–11367  DOI: https://doi.org/10.1016/j.ifacol.2020.12.546
  * [2] Timothy P Lillicrap et al.  “Continuous control with deep reinforcement learning”  In _arXiv preprint arXiv:1509.02971_ , 2015
  * [3] K.J. Hunt, D. Sbarbaro, R. Żbikowski and P.J. Gawthrop  “Neural networks for control systems—A survey”  In _Automatica_ 28.6, 1992, pp. 1083–1112  DOI: https://doi.org/10.1016/0005-1098(92)90053-I
  * [4] Chao Huang et al.  “Reachnn: Reachability analysis of neural-network controlled systems”  In _ACM Transactions on Embedded Computing Systems (TECS)_ 18.5s ACM New York, NY, USA, 2019, pp. 1–22
  * [5] Kyle D. Julian and Mykel J. Kochenderfer  “A Reachability Method for Verifying Dynamical Systems with Deep Neural Network Controllers”, 2019  arXiv:1903.00520 [cs.SY]
  * [6] Joseph A. Vincent and Mac Schwager  “Reachable Polyhedral Marching (RPM): A Safety Verification Algorithm for Robotic Systems with Deep Neural Network Components”  In _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , 2021, pp. 9029–9035  DOI: 10.1109/ICRA48506.2021.9561956
  * [7] Hengjun Zhao, Xia Zeng, Taolue Chen and Zhiming Liu  “Synthesizing barrier certificates using neural networks”  In _Proceedings of the 23rd international conference on hybrid systems: computation and control_ , 2020, pp. 1–11
  * [8] Qingye Zhao et al.  “Verifying neural network controlled systems using neural networks”  In _Proceedings of the 25th ACM International Conference on Hybrid Systems: Computation and Control_ , 2022, pp. 1–11
  * [9] S. F.  “Set-Theoretic Methods in Control”  Birkhäuser, 2008, pp. 103–106
  * [10] Guido F Montufar, Razvan Pascanu, Kyunghyun Cho and Yoshua Bengio  “On the number of linear regions of deep neural networks”  In _Advances in neural information processing systems_ 27, 2014
  * [11] Atilim Gunes Baydin, Barak A Pearlmutter, Alexey Andreyevich Radul and Jeffrey Mark Siskind  “Automatic differentiation in machine learning: a survey”  In _Journal of Marchine Learning Research_ 18 Microtome Publishing, 2018, pp. 1–43
  * [12] Adam Paszke et al.  “Pytorch: An imperative style, high-performance deep learning library”  In _Advances in neural information processing systems_ 32, 2019
  * [13] Kurt Hornik, Maxwell Stinchcombe and Halbert White  “Multilayer feedforward networks are universal approximators”  In _Neural Networks_ 2.5, 1989, pp. 359–366  DOI: https://doi.org/10.1016/0893-6080(89)90020-8
  * [14] Arne Brondsted  “An introduction to convex polytopes”  Springer Science & Business Media, 2012
  * [15] Lingyang Chu et al.  “Exact and consistent interpretation for piecewise linear neural networks: A closed form solution”  In _Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, 2018, pp. 1244–1253
  * [16] Boris Hanin and David Rolnick  “Complexity of linear regions in deep networks”  In _International Conference on Machine Learning_ , 2019, pp. 2596–2604  PMLR