On the stability analysis of deep neural network representations of an optimal state-feedback

TL;DR

This paper investigates the stability of nonlinear systems controlled by neural networks trained to approximate optimal feedback laws, introducing differential algebraic techniques and Taylor maps to analyze stability and robustness.

Contribution

It presents a novel method combining differential algebra and Taylor maps to analyze stability of neural network controllers derived from optimal control theory.

Findings

The proposed method effectively assesses stability margins and robustness.

Different neural network architectures influence stability and delay margins.

Application to a quadcopter demonstrates practical utility.

Abstract

Recent work have shown how the optimal state-feedback, obtained as the solution to the Hamilton-Jacobi-Bellman equations, can be approximated for several nonlinear, deterministic systems by deep neural networks. When imitation (supervised) learning is used to train the neural network on optimal state-action pairs, for instance as derived by applying Pontryagin's theory of optimal processes, the resulting model is referred here as the guidance and control network. In this work, we analyze the stability of nonlinear and deterministic systems controlled by such networks. We then propose a method utilising differential algebraic techniques and high-order Taylor maps to gain information on the stability of the neurocontrolled state trajectories. We exemplify the proposed methods in the case of the two-dimensional dynamics of a quadcopter controlled to reach the origin and we study how different architectures of the guidance and control network affect the stability of the target equilibrium point and the stability margins to time delay. Moreover, we show how to study the robustness to initial conditions of a nominal trajectory, using a Taylor representation of the neurocontrolled neighbouring trajectories.

Tables

Table 1. TABLE I: Bounds for random sampling of initial conditions

State variable	Interval
$y$	$\text{interval}-1010m$
$z$	$\text{interval}-1010m$
$v_y$	$\text{interval}-55m{s}^{-1}$
$v_z$	$\text{interval}-55m{s}^{-1}$
$\theta$	$\text{interval}-\frac{\pi }{4}\frac{\pi }{4}\mathrm{rad}$

Table 2. TABLE III: Stability margins for various G&CNETs

	${\zeta }_{10\%}$	$T$	${\tau }^{\ast }$
${𝒩}_1^{50}$	2.27s	2.41s	0.137s
${𝒩}_2^{50}$	1.66s	–	0.034s
${𝒩}_3^{50}$	1.73s	2.63s	0.023s
${𝒩}_4^{50}$	1.26s	2.52s	0.026s
${𝒩}_5^{50}$	1.19s	3.32s	0.037s
${𝒩}_6^{50}$	0.89s	2.07s	0.039s
${𝒩}_7^{50}$	0.59s	2.00s	0.056s
${𝒩}_8^{50}$	0.71s	8.34s	0.069s
${𝒩}_9^{50}$	0.78s	2.68s	0.091s

Equations

$${\mathbf{u}}^{\ast }(\mathbf{x})=\underset{\mathbf{u}\in \mathcal{U}}{\mathrm{argmin}}\left\{\ell (\mathbf{x},\mathbf{u})+\mathbf{f}(\mathbf{x},\mathbf{u})\cdot {\nabla }_{\mathbf{x}}v(\mathbf{x})\right\}.$$

$$\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\mathcal{N}(\mathbf{x}))$$

$$\mathcal{N}(\mathbf{x}):\enskip\left\{\begin{array}[]{l}L^{(0)}=\mathbf{g}_{\textrm{pre}}(\mathbf{x})\\ L^{(i+1)}=\sigma_{i}(\boldsymbol{W}^{(i)}L^{(i)}+\mathbf{b}^{(i)}),\enskip\forall i=0..l\\ \mathcal{N}(\mathbf{x})=\mathbf{g}_{\textrm{post}}(L^{(l+1)})\end{array}\right.$$

$$a_{ij}=\frac{d{f}_i}{d{x}_j}=\underset{a_{s,ij}}{\underbrace{\frac{\partial f_i}{\partial x_j}}}+\underset{a_{\mathcal{N},ij}}{\underbrace{\sum _{k=1}^m\frac{\partial f_i}{\partial u_k}\text{leavevmode}to27.81pt\text{vbox}to27.81pt\text{pgfpicture}\text{makeatletter}\text{lower}-13.90312ptto0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{definecolor}pgfstrokecolorrgb\mathrm{0,0,0}\text{pgfsys@color@rgb@stroke}000\text{pgfsys@invoke}\text{pgfsys@color@rgb@fill}000\text{pgfsys@invoke}\text{pgfsys@setlinewidth}0.4pt\text{pgfsys@invoke}\text{nullfont}to0.0pt\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@moveto}13.70313pt0.0pt\text{pgfsys@curveto}13.70313pt7.56812pt7.56812pt13.70313pt0.0pt13.70313pt\text{pgfsys@curveto}-7.56812pt13.70313pt-13.70313pt7.56812pt-13.70313pt0.0pt\text{pgfsys@curveto}-13.70313pt-7.56812pt-7.56812pt-13.70313pt0.0pt-13.70313pt\text{pgfsys@curveto}7.56812pt-13.70313pt13.70313pt-7.56812pt13.70313pt0.0pt\text{pgfsys@closepath}\text{pgfsys@moveto}0.0pt0.0pt\text{pgfsys@stroke}\text{pgfsys@invoke}\text{pgfsys@beginscope}\text{pgfsys@invoke}\text{pgfsys@transformcm}1.00.00.01.0-7.95035pt-4.44444pt\text{pgfsys@invoke}\text{definecolor}pgfstrokecolorrgb\mathrm{0,0,0}\text{pgfsys@color@rgb@stroke}000\text{pgfsys@invoke}\text{pgfsys@color@rgb@fill}000\text{pgfsys@invoke}\frac{\partial {\mathcal{N}}_k}{\partial x_j}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{pgfsys@discardpath}\text{pgfsys@invoke}\text{lxSVG@closescope}\text{pgfsys@endscope}\text{hss}\text{lxSVG@closescope}\text{endpgfpicture}}}.$$

