proofProof potProof of Theorem LABEL:thm

[ orcid=0000-0002-7724-8878]

[orcid=0000-0003-4201-4745] [1] 1]organization=University of Kentucky,
city=Lexington, state=Kentucky, country=USA

[cor1]Corresponding author

# Invariant Set Estimation for Piecewise Affine Dynamical Systems Using
Piecewise Affine Barrier Function

Pouya Samanipour  Hasan A.Poonawala  [<EMAIL_ADDRESS>
###### Abstract

This paper introduces an algorithm for approximating the invariant set of
closed-loop controlled dynamical systems identified using ReLU neural networks
or piecewise affine ($\mathrm{PWA}$) functions, particularly addressing the
challenge of providing safety guarantees for ReLU networks commonly used in
safety-critical applications. The invariant set of $\mathrm{PWA}$ dynamical
system is estimated using ReLU networks or its equivalent $\mathrm{PWA}$
function. This method entails formulating the barrier function as a
$\mathrm{PWA}$ function and converting the search process into a linear
optimization problem using vertices. We incorporate a domain refinement
strategy to increase flexibility in case the optimization does not find a
valid barrier function. Moreover, the objective of optimization is to maximize
the invariant set based on the current partition. Our experimental results
demonstrate the effectiveness and efficiency of our approach, demonstrating
its potential for ensuring the safety of $\mathrm{PWA}$ dynamical systems.

###### keywords:

Barrier Function Piecewise Affine Dynamics ReLU Neural Networks Forward
Invariant Set

##  1 Introduction

Safety is paramount in various applications, such as robotics and autonomous
driving. It has become increasingly common for these applications to be
deployed in diverse environments as they gain more advanced features and
autonomy. It is imperative to guarantee safety for both robots or autonomous
vehicles and the humans engaging with them [20]. Conventional control methods
often rely on pre-established models and manually designed controllers, a
paradigm that might not be optimal in the face of the complex and dynamic
environment [3].

One common approach in the literature that can guarantee safety automatically
is the model predictive controller (MPC) [24]. MPC consists of making
decisions based on a predictive model of the system [24]. Control sequences
must be optimized over a finite time horizon while considering the system’s
dynamics and constraints. In the context of the MPC controller, due to the
fact that it heavily depends on the accuracy of the predictive model,
learning-based approaches are used in order to learn the predictive model [22,
35, 21].

The other common approach for guaranteeing safety automatically is using
barrier functions (BF). To verify safety in dynamical systems, a barrier
function synthesized by the sum of squares was initially proposed by [27]. In
this case, the BF is dealing with a closed-loop dynamic with a known safe
controller. For dealing with dynamical systems with unknown safe controllers,
as described in [41], the control barrier function (CBF) was initially
introduced, followed by reciprocal and zeroing less conservative barrier
functions described in [2, 1]. CBFs are commonly employed as safety filters in
control systems. Since these methods are highly dependent on the model, the
barrier function should be robust against uncertainty [10].

Accordingly, the use of machine learning approaches to learn the model and the
control policy using data has gained significant attention in recent years [3,
14, 36, 30, 16, 11]. Learning methods have led to an increased interest in
rectified linear unit (ReLU) neural networks (NN). This is due to the
universal approximation theorem of ReLU neural networks [19, 13]. ReLU NNs are
widely used to learn dynamics and control policy [7, 34, 38]. Learning
approaches are challenged by the importance of ensuring safety during training
and testing [23]. Thus, it is necessary to develop an automated approach for
determining the forward invariance set for learned dynamics and control
policy.

Supervised learning and reinforcement learning (RL) can be used to ensure
safety. A model-free RL controller is combined with a model-based controller
using CBF in [6]. The method utilizes a Gaussian Process (GP) model to learn
unknown system dynamics and derives CBFs from the GP model. An off-policy
actor-critic RL algorithm is employed in [23] to maintain a safe policy
without prior knowledge of the dynamics of the system. The CBF formulation of
[39] accommodates a variety of uncertainties learned episodically by NN by
placing uncertainty directly into its formulation. An adaptive sampling method
is used to expand the safe region of a dynamical system with an unknown
disturbance term in [40]. Through chance-constrained optimization problems,
[17] provides a more general approach to optimizing system behavior and
ensuring system safety.

The issue of effectively addressing the synthesis of control barrier functions
for broad categories of systems remains a gap in prior research endeavors.
Nonetheless, recent investigations have introduced various methodologies [9,
31, 28, 37, 43, 32] to confront this challenge. The reachability analysis has
become a popular method for determining the safe set [9], despite potential
computational difficulties or a tendency to produce overly conservative
estimates. In [37], a support vector machine is leveraged to parameterize the
CBF in real-time, while [39] introduces an optimization-based approach to
acquiring knowledge about uncertainty in CBF constraints through expert
demonstrations. Another approach, presented in [28, 29], involves employing a
soft minimum and soft max barrier function over a finite horizon to synthesize
safe control for nonlinear systems subject to actuator constraints. The
utilization of a feedforward neural network with ReLU activation function is
proposed in [43] for barrier certificates. The verification of the barrier
certificate is addressed using Mixed Integer Linear Programming (MILP) and
Mixed Integer Quadratically Constrained Problems (MIQCP) with non-linear
terms. Furthermore, [42] discusses verifying the ReLU neural network in the
context of the barrier function. The computational complexity of the verifier,
however, remains a challenge.

This paper introduces a novel framework for determining the forward invariant
set of $\mathrm{PWA}$ dynamical systems or their equivalent ReLU neural
networks, bypassing the need for verification procedures. The contributions
can be described as follows:

  1. 1.

Introduction of a framework for determining the forward invariant set of PWA
dynamical systems or ReLU neural networks without relying on traditional
verification procedures.

  2. 2.

The BF is represented as a $\mathrm{PWA}$ function with a structure similar to
dynamical systems, with constraints expressed as linear equations.

  3. 3.

Adding a refinement step to the algorithm to enhance flexibility through the
use of higher capacity $\mathrm{PWA}$ Barrier Functions.

The paper’s organization is as follows: Section 2 provides an overview of
relevant concepts, Section 3 details the parameterization of the barrier
function and assumptions for the forward invariant set, Section 4 presents the
main algorithm, and Section 5 showcases simulation results.

##  2 Preliminaries

In this paper, we strive to estimate the invariant set for dynamical systems
described by $\mathrm{PWA}$ functions. We first introduce the concept of
positive invariant set, barrier function, and then describe the $\mathrm{PWA}$
function and required notation for the remainder of the paper.

###  2.1 Forward Invariant Set and Barrier Function Description

Let us consider the following nonlinear dynamics:

$\dot{x}=f_{cl}(x).$ (1)

We assume that the closed-loop dynamics, represented by $f_{cl}$, is locally
Lipschitz continuous. Under this assumption, for any given initial condition
$x_{0}\in\mathbb{R}^{n}$, there exists a time interval denoted as
$I(x_{0})=[0,\tau_{\text{max}})$, within which a unique solution
$x:I(x_{0})\rightarrow\mathbb{R}^{n}$ exists. This solution satisfies the
differential equation (1) with the initial condition $x_{0}$ [1].

######  Definition 1 (Forward invariant set [1])

In the case of general controlled nonlinear dynamics described by equation
(1), let us consider a super level set $C$ defined by a continuous
differentiable function $h:D\subset\mathcal{R}^{n}\rightarrow\mathcal{R}$ as
follows:

$C=\\{x\in D:h(x)\geq 0\\},$ (2)

where $D$ represents the domain of $h$. The set $C$ is said to be forward
invariant if, for every $x_{0}\in C$, the solution $x(t)$ satisfies $x(t)\in
C$ for all $x_{0}$ and all $t\in I(x_{0})$. If the set $C$ is forward
invariant, it implies that the system described by equation (1) is safe with
respect to $C$.

######  Definition 2 (Barrier function [1])

Let $C\subset D\subseteq\mathbb{R}^{n}$ be the superlevel set of a
continuously differentiable function $h(x)$. We define $h(x)$ as a Barrier
function if there exists an extended class $K_{\infty}$ function $\alpha$ such
that:

$L_{f_{cl}}h(x)\geq-\alpha(h(x)),\quad\text{for all }x\in D,$ (3)

where $L_{f_{cl}}h(x)$ represents the Lie derivative of $h$ along the closed-
loop dynamics $f_{cl}$. The equation (3) is the generalized form of Nagumo’s
theorem [25]. Equation (3) implies that $h$ can increase or decrease along the
trajectory of the system described by equation (1) within $Int(C)$.
Additionally, it ensures that $h$ remains increasing along the trajectories of
(1) on $\partial C$ as well as on the complement of $C$ within $D$.

###  2.2 Representations of Piecewise Affine Functions

First, we introduce the notation and definitions for $\mathrm{PWA}$ functions
that will be used throughout this paper, and then we describe the standard
representation of $\mathrm{PWA}$ functions 2.3. These definitions and
notations establish a consistent framework for the subsequent discussions and
analyses in the paper.

#### Notation

Given a set $S$, the set $I(S)$ is defined as the index set corresponding to
each element in $S$. The convex hull of $S$ is represented by
$\mathrm{conv}\left(S\right)$. The interior of $S$ is represented as
$\text{Int}(S)$. The boundary of $S$ is expressed as $\partial S$. The closure
of $S$ is denoted as $\overline{S}$, encompassing both $\text{Int}(S)$ and its
boundary. In the context of vectors, $A^{T}$ denotes the transpose of $A$. The
inner product of two vectors is denoted by $\langle\cdot,\cdot\rangle$.
Additionally, it is important to note that the symbols $\preceq$, $\succ$, and
$\prec$ carry the same element-wise interpretation as the symbols $\leq$, $>$,
and $<$, respectively.

###### Definition 3

This paper introduces the concept of a partition $\mathcal{P}$, which is
defined as a collection of subsets $\\{X_{i}\\}_{i\in I(\mathcal{P})}$, where
each $X_{i}$ represents a closed subset of $\mathcal{R}^{n}$ for all $i\in
I(\mathcal{P})$. In partition $\mathcal{P}$,
$\text{Dom}(\mathcal{P})=\cup_{i\in I(\mathcal{P})}X_{i}$ and
$\mathrm{Int}\left(X_{i}\right)\cap\mathrm{Int}\left(X_{j}\right)=\emptyset$
for $i\neq j$. As a result of these notations and definitions, we can proceed
with the introduction of the $\mathrm{PWA}$ functions.

###  2.3 Piecewise Affine Functions

We represent a $\mathrm{PWA}(x)$ using an explicit parameterization with
respect to a partition $\mathcal{P}=\\{X_{i}\\}_{i\in I(\mathcal{P})}$, a
collection of matrices $\mathbf{A}_{\mathcal{P}}=\\{A_{i}\\}_{i\in
I(\mathcal{P})}$, and vectors $\mathbf{a}_{\mathcal{P}}=\\{a_{i}\\}_{i\in
I(\mathcal{P})}$. The parameterization is defined as follows:

$\displaystyle\mathrm{PWA}(x)$ $\displaystyle=A_{i}x+a_{i},\textrm{ if }x\in
X_{i},\text{ where}$ (4) $\displaystyle X_{i}$
$\displaystyle=\\{x\in\mathbb{R}^{n}\colon E_{i}x+e_{i}\succeq 0\\}.$

This parameterized formulation concisely represents the piecewise affine
function $\mathrm{PWA}(x)$ based on the given partition $\mathcal{P}$ and the
associated matrices and vectors.

###### Assumption 1

In this paper, we specifically assume that all $\mathrm{PWA}$ functions
presented herein, with their explicit form, satisfy the continuity constraints
described in [26, 15]).

