# Unconstrained Parametrization of Dissipative and Contracting Neural Ordinary
Differential Equations

Daniele Martinelli, Clara Lucía Galimberti, Ian R. Manchester, Luca Furieri,
and Giancarlo Ferrari-Trecate D. Martinelli, C. L. Galimberti, L. Furieri, and
G. Ferrari-Trecate are with the Institute of Mechanical Engineering, EPFL,
Switzerland. E-mail addresses: {daniele.martinelli, clara.galimberti,
luca.furieri<EMAIL_ADDRESS>Manchester is with the
Australian Centre for Robotics (ACFR) and the School of Aerospace, Mechanical
and Mechatronic Engineering, University of Sydney, Australia. E-mail address:
<EMAIL_ADDRESS>supported by the Swiss National Science
Foundation (SNSF) under the NCCR Automation (grant agreement 51NF40_80545).
Luca Furieri is also grateful to the SNSF for the Ambizione grant
PZ00P2_208951.

###### Abstract

In this work, we introduce and study a class of Deep Neural Networks (DNNs) in
continuous-time. The proposed architecture stems from the combination of
Neural Ordinary Differential Equations (Neural ODEs) with the model structure
of recently introduced Recurrent Equilibrium Networks (RENs). We show how to
endow our proposed NodeRENs with contractivity and dissipativity — crucial
properties for robust learning and control. Most importantly, as for RENs, we
derive parametrizations of contractive and dissipative NodeRENs which are
unconstrained, hence enabling their learning for a large number of parameters.
We validate the properties of NodeRENs, including the possibility of handling
irregularly sampled data, in a case study in nonlinear system identification.

##  I Introduction

Learning complex nonlinear mappings from data is a fundamental challenge in
various engineering applications, including computer vision, healthcare,
internet of things, and smart cities [1]. Deep Neural Networks (DNNs) have
emerged as powerful tools for this task, thanks to their adaptability and
ability to generalize their predictions from large amounts of data. However,
DNNs often lack robustness: a slight change in the input data may yield highly
different outputs, eventually leading to poor generalization capabilities [2,
3, 4, 5]. This deficiency becomes especially critical when DNNs are
implemented in real-world systems like safety-critical power grids or human-
interacting robots. In these applications, sensor measurements are inevitably
affected by multiple sources of noise and uncertainty, and a lack of
robustness may result in substantial losses.

To endow DNN models with formal stability and robustness guarantees, [3]
proposes to equip DNN layers with a dynamical system interpretation. The work
[6] establishes dynamical DNN models which universally approximate all
nonlinear dynamical systems defined in discrete-time and proposes a set of
convex constraints that enforce the stability of the DNN model during
training. In [7], the authors have developed discrete-time Recurrent
Equilibrium Networks (RENs) that result from the closed-loop interconnection
of a discrete-time linear dynamical system with a static nonlinearity. A main
contribution of [7] is to provide an unconstrained parametrization (also known
as direct parametrization) of a class of RENs with properties of stability and
dissipativity that are _built-in_ , i.e., that holds for any choice of the
parameters, and without the need to constrain them to a subset. This property
enables parameter optimization for very deep models through unconstrained-
gradient-descent-based algorithms, while ensuring stability and dissipativity
at any iteration. The usefulness of RENs for system identification and optimal
control has been demonstrated in [7] and [8], respectively. However, the
corresponding stability and dissipative guarantees are only compatible with
discrete-time or sampled-data systems. Studying how the properties of
nonlinear discrete-time dynamics port to their continuous-time counterpart is
usually a challenging problem [9].

The recently proposed Neural Ordinary Differential Equations (Neural ODEs)
[10] offer a bridge between DNNs and continuous-time dynamics. By unfolding
infinitely many DNN layers using a parametrized ODE and employing the adjoint
sensitivity method for training [10] , Neural ODEs offer several advantages
over discrete DNN models. These include memory efficiency, adaptive
computations beyond Euler-like integration methods, and the ability to handle
data arriving at arbitrary time instants.

Expanding on Neural ODEs, [11] introduced a continuous-time analogous of
Recurrent Neural Networks, albeit without stability or dissipativity
guarantees. The work [12] has proposed stable Hamiltonian Neural ODEs which
further exhibit non-vanishing gradient flows. Additionally, [5] developed a
class of Hamiltonian ODEs with contractivity properties. Furthermore, [13]
derived architectures ensuring convergence to a stable set. However, an
architecture combining the robustness and expressiveness of [6, 7] with the
benefits of Neural ODEs is not available to the best of our knowledge.

###  I-A Contributions

In this paper, we establish an alternative version of the REN architecture of
[7] that is compatible with Neural ODEs in continuous-time. Consequently, we
call our architectures NodeRENs. After showing that NodeRENs induce well-posed
dynamical systems in continuous-time, we establish that our models enjoy
stability and dissipativity properties by design, i.e., without any
constraints on the space of the parameters that describe the DNNs — akin to
their discrete-time counterpart [7]. NodeRENs, also inherit the advantages of
Neural ODEs in continuous-time; they are compatible with any integration
scheme that can be chosen based on a trade-off between accuracy and
computational resources, and they can be evaluated at any chosen time instant,
without the requirement to be uniformly sampled in time.

The paper is structured as follows. In Section II, we describe the problem
setting, including the definition of contractivity and Integral Quadratic
Constraints (IQCs). In Section III, we introduce the NodeREN model and detail
the steps to obtain NodeRENs that are contracting and dissipative by design.
Moreover, in Section IV the properties and performance of NodeRENs are
validated on a nonlinear system identification problem. Finally, in Section V
we summarize the conclusions of this work, outlining potential future
directions for research.

###  I-B Notation

We denote the set of non-negative real numbers as $+_{0}$. For $T>0$, let
$PC([0,T],n)$ be the space of piecewise-continuous functions in the time
interval $[0,T]$. We represent the Euclidean norm of $v\in n$ with $|v|$. For
a square matrix $X$, we use the notation $X\succ 0$ $(X\succeq 0)$ and $X\prec
0$ $(X\preceq 0)$ to denote positive (semi-) definiteness and negative (semi-)
definiteness, respectively. The minimum and maximum eigenvalues of the square
matrix $X$ are denoted as $\lambda_{min}(X)$ and $\lambda_{max}(X)$,
respectively. We use $(*)$ to represent a symmetric term in a quadratic
expression, e.g., $(X+Y)Q(*)^{\top}=(X+Y)Q(X+Y)^{\top}$, for some $Q\in
n\times n$ and $X,Y\in q\times n$. With $(*)$, we also indicate elements in
symmetric matrices that can be obtained by symmetry. We use $[{M}]_{p\times
m}$ to indicate the block matrix with the first $p$ rows and $m$ columns of
${M}\in n\times n$ with $n\geq p$ and $n\geq m$.

##  II Problem Formulation

Consider a nonlinear system $\Sigma_{\theta}$ in continuous-time

$\Sigma_{\theta}\>=\>\begin{cases}\>\dot{x}(t)=f_{\theta}(x(t),u(t))\\\
\>y(t)=g_{\theta}(x(t),u(t))\end{cases}\,,$ (1)

where $x(t)\in n$, $y(t)\in p$ and $u(t)\in m$ denote the state, output and
input of the system at any time $t\in+_{0}$, respectively, and $x(0)=x_{0}$.
In (1), $\theta\in n_{\theta}$ denotes a vector of parameters affecting the
behavior of the system. The idea behind Neural ODEs [10] is to interpret some
specific DNNs (e.g., ResNets [14]) as discretized ODEs. Then, the family of
Neural ODEs can be obtained when the discretization step converges to zero,
effectively passing from discrete-time DNNs to continuous-time models as
(1).111For instance, standard ResNet models whose hidden states evolve as
$h_{k+1}=h_{k}+\Delta t\,f(h_{k},\theta_{k})$ can be interpreted as the ODE
$\dot{h}(t)=f(h(t),\theta(t))$ as the discretization step $\Delta t$ tends to
zero. Accordingly, in this paper, we consider the following learning problem:

$\displaystyle\min_{\theta}\quad{L}(\theta,{\mathcal{Z}})$ (2)
$\displaystyle\textrm{subject to}\;\;\eqref{eq:diff_system}\,.$

$L$ is a scalar loss function that can depend on the trajectories ${x}(t)$ and
${y}(t)$ of (1), for $t\in[0,T]$ with $T>0$, and ${\mathcal{Z}}$ is a given
training dataset. For instance, in a classification task, $y(T)$ can play the
role of the DNN output layer which is compared to the label $y^{j}_{label}$ of
the point $x^{j}(0)$ contained in a dataset
${\mathcal{Z}}=\left\\{(x^{j}(0),y^{j}_{label})\right\\}_{j=1}^{N}$ for
$N\in\mathbb{N}$. Neural ODEs can exploit any chosen numerical method to
perform the forward propagation and the gradient computations through (1),
including those based on adaptive sampling to guarantee a desired level of
precision (e.g., Dormand-Prince, Bogacki–Shampine [15]). Although possible,
back-propagating through numerical solver’s operations can be highly expensive
in terms of memory and introduce numerical errors. As an alternative, to solve
(2) one can use the adjoint method — we refer to [10] for full details.

Neural-ODEs in their general form are not guaranteed to yield stable or
dissipative dynamical flows. Such properties are fundamental in optimal
control and system identification, as well as in robust learning problems
dealing with noisy features [5] and adversarial attacks [16]. Specifically,
motivated by [7], in this paper, we focus on continuous-time systems that are
contracting and systems that satisfy incremental IQCs. We proceed with
formally defining both.

Firstly, given a system $\Sigma_{\theta}$ with initial condition $a\in n$ and
an input function $u^{[1]}\in PC([0,\infty],m)$, let $x_{a}^{[1]}$ and
$y_{a}^{[1]}$ denote the corresponding state and output trajectories,
respectively.

######  Definition 1.

A system $\Sigma_{\theta}$ in the form (1) is said to be contracting if for
any two initial conditions $a,b\in n$ and given the same input trajectory
${u}^{[1]}\in PC([0,\infty],m)$, the corresponding state trajectories
$x_{a}^{[1]}$ and $x_{b}^{[1]}$ satisfy:

$|x_{a}^{[1]}(t)-x_{b}^{[1]}(t)|\leq\kappa e^{-ct}|a-b|\,,$ (3)

for all $t\in+_{0}$ and for some $c>0$, $\kappa>0$.

The Definition 1 can be interpreted as follows: a contracting system ‘forgets’
the initial condition exponentially fast as time progresses. Hence, all
trajectories converge to each other, independently of the initial state.

Next, we define dissipative systems that satisfy incremental IQCs. Dissipative
systems cannot increase their internal energy despite external inputs (e.g.,
feedback control actions or disturbances). As such, they can be designed to
possess finite input-output gains, a crucial property in nonlinear control
theory [17], as well as robust learning [18]. In order to formally define
these properties, let $a,b\in n$ be two initial conditions and
$u^{[1]},u^{[2]}\in PC([0,\infty],m)$ be two input trajectories. Define the
relative displacements as

$\begin{gathered}\Delta y(t)=y_{a}^{[1]}(t)-y_{b}^{[2]}(t)\>,\>\Delta
u(t)=u^{[1]}(t)-u^{[2]}(t)\>,\\\ \Delta
x(t)=x_{a}^{[1]}(t)-x_{b}^{[2]}(t)\>.\end{gathered}$ (4)

Let

$s_{\Delta}(\Delta u(t),\Delta y(t))=\begin{bmatrix}\Delta y(t)\\\ \Delta
u(t)\end{bmatrix}^{\top}\begin{bmatrix}Q&S^{\top}\\\
S&R\end{bmatrix}\begin{bmatrix}\Delta y(t)\\\ \Delta u(t)\end{bmatrix}\,,$ (5)

be a quadratic function $s_{\Delta}:m\times p\rightarrow$ parametrized by
$Q\in p\times p$, $S\in m\times p$, $R\in m\times m$. We recall the definition
of dissipative systems that satisfy IQCs according to the supply rate (5).

######  Definition 2.

A system $\Sigma_{\theta}$ in form (1) satisfies the incremental IQC defined
by the matrices $(Q,S,R)$, with $Q\preceq 0$ and $R=R^{\top}$, if there exists
a function $\mathcal{S}:n\to\mathbb{R}^{+}_{0}$ such that, for any two initial
conditions $a,b\in n$ and any two possible input functions $u^{[1]},u^{[2]}\in
PC([0,\infty],m)$,

$\mathcal{S}(\Delta x(t_{1}))\leq\mathcal{S}(\Delta
x(t_{0}))+\int_{t_{0}}^{t_{1}}s_{\Delta}(\Delta u(t),\Delta y(t))dt\,,$ (6)

for every $t_{1}\geq t_{0}$, where $\Delta x$, $\Delta y$ and $\Delta u$ are
defined in (4) and the supply rate $s_{\Delta}$ is defined in (5).

Despite being restricted to quadratic supply functions according to (5), IQCs
can certify many incremental properties by appropriately selecting the values
of $(Q,S,R)$. It is worth noting that the assumptions in Definition 2

$Q\preceq 0\,,\quad R=R^{\top}\,,$ (7)

are fulfilled for several incremental properties of interest, see Table I.

Incremental Property | $Q$ | $R$ | $S$ | $s_{\Delta}(\Delta u,\Delta y)$
---|---|---|---|---
$L_{2}$\- gain bound $(\gamma\geq 0)$ | $-\frac{1}{\gamma}I$ | $\gamma I$ | $0$ | $\gamma|\Delta u|^{2}-\frac{1}{\gamma}|\Delta y|^{2}$
Passivity | $0$ | $0$ | $\frac{1}{2}I$ | $\Delta u^{\top}\Delta y$
Input Passivity $(\nu\geq 0)$ | $0$ | $-2\nu I$ | $I$ | $\Delta u^{\top}\Delta y-\nu|\Delta u|^{2}$
Output Passivity $(\varepsilon\geq 0)$ | $-2\varepsilon I$ | $0$ | $I$ | $\Delta u^{\top}\Delta y-\varepsilon|\Delta y|^{2}$
TABLE I: Choices of $(Q,S,R)$ to verify different incremental properties.

To summarize, our goal is to learn the parameters $\theta$ of a dynamical
system (1) — i.e., a Neural ODE — that optimize a given cost as per (2), with
the hard constraint that either (3) or (6) (or both) hold. This new
requirement (i.e., the hard constraint) must be satisfied for all the
parameters $\theta$ we optimize over. In other words, we consider a form of
fail-safe learning, in the sense that the property of the model to be
contracting or dissipative must be guaranteed during and after parameter
optimization.

######  Remark 1.

In solving (2), one must select a numerical integration scheme for simulating
the trajectories of (1) to be used in forward propagation and gradient
computations. However, the properties of the continuous-time model (1) may not
be preserved in the discretized one, in general. While our focus in this paper
is on the contractivity and dissipativity properties of the continuous-time
model (1), it is worth noting that integration methods preserving
contractivity are discussed in [19]. Adapting these methods to IQC properties,
or developing new ones, represents an interesting future research direction.

##  III Contractivity and Dissipativity of NodeRENs

In this section, we present and analyze novel DNN structures that arise from
appropriately combining Neural ODEs [10] with RENs [7]. The starting idea is
to define the functions $f_{\theta}$ and $g_{\theta}$ in (1) through the model
utilized for the layer equation in [7]. Specifically, we consider

$\displaystyle\begin{bmatrix}\dot{x}(t)\\\ v(t)\\\
y(t)\end{bmatrix}=\overbrace{{\begin{bmatrix}A&B_{1}&B_{2}\\\
C_{1}&D_{11}&D_{12}\\\
C_{2}&D_{21}&D_{22}\end{bmatrix}}}^{\tilde{A}}\begin{bmatrix}x(t)\\\ w(t)\\\
u(t)\end{bmatrix}+\overbrace{\begin{bmatrix}b_{x}\\\ b_{v}\\\
b_{y}\end{bmatrix}}^{\tilde{b}},$ (8) $\displaystyle w(t)=\sigma(v(t)),$ (9)

where $x(t)\in n$, $u(t)\in m$ and $y(t)\in p$ are respectively the state, the
input, and output at time $t$. The function $\sigma(\cdot)$ represents a
nonlinear map and it is applied entry-wise. The input and output of
$\sigma(\cdot)$ are $v(t),w(t)\in\mathbb{R}^{q}$, respectively. We denote
system (8)-(9) as a NodeREN, and it can be interpreted as an affine time-
invariant system in closed-loop with a static nonlinearity $\sigma(\cdot)$. In
NodeRENs, the set of trainable parameters $\theta\in\mathbb{R}^{n_{\theta}}$
consists of the set $\tilde{A}$ of matrices $(A$, $B_{1}$, $B_{2}$, $C_{1}$,
$C_{2}$, $D_{11}$, $D_{12}$, $D_{21}$ and $D_{22})$ and the set $\tilde{b}$ of
vectors $(b_{x},b_{v},b_{y})$ in (8). It is worth noting that the model
(8)-(9) is highly flexible, as it encompasses many important neural network
architectures. For more details, see [7].

The work [7] has established conditions that ensure contractivity and IQC
properties of the _discrete-time_ dynamics induced by $\theta$. However, the
same conditions on $\theta$ fail to ensure these properties for the
continuous-time solutions of (8)-(9), in general.

Towards establishing contractivity and IQC-based dissipativity in continuous-
time for the model (8)-(9), consider two different possible trajectories of
the system, starting from two initial conditions $a$, $b\in n$ and two input
functions $u^{[1]}$, $u^{[2]}\in PC([0,\infty],m)$. Then, define the
incremental form of the system (8)-(9)

$\displaystyle\begin{bmatrix}\Delta\dot{x}(t)\\\ \Delta v(t)\\\ \Delta
y(t)\end{bmatrix}=\begin{bmatrix}A&B_{1}&B_{2}\\\ C_{1}&D_{11}&D_{12}\\\
C_{2}&D_{21}&D_{22}\end{bmatrix}\begin{bmatrix}\Delta x(t)\\\ \Delta w(t)\\\
\Delta u(t)\end{bmatrix}\,,$ (10) $\displaystyle\Delta
w(t)=\sigma(v_{b}^{[2]}(t)+\Delta v(t))-\sigma(v_{b}^{[2]}(t))\,,$ (11)

where $\Delta x$, $\Delta y$ and $\Delta u$ are defined in (4). Moreover,
$\Delta v(t)=v_{a}^{[1]}(t)-v_{b}^{[2]}(t)$, where $v_{a}^{[1]}$ and
$v_{b}^{[2]}$ are the inputs of $\sigma(\cdot)$ for each trajectory. Next, we
introduce two technical assumptions.

######  Assumption 1.

The function $\sigma(\cdot)$ belongs to $PC([0,\infty],)$ and its slope is
restricted to the interval $[0,1]$, that is

$0\leq\frac{\sigma(y)-\sigma(x)}{y-x}\leq 1\>,\quad\forall
x,y\in\mathbb{R}\>,\quad x\neq y.$

It is important to notice that, under Assumption 1, $\Delta v(t)$ and $\Delta
w(t)$ verify the following inequality

$\Gamma(t)=\begin{bmatrix}\Delta v(t)\\\ \Delta
w(t)\end{bmatrix}^{\top}\begin{bmatrix}0&\Lambda\\\
\Lambda&-2\Lambda\end{bmatrix}\begin{bmatrix}*\end{bmatrix}\geq 0,\quad\forall
t\in\mathbb{R},$ (12)

for any diagonal matrix $\Lambda\succ 0$. Note that most of the popular
activation functions used in the literature, such as the logistic function,
$ReLU(\cdot)$ and $\tanh{(\cdot)}$, satisfy Assumption 1.

######  Assumption 2.

$D_{11}$ in (8) is strictly lower-triangular.

Assumption 2 enforces that each scalar entry of $\Delta v(t)$ only depends on
the ones above it through (10). Hence, it becomes possible to explicitly
compute $\Delta v(t)$. This simplifies the numerical calculation of the
solutions to (10)-(11), while guaranteeing that the model is very expressive
thanks to the recursive application of the nonlinearity $\sigma$ on successive
entries of $\Delta v(t)$. The case where $D_{11}$ is not lower-triangular is
left for future work. We refer the interested reader to [7] for a discussion
on implicit RENs in discrete-time.

In general, an ODE model may not admit a unique solution for a given initial
condition and input trajectory [20]. We now proceed to show that (8)-(9)
admits a unique solution for any choice of the parameters $\theta$.

######  Lemma 1.

Let Assumptions 1 and 2 hold. Then, the model (8)-(9) admits a unique solution
in time for any $u(t)\in PC([0,\infty],m)$, $x(0)\in n$ and for any choice of
the parameters $\theta\in\mathbb{R}^{n_{\theta}}$.

The proof, reported in Section -A for completeness, is based on deriving a
global Lipschitz constant by iterating through the nonlinearities defining the
vector $v(t)$.

###  III-A Characterizing Contracting and Dissipative NodeRENs

Here, we derive sufficient conditions over the parameters $\theta$ to
guarantee contractivity and dissipativity of NodeRENs (8)-(9) according to a
specified supply rate (5).

######  Theorem 1.

A NodeREN (8)-(9) is contracting according to (3), if there exists a matrix
$P\succ 0$ and a diagonal matrix $\Lambda\succ 0$ such that

$\displaystyle\begin{bmatrix}-A^{\top}P-PA&-C_{1}^{\top}\Lambda-PB_{1}\\\
*&W\end{bmatrix}\succ 0\,,$ (13) $\displaystyle W=2\Lambda-\Lambda
D_{11}-D_{11}^{\top}\Lambda.$ (14)

The proof is reported in Section -B. For the rest of the paper, we refer to
NodeRENs that comply with (13)-(14) as C-NodeRENs.

######  Remark 2.

One may want to impose a certain convergence rate for a C-NodeREN. For this
purpose, (13) can be modified as follows:

$\begin{bmatrix}-A^{\top}P-PA-2\tilde{\omega}P&-C_{1}^{\top}\Lambda-PB_{1}\\\
*&W\end{bmatrix}\succ 0\,,$

where the term $-2\tilde{\omega}P$ enforces that (3) holds with
$c=\tilde{\omega}$. This can be directly seen by following the proof of
Theorem 1 starting from (31).

Next, we characterize NodeRENs that comply with incremental IQCs.

######  Theorem 2.

A NodeREN (8)-(9) satisfies the incremental IQC described by the triple
$(Q,S,R)$ fulfilling (7), if there exists a matrix $P\succ 0$ and a diagonal
matrix $\Lambda\succ 0$ such that

$\scalebox{0.93}{\mbox{$\displaystyle\begin{bmatrix}-A^{\top}P-PA&-C_{1}^{\top}\Lambda-
PB_{1}&-PB_{2}+C_{2}^{\top}S^{\top}\\\ *&W&-\Lambda
D_{12}+D_{21}^{\top}S^{\top}\\\ *&*&\mathcal{R}\end{bmatrix}$}}\\\
\scalebox{0.93}{\mbox{$\displaystyle+\begin{bmatrix}C_{2}^{\top}\\\
D_{21}^{\top}\\\ D_{22}^{\top}\end{bmatrix}Q\begin{bmatrix}C_{2}^{\top}\\\
D_{21}^{\top}\\\ D_{22}^{\top}\end{bmatrix}^{\top}\succ 0\,,$}}$ (15)

with $W$ given by (14) and

$\mathcal{R}=R+SD_{22}+D_{22}^{\top}S^{\top}\,.$ (16)

The proof is reported in Section -C. For the rest of the paper, we refer to
NodeRENs that comply with (15)-(16) as IQC-NodeRENs.

######  Remark 3.

Note that all IQC-NodeRENs are also contracting for any choice of $(Q,S,R)$
such that (7) holds. Indeed, if there exists $P\succ 0$ and a diagonal matrix
$\Lambda\succ 0$ such that (15) holds, then (13) also holds because $Q\preceq
0$ and the left-hand-side of (13) is a principal minor of that of (15).

It should be noted that solving the minimization problem (2) using the above
parametrizations of C-NodeRENs or IQC-NodeRENs requires solving a semidefinite
program (SDP) multiple times during the optimization process. This can be
computationally intractable for DNNs with several parameters. In the following
section, we propose a parametrization of a rich class of C-NodeRENs and IQC-
NodeRENs that circumvents this problem and allows for unconstrained
optimization. It is worth noting that this approach is similar to [7], but the
techniques differ due to the continuous-time nature of the proposed NodeRENs.

###  III-B Direct Parametrization of NodeRENs

Our goal is to define new free parameters $\theta_{C}\in\mathbb{R}^{n_{C}}$
and $\theta_{IQC}\in\mathbb{R}^{n_{IQC}}$ that can be mapped onto the
parameters $\theta$ of the NodeREN model, given by equations (8)-(9), and such
that either (13) or (15) are verified. First, we focus on constructing
C-NodeRENs from the parameters $\theta_{C}\in\mathbb{R}^{n_{C}}$ given by

$\theta_{C}=\\{X,{B}_{2},C_{2},{D}_{12},D_{21},D_{22},\tilde{b},U,Y_{1},X_{P}\\}\,,$
(17)

where $X\in(n+q)\times(n+q)$, $B_{2}\in n\times m$, $C_{2}\in p\times n$,
$D_{12}\in q\times m$, $D_{21}\in p\times q$, $D_{22}\in p\times m$,
$\tilde{b}\in(n+q+p)$ $U\in n\times q$, $Y_{1}\in n\times n$ and $X_{P}\in
n\times n$.

######  Theorem 3.

For any $\theta_{C}\in n_{C}$ defined in (17), and $\epsilon,\epsilon_{P}>0$
the two following statements hold.

  1. 1.

There are matrices $(Y,W,Z,P)$ of appropriate dimensions such that

$\displaystyle\begin{bmatrix}-Y^{\top}-Y&-U-Z\\\
*&W\end{bmatrix}=X^{\top}X+\epsilon I\,,$ (18) $\displaystyle
P=X_{P}^{\top}X_{P}+\epsilon_{P}I\,.$ (19)

  2. 2.

There are matrices $\\{\tilde{A},\tilde{b}\\}$ defined in terms of
$\theta_{C}$ and the matrices $(Y,W,Z,P)$ defined in point $1$, such that the
corresponding NodeREN (8)-(9) is contracting.

Moreover, matrices $(Y,W,Z,P)$ and $\\{\tilde{A},\tilde{b}\\}$ can be computed
as described in the proof of the Theorem.

The proof is reported in -D. In conclusion, an unconstrained parametrization
of a class of C-NodeRENs is obtained by first freely choosing $\theta_{C}$
defined in (17), and then recovering the matrix $\tilde{A}$ as per (34)-(37)
in Section -D.

Next, we focus on NodeRENs that satisfy IQC constraints by design, i.e. for
any choice of parameters $\theta_{IQC}\in n_{IQC}$ given by

$\theta_{IQC}=\\{X_{R},B_{2},C_{2},D_{21},\tilde{b},X_{3},T,U,Y_{1},X_{P}\\}\,,$
(20)

where $X_{R}\in(n+q)\times(n+q)$, $B_{2}\in n\times m$, $C_{2}\in p\times n$,
$D_{21}\in p\times q$, $\tilde{b}\in(n+p+q)$, $X_{3}\in s\times s$, $T\in
q\times m$, $U\in n\times q$, $Y_{1}\in n\times n$ and $X_{P}\in n\times n$.
In the next theorem, we provide a procedure to construct an IQC-NodeREN for
any choice of $\theta_{IQC}$.

######  Theorem 4.

Let $(Q,S,R)$ be such that (7) holds. Assume also that there exists $\delta>0$
satisfying $R-S(Q-\delta I)^{-1}S^{\top}\succ 0$. For any $\theta_{IQC}\in
n_{IQC}$ defined in (20), and $\epsilon,\epsilon_{P}>0$ the two following
statements hold.

  1. 1.

There are matrices $(Y,W,Z,P,D_{22})$ of appropriate dimensions such that

$\displaystyle\tilde{\mathcal{R}}=R+SD_{22}+D_{22}^{\top}S^{\top}+D_{22}^{\top}QD_{22}\succ
0\,,$ (21) $\displaystyle\begin{bmatrix}-Y^{\top}-Y&-U-Z\\\
*&W\end{bmatrix}-\Psi=X_{R}^{\top}X_{R}+\epsilon I\,,$ (22)

where $P$ is constructed as in (19), and

$\Psi=\begin{bmatrix}V\\\
\tilde{T}\end{bmatrix}\tilde{\mathcal{R}}^{-1}\begin{bmatrix}*\end{bmatrix}^{\top}-\begin{bmatrix}C_{2}^{\top}\\\
D_{21}^{\top}\end{bmatrix}Q\begin{bmatrix}*\end{bmatrix}^{\top}\,,$ (23)
$\displaystyle\tilde{T}=-T+D_{21}^{\top}S^{\top}+D_{21}^{\top}QD_{22}\>,$ (24)
$\displaystyle V=-PB_{2}+C_{2}^{\top}S^{\top}+C_{2}^{\top}QD_{22}\,.$ (25)

  2. 2.

There are matrices $\\{\tilde{A},\tilde{b}\\}$ defined in terms of
$\theta_{IQC}$ and the matrices $(Y,Z,W,P,D_{22})$ defined in point 1, such
that the corresponding NodeREN (8)-(9) satisfies the incremental IQCs
parametrized by $(Q,S,R)$.

Moreover, matrices $(Y,W,Z,P)$ and $\\{\tilde{A},\tilde{b}\\}$ can be computed
as described in the proof of the Theorem.

The proof of Theorem 4 shares similarities with the one of Proposition 2 in
[7] and it is reported in Section -E. Note that the assumption that there
exists a $\delta$ such that $R-S(Q-\delta I)^{-1}S^{\top}\succ 0$ is not
restrictive. Indeed, for all the most relevant cases, reported in Table I,
finding an appropriate value of $\delta$ is straightforward. In conclusion, an
unconstrained parametrization of IQC-NodeRENs is obtained by first freely
choosing the parameters $\theta_{IQC}\in n_{IQC}$, then constructing $D_{22}$
as per (42) in Appendix -E, and last recovering the matrices of $\tilde{A}$
using (34)-(37), (42) and (43).

##  IV Simulations & Results

In this section, we showcase the application of NodeRENs for identifying
stable nonlinear systems. We begin by assessing the performance of NodeRENs
when using various integration methods during training. Next, we highlight the
advantages of contractivity-by-design by training a general non-contractive
NodeREN (10)-(11), which fails to learn stable behavior for the same task.
Finally, we demonstrate the ability of NodeRENs to learn from irregularly
sampled time-series data while maintaining the guarantees by design. Our code
is accessible at https://github.com/DecodEPFL/NodeREN.

We remark that NodeREN trajectories are time-reversible, as established in
Lemma 1, which implies that they cannot intersect for distinct initial
conditions. To enhance model expressiveness, we employ augmented state
vectors, initializing additional scalar states to zero. This technique, known
as feature augmentation, was introduced in [21] for general Neural ODEs.

###  IV-A Continuous-Time System Identification

Here, we consider the system identification of a black-box system. For this
experiment, we assume that the unknown system dynamics are those of a
nonlinear pendulum

$\ell\ddot{\alpha}(t)+\beta\dot{\alpha}(t)+g\sin{\alpha(t)}=0\;,$ (26)

where $\alpha(t)$ is the angle position of the pendulum at time $t$ with
respect to the vertical axis, $\beta>0$ is the viscous damping coefficient,
$g$ is the gravitational acceleration and $\ell>0$ is the length of the
pendulum. For our experiments, we have chosen $\ell=$0.5\text{\,}\mathrm{m}$$
and $\beta=$1.5\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$$.

We consider the scenario where the system dynamics (26) to identify are
completely unknown and the only prior knowledge we have is that the system is
stable around the origin. Hence, we choose to train a C-NodeREN, which is
guaranteed to be contracting (and so, stable around the origin if $\tilde{b}$
is null) for any choice of the trainable parameters $\theta_{C}$ — i.e. even
if we stop the optimization prematurely.

The training data are given by noisy measurements of $\alpha(t)$ and
$\dot{\alpha}(t)$ across $N\in\mathbb{N}$ different experiments performed in
the time interval $[0,T_{end}]$. For each experiment $i=1,\ldots,N$, the
system (26) starts from a different and known initial condition
$\alpha_{i}(0)$ and $\dot{\alpha}_{i}(0)$. Due to the possibility of random
time delays during data acquisition, we assume that trajectory measurements
are taken at random time instants across $[0,T_{end}]$ for each experiment. It
is worth noting that inconsistent measurement times are common in control
applications; however, incorporating this feature into discrete-time systems
would pose a significant modeling challenge.

Precisely, the training data consist of
$\mathcal{Z}=\\{\mathcal{Z}_{i}\\}_{i=1}^{N}$, where

$\mathcal{Z}_{i}=\left\\{(t_{j},z_{j}),\text{ for all
}t_{j}\in\mathcal{T}_{i}\subset[0,T_{end}]\subset\mathbb{R}\right\\}\,,$ (27)

where $T_{end}=$3\text{\,}\mathrm{s}$$, $\mathcal{T}_{i}$ contains
$n_{i}\in\mathbb{N}$ time instants, and for each experiment $i=1,\ldots,N$,
the vector $z_{j}\in\mathbb{R}^{2}$ is measured as
$z_{j}=\begin{bmatrix}\alpha(t_{j})&\dot{\alpha}(t_{j})\end{bmatrix}^{\top}+w(t_{j})$
with $w(t_{j})\sim\mathcal{N}(0,0.01I_{2})$ being drawn according to a
Gaussian distribution. We minimize over $\theta_{C}$ the loss function

$L(\theta_{C},\mathcal{Z})=\frac{1}{N}\sum_{i=1}^{N}\frac{1}{n_{i}}\sum_{(t_{j},z_{j})\in\mathcal{Z}_{i}}|y(t_{j})-z_{j}|^{2}\,,$

where, with slight abuse of notation, $y(t_{j})\in\mathbb{R}^{2}$ denotes the
output at time $t_{j}$ of the C-NodeREN with parameters $\theta_{C}$.

The chosen C-NodeREN architecture has $n=4$ and $q=5$. It comes with
$\theta_{C}\in\mathbb{R}^{\sim 150}$, and it is trained using one of the
following integration methods: forward Euler (euler), Runge-Kutta of order 4
(rk4), and Dormand-Prince-Shampine of order 5 (dopri5). The training is
performed with $N=200$ experiments and using Adam [22].

To test the prediction performance of our trained models, we construct a test
dataset $\mathcal{Z}_{test}$ in the same way as (27), using a longer time
window $[0,T_{end}^{test}]$ with $T_{end}^{test}=$8\text{\,}\mathrm{s}$$ and
new sets of time instants $\mathcal{T}_{i}^{test}$. We then compute the test
loss $L(\theta_{C}^{\star},\mathcal{Z}_{test})$, where $\theta_{C}^{\star}$
are the parameters where the training converged to. In Figure 1, we report the
resulting test losses for different methods, along with their Number of
Function Evaluations (NFE) used for the prediction in the time window
$[0,T_{end}^{test}].$ Both euler and rk4 have been applied using
$\\{100,200\\}$ and $\\{50,200\\}$ integration steps, respectively, in
$[0,T_{end}^{test}]$. One of the advantages of variable-steps methods such as
dopri5 is the possibility to choose a tolerance value _tol_ to fix the
computational error during the ODE simulation. Indeed, it is possible to
obtain a trade-off between accuracy and the required number of function
evaluations just by properly setting this parameter. It is important to
emphasize how for fixed-step methods (e.g., euler, rk4), the chosen number of
steps may be insufficient to fully capture the dynamics of the trained model,
which can result in numerical instability. In contrast, adaptive methods
(e.g., dopri5) dynamically adjust the step size to ensure bounded integration
errors, hence promoting numerical stability in the trajectory simulations.
Such freedom enables us to either control the NFEs with a fixed step approach,
or to achieve more accurate simulations with an adaptive step size method. As
anticipated in Remark 1, it would be interesting to endow C-NodeRENs with
specific integration methods that preserve contractivity, e.g. [19].

Then, for the same identification task, we have trained a general NodeREN
(10)-(11) (that we denote as G-NodeREN), with trainable parameters
$\theta=\\{\tilde{A},\tilde{b}\\}$. Given the higher flexibility of G-NodeRENs
with respect to C-NodeRENs, one could expect a better performance in the time
window $[0,T_{end}]$ considered for the training. However, a trained G-NodeREN
can exhibit unstable dynamics when simulating for longer horizons; this
phenomenon is illustrated in Figure 2. Instead, C-NodeRENs only limit the
search to stable models, resulting in better identification performance for
the considered example. In Figure 2, we also display trajectory tubes from
different perturbed initial conditions. Note that the diameter of the tube
associated with the C-NodeREN decreases over time, showcasing its
contractivity.

Figure 1: Comparison between the number of function evaluations (NFE) and loss
on testing dataset for three integration methods.

###  IV-B Irregularly sampled-data

To show the robustness of C-NodeRENs concerning the irregular sample data, we
compare the test losses of $10$ C-NodeRENs trained on $10$ different training
datasets having the same initial conditions, but different sampling times. We
observe that, despite using differently sampled data, all the models have
similar test losses: they all lay in the interval $[3.7\times
10^{-4};9.1\times 10^{-3}]$, depending on how well the samples were
distributed.

Figure 2: Noisy test trajectory of the angular position $\alpha(t)$ of the
real system (in gray). Predicted trajectories of the trained C-NodeREN (in
orange) and G-NodeREN (in blue). Shaded areas represent tubes of additional
trajectories resulting from perturbations on the initial state of the systems.
The black dashed line indicates the end of the training horizon
$T_{end}=$3\text{\,}\mathrm{s}$$.

##  V Conclusions

In this work, we have established a class of Neural ODEs that generalize RENs
for continuous-time scenarios. Specifically, the proposed models guarantee
relevant continuous-time system-theoretic properties such as contractivity and
dissipativity by design. The resulting DNNs can be trained using unconstrained
optimization — which enables the use of large networks — and the resulting
architectures are flexible in the choice of an ODE integration scheme that is
suitable for the learning problem at hand. We have showcased the performance
of NodeRENs through the task of identifying the behavior of a nonlinear
pendulum, even with irregularly sampled data.

Future research should focus on integration schemes preserving contractivity
and IQC properties, exploring distributed NodeREN architectures, and analyzing
their generalization capabilities.

## References

  * [1] I. H. Sarker, “Machine learning: Algorithms, real-world applications and research directions,” SN computer science, vol. 2, no. 3, p. 160, 2021\.
  * [2] C. Szegedy, W. Zaremba, I. Sutskever, J. Bruna, D. Erhan, I. Goodfellow, and R. Fergus, “Intriguing properties of neural networks,” International Conference on Learning Representations, 2014.
  * [3] E. Haber and L. Ruthotto, “Stable architectures for deep neural networks,” Inverse Problems, vol. 34, p. 014004, dec 2017.
  * [4] M. Cheng, J. Yi, P.-Y. Chen, H. Zhang, and C.-J. Hsieh, “Seq2Sick: Evaluating the robustness of sequence-to-sequence models with adversarial examples,” Proceedings of the AAAI Conference on Artificial Intelligence, vol. 34, pp. 3601–3608, Apr. 2020.
  * [5] M. Zakwan, L. Xu, and G. Ferrari-Trecate, “Robust classification using contractive Hamiltonian neural ODEs,” IEEE Control Systems Letters, vol. 7, pp. 145–150, 2023.
  * [6] K.-K. K. Kim, E. R. Patrón, and R. D. Braatz, “Standard representation and unified stability analysis for dynamic artificial neural network models,” Neural Networks, vol. 98, pp. 251–262, 2018.
  * [7] M. Revay, R. Wang, and I. R. Manchester, “Recurrent equilibrium networks: Flexible dynamic models with guaranteed stability and robustness,” IEEE Transactions on Automatic Control, 2023.
  * [8] R. Wang and I. R. Manchester, “Youla-REN: Learning nonlinear feedback policies with robust stability guarantees,” 2022 American Control Conference (ACC), pp. 2116–2123, 2021.
  * [9] K. J. Åström, P. Hagander, and J. Sternby, “Zeros of sampled systems,” Automatica, vol. 20, no. 1, pp. 31–38, 1984.
  * [10] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud, “Neural ordinary differential equations,” Advances in Neural Information Processing Systems, 2018.
  * [11] P. Kidger, J. Morrill, J. Foster, and T. Lyons, “Neural controlled differential equations for irregular time series,” Advances in Neural Information Processing Systems, vol. 33, pp. 6696–6707, 2020.
  * [12] C. L. Galimberti, L. Furieri, L. Xu, and G. Ferrari-Trecate, “Hamiltonian deep neural networks guaranteeing non-vanishing gradients by design,” IEEE Transactions on Automatic Control, 2023.
  * [13] S. Massaroli, M. Poli, M. Bin, J. Park, A. Yamashita, and H. Asama, “Stable neural flows,” arXiv preprint arXiv:2003.08063, 2020.
  * [14] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 770–778, 2016.
  * [15] E. Hairer, G. Wanner, and S. P. Nørsett, Solving Ordinary Differential Equations I.  Springer Berlin Heidelberg, 1993.
  * [16] Q. Kang, Y. Song, Q. Ding, and W. P. Tay, “Stable neural ODE with Lyapunov-stable equilibrium points for defending against adversarial attacks,” Advances in Neural Information Processing Systems, vol. 34, pp. 14925–14937, 2021.
  * [17] A. Van der Schaft, L2-gain and passivity techniques in nonlinear control.  Springer, 2000.
  * [18] P. Pauli, A. Koch, J. Berberich, P. Kohler, and F. Allgöwer, “Training robust neural networks using Lipschitz bounds,” IEEE Control Systems Letters, vol. 6, pp. 121–126, 2021.
  * [19] I. R. Manchester, “Contracting nonlinear observers: Convex optimization and learning from data,” in 2018 Annual American Control Conference (ACC), pp. 1873–1880, 2018.
  * [20] J. Lygeros and F. Ramponi, “Lecture notes on linear system theory,” Automatic Control Laboratory, ETH Zurich, 2015.
  * [21] E. Dupont, A. Doucet, and Y. W. Teh, “Augmented neural ODEs,” Advances in neural information processing systems, vol. 32, 2019.
  * [22] S. Ruder, “An overview of gradient descent optimization algorithms,” arXiv preprint arXiv:1609.04747, 2017.
  * [23] H. K. Khalil, Nonlinear systems; 3rd ed. Upper Saddle River, NJ: Prentice-Hall, 2002.

###  -A Proof of Lemma 1

Given $u(t)\in PC([0,\infty],m)$, we can rewrite (8) as:

$\dot{x}(t)=p(x(t),t)=Ax(t)+B_{1}w(t)+B_{2}u(t)+b_{x},$ (28)

where we define $p(x(t),t)=f(x(t),u(t),t)$. It is well-known that a solution
exists and is unique if $p(\cdot,\cdot)$ is globally Lipschitz in its first
argument and piece-wise continuous in its second (e.g., [20, Theorem 3.6]).
Therefore, we prove that these properties hold for the model in (8) under
Assumptions 1 and 2. First, $p(\cdot,\cdot)$ is piece-wise continuous in its
second argument because it is the composition of piece-wise continuous
functions under Assumption 1. Then, we prove that $p(\cdot,\cdot)$ is globally
Lipschitz in its first variable, that is

$\exists\kappa>0:|p(x_{1},t)-p(x_{2},t)|\leq\kappa|x_{1}-x_{2}|\,,$

for every $x_{1},x_{2}\in n$ and $t\in\mathbb{R}$. For $j\in\\{1,2\\}$, denote
$w_{j}\in\mathbb{R}^{q}$ as
$w_{j}=\left[w_{j}^{1},\ldots,w_{j}^{q}\right]^{\top}=\sigma(v_{j})=\left[\sigma\left(v_{j}^{1}\right),\ldots,\sigma\left(v_{j}^{q}\right)\right]^{\top}$,
where $v_{j}=C_{1}x_{j}+D_{11}w_{j}+D_{12}\hat{u}$ for any possible value
$\hat{u}$ of the function $u(t)$ in time. By using the triangular inequality,
we obtain

$\displaystyle|p(x_{1},t)-p(x_{2},$ $\displaystyle
t)|=|A(x_{1}-x_{2})+B_{1}(w_{1}-w_{2})|$
$\displaystyle\leq|A|\,|x_{1}-x_{2}|+|B_{1}|\,|w_{1}-w_{2}|\,.$ (29)

Using Assumption 1 and Assumption 2 we verify that
$|w_{1}^{1}-w_{2}^{1}|\leq|C_{1}^{1:}||x_{1}-x_{2}|$, where $C_{1}^{i:}$
denotes the $i^{th}$ row of $C_{1}$, and we define $\kappa_{1}=|C_{1}^{1:}|$.
Proceeding similarly for each entry of $w_{1}-w_{2}$ we derive that, for every
$r=2,3,\ldots,q$, it holds that
$|w_{1}^{r}-w_{2}^{r}|\leq\kappa_{r}|x_{1}-x_{2}|$, with Lipschitz constants
$\kappa_{r}=|C_{1}^{r:}|+\sum_{\ell=1}^{r-1}|D_{11}^{r,\ell}|\kappa_{\ell}$,
where $D_{11}^{{i},j}$ denotes the element in the ${i}^{th}$ row and $j^{th}$
column of $D_{11}$. Finally, from (29) we conclude
$|p(x_{1},t)-p(x_{2},t)|\leq\kappa_{tot}|x_{1}-x_{2}|$, where
$\kappa_{tot}=|A|+\sum_{i=1}^{q}|B_{1}^{i:}|\kappa_{i}$. Thus, the
differential equation $\dot{x}(t)=f(x{(t)},u(t),t)=p(x{(t)},t)$ is globally
Lipschitz in $x$ with constant $\kappa_{tot}\in+_{0}$ and its solution exists
and it is unique.

###  -B Proof of Theorem 1

Define the function

$V(\Delta x(t))=\Delta x(t)^{\top}P\Delta x(t)\,.$ (30)

Let

$\Phi=(-A^{\top}P-PA)-(C_{1}^{\top}\Lambda+PB_{1})W^{-1}(*)^{\top}\succ 0,$

be the Schur complement of (13). Then, using properties of matrix
inequalities, it is always possible to find an
$\omega\in(0,{\lambda_{min}(\Phi)}\>{\lambda_{max}(P)}^{-1}]$ such that
$(\Phi-\omega P)\succ 0.$ Hence,

$\begin{bmatrix}-A^{\top}P-PA-\omega P&-C_{1}^{\top}\Lambda-PB_{1}\\\
*&W\end{bmatrix}\succ 0.$ (31)

By left- and right-multiplying the inequality (31) with $\begin{bmatrix}\Delta
x(t)^{\top}&\Delta w(t)^{\top}\end{bmatrix}$ and $\begin{bmatrix}\Delta
x(t)^{\top}&\Delta w(t)^{\top}\end{bmatrix}^{\top}$, respectively, and using
(12), we obtain by substitution that

$\dot{V}(\Delta x(t))\>+\omega V(\Delta x(t))<-\Gamma(t)\leq 0\;.$

Thus, according to Lyapunov’s exponential stability theorem [23], the
incremental system is globally exponentially stable, that is, it satisfies (3)
with $c=\nicefrac{{\omega}}{{2}}$.

###  -C Proof of Theorem 2

From Definition 2, if the function $\mathcal{S}(\Delta x(t))$ is continuous
and differentiable, (6) can be rewritten as:

$\frac{d}{dt}\Bigl{(}\mathcal{S}(\Delta
x(t))\Bigl{)}\>\leq\>s_{\Delta}\bigl{(}\Delta u(t),\Delta y(t)\bigl{)}\;.$
(32)

Consider as $\mathcal{S}(\Delta x(t))$ the function $V(\Delta x(t))$ as in
(30), which is continuous and differentiable in $\Delta x(t)$. By left- and
right-multiplying the inequality (15) with $\begin{bmatrix}\Delta
x(t)^{\top}&\Delta w(t)^{\top}&\Delta u(t)^{\top}\end{bmatrix}$ and
$\begin{bmatrix}\Delta x(t)^{\top}&\Delta w(t)^{\top}&\Delta
u(t)^{\top}\end{bmatrix}^{\top}$ respectively, we obtain by substitution

$\dot{V}(\Delta x(t))-\begin{bmatrix}\Delta y(t)\\\ \Delta
u(t)\end{bmatrix}^{\top}\begin{bmatrix}Q&S^{\top}\\\
S&R\end{bmatrix}\begin{bmatrix}*\end{bmatrix}<-\Gamma(t)\leq 0\,.$

Thus, the inequality (32) holds.

###  -D Proof of Theorem 3

_Part 1)_ : Define

