# Motion by mean curvature and Dyson Brownian Motion

Ching-Peng Huang  ,  Dominik Inauen  and  Govind Menon  Division of Applied
Mathematics, Brown University, 182 George St., Providence, RI 02912.  Institut
für mathematik, Universität Leipzig, D-04109, Leipzig, Germany.
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract.

We construct Dyson Brownian motion for $\beta\in(0,\infty]$ by adapting the
extrinsic construction of Brownian motion on Riemannian manifolds to the
geometry of group orbits within the space of Hermitian matrices. When $\beta$
is infinite, the eigenvalues evolve by Coulombic repulsion and the group
orbits evolve by motion by (minus one half times) mean curvature.

##  1\. Introduction

Fix $\beta\in(0,\infty]$ and $n$ standard independent Wiener processes
$\\{B_{j}(t)\\}_{j=1}^{n}$, $t\in[0,\infty)$. Dyson Brownian motion refers to
the unique weak solution to the Itô equation

$d\lambda_{j}=\sqrt{\frac{2}{\beta}}dB_{j}+\sum_{k\neq
j}\frac{dt}{\lambda_{j}-\lambda_{k}}\,,$ (1.1)

within the Weyl chamber

$\mathcal{W}_{n}=\\{(\lambda_{1},\ldots,\lambda_{n})\in\mathbb{R}^{n}\left|\lambda_{1}<\lambda_{2}\ldots<\lambda_{n}\right.\\}.$
(1.2)

In this paper, we introduce a new stochastic model for equation (1.1). To this
end, we first review the most relevant past constructions to provide some
context for our work.

Let $\mu_{n,\beta}$ denote the probability measure on $\mathcal{W}_{n}$ with
density proportional to the weight $\prod_{1\leq j<k\leq
n}(\lambda_{k}-\lambda_{j})^{\beta}$. There are two fundamentally different
classes of random matrices whose eigenvalues have law $\mu_{n,\beta}$. For
$\beta=1$,$2$ and $4$ these are the self-dual Gaussian ensembles (GOE, GUE and
GSE) of real symmetric, complex Hermitian and quaternionic matrices introduced
in the 1960s by Dyson, Gaudin and Mehta [15]. For $\beta\in(0,\infty]$, these
are the Gaussian $\beta$ ensembles (G$\beta$E) of real, symmetric tridiagonal
matrices, introduced by Dumitriu and Edelman [5].

Orthogonal polynomials played an important role in the first studies of these
ensembles (see [4, 5, 15]). However, dynamic models, especially equation
(1.1), play an important role in understanding universality [7]. Once a random
matrix ensemble has been chosen, natural time dynamics can be obtained by
replacing a single matrix drawn from the ensemble with a matrix-valued process
whose equilibrium measure is the given ensemble. Dyson obtained equation (1.1)
in this way for $\beta=1$, $2$ and $4$, replacing each Gaussian self-dual
ensemble with its associated Ornstein-Uhlenbeck process [6].

For $\beta\in(0,\infty]$ the extension of Dyson’s approach to the G$\beta$E
ensemble is a subtle problem. Holcomb and Paquette have constructed diffusions
of tridiagonal matrices whose eigenvalues satisfy equation (1.1) using
orthogonal polynomials and the Lanczos algorithm [9]. On the other hand,
Yabuoku has studied the eigenvalue process for the G$\beta$E diffusion
obtained by choosing independent Ornstein-Uhlenbeck processes on the diagonal
and independent Bessel processes on the off-diagonal (see [21, (2.1)]). He
shows that the eigenvalue process of a G$\beta$E diffusion depends on
additional minors and does not satisfy Dyson Brownian motion. The gap between
these results arises because the diffusions of tridiagonal matrices introduced
by Holcomb and Paquette are described somewhat implicitly in terms of a
separation between eigenvalue and eigenvector dynamics ([9, Section 3]); these
conditions do not hold for the G$\beta$E diffusion.

In the range $\beta\in[0,2]$ a model for equation (1.1) inspired by free
probability has been constructed by Allez, Bouchaud and Guionnet [1, 2]. They
construct a stochastic process $S_{t}$ of real, symmetric matrices whose
eigenvalues satisfy (1.1). Roughly, the process $S_{t}$ is a scaling limit
that interpolates between free convolution and standard convolution steps.
Despite the narrower range of $\beta$, and a different matrix model, this work
contains certain observations that reappear in [9].

Thus, the existence of natural time-dependent matrix models whose eigenvalues
satisfy (1.1) for arbitrary $\beta\in(0,\infty]$ is not fully settled. This is
the question we address.

The main contribution in this work is a geometric interpretation of equation
(1.1). For each $\beta\in(0,\infty]$, we construct a stochastic process
$M_{t}$ in the space of Hermitian matrices (equation (2.1) below) whose
eigenvalue process has the same law as the solutions to (1.1). The main new
tool in our approach is Riemannian geometry. Specifically, we use Riemannian
submersions of group orbits and a probabilistic interpretation of mean
curvature to obtain equation (1.1). The use of Riemannian submersion allows us
to view $\beta$ as a parameter that describes an anisotropic splitting between
noise in the tangent and normal directions (not an inverse temperature, as in
Dyson’s work). A similar role for $\beta$ has been observed by Holcomb and
Paquette [9, Thm.7]; our approach provides a systematic geometric explanation
for its importance. Second, we show that the Coulombic repulsion in equation
(1.1) corresponds to the mean curvature of group orbits. This is not a lucky
accident: it is a general principle corresponding to the gradient descent of
Boltzmann entropy for group orbits.

In order to explain the main new ideas in the simplest terms, we focus on the
explanation of the model, relying on previous work on well-posedness for Dyson
Brownian motion and standard calculations in random matrix theory to minimize
technicalities. The result in this paper is part of an effort by the authors
to develop previously unnoticed connections between three well-studied
problems: the construction of Brownian motion on Riemannian manifolds, Dyson
Brownian motion, and the isometric embedding problem for Riemannian manifolds.
At present, this interplay provides a new formulation of the embedding problem
for Riemannian manifolds [13], new interacting particle systems akin to Dyson
Brownian motion [16], and a systematic derivation of SDE for eigenvalue
processes of other classes of random matrices using Riemannian submersion
[11].

##  2\. Statement of results

###  2.1. The model

Let $\mathbb{H}(n)$ and $\mathbb{A}(n)$ denote the spaces of Hermitian and
anti-Hermitian matrices respectively. We equip these spaces with the Frobenius
norm $\|M\|^{2}=\mathrm{Tr}(MM^{*})$. Consider the smooth group action of the
unitary group $U(n)$ on the open, dense submanifold $V\subset\mathbb{H}(n)$ of
Hermitian matrices with simple spectrum given by

$\theta:U(n)\times V\to V,\,(Q,M)\mapsto\theta(Q,M)=\theta_{M}(Q)=QMQ^{*}\,.$

Since $U(n)$ is compact the isospectral orbits $\Sigma_{M}:=\theta(U(n),M)$
are embedded smooth submanifolds of $V$ (see Lemma 1 below). Consequently, at
any point $M\in V$ the tangent space $T_{M}V=\mathbb{H}(n)$ splits into the
$n^{2}-n$ dimensional tangent space $T_{M}\Sigma_{M}$ to the isospectral orbit
through $M$ and the $n$-dimensional normal space $N_{M}\Sigma_{M}$. Denote by
$P_{M}:\mathbb{H}(n)\to T_{M}\Sigma_{M}$ and $P^{\perp}_{M}:\mathbb{H}(n)\to
N_{M}\Sigma_{M}$ the respective orthogonal projections.

We construct a process $\\{M_{t}\\}_{t\geq 0}$ by a suitable projection of
standard Brownian motion on $\mathbb{H}(n)$ onto the tangent and normal spaces
to isospectral orbits. More precisely, assume given $M_{0}\in V$, let
$\\{X_{t}\\}_{t\geq 0}$ be a standard Brownian motion on $\mathbb{H}(n)$
starting at $X_{0}\in V$, and consider the Itô SDE

$dM_{t}=P_{M_{t}}dX_{t}+\sqrt{\frac{2}{\beta}}P^{\perp}_{M_{t}}dX_{t}\,.$
(2.1)

Using an explicit description of the projection operators and standard SDE
theory one can show (Lemma 3) that for every $\beta>0$ there exists a stopping
time $\tau_{\beta}$ and a solution $M_{t}$ of (2.1) on $[0,\tau_{\beta})$. We
then have

######  Theorem 1.

The eigenvalues of $M_{t}$ solve the equation (1.1) for
$t\in[0,\tau_{\beta})$.

The projection operators are smooth when the spectrum is simple. Since the
eigenvalues do not collide when $\beta\geq 1$, a simple bootstrap argument
shows that when $\beta\geq 1$, the stopping time $\tau_{\beta}=\infty$ (see
Lemma 3).

Theorem (1) may be established in a direct way once one has identified the Itô
equation (2.1). The main insight underlying equation (2.1), and thus Theorem
1, is a probabilistic interpretation of mean curvature. Let us now explain
this idea.

###  2.2. Brownian motion on Riemannian manifolds and mean curvature

There are two standard constructions of Brownian motion on Riemannian
manifolds using SDEs, referred to as the intrinsic and extrinsic constructions
respectively [10, 12]. The extrinsic construction goes as follows. Assume
$\mathcal{M}$ is a smooth $d$-dimensional manifold and assume given a smooth
embedding $u:\mathcal{M}\to\mathbb{R}^{q}$. Then Brownian motion on the
embedded submanifold $\Sigma=u(\mathcal{M})$ may be constructed as the
solution to the Stratonovich equation

$dZ_{t}=P_{Z_{t}}\circ dW_{t}$ (2.2)

where $P_{Z}$ is the orthonormal projection onto $T_{Z}\Sigma$ in
$\mathbb{R}^{q}$ and $W_{t}$ is a standard Wiener process in $\mathbb{R}^{q}$
[10].

The use of the Stratonovich formulation is crucial when one studies stochastic
processes on manifolds, since it accounts naturally for invariance under
coordinate transformations. But equation (2.2) also admits the equivalent Itô
formulation

$dZ_{t}=P_{Z_{t}}dW_{t}+\frac{1}{2}H(Z_{t})\,dt$ (2.3)

where $H(Z)$ is the mean curvature vector of the embedding $u$ at the point
$Z$. This identity is due to Stroock [20, Thm 4.4.2]; it was rediscoved by two
of the authors in their work on the isometric embedding problem [13, Thm.2].

The mean curvature vector of an embedding is defined111In differential
geometry , one frequently finds an additional normalization
$H=\frac{1}{n}\mathrm{tr}{II}$. We do not use this normalization factor,
adopting the convention of most texts in geometric analysis. as the trace of
the second fundamental form, but equations (2.2) and (2.3) show that it may be
approached directly from SDE theory. We obtain equation (2.3) by beginning
with (2.2), using the conversion rule between the Itô and Stratonovich
formulations to compute the Itô correction, recognizing finally that the Itô
correction has a fundamental geometric meaning. A related identity involving
mean curvature in the case of a Riemannian submersion was obtained by Pauwels
[18].

The intuitive content of equation (2.3) is that stochastic fluctuations in the
tangent space give rise to a ‘centrifugal force’ given by the mean curvature.
Let us illustrate this idea with an example. Let $\Sigma$ be a sphere of
radius $r$ in $\mathbb{R}^{q}$. We compute the projections explicitly, to see
that equations (2.2) and (2.3) take the form

$dZ_{t}=\left(I-\frac{Z_{t}Z_{t}^{T}}{|Z_{t}|^{2}}\right)\circ dW_{t},\quad
dZ_{t}=\left(I-\frac{Z_{t}Z_{t}^{T}}{|Z_{t}|^{2}}\right)\,dW_{t}-\frac{q-1}{2r}\frac{Z_{t}}{|Z_{t}|}\,dt.$
(2.4)

The Stratonovich form ensures that the constraint $|Z_{t}|=r$ holds for all
$t$. The ‘centrifugal force’ is $(q-1)/2r$ and it arises as follows. If we had
naively attempted to construct Brownian motion on the sphere with the Itô SDE

$dZ_{t}=\left(I-\frac{Z_{t}Z_{t}^{T}}{|Z_{t}|^{2}}\right)\,dW_{t},\quad\text{we
would find}\quad d|Z_{t}|=\frac{q-1}{2r|Z_{t}|}\,dt.$ (2.5)

Thus, the radial process $|Z_{t}|$ solves a deterministic equation, even
though the evolution of $Z_{t}$ is purely stochastic. Further, while each
point $Z_{t}$ moves tangentially to the sphere of radius $|Z_{t}|$, this
evolution has the effect of pushing spheres outwards normally by minus a half
times the mean curvature. (Observe that the mean curvature vector for the
sphere points inward).

###  2.3. Motion by mean curvature

This observation acquires greater depth when we recall the concept of motion
by mean curvature. A family of immersions
$u:\mathcal{M}\times[0,T]\to\mathbb{R}^{q}$ is said to evolve by motion by
mean curvature on the time interval $[0,T]$ if the velocity of each point
$u(x,t)$ is $H(u(x,t))$ for $t\in[0,T]$. Motion by mean curvature has been
extensively studied in geometric analysis [3, 8]. It is related to Dyson
Brownian motion as follows.

######  Theorem 2.

Assume $\beta=\infty$ in equation (2.1). The eigenvalues of the process
solving $dM_{t}=P_{M_{t}}dX_{t}$ evolve deterministically by Coulombic
repulsion for $t\in[0,\infty)$. Moreover, the corresponding isospectral orbits
$\Sigma_{M_{t}}$ move by minus a half times the mean curvature.

