# Protein-protein docking using a tensor train black-box optimization method

Dmitry Morozov<EMAIL_ADDRESS>Terra Quantum AG, Kornhausstrasse 25,
9000 St. Gallen, Switzerland  Artem Melnikov  Terra Quantum AG,
Kornhausstrasse 25, 9000 St. Gallen, Switzerland  Vishal Shete  Terra Quantum
AG, Kornhausstrasse 25, 9000 St. Gallen, Switzerland  Michael Perelshtein
<EMAIL_ADDRESS>Terra Quantum AG, Kornhausstrasse 25, 9000 St. Gallen,
Switzerland

###### Abstract

Black-box optimization methods play an important role in many fields of
computational simulation. In particular, such methods are often used in the
design and modelling of biological systems, including proteins and their
complexes with various ligands. This work is mainly focused on the protein-
protein docking that plays a key role in modern drug-design workflows. We
develop a black-box approach for such docking problems using a novel technique
based on the tensor-train decomposition of high-dimensional interaction
functions. Our method shows an advantage in terms of the discovered global
minima and has a high potential for further implementation on a wide range of
devices, including graphical processing units and quantum processing units.

##  I Introduction

Protein-protein interactions play a crucial role in organisms [1].
Understanding how and why individual molecules could form complexes of a very
distinct structure is challenging. There are experimental methods to determine
the structures of complexes, like crystallography and CryoEM techniques, which
are typically deposited into the Protein Data Bank (PDB) [2]. However, they
can cover only a very tiny fraction of all biological complexes. To overcome
these limitations, computational methods for predicting and building protein
complexes are becoming more popular, as a complementary addition to the
experimental techniques.

Molecular docking is one of the major methods in studying protein-protein
interactions at the atomistic scale [3, 4]. Essentially, the method aims to
sample the geometrical space of two or more proteins and search for the best
poses of one protein (the ligand) with respect to the other protein (the
receptor) using scoring functions. The main challenges are: (i) to design the
scoring function to give an acceptable precision; and (ii) to efficiently
sample the conformational space to find the best possible structures.

Today’s most widely used sampling methods are based on the Fast Fourier
Transform (FFT). They are capable of generating geometrically complementary
rigid-body docking poses [5]. However, they are typically limited in the use
of very specific scoring functions and thus, in many cases require an a
posteriori recalculation of the scores for the obtained structures with
additional external software [6]. On the other hand, there are existing black-
box approaches for searching minimas in high-dimensional spaces. As an
example, swarm intelligence algorithms, such as Particle Swarm Optimization
(PSO) [7] and Glowworm Swarm Optimization (GSO) methods [8], have been
successfully applied to the protein docking problem using a large variety of
scoring functions [9].

In this work, we are aiming to apply another class of black-box global
optimization algorithms based on the tensor train (TT) decomposition of the
high-dimensional scoring (or energy) function [10] that were developed over
the last years [11]. The expected advantages of such methods is their
agnosticism to the initial conditions and high-efficiency in sampling even
very heterogeneous coordinate spaces that have a large number of local minima
in their structure. These properties make such methods promising candidates
for exploring interactions in biological complexes.

Another advantage of the tensor train methods developed by Terra Quantum is
their capability for efficient implementation on both CPUs and GPUs. As can be
seen from the results in this paper, usage of these algorithms already
produces highly promising results. Moreover, the tensor train algorithms
possess a structure similar to quantum circuits [12, 13], which paves the way
towards their future implementation on Quantum Processing Units (QPUs) once
such QPUs are of sufficient scale and fidelity [14]. Thus, we are able to
harness some advantages of our Quantum Software approach [15], which will grow
once quantum hardware matures and numerical methods of quantum information
processing via tensor networks advances.

##  II Methods

Three different algorithms for solving docking problem were employed in the
study: the glowworm swarm optimization (GSO) algorithm as implemented in the
LightDock [16] package; the tensor train based minima search implemented in
the ttpy library (TTPY) [17]; and our own tensor train based black-box
optimization method packaged in a Python library (TTOPT), which will be
described in the following section.

###  II.1 Tensor Train assisted docking

Here we will describe our workflow for performing docking simulations using
the optimization algorithm based on the Tensor Train (TT) decomposition [10,
11, 18]. The overall scheme is shown in Fig. 1. Initially, as in the other
docking protocols, the preparation and cleanup of the ligand and receptor
structures should be performed. At that stage, we also decide on the target
(scoring) function used to evaluate the docking poses of the ligand molecule
with respect to the receptor. This stage could also involve fixing the missing
residues and side-chains, adding hydrogenes, building molecular topologies and
partially relaxing the structures.

Figure 1: Workflow for performing docking simulations with the TTOPT approach.
(a) - structure preparation, (b) - grid and target function definition, (c) -
decomposition of the target function into a tensor train and search for optima
on a grid, (d) - clusterization of the obtained minima geometries.

For the next stage, one would need to define a grid spanned around the
receptor Fig. 1b in which optimization of the target function with respect to
the ligand position will be performed. It could be both minimization of the
energy or maximization of the scoring. In general, our TT algorithm requires
discretization of the search on a uniform grid. We define $d$ as the number of
variables (dimensions) and $ngrid$ as the grid size in one dimension. In the
current work, we use rigid-body models of the receptor and ligand proteins
meaning that the total search space has 7 coordinates in total: 3 coordinates
for the ligand’s center of masses and 4 coordinates for its rotation
quaternion. It should be noted that while here we are focused mainly on rigid-
body protein-protein docking, the same approach is also applicable to other
systems such as small molecules and nucleic acids, where the flexibility of
these systems could be straightforwardly treated as additional dimensions in
the TT decomposition.

The TTOPT algorithm would be used to search for the possible optima of the
target function in the defined grid space. Unlike brute-force grid search,
this algorithm does not estimate the cost function at all points of the grid
but dynamically provides the next set of evaluating points in the search space
based on the knowledge accumulated during all previous evaluations. The TTOPT
algorithm applies the cross-approximation technique to the tensors [19], which
in turn is based on the MaxVol routine [20]. Essentially, it performs a
decomposition of the multi-dimensional target function into a compressed form
as the product of small tensors connected into a chain (Fig. 1c) and expanded
on a grid.

In addition, the iterative power algorithm by Soley et.al.[21] is used for the
global minima search within the TT representation. The total number of the
target function evaluations required by our approach to find global minima
could be approximated with $O(nswp*d*ngrid*rmax^{2})$, where $d$ is the number
of dimensions, $ngrid$ is the grid-size, $rmax$ is the maximum rank of the
tensors with which the algorithm tries to approximate the target function via
tensor train decomposition and $nswp$ is the number of iterations which the
algorithm performs to improve the decomposition. Overall, our TTOPT approach
scales much better with the system size than the $ngrid^{d}$ brute force grid
search and shows much better parallelization capabilities than other black-box
methods, such as Bayesian optimization or genetic algorithms. We implement the
TTOPT core functional using C++ and package the whole framework in a Python
library.

