#  $\beta$-Variational autoencoders and transformers for reduced-order
modelling of fluid flows

Alberto Solera-Rico${}^{\text{1,2}}$, Carlos Sanmiguel Vila${}^{\text{1,2}}$,
M. Á. Gómez${}^{\text{2}}$, Yuning Wang${}^{\text{4}}$, Abdulrahman
Almashjary${}^{\text{3}}$, Scott T. M. Dawson${}^{\text{3}}$, Ricardo
Vinuesa${}^{\text{4}}$ 1: Aerospace Engineering Research Group, Universidad
Carlos III de Madrid, Leganés, Spain
2: Subdirectorate General of Terrestrial Systems, Spanish National Institute
for Aerospace Technology (INTA), San Martín de la Vega, Spain
3: Mechanical, Materials, and Aerospace Engineering Department, Illinois
Institute of Technology, Chicago, IL 60616
4: FLOW, Engineering Mechanics, KTH Royal Institute of Technology, SE-100 44
Stockholm, Sweden

###### Abstract

Variational autoencoder (VAE) architectures have the potential to develop
reduced-order models (ROMs) for chaotic fluid flows. We propose a method for
learning compact and near-orthogonal ROMs using a combination of a $\beta$-VAE
and a transformer, tested on numerical data from a two-dimensional viscous
flow in both periodic and chaotic regimes. The $\beta$-VAE is trained to learn
a compact latent representation of the flow velocity, and the transformer is
trained to predict the temporal dynamics in latent space. Using the
$\beta$-VAE to learn disentangled representations in latent-space, we obtain a
more interpretable flow model with features that resemble those observed in
the proper orthogonal decomposition, but with a more efficient representation.
Using Poincaré maps, the results show that our method can capture the
underlying dynamics of the flow outperforming other prediction models. The
proposed method has potential applications in other fields such as weather
forecasting, structural dynamics or biomedical engineering.

## Introduction

Turbulent flows are an important and ubiquitous phenomenon in nature and
engineering, with applications ranging from aircraft design to weather
forecasting. Understanding the behaviour of fluid flows is often challenging
due to their complex spatio-temporal dynamics involving a large number of
degrees of freedom and complex non-linear interactions [1]. As a result, there
is a growing interest in developing reduced-order models (ROMs) of fluid flow
dynamics that can capture the key underlying dynamics of the flow while
reducing the problem dimensionality [2, 3]. Developing ROMs is one of the most
prominent research fields since they facilitate finding low-dimensional
representations that can be applied to perform flow control applications [4]
or reduce the computational cost of numerical simulations [5, 6]. One of the
most used techniques for dimensionality reduction in the fluid-dynamics
community is proper orthogonal decomposition (POD), which involves finding the
dominant modes of variation in a given dataset and projecting the data onto a
lower-dimensional subspace spanned by these modes. Another famous linear
approach is the dynamic-mode decomposition (DMD), which identifies dynamic
modes that govern the evolution of the system over time [7]. While POD and
DMD, as well as the extensions of these methods such as the spectral POD [8]
or the higher-order DMD [9], have successfully reduced the dimensionality of
some flows [10], their optimal linear bases exhibit limitations when working
with turbulent flows, which typically involve complex non-linear interactions
[11].

In recent years, machine-learning (ML) techniques have emerged as promising
approaches for developing ROMs of fluid-flow dynamics [6, 12, 13, 14, 15, 16].
One of the potential ML techniques adopted to create non-linear ROMs is the
neural networks with convolutional autoencoder architectures. These
architectures comprise of both an encoder and a decoder trained to minimize
the reconstruction error between the encoded-decoded data and the initial
data. The resulting encoder allows obtaining a latent-space composed of non-
linear representations. Then, neural-network architectures suitable for
temporal predictions [17, 12], such as long short-term memory (LSTM) networks,
can be used to model the dynamics of the non-linear latent-space, resulting in
a fast surrogate model for fluid-flow predictions.

Among the different autoencoder architectures, variational autoencoders (VAEs)
have proven to be effective for encoding the spatial information of fluid
flows in non-linear low-dimensional latent-spaces [12]. Unlike a standard
autoencoder, a VAE architecture is based on a probabilistic framework for
describing an observation in latent-space by including an additional loss term
on the latent-space variables. However, the low-dimensional representations
obtained using these architectures lack the orthogonality of the classical
linear-decomposition techniques. To overcome this issue, the $\beta$-VAE
architecture introduced in Refs. [18, 19] modifies the loss function of the
VAE by adding a regularisation parameter to balance the reconstruction
accuracy with regularisation and mode disentanglement. The potential of these
architectures to develop a compact and near-orthogonal ROM was reported in
Ref. [20], where they tested this architecture in a high-fidelity simulation
of a turbulent flow through a simplified urban environment. The results showed
a ROM that was able to capture up to $87.36\%$ of the original energy with
only five modes, compared to the $32.41\%$ obtained with five POD modes.

For the temporal predictions, the LSTM has been shown to be an effective
architecture in turbulent flows [17, 21, 22]. However, in the latest years,
another neural-network architecture known as transformer [23] appears to have
the potential to outperform the LSTM and allow the development of more complex
ROMs. The transformer is a deep neural-network architecture that has gained
widespread attention in recent years due to their state-of-the-art performance
in natural-language-processing (NLP) tasks such as language translation [24]
or text generation [25]. Unlike traditional recurrent neural networks like
LSTM, which process sequential data one element at a time, transformers are
designed to capture long-range dependencies between elements in a sequence.
This capability is achieved using self-attention mechanisms which allow the
network to attend to different parts of the input sequence at each network
layer. As a result, transformers have been shown to outperform previous state-
of-the-art methods in several NLP benchmarks, often by large margins.

The success of transformers in NLP has led to their application in other
domains, including computer vision [26, 27], audio processing [28], and
robotics [29]. In particular, due to their ability to capture long-range
dependencies, transformers are particularly well suited to model dynamic
systems [30]. In fact, transformers are able to represent the multi-scale
character of turbulence in long temporal sequences [31]; this can only be
captured by LSTMs when separately predicting modes of different ranges of
frequencies [32]. In these applications, the goal is to learn a low-
dimensional representation of the system that captures the underlying
dynamics, which can then be used to make predictions or generate new
trajectories. Transformers are a promising tool for this task, as they can
learn complex temporal dependencies and capture long-term trends in the data
while allowing efficient parallel processing.

The potential combination of autoencoder architectures, which enable obtaining
near-orthogonal non-linear latent-spaces, with transformer architectures for
the dynamics of the temporal predictions, is a powerful tool that can be
employed to model complex flows with a higher level of accuracy. For this
reason, in this paper, we propose a $\beta$-VAE and transformer-based model
for encoding the fluid-flow velocity fields and learning a ROM of its spatio-
temporal dynamics. Two flow cases are analysed, namely a periodic and a
chaotic configuration of a two-dimensional, viscous flow over two collinear
flat plates, obtained by numerical simulation. Flow past multiple bodies in
close configuration is relevant for a range of applications, such as buildings
or chimneys in urban and industrial environments, power lines, offshore
structures, and heat exchangers. Even at relatively low Reynolds numbers, such
flows can exhibit substantially more complexity than flow over a single,
isolated body [33, 34, 35].

The resulting latent-space of $\beta$-VAEs is analysed using POD as a
reference case to analyse the resulting spatial mode features, and different
latent-spaces are tested. The temporal predictions of the latent-space
dynamics performed using transformer-based architectures are compared with
other ML temporal models, including LSTMs and Koopman with Non-linear Forcing
[36]. Finally, the predictions are assessed using the reconstructed predicted
fields and Poincaré maps to assess the dynamic behaviour of the resulting
ROMs.

## Results

We compare the behaviour of $\beta$-VAE and POD as defined in Methods using
the following metrics to compare their energy reconstruction and orthogonality
performance.

For the $\beta$-VAE, following Ref. [20], the energy percentage $E_{k}$ that
is captured by the model reconstructions is defined as:

$E_{k}=\left(1-\left\langle\frac{\sum_{i=1}^{N_{p}}\sum_{j=1}^{2}(u_{j}-\tilde{u}_{j})^{2}}{\sum_{i=1}^{N_{p}}\sum_{j=1}^{2}u_{j}^{2}}\right\rangle\right)\times
100,$ (1)

where $\langle\cdot\rangle$ indicates ensemble averaging in time, $N_{p}$ is
the number grid points along the spatial domain $u_{j}$ denotes the $j^{th}$
reference value of fluctuating component velocity and $\tilde{u_{j}}$ its
reconstruction, respectively. To measure the orthogonality of the latent
variables, we compute the determinant of the correlation matrix multiplied by
100 and refer to it as $\text{det}_{\textbf{R}}$, where
$\textbf{R}=(\textbf{R})_{dxd}$ is the correlation matrix defined by:

$R_{ii}=1,\quad R_{ij}=\frac{C_{ij}}{\sqrt{C_{ii}C_{jj}}},$ (2)

for all $1\leq i\neq j\leq d$ where $C_{ij}$ denotes the components $i$, $j$
of the covariance matrix C and $d$ is the dimension of the latent-space. A
value of $1$ is reached when all the variables are completely uncorrelated
($R_{ij}=0$) and zero when they are completely correlated ($R_{ij}=1$), this
metric is used to report the independence of the latent variables.

For the POD, the modes are orthogonal by definition and the energy percentage
captured by the POD reconstruction is determined using the cumulative energy
recovery defined as:

$E_{k}=\frac{\sum_{i=1}^{N_{m}}\lambda_{i}}{\sum_{i}\lambda_{i}}\times 100,$
(3)

where $\lambda_{i}$ are the eigenvalues representative of the energy
contribution of each mode.

The modes have been ordered according to their cumulative contribution to the
reconstructed energy, $E_{k}$, following Ref. [20]. The first mode is chosen
as the one with the largest individual contribution to $E_{k}$, and the next
are those that have the maximum contribution when added to the previous modes.