###### Assumption 2

In our study, all cells are assumed to be bounded. The vertex representation
of the $\mathrm{PWA}$ dynamics can therefore be used. In a cell $X_{i}$, the
vertices represent the faces of 0 dimensions. Hence, $X_{i}$ can be
represented mathematically as the convex hull of its vertex set:

$X_{i}=\mathrm{conv}\left(\mathcal{F}_{0}(X_{i})\right),$ (5)

where $\mathcal{F}_{0}(X_{i})$ refers to the set of vertices that belong to
the cell $X_{i}$ [12].

###### Assumption 3

This paper assumes that the

$\textbf{0}\in X_{i}\Rightarrow\textbf{0}\in\mathcal{F}_{0}(X_{i}).$

##  3 Problem formulation

The objective of this paper is to identify a forward invariant set for
dynamical systems characterized by the equation:

$\dot{x}=\mathrm{PWA}(x)$ (6)

where $x$ belongs to the $\mathcal{D}\subset\mathcal{R}^{n}$ representing the
state variables and the term $\mathrm{PWA}$ denotes a piecewise affine
function, which was elaborated in Section 2.3. To align with the notation used
for the concepts of forward invariance and barrier functions, we define
$\mathcal{D}$ as:

$\mathcal{D}\triangleq\\{x\in\mathcal{R}^{n}:x\in\mathrm{Dom}(\mathcal{P})\\}.$
(7)