$$\dot{\mathbf{x}}(\mathbf{t})=\mathbf{f}(\mathbf{x}(t),\mathcal{N}({\mathbf{x}}^{\tau }))$$

$$\dot{\delta \mathbf{x}}={\mathbf{A}}_s\delta \mathbf{x}+{\mathbf{A}}_{\mathcal{N}}\delta {\mathbf{x}}^{\tau }.$$

$$\left|{\mathbf{A}}_s+{\mathbf{A}}_N{e}^{-\lambda \tau }-\lambda \mathbf{I}\right|=0.$$

$$\delta {\mathbf{x}}^{\tau }$$

$${\dot{\mathbf{x}}}_p$$

$$\left[\begin{array}{c}\delta \dot{\mathbf{x}}\\ {\dot{\mathbf{x}}}_p\end{array}\right]=\left[\begin{array}{cc}{\mathbf{A}}_s+{\mathbf{A}}_{\mathcal{N}}{\mathbf{D}}_p& {\mathbf{A}}_{\mathcal{N}}{\mathbf{C}}_p\\ {\mathbf{B}}_p& {\mathbf{A}}_p\end{array}\right]\left[\begin{array}{c}\delta \mathbf{x}\\ {\mathbf{x}}_p\end{array}\right]$$

$$\mathbf{x}(t)={\Phi }_{\mathcal{N}}(t,{\mathbf{x}}_0)$$

$$\mathbf{x}(t)\approx {\mathcal{M}}_{{\mathbf{x}}_0}^{k,t}(\delta {\mathbf{x}}_0).$$

$$\Vert \delta \mathbf{x}\Vert \le \Vert {\mathbf{A}}_1\Vert \Vert \delta {\mathbf{x}}_0\Vert +\Vert {\mathbf{A}}_2\Vert \Vert \delta {\mathbf{x}}_0{\Vert }^2+\Vert {\mathbf{A}}_3\Vert \Vert \delta {\mathbf{x}}_0{\Vert }^3+\dots$$

$$\Vert \delta \mathbf{x}\Vert \le \sum _{i=0}^k\Vert {\mathbf{A}}_i\Vert {ϵ}^i=\sum _{i=0}^k{b}_i{ϵ}^i$$

$$ϵ=\underset{k\to \infty }{\lim }\frac{|b_k|}{|b_{k+1}|}$$

$$\left\{\begin{array}[]{l}\dot{\mathbf{r}}=\mathbf{v}\\ \dot{\mathbf{v}}=c_{1}\frac{u_{1}}{m}\hat{\mathbf{i}}_{\theta}+\mathbf{g}-\frac{1}{2}\mathbf{v}\\ \dot{\theta}=c_{2}u_{2}\end{array}.\right.$$

$$\begin{array}[]{l}J=\int_{0}^{t_{f}}\left(c_{1}^{2}u_{1}^{2}+c_{2}^{2}u_{2}^{2}\right)dt\\ \end{array}$$

$$\begin{gathered} \delta x_i=a_{1,ij}\delta x_{0,j}+a_{2,ijk}\delta x_{0,j}\delta x_{0,k}+ \\ {a}_{3,ijkl}\delta x_{0,j}\delta x_{0,k}\delta x_{0,l}+\dots \end{gathered}$$

$$\Vert \delta \mathbf{x}\Vert \le \Vert {\mathbf{A}}_1\Vert \Vert \delta {\mathbf{x}}_0\Vert +\Vert {\mathbf{A}}_2\Vert \Vert \delta {\mathbf{x}}_{02}\Vert +\Vert {\mathbf{A}}_3\Vert \Vert \delta {\mathbf{x}}_{03}\Vert +\dots$$

Full text

I Introduction

The optimal feedback of several deterministic, non-linear systems of interest in aerospace applications has been, recently, directly represented by deep neural networks trained using techniques such as imitation learning [1, 2, 3, 4] or reinforcement learning [5, 6, 7]. Regardless of the training details, the neural network is approximating the solution to the Hamilton-Jacobi-Bellman (HJB) equations. The revival of interest in such methods is due to recent advances in deep learning, not limited to: learning algorithms, regularisation techniques, exploitation of GPUs for faster training and large datasets.

The term G&CNET (Guidance and Control Network) is here used to refer to one such representation, in particular to a feedforward, fully-connected neural network trained using supervised learning to approximate the optimal feedback—a function relating the state to the optimal action—of an autonomous, deterministic system. The optimal state-action pairs that constitute the training set are computed by applying Pontryagin’s theory of optimal processes [8], or by solving the HJB equations. In other words, G&CNETs are neural networks imitating the optimal feedback of a nonlinear, autonomous and deterministic system, where the word “imitating” has to be taken in the sense of the machine learning technique of imitation learning [9] of expert’s actions. G&CNETs can be viewed as an alternative to the widely used two-degrees-of-freedom approach to optimal control, based on tracking a pre-computed optimal guidance profile. In the two-degrees-of-freedom approach, the guidance and control problem is decoupled into two tasks: trajectory generation, which is done offline, and the trajectory tracking, which is taken care by an on-board controller or a deep neural network. On the other hand, in the case of G&CNETs both tasks are performed by the neural network in real-time; the remainder of the trajectory is uniquely defined from the current state and perturbations from a nominal trajectory result in alternative profiles. The idea behind G&CNETs is not new and similar schemes were studied in the 90s [3]. It is, though, only in recent years, thanks to advances in the use of neural networks and computer hardware, that the idea was successfully deployed on higher dimensional and more complex systems and thus received renewed attention [4, 1, 2].

Neurocontrollers have been studied in the past and results on their stability obtained in several cases [7, 10, 11, 12]. In particular, stable neurocontrollers can be designed by applying techniques from the adaptive control framework [13]; the adaption of the neurocontroller weights is performed in real-time satisfying a predefined Lyapunov function [10, 11]. Although training of the neural network is entirely avoided, the technique is only applicable to shallow neural networks. Separately, the notion of stability around a nominal trajectory used to analyze controllers which are developed to perform the tracking task in a two-degrees-of-freedom design, cannot be transferred to G&CNETs. Perturbations along a nominal trajectory are not driven to zero, but drift as the neurocontroller applies the approximate learned optimal feedback. Furthermore, the generally perceived black-box nature of neural networks leads to skepticism on the use of such controllers. Especially in cases where safety and validation are of paramount importance, such as the automotive or space industry, other solutions are employed even if neural networks would provide competitive performance.

In this paper we first analyze the linear stability of a system controlled by a G&CNET and then we propose a new method, based on the use of differential algebra and high-order Taylor maps (HOTM) [14], to gain information on the stability of neurocontrolled state trajectories with respect to their initial conditions. After presenting a general methodology, we introduce the case study of a two-dimensional quadcopter dynamics. We train several G&CNETs with different architectures to approximate the power-optimal control of the system (to reach a hovering equilibrium at the origin) and we study their linear behaviour in proximity to the equilibrium point. We then show how to use Taylor models to describe the system behaviour around any neurocontrolled trajectory and obtain stability estimates

111The code and data to reproduce the results in this paper is made available via the github repository https://github.com/darioizzo/neurostability.

II Methodology

In this section we state formally the optimal control problem considered and introduce the formalism used to construct neural network representations of the optimal state-feedback or G&CNETs

Later we present the linear stability study of a neurocontrolled system highlighting the explicit relation between the linearized system and neural network gradients. We then introduce the use of a high-order Taylor model of a neurocontrolled trajectory in proximity of a reference trajectory to study its stability with respect to perturbations of the initial conditions.

II-A Optimal Control problem statement

We study the non-linear, autonomous system $\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\mathbf{u})$, where $\mathbf{x}(\cdot):\mathbb{R}\rightarrow\mathbb{R}^{n}$, $\mathbf{u}(\cdot):\mathbb{R}\rightarrow\mathbb{R}^{m}$ and $\mathbf{f}(\cdot,\cdot):\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$. The free terminal time optimal control problem is then introduced as the problem of finding a control $\mathbf{u}(\cdot)\in\mathcal{U}$ able to steer in the time $t_{f}$ such a system from an initial state $\mathbf{x}_{0}$ to a target state $\mathbf{x}_{f}$ minimizing the cost functional $J(\mathbf{x},\mathbf{u},t_{f})=\int_{t_{0}}^{t_{f}}\ell(\mathbf{x},\mathbf{u})dt$. The functional space $\mathcal{U}$ is the space of all piecewise continuous functions assuming values in some closed region $U\subset\mathbb{R}^{m}$. The relationship between the optimal value attained by $J$ and the initial state $\mathbf{x}_{0}$ is called the value function and is here indicated with $v(\mathbf{x}_{0})$. We restrict our attention to problems in which the cost rate $\ell$ does not depend on time. This means the HJB equations have a time-invariant form. In all points $\mathbf{x}$ where the value function is differentiable, the optimal control is given by [15]:

[TABLE]

Beyond classic solutions, the optimal value function can be given uniquely by considering the viscosity solution to the HJB equations [16]. Furthermore, if equation (1) has a single minimizer, the optimal control problem also admits a unique solution $\mathbf{u}^{*}(\mathbf{x})$ called the optimal feedback. The relation above reveals how the optimal control for such systems is purely reactive, as it depends only on the current system state and not on its history. When controlled by the optimal feedback the system will acquire the desired target state at $t^{*}_{f}$. The resulting optimal trajectories are denoted by $\mathbf{x}^{*}(t)$.

II-B Neural network representation of the optimal feedback

Fundamental problems in aerospace engineering, such as low-thrust spacecraft interplanetary transfers, spacecraft landing, unmanned aerial vehicle control and rocket guidance, all fall into the description above. Finding a solution to these problems is known to require significant computational resources. This is due to the complex structure of the resulting control problem as well as the mathematical and numerical issues connected with the HJB equation or with the application of Pontryagin’s maximum principle. As a consequence, approximating $\mathbf{u}^{*}(\mathbf{x})$, or the value function $v(\mathbf{x})$, is desirable and proposals have been widely researched in the past both non-neural network based [17] and neural network based [10, 5, 18, 1]. We take the approach discussed in [1] and study neural networks, with feedforward, fully-connected architectures, trained on a database of optimal state-action pairs computed by solving, from different initial conditions $\mathbf{x}_{0}$, the two-point boundary value problem arising from the application of Pontryagin’s maximum principle [8]. The resulting trained networks are here indicated with the term G&CNET as they can be used on-board in real time as a substitute for more classical guidance and control systems.