### Analysis of the latent-spaces

In this section, the $\beta$-VAE architecture and the POD are applied to two
flow cases generated from a numerical simulation of a two-dimensional, viscous
flow around two collinear flat plates. The characteristics of both cases are
discussed next:

  * •

Periodic flow with Reynolds number based on the freestream velocity
$U_{\infty}$ and (single-plate) chord length $c$ of $Re=40$, where the two
collinear plates are arranged at an angle of $90^{\circ}$ with respect to the
incoming flow. Total of 1,000 instantaneous flow fields.

  * •

Chaotic flow with Reynolds number of $Re=100$, where the two collinear plates
are arranged at an angle of $80^{\circ}$ with respect to the incoming flow.
Total of 9,095 instantaneous flow fields.

In both cases the domain where data is collected has dimensions $96c\times
28c$, where the spatial resolution is $300\times 98$ grid points with a
uniform grid spacing and a convective time $\Delta t=c/U_{\infty}=t_{c}$
separation between snapshots (downsampled in both space and time from the
original simulations). A more detailed dataset description is given in the
Methods section.

(a)

(b)

Figure 1: Fraction of energy reconstructed by POD as a function of number of
modes. (a) $Re=40$, $\alpha=90^{\circ}$, (b) $Re=100$, $\alpha=80^{\circ}$.

Figure 2: Resulting modes for the $Re=40$, $\alpha=90$ case: a) POD, b)
$\beta$-VAE. The first column contains the streamwise-velocity component modes
sampled with a unit value, and the second one the crosswise-velocity
component. The third column represents the frequency content in the temporal
coefficient associated with each of the modes ($f_{c}=U_{\infty}/c$); note
that in b) this frequency content is evaluated in the latent-space.

The periodic case is used as a benchmark to validate the physical soundness of
the $\beta$-VAE model representations. This flow case can be adequately
represented only using two POD modes that can represent $E_{k}=98.4\%$ of the
kinetic energy of the fluctuations with respect to the base flow, as observed
in Figure 1(a). In Figure 2, the first two rows show the two most energetic
POD modes, which are identified as a shedding wake. For the $\beta$-VAE with
$\beta=0.001$, the results from the last two rows in Figure 2 show the spatial
modes which, as explained previously, are obtained only using the decoder
network with an input vector $\bm{s}_{i}$ containing a unit value in the
desired $i^{th}$ position and zero elsewhere. The reconstructed energy in this
case is equal to $E_{k}=99.4\%$ and $\text{det}_{\textbf{R}}>0.999$. It is
observed from Figure 2 that the spatial modes obtained from both methods
exhibit the same pattern of shedding flow. The spectrum of the temporal
coefficients for both methods is shown in the last column of Figure 2. The
spectra analysis is from the resulting $\textbf{r}_{i}$ and $\textbf{a}_{i}$
time coefficients for the $\beta$-VAE and POD, respectively. This further
confirms that the dynamics associated with the spatial modes are also in good
agreement: both methods can capture the same characteristic frequency. This
result shows that the latent-space also exhibits meaningful physical phenomena
of the flows.

Since the $\beta$-VAE architecture requires the user to set a latent-space
dimension $d$, it can be argued that $d$ can be set to a value larger than 2.
However, we observed that models with larger latent-spaces produce only two
meaningful modes, as the remaining modes have negligible values. This
behaviour shows that the $\beta$-VAE regularisation effectively avoids the
artificial creation of more modes than necessary to represent the solution. In
the work by Eivazi et al. [20] it was shown that the $\beta$-VAE produces
compact representations of latent-spaces, which suggests that these
architectures may be a good framework in cases that can be represented with
few energetic phenomena.

After this first assessment of the latent-spaces created by the $\beta$-VAE
architectures, a non-linear and higher-dimensional chaotic case at $Re=100$
with the two collinear plates arranged at an angle of $80^{\circ}$ with
respect to the incoming flow is tested. As well as the higher Reynolds number,
the change in angle adds additional complexity to the configuration as the
geometry is no longer symmetric with respect to the freestream. In this case,
the size of the latent-space is set to be as small as possible but large
enough to allow the separation of different physical effects in the resulting
modes. In this relatively complex case, it is expected to find more variety of
effects in the fluid flow, requiring a more nuanced mapping to the low-
dimensional latent-space, and an evident lack of performance of linear methods
such as POD.

To define an appropriate $\beta$ value, different values were tested between
0.001 and 0.05 and the $\text{det}_{\textbf{R}}$ and $E_{k}$ were calculated
to evaluate the performance in both metrics with a fixed latent-space. The
value of $\beta$ is chosen high enough to produce a near-orthogonal latent-
space representation but as small as possible to avoid increasing the
reconstruction error. After this study, a $\beta=0.05$ is chosen, and an
analysis of the appropriate latent-space is followed. For this analysis, the
loss terms from the $\beta$-VAE are analysed during the training process as
reported in Figure 3, with a particular interest in the Kullback–Leibler (KL)
divergence loss. The figure shows that, as expected, the reconstruction loss
is lower for the $\beta$-VAE model with larger latent-space $d=20$, as more
information is allowed to flow through the autoencoder bottleneck. The less
constrained latent-space also allows for a lower KL divergence loss, meaning
the latent-space distributions are closer to standard normal distributions.
This lower KL loss also reflects the better ability of the model to produce
disentangled representations in the latent-space. In this case, $E_{k}=98.0\%$
and $\text{det}_{\textbf{R}}=0.944$, indicating both a good reconstruction
capability and low mode correlation.

Figure 3: Evolution of $\beta$-VAE losses during training: (dotted) 5 modes,
(solid) 10 modes, (dashed) 20 modes. Losses are defined in the Methods
section.

On the other hand, the model with a smaller latent-space $d=5$ converges to a
relatively higher reconstruction loss. Even the KL loss term appears unable to
converge during the training process because the bottleneck is too strict. For
$d=5$, the reconstructed energy $E_{k}=86.6\%$ and the determinant of the
correlation matrix $\text{det}_{\textbf{R}}=0.993$, reflect the loss of
reconstruction capability due to the constrained bottleneck. The information
is compressed into a few modes with less freedom for disentangled
representation, a fact that implies that even for these architectures, a
minimum number of modes is required to obtain an appropriate latent-space. The
latent-space of size $d=10$ leads to a slight increase in reconstruction error
compared to the model with $d=20$, with an increased KL loss. However, in this
case, the behaviour of the KL loss with a maximum and a later decay is
indicative of the fact that the latent-space is sufficiently large to obtain
convergence of the KL loss, and the gains in terms of reconstruction loss, in
that case, are modest to justify a double number of modes. For these reasons,
the $\beta$-VAE model with a latent-space of size $d=10$ was chosen as a
reference for the present study. The values for the final model are
$E_{k}=96.4\%$ and $\text{det}_{\textbf{R}}=0.982$, indicating a good
compromise between the reconstruction capability and the expected low
correlation between modes.

Figure 4: Resulting $\beta$-VAE spatial modes for $Re=100$,
$\alpha=80^{\circ}$ ranked according to their contribution to $E_{k}$. The
first column contains $u$ values for modes sampled $\textbf{r}_{i}$. The
second column shows the $v$ values for the same inputs. The third column
represents the frequency content of each individual mode
($f_{c}=U_{\infty}/c$). Figure 5: POD spatial modes for $Re=100$,
$\alpha=80^{\circ}$. The first column contains $u$ values for modes sampled
$r_{i}=1$. The second column shows the $v$ values for the same inputs. The
third column represents the frequency content of each individual mode
($f_{c}=U_{\infty}/c$).

Figure 4 shows, for each mode, the $u$ and $v$ components and the spectrum
obtained from the temporal modes associated with each mode, while Figure 5
shows the same quantities obtained using POD. The spectral analysis of the
temporal dynamics of the modes is used to identify different modes
representing the same physical effect and provides helpful information for
determining the size of the latent-space, since it allows more meaningful
comparison with the POD. In this case, the $\beta$-VAE spatial modes are a
non-linear combination of the most relevant flow features. Comparing the
$\beta$-VAE and POD spectra, it can be seen that the $\beta$-VAE model can
better separate different phenomena, each of which has its characteristic
frequency associated with it, which appears as a peak in the spectrum. The
different effects can be further isolated as individual modes by using larger
latent-spaces. However, it is observed that the frequencies associated with
the most energetic POD modes appear in the $\beta$-VAE dynamics, suggesting
the idea that the $\beta$-VAE latent-spaces can find non-linear modes that
effectively represent the dynamics of the system. The effect of the mode
disentanglement produced by $\beta$-VAE models can be observed in the
correlation matrix between the time series of each mode, Figure 6. As
reported, the correlation between different modes of the time series is almost
negligible due to the near-orthogonal representation in the latent-space. It
is interesting to mention that, as in POD, for the case with $d=20$ analysing
the temporal dynamics of the newly included modes more characteristic
frequencies can be observed, which also appear in the range of the 10 to 20
POD modes. This result suggests that the latent-spaces created with these
methods focus on obtaining the most energetic phenomena.

Figure 6: Correlation matrix corresponding to the $\beta$-VAE mode
coefficients for the case with $Re=100$, $\alpha=80^{\circ}$
$\text{det}_{\textbf{R}}=0.982$. Figure 7: Visualisation of the first
$\beta$-VAE mode with different latent input values. The first column contains
the streamwise-velocity component modes and the second one the crosswise-
velocity component. Case with $Re=100$, $\alpha=80^{\circ}$