Finally, after the global optima search, we perform an analysis of the
obtained receptor-ligand complex structures in terms of the target function
values and clusterize them according to the root mean square deviation (RMSD)
of the ligand C$\alpha$ protein backbone atoms (Fig. 1d). Obtained in such a
way, the best complexes from all found clusters could be further used to
perform a comparison with the experimental data or for conducting simulations
such as molecular dynamics, binding free-energy evaluations or direct
mutagenesis using more elaborate approaches.

###  II.2 Systems for docking

To perform docking simulations, we have used 10 systems from the CAPRI scoring
benchmark [22]. The Protein Data Bank (PDB) structures along with the chain
identifiers, which have been used as the receptor and ligand parts of the
respecting complexes to perform the simulations, are summarized in Table 1.
The fast implementation of the DFIRE energy function [23] provided in the
LightDock software package was used as a target function for performing global
minima search.

Table 1: Systems considered in this work. Chain - ID taken from the original PDB, Length - number of the residues in the Chain, Weight - molecular weight of the Chain in kDa  PDBID | Receptor | Ligand
---|---|---
Chain | Length | Weight | Chain | Length | Weight
3BX1 | A | 269 | 26.7 | C | 181 | 19.9
3U43 | B | 132 | 15.3 | A | 92 | 10.0
3R2X | AB | 492 | 56.8 | C | 82 | 10.6
3Q87 | B | 164 | 19.2 | A | 122 | 14.2
2WPT | B | 113 | 15.2 | A | 82 | 9.9
3E8L | A | 223 | 23.3 | C | 176 | 20.2
3FM8 | C | 364 | 45.6 | B | 95 | 14.0
2W83 | A | 161 | 19.0 | CD | 129 | 17.7
2REX | B | 171 | 21.9 | A | 107 | 13.3
2VDU | D | 376 | 51.3 | F | 201 | 29.7

In the GSO optimization, 128 swarms have been simulated. Each swarm consisted
of the 600 glowworms randomly placed at the beginning of the simulation. 200
algorithm iterations have been performed for each swarm. In addition, the best
glowworm in each swarm has been minimized for 100 steps using the Powell
minimization method [24] at each GSO iteration. That resulted in $\sim 140000$
(i.e., $600*200+200*100)$ total evaluations of the target function.

The tensor train based methods used the same parameters for the generation of
the target function approximation with: the maximum rank of the tensor-train
expansion $rmax=3$; the number of grid points $ngrid=64$ along each of the 7
dimensions (3 for the center of masses translation and 4 for the rotation
quaternion); and the total number of iterations $nswp=10$. The parameters were
chosen such that the methods resulted in $\sim 120000$ target function
evaluations. To be consistent with the GSO algorithm, we performed the TT-
based minima search 128 times as well, starting with different random number
seeds. In addition, for the TTPY method, we performed local energy
minimization of the final structures for 100 steps using the same Powell
approach. In TTOPT, we performed additional local minimizations for 20 steps
for the best found points at the end of each iteration.

Figure 2: Scores evaluated with the DFIRE energy function for the five top
clusters of structures obtained for docking complexes with Terra Quantum’s
TTOPT (blue) and the TTPY (orange) and GSO (gray) approaches. Higher scores
indicate more optimal solutions.

All simulations were performed in the modified version of the LightDock
package that was changed in order to enable usage of the all three global
minima search methods.

##  III Results and Discussion

The simulations were carried out using 32 CPU cores on the Google Cloud
Platform. Afterwards, the simulations results were post-processed and
clusterized. For tensor-train based methods, the positions of the C$\alpha$
protein backbone atoms for the best ligand positions found in each of the 128
simulations were generated and clusterized using RMSD criteria of 4Å. For the
GSO approach, the five best found structures of the receptor-ligand complex
for each swarm were generated and further clusterized using the same 4Å RMSD
criteria for the C$\alpha$ atoms. The five best-scoring clusters were used for
comparisons of each complex. The results are shown in Fig. II.2. In the
majority of cases, TTOPT outperforms GSO quite significantly in terms of the
found global minima. It also outperforms the TTPY minimizer for most cases.
However, in complex 3Q87, we observed almost the same results for all three
methods. Inspection of the found geometries for the best-scoring complexes
shows that they are nearly identical as presented for the top cluster on Fig.
3a. But in the other studied systems, the tensor-train based approaches tend
to find better minima than GSO as shown for the 2REX system in Fig. 3b.

Figure 3: Overlap of the best structures for the docking ligand protein found
by employing the GSO (gray), TTOPT (blue) and TTPY (orange) algorithms. The
receptor protein is shown in red. For the 3Q87 panel (a), the complexes are
nearly identical. For the 2REX panel (b), the tensor-train based approaches
found better minima.

##  IV Conclusion

We have presented a novel protein–protein docking approach based on the
tensor-train decomposition of multi-dimensional scoring function landscapes.
Our own TTOPT algorithm was implemented in the LightDock package. It
demonstrated promising capabilities for sampling the translational and
rotational space of rigid body protein–protein docking simulations. We
demonstrated that our method outperforms LightDock’s native GSO algorithm as
well as the tensor-train based black-box global optimizer implemented in the
ttpy library. Due to the nature of our optimization approach, it could also be
employed for other chemical and biological modeling problems, such as the
folding of proteins and nucleic acids. More broadly, it could be used in
machine learning pipelines and for the optimization of hyperparameters [25],
for instance in the prediction of drug responses [26].

##  V Acknowledgments

We are grateful to Prof. V. Vinokur and Dr. Y. Alexeev from Argonne National
Laboratory for valuable discussions and suggestions.