It is important to note that no assumption is made regarding the forward
invariance of $\mathcal{D}$. Specifically, this study is concerned with
continuous $\mathrm{PWA}$ functions based on polytopes. The objective of this
section is to identify the invariant set for (6) which represents a closed-
loop system under a designed controller for which an invariant set exists, but
is initially unknown.Various methods can be employed to design the controller,
such as traditional control techniques and explicit model predictive control
(eMPC).

###  3.1 Piecewise affine barrier function:

This section introduces the formulation of the BF as a $\mathrm{PWA}$ function
for the dynamical system described by equation (6), considering a given
partition $\mathcal{P}$. We consider the following index set.

$\displaystyle I_{\partial\mathcal{D}}$ $\displaystyle=\\{i\in I\colon
X_{i}\cap\partial\mathcal{D}\neq\emptyset\\}.$ (8)

The set $I_{\partial\mathcal{D}}$ represents the indices of cells containing
the boundary of $\mathcal{D}$. These cells have points in common with the
boundary of the domain $\mathcal{D}$. In order to determine which vertices are
in the $\partial\mathcal{D}$ and which vertices are in the $Int(\mathcal{D})$,
we can use the sets $I_{b}$ and $I_{int}$ as follows:

$\displaystyle I_{b}$ $\displaystyle=\\{(i,k):v_{k}\in\partial\mathcal{D},i\in
I_{\partial D},v_{k}\in\mathcal{F}_{0}(X_{i})\\}$ (9) $\displaystyle I_{int}$
$\displaystyle=\\{(i,k):v_{k}\notin\partial\mathcal{D},i\in
I(\mathcal{P}),v_{k}\in\mathcal{F}_{0}(X_{i})\\}.$ (10)

The set $I_{b}$ and $I_{int}$ consists of pairs of indices, where the first
element denotes the cell index $i$, and the second element represents the
vertex index $k$. $I_{b}$ represents the indices of the vertices on the
boundary and their associated cells, while $I_{int}$ represents the indices of
vertices in the interior of the partition and their associated cells. In
Figure 1, samples for these index cells are provided.

Figure 1: A sample partition of a $\mathrm{PWA}$ dynamics. Cells $X_{1}$
,$X_{2}$, $X_{3}$ and $X_{4}$ share points with the boundary, hence,
$\\{1,2,3,4\\}\in I_{\partial\mathcal{D}}$. While cells $X_{5}$,$X_{6}$, and
$X_{7}$ do not share any points with the boundary. As a result
$\\{5,6,7\\}\notin I_{\partial\mathcal{D}}$. Furthermore, $v_{1},v_{2}$ are
vertices on the boundary, while $v_{3}$ and $v_{4}$ are not on the boundary.
Hence, $\\{(1,1),(1,2),(2,2)\\}\in I_{b}$ and $\\{(6,3),(6,4)\\}\in I_{int}$.

The main idea in this paper is to formulate the barrier function using a
$\mathrm{PWA}$ formulation for the cell $X_{i}$ as follows:

$h_{i}(x)=p_{i}^{T}x+q_{i}\quad\text{for}\quad i\in I(\mathcal{P})$ (11)

where $p_{i}\in\mathcal{R}^{n}$ and $q_{i}$ is a scalar. The function $h(x)$
is continuous and differentiable within the interior of the cell. This allows
us to compute the derivative of the candidate Barrier function along the
dynamics described by (1) within the interior of the cell $X-\\{0\\}$. The
derivative is expressed as:

$\mathcal{L}_{f}h=\langle\nabla h,f(x)\rangle,$ (12)

where $\nabla h$ represents the gradient of $h(x)$, and $\mathcal{L}$
represents the Lie derivative.

###### Assumption 4

Let’s assume the candidate barrier function $h_{i}(x)=p_{i}^{T}x+q_{i}$ and
the dynamics equation given by $\dot{x}=A_{i}x+a_{i}$ for all $x\in X_{i}$.
Then the derivative of the barrier function along the trajectories at
$\mathbf{x}^{\prime}$, when $h(\mathbf{x}^{\prime})$ is differentiable can be
computed as:

$\dot{h}=p_{i}^{T}(A\mathbf{x}^{\prime}+a).$ (13)

This paper has the following objectives:

  * •

develop an optimization problem to find the invariant set for the (6).

  * •

Maximize the forward invariant set obtained from $\mathrm{PWA}$ barrier
function.

To achieve these goals, we must form the constraints based on the requirements
of the barrier function as described in (3). In order to meet the goal of
finding a maximized invariant set using a barrier function formulated as a
$\mathrm{PWA}$ function, the following propositions and theorems are
presented.

###### Proposition 1

The following are equivalent:

  1. (a)

$C\subseteq\mathcal{D}$

  2. (b)

$p_{i}^{T}v_{k}+q_{i}\leq 0,\quad\forall(i,k)\in I_{b}$

######  Proof 3.1.

Regarding the definition $C$ in (2) and (a), if $C=D$ we can conclude that

$p_{i}^{T}v_{k}+q_{i}=0,\quad(i,k)\in I_{b}.$

Furthermore, if $C\subset D$, we may find points $x_{d}\in D$ where
$x_{d}\notin C$. In light of Assumption 2 and parametrization of the barrier
function, (11), it can be concluded that there is a vertex $v_{d}$ at which
$p_{i}^{T}v_{d}+q_{i}<0$. Accordingly, we can conclude (b).

It is essential that the obtained invariant set must be a subset of
$\mathcal{D}$. Using Proposition 1, we can satisfy this requirement by
maintaining the barrier function non-positive on the $\partial\mathcal{D}$ of
the $\mathcal{D}$. However, there is no guarantee that such $p_{i}$ and
$q_{i}$ values exist. In order to ensure the feasibility of the optimization
problem, b must be relaxed to:

$p_{i}^{T}v_{k}+q_{i}-\tau_{b_{i}}\leq 0,\hskip 1.99997pt(i,k)\in
I_{b},v_{k}\in\mathcal{F}_{0}(X_{i}),$ (14)

where $\tau_{b_{i}}$ represents a positive slack variable introduced to ease
the constraint. This slack variable allows us to identify cells $X_{i}$ that
might necessitate a higher capacity $\mathrm{PWA}$ function. The process for
increasing the capacity of the $\mathrm{PWA}$ function is explained in section
4.2.

######  Proposition 1.

The candidate $\mathrm{PWA}$ barrier function is continuous if and only if the
following condition is met:

$\displaystyle h_{i}(v_{k})=$ $\displaystyle h_{j}(v_{k}),\quad i\neq j\in
I(\mathcal{P}),$ $\displaystyle
v_{k}\in\mathcal{F}_{0}(X_{i})\cap\mathcal{F}_{0}(X_{j}).$ (15)

######  Proof 3.2.

Clearly, if the barrier function is continuous, then the barrier function
value must be the same at all common points, such as vertices. Conversely, if
(1) is satisfied due to the parametrization of the barrier function and
convexity of the common faces between any two adjacent cells, we can conclude
that the candidate barrier function is continuous on the boundary of the
cells. As stated in (11), the candidate barrier function is continuous
throughout the cell’s interior. Therefore, the candidate $\mathrm{PWA}$
barrier function is continuous.

######  Proposition 2.

The following are equivalent if we assume $\alpha(x)=\tilde{\alpha}x$, where
$\tilde{\alpha}>0$:

$\displaystyle(a)$ $\displaystyle\quad\frac{\partial h}{\partial
x}\dot{x}+\alpha(h(x))\geq 0\quad\forall x\in\mathcal{D}$ $\displaystyle(b)$
$\displaystyle\quad
p_{i}^{T}(A_{i}v_{k}+a_{i})+\tilde{\alpha}(p_{i}^{T}v_{k}+q_{j})\geq 0,\forall
i\in I(\mathcal{P}),$ $\displaystyle\quad k\in\mathcal{F}_{0}(X_{i}).$ (16)

######  Proof 3.3.

Since the cells are convex polytopes, we need only check (a) at all vertices
to ensure that all points within the interior of the cells satisfy (a).
Therefore, (b) can be derived directly from (a) by using (6), (13) and
assumption $\alpha(x)=\tilde{\alpha}x$ for all vertices. These conditions also
account for the non-differentiable nature of $h$ at boundaries of cells,
however we suppress the fairly involved details, resulting in an abuse of
notation. Our prior work in [34, 33] discusses how the theory of set-valued
Lie derivatives and continuous PWA functions justifies this claim. If the PWA
dynamics involved switching, this claim would no longer hold.

######  Theorem 3.4.

A barrier function (11) is considered to be a valid barrier function if
Proposition 1, Proposition 1 and Proposition 2 hold.

######  Proof 3.5.

It is guaranteed by the first proposition that the obtained invariant set is a
subset of the partition domain $\mathcal{D}$. The barrier function will be
continuous as we used the proposition 1, and it will be differentiable within
the cell’s interior. Furthermore, Proposition 2 guarantees the generalized
Nagumo’s condition.

Using the propositions presented so far, the first objective has been
achieved. In order to achieve the second goal, maximizing the positive
invariant set with respect to the current partition, we need to ensure that as
many vertices in the interior of the domain as possible are also inside the
invariant set or, in other words, $h_{i}(v_{k})\geq 0$ where $(i,k)\in
I_{int}$.

###### Assumption 5

To achieve the second objective, the following inequality can be utilized:

$p_{i}^{T}v_{k}+q_{i}+\tau_{int_{i}}>0,\hskip 1.99997pt(i,k)\in I_{int}.$ (17)

In (17), $\tau_{int_{i}}$ represents a slack variable introduced to relax the
constraint. By incorporating this slack variable, some flexibility is
provided. The importance of this flexibility arises from the fact that not all
interior cells may be considered part of the invariant set. Including the
slack variable makes it possible to identify interior cells that need to be
refined to achieve the desired forward invariant set. In 4.1, it will be shown
how this constraint will help us to maximize the invariant set.

By integrating these three Propositions and (17), we are able to construct a
systematic approach to searching for the invariant set using the barrier
function.

##  4 Main Algorithm

According to the assumptions and propositions made in the previous section,
the main algorithm will be constructed. As a first step, the optimization
problem is formulated, and then the refinement process is described to
increase the capacity of the $\mathrm{PWA}$ function. In the end, we will
discuss the search process.

###  4.1 Forming optimization problem

We can create an optimization problem to derive the barrier function using
equations (14) through (17). The cost function can be defined as follows:

$\mathcal{J}=\sum_{i=1}^{M}\tau_{b_{i}}+\sum_{i=1}^{N}\tau_{int_{i}}$ (18)

$M$ and $N$ are the number of elements in $I_{b}$ and $I_{int}$ respectively.

The following represents the ultimate form of this optimization problem:

$\displaystyle\min_{p_{i},q_{i},\tau_{int_{i}},\tau_{b_{i}}}\quad\mathcal{J}$
(19) Subject to: $\displaystyle
p_{i}^{T}v_{k}+q_{i}-\tau_{b_{i}}\leq-\epsilon_{1},\hskip
1.99997pt\forall(i,k)\in I_{b}$ $\displaystyle
p_{i}^{T}v_{k}+q_{i}+\tau_{int_{i}}\geq\epsilon_{2},\hskip
1.99997pt\forall(i,k)\in I_{int}$ $\displaystyle
p_{i}^{T}(A_{i}v_{k}+a_{i})+\alpha(p_{i}^{T}v_{k}+q_{i})\geq\epsilon_{3},\hskip
1.99997pt\forall i\in I(\mathcal{P})$ $\displaystyle
h_{i}(v_{k})=h_{j}(v_{k}),\forall
v_{k}\in\mathcal{F}_{0}(X_{i})\cap\mathcal{F}_{0}(X_{j})$
$\displaystyle\tau_{b_{i}}\geq 0,\forall i\in I_{\partial\mathcal{D}}$

where $\epsilon_{1},\epsilon_{2},\epsilon_{3}>0$.

######  Lemma 4.6.

Optimization problem (19) is always feasible.

######  Proof 4.7.

As a result of the construction of the optimization problem, this result is
obtained.

Although the optimization problem is always feasible, The solution will be the
valid control barrier function if and only if:

$\sum_{i=1}^{M}\tau_{b_{i}}=0.$ (20)

If a non-zero $\tau_{b_{i}}$ exists, it indicates that the obtained solution
is not a valid barrier function. In other words, based on the current
solution, no forward invariant set $C$ satisfies the condition
$C\subseteq\mathcal{D}$. In order to address this issue, a refinement process
is proposed. Refinement involves dividing each cell into smaller subcells,
allowing for the construction of higher-capacity $\mathrm{PWA}$ barrier
functions. The details of the refinement procedure will be discussed in the
following sections.

######  Lemma 4.8.

The valid invariant set obtained from the optimization problem (19) is the
largest possible invariant set with the current partition and $\mathrm{PWA}$
barrier function.