In POD orthogonality is ensured, being the mapping between physical and latent
spaces linear. It is well known that chaotic flows often involve complex, non-
linear interactions between fluid particles, which can lead to non-linear
relationships between variables. As a result, the dominant modes of
variability may not capture the complex, non-linear behaviour of the flow.
This effect can be seen in Figure 7, where the first input to the $\beta$-VAE
decoder network is sampled with different scalar values while the remaining
inputs are kept to zero. The mode is sampled with input values in the range
$r_{i}\in[-2,2]$ because during VAE training the latent-space is sampled from
a near-standard normal distribution, driven by the KL-loss regularisation
during training. The non-linear representation is evident as the shedding wake
patterns change with the latent input value. This non-linear representation
allows the $\beta$-VAE architecture to reproduce $E_{k}=96.42\%$ with only ten
modes. In contrast, POD only yields $E_{k}=48.79\%$ and would require using
over 100 modes to obtain a reconstruction above $E_{k}=95\%$.

### Latent-space predictor models

(a)

(b)

Figure 8: (a) Reference values and predictions of latent-space modes by
different models. (b) Average prediction error over the prediction horizon.
This is shown for the case with $Re=100$, $\alpha=80^{\circ}$.

In this section, the temporal dynamics of the latent-space $\bm{r}_{i}$ for
the case $d=10$ are combined with different ML models to implement a framework
able to predict the temporal dynamics in the latent-space that can be used
with the decoder from the $\beta$-VAE to obtain flow-field temporal
predictions. With this purpose, a transformer model [23], a Koopman with Non-
linear Forcing (KNF) model [36], and an LSTM network [37] are implemented and
compared. The KNF and LSTM models have previously been analysed and compared
in Ref. [22] and applied to predict the temporal dynamics of a low-order model
of near-wall turbulence, showing that both approaches can reproduce the
temporal dynamics of this system. Furthermore, the transformer has been used
in the context of temporal predictions of turbulent flows in Ref. [31]. The
three model architectures have been tuned to obtain the lowest mean-squared
error over the validation data.

All three predictor models are trained to predict the next time step of a
temporal sequence of previous latent-space $\bm{r}_{i}$ vectors. Further
details of the process and hyperparameter choice are detailed in the Methods
section. Once trained, the models are used recursively to predict the next
time steps. By predicting long time series, we can determine whether the model
has been able to learn the system dynamics correctly. Figure 8(a) shows the
true $\bm{r}_{i}$ values for the test data and the corresponding $\bm{p}_{i}$
predictions from each model. As expected for a chaotic dynamical system, all
predictor models diverge from the original trajectory after several time steps
and appear to exhibit similar quantitative performance. To further evaluate
the deviation error after several time steps, the ensemble average of the L2
error norm to the prediction horizon is defined as:

$\varepsilon(\Delta
t)=\left\langle\left(\sum_{i=0}^{d}\Big{(}p_{i}(t)-r_{i}(t)\Big{)}^{2}\right)^{1/2}\right\rangle_{\rm{ensemble}}.$
(4)

This quantity is shown in Figure 8(b) and the ensemble is computed over all
the possible windows in the test data to capture the average trend. This plot
shows how the error $\varepsilon(\Delta t)$ increases as the prediction
horizon $\Delta t$ gets longer. It can be seen that the error growth of the
LSTM and KNF models is slower, which at first sight would suggest a better
prediction by these models. However, since all predictions tend to diverge
after several time steps, this result alone does not imply that the LSTM or
KNF models would properly capture the system dynamics. Therefore, this aspect
is further evaluated in terms of the dynamic behaviour of the system through
the Poincaré-map analysis using the predictions for 1,024 additional time
steps not seen during the training process.

(a) LSTM

(b) KNF

(c) Transformer

Figure 9: Poincaré sections with plane $r_{0}=0$ and
$\text{d}r_{0}/\text{d}t>0$. Black lines correspond to the original data and
blue lines to the predicted data.

The Poincaré map is constructed as the intersection of the latent vectors with
the hyperplane $r_{1}=0$ on the $r_{6}-r_{10}$ space with direction
$\text{d}r_{1}/\text{d}t>0$, using then a probability density function (PDF)
to fit the resulting intersection points. These distributions are plotted for
the true data and the model predictions in Figure 9. This figure shows the
correlation between the amplitudes of the $r_{5}$ and $r_{9}$ latent variables
at the intersection. It can be seen from these figures that the LSTM and, to a
lesser extent, the KNF models do not adequately reproduce the variability of
this correlation. This result can be explained by the tendency of the LSTM
model to converge towards a harmonic behaviour as the prediction progresses,
which does not capture the variability in the evolution of the temporal
coefficients that represent this chaotic flow. The transformer produces
accurate predictions of the dynamics of the system using the patterns and
correlations between the different $r_{i}$ in the latent-space learned by the
$\beta$-VAE. This comparison between models suggests that focusing only on
instantaneous predictions may not be the correct approach to develop ROM
models based on ML techniques.

Figure 10: True, reconstructed and predicted fields for $u$ and $v$ velocity
components. Top two rows show the prediction of the next time step
$(t+t_{c})$, while the bottom rows show the prediction of the snapshot at
$t+5t_{c}$ time step.

Using the latent vectors predicted with the transformer, we can develop a ROM
of the flow in the latent-space that captures the underlying dynamics, as
shown by the Poincaré maps, and then using these predicted $r_{i}$ values the
$\beta$-VAE decoder network can be used to project the predicted latent-space
vector back into physical space, completing the flow field prediction model.
The results of this prediction are shown in Figure 10. In the figure, the
first column shows the actual $t+n$ flow field in the dataset, the second
column shows the same snapshot compressed into latent-space and decompressed
by the $\beta$-VAE network, and the last column shows the field predicted by
the ROM created using the $\beta$-VAE and transformer. The first two rows
correspond to a prediction of the next time step, while the bottom two rows
correspond to the result of the recurrent application of the transformer model
to predict the next 5 time steps, being the time step equal to a one
convective time $t_{c}=c/U_{\infty}$. The calculation of the reconstructed
energy of the predicted fields for time steps $t+t_{c}$ and $t+5t_{c}$ yields
$E_{k}$ values of 95.86% and 94.96%, respectively. These results support the
quality of the reconstruction observed in Figure 10, even after several time
steps equal to the convective time, and confirm the robustness of the
methodology proposed in this work.

## Discussion

This study presents and evaluates a ROM framework based on $\beta$-VAE
architectures to produce robust non-linear latent-spaces, combined with the
time-prediction model obtained by means of a transformer architecture. Using a
two-dimensional viscous flow around two collinear flat plates in periodic and
chaotic regimes, a first analysis of the $\beta$-VAE capabilities and the
latent-spaces generated using these techniques is assessed and compared with
the ROM obtained through POD modes. For the periodic case, it is observed that
the same modal features are obtained, representing a vortex shedding with an
equivalent number of energetic modes. For the chaotic case, the $\beta$-VAE
learns a compact near-orthogonal latent-space that significantly improves the
energy reconstruction obtained by the POD using ten modes, $E_{k}=95.3\%$
against $E_{k}=48.79\%$. While the $\beta$-VAE ROM comprises 10 modes, a total
of 100 POD modes would be required to match the result obtained by the
$\beta$-VAE (in terms of flow reconstruction). The resulting modes from the
$\beta$-VAE analysis exhibit temporal dynamics with shedding frequencies
equivalent to those observed in the most energetic POD modes suggesting that
the $\beta$-VAE is able to obtain non-linear modes that represent the most
energetic or representative flow features.

The latent-space generated by the $\beta$-VAE is combined with a transformer
architecture to predict the temporal dynamics. Its performance is compared
against other models used in previous studies, such as LSTM or KNF. The
results show that despite the fact that the prediction error is similar for
all the methods, only the transformer can learn the correlations among the
various temporal coefficients. Combining the $\beta$-VAE and the transformer
models, we obtain a ROM model that can produce predictions with $E_{k}$ equal
to 95.86% and 94.96% at $t+t_{c}$ and $t+5t_{c}$. This study constitutes a
proof of concept of the capabilities of these novel ML techniques to generate
compact and robust non-linear ROMs that can be employed to generate
predictions in chaotic flows while capturing the most relevant flow features.
The combination of various VAE and transformer architectures exhibits great
potential for the future of ROM development in fluid mechanics, leveraging the
inherent nonlinearities of these techniques.

## Methods

###  $\beta$-VAE implementation details

The variational autoencoder (VAE) described in Kingma and Welling [38] is one
of the most common architectures used in generative models. In the basic
autoencoder architecture, the input x is directly encoded as a vector r using
the encoder $\mathcal{E}$, $\textbf{r}=\mathcal{E}(\textbf{x})$ in a lower-
dimensional space of size $d$, which can be mapped back to the original space
by the decoder $\tilde{\textbf{x}}=\mathcal{D}(\textbf{r})$, being
$\tilde{\textbf{x}}$ the reconstructed output. A disadvantage of conventional
autoencoders is the fact that there are no constraints imposed on the learned
latent-space. The lack of constraints in the generated latent-space can lead
to overfitting and poor generalisation performance, especially for high-
dimensional input data. In contrast, in VAE architectures the encoder maps
input data to a parameterised prior distribution, usually a Gaussian
distribution, in the latent-space
$\bm{\mu},\bm{\sigma}=\mathcal{E}(\textbf{x})$, where $\bm{\mu}$ and
$\bm{\sigma}$ denote the mean and the standard deviation, respectively. The
prior Gaussian distribution encourages the model to learn a compact and smooth
representation of the input data. The distribution of the latent-space is then
randomly sampled, and this sample is decoded back into the original space by
the decoder $\tilde{\textbf{x}}=\mathcal{D}(\textbf{s})$, with s following a
distribution $\textbf{s}\sim N(\bm{\mu},\bm{\sigma})$. The architecture can
therefore produce different outputs for the same input data, each sampled from
the corresponding latent distribution. The encoder and decoder neural networks
are trained by gradient descent and backpropagation, enabled by a re-
parameterisation of the distribution sampling [38]. The training process for
the VAE involves simultaneous training of the encoder and decoder networks
with a compound loss function $\mathcal{L}$,