## References

  * [1] David Eisenberg, Edward M Marcotte, Ioannis Xenarios, and Todd O Yeates.  Protein function in the post-genomic era.  Nature, 405(6788):823–826, 2000.
  * [2] Helen M Berman, John Westbrook, Zukang Feng, Gary Gilliland, Talapady N Bhat, Helge Weissig, Ilya N Shindyalov, and Philip E Bourne.  The protein data bank.  Nucleic acids research, 28(1):235–242, 2000.
  * [3] Xuan-Yu Meng, Hong-Xing Zhang, Mihaly Mezei, and Meng Cui.  Molecular docking: a powerful approach for structure-based drug discovery.  Current computer-aided drug design, 7(2):146–157, 2011.
  * [4] Jerome De Ruyck, Guillaume Brysbaert, Ralf Blossey, and Marc F Lensink.  Molecular docking as a popular tool in drug design, an in silico travel.  Advances and applications in bioinformatics and chemistry: AABC, 9:1, 2016.
  * [5] Nataraj S Pagadala, Khajamohiddin Syed, and Jack Tuszynski.  Software for molecular docking: a review.  Biophysical reviews, 9(2):91–102, 2017.
  * [6] Tammy Man-Kuang Cheng, Tom L Blundell, and Juan Fernandez-Recio.  pydock: electrostatics and desolvation for effective scoring of rigid-body protein–protein docking.  Proteins: Structure, Function, and Bioinformatics, 68(2):503–515, 2007.
  * [7] Xiaofan Li, Iain H Moal, and Paul A Bates.  Detection and refinement of encounter complexes for protein–protein docking: taking account of macromolecular crowding.  Proteins: Structure, Function, and Bioinformatics, 78(15):3189–3196, 2010.
  * [8] KN Krishnanand and Debasish Ghose.  Glowworm swarm optimization for simultaneous capture of multiple local optima of multimodal functions.  Swarm intelligence, 3(2):87–124, 2009.
  * [9] Brian Jiménez-García, Jorge Roel-Touris, Miguel Romero-Durana, Miquel Vidal, Daniel Jiménez-González, and Juan Fernández-Recio.  Lightdock: a new multi-scale approach to protein–protein docking.  Bioinformatics, 34(1):49–55, 2018.
  * [10] Ivan Oseledets.  Tensor-train decomposition.  SIAM J. Scientific Computing, 33:2295–2317, 01 2011.
  * [11] Konstantin Sozykin, Andrei Chertkov, Roman Schutski, Anh-Huy Phan, Andrzej Cichocki, and Ivan Oseledets.  Ttopt: A maximum volume quantized tensor train-based optimization and its application to reinforcement learning.  arXiv preprint arXiv:2205.00293, 2022.
  * [12] William Huggins, Piyush Patil, Bradley Mitchell, K Birgitta Whaley, and E Miles Stoudenmire.  Towards quantum machine learning with tensor networks.  Quantum Science and technology, 4(2):024001, 2019.
  * [13] Manuel S Rudolph, Jing Chen, Jacob Miller, Atithi Acharya, and Alejandro Perdomo-Ortiz.  Decomposition of matrix product states into shallow quantum circuits.  arXiv preprint arXiv:2209.00595, 2022.
  * [14] Mohammad Kordzanganeh, Markus Buchberger, Maxim Povolotskii, Wilhelm Fischer, Andrii Kurkin, Wilfrid Somogyi, Asel Sagingalieva, Markus Pflitsch, and Alexey Melnikov.  Benchmarking simulated and physical quantum processing units using quantum and hybrid algorithms.  arXiv preprint arXiv:2211.15631, 2022.
  * [15] Michael Perelshtein, Asel Sagingalieva, Karan Pinto, Vishal Shete, Alexey Pakhomchik, Artem Melnikov, Florian Neukart, Georg Gesek, Alexey Melnikov, and Valerii Vinokur.  Practical application-specific advantage through hybrid quantum computing.  arXiv preprint arXiv:2205.04858, 2022.
  * [16] Jorge Roel-Touris, Alexandre MJJ Bonvin, and Brian Jiménez-García.  Lightdock goes information-driven.  Bioinformatics, 36(3):950–952, 2020.
  * [17] Ivan Oseledets.  ttpy v.1.2.1, 2022.
  * [18] Dmitry Zheltkov and Eugene Tyrtyshnikov.  Global optimization based on TT-decomposition.  Russian Journal of Numerical Analysis and Mathematical Modelling, 35(4):247–261, 2020.
  * [19] Ivan Oseledets and Eugene Tyrtyshnikov.  TT-cross approximation for multidimensional arrays.  Linear Algebra and its Applications, 432(1):70–88, 2010.
  * [20] Sergei Goreinov, Ivan Oseledets, Dimitry Savostyanov, Eugene Tyrtyshnikov, and Nikolay Zamarashkin.  How to find a good submatrix.  In Matrix Methods: Theory, Algorithms And Applications: Dedicated to the Memory of Gene Golub, pages 247–256. World Scientific, 2010\.
  * [21] Micheline B Soley, Paul Bergold, and Victor S Batista.  Iterative power algorithm for global optimization with quantics tensor trains.  Journal of Chemical Theory and Computation, 17(6):3280–3291, 2021\.
  * [22] Marc F Lensink and Shoshana J Wodak.  Score_set: a capri benchmark for scoring protein complexes.  Proteins: Structure, Function, and Bioinformatics, 82(11):3163–3169, 2014.
  * [23] CHI Zhang, Song Liu, and Yaoqi Zhou.  Accurate and efficient loop selections by the dfire-based all-atom statistical potential.  Protein science, 13(2):391–399, 2004.
  * [24] Michael JD Powell.  An efficient method for finding the minimum of a function of several variables without calculating derivatives.  The computer journal, 7(2):155–162, 1964.
  * [25] Asel Sagingalieva, Andrii Kurkin, Artem Melnikov, Daniil Kuhmistrov, Michael Perelshtein, Alexey Melnikov, Andrea Skolik, and David Von Dollen.  Hyperparameter optimization of hybrid quantum neural networks for car classification.  arXiv preprint arXiv:2205.04878, 2022.
  * [26] Asel Sagingalieva, Mohammad Kordzanganeh, Nurbolat Kenbayev, Daria Kosichkina, Tatiana Tomashuk, and Alexey Melnikov.  Hybrid quantum neural network for drug response prediction.  arXiv preprint arXiv:2211.05777, 2022.


# KNODE-MPC: A Knowledge-based Data-driven Predictive Control Framework for
Aerial Robots