######  Proof 4.9.

Let’s consider that there exists an invariant set $C^{\prime}$ where $C\subset
C^{\prime}$. As a result, there are points in a number of cells that belong to
$C^{\prime}$ but are not included in $C$. For simplicity, let’s assume the
difference of invariant set $C^{\prime}$ and $C$ is only at cell $X_{1}$ at
some points such as $x_{1}$. This assumption does not pose any problems for
generalization.At $x_{1}$, the barrier function corresponding to $C^{\prime}$,
$h^{\prime}(x)$, is positive, whereas the barrier function corresponding to
$C$, $h(x)$, is negative. Therefore, we can conclude that
$h_{1}^{\prime}(x_{1})>h_{1}(x_{1})$. Due to the parametrization of the
barrier function as a $\mathrm{PWA}$ function, and
$h_{1}^{\prime}(x_{1})>h_{1}(x_{1})$, there exists a
$v_{k}\in\mathcal{F}_{0}(X_{i})$ where $h_{1}^{\prime}(v_{k})>h_{1}(v_{k})$.
Regarding (17), there exist a $\tau^{\prime}_{int_{1}}<\tau_{int_{1}}$, where
$\tau^{\prime}_{int_{1}}$ and $\tau_{int_{1}}$ are the slack variables in (17)
for $h_{1}^{\prime}(x)$ and $h_{1}(x)$ respectively. Moreover, according to
the assumption $C\subset C^{\prime}$ we can conclude that the slack variable
at all the other interior cells is

$\tau^{\prime}_{int_{i}}\leq\tau_{int_{i}}\quad\forall i\in I_{int}-\\{1\\}.$

Therfore, $\sum_{i=1}^{N}\tau_{int_{i}}>\sum_{i=1}^{N}\tau^{\prime}_{int_{i}}$
which is contrary to the minimization of the slack variable $\tau_{int_{i}}$
in the optimization problem (19). As a result, there is no invariant set
associated with the current partition where $C\subset C^{\prime}$.

0: $\mathrm{PWA}(x)$, Vertices

Solve Optimization Problem (19) with the initial $\mathrm{PWA}$ dynamics.

while $\Sigma_{i=1}^{N}\tau_{b_{i}}\neq 0$ and $\textit{computational-
time}<3600s$ do

for $\tau_{b_{i}}\neq 0$ and $\tau_{int_{i}}\neq 0$ do

1- Refine $X_{i}$ using vector field refinement Section 4.2.

end for

Solve Optimization Problem (19) with the refined $\mathrm{PWA}$ dynamics

end while

return $h(x)=p_{i}^{T}x+q_{i}\quad\forall i\in I(\mathcal{P})$.

Algorithm 1 Barrier function search using vertices with known control

###  4.2 Refinement

A refinement of the current partition aims to increase the flexibility of the
Barrier function search process. In mathematical terms, given two partitions
$\mathcal{P}=\\{Y_{i}\\}_{i\in I}$ and $\mathcal{R}=\\{Z_{j}\\}_{j\in J}$ of a
set
$S=\mathrm{Dom}\left(\mathcal{P}\right)=\mathrm{Dom}\left(\mathcal{R}\right)$,
we say that $\mathcal{R}$ is a refinement of $\mathcal{P}$ if $Z_{j}\cap
Y_{i}\neq\emptyset$ implies that $Z_{j}\subseteq Y_{i}$. To increase
flexibility, we introduce a refinement on the cell $X_{i}$ with a nonzero
slack variable, $\tau_{b_{i}}$. As a result of the refinement of the cell
$X_{i}$, consider that two new sub-cells will be created, $X_{i_{1}}$ and
$X_{i_{2}}$. Each sub-cell can be parameterized with a $\mathrm{PWA}$ Barrier
function, denoted as $h_{i_{1}}=p_{i_{1}}^{T}x+q_{i_{1}}$ and
$h_{i_{2}}=p_{i_{2}}^{T}x+q_{i_{2}}$, respectively. By dividing the cell into
smaller sub-cells and introducing separate parameterizations, the candidate
Barrier function for cell $X_{i}$ achieves a higher capacity $\mathrm{PWA}$
function, thereby increasing flexibility in the Barrier function search
process. This paper’s refinement technique includes two steps: choosing a new
vertex and using Delaunay triangulation to form subcells. The first step is
choosing a new vertex on the boundary of the considered cell with a nonzero
slack variable using the vector field approach [33]. This method is designed
to decrease the variation of the vector field of the new subcells. Therefore,
the new vertex is chosen in a way where its vector field will be the angle
bisector of the edge, the faces of 1 dimensions [12], with the biggest
variation of the vector field. Based on the newly appointed vertex, subcells
will be formed using Delaunay triangulation. A Delaunay triangulation will
also be helpful in maintaining continuity. For more detail, please see [33].

###  4.3 Search Algorithm

The pseudo-code for the search algorithm is presented in Algorithm 1. The
objective is to identify a valid Barrier function on the initial partition by
performing an optimization, as given in (19). In cases where a valid solution
is not obtained initially, the algorithm proceeds to refine the partition
defining the barrier function. This refinement is directed towards cells
exhibiting $\tau_{b_{i}}\neq 0$ or $\tau_{int_{i}}\neq 0$. Subsequently, the
selected cells are subdivided into smaller sub-cells, and the optimization
process is repeated on the refined partition. The iterative refinement process
continues until a valid solution is found or the maximum computation time is
reached.

As a result of the systematic approach outlined by the algorithm, the
partition can be refined in a methodical manner while actively searching for a
valid Barrier function. Using this iterative approach increases the capacity
of the Barrier function, although it may significantly increase the size the
optimization problem.

In this search algorithm, the next important issue is how to select
$\tilde{\alpha}$ in (2). The literature does not provide a systematic method
for selecting the $\kappa_{\infty}$ function in (2). Choosing
$\tilde{\alpha}=0$ will force all level sets of the barrier function to be
invariant, resulting in a conservative result. This paper uses the bisection
method to find the least conservative $\tilde{\alpha}$. Bisection is performed
on a predefined interval $[\tilde{\alpha}_{1},\tilde{\alpha}_{2}]$. The
Algorithm 1 will be tested by the $\tilde{\alpha}_{1}$ and
$\tilde{\alpha}_{2}$. Upon bisection, the greatest $\tilde{\alpha}$ that
satisfies the search algorithm for finding the valid barrier function is
selected. This is done by selecting $\tilde{\alpha}_{2}$ so that it violates
the algorithm’s computation time limit and selecting
$\tilde{\alpha}_{1}=10^{-3}$ so that it is conservative but still within the
algorithm’s operational time limit. Please note that the search interval will
depend on the nature of the problem, and we may need to extend the search
interval in order to find the least conservative invariant set.