$H=\begin{bmatrix}H_{11}&H_{12}\\\
H_{21}&H_{22}\end{bmatrix}=X^{\top}X+\epsilon I\,,$ (33)

where $H_{11}\in n\times n$, $H_{12}\in n\times q$, $H_{21}\in q\times n$ and
$H_{22}\in q\times q$. Then, one can set

$Y=-\frac{1}{2}(H_{11}+Y_{1}-Y_{1}^{\top})\,,$ (34) $W=H_{22}\,,\quad
Z=-H_{12}-U\,.$ (35)

Thus, starting from $H$ as in (33), $(Y,W,Z)$ can be built using (34)-(35) and
$P$ using (19).

_Part 2)_ : First, note that $H\succ 0$ and $P\succ 0$ are guaranteed by
construction because $X^{\top}X+\epsilon I$ and
$X_{P}^{\top}X_{P}+\epsilon_{P}I$ are positive definite for any $X$ and
$X_{P}$. From Assumption 2, ${D}_{11}$ is strictly lower triangular. Thus,
using (14) and being $W=W^{\top}$ by construction, it is possible to split $W$
as: $W=W_{diag}+W_{low}+W_{low}^{\top}$, where $W_{diag}$ is a diagonal
matrix, and $W_{low}$ is a strictly lower triangular matrix. By setting