Even before training commences, we can infer several properties of the G&CNET by looking at the state-control pairs comprising the database. One such property is their behaviour after reaching the target state $\mathbf{x}_{f}$; critical for the use of these controllers on-board. This is strongly influenced by the controls associated with $\mathbf{x}_{f}$, which is by far the majority state in the database as it appears once for each optimal trajectory. In particular, the value function gradient $\nabla_{\mathbf{x}}v(\mathbf{x})$ is not uniquely defined at the target point $\mathbf{x}_{f}$ and neither is the corresponding optimal control $\mathbf{u}^{*}(\mathbf{x}_{f})$. As a consequence, depending on the initial conditions chosen, the training database will contain contradictory state-action pairs defining the policy in $\mathbf{x}_{f}$. In most cases the neural network will learn to associate a value close to the average control at $\mathbf{x}_{f}$. Whilst manipulation of the initial conditions used in trajectory generation is a viable approach as done in [2], we propose to overwrite $\mathbf{u}^{*}(\mathbf{x}_{f})$ as a post-processing step following database creation. This is in the spirit of [3] which proposed the exploitation of a priori knowledge (i.e. the specific control problem) when designing the controller. Proposing a suitable $\mathbf{u}^{*}(\mathbf{x}_{f})$ is made easier when we restrict our attention to optimal control problems where $\mathbf{x}_{f}$ is also an equilibrium point for the system under some control $\mathbf{u}_{e}$, that is we assume $\exists!\mathbf{u}_{e}$ such that $\mathbf{f}(\mathbf{x}_{f},\mathbf{u}_{e})=\mathbf{0}$. Consequently, we extend the optimal feedback definition by setting $\mathbf{u}^{*}(\mathbf{x}_{f})=\mathbf{u}_{e}$. From this point onwards we will use $\mathbf{x}_{e}$ to denote the final target state. Note that $\mathbf{x}_{e}$ is a globally asymptotically stable equilibrium point for the optimally controlled system $\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\mathbf{u}^{*})$, as can be promptly shown realizing that the value function $v(\mathbf{x})$ exists unique and is a Lyapunov function for the system (every meaningful value function is a Lyapunov function and every Lyapunov function is a meaningful value function [19]).

We indicate the optimal state feedback as approximated by a G&CNET with $\mathcal{N}(\mathbf{x})\approx\mathbf{u}^{*}(\mathbf{x})$. Substituting this for $\mathbf{u}$ in the system dynamics gives:

[TABLE]

$\mathcal{N}(\mathbf{x})$ has a functional form that is hierarchical in its transformation from state to control. Intermediate computations can be considered on a layer-by-layer basis as we show below:

[TABLE]

The number of layers (network depth) $l$ as well as the weight matrix $\boldsymbol{W}^{(i)}$ and bias vector $\mathbf{b}^{(i)}$ dimensionality across all layers (network width) constitute the network architecture. $\sigma_{i}$ is a non-linear function termed activation function selected for each layer. $\mathbf{g}_{\textrm{pre}}(\cdot)$ and $\mathbf{g}_{\textrm{post}}(\cdot)$ are transformation functions resulting from the data preprocessing step of the standard machine learning pipeline. They are necessary for two reasons: (1) many machine learning algorithms have an assumption on the distribution of the training data; (2) to re-scale neural network outputs allowing for the use of arbitrary activation functions in the output layer. It is worth stating that these functions are unchanged during learning and therefore the same transformation functions are used when we consider different network architectures in Sec. III. The neural network parameters ($\{\boldsymbol{W}^{(i)},\mathbf{b}^{(i)}\}_{i}$) are found during training, typically using some variant of stochastic gradient descent. This stochasticity in the optimisation algorithm means that in repeated training runs the fitted parameters could be different and therefore properties of the neurocontrolled system (2) could change. We note that even for networks with a single hidden layer, the expanded form shown in equation (3) provides little mathematical insight on the resulting control structure. Furthermore, its highly non-linear form gives rise to skepticism on the use of such controllers in applications where the satisfaction of requirements must be proved using formal mathematical tools. Despite this, it is worth stating that $\mathcal{N}(\mathbf{x})$ is a function in the differentiability class $\mathrm{C}^{k}$ where $k$ is the lowest amongst all the differentiability classes of its activation functions $\sigma_{i}$. This suggests that the straightforward deployment of techniques such as linearisation and even high-order Taylor map approximations is possible.

In this work, we restrict our attention to neural networks that use softplus activation functions ($\sigma_{\textrm{soft}}(x)=\log(1+\exp(x))$) in the hidden layers and tanh activation functions in the output layer (see Fig. 1). Consequently, the network outputs are continuous and differentiable everywhere since they belong to differentiability class $C^{\infty}$. This avoids issues that could arise when the popular ReLU activation function ($\sigma_{\textrm{relu}}(x)=\textrm{max}(0,x)$) is used. We also note that for the dynamical model considered in section III, comparable training performance was observed when softplus activation functions were used instead of ReLUs [2]. Since a G&CNET is expected to approximate the optimal feedback $\mathbf{u}^{*}$, a good requirement on the neurocontrolled system described by equation (2) is to also produce in $\mathbf{x}_{e}$ a globally stable equilibrium point.

II-C Neurocontroller dynamics linearisation

II-C1 Linear stability of the equilibrium

We stated earlier that $\mathbf{x}_{e}$ is a globally asymptotically stable equilibrium point for the optimally controlled system. Here we demonstrate that the asymptotic stability of $\mathbf{x}_{e}$, at least in the local sense, for the system controlled by a G&CNET can be verified. We also emphasise that the representation of the optimal feedback by a neural network allows for the evaluation of linear stability margins. This is not available for $\mathbf{u}^{*}(\mathbf{x})$ as the extended definition has a discontinuity point at $\mathbf{x}_{e}$ where the control instantaneously takes the value $\mathbf{u}_{e}$. This means the optimal dynamics in a neighbourhood of $\mathbf{x}_{e}$ does not follow laws of type $\exp(-ct)$ and thus cannot be described by a linear system.

Consider the neurocontrolled system in equation (2) and derive, around its equilibrium point $\mathbf{x}_{e}$, the linearized form $\dot{\mathbf{x}}=\mathbf{A}\mathbf{x}$ where $\mathbf{A}=\nabla\mathbf{f}\left|{}_{\mathbf{x}=\mathbf{x}_{e}}\right.$. Let us expand this expression further to highlight the relationship between the linearized dynamics and the neural network gradients. Indicating the matrix $\mathbf{A}$ through its components $a_{ij}$ we have:

[TABLE]