$\mathcal{L}(\bm{x})=\underbrace{\frac{1}{N_{t}}\sum_{i=1}^{N_{t}}(\bm{x}-\tilde{\bm{x}})^{2}}_{\text{Reconstruction
loss}}-\underbrace{\frac{1}{2}\sum^{d}_{i=1}\big{(}1+\log(\sigma_{i}^{2})-\mu_{i}^{2}-\sigma_{i}^{2}\big{)}}_{\text{KL
loss}},$ (5)

where $N_{t}$ denotes the total number of points the in x. The first term in
the loss function is the reconstruction loss, $\mathcal{L}_{rec}$, which
measures the model accuracy in reconstructing the input data from the reduced
latent-space representation. The second term in the loss function is the
Kullback–Leibler (KL) divergence loss [39], $\mathcal{L}_{KL}$, which measures
the difference between the generated probability distribution and a prior
probability distribution, typically Gaussian. The overall training goal is to
optimise the model to produce accurate reconstructions while keeping the
latent distributions close to a standard normal distribution.

The original VAE architecture was extended by Higgins et al. [40] to the
$\beta$-variational autoencoder ($\beta$-VAE) architecture, which aims to
promote disentangled representations in latent-space [18]. The $\beta$-VAE
loss function, defined in equation (6), includes a scalar hyperparameter
$\beta\geq 0$ that modulates the trade-off between reconstruction accuracy and
latent-space disentanglement. A higher value of $\beta$ results in a more
disentangled representation but may decrease reconstruction accuracy, whereas
lower value of $\beta$ may result in a more accurate reconstruction but a less
disentangled latent-space. As a result, by incorporating the $\beta$ parameter
into the loss function, the $\beta$-VAE architecture can achieve a more
interpretable and disentangled latent-space, which may improve the model
generalisation ability.

$\mathcal{L}(\bm{x})=\mathcal{L}_{rec}-\frac{\beta}{2}\mathcal{L}_{KL}.$ (6)

The $\beta$-VAE architecture is adapted from Ref. [20]. The encoder is used to
generate the temporal modes for the latent-space representation using the mean
vector for each time step $\textbf{r}_{t}=\bm{\mu}_{t}$. The mean vectors are
generated by applying the encoder $\mathcal{E}$ to each input timestep:
$\bm{\mu}_{t},\bm{\sigma}_{t}=\mathcal{E}(\bm{x}_{t})$. The spatial modes are
reconstructed one by one: $Y_{i}=\mathcal{D}(\bm{s}_{i})$. Here $Y_{i}$ is the
$i^{th}$ spatial mode, $\mathcal{D}$ is the decoder and $\bm{s}_{i}$ is the
vector that selects the mode to be reconstructed
($\bm{s}_{i}=\delta_{i}^{j}=(0,\cdots,1,\cdots,0)$). The $\beta$-VAE
architecture is sketched in Fig. 11.

Figure 11: Model architectures: (top) $\beta$-VAE encoder, (middle)
$\beta$-VAE decoder and (bottom) transformer.

The nature of the data determines the choice of a convolutional neural network
(CNN) to build the encoder, since the patterns in the flow can be better
captured by convolutional layers that preserve the spatial relationships
between points in the input. Figure 11 shows a representation of the model. In
the encoder, we use six convolutional layers with a stride of two so that each
layer halves the spatial dimension. The spatial reduction allows the
subsequent layers to capture information at larger scales in the input flow
data, which also has similarities to the multiple scales of the studied
chaotic flow and can help to represent it. The number of filters in each layer
progressively increases to preserve the information flow while reducing the
spatial dimension. After the sixth convolutional layer, the spatial
information is discarded, and a fully-connected layer is added to combine the
information. Finally, two parallel layers with the same number of units as the
latent-space dimension are used to output the mean and variance of the latent
statistical distributions. During training, the distributions are sampled to
generate inputs to the decoder network, while only the $\bm{\mu}$ values are
used to encode the time series used later by the predictor.

The decoder model is designed as an almost symmetric network to the encoder.
Latent-space samples are fed into a fully-connected layer, and its output is
reshaped to the same shape as the last convolutional layer in the encoder. Six
transposed convolution layers are then used to increase the spatial dimension
with decreasing filters. A final transposed convolution layer with two filters
produces the two output channels. The chosen activation function is the
exponential linear unit (ELU) [41] for all layers except the last, where the
activation is linear.

### Proper orthogonal decomposition

In this section POD is employed as a reference to compare the resulting modes
from the $\beta$-VAE latent-space with those obtained using this classical
method. In particular, the snapshot method [42] has been used for the present
study. Considering the streamwise and crosswise velocity components defined as
$U(\textbf{x},t)$ and $V(\textbf{x},t)$, respectively, being
$\textbf{x}=(x,y)$ with $x$ and $y$ the streamwise and crosswise coordinates
and $t$ the time. We decompose the velocity components as:

$U(\textbf{x},t)=\overline{U}(\textbf{x})+u(\textbf{x},t),\quad
V(\textbf{x},t)=\overline{V}(\textbf{x})+v(\textbf{x},t)$ (7)

where $\overline{U}(\textbf{x})$ and $\overline{V}(\textbf{x})$ are the
streamwise and crosswise time-averaged velocity components and
$u(\textbf{x},t)$ and $v(\textbf{x},t)$ are the streamwise and crosswise
fluctuating velocity components. The fluctuating quantities can be
approximated as a linear combination of basis functions $\phi_{i}(\textbf{x})$
as:

$u(\textbf{x},t)\approx\sum_{i=1}^{N_{m}}a_{i}^{u}(t)\phi_{i}^{u}(\textbf{x}),\quad
v(\textbf{x},t)\approx\sum_{i=1}^{N_{m}}a_{i}^{v}(t)\phi_{i}^{v}(\textbf{x}),$
(8)

where $a_{i}(t)$ are time-dependent coefficients and $N_{m}$ is the number of
basis functions used. Here we assume a number of images $N_{t}$, each one
consisting of $N_{p}$ grid points along the spatial domain x, with
$N_{t}<N_{p}$. Following the snapshot method [42], each image can be treated
as an $N_{p}$-dimensional vector and the data can be arranged into an
$N_{t}\times N_{p}$ snapshot matrix:

$\textbf{u}=\left[\begin{array}[]{ccccc}u(x_{1},t_{1})&\cdots&u(x_{N_{p}},t_{1})\\\
\vdots&\ddots&\vdots\\\ u(x_{1},t_{N_{t}})&\cdots&u(x_{N_{p}},t_{N_{t}})\\\
\end{array}\right];\
\textbf{v}=\left[\begin{array}[]{ccccc}v(x_{1},t_{1})&\cdots&v(x_{N_{p}},t_{1})\\\
\vdots&\ddots&\vdots\\\ v(x_{1},t_{N_{t}})&\cdots&v(x_{N_{p}},t_{N_{t}})\\\
\end{array}\right].$ (9)

Using the snapshot matrix, the two-point correlation matrix can be written as
$\textbf{G}=\textbf{u}\ \textbf{u}^{T}+\textbf{v}\ \textbf{v}^{T}$, where the
superscript T refers to the matrix transpose. Solving the eigenvalue problem
of G returns the eigenvalues $\lambda_{i}$ and the left and right eigenvector
matrices. The left and right eigenvector matrices are respectively the matrix
$\bm{\Psi}$ containing in its columns the temporal modes $\textbf{a}_{i}$
(which are orthonormal vectors of length $N_{t}$ and unitary norm) and its
inverse (i.e. its transpose). Note that the columns of $\bm{\Psi}$ form a
basis of rank $N_{t}$ and that the eigenvalues $\lambda_{i}$ are
representative of the energy contribution of each mode. The orthonormal
spatial modes $\phi_{n}(\textbf{x})$ can then easily be computed as
$\bm{\Sigma_{u}}\ \bm{\phi_{u}}=\bm{\Psi}^{T}\ \textbf{u}$ and
$\bm{\Sigma_{v}}\ \bm{\phi_{v}}=\bm{\Psi}^{T}\ \textbf{v}$, where
$\bm{\Sigma_{u}}$ and $\bm{\Sigma_{v}}$ are diagonal matrices which, in each
$i^{th}$ diagonal element, contain the streamwise and wall-normal Reynolds-
stress contributions of the $i^{th}$ mode.

### Time-series prediction models

The $\beta$-VAE encoder network is used to generate the entire dataset latent-
space time series. The last 1024 time steps are then used as test data and the
rest is used for training. The transformer model predicts the next latent-
space vector in a sequence of latent vectors generated by the $\beta$-VAE. The
input to the transformer at each time step is a sequence of 32 time steps, and
the output is the prediction of the next time step for the latent vector in
the sequence. The transformer is trained to minimise the difference between
the prediction of the next time step and the true data using a mean-squared-
error loss function.

A time-embedding vector [43] is added to each input latent-space vector to
incorporate temporal information into the model. The time-embedding vector
represents the time step at which the latent vector was generated, allowing
the model to distinguish between latent vectors generated at different time
steps and to capture the temporal dynamics of the physical system. The time
embedding is a vector which is learned together with the other parameters of
the transformer model during training. It is randomly initialised at the start
of training and updated during backpropagation in the same way as the other
parameters of the model.