$\Lambda=\frac{1}{2}W_{diag}\,,\quad D_{11}=-\Lambda^{-1}W_{low}\,,$ (36)

and the matrices $A,B_{1},C_{1}$ can be retrieved as

$C_{1}=\Lambda^{-1}U^{\top}\,,\quad A=P^{-1}Y\,,\quad B_{1}=P^{-1}Z\,,$ (37)

where $P^{-1}$ and $\Lambda^{-1}$ exist since $P\succ 0$ and $\Lambda\succ 0$.
Finally, we can show that (13) holds by substituting $A$, $B_{1}$, $C_{1}$,
$D_{11}$, $\Lambda$, and $P$ with their definitions in (33)-(37). Thus, the
corresponding NodeREN (8)-(9) specified by $\\{\tilde{A},\tilde{b}\\}$ is
contracting.

###  -E Proof of Theorem 4

We first need a technical lemma.

######  Lemma 2.

Let $p,m\in\mathbb{N}$, and $s=\max(p,m)$. Let $M\in s\times s$, and
$M=M^{\top}\succ 0$. Define

$F=(I-M)(I+M)^{-1}\,,$ (38)

and $\tilde{F}=[F]_{p\times m}\,.$ Then it follows that
$I-\tilde{F}^{\top}\tilde{F}\succ 0$.