A minor difference with immersions in $\mathbb{R}^{q}$ is that the mean
curvature vector is the trace of the second fundamental form with respect to
the Frobenius metric on $\mathbb{H}(n)$. The flow is described precisely in
equation (3.6) below.

Theorem 2 corresponds to a gradient descent of Boltzmann entropy in the
following sense. As in the example above, we see that each matrix $M_{t}$ on
the group orbit moves tangentially (and stochastically), whereas the group
orbit as a whole evolves normally (and deterministically) by minus a half
times the mean curvature. The group orbits foliate the space $\mathbb{H}(n)$
and the group action is an isometry. In this setting, it is known that the
mean curvature at each point on the group orbit $\Sigma_{M_{t}}$ is the
gradient (with respect to the Frobenius norm) of
$-\log\mathrm{vol}(\Sigma_{M_{t}}))$ [17, p.3350]. The volume of the group
orbit depends only on the eigenvalues $\Lambda_{t}$ of $M_{t}$. By
interpreting $\Lambda_{t}$ as a macrostate, and each point $M_{t}$ as a
microstate, we see that $\log\mathrm{vol}(\Sigma_{\Lambda_{t}}))$ may be
interpreted as a Boltzmann entropy obtained using the theory of Brownian
motion on Riemannian manifolds. Further analysis from this viewpoint may be
found in [11, 13, 16].

##  3\. Proofs

We provide self-contained proofs of Theorem 1 and Theorem 2. An alternative
approach, which uses Pauwel’s theorem on the relationship between Riemannian
submersion and Brownian motion on Riemannian manifolds, and applies to other
eigenvalue processes, has been pursued by the first author [11].

######  Lemma 1.

Assume $M\in V$. The group orbit $\Sigma_{M}$ is an embedded submanifold in
$\mathbb{H}(n)$ of dimension $n^{2}-n$.

###### Proof.

Fix a diagonal matrix $\Lambda\in V$ with the same spectrum, so that
$\Sigma_{M}=\Sigma_{\Lambda}$. As can be checked, the isotropy group
$G_{\Lambda}=\\{Q\in U(n):Q\Lambda Q^{*}=\Lambda\\}$ is given by diagonal
matrices with entries $\mu_{i}\in\mathbb{S}^{1}$ and is hence isomorphic to
$\mathbb{T}^{n}$. Since $G_{\Lambda}$ is a closed subgroup of $U(n)$, the left
coset space $U(n)/G_{\Lambda}$ is a smooth manifold of dimension $n^{2}-n$.
The orbit map $\theta_{\Lambda}$ then descends to the quotient and gives a
smooth embedding $\theta_{\Lambda}:U(n)/G_{\Lambda}\to V$. Hence
$\dim\Sigma_{M}=\dim U(n)/G_{\Lambda}=n^{2}-n$. ∎

In the proof of Theorems 1 and 2 it is convenient to have an explicit
description of the tangent and normal spaces to an orbit.

######  Lemma 2.

Fix a diagonal matrix $\Lambda\in V$ and $M=Q\Lambda
Q^{*}\in\Sigma_{\Lambda}$. It then holds

$T_{M}\Sigma_{\Lambda}=\\{[QAQ^{*},M]:A\in\mathbb{A}(n)\text{ with
}A_{jj}=0\text{ for }j=1,\ldots,n\\}\,,$ (3.1)

and

$N_{M}\Sigma_{\Lambda}=\\{QDQ^{*}:D=\mathrm{diag}(\mu_{1},\ldots,\mu_{n})\text{
with }\mu_{j}\in\mathbb{R}\text{ for }j=1,\ldots,n\\}\,.$ (3.2)

###### Proof.

Fix any skew-hermitian $A\in\mathbb{A}(n)$ and consider the curve

$\gamma(t)=\exp(tA)\Lambda\exp(-tA)\,.$

Since $\gamma(t)\in\Sigma_{\Lambda}$ for any $t\in\mathbb{R}$ and
$\gamma(0)=\Lambda$, we find that

$[A,\Lambda]=\gamma^{\prime}(0)\in T_{\Lambda}\Sigma_{\Lambda}\,.$

Since $[iD,\Lambda]=0$ for any real diagonal matrix $D$ it follows by counting
dimensions that

$T_{\Lambda}\Sigma_{\Lambda}=\\{[A,\Lambda]:A\in\mathbb{A}(n)\text{ with
}A_{jj}=0\text{ for }j=1,\ldots,n\\}.\,.$

Observing that $[A,\Lambda]$ has empty diagonal for any $A\in\mathbb{A}(n)$ we
find that

$N_{\Lambda}\Sigma_{\Lambda}=\\{D=\mathrm{diag}(\mu_{1},\ldots,\mu_{n}):\mu_{j}\in\mathbb{R}\text{
for all }j=1,\ldots,n\\}\,.$

Since the action $\theta$ restricts to a transitive action on the orbits and
$\theta_{Q}:\Sigma_{\Lambda}\to\Sigma_{\Lambda}$ is a local diffeomorphism for
any $Q\in U(n)$ the claim follows. ∎

######  Lemma 3.

For every $\beta>0$ there exists a stopping time $\tau_{\beta}$ and a solution
$M_{t}$ of (2.1) on $[0,\tau_{\beta})$. Moreover, for $\beta\geq 1$ it holds
$\tau_{\beta}=+\infty$ almost surely.

###### Proof.

We claim that the maps $M\mapsto P_{M}$ and $M\mapsto P_{M}^{\perp}$ are
smooth on the open, dense subset $V\subset\mathbb{H}(n)$ of hermitian matrices
with simple spectrum. The local existence then follows from standard SDE
theory with $\tau_{\beta}$ being (bounded from below by) the first collision
time of the eigenvalues. The assertion that $\tau_{\beta}=+\infty$ almost
surely for $\beta\geq 1$ then follows from Theorem 1 and known non-collision
results of Dyson Brownian Motion for $\beta\geq 1$ (which are consequences of
”McKean’s argument”, see e.g. [14, Proposition 4.3]).

To show that the projections are smooth we fix $M=Q\Lambda Q^{*}\in V$ and
$N\in\mathbb{H}(n)$. From the description in Lemma 2 we find that

$N=P_{M}(N)+P^{\perp}_{M}(N)=[QAQ^{*},M]+QDQ^{*}$

for some $A\in\mathbb{A}(n)$ with empty diagonal and $D$ diagonal matrix with
real entries. Since $Q^{*}[QAQ^{*},M]Q=[A,\Lambda]$ has empty diagonal it
follows that

$D=\sum_{j=1}^{n}(Q^{*}NQ)_{jj}e_{j}\otimes e_{j}\,,$

where $e_{j}$ is the $j$-th standard basis vector of $\mathbb{R}^{n}$.
Therefore we have

$P_{M}^{\perp}(N)=\sum_{j=1}^{n}(Q^{*}NQ)_{jj}Qe_{j}\otimes
e_{j}Q^{*}=\sum_{j=1}^{n}\mathrm{Tr}(Q^{*}NQe_{j}\otimes
e_{j})P_{j}=\sum_{j=1}^{n}\mathrm{Tr}(NP_{j})P_{j}\,,$

where $P_{j}=P_{j}(M):=Qe_{j}(Qe_{j})^{*}$ is the (spectral) projection onto
the eigenspace associated to eigenvalue $\lambda_{j}(M)$. Observe that the
ordered spectrum $\lambda_{1},\ldots,\lambda_{n}:V\to\mathbb{R}$ is a family
of smooth functions on $V$ thanks to the implicit function theorem.
Consequently, also the projection operators $P_{j}$ are smooth on V, which can
be seen, for example, from the formula

$P_{j}(M)=-\frac{1}{2\pi
i}\oint_{\gamma}(M-\lambda\mathrm{Id})^{-1}\,d\lambda\,,$

where $\gamma\subset\mathbb{C}$ is a Jordan curve such that $\lambda_{j}(M)$
is the only eigenvalue of $M$ contained in its interior (see e.g. [19]). ∎

###### Proof (of Theorem 1):.

We get the equation for the eigenvalues by Itô’s formula. Consider the ordered
spectrum $\lambda_{1},\ldots,\lambda_{n}:V\to\mathbb{R}$. These are smooth
functions by the implicit function theorem and the Hadamard variation formulae
show that for $A,B\in\mathbb{H}(n)$

$D\lambda_{j}|_{Q\Lambda Q^{*}}(A)=\left(Q^{*}AQ\right)_{jj}\,,$ (3.3)

and

$D^{2}\lambda_{j}|_{Q\Lambda Q^{*}}(A,B)=2\sum_{k\neq
j}\frac{\left(Q^{*}AQ\right)_{kj}\left(Q^{*}BQ\right)_{kj}}{\lambda_{j}-\lambda_{k}}\,.$
(3.4)

Let then $\\{E_{\alpha}\\}_{\alpha=1}^{n^{2}}$ be the standard basis
(orthonormal with respect to the Frobenius metric) of $\mathbb{H}(n)$, i.e.,

$\\{E_{\alpha}\\}=\\{e_{j}\otimes e_{j},\frac{1}{\sqrt{2}}(e_{k}\otimes
e_{l}+e_{l}\otimes e_{k}),\frac{i}{\sqrt{2}}(e_{k}\otimes e_{l}-e_{l}\otimes
e_{k})\\}\,,$

where $e_{j}$ is the $j$-th standard basis vector of $\mathbb{R}^{n}$ and
$j=1,\ldots,n$, $1\leq k<l\leq n$. Then (2.1) reads

$dM_{t}=\sum_{\alpha=1}^{n^{2}}\left(P_{M_{t}}(E_{\alpha})+\sqrt{\frac{2}{\beta}}P^{\perp}_{M_{t}}(E_{\alpha})\right)dX_{t}^{\alpha}\,,$

where now $X_{t}^{\alpha}$ are jointly independent standard Wiener processes
on $\mathbb{R}$ starting at $X_{0}^{\alpha}$. By Itô’s formula it follows

$\displaystyle d\lambda_{j}$
$\displaystyle=\sum_{\alpha=1}^{n^{2}}D\lambda_{j}|_{M_{t}}\left(P_{M_{t}}(E_{\alpha})+\sqrt{\frac{2}{\beta}}P^{\perp}_{M_{t}}(E_{\alpha})\right)dX_{t}^{\alpha}$
$\displaystyle\quad+\frac{1}{2}\sum_{\alpha=1}^{n^{2}}D^{2}\lambda_{j}|_{M_{t}}\left(P_{M_{t}}(E_{\alpha})+\sqrt{\frac{2}{\beta}}P^{\perp}_{M_{t}}(E_{\alpha}),P_{M_{t}}(E_{\alpha})+\sqrt{\frac{2}{\beta}}P^{\perp}_{M_{t}}(E_{\alpha})\right)dt\,.$

Since $P_{M_{t}}(E_{\alpha})$ is tangent to the isospectral manifold we have
$D\lambda_{j}|_{M_{t}}\left(P_{M_{t}}(E_{\alpha})\right)=0$ for any $\alpha$.
Moreover, writing $M_{t}=Q_{t}\Lambda_{t}Q^{*}_{t}$ for a fixed $t>0$ we infer
from the description (3.2) that $Q^{*}_{t}P^{\perp}_{M_{t}}(E_{\alpha})Q_{t}$
is diagonal and hence

$D^{2}\lambda_{j}|_{M_{t}}\left(P_{M_{t}}(E_{\alpha}),P^{\perp}_{M_{t}}(E_{\alpha})\right)=D^{2}\lambda_{j}|_{M_{t}}\left(P^{\perp}_{M_{t}}(E_{\alpha}),P^{\perp}_{M_{t}}(E_{\alpha})\right)=0$

for any $\alpha$ thanks to (3.4). Consequently,

$d\lambda_{j}=\sqrt{\frac{2}{\beta}}\sum_{\alpha=1}^{n^{2}}D\lambda_{j}|_{M_{t}}\left(P^{\perp}_{M_{t}}(E_{\alpha})\right)dX_{t}^{\alpha}+\frac{1}{2}\sum_{\alpha=1}^{n^{2}}D^{2}\lambda_{j}|_{M_{t}}\left(P_{M_{t}}(E_{\alpha}),P_{M_{t}}(E_{\alpha})\right)dt\,.$

Observe now that

$Z_{t}:=\sum_{\alpha=1}^{n^{2}}\int_{0}^{t}D\lambda_{j}|_{M_{s}}\left(P^{\perp}_{M_{s}}(E_{\alpha})\right)dX_{s}^{\alpha}$

is a real-valued martingale with quadratic variation

$\displaystyle\langle Z\rangle_{t}$
$\displaystyle=\sum_{\alpha=1}^{n^{2}}\int_{0}^{t}\left(D\lambda_{j}|_{M_{s}}\left(P^{\perp}_{M_{s}}(E_{\alpha})\right)\right)^{2}ds\
=\sum_{\alpha=1}^{n^{2}}\int_{0}^{t}\left(D\lambda_{j}|_{M_{s}}\left(E_{\alpha}\right)\right)^{2}ds$
$\displaystyle=\sum_{\alpha=1}^{n^{2}}\int_{0}^{t}\left(\left(Q^{*}_{s}E_{\alpha}Q_{s}\right)_{jj}\right)^{2}ds=\sum_{\alpha=1}^{n^{2}}\int_{0}^{t}\langle
E_{\alpha},Q_{s}e_{j}\otimes e_{j}Q^{*}_{s}\rangle^{2}ds$
$\displaystyle=\int_{0}^{t}|Q_{s}e_{j}\otimes e_{j}Q^{*}_{s}|^{2}\,ds=t\,,$

since $\\{E_{\alpha}\\}$ is an orthonormal basis. By Lévy’s characterisation
it therefore follows that $dZ=dB_{j}$ for a real-valued Brownian motion
$B_{j}$, i.e. ,

