# Validation of Neural Network Controllers for Uncertain Systems Through Keep-
Close Approach: Robustness Analysis and Safety Verification

Abdelhafid Zenati
Department of Engineering,
School of Mathematics, Computer Science and Engineering,
City, University of London Northampton Square,
London, UK, EC1V 0HB.
<EMAIL_ADDRESS>
& Nabil Aouf
School of Mathematics, Computer Science and Engineering,
City, University of London Northampton Square,
London, UK, EC1V 0HB.
<EMAIL_ADDRESS>
& David Sanchez de la Llana
European Space Research and Technology Centre,
ESA, 2201 AZ Noordwijk, The Netherlands.
<EMAIL_ADDRESS>
& Samir Bennani
European Space Research and Technology Centre,
ESA, 2201 AZ Noordwijk, The Netherlands.
<EMAIL_ADDRESS>
The main author

###### Abstract

Among the major challenges in neural control system technology is the
validation and certification of the safety and robustness of neural network
(NN) controllers against various uncertainties including unmodelled dynamics,
nonlinearities, and time delays. One way in providing such validation
guarantees is to maintain the closed-loop system output with a NN controller
when its input changes within a bounded set, close to the output of a robustly
performing closed-loop reference model. This paper presents a novel approach
to analysing the performance and robustness of uncertain feedback systems with
NN controllers. Due to the complexity of analysing such systems, the problem
is reformulated as the problem of dynamical tracking errors between the
closed-loop system with a neural controller and an ideal closed-loop reference
model. Then, the approximation of the controller error is characterised by
adopting the differential mean value theorem (DMV) and the Integral Quadratic
Constraints (IQCs) technique. Moreover, the Relative Integral Square Error
(RISE) and the Supreme Square Error (SSE) bounded set are derived for the
output of the error dynamical system. The analysis is then performed by
integrating Lyapunov theory with the IQCs-based technique. The resulting
worst-case analysis provides the user a prior knowledge about the worst case
of RISE and SSE between the reference closed-loop model and the uncertain
system controlled by the neural controller. By applying the novel approach
presented in this paper to two different problems, the theoretical results
have been validated. The first problem, the Single-Link Robot Arm Control
Problem, focused on a single output, while the second problem, the Deep
Guidance and Control of Apollo Lander, dealt with multiple outputs. As a
result of these practical tests, the proposed approach has proven to be
effective and applicable to various scenarios.

_Keywords_ NN Controller $\cdot$ Lyapunov $\cdot$ IQCs $\cdot$ Robustness
Analysis $\cdot$ Safety and Verification.

##  1 Introduction

Neural networks (NNs) show great performances when adopted for tasks where
Artificial Intelligence (AI) based solutions are considered. Due to the latest
advancements in AI, stabilising dynamical systems with NN controllers is
invigorated (Schulman et al., 2015; Pauli et al., 2022). However, feedback
systems with NN controllers still suffer from a lack of robustness, safety
certificates and, validation due to their complex, opaque, and large-scale
nature structures (Fazlyab et al., 2020; Deka et al., 2022). Particularly, NN
controllers have many types of nonlinear activation functions (Yin et al.,
2021b) in addition to their vulnerability to adversarial attacks and system
uncertainties (Szegedy et al., 2013; Zheng et al., 2016). These drawbacks make
the application of classical stability analysis methods such as Lyapunov
theory to NN controller-based closed-loop systems difficult. This leads to
limited adoption of NN controllers in safety-critical applications (Khalil,
2015; Xiang et al., 2018a).

Over the course of the last three decades, the framework of uncertain system
systematic analysis via Integral Quadratic Constraints (IQCs) was developed
and used (Megretski and Rantzer, 1997; Veenman and Scherer, 2013). The IQCs
approach is proven to be efficient in capturing the properties of various
classes of uncertainties such as slope-restricted nonlinearities, dynamics
that can not be modeled,and time delays in addition to systems’ nonlinearities
(Lessard et al., 2016; Veenman et al., 2016). The flexibility of the approach
and its generality made it one of the most important techniques utilized to
perform and assess an uncertain dynamical system’s stability analysis and
performance, respectively (Wang and Manchester, 2019). In (Iannelli et al.,
2019), the problem of stability analysis of nonlinear uncertain systems is
formulated using the definition of the augmented plant via the state-space
factorization of the IQC. The certificates for the region of attraction are
established with both soft and hard IQC factorizations. Then, the formulated
problem is solved using a numerical approach such as Sum of Squares (SOS).
Based on solving Linear Matrix Inequalities (LMIs), (Fetzer et al., 2017)
proposes a novel dynamic multipliers-dependent shift of the Kalman-Yakubovich-
Popov (KYP) (Rantzer, 2015) certificate that is expressed using constraints.
The novel multipliers allow the formulation of the regional stability
properties within IQCs theory, the verification of output constraints and the
hard time-domain state. (Summers and Packard, 2010) derives the certification
conditions for local $\mathcal{L}_{2}$ gain of locally stable interconnected
systems. Using estimation of reachability, the approach enhances the local
$\mathcal{L}_{2}$ gain calculation. Therefore, a set of local (IQCs) can be
established for a fixed nonlinear dynamical system. Then by the local IQCs, a
local $\mathcal{L}_{2}$ gain is obtained. The possibility of using IQCs in the
analysis systems subject to saturated linear feedback and the design of its
control is the objective of the work in (Fang et al., 2008). The research
there focuses on establishing a IQCs-based conditions set under which an
ellipsoid is contractively invariant for such a system. The paper results show
that these set invariant conditions are necessary for such a study.

Several formal verification methods recently propose to investigate the
stability of closed-loop systems interconnected with NN controllers. In (Xiang
et al., 2018b), the authors study the verification problems for a class of
piece-wise linear systems with NN controllers. By making an extension to the
reachable set estimation, a reduction of the safety verification is
established to check for empty junctions between unsafe regions and the
reachable set. In (Ivanov et al., 2020), converting the NN controller into a
hybrid system equivalent to the original one, the work addresses the issue of
safety guarantees for autonomous systems with NN controllers. Moreover, the
scalability of the proposed approach using Taylor series approximation with
worst-case error bounds is looked at. Using a similar approach as (Ivanov et
al., 2020), the work in (Ivanov et al., 2021) emphasizes on NNs with both tanh
and sigmoid activation functions. Through the Taylor model preconditioning and
shrink wrapping, a Taylor-model-based reachability algorithm is developed and
a parallel implementation is provided. (Yin et al., 2021a) presents an
approach, which permits the closed loop system with NN controllers to involve
different types of perturbations by using IQCs technique and to capture their
input/output dynamical behavior.

In this paper, we aim at certifying the safety and providing robustness
guarantees for NN controllers against various systems uncertainties including
unmodelled dynamics, slope-restricted non-linearities, and time delays.
Inspiring from comparison principle (Nam et al., 2021; et all, 2022), the core
idea of this work is to provide guarantees in maintaining the output of the
uncertain system with a trained NN controller close to an ideal reference
closed-loop model when its inputs change within a bounded set. To this end, a
novel approach to analyse the dynamical error stability of the feedback system
with a NN controller and a reference closed-loop model is presented. Since
analysing such systems is complicated, the problem is reformulated as the
dynamical tracking errors between the uncertain interconnected system with a
NN and the reference closed-loop model. Approximation of the neural controller
error is then performed using the differential mean value theorem (DMV) and
integral quadratic constraints (IQCs) in order to eliminate the neural
controller’s nonlinearity. As the time integral of the square of the error
(ISE) is used to evaluate the system’s performance (Atkins et al., 2013), the
RISE and SSE bounded sets are derived based on generating dissipation-
inequality conditions. These conditions of the RISE and SSE bound the worst-
case $\mathcal{L}_{2}\rightarrow\mathcal{L}_{2}$ and
$\mathcal{L}_{2}\rightarrow\mathcal{L}_{\infty}$ induced gain of the dynamical
error system, respectively. In addition, Lyapunov theory is integrated with
IQCs-based techniques to accomplish the analysis of the worst-case gains. This
analysis determines the worst-case ISE between the closed-loop reference model
and the neural-controlled uncertain system. Using the difference between the
closed-loop reference system and interconnected system with a NN is an
appealing approach for the following reasons: (i) it can eliminate the non-
linearity of the NN controller by applying the Differential Mean Value theorem
which is more suitable and easier to analyse. (ii) The work in (Pfifer and
Seiler, 2015) considered the trivial initial conditions to facilitate the
analysis. However, these conditions limit the application of their results to
real challenging problems. Indeed, in the real world, the initial condition
has a big effect on the stability of the system and generally, it is not null.
Our keep-close technique can function on any initial condition of the system.
Using the difference between the closed-loop reference system and the
interconnected system with a NN controller, the initial conditions become null
which makes the analysis easier, plausible, and more general. In (Fazlyab et
al., 2020), the authors analysed the stability of feedback systems with NN
controllers using IQCs approach to capture their input/output behavior. Their
theoretical result relies on semi-definite programming to estimate the
Attraction Region (ROA) of the equilibrium point. Our keep-close approach is
more general as it covers the tracking in general where the convergence to an
equilibrium point is a particular case. Moreover, from verification and safety
point of view, the keep-close approach is more practical than the result in
(Fazlyab et al., 2020) because it can indicate the occurrence likelihood of
any undesirable behavior before the convergence to the equilibrium point. The
novel approach outlined in this paper has been implemented and tested in
practice. Validation was carried out on two specific problems: the Single-Link
Robot Arm Control Problem, which focused on a single output scenario, and the
Deep Guidance and Control of Apollo Lander, which dealt with multiple outputs.
These tests provide empirical support for the efficacy and adaptability of the
proposed approach across different application domains.

The rest of this paper is organized as follows: The problem formulation and
preliminaries are given in Section 2. Section 3 is devoted to the robustness
analysis when the system is submitted to uncertainties. In Section 3, a
numerical example is provided to illustrate the results of our proposed
solution. Conclusions are given in Section 4.

###  1.1 Notation

  * •

$\mathbb{R}$, $\mathbb{R}_{+}$, and $\mathbb{C}$ refer to the set of real,
non-negative real, and complex numbers, respectively. Also, in the complex
plane, the notation $\mathbb{RH}$ denote to all analytic functions within the
closed exterior of the unit disk. $M^{T}$ and $M^{*}$ are the transpose and
complex conjugate transpose of the matrix $M$, respectively.

  * •

Consider the LPV system $\Sigma$ where its matrices of state space depend on a
time-varying parameter vector $s(t):\mathbb{R}_{+}\rightarrow\mathcal{S}$ and
$\mathcal{S}\subset\mathbb{R}^{n}$. $s(t)$ a continuously differentiable
function of time where the parameter rates of variation
$\dot{s}(t):\mathbb{R}_{+}\rightarrow\dot{\mathcal{S}}$. $\dot{\mathcal{S}}$
is a hyperrectangle given by:

$\dot{\mathcal{S}}=\left\\{\rho\in\mathbb{R}^{n}~{}|~{}\underline{\rho}_{i}\leq\rho_{i}\leq\overline{\rho}_{i},~{}i=1\cdots
n\right\\}$ (1)

Let a parameter-dependent matrix $P$ a continuously differentiable function of
the parameter $s$ gievn by: $P:\mathcal{S}\rightarrow\mathbb{S}^{n}$,
$\mathbb{S}^{n}$ refers to the set of $n\times n$ symmetric matrices. The
differential operator $\partial
P:\mathcal{S}\times\dot{\mathcal{S}}\rightarrow\mathbb{S}^{n}$ is given by:

$\partial P(s,\dot{s})=\sum\limits_{k=1}^{n}\dfrac{\partial P(s)}{\partial
s_{k}}\dot{s}_{k}$ (2)

  * •

The space $\mathcal{L}_{2}$ denotes the set of functions
$\varphi:\mathbb{R}_{+}\rightarrow\mathbb{R}^{n}$ satisfying
$||\varphi||_{2}<\infty$

$\displaystyle||\varphi||_{2}=\left[\int\limits_{0}^{+\infty}\ \
\varphi(\tau)^{T}\varphi(\tau)\,\mathrm{\textit{d}\tau}\right]^{\frac{1}{2}}$
(3)

  * •

The space $\mathcal{L}_{\infty}$ refers to the set of functions
$\varphi:\mathbb{R}_{+}\rightarrow\mathbb{R}^{n}$ satisfying
$||\varphi||_{\infty}<\infty$

$\displaystyle||\varphi||_{\infty}=\sup_{t\in[0,\infty)}\left[\varphi(t)^{T}\varphi(t)\right]^{\frac{1}{2}}$
(4)

  * •

We recall here the Differential Mean Value (DMV) Theorem:

######  Lemma 1 (DMV Theorem (Zemouche et al., 2005; Sharma et al., 2016))