###### Proof.

Notice that $F$ as defined in (38) is the Cayley transformation of $M$. First,
we list several properties of the Cayley transformation.

  1. P1

$M=M^{\top}\Rightarrow F=F^{\top}$,

  2. P2

$I+F=2(I+M)^{-1}$,

  3. P3

$I-F=2(I+M)^{-1}M$.

Moreover, note that the following properties hold for any symmetric matrices
$M,N$:

  1. P4

If $M\succ 0$, $N\succ 0$ and $MN=NM$, then $MN\succ 0$,

  2. P5

$(I+M)^{-1}M=M(I+M)^{-1}$.

We are now ready to prove Lemma 2 based on these properties.

Let $M$ be such that $M=M^{\top}\succ 0$. First, we prove that

$I-F^{\top}F\succ 0\,.$ (39)

This is equivalent to proving

$(I+F)^{\top}(I-F)\succ 0\,,$ (40)

since $F=F^{\top}$. Using P4, we prove (40) in three steps.

  1. 1.

From property P2, $(I+F)\succ 0$ since $2(I+M)^{-1}\succ 0$ thanks to $M\succ
0$.

  2. 2.

From property P3, $(I-F)\succ 0$ since $2(I+M)^{-1}M\succ 0$ from properties
P4 and P5.

  3. 3.