The transformer model comprises a set of encoder and decoder layers, each
consisting of a multi-head attention block and a feed-forward neural network.
The attention blocks allow the model to weight the importance of different
parts of the input sequence when making predictions [44], while the feed-
forward network allows it to learn complex non-linear relationships between
the input and output sequences. Two transformer blocks are concatenated with
skip connections around them to enhance gradient propagation towards the input
during the backpropagation phase of training [45]. After the transformer
blocks, feed-forward layers learn the complex non-linear relationships in the
transformer output and form the final latent-space vector prediction. The
architecture is illustrated in Figure 11. To tune the number of blocks, layers
and neurons per layer among other hyperparameters, the KerasTuner framework
[46] was employed using a Bayesian-optimisation algorithm to evaluate
different models, using the validation error as the reference metric. The LSTM
model architecture was also tuned searching the lowest prediction error in the
test data. The resulting architecture includes two layers of 128 LSTM
elements, followed by a fully-connected layer of 128 neurons and a final
output layer matching the latent-space size. The input for the LSTM includes
the previous 16 time steps, the same as for the transformer model. The KNF-
model implementation is based in the code from Ref. [22]. After the
hyperparameter tuning, the number of previous time steps used to predict is 5,
and the maximum order of the functions for construction of the forcing term is
3 for polynomial functions and 4 for trigonometric functions.

### Training setup

The TensorFlow 2.10 deep-learning framework [47] was used to implement the
models and the training pipeline. An NVIDIA A100-SXM GPU was used to train the
models within the Google Colaboratory framework. Training the $\beta$-VAE
model took about 1 hour while training the transformer model took only a few
minutes. The $\beta$-VAE model is trained using the Adam algorithm [48] with a
lookahead mechanism [49] and gradient clipping of 0.5 for stability
requirements [50]. The synchronisation period for the lookahead mechanism is
6, and the slow step size is 0.4. The learning rate is variable, starting at
5$\times$10-4 and decreasing with a piecewise constant decay to 3$\times$10-5.
The model is trained over 1000 epochs. The encoder network has 3.61$\times$106
trainable parameters, while the decoder network has 3.48$\times$106. The
entire dataset is used for training the $\beta$-VAE model. The resulting
temporal dynamics from the $\beta$-VAE are then used to train a transformer
architecture using the NovoGrad algorithm [51] with a gradient-clipping value
of 0.5 for stability reasons. The learning rate was set to 0.001, and the
model was trained for 256 epochs. The number of trainable parameters for this
model is 2.41$\times$105. In this case, the last 1024 time steps are not used
during training and are used to test the prediction models.

### Dataset description and pre-processing

The source data set is an incompressible, two-dimensional, viscous flow over
two collinear flat plates. Numerical simulation is used to solve the governing
Navier–Stokes equations for this flow using an immersed-boundary projection
method (IBPM) [52, 53]. The two plates have a chord length $c$ and are
separated by a gap of the same size. The free stream velocity is $U_{\infty}$,
and the Reynolds number, $Re$, is defined based on the free stream velocity
and single-plate chord length.

The simulations are performed on a series of nested grids of increasing size
and decreasing resolution. The finest resolution for the simulation has a grid
spacing of $\Delta x=\Delta y=0.02c$, and the total computational domain has
dimensions of $96c\times 28c$, with the upstream boundary at $x=-9c$. The
dataset generation and its characteristics are described in more detail in
Ref. [54]. The time step in the dataset (downsampled from the simulations) was
chosen to be equal to the convective time $\Delta t=c/U_{\infty}=t_{c}$, which
is sufficient to resolve the dynamics of coherent vortical structures in the
wake.

The dataset employed for this study is downsampled from the original mesh to a
spatial resolution of $300\times 98$ with a uniform grid spacing, to reduce
the GPU requirements for training. The fluctuation components of the
streamwise $u$ and crosswise $v$ velocity were stacked as separate channels in
the dataset. The values in the dataset were standardised by subtracting the
pixel-wise average and dividing by the standard deviation for each component
over the entire dataset. The resulting dataset size was reduced to 1.8 GB and
fully loaded into memory during training. This approach allowed the training
process to be efficient while still capturing the essential features of the
fluid flow.

## Data availability

The dataset will be made publicly available after publication.

## Code availability

The codes will be made publicly available after publication.

## Acknowledgements

The authors would like to thank to Dr. Hamidreza Eivazi and Mr. Samuel Molina
for their technical assistance in the first part of the study and Dr. Steve
Brunton for his insightful suggestions. Some of the deep-learning models were
trained by means of resources provided by the National Academic Infrastructure
for Supercomputing in Sweden (NAISS). RV acknowledges the financial support
from ERC grant no. ‘2021-CoG-101043998, DEEPCONTROL’. Views and opinions
expressed are however those of the author(s) only and do not necessarily
reflect those of the European Union or the European Research Council. Neither
the European Union nor the granting authority can be held responsible for
them.

## Author contributions

ASR: Methodology, Software, Validation, Investigation, Data Curation, Writing
- Original Draft, Writing - Review & Editing, Visualization. CSV: Methodology,
Software, Validation, Investigation, Data Curation, Writing - Original Draft,
Writing - Review & Editing, Visualization. MAG: Methodology, Writing - Review
& Editing. YW: Methodology, Software, Validation, Writing - Review & Editing.
AA: Data generation. SD: Data generation, Writing - Original Draft, Writing -
Review & Editing RV: Ideation, Methodology, Writing - Review & Editing,
Resources, Funding acquisition.

## Competing interests

The authors declare that they have no competing interests.

## References