Kong Yao Chee∗, Tom Z. Jiahao∗ and M. Ani Hsieh The authors are with the GRASP
Laboratory, University of Pennsylvania, Philadelphia, PA 19104, USA.
{ckongyao, zjh<EMAIL_ADDRESS>Equal contribution.

###### Abstract

In this work, we consider the problem of deriving and incorporating accurate
dynamic models for model predictive control (MPC) with an application to
quadrotor control. MPC relies on precise dynamic models to achieve the desired
closed-loop performance. However, the presence of uncertainties in complex
systems and the environments they operate in poses a challenge in obtaining
sufficiently accurate representations of the system dynamics. In this work, we
make use of a deep learning tool, knowledge-based neural ordinary differential
equations (KNODE), to augment a model obtained from first principles. The
resulting hybrid model encompasses both a nominal first-principle model and a
neural network learnt from simulated or real-world experimental data. Using a
quadrotor, we benchmark our hybrid model against a state-of-the-art Gaussian
Process (GP) model and show that the hybrid model provides more accurate
predictions of the quadrotor dynamics and is able to generalize beyond the
training data. To improve closed-loop performance, the hybrid model is
integrated into a novel MPC framework, known as KNODE-MPC. Results show that
the integrated framework achieves 73% improvement in simulations and more than
14% in physical experiments, in terms of trajectory tracking performance.

##  I Introduction

The advent of model predictive control (MPC) has enabled control processes to
take advantage of the rapid advancement in optimization techniques and the
computational power of modern hardware [1]. The ability to optimize over
dynamic models and incorporate various constraints has bestowed model
predictive controllers with feasibility and stability guarantees for complex
control tasks [2]. However, identifying accurate dynamic models remains a
central challenge to MPC. Building accurate first-principle models often
requires domain expertise and deep physical insights into the systems. Even
for experts, modeling can be a daunting and tedious task when the system or
its operating environment is complex.

Recent advances in machine learning algorithms have enabled efficient
discovery of coherent patterns in complex data. Neural networks (NNs) have
been used for various tasks in computer vision, natural language processing,
and recommender systems. Most recently, a new family of neural networks –
neural ordinary differential equations (NODE), has been shown as an effective
tool for extracting dynamic models from data. NODE approximates a continuous-
depth neural network which is used to directly model differential equations.
Knowledge-based NODE (KNODE) was proposed to leverage NODE’s compatibility
with first principles. It combines first-principle models and NNs into hybrid
models to improve the open-loop prediction accuracy of nonlinear dynamic
models [3]. However, it is unclear whether such hybrid models could improve
the closed-loop performance of robotic systems in the real world.

In this work, we use KNODE to model the dynamics of a quadrotor. In
particular, we utilize the neural network in the hybrid model to account for
residual and uncertain dynamics within the system. This hybrid model is then
incorporated into a model predictive control framework, known as KNODE-MPC. A
schematic of the proposed framework is depicted in Fig. 1. We demonstrate that
the integrated framework not only provides better predictions of future
states, but also improves the closed-loop performance in both simulations and
real-world experiments.

Figure 1: Schematic of the KNODE-MPC framework, applied to a quadrotor system.
The NN accounts for residual or uncertain dynamics in the system, and is
combined with a model derived from first principles to form a KNODE model. The
KNODE model is then incorporated into a MPC framework, where optimal control
commands are generated and applied to the quadrotor.

##  II Related Work

In recent years, increasing data availability has fueled rapid development in
data-driven system identification. One body of work aims to extract the
governing equations of dynamic models directly from data [4, 5, 6]. A major
contribution uses sparse regression to reconstruct dynamic models from a
predetermined library of basis functions [4]. This method has been shown as an
effective alternative to system identification [7]. However, a fundamental
challenge to this method lies in coming up with the library of functions. If
the correct basis functions are missing from the library, the sparse
regression fails to identify the correct dynamic model. Such methods also do
not scale well to high-dimensional systems as the library of functions can
become unreasonably large. This prevents deployment on real-world systems,
which are very often high-dimensional.

On the other hand, neural networks have recently been shown as an effective
parameterization of dynamic models. Stochastic neural networks were used to
model residual physics in contact models [8], and subsequently ContactNets was
proposed to model frictional contact behaviors [9]. Recurrent neural networks
(RNNs) have been used to model the dynamics of high-dimensional flows [10] and
spatio-temporally chaotic systems [11]. NODE was developed for similar
purposes but differs from RNNs by explicitly modeling differential equations
from data [12]. Notably, NODE has enabled various deep learning tools to
facilitate system identification. For example, convolutional neural networks
have been used to reduce the dimensionality of high-dimensional systems [3].
In particular, KNODE has demonstrated neural network’s generality by modeling
a variety of systems with nonlinear and chaotic dynamics, and has also shown
NODE’s compatibility with first-principle dynamic models [3]. These
developments have opened up new avenues for system identification for MPC.

There are a number of works in the literature that integrate data-driven
models into a predictive control framework. One approach to synthesizing
models for MPC is through Gaussian processes (GP) regression. In [13], GP
regression is used to learn the residual errors between the true and nominal
dynamics of a race car, accounting for prior knowledge of the system. The
authors in [14] extend this framework to model residual dynamics in
quadrotors. However, these methods assume that the relationship between the
residual and true dynamics is known and this limits the nature of
uncertainties that these models are able to accommodate. Another drawback of
GP regression is computational complexity. This implies the need to select a
small number of training data points that best represent the dynamics. Since
the dynamics are either unknown or uncertain, it is difficult to select these
points in practice. Another increasingly popular approach that incorporates
learning into MPC is model-based reinforcement learning (MBRL). In [15],
neural networks are used to learn a dynamics model, which is then applied to a
MPC framework to solve iterative tasks. The authors in [16] use MBRL to design
low-level controllers using sensors on-board the quadrotor. Compared to the
approaches that use GP regression and our proposed framework, these methods do
not account for prior knowledge of the system dynamics and are therefore less
sample-efficient.

Our contributions in this work are three-fold. First, we employ KNODE to
develop a high-fidelity quadrotor model capable of capturing poorly understood
uncertainties and residual dynamics. The KNODE model combines prior knowledge
of the system dynamics with a neural ODE, and provides a better representation
of the quadrotor dynamics. Second, we show that the KNODE model improves the
accuracy of state predictions significantly, over both the nominal model, as
well as a state-of-the-art GP-based prediction model, which is similar to
those described in [13] and [14]. In particular, we compare our framework
against a more general variant of these GP models that does not impose any
structural assumptions between the system uncertainties and the true dynamics.
Third, the hybrid model is integrated into a novel model predictive control
framework, known as KNODE-MPC. Simulations and physical experiments are
conducted with KNODE-MPC to evaluate the improvement in closed-loop
performance. Empirical results show that the framework outperforms a nominal
model predictive controller with a first-principle model.

##  III Methodology

###  III-A Quadrotor Dynamics

In this work, we use a six degrees-of-freedom quadrotor system to illustrate
our proposed KNODE-MPC framework. The model used to describe this system is
similar to that in [17]. Denoting the mass and moments of inertia as $m$ and
$J:=\text{diag}(J_{xx},\,J_{yy},\,J_{zz})$, the dynamics of the quadrotor,
governed by its equations of motion, can be described as:

$\begin{bmatrix}\ddot{x}\\\ \ddot{y}\\\
\ddot{z}\end{bmatrix}:=\ddot{\mathbf{r}}=\frac{1}{m}\left(\begin{bmatrix}0\\\
0\\\ -mg\end{bmatrix}+R_{B}^{W}\begin{bmatrix}0\\\ 0\\\
u_{1}\end{bmatrix}\right),\\\ $ (1) $\begin{bmatrix}\dot{p}\\\ \dot{q}\\\
\dot{r}\end{bmatrix}:=\dot{\boldsymbol{\omega}}=J^{-1}\left(\mathbf{u}_{2}-\begin{bmatrix}p\\\
q\\\ r\end{bmatrix}\times J\begin{bmatrix}p\\\ q\\\ r\end{bmatrix}\right),$
(2)

where $R_{B}^{W}$ is the transformation matrix between the quadrotor body
frame and the world frame. $u_{1}$ is the summation of the motor forces, i.e.,
$u_{1}:=\sum_{i=1}^{4}F_{i}$ and $\mathbf{u}_{2}$ is the moment vector that is
related to the motor forces by

$\mathbf{u_{2}}=\begin{bmatrix}0&L&0&-L\\\ -L&0&L&0\\\
\gamma&-\gamma&\gamma&-\gamma\end{bmatrix}\begin{bmatrix}F_{1}\\\ F_{2}\\\
F_{3}\\\ F_{4}\end{bmatrix},$ (3)

where $L$ is the arm length of the quadrotor and $\gamma$ is the ratio between
the moment and thrust coefficient of the motors. The orientation in terms of
Euler angles velocities is related to the angular velocities by

$\boldsymbol{\omega}=\begin{bmatrix}cos(\theta)&0&-\cos(\phi)\sin(\theta)\\\
0&1&\sin(\phi)\\\
\sin(\theta)&0&\cos(\phi)\cos(\theta)\end{bmatrix}\begin{bmatrix}\dot{\phi}\\\
\dot{\theta}\\\ \dot{\psi}\end{bmatrix}.$ (4)

By defining the state as
$\mathbf{x}:=[\mathbf{r}\;\dot{\mathbf{r}}\;\phi\;\theta\;\psi\;\boldsymbol{\omega}]\in\mathbb{R}^{n}$,
we have the following compact representation of the quadrotor system,

$\dot{\mathbf{x}}=f(\mathbf{x},\mathbf{u}),$ (5)

with $\mathbf{u}=[u_{1}\;\mathbf{u}_{2}]^{T}\in\mathbb{R}^{m}$ and
$f:\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$ denotes a
mapping given by (1), (2) and (4). In this work, both quaternions and Euler
angles are used interchangeably to represent orientation.

###  III-B Knowledge-based Neural ODEs

Although the nominal quadrotor model in (5) represents the system with certain
fidelity, it is generally not sufficient for applications where a more
accurate model is required. For instance, the quadrotor experiences
aerodynamic and environmental disturbances in flight, and these are not
captured by the nominal model. In addition, system parameters such as the
moments of inertia are inherently difficult to measure or estimate. These
uncertainties cause deviations between the true dynamics of the system and
those given by the nominal model. The KNODE framework tackles this challenge
by learning these uncertainties and residual dynamics using a data-driven
approach. By combining prior knowledge of the system with the learned
dynamics, this hybrid model not only improves accuracy, but is also more
sample efficient.

Consider the following representation of the true system dynamics

$\dot{\mathbf{x}}=f_{t}(\mathbf{x},\mathbf{u},\,\Delta(\mathbf{x},\mathbf{u})),$
(6)

where
$f_{t}:\mathbb{R}^{n}\times\mathbb{R}^{m}\times\mathbb{R}^{p}\rightarrow\mathbb{R}^{n}$
denotes the true dynamics. The subscript $t$ denotes the true system.
$\Delta:\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{p}$ is a
mapping from the state-input space to the uncertainty space. Next, we
parameterize the uncertainty $\Delta$ using an NN with parameters
$\boldsymbol{\theta}$. This represents part of the true dynamics that are
absent in the nominal model. With this parameterization, the hybrid KNODE
model (subscript $h$) can be written as

$\dot{\mathbf{x}}=f_{h}(\mathbf{x},\mathbf{u},\,\Delta_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{u})),$
(7)