We prove that $(I+F)^{\top}$ and $(I-F)$ are commutative:

$\displaystyle(I+F)^{\top}(I-F)$ $\displaystyle=I-F+F^{\top}+F^{\top}F$
$\displaystyle=I+F^{\top}-F+FF^{\top}$ $\displaystyle=(I-F)(I+F)^{\top}$

where the second equality is valid from P1.

Then, (40) holds and this proves that $I-F^{\top}F\succ 0$.

Last, we prove that (39) implies $I-\tilde{F}^{\top}\tilde{F}\succ 0$. Note
that if $p=m$ then $\tilde{F}=F$ and this fact is trivial. We consider the two
remaining cases $p>m$ and $m>p$ separately.

Assume that $p>m$. Then, $F=\begin{bmatrix}\tilde{F}\\\ F_{b}\end{bmatrix}$.
Hence, $I-F^{\top}F=I-\tilde{F}^{\top}\tilde{F}-F_{b}^{\top}F_{b}\succ 0$,
which implies that $I-\tilde{F}^{\top}\tilde{F}\succ F_{b}^{\top}F_{b}\succeq
0$.

Assume that $m>p$. Then, $F=\begin{bmatrix}\tilde{F}&F_{r}\end{bmatrix}$.
Hence,

$I-F^{\top}F=\begin{bmatrix}I-\tilde{F}^{\top}\tilde{F}&-\tilde{F}^{\top}F_{r}\\\
-F_{r}^{\top}\tilde{F}&I-F_{r}^{\top}F_{r}\end{bmatrix}\succ 0\,,$

and by applying the Schur complement, the above equation implies
$I-\tilde{F}^{\top}\tilde{F}\succ 0$. ∎

Note that the property (38) has also been used in [7]. In this work, we have
derived a full proof for completeness. We are now ready to present the proof
of Theorem 4.

###### Proof.

_Part 1)_ : Define $\mathcal{Q}=Q-\delta I$. Since $\mathcal{Q}\prec 0$ and
$R-S(Q-\delta I)^{-1}S^{\top}\succ 0$, there exist $L_{Q}\succ 0$ and
$L_{R}\succ 0$ with

$L_{Q}^{\top}L_{Q}=-\mathcal{Q}\>,\quad
L_{R}^{\top}L_{R}=R-S\mathcal{Q}^{-1}S^{\top}\,.$ (41)

Next, define $M=X_{3}^{\top}X_{3}+\epsilon I$ and
$\tilde{F}=\left[(I-M)(I+M)^{-1}\right]_{p\times m}$. Construct $D_{22}$ as

$D_{22}=-\mathcal{Q}^{-1}S^{\top}+L_{Q}^{-1}\tilde{F}L_{R}.$ (42)

We now show that this choice of $D_{22}$ leads to $\tilde{\mathcal{R}}\succ
0$. Let $\Xi=R+SD_{22}+D_{22}^{\top}S^{\top}+D_{22}^{\top}\mathcal{Q}D_{22}$.
Since $\tilde{\mathcal{R}}=\Xi+\delta D_{22}^{\top}D_{22}$, it suffices to
prove that $\Xi\succ 0$. After plugging (42) in the definition of $\Xi$, we
obtain $\Xi=L_{R}^{\top}(I-\tilde{F}^{\top}\tilde{F})L_{R}$ by cancellation of
addends and by using (41). By Lemma 2 and since $L_{R}\succ 0$, we conclude
that $\Xi\succ 0$.

Finally, we indicate how to construct $(Y,W,Z)$ so that (22) holds. Compute
$\Psi$ from (23) and define

$H=\begin{bmatrix}H_{11}&H_{12}\\\
H_{21}&H_{22}\end{bmatrix}=X_{R}^{\top}X_{R}+\epsilon I+\Psi\,.$

Then, using $H$, it is possible to use the same steps of Theorem 3 to compute
$(Y,W,Z)$ from (34)-(35), and $P$ from (19). This choice satisfies the
equality (22) by construction.

_Part 2)_ : We show how to parametrize the matrices in $\tilde{A}$. First,
note that $P\succ 0$ is valid due to (19). Construct $\Lambda$ and $D_{11}$
using (36), and define the matrices $(A,B_{1},C_{1},D_{11})$ according to
(37). Set

$D_{12}=\Lambda^{-1}T\,.$ (43)

Then, when choosing the parameters according to part 1 of this Theorem and
plugging in the matrices just defined, $\Psi$ satisfies
$\Psi=\Psi^{\top}\succeq 0$ since $\tilde{\mathcal{R}}\succ 0$ and $Q\preceq
0$. Moreover, $H-\Psi=X_{R}^{\top}X_{R}+\epsilon I\succ 0$. Thus,
$H=H^{\top}\succ 0$, and we have $W=W^{\top}\succ 0$ and $\Lambda\succ 0$.
Thus, the NodeREN defined by $\\{\tilde{A},\tilde{b}\\}$ complies with (15),
and thus it satisfies the incremental IQCs parametrized by $(Q,S,R)$. ∎


# Verifying Safety of Behaviour Trees in Event-B

Matteo Tadiello Elena Troubitsyna KTH Royal Institute of Technology
Stockholm, Sweden<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Behavior Trees (BT) are becoming increasingly popular in the robotics
community. The BT tool is well suited for decision-making applications
allowing a robot to perform complex behavior while being explainable to humans
as well. Verifying that BTs used are well constructed with respect to safety
and reliability requirements is essential, especially for robots operating in
critical environments. In this work, we propose a formal specification of
Behavior Trees and a methodology to prove invariants of already used trees,
while keeping the complexity of the formalization of the tree simple for the
final user. Allowing the possibility to test the particular instance of the
behavior tree without the necessity to know the more abstract levels of the
formalization.

##  1 Introduction

Autonomous Systems (AS) like Humanoid Robots, Autonomous Vehicles, or Unmanned
Aerial Vehicles are becoming increasingly complex and need to interact with
dynamic environments and with each other. For this reason, robots require
tools to enable advanced perception and understanding of the environment, or
capabilities to operate in complex situations. Artificial Intelligence is
extending the capability of perception and action of the agents and allows
robots to operate in environments not suitable for robots just a few years
ago.

In most common scenarios the complexity of the environment requires to the
robot to have different skills, the capability of different actions, and hence
also a certain degree of reasoning and understanding of which action to take
and when. A relevant example could be an urban road, with car, pedestrian, and
signals. That’s why planning becomes more important with the increasing
complexity of the Operational Design Domain (ODD). Various techniques have
been developed to solve this challenge. They include a variety of machine
learning approaches as well as Behaviour Trees. Behaviour Trees started as a
tool for video-game development of Non-Playable Character [20], and are now
being applied also to robotics due to their flexibility and high usability.

Usually, robotics systems should be dependable, i.e., being able to operate
without putting at risk others, or themselves, while continuing to provide
their service in a reliable way. Formal approaches can be used to provide
these characteristics and ensure the safety and reliability of the systems.
Applying formal methods to AI techniques, however, is still challenging given
the enormous search space and states possibility, and became particularly
relevant when dealing with Neural Networks with billions of parameters. The
importance of safety in robots that need to operate in a day-to-day basis in a
critical environment is extremely important, as well as the capability of the
system to explain why certain decisions have been made. For this reason
Behaviour Trees (BT) have become particularly interesting in the robotics
industry. This mathematical model follows simple principles and allow
readability, modularity and re-usability. And for this reason, it has gain
popularity especially in decision making tasks.

However, even with these properties BT can become quite complex and easy to
contain design mistakes that could mine the safety and reliability of the
system. For these reasons, it is important to ensure that certain properties
are kept during the execution of the BT that ensure the dependability of the
system built. The usage of formal methods to ensure the correctness of BT is
still scarce, and difficult to build.

In this work, we try to advance in this direction by providing a formal
specification of the Behaviour Tree framework that will allow the verification
of correctness of specific BT instances and holding of interesting invariant
properties that the system should maintain. In particular, we are able to
define safety invariants and prove that these hold for specific BTs, while
maintaining the re-usability and modularity of the tree structure. The system
modeling has been built using the Event-B modeling method [3] and the Rodin
[2] platform as a support for the refinements and mathematical proofs.

The rest of this paper is structured in the following way: we first provide an
introductory explanation of Behaviour Trees and Event-B in the Background
section, we continue with the Method section presenting the BT formal
specification, followed by a case study section where we show a simple BT
instance and the method to convert a BT to a formal specification that can be
checked. We will then move to the related work and end with the conclusion
section.

##  2 Background

###  2.1 Behaviour Trees

At a high level, a Behaviour Tree is a mathematical model that structures the
switching of different tasks, which a virtual agent or a robot can perform
[10]. This model uses a rooted tree where every node has a single parent and
no loops. The nodes in a BT can be subdivided into three categories:

  * •

root: the root node. It does not have parents, and it normally has just one
child.

  * •

internal node: also known as control flow nodes, those are nodes that are at
the intermediary level of the tree. They all have a parent and at least a
child. They represent how to navigate the tree.

  * •

leaf node: also known as execution nodes, these nodes have a parent but do not
have a child. They represent the tasks that the agent needs to perform (like
”go to point A”) or conditions (like ”battery ¿ than 30%) that needs to be
true before executing an action.

In practice, the BT is navigated from the root to the leaves to decide which
task needs to be executed at every time. To navigate the tree a process called
_ticking_ is used. An example of a BT and its ticking process can be seen in
Fig 1. The navigation starts by ticking the root and then recursively ticking
the left-most child of the node. When a leaf node is ticked it immediately
replies with the status and sends the signal back to its parents, then the
continuation of the navigation is pursued by following the rules of the
internal nodes. A node is executed only if it has been ticked. At each time
step, the tree is explored and returns the execution status of the agent. The
output value can be $SUCCESS$, $RUNNING$, or $FAILURE$ which, respectively,
represent if the behavior of the agent, is in a success state, in a failure,
or if it is still running. We call this the status of the tree.

The status of the tree is due to the composition of the statuses of the nodes
in each subtree analyzed, and each node follows different rules to provide its
status.

More specifically, for the leaf nodes we can distinguish two kinds:

  * •