This shows that the linearized dynamics is the sum of two matrices: $\mathbf{A}_{s}$ with components $a_{s,ij}$ representing the system dynamics, and $\mathbf{A}_{\mathcal{N}}$ with components $a_{\mathcal{N},ij}$ representing the network feedback. It is then clear that to compute the linear dynamics one needs to compute the derivatives of the neural network outputs with respect to its inputs. This task is most efficiently solved using automatic differentiation, a technique of widespread use in deep learning as it forms the core of the backpropagation algorithm used to train neural networks. As a consequence, there are an abundance of software frameworks that implement this (e.g. Tensorflow, Pytorch).

The dynamics of the neurocontrolled system in a local regime about the equilibrium point can be determined by the eigenvalues $\lambda_{i}=\alpha_{i}+j\omega_{i}$ of $\mathbf{A}$. Briefly, an initial perturbation $\delta\mathbf{x}_{0}$ will evolve as $\delta\mathbf{x}=e^{\mathbf{A}t}\delta\mathbf{x}_{0}$, where we have introduced the variable $\mathbf{x}=\mathbf{x}_{e}+\delta\mathbf{x}$. The solution $\delta\mathbf{x}=e^{\mathbf{A}t}\delta\mathbf{x}_{0}$ is a linear combination of terms in the form $e^{-\lambda_{i}t}$. Therefore, $\alpha_{i}$ defines the asymptotic behavior of the system and $\omega_{i}$ determines the frequency $f=\frac{\omega_{i}}{2\pi}$ of its oscillatory behaviour.

II-C2 Stability to time delay

The principal reason for the use of G&CNETs is for on-board, real-time optimal control [1]. In such contexts, the time scale to apply the optimal action is typically short. This makes studying the neurocontroller behaviour with respect to feedback delays important for ensuring the robustness of the proposed controller. To this effect, we introduce the time delay $\tau$ in the system dynamics as follows:

[TABLE]

where $\mathbf{x}^{\tau}=\mathbf{x}(t-\tau)$. Essentially, we assume that the current action is calculated based on the system state delayed by $\tau$. The linearized state space representation around $\mathbf{x}_{e}$ is:

[TABLE]

This linear system with time delay has the general solution [20]: $\delta\mathbf{x}=e^{\mathbf{\Phi}t}\mathbf{\delta}\mathbf{x}_{0}$ where the eigenvalues of $\mathbf{\Phi}$ are obtained from the time-delayed characteristic equation:

[TABLE]

We note that for $\tau=0$ (no delay) we recover the stability study of the non-delayed system. The solution to the above equation forms the root locus of the system for different delay values. Detecting the first occurrence of the imaginary axis crossing by one of the roots allows us to determine the critical time delay, $\tau^{*}$, that destabilizes the system. In practice, however, the solution to equation (7) comes with quite some issues associated to the initial guess and the non-linear nature of the equation. While eigenvalue tracking approaches can help to alleviate such issues, the resulting numerical methods are cumbersome and inefficient. As a consequence, in this work we use the Padé approximation of order 5 for the time delay, to obtain the initial guess for solving equation (7). The time delay in the Laplace domain can be approximated as a rational function of two polynomials of equal degree, $e^{-s\tau}\approx P(s)/Q(s)$, where $s$ is the Laplace variable. Transforming the approximation in the time domain and applying it separately on each state of the original system results in:

[TABLE]

where the matrices $\mathbf{A}_{p}$, $\mathbf{B}_{p}$, $\mathbf{C}_{p}$ and $\mathbf{D}_{p}$ depend on $\tau$, and $\mathbf{x}_{p}$ are internal states of the system. The input to this linear system is the current state $\delta\mathbf{x}$ and the output is the approximation of the delayed state $\delta\mathbf{x}^{\tau}$. The interconnection of the delayed system (6) with the approximation (8-9)

[TABLE]

gives a linear system with an augmented state. An estimate for $\tau^{*}$ can be obtained by calculating the eigenvalues of the system matrix.

II-D High-order Taylor maps of neurocontrolled trajectories

So far, we have only applied standard tools of control theory to the study of the system described by equation (2). By doing so, we were able to look into the behaviour of a G&CNET controlled system in a neighbourhood of the equilibrium point. We now propose a method based on differential algebra and high-order Taylor maps [21] to study the solutions of (2) in a neighbourhood of a reference solution (or nominal trajectory). The basic idea is to use differential algebraic techniques to compute a high-order Taylor representation of such solutions. These expansions can then be be used to analyze the global behaviour of the neurocontrolled system to parametric uncertainties. In contrast to the standard approach, we do not seek to investigate if small perturbations along a given tracked trajectory die out. Instead, we look into the adjustments the neurocontroller makes as it is perturbed from an initial nominal trajectory defined by $\mathbf{x}_{0}$. We note that this same technique can also be applied to the analysis of perturbations with respect to other model parameters.

Let us indicate the solution to the initial value problem of (2) with:

[TABLE]

where the subscript reminds us that the system dynamics is controlled via a neural network $\mathcal{N}$. We may then compute its Taylor representation of order $k$ and write:

[TABLE]

The map $\mathcal{M}^{k,T_{j}}_{\mathbf{x}_{0}}$ describes, to order $k$ and at time $t$, the solution to the initial value problem with initial condition $\mathbf{x}_{0}+\delta\mathbf{x}_{0}$. Such representations are called Taylor models, or high-order Taylor maps, and are computed efficiently using differential algebraic techniques such as those implemented in the open source python library pyaudi [22] and used here. We compute such maps by integrating numerically, with a variable step Runge-Kutta-Fehlberg scheme, the dynamics of the system from the initial condition $\mathbf{x}_{0}+\delta\mathbf{x}_{0}$ up to the time $T$. As the numerical integration is defined over generalized dual numbers [22] rather than floating point numbers, the resulting $\mathcal{M}^{k,T_{j}}_{\mathbf{x}_{0}}$ are computed at times $T_{j}$ defined by the adaptive step size procedure of the numerical integrator (see [23] for a recent overview of differential algebra and its use to compute high-order Taylor maps, or [14] for an earlier work).