where $\Delta_{\boldsymbol{\theta}}$ is a neural network, and $f_{h}$ includes
both the nominal model and its coupling with the neural network. This
framework is general enough to allow coupling between different components of
the state and input, and it accommodates uncertain forces and moments readily.
Furthermore, compared to non-parametric methods like Gaussian process (GP)
regression, this parametric approach allows training on large amounts of data
while keeping the size of the uncertainty model constant, and can potentially
be computationally faster during inference. Hence, by modeling uncertainty in
the dynamics, this hybrid model gives a more accurate representation of the
true system.

###  III-C Model Predictive Control

An MPC framework uses an optimization-based approach to generate a sequence of
optimal control inputs in a receding horizon manner. In its formulation, a
model of the system, together with state and input constraints, are considered
within an optimization problem of the following form,

$\begin{split}\underset{\mathbf{u}}{\textnormal{minimize}}\quad\;&\sum_{i=1}^{N-1}\left(x_{i}^{T}Qx_{i}+u_{i}^{T}Ru_{i}\right)+x_{N}^{T}Px_{N}\\\
\text{subject
to}\quad\;&x_{i+1}=f_{d}(x_{i},u_{i}),\quad\forall\,i=0,\dotsc,N-1,\\\
\;&x_{i}\in\mathcal{X},\quad
u_{i}\in\mathcal{U},\quad\forall\,i=0,\dotsc,N-1,\\\
\;&x_{0}=x(k),\end{split}$ (8)

where $\mathcal{X}$ and $\mathcal{U}$ are sets defining the state and input
constraints and $x(k)$ is the initial state to the optimization problem at
time step $k$. $Q$, $R$ and $P$ are weighting matrices for the stage, input
and terminal costs respectively. At each time step, (8) is solved to obtain a
sequence of control inputs $\mathbf{u}:=[u_{0},\dotsc,u_{N-1}]^{T}$ and
$u_{0}$ is applied to the system. The discrete-time model
$f_{d}:\mathbb{R}^{n}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$ is derived
from a continuous-time counterpart using a numerical solver. An explicit
fourth-order Runge-Kutta method is applied to obtain the predicted states
$x_{i+1},\,i=0,\dotsc,N-1$. In the ideal case, the true dynamics in (6) are
discretized, but since the uncertainty is unknown, a discrete-time version of
the nominal model described in (5) is commonly used in this framework.

There exist robust and stochastic MPC methods in the literature that account
for uncertainty by assuming that the uncertainty is either bounded [18, 19] or
characterized by a probability distribution [20, 21]. However, these
approaches require additional assumptions on the nature of the uncertainty,
which can be challenging to ascertain in practice. In most cases, the closed-
loop performance is sensitive to the accuracy of the model and hence, using
the KNODE model in (7) is expected to yield better closed-loop system
performance than a nominal model. This motivates the KNODE-MPC framework,
where a hybrid knowledge-based NODE model is integrated into the MPC
framework.

###  III-D KNODE Training

While the KNODE framework is applicable to the general case shown in (6), we
shall focus on a more practical scenario. Assuming that the uncertainty is
separable from the nominal dynamics, the true dynamics can then be expressed
as

$\dot{\mathbf{x}}=f_{t}(\mathbf{x},\mathbf{u},\,\Delta(\mathbf{x},\mathbf{u}))=f(\mathbf{x},\mathbf{u})+\Delta(\mathbf{x},\mathbf{u}),$
(9)

with the first component given by the nominal dynamics in (5). Following the
parameterization described in III-B, the hybrid model that consists of both
the nominal model and neural network is given by

$\dot{\mathbf{x}}=f_{h}(\mathbf{x},\mathbf{u},\Delta_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{u}))=f(\mathbf{x},\mathbf{u})+\Delta_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{u}),$
(10)