$d\lambda_{j}=\sqrt{\frac{2}{\beta}}dB_{j}+\frac{1}{2}\sum_{\alpha=1}^{n^{2}}D^{2}\lambda_{j}|_{M_{t}}\left(P_{M_{t}}(E_{\alpha}),P_{M_{t}}(E_{\alpha})\right)dt\,.$

Now

$\displaystyle\frac{1}{2}\sum_{\alpha=1}^{n^{2}}D^{2}\lambda_{j}|_{M_{t}}\left(P_{M_{t}}(E_{\alpha}),P_{M_{t}}(E_{\alpha})\right)$
$\displaystyle=\frac{1}{2}\sum_{\alpha=1}^{n^{2}}D^{2}\lambda_{j}|_{M_{t}}\left(E_{\alpha},E_{\alpha}\right)$
$\displaystyle=\sum_{\alpha=1}^{n^{2}}\sum_{k\neq
j}\frac{\left(\left(Q^{*}_{t}E_{\alpha}Q_{t}\right)_{kj}\right)^{2}}{\lambda_{j}-\lambda_{k}}$
$\displaystyle=\sum_{\alpha=1}^{n^{2}}\sum_{k\neq j}\frac{\langle
E_{\alpha},Q_{t}e_{j}\otimes
e_{k}Q^{*}_{t}\rangle^{2}}{\lambda_{j}-\lambda_{k}}$
$\displaystyle=\sum_{k\neq j}\frac{|Q_{t}e_{j}\otimes
e_{k}Q^{*}_{t}|^{2}}{\lambda_{j}-\lambda_{k}}=\sum_{k\neq
j}\frac{1}{\lambda_{j}-\lambda_{k}}\,,$

which shows the claim. ∎

###### Proof (of Theorem 2):.

Suppose that $dM_{t}=P_{M_{t}}dX_{t}$. As in the proof of Theorem 1 it follows
that

$\displaystyle d\lambda_{j}$
$\displaystyle=\sum_{\alpha=1}^{n^{2}}D\lambda_{j}|_{M_{t}}\left(P_{M_{t}}(E_{\alpha})\right)dX_{t}^{\alpha}+\frac{1}{2}\sum_{\alpha=1}^{n^{2}}D^{2}\lambda_{j}|_{M_{t}}\left(P_{M_{t}}(E_{\alpha}),P_{M_{t}}(E_{\alpha})\right)dt$
$\displaystyle=\sum_{k\neq j}\frac{dt}{\lambda_{j}-\lambda_{k}}\,,$ (3.5)

which shows the deterministic Coulombic repulsion of the eigenvalues.

We moreover claim that the orbits $\Sigma_{M_{t}}$ evolve by time-reversed,
scaled mean curvature flow. More precisely, let $\Lambda_{t}$ be a diagonal
matrix with the same spectrum as $M_{t}$, so that
$\Sigma_{M_{t}}=\Sigma_{\Lambda_{t}}$. As in Lemma 2 we consider the closed
subgroup $T\subset U(n)$ of diagonal matrices with entries in
$\mathbb{S}^{1}$. We then define the family of embeddings

