a b d e k l H L p r s S F t w th div tr Re Im AA SS Si gg ss

# Gradient flow encoding with distance optimization adaptive step size

Kyriakos Flouris
Computer Vision Laboratory, ETH Zurich
<EMAIL_ADDRESS>&Anna Volokitin
Computer Vision Laboratory, ETH Zurich
<EMAIL_ADDRESS>Gustav Bredell
Computer Vision Laboratory, ETH Zurich
<EMAIL_ADDRESS>Ender Konukoglu
Computer Vision Laboratory, ETH Zurich
<EMAIL_ADDRESS>
###### Abstract

The autoencoder model uses an encoder to map data samples to a lower
dimensional latent space and then a decoder to map the latent space
representations back to the data space. Implicitly, it relies on the encoder
to approximate the inverse of the decoder network, so that samples can be
mapped to and back from the latent space faithfully. This approximation may
lead to sub-optimal latent space representations. In this work, we investigate
a decoder-only method that uses gradient flow to encode data samples in the
latent space. The gradient flow is defined based on a given decoder and aims
to find the optimal latent space representation for any given sample through
optimisation, eliminating the need of an approximate inversion through an
encoder. Implementing gradient flow through ordinary differential equations
(ODE), we leverage the adjoint method to train a given decoder. We further
show empirically that the costly integrals in the adjoint method may not be
entirely necessary. Additionally, we propose a $2^{nd}$ order ODE variant to
the method, which approximates Nesterov’s accelerated gradient descent, with
faster convergence per iteration. Commonly used ODE solvers can be quite
sensitive to the integration step-size depending on the stiffness of the ODE.
To overcome the sensitivity for gradient flow encoding, we use an adaptive
solver that prioritises minimising loss at each integration step. We assess
the proposed method in comparison to the autoencoding model. In our
experiments, GFE showed a much higher data-efficiency than the autoencoding
model, which can be crucial for data scarce applications.

##  1 Introduction

Auto-encoders are wildly successful as artificial neural network architectures
for unsupervised learning [10]. The appeal is centred on learning a lower
dimensional representation of the input data allowing extremely efficient
computation. The idea is simple as it is brilliant, provided there is
correlation between the input data the latent dimension can be leveraged to
output a model of the input. Nevertheless, the encoding process is semi-
arbitrary as there is no direct learning of the encoder, rather than the
decoder is optimised to rectify the encoding process. This sub-optimal latent
space representations can lead to inefficient learning, [6]. In part to
compensate for this learning process extensive work has been done to
regularise the latent space more directly, [8].

Auto-encoders depend on the encoder to approximate the inverse of the decoder
network, this approximate inversion requires to be learned with additional
parameters and therefore may get in the way of learning with fewer data and
adversely affect the latent space structure. Some research in flow models,
[4], can resolve this issue by using invertible maps however they are
constrained to equi-dimensional latent representations. What if we can
eliminate the need of an encoder neural network altogether, retain the
advantage of lower-dimensional latent space, while obtaining a directly
optimised representation? This may allow us to create a model that can map
images to lower dimensional latent space and reconstruct faithfully with fewer
images and less iterations.

In this work we propose a novel encoding method namely a gradient flow
encoding (GFE). This decoder-only method at each training step primarily
determines the optimal latent space representation for each sample via a
gradient flow, namely, an ordinary differential equation (ODE) solver. The
decoder is updated as usual by minimising the loss between input image and its
reconstruction retrieved from the optimal latent space representation for the
image. The method, albeit being slower, is considerably superior in data
efficiency and ultimately obtains better reconstructions with very few number
of training samples when compared to a conventional auto-encoder (AE). A side-
advantage of GFE is ’halving’ the size of the AE neural network, no need for
an encoder.

Using a gradient flow for encoding can be a computationally challenging task.
Traditionally, ODE solvers with adaptive step size are propagating the
variables in such a manner as to minimise the error of integration. These
solvers are crucial for the accuracy of the ODE solution. However, they also
impede training of ODE based neural networks if the underlying ODE defined by
the network becomes stiff, [3, 5]. Adaptive step size becomes too small and
integration takes long. When using GFE at each training step, using an ODE
solver with adaptive step size can be a debilitating factor.

For integrating gradient flow in an encoding scheme, we observe that the exact
path of the gradient flow may be less important than the final point of
convergence. Therefore, adapting step size to minimise error in integration is
not necessarily the best approach. Using a fixed step size is also not an
option as it may result in poor convergence and stability issues in training
the decoder network. Here, we develop an adaptive solver for the gradient
flow, which adapts each optimisation (’time’) step such as to directly
minimise the distance, i.e. loss between input and output. Furthermore, this
adaptive minimise distance (AMD) solver is modified to include a loss
convergence assertion to improve performance. The AMD method can possibly make
gradient flow a computationally feasible module to be used in neural network
architectures.

Here, an adjoint method, [2], as also used in other recent works [3], is
implemented to properly optimise the latent space and decoder of the GFE.
Consequently, for efficiency considerations we show that a full adjoint method
is not necessarily needed and an approximation is utilised. Furthermore, we
present the implementation of a Nesterov $2^{nd}$ order ODE solver with
accelerated convergence per training data size. Ultimately, the approximate
GFE utilising AMD (GFE-amd) is employed for testing and comparison with a
traditional AE solver.

##  2 Relevant works