where $\Delta_{\boldsymbol{\theta}}$ denotes the neural network. Training of
$\Delta_{\boldsymbol{\theta}}$ is done with the trajectory data collected from
closed-loop simulations or experiments. The training process in our work
differs from [3] in that control inputs in the training data need to be
injected into the model to make the one-step-ahead predictions. Similar to
[3], we first divide up trajectory data from each time step, and then use the
state at each of these time steps as initial conditions to make one-step-ahead
predictions using (10). Specifically, we consider a set of state observations
and control inputs
$\mathcal{O}=[(\mathbf{x}(t_{1}),\mathbf{u}(t_{1})),(\mathbf{x}(t_{2}),\mathbf{u}(t_{2})),\cdots,(\mathbf{x}(t_{N}),\mathbf{u}(t_{N}))]^{T}$,
sampled at the times $T=\\{t_{1},t_{2},\cdots,t_{N}\\}$. For a time step
$t_{i}$, we compute the one-step-ahead prediction given by

$\hat{\mathbf{x}}(t_{i+1})=\mathbf{x}(t_{i})+\int^{t_{i+1}}_{t_{i}}f_{h}(\mathbf{x}(t_{i}),\mathbf{u}(t_{i}))dt,$
(11)

where $\hat{\mathbf{x}}(t_{i+1})$ denotes the state prediction at time
$t_{i+1}$. The learning framework in this work performs the integration in
(11) numerically using the explicit fourth order Runge-Kutta method. Next, to
compute the mean squared error (MSE) between the predictions and the ground
truth, we define a loss function as

$L(\boldsymbol{\theta}):=\frac{1}{N-1}\sum^{N}_{i=2}\left(\hat{\mathbf{x}}(t_{i})-\mathbf{x}(t_{i})\right),$
(12)

and the optimization problem is therefore given by

$\displaystyle\min_{\boldsymbol{\theta}}$ $\displaystyle
L(\boldsymbol{\theta})$ (13) s.t.
$\displaystyle\dot{\mathbf{x}}=f_{h}(\mathbf{x},\mathbf{u},\Delta_{\boldsymbol{\theta}}(\mathbf{x},\mathbf{u})),$

where the constraint comprises of the dynamic model. To train the hybrid model
on multiple trajectories, we can sum up the losses computed on each
trajectory. Similar to [3, 12], we solve the optimization task in (13) using
the adjoint sensitivity method which propagates gradients from the loss
function to the neural network parameters. The parameters are then iteratively
updated based on the gradients using readily available optimizers like Adam
[22]. The adjoint sensitivity method has been noted as a memory-efficient
alternative to the conventional backpropagation [12], and its proof can be
found in [23].

###  III-E Evaluation and Verification

In our experimental design and evaluation process, we set out to answer the
following questions about our proposed KNODE-MPC framework. (a) How accurate
are the state predictions given by the KNODE model, as compared to state-of-
the-art data-driven models in [13] and [14]? (b) How well does the KNODE model
generalize beyond training data? (c) How much improvement does the KNODE-MPC
framework provide, in terms of closed-loop trajectory tracking, compared to a
nominal MPC framework?

To evaluate the accuracy of state predictions and to verify the generalization
capability of the KNODE model, a generalized variant of the model in [14] is
used as a benchmark. In [14], it is assumed that the relationship between the
uncertainties and equations of motion are known and this is incorporated as
user-defined selection matrices within the model (see Section III.C of [14]).
The model also assumes that the uncertainties are decoupled and only allows
for a mapping between the body velocities and acceleration disturbances (see
Section III.E of [14]), without accounting for any rotational disturbances.
These assumptions are often difficult to ascertain in practice.

In this work, we implement a general version of this GP model that removes
these assumptions, which can be compactly represented as

$\dot{\mathbf{x}}_{\textbf{gp}}=f_{\textbf{gp}}(\mathbf{x},\mathbf{u})=f(\mathbf{x},\mathbf{u})+\boldsymbol{\mu}\left(\begin{bmatrix}\mathbf{x}\\\
\mathbf{u}\end{bmatrix}\right),$ (14)

where $f(\mathbf{x},\mathbf{u})$ represents the nominal model in (5) and
$\boldsymbol{\mu}$ denotes the mean posterior of a GP. A nominal model that
does not account for uncertainty is also included in our tests as a baseline.
To demonstrate the effectiveness of both data-driven models, we incorporate
nonlinear components into the true dynamics and compare the prediction
accuracy across three models; the nominal, KNODE and GP models from (5), (10)
and (14) respectively. Quantitatively, we want to measure the spatial
similarity between the true and predicted trajectories and this can be
evaluated using a distance metric generated by the dynamic time warping (DTW)
algorithm [24]. Next, for evaluation of the trajectory tracking performance of
KNODE-MPC, we compare it against a nominal MPC framework, both in simulations
and in physical experiments. The nominal MPC framework uses a model that does
not consider uncertainty such as aerodynamic forces or environmental
disturbances in its formulation. The closed-loop trajectories obtained from
both frameworks are compared against the planned trajectory using the same
distance metric.

##  IV Experiments and Results

###  IV-A Setup

Simulations: A nominal quadrotor model is constructed using the equations of
motion given in (5). An explicit 5${}^{\text{th}}$ order Runge-Kutta method
(RK45) with a sampling interval of 2 milliseconds is used for numerical
integration to generate dynamic responses of the quadrotor. We assume that the
model predictive controller has access to perfect measurements of the
quadrotor dynamics. The MPC architecture is implemented in CasADi [25] and the
optimization routine generates optimal thrust and moment commands. These
commands act as inputs to the quadrotor model, which simulates the closed-loop
responses. Two classes of trajectories are considered in the simulations;
circular and lemniscate trajectories. Both the GP and KNODE models are trained
on circular trajectories of radii 3m and 6m, and a circular trajectory of
radius 4m is used as validation data. To model uncertainty within the system,
nonlinear aerodynamic drag effects from the fuselage and rotors are
incorporated into the true quadrotor dynamics. These uncertain dynamics are
absent in the nominal model and are to be learnt by the GP and KNODE models.

The KNODE model uses a two layer neural network architecture with 64 and 16
hidden units each. A hyperbolic tangent function is used as the activation
function. Each of the two trajectories in the training data has a duration of
8s and 4000 data points. The GP model is trained using the scikit-learn
package [26] and uses a kernel consisting of the product of constant and
radial basis function (RBF) kernels. Eighty data points, sampled at regular
intervals, are used for training. The inputs for the training data are data
samples containing the state-input vector $[\mathbf{x}\;\mathbf{u}]^{T}$,
while the training labels are the differences between the true and nominal
state derivatives, i.e.,
$\Delta\dot{x}:=f_{\textbf{gp}}(\mathbf{x},\mathbf{u})-f(\mathbf{x},\mathbf{u})$.
With the trained models, simulations are conducted along both types of
trajectories to evaluate the prediction accuracy of the KNODE model against
the nominal and GP models. To assess generalization performance, the models
are tested on trajectories beyond the training data. These include circular
trajectories of different radii, as well as lemniscate trajectories. Lastly,
to test closed-loop performance, the KNODE-MPC framework is simulated on
circular trajectories of various radii and compared against a nominal model
predictive controller.

Figure 2: Schematic of the experimental setup: The Crazyflie experimental
platform operating with a Vicon motion capture system. KNODE-MPC publishes
control commands through the Crazyflie ROS server and into the commander
module on-board Crazyflie.