Each map $\mathcal{M}^{k,T_{j}}_{\mathbf{x}_{0}}$ is a $k$-th order Taylor polynomial in $\delta\mathbf{x}_{0}$ and can be used to study the neurocontroller behaviour around a nominal trajectory since it represents an explicit form of the solutions of (2), at least in a ball $B_{\epsilon}$ of radius $\epsilon$ around the nominal trajectory initial conditions $\mathbf{x}_{0}$. Consequently, starting from an initial perturbation in $B_{\epsilon}$, the asymptotic behaviour of the system will then be described by $\mathcal{M}^{\infty,t}_{\mathbf{x}_{0}}$ when $t\rightarrow\infty$. This fact can be used to verify numerically if the neurocontrolled system will eventually reach a target point under the effect of external disturbances. In the following we discuss the problem of providing an estimate for $\epsilon$.

II-D1 Estimating the convergence radius of the Taylor model

Let us introduce $\delta\mathbf{x}=\mathcal{M}^{k,t}_{\mathbf{x}_{0}}-\Phi_{\mathcal{N}}(t,\mathbf{x}_{0})$, the evolution of the initial perturbation $\delta\mathbf{x}_{0}$. In order to obtain an estimate for the convergence radius of this Taylor series, we bound the magnitude of $\delta\mathbf{x}$ as:

[TABLE]

where $\mathbf{A}_{i}$ are the matrices resulting from unfolding the tensor of the $i$-th Taylor term, i.e. $i=1$ represents the gradient, $i=2$ the Hessian and so forth. A detailed derivation of (10) is given in the Appendix. As a result for an expansion of order $k$, we can write:

[TABLE]

where we introduced, for shortness, $b_{i}$ as the corresponding matrix norms and $\lVert\delta\mathbf{x}_{0}\rVert=\epsilon$. We may then study the convergence radius of the series $\sum_{i=0}^{\infty}b_{i}\epsilon^{i}$ to obtain a conservative estimate of the radius of convergence of $\mathcal{M}^{\infty,t}_{\mathbf{x}_{0}}$. Applying the ratio test (D’Alembert’s criterion) we have: $\epsilon\lim_{k\rightarrow\infty}\frac{|b_{k+1}|}{|b_{k}|}<1$ and hence we can define:

[TABLE]

and guarantee that the series $\mathcal{M}^{\infty,t}_{\mathbf{x}_{0}}$ converges in the ball $B_{\epsilon}$. Therefore, starting from initial conditions perturbed in $B_{\epsilon}$, $\mathbf{x}_{e}$ will be an asymptotically stable point in $B_{\epsilon}$ if and only if $\lVert\mathbf{A}_{j}\rVert\rightarrow 0\quad\forall j$, as $t\rightarrow\infty$. Such a condition can easily be verified numerically up to $j=k$ once $\mathcal{M}^{k,t}_{\mathbf{x}_{0}}$ is computed for a sufficiently large $t$.

To summarize, the use of the techniques described in this section allow us to study with respect to G&CNETs: a) the local stability of the equilibrium point $\mathbf{x}_{e}$; b) the linear time delay margin of the system; c) the stability of a nominal neurocontrolled trajectory originating from $\mathbf{x}_{0}$. In the next section we will use these techniques to study a practical test case: the two-dimensional dynamics of a quadcopter.

III Numerical experiments

III-A The quadcopter model

Consider the system shown in Fig. 2 representing a quadcopter whose two-dimensional dynamics is defined by the following set of ordinary differential equations:

[TABLE]

The quadcopter state includes its position $\mathbf{r}=(y,z)$, its velocity $\mathbf{v}=(v_{y},v_{z})$ and its orientation $\theta$. We also refer to the state using the variable $\mathbf{x}=[\mathbf{r},\mathbf{v},\theta]$ and to the controls using the variable $\mathbf{u}=[u_{1},u_{2}]$. In the rest of this paper, we use data from the ®Parrot Bebop drone. The mass of the quadcopter is set to be $m=$0.389,05\text{,}\mathrm{kg} and the acceleration due to the Earth’s gravity is $\mathbf{g}=(0,-g)$ where $g=$9.81\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-2}. Note that the dynamics of the quadcopter is also perturbed by a drag term and its coefficient is set to $\frac{1}{2}$. The control $u_{1}\in[0,1]$ corresponds to a thrust action applied along the direction ${\mathbf{i}}_{\theta}=[\sin\theta,\cos\theta]$ bounded by a maximum magnitude $c_{1}=$9.1\text{,}\mathrm{N}. The control $u_{2}\in[-1,1]$ models a pitch rate bounded by $c_{2}=$35\text{\,}\mathrm{rad}\text{\,}{\mathrm{s}}^{-1}.

For the above non-linear dynamical system we consider the optimal control problem of steering the state from any initial state to the target state, $\{\mathbf{r}=(0,0),\mathbf{v}=(0,0),\theta=0\}$ (i.e. $\mathbf{x}_{f}=\mathbf{O}$), minimizing the cost function:

[TABLE]

We solve the resulting control problem from 200,000 initial conditions randomly sampled from the intervals reported in Table I. For each optimal trajectory we insert in a database the state-action pairs at 59 equally spaced points, taking care to eliminate the last point and substitute it with $(\mathbf{x}_{f},\mathbf{u}^{*}(\mathbf{x}_{f}))$ where $\mathbf{u}^{*}(\mathbf{x}_{f})=[\frac{mg}{c_{1}},0]$ thereby ensuring the origin is an equilibrium point for the system (12). We then train neural networks with varying architectures $\mathcal{N}_{l}^{\eta}$, where $l$ is the network depth and $\eta$ is the network width held constant for all hidden layers, as a regression task on the database containing 11,800,000 optimal state-action pairs.