DeepSDF: Learning Continuous Signed Distance Functions for Shape
Representation [6]: The authors replace the auto-encoder network with an auto-
decoder where, a similar to here, latent vector $z$ is introduced. $z$
represents the encoding of a desired shape. They map this latent vector to a
3D shape represented by a continuous signed Distance Function. I.e. For shape
$i$ with function $f_{\t{$}}andcoordinate$\mathbf{x}$,$f_\t{(}z_i,\mathbf{x})
\approx
SDF^i(\mathbf{x})$.Byconditioningofthenetworkoutputonalatentvector,theymodelmultipleSDFswithasingleneuralnetwork.\par\par\par\par\par\par\par\par\par\par\par\@@numbered@section{section}{toc}{Method}
An auto-encoder funnels an input $y$ into a lower dimensional latent space
representation $z$ using an encoder network $E$ and reconstructs it back using
a decoder network $D$. Here, the encoder $E$ and decoder $D$ networks can be
thought as approximate inverses of each other. During training, each sample is
mapped to the latent space using $E$, mapped back using $D$ and the
reconstruction error is minimised with respect to the parameters of both of
the networks. \par At any point in this training process, given the decoder an
optimal latent space representation for each sample can be defined as the
$z^{*}$ that minimises the reconstruction error. The encoder however, does not
map the sample to that optimal $z^{*}$. So, instead of trying to get the
reconstruction of $z^{*}$ to get closer to the sample, its parameters are
updated to get the reconstruction of $E(z)$ to get closer to the sample. This
may not be the most efficient use of data samples. Alternatively, one can
determine $z^{*}$ for each sample and update the decoder's parameters
according to this optimal latent representation. \par\par Determination of
$z^{*}$ for each sample $y$ can be formulated as an optimisation problem
\begin{equation*}z^{*}=\arg_{z}\min l(y,D(z,\th)),\end{equation*} where
$\theta$ represent the parameters of the decoder network and $l(\cdot,\cdot)$
is a distance function between the sample and its reconstruction by the
decoder, which can be defined based on the application. One obvious form can
be the $L_{2}$ distance $\lVert D(z,\th)-y\rVert^{2}_{2}$. The optimisation
can be achieved by a gradient decent minimisation. In order to integrate the
minimisation in the training of the decoder network, a continuous gradient
decent algorithm is implemented via a solution to an ordinary differential
equation $dz/dt=-\alpha(t)\nabla_{z}l(y,D(z(t),\theta))$, where time $t$ is
the continuous parameter describing the extent of the minimisation and
$\alpha(t)$ is a scaling factor that can vary with time. When the extremum is
reached the $z$ comes to a steady state. In practice we compute the optimal
$z^{*}$ by integrating the ODE
\begin{equation*}z^{*}=z(\t{)}=\int_{0}^{\t{}}-\alpha(t)\nabla_{z}l(y,D(z(t),\th))dt,\
z(0)=z_{0},\end{equation*} where $z_{0}$ is the initialisation of the
optimisation, which is set as 0 vector in our experiments. Consequently to the
minimisation of $z\rightarrow z^{*}\equiv z(t=\t{)}$ for a given $D$ ('forward
model'), the decoder is trained with a total loss function for a given
training set $M$ ('backward'),
\begin{equation}\mathcal{L}(\theta)=\sum_{m=1}^{M}l(y_{m},D(z_{m}^{*},\theta)).\end{equation}
At each iteration while searching for the $\arg_{\th}\min\mathcal{L}(\theta)$,
a new $z^{*}$ is recalculated for each sample.
$\arg_{\th}\min\mathcal{L}(\theta)$ is computed via the adjoint method as
explained in Section \ref{sec:adjoint}.
\par\par\par\@@numbered@section{subsection}{toc}{The adjoint method for the
gradient flow} As described above, after finding $z^{*}\equiv
z(\t{)}=\arg_{z}\min l(y,D(z,\th))$ the total loss is minimised with respect
to the model parameters. The dependence of $z^{*}$ to $\th$ creates an
additional dependence of $\mathcal{L}(\th)$ to $\th$ via $z^{*}$. For
simplicity, let us consider the cost of only one sample $y$, effectively
$l(y,D(z^{*},\th))$. We will compute the total derivative
$d_{\th}l(y,D(z^{*},\th))$ for this sample and the derivative of the total
cost for a batch of samples can be computed as the sum of the sample
derivatives in the batch. The total derivative $d_{\th}l(y,D(z^{*},\th))$ is
computed as
\begin{equation*}d_{\theta}l(y,D(z^{*},\th))=\partial_{\theta}l(y,D(z^{*},\th))+\partial_{z^{*}}l(y,D(z^{*},\th))\partial_{\theta}z^{*}.\end{equation*}
The derivative $\partial_{z^{*}}l(y,D(z^{*},\th))d_{\theta}z^{*}$ can be
computed using the adjoint method and leads to the following set of equations
\@@amsalign dz/dt=-\alpha(t)\nabla_{z}l(y,D(z(t),\th)),\text{ with }z(0)=0\\\
d\l/dt=-\alpha(t)\l^{T}\nabla^{2}_{z}l(y,D(z(t),\th)),\text{ with
}\l(\tau)=-\nabla_{z}l(y,D(z(\t{)},\th))\\\
d_{\th}l(y,D(z^{*},\th))=\partial_{\th}l(y,D(z^{*},\th))-\int_{0}^{\t{}}\alpha(t)\lambda^{T}\partial_{\th}\nabla_{z}l(z(t),\th)dt,
where we used $z^{*}=z(\t{)}$. Equations [\ref{eq:adjoint1}-\ref{eq:adjoint3}]
define the so called adjoint method for gradient flow optimisation of the
loss. Due to the cost of solving all three equations we empirically find that
for this work sufficient and efficient optimisation can be accomplished by
ignoring the integral (``adjoint function'') part of the method.
Theoretically, this is equivalent to ignoring the higher order term of the
total differential
$d_{\th}l(y,D(z^{*},\th))=\partial_{\th}l(y,D(z^{*},\th))+\partial_{z}l(y,D(z^{*},\th))\partial_{\th}z\approx\partial_{\th}l(y,D(z^{*},\th))$.
Reducing Equations [\ref{eq:adjoint1}-\ref{eq:adjoint3}] to \@@amsalign
dz/dt=-\alpha(t)\nabla_{z}l(y,D(z(t),\th)),\text{ with }z(0)=0\\\
d_{\th}\mathcal{L}=\partial_{\th}l(y,D(z^{*},\th)), i.e. optimise the latent
space via solving an ordinary differential equation and minimise the loss
``naively'' with respect to the parameters ignoring the dependence of
$z(\t{)}$ to $\th$. \par\@@numbered@section{subsection}{toc}{Nesterov 2-nd
order accelerated gradient flow} \par The gradient flow described above is
based on naive gradient descent. This method may be slow in convergence. The
convergence per iteration can be further increased by considering Nesterov's
accelerated gradient descent. A second differential equation approximating
Nesterov's accelerated gradient method has been developed in
\cite[cite]{[\@@bibref{}{nesterov2016}{}{}]}. This $2^{nd}$ order ODE equation
for $z$ reads
\begin{equation}\frac{d^{2}z}{dt^{2}}+\frac{3}{t}\frac{dz}{dt}+\nabla_{z}l(y,D(z,\th))=0\end{equation}
for $dz/dt|_{0}=z(0)=0$. To be able to use this in the framework of the
gradient flow encoding we split the $2^{nd}$ order ODE into two interacting
$1^{st}$ order equations and solve the simultaneously. Specifically, solving
\begin{equation}\frac{dv}{dt}=\frac{3}{t+\epsilon}v+\nabla_{z}l(y,D(z,\th)),\
\ \ \ \ \ \ \frac{dz}{dt}=v,\end{equation} where $\epsilon$ ensures stability
at small $t$. \par\par\@@numbered@section{subsection}{toc}{Fixed grid ODE
solver} Being an ODE, the gradient flow can be solved with general ODE
solvers. Unfortunately, generic adaptive step size solvers are not useful
because the underlying ODE becomes stiff quickly during training, the step
size is reduced to extremely small values and the time it takes to solve
gradient flow ODEs at each iteration takes exorbitant amount of time. Fixed
time-step or time grid solvers can be used, despite the stiffness. However, we
empirically observed that these schemes can lead to instabilities in the
training, see Figure\leavevmode\nobreak\ \ref{fig:all_methods}. To demonstrate
this, we experimented with a $4^{th}$ order Runge-Kutta method with fixed step
size. The $\d{t}$ slices are predefined in logarithmic series such as $\d{t}$
is smaller closer to $t=0$, where integrands, $-\nabla_{z}l(y,D(z,\th))$, are
more rapidly changing. Similarly, $\alpha$ is empirically set to
$e^{(-2t/\t{)}}$ to facilitate faster convergence of $z$, see Figure
\ref{fig:conv}. For the GFE full adjoint method the integrands for each time
slice are saved during the forward pass so they can be use for the calculation
of the adjoin variable $\l$ in the backward pass. We used the same strategy
for both basic gradient flow and the $2^{nd}$ order model.
\par\begin{figure}[ht] \begin{center}
\includegraphics[width=216.81pt]{mnist_GFE_conv.pdf} \@@toccaption{{\lx@tag[
]{{1}}{ Typical convergence curve of $z$ during an ODE integration for GFE and
GFE $2^{nd}$ order using the fixed grid ODE solver. Experiments were performed
on MNIST. Mean cross-entropy loss plotted against time for three training
iterations, darker corresponds to earlier iteration. Curves are obtained using
the same batch of training images. The plot shows how convergence curves of
$z$ changes with training iterations. Darker curves are for earlier
iterations. In later iterations, the convergence curves reach lower loss,
meaning updates of the decoder network's parameters lead to more faithful
latent space representations, as expected.}}}\@@caption{{\lx@tag[: ]{{Figure
1}}{ Typical convergence curve of $z$ during an ODE integration for GFE and
GFE $2^{nd}$ order using the fixed grid ODE solver. Experiments were performed
on MNIST. Mean cross-entropy loss plotted against time for three training
iterations, darker corresponds to earlier iteration. Curves are obtained using
the same batch of training images. The plot shows how convergence curves of
$z$ changes with training iterations. Darker curves are for earlier
iterations. In later iterations, the convergence curves reach lower loss,
meaning updates of the decoder network's parameters lead to more faithful
latent space representations, as expected.}}}\end{center} \end{figure}
\par\par\@@numbered@section{subsection}{toc}{Adaptive minimise distance (AMD)
solver} The step size in any optimisation algorithm is of fundamental
importance for reaching a local extremum in any optimisation problem. The same
applies to the search for optimal latent space representation for each sample
by solving gradient flow equation through its time discretisation. As
discussed, fixed time step or time grid solvers empirically lead to
instabilities during training of the decoder. Adaptive step size ODE solvers
can theoretically solve the instabilities but due to the possible stiffness of
the gradient flow equation, their use is impeded. Note that, the adaptive step
size solvers focus on an accurate integration, which is not necessarily an
advantage for integrating gradient flow in the training of a decoder. A solver
that focuses on reaching a local extremum while searching for the optimal $z$
would be more advantageous. The path it takes to reach there may be less
crucial as long as it reaches the same local extremum. This is in contrast to
generic ODE solvers where the path can substantially change the end point,
thus an focusing on accurate integration is more important. To this end, we
develop an adaptive step size method, which guarantees a reduction of the loss
at each step. \par The method follows a similar structure to an explicit ODE
solver, such as the feed-forward Euler method, but without a fixed grid. The
problem lies in solving Equation \ref{eq:approx1} while taking time-steps of
appropriate size that reduce $l(y,D(z(t),\th))$ at each $t$. This approach is
in essence a gradient descent method that uses step-size selection
mechanisms\leavevmode\nobreak\ \cite[cite]{[\@@bibref{}{bertsekas_2016}{}{}]}.
Lastly, viewing the time-step in solving
$d{z}/dt=-\alpha(t)\nabla_{z}l(y,D(z(t),\theta))$ as a tool to minimise the
loss makes $\alpha(t)$ obsolete. Its role is now overtaken by the time-step
$\d{t}$ and it can be set to 1 for all t. \par In the AMD method, at each time
$t$ the time step is chosen based on finding the smallest $m=0,1,\dots,$ that
satisfies
\begin{equation}l\left(y,D(z(t_{n})-\beta^{m}s_{n}\nabla_{z}l(y,D(z(t_{n}),\theta)),\th)\right)<l\left(y,D(z(t_{n}),\th)\right)\end{equation}
with $\beta\in(0,1)$ (set as 0.75 in our experiments) and $s_{n}$ is a scaling
factor. At each time point $t_{n}$ the time step is chosen as
$\d{t}_{n}=\beta^{m}s_{n}$. The scaling factor is updated at each iteration as
\begin{equation*}\hat{s}_{n}=\max(\kappa s_{n-1},s_{0}),\
s_{n}=\min(\hat{s}_{n},s_{max}),\ s_{max}=10,\ s_{0}=1,\
\kappa=1.1\end{equation*} Based on this $t_{n+1}=t_{n}+\d{t}_{n}$ and
\begin{equation}z(t_{n+1})=z(t_{n})-\d{t}_{n}\nabla_{z}l(y,D(z(t_{n}),\theta)\end{equation}
\par\par\par At the end, if time step chosen goes beyond $\tau$, a smaller
time step is used to reach $\tau$ exactly. The solution of the integral
\ref{eq:approx1} is then $z(\tau)$. Furthermore, the AMD solver is using the
gradient of the convergence curve (see Figure \ref{fig:conv}) to assert if the
loss function is sufficiently optimised to assign a new final $\tau^{\prime}$
and stop in order to avoid unnecessary integration.
\par\par\begin{algorithm}[tb] \@@toccaption{{\lx@tag[ ]{{1}}{training GFE-amd
for one sample per batch}}}\@@caption{{\lx@tag[ ]{{\bf Algorithm~1}}{training
GFE-amd for one sample per batch}}}
\lx@algorithmicx@beginlist@{}\lx@algorithmicx@@item\noindent\hskip
0.0pt{Input: } z(0) \lx@algorithmicx@@item\noindent\hskip 0.0pt{Output:}
{$\mathbf{y}_{b}$, $\mathbf{\theta}_{D}$, $\mathbf{z}^{*}_{b}$} {{for}}\
{\emph{$b\in\\{1,2,...,B\\}$}}\ {{do}}
\lx@prepend@indentation@{}\algocf@printnl{1}\lx@algo@@endline\lx@algo@@startline
{{while}}\ {\emph{$\frac{d}{db}{L}_{CE}(b)<0.01$}}\ {{do}}
\lx@prepend@indentation@{\hskip 5.0pt\lx@algo@rule\hskip
10.00002pt}\algocf@printnl{2}\lx@algo@@endline\lx@algo@@startline {{while}}\
{\emph{$l(z({t}^{n}_{b}+\delta{t}^{n}_{b})<l(z({t}^{n}_{b}$))}}\ {{do}}
\lx@prepend@indentation@{\hskip 5.0pt\lx@algo@rule\hskip 10.00002pt\hskip
5.0pt\lx@algo@rule\hskip
10.00002pt}\algocf@printnl{3}\lx@algo@@endline\lx@algo@@startline
\lx@algorithmicx@@item\noindent\hskip 0.0pt$z({t}^{n+1}_{b})\leftarrow
z({t}^{n}_{b})+[-\nabla_{z}l(z({t}^{n}_{b}))]\delta{t}^{n}_{b}$
\lx@algorithmicx@@item\noindent\hskip
0.0pt{$\delta{t}^{n+1}_{b}\leftarrow\beta^{r}\delta{t}^{n}_{b}$}
\lx@algorithmicx@@item\noindent\hskip 0.0pt{$r\leftarrow r+1$}
\lx@prepend@indentation@{\hskip 5.0pt\lx@algo@rule\hskip 10.00002pt\hskip
5.0pt\lx@algo@rule\hskip 10.00002pt\hskip 5.0pt\lx@algo@rule\hskip
10.00002pt}\algocf@printnl{4}\lx@algo@@endline\lx@algo@@startline end
while\lx@prepend@indentation@{\hskip 5.0pt\lx@algo@rule\hskip 10.00002pt\hskip
5.0pt\lx@algo@rule\hskip
10.00002pt}\algocf@printnl{5}\lx@algo@@endline\lx@algo@@startline{}\lx@algorithmicx@@item\noindent\hskip
0.0pt{$n\leftarrow n+1$} \lx@prepend@indentation@{\hskip
5.0pt\lx@algo@rule\hskip 10.00002pt\hskip 5.0pt\lx@algo@rule\hskip
10.00002pt}\algocf@printnl{6}\lx@algo@@endline\lx@algo@@startline end
while\lx@prepend@indentation@{\hskip 5.0pt\lx@algo@rule\hskip
10.00002pt}\algocf@printnl{7}\lx@algo@@endline\lx@algo@@startline{}\lx@algorithmicx@@item\noindent\hskip
0.0pt{$z^{*}_{b}\leftarrow z({t}^{n}_{b})$}
\lx@algorithmicx@@item\noindent\hskip 0.0pt{$L(\theta_{D})\leftarrow
l(y,D(z^{*},\theta_{D})$} \lx@algorithmicx@@item\noindent\hskip
0.0pt{$\theta_{D}\leftarrow$ backpropagate $L(\theta_{D})$}
\lx@prepend@indentation@{\hskip 5.0pt\lx@algo@rule\hskip
10.00002pt}\algocf@printnl{8}\lx@algo@@endline\lx@algo@@startline end
for\lx@prepend@indentation@{}\algocf@printnl{9}\lx@algo@@endline\lx@algo@@startline{}\lx@algorithmicx@endlist
\end{algorithm} \par\par\par\@@numbered@section{section}{toc}{Experimental
setup} \par For training with MNIST and FashionMNIST datasets we implement a
sequential linear neural network. The decoder network architecture corresponds
to four gradually increasing linear layers with ELU non-linearities in-
between. The exact reverse is used for the encoder of the AE. A schematic
diagram of the GFE method is shown in Figure \ref{fig:schematic}.
\par\par\begin{figure}[ht] \begin{center}
\includegraphics[width=173.44534pt]{schematic_v1.pdf} \@@toccaption{{\lx@tag[
]{{2}}{ Schematic diagram of the GFE implementation. The two optimisation
correspond to (Step 1) the latent space optimisation i.e. 'encoding' and (Step
2) the parameter update of the neural network.}}}\@@caption{{\lx@tag[:
]{{Figure 2}}{ Schematic diagram of the GFE implementation. The two
optimisation correspond to (Step 1) the latent space optimisation i.e.
'encoding' and (Step 2) the parameter update of the neural
network.}}}\end{center} \end{figure} \par The network training is carried out
with a momentum gradient decent optimiser (RMSprop), learning rate $0.0005$,
$\varepsilon=1\times 10^{-6}$ and $\alpha=0.9$. The GFE and AE are considered
trained after one epoch and twelve epochs respectively.
\par\par\@@numbered@section{section}{toc}{Results and Discussion}
\par\par\par\par\par Initially a relative comparison between the full adjoint
and the approximate fixed grid GFE methods is carried out to assess the
relevance of the higher order term. Specifically, we carry out MNIST
experiments for a fixed network random seed, where we trained the Decoder
using the different GFE methods and computed cross entropy loss over the
validation set. The proper adjoint solution requires Equations
\ref{eq:adjoint1} and \ref{eq:adjoint2} to be solved for each slice of the
integral in Equation \ref{eq:adjoint3}. Given $N$ time-slices (for sufficient
accuracy $N\approx 100$), this requires $\mathcal{O}(5N)$ calls to the model
$D$ for each training image. The approximate method as in Equations
\ref{eq:approx1} and \ref{eq:approx2} requires only $\mathcal{O}(N)$ passes.
From Figure \ref{fig:all_methods}\leavevmode\nobreak\ (left) it is evident
that the 5-fold increase in computational time is not cost-effective as the
relative reduction in loss convergence per iteration is not significant. \par
Furthermore, to increase convergence with respect to training data the
accelerated gradient flow $2^{nd}$ order GFE is implemented in Section
\ref{sec:nesterov}. From Figure \ref{fig:all_methods}\leavevmode\nobreak\
(right), the accelerated gradient method increases initially the convergence
per iteration relative to GFE, nevertheless it is slightly more
computationally expensive due to solving a coupled system. Additionally, from
the same Figure certain stability issues are observed for both GFE and second
order GFE methods later on despite the initial efficient learning. In order to
guarantee stability the GFE-amd method is implemented as explained in Section
\ref{sec:amd}. The black curve in Figure
\ref{fig:all_methods}\leavevmode\nobreak\ (right) shows a clear improvement of
the GFE-amd over the later methods. Importantly, this result is robust to O of
the experiment. \par\par\par\begin{figure}[ht] \begin{center}
\includegraphics[width=433.62pt]{mnist_all_methods.pdf}
\@@toccaption{{\lx@tag[ ]{{3}}{ {Left} Validation mean cross-entropy loss
plotted against MNIST training iterations for the approximate and full adjoint
GFE methods. The full adjoint has a slight advantage over the approximate.
{Right} Validation mean cross-entropy loss plotted against MNIST training
iterations for the GFE, $2^{nd}$ order GFE and GFE-amd methods. The GFE-amd is
both more stable and approaches a better convergence relative to the other
methods}}}\@@caption{{\lx@tag[: ]{{Figure 3}}{ {Left} Validation mean cross-
entropy loss plotted against MNIST training iterations for the approximate and
full adjoint GFE methods. The full adjoint has a slight advantage over the
approximate. {Right} Validation mean cross-entropy loss plotted against MNIST
training iterations for the GFE, $2^{nd}$ order GFE and GFE-amd methods. The
GFE-amd is both more stable and approaches a better convergence relative to
the other methods}}}\end{center} \end{figure} \par\par\par\par\par A direct
comparison of the GFE-amd to a conventional AE for an MNIST training can be
seen in Figure \ref{fig:GFE_AE_step_relative}. The x-axis shows the number of
training images (instead of iterations) the algorithm sees until that point in
the training. The training is based on mini-batch training using the data with
replacement, going over the training data multiple times. The GFE-amd is
substantially superior in learning per training image, reaching near
convergence with at 800000 images, see Figure
\ref{fig:GFE_AE_step_relative}\leavevmode\nobreak\ (left). This is a
consequence of the efficiently optimised latent space. Nevertheless, this
comes at a higher computational cost for each iteration due to the ODE solver
as seen from Figure \ref{fig:GFE_AE_step_relative}\leavevmode\nobreak\
(right). Importantly, the optimisation of the network parameters is performed
using Adam optimiser for both AE and GFE models. So the difference we see can
be attributed to better gradients GFE model generates to update the decoder
network at each training iteration. \par\begin{figure}[ht] \begin{center}
\includegraphics[width=433.62pt]{mnist_AE_AMD_val.pdf} \end{center}
\@@toccaption{{\lx@tag[ ]{{4}}{ {Left} Validation mean cross-entropy loss
plotted against number of MNIST training images for the GFE-amd and AE
methods. The former shows a significant convergence with a tiny amound of
training images. {Right} Validation mean cross-entropy loss plotted against
time for the GFE-amd and AE methods. The latter is significantly faster to the
former as much more iterations are carried out in the same time-
span}}}\@@caption{{\lx@tag[: ]{{Figure 4}}{ {Left} Validation mean cross-
entropy loss plotted against number of MNIST training images for the GFE-amd
and AE methods. The former shows a significant convergence with a tiny amound
of training images. {Right} Validation mean cross-entropy loss plotted against
time for the GFE-amd and AE methods. The latter is significantly faster to the
former as much more iterations are carried out in the same time-
span}}}\end{figure} \par This increase in computation is not necessarily a
disadvantage considering the efficient learning of the GFE-amd method. In
Table \ref{tb:trainingimages}, the average cross entropy loss for a complete
test-set is recorded for both methods for some small number of training
images. The GFE-amd is able to learn quite well even after seeing a tiny
fraction of the total training data. Furthermore, the GFE-amd method
noticeably improves an AE trained decoder when it is used to test, the result
of an optimised latent space even without a network parameters update.
\par\par\par\par\par\begin{table}[ht] \centering\begin{tabular}[]{||c c c
c||}\hline\cr Number of Training images&AE&GFE-amd&train:AE test:GFE-
amd\\\\[1.50694pt] \hline\cr\hline\cr 480 (0.24$\%$)&0.2660&0.2098&0.2634\\\
\hline\cr 960 (0.49$\%$)&0.2618&0.1987&0.2525\\\ \hline\cr 1920
(0.98$\%$)&0.2488&0.1558&0.2323\\\ \hline\cr 3840
(1.95$\%$)&0.2195&0.1336&0.2038\\\ \hline\cr 5760
(2.9$\%$)&0.1954&0.1136&0.1829\\\ \hline\cr\end{tabular}
\@@toccaption{{\lx@tag[ ]{{1}}{ Test-set average cross entropy loss for
different number of training data, $\%$ reflects the percentage relative to
total training data needed for convergence of the method. Carried out for the
AE GFE-amd methods. The forth column represents testing the decoder of an AE
trained network with GFE-amd. The GFE-amd is far superior here in learning
with a limited image sample.}}}\@@caption{{\lx@tag[: ]{{Table 1}}{ Test-set
average cross entropy loss for different number of training data, $\%$
reflects the percentage relative to total training data needed for convergence
of the method. Carried out for the AE GFE-amd methods. The forth column
represents testing the decoder of an AE trained network with GFE-amd. The GFE-
amd is far superior here in learning with a limited image
sample.}}}\@add@centering\end{table} \par To verify the overall quality of the
method both the AE and GFE-amd are tested when converged as shown in Table
\ref{tb:loss}. The GFE-amd performs very similar to AE both for MNIST,
SegmentedMNIST and FMNIST. It is worth noting that the GFE-amd trainings are
on average converged at $1/12^{th}$ of the number of iterations relative to
the AE. For the segmented MNIST the networks are fully trained while seeing
only the first half (0-4) of the MNIST labels and they are tested with the
second half (5-9) of the labels. The GFE-amd shows a clear advantage over the
AE emphasizing the versatility of a GFE-amd trained neural network.
\par\par\begin{table}[ht] \centering\begin{tabular}[]{||c c c||}\hline\cr
Dataset&AE (Test-set)&GFE-amd (Test-set)\\\\[1.50694pt] \hline\cr\hline\cr
MNIST&0.0843&0.0830\\\ \hline\cr SegmentedMNIST&0.1205&0.1135\\\ \hline\cr
FMNIST&0.2752&0.2764\\\ \hline\cr\end{tabular} \@@toccaption{{\lx@tag[ ]{{2}}{
Test-set average cross entropy loss for trained networks with the AE and GFE-
amd methods. Segmented MNIST represents training for half the labels and
testing with the other half.}}}\@@caption{{\lx@tag[: ]{{Table 2}}{ Test-set
average cross entropy loss for trained networks with the AE and GFE-amd
methods. Segmented MNIST represents training for half the labels and testing
with the other half.}}}\@add@centering\end{table} \par Sample test-set
reconstructions with a fixed network random seed for GFE-amd and AE methods
are shown in Figure \ref{fig:images}. From Figure
\ref{fig:images}\leavevmode\nobreak\ (a) it is evident that the GFE-amd is
superior in producing accurate reconstructions with the very limited amount of
data. Figure \ref{fig:images}\leavevmode\nobreak\ (b) indicates that both GFE-
amd and AE generate similar reconstructions when properly trained.
\par\begin{figure}[ht] \begin{center}
\includegraphics[width=433.62pt]{training_images.pdf} \end{center}
\@@toccaption{{\lx@tag[ ]{{5}}{ {(a)} Test-set reconstructions for trained
GFE-amd (left) and AE (right) that only see $1\%$ of MNIST (top) and
FashionMNIST (bottom) training images. {(b)} Test-set reconstructions for
fully trained GFE-amd (left) and AE (right) with MNIST (top) and FashionMNIST
(bottom) training images. Note: The labels are identical in the respective
reconstructions.}}}\@@caption{{\lx@tag[: ]{{Figure 5}}{ {(a)} Test-set
reconstructions for trained GFE-amd (left) and AE (right) that only see $1\%$
of MNIST (top) and FashionMNIST (bottom) training images. {(b)} Test-set
reconstructions for fully trained GFE-amd (left) and AE (right) with MNIST
(top) and FashionMNIST (bottom) training images. Note: The labels are
identical in the respective reconstructions.}}}\end{figure} \par\par\par
Finally, to further compare latent space representation, we visualize the
samples in the latent space using the t-distributed stochastic neighbour
embedding (t-SNE) map, \cite[cite]{[\@@bibref{}{tsne2008}{}{}]}. This is
calculated for the GFE optimised $z\rightarrow z*$. This is shown for AE and
GFE, MNIST trained neural networks in Figure \ref{fig:tsne}. The latent space
representations are similar when both models see the entire MNIST dataset
multiple times. Similar t-SNE plots for models that only see 1\% of the data
during training are given in Figure \ref{fig:tsne_training}. Latent space
structure of GFE is very similar for both cases, while AE's latent space
structure is very different, not clustering the different numbers. This result
is inline with Table\leavevmode\nobreak\ \ref{tb:trainingimages}. GFE uses
training images more efficiently at each iteration thanks to latent space
optimisation to invert the decoder rather than using an approximation through
an encoder network. \par\par\begin{figure}[ht] \begin{center}
\includegraphics[width=346.89731pt]{tsne_GFE_AE.pdf} \end{center}
\par\@@toccaption{{\lx@tag[ ]{{6}}{ t-SNE map of the latent space plotted for
MNIST trained {Left} GFE and {Right} AE methods.}}}\@@caption{{\lx@tag[:
]{{Figure 6}}{ t-SNE map of the latent space plotted for MNIST trained {Left}
GFE and {Right} AE methods.}}}\par\begin{center}
\includegraphics[width=346.89731pt]{tsne_training.pdf} \end{center}
\par\@@toccaption{{\lx@tag[ ]{{7}}{ t-SNE map of the latent space plotted for
MNIST trained with $1\%$ training data {Left} GFE and {Right} AE methods. The
GFE latent space is well optimised even with a fraction of the data.
}}}\@@caption{{\lx@tag[: ]{{Figure 7}}{ t-SNE map of the latent space plotted
for MNIST trained with $1\%$ training data {Left} GFE and {Right} AE methods.
The GFE latent space is well optimised even with a fraction of the data.
}}}\end{figure} \par\par\par\@@numbered@section{section}{toc}{Conclusions} To
this end, a gradient flow encoding, decoder-only method was investigated. The
decoder depended gradient flow searches for the optimal latent space
representation, which eliminates the need of an approximate inversion. The
full adjoint solution and its approximation or leveraged for training and
compared. Furthermore, we present a $2^{nd}$ order ODE variant to the method,
which approximates Nesterov's accelerated gradient descent, with faster
convergence per iteration. Additionally, an adaptive solver that prioritises
minimising loss at each integration step is described and utilised for
comparative tests to the autoencoding model. The gradient flow encoding shows
a much higher data-efficiency than the autoencoding model.
\par\bibstyle{abbrv} \thebibliography\lx@bibitem{bertsekas_2016}
D.\leavevmode\nobreak\ P. Bertsekas. \lx@bibnewblock{\em Nonlinear
programming}. \lx@bibnewblock Athena Scientific, 2016.
\par\lx@bibitem{butcher2008numerical} J.\leavevmode\nobreak\ C. Butcher and
N.\leavevmode\nobreak\ Goodwin. \lx@bibnewblock{\em Numerical methods for
ordinary differential equations}, volume\leavevmode\nobreak\ 2.
\lx@bibnewblock Wiley Online Library, 2008. \par\lx@bibitem{chen2018neural}
R.\leavevmode\nobreak\ T. Chen, Y.\leavevmode\nobreak\ Rubanova,
J.\leavevmode\nobreak\ Bettencourt, and D.\leavevmode\nobreak\ K. Duvenaud.
\lx@bibnewblock Neural ordinary differential equations. \lx@bibnewblock In
{\em Advances in neural information processing systems}, pages 6571--6583,
2018. \par\lx@bibitem{DBLP:journals/corr/DinhKB14} L.\leavevmode\nobreak\
Dinh, D.\leavevmode\nobreak\ Krueger, and Y.\leavevmode\nobreak\ Bengio.
\lx@bibnewblock{NICE:} non-linear independent components estimation.
\lx@bibnewblock In Y.\leavevmode\nobreak\ Bengio and Y.\leavevmode\nobreak\
LeCun, editors, {\em 3rd International Conference on Learning Representations,
{ICLR} 2015, San Diego, CA, USA, May 7-9, 2015, Workshop Track Proceedings},
2015. \par\lx@bibitem{grathwohl2019ffjord} W.\leavevmode\nobreak\ Grathwohl,
R.\leavevmode\nobreak\ T.\leavevmode\nobreak\ Q. Chen, J.\leavevmode\nobreak\
Bettencourt, I.\leavevmode\nobreak\ Sutskever, and D.\leavevmode\nobreak\
Duvenaud. \lx@bibnewblock Ffjord: Free-form continuous dynamics for scalable
reversible generative models. \lx@bibnewblock{\em International Conference on
Learning Representations}, 2019. \par\lx@bibitem{deepsdf}
J.\leavevmode\nobreak\ Park, P.\leavevmode\nobreak\ Florence,
J.\leavevmode\nobreak\ Straub, R.\leavevmode\nobreak\ Newcombe, and
S.\leavevmode\nobreak\ Lovegrove. \lx@bibnewblock Deepsdf: Learning continuous
signed distance functions for shape representation. \lx@bibnewblock In {\em
Proceedings of the IEEE Conference on Computer Vision and Pattern
Recognition}, pages 165--174, 06 2019. \par\lx@bibitem{nesterov2016}
W.\leavevmode\nobreak\ Su, S.\leavevmode\nobreak\ Boyd, and
E.\leavevmode\nobreak\ J. Cand{{\\`{e}}}s. \lx@bibnewblock A differential
equation for modeling nesterov's accelerated gradient method: Theory and
insights. \lx@bibnewblock{\em Journal of Machine Learning Research}, 17(153):1
--43, 2016. \par\lx@bibitem{Tschannen2018RecentAI} M.\leavevmode\nobreak\
Tschannen, O.\leavevmode\nobreak\ Bachem, and M.\leavevmode\nobreak\ Lucic.
\lx@bibnewblock Recent advances in autoencoder-based representation learning.
\lx@bibnewblock{\em ArXiv}, abs/1812.05069, 2018. \par\lx@bibitem{tsne2008}
L.\leavevmode\nobreak\ van\leavevmode\nobreak\ der Maaten and
G.\leavevmode\nobreak\ Hinton. \lx@bibnewblock Visualizing data using {t-SNE}.
\lx@bibnewblock{\em Journal of Machine Learning Research}, 9:2579--2605, 2008.
\par\lx@bibitem{stackedAE} P.\leavevmode\nobreak\ Vincent,
H.\leavevmode\nobreak\ Larochelle, I.\leavevmode\nobreak\ Lajoie,
Y.\leavevmode\nobreak\ Bengio, and P.-A. Manzagol. \lx@bibnewblock Stacked
denoising autoencoders: Learning useful representations in a deep network with
a local denoising criterion. \lx@bibnewblock{\em J. Mach. Learn. Res.},
11:3371\textendash 3408, Dec. 2010. \par\endthebibliography
\par\par\@add@PDF@RDFa@triples\LTX@newpage\par\end{document}$


# Towards enabling reliable immersive teleoperation through Digital Twin: A
UAV command and control use case

Nassim Sehad1, Xinyi Tu2, Akash Rajasekaran1, Hamed Hellaoui1, Riku Jäntti1,
and Mérouane Debbah3 1Aalto University, School of Electrical Engineering,
Information and Communications Engineering, Espoo, Finland. 2Aalto University,
School of Engineering, Mechanical Engineering, Espoo, Finland. 3Khalifa
University of Science and Technology, Abu Dhabi, United Arab Emirates. Emails:
{nassim.sehad, xinyi.tu, akash.rajasekaran, hamed.hallaoui,
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

This paper addresses the challenging problem of enabling reliable immersive
teleoperation in scenarios where an Unmanned Aerial Vehicle (UAV) is remotely
controlled by an operator via a cellular network. Such scenarios can be quite
critical particularly when the UAV lacks advanced equipment (e.g., Lidar-based
auto stop) or when the network is subject to some performance constraints
(e.g., delay). To tackle these challenges, we propose a novel architecture
leveraging Digital Twin (DT) technology to create a virtual representation of
the physical environment. This virtual environment accurately mirrors the
physical world, accounting for 3D surroundings, weather constraints, and
network limitations. To enhance teleoperation, the UAV in the virtual
environment is equipped with advanced features that may be absent in the real
UAV. Furthermore, the proposed architecture introduces an intelligent logic
that utilizes information from both virtual and physical environments to
approve, deny, or correct actions initiated by the UAV operator. This
anticipatory approach helps to mitigate potential risks. Through a series of
field trials, we demonstrate the effectiveness of the proposed architecture in
significantly improving the reliability of UAV teleoperation.

###### Index Terms:

Teleoperation, Digital Twin, 5G and beyond, Unmanned Aerial Vehicles, Internet
of Things, Virtual Reality.

##  I Introduction

In recent times, the fusion of teleoperation and Internet of Things offers an
enhanced paradigm for controlling Unmanned Aerial Vehicles (UAVs). At the
nexus of this evolution are Extended Reality (XR) technologies, including
Virtual Reality (VR), Augmented Reality (AR), and Mixed Reality (MR), which
gained considerable attention thanks to their ability to provide a more
intuitive and immersive control experience. Another pivotal advancement that
complements this experience is haptic feedback technology that delivers
tactile information to operators [1]. [2] proposed VRHapticDrones, a system
providing haptic feedback in VR without wearable devices, although limited by
collision risks and trajectory planning complexities. Similarly, [3]
introduced HoverHaptics, autonomous quadcopters serving as haptic devices in
VR. This work addressed challenges like speed limitations, control accuracy,
and safety concerns using display techniques and collision avoidance
strategies.

Despite these advancements, the use of UAVs for teleoperation is still facing
challenges such as Network Latency (NL), environmental disturbances, and long-
distance operation, which can lead to an unstable control of the remote
system. This becomes more complex as UAVs could lack advanced equipment, such
as Lidar-based auto stop. Furthermore, novel user interfaces, such as Brain-
Computer Interfaces (BCIs) combined with XR for Human-Computer Interaction
(HCI), have been shown to introduce significant reaction delays exceeding 400
ms [4][5][6]. However, these studies do not explore methods for adapting to
network challenges or establishing autonomous control mechanisms in cases of
incidents caused by high delays. To address these challenges, we propose a
novel architecture in this paper that enables reliable immersive teleoperation
through the use of Digital Twins (DTs). The DT creates a virtual
representation of a physical UAV, rendered in the VR environment.

The main contributions of this paper are as follows:

  1. 1.

Proposing a novel architecture for DT-enabled UAV immersive teleoperation in a
VR environment, incorporating detailed 3D surroundings, weather, and network
constraints;

  2. 2.

Enhancing the virtual UAV with advanced features absent in the real UAV,
improving teleoperation reliability;

  3. 3.

Introducing an intelligent decision-making logic that leverages information
from both the virtual and physical environments to approve, deny, or modify
actions initiated by the UAV operator, thereby anticipating and mitigating
potential risks;

  4. 4.

Implementing the proposed architecture and validating its effectiveness
through a series of field trials, demonstrating its practical applicability.

The remainder of this paper is structured as follows: Section II presents an
overview of DT technology, discussing its various challenges when applied to
teleoperation. Section III introduces a motivating use case scenario that
illustrates the benefits of employing DTs to enhance reliability in immersive
teleoperation for UAV operators. In Section IV, we describe our proposed
architecture for DT-enabled UAV teleoperation. Section V provides a summary of
the field trials conducted and the results obtained. Finally, Section VI
concludes the paper and highlights potential avenues for future research.

##  II Digital twins for teleoperation

###  II-A Overview

DTs are digital replicas of physical assets, processes, and systems that can
be used for simulation, analysis, and optimization in various fields including
manufacturing, healthcare, and transportation [7].

In the context of teleoperation, DTs have shown great potential to improve the
efficiency and safety of remote systems [8]. DTs can be used to monitor the
health and performance of a remote system in real-time, enabling proactive
maintenance and minimizing downtime. This is particularly important for UAV
teleoperation, where any downtime can cause significant delays and
disruptions. DTs can also be used to identify potential issues before they
become critical, allowing for early intervention and prevention of more severe
problems. By simulating the behavior of a UAV in different environments, DTs
can identify potential hazards and risks before a UAV is deployed. This allows
better risk assessment and mitigation strategies, reducing the likelihood of
accidents or incidents.

###  II-B Challenges

Despite the promising potential of DTs in teleoperation, there are still
several challenges and research gaps that need to be addressed. One of the
main challenges is the availability and quality of data. DTs require a large
amount of data to accurately simulate the behavior of a remote system [9].
However, the data is often fragmented and inconsistent. This can result in
inaccuracies and errors of DTs, which further lead to suboptimal performance
or even failures in real-world operations. Therefore, it is crucial to ensure
that the data used for developing DTs is accurate, reliable, and up-to-date.
This requires access to a variety of data sources, which, for UAV
teleoperation, include sensors, telemetry data, and environment data.

Another challenge is the complexity and flexibility of the models used to
create DTs [10]. Models should be accurate to ensure DTs reliably reflect the
behaviors of the remote system. Models should also be flexible enough to
account for dynamic environments, where UAV teleoperation faces significant
challenges such as changing weather conditions, interference from other
wireless systems, and safety concerns for people and property. This requires
sophisticated modeling techniques to capture the complex interactions between
a UAV and its environment.

In addition, the integration of DTs with other systems and processes in the
teleoperation ecosystem can also be challenging, especially the synergies
between DTs and XR [11], which is critical for immersive teleoperation. In the
UAV case, data from DTs may need to be integrated with the UAV control system,
the teleoperation platform, or its VR interface, to ensure seamless operation
[12]. This requires a deep understanding of the technical requirements and
constraints of each system, as well as the ability to design and implement
integration solutions.

Finally, there are also technical challenges associated with the deployment of
DTs for teleoperation in real-world environments [13]. For instance, the
computing resources required to run a DT in real-time may be significant,
demanding high-performance computing infrastructure [14]. This can lead to
high power consumption, posing a challenge for the onboard computing of DTs in
head-mounted devices (HMDs) [15]. In addition, the communication and
networking infrastructure used to transmit data between DTs and the remote
systems should be robust and reliable [16], particularly in challenging
environments for UAV remote control. Addressing these technical challenges
requires a comprehensive understanding of the underlying technologies and
infrastructures required to support DT-enabled UAV teleoperation.

Overall, a standard architecture for integrating and developing DTs in UAV
immersive teleoperation is still missing, which should cover data acquisition,
modeling, simulation, and validation to ensure the accuracy and reliability of
DT-enabled UAV command and control.

Figure 1: Our proposed architecture for digital twin-enabled UAV
teleoperation.

##  III A use case scenario

Thomas, an operational engineer for a power line inspection firm, pilots UAVs
through challenging terrains like marshlands, forests, and mountains, using a
5G network. Using First Person View (FPV) mode and a HMD with joysticks, he
captures images and videos of power lines. Despite being a professional pilot,
ensuring the safety of the operation is challenging. This is primarily due to
long-distance operations, introducing significant delays, diminishing video
quality and making safe piloting arduous. Misjudgments not only endanger the
UAV but also the power grid.

To prevent the UAV from crashing, it is important to consider factors
exacerbating the delays, primarily unstable connections between the UAV and
Base Stations (BS) due to the mobility of the UAV, especially at high
altitudes. Here, the reduced gains of the down-tilted BS antennas can lead to
suboptimal handover and the selection of farther BSs [17]. Such issues reduce
bandwidth and may degrade or interrupt video transmission. Incorporating a DT
to pre-emptively execute commands and monitor the state of the UAV can help
mitigate these issues, allowing real-time adjustments to Thomas’s piloting
instructions.

##  IV Our proposed architecture for digital twin-enabled UAV teleoperaton

This section introduces a novel architecture for DT-enabled UAV teleoperation.
The proposed architecture enables a user wearing an HMD to reliably control a
remote UAV located anywhere in the world with 5G connectivity. Furthermore,
the proposed architecture also incorporates intelligent logic to monitor
different parameters from the virtual and the physical environments to
approve/deny/correct actions initiated by the operator in a way to anticipate
some risks. The components of the proposed architecture are illustrated in
Figure 1 and are explained in detail below.

###  IV-A Real UAV

The UAV is a quadcopter equipped with a 360∘ camera, an onboard computer, a
flight controller, and a 5G modem. The onboard computer receives user commands
via telemetry from the flight controller API module and sends captured video
to the video module for encoding and streaming. Both modules are located on
the edge server.

###  IV-B Virtual UAV/Digital twin

The DT of the UAV and the world is a virtual representation that runs on a
local computer situated near the edge server. Its primary responsibility is to
simulate real-world physics and the UAV’s behavior, thus providing a highly
accurate and synchronized virtual model of the physical UAV.

The DT consists of two main components: a Webots simulation environment [18]
and a UAV model from a Software In The Loop (SITL) simulator. The Webots
environment includes the world simulation and a virtual LiDAR sensor that
reports distance to facing objects, while the SITL simulator encompasses the
dynamics and motion of the UAV. The Webots simulation and the UAV model are
connected through socket communication. The world simulation in the DT
replicates real-world physics within the virtual environment, taking into
account various factors that a UAV may encounter in the field, such as weight,
dimensions, and weather conditions, namely, wind speed, direction, and
turbulence provided by the weather forecast API. A virtual LiDAR sensor is
integrated within the Webots simulation to accurately measure distances to
nearby obstacles, enabling the virtual UAV to detect and avoid potential
collisions in the simulated environment. Figure 2 illustrates the UAV in the
real world (without LiDAR) alongside its DT equivalent within a 3D map, which
is equipped with a LiDAR to prevent collisions with approaching obstacles. The
figure also represents the experimental scenario utilized for validating the
proposed architecture, as described in Section V Results and Analysis.

Figure 2: Comparison of the DT and the real-world UAV.

###  IV-C VR user

The VR user operates the UAV through a HMD with 5G connectivity, which
displays the live omnidirectional video feed from the UAV’s 360∘ camera in a
virtual environment as if they were physically present at the UAV’s location.
The HMD is equipped with joysticks that allow the user to intuitively control
the UAV’s movements, including altitude, yaw, pitch, and roll. To ensure low
latency and smooth control of the UAV, the HMD continuously sends consecutive
commands through the MQ Telemetry Transport (MQTT) protocol to the monitoring
and decision module on the edge server over a 5G network.

###  IV-D Edge server

The edge server is a crucial component of the DT-enabled UAV teleoperation
system, as it serves as an intermediary between various entities, such as the
real and virtual UAVs, the VR user, and external services. The edge server is
designed to be modular, with five main modules handling distinct
functionalities: video encoding/streaming, monitoring and decision, flight
control API, weather forecast API, and a measurement database. These modules
work in tandem to provide a seamless and reliable teleoperation experience.

The video encoding/streaming module is responsible for handling the 360∘ video
stream received from the real UAV’s camera. Upon receiving the video feed, it
encodes and compresses the stream for efficient transmission, before sending
it to the VR user’s HMD. This ensures that the user has a real-time view of
the UAV’s surroundings, enabling them to make informed decisions while
controlling the UAV.

The monitoring and decision module serves as the central hub for processing
and managing data from various sources. It receives the user command flow from
the VR user, UAV state and LiDAR distance data from the virtual UAV, wind
speed data from the weather forecast API module, and sensorial data from the
real UAV. Taking into account factors such as the current UAV state,
environmental conditions, and potential obstacles, the module makes command
decisions based on this information. It is important to note that the real UAV
and its digital twin have identical starting positions for any given mission.
This module incorporates a logic allowing to approve/deny commands initiated
by the operator. The decision to forward the user’s command to the flight
control API module is based on whether the condition in Equation (1) is met.
This equation is provided as

$\small LD\leq TH+EDL,$ (1)

where the LiDAR Distance ($LD$) refers to the distance to a facing object
reported by the DT UAV’s LiDAR, and the ThresHold ($TH$) is a fixed predefined
distance at which the UAV is programmed to stop. $EDL$ stands for the Error
Distance related to the Latency and is defined in Equation (2) as

$\small EDL=(CL+ENL)\times US,$ (2)

where $CL$ stands for the average of a measured Command Latency between the DT
and the real UAV which variates according to the Estimated Network Latency
($ENL$). As for $US$, it refers to the UAV’s Speed in the DT/Real-world at
position coordinates (${U_{x},U_{y},U_{z}}$).

This allows Delay Compensation (DC) and collision avoidance following
Algorithm 1, facilitating effective decision-making.

Algorithm 1 DC-based Monitoring and Decision Algorithm

$\vec{U}=(U_{x},U_{y},U_{z})$, $LD(t)$, $TH=1m$, $CL=146ms$

function Monitoring_Decision

$EDL\leftarrow US\times[CL+\operatorname{GET\\_ENL}(\vec{U})]$

if $LD(t)\leq TH+EDL$ and $command\\_direction=forward$ then

Send stop command to the flight control API

else

Forward command to the flight control API

end if

end function

function get_ENL($\vec{U}$)

$closest\\_NL\leftarrow\text{null}$

$min\\_distance\leftarrow\infty$

$NL\leftarrow\operatorname{read\\_NL\\_from\\_file}("measurement")$

for each $NL$ in $measurement$ do

$Dist\leftarrow\operatorname{compute\\_distance}(\vec{U},NL.lat,NL.lon,NL.alt)$

if $Dist<min\\_distance$ then

$ENL\leftarrow NL$

$min\\_distance\leftarrow distance$

end if

end for

return $ENL$

end function

The flight control API module is responsible for relaying command decisions
from the monitoring and decision module to both the real and virtual UAVs. By
ensuring that the command decisions are executed on both UAVs, this module
maintains synchronization between the real and DT UAVs. This allows for
accurate testing and validation of command decisions in the virtual
environment before applying them to the physical UAV, thus reducing the risk
of accidents and improving operational efficiency.

The weather forecast API module communicates with an external weather service
to obtain real-time wind speed data for the UAV’s location. It receives the
UAV’s position from the monitoring and decision module and sends the wind
speed information in return. By considering the impact of wind on UAV
performance, this module helps improve the accuracy and safety of UAV
operations.

The measurement module maintains a database containing information about the
average NL between the UAV and the edge server, based on the specific BS to
which the UAV is connected. This database is constructed from measurements, as
elaborated in Section IV-E Measurement Process. The database is indexed by the
current cell ID that the UAV is connected to and the UAV’s position, which are
obtained from the UAV’s cellular modem and Global Positioning System (GPS)
respectively. These indices are used to provide delay estimates for the
communication between the UAV and the edge server.

###  IV-E Measurement process

The network latency estimation process in our proposed architecture relies on
real-world measurements to account for network variability. In particular, we
concentrate on the non-reliable part of the network, which contributes to
fluctuations in delay times. This unreliable segment is related to the link
between the UAV and the edge server, where the elevated altitude and high
speed of the UAV have a significant impact on handover and throughput. This
impact is evident in the scenarios depicted in Figure 3 and the measurements
presented in Figures 4, 5, and 6. These figures illustrate the variations of
handover measurement compared to throughput, altitude, and latency,
respectively.

To obtain these measurements, we conducted a series of UAV flights,
representing two scenarios. In the first scenario, the UAV maintained a fixed
position while its altitude varied between $0-130m$ as depicted in Figure 3
(1) of the altitude trajectory. In the second scenario, the UAV maintained a
fixed altitude of $10m$ and performed a linear back-and-forth trajectory of
$100m$ at a maximum speed of $5m/s$ as shown in Figure 3 (2) of the horizontal
trajectory. During these UAV flights, UAV position, throughput, cell ID,
speed, and NL were recorded. This allowed us to construct the measurement
database and closely examine how changes in altitude and position within a
small area influenced the network’s overall performance. By incorporating
these real-world measurements into the DT, we can accurately estimate delays
and make more informed decisions when relaying user commands to the UAV.

Figure 3: UAV flight paths along serving base station IDs and locations with
(1) altitude trajectory and (2) horizontal trajectory.  Figure 4: Handover
compared to downlink throughput for trajectory 1. Figure 5: Handover rate in
function of altitude range for trajectory 1. Figure 6: Handover compared to
latency for trajectory 2.

##  V Results and Analysis

We present the results of two experimental scenarios depicted in Figure 3.
Figures 4 and 5 depict two critical performance metrics for a UAV
communication system during a handover between two base stations as altitude
changes. Figure 4 shows the measured throughput as a function of the BS’s cell
IDs during the handover, while Figure 5 displays the handover rate as a
function of altitude. Both figures relate to the first scenario shown in
Figure 3 (1), which illustrates the UAV’s altitude trajectory. Our
observations indicate a considerable reduction in throughput from an average
of $60Mbps$ to $10Mbps$ when the UAV ascends beyond $20m$ altitude due to high
interference from neighboring base stations in line of sight (LoS) with the
UAV, which leads to increased handover rates. Figure 6 presents latency versus
handover data from the second scenario shown in Figure 3 (2). We observed 6
base stations attachment with an increase in handover frequency near the lake
area, representing the edges of different cells, while the average delay
remained relatively constant with the distances to non-LoS BSs.

To validate our system and demonstrate the effectiveness of the DC method, we
first measured the average End-to-End (E2E) CL between the DT and the real
UAV. This measured CL showed an average delay of $146.04ms$ with a standard
deviation of $27.23ms$. We then conducted 20 UAV flight experiments with a VR
operator sending a command flow to both the DT and the real UAV that was
moving towards an obstacle with a speed ranging from $1-5m/s$, without
employing delay compensation. Theoretically, the difference between the
traveled distance of the real UAV and the digital one, denoted as
$D_{\Delta}$, would be proportional to the theoretical CL (TCL) described by
Equation (3) as

$\small D_{\Delta}=(CL+NL)\cdot US\implies TCL=D_{\Delta}\cdot\frac{1}{US}.$
(3)

The average TCL measured was $169.01ms$ with a standard deviation of $19.5ms$,
highlighting the effectiveness of our DT physics in accurately mimicking the
real world, as otherwise there would have been a significant difference
between the two distances due to the inertia and frictional force with the
air. The difference in the traveled distances matches the measured delay, with
an average error of $22.96ms$ ($14.57\%$) due to the NL variation,
demonstrating the reliability of our measurements and the realism of our
implementation, thanks to considering weather variations and other real-world
parameters. Upon applying our proposed DC method, the TCL ($TCL_{DC}$) showed
a significant reduction, with an average delay of $89.64ms$ and a standard
deviation of $11.86ms$, achieving a delay reduction of $53.04\%$ compared to
the previous measurements. Figure 7 shows the comparison of the measured CL,
deducted TCL, and deducted $TCL_{DC}$.

Figure 7: Measured latency and corresponding latency deducted from distances.

##  VI Conclusion

In conclusion, this paper presents a novel architecture that significantly
enhances the reliability of UAV teleoperation in critical scenarios. By
leveraging DT technology, we create an immersive virtual environment mirroring
the physical world, incorporating 3D surroundings, weather constraints, and
network limitations. Our architecture employs an intelligent logic that uses
data from both the virtual and physical worlds to make informed decisions
regarding operator-initiated actions. Furthermore, the proposed delay
compensation method has shown significant improvement in reducing network
latency, which makes this architecture even more promising for real-world
deployment. Field trials have substantiated the efficacy of this architecture,
showcasing its potential in mitigating risks and enhancing UAV teleoperation
reliability.

Looking ahead, the principles and mechanisms underscored herein can find
meaningful integration within healthcare, facilitating remote medical
interventions and surveillance. Furthermore, as the industry propels towards
Industry 5.0, our architecture can play an instrumental role in fostering
seamless human-machine collaboration, driving operational efficiencies and
ushering in a new era of industrial symbiosis. Future research can explore the
integration of advanced artificial intelligence algorithms for decision-
making, and autonomous flight capabilities.

## References

  * [1] Hassam Ahmed Malik, Shahzad Rasool, Adnan Maqsood and Rizwan Riaz  “Effect of haptic feedback on pilot/operator performance during flight simulation”  In _Applied Sciences_ 10.11 MDPI, 2020, pp. 3877
  * [2] Matthias Hoppe et al.  “VRHapticDrones: Providing haptics in virtual reality through quadcopters”  In _Proceedings of the 17th International Conference on Mobile and Ubiquitous Multimedia_ , 2018, pp. 7–18
  * [3] Parastoo Abtahi et al.  “Beyond the force: Using quadcopters to appropriate objects and the environment for haptics in virtual reality”  In _Proceedings of the 2019 CHI Conference on Human Factors in Computing Systems_ , 2019, pp. 1–13
  * [4] Soram Kim et al.  “P300 brain-computer interface-based drone control in virtual and augmented reality”  In _Sensors_ 21.17 MDPI, 2021, pp. 5765
  * [5] Tarik Taleb, Nassim Sehad, Zinelaabidine Nadir and Jaeseung Song  “VR-based immersive service management in B5G mobile systems: A UAV command and control use case”  In _IEEE Internet of Things Journal_ 10.6 IEEE, 2022, pp. 5349–5363
  * [6] Nassim Sehad et al.  “Locomotion-Based UAV Control Towards the Internet of Senses”  In _IEEE Transactions on Circuits and Systems II: Express Briefs_ 70.5 IEEE, 2023, pp. 1804–1808
  * [7] Mengnan Liu, Shuiliang Fang, Huiyue Dong and Cunzhi Xu  “Review of digital twin about concepts, technologies, and industrial applications”  In _Journal of Manufacturing Systems_ 58 Elsevier, 2021, pp. 346–361
  * [8] Roland Rosen, Georg Von Wichert, George Lo and Kurt D Bettenhausen  “About the importance of autonomy and digital twins for the future of manufacturing”  In _Ifac-Papersonline_ 48.3 Elsevier, 2015, pp. 567–572
  * [9] Xinyi Tu et al.  “A mixed reality interface for a digital twin based crane”  In _Applied Sciences_ 11.20 MDPI, 2021, pp. 9480
  * [10] Fei Tao et al.  “Digital twin modeling”  In _Journal of Manufacturing Systems_ 64 Elsevier, 2022, pp. 372–389
  * [11] Xinyi Tu et al.  “TwinXR: Method for using digital twin descriptions in industrial eXtended reality applications”  In _Frontiers in Virtual Reality_ 4 Frontiers, 2023, pp. 1019080
  * [12] Naoki Kikuchi, Tomohiro Fukuda and Nobuyoshi Yabuki  “Future landscape visualization using a city digital twin: Integration of augmented reality and drones with implementation of 3D model-based occlusion handling”  In _Journal of Computational Design and Engineering_ 9.2 Oxford University Press, 2022, pp. 837–856
  * [13] Wen Sun, Haibin Zhang, Rong Wang and Yan Zhang  “Reducing offloading latency for digital twin edge networks in 6G”  In _IEEE Transactions on Vehicular Technology_ 69.10 IEEE, 2020, pp. 12240–12251
  * [14] David Concannon, Ronan Flynn and Niall Murray  “A quality of experience evaluation system and research challenges for networked virtual reality-based teleoperation applications”  In _Proceedings of the 11th ACM workshop on immersive mixed and virtual environment systems_ , 2019, pp. 10–12
  * [15] Saeed Hamood Alsamhi et al.  “Computing in the sky: A survey on intelligent ubiquitous computing for uav-assisted 6g networks and industry 4.0/5.0”  In _Drones_ 6.7 MDPI, 2022, pp. 177
  * [16] Chao Yang et al.  “Extended Reality Application Framework for a Digital-Twin-Based Smart Crane”  In _Applied Sciences_ 12.12 MDPI, 2022, pp. 6030
  * [17] Samira Hayat et al.  “An experimental evaluation of LTE-A throughput for drones”  In _Proceedings of the 5th Workshop on Micro Aerial Vehicle Networks, Systems, and Applications_ , 2019, pp. 3–8
  * [18] Olivier Michel  “Cyberbotics ltd. webots™: professional mobile robot simulation”  In _International Journal of Advanced Robotic Systems_ 1.1 SAGE Publications Sage UK: London, England, 2004, pp. 5

  *[BSs]: base station