Conditions: these nodes represent guards check, they are normally placed
before tasks nodes, they can provide just $SUCCESS$ or $FAILURE$, but not
$RUNNING$)

  * •

Tasks or Actions: these nodes represent the actual actions executed by the
agent. Once ticked, they immediately return SUCCESS, RUNNING, or FAILURE.

The leaf nodes are the ones that check the status of the agent and execute
actions. The internal nodes instead define how leaf nodes are, or are not,
executed.

The internal nodes, analyzed, are divided into two categories:

  * •

Sequence: these nodes execute a series of children in sequence. The next child
is executed only if the previous child has returned SUCCESS. The sequence node
return to their parent SUCCESS if and only if all the children have returned
SUCCESS. It returns RUNNING if one of the children is RUNNING, and FAILURE
otherwise (hence one of the children has returned FAILURE). The sequence
symbol is a box with the ”$\rightarrow$” label.

  * •

Fallback: these nodes tick a series of children following a left-to-right
order, starting from the leftmost child. The next child is executed only if
the previous child has returned FAILURE. The fallback node return to their
parent FAILURE if and only if all the children have returned FAILURE. It
returns RUNNING if one of the children is RUNNING, and SUCCESS otherwise
(hence one of the children has returned SUCCESS). The fallback symbol is a box
with the ”$?$” label.

Figure 1: A sample of a BT and its ticking process. The execution of the tree
returns RUNNING while executing the action of opening the door.

BT can be extended also with other two internal nodes categories called
PARALLEL and DECORATORS. Those nodes, less frequent, are not treated in this
work.

The main strengths of Behaviour Tree are that they are simple and readable,
and each subtree is reusable given the modularity property of this structure.
But they can represent also complex behaviors for robots operating in
challenging environments. Moreover, BT can also be extended with different
properties like Stochastic BTs [6] or automatically generated as in [8].

###  2.2 Event-B and Rodin

The formal development framework used in this work is the _Event-B_ formalism.
In Event-B a system specification is defined using the concept of _abstract-
state machine_. This model encapsulates a collection of variables and contains
the set of operations to interact with the states of the machine called
Events. We can thus define the behavior of the machine by describing the
dynamic part of the state machine. Another component of the Event-B method is
the so-called context, these components contain the static part of the system.
For instance, in the context user-defined types are defined or constant
components of the system. The context can be seen by the abstract machines as
shown in 2.

Machine $\;M$ Variables $\;v$ Invariants $\;I$ Events $Initialisation$
$evt_{1}$ $\cdots$ $evt_{N}$ $\longrightarrow$ Context $\;C$ Carrier Sets
$\;d$ Constants $\;c$ Axioms $\;A$

Figure 2: Event-B machine and context

The dynamic part of the system is represented by the variables of the
abstract-state machines and the Cartesian product of their range represents
the state space of the machine. The variables are strongly typed by invariant
statements and are initialized using a special initialization event. The
invariant statements, not only define the type of the variables but most
importantly, are used to define properties that the system should hold for the
entire time of execution. For instance, a possible invariant could be a safety
property like: ”the distance to the obstacles needs to be at least 5 meters
from the robot”.

The changing of the variables, and hence the state of the machine, is
represented through a set of Event clauses defined in the machine. An event is
normally defined as follows:

$\mathbf{evt}\;\widehat{=}\;\mathbf{any}\;variables\;\mathbf{where}\;guards\;\mathbf{then}\;S\;\mathbf{end}$

where variables is a new list of local variables, guards are predicates
defining when the event can be triggered, and $S$ is a set of statements or
assignment that changes the state of the machine. However, both the variables
and the guards are optional in the definition of an event. Multiple events can
be enabled at the same time, in this case, any of them can be chosen to be
executed in a non-deterministic way.

One of the strong points of Event-B is that it allows the possibility of
machine refinement. This allows us to start the definition from a very high
and abstract level, and include more details, variables, and events at every
refinement. We can then introduce implementation details in steps while
preserving functional correctness.

The semantics of an Event-B model is formulated as a collection of proof
obligations [3], and the process of verification and generation and proving of
proof obligations, has been facilitated by the Rodin platform [2].

##  3 Specification of Behaviour Trees in Event-B

The formal specification of a generic Behaviour Tree has been implemented and
proven using Event-B and Rodin111The source code can be found here:
https://github.com/tadteo/BT_eventb. To do that, we have used 4 different
refinements. We built initially the abstract figure of the nodes, till
arriving at the instances of the BT, which can be found as a refinement of the
4th level. Each abstract machine is coupled with a corresponding context. The
four abstraction levels are described in detail in the following subsections,
and are:

  * •

Nodes

  * •

Tree

  * •

Behavior Tree

  * •

Specific Instance of a BT.

###  3.1 Specifying a node

We start defining the nodes and certain properties that the nodes composing
the Tree need to have. At this level, we are not interested in the connection
and relationship between the nodes but just in their properties.

First of all, we define two new types to define the kind of node and the
result that each node can return.

The $TYPE$ type represents the category that each node can be and it is
defined as follow:

$TYPE=\\{ROOT,SEQUENCE,FALLBACK,CONDITION,ACTION\\}$

It contains the types for the root, control flow nodes, and execution nodes,
this will be used in a future refinement to define the different navigation
behavior of the nodes.

The RESULT type, instead, represents the value of return that each node can
contain. We define it as follows:

$RESULT=\\{SUCCESS,RUNNING,FAILURE,UNKNOWN\\}$

The first three elements are the normal return value of a node in a BT, the
$UNKNOWN$ element instead is used as the default value when a node has not
been ticked yet, and hence the return value is unknown.

We wanted to define the single nodes as a data structure containing the
following members:

        node {
          int n_id; //constant value in context
          TYPE n_type; //constant value in context
          bool n_tick; //dynamic value in machine
          RESULT n_result; //dynamic value in machine
        };

Where n_id represents a unique identifier of the node that starts from 0 with
the root. This value is used to define an order on the nodes used in a future
refinement to identify the order of the children, used to navigate the tree
from left to right. And n_tick is a boolean variable to see if the node has
been ticked or not.

Since in Event-B there is not a standard definition of a structure, the
members of the nodes have been defined as relationships between the nodes and
the corresponding type as follows:

$\textnormal{{n_type}}\in NODES\rightarrow TYPE$ $\textnormal{{n_id}}\in
NODES\rightarrow\mathbb{N}$ $\textnormal{{n_tick}}\in NODES\rightarrow BOOL$
$\textnormal{{n_result}}\in NODES\rightarrow RESULT$

and initialized as:

$\textnormal{{n_tick}}\in NODES\times\\{FALSE\\}$ $\textnormal{{n_result}}\in
NODES\times\\{UNKNOWN\\}$

where $NODES$ is the set of all possible nodes that can be used to build the
tree, and that has been defined in the node context together with the types,
and the constant values. n_id and n_type does not require initialization since
they are defined in the context and remain constant.

###  3.2 Formal model of tree

In the first refinement, we give the nodes a structure and build the rules to
construct a tree. We used the following four requirements to define the tree
structure:

  * •

Req1: There is a single root object

  * •

Req2: Each object other than the root has a parent

  * •

Req3: There are no loops in the parent structure

  * •

Req4: Each object is reachable from the root

First of all, we define the nodes constant which represents the set of nodes
(a subset of NODES) in the tree, with root as the ROOT node.

For Req1 we define an axiom that states that for every pair of nodes, if one
is a ROOT then the other is not ROOT:

$\forall n1,n2\cdot(n1\in nodes\land n2\in nodes\land n1\neq n2\land
n\\_type(n1)=ROOT\implies n\\_type(n2)\neq ROOT)$ (1)

For Req2 we define the relation parent as following:

$parent\in nodes\backslash\\{root\\}\rightarrow nodes$ (2)

For the no loop property we use the inverse of the parent function (i.e. the
children)222In Event-B in Rodin, the inverse is written using the $\sim$
symbol.In this case, the parent function would be written as:
$parent\sim[\\{node\\}]$ in the following way:

$\forall n\cdot(n\subseteq parent^{-1}[n]\implies n=\emptyset)$ (3)

For Req4 we use a transitive closure on a node relation, that in a graph
represents a relation to any other node in the graph [15].

To achieve that, we first define a generic relation between two nodes, and
then we define the transitive closure of it.