## References

  * Jiménez [2018] J. Jiménez, Coherent structures in wall-bounded turbulence, Journal of Fluid Mechanics 842, P1 (2018).
  * Taira  _et al._ [2017] K. Taira, S. L. Brunton, S. T. Dawson, C. W. Rowley, T. Colonius, B. J. McKeon, O. T. Schmidt, S. Gordeyev, V. Theofilis, and L. S. Ukeiley, Modal analysis of fluid flows: An overview, Aiaa Journal 55, 4013 (2017).
  * Taira  _et al._ [2020] K. Taira, M. S. Hemati, S. L. Brunton, Y. Sun, K. Duraisamy, S. Bagheri, S. T. Dawson, and C.-A. Yeh, Modal analysis of fluid flows: Applications and outlook, AIAA journal 58, 998 (2020).
  * Rowley and Dawson [2017] C. W. Rowley and S. T. Dawson, Model reduction for flow analysis and control, Annual Review of Fluid Mechanics 49, 387 (2017).
  * Lucia  _et al._ [2004] D. J. Lucia, P. S. Beran, and W. A. Silva, Reduced-order modeling: new approaches for computational physics, Progress in Aerospace Sciences 40, 51 (2004).
  * Brunton  _et al._ [2020] S. L. Brunton, B. R. Noack, and P. Koumoutsakos, Machine learning for fluid mechanics, Annu. Rev. Fluid Mech. 52, 477 (2020).
  * Schmid [2010] P. J. Schmid, Dynamic mode decomposition of numerical and experimental data, Journal of fluid mechanics 656, 5 (2010).
  * Sieber  _et al._ [2016] M. Sieber, C. O. Paschereit, and K. Oberleithner, Spectral proper orthogonal decomposition, Journal of Fluid Mechanics 792, 798 (2016).
  * Le Clainche and Vega [2017] S. Le Clainche and J. M. Vega, Higher order dynamic mode decomposition, SIAM Journal on Applied Dynamical Systems 16, 882 (2017).
  * Abreu  _et al._ [2020] L. I. Abreu, A. V. Cavalieri, P. Schlatter, R. Vinuesa, and D. S. Henningson, Spectral proper orthogonal decomposition and resolvent analysis of near-wall coherent structures in turbulent pipe flows, Journal of Fluid Mechanics 900, A11 (2020).
  * Berkooz  _et al._ [1993] G. Berkooz, P. Holmes, and J. L. Lumley, The proper orthogonal decomposition in the analysis of turbulent flows, Annual review of fluid mechanics 25, 539 (1993).
  * Eivazi  _et al._ [2020] H. Eivazi, H. Veisi, M. H. Naderi, and V. Esfahanian, Deep neural networks for nonlinear model order reduction of unsteady flows, Physics of Fluids 32, 105104 (2020).
  * Hijazi  _et al._ [2020] S. Hijazi, G. Stabile, A. Mola, and G. Rozza, Data-driven pod-galerkin reduced order model for turbulent flows, Journal of Computational Physics 416, 109513 (2020).
  * Murata  _et al._ [2020] T. Murata, K. Fukami, and K. Fukagata, Nonlinear mode decomposition with convolutional neural networks for fluid dynamics, Journal of Fluid Mechanics 882, A13 (2020).
  * Vinuesa and Brunton [2022] R. Vinuesa and S. L. Brunton, Enhancing computational fluid dynamics with machine learning, Nature Computational Science 2, 358 (2022).
  * Fukami  _et al._ [2020] K. Fukami, T. Nakamura, and K. Fukagata, Convolutional neural network based hierarchical autoencoder for nonlinear mode decomposition of fluid field data, Physics of Fluids 32, 095110 (2020).
  * Srinivasan  _et al._ [2019] P. A. Srinivasan, L. Guastoni, H. Azizpour, P. Schlatter, and R. Vinuesa, Predictions of turbulent shear flows using deep neural networks, Physical Review Fluids 4, 054603 (2019).
  * Burgess  _et al._ [2018] C. P. Burgess, I. Higgins, A. Pal, L. Matthey, N. Watters, G. Desjardins, and A. Lerchner, Understanding disentangling in $\beta$-VAE, arXiv preprint arXiv:1804.03599 (2018).
  * Higgins  _et al._ [2017a] I. Higgins, L. Matthey, A. Pal, C. Burgess, X. Glorot, M. Botvinick, S. Mohamed, and A. Lerchner, $\beta$-vae: Learning basic visual concepts with a constrained variational framework, in  _International conference on learning representations_ (2017).
  * Eivazi  _et al._ [2022] H. Eivazi, S. Le Clainche, S. Hoyas, and R. Vinuesa, Towards extraction of orthogonal and parsimonious non-linear modes from turbulent flows, Expert Systems with Applications 202, 117038 (2022).
  * Nakamura  _et al._ [2021] T. Nakamura, K. Fukami, K. Hasegawa, Y. Nabae, and K. Fukagata, Convolutional neural network and long short-term memory based reduced order surrogate for minimal turbulent channel flow, Physics of Fluids 33, 025116 (2021).
  * Eivazi  _et al._ [2021] H. Eivazi, L. Guastoni, P. Schlatter, H. Azizpour, and R. Vinuesa, Recurrent neural networks and Koopman-based frameworks for temporal predictions in a low-order model of turbulence, International Journal of Heat and Fluid Flow 90, 108816 (2021).
  * Vaswani  _et al._ [2017] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. u. Kaiser, and I. Polosukhin, Attention is all you need, in  _Advances in Neural Information Processing Systems_, Vol. 30, edited by I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett (Curran Associates, Inc., 2017).
  * Yang  _et al._ [2020] S. Yang, Y. Wang, and X. Chu, A survey of deep learning techniques for neural machine translation, arXiv preprint arXiv:2002.07526 (2020).
  * Brown  _et al._ [2020] T. Brown, B. Mann, N. Ryder, M. Subbiah, J. D. Kaplan, P. Dhariwal, A. Neelakantan, P. Shyam, G. Sastry, A. Askell, _et al._ , Language models are few-shot learners, Advances in neural information processing systems 33, 1877 (2020).
  * Khan  _et al._ [2022] S. Khan, M. Naseer, M. Hayat, S. W. Zamir, F. S. Khan, and M. Shah, Transformers in vision: A survey, ACM computing surveys (CSUR) 54, 1 (2022).
  * He  _et al._ [2020] S. He, W. Liao, H. R. Tavakoli, M. Yang, B. Rosenhahn, and N. Pugeault, Image captioning through image transformer, in  _Proceedings of the Asian Conference on Computer Vision (ACCV)_ (2020).
  * Gong  _et al._ [2021] Y. Gong, Y.-A. Chung, and J. Glass, Ast: Audio spectrogram transformer, arXiv preprint arXiv:2104.01778 (2021).
  * Prakash  _et al._ [2021] A. Prakash, K. Chitta, and A. Geiger, Multi-modal fusion transformer for end-to-end autonomous driving, in  _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ (2021) pp. 7077–7087.
  * Geneva and Zabaras [2022] N. Geneva and N. Zabaras, Transformers for modeling physical systems, Neural Networks 146, 272 (2022).
  * Yousif  _et al._ [2023] M. Z. Yousif, M. Zhang, L. Yu, R. Vinuesa, and H. Lim, A transformer-based synthetic-inflow generator for spatially developing turbulent boundary layers, Journal of Fluid Mechanics 957, A6 (2023).
  * Borrelli  _et al._ [2022] G. Borrelli, L. Guastoni, H. Eivazi, P. Schlatter, and R. Vinuesa, Predicting the temporal dynamics of turbulent channels through deep learning, International Journal of Heat and Fluid Flow 96, 109010 (2022).
  * Zhou and Alam [2016] Y. Zhou and M. M. Alam, Wake of two interacting circular cylinders: A review, International Journal of Heat and Fluid Flow 62, 510 (2016).
  * Deng  _et al._ [2020] N. Deng, B. R. Noack, M. Morzyński, and L. R. Pastur, Low-order model for successive bifurcations of the fluidic pinball, Journal of fluid mechanics 884, A37 (2020).
  * Ren  _et al._ [2021] C. Ren, L. Cheng, C. Xiong, F. Tong, and T. Chen, Bistabilities in two parallel kármán wakes, Journal of Fluid Mechanics 929, A5 (2021).
  * Khodkar and Hassanzadeh [2021] M. Khodkar and P. Hassanzadeh, A data-driven, physics-informed framework for forecasting the spatiotemporal evolution of chaotic dynamics with nonlinearities modeled as exogenous forcings, Journal of Computational Physics 440, 110412 (2021).
  * Hochreiter and Schmidhuber [1997] S. Hochreiter and J. Schmidhuber, Long Short-Term Memory, Neural Computation 9, 1735 (1997), https://direct.mit.edu/neco/article-pdf/9/8/1735/813796/neco.1997.9.8.1735.pdf .
  * Kingma and Welling [2013] D. P. Kingma and M. Welling, Auto-encoding variational bayes, arXiv preprint arXiv:1312.6114 (2013).
  * Kullback and Leibler [1951] S. Kullback and R. A. Leibler, On information and sufficiency, The Annals of Mathematical Statistics 22, 79 (1951).
  * Higgins  _et al._ [2017b] I. Higgins, L. Matthey, A. Pal, C. Burgess, X. Glorot, M. Botvinick, S. Mohamed, and A. Lerchner, beta-VAE: Learning basic visual concepts with a constrained variational framework, in  _International Conference on Learning Representations_ (2017).
  * Clevert  _et al._ [2015] D.-A. Clevert, T. Unterthiner, and S. Hochreiter, Fast and accurate deep network learning by exponential linear units (elus), arXiv preprint arXiv:1511.07289 (2015).
  * Sirovich [1987] L. Sirovich, Turbulence and the dynamics of coherent structures. i. coherent structures, Quarterly of applied mathematics 45, 561 (1987).
  * Kazemi  _et al._ [2019] S. M. Kazemi, R. Goel, S. Eghbali, J. Ramanan, J. Sahota, S. Thakur, S. Wu, C. Smyth, P. Poupart, and M. Brubaker, Time2vec: Learning a vector representation of time, arXiv preprint arXiv:1907.05321 (2019).
  * Bahdanau  _et al._ [2014] D. Bahdanau, K. Cho, and Y. Bengio, Neural machine translation by jointly learning to align and translate, arXiv preprint arXiv:1409.0473 (2014).
  * He  _et al._ [2016] K. He, X. Zhang, S. Ren, and J. Sun, Deep residual learning for image recognition, in  _Proceedings of the IEEE conference on computer vision and pattern recognition_ (2016) pp. 770–778.
  * O’Malley  _et al._ [2019] T. O’Malley, E. Bursztein, J. Long, F. Chollet, H. Jin, L. Invernizzi, _et al._ , Kerastuner, https://github.com/keras-team/keras-tuner (2019).
  * Abadi  _et al._ [2015] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, J. Shlens, B. Steiner, I. Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng, TensorFlow: Large-scale machine learning on heterogeneous systems (2015), software available from tensorflow.org.
  * Kingma and Ba [2014] D. Kingma and J. Ba, Adam: A method for stochastic optimization, International Conference on Learning Representations (2014).
  * Zhang  _et al._ [2019] M. Zhang, J. Lucas, J. Ba, and G. E. Hinton, Lookahead optimizer: k steps forward, 1 step back, Advances in neural information processing systems 32 (2019).
  * Pascanu  _et al._ [2013] R. Pascanu, T. Mikolov, and Y. Bengio, On the difficulty of training recurrent neural networks, in  _International conference on machine learning_ (Pmlr, 2013) pp. 1310–1318.
  * Ginsburg  _et al._ [2019] B. Ginsburg, P. Castonguay, O. Hrinchuk, O. Kuchaiev, V. Lavrukhin, R. Leary, J. Li, H. Nguyen, Y. Zhang, and J. M. Cohen, Stochastic gradient methods with layer-wise adaptive moments for training of deep networks, arXiv preprint arXiv:1905.11286 (2019).
  * Taira and Colonius [2007] K. Taira and T. Colonius, The immersed boundary method: a projection approach, Journal of Computational Physics 225, 2118 (2007).
  * Colonius and Taira [2008] T. Colonius and K. Taira, A fast immersed boundary method using a nullspace approach and multi-domain far-field boundary conditions, Computer Methods in Applied Mechanics and Engineering 197, 2131 (2008).
  * Asztalos  _et al._ [2023] K. J. Asztalos, A. Almashjary, and S. Dawson, Galerkin spectral estimation of vortex-dominated wake flows, arXiv preprint arXiv:2302.06412 (2023).


# A Variable Autonomy approach for an Automated Weeding Platform

Ionut Moraru<EMAIL_ADDRESS>0000-0002-8314-037X School of Computer
Science, University of LincolnLincolnUK ,  Tsvetan Zhivkov
<EMAIL_ADDRESS>Lincoln Institute for Agri-food Technology, University
of LincolnLincolnUK ,  Shaun Coutts<EMAIL_ADDRESS>Lincoln Institute
for Agri-food Technology, University of LincolnLincolnUK ,  Dom Li
<EMAIL_ADDRESS>Lincoln Institute for Agri-food Technology, University of
LincolnLincolnUK and  Elizabeth I Sklar<EMAIL_ADDRESS>Lincoln
Institute for Agri-food Technology, University of LincolnLincolnUK

(2023)

###### Abstract.

Climate change, increase in world population and the war in Ukraine have led
nations such as the UK to put a larger focus on food security, while
simultaneously trying to halt declines in biodiversity and reduce risks to
human health posed by chemically-reliant farming practices. Achieving these
goals simultaneously will require novel approaches and accelerating the
deployment of Agri-Robotics from the lab and into the field. In this paper we
describe the ARWAC robot platform for mechanical weeding. We explain why the
mechanical weeding approach is beneficial compared to the use of pesticides
for removing weeds from crop fields. Thereafter, we present the system design
and processing pipeline for generating a course of action for the robot to
follow, such that it removes as many weeds as possible. Finally, we end by
proposing a trust-based ladder of autonomy that will be used, based on the
users’ confidence in the robot system.

Robotics, Human-Robotic-Interaction, Autonomous Planning

††copyright: acmcopyright††journalyear: 2023††conference: Variable Autonomy
for human-robot Teaming (VAT) workshop at the 18th ACM/IEEE International
Conference on Human-Robot Interaction (HRI); March 13, 2023; Stockholm,
Sweden††booktitle: Variable Autonomy for human-robot Teaming (VAT) workshop at
the 18th ACM/IEEE International Conference on Human-Robot Interaction (HRI),
March 13, 2023, Stockholm, Sweden