##  5 Results and simulations

We provide three examples to illustrate the performance of the Barrier
function search algorithm presented in Algorithm 1. All computations are
implemented using Python 3.9 and the Mosek optimization package [4] on a
computer with a 2.1 GHz processor and 8 GB RAM. A tolerance of $10^{-8}$ is
used to determine if a number is nonzero. Moreover, the examples employ a
tolerance of ${\epsilon_{1}=\epsilon_{2}=\epsilon_{3}=10^{-4}}$. These
examples demonstrate the proposed algorithm’s effectiveness and efficiency in
finding valid $\mathrm{PWA}$ Barrier functions or the forward invariant set,
highlighting the algorithm’s performance in different scenarios.

######  Example 5.10 (Inverted Pendulum [28]).

The inverted pendulum system can be described by the following equations [28]:

$\displaystyle\dot{x}_{1}$ $\displaystyle=x_{2}$ $\displaystyle\dot{x}_{2}$
$\displaystyle=\sin(x_{1})+u$

where $x_{1}$ represents the angle of the pendulum and $x_{2}$ is its angular
velocity. The control input $u$ is given by $u=3x_{1}+3x_{2}$ and is subject
to saturation between the lower bound $-1.5$ and the upper bound $1.5$. The
domain $\mathcal{D}$ of the system is defined as
$\mathcal{D}=\pi-||x||_{\infty}\geq 0$, where $||x||_{\infty}$ represents the
infinity norm of the state vector $x$. We approximate the right hand side of
the dynamics using a ReLU neural network which is equivalent to $PWA(x)$. To
train a neural network for this system, a dataset of $10000$ samples is used.
The neural network is a single-hidden layer ReLU network with eight neurons.
The training is performed using PyTorch 2.0.0 with ADAM optimization.

For this problem, we have set $\tilde{\alpha}=0.2$. Algorithm 1 started with
34 cells and ended up with 680 due to the refinement. The forward invariant
set is approximated, as shown in Fig.2. Moreover, A sample trajectory with the
initial conditions $x_{0}=[-1,3]^{T}$ is shown in Figure 3. Additionally, the
value of the barrier function for the specified trajectory is shown in Fig.4.
Details about this example can be found in Table.1. The computation time is
based on the average time it takes to find the solution after it has been run
ten times in Table.1.

Figure 2: Vector fields for the Inverted Pendulum identified using ReLU NN and
its estimated Invariant set Figure 3: $x_{1}$ and $x_{2}$ for a sample
trajectory with initial condition $x_{0}=[-1,3]^{T}$ Figure 4: $h(x)$ along
the sample trajecty. Since the trajectory starts from outside the invariant
set, the value of $h(x)$ is negative. As soon as they enter the invariant set,
the value of $h(x)$ becomes positive.

######  Example 5.11 (Explicit model-predictive controller [33]).

To formulate a control policy for a discrete-time linear dynamical system, we
employ Model Predictive Control (MPC) following the methodology outlined in
[33]. The dynamics of the considered system are governed by the following
equations:

$\displaystyle x_{t+1}=\begin{bmatrix}1&1\\\
0&1\end{bmatrix}x_{t}+\begin{bmatrix}1\\\ 0.5\end{bmatrix}u_{t}$ (21)

To address the MPC problem, we adhere to the specifications presented in [33],
encompassing the stage cost, actuator constraints, and state constraints. The
explicit MPC controller is implemented using the MPT3 toolbox [18] in the
Matlab environment. A discretization step of 0.01 seconds is employed to
convert the discrete-time to the continuous-time dynamics.

The effectiveness of the proposed algorithm for the continuous form of (21)
with the explicit MPC controller is validated. The trail and error method is
used in this example in order to obtain $\tilde{\alpha}=4.5$ as the least
conservative value. In Fig.5, the invariant set of the 2-D MPC is depicted,
providing a visual representation of the system’s behavior. Additionally,
Fig.7 illustrates the value of the barrier function along sample trajectories
for the 2-D MPC. Notably, the algorithm initiates with 26 cells and undergoes
four iterations of refinement, resulting in an increased number of cells
(136).

Figure 5: 2-D mpc vector fields, invariant set Figure 6: The states for the
sample trajectory with the initial condition $x_{0}=[4,-3.5]^{T}$. Figure 7:
$h(x)$ along the sample trajectory for the 2-D MPC.

######  Example 5.12 (4-D Example [5]).

In this example, we explore a 4-D Model Predictive Control (MPC) scenario
based on the formulation presented in [5]. The system dynamics are given by:

$\displaystyle x_{t+1}=$
$\displaystyle\begin{bmatrix}0.4346&-0.2313&-0.6404&0.3405\\\
-0.6731&0.1045&-0.0613&0.3400\\\ -0.0568&0.7065&-0.086&0.0159\\\
0.3511&0.1404&0.2980&1.0416\end{bmatrix}x_{t}+$ (22)
$\displaystyle\begin{bmatrix}0.4346,-0.6731,-0.0568,0.3511\end{bmatrix}^{T}u_{t}.$

This example inherits the key parameters from [5], such as a state constraint
of $\lVert x\rVert_{\infty}\leq 4$, an input constraint of $\lVert
u\rVert_{\infty}\leq 1$, a prediction horizon of $T=10$, a stage cost of
$Q=10I$, and $R=1$. The design of the MPC controller incorporates a state
constraint $\lVert x\rVert_{\infty}\leq 4$, an input constraint $\lVert
u\rVert_{\infty}\leq 1$. The terminal cost utilizes $P_{\infty}$, a solution
to the Riccati equation defined by $(A,B,Q,R)$, and the terminal set is
selected to maximize the positive invariant set.

We find the invariant set of a continuous-time linear system corresponding to
(22) when using an explicit MPC controller for feedback. We assume that (22)
corresponds to a sampling time of $0.01$s. The closed loop is a $\mathrm{PWA}$
dynamical system with 193 cells. We employ vector field refinement techniques
to search for the continuous $\mathrm{PWA}$ Barrier function using Algorithm
1. In this example, we found $\tilde{\alpha}=9.7$ by trial and error. With the
vector field refinement, the algorithm successfully converges, finding the
solution in 2886 seconds and generating 17689 cells.