We report in Table II the performance (mean absolute error evaluated on a held-out partition of the training database) of networks with depth $l\in\{1..9\}$ and width $\eta\in\{50,100,200\}$. For further details on the training procedure see [2].

A consequence of approximating the optimal feedback using neural networks is that the equilibrium point of the system (12) as controlled by a G&CNET $\mathbf{u}=\mathcal{N}_{l}^{\eta}(\mathbf{x})$ may not be exactly at the origin. In our experiments we find that the offset of the equilibrium from the origin is rather small, between $10^{-3}$ and $10^{-2}$, particularly when compared to the intervals used for sampling of the initial conditions. In our case, we can make the origin the equilibrium point for all G&CNETs simply by evaluating the true equilibrium point $\hat{\mathbf{x}}_{l}^{\eta}$ followed by shifting the axis i.e. $\mathbf{x}\rightarrow\mathbf{x}+\hat{\mathbf{x}}_{l}^{\eta}$.

We show in Fig. 3 a few trajectories resulting from the use of one particular G&CNET $\mathcal{N}_{3}^{100}$ to control the quadcopter dynamics. The optimal trajectories, as evaluated by solving the two-point boundary value problem resulting from the application of Pontryagin’s maximum principle, are indistinguishable from the ones produced by the network.

III-A1 Hovering stability and stability to time delay

After shifting our axis, the quadcopter dynamics (12) controlled by each G&CNET $\mathbf{u}=\mathcal{N}_{l}^{\eta}(\mathbf{x})$ has in the origin $\mathbf{O}$ its equilibrium point and thus $\mathcal{N}_{l}^{\eta}(\mathbf{O})=[\frac{mg}{c_{1}},0]$. This corresponds to the quadcopter hovering over the origin and cancelling exactly the gravitational pull with its upward thrust. We thus proceed to study the behavior of the quadcopter close to the hovering conditions by computing all the eigenvalues $\lambda_{i}$ of $\mathbf{A}$ resulting from (4), for each of the trained networks. In all cases, the modes computed resulted to be stable, that is $\alpha_{i}<0\;\forall i=1..5$, proving that all trained networks result in a stable dynamics near the hovering position. For each eigenvalue, we then proceeded to compute the decay time, $\zeta_{{10\%}_{i}}=\frac{\ln(0.1)}{\alpha_{i}}$, which is defined as the time the envelope associated to the $i$-th mode takes to decay to 10% of its starting value, and the oscillating period, $T_{i}=\frac{2\pi}{\omega_{i}}$, associated to the $i$-th mode. In Table III we report the largest values for $\zeta_{{10\%}}$ and $T$ for each network $\mathcal{N}_{l}^{\eta}$. At most one pair of complex conjugate eigenvalues was present in all cases tested. In three cases no oscillatory behaviour is present as all eigenvalues are real.

We then perform the analysis of the time-delayed system solving (7) for different values of the time delay $\tau$ and thus producing a root locus which we show in Fig. 4 for two G&CNETs: $\mathcal{N}_{100}^{3}$ and $\mathcal{N}_{200}^{4}$. The eigenvalues of the non-delayed system are also indicated in the figure. This shows that time delays initially have a stabilizing effect on the most stable roots but later give raise to instability represented by the crossing of the imaginary axis. The computed values for the critical time delay $\tau^{*}$ are reported in Table III.

From previous work [1, 2], we know that neural networks with greater depths and widths (i.e. increased parameterisation) are better representations of the optimal feedback. This is repeated in Table II for which we observe the largest reduction in error when the network goes from 1 to 2 hidden layers. However, the stability margins characterising the neurocontroller behaviour in a neighbourhood of the equilibrium point (columns 1 and 2 of Table III) indicate no clear relationship to the network architecture (and thus by extension the optimal feedback approximation performance). One could argue that this is not surprising for two reasons, firstly the metric in Table II is a measure of global performance, and secondly, as noted in section II, the neurocontroller behaviour in a neighbourhood of the equilibrium point is not expected to match the optimal behaviour. Despite this, we observe the critical time delay (column 3 of Table III) for G&CNETs with depth 1 to be considerably higher than the others. In the context of the on-board implementation of G&CNETs, this presents an interesting choice between greater optimality or improved robustness to time delays.

III-A2 High-order Taylor maps for the neurocontrolled quadcopter.