Let $y_{1},y_{2},x\in\mathbb{R}^{n}$ and
$\Upsilon:~{}\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ be a continuous and
differentiable function on convex hull of the set $\mathbf{Co}(y_{1},y_{2})$.
Then, there are constant vectors
$s_{1},\cdots,s_{m}\in\mathbf{Co}(y_{1},y_{2})$, $s_{i}\neq y_{1},~{}s_{i}\neq
y_{2}$ , $i=1,\cdots,n$ such that:

$\Upsilon(y_{1})-\Upsilon(y_{2})=\left(\sum\limits_{i=1}^{m}\sum\limits_{j=1}^{n}e_{m}(i)e_{n}^{T}(j)\dfrac{\partial\Upsilon_{i}}{\partial
x_{j}}(s_{i})\right)(y_{1}-y_{2})$ (5)

where $\mathbf{Co}(y_{1},y_{2})={\alpha y_{1}+(1-\alpha)y_{2},0\leq\alpha\leq
1}$ is the convex domain of $y_{1},y_{2}$ and $e_{m}(i)$ is a the canonical
basis vector of $\mathbb{R}^{m}$ with 1 at $i-th$ component with

$e_{m}(i)=(0,\cdots,0,1,0,\cdots,0)^{T}\in\mathbb{R}^{m}$ (6)

  * •

As the aim in this work is to estimate the Integral Square Error (ISE) and
worst expected Supreme Square Error (SSE), it is worth to define those metrics
as following:

######  Definition 1 ( RISE & SSE)

Let us consider $y$ and $\hat{y}$ are the outputs of the closed-loop
interconnected-based NN system and the closed-loop reference model,
respectively. Then the Relative Integral Square Error (RISE) and the Supreme
Square Error (SSE) are given by :

$RISE\triangleq\dfrac{||y-\hat{y}||_{2}}{\sqrt{||y||_{2}^{2}+||\hat{y}||_{2}^{2}}},~{}SSE\triangleq\dfrac{||y-\hat{y}||_{\infty}}{\sqrt{||y||_{2}^{2}+||\hat{y}||_{2}^{2}}}$
(7)

where the worst case of the Relative Integral Square Error $\gamma$ and the
worst case of the Supreme Square Error $\sigma$ satisfy:

$\dfrac{||y-\hat{y}||_{2}}{\sqrt{||y||_{2}^{2}+||\hat{y}||_{2}^{2}}}\leq\gamma,~{}~{}\dfrac{||y-\hat{y}||_{\infty}}{\sqrt{||y||_{2}^{2}+||\hat{y}||_{2}^{2}}}\leq\sigma$
(8)

##  2 Problem formulation

###  2.1 Keep-Close Approach

A robustness analysis of feedback systems based on NN controllers against
potential uncertainties is presented here. To achieve this, the output of the
closed-loop system with a NN controller is kept close to the output of an
ideal and trusted closed-loop reference model when its input changes within a
bounded set. Consider the closed-loop interconnected- based NN system with the
uncertainty $\Delta$ in Figure 1, (i.e. the uncertain plant $F(P,\Delta)$ and
NN $\pi$). The perturbed plant $F(P,\Delta)$ is an interconnection of a
nominal plant $P$ and an uncertainty $\Delta$ which can be described
mathematically by the following dynamical equations:

$\Sigma_{F(P,\Delta)}^{\pi}:=\left\\{\begin{array}[]{l}\dot{x}(t)=Ax(t)+Bu(t)+\tilde{B}q_{\delta}(t)\\\
y(t)=Cx(t)+Du(t)+\tilde{D}q_{\delta}(t)\\\
q_{\delta}(t)=\Delta(p_{\delta}(t)=y(t))\\\
u(t)=\pi(y(t))\in\mathcal{L}_{2}\\\ x(0)=x_{0}\end{array}\right.$ (9)

where $x\in\mathbb{R}^{n}$ is the system’s state vector, control inputs of the
neural controller $u\in\mathbb{R}^{m}$ and $(A,B,C,D,\tilde{B},\tilde{D})$ are
known, and the perturbation operator $\Delta$ is a bounded, causal. In this
study, we aim to ensure that the closed-loop system with a NN controller
$\Sigma_{F(P,\Delta)}^{\pi}$ keeps its outputs close to a given reference
model $\Sigma_{Ref}^{\pi^{*}}$ (a closed-loop system with an ideal, and
classically proven controller). This reference model can be given by:

$\Sigma_{Ref}^{\pi^{*}}:=\left\\{\begin{array}[]{l}\dot{\hat{x}}(t)=A\hat{x}(t)+B\hat{u}(t)\\\
\hat{y}(t)=C\hat{x}(t)+D\hat{u}(t)\\\
\hat{u}(t)=\pi^{*}(\hat{y}(t),t)\in\mathcal{L}_{2}\\\
\hat{x}(0)=x_{0}\end{array}\right.$ (10)

with reference system’s state vector $\hat{x}\in\mathbb{R}^{n}$, control
inputs $\hat{u}\in\mathbb{R}^{m}$ and $(A,B,C,D)$ known as it is mentioned in
(9). To accomplish our aim stated above, an analysis of the dynamical error
system between the two closed-loop systems (i.e. the reference model and the
system closed with a NN controller) is presented in the following section. We
start by mathematically describing the dynamical error system. Let us assume
that the state of dynamical error system is given by
$\zeta(t)=x(t)-\hat{x}(t)$ and $z(t)=p_{\delta}(t)-\hat{y}(t)$, then,
$\dot{\zeta}(t)=\dot{x}(t)-\dot{\hat{x}}(t)$, where

$\Sigma_{\tilde{F}(\tilde{P},\Delta)}^{\mu}:=\left\\{\begin{array}[]{l}\dot{\zeta}(t)=A\zeta(t)+B\mu(t)+\tilde{B}q_{\delta}(t)\\\
z(t)=C\zeta(t)+D\mu(t)+\tilde{D}q_{\delta}(t)\\\
\mu(t)=\pi(y(t))-\pi^{*}(\hat{y}(t),t)\\\ \zeta(0)=0\end{array}\right.$ (11)

where $\mu(t)$ is the approximation of the controller error (between the NN
controller and the ideal controller adopted in the closed loop reference
model). Also, the initial condition of the dynamical error state $\zeta(t)$ is
$\zeta(0)=0$ because $x(0)=\hat{x}(0)=x_{0}$.

Figure 1: Feedback system with plant $F(P,\Delta)$ and NN $\pi$ with closed
loop Reference Model.

###  2.2 Integral Quadratic Constraints (IQCs) Concept

The key concept of IQCs is instead of investigating the dynamical behaviour of
the system that contains uncertainty $\Delta$, we analyze the system where
$\Delta$ is removed and the signals $(p(t),~{}q(t))$ are utilized to enforce
the following constraints:

$\displaystyle\int\limits_{-\infty}^{+\infty}\ \
\begin{bmatrix}\hat{p}(j\omega)\\\
\hat{q}(j\omega)\end{bmatrix}^{*}\Pi(j\omega)\begin{bmatrix}\hat{p}(j\omega)\\\
\hat{q}(j\omega)\end{bmatrix}\,\mathrm{\textit{d}\tau}\geq 0$ (12)

where $\Pi(j\omega)$ is named an “IQC multiplier” or can simply called a
“multiplier”, $\hat{p}(j\omega)$ and $\hat{q}(j\omega)$ are Fourier transforms
of $p(t)$ and $q(t)$, respectively. Since (12) is valid for all admissible
choices of $\Delta$, then all properties can be proved for the constrained
system are hold also for the original system. Therefore, the mathematical
relationship between input-output quantities of system components can be
described via IQCs. The IQC in (12) can be described in the time domain based
on Parseval Theorem (Kelkar et al., 1983). Assuming that $\Pi$ is function
that is uniformly bounded and rational. Then, $\Pi(j\omega)$ can be factorized
as:

$\Pi(j\omega)=\Psi(j\omega)^{*}\mathcal{M}\Psi(j\omega),$ (13)

where $\Psi(j\omega)\in\mathbb{RH}$ and $\mathcal{M}$ represents a constant
matrix (Veenman et al., 2016). $\Psi(j\omega)$ is expressed as:

$\Psi(j\omega):=C_{\Psi}j[\omega I-A_{\Psi}][B_{\Psi 1}~{}B_{\Psi 2}]+[D_{\Psi
1}~{}D_{\Psi 2}]$ (14)

Then, the IQC of (12) can be written as:

$\displaystyle\int\limits_{0}^{+\infty}\ \
r(\tau)^{T}\mathcal{M}r(\tau)\,\mathrm{\textit{d}\tau}\geq 0$ (15)

where $r(t)$ is the output of the following linear system:

$\Psi:=\left\\{\begin{array}[]{ll}\dot{\psi}(t)&=A_{\Psi}\psi(t)+B_{\Psi
1}p(t)+B_{\Psi 2}q(t)\\\ r(t)&=C_{\Psi}\psi(t)+D_{\Psi 1}p(t)+D_{\Psi
2}q(t)\\\ \psi(0)&=0\end{array}\right.$ (16)

The time-domain constraint (17) is generally satisfied over infinite time
intervals. A hard IQC satisfies the following more restrictive conditions: If
$\Delta$ is any causal and bounded operator that fulfills (12) then:

$\displaystyle\int\limits_{0}^{T}\ \
r(\tau)^{T}\mathcal{M}r(\tau)\,\mathrm{\textit{d}\tau}\geq 0$ (17)

holds for all $T\geq 0$. Conversely, the soft IQCs in the time-domain
constraint do not require to hold overall finite time intervals. This property
is critical because the major technical and pedagogical problems in the IQCs
framework occur when we utilize soft IQCs. Unfortunately, there is an
ambiguity surrounding both terms soft IQC and hard IQC. Particularly, the
unicity of factorization of $\Pi(j\omega)$ as
$\Psi^{*}(j\omega)\mathcal{M}\Psi(j\omega)$, is not unique. For clarity, the
following definition will be adopted in the rest of the paper.

######  Definition 2 ( Hard IQC (Pfifer and Seiler, 2015))

A uniformly bounded and rational function
$\Pi:j\mathbb{R}\rightarrow\mathbb{C}$ admits a hard IQC factorization if
there exists $\Psi(j\omega)\in\mathbb{RH}$ and matrix $\mathcal{M}$ with
$\mathcal{M}_{ij}\in\mathbb{C}$ , such that
$\Pi(j\omega)=\Psi^{*}(j\omega)\mathcal{M}\Psi(j\omega)$, and any bounded,
causal operator $\Delta$ which satisfies the IQC defined by $\Pi(j\omega)$
also satisfies:

$\displaystyle\int\limits_{0}^{T}\ \
r(\tau)^{T}\mathcal{M}r(\tau)\,\mathrm{\textit{d}\tau}\geq 0$ (18)

for all $T\geq 0$ and for all $p(t)\in\mathcal{L}_{2}([0~{}\infty))$,
$q=\Delta(p)$. $(\Psi(j\omega),\mathcal{M})$ is a hard IQCs factorization of
$\Pi$.

###  2.3 Controller Error Analysis

This subsection analyses the approximation of the controller error $\mu(t)$
(i.e. the difference between the NN controller and the ideal controller
adopted in the reference closed-loop model). The objective of this analysis is
to find a less conservative and better representative expression of $\mu(t)$
where we can extract and dissociate certain and uncertain quantities from
$\mu(t)$. Then, this expression can be directly used in the dynamic error
system (11). To reach the objective mentioned above let us start by
introducing the following lemma:

######  Lemma 2 ( Controller Error)

Consider the dynamics of the error in (11) and the bounded NN controller
$\pi$, the approximation error $\mu(t)$ given in (11) can be expressed as:

$\mu(t)=\left[I-\Lambda(s)D\right]^{-1}[\Lambda(s)C\zeta(t)+\Lambda(s)\tilde{D}\Delta(y(t))+\epsilon(\hat{y}(t))]$
(19)

if $D=0$, therefore

$\mu(t)=\Lambda(s)C\zeta(t)+\Lambda(s)\tilde{D}\Delta(y(t))+\epsilon(\hat{y}(t))$
(20)

and the vectors $s_{1},\cdots,s_{n}\in\mathbf{Co}(y(t),\hat{y}(t))$,
$s_{i}\neq y_{1},~{}s_{i}\neq y_{2}$ , $i=1,\cdots,n$ such that

$\displaystyle\Lambda(s)=\sum\limits_{i=1}^{m}\sum\limits_{j=1}^{n}e_{m}(i)e_{n}^{T}(j)\dfrac{\partial\pi_{i}}{\partial
y_{j}}(s_{i})$ (21)

where $e_{m}(i)=(0,\cdots,0,1,0,\cdots,0)^{T}\in\mathbf{R}^{m}$ is a canonical
basis vector of $\mathbb{R}^{m}$ with 1 at $i-th$ component.

##### Proof:

From (11), the error approximation $\mu(t)$ is given by:

$\mu(t)=\pi(y(t))-\pi^{*}(\hat{y}(t),t)$ (22)

Adding and subtracting $\pi(\hat{y}(t))$, the error approximation $\mu(t)$ can
be expressed by:

$\begin{array}[]{ll}\mu(t)&=\pi(y(t))-\pi^{*}(\hat{y}(t),t)+\pi(\hat{y}(t))-\pi(\hat{y}(t))\\\
&=\pi(y(t))-\pi(\hat{y}(t))+\epsilon(\hat{y}(t))\end{array}$ (23)

where the trained NN controller error $\epsilon(y(t))$ given by:

$\epsilon(\hat{y}(t))=\pi(\hat{y}(t))-\pi^{*}(\hat{y}(t),t)$ (24)

On the other hand, using the Differential Mean Value Theorem in lemma 1 and
the fact that the function $\pi$ is continuous and differentiable on convex
hull of the set $\mathbf{Co}(y(t),\hat{y}(t))$, then, there are vectors
$s_{1},\cdots,s_{n}\in\mathbf{Co}(y(t),\hat{y}(t))$, $s_{i}\neq
y,~{}s_{i}\neq\hat{y}$ , $i=1,\cdots,n$. Because
$z(t)=C\zeta(t)+D\mu(t)+\tilde{D}q_{\delta}(t)$, then

$\left.\begin{array}[]{l}\pi(y(t))-\pi(\hat{y}(t))=\Lambda(s)(y(t)-\hat{y}(t))=\Lambda(s)z(t)\\\
\\\
=\Lambda(s)\left[C\zeta(t)+D\mu(t)+\tilde{D}q_{\delta}(t)\right]\end{array}\right.$
(25)

where

$\displaystyle\Lambda(s)=\sum\limits_{i=1}^{m}\sum\limits_{j=1}^{n}e_{m}(i)e_{n}^{T}(j)\dfrac{\partial\pi_{i}}{\partial
x_{j}}(s_{i}),~{}q_{\delta}(t)=\Delta(y(t))$

The equations (23) and (25) allow us to conclude that:

$\mu(t)=\left[I-\Lambda(s)D\right]^{-1}[\Lambda(s)C\zeta(t)+\Lambda(s)\tilde{D}q_{\delta}(t)+\epsilon(\hat{y}(t))]$
(26)

where $\epsilon(\hat{y}(t))=\pi(\hat{y}(t))-\pi^{*}(\hat{y}(t),t)$.

$\square$

During the design of a control system, the feed-through matrix $D$ can be
adjusted to achieve desired control objectives and ensure that $I-\Lambda(s)D$
is invertible. Control strategies and the desired impact of inputs on outputs,
in conjunction with state feedback control, determine whether the feed-through
matrix is nonzero or zero. The feed-through matrix $D$ is usually a zero
matrix, i.e. $D=0$.

The result of (19) represents less conservative and better representative
expression of $\mu(t)$. The certain quantity $\Lambda(s)C\zeta(t)$ depends on
the weights of the neural controller and the state of the dynamical error
system $\zeta(t)$ which is similar to what we call within the control theory
community as state-feedback control. For the uncertain quantities, we have two
terms $\Lambda(s)D_{p}\Delta(y(t))$ and $\epsilon(\hat{y}(t))$ depending on
the outputs of the interconnected system with a NN controller and the closed
loop reference model, respectively. The term $\Lambda(s)D_{p}\Delta(y(t))$
represents the effect of the uncertainty of the plant on the controller
whereas the term $\epsilon(\hat{y}(t))$ illustrates the training error effect
on the dynamical error system which is very small and bounded because the
purpose of the training is to minimize this quantity $|\epsilon(\hat{y}(t))|$.

###  2.4 IQCs Characteristics of Uncertainties

The perturbation acting on the feedback system can represent diverse types of
uncertainty, such as unmodeled dynamics, saturation, time delay, and slope-
restricted nonlinearities. Knowing that the term
$\epsilon(\hat{y}(t))=\pi(\hat{y}(t))-\pi^{*}(\hat{y}(t),t)$ represents the
training error which is very small and bounded because the purpose of the
training is to minimize this quantity $\min|\epsilon(\hat{y}(t))|$. Also, the
input-output relationship of the uncertainties $\Delta(.)$ can be
mathematically described, as (18), by integral quadratic constraints (IQCs).
Accordingly, the IQCs technique can be adopted to characterise the
uncertainties $\epsilon(.)$ and $\Delta(.)$. This later leads to the following
assumption:

###### Assumption 1

Assume that the outputs of the closed-loop reference model and the closed-loop
reference system and the interconnected system with a NN controller
$\hat{y}(t)=p_{\epsilon}(t)\in\mathcal{L}_{2}$ and
$y(t)=p_{\delta}(t)\in\mathcal{L}_{2}$, having bounded, causal operators
$\Delta(.)$ and $\epsilon(.)$ which satisfy the IQCs defined by
$\Pi_{\epsilon}(j\omega)$ and $\Pi_{\delta}(j\omega)$, respectively. The
input-output relationship of the uncertainties $\Delta(.)$ and $\epsilon(.)$
can be mathematically described, as (18), by integral quadratic constraints
(IQCs) as follows: :

$\displaystyle\left\\{\begin{array}[]{ll}\int\limits_{0}^{T}\ \
r_{\delta}(\tau)^{T}\mathcal{M}_{\delta}r_{\delta}(\tau)\,\mathrm{\textit{d}\tau}\geq
0\\\ \int\limits_{0}^{T}\ \
r_{\epsilon}(\tau)^{T}\mathcal{M}_{\epsilon}r_{\epsilon}(\tau)\,\mathrm{\textit{d}\tau}\geq
0\end{array}\right.$ (29)

for all $T\geq 0$ and for all
$\hat{y}(t),~{}y(t)\in\mathcal{L}_{2}([0~{}\infty))$.
$(\Psi_{\epsilon}(j\omega),\mathcal{M}_{\epsilon})$ and
$(\Psi_{\delta}(j\omega),\mathcal{M}_{\delta})$ are a hard IQCs factorization
of $\Pi_{\epsilon}$ and $\Pi_{\delta}$, respectively.

###### Corollary 1

Consider Assumption 1, the signals $r_{\delta}(t)$ and $r_{\epsilon}(t)$ in
(29) are the output of the following LTI systems, respectively:

$\Psi_{\delta}:=\left\\{\begin{array}[]{ll}\dot{\psi}_{\delta}(t)&=A_{\delta}\psi_{\delta}(t)+B_{\delta
1}q_{\delta}(t)+B_{\delta 2}p_{\delta}(t)\\\
r_{\delta}(t)&=C_{\delta}\phi(t)+D_{\delta 1}q_{\delta}(t)+D_{\delta
2}p_{\delta}(t)\\\ \psi_{\delta}(0)&=0\end{array}\right.$ (30)

and

$\Psi_{\epsilon}:=\left\\{\begin{array}[]{ll}\dot{\psi}_{\epsilon}(t)&=A_{\epsilon}\psi_{\epsilon}(t)+B_{\epsilon
1}q_{\epsilon}(t)+B_{\epsilon 2}p_{\epsilon}(t)\\\
r_{\epsilon}(t)&=C_{\epsilon}\phi(t)+D_{\epsilon 1}q_{\epsilon}(t)+D_{\epsilon
2}p_{\epsilon}(t)\\\ \psi_{\epsilon}(0)&=0\end{array}\right.$ (31)

where $\Psi_{\delta}$ and $\Psi_{\epsilon}$ represent the ‘virtual’ filters as
in (16) applied to the inputs $p_{\delta}$ and $p_{\epsilon}$ and outputs
$q_{\delta}$ and $q_{\epsilon}$ of $\Delta$ and $\epsilon$ and a constraint on
the outputs $r_{\delta}$ and $r_{\epsilon}$ of (16). Thus, an extended filter
can be written as follows:

$\Xi:=\left\\{\begin{array}[]{ll}\dot{\xi}(t)&=A_{\xi}\xi(t)+B_{\xi
1}q(t)+B_{\xi 2}p(t)\\\ r(t)&=C_{\xi}\xi(t)+D_{\xi 1}q(t)+D_{\xi 2}p(t)\\\
\xi(0)&=0\end{array}\right.$ (32)

where

$\displaystyle\xi(t)=\begin{bmatrix}\psi_{\delta}(t)\\\
\psi_{\epsilon}(t)\end{bmatrix},~{}~{}p(t)=\begin{bmatrix}p_{\delta}(t)\\\
p_{\epsilon}(t)\end{bmatrix},~{}q(t)=\begin{bmatrix}q_{\delta}(t)\\\
q_{\epsilon}(t)\end{bmatrix}$

and the matrices of (32) are defined as follow:

$\displaystyle A_{\xi}=\begin{bmatrix}A_{\delta}&0\\\
0&A_{\epsilon}\end{bmatrix},B_{\xi 1}=\begin{bmatrix}B_{\delta 1}&0\\\
0&B_{\epsilon 1}\end{bmatrix},~{}B_{\xi 2}=\begin{bmatrix}B_{\delta 2}&0\\\
0&B_{\epsilon 2}\end{bmatrix}$ $\displaystyle
C_{\xi}=\begin{bmatrix}C_{\delta}&0\\\ 0&C_{\epsilon}\end{bmatrix},D_{\xi
1}=\begin{bmatrix}D_{\delta 1}&0\\\ 0&D_{\epsilon 1}\end{bmatrix},~{}D_{\xi
2}=\begin{bmatrix}D_{\delta 2}&0\\\ 0&D_{\epsilon 2}\end{bmatrix}$

##  3 Robust Tracking Analysis via Standard IQCs

Considering the error dynamics system of (11), the conservative expression of
the $\mu$ in (20), and let
$\chi(t)=\begin{bmatrix}\zeta(t)~{}\xi(t)\end{bmatrix}^{T}$ be the extended
average state vector with dynamics given by the following (omitting $s$ for a
shorter notation):

$\Sigma_{ex}\left\\{\begin{array}[]{ll}\dot{\chi}(t)&=\mathcal{A}\chi(t)+\mathcal{B}_{1}q(t)+\mathcal{B}_{2}\eta(t)\\\
r(t)&=\mathcal{C}_{1}\chi(t)+\mathcal{D}_{11}q(t)+\mathcal{D}_{12}\eta(t)\\\
z(t)&=\mathcal{C}_{2}\chi(t)+\mathcal{D}_{21}q(t)+\mathcal{D}_{22}\eta(t)\\\
\eta(t)&=\begin{bmatrix}y^{T}(t)&\hat{y}^{T}(t)\end{bmatrix}^{T},~{}\chi(0)=0\end{array}\right.$
(33)

where the matrices of the Linear Parameter Varying (LPV) extended system (33)
have the appropriate dimensions.

$\begin{array}[]{l}\mathcal{A}(s)=\begin{bmatrix}A+B\left[I-\Lambda(s)D\right]^{-1}\Lambda(s)C&0\\\
0&A_{\xi}\end{bmatrix}\\\
\mathcal{B}_{1}(s)=\begin{bmatrix}B\left[I-\Lambda(s)D\right]^{-1}\Lambda(s)\tilde{D}+\tilde{B}&B\\\
B_{\delta 2}&0\\\ 0&B_{\epsilon
2}\end{bmatrix},~{}\mathcal{B}_{2}=\begin{bmatrix}0\\\ B_{\xi}\end{bmatrix}\\\
\mathcal{C}_{1}=\begin{bmatrix}0\\\
C_{\xi}\end{bmatrix},~{}\mathcal{C}_{2}(s)=\begin{bmatrix}C+D\left[I-\Lambda(s)D\right]^{-1}\Lambda(s)C&0\\\
0&0\end{bmatrix}\\\ \mathcal{D}_{11}=\begin{bmatrix}0\\\ D_{\xi
2}\end{bmatrix},~{}\mathcal{D}_{12}=\begin{bmatrix}0\\\ D_{\xi
2}\end{bmatrix},~{}\mathcal{D}_{22}=0\\\
~{}\mathcal{D}_{21}(s)=\begin{bmatrix}\tilde{D}+D\left[I-\Lambda(s)D\right]^{-1}\Lambda(s)\tilde{D}&D\left[I-\Lambda(s)D\right]^{-1}\\\
0&0\end{bmatrix},\\\ \end{array}$ (34)

This extended system is obtained using the dynamics error system of (11).
According to (Pfifer and Seiler, 2015), the approximation error $\mu(t)$
expression in (19), and the extended virtual filter (32). The new system (33)
is equivalent to Fig. 2, where the system has an input of the state of the
uncertain system $\eta(t)$ and its output is the output of the error $z(t)$
with the uncertainties $\Delta$ and $\epsilon$, see Fig. 2. In the other hand,
a dissipation inequality can be formulated to upper bound the worst-case
$\mathcal{L}_{2}$ gain of $\tilde{F}(\tilde{P},\Delta,\epsilon)$ of Fig. 2
utilizing the extended system (33) and the time domain IQCs (29).

Figure 2:  Graphical interpretation of the IQCs with the equivalent system.

###  3.1  Relative Integral Square Error (RISE)

The following theorem derives the conditions for the worst case Relative
Integral of the Square of the Error (RISE) metric. In control theory, RISE is
used to assess the performance of a system (Atkins et al., 2013). RISE
provides valuable information about the behavior of the uncertain closed-loop
system with an NN controller compared to the closed-loop reference model.

######  Theorem 1 (Worst-case of RISE)

Let $\Delta$ and $\epsilon$ satisfy $IQC(\Xi,M)$ and assume that the extended
system $\Sigma_{ex}$ of (33) is well-posed. If there exists a scalar
$\lambda>0$ and a continuously differentiable matrix $P=P^{T}$ with
appropriate dimensions such that $P\geq 0$ and for all
$(s,\dot{s})\in\Omega\times\dot{\mathcal{S}}$ (Omitting $s$ and $\dot{s}$ for
a shorter notation):

$\displaystyle\begin{array}[]{c}\begin{bmatrix}P\mathcal{A}+\mathcal{A}P^{T}+\partial
P&P\mathcal{B}_{1}&P\mathcal{B}_{2}\\\ \mathcal{B}_{1}^{T}P&0&0\\\
\mathcal{B}_{2}^{T}P&0&-I\end{bmatrix}\\\
+\dfrac{1}{\gamma^{2}}\begin{bmatrix}\mathcal{C}_{2}^{T}\\\
\mathcal{D}_{21}^{T}\\\
\mathcal{D}_{22}^{T}\end{bmatrix}\begin{bmatrix}\mathcal{C}_{2}&\mathcal{D}_{21}&\mathcal{D}_{22}\end{bmatrix}\\\
+\lambda\begin{bmatrix}\mathcal{C}_{1}^{T}\\\ \mathcal{D}_{11}^{T}\\\
\mathcal{D}_{12}^{T}\end{bmatrix}M\begin{bmatrix}\mathcal{C}_{1}&\mathcal{D}_{11}&\mathcal{D}_{12}\end{bmatrix}<0\end{array}$
(38)

Then,

  * •

the Relative Integral Square Error (RISE) satisfies:

$\sup_{s\in\Omega.}\sup_{0\neq\eta\in\mathcal{L}_{2},\zeta(0)=0}\left[\dfrac{||y-\hat{y}||_{2}}{\sqrt{||y||_{2}^{2}+||\hat{y}||_{2}^{2}}}\right]\leq\gamma$
(39)

and

$||y-\hat{y}||_{2}<\
\frac{2\gamma}{1-\gamma}\left|\left|\hat{y}\left(t\right)\right|\right|_{2}$
(40)

##### Proof:

Consider the different signals $(\chi,z,p,\eta,r)$ of the error dynamics
system in with the input $\eta\in\mathcal{L}_{2}$ and solution trajectory
$\chi$ with zero initial conditions $\chi(0)=0$. Assuming the well-posedness
of (33) means that all signals of this system are well-defined. By assumption,
the uncertainties $\Delta$ and $\epsilon$ satisfy the IQC defined by
$(\Xi,M)$, and hence the quantity $r$ must satisfy the IQC time domain (29)
for any $T>0$. Also, using the state space representation of the virtual
filter in (16), the quantity $r^{T}Mr$ is expressed as follow:

$rMr^{T}=\begin{bmatrix}\underline{\xi}\\\ p\\\
q\end{bmatrix}^{T}\begin{bmatrix}C_{\xi}^{T}MC_{\xi}&C_{\xi}^{T}MD_{\xi}\\\
D_{\xi}^{T}MC_{\xi}&D_{\xi}^{T}MD_{\xi}\end{bmatrix}\begin{bmatrix}\underline{\xi}\\\
p\\\ q\end{bmatrix}$ (41)

The following proof is based on defining a parameter-dependent storage
function or Lyapunov function
$V:\mathbb{R}^{\chi}\times\Omega\rightarrow\mathbb{R}_{+}$ given by
$V(\chi,s)=\chi^{T}P(s)\chi$. Clearly, the condition (38) is a strict
inequality which implies that there exists very small quantity $\varepsilon>0$
such that the following perturbed matrix inequality holds for all
$s\in\Omega$:

$\displaystyle\begin{array}[]{l}\begin{bmatrix}P\mathcal{A}^{T}+\mathcal{A}P^{T}+\partial
P&P\mathcal{B}_{1}&P\mathcal{B}_{2}\\\ \mathcal{B}_{1}^{T}P&0&0\\\
\mathcal{B}_{2}^{T}P&0&-(1-\varepsilon)I\end{bmatrix}\\\
+\dfrac{1}{\gamma^{2}}\begin{bmatrix}\mathcal{C}_{2}^{T}\\\
\mathcal{D}_{21}^{T}\\\
\mathcal{D}_{22}^{T}\end{bmatrix}\begin{bmatrix}\mathcal{C}_{2}&\mathcal{D}_{21}&\mathcal{D}_{22}\end{bmatrix}\\\
+\lambda\begin{bmatrix}\mathcal{C}_{1}^{T}\\\ \mathcal{D}_{11}^{T}\\\
\mathcal{D}_{12}^{T}\end{bmatrix}M\begin{bmatrix}\mathcal{C}_{1}&\mathcal{D}_{11}&\mathcal{D}_{12}\end{bmatrix}\leq
0\end{array}$ (45)

Taking into account (41) and multiplying left and right (45) by
$[\chi^{T},q^{T},\eta^{T}]$ and its transpose $[\chi^{T},q^{T},\eta^{T}]^{T}$
where we can see that $V$ satisfies the following dissipation inequality:

$\displaystyle\nabla_{\chi}V\dot{\chi}+\nabla_{s}V\dot{s}\leq-\dfrac{1}{\gamma^{2}}z^{T}z+(1-\varepsilon)\eta^{T}\eta-\lambda
rMr^{T}$ (46)

This yields that the dissipation inequality (46) is equivalent to the
perturbed linear matrix inequality (38). Clearly, integrating the dissipation
inequality gotten in (46) from lower bound time $t=0$ to upper bound time
$t=T$ with the initial condition satisfies $\chi(0)=0$ to yield:

$\displaystyle V(\chi(T))+\lambda\int\limits_{0}^{T}\ \
r(\tau)^{T}\mathcal{M}r(\tau)\,\mathrm{\textit{d}\tau}\leq(1-\varepsilon)$
$\displaystyle\times\int\limits_{0}^{T}\ \
\eta(\tau)\eta(\tau)^{T}\,\mathrm{\textit{d}\tau}-\dfrac{1}{\gamma^{2}}\int\limits_{0}^{T}\
\ z(\tau)z(\tau)^{T}\,\mathrm{\textit{d}\tau}$ (47)

It follows from the positivity condition of IQCs as in (29), $\lambda\geq 0$,
and also due to the non-negativity of the storage function $V$ leads to:

$\dfrac{1}{\gamma^{2}}\ \int\limits_{0}^{T}\ \
z(\tau)z(\tau)^{T}\,\mathrm{\textit{d}\tau}\leq(1-\varepsilon)\int\limits_{0}^{T}\
\ \eta(\tau)\eta(\tau)^{T}\,\mathrm{\textit{d}\tau}$ (48)

As we assume that $\varepsilon>0$ , thus

$\displaystyle\lim_{T\rightarrow+\infty}\left[\dfrac{1}{\gamma^{2}}\
\int\limits_{0}^{T}\ \
z(\tau)z(\tau)^{T}\,\mathrm{\textit{d}\tau}<\int\limits_{0}^{T}\ \
\eta(\tau)\eta(\tau)^{T}\,\mathrm{\textit{d}\tau}\right]$ (49)

Knowing that

$\begin{array}[]{l}\left|\left|z(t)\right|\right|_{2}=||y(t)-\hat{y}(t)||_{2}\\\
\left|\left|\eta(t)\right|\right|_{2}=\sqrt{||y(t)||_{2}^{2}+||\hat{y}(t)||_{2}^{2}}\end{array}$
(50)

this allows us to deduce that:

$\sup_{s\in\Omega.}\sup_{0\neq
y,\hat{y}\in\mathcal{L}_{2},\zeta(0)=0}\dfrac{||y-\hat{y}||_{2}}{\sqrt{||y||_{2}^{2}+||\hat{y}||_{2}^{2}}}\leq\gamma$
(51)

Additionally, the next part of the proof will be based on the following
inequalities

$\displaystyle\begin{array}[]{l}\sqrt{||y||_{2}^{2}+||\hat{y}||_{2}^{2}}\leq||y||_{2}+||\hat{y}||_{2}\\\
||y||_{2}=||y-\hat{y}+\hat{y}||_{2}\leq||y-\hat{y}||_{2}+||\hat{y}||_{2}\end{array}$
(54)

Using (54) and (51), we get

$||y-\hat{y}||_{2}<\
\frac{2\gamma}{1-\gamma}\left|\left|\hat{y}\left(t\right)\right|\right|_{2}$
(55)

$\square$ The value $\gamma$ in Theorem 1 generally represents the largest
$\mathcal{L}_{2}$ gain overall allowable parameter trajectories. In this case,
it gives information about the worst case of the worst case of RISE between
the closed loop reference model and the closed-loop system with neural
controller under neural controller parameters $\Lambda$ and its training error
$\epsilon$. This information is valuable as it provides prior knowledge about
the behavior of the interconnected plant submitted to the uncertainty
$\Delta$. In other words, the result of this theorem can be used to improve
the robustness of the NN controller by making the dynamical error smaller
which will guarantee a better performance of the NN controller when it is
implemented. This worst-case measure is also an indication of the appearance
likelihood of any undesirable behavior in the closed-loop system when adopting
the designed NN controller.

On the other hand, the equation (40) holds substantial significance as it
delineates the correlation between the expected error in a system with neural
networks (NN) and the reference model, with the upper limit set by term
$\frac{2\gamma}{1-\gamma}$ multiplied by the reference model norm. As a result
of this crucial relationship, it is possible to understand how the NN operates
and predict the worst-case scenario for the system. As a result of
understanding this equation, we can improve risk management as well as fine-
tune the system so that deviations from the reference model are minimized.

###  3.2 Supreme Square Error (SSE)

The Linear Parameter Varying (LPV) model approximation of the robust tracking
problem developed above can be further analysed for more conservative
conditions. Indeed, LPV robust performance analysis of (33) can be extended to
cover less conservative performance measurements than the Relative Integral
Square Error (RISE). The following theorem is to derive the condition for
worst case of Supreme Square Error (SSE) $\sigma$ that less conservative
conditions as $||y-\hat{y}||_{2}\geq||y-\hat{y}||_{\infty}$. More precisely
than RISE, SSE provides significant information about the dynamical behavior
of the interconnected system with a NNs controller compared to the closed-loop
reference model. Based on the definition 1, SSE is a good estimation for the
maximum error between the two systems’ outputs. If the user can find the bound
set of the SSE which leads to increase the level of the truest for the user to
adopt these schemes.

######  Theorem 2 (Worst-case of SSE)

Let $\Delta$ and $\epsilon$ satisfy $IQC(\Psi,M)$ and the extended system
$\Sigma_{ex}$ of (33) is well-posed, if there exists a scalar $\lambda>0$ and
a continuously differentiable matrix $P=P^{T}$ with appropriate dimension such
that $P\geq 0$ and for all
$(s,\dot{s})\in\Omega\times\dot{\mathcal{S}}$(Omitting $s$ and $\dot{s}$ for a
shorter notation):

$\displaystyle\begin{array}[]{c}\begin{bmatrix}P\mathcal{A}+\mathcal{A}P^{T}+\partial
P&P\mathcal{B}_{1}&P\mathcal{B}_{2}\\\ \mathcal{B}_{1}^{T}P&0&0\\\
\mathcal{B}_{2}^{T}P&0&-I\end{bmatrix}\\\
+\lambda\begin{bmatrix}\mathcal{C}_{1}^{T}\\\ \mathcal{D}_{11}^{T}\\\
\mathcal{D}_{12}^{T}\end{bmatrix}M\begin{bmatrix}\mathcal{C}_{1}&\mathcal{D}_{11}&\mathcal{D}_{12}\end{bmatrix}<0\end{array}$
(58)

and

$\begin{bmatrix}P&\mathcal{C}_{2}^{T}\\\
\mathcal{C}_{2}&\sigma^{2}I\end{bmatrix}>0$ (59)

Then, the Supreme Square Error (SSE) satisfies:

$\sup_{s\in\Omega.}\sup_{0\neq\eta\in\mathcal{L}_{2},\zeta(0)=0}\left[\dfrac{||y-\hat{y}||_{\infty}}{\sqrt{||y||_{2}^{2}+||\hat{y}||_{2}^{2}}}\right]<\sigma$
(60)

##### Proof:

Having the signals $(\chi,z,\eta,q,r)$ generated for the extended system given
in 33 where we consider the input $\eta\in\mathcal{L}_{2}$ and parameter
trajectory $s\in\Omega$ assuming zero initial conditions. Well-posedness
assumption of the extended LPV system implies that all signals
$(\chi,z,\eta,q,r)$ are well-defined. By assumption as well, the uncertainties
$\Delta$ and $\epsilon$ satisfy the IQCs defined by $(\Xi,M)$ and hence for
any $T>0$, the time domain expression of the IQCs in (18) must be satisfied by
the signal $r$. Moreover, defining a storage function
$V:\mathbb{R}^{n_{\chi}}×\Omega\rightarrow\mathbb{R}_{+}$ by
$V(\chi,s)=\chi^{T}P(s)\chi$. Also, the state space representation of the
virtual filter (16) gives (41). Evaluating of the perturbed LMI (58) at
$(s(t),\dot{s}(t))$ after multiplying its left/right by
$[\chi^{T},q^{T},\eta^{T}]$ and $[\chi^{T},q^{T},x^{T}]^{T}$, respectively,
yields the following dissipation inequality:

$\nabla_{\chi}V\dot{\chi}+\nabla_{s}V\dot{s}\leq(1-\varepsilon)\eta\eta^{T}-\lambda
rMr^{T}$ (61)

This means that the dissipation inequality (61) is equivalent to the perturbed
linear matrix inequality (58). Then, integrating along the state/parameter
trajectory the obtained dissipation inequality (61) from lower bound $t=0$ to
upper bound $t=T$. Recalling that $\chi(0)=0$, thus:

$\begin{array}[]{ll}V(\chi(T),s(T))\leq&(1-\varepsilon)\int\limits_{0}^{T}\ \
\eta(\tau)\eta(\tau)^{T}\,\mathrm{\textit{d}\tau}\\\
&-\lambda\int\limits_{0}^{T}\ \
r(\tau)Mr(\tau)^{T}\,\mathrm{\textit{d}\tau}\end{array}$ (62)

It follows from the IQCs condition (18), that

$V(\chi(T),s(T))\leq(1-\varepsilon)\int\limits_{0}^{T}\ \
\eta(\tau)\eta(\tau)^{T}\,\mathrm{\textit{d}\tau}$ (63)

Now, applying Schur complement on (59) and multiplying left and right by
$\chi,~{}z$ and its transpose $\chi^{T}$, respectively, results:

$\dfrac{1}{\sigma^{2}}z(t)z^{T}(t)\leq\chi^{T}(t)P(s)\chi(t)~{}~{}\forall t>0$
(64)

Evaluating (64) at $t=T$, and applying (63) yields:

$\dfrac{1}{\sigma^{2}}z(T)z^{T}(T)\leq(1-\varepsilon)\int\limits_{0}^{T}\ \
\eta(\tau)\eta(\tau)^{T}\,\mathrm{\textit{d}\tau}$ (65)

Knowing the $\varepsilon>0$, therefore

$\dfrac{1}{\sigma^{2}}z(T)z^{T}(T)\leq\int\limits_{0}^{T}\ \
\eta(\tau)\eta(\tau)^{T}\,\mathrm{\textit{d}\tau}$ (66)

Considering the supremum over $T$ to demonstrate that
$||z||_{\infty}\leq\sigma||\eta||_{2}$. Knowing that this result holds for any
given input $\eta\in\mathcal{L}_{2}$, admissible trajectory solution
$s\in\Omega$, and uncertainty $\Delta,\epsilon\in IQC(\Xi,M)$. Thus,

$\sup_{s\in\Omega.}\sup_{0\neq\eta\in\mathcal{L}_{2},\zeta(0)=0}\dfrac{\left|\left|z(t)\right|\right|_{\infty}}{\left|\left|\eta(t)\right|\right|_{2}}<\sigma$
(67)

Having $\left|\left|z(t)\right|\right|_{\infty}=||y(t)-\hat{y}(t)||_{\infty}$
and
$\left|\left|\eta(t)\right|\right|_{2}=\sqrt{||y(t)||_{2}^{2}+||\hat{y}(t)||_{2}^{2}}$,
this implies

$\dfrac{\left|\left|z(t)\right|\right|_{\infty}}{\left|\left|\eta(t)\right|\right|_{2}}=\dfrac{||y-\hat{y}||_{\infty}}{\sqrt{||y||_{2}^{2}+||\hat{y}||_{2}^{2}}}$
(68)

$\square$ Theorem 2 is very important as it provides the maximum deviation
expected of the closed-loop system with a NN controller from the closed-loop
reference model. Having such information will give the user confidence to
implement the NN controller by making the dynamical error small to insure the
best performance possible. This worst case of SSE is a good indicator of the
occurrence likelihood of any undesirable behavior in the closed-loop system
when adopting the designed NN controller.

##  4  Validation and Numerical Simulation

In order to validate the theoretical results we have proved in the previous
section, we will conduct tests using two different problems. Our first problem
involves controlling a robot arm with only a single output, called the Single-
Link Robot Arm Control Problem. A simpler and more focused scenario will allow
us to evaluate the effectiveness of our approach.

Further, we will also test our approach on the more complex problem of Deep
Guidance and Control of Apollo Lander, which involves multiple outputs. During
the descent and landing phase of the Apollo Lander, precise control and
guidance is required. We can further evaluate the robustness and adaptability
of our approach by examining this more complex scenario with multiple outputs.

We notice here that the property of the activation function of the controller
in quantity $\Lambda(s)$ allows us to express the set of a bounded convex
domain $\mathcal{G}_{m}^{n}$ of the matrix $\Lambda(s)$ as follows:

$\mathcal{G}_{m}^{n}=\left\\{a_{\text{ij}}\leq\frac{\partial\pi_{i}}{\partial
y_{j}}\leq b_{\text{ij}},\ i=1,\cdots,n,\ j=1,\cdots,m\right\\}$ and the set
of vertices is defined by

$\mathcal{V}_{\mathcal{G}_{m}^{n}}=\left\\{\mathcal{V=}\left[\mathcal{V}_{11},\cdots,\mathcal{V}_{1m},\cdots,\mathcal{V}_{\text{nm}}\right]:\mathcal{V}_{\text{ij}}\in\left\\{a_{\text{ij}},\
b_{\text{ij}}\right\\}\right\\}$ (69)

Also, IQCs provides a tool to study the stability of feedback connections of
stable LTI systems and perturbation blocks as shown in Fig. 2. IQCs
corresponding to many perturbation blocks is already available in the
literature. Furthermore, a MATLAB toolbox for IQC-based stability and
performance analysis exists and is called IQC-$\beta$. For different
perturbations, different multipliers, e.g. constant or frequency dependent,
are necessary to lead to non-conservative analysis and design results. Within
IQC-$\beta$ toolbox the different multipliers are chosen automatically for
user’s convenience. As the dissipation inequality (46) is equivalent to the
perturbed linear matrix inequality (38), the user, in IQC-$\beta$ toolbox, has
only to choose the type of IQC, which is more physically related, e.g. LTI
uncertain parameter (Jonsson et al., 2004).

Since the modelling open source toolbox IQC-$\beta$ we are using in this work
does not support LPV formulation, the source of parameters variation,
represented by the matrix $\Lambda(s)$, will therefore be modelled using the
vertices defined in (69). More information can be found in (Boyd and
Vandenberghe, 2001).

###  4.1  Single-Link Robot Arm Control Problem

####  4.1.1  Step 1: Problem Setting

The simulation validation that we conduct in this section uses the worst-case
of RISE metric resulting from Theorem 1 since the open source tool we adopt
for simulation calculations can only provide $\mathcal{L}_{2}$ to
$\mathcal{L}_{2}$ gain.

Consider the nonlinear single-link robot arm example, shown in Fig. 3, with
parameters $m=0.15kg$, length $l=0.5m,$ and friction coefficient
$\mu=0.5Nms/rad$. The state space representation of the arm’s motion is given
by:

$\left\\{\begin{array}[]{ll}\dot{\omega}(t)&=-10\sin(\theta(t))-2\omega(t)+\tau(t)\\\
\dot{\theta}(t)&=\omega(t)\end{array}\right.$ (70)

Figure 3:  Single-link robot arm.

####  4.1.2  Step 2: Linear Reference Model

We start by training a neural controller with the aim that the arm with this
neural controller tracks the closed-loop reference model given by:

$\left\\{\begin{array}[]{ll}\dot{\omega}_{r}(t)&=-9\theta_{r}(t)-6\omega_{r}(t)+9r(t)\\\
\dot{\theta}_{r}(t)&=\omega_{r}(t)\end{array}\right.$ (71)

This means that the dynamic behavior of the system (70) in the closed-loop
system after the training of the neural controller is close to the desired
closed loop dynamics (71). The ideal case is when the system with a NN
controller (70) matches the closed loop reference model (71) at all time $t$.
However, this is nearly impossible in practice. The objective of this example
is to determine the worst relative squared difference, $\gamma$, between the
robot arm with neural control system closing the loop, while subject to its
nonlinear functions, and the ideal closed loop reference model, where

$\sup_{s\in\Omega.}\sup_{0\neq[\theta(t)~{}\theta_{r}(t)]^{T}\in\mathcal{L}_{2},\zeta(0)=0}\frac{\left|\left|\theta\left(t\right)-\theta_{r}(t)\right|\right|_{2}}{\left|\left|[\theta^{T}(t)~{}\theta_{r}^{T}(t)]^{T}\right|\right|_{2}}<\gamma$
(72)

To this end, let assume the static non-linearity
$\Delta(\theta)=\theta-\sin(\theta)$ which is slope-restricted and sector
bounded. We rearrange (70) then into the form:

$\left\\{\begin{array}[]{ll}\dot{\omega}(t)&=-9\theta(t)-2\omega(t)+\tau(t)+10q(t)\\\
\dot{\theta}(t)&=\omega(t),~{}p(t)=\theta(t)\\\
\Delta(\theta)&=q(t)=\theta-\sin(\theta)\end{array}\right.$ (73)

####  4.1.3  Step 3: Error Dynamical System

Now, we generate the error dynamical system using the closed loop reference
model (71) and the new mathematical form of the arm’s motion (73). Let
$\delta\omega=\omega-\omega_{r}$, $\delta\theta=\theta-\theta_{r}$ and
$\delta\tau=\tau-\tau_{r}$ then

$\displaystyle\dfrac{d}{dt}\begin{bmatrix}\delta\omega\\\
\delta\theta\end{bmatrix}=\begin{bmatrix}-6&-9\\\
1&0\end{bmatrix}\begin{bmatrix}\delta\omega\\\
\delta\theta\end{bmatrix}+\begin{bmatrix}10\\\
0\end{bmatrix}q(t)+\begin{bmatrix}1\\\ 0\end{bmatrix}\delta\tau(t)$ (74)

where $\tau_{r}(t)=\pi^{*}\left(\theta,t\right)$ and

$\pi^{*}\left(\theta,t\right)=10\theta_{r}(t)+4\omega_{r}(t)-9r(t)-10\sin(\theta_{r}(t))$
(75)

The Neural controller architecture adopted in this example is given in Fig. 4.
The output of the neural controller is given by :

$\pi\left(\theta,t\right)=w_{4}\tanh\left(w_{1}\theta\left(t\right)+w_{2}r\left(t\right)+w_{3}\pi\left(\theta\left(t\right)\right)+b_{1}\right)+b_{2}$
(76)

where $\tau(t)=\pi\left(\theta,t\right)$.

Figure 4:  Neural controller $\pi^{*}$ architecture.

Differential Mean Value Theorem (Zemouche et al., 2005) introduced in lemma 1,
can convert the dynamical error system to an LPV system without nonlinearities
of the neural controller. According to lemma 2, the approximation error
$\delta\tau(t)=\mu(t)$ given in (11) can be expressed as :

$\delta\tau(t)=\mu(t)=\Lambda(s)C\zeta(t)+\Lambda(s)D_{p}\Delta(y(t))+\epsilon(\hat{y}(t))$
(77)

and the vectors $s_{1},\cdots,s_{n}\in\mathbf{Co}(y(t),\hat{y}(t))$,
$s_{i}\neq y_{1},~{}s_{i}\neq y_{2}$ , $i=1,\cdots,n$ such that

$\displaystyle\Lambda(s)=\sum\limits_{i=1}^{m}\sum\limits_{j=1}^{n}e_{m}(i)e_{n}^{T}(j)\dfrac{\partial\pi_{i}}{\partial
y_{j}}(s_{i})$

where

$\epsilon(\hat{y}(t))=\pi(\hat{y}(t))-\pi^{*}(\hat{y}(t),t)$

In this case, we have $C=\begin{bmatrix}0&1\end{bmatrix}$,
$\zeta(t)=\begin{bmatrix}\delta\omega(t)&\delta\theta(t)\end{bmatrix}^{T}$,
$D_{p}=0$, $\hat{y}(t)=\theta_{r}$ and $y(t)=\theta(t)$. Therefore,

$\delta\tau(t)=\Lambda(s)\delta\theta(t)+\epsilon(\theta_{r}(t))$ (78)

####  4.1.4  Step 4: Characterisation of the nonlinearities $\Delta$ and
$\epsilon$ using IQCs

To choose the type of IQCs of the uncertainties considered for this simulation
example, let us start with the static non-linearity
$\Delta(\theta)=\theta-\sin(\theta)$. It is clear from the drawn figure Fig. 5
that the uncertainty $\Delta(\theta)$ is slope-restricted and sector bounded.
Fig. 5 provides the IQCs characterisation of this non-linearity and plots the
graph $\Delta=f(\theta)$ where $\Delta(\theta(t))=\varphi(\theta,t)$ and
$\alpha\theta(t)\leq\varphi(\theta,t)\leq\beta\theta(t)$ for all
$\theta\in[-\dfrac{\pi}{2},~{}\dfrac{\pi}{2}]$. Simple calculation of the
slope gives $\alpha=0$ and $\beta=0.364$ representing the bounds of the
uncertainty $\Delta$.

Figure 5:  IQCs characterisation of the non-linearity $\Delta$.

Similarly, we have the nonlinearity
$\epsilon(\theta_{r})=\pi(\theta_{r})-\pi^{*}(\theta_{r})$, where $\pi^{*}$
and $\pi$ are given in (75) and (76), respectively. The graph of
$\epsilon=\hat{f}(\theta_{r})$ is shown in Fig 6 where
$\epsilon(\theta_{r}(t))=\varrho(\theta_{r},t)$ and
$|\varrho(\theta_{r},t)|\leq 0.2$ for all
$\theta_{r}\in[-\dfrac{\pi}{2},~{}\dfrac{\pi}{2}]$. The IQCs characterisation
of this nonlinearity is thus conducted as we show in Fig 6 that it is bounded
by 0.2 and -0.2.

Figure 6:  IQCs characterisation of the non-linearity $\epsilon$ .

####  4.1.5  Step 5: IQCs based analysis results

Now, we would like to compute an upper bound on the $\mathcal{L}_{2}$-gain
from the output error $\theta-\theta_{r}$ to the outputs vector
$[\theta,~{}\theta_{r}]$. We use IQC-$\beta$ toolbox to express and formulate
the $\mathcal{L}_{2}$ gain estimation problem. Fig. 7 explains how to adapt
our case to MATLAB IQCs-$\beta$ Toolbox. We only need the two linear systems,
the closed loop reference model and the dynamical error system, since the
output of the closed-loop system with the NN controller is the sum of the
outputs of the those two systems. The bounds of the uncertainty $\Delta$ and
the nonlinearity $\epsilon$ derived in the previous sub-section are also
provided as inputs to the IQCs-$\beta$ Toolbox.

Figure 7:  Adaption of IQCs-$\beta$ MATLAB tool for the problem.

According to the results given by IQCs-$\beta$ for the model in Fig.8
corresponding to the original Single-Link Robot Arm Control Problem:

$\sup_{s\in\Omega.}\sup_{0\neq[\theta(t)~{}\theta_{r}(t)]^{T}\in\mathcal{L}_{2},\zeta(0)=0}\frac{\left|\left|\theta\left(t\right)-\theta_{r}(t)\right|\right|_{2}}{\left|\left|[\theta^{T}(t)~{}\theta_{r}^{T}(t)]^{T}\right|\right|_{2}}<\gamma=0.05669$
(79)

$\gamma$ is the worst case of the mean integrated squared error between the
closed loop system with the trained neural controller and the closed loop
reference system and is less than $5.669\%$. According to (Atkins et al.,
2013), the time Integral of the Square of the Error (ISE) in a control system
is used as a measure of the performance of the system.

The value of $\gamma$ represents the worst case of the relative error between
the closed loop reference model, used for collecting data to train the NN
controller, and the uncertain system with the designed neural controller.

Figure 8:  Case Study A: (a) Plot for the closed-loop system with the NN
within the region of expectation (b) Plot of the output of reference model
response $\theta_{r}$ indicated in blue and the one of the closed-loop
uncertain robotic arm system with neural controller $\theta$ indicated in red.
(c) The controllers $\pi$ and $\pi^{*}$ outputs. (d) The relative squared
difference between the robot arm neural control system is subject to
uncertainty in nonlinear functions and the reference model.  Figure 9: Case
Study B: (a) Plot for the closed-loop system with the NN within the region of
expectation (b) Plot of the output of reference model response $\theta_{r}$
indicated in blue and the one of the closed-loop uncertain robotic arm system
with neural controller $\theta$ indicated in red. (c) The controllers $\pi$
and $\pi^{*}$ outputs. (d) Relative squared difference between the robot arm
neural control system subject to uncertainty in nonlinear functions and the
reference model.

This worst-case result can be used to improve the robustness of the NN
controller by making the dynamical error smaller which in turn will guarantee
a better performance of the NN controller when it is implemented. This worst-
case measure is also an indication of the appearance likelihood of any
undesirable behavior in the closed loop system when adopting the designed NN
controller.

To illustrate further the results of our controlled Single-Link Robot Arm with
our designed neural controller, we propose two study cases: A and B. The
reference of the first case study is periodically applied to the uncertain
system with a neural controller and the reference closed-loop model. Using the
equation given in (40), we can plot the expected region by utilise the
information about the worst-case of the error and the known linear reference
model as follow:

$\left|\left|\theta\left(t\right)-\theta_{r}(t)\right|\right|_{2}<\
\frac{2\gamma}{1-\gamma}\left|\left|\theta_{r}(t)\right|\right|_{2}=0.1202\left|\left|\theta_{r}(t)\right|\right|_{2}$
(80)

Fig. 8 shows in sub-figure (a) the response of the system with the trained NN
controller and clearly, the output remains inside the expected bound region.
In sub-figure (b) of Fig. 8, the two responses (i.e. closed loop reference
model and the uncertain plant with NN controller) for this first case study
are shown. We can see from this figure that the uncertain plant with NN
controller keeps close to the response of the closed-loop reference model. The
neural controller output is also very close to the output of the reference
controller used in the reference closed-loop model as shown in sub-figure (c)
of Fig. 8. Moreover, we can see in sub-figure (d) of Fig. 8 that the value of
$\frac{\left|\left|\theta\left(t\right)-\theta_{r}(t)\right|\right|_{2}}{\left|\left|[\theta^{T}(t)~{}\theta_{r}^{T}(t)]^{T}\right|\right|_{2}}$
reaches the maximum of $0.045$ approximately which is less than the worst-case
$\gamma=0.05669$. The same thing said above for the sub-figures of Fig. 8 can
be said regarding the non-periodic case study B illustrated in the sub-figures
of Fig 9 where the value
$\frac{\left|\left|\theta\left(t\right)-\theta_{r}(t)\right|\right|_{2}}{\left|\left|[\theta^{T}(t)~{}\theta_{r}^{T}(t)]^{T}\right|\right|_{2}}$
has a maximum pic of around $0.051$ which is less than the worst case $\gamma$
as well and the response of the system with NN controller is always inside the
region of performance expectation.

###  4.2 Deep Guidance and Control of Apollo Lander

In the classical Guidance and Control of Apollo Lander described in (Amaral,
2019) , GC system generates trajectory and thrust commands to manoeuvre the
lander for a given initial position, velocity, and desired landing target. The
control part is compensated for by the Martian atmosphere and converted into
inertial coordinates before the position commander uses the information
(Millour et al., 2015).

####  4.2.1  Step 1: Problem Setting

We perform the test of the trained network in closed-loop using the parameters
listed in Table 1 with the aim of validating its performance. In terms of
applying aerodynamic forces and torques which are external disturbances on the
lander, coefficient of drag $C_{D}$ is $2.0$, and it is assumed to be
constant, whereas surface-relative head wind is in the range of $20m/s$.
Furthermore, the air density $\rho$ is assumed to be $0.023kg/m^{3}$
throughout. Under this scenario, the worst-case aerodynamic resistance during
the descent phase is less than 2% of the vehicle’s maximum thrust. Due to the
close-loop control action, aerodynamic forces and torques are automatically
compensated.

Table 1: Simulation Parameters  Parameter | value
---|---
Initial Position $p_{0}$ $[m]$ | $[-5632.2,709.25,6190.5]^{T}$
Initial Velocity $v_{0}$ $[m/s]$ | $[206.48,-26.006,-103.51]^{T}$
Initial Acceleration $a_{0}$ $[m/s^{2}]$ | $[-2.1124,0.24947,-2.6404]^{T}$
Intial Mass $m_{0}$ $[Kg]$ | 600.3
The capsule reference area $S$ | 5.137426149499100
Mars Gravity $g$ $[m/s^{2}]$ | $[0,0,3.725258]^{T}$
Engine Maximum Thrust $T$ $[N]$ | 3600
Maximum Throttle Level $T_{max}$ | 100%
Minimum Throttle Level $T_{min}$ | 30%

##### Atmosphere

As part of the Mars environment module, the horizontal wind is included in the
atmospheric model. Based on lookup tables, an atmospheric model specifies the
properties of the Martin atmosphere at a specific location. The atmospheric
properties are obtained directly from the Mars Climate Database (MCD) v.5.2 at
the mission’s specified landing coordinates. MCD is derived from numerical
simulations of the Martian atmosphere using General Circulation Models (GCM)
and validated by observations of the Martian atmosphere. For more details see
(Millour et al., 2015).

##### Lander Dynamics

In closed-loop simulations, Newton’s second law is used to model a point
mass’s robotic lander motion, where the lander motion is expressed in the LENU
frame as follow:

$\frac{d\left[r\right]_{LENU}}{dt}=\left[v\right]_{LENU},$ (81)

where ($\left[v\right]_{LENU}$) refer to the relative velocity of the robotic
lander. In in the LENU frame and the lander mass ($m$), The dynamics equation
is given based on the force resultant ($\left[F\right]_{LENU}$) expressed as:

$\frac{d\left[v\right]_{LENU}}{dt}=\left[a\right]_{LENU}=\frac{\left[F\right]_{LENU}}{m}.$
(82)

##### Aerodynamics

In the simulations, a lookup table is used as the aerodynamics model of the
Apollo capsule. In this module, only aerodynamic drag forces are considered
following equation:

$F_{D}=\frac{1}{2}\rho v^{2}C_{D}\left(\text{Ma}\right)S,$ (83)

where $\rho$ represents the atmosphere density, $C_{D}$ denotes the capsule
drag coefficient as a function of the Mach number (Ma), the capsule reference
area represented by the $S$, and $v$ is the capsule velocity relative to the
fluid velocity expressed in the landing site East-North-UP frame.

Figure 10: Plot of the 1537 trajectories (Amaral, 2019). The marked dots
representing the trigger locations. An optimal altitude to trigger is in the
range of 2 km to 3 km, and in some trajectories, overshooting the target is
the best behavior.

##### Actuators

In the actuator model, four thrusters engines are used to generate forces
concerning the center of mass of the vehicle. The main engine responsible on
the force $F_{Thr}$ is located along with the vehicle $z$-axis, while the
remaining are the sources of saturated forces $F_{Trq}$ are placed at a
distance $l$ from the mass center to generate pure torques along positive and
negative axes. To calculate the fuel consumption, the mass flow rate is
integrated:

${\dot{m}}_{fuel}={V_{e}}^{-1}\left(F_{Thr}+\sum F_{Trq}\right),$ (84)

where $V_{e}=g_{0}I_{sp}$.

##### Apollo Guidance

The algorithm is based on the idea of following the reference trajectory
backwards in time from the target point to the present moment. The time-to-go
($t_{go}$), represents the time left until the desired landing is reached, so
it tends to zero as time gets closer to the goal. In (Klumpp, 1974), the
reference trajectory satisfies a two-boundary problem with five degrees of
freedom. Accordingly, the reference trajectory must be a polynomial function
of fourth or higher order. Therefore, the reference position ($r_{d}$ given at
$t_{go}$ as a function of the target state as:

$r_{d}=\
r_{t}+v_{t}t_{go}+a_{t}\frac{{t_{go}}^{2}}{2}+j_{t}\frac{{t_{go}}^{3}}{6}+s_{t}\frac{{t_{go}}^{4}}{24},$
(85)

where $r_{t}$ refers to the target position, $v_{t}$ represents the target
velocity and $a_{t}$, $j_{t}$, $s_{t}$ denote the target acceleration, jerk
and snap, respectively.

According to (85), the work of (Huang et al., 2022) describes the acceleration
command ($a_{cmd}$) as a quadratic polynomial i.e.

$a_{\text{cmd}}=\ C_{0}+C_{1}t_{go}+{{C_{2}t}_{go}}^{2}.$ (86)

The coefficients of the polynomial are derived by solving a system of three
equations comprising (86). The two other equations related to the target
velocity and position equations are obtained by by integrating the target
acceleration equation:

$\left\\{\begin{array}[]{l}a_{t}=\ C_{0}+C_{1}t_{go}+{{C_{2}t}_{go}}^{2}\\\
v_{t}=C_{0}t_{go}+\frac{1}{2}C_{1}{t_{go}}^{2}+\frac{1}{3}C_{2}{t_{go}}^{3}+v\\\
r_{t}=vt_{go}+{\frac{1}{2}C}_{0}{t_{go}}^{2}+\frac{1}{6}C_{1}{t_{go}}^{3}+\frac{1}{12}C_{2}{t_{go}}^{4}+r\\\
\end{array}\right.$ (87)

where $a_{t}$ represents the reference acceleration vector, $r$ and $v$ denote
the vehicle current position and velocity, respectively. The following values
for the polynomial coefficients $C_{0}$, $C_{1}$ and $C_{2}$) are obtained by
solving the following matrix system:

$\begin{bmatrix}a_{t}\\\ v_{t}\\\
r_{t}\end{bmatrix}=\begin{bmatrix}1&t_{go}&t_{go}^{2}\\\
t_{go}&\frac{1}{2}t_{go}^{2}&\frac{1}{3}t_{go}^{3}\\\
\frac{1}{2}t_{go}^{2}&\frac{1}{6}{t_{go}}^{3}&\frac{1}{12}{t_{go}}^{4}\end{bmatrix}\begin{bmatrix}C_{0}\\\
C_{1}\\\ C_{2}\end{bmatrix}+\begin{bmatrix}0\\\ v\\\ v_{t_{go}}\end{bmatrix}.$
(88)

Therefore,

$\left\\{\begin{array}[]{l}C_{0}=\
a_{t}-6\left(\frac{v_{t}+v}{t_{go}}\right)+12\left(\frac{r_{t}-r}{{t_{go}}^{2}}\right)\\\
C_{1}=\
-6\left(\frac{a_{t}}{t_{go}}\right)+6\left(\frac{{5v}_{t}+3v}{{t_{go}}^{2}}\right)-48\left(\frac{r_{t}-r}{{t_{go}}^{3}}\right)\\\
C_{2}=\
6\left(\frac{a_{t}}{{t_{go}}^{2}}\right)-12\left(\frac{{2v}_{t}+v}{{t_{go}}^{3}}\right)+36\left(\frac{r_{t}-r}{{t_{go}}^{4}}\right)\end{array}\right.$
(89)

By defining the vertical component of the acceleration profile as a linear
function of $t_{go}$. Hence, solving (86) for $t_{go}$ two solutions are
obtained:

$\begin{array}[]{l}\mathbf{If}~{}\left|\left(a_{t}\right)_{z}\right|>0\\\
t_{go}=\frac{2\left(v_{t}\right)_{z}+\left(v\right)_{z}}{\left(a_{z}\right)_{t}}+\sqrt{\left[\frac{2\left(v_{t}\right)_{z}+\left(v\right)_{z}}{\left(a_{t}\right)_{z}}\right]^{2}+\frac{6\left[\left(\mathbf{r}\right)_{z}-\left(r_{t}\right)_{z}\right]}{\left(a_{t}\right)_{z}}},\\\
\mathbf{else}~{}\mathbf{if}\left(a_{t}\right)_{z}=0\\\
t_{go}=\frac{3\left[\left(r_{t}\right)_{z}-\left(r\right)_{z}\right]}{\left(v\right)_{z}+2\left(v_{t}\right)_{z}},\end{array}$
(90)

where $(.)_{z}$ denotes the vertical component of the variable $(.)$

A robustness analysis of feedback systems Apollo lander based on NN
controllers against potential uncertainties is presented here. To achieve
this, we need to follow the below pattern where we aiming at keeping the
output of the closed-loop system with a NN controller close to the output of
an ideal and trusted closed-loop reference model when its input changes within
a bounded set.

####  4.2.2  Step 2: Linear Reference Model

We start by training a neural controller with the aim that the lander system
with this neural controller tracks classical Apollo controller provided in
(Amaral, 2019). A linear reference model is what we want. To generate the
linear closed-loop reference model, we use the Five Tau rule, which states
that transient events die down after five Tau seconds. In other words, a
transition from a (steady) state to another (steady) state lasts only for five
Tau seconds. Let consider the following form of the reference model:

$\left[\begin{matrix}\frac{d\hat{r}}{dt}\\\ \frac{d\hat{v}}{dt}\\\
\end{matrix}\right]=\left[\begin{matrix}\hat{v}(t)\\\ \
a\hat{v}(t)+\hat{u}(t)\\\ \end{matrix}\right]$ (91)

where

$a=\left(\begin{matrix}a_{x}&0&0\\\ 0&a_{y}&0\\\ 0&0&a_{z}\\\
\end{matrix}\right)$

MATLAB is used to find the matrix $a$ using regression linear technique:

$a=\arg\min_{v}(\frac{1}{2}\frac{\rho}{m}v^{2}C_{D}\left(\mathrm{Ma}\right)S\
-av(t))$ (92)

As a result, we obtain the matrix $a$ as follows:

$a=\left(\begin{matrix}-0.0087&0&0\\\ 0&-0.0075&0\\\ 0&0&-0.0077\\\
\end{matrix}\right)$ (93)

The linear reference model is now complete, and we can move on to finding the
best controller which has dynamical behaviour too close to the nonlinear
classical one. MATLAB is used to implement the poles placement technique
(state feedback control) along with Five Tau rule. The linear closed-loop
reference model and the nonlinear classical Apollo are compared via simulation
for some scenarios to illustrate their close responses. Both systems have
closed responses, as shown in Figure 11 and Figure 12.

Figure 11:  Position: linear reference model vs. the nonlinear Apollo Model.
Figure 12:  Velocity: linear reference model vs. the nonlinear Apollo .

The neural controller in TABLE 2 was trained with the aim that the lander
dynamics with this neural controller tracks the closed-loop reference Apollo
model. This means that the dynamic behaviour of the system in the closed-loop
system after the training of the neural controller is close to the linear
reference closed loop dynamics. In an ideal world, the closed loop Apollo
model would be matched by the system with a NN controller at all times t. In
reality, this is almost impossible. In this example, we are calculating the
worst relative squared difference, $\gamma$, between the classical Apollo
system provided in (Amaral, 2019) with neural control system closing the loop,
while subject to its nonlinear functions, and the ideal closed loop reference
model. For this purpose, let assume the static non-linearity
$\Delta\left(v(t)\right)=\frac{1}{2}\frac{\rho}{m}v^{2}(t)C_{D}\left(\mathrm{Ma}\right)S\
-av(t)$ which is slope-restricted, and sector bounded. The Lander dynamics are
then rearranged into the following form:

$\left\\{\begin{matrix}\left(\begin{matrix}\frac{dr}{dt}\\\ \frac{dv}{dt}\\\
\end{matrix}\right)=\left[\begin{matrix}v\\\
av\left(t\right)+u\left(t\right)+\Delta\left(v\left(t\right)\right)\\\
\end{matrix}\right]\\\ \
\Delta\left(v(t)\right)=\frac{1}{2}\frac{\rho}{m}v^{2}(t)C_{D}\left(\mathrm{Ma}\right)S\
-av(t)\\\ y\left(t\right)=\left(\begin{matrix}r\left(t\right)\\\
v\left(t\right)\\\ \end{matrix}\right)\\\ \end{matrix}\right.$ (94)

Table 2: Neural Controller Architecture Layer | No. Neurons/Activation Function
---|---
Input layer | 6/ Linear $[r^{T}(t),v^{T}(t)]^{T}$
$1^{st}$ layer | 40/Tanh
$2^{nd}$ layer | 40/Sigmoid
$3^{rd}$layer | 40/Tanh
Output layer | 3/ Linear

####  4.2.3  Step 3: Error Dynamical System

By using the closed loop reference model and the new mathematical form of the
Lander dynamical motion, we can now construct the error dynamical system. Let
$\delta r=r-\hat{r}$, $\delta v=v-\hat{v}$ and $\delta u=u-\hat{u}$ then

$\left\\{\begin{matrix}\left(\begin{matrix}\frac{d\delta r}{dt}\\\
\frac{d\delta v}{dt}\\\ \end{matrix}\right)=\left[\begin{matrix}\delta v\\\
a\delta v\left(t\right)+\delta
u\left(t\right)+\Delta\left(v\left(t\right)\right)\\\ \end{matrix}\right]\\\ \
\Delta\left(v(t)\right)=\frac{1}{2}\frac{\rho}{m}v^{2}(t)C_{D}\left(\mathrm{Ma}\right)S\
-av(t)\\\ z\left(t\right)=\left(\begin{matrix}\delta r\left(t\right)\\\ \delta
v\left(t\right)\\\ \end{matrix}\right)\\\ \end{matrix}\right.$ (95)

According to lemma 2, the approximation error $\delta u=u-\hat{u}$ given in
the equation above can be expressed as:

$\delta u\left(t\right)=\Lambda\left(s\right)\left(\begin{matrix}\delta
r\left(t\right)\\\ \delta v\left(t\right)\\\
\end{matrix}\right)+\epsilon\left(\hat{v}\left(t\right)\right)$ (96)

where
$\epsilon\left(\hat{v}\left(t\right)\right)=u\left(\hat{v}\left(t\right)\right)-\hat{u}\
\left(\hat{v}\left(t\right)\right)$. Then, the dynamical error became:

$\displaystyle\left\\{\begin{array}[]{ll}\begin{matrix}\left[\begin{matrix}\frac{d\delta
r}{dt}(t)\\\ \frac{d\delta v}{dt}(t)\\\
\end{matrix}\right]&=\left[\begin{matrix}\delta v(t)\\\ a\delta
v\left(t\right)+\Delta\left(v\left(t\right)\right)\\\
\end{matrix}\right]+B\Lambda\left(s\right)\left[\begin{matrix}\delta
r\left(t\right)\\\ \delta v\left(t\right)\\\ \end{matrix}\right]\\\
&+B\epsilon\left(\hat{v}\left(t\right)\right)\\\ \
\Delta\left(v(t)\right)&=\frac{1}{2}\frac{\rho}{m}v^{2}(t)C_{D}\left(\mathrm{Ma}\right)S\
-av(t)\\\ z\left(t\right)&=\left(\begin{matrix}\delta r\left(t\right)\\\
\delta v\left(t\right)\\\ \end{matrix}\right),\ \
B=\left[\begin{matrix}\begin{matrix}\begin{matrix}0&0&0\\\ \end{matrix}\\\
\begin{matrix}0&0&0\\\ \end{matrix}\\\ \begin{matrix}0&0&0\\\ \end{matrix}\\\
\end{matrix}\\\ \begin{matrix}\begin{matrix}1&0&0\\\ \end{matrix}\\\
\begin{matrix}0&1&0\\\ \end{matrix}\\\ \begin{matrix}0&0&1\\\ \end{matrix}\\\
\end{matrix}\\\ \end{matrix}\right]\\\ \end{matrix}\end{array}\right.$ (98)

####  4.2.4  Step 4: IQCs Characterisation of the nonlinearities $\Delta$ and
$\epsilon$

To choose the type of IQCs of the uncertainties considered for this example,
let start with the static non-linearity
$\Delta\left(v(t)\right)=\frac{1}{2}\frac{\rho}{m}v^{2}(t)C_{D}\left(\mathrm{Ma}\right)S\
-av(t)$. It is clear from the drawn Figure 13 that the uncertainty
$\Delta\left(v(t)\right)$ is slope-restricted and sector bounded. Figure 13
provides the IQCs characterisation of this non-linearity and plots the graph
$\Delta\left(v(t)\right)$ where
$\Delta\left(v\left(t\right)\right)=\varphi\left(v,t\right)and\alpha
v\left(t\right)\leq\varphi\left(v,t\right)\leq\beta v\left(t\right)$ for all
$v$. Simple calculation of the slope gives
$\alpha=\left[\begin{matrix}-0.003&0.013&0.006\\\ \end{matrix}\right]^{T}$ and
$\beta=\left[\begin{matrix}0.008&-0.004&-0.004\\\ \end{matrix}\right]^{T}$
representing the bounds of the uncertainty $\Delta.$

Similarly, we have the nonlinearity
$\epsilon\left(\hat{v}\left(t\right)\right)=u\left(\hat{v}\left(t\right)\right)-\hat{u}\
\left(\hat{v}(t)\right)$. The graph of
$\epsilon=\hat{f}\left(\hat{v}\left(t\right)\right)$ is shown in Figure 14
where
$\hat{\alpha}\hat{v}\left(t\right)\leq\hat{f}\left(\hat{v}\left(t\right)\right)\leq\hat{\beta}\hat{v}\left(t\right).$
The IQCs characterisation of this nonlinearity is thus conducted to in
$\hat{\beta}=\left[\begin{matrix}-0.052&0&0.0053\\\ \end{matrix}\right]^{T}$
and $\widehat{\alpha}=\left[\begin{matrix}0.022&-0.024&-0.0043\\\
\end{matrix}\right]^{T}$.

Figure 13:  IQCs characterisation of the non-linearity $\Delta$ . Figure 14:
IQCs characterisation of the non-linearity $\epsilon$.

####  4.2.5  Step 5: IQCs based Robustness Analysis Results

According to the results given by IQC-$\beta$ for the model corresponding to
the original Apollo Robustness Control Problem:

$\displaystyle\left[\begin{matrix}\begin{matrix}\frac{2\gamma_{x}}{1-\gamma_{x}}&\frac{2\gamma_{y}}{1-\gamma_{y}}&\frac{2\gamma_{z}}{1-\gamma_{z}}\end{matrix}&\begin{matrix}\frac{2\gamma_{v_{x}}}{1-\gamma_{v_{x}}}&\frac{2\gamma_{v_{y}}}{1-\gamma_{v_{y}}}&\frac{2\gamma_{v_{z}}}{1-\gamma_{v_{z}}}\\\
\end{matrix}\\\ \end{matrix}\right]^{T}$
$\displaystyle=\left[\begin{matrix}\begin{matrix}0.0988&0.1992&0.1581\\\
\end{matrix}&\begin{matrix}0.0119&0.0244&0.0158\\\ \end{matrix}\\\
\end{matrix}\right]^{T}$

The worst relative error between the closed loop system with the trained
neural controller and the closed loop reference system and is less than
$\left[\begin{matrix}\begin{matrix}9.88\%&19.92\%&15.81\%\end{matrix}&\begin{matrix}1.19\%&2.44\%&1.58\%\\\
\end{matrix}\\\ \end{matrix}\right]^{T}$ of the reference model outputs. To
illustrate further the results of our controlled Lander with our designed
neural controller, we perform some simulation results. Figure 15 and Figure 16
show the response of the system with the trained NN controller for different
two scenarios with different initial position and velocity. Clearly, the
output remains inside the expected bound region.

Figure 15: Simulation Results, Successful Landing Scenario. Figure 16:
Simulation Results, Failed Landing Scenario.

Clearly, the novel approach results can be used to improve the robustness of
the NN controller by making the dynamical error smaller which in turn will
guarantee a better performance of the NN controller when it is implemented.
The worst-case measure estimated here is also an indication of the appearance
likelihood of any undesirable behaviour in the closed loop system when
adopting the designed NN controller. The envelop can be reduced by using IQC
toolbox that is compatible with LPV system also hybrid systems as we have
discontinuity in the activation functions such Relu and Leaky-Relu which
result in hybrid system case when we use the Jacobian of the system.

##  5 Conclusion

A novel technique named keep-close was presented to guarantee the safety and
robustness of feedback systems with NN controllers against various uncertainty
types. The theoretical results we reached for the worst-case analysis
presented provides the user with a prior knowledge about the worst case of the
error between the reference model and the state of the system with a neural
controller including perturbations. The key idea of the new approach is to
maintain the closed-loop system output with a NN controller when its input
changes within a bounded set, close to the output of a robustly performing
closed loop reference model. The keep-close approach provides the user a prior
knowledge about the worst-case RISE and SSE between the reference closed-loop
model and the uncertain system controlled by the neural controller. This
information is very important to provide some levels of validation to AI-based
ML techniques by using robust control theory and increasing the level of trust
for users to adopt these schemes.Through the practical implementation and
testing of two distinct problems, namely the Single-Link Robot Arm Control
Problem and the Deep Guidance and Control of Apollo Lander, the theoretical
results of the innovative approach presented in this paper have been
successfully validated. These real-world experiments have substantiated the
efficacy and adaptability of the proposed approach across diverse application
domains, providing empirical evidence to support its effectiveness and
versatility.

## Acknowledgment

The authors would like to thank Spin-Works as a company and Tiago Amaral
specifically for their software to get generate data.

## References

  * Amaral [2019] Tiago Amaral.  Guidance optimization for mars pinpoint landing with optimal trigger and re-optimization.  2019\.
  * Atkins et al. [2013] Anthony G Atkins, Tony Atkins, and Marcel Escudier.  _A dictionary of mechanical engineering_.  Oxford University Press, 2013.
  * Boyd and Vandenberghe [2001] S Boyd and L Vandenberghe.  Convex optimization with engineering applications.  _Lecture Notes, Stanford University, Stanford_ , 2001.
  * Deka et al. [2022] Shankar A Deka, Dušan M Stipanović, and Claire J Tomlin.  Dynamically computing adversarial perturbations for recurrent neural networks.  _IEEE Transactions on Control Systems Technology_ , 2022.
  * et all [2022] Abdelhafid Zenati et all.  On the admissibility and stability of multi-agent nonlinear interconnected positive systems with heterogeneous delays.  _IEEE Transactions on Automatic Control_ , 2022.
  * Fang et al. [2008] Haijun Fang, Zongli Lin, and Mario Rotea.  On iqc approach to the analysis and design of linear systems subject to actuator saturation.  _Systems & control letters_, 57(8):611–619, 2008.
  * Fazlyab et al. [2020] Mahyar Fazlyab, Manfred Morari, and George J Pappas.  Safety verification and robustness analysis of neural networks via quadratic constraints and semidefinite programming.  _IEEE Transactions on Automatic Control_ , 2020.
  * Fetzer et al. [2017] Matthias Fetzer, Carsten W Scherer, and Joost Veenman.  Invariance with dynamic multipliers.  _IEEE Transactions on Automatic Control_ , 63(7):1929–1942, 2017.
  * Huang et al. [2022] Xiangyu Huang, Chao Xu, Jinchang Hu, Maodeng Li, Minwen Guo, Xiaolei Wang, Yu Zhao, Baocheng Hua, and Yunpeng Wang.  Powered-descent landing gnc system design and flight results for tianwen-1 mission.  _Astrodynamics_ , 6:3–16, 2022.
  * Iannelli et al. [2019] Andrea Iannelli, Peter Seiler, and Andrés Marcos.  Region of attraction analysis with integral quadratic constraints.  _Automatica_ , 109:108543, 2019.
  * Ivanov et al. [2020] Radoslav Ivanov, Taylor J Carpenter, James Weimer, Rajeev Alur, George J Pappas, and Insup Lee.  Verifying the safety of autonomous systems with neural network controllers.  _ACM Transactions on Embedded Computing Systems (TECS)_ , 20(1):1–26, 2020.
  * Ivanov et al. [2021] Radoslav Ivanov, Taylor Carpenter, James Weimer, Rajeev Alur, George Pappas, and Insup Lee.  Verisig 2.0: Verification of neural network controllers using taylor model preconditioning.  In _International Conference on Computer Aided Verification_ , pages 249–262. Springer, 2021.
  * Jonsson et al. [2004] Ulf Jonsson, Chung-Yao Kao, Alexandre Megretski, and Anders Rantzer.  A guide to iqc$\beta$: A matlab toolbox for robust stability and performance analysis.  Technical report, Technical report, 2004. URL http://www. ee. unimelb. edu. au/staff/cykao …, 2004.
  * Kelkar et al. [1983] SS Kelkar, LL Grigsby, and J Langsner.  An extension of parseval’s theorem and its use in calculating transient energy in the frequency domain.  _IEEE Transactions on Industrial Electronics_ , 1(1):42–45, 1983.
  * Khalil [2015] Hassan K Khalil.  _Nonlinear control_ , volume 406.  Pearson New York, 2015.
  * Klumpp [1974] Allan R Klumpp.  Apollo lunar descent guidance.  _Automatica_ , 10(2):133–146, 1974.
  * Lessard et al. [2016] Laurent Lessard, Benjamin Recht, and Andrew Packard.  Analysis and design of optimization algorithms via integral quadratic constraints.  _SIAM Journal on Optimization_ , 26(1):57–95, 2016.
  * Megretski and Rantzer [1997] Alexandre Megretski and Anders Rantzer.  System analysis via integral quadratic constraints.  _IEEE Transactions on Automatic Control_ , 42(6):819–830, 1997.
  * Millour et al. [2015] Ehouarn Millour, Francois Forget, Aymeric Spiga, T Navarro, Jean-Baptiste Madeleine, Luca Montabone, Alizée Pottier, Franck Lefevre, Franck Montmessin, Jean-Yves Chaufray, et al.  The mars climate database (mcd version 5.2).  In _European planetary science congress_ , volume 10, pages 2015–2438, 2015.
  * Nam et al. [2021] Phan Thanh Nam, Tran Ngoc Nguyen, Hieu Trinh, et al.  Comparison principle for positive time-delay systems: An extension and its application.  _Journal of the Franklin Institute_ , 358(13):6818–6834, 2021.
  * Pauli et al. [2022] Patricia Pauli, Julian Berberich, and Frank Allgöwer.  Robustness analysis and training of recurrent neural networks using dissipativity theory.  _at-Automatisierungstechnik_ , 70(8):730–739, 2022.
  * Pfifer and Seiler [2015] Harald Pfifer and Peter Seiler.  Robustness analysis of linear parameter varying systems using integral quadratic constraints.  _International Journal of Robust and Nonlinear Control_ , 25(15):2843–2864, 2015.
  * Rantzer [2015] Anders Rantzer.  On the kalman-yakubovich-popov lemma for positive systems.  _IEEE Transactions on Automatic Control_ , 61(5):1346–1349, 2015.
  * Schulman et al. [2015] John Schulman, Sergey Levine, Pieter Abbeel, Michael Jordan, and Philipp Moritz.  Trust region policy optimization.  In _International conference on machine learning_ , pages 1889–1897. PMLR, 2015.
  * Sharma et al. [2016] V Sharma, V Agrawal, BB Sharma, and R Nath.  Unknown input nonlinear observer design for continuous and discrete time systems with input recovery scheme.  _Nonlinear Dynamics_ , 85:645–658, 2016.
  * Summers and Packard [2010] Erin Summers and Andrew Packard.  L 2 gain verification for interconnections of locally stable systems using integral quadratic constraints.  In _49th IEEE Conference on Decision and Control (CDC)_ , pages 1460–1465. IEEE, 2010.
  * Szegedy et al. [2013] Christian Szegedy, Wojciech Zaremba, Ilya Sutskever, Joan Bruna, Dumitru Erhan, Ian Goodfellow, and Rob Fergus.  Intriguing properties of neural networks.  _arXiv preprint arXiv:1312.6199_ , 2013.
  * Veenman and Scherer [2013] Joost Veenman and Carsten W Scherer.  Stability analysis with integral quadratic constraints: A dissipativity based proof.  In _52nd IEEE Conference on Decision and Control_ , pages 3770–3775. IEEE, 2013.
  * Veenman et al. [2016] Joost Veenman, Carsten W Scherer, and Hakan Köroğlu.  Robust stability and performance analysis based on integral quadratic constraints.  _European Journal of Control_ , 31:1–32, 2016.
  * Wang and Manchester [2019] Ruigang Wang and Ian R Manchester.  Robust contraction analysis of nonlinear systems via differential iqc.  In _2019 IEEE 58th Conference on Decision and Control (CDC)_ , pages 6766–6771. IEEE, 2019.
  * Xiang et al. [2018a] Weiming Xiang, Patrick Musau, Ayana A Wild, Diego Manzanas Lopez, Nathaniel Hamilton, Xiaodong Yang, Joel Rosenfeld, and Taylor T Johnson.  Verification for machine learning, autonomy, and neural networks survey.  _arXiv preprint arXiv:1810.01989_ , 2018a.
  * Xiang et al. [2018b] Weiming Xiang, Hoang-Dung Tran, Joel A Rosenfeld, and Taylor T Johnson.  Reachable set estimation and safety verification for piecewise linear systems with neural network controllers.  In _2018 Annual American Control Conference (ACC)_ , pages 1574–1579. IEEE, 2018b.
  * Yin et al. [2021a] He Yin, Peter Seiler, and Murat Arcak.  Stability analysis using quadratic constraints for systems with neural network controllers.  _IEEE Transactions on Automatic Control_ , 67(4):1980–1987, 2021a.
  * Yin et al. [2021b] He Yin, Peter Seiler, and Murat Arcak.  Stability analysis using quadratic constraints for systems with neural network controllers.  _IEEE Transactions on Automatic Control_ , 67(4):1980–1987, 2021b.
  * Zemouche et al. [2005] Ali Zemouche, Mohamed Boutayeb, and G Iulia Bara.  Observer design for nonlinear systems: An approach based on the differential mean value theorem.  In _Proceedings of the 44th IEEE Conference on Decision and Control_ , pages 6353–6358. IEEE, 2005.
  * Zheng et al. [2016] Stephan Zheng, Yang Song, Thomas Leung, and Ian Goodfellow.  Improving the robustness of deep neural networks via stability training.  In _Proceedings of the ieee conference on computer vision and pattern recognition_ , pages 4480–4488, 2016.