Example |  Initial number of cells |  final number of cells |  Computational time(s)
---|---|---|---
Inverted Pendulum |  34 |  680 |  10.3
2-D MPC |  26 |  136 |  1.21
4-D MPC |  192 |  17689 |  2886
Table 1: Summary of results for finding the valid barrier function using the
proposed methods Algorithm. 1.

#### Limitations

In order to find the $\mathrm{PWA}$ barrier function for the $\mathrm{PWA}$
dynamical system with $n$ state variables, we require $(n+1)n_{r}$
optimization variables, where $n_{r}$ represents the number of cells. In this
paper, we added the $M+N$ slack variable as described in 4.1. It has been
demonstrated in [33] that $n$ influences the number of cells $n_{r}$. As
described in [8], the computational complexity of a linear optimization
problem is $O(n_{p}^{3/4}log(\frac{n_{p}}{\epsilon}))$ where $n_{p}$ is the
number of variables and $\epsilon$ is the accuracy. As a result, $n$ as the
dimension of the state space directly impacts the optimization process’s
computational time. When applying our algorithm to various systems, it is
essential to consider these complexities and challenges.

##  6 Conclusion

This paper presents a novel computational framework to estimate the invariant
set for $\mathrm{PWA}$ dynamical systems or their equivalent ReLU neural
networks, addressing the challenge of ensuring safety guarantees. The key
contribution is formulating the barrier function as a $\mathrm{PWA}$ function
and casting the search process as a linear optimization problem. The method
intelligently refines the domain when the optimization fails, enhancing the
search process’s flexibility. The optimization problem aims to maximize the
invariant set, improving safety guarantees. The proposed method effectively
estimates the invariant set through non-trivial examples while maintaining
computational efficiency. The results demonstrate the approach’s capability
and applicability in providing safety guarantees for $\mathrm{PWA}$ dynamics
and ReLU neural networks in safety-critical systems. Further exploration and
extension of this approach hold potential for more sophisticated techniques in
addressing safety concerns in neural network-based systems.