Let us consider, as an example, the nominal trajectory of the quadcopter resulting from the neurocontroller $\mathcal{N}_{100}^{3}$ and the initial condition $\mathbf{x}_{0}=[-4,0,0,0,0]\,$\mathrm{m}. The trajectory, also shown in Fig. [3](#S2.F3)b, is a manoeuvre where the quadcopter changes its position horizontally by $4\text{\,}\mathrm{m}$. We proceed to study the stability of such a nominal trajectory with respect to perturbations of the initial conditions $\mathbf{x}_{0}+\delta\mathbf{x}_{0}$. In particular we would like to have an indication on the neurocontroller being able to bring the quadcopter to the final target state for all perturbations in some ball $B_{\epsilon}$. Note that, for this manouvre, the corresponding optimal time is $t^{*}_{f}=$2.06\text{\,}\mathrm{s}. Following the approach detailed in section II, we use differential algebraic techniques and a Runge–Kutta–Fehlberg adaptive step numerical integration scheme to compute the high-order Taylor maps $\mathbf{x}_{T}\approx\mathcal{M}^{k,T_{j}}_{\mathbf{x}_{0}}(\delta\mathbf{x}_{0})$ representing the Taylor expansion of order $k$ of the quadcopter state at the time grid points $T_{j}$ defined by the adaptive stepper of the numerical integration scheme. The Taylor expansion is taken with respect to perturbations $\delta\mathbf{x}$ of the nominal initial conditions $\mathbf{x}_{0}$. Since $n=5$, i.e. the state has a dimension of five, each resulting Taylor polynomial contains 5, 20, 55, 125, 251, 461, 791 terms as the order grows from linear. Clearly, if: $\lim_{k\rightarrow\infty,t\rightarrow\infty}\mathcal{M}^{k,t}_{\mathbf{x}_{0}}(\delta\mathbf{x}_{0})=\mathbf{x}_{e}$ then the neurocontroller is guaranteed to be able to drive the quadcopter to the final desired hovering position $\forall\delta\mathbf{x}_{0}\in B_{\epsilon}$. When this is the case, a formal guarantee on the stability of the nominal trajectory is also obtained. Taking a numercial approach, since we cannot compute the Taylor model for $t\rightarrow\infty$, nor can we compute the map at order $k\rightarrow\infty$, we stop the numerical integration at $T=\frac{3}{2}t_{f}^{*}$ (with $t^{*}_{f}$ representing the time of the optimal manoeuvre) giving the system enough time to reach the equilibrium position and stabilize itself. Inspecting the Taylor map computed for $t=T$ reveals that $\lVert\mathbf{A}_{j}\rVert<1e^{-3}\quad\forall j\in[1..7]$, and allows us to conclude that the equilibrium position is indeed acquired, with the said precision, from all initial conditions in an $\epsilon$ ball around the nominal $\mathbf{x}_{0}$. The ball size $\epsilon$ is then computed from (11) using $k=7$.

We can actually estimate the convergence radii for the map at all $T_{j}$. The result is shown in Fig. 6 where $\epsilon\approx\frac{|b_{6}|}{|b_{7}|}$ is plotted for the map $\mathcal{M}_{\mathbf{x}_{0}}^{7,t}$ against $t$ (the times $t$ are sampled in the grid $T_{j}$ defined by the adaptive stepper of the numerical integration scheme). The convergence radius at $T=\frac{3}{2}t_{f}^{*}$ is approximately 0.1. Around the optimal time $t^{*}_{f}$ we also observe a dramatic decrease of the estimated convergence radius of the map. In order to confirm this observation we simulate in Fig. 5 the quadcopter state during the nominal trajectory and the state predicted by the various maps $\mathcal{M}^{7,T_{j}}_{\mathbf{x}_{0}}(\delta\mathbf{x}_{0})$ for a rather large perturbation $\delta\mathbf{x}_{0}=[0.1,0.1,0.1,0.1,0.1]\,$\mathrm{m}. We also simulate the ground truth by numerical integration of ([2](#S2.E2)) using the initial conditions $\mathbf{x}_{0}+\delta\mathbf{x}_{0}$. As expected, the various maps represent the new optimal trajectory, corresponding to the perturbed initial conditions, quite accurately for most time instants, failing catastrophically for $t\approx$2\text{\,}\mathrm{s}, where the selected perturbation is evidently larger than the corresponding map convergence radius.

IV Conclusions

We show a general methodology to analyse the behaviour of neurocontrollers. Our method enables the study of neurocontrolled feedback dynamics in proximity of the equilibrium point by deriving linear stability margins as well as time delay margins. We then propose the use of high-order Taylor models to study the robustness of a nominal trajectory to perturbations on its initial conditions. The methodology is successfully applied to the test case of a two-dimensional quadcopter dynamical model where we show, for several G&CNETs with varying architectures, the formal stability guarantees for the neurocontrolled hovering behaviour as well as for the stability of a nominal manoeuvre. We found that as soon as the network is deep, i.e. has more than one hidden layer, there seem to be no relation between the linear stability margins obtained and the network architecture. We also found that the computed Taylor model of the state describes well the neurocontrolled trajectories in a neighbourhood of a nominal trajectory and we propose a method to estimate the size of said neighbourhood. Our results constitute a first step to increasing the trust on the on-board use of an optimal feedback controller represented by a neural network, thus narrowing the gap between control theory and machine learning.

[Tensor unfolding] Let us write the $i$-th component of $\delta\mathbf{x}=\mathcal{M}^{k,t}_{\mathbf{x}_{0}}-\Phi_{\mathcal{N}}(t,\mathbf{x}_{0})$, namely the time evolution of an initial perturbation $\delta\mathbf{x}_{0}$, using the Einstein notation

[TABLE]

where $\delta x_{0,i}$ denotes the $i$-th component of $\delta\mathbf{x}_{0}$. The terms $a_{1}$, $a_{2}$ and $a_{3}$ are components of the gradient, the Hessian and the third-order partial derivatives tensors. We take the norm of the previous relation and we apply the triangle inequality and the properties of the induced matrix norms:

[TABLE]

To write the previous equation in a matrix-vector notation, we have unfolded the tensors such that all the partial derivatives corresponding to the same $\Phi_{\mathcal{N},i}$ are in the same row, e.g. for $\mathbf{A}_{2,ij}=a_{2,ikl}$ and $\mathbf{A}_{3,ij}=a_{3,iklq}$ with $k=(j-1)\mod n+1$, $l=(j-1)\div n+1$ and $q=(j-1)\div n^{2}+1$, $j=1,2,\ldots,n^{2}$ and $j=1,2,\ldots,n^{3}$ for each matrix respectively. We apply the same unfolding to the vectors $\delta\mathbf{x}_{02}$ and $\delta\mathbf{x}_{03}$ as $\delta\mathbf{x}_{02,j}=\delta x_{0,k}\delta x_{0,l}$ and $\delta\mathbf{x}_{03,j}=\delta x_{0,k}\delta x_{0,l}\delta x_{0,q}$ where the indices $k$, $l$ and $q$ depend on $j$ as before. Therefore, $\lVert\delta\mathbf{x}_{02}\rVert=\lVert\delta\mathbf{x}_{0}\rVert^{2}$ and $\lVert\delta\mathbf{x}_{03}\rVert=\lVert\delta\mathbf{x}_{0}\rVert^{3}$ under any $p$-norm.