Figure 1. ARWAC Weeding Platform at Riseholme Campus, University of Lincoln.
October 2022.

ARWAC Weeding Platform at Riseholme Campus, University of Lincoln.

##  1\. Introduction

Figure 2. System design diagram of the ARWAC platform.

System design with regards to data flow.

As the fields of AI and Robotics mature, their impact will start being felt
outside of well-structured environments that have little to no outside factors
modifying them. We take inspiration from domains such as logistics and
warehouse management (Bogue, 2016; Mahroof, 2019), and look at how this
evolution in robotics will change agriculture.

The need to feed the world’s growing population as the effects of climate
change accelerate, while simultaneously reducing the environmental impacts of
agriculture has driven a long term focus on agri-robotics as a technology that
could potentially deliver these twin goals (Duckett et al., 2018). This
ambition has evolved from a distant goal to a necessity in recent years, as
nations feel the effects of global supply chain disruptions from the pandemic
and the war in Ukraine. The subjects of interest in agri-robotics have been
focusing on either picking soft fruits and vegetables (Huang et al., 2020;
Navas et al., 2021), precision herbicide spraying (Salazar-Gomez et al.,
2022), crop monitoring and optimisation (Harman et al., 2022), or advancing
auxiliary technology for aiding robotics in remote fields (Zhivkov and Sklar,
2022; Zhivkov et al., 2021; Moraru et al., 2021).

Our work focuses on non-chemical weed control. The use of critical chemical
tools such as herbicides will continue to come under pressure from evolved
weed resistance and tightening regulations to protect human and environmental
health (Hicks et al., 2018). In order to maintain global food security, both
more efficient use of herbicide and non-chemical weed and pest control will be
needed before we lose these valuable tools. Robotics and autonomy will be
critical technologies to enable this change to occur.

The rest of the paper is comprised of three sections. Section 2 describes the
ARWAC weeding robot platform, as shown in Figure 1. Section 3 proposes three
levels of autonomy that aim at increasing the users’ trust in the robot
platform. Finally, Section 4 ends with a presentation of related works from
the field of Agri-Robotics, focusing on systems designed for weeding and pest
removal.

##  2\. ARWAC Weeding Robot Platform

The ARWAC111ARWAC website - http://arwac.uk/ platform is a rover-type robot
built specifically with the intention of creating an autonomous and carbon-
neutral approach to identifying and removing blackgrass from fields without
harming the cereal-type crops. The robot is bespoke and unique in its approach
to tackle the challenge of weeding in agriculture. The common approach to
weeding is to improve on simpler and lower probability of failure methods,
such as spot-spraying. However, mechanical weeding is more complex as it
requires not only vision systems to identify weeds, but also precision tools
to target the weeds. Moreover, errors in tool displacement can cause far more
damage than precision spot-spraying, therefore requiring accurate crop row
detection and navigation systems. The ARWAC robot uses a mechanical weeding
tool, which removes the need for aerosols and chemicals that are bad for the
environment and human health, and can damage and kill crops.

The novelty behind the ARWAC robot is that it was designed to be relatively
light for agricultural tasks, with a weight under 250 kg, as this will reduce
soil compaction, increase battery life and can lead to the robot being able to
do more passes in a season under a wider range of ground conditions,
increasing the effectiveness of mechanical weeding. Moreover, it is a low cost
platform utilising a low power micro-processor for driving the motors and
Ackermann steering, and an onboard computer to perform sensor fusion and
autonomous planning.

###  2.1. System Design

The system design of the robot consists of a central computer that receives
data from an ESP32222ESP32 - https://docs.espressif.com/projects/esp-
idf/en/latest/esp32/ micro-processor, Ouster OS1333Ouster OS1 -
https://data.ouster.io/downloads/hardware-user-manual/hardware-user-manual-
revd-os1.pdf LiDAR, RGB camera and Ardusimple RTK-GPS (sensors), this is
illustrated in Figure 2. It uses the data from the micro-processor to
determine if the motors are functional and if the robot is moving, whereas
inputs from the LiDAR, RGB and RTK-GPS are used to determine a safe route to
navigate in the environment and identify weeds in the field. The central
computer uses data from the sensors to send control commands back to the
micro-processor to stop, move or avoid collisions. Additionally, the sensor
data is also used to send control commands to the micro-processor to control
the weeding tool. It should be noted that the central computer is using ROS1
(Robot Operating System (Quigley et al., 2009)) for communication and control
of sensors and robot motion.

Figure 3. Three part pipeline – mapping, processing and planning – for
generating a sequence of actions that the ARWAC Platform uses for mechanically
removing the blackgrass from the field of crops.

The control code is abstracted from the planner, allowing the planner to be
used as the main interface to control the robot’s behaviour, i.e., set,
remove, stop and edit robot tasks. The user (farmer) will be able to engage
with the robot via a simplified mobile app or GUI, made to directly
communicate with the planner. Rosbridge444rosbridge_suite package -
http://wiki.ros.org/rosbridge_suite is a ROS1 package that provides an easy to
use interface for any application, specifically web-based, to interact and
control ROS-based applications. This is what we will use to allow the user to
interact with the robot and configure the desired level of autonomy, with
basic descriptions of each level.

###  2.2. System Processing Pipeline

The project’s proposed processing pipeline can be seen in Figure 3 and is
split in three parts: mapping, pre-processing and planning. The first one
consists of a drone being flown over the field of crops, taking pictures which
will then be stitched together as a map. This will then be uploaded to a
central computer which will generate a probability weed map, which will inform
the system of how likely a spot is to have blackgrass growing in it. This step
could also be done with the input of farmers with a suitable interface to
allow them to input their knowledge of weed distribution on their fields to a
digital map.

The final step consists of synthesising a plan based on the areas identified
as most likely to have blackgrass growing around the arable crop. This part is
vital, as it will offer the agent high-level reasoning process. This is based
on PDDL planning (Fox and Long, 2003; Moraru et al., 2019), where solvers
create a plan based on the actions the agent poses to modify the agent such
that at the end, the environment will look as desired by the user. For this
implementation, there will be the possibility of adapting the plan based on
the users input, where we have three approaches based on users trust in the
systems reasoning process.

##  3\. Trust-based Autonomy Levels

Safe and trusted/explainable AI is one of the fastest growing sub-fields in
computer science research, due to many more AI approaches being incorporated
into systems currently being in production. In Ref. (Jacovi et al., 2021), the
notion of trust in AI is formalised from many different perspectives, and we
choose to work with the view of ”trust in model correctness”, which they
describe as being not how well the approach works, but how easy it is for the
user to distinguish where the AI succeeds or fails.

In a perfect world, where the user trusts the robot platform, it would be able
to act completely independent from the user and carry out its plan without any
need of validation or verification from outside. However, agri-robotic
platforms are still far from perfect and cause considerable economical and
physical harm to the crops or the user, with the ARWAC robot having the same
constraints. For such a system to be broadly adopted by the agricultural
sector, having the users’ trust will be vital to be earned, as shown in (Bedué
and Fritzsche, 2021; Canal et al., 2020).

For users to better understand and learn to trust the capabilities of the
ARWAC weeding platform, we have encorporated the robot with three levels of
autonomy for high-level reasoning. Their intended aim is so they would lead to
gaining the users trust in the platforms’s ability to complete its mission in
a safe and efficient manner, while also being able to verify that there are no
errors in the processing pipeline for their specific farm.

The three levels of autonomy are:

  * •

Low-Trust Autonomy – the user will come up with a plan for the platform to
follow, action by action, ;

  * •

Partial-Trust Autonomy – the robot comes up with a plan which it proposes to
the user (Canal et al., 2022), who will then challenge it until they are
satisfied with it, or reject and come up with a new one (Krarup et al., 2021);

  * •

Full-Trust Autonomy – it synthesises a plan which it will follow without the
need for human validation.

These three levels would allow time for the users to fully understand the
capabilities and limitations of the ARWAC weeding platform. It will do this by
first by first operating in a fully manual mode entitled ”Low-Trust autonomy”,
made with the intention that the user would pay closer attention to the
platform and which would lead to the user learning how the platform would
operate specifically on their farm. This mode is intended to be entered
initially, and used until the user is familiar with the robot and has used the
platform on a representative number of their fields.

Following this mode, the second level of autonomy would be the ”Partial-Trust
autonomy”, is designed for familiarising the user with platforms planning
approach to weeding a farm. This would be done by having the robot propose a
plan of action to the user, who would then be able to challenge it by asking
for contrastive explanations (Krarup et al., 2021).

Finally, once the user is accustomed to how the platform would tackle weeding,
from their individual actions to their own plans, the system has a ”Full-Trust
Autonomy” mode, which would allow for the ARWAC weeding platform to fully
autonomously complete a user selected field without any human intervention.

Our current schedule is to evaluate the capabilities of the ARWAC robot during
the spring season (late February to early May 2023), leading into studies of
users’ direct interaction with it. Once the platform is fully validated, we
will start evaluating our proposed autonomy levels such that we ensure obvious
bugs will not play the main role in how users will build trust in the system.

##  4\. Related Works

Outdoor robotic activities range from agricultural555Agri-robotics -
https://sagarobotics.com/, https://www.smallrobotcompany.com/,
delivery666Robot Delivery - https://www.starship.xyz/ to autonomous
shuttling777Robot Shuttle - https://www.navya.tech/en. Generally, outdoor
environments in this field are commonly referred to as unstructured or semi
structured. There are two main reasons for this. Firstly, in outdoor
environments there are few or no permanent obstacles/walls that can be used to
memorise (build) a map, which can be used to localise a robot with. Secondly,
communication and control of dynamic events cannot be guaranteed. For example,
in an indoor environment, navigating around people is the dynamic event, and
communication can be guaranteed to a certain extent, with predicted behaviour
and fallback methods. Comparatively, in an outdoor environment, dynamic events
are unbounded (weather, animals, vehicles, etc.,) and communication cannot be
guaranteed. This makes outdoor navigation and localisation more complex than
it is for indoor, and GPS data is a necessity (Ohno et al., 2004). However,
there are multiple methods that can be used to improve sensor fusion and
thereby improve outdoor localisation (O’Mahony et al., 2019); e.g.
recognition/horizon detection.