## References

  * Ames et al. [2019] Ames, A.D., Coogan, S., Egerstedt, M., Notomista, G., Sreenath, K., Tabuada, P., 2019\.  Control barrier functions: Theory and applications, in: 2019 18th European control conference (ECC), IEEE. pp. 3420–3431.
  * Ames et al. [2014] Ames, A.D., Grizzle, J.W., Tabuada, P., 2014.  Control barrier function based quadratic programs with application to adaptive cruise control, in: 53rd IEEE Conference on Decision and Control, IEEE. pp. 6271–6278.
  * Anand et al. [2021] Anand, A., Seel, K., Gjærum, V., Håkansson, A., Robinson, H., Saad, A., 2021\.  Safe learning for control using control lyapunov functions and control barrier functions: A review.  Procedia Computer Science 192, 3987–3997.
  * ApS [2022] ApS, M., 2022.  Mosek optimizer api for python.  Version 9, 6–4.
  * Chen et al. [2021] Chen, S., Fazlyab, M., Morari, M., Pappas, G.J., Preciado, V.M., 2021.  Learning lyapunov functions for hybrid systems, in: Proceedings of the 24th International Conference on Hybrid Systems: Computation and Control, pp. 1–11.
  * Cheng et al. [2019] Cheng, R., Orosz, G., Murray, R.M., Burdick, J.W., 2019\.  End-to-end safe reinforcement learning through barrier functions for safety-critical continuous control tasks, in: Proceedings of the AAAI conference on artificial intelligence, pp. 3387–3395.
  * Dai et al. [2021] Dai, H., Landry, B., Yang, L., Pavone, M., Tedrake, R., 2021\.  Lyapunov-stable neural-network control.  arXiv preprint arXiv:2109.14152 .
  * El Ghami et al. [2012] El Ghami, M., Guennoun, Z., Bouali, S., Steihaug, T., 2012\.  Interior-point methods for linear optimization based on a kernel function with a trigonometric barrier term.  Journal of Computational and Applied Mathematics 236, 3613–3623.
  * Gillula and Tomlin [2012] Gillula, J.H., Tomlin, C.J., 2012\.  Guaranteed safe online learning via reachability: tracking a ground target using a quadrotor, in: 2012 IEEE International Conference on Robotics and Automation, IEEE. pp. 2723–2730.
  * Hamdipoor et al. [2023] Hamdipoor, V., Meskin, N., Cassandras, C.G., 2023.  Safe control synthesis using environmentally robust control barrier functions.  European Journal of Control 74, 100840\.  URL: https://www.sciencedirect.com/science/article/pii/S0947358023000699, doi:https://doi.org/10.1016/j.ejcon.2023.100840. 2023 European Control Conference Special Issue.
  * Harapanahalli et al. [2023] Harapanahalli, A., Jafarpour, S., Coogan, S., 2023.  Forward invariance in neural network controlled systems.  IEEE Control Systems Letters 7, 3962–3967.  doi:10.1109/LCSYS.2023.3341980.
  * Henk et al. [2017] Henk, M., Richter-Gebert, J., Ziegler, G.M., 2017.  Basic properties of convex polytopes, in: Handbook of discrete and computational geometry. Chapman and Hall/CRC, pp. 383–413.
  * Huang [2020] Huang, C., 2020.  Relu networks are universal approximators via piecewise linear or constant functions.  Neural Computation 32, 2249–2278.
  * Jebellat et al. [2021] Jebellat, I., Pishkenari, H.N., Jebellat, E., 2021.  Training microrobots via reinforcement learning and a novel coding method, in: 2021 9th RSI International Conference on Robotics and Mechatronics (ICRoM), IEEE. pp. 105–111.
  * Johansson [1999] Johansson, M., 1999.  Piecewise linear control systems.  Ph.D. thesis. Ph. D. Thesis, Lund Institute of Technology, Sweden.
  * Kargar et al. [2023] Kargar, M., Zhang, C., Song, X., 2023.  Integrated optimization of powertrain energy management and vehicle motion control for autonomous hybrid electric vehicles.  IEEE Transactions on Vehicular Technology .
  * Khojasteh et al. [2020] Khojasteh, M.J., Dhiman, V., Franceschetti, M., Atanasov, N., 2020\.  Probabilistic safety constraints for learned high relative degree system dynamics, in: Learning for Dynamics and Control, PMLR. pp. 781–792.
  * Kvasnica et al. [2004] Kvasnica, M., Grieder, P., Baotić, M., Morari, M., 2004\.  Multi-parametric toolbox (mpt), in: Hybrid Systems: Computation and Control: 7th International Workshop, HSCC 2004, Philadelphia, PA, USA, March 25-27, 2004. Proceedings 7, Springer. pp. 448–462.
  * Leshno et al. [1993] Leshno, M., Lin, V.Y., Pinkus, A., Schocken, S., 1993\.  Multilayer feedforward networks with a nonpolynomial activation function can approximate any function.  Neural networks 6, 861–867.
  * Liu et al. [2023] Liu, S., Liu, C., Dolan, J., 2023.  Safe control under input limits with neural control barrier functions, in: Conference on Robot Learning, PMLR. pp. 1970–1980.
  * Makdesi et al. [2023] Makdesi, A., Girard, A., Fribourg, L., 2023.  Safe learning-based model predictive control using the compatible models approach.  European Journal of Control 74, 100849\.  URL: https://www.sciencedirect.com/science/article/pii/S094735802300078X, doi:https://doi.org/10.1016/j.ejcon.2023.100849. 2023 European Control Conference Special Issue.
  * Manzano et al. [2020] Manzano, J.M., Limon, D., de la Peña, D.M., Calliess, J.P., 2020\.  Robust learning-based mpc for nonlinear constrained systems.  Automatica 117, 108948\.
  * Marvi and Kiumarsi [2021] Marvi, Z., Kiumarsi, B., 2021\.  Safe reinforcement learning: A control barrier function optimization approach.  International Journal of Robust and Nonlinear Control 31, 1923–1940.
  * Morari and Lee [1999] Morari, M., Lee, J.H., 1999\.  Model predictive control: past, present and future.  Computers & chemical engineering 23, 667–682.
  * Nagumo [1942] Nagumo, M., 1942.  Über die lage der integralkurven gewöhnlicher differentialgleichungen.  Proceedings of the Physico-Mathematical Society of Japan. 3rd Series 24, 551–559.
  * Poonawala et al. [2021] Poonawala, H.A., Lauffer, N., Topcu, U., 2021.  Training classifiers for feedback control with safety in mind.  Automatica 128, 109509\.
  * Prajna et al. [2007] Prajna, S., Jadbabaie, A., Pappas, G.J., 2007.  A framework for worst-case and stochastic safety verification using barrier certificates.  IEEE Transactions on Automatic Control 52, 1415–1428.
  * Rabiee and Hoagg [2023a] Rabiee, P., Hoagg, J.B., 2023a.  Soft-minimum and soft-maximum barrier functions for safety with actuation constraints.  arXiv preprint arXiv:2305.10620 .
  * Rabiee and Hoagg [2023b] Rabiee, P., Hoagg, J.B., 2023b.  Soft-minimum barrier functions for safety-critical control subject to actuation constraints, in: 2023 American Control Conference (ACC), pp. 2646–2651.  doi:10.23919/ACC55779.2023.10156245.
  * Rajoli et al. [2023] Rajoli, H., Lotfi, F., Atyabi, A., Afghah, F., 2023\.  Triplet loss-less center loss sampling strategies in facial expression recognition scenarios, in: 2023 57th Annual Conference on Information Sciences and Systems (CISS), IEEE. pp. 1–6.
  * Robey et al. [2020] Robey, A., Hu, H., Lindemann, L., Zhang, H., Dimarogonas, D.V., Tu, S., Matni, N., 2020.  Learning control barrier functions from expert demonstrations, in: 2020 59th IEEE Conference on Decision and Control (CDC), IEEE. pp. 3717–3724.
  * Safari and Hoagg [2023] Safari, A., Hoagg, J.B., 2023\.  Time-varying soft-maximum control barrier functions for safety in an a priori unknown environment.  arXiv preprint arXiv:2310.05261 .
  * Samanipour and Poonawala [2023a] Samanipour, P., Poonawala, H.A., 2023a.  Automated stability analysis of piecewise affine dynamics using vertices, in: 2023 59th Annual Allerton Conference on Communication, Control, and Computing (Allerton), pp. 1–8.  doi:10.1109/Allerton58177.2023.10313502.
  * Samanipour and Poonawala [2023b] Samanipour, P., Poonawala, H.A., 2023b.  Stability analysis and controller synthesis using single-hidden-layer relu neural networks.  IEEE Transactions on Automatic Control , 1–12doi:10.1109/TAC.2023.3270219.
  * Seel et al. [2021] Seel, K., Grøtli, E.I., Moe, S., Gravdahl, J.T., Pettersen, K.Y., 2021.  Neural network-based model predictive control with input-to-state stability, in: 2021 American Control Conference (ACC), IEEE. pp. 3556–3563.
  * Shamsoshoara et al. [2023] Shamsoshoara, A., Lotfi, F., Mousavi, S., Afghah, F., Guvenc, I., 2023.  Joint path planning and power allocation of a cellular-connected uav using apprenticeship learning via deep inverse reinforcement learning.  arXiv preprint arXiv:2306.10071 .
  * Srinivasan et al. [2020] Srinivasan, M., Dabholkar, A., Coogan, S., Vela, P.A., 2020\.  Synthesis of control barrier functions using a supervised machine learning approach, in: 2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), IEEE. pp. 7139–7145.
  * Sun et al. [2022] Sun, S., Zhang, Y., Luo, X., Vlantis, P., Pajic, M., Zavlanos, M.M., 2022.  Formal verification of stochastic systems with relu neural network controllers, in: 2022 International Conference on Robotics and Automation (ICRA), IEEE. pp. 6800–6806.
  * Taylor et al. [2020] Taylor, A., Singletary, A., Yue, Y., Ames, A., 2020.  Learning for safety-critical control with control barrier functions, in: Learning for Dynamics and Control, PMLR. pp. 708–717.
  * Wang et al. [2018] Wang, L., Theodorou, E.A., Egerstedt, M., 2018.  Safe learning of quadrotor dynamics using barrier certificates, in: 2018 IEEE International Conference on Robotics and Automation (ICRA), IEEE. pp. 2460–2465.
  * Wieland and Allgöwer [2007] Wieland, P., Allgöwer, F., 2007\.  Constructive safety using control barrier functions.  IFAC Proceedings Volumes 40, 462–467.
  * Zhang et al. [2023] Zhang, H., Wu, J., Vorobeychik, Y., Clark, A., 2023\.  Exact verification of relu neural control barrier functions.  arXiv preprint arXiv:2310.09360 .
  * Zhao et al. [2021] Zhao, Q., Chen, X., Zhang, Y., Sha, M., Yang, Z., Lin, W., Tang, E., Chen, Q., Li, X., 2021.  Synthesizing relu neural networks with two hidden layers as barrier certificates for hybrid systems, in: Proceedings of the 24th International Conference on Hybrid Systems: Computation and Control, pp. 1–11.