$node\\_rel=NODES\leftrightarrow NODES$ $tcl\in node\\_rel\rightarrow
node\\_rel$ $\forall r\cdot(r\in node\\_rel\implies r\subseteq tcl(r)$
$\forall r\cdot(r\in node\\_rel\implies r;tcl(r)\subseteq tcl(r)$ // The
unfolding of tcl is also part of tcl $\forall r,t\cdot(r\in node\\_rel\land
r\subseteq t\land r;t\subseteq t\implies tcl(r)\subseteq t$ //tcl(r) is least
$\forall r\cdot(r\in node\\_rel\implies r\cup(r;tcl(r))$

###  3.3 Behavior Tree level

The Behavior Tree level starts from the definition of the tree present in the
level above and adds all the rules to allow the correct navigation and return
of the result.

In this level, we simulate the ticking process, which starts from a completely
unvisited and unticked tree. As the first step, the root needs to be ticked.
After this, a chain reaction of events brings to the exploration of the tree
which eventually returns a result. The exploration ends whenever the result
arrives at the root node. At this point, the tree got reinitialized and the
root get ticked again, restarting the cycle.

To achieve correct navigation of the tree, each type of node needs to be
treated separately with specific events.

####  3.3.1 Root related events

To deal with the root we have 3 distinct events:

The first tick_root is used to tick the root node. It is used whenever the
tick status of the root is FALSE and the result is $UNKNOWN$:

$\mathsf{tick\\_root}\;\widehat{=}$ $\textbf{any}\;\;node$ where
$n\\_tick(node)=FALSE$ $n\\_type(node)=ROOT$ $n\\_result(node)=UNKNOWN$
$node\in NODES$ $node\in nodes$ then $n\\_tick(node):=TRUE$

We then have an event that is used to tick the first unticked child of the
root as soon as the root is ticked. In our specification, the root has just
one child and is treated as a special node. This allows a simpler way to
manage the start and end of an exploration of the tree.

$\mathsf{root\\_ticked}\;\widehat{=}$ any $\;\;node$ $\;\;child$ where
$n\\_tick(node)=TRUE$ $n\\_type(node)=ROOT$ $n\\_tick(child)=FALSE$ $node\in
nodes$ $child\in\\{x|x\in parent^{-1}[\\{node\\}]\land
n\\_tick(x)=FALSE\land\\\ \forall y\cdot(y\in parent^{-1}[\\{node\\}]\land
n\\_tick(y)=FALSE\implies n\\_id(y)\geq n\\_id(x)\\}$ then
$n\\_tick(child):=TRUE$ $analyzing\\_subtree(child):=TRUE$

Figure 3: The event $\mathsf{root\\_ticked}$

To identify the first unticked child we make use of the id of the node. The id
of the nodes has been given in a way that respects the order of visiting of
the tree, following a left-to-right order, with siblings having the leftmost
with a smaller id and the rightmost with the highest id. When choosing which
node to visit, we select as a child the node that has the n_id smaller than
all the other nodes, and that is not yet ticked. Finally, when the ticking
process return to the root an event called result_arrived is used to copy the
result of the child of the root to the root.

$\mathsf{result\\_arrived}\;\widehat{=}$ any $\;\;node$ $\;\;child$ where
$node\in NODES$ $n\\_tick(node)=TRUE$ $n\\_type(node)=ROOT$
$n\\_result(node)=UNKNOWN$ $child\in parent^{-1}[\\{node\\}]$
$child\in\\{x|x\in parent^{-1}[\\{node\\}]\land n\\_result(x)\neq UNKNOWN\\}$
then $n\\_result(node):=n\\_result(child)$

Before restarting the ticking process an event called root_reinitialize is
used to reinitialize the tree nodes as unexplored and with unknown results.

$\mathsf{root\\_reinitialize}\;\widehat{=}$ any $\;\;node$ where
$n\\_tick(node)=TRUE$ $n\\_type(node)=ROOT$ $n\\_result(node)\neq UNKNOWN$
$node\in NODES$ $node\in nodes$ then $n\\_result:=NODES\times\\{UNKNOWN\\}$
$n\\_tick:=NODES\times\\{FALSE\\}$

Figure 4: The event $\mathsf{root\\_reinitialize}$

####  3.3.2 Leaf nodes related events

The leaf nodes can be of two types: $CONDITION$ or $ACTION$. The leaf nodes
are the only ones that directly provide a result with the ticking of the
nodes. However, these events need to be refined at the next level for each
case of result ($SUCCESS$, $FAILURE$ for the $CONDITION$ nodes and $SUCCESS$,
$FAILURE$ and $RUNNING$ for the $ACTION$ nodes). And providing the correct
logic to return the result wanted, since for each BT the logic to return the
result can be different.

At this level of abstraction, the only action executed by these nodes is
setting the analyzing_subtree variable of the parent node to FALSE. This is
done as follows:

$analyzing\\_subtree(parent(node)):=FALSE$ (4)

####  3.3.3 Internal nodes related events

The internal nodes are the ones in charge of controlling the flow of execution
and analysis of the tree. In this work, we have analyzed the sequence node and
the fallback node. To formalize the behavior of this kind of node, we have
used 5 different events for each kind covering all the possible ways in which
the node can be found while analyzing the tree. At each step of the tree
exploration, the control flow nodes either have enough information to return a
value to the parents or need to continue the analysis of the children and tick
the next child on the list. The five events have the following form:  1.
*_ticked_initial,  2. *_ticked_success,  3. *_ticked_failure,  4.
*_ticked_running,  5. *_ticked_continue,  with * used as a wildcard that can
be substituted with the name of the node (like fallback_ticked_continue).

The *_ticked_initial event is used to start the ticking process of the
children and it’s only purpose is to tick the first children of these nodes.

Let’s analyze now the events of the fallback node.
The definition is ”The Fallback node […] corresponds to routing the ticks to
its children from the left until it finds a child that returns either SUCCESS
or RUNNING, then it returns SUCCESS or RUNNING accordingly to its own parent.
It returns FAILURE if and only if all its children return FAILURE. Note that
when a child returns RUNNING or SUCCESS, the Fallback node does not route the
ticks to the next child (if any)” [10].
To be able to route the ticks to the correct child (from left to right), we
first need to correctly identify the first unticked child and last analyzed
child. To do that we use a similar approach to the one used in the root_ticked
event. We use the inverse of parent to identify the children of a node and
then we use the n_id of the children to analyze the order, finally, we use the
n_result and n_tick variables to check the correct value we are looking for.

For instance, to identify if the last child analyzed has returned the
$RUNNING$ result (used in fallback_ticked_running, Figure 5, we use the
following guard:

$\exists x\cdot(x\in parent^{-1}[\\{node\\}]\land n\\_tick(x)=TRUE\land\\\
\forall y\cdot(y\in parent^{-1}[\\{node\\}]\land n\\_tick(y)=TRUE\implies
n\\_id(y)\leq n\\_id(x))\\\ \land n\\_result(x)=RUNNING)$ (5)

In this code, we check if exists a node of the children set, which has been
ticked, that has the greatest id of all the other ticked children, and that
has the result set to RUNNING. To find the first unticked children instead we
use the same approach but instead check that the n_tick variable is set to
false and that the n_id value is the minimum of the unticked children. This
can be seen in the events fallback_ticked_continue Fig. 5, and
sequence_ticked_continue Fig. 6,

In some events, we want also to be able to identify whenever all the children
have been analyzed (ticked). In these cases, we just check that the set of
children of a node that has n_tick set to FALSE is an empty set. The guard
used to check this is:

$\\{x|\;x\in parent^{-1}[\\{node\\}]\;\land\;n\\_tick(x)=FALSE\\}=\emptyset$
(6)

With these guards, we are now able to define the 4 events that allow us to
correctly manage the Fallback node (Figure 5).

Events $\mathsf{fallback\\_ticked\\_running}\;\widehat{=}$
$\textbf{any}\;\;node$ where $n\\_type(node)=FALLBACK$ $n\\_tick(node)=TRUE$
$n\\_result(node)=UNKNOWN$ $\exists x\cdot(x\in parent^{-1}[\\{node\\}]\land
n\\_tick(x)=TRUE\land\\\ \forall y\cdot(y\in parent^{-1}[\\{node\\}]\land
n\\_tick(y)=TRUE\implies\\\ n\\_id(y)\leq n\\_id(x))\land
n\\_result(x)=RUNNING)$ then $n\\_result(node):=RUNNING$
$analyzing\\_subtree(parent(node)):=FALSE$ end
$\mathsf{fallback\\_ticked\\_failure}\;\widehat{=}$ $\textbf{any}\;\;node$
where $n\\_type(node)=FALLBACK$ $n\\_tick(node)=TRUE$
$n\\_result(node)=UNKNOWN$ $\exists x\cdot(x\in parent^{-1}[\\{node\\}]\land
n\\_tick(x)=TRUE\land\\\ \forall y\cdot(y\in parent^{-1}[\\{node\\}]\land
n\\_tick(y)=TRUE\implies\\\ n\\_id(y)\leq n\\_id(x))\land
n\\_result(x)=FAILURE)$ $\\{x|\;x\in
parent^{-1}[\\{node\\}]\;\land\;n\\_tick(x)=FALSE\\}=\emptyset$ then
$n\\_result(node):=FAILURE$ $analyzing\\_subtree(parent(node)):=FALSE$ end
$\mathsf{fallback\\_ticked\\_success}\;\widehat{=}$ $\textbf{any}\;\;node$
where $n\\_type(node)=FALLBACK$ $n\\_tick(node)=TRUE$
$n\\_result(node)=UNKNOWN$ $\exists x\cdot(x\in parent^{-1}[\\{node\\}]\land
n\\_tick(x)=TRUE\land\\\ \forall y\cdot(y\in parent^{-1}[\\{node\\}]\land
n\\_tick(y)=TRUE\implies\\\ n\\_id(y)\leq n\\_id(x))\land
n\\_result(x)=SUCCESS)$ then $n\\_result(node):=SUCCESS$
$analyzing\\_subtree(parent(node)):=FALSE$ end
$\mathsf{fallback\\_ticked\\_continue}\;\widehat{=}$
$\textbf{any}\;\;node,child$ where $n\\_type(node)=FALLBACK$
$n\\_tick(node)=TRUE$ $n\\_result(node)=UNKNOWN$ $\exists x\cdot(x\in
parent^{-1}[\\{node\\}]\land n\\_tick(x)=TRUE\land\\\ \forall y\cdot(y\in
parent^{-1}[\\{node\\}]\land n\\_tick(y)=TRUE\implies\\\ n\\_id(y)\leq
n\\_id(x))\land n\\_result(x)=FAILURE)$ $\\{x|\;x\in
parent^{-1}[\\{node\\}]\;\\\ \hskip
10.0pt\land\;n\\_tick(x)=FALSE\\}=\emptyset$ $node\in nodes$ $child\in\\{x\in
parent^{-1}[\\{node\\}]\\\ \hskip 10.0pt\land n\\_tick(x)=FALSE\land\\\
\forall y\cdot(y\in parent^{-1}[\\{node\\}]\land n\\_tick(y)=FALSE\implies\\\
n\\_id(y)\geq n\\_id(x))$ $analyzing\\_subtree(node):=FALSE$ then
$n\\_tick(child):=TRUE$ $analyzing\\_subtree(node):=TRUE$ end end

Figure 5: The events to manage the Fallback node

At every event either we return a result or we tick a child node. And we set
the analyzing_subtree label true or false accordingly.

For the Sequence node, we follow a similar approach we the modified event to
control it according to its definition.
The definition is ”The Sequence node […], corresponds to routing the ticks to
its children from the left until it finds a child that returns either Failure
or Running, then it returns Failure or Running accordingly to its own parent.
It returns success if and only if all its children return Success. Note that
when a child returns Running or Failure, the Sequence node does not route the
ticks to the next child (if any).” [10].

The guards and events are similar to the Fallback node but with different
logic. The events of the Sequence node can be seen in Figure 6.

Events $\mathsf{sequence\\_ticked\\_running}\;\widehat{=}$
$\textbf{any}\;\;node$ where $n\\_type(node)=SEQUENCE$ $n\\_tick(node)=TRUE$
$n\\_result(node)=UNKNOWN$ $\exists x\cdot(x\in parent^{-1}[\\{node\\}]\land
n\\_tick(x)=TRUE\land\\\ \forall y\cdot(y\in parent^{-1}[\\{node\\}]\land
n\\_tick(y)=TRUE\implies\\\ n\\_id(y)\leq n\\_id(x))\land
n\\_result(x)=RUNNING)$ then $n\\_result(node):=RUNNING$
$analyzing\\_subtree(parent(node)):=FALSE$ end
$\mathsf{sequence\\_ticked\\_failure}\;\widehat{=}$ $\textbf{any}\;\;node$
where $n\\_type(node)=SEQUENCE$ $n\\_tick(node)=TRUE$
$n\\_result(node)=UNKNOWN$ $\exists x\cdot(x\in parent^{-1}[\\{node\\}]\land
n\\_tick(x)=TRUE\land\\\ \forall y\cdot(y\in parent^{-1}[\\{node\\}]\land
n\\_tick(y)=TRUE\implies\\\ n\\_id(y)\leq n\\_id(x))\land
n\\_result(x)=FAILURE)$ then $n\\_result(node):=FAILURE$
$analyzing\\_subtree(parent(node)):=FALSE$ end
$\mathsf{sequence\\_ticked\\_success}\;\widehat{=}$ $\textbf{any}\;\;node$
where $n\\_type(node)=SEQUENCE$ $n\\_tick(node)=TRUE$
$n\\_result(node)=UNKNOWN$ $\exists x\cdot(x\in parent^{-1}[\\{node\\}]\land
n\\_tick(x)=TRUE\land\\\ \forall y\cdot(y\in parent^{-1}[\\{node\\}]\land
n\\_tick(y)=TRUE\implies\\\ n\\_id(y)\leq n\\_id(x))\land
n\\_result(x)=SUCCESS)$ $\\{x|\;x\in parent^{-1}[\\{node\\}]\;\\\ \hskip
10.0pt\land\;n\\_tick(x)=FALSE\\}=\emptyset$ $node\in nodes$ then
$n\\_result(node):=SUCCESS$ $analyzing\\_subtree(parent(node)):=FALSE$ end
$\mathsf{sequence\\_ticked\\_continue}\;\widehat{=}$
$\textbf{any}\;\;node,child$ where $n\\_type(node)=SEQUENCE$
$n\\_tick(node)=TRUE$ $n\\_result(node)=UNKNOWN$ $\exists x\cdot(x\in
parent^{-1}[\\{node\\}]\land n\\_tick(x)=TRUE\land\\\ \forall y\cdot(y\in
parent^{-1}[\\{node\\}]\land n\\_tick(y)=TRUE\implies\\\ n\\_id(y)\leq
n\\_id(x))\land n\\_result(x)=SUCCESS)$ $\\{x|\;x\in
parent^{-1}[\\{node\\}]\;\\\ \hskip
10.0pt\land\;n\\_tick(x)=FALSE\\}=\emptyset$ $node\in nodes$ $child\in\\{x\in
parent^{-1}[\\{node\\}]\\\ \hskip 10.0pt\land n\\_tick(x)=FALSE\land\\\
\forall y\cdot(y\in parent^{-1}[\\{node\\}]\land n\\_tick(y)=FALSE\implies\\\
n\\_id(y)\geq n\\_id(x))$ $analyzing\\_subtree(node):=FALSE$ then
$n\\_tick(child):=TRUE$ $analyzing\\_subtree(node):=TRUE$ end end

Figure 6: The events to manage the Sequence node

###  3.4 Instantiating generic specification

The fourth refinement finally allows us to deal with real instances of BTs. At
this level, actual Behaviour Trees can be verified and checked for invariants
of interest for the agent. We now explain how the context and the machine need
to be built to properly check an existing tree.

The context contains the static part of the machine, the BT in itself is a
static part of the agent, and for this reason, the definition of the topology
and of the nodes type of the BT needs to be defined here.

To define a BT instance we start by extending the BT_level and listing all the
nodes contained in the tree (excluding the root, already defined in the
previous refinement), as constants.

We then need to define the axioms regarding the type: that every node defined
in the constant part is of type nodes. After this we need to initialize all
the nodes giving to them a n_type, and a n_id. Finally, we construct the
parent function, defining for each node which node is its parent, defining in
this way the topology of the tree.

Numbering the nodes correctly is essential to correctly navigate a tree. To
correctly give a n_id to the nodes, is required to follow the same principle
as a breadth-first search: firstly considering the depth, and then starting
counting the children from left to right. Therefore, The root, receive the
number $0$. We will then give to all its children the next numbers starting by
counting from the left-most node (i.e. the first to be ticked). When all the
children have been assigned a number we start to analyze the children of the
first child, numbering them left-to-right in the same way as before. And so
on, until all the nodes have been numbered.

Moving on the machine of the BT instance, we start as a refinement of the
BT_machine of the previous level. The only effort required here is the
extension and refinement of the events of the leaf nodes. Since the leaf nodes
(conditions and actions) are very specific to the application, we need to
extend the corresponding events for any action and condition we want to check.
For instance, let’s suppose that our agent has two condition nodes:
condition_1 and condition_2, the first condition checks if _the agent is more
than 3 meters from an obstacle_ , and the second that _the battery level is
more than 30%_. To cover the two conditions we would need to extend the
condition_ticked event of the BT level four times:  1.
condition_1_ticked_success,  2. condition_1_ticked_failure,  3.
condition_2_ticked_success,  4. condition_2_ticked_failure,  covering all the
possible results that the nodes can output. The modifications to be added
would be minimal, in the condition_1 case we would just need a guard to check
the condition (something like $distance\\_to\\_obstacle\geq 3$) and an action
that explain which value to return ($n\\_result(node):=SUCCESS$). For the
failure case, we would add the condition for failure
($distance\\_to\\_obstacle<3$) and the corresponding returning value
($n\\_result(node):=FAILURE$)

Moreover, if we want to check some safety invariant or other parameters, we
would need to model the part of the environment we are interested to check and
how the various actions would change the environment itself. For instance, if
a robot is moving in an environment with some obstacles we could define two
variables (pos_x, pos_y to define the position of the robot and another
variable (distance_to_obstacle), to keep track of the distance to the closest
obstacle. These would then be modified by actions. For instance, if we have an
action ”move one step up”, in the event of the corresponding action we would
also need to add something like: $pos\\_y:=pos\\_y+1$.

##  4 Case Study – A moving robot

As a case study, we analyzed a very simple case, considering a robot that is
moving toward a wall. The goal is to avoid the robot from going too close to a
certain distance to the obstacle Fig: 7.

Figure 7: Case study of a robot moving towards a wall, without going too close
to it. The robot can move closer to the wall, but it is not allowed to enter
the unsafe area.

To model this simple environment we have defined in the machine different
variables:  1. distance_to_object,  2. time,  3. prev_time . We use the
distance_to_object variable to identify the position of the robot and hence
its distance to the wall, while a _time_ variable (together with the previous
time step) is used to control the real-time passed during the execution.

To define the behavior of the robot we have defined BT represented in Figure
8.

Figure 8: Behaviour tree used to model the decision-making process of the
agent.

As invariant to be checked we set:

$distance\\_to\\_object\geq 3$

In the initialization, we set the robot at a distance of 10 meters from the
wall.

The first leaf node of the tree is a condition that the robot needs to be at
least 5 meters from the wall. If this condition holds then the robot can
perform the action of moving one step closer to the wall.

To model this environment, we followed the same approach described at the
instance level above. We first defined the tree in the context described
below. In the machine, we define the variables described above to define the
environment and we refine the leaf node events.

For the condition, we create the two refinements of the event for the
$SUCCESS$ and $FAILURE$ case as described in the methodology section. For the
action, we simplify the scenario enforcing that, once selected, the action
will take effect instantly and that can never fail. Allowing us to refine the
action event just for the $SUCCESS$ case.

$\textbf{context}\;robot\\_wall\\_context\;$
$\textbf{extends}\;03\\_BT\\_context\;$ constants $sequence\\_1$
$condition\\_1$ $action\\_1$ axioms $sequence\\_1\in nodes$ $condition\\_1\in
nodes$ $action\\_1\in nodes$ $n\\_type=\\{root\mapsto
ROOT,\;sequence\\_1\mapsto SEQUENCE,\;condition\\_1\mapsto
CONDITION,\;action\\_1\mapsto ACTION\\}$ $n\\_id=\\{root\mapsto
0,\;sequence\\_1\mapsto 1,\;condition\\_1\mapsto 2,\;action\\_1\mapsto 3\\}$
$parent=\\{sequence\\_1\mapsto root,\;condition\\_1\mapsto
sequence\\_1,\;action\\_1\mapsto sequence\\_1\\}$

Figure 9: The context of the BT for the robot moving towards the wall

We can note that the action ticked event is very similar to the previous
definition, with the only difference due to the addition of the last two
propositions.

$\mathsf{action\\_ticked\\_success}\;\widehat{=}$
$\textbf{REFINES}\;action\\_ticked$ $\textbf{any}\;\;node$ where
$n\\_type(node)=ACTION\;\text{// not changed}$ $n\\_tick(node)=TRUE\;\text{//
not changed}$ $n\\_result(node)=UNKNOWN\;\text{// not changed}$ $node\in
nodes\;\text{// not changed}$ then
$analyzing\\_subtree(parent(node)):=FALSE\;\text{// not changed}$
$n\\_result(node):=SUCCESS\;\text{// added in refinement}$
$distance\\_to\\_object:=distance\\_to\\_object-1\;\text{// added in
refinement}$ end

Figure 10: The context of the BT for the robot moving towards the wall

We then added a statement to increase the timestep whenever the result reaches
the root node, considering an entire exploration of the tree as a timestep in
the real world since the changes in the position of the robot have been
applied to the action event.

In general, just with these modifications at the machine from the BT level, we
are able to ensure the safety properties of the system by proving the safety
invariant that we have defined for this BT. This ensures us that, after the
initialization, the safety statement that we have described will hold for each
possible state that the BT can reach. This can be expanded, with multiple
safety statements and multiple requirements, allowing us to formally prove the
safety of the system that will implement this specific BT.

##  5 Related Works

Although Behavior Trees have started to become widely used in various robotics
applications, still little research has been done to increase the
dependability of robots using them. The importance of providing safe and
reliable agents becomes more and more critical, especially given the fact that
these tools are starting to be used at a larger level, especially with the
support of open source software like the Nav2 package for ROS [14]. Work in
the safety assessment of Behavior Trees has been done, for instance, [13], but
these studies normally refer to specific cases and become difficult to
generalize for other BT. BT can also be safely generated as in [18, 7] and
that BT is an optimal tool to build safe and reliable systems [9]. However,
formal verification of safety properties is rarely done on new trees
especially in manually constructed ones, risking in this way the overall
dependability of the systems using them. Other works propose the formal
verification of BT using Linear Temporal Logic [5], but they do not provide a
tool to support the verification process.

Event-B was used for modelling robotic systems by Troubitsyna et al. [12, 19,
16]. These works adopted goal-oriented model as a basis of robot’s decision
making and demonstrated goal reachability despite robot’s failures. However,
these works rely on an assumption that the robot’s operating environment
remains static, i.e., once a certain goal is accomplished, its status remains
unchanged throughout the entire system execution. In our work, reliance on BT
allows us to adopt dynamic model of system environment, because at each tick
the status of all tasks is getting checked, i.e., this work offers a more
flexible basis for modeling robot’s decision making.

##  6 Conclusion

In this work, we present a formal specification of Behaviour Trees in Event-B.
We also provide a methodology to test new instances of BT in a fairly simple
manner for newcomers hiding the definition of the BT at more abstract levels,
allowing a user to focus just on test and checking the properties the specific
BT needs to hold. Our approach aimed to facilitate the formal verification of
Behaviour Trees without the need to build the formalization of the tree and
the control flow node from scratch and simplifying also the modeling of the
leaf nodes.

In future works, we propose to improve the actual specification by adding the
support for parallel nodes and the most common decorator used in practice.
Moreover, we would like to simplify further the possibility to check already
used BT, creating a tool to automatically import trees built using already
used libraries like py_tree [17] or Behaviour_Tree.CPP [4]. Of particular
interest is also the analysis of already used BT in robotic frameworks like
the ones used in the Navigation 2 [14] package in ROS.

## References

  * [1]
  * [2] _Event-B platform,_.  Available at http://www.event-b.org/.
  * [3] Jean-Raymond Abrial (2010): _Modeling in event-B: System and software engineering_.  Cambridge University Press, 10.1017/CBO9781139195881.
  * [4] BehaviorTree: _Behaviortree/BehaviorTree.CPP: Behavior trees library in C++. batteries included._ Available at https://github.com/BehaviorTree/BehaviorTree.CPP.
  * [5] Oliver Biggar & Mohammad Zamani (2020): _A Framework for Formal Verification of Behavior Trees With Linear Temporal Logic_.  IEEE Robotics and Automation Letters 5(2), pp. 2341–2348, 10.1109/LRA.2020.2970634.  Conference Name: IEEE Robotics and Automation Letters.
  * [6] Michele Colledanchise, Alejandro Marzinotto & Petter Ogren (2014): _Performance analysis of stochastic behavior trees_.  pp. 3265–3272, 10.1109/ICRA.2014.6907328.
  * [7] Michele Colledanchise, Richard M. Murray & Petter Ögren (2017): _Synthesis of Correct-by-Construction Behavior Trees_.  In: 2017 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), IEEE Press, p. 6039–6046, 10.1109/IROS.2017.8206502.
  * [8] Michele Colledanchise, Ramviyas Parasuraman & Petter Ögren (2015): _Learning of Behavior Trees for Autonomous Agents_.  CoRR abs/1504.05811, 10.1109/TG.2018.2816806.
  * [9] Michele Colledanchise & Petter Ögren (2014): _How Behavior Trees modularize robustness and safety in hybrid systems_.  In: 2014 IEEE/RSJ International Conference on Intelligent Robots and Systems, pp. 1482–1488, 10.1109/IROS.2014.6942752.
  * [10] Michele Colledanchise & Petter Ögren (2018): _Behavior Trees in Robotics and AI: An Introduction_.  CRC Press, 10.1201/9780429489105.  Available at http://arxiv.org/abs/1709.00084.
  * [11] Razan Ghzouli, Thorsten Berger, Einar Broch Johnsen, Swaib Dragule & Andrzej Wąsowski (2020): _Behavior trees in action: a study of robotics applications_.  In: Proceedings of the 13th ACM SIGPLAN International Conference on Software Language Engineering, SLE 2020, Association for Computing Machinery, New York, NY, USA, pp. 196–209, 10.1145/3426425.3426942.
  * [12] Linas Laibinis, Inna Pereverzeva & Elena Troubitsyna (2017): _Formal reasoning about resilient goal-oriented multi-agent systems_.  Science of Computer Programming 148, pp. 66–87, 10.1016/j.scico.2017.05.008.  Available at https://www.sciencedirect.com/science/article/pii/S0167642317301181.  Special issue on Automated Verification of Critical Systems (AVoCS 2015).
  * [13] Peter A. Lindsay, Kirsten Winter & Nisansala Yatapanage (2010): _Safety Assessment Using Behavior Trees and Model Checking_.  In: 2010 8th IEEE International Conference on Software Engineering and Formal Methods, IEEE, Pisa, Italy, pp. 181–190, 10.1109/SEFM.2010.23.  Available at http://ieeexplore.ieee.org/document/5637427/.
  * [14] Steve Macenski, Francisco Martín, Ruffin White & Jonatan Ginés Clavero (2020): _The Marathon 2: A Navigation System_.  In: 2020 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 2718–2725, 10.1109/IROS45743.2020.9341207.  Available at http://arxiv.org/abs/2003.00368.
  * [15] W.F. McColl & K. Noshita (1986): _On the number of edges in the transitive closure of a graph_.  Discrete Applied Mathematics 15(1), pp. 67–73, 10.1016/0166-218X(86)90020-X.  Available at https://www.sciencedirect.com/science/article/pii/0166218X8690020X.
  * [16] Inna Pereverzeva, Elena Troubitsyna & Linas Laibinis (2012): _Formal Goal-Oriented Development of Resilient MAS in Event-B_.  In Mats Brorsson & Luís Miguel Pinho, editors: Reliable Software Technologies – Ada-Europe 2012, Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 147–161. 10.1007/978-3-642-30598-6_11.
  * [17] Splintered-Reality: _Splintered-reality/py_trees: Python implementation of behaviour trees._ Available at https://github.com/splintered-reality/py_trees.
  * [18] Tadewos G. Tadewos, Laya Shamgah & Ali Karimoddini (2019): _Automatic Safe Behaviour Tree Synthesis for Autonomous Agents_.  In: 2019 IEEE 58th Conference on Decision and Control (CDC), pp. 2776–2781, 10.1109/CDC40024.2019.9030183.
  * [19] Inna Vistbakka & Elena Troubitsyna (2021): _Modelling resilient collaborative multi-agent systems_.  Computing 103, pp. 535–557, 10.1007/s00607-020-00861-2.  Available at https://link.springer.com/article/10.1007/s00607-020-00861-2.
  * [20] Michael Waltham & Deshen Moodley (2016): _An Analysis of Artificial Intelligence Techniques in Multiplayer Online Battle Arena Game Environments_.  In: Proceedings of the Annual Conference of the South African Institute of Computer Scientists and Information Technologists, SAICSIT ’16, Association for Computing Machinery, New York, NY, USA, 10.1145/2987491.2987513.