Agri-robotic weeding platforms exist in two main forms, herbicide sprayers and
mechanical weeders. Both forms fall under the umbrella term of precision
agriculture (PA). Moreover, both forms generally make use of machine vision
systems to identify the exact location of weeds (Li et al., 2022; Salazar-
Gomez et al., 2021, 2022), which allow the sprayer or weeding tool to carry
out its function. Mechanical weeders require combinations of weeding and
cultivation strategies to achieve acceptable weed control, compared to their
sprayer counterpart (Machleb et al., 2020). There are multiple tool mechanism
that are used by mechanical weeders such as harrow (Bakker et al., 2010) and
rotating (Nørremark et al., 2008) tools. The ARWAC robot makes use of a
harrow-like mechanism with blades that cut the identified weeds.

###### Acknowledgements.

We would like to express our gratitude toward ARWAC for providing their robot
platform and their continued help toward the project. Also to the farmers
whose continued help and input has been and will be critical to the success of
this project. This work is supported by Innovate UK (UKRI) under grant no.
#10027958 and project titled “ARWAC attack Blackgrass in Farming”.

## References

  * (1)
  * Bakker et al. (2010) Tijmen Bakker, Kees Asselt, Jan Bontsema, Joachim Müller, and Gerrit Straten. 2010.  Systematic design of an autonomous platform for robotic weeding.  _Journal of Terramechanics_ 47, 2 (2010), 63–73.  https://doi.org/10.1016/j.jterra.2009.06.002
  * Bedué and Fritzsche (2021) Patrick Bedué and Albrecht Fritzsche. 2021.  Can we trust AI? An empirical investigation of trust requirements and guide to successful AI adoption.  _Journal of Enterprise Information Management_ (2021).
  * Bogue (2016) Robert Bogue. 2016\.  Growth in e-commerce boosts innovation in the warehouse robot market.  _Industrial Robot: An International Journal_ 43, 6 (2016), 583–587.  https://doi.org/10.1108/IR-07-2016-0194
  * Canal et al. (2020) Gerard Canal, Rita Borgo, Andrew Coles, Archie Drake, Dong Huynh, Perry Keller, Senka Krivić, Paul Luff, Quratul-ain Mahesar, Luc Moreau, et al. 2020\.  Building trust in human-machine partnerships.  _Computer Law & Security Review_ 39 (2020), 105489.
  * Canal et al. (2022) Gerard Canal, Senka Krivić, Paul Luff, and Andrew Coles. 2022.  PlanVerb: Domain-Independent Verbalization and Summary of Task Plans. In _Proceedings of the AAAI Conference on Artificial Intelligence_ , Vol. 36. 9698–9706.
  * Duckett et al. (2018) Tom Duckett, Simon Pearson, Simon Blackmore, Bruce Grieve, Wen-Hua Chen, Grzegorz Cielniak, Jason Cleaversmith, Jian Dai, Steve Davis, Charles Fox, et al. 2018\.  Agricultural robotics: the future of robotic agriculture.  _arXiv preprint arXiv:1806.06762_ (2018).
  * Fox and Long (2003) Maria Fox and Derek Long. 2003.  PDDL2. 1: An extension to PDDL for expressing temporal planning domains.  _Journal of artificial intelligence research_ 20 (2003), 61–124.
  * Harman et al. (2022) Helen Harman, Elizabeth Sklar, et al. 2022\.  Multi-agent task allocation for harvest management.  _Frontiers in Robotics and AI_ (2022).
  * Hicks et al. (2018) Helen L Hicks, David Comont, Shaun R Coutts, Laura Crook, Richard Hull, Ken Norris, Paul Neve, Dylan Z Childs, and Robert P Freckleton. 2018.  The factors driving evolved herbicide resistance at a national scale.  _Nature ecology & evolution_ 2, 3 (2018), 529–536.
  * Huang et al. (2020) Zhuoling Huang, Elizabeth Sklar, and Simon Parsons. 2020\.  Design of Automatic Strawberry Harvest Robot Suitable in Complex Environments _(HRI ’20)_. Association for Computing Machinery, New York, NY, USA, 567–569.  https://doi.org/10.1145/3371382.3377443
  * Jacovi et al. (2021) Alon Jacovi, Ana Marasović, Tim Miller, and Yoav Goldberg. 2021.  Formalizing trust in artificial intelligence: Prerequisites, causes and goals of human trust in AI. In _Proceedings of the 2021 ACM conference on fairness, accountability, and transparency_. 624–635.
  * Krarup et al. (2021) Benjamin Krarup, Senka Krivic, Daniele Magazzeni, Derek Long, Michael Cashmore, and David E Smith. 2021\.  Contrastive explanations of plans through model restrictions.  _Journal of Artificial Intelligence Research_ 72 (2021), 533–612.
  * Li et al. (2022) Yong Li, Zhiqiang Guo, Feng Shuang, Man Zhang, and Xiuhua Li. 2022.  Key technologies of machine vision for weeding robots: A review and benchmark.  _Computers and Electronics in Agriculture_ 196 (2022), 106880.  https://doi.org/10.1016/j.compag.2022.106880
  * Machleb et al. (2020) Jannis Machleb, Gerassimos G. Peteinatos, Benjamin L. Kollenda, Dionisio Andújar, and Roland Gerhards. 2020.  Sensor-based mechanical weed control: Present state and prospects.  _Computers and Electronics in Agriculture_ 176 (2020), 105638.  https://doi.org/10.1016/j.compag.2020.105638
  * Mahroof (2019) Kamran Mahroof. 2019\.  A human-centric perspective exploring the readiness towards smart warehousing: The case of a large retail distribution warehouse.  _International Journal of Information Management_ 45 (2019), 176–190.  https://doi.org/10.1016/j.ijinfomgt.2018.11.008
  * Moraru et al. (2021) Ionut Moraru, Gerard Canal, and Simon Parsons. 2021\.  Using Plan Libraries for Improved Plan Execution. In _UKRAS21 Conference:“Robotics at home” Proceedings_.
  * Moraru et al. (2019) Ionut Moraru, Stefan Edelkamp, Santiago Franco, and Moises Martinez. 2019.  Simplifying automated pattern selection for planning with symbolic pattern databases. In _Joint German/Austrian Conference on Artificial Intelligence (Künstliche Intelligenz)_. Springer, 249–263.
  * Navas et al. (2021) Eduardo Navas, Roemi Fernández, Delia Sepúlveda, Manuel Armada, and Pablo Gonzalez-de Santos. 2021.  Soft grippers for automatic crop harvesting: A review.  _Sensors_ 21, 8 (2021), 2689.
  * Nørremark et al. (2008) M. Nørremark, H.W. Griepentrog, J. Nielsen, and H.T. Søgaard. 2008.  The development and assessment of the accuracy of an autonomous GPS-based system for intra-row mechanical weed control in row crops.  _Biosystems Engineering_ 101, 4 (2008), 396–410.  https://doi.org/10.1016/j.biosystemseng.2008.09.007
  * Ohno et al. (2004) Kazunori Ohno, Takashi Tsubouchi, Bunji Shigematsu, and Shin’ichi Yuta. 2004.  Differential GPS and odometry-based outdoor navigation of a mobile robot.  _Advanced Robotics_ 18, 6 (2004), 611–635.
  * O’Mahony et al. (2019) Niall O’Mahony, Sean Campbell, Anderson Carvalho, Suman Harapanahalli, Gustavo Adolfo Velasco-Hernández, Daniel Riordan, and Joseph Walsh. 2019\.  Adaptive Multimodal Localisation Techniques for Mobile Robots in Unstructured Environments : A Review. In _IEEE 5th World Forum on Internet of Things (WF-IoT)_. 799–804.  https://doi.org/10.1109/WF-IoT.2019.8767330
  * Quigley et al. (2009) Morgan Quigley, Ken Conley, Brian Gerkey, Josh Faust, Tully Foote, Jeremy Leibs, Rob Wheeler, Andrew Y Ng, et al. 2009\.  ROS: an open-source Robot Operating System. In _ICRA workshop on open source software_ , Vol. 3. Kobe, Japan, 5.
  * Salazar-Gomez et al. (2021) Adrian Salazar-Gomez, Madeleine Darbyshire, Junfeng Gao, Elizabeth I Sklar, and Simon Parsons. 2021\.  Towards practical object detection for weed spraying in precision agriculture.  https://doi.org/10.48550/ARXIV.2109.11048
  * Salazar-Gomez et al. (2022) Adrian Salazar-Gomez, Madeleine Darbyshire, Junfeng Gao, Elizabeth I Sklar, and Simon Parsons. 2022\.  Beyond mAP: Towards practical object detection for weed spraying in precision agriculture. In _IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 9232–9238.
  * Zhivkov et al. (2021) Tsvetan Zhivkov, Adrian Gomez, Junfeng Gao, Elizabeth Sklar, and Simon Parsons. 2021.  The need for speed: How 5G communication can support AI in the field. In _Proceedings of UKRAS21 Conference: Robotics at home_. UK-RAS.  https://doi.org/10.31256/On8Hj9U
  * Zhivkov and Sklar (2022) Tsvetan Zhivkov and Elizabeth I Sklar. 2022.  5G on the Farm: Evaluating Wireless Network Capabilities for Agricultural Robotics.  _arXiv preprint arXiv:2301.01600_ (2022).