Physical Experiments: We use the open-source Crazyflie 2.1 quadrotor [27] as
the experimental platform. An image of a Crazyflie is shown in Fig. 2. The
Crazyflie, together with motion trackers, weighs 32 g and has a size of 9 cm2.
The software architecture is developed using the CrazyROS wrapper [28]. A
laptop running on Intel i7 CPU acts as the base station and communication with
the Crazyflie is done via the Crazyradio PA at a rate of 250 Hz. To obtain
pose information of the quadrotor, we use a Vicon motion capture system which
communicates with the base station via a Vicon bridge and publishes at an
approximate rate of 88 Hz. Linear velocities are estimated from the positions
obtained from the Vicon, while angular velocities are measured from the
gyroscope sensors on-board the quadrotor. This experimental setup is depicted
in Fig. 2.

To train and validate the KNODE model, data is collected by flying the
Crazyflie with the nominal MPC framework along circular trajectories of radii
0.5 and 1m, each for a duration of 76 seconds and 40000 data points. We
preprocess the velocity data using a fifth order Butterworth low pass filter
before training. Similar to the simulation setup, both the nominal and KNODE-
MPC frameworks are implemented in CasADi and run on the base station. The
KNODE in physical experiments uses a one-layer neural network with 32 hidden
units and the hyperbolic tangent activation function. The frameworks generate
three-dimensional acceleration commands and run on top of a geometric
controller [17], as well as attitude and thrust controllers within the
Crazyflie firmware. For each of the two frameworks, we conduct five test runs
along circular trajectories of radii 0.5m and 1.25m to evaluate closed-loop
trajectory tracking performance.

We refer the reader to the video in the supplemental materials for a better
understanding of the physical experiments.

###  IV-B Simulation Results

To illustrate the accuracy of state predictions, the 3D DTW position errors
between the predicted trajectories and the true trajectory are plotted in Fig.
3. Both the KNODE and GP models are able to account for the nonlinear
aerodynamic drag effects included in the true dynamics and provide accurate
predictions, outperforming the nominal model by more than a factor of 10.
Moreover, the accuracy of the KNODE model is observed to be substantially
higher than that of the GP model across all trajectories. In particular, by
contrasting the DTW errors of both models, an average improvement of 82% was
observed for the KNODE model under circular trajectories. It is also noted
that although the GP model provide reasonable predictions on the training data
(circular trajectories of radii 3m and 6m), there are cases in which it does
not generalize. For instance, the accuracy of the GP model under a circular
trajectory of radius 2m is comparable to that of a nominal model, as shown in
Fig. 3. On the other hand, the accuracy of the KNODE model is consistent
across all trajectories, which implies that the KNODE model is capable of both
interpolating within and extrapolating outside of its training domain. For
lemniscate trajectories, KNODE outperforms the GP model by an average of 63%
and this is shown in Fig. 4. These results for lemniscate trajectories
ascertain the generalization capability of the KNODE model.

Figure 3: DTW prediction errors for the nominal, GP and KNODE models, tested
with circular trajectories of different radii. Figure 4: DTW prediction errors
for the nominal, GP and KNODE models, tested with lemniscate trajectories of
different radii.

Closed-loop simulation results for both the nominal and KNODE-MPC frameworks
are plotted in Fig. 5. Compared to nominal MPC, KNODE-MPC reduces the 3D
trajectory tracking DTW errors by 73% on average. This demonstrates the
effectiveness of the integrated framework in the presence of uncertain
dynamics. Fig. 6 provides a comparison between the closed-loop trajectories
under the two frameworks. Since the model in the nominal MPC framework does
not account for residual or uncertain dynamics, the quadrotor trajectory
deviates from the planned trajectory, with position errors accumulating over
time. On the other hand, the KNODE model compensates for the residual dynamics
and allows the quadrotor to follow the planned trajectory more closely.

Figure 5: Trajectory tracking errors for the nominal and KNODE models, under
the nominal and KNODE MPC frameworks, along circular trajectories of different
radii. Figure 6: 3D trajectories for the nominal and KNODE models, under the
nominal and KNODE MPC frameworks, under a circular trajectory of radius = 2m.

###  IV-C Physical Experimental Results

The system and planned trajectories of the quadrotor with both nominal MPC and
KNODE-MPC under a circular trajectory of radius 0.5m are shown in Fig 7. The
boundaries of the red annulus indicate the maximum deviation of the system
trajectories under the nominal MPC framework. The same annulus is overlaid on
the right subplot and it provides a qualitative assessment of the improvement
in accuracy with the KNODE-MPC framework. Notably, the trajectories under
KNODE-MPC do not exceed the boundaries of the annulus, which implies there are
smaller deviations from the planned trajectory. By taking the average across
five runs, the nominal MPC yields an DTW x-y position error of 4.69, while an
error of 2.54 is obtained with KNODE-MPC. These indicate an improvement of
over 45%. In the z axis, the average DTW errors between the frameworks are
comparable, valued at 11.59 and 12.11 respectively.

Next, the generalization ability of KNODE-MPC is verified using trajectories
beyond its training dataset. In this scenario, the planned trajectory has a
radius of 1.25m, which differs from training trajectories, which have radii of
0.5 and 1m. Fig. 8 depicts the system and planned trajectories under both
frameworks. Here, the boundaries of the red annulus also indicate the maximum
deviation of the system trajectories under the nominal MPC and it is overlaid
on the right subplot. The trajectories generated by KNODE-MPC lie strictly
within the red annulus, which implies the closed-loop tracking performance of
KNODE-MPC is more accurate than nominal MPC. Quantitatively, the average DTW
x-y position error across five test runs for both the nominal and KNODE-MPC
are 10.3 and 8.81 respectively, which yields a 14.5% improvement in trajectory
tracking for KNODE-MPC. The differences in performance improvement between
simulation and physical experimental results can be attributed to the
different control architectures between the two setups. In particular, the
KNODE-MPC framework in simulations generate control commands to actuate the
quadrotor directly, while in physical experiments, it runs on top of the
geometric controller and low-level controllers within the firmware, and it
does not have direct control over the quadrotor dynamics.

Figure 7: Top view: Quadrotor trajectories under both nominal MPC and KNODE-
MPC, with a planned trajectory of radius 0.5m. Figure 8: Top view: Quadrotor
trajectories under both nominal MPC and KNODE-MPC, with a planned trajectory
of radius 1.25m.

##  V Conclusion and Future Work

In this work, we presented the KNODE-MPC framework with an application to
quadrotor control. The knowledge-based, data-driven model, KNODE, provides a
more accurate representation of the system dynamics, as compared to a state-
of-the-art GP model. By incorporating KNODE into the MPC framework, the
predictive controller generates optimal control commands that yield better
closed-loop performance. This framework is tested extensively in simulations
and physical experiments. Results show significant improvements in terms of
state predictions, as well as closed-loop trajectory tracking. Given its
flexibility, the framework can be applied to robotic systems operating in
complex and uncertain environments such as marine robots or robotic teams. In
the near future, we would like to extend this framework to an online learning
setting, where the KNODE model can be updated and refined in real-time, which
can further improve the performance of closed-loop control tasks.