$F:U(n)/T\times[0,+\infty[\to\mathbb{H}(n),F(QT,t)=Q\Lambda_{t}Q^{*}\,,$

so that $\Sigma_{M_{t}}=F(U(n)/T,t)$. We claim that

$\frac{\partial}{\partial t}F(QT,t)=-\frac{1}{2}H_{Q\Lambda_{t}Q^{*}}\,,$
(3.6)

where $H_{M}$ is the mean curvature vector of $\Sigma_{M}$ at $M$, i.e.,

$H_{M}=\mathrm{trace}\,\mathrm{II}_{M}$

for $\mathrm{II}_{M}:T_{M}\Sigma_{M}\times T_{M}\Sigma_{M}\to N_{M}\Sigma_{M}$
the second fundamental form of the embedding $\Sigma_{M}\subset\mathbb{H}(n)$.

First, observe that $H_{Q\Lambda Q^{*}}=QH_{\Lambda}Q^{*}$, so that it
suffices to show (3.6) for $Q=Id$. Indeed, if $\\{X_{\alpha}\\}$ is an
orthonormal basis of $T_{\Lambda}\Sigma_{\Lambda}$, then
$Y_{\alpha}=QX_{\alpha}Q^{*}=d\theta_{Q}|_{\Lambda}(X_{\alpha})$ is an
orthonormal basis of $T_{Q\Lambda Q^{*}}\Sigma_{\Lambda}$, since
$d\theta_{Q}|_{\Lambda}$ is an isometry. Thus, setting $M=Q\Lambda Q^{*}$,

$H_{M}=\sum_{\alpha=1}^{n^{2}-n}\mathrm{II}_{M}(Y_{\alpha},Y_{\alpha})\,.$

Since $\mathbb{H}(n)$ is flat, $\mathrm{II}_{M}(Y_{\alpha},Y_{\alpha})$ is
given by

$\mathrm{II}_{M}(Y_{\alpha},Y_{\alpha})=P^{\perp}_{M}\left(\frac{d^{2}}{dt^{2}}|_{t=0}\gamma(t)\right)\,,$

for $\gamma:I\to\Sigma_{M}$ with $\gamma(0)=M$ and
$\gamma^{\prime}(0)=Y_{\alpha}$. We can take
$\gamma(t)=Q\tilde{\gamma}(t)Q^{*}$, where $\tilde{\gamma}:I\to\Sigma_{M}$
satisfies $\tilde{\gamma}(0)=\Lambda$ and
$\tilde{\gamma}^{\prime}(0)=X_{\alpha}$, so that

$\mathrm{II}_{M}(Y_{\alpha},Y_{\alpha})=P^{\perp}_{M}\left(Q\frac{d^{2}}{dt^{2}}|_{t=0}\tilde{\gamma}(t)Q^{*}\right)\,.$

However, from the description (3.2) it follows

$P^{\perp}_{Q\Lambda Q^{*}}(QHQ^{*})=QP^{\perp}_{\Lambda}(H)Q^{*}$

for any $H\in\mathbb{H}(n)$, from which we infer $H_{Q\Lambda
Q^{*}}=QH_{\Lambda}Q^{*}$.

Assume therefore $Q=Id$ and fix $X=[A,\Lambda]\in T_{\Lambda}\Sigma_{\Lambda}$
(recall the description from Lemma 2). In the following, for a given
$H\in\mathbb{H}(n)$ we denote by $D(H)$ the diagonal matrix with the same
diagonal as $H$. Defining $\gamma(t)=\exp(tA)\Lambda\exp(-tA)$ for
$t\in\mathbb{R}$ we find

$\mathrm{II}_{\Lambda}(X,X)=P^{\perp}_{\Lambda}\left(\frac{d^{2}}{dt^{2}}|_{t=0}\gamma(t)\right)=D\left([A,[A,\Lambda]]\right)=D\left([A,X]\right)\,.$
(3.7)

To compute $H_{\Lambda}$ we pick the explicit orthonormal basis of
$T_{\Lambda}\Sigma_{\Lambda}$ given by

$e_{kl}:=\frac{1}{\sqrt{2}}\left(e_{k}\otimes e_{l}+e_{l}\otimes
e_{k}\right)\,,iE_{kl}:=\frac{i}{\sqrt{2}}\left(e_{k}\otimes
e_{l}-e_{l}\otimes e_{k}\right)\,,$

for $1\leq k<l\leq n$. Observe that

$e_{kl}=\left[\frac{1}{\lambda_{l}-\lambda_{k}}E_{kl},\Lambda\right],iE_{kl}=\left[\frac{i}{\lambda_{l}-\lambda_{k}}e_{kl},\Lambda\right]\in
T_{\Lambda}\Sigma_{\Lambda}\,,$

so that the collection is indeed an orthonormal basis of
$T_{\Lambda}\Sigma_{\Lambda}$. It follows by (3.7) that

$\mathrm{II}_{\Lambda}(e_{kl},e_{kl})=\frac{1}{\lambda_{l}-\lambda_{k}}D([E_{kl},e_{kl}])=\mathrm{II}_{\Lambda}(iE_{kl},iE_{kl})\,.$

But $[E_{kl},e_{kl}]=e_{k}\otimes e_{k}-e_{l}\otimes e_{l}$, so

$\displaystyle H_{\Lambda}$
$\displaystyle=2\sum_{k<l}\mathrm{II}_{\Lambda}(e_{kl},e_{kl})=2\sum_{k<l}\frac{1}{\lambda_{l}-\lambda_{k}}\left(e_{k}\otimes
e_{k}-e_{l}\otimes e_{l}\right)$
$\displaystyle=2\sum_{k=1}^{n}\left(\sum_{l\neq
k}\frac{1}{\lambda_{l}-\lambda_{k}}\right)e_{k}\otimes e_{k}\,.$

Consequently, we find from (3.5)

$\frac{\partial}{\partial
t}F(T,t)=\sum_{j=1}^{n}\lambda_{j}^{\prime}e_{j}\otimes
e_{j}=\sum_{j=1}^{n}\left(\sum_{k\neq
j}\frac{1}{\lambda_{j}-\lambda_{k}}\right)e_{j}\otimes
e_{j}=-\frac{1}{2}H_{\Lambda_{t}}\,,$

as claimed. ∎

## References

  * [1] Romain Allez, Jean-Philippe Bouchaud, and Alice Guionnet, _Invariant beta ensembles and the Gauss-Wigner crossover_ , Physical Review Letters 109 (2012), no. 9, 094102.
  * [2] Romain Allez and Alice Guionnet, _A diffusive matrix model for invariant $\beta$-ensembles_, Electron. J. Probab. 18 (2013), no. 62, 30. MR 3078021
  * [3] Tobias Holck Colding, William P. Minicozzi, II, and Erik Kjæ r Pedersen, _Mean curvature flow_ , Bull. Amer. Math. Soc. (N.S.) 52 (2015), no. 2, 297–333. MR 3312634
  * [4] Percy. A. Deift, _Orthogonal polynomials and random matrices: a Riemann-Hilbert approach_ , Courant Lecture Notes in Mathematics, vol. 3, New York University, Courant Institute of Mathematical Sciences, New York; American Mathematical Society, Providence, RI, 1999. MR 1677884
  * [5] Ioana Dumitriu and Alan Edelman, _Matrix models for beta ensembles_ , J. Math. Phys. 43 (2002), no. 11, 5830–5847. MR 1936554
  * [6] Freeman J. Dyson, _A Brownian-motion model for the eigenvalues of a random matrix_ , J. Mathematical Phys. 3 (1962), 1191–1198. MR 148397
  * [7] László Erdős and Horng-Tzer Yau, _A dynamical approach to random matrix theory_ , Courant Lecture Notes in Mathematics, vol. 28, Courant Institute of Mathematical Sciences, New York; American Mathematical Society, Providence, RI, 2017. MR 3699468
  * [8] Yoshikazu Giga, _Surface evolution equations_ , Monographs in Mathematics, vol. 99, Birkhäuser Verlag, Basel, 2006, A level set approach. MR 2238463
  * [9] Diane Holcomb and Elliot Paquette, _Tridiagonal models for Dyson Brownian motion_ , arXiv preprint arXiv:1707.02700 (2017).
  * [10] Elton P. Hsu, _Stochastic analysis on manifolds_ , Graduate Studies in Mathematics, vol. 38, American Mathematical Society, Providence, RI, 2002. MR 1882015
  * [11] Ching-Peng Huang, _Eigenvalue processes in light of Riemannian submersion and gradient flow of isospectral orbits_ , arXiv preprint arXiv:2210.04858 (2022).
  * [12] Nobuyuki Ikeda and Shinzo Watanabe, _Stochastic differential equations and diffusion processes_ , second ed., North-Holland Mathematical Library, vol. 24, North-Holland Publishing Co., Amsterdam; Kodansha, Ltd., Tokyo, 1989\. MR 1011252
  * [13] Dominik Inauen and Govind Menon, _Stochastic Nash evolution_ , Preprint (2023).
  * [14] Eberhard Mayerhofer, Oliver Pfaffel, and Robert Stelzer, _On strong solutions for positive definite jump diffusions_ , Stochastic Process. Appl. 121 (2011), no. 9, 2072–2086. MR 2819242
  * [15] Madan Lal Mehta, _Random matrices_ , third ed., Pure and Applied Mathematics (Amsterdam), vol. 142, Elsevier/Academic Press, Amsterdam, 2004. MR 2129906
  * [16] Govind Menon and Tianmin Yu, _Stochastic flows in the Siegel half-space_ , Preprint (2023).
  * [17] Tommaso Pacini, _Mean curvature flow, orbits, moment maps_ , Trans. Amer. Math. Soc. 355 (2003), no. 8, 3343–3357. MR 1974691
  * [18] E. J. Pauwels, _Riemannian submersions of Brownian motions_ , Stochastics Stochastics Rep. 29 (1990), no. 4, 425–436. MR 1124160
  * [19] Michael Reed and Barry Simon, _Methods of modern mathematical physics. IV. Analysis of operators_ , Academic Press [Harcourt Brace Jovanovich, Publishers], New York-London, 1978. MR 0493421
  * [20] Daniel W. Stroock, _An introduction to the analysis of paths on a Riemannian manifold_ , Mathematical Surveys and Monographs, vol. 74, American Mathematical Society, Providence, RI, 2000. MR 1715265
  * [21] Satoshi Yabuoku, _Eigenvalue processes of symmetric tridiagonal matrix-valued processes associated with Gaussian beta ensemble_ , arXiv preprint arXiv:2204.00736 (2022).


# Physics-Informed LSTM-Based Delay Compensation Framework for Teleoperated
UGVs

Ahmad Abubakar, Yahya Zweiri, AbdelGafoor Haddad, Mubarak Yakubu,
Ruqayya Alhammadi, and Lakmal Seneviratne Manuscript submitted February, 2024.
This work was supported by Khalifa University under Award No. RC1-2018-KUCARS
(Corresponding author: Ahmad Abubakar.)A. Abubakar, A. Haddad, M. Yakubu, and
R. Alhammadi are with the Department of Mechanical Engineering, Khalifa
University of Science and Technology, Abu Dhabi, United Arab Emirates (e-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS>100041348@ku.ac.ae).Y. Zweiri and L.
Seneviratne are with the Advanced Research and Innovation Center (ARIC) and
Khalifa University Center for Autonomous Robotic Systems (KUCARS). They are
also with the Department of Mechanical Engineering and the Department of
Aerospace Engineering, respectively, Khalifa University of Science and
Technology, Abu Dhabi, United Arab Emirates (e-mail:
<EMAIL_ADDRESS>yahya.zweiri@ku.ac.ae).

###### Abstract

Bilateral teleoperation of low-speed Unmanned Ground Vehicles (UGVs) on soft
terrains is crucial for applications like lunar exploration, offering
effective control of terrain-induced longitudinal slippage. However, latency
arising from transmission delays over a network presents a challenge in
maintaining high-fidelity closed-loop integration, potentially hindering UGV
controls and leading to poor command-tracking performance. To address this
challenge, this paper proposes a novel predictor framework that employs a
Physics-informed Long Short-Term Memory (PiLSTM) network for designing
bilateral teleoperator controls that effectively compensate for large delays.
Contrasting with conventional model-free predictor frameworks, which are
limited by their linear nature in capturing nonlinear and temporal dynamic
behaviors, our approach integrates the LSTM structure with physical
constraints for enhanced performance and better generalization across varied
scenarios. Specifically, four distinct predictors were employed in the
framework: two compensate for forward delays, while the other two compensate
for backward delays. Due to their effectiveness in learning from temporal
data, the proposed PiLSTM framework demonstrates a 26.1% improvement in delay
compensation over the conventional model-free predictors for large delays in
open-loop case studies. Subsequently, experiments were conducted to validate
the efficacy of the framework in close-loop scenarios, particularly to
compensate for the real-network delays experienced by teleoperated UGVs
coupled with longitudinal slippage. The results confirm the proposed framework
is effective in restoring the fidelity of the closed-loop integration. This
improvement is showcased through improved performance and transparency, which
leads to excellent command-tracking performance.

###### Index Terms:

Bilateral teleoperation, delay compensation, model-free predictor, physics-
informed LSTM, UGVs.

##  I Introduction

The low-speed teleoperated Unmanned Ground Vehicles (UGVs) have gained
significant attention in various fields, including industrial automation,
military operations, and space exploration [1, 2, 3]. These UGVs are typically
used in scenarios that require careful navigation and interaction with the
remote environment [4]. A good example is the navigation of space rovers on
soft terrains in planetary exploration [5].

Terrain traversability is one of the key challenges for teleoperated UGVs. In
particular, soft terrains with some degree of roughness, such as regolith [6],
pose unique challenges for UGV navigation due to the longitudinal slippage
induced by wheel-terrain interaction. To enhance the navigation and ensure
precise UGV control, the operator should be aware of the induced slippage and
apply appropriate commands through a bilateral teleoperation system [7, 8].
Telepresence, which encompasses the operator’s situational awareness, enables
proactive and timely modification of control inputs, to optimize the terrain
traversability and the UGV’s stability and excellent command-tracking
performance. Thus, haptic feedback could be a good approach to render slippage
awareness in teleoperated UGV as it is difficult to detect slippage through
video feedback alone [9].

However, latency (network delays) within the communication channel between the
operator and the UGV can significantly affect telepresence by compromising the
fidelity of the closed-loop integration necessary for timely responses to
environmental changes [3]. This delay may lead to degraded performance and
transparency, with the potential to destabilize the system, particularly as
the communication channel becomes nonpassive, even with small time delays [10,
11]. Meanwhile, the threshold for small delays that cause instability varies,
depending on the specific system dynamics and application of the system [12].
To address the challenges posed by latency in low-speed teleoperated UGVs
subjected to induced longitudinal slippage, researchers have focused on
various approaches including passivity-based, energy-based, and control-based
strategies [3]. Fig. 1 illustrates how varying degrees of telepresence and
terrain traversability correlate with the overall performance of the
teleoperated UGV.

Figure 1: The conceptual relationship illustrates how the command-tracking
performance of teleoperated UGV is affected by both the terrain properties and
the latency.

In this paper, we examine control-based delay compensation techniques for low-
speed teleoperated UGVs subject to longitudinal slippage in the presence of
significant delays.

###  I-A Related Work

Luz et al. [13] highlight the importance of slippage awareness via haptic
feedback in teleoperated UGVs operating in low-light conditions, showing
improved control over induced slippage and faster task completion compared to
relying on virtual feedback alone. Nonetheless, such discrete haptic feedback
does not precisely convey precise slippage awareness to the operator [4]. In
[14], it is demonstrated that continuous haptic feedback, particularly based
on the induced velocity loss of the UGV, can render slippage awareness more
effectively than discrete feedback. Alternatively, the torque at the wheel
contact point has been adopted to provide continuous haptic feedback for
slippage awareness [15], though it is not as effective as the velocity loss-
based feedback. Consequently, these feedbacks significantly improve the
control of induced longitudinal slippage. Building on [14], the authors
further extended their approach in [7] by incorporating control of both
lateral slide and longitudinal slippage simultaneously. However, a trade-off
between command-tracking performance and system transparency was observed;
specifically, enhancing performance sometimes introduced non-passivity issues
due to negative slip conditions. To address these issues, a slippage
compensator based on a feedforward controller was proposed [7]. Moreover, Li
et al. [16] addressed the non-passivity issue by modifying the slide angle in
the UGV’s kinematics and employing the Time Domain Passivity Approach (TDPA).
Additionally, they proposed a slippage-dependent local controller in [9].
Their strategies resulted in more effective control of both longitudinal
slippage and lateral sliding, significantly enhancing the UGV’s command-
tracking performance on soft terrains. Notwithstanding, these studies have
overlooked the impact of latency (communication delays) on system performance.
Subsequent research addressed this by investigating delay effects and revealed
that even small delays could not only degrade the performance but introduce
additional non-passivity issues [10]. Additionally, they propose a Wave
Variable (WV) passivity approach with derivative control to stabilize the
system and enhance the command-tracking performance.

However, the proposed WV transformation with the derivative controller in [10]
leads to a more conservative nature and poor transparency. While various
studies have leveraged the passivity approaches of WV and TDPA [3, 4]. Others
have utilized control techniques that include model-based approaches, model-
mediated approaches, and model-free approaches [3, 17]. Among these, model-
free predictors have been widely employed in bilateral teleoperated UGV to
address delay issues, due to the complexities of modeling driver behavior and
the complex nonlinear dynamics of environment interactions [12].

Zheng et al. [11] introduced the first model-free predictor designed for
closed network systems operating below 1 Hz, demonstrating its capacity to
compensate for small delays (135 ms) and addressing the non-passivity issue in
the system. It was later validated on teleoperated UGVs to show its efficacy
in compensating large delays of 600 ms in real-world scenarios [18, 19]
Building on this, they further combined the model-free predictor with a model-
based approach to compensate for large delays of high-speed teleoperated UGVs
[12]. The results showed that this approach yielded more accurate predictions
than either approach alone with a small delay. Nevertheless, challenges remain
in dealing with the uncertainty of real network delays. Sridhar et al. [20]
modified the model-free predictor to better accommodate uncertain network
delays by introducing an adaptively varying parameter that depends on delay
fluctuation. Despite the enhanced robustness, the approach encountered issues
with oscillatory error gain in higher frequency coupling variables [3]. In a
recent study [21], the model-free predictor was modified by introducing an
additional degree of freedom as a potential solution to address both the
uncertainties in delays and the challenges of higher oscillatory error gain.
Results from a connected testbed with a delay of 175 ms demonstrate robustness
to varying delays and improved performance, leading to a 30% increase in the
fidelity of the closed-loop integration compared to the original predictor.
Nevertheless, despite these improvements, the technique falls short in
modeling the complex nonlinear, and temporal dynamic behaviors in teleoperated
UGV systems coupled with longitudinal slippage, particularly in handling large
delays, such as those encountered in long-distance Earth-to-Lunar networks.

###  I-B Contributions

In this work, we advance the recent delay compensation framework [21] using
the Recurrent Neural Network (RNN) to develop a high-fidelity closed-loop
integration for teleoperated UGVs subjected to longitudinal slippage. The
original contributions of this work are outlined as follows.

  1. 1.

We developed a novel physics-informed LSTM (PiLSTM)-based predictor framework
to improve performance and effectively compensate for large network delays.

  2. 2.

A comprehensive performance comparison of delay compensation is conducted
between the conventional predictor framework [21] and our PiLSTM-based
predictor framework in open-loop case studies.

  3. 3.

The developed framework is validated for its effectiveness in mitigating the
negative effects of significant real-network delays on the teleoperated UGVs
in closed-loop scenarios.

The structure of the remainder of this article is as follows. In Section II,
we delve into integrated system dynamics and controls for UGV bilateral
teleoperation on soft terrains. Section III discusses the detailed design of
the physics-informed LSTM-based predictor, including problem definition, data
generation and preprocessing, model development, and performance evaluation.
Section IV presents the human-in-the-loop experiment test platform and the
experimental design protocol. Section V details the comparative analysis of
delay compensation performance and evaluates overall teleoperation
effectiveness. Finally, Section VI summarizes the findings and their
significance, and outlines future research directions.

##  II System Dynamics and Controls

The core of this research lies in delay compensation for the bilateral
teleoperation of low-speed UGVs subjected to induced longitudinal slippage,
integrating non-holonomic motion commands and environmental force feedback to
improve operator awareness of terrain-induced slippage. The architecture of
the proposed predictor framework, as illustrated in Fig. 2, comprises several
key components that include the haptic device, the communication channel, UGV
operating on soft terrain, and both forward and backward predictors.

Figure 2: Proposed predictor framework, comprising two model-free predictors
namely; forward predictor and backward predictor, and a blended architecture
for the teleoperated UGV subjected to terrain-induced slippage. The Red lines
represent the Delayed case, and the Green Lines represent the Predicted case.

###  II-A Haptic device (master) model

The Phantom haptic device was utilized as the haptic interface. The two
Degrees of Freedom (DOF) of the haptic device (joint 1 and joint 2) were
employed to encode the linear and angular velocities of the UGV, as shown in
the operator side of Fig. 2. This approach is commonly used in tele-driving
UGV with non-holomorphic behavior [7]. The nonlinear dynamic equation of the
Phantom haptic device was derived in detail in [22], which was eventually
reduced to linear dynamics [16] in (1).

$\mathbf{M}_{m}\ddot{q}_{m}(t)+\mathbf{C}_{m}\dot{q}_{m}(t)=u_{m}(t)+f_{h}(t)$
(1)

where $\mathbf{M}_{m}$ represents the mass matrix of the robot,
$\mathbf{C}_{m}$ denotes the centrifugal matrix,
$q_{m}=\begin{bmatrix}q_{mv}&q_{m\omega}\end{bmatrix}^{T}$ is the joints
position variable, $u_{m}=\begin{bmatrix}u_{mv}&u_{m\omega}\end{bmatrix}^{T}$
is the joints control input, and
$f_{h}=\begin{bmatrix}f_{hv}&f_{h\omega}\end{bmatrix}^{T}$ represents the
external force applied by the human operator.

To address potential instability resulting from workspace mismatch between the
master and the slave, a new dynamic variable was introduced,
$x_{m}=\lambda\dot{q}_{m}+q_{m}\quad(0<\lambda<1)$. As a result, the robot
controller can be expressed as $u_{m}=\bar{u}_{m}+u_{m}^{*}$, where
$\bar{u}_{m}$ is the teleoperation controller, and $u_{m}^{*}$ is the local
controller, given by
$u_{m}^{*}=\mathbf{B}_{vi}\dot{q}_{m}+\mathbf{B}_{p}q_{m}$. The modified
linear dynamics result in a first-order system as in (2).

$\mathbf{\bar{M}}_{m}\dot{x}_{m}(t)+\mathbf{\bar{C}}_{m}x_{m}(t)=u_{m}(t)+f_{h}(t)$
(2)

where $\mathbf{\bar{M}}_{m}=\mathbf{M}_{m}/\lambda$ and
$\mathbf{\bar{C}}_{m}=\mathbf{C}_{m}/\lambda$ are the equivalent mass and
centrifugal force coefficients, respectively.
$x_{m}=\begin{bmatrix}x_{mv}&x_{m\omega}\end{bmatrix}^{T}$ is the new joints
variable for motion commands.

###  II-B UGV (slave) Model

For this research, the low-speed teleoperated UGV utilized was the Khalifa
University Lunar rover, which is an exact replica of the UAE Rashid rover [6],
as displayed in the remote environment of Fig. 2. The UGV’s mobility concept
relied on its kinematic models that dictate the non-holomorphic motion,
especially at low speeds. Inspired by the Boogie rover, the UGV locomotion
principle assumes that the speeds of the front and rear wheels are nearly
equal when operating on flat terrain. The kinematic model in (3) corresponds
to the assumption of pure rolling, is formulated in the UGV frame [7].

$\left[\begin{array}[]{c}v_{s}(t)\\\
\omega_{s}(t)\end{array}\right]=\left[\begin{array}[]{cc}1/2&1/2\\\
-1/(2b)&1/(2b)\end{array}\right]\left[\begin{array}[]{c}v_{r}(t)\\\
v_{l}(t)\end{array}\right]=\left[\begin{array}[]{c}u_{sv}(t)\\\
u_{s\omega}(t)\end{array}\right]$ (3)

where $v_{s}(t)$ and $\omega_{s}(t)$ represent the measured linear and angular
velocity of the UGV, respectively. $u_{sv}(t)$, and $u_{s\omega}(t)$ represent
the desired linear and angular velocity of the UGV, respectively. $b$ is the
distance between the right and the left wheels, $v_{r}(t)$ and $v_{l}(t)$ are
the right and the left measured wheel velocities, respectively. The
transformation matrix, denoted as $\mathbf{E}(b)$, is:

$\mathbf{E}(b)=\left[\begin{array}[]{cc}1/2&1/2\\\
-1/(2b)&1/(2b)\end{array}\right]$ (4)

On soft terrains, like loose soil, the relationship between the linear
velocity of each wheel and the product of the wheel’s angular velocity and
radius may not be equal due to wheel longitudinal slippage caused by
interaction forces from the terrain [9]. Here we use (5) to describe the level
of wheel slippage. This slippage deviation from the original kinematic model
(in (3)) was recognized by Li et al. [16]. They proposed a modified kinematic
model for the UGV on soft terrains, expressed by (6).

$s_{r}(t)=\frac{v_{rd}(t)-v_{r}(t)}{v_{r}(t)},\
s_{l}(t)=\frac{v_{ld}(t)-v_{l}(t)}{v_{l}(t)}$ (5)

$\left[\begin{array}[]{c}v_{s}(t)\\\
\omega_{s}(t)\end{array}\right]=\left[\begin{array}[]{c}u_{sv}(t)\\\
u_{s\omega}(t)\end{array}\right]-\mathbf{E}(b)\left[\begin{array}[]{c}v_{rd}(t)-v_{r}(t)\\\
v_{ld}(t)-v_{l}(t)\end{array}\right]$ (6)

Where $v_{rd}(t)$ and $v_{ld}(t)$ are the right and the left desired wheel
velocities, respectively. $s_{r}(t)$ and $s_{l}(t)$ are the right and the left
desired wheel slippage, respectively. Besides, we only care about the case of
$v_{rd}\geq v_{r}$ in (5), corresponding to $s_{r}\geq 0$, to avoid the non-
passivity issues due to negative slippage [7]. To solve this issue, their
feedforward controller (FFC) was adopted to compensate for the induced
slippage to improve the command-tracking performance in (6). Now the terrain-
induced longitudinal slippage $f_{e}(t)$, the functional relationship between
its linear and angular velocity, is defined based on the induced velocity loss
of the UGV, as in (7):

$\left[\begin{array}[]{c}f_{ev}(t)\\\
f_{e\omega}(t)\end{array}\right]=\left[\begin{array}[]{c}u_{sv}(t)-v_{s}(t)\\\
u_{s\omega}(t)-\omega_{s}(t)\end{array}\right]=\mathbf{E}(b)\left[\begin{array}[]{c}v_{rd}(t)-v_{r}(t)\\\
v_{ld}(t)-v_{l}(t)\end{array}\right]$ (7)

In this study, $f_{e}$ is the environmental force, which is assumed the
friction force exerted by the soft terrains on the wheel and the UGV has a
high-precision built-in velocity controller for each wheel motor, thus the
control input is $u_{s}$. As the induced slippage is well compensated by the
FFC, the command-tracking performance in (6) becomes ideal in (8). Refer to
[10, 16] for more detailed.

$\begin{bmatrix}v_{s}(t)\\\
\omega_{s}(t)\end{bmatrix}\overset{f_{e}\rightarrow
0}{\longrightarrow}\begin{bmatrix}u_{sv}(t)\\\ u_{s\omega}(t)\end{bmatrix}$
(8)

###  II-C Proposed Bilateral Teleoperator Controls

While the existing bilateral teleoperator controllers have served their
purpose in mapping coordination between the haptic device (master) and the UGV
(slave), as demonstrated by the delayed case shown with red lines in Fig. 2,
their performance and transparency are poor in the presence of large delays
$T$ [10]. To address this shortcoming, we propose new bilateral teleoperator
controllers designed to enhance performance and transparency while maintaining
overall system stability. These controllers ensure effective coordination
between the haptic device (master) and the UGV (slave), as demonstrated by the
predicted case with green lines in Fig. 2, even in the presence of large
delays $T$. The control designs incorporate simple proportional derivative
controllers, with the mathematical description of the transition outlined as
follows:

Master Controller: To achieve perfect coordination between $(x_{mv},v_{s})$
and $(x_{m\omega},\omega_{s})$ at the operator side. The joint control input
$u_{m}$ in (2) is switched from (9) to (10) by introducing predicted feedback
($\hat{f}_{ev}(t)$, $\hat{f}_{e\omega}(t)$):

$\left[\begin{array}[]{l}u_{mv}(t)\\\
u_{m\omega}(t)\end{array}\right]=\left[\begin{array}[]{l}-k_{mv}\
f_{ev}(t-T)\\\ -k_{m\omega}\ f_{e\omega}(t-T)\end{array}\right]$ (9)

$\left[\begin{array}[]{l}u_{mv}(t)\\\
u_{m\omega}(t)\end{array}\right]=\left[\begin{array}[]{l}-k_{mm}\
\hat{f}_{ev}(t)\\\ -k_{m\omega}\ \hat{f}_{e\omega}(t)\end{array}\right]$ (10)

Slave Controller: To achieve perfect coordination between $(x_{mv},v_{s})$ and
$(x_{m\omega},\omega_{s})$ at the remote side. The control input $u_{s}$ in
(6) is switched from (11) to (12) by introducing predicted motion commands
($\hat{x}_{mv}(t)$, $\hat{x}_{m\omega}(t)$):

$\left[\begin{array}[]{l}u_{sv}(t)\\\
u_{s\omega}(t)\end{array}\right]=\left[\begin{array}[]{l}k_{sv}\
x_{mv}(t-T)\\\ k_{s\omega}\ x_{m\omega}(t-T)\end{array}\right]$ (11)

$\left[\begin{array}[]{l}u_{sv}(t)\\\
u_{s\omega}(t)\end{array}\right]=\left[\begin{array}[]{l}k_{sv}\
\hat{x}_{m1}(t)\\\ k_{s\omega}\ \hat{x}_{m2}(t)\end{array}\right]$ (12)

The proportional control gains $k_{mv}$, $k_{m\omega}$, $k_{sv}$, and
$k_{s\omega}$ are maintained positive [14]. Therefore, in the designed
teleoperation system above, if the desired linear velocity $x_{mv}(t)$ is
greater than the actual linear velocity $v_{s}(t)$, a backward force will be
felt by the human operator that pushes back on the master robot joint 1.
Similarly, for any difference in angular velocities $x_{m\omega}(t)$ and
$\omega_{s}(t)$, a backward force will be felt on joint 2. Conversely, if the
desired linear velocity $x_{mv}(t)$ is smaller than the actual linear velocity
$v_{s}(t)$, a forward force will be felt by the human operator that pulls the
master robot joint 1 forward. Similarly, for any difference in angular
velocities $x_{m\omega}(t)$ and $\omega_{s}(t)$, a backward force will be felt
on joint 2. Thus, the force feedback serves as a slippage awareness guide to
the human operator, enabling more effective command execution, as demonstrated
in [10, 14, 16].

##  III Predictor framework Design

The fundamental form of the conventional model-free predictors adopted in this
paper was initially developed to compensate for delays in general networked
closed-loop systems [11], then it was later expanded for the high-speed
teleoperated UGV in [12]. In [21], the modified predictor has been developed
based on the first-order plus time delay system dynamics, as described in
(13).

$\displaystyle\dot{x}_{p}(t)$
$\displaystyle=\dot{x}(t-T)+\beta\left[x(t-T)-x_{p}(t-T)\right]$ (13)
$\displaystyle\quad+\alpha\left[\dot{x}(t-T)-\dot{x}_{p}(t-T)\right],$
$\displaystyle\hat{x}(t)$ $\displaystyle=x_{p}(t).$

The predictor equation, (13), incorporates the actual state $x(t-T)$ and its
derivative $\dot{x}(t-T)$ at a previous time step $t-T$; the predicted state
$x_{p}(t-T)$ and its derivative $\dot{x}_{p}(t-T)$ at the same previous time
step $t-T$; and the predicted state $x_{p}(t)$ at the current time step $t$,
to estimate the predicted undelayed state $\hat{x}(t)$, which is an
approximation of the actual undelayed state $x(t)$. The performance of the
predictor is adjusted using the regularization parameters $\beta$ and $\alpha$
for a specific delay $T$. Details of the predictor stability analysis and
performance can be found in [21].

The primary limitation of the above conventional predictor dynamics lies in
their failure to capture the complex, temporal, and unknown dynamics behaviors
of the coupling variable, significantly compromising the prediction
performance [21]. A promising solution to address this limitation is the
implementation of Recurrent Neural Network (RNN)-based predictors. These
predictors, as proven in [23, 24], are effective at learning nonlinear and
temporal dynamics behaviors, thus making them ideal tools for modeling highly
nonlinear and time-dependent systems. However, they often struggle with
generalization in complex dynamic systems [26, 28]. This challenge paves the
way for hybrid models that combine data-driven insights with some physical
constraints [29, 30]. Capitalizing on this, our study introduces a new RNN-
based predictor framework, fundamentally structured based on the Long Short-
Term Memory (LSTM) networks, augmented with some physical constraints, similar
to the architecture in [30]. The proposed framework consists of two distinct
predictors: a forward predictor tailored to predict motion commands and a
backward predictor tailored to force feedback, as in Fig. 2.

###  III-A Problem Definition

The proposed predictor is utilized to estimate the undelayed variables from
delayed variables. The focus of the models is to perform a single-step
prediction of the coupling variables, particularly in the context of a track-
following task. Given the nature of the task, it is imperative to accurately
predict future variables at a time step.

The problem can be defined as:

$\mathcal{H}_{\Theta}(X_{t},X_{t-1},X_{t-2},\ldots,X_{t-n})=Y_{t+1},$ (14)

where $\mathcal{H}_{\Theta}$ denotes the PiLSTM network with parameters
$\Theta$, mapping the sequence inputs $(X_{t},X_{t-1},X_{t-2},\ldots,X_{t-n})$
to a single output $Y_{t+1}$ at time step $t$ during the learning process.
Each input vector is defined as
$X=[x(t-T),x_{p}(t-T),\dot{x}(t-T),\dot{x}_{p}(t-T)]$ and corresponds to the
one-step ahead actual output $Y=x(t)$.

###  III-B Data and Preproccesing

The dataset was generated from the human-in-the-loop experiment of the
teleoperated UGV on soft terrains, which is discussed in detail in Section IV.
The dataset acquired different path patterns. We initiated the process by
gathering comprehensive data of motion commands $x_{m}(t)$ and the force
feedback signals $f_{e}(t)$, saved at a frequency of 10 Hz from the
experiments. The data comprise five essential variables for each of the
coupling variables, including the undelayed actual value $x(t)$, delayed
actual value $x(t-T)$, the delayed predicted value $x_{p}(t-T)$, the
derivative of the delayed actual value $\dot{x}(t-T)$, and the derivative of
the delayed predicted value $\dot{x}_{p}(t-T)$. These variables are assumed to
capture the essential information necessary for accurately modeling the
complexity of the empirical data.

The training dataset was recorded from three randomly selected trained
operators with a sampling rate of 10 Hz for a total simulation time of 300 s
(comprising about 3000 acquisitions). For the testing datasets, three sets
were collected from different trained operators, with each testing dataset
maintaining the same sampling rate as the training set but spanning a shorter
simulation duration of 30 s. Furthermore, the training dataset is further
partitioned into two distinct segments for training and validation purposes,
in the ratio of 70% to 30%, respectively. Each segment is structured with the
input features $x(t-T)$, $x_{p}(t-T)$, $\dot{x}(t-T)$, and $\dot{x}_{p}(t-T)$;
and the target output feature $x(t)$, which represents the ground truth values
that the model aims to learn and predict. Additionally, the input features
were preprocessed by normalizing and scaling the features to ensure
homogeneity among the features, similar to the approach described in [24].

###  III-C Model Development

Here, we discuss the detailed design of the four distinct models, each
tailored to one of the four coupling variables of the bilateral teleoperation
system. These include the desired linear velocity $x_{mv}(t)$, desired angular
velocities $x_{m\omega}(t)$, longitudinal force feedback $f_{ev}(t)$, and
lateral force feedback $f_{e\omega}$.

####  III-C1 Network Architecture

The architecture of our physics-informed LSTM network is designed to handle
sequential data with complex nonlinearity and temporal dependencies [30]. The
network is composed of an input layer, a dense layer, an activation layer,
LSTM layers, and an output layer, along with a physics constraint, as depicted
in Fig. 3. Each LSTM layer captures different levels of abstraction within the
sequence data.

The input layer receives the input sequence, composed of a feature vector of
fixed length, denoted as $(X_{1},X_{2},\ldots,X_{n})$, where each $X_{i}$ is
defined as
$X_{i}=[x(t_{i}-T),x_{p}(t_{i}-T),\dot{x}(t_{i}-T),\dot{x}_{p}(t_{i}-T)]$ for
$i=1,\ldots,n$. An LSTM layer is comprised of a series of LSTM cells, each
employing a sigmoid activation function $\sigma$ for gate operations and a
hyperbolic tangent $\tanh$ for updating cell states. Its unique feedback
mechanism, where the output of each cell is fed back into itself at the next
time step, endows the network with the capability to capture long-term
temporal dependencies from the input sequence. The output layer produces
single-step predicted values, represented by $\hat{Y}_{n}$, where
$\hat{Y}=\hat{x}(t_{i})=x_{p}$. Thus, each of our four models shares the same
architectural structure but may vary in the number of layers and units
according to their complexity.

Figure 3: Physics-informed LSTM architecture for the proposed predictor:
Number of fixed-length input features $n$, number of dense layer hidden units
$m$, number of LSTM layers $k$, number of LSTM layer hidden units $l$, and a
single-step target output at $n$.

####  III-C2 Loss function formulation

Having discussed how to design the network architecture, we next address how
to obtain the network parameters $\Theta$. For this, the optimization problem
is defined as follows:

$\Theta^{*}=\min\mathcal{L}(\Theta)\vspace{-5pt}$ (15)

where $\mathcal{L}$ represents a custom loss function, which must be
differentiable with respect to its parameters. In neural networks addressing
this kind of regression problem, the mean absolute error is the most commonly
used for data loss as given by:

$\mathcal{L}_{d}(\Theta)=\frac{1}{N}\sum_{i=1}^{N}\left|x_{p}(t_{i})-x(t_{i})\right|\vspace{-5pt}$
(16)

where $N$ represents the total number of samples in a run. This loss function
is selected due to its inherent benefit of reducing the impact of outliers. To
further enhance model generalization to unseen data [30], some physical
constraints are enforced on the derivative of the variable. By employing a
method in [25] to approximate the temporal derivative of past states on the
right-hand side of the first-order dynamics in (13), it is transformed into a
standard Delay Differential Equation (DDE). The compacted form of this DDE is
described in (17), with $\psi(t)=t-T$. In turn, the loss function assimilates
the DDE constraint, which introduces a soft regularization in the optimization
process to ensure compliance with fundamental principles. The physics-informed
loss is defined based on the L2-norm of the DDE residual, as presented in
(18).

$\displaystyle\dot{x}_{p}(t)$
$\displaystyle=\mathcal{F}(t,x(\psi(t)),x_{p}(\psi(t))$ (17)
$\displaystyle\hat{x}(t)$ $\displaystyle=x_{p}(t)$
$\mathcal{L}_{p}(\Theta)=\frac{1}{M}\sum_{j=1}^{M}\left\|\dot{x}_{p}(t_{j})-\mathcal{F}(t_{j},x(\psi(t_{j})),x_{p}(\psi(t_{j}))\right\|^{2}\vspace{-5pt}$
(18)

where $M$ represents the number of collocation points where the DDE constraint
is enforced. The overall cost function $\mathcal{L}(\Theta)$ is a weighted sum
of these two loss terms and is minimized with respect to the network
parameters $\Theta$ during the training process. By optimizing this loss, the
network learns a model $\mathcal{H}_{\Theta}$ that not only fits the observed
data but also aligns with the specified dynamics of the system.

$\mathcal{L}(\Theta)=\mathcal{L}_{d}(\Theta)+\lambda\mathcal{L}_{p}(\Theta)$
(19)

The coefficient $\lambda$ is a hyperparameter that balances the two loss
components [30], allowing you to adjust the relative importance of data
fitting and physical constraint adherence during training as depicted in Fig.
3.

####  III-C3 Training

The training protocol for each of the models is initiated by inputting a
sequence of 50 previous data points for each model (5 s of historical data).
This predetermined sequence length, denoted as $n$, serves as the input vector
for the network. The network parameters undergo iterative refinement through
an optimization routine predicated on a composite loss function, as delineated
in (19). The backpropagation through time (BPTT) method is applied, and
gradient descent algorithms are also employed using the Adam optimizer for its
proficiency in stable convergence and adaptive learning rate adjustments. The
pseudo-code of the training process, encapsulated in Algorithm 1, delineates
the computational steps undertaken during the model learning phase.

Algorithm 1 Physics-Informed LSTM Training Pseudo Code

INPUT: $X=[X_{1},X_{2},\ldots,X_{n}]$, each $X_{i}\in\mathbb{R}^{n}$

OUTPUT: $Y=Y_{n}$, where $Y_{n}\in\mathbb{R}^{n}$

PARAMETERS: $\Theta$ = [$W_{in}$, $b_{in}$, $W_{f}$, $U_{f}$, $b_{f}$,
$W_{\tilde{c}}$, $U_{\tilde{c}}$, $b_{\tilde{c}}$, $W_{i}$, $U_{i}$, $b_{i}$,
$W_{o}$, $U_{o}$, $b_{o}$, $W_{out}$, $b_{out}$]

INITIALIZE: hidden state $h_{0}=\vec{0}$, cell state $c_{0}=\vec{0}$, learning
rate $\alpha$, physics loss weight $\lambda$

DEFINE: differential operator $\mathcal{D}$, DDE constraint $\mathcal{R}$

TRAIN and
VALIDATE$(X_{train},Y_{train},X_{val},Y_{val},\mathcal{D},\mathcal{R})$

for each epoch do

for each batch $(X_{batch},Y_{batch})$ in $(X_{train},Y_{train})$ do

Compute LSTM outputs for batch: $H_{batch}$

Compute Dense output for batch: $\hat{Y}_{batch}$

Apply $\mathcal{D}$ to get $\dot{\hat{Y}}_{batch}$

Compute data loss: $\mathcal{L}_{d}(\Theta)=f(\hat{Y}_{batch},Y_{batch})$

Compute phy. loss:
$\mathcal{L}_{p}(\Theta)=f(\dot{\hat{Y}}_{batch},\mathcal{R}(Y_{batch}))$

Compute total loss:
$\mathcal{L}(\Theta)=\mathcal{L}_{d}(\Theta)+\lambda\mathcal{L}_{p}(\Theta)$

Update $\Theta$ using gradient descent:

$\Theta=\Theta-\alpha\nabla_{\Theta}\mathcal{L}(\Theta)$

end for

Take $(X_{val},Y_{val})$

Compute output for validation: $\hat{Y}_{val}$

Compute validation loss: $\left\lVert\hat{Y}_{val}-Y_{val}\right\rVert_{2}$

end for

CHECK POINT: $\Theta^{*}=\Theta_{best}$

####  III-C4 Hyperparameter tuning

Bayesian optimization is employed to obtain the optimal hyperparameters for
the network models. Table I summarizes the optimal hyperparameters as obtained
for each of the four distinct models.

TABLE I: Optimal Hyperparameters for Predictors Coupling Variables |  Initial Learn Rate |  Dense Hidden Units, $m$ |  LSTM Depth, $k$ |  LSTM Hidden Units, $l$ |  Threshold
---|---|---|---|---|---
$x_{mv}$ |  0.00230 |  118 |  2 |  162 |  0.07
$x_{m\omega}$ |  0.00016 |  174 |  3 |  92 |  1.12
$f_{ev}$ |  0.00095 |  151 |  2 |  188 |  0.80
$f_{e\omega}$ |  0.00076 |  207 |  4 |  75 |  0.92

###  III-D Prediction Performance

During the training process, our methodology entailed a comprehensive
evaluation of various RNN layers that include LSTM, GRU, and hybrid LSTM-GRU
to have the best models. For rigorous model validation, we utilized testing
datasets that were distinct from those employed during the training phase, to
prove the model generalization and their efficacy in extrapolating insights
from new, unseen data. Thus, to prevent overfitting the model with the
training data, the performance of the developed models is quantified using the
Root Mean Square Error (RMSE), mathematically described as follows:

$RMSE=\sqrt{\frac{1}{N}\sum_{i=1}^{N}(x_{p}(t_{i})-x(t_{i}))^{2}}$ (20)

where $N$ represents the total number of samples within the dataset. The
evaluation was conducted in an offline mode, enabling a thorough analysis of
the models’ capabilities without the constraints of real-time processing. This
approach allowed for a quick and detailed evaluation of the models’
performance. Table II provides the performance analysis of the final developed
models.

TABLE II: Performance Metrics for Different RNN Network Architectures Coupling Variables |  Network Arch. |  Train RMSE |  Test #1 RMSE |  Test #2 RMSE |  Test #3 RMSE
---|---|---|---|---|---
|  |  x$10^{-3}$ |  x$10^{-3}$ |  x$10^{-3}$ |  x$10^{-3}$
$x_{mv}$ |  LSTM |  2.90 |  3.71 |  3.26 |  4.05
|  GRU |  3.31 |  4.49 |  4.04 |  4.92
|  LSTM+GRU |  3.31 |  4.36 |  4.00 |  4.84
$x_{m\omega}$ |  LSTM |  50.7 |  67.9 |  62.1 |  72.3
|  GRU |  52.2 |  76.1 |  75.2 |  81.1
|  LSTM+GRU |  51.9 |  75.0 |  70.9 |  72.6
$f_{ev}$ |  LSTM |  15.7 |  21.4 |  19.9 |  31.2
|  GRU |  15.2 |  22.3 |  20.8 |  25.6
|  LSTM+GRU |  15.3 |  21.2 |  24.9 |  29.0
$f_{e\omega}$ |  LSTM |  36.1 |  57.9 |  52.1 |  61.3
|  GRU |  37.4 |  59.1 |  55.2 |  68.1
|  LSTM+GRU |  37.0 |  59.0 |  54.9 |  66.7

The analysis of the error values from Table II indicates that the average
generalization error for $x_{m}$ (i.e., 26% for the LSTM) is consistently
lower than those for $f_{e}$ (i.e., 52% for the LSTM), which might suggest
that the models are better at predicting motion commands than force feedback.
This disparity in prediction accuracy may be attributed to the inherent nature
of the data, where motion commands exhibit a more regular and predictable
pattern compared to the more complex and less predictable patterns of force
feedback. In terms of network architecture, the LSTM consistently outperforms
the GRU and the combined LSTM+GRU in predicting both $x_{m}$ and $f_{e}$, as
evidenced by the lowest error values highlighted in green. Moreover, the LSTM
seems to show the least variation in error rates across the three test
datasets, suggesting better robustness to varying conditions compared to the
other two. Conversely, the GRU exhibits the highest error values, which are
highlighted in brown. Furthermore, Test 2 yields the most accurate predictions
for all variables with the least average generalization error, while Test 3
consistently shows the least accurate performance. This pattern suggests that
the closeness of the test data to the training data could be a determining
factor in the prediction performance of the models. Finally, the best models
(least errors) of the four coupling variables are chosen.

##  IV Human-in-the-loop Experiments

Human-in-the-loop experiments as in [10, 19, 31] were performed on a test
platform to evaluate the performance of teleoperated UGVs on soft terrain with
and without the predictor framework in the presence of large delays.

###  IV-A Test Platform

A real-time human-in-the-loop test platform has been developed, integrating a
Torch-X haptic device with MATLAB/Simulink and the Vortex Studio simulator to
emulate a low-speed teleoperated UGV, as shown in Fig. 4. The choice of Vortex
Studio was based on previous work [14, 16], which demonstrated its high-
fidelity and realistic simulation capabilities for UGVs operating on soft
terrains. Additionally, the dynamic wheel-terrain interactions were
experimentally validated in [27]. This setup can reproduce real-world
longitudinal slippage of the UGVs on soft terrains.

Figure 4: Human-in-the-loop test platform for teleoperated UGV traversing soft
terrain with slippage awareness.

The operator station is equipped with a Torch-X haptic device to generate
motion commands and perceive force feedback signals, as well as a monitor to
provide visual feedback from the virtual environment. The device, which has
three joints, activated the first two for operation while the third is looked
with high gain to encode non-holomorphic driving commands. The haptic device
is connected as a control interface to the simulator via Simulink. To regulate
the observed high-frequency commands and feedback, a first-order low-pass
filter with a cutoff frequency of 0.8 Hz is applied to the coupling variables
to prevent exceeding the predictor framework’s bandwidth of 1.0 Hz.

The virtual environment, assumed as the remote environment, receives motion
commands from the operator station to control the UGV. This UGV is simulated
in Vortex Studio with a comprehensive model that includes a kinematic model,
sensor suite, a Reece tire model, and a deformable terrain model, providing a
dynamic representation of the UGV within the virtual environment. The
deformable terrain property is fine-tuned using empirical lunar soil
parameters from the Mohammed Bin Rashid Space Centre in [27], ensuring the
simulated terrain closely mirrors the mechanical properties of lunar soil.
Moreover, sensors measure the environmental forces in (7), which are
transmitted back to the operator station to provide real-time slippage
awareness.

###  IV-B Experiment Design

A roadmap was established to ensure a diverse array of test tracks, generated
in a virtual environment with the lunar environment in mind, as depicted in
Fig. 5. Three distinct track patterns, labeled A, B, and C (in gray), are
designed as soft terrains, starting from the green mark and ending at the red
mark. Track A is a straight pattern, measuring $10$ meters in length and $2$
meters in width. Track B features right turns, while Track C includes both
right and left turns. Each track is designed to exhibit position-dependent
variable slippage, characterized by an internal friction angle, $\phi(z)$,
which varies as a function of the track’s position, $z$. This function is
defined within the range $0.50\,\text{rad}\leq\phi(z)\leq 0.95\,\text{rad}$,
to introduce potential sand traps [10]. The non-track area (in brown) serves
as a landmark, giving operators a sense of the distance to the turns and the
speed. The safe speed on the tracks is predetermined as a maximum allowable
speed of $0.1\,\text{m/s}$, allowing the UGV to navigate soft terrain without
tipping over. The UGV parameters can be found in [27].

Figure 5: Designated soft terrain tracks A, B, and C; UGV, and non-track zone
in the virtual environment.

Furthermore, the one-way communication delay from the operator station to the
virtual environment is assumed to be 1.25 s, representing the latency for the
long earth-to-lunar distance (380,000 km) network. We employ varying delays
obtained on a real network of $1.0\sim\mathcal{U}(-0.25,0.25)$ using the User
Datagram Protocol to simulate both forward and backward delays for such long
distances. The visual interface simulates lunar low-light conditions, about
0.5 lux, while the simulation editor is set to a lunar gravity of
$1.633\,\text{m/s}^{2}$.

Three distinct cases are introduced under different conditions: (I) Ideal
case, without any delays, serving as a baseline for performance; (II) Delayed
case, with round-trip delays and without a predictor-based framework, to
assess the impact of delays on performance; and (III) Predicted case, with the
same amount of delays and the predictor framework, to assess performance
improvement due to delay compensation. To conduct these experiments, five
operators were selected to perform a track-following task. Before data
collection, these operators received training to familiarize them with the
test platform and practice operating under different conditions—with and
without delays, and with and without the predictor framework. After the
training, each operator repeated each of the three experimental cases, aiming
to complete the track as fast as possible, with excellent command-tracking
performance and without being entrapped. A portion of the collected data on
coupling variables has been allocated for the training and testing phases of
the predictor models in Section III.

###  IV-C Predictor Framework implementation

Each PiLSTM predictor, outlined in Section III, is tailored to a specific
coupling variable. The predictors $\hat{x}_{mv}$ and $\hat{x}_{m\omega}$ are
implemented on the remote environment as the forward predictor, and
$\hat{f}_{ev}$ and $\hat{f}_{e\omega}$ are on the operator station as the
backward predictor, integrated to form the framework shown in Fig. 2. This
framework is designed for optimal delay compensation in bilateral teleoperated
UGVs subjected to terrain-induced slippage [10]. However, seeking such optimal
performance is beyond the scope of this particular application. Finally, the
developed predictor framework and the blended architecture for the
teleoperated UGV are evaluated in two case studies discussed in subsequent
sections.

###  IV-D Evaluation Metrics

The experiment, constrained by operating the teleoperated UGV under ideal,
predicted, and delayed conditions, is evaluated through two case studies, each
with its evaluation metrics. The first case study, an open-loop scenario,
enables comparative analysis of all conditions simultaneously [12], while the
second, a closed-loop scenario, allows for their individual examination [16].
Evaluation metrics are as follows.

####  IV-D1 Delay compensation metrics

In the open-loop case study, the predictors’ delay compensation performance is
evaluated. The output of the ideal case is considered the reference for the
other two cases. The compensation performance in the time domain is quantified
by normalized performance metrics, in percentage, denoted as $\Delta_{n}$
[21]. This metric is computed as the ratio of the L2 norm of the difference
between the predicted output and the ideal output to the L2 norm of the
difference between the delayed output and the ideal output as:

$\Delta_{n}:=\frac{\lVert\hat{x}_{pred}(t)-x_{ideal}(t)\rVert_{2}}{\lVert
x_{del}(t)-x_{ideal}(t)\rVert_{2}}\times 100\%$ (21)

where $x$ symbolizes the coupling variable of interest. A $\Delta_{n}=0\%$
indicates that the predictor perfectly compensates for the delays, achieving
an ideal performance. A $\Delta_{n}<100\%$ suggests that the predictor
compensates the delay to some meaningful amount, with its proximity to zero
reflecting the extent of compensation. Conversely, $\Delta_{n}>100\%$
indicates that the predictor amplifies the negative effect of the delay.

####  IV-D2 Teleoperation metrics

In the closed-loop case study, the teleoperation performance of the predictor
framework is evaluated. The ideal case is considered as the baseline for the
other two cases. Three metrics quantify the teleoperation performance [16]:
performance, transparency, and track completion time. The performance-velocity
tracking $\Omega$ is quantified by the deviations between the commanded motion
and the UGV’s actual measured states. Transparency-force feedback tracking
$\Gamma$, in contrast, is quantified by the deviation between what the
operator perceives and the environmental force-induced slippage. Track
completion time captures how quickly the UGV is operated along the track from
the start point to the end point. The first two metrics directly reflect the
fidelity of the closed-loop integration, where smaller values lead to higher
fidelity, and larger values lead to lower fidelity. The mathematical
expressions for these metrics are [10]:

$\Omega:=\|x_{m}(t)-v_{s}(t)\|_{2},\quad\Gamma:=\|f_{h}(t)-f_{e}(t)\|_{2}$
(22)

However, the force applied by the operator $f_{h}$ cannot be directly
measured. Thus, an estimation technique based on the system dynamics [9, 10,
14, 16] is utilized as:

$f_{h}=\bar{M}_{m}\dot{x}_{m}+\bar{C}_{m}x_{m}-u_{m}$ (23)

##  V Result and Discussion

The proposed teleoperator controller gains in (10) and (12) were set to
$k_{m}=5$ and $k_{s}=1$, to achieve a balanced trade-off between performance
and transparency. The optimal parameter values for the conventional predictors
in (13) are $\alpha=0.57,\beta=1.12$ for the forward predictor, and
$\alpha=0.64,\beta=0.91$ for the backward predictor, as they exhibit different
bandwidths.

###  V-A Case Study 1

This study compares the delay compensation performance between the developed
PiLSTM predictor and the conventional (Conv) predictor frameworks in open-loop
scenarios [21].

A comparative analysis of the coupling variables for Operator #1 is depicted
in Fig. 6. The predicted outputs from both the PiLSTM and the Conv are aligned
with the Ideal output to some extent. However, the PiLSTM outputs are closer
to the ideal outputs than those of the Conv, indicating better delay
compensation by PiLSTM. Notably, the PiLSTM demonstrated a smoother curve,
while the Conv shows more fluctuation, as evidenced by the curves. These
performance differences are clearly observed in the corresponding error
metrics presented in Fig. 7. It can be observed that the PiLSTM maintains a
lower error margin compared to the Conv, with fewer and lower spikes
throughout the observed period. In contrast, the Conv exhibits higher peaks at
specific instances, which could translate to less accurate force feedback to
the operator. This results in errors higher than the delayed outputs, for
instance, in the Conv predicted motion commands $\hat{x}_{mv}$ and
$\hat{x}_{m\omega}$ at 8.2 s and at 23.5 s, respectively, and for the force
feedbacks $\hat{f}_{mv}$ and $\hat{f}_{m\omega}$ at 5.1 s and 15.8 s,
respectively. These peaks typically occur in response to abrupt changes in the
data, as discussed in [21]. Further encapsulating the delay compensation
performance analysis, the normalized performance metrics are summarized in
Table III. Evidently, the PiLSTM predicted output $\hat{x}_{mv}$ demonstrates
a significantly lower error of 38.3% compared to the Conv error of 70.5%.
Similarly, the predicted output $\hat{f}_{mv}$ from PiLSTM shows a lower error
of 39.3% compared to the 78.9% of Conv output.

Similarly, performance comparisons for Operators 2, 3, 4, and 5 were conducted
in the remaining experiments. The comparative analysis of the coupling
variables for Operator #2, along with the corresponding error metrics, is
displayed in Fig. 8 and Fig. 9, respectively. Complementing these visual
representations, the normalized performance metrics for this operator and all
the remaining operators, which encapsulate the delay compensation performance
analysis, are summarized in Table III. The best prediction accuracy for PiLSTM
is 23.6%, compared to the Conv best of 49.4% for $\hat{x}_{mv}$, as
highlighted in green, a difference of 25.8%. Furthermore, the average
performance metric for all variables (framework) for Operator #4 is 32.3% for
PiLSTM, in contrast to Conv’s 58.4%, indicating an improvement of around
26.1%. However, the Conv has the highest error of 116.9% shaded in brown,
indicating that the framework amplifies the negative effect of the delay at
some instances instead of compensating, as highlighted in red.

In conclusion, the analysis conducted across multiple operators in open-loop
scenarios consistently showed that the developed PiLSTM framework outperforms
the Conv framework in delay compensation, with an average improvement of
26.1%, thus demonstrating its enhanced capability for accurate and reliable
bilateral teleoperator controls.

Figure 6: Operator #1: Performance comparison of motion commands (top) and force feedbacks (bottom) among the ideal case, delayed case (1.25 s), conventional predicted case, and PiLSTM predicted case in an open-loop scenario. Our PiLSTM Predictor leads to better delay compensation performance compared to the conventional predictor. Figure 7: Operator #1: Absolute prediction error comparison of motion commands (top) and force feedbacks (bottom) among the ideal case, delayed case (1.25 s), conventional predicted case, and PiLSTM predicted case. Limitations on the delay compensation performance of both predictors occur at specific instances, still, the conventional predictor yields relatively poor delay compensation performance. Figure 8: Operator #2: Performance comparison of motion commands (top) and force feedbacks (bottom) among the delayed case (1.25 s), conventional predicted case, and PiLSTM predicted case in an open-loop scenario. Figure 9: Operator #2: Absolute prediction error comparison of motion commands (top) and force feedbacks (bottom) among the ideal case, delayed case (1.25 s), conventional predicted case, and PiLSTM predicted case. TABLE III: Normalized performance metrics in percentage for open-loop case studies Operators | $\Delta_{n}-x_{mv}(\%)$ | $\Delta_{n}-x_{m\omega}(\%)$ | $\Delta_{n}-f_{ev}(\%)$ | $\Delta_{n}-f_{e\omega}(\%)$
---|---|---|---|---
|  PiLSTM |  Conv |  PiLSTM |  Conv |  PiLSTM |  Conv |  PiLSTM |  Conv
Operator #1 |  38.3 |  70.5 |  42.5 |  68.7 |  39.2 |  78.9 |  44.5 |  89.9
Operator #2 |  28.7 |  58.0 |  37.9 |  72.1 |  33.7 |  97.8 |  35.3 |  102.7
Operator #3 |  54.2 |  86.3 |  50.8 |  116.9 |  49.7 |  88.1 |  58.1 |  98.0
Operator #4 |  23.6 |  49.4 |  34.5 |  53.0 |  33.4 |  69.2 |  41.7 |  78.1
Operator #5 |  52.0 |  108.3 |  72.5 |  82.1 |  57.4 |  79.6 |  74.2 |  91.6

###  V-B Case Study 2

This study evaluates the teleoperation performance of the developed PiLSTM
predictor frameworks, along with the blended architecture of teleoperated UGVs
in closed-loop scenarios [10].

The performance and transparency of the ideal case for operator #1 in the
track-following task are depicted in Fig. 13. The natural deviation of the UGV
states from the given motion command due to induced longitudinal slippage is
well-regulated, as the operator can perceive the induced slippage and apply
the necessary correction through high-fidelity closed-loop integration. As a
result, the command-tracking performance becomes nearly perfect, as evidenced
by the velocity tracking and force feedback tracking metrics, which is
consistent with the results obtained in [9, 16]. However, after introducing a
forward and backward delay of $1.0\sim\mathcal{U}(-0.25,0.25)$ s each in this
ideal case for the same operator—Delayed case, the system’s performance and
transparency significantly deteriorate, as shown by the velocity and force
feedback tracking metrics in Fig. 14. This decline is attributed to the
degradation in the fidelity of the closed-loop integration because the
operator perceives delayed force feedback-induced slippage, which in turn
delays his responses to the UGV. As this effect accumulates, the system
becomes less stable, leading to fast fluctuations, until the operator reduces
their pace and initiates the process again, as evident at around 32 s. Thus
making it difficult for the operator to control the UGV to maintain good
command-tracking performance. This finding is consistent with the results in
[10]. Furthermore, the time taken to accomplish the task increased from 148 s
to 206 s.

Building upon the observed performance deterioration due to the introduced
delays, we subsequently integrated the developed PiLSTM predictor framework to
compensate for the delays—Predicted case. The implementation of this framework
marked a significant improvement in both the performance and transparency of
the system for operator #1, as evident by the velocity and force feedback
tracking metrics presented in Fig. 15. The framework effectively reduced the
discrepancy between the motion commands and the actual state of the UGV, as
well as the environmental force and the force perceived by the operator, by
accurately predicting the undelayed commands and force feedback induced
slippage. The command-tracking performance approached more stable levels
observed in the ideal case, indicating a near recovery to the system’s optimal
state and the restoration of high-fidelity closed-loop integration. This
substantial improvement underscores the efficacy of our delay compensation
framework in restoring system performance. The enhanced time required to
complete the task is now significantly reduced from the delayed scenario to
182 s. Refer to supplementary material for zoomed plots.

Table IV illustrates the effectiveness of the PiLSTM framework in reducing
both velocity (in m/s for $x_{mv}$ and rad/s for $x_{m\omega}$) and force
feedback tracking errors (in N). The Ideal case consistently shows the lowest
error values across all operators, signifying perfect performance. In
contrast, the Delayed case presents relatively higher error values,
highlighting the compromised fidelity of the closed-loop integration. Notably,
the Predicted case shows reduced tracking error in both $x_{mv}$ and
$x_{m\omega}$ compared to the Delayed case, closely approaching the
performance of the Ideal case. This trend of improvement is similarly observed
in the metrics $f_{ev}$ and $f_{v\omega}$, and is consistent across all
operators. However, for certain cases, the PiLSTM framework exhibits large
errors for $f_{v\omega}$ closer to the Delayed case, indicating some
challenges in accurately predicting force feedback. This is likely due to the
less predictable nature of induced slippage. For instance, $f_{v\omega}$ for
operators #2 and #3 are approximately 18.71N and 13.23N, respectively, shaded
in red. Additionally, the prediction error variation for $x_{m}$ is
considerably smaller than for $f_{e}$, showing the framework is more robust in
predicting motion commands than force feedback. Overall, this analysis
suggests the efficacy of the developed framework in compensating for large
delays and subsequently enhancing the fidelity of closed-loop integration.

The bar chart depicted in Fig. 16 compares the track completion times across
the three conditions for the five operators. Notably, the Predicted case
showcases a significant reduction in completion times compared to the Delayed
case, highlighting the practical benefits of the framework in delayed
teleoperated UGV operations. Although the Ideal case remains the best-
performing scenario, the Predicted case narrows this gap, particularly for
operators 1 and 4.

Figure 10: Ideal case - Operator #1: Performance and transparency comparison of motion commands (top) and force feedbacks (bottom), respectively. Baseline performance Figure 11: Delayed case - Operator #1: Performance and transparency comparison of motion commands (top) and force feedbacks (bottom), respectively. With the introduced communication delay of 1.25 s. Figure 12: Predicted case - Operator #1: Performance and transparency comparison of motion commands (top) and force feedbacks (bottom), respectively. With the delays and the PiLSTM predictor framework. TABLE IV: Performance and Transparency Metrics for Closed-Loop Case Studies Operators | Velocity Tracking Error - $\Omega$ | Force Feedback Tracking Error - $\Gamma$
---|---|---
| Ideal Case | Delayed Case | Predicted Case | Ideal Case | Delayed Case | Predicted Case
| $x_{mv}$ | $x_{m\omega}$ | $x_{mv}$ | $x_{m\omega}$ | $x_{mv}$ | $x_{m\omega}$ | $f_{ev}$ | $f_{e\omega}$ | $f_{ev}$ | $f_{e\omega}$ | $f_{ev}$ | $f_{e\omega}$
Operator #1 | 0.023 | 1.66 | 0.083 | 5.17 | 0.030 | 3.13 | 0.89 | 2.84 | 4.09 | 15.41 | 1.94 | 8.44
Operator #2 | 0.020 | 1.89 | 0.077 | 6.52 | 0.054 | 4.91 | 0.91 | 3.11 | 4.17 | 16.21 | 1.08 | 13.23
Operator #3 | 0.031 | 2.01 | 0.11 | 7.61 | 0.074 | 3.96 | 1.12 | 3.17 | 5.53 | 20.72 | 2.79 | 9.65
Operator #4 | 0.019 | 1.05 | 0.072 | 5.09 | 0.039 | 4.42 | 1.09 | 2.56 | 3.89 | 14.56 | 1.26 | 9.31
Operator #5 | 0.030 | 1.96 | 0.12 | 7.52 | 0.081 | 6.17 | 1.14 | 4.09 | 5.99 | 21.62 | 2.84 | 18.71
Figure 13: Ideal case zoomed plot - Operator #1: Performance and transparency
comparison of motion commands (top) and force feedbacks (bottom),
respectively. With perfect velocity and force feedback tracking performance,
the operator was able to well perceive induced slippage through force feedback
and promptly make the necessary corrections to achieve excellent command-
tracking performance. Thus making the system attain high-fidelity closed-loop
integration.  Figure 14: Delayed case zoomed plot - Operator #1: Performance
and transparency comparison of motion commands (top) and force feedbacks
(bottom), respectively. After introducing a communication delay of 1.25 s in
both forward and backward channels, both the velocity and force feedback
tracking become relatively poor, this is because the operator perceived
delayed force feedback, which in turn delayed his responses. As this effect
accumulates, the system becomes less stable, leading to fast fluctuations,
until the operator reduces their pace and initiates the process again. Thus,
poor fidelity closed-loop integration. Figure 15: Predicted case zoomed plot -
Operator #1: Performance and transparency comparison of motion commands (top)
and force feedbacks (bottom), respectively, with the PiLSTM predictor
framework. Integrating the developed framework leads to relatively better
velocity and force feedback tracking, this is because the operator perceived
predicted force feedback, which helped him to respond quickly. Thus, restoring
the fidelity of the closed-loop integration. Figure 16: Track completion time
in delayed, predictor, and ideal cases for all operators. Delays significantly
increase the completion time for all operators, and the PiLSTM predictor
framework significantly improves it compared to the delayed case for operators
1 and 4.

##  VI Conclusion

In this work, a PiLSTM delay compensation framework is presented for low-speed
teleoperated UGVs subjected to terrain-induced slippage and applied in a
human-in-the-loop experiment. Large network delays degrade the fidelity of the
closed-loop integration, making the system less stable and resulting in poor
command tracking performance. We leverage the capabilities of RNNs to learn
complex dynamics from data, thereby improving the prediction performance over
the conventional predictor framework. Due to the inability of purely data-
driven models to extrapolate complex dynamics, we integrate physical
constraints to enhance the predictor’s generalizations. Our results reveal
that the developed PiLSTM predictor framework outperforms the conventional
framework, with an improvement of around 26.1% in open-loop studies. We
further validate the framework in closed-loop studies through human-in-the-
loop experiments, which proves the framework’s efficacy in restoring the
fidelity of the closed-loop integration, leading to a more stable system with
improved command-tracking performance. Not only this, it also succeeded in
reducing task completion times underlining its potential to improve the
traversability and safety of teleoperated UGV for planetary exploration
despite communication delays. Future work may explore a hybrid PiLSTM and
model-based approach for lower prediction error in large-delay compensation.

## References

  * [1] C. Li, Y. Tang, Y. Zheng, P. Jayakumar, and T. Ersal, ”Modeling Human Steering Behavior in Teleoperation of Unmanned Ground Vehicles With Varying Speed,” Hum. Factors, vol. 64, no. 3, 2022.
  * [2] A. Abubakar, K. I. Dahiru, S. H. Sulaiman, and H. Mustapha, “Robust polytopic LPV based Adaptive Cruise Control design for Autonomous Vehicle System,” Int. J. Sci. Eng. Res., vol. 10, no. 12, pp. 1367–1373, 2019.
  * [3] S. N. F. Nahri, S. Du, and B. J. Van Wyk, “A Review on Haptic Bilateral Teleoperation Systems,” J. Intell. Robot. Syst. Theory Appl., vol. 104, no. 1, 2022.
  * [4] S. Opiyo, J. Zhou, E. Mwangi, W. Kai, and I. Sunusi, ”A Review on Teleoperation of Mobile Ground Robots: Architecture and Situation Awareness,” Int. J. Control Autom. Syst., vol. 19, no. 3, pp. 1384-1407, 2021.
  * [5] K. Althoefer, C. P. Tan, Y. H. Zweiri, and L. D. Seneviratne, ”Hybrid Soil Parameter Measurement and Estimation Scheme for Excavation Automation,” IEEE Trans. Instrum. Meas., vol. 58, no. 10, pp. 3633-3641, 2009.
  * [6] R. Almehisni, P. Garg, and M. Wali, ”Thermal Management System Design and Analysis of Rashid Rover - Emirates Lunar Mission,” in Proc. Int. Astronaut. Congr., 2021.
  * [7] W. Li, N. Yang, J. Wang, and D. Ren, ”Kinematic Teleoperation of Wheeled Mobile Robot with Slippage Compensation on Soft Terrains,” IEEE Access, vol. 7, pp. 110982-110991, 2019.
  * [8] A. Abubakar, I. K. Dahiru, S. H. Sulaiman and A. B. Kunya, ”Robust H-Infinity Control for Magnetic Levitation System,” 2019 2nd International Conference of the IEEE Nigeria Computer Chapter (NigeriaComputConf), Zaria, Nigeria, 2019, pp. 1-6.
  * [9] W. Li, J. Guo, L. Ding, J. Wang, H. Gao, and Z. Deng, ”Teleoperation of Wheeled Mobile Robot With Dynamic Longitudinal Slippage,” IEEE Trans. Control Syst. Technol., vol. 31, no. 1, pp. 99-113, 2023.
  * [10] W. Li, L. Ding, H. Gao, and M. Tavakoli, ”Haptic Tele-Driving of Wheeled Mobile Robots under Nonideal Wheel Rolling, Kinematic Control and Communication Time Delay,” IEEE Trans. Syst. Man, Cybern. Syst., vol. 50, no. 1, pp. 336-347, 2020.
  * [11] Y. Zheng, M. J. Brudnak, P. Jayakumar, J. L. Stein, and T. Ersal, ”A Predictor-Based Framework for Delay Compensation in Networked Closed-Loop Systems,” IEEE/ASME Trans. Mechatronics, vol. 23, no. 5, pp. 2482-2493, 2018.
  * [12] Y. Zheng, M. J. Brudnak, P. Jayakumar, J. L. Stein, and T. Ersal, ”A Delay Compensation Framework for Predicting Heading in Teleoperated Ground Vehicles,” IEEE/ASME Trans. Mechatronics, vol. 24, no. 5, pp. 2365-2376, 2019.
  * [13] R. Luz, J. Corujeira, J. L. Silva, and R. Ventura, ”Traction Awareness Through Haptic Feedback for the Teleoperation of UGVs,” in Proc. 27th IEEE Int. Symp. on Robot and Human Interactive Communication, Nanjing, China, pp. 313-318, 2018.
  * [14] W. Li, L. Ding, Z. Liu, W. Wang, H. Gao, and M. Tavakoli, ”Kinematic Bilateral Teledriving of Wheeled Mobile Robots Coupled with Slippage,” IEEE Trans. Ind. Electron., vol. 64, no. 3, pp. 2147-2157, 2017.
  * [15] Y.-J. Pan, C. Canudas-de-Wit, and O. Sename, ”A New Predictive Approach for Bilateral Teleoperation With Applications to Drive-by-Wire Systems,” IEEE Trans. Robot., vol. 22, no. 6, pp. 1146-1162, 2006.
  * [16] W. Li, N. Yang, J. Wang, and Z. Deng, ”Teleoperation of Wheeled Mobile Robots Subject to Longitudinal Slipping and Lateral Sliding by Time-Domain Passivity Controller,” Mechatronics, vol. 81, 2022.
  * [17] A. Abubakar and I. Janajreh, “Comparative study for Active Noise Cancellation using Adaptive filtering and Standing wave pattern,” Int. J. Therm. Environ. Eng., vol. 17, no. 02, pp. 1–11, 2020.
  * [18] D. J. Gorsich, P. Jayakumar, M. P. Cole, C. M. Crean, A. Jain, and T. Ersal, ”Evaluating mobility vs. latency in unmanned ground vehicles,” J. Terramechanics, vol. 80, pp. 11–19, 2018.
  * [19] Y. Zheng, M. J. Brudnak, P. Jayakumar, J. L. Stein, and T. Ersal, ”Evaluation of a Predictor-Based Framework in High-Speed Teleoperated Military UGVs,” IEEE Trans. Human-Machine Syst., vol. 50, no. 6, pp. 561–572, 2020.
  * [20] N. Sridhar, R. Lima, U. Rai, V. Vakharia, K. Das, and B. P., ”Transparency Enhancement in Teleoperation: An Improved Model-Free Predictor for Varying Network Delay in Telerobotic Application,” in 2021 European Control Conference, Delft, Netherlands, 2021.
  * [21] S. Guo, Y. Liu, Y. Zheng, and T. Ersal, ”A Delay Compensation Framework for Connected Testbeds,” IEEE Trans. Syst., Man, Cybern.: Syst., vol. 52, no. 7, pp. 4163-4176, 2022.
  * [22] C. Cruz-Valverde, O. A. Domínguez-Ramírez, E. R. Ponce-de-León-Sánchez, I. Trejo-Mota, and G. Sepúlveda-Cervantes, ”Kinematic and Dynamic Modeling of the PHANToM Premium 1.0 Haptic Device: Experimental Validation,” in 2010 IEEE Electronics, Robotics and Automotive Mechanics Conference, Cuernavaca, Mexico, 2010.
  * [23] M. H. Jeong, T. Y. Lee, S. B. Jeon, and M. Youm, ”Highway Speed Prediction Using Gated Recurrent Unit Neural Networks,” Appl. Sci., vol. 11, no. 7, 2021.
  * [24] Y. Zou, L. Ding, H. Zhang, T. Zhu, and L. Wu, ”Vehicle Acceleration Prediction Based on Machine Learning Models and Driving Behavior Analysis,” MDPI Appl. Sci., vol. 2022.
  * [25] X. A. Ji and G. Orosz, ”Learning Time Delay Systems with Neural Ordinary Differential Equations,” IFAC-PapersOnLine, vol. 55, no. 36, pp. 79–84, 2022.
  * [26] A. Abubakar, M. B. Mohiuddin, O. A. Hay, M. Yakubu, and R. A. Alhammadi, “Robust Control for Autonomous Vehicle Lateral Dynamics: A Comparative Study of Gain-Schedule LPV and Non-linear MPC”.in 8th International Conference on Robotics, Control and Automation. Shanghai, China — January 12-14, 2024.
  * [27] A. Abubakar, R. Alhammadi, Y. Zweiri, and L. Seneviratne, ”Advance Simulation Method for Wheel-Terrain Interactions of Space Rovers: A Case Study on the UAE Rashid Rover,” in 21st International Conference on Advanced Robotics, no. December, pp. 526–532, 2023 Abu Dhabi, UAE, 2023.
  * [28] A. Thomas, A. Manuel, A. Roehrl, and V. Brandtstetter, “Modeling System Dynamics with Physics-Informed Neural Networks Based on Lagrangian Mechanics,” IFAC-PapersOnLine, vol. 53, no. 2, pp. 9195–9200, 2020.
  * [29] F. Djeumou, C. Neary, and E. Goubault, “Neural Networks with Physics-Informed Architectures and Constraints for Dynamical Systems Modeling,” in Learning for Dynamics and Control Conference, pp. 263-277, PMLR, 2022.
  * [30] R. Zhang, Y. Liu, and H. Sun, “Physics-Informed Multi-LSTM Networks for Metamodeling of Nonlinear Structures,” Comput. Methods Appl. Mech. Eng., vol. 369, p. 113226, 2020.
  * [31] P. M. Kebria, A. Khosravi, S. Nahavandi, P. Shi, and R. Alizadehsani, ”Robust Adaptive Control Scheme for Teleoperation Systems With Delay and Uncertainties,” IEEE Trans. Cybern., vol. 50, no. 7, pp. 3243-3253, 2020.