## References

  * [1] S. J. Qin and T. A. Badgwell, “A survey of industrial model predictive control technology,” _Control Eng. Pract._ , vol. 11, no. 7, pp. 733–764, 2003.
  * [2] F. Borrelli, A. Bemporad, and M. Morari, _Predictive control for linear and hybrid systems_. Cambridge Univ. Press, 2017.
  * [3] T. Z. Jiahao, M. A. Hsieh, and E. Forgoston, “Knowledge-based learning of nonlinear dynamics and chaos,” 2021.
  * [4] S. L. Brunton, J. L. Proctor, and J. N. Kutz, “Discovering governing equations from data by sparse identification of nonlinear dynamical systems,” _Proc. of the Nat. Acad. of Sciences_ , vol. 113, no. 15, pp. 3932–3937, 2016\.
  * [5] A. A. R. AlMomani, J. Sun, and E. Bollt, “How entropic regression beats the outliers problem in nonlinear system identification,” _Chaos_ , vol. 30, no. 1, p. 013107, 2020.
  * [6] K. Champion, B. Lusch, J. N. Kutz, and S. L. Brunton, “Data-driven discovery of coordinates and governing equations,” _Proc. of the Nat. Acad. of Sciences_ , vol. 116, no. 45, pp. 22 445–22 451, 2019. [Online]. Available: https://www.pnas.org/content/116/45/22445
  * [7] E. Kaiser, J. N. Kutz, and S. L. Brunton, “Sparse identification of nonlinear dynamics for model predictive control in the low-data limit,” _Proc. of the Roy. Soc. A: Math., Phys. and Eng. Sciences_ , vol. 474, no. 2219, p. 20180335, Nov 2018. [Online]. Available: http://dx.doi.org/10.1098/rspa.2018.0335
  * [8] A. Ajay, J. Wu, N. Fazeli, M. Bauzá, L. Kaelbling, J. Tenenbaum, and A. Rodriguez, “Augmenting physical simulators with stochastic neural networks: Case study of planar pushing and bouncing,” _2018 IEEE/RSJ Int. Conf. on Intell. Robots and Systems (IROS)_ , pp. 3066–3073, 2018.
  * [9] S. Pfrommer, M. Halm, and M. Posa, “Contactnets: Learning discontinuous contact dynamics with smooth, implicit representations,” 2020.
  * [10] M. Qraitem, D. Kularatne, E. Forgoston, and M. A. Hsieh, “Bridging the gap: Machine learning to resolve improperly modeled dynamics,” _Physica D: Nonlinear Phenomena_ , vol. 414, 12 2020.
  * [11] A. Wikner, J. Pathak, B. Hunt, M. Girvan, T. Arcomano, I. Szunyogh, A. Pomerance, and E. Ott, “Combining machine learning with knowledge-based modeling for scalable forecasting and subgrid-scale closure of large, complex, spatiotemporal systems,” _Chaos_ , vol. 30, p. 053111, 05 2020.
  * [12] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud, “Neural ordinary differential equations,” in _Advances in Neural Inf. Process. Syst._ Curran Associates, Inc., 2018, pp. 6571–6583. [Online]. Available: http://papers.nips.cc/paper/7892-neural-ordinary-differential-equations.pdf
  * [13] J. Kabzan, L. Hewing, A. Liniger, and M. N. Zeilinger, “Learning-based model predictive control for autonomous racing,” _IEEE Robot. Autom. Lett._ , vol. 4, no. 4, pp. 3363–3370, 2019.
  * [14] G. Torrente, E. Kaufmann, P. Föhn, and D. Scaramuzza, “Data-driven mpc for quadrotors,” _IEEE Robot. Autom. Lett._ , vol. 6, no. 2, pp. 3769–3776, 2021.
  * [15] G. Williams, N. Wagener, B. Goldfain, P. Drews, J. M. Rehg, B. Boots, and E. A. Theodorou, “Information theoretic mpc for model-based reinforcement learning,” in _2017 IEEE Int. Conf. on Robot. and Automat._ IEEE, 2017, pp. 1714–1721.
  * [16] N. O. Lambert, D. S. Drew, J. Yaconelli, S. Levine, R. Calandra, and K. S. J. Pister, “Low-level control of a quadrotor with deep model-based reinforcement learning,” _IEEE Robot. Autom. Lett._ , vol. 4, no. 4, pp. 4224–4230, 2019.
  * [17] D. Mellinger and V. Kumar, “Minimum snap trajectory generation and control for quadrotors,” in _2011 IEEE Int. Conf. on Robot. and Automat._ , 2011, pp. 2520–2525.
  * [18] A. Bemporad and M. Morari, “Robust model predictive control: A survey,” in _Robustness in identification and control_. Springer, 1999, pp. 207–226.
  * [19] D. Q. Mayne, E. C. Kerrigan, E. Van Wyk, and P. Falugi, “Tube-based robust nonlinear model predictive control,” _Int. J. of Robust and Nonlinear Control_ , vol. 21, no. 11, pp. 1341–1353, 2011.
  * [20] A. Mesbah, “Stochastic model predictive control: An overview and perspectives for future research,” _IEEE Control Syst. Mag._ , vol. 36, no. 6, pp. 30–44, 2016.
  * [21] A. Mesbah, S. Streif, R. Findeisen, and R. D. Braatz, “Stochastic nonlinear model predictive control with probabilistic constraints,” in _2014 Amer. Control Conf._ , 2014, pp. 2413–2419.
  * [22] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _CoRR_ , vol. abs/1412.6980, 2015.
  * [23] Y. Cao, S. Li, L. Petzold, and R. Serban, “Adjoint sensitivity analysis for differential-algebraic equations: The adjoint DAE system and its numerical solution,” _SIAM J. Sci. Comput._ , vol. 24, pp. 1076–1089, 2003.
  * [24] H. Sakoe and S. Chiba, “Dynamic programming algorithm optimization for spoken word recognition,” _IEEE Trans. on Acoust., Speech, and Signal Process._ , vol. 26, no. 1, pp. 43–49, 1978.
  * [25] J. A. E. Andersson, J. Gillis, G. Horn, J. B. Rawlings, and M. Diehl, “CasADi – A software framework for nonlinear optimization and optimal control,” _Math. Program. Computation_ , vol. 11, no. 1, pp. 1–36, 2019\.
  * [26] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, “Scikit-learn: Machine learning in Python,” _J. of Mach. Learn. Res._ , vol. 12, pp. 2825–2830, 2011.
  * [27] Bitcraze, “Crazyflie 2.1.” [Online]. Available: https://www.bitcraze.io/products/crazyflie-2-1/
  * [28] W. Hönig and N. Ayanian, _Flying Multiple UAVs Using ROS_. Springer Int. Publishing, 2017, pp. 83–118.