ilesRGB.icc

# QMKPy: A Python Testbed for the Quadratic Multiple Knapsack Problem

Karl-Ludwig Besser  0000-0002-1597-8963 Eduard A. Jorswieck
0000-0001-7893-8435 0000-0002-9402-1982

††margin:  DOI: 10.21105/joss.04718 Software • Review • Repository • Archive
Editor: Mehmet Hakan Satman
Reviewers: • @ulf1 • @max-little Submitted: 08 August 2022
Published: 02 November 2022 License
Authors of papers retain copyright and release the work under a Creative
Commons Attribution 4.0 International License (CC BY 4.0).

## Summary

QMKPy provides a Python framework for modeling and solving the quadratic
multiple knapsack problem (QMKP). It is primarily aimed at researchers who
develop new solution algorithms for the QMKP. QMKPy therefore mostly functions
as a testbed to quickly implement novel algorithms and compare their results
with existing ones. However, the package also already includes implementations
of established algorithms for those who only need to solve a QMKP as part of
their research.

The QMKP is a type of knapsack problem which has first been analyzed by Hiley
& Julstrom (2006). For a basic overview of other types of knapsack problems
see, e.g., Kellerer et al. (2004). As in the regular multiple knapsack
problem, the goal is to assign $N\in\mathbb{N}$ items with given weights
$w_{i}\in\mathbb{R}_{+}$ and (non-negative) profit values
$p_{i}\in\mathbb{R}_{+}$ to $K\in\mathbb{N}$ knapsacks with given weight
capacities $c_{u}\in\mathbb{R}_{+}$, such that a total profit is maximized. In
the QMKP, there exist additional joint profits $p_{ij}\in\mathbb{R}_{+}$ which
are yielded when items $i$ and $j$ are packed into the same knapsack.

Mathematically, the QMKP is described by the following optimization problem

$\displaystyle\max\quad$
$\displaystyle\sum_{u\in\mathcal{K}}\Bigg{(}\sum_{i\in\mathcal{A}_{u}}p_{i}$
$\displaystyle+$ $\displaystyle\sum_{\begin{subarray}{c}j\in\mathcal{A}_{u}\\\
j\neq i\end{subarray}}p_{ij}\Bigg{)}$ (1a) $\displaystyle\mathrm{s.t.}\quad$
$\displaystyle\sum_{i\in\mathcal{A}_{u}}w_{i}\leq c_{u}$ $\displaystyle\forall
u\in\mathcal{K}$ (1b) $\displaystyle\sum_{u=1}^{K}a_{iu}\leq 1$
$\displaystyle\forall i\in\\{1,2,\dots,N\\}$ (1c)

where $\mathcal{K}=\\{1,2,\dots,K\\}$ describes the set of knapsacks,
$\mathcal{A}_{u}\subseteq\\{1,2,\dots,N\\}$ is the set of items that are
assigned to knapsack $u$ and $a_{iu}\in\\{0,1\\}$ is a binary variable
indicating whether item $i$ is assigned to knapsack $u$.

## Statement of need

The QMKP is an NP-hard optimization problem and therefore, there exists a
variety of (heuristic) algorithms to find good solutions for it. While Python
frameworks already exist for the standard (multiple) knapsack problem and the
quadratic knapsack problem, they do not consider the _quadratic multiple_
knapsack problem. However, this type of problem arises in many areas of
research. In addition to the typical problems in Operations Research, it also
occurs in distributed computing (Rust et al., 2020) and in the area of
wireless communications (Besser et al., 2022).

For the classic knapsack problem and the quadratic _single_ knapsack problem,
many well-known optimization frameworks like Gurobi (Gurobi Optimization, LLC,
2022) and OR-Tools (Perron & Furnon, 2022) provide solution routines. However,
they are not directly applicable to the QMKP. Furthermore, it can be difficult
for researchers to reproduce results that rely on commerical software.

Therefore, QMKPy aims to close that gap by providing an open source testbed to
easily implement and compare solution algorithms. Additionally, Python is
widely used among researchers and enables easy to read implementations. This
further supports the goal of QMKPy to promote sharable and reproducible
solution algorithms.

For initial comparisons, the software already implements multiple solution
algorithms for the QMKP, including a _constructive procedure (CP)_ based on
Algorithm 1 from Aïder et al. (2022) and the greedy heuristic from Hiley &
Julstrom (2006). A second algorithm that is included is the _fix and complete
solution (FCS) procedure_ from Algorithm 2 in Aïder et al. (2022).
Additionally, a collection of reference QMKP instances that can be used with
QMKPy is provided at (Besser, 2022). This dataset includes the well-known
reference problems from Hiley & Julstrom (2006), which in turn are based on
the quadratic single knapsack problems from Billionnet & Soutif (2004).

The open source nature of QMKPy and its aim at researchers encourages the
implementation of more algorithms for solving the QMKP that can become part of
the QMKPy framework.

The most notable benefits when implementing an algorithm using QMKPy are the
following:

  * •

No additional overhead is required. Only a single function with the novel
solution algorithm needs to be implemented.

  * •

Generic unit tests are available to make sure that the novel algorithm
fulfills the set of basic criteria.

  * •

The ability of loading and saving problem instances allows for quick and easy
testing of any algorithm against reference datasets. This enables reproducible
research and creates a high degree of comparability between different
algorithms.

## References

reAïder, M., Gacem, O., & Hifi, M. (2022). Branch and solve strategies-based
algorithm for the quadratic multiple knapsack problem. _Journal of the
Operational Research Society_ , _73_(3), 540–557.
https://doi.org/10.1080/01605682.2020.1843982

preBesser, K.-L. (2022). _QMKPy datasets_ (Version v1.0) [Computer software].
Zenodo. https://doi.org/10.5281/ZENODO.7157144

preBesser, K.-L., Jorswieck, E. A., & Coon, J. P. (2022, September). Multi-
user frequency assignment for ultra-reliable mmWave two-ray channels. _20th
International Symposium on Modeling and Optimization in Mobile, Ad Hoc, and
Wireless Networks (WiOpt)_.

preBillionnet, A., & Soutif, Éric. (2004). Using a mixed integer programming
tool for solving the 0–1 quadratic knapsack problem. _INFORMS Journal on
Computing_ , _16_(2), 188–197. https://doi.org/10.1287/ijoc.1030.0029

preGurobi Optimization, LLC. (2022). _Gurobi Optimizer Reference Manual_.
https://www.gurobi.com

preHiley, A., & Julstrom, B. A. (2006). The quadratic multiple knapsack
problem and three heuristic approaches to it. _Proceedings of the 8th Annual
Conference on Genetic and Evolutionary Computation - GECCO ’06_ , 547–552.
https://doi.org/10.1145/1143997.1144096

preKellerer, H., Pferschy, U., & Pisinger, D. (2004). _Knapsack problems_.
Springer Berlin Heidelberg. https://doi.org/10.1007/978-3-540-24777-7

prePerron, L., & Furnon, V. (2022). _OR-tools_ (Version v9.4) [Computer
software]. Google. https://developers.google.com/optimization/

preRust, P., Picard, G., & Ramparany, F. (2020). Resilient distributed
constraint optimization in physical multi-agent systems. _24th European
Conference on Artificial Intelligence (ECAI 2020)_ , 195–202.
https://doi.org/10.3233/FAIA200093

p


# Visual-Geometry GP-based Navigable Space for Autonomous Navigation

Mahmoud Ali, Durgkant Pushp, Zheng Chen, and Lantao Liu All authors are with
the Luddy School of Informatics, Computing, and Engineering, Indiana
University, Bloomington, IN 47408 USA. {alimaa, dpushp, zc11<EMAIL_ADDRESS>
###### Abstract

Autonomous navigation in unknown environments is challenging and demands the
consideration of both geometric and semantic information in order to parse the
navigability of the environment. In this work, we propose a novel space
modeling framework, Visual-Geometry Sparse Gaussian Process (VG-SGP), that
simultaneously considers semantics and geometry of the scene. Our proposed
approach can overcome the limitation of visual planners that fail to recognize
geometry associated with the semantic and the geometric planners that
completely overlook the semantic information which is very critical in real-
world navigation. The proposed method leverages dual Sparse Gaussian Processes
in an integrated manner; the first is trained to forecast geometrically
navigable spaces while the second predicts the semantically navigable areas.
This integrated model is able to pinpoint the overlapping (geometric and
semantic) navigable space. The simulation and real-world experiments
demonstrate that the ability of the proposed VG-SGP model, coupled with our
innovative navigation strategy, outperforms models solely reliant on visual or
geometric navigation algorithms, highlighting a superior adaptive behavior. We
provided a demonstration video111Video: https://youtu.be/0s6VSj5Z1dg and open-
sourced our code 222Code: https://github.com/mahmoud-a-ali/vg-nav.

##  I Introduction and Related Work

It is a challenging task for autonomous robots to navigate in real-world,
unpredictable field environments which typically contain considerable
complexity, with obstacles that may include large rocks, shrubs, and tree
stumps, as well as terrains of mixed composition such as asphalt, sand, and
mud. When deploying a robot in such a setting, it must determine the navigable
portions of the scene it captures. One approach to identifying navigable space
involves utilizing the geometry or structure of surrounding objects, without
considering their semantic interpretation. This geometric method typically
leverages a 3D representation of the environment, such as point clouds or
depth observations [1, 2]. Geometric methodologies commonly employ global or
local representation of the navigable space of the environment. Global
representations include the use of an occupancy grid map, which encodes the
likelihood of occupancy within a 2D spatial domain [3], a 2.5D elevation
map[4] that outlines the terrain elevation, or 3D maps that uses voxels to
model the 3D space, such as Octomaps [5]. On the other hand, local
representations such as cost and traversability maps are used to aid local
planners in devising navigational strategies [6]. Another approach to
identifying the navigable spaces is the semantic approach which performs
semantic segmentation of the scene to predict its traversability using visual
input [7, 8, 9]. Traversability may be assigned manually based on terrain
types, for example attributing low traversability to mud areas and high
traversability to grass areas [10, 11, 12]. In contrast, recent learning-based
approaches [13, 14, 15, 16, 17] circumvent heuristic methods by leveraging
data to learn traversability costs.

(a)

(b)

Figure 1: Visual-Geometry Navigable Space: (a) shows LiDAR’s pointcloud
(yellow) and camera’s colored-pointcloud; (b) shows the geometry navigable
space (grey-coded circular surface, white regions represent free space), the
visual navigable space (grey-coded vertical plane, white regions represents
the navigable classes in the camera field of view), and the local navigation
points (colored-circles, where color represents the go-to-goal cost).

Moreover, visual navigation has been approached through various methods like
utilizing motion information from targets or optical flow features [18, 19],
storing environment images with corresponding control actions [20], using
image signal entropy for non-holonomic robot control [21], feature tracking
[22, 23], and path-finding through image segmentation based on control
parameters from path boundaries [24]. However, visual navigation may fall
short in challenging environments where visually appealing navigable scenes
(e.g. high-slope grassland) are non-navigable geometrically. Gaussian Process
(GP) has been used for modeling continuous spatial phenomena and creating
navigability maps based on geometric [25] or semantic interpretations [26].
While standard GPs are limited in real-time applications due to their high
time complexity of $\mathcal{O}(n^{3})$. Sparse GP (SGP) methods reduce GP’s
complexity to $\mathcal{O}(nm^{2})$ by selecting $m$ inducing points to
approximate the full dataset under Bayesian rules [27, 28, 29].

LiDAR-based geometric navigation is more robust than visual-based approaches
however, they may incorrectly identify navigable spaces due to a lack of
semantic consideration. While visual methods offer flexibility through pixel
analysis but are limited by weather due to visual dependency. Recent
navigation strategies merge geometric and visual data, employing end-to-end
learning to analyze commands and trajectories across terrains [13, 15], but
struggle with 3D terrain complexity due to the high cost of data analysis.
Other methods use global maps to integrate semantic and geometric costs into a
combined navigable space [30, 10]. In contrast, we propose a local mapless
navigation approach that integrates visual and geometric navigable spaces
using two local SGP models.

##  II Methodology

We propose a new framework that combines the geometry and the semantics of the
robot’s surroundings to identify the navigable spaces around the robot.
Briefly, the output of the RGBD camera and LiDAR are processed in a parallel
way where the RGB image is segmented by labeling each pixel with a unique
class (grass, tree, asphalt, etc), and the pointcloud is represented by an
occupancy surface around the robot. Consequentially, the segmented image is
converted into a binary navigability image (navigable and non-navigable
pixels) based on a defined set of navigable classes. The navigability is then
projected on the camera’s depth pointcloud which is represented as a Visual
SGP model (V-SGP) to find the (visual) navigable spaces in front of the robot.
On the other hand, the occupancy surface is represented by a Geometry SGP
model (G-SGP) to assess the free and occupied spaces and to identify the
(geometric) navigable space around the robot. The visual-based and the
geometry-based navigable spaces are coupled to calculate a more accurate
navigable space, based on which the Local Navigation Points (LNPs) are
generated to drive the robot to its destination, see Fig. 1.

(a)

(b)

(c)

(d)

Figure 2: Geometry-Navigable Space: (b) shows the raw pointcloud in yellow and
the original occupancy surface , where warmer colors indicate less occupancy;
(c) shows the predicted occupancy surface using the G-SGP model; (d) shows the
variance surface, i.e., the uncertainty associated with the predicted
occupancy, where white color indicates highly uncertain (free) points. The
Geometry-LNPs (G-LNPs) are shown as colored circles, where green indicates
less go-to-goal cost associated with each G-LNP.

###  II-A Geometry Navigable Space: The G-SGP Model

The LiDAR’s pointcloud,
$\mathbf{P}_{l}=\\{(x_{i},y_{i},z_{i})\\}_{i=1}^{n_{l}}$, is converted to the
occupancy surface, $\mathcal{S}_{g}$, representation [31], where each
(geometry) point, $\mathbf{p_{g}}_{i}$, on $\mathcal{S}_{g}$ is defined by its
azimuth, $\alpha_{i}$, and elevation, $\beta_{i}$, angles, and given an
occupancy value, $\Omega_{i}$, equals to the difference between the surface
radius, $\rho_{g}$, and the point radius, $\rho_{i}$, as follows
$\Omega_{i}=\rho_{g}-\rho_{i}$. The surface regions with projected points
represent the occupied space of $\mathcal{S}_{g}$. In contrast, other regions
with no points represent the free space of $\mathcal{S}_{g}$, see Fig.
LABEL:fig_sgp_oc_b. The $n_{g}$ projected points on $\mathcal{S}_{g}$ form the
geometric data set
$\mathcal{D}_{g}=\left\\{\left(\mathbf{p_{g}}_{i},\Omega_{i}\right)\right\\}_{i=1}^{n_{g}}$,
where $\mathbf{p_{g}}_{i}=(\alpha_{i},\beta_{i})$, and $\Omega_{i}$ is the
occupancy of $\mathbf{p_{g}}_{i}$. Subsequently, the geometric data set
$\mathcal{D}_{g}$ is employed to train a 2D variational SGP (G-SGP), to model
the probability of occupancy, $f_{g}(\mathbf{p_{g}}_{i})$, over
$\mathcal{S}_{g}$ as follows:

$\begin{gathered}f_{g}(\mathbf{p_{g}})\sim{SGP_{g}}\left(m_{g}(\mathbf{p_{g}}),k_{g}\left(\mathbf{p_{g}},\mathbf{p_{g}}^{\prime}\right)\right),\\\
k_{\mathrm{g}}\left(\mathbf{p_{g}},\mathbf{p_{g}}^{\prime}\right)=\sigma_{1}^{2}\left(1+\frac{\left(\mathbf{p_{g}}-\mathbf{p_{g}}^{\prime}\right)^{2}}{2\alpha_{1}\ell_{1}^{2}}\right)^{-\alpha_{1}},\end{gathered}$
(1)

where $m_{g}(\mathbf{p_{g}})$ is the zero mean function, and
$k_{g}\left(\mathbf{p_{g}},\mathbf{p_{g}}^{\prime}\right)$ is a Rational
Quadratics (RQ) kernel with a length-scale $\ell_{1}$, a signal variance
$\sigma_{1}^{2}$, and a relative weighting factor $\alpha_{1}$. A Gaussian
noise $\epsilon_{g}$ is added to the predicted occupancy to reflect the
measurement noise. The probability of occupancy $\Omega_{g}^{*}$ for any query
point $\mathbf{p_{g}}^{*}$ on $\mathcal{S}_{g}$ is calculated by the GP
prediction as follows,

$\begin{gathered}p_{g}(\Omega_{g}^{*}|\bm{\Omega}_{g})=\mathcal{N}_{g}(\Omega_{g}^{*}|m_{\bm{\Omega_{g}}}(\bm{p_{g}}^{*}),k_{\bm{\Omega_{g}}}(\bm{p_{g}}^{*},\bm{p_{g}}^{*})+\sigma_{n_{g}}^{2}),\\\
m_{\bm{\Omega_{g}}}(\bm{p_{g}})=K_{\bm{p_{g}}n_{g}}\left(\sigma_{g_{n}}^{2}I+K_{n_{g}n_{g}}\right)^{-1}\bm{\Omega_{g}},\\\
k_{\bm{\Omega_{g}}}\left(\bm{p_{g}},\mathbf{p_{g}}^{\prime}\right)=k\left(\bm{p_{g}},\bm{p_{g}}^{\prime}\right)-K_{\mathbf{p_{g}}n_{g}}\left(\sigma_{n_{g}}^{2}I+K_{n_{g}n_{g}}\right)^{-1}K_{n_{g}\bm{p_{g}}^{\prime}},\end{gathered}$
(2)

where $m_{\bm{\Omega_{g}}}(\bm{p_{g}})$ and
$k_{\bm{\Omega_{g}}}\left(\bm{p_{g}},\bm{p_{g}}^{\prime}\right)$ are the
posterior mean and covariance functions [29], $K_{n_{g}n_{g}}$ is $n_{g}\times
n_{g}$ co-variance matrix of the inputs, $K_{p_{g}n_{g}}$ is
$n_{g}$-dimensional row vector of kernel function values between $\bm{p_{g}}$
and the inputs, and $K_{n_{g}p_{g}}=K_{p_{g}n_{g}}^{T}$. We leverage the
variational SGP approach [29] to estimate the kernel hyperparameters $\Theta$
and to select the inducing points $X_{m}$, more details about the
implementation of the G-SGP model can be found in our previous work [31]. Fig.
LABEL:fig_sgp_oc_c shows the predicted occupancy $\mu_{g}$ on predicted
occupancy surface $\mathcal{S}_{\mu_{g}}$, where the prediction uncertainty
$\sigma_{g}$ is shown as the variance surface $\mathcal{S}_{\sigma_{g}}$ in
Fig. LABEL:fig_sgp_oc_d. Regarding the accuracy of the SGP occupancy model,
the reconstructed pointcloud from a G-SGP model with 400 inducing points has
an average error of $12\text{\,}\mathrm{cm}$ [31].

The variance surface $\mathcal{S}_{\sigma_{g}}$ discriminates efficiently
between the free space (white regions with a variance higher than a threshold
$V_{g_{th}}$) and the occupied space (dark regions with a variance less than
$V_{g_{th}}$) around the robot [32, 33], and reflects the terrain elevation
(the boundary between free and occupied space) in the local observation [34].
Therefore, $\mathcal{S}_{\sigma_{g}}$ is used to define a set of geometrical-
feasible LNPs (G-LNPs) around the robot in the free space that are considered
navigable based on the robot’s maximum roll and pitch angles. G-LNPs are the
lowest free points on the variance surface whose elevation angles are bounded
by the safe elevations that the robot can climb, ($\beta_{min_{s}}$ and
$\beta_{max_{s}}$);
$\text{G-LNPs}={(\alpha_{i},\beta_{j^{*}})}\,|\,-\pi<\alpha_{i}<\pi\\}$; where
$\beta_{min_{s}}<\beta_{j}^{*}<\beta_{max_{s}}$.

Formally, a G-LNP is defined as
$\mathbf{g}_{lnp_{i}}=(\alpha_{i},\beta_{i},\rho_{i})$, where $\alpha_{i}$
defines the direction of $\mathbf{g}_{lnp_{i}}$ with respect to the robot
heading, $\beta_{i}$ is the elevation of $\mathbf{g}_{lnp_{i}}$ with respect
to the robot, and $\rho_{i}$ is the distance between $\mathbf{g}_{lnp_{i}}$
and the robot predicted by the G-SGP model;
$\hat{\rho_{i}}=\rho_{g}-\hat{\Omega_{i}}$, where
$\hat{\Omega_{i}}=\mu_{g_{i}}$ and
$(\mu_{g_{i}},\sigma_{g_{i}})=\mathcal{SGP}_{g}((\alpha_{i},\beta_{i}))$. The
Cartesian coordinates of $\mathbf{g}_{lnp_{i}}$ within the global world frame
$\mathcal{W}$ are derived as
$(x^{w}_{i},y^{w}_{i},z^{w}i)=\prescript{W}{}{\mathbf{T}}{R}\cdot(x^{R}_{i},y^{R}_{i},z^{R}i)$,
where $\prescript{W}{}{\mathbf{T}}{R}$ represents the robot’s localization.
The coordinates $(x^{R}_{i},y^{R}_{i},z^{R}i)$ denote the position of
$\mathbf{g}{lnp_{i}}$ in the robot frame $\mathcal{R}$, calculated from its
spherical coordinates $(\alpha_{i},\beta_{i},\rho_{i})$.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 3: Visual-Navigable Space: V-LNPs are shown as colored-circles in (h),
where the color shows the cost assigned to each V-LNP.

###  II-B Visual Navigable Space: V-SGP Model

The RGB image contains crucial information contributing to obstacle
identification, which may not be revealed through LiDAR sensing (geometry-
navigable space). The semantic image $I_{seg}$ and the associated class labels
$I_{cls}$ can be obtained by $I_{seg},I_{cls}=g(I_{rgb},\Theta)$ where $g(.)$
represents any existing image segmentation model with parameter $\Theta$, see
Fig. LABEL:fig_vis_mdl_c. In this paper we use the mask2former segmentation
model[9]. Utilizing our domain knowledge regarding the physical properties of
objects within the scene, we construct the navigability image, denoted as
$I_{nav}$, by assigning a value $0$ to all pixels whose class labels are
categorized as non-navigable and a value $255$ to all other pixels, see Fig.
LABEL:fig_vis_mdl_d. The navigability image provides a dynamic categorization
of navigable classes which is crucial in defining flexible visually navigable
spaces. Depending on the robot’s capabilities and the required behavior,
various terrains such as grass, snow, or dirt can be designated as navigable
or non-navigable, enhancing the robot’s adaptive behavior across different
scenarios.

The navigability image ${I}_{nav}$ is projected on the visual depth pointcloud
$\mathbf{P}_{c}=\\{(x_{i},y_{i},z_{i},rgb_{i})\\}_{i=1}^{n_{v}}$ by replacing
the $rgb_{i}$ value with the binary navigability value, $\iota_{i}$, where
each point is set to either navigable or non-navigable, resulting in the
navigability pointcloud
$\mathbf{P}_{v}=\\{(x_{i},y_{i},z_{i},\iota_{i})\\}_{i=1}^{n_{v}}$, see Fig.
LABEL:fig_vis_mdl_e. Similar to $\mathcal{S}_{g}$ in Sec. II-A, the
navigability points are projected on a curved visual surface $\mathcal{S}_{v}$
with a predefined radius $\rho_{v}$, where the visual surface span is aligned
with the camera field of view (FoV). Each point, $\mathbf{p_{v}}_{i}$, on
$\mathcal{S}_{v}=\\{(\alpha_{i},\beta_{i},\rho_{i},\iota_{i})\\}_{i=1}^{n_{v}}$,
is represented by its azimuth, elevation, radius, and navigability values, see
Fig. LABEL:fig_vis_mdl_e. $\mathcal{S}_{v}$ is then decomposed into two data
sets, the visual navigability dataset
$\mathcal{D}_{nav}=\\{\mathbf{p_{v}}_{i},{\iota}_{i}\\}_{i=1}^{n_{v}}$ and the
visual depth data set
$\mathcal{D}_{dpth}=\\{\mathbf{p_{d}}_{i},{\rho}_{i}\\}_{i=1}^{n_{v}}$, where
$\mathbf{p_{d}}_{i}=(\alpha_{i},\beta_{i})$. $\mathcal{D}_{dpth}$ is then
filtered to include only the navigability points whose $\rho$ less than a
constant $\rho_{d}$, i.e.,
$\mathcal{D}_{dpth}=\\{\mathbf{p_{d}}_{i},{\rho}_{i}\\}_{i=1}^{n_{v}}:\mathbf{p_{d}}_{i}=\mathbf{p_{v}}_{i}\mid\rho_{i}<\rho_{d}$.
Thereafter, the visual surface is modeled using two SGP models. The first one
is the SGP Depth (D-SGP) model , which is an SGP regression model trained on
$\mathcal{D}_{dpth}$ in a similar way as the G-SGP model to predict the
occupancy of $\mathbf{p_{d}}_{i}$ as
${\Omega_{d}}_{i}=f_{d}(\mathbf{p_{d}}_{i})+\epsilon_{d}$:

$\begin{gathered}f_{d}(\mathbf{p_{d}})\sim{SGP_{d}}\left(m_{d}(\mathbf{p_{d}}),k_{d}\left(\mathbf{p_{d}},\mathbf{p_{d}}^{\prime}\right)\right),\\\
k_{\mathrm{d}}\left(\mathbf{p_{d}},\mathbf{p_{d}}^{\prime}\right)=\sigma_{2}^{2}\left(1+\frac{\left(\mathbf{p_{d}}-\mathbf{p_{d}}^{\prime}\right)^{2}}{2\alpha_{2}\ell_{2}^{2}}\right)^{-\alpha_{2}},\end{gathered}$
(3)

where $m_{d}(\mathbf{p_{d}})$,
$k_{d}\left(\mathbf{p_{d}},\mathbf{p_{d}}^{\prime}\right)$, and $\epsilon_{d}$
are similar to Eq. (1) On the other hand, the second model is the SGP
Navigability (N-SGP) model, which is an SGP classification model trained on
$\mathcal{D}_{nav}$ to predict the navigability status $\iota$ of any point
$\mathbf{p_{v}}_{i}$. The N-SGP classification model uses the same variational
SGP framework as the G-SGP model where
${\Omega_{v}}_{i}=f_{v}(\mathbf{p_{v}}_{i})+\epsilon_{v}$, however, its output
probability is thresholded to decide whether $\mathbf{p_{v}}_{i}$ is navigable
or not, see Fig. LABEL:fig_vis_mdl_g,

$\begin{gathered}f_{n}(\mathbf{p_{v}})\sim{SGP_{v}}\left(m_{v}(\mathbf{p_{v}}),k_{v}\left(\mathbf{p_{v}},\mathbf{p_{v}}^{\prime}\right)\right),\\\
k_{\mathrm{v}}\left(\mathbf{p_{v}},\mathbf{p_{v}}^{\prime}\right)=\sigma_{3}^{2}\left(1+\frac{\left(\mathbf{p_{v}}-\mathbf{p_{v}}^{\prime}\right)^{2}}{2\alpha_{3}\ell_{3}^{2}}\right)^{-\alpha_{3}},\end{gathered}$
(4) $\begin{gathered}\iota_{\mathbf{v}_{i}}=\begin{cases}\text{255
"navigable"}&\hat{\Omega_{v_{i}}}>\Omega_{v_{th}}\\\ \text{0 "non-
navigable"}&\text{otherwise.}\end{cases}\end{gathered}$ (5)

where $m_{d}(\mathbf{p_{d}})$,
$k_{d}\left(\mathbf{p_{d}},\mathbf{p_{d}}^{\prime}\right)$ are similar to Eq.
(1), and $\Omega_{v_{th}}$ is a predefined threshold. Fig.
LABEL:fig_vis_dpth_mdl_a shows the pointcloud predicted by the D-SGP model,
$\hat{\mathbf{P}_{c}}=\\{(x_{i},y_{i},z_{i},\sigma_{d_{i}})\\}_{i=1}^{N_{c}}$,
where the grey-color intensity indicates the uncertainty $\sigma_{d}$ of the
predicted radius $\hat{\rho}$, i.e., darker points are more certain. Since the
model is trained only on the points whose radius values $\rho_{i}$ are less
than $\rho_{d}$, the D-SGP prediction in the range of $\rho_{d}$ is certain.
That is why the dark (certain) points in Fig. LABEL:fig_vis_dpth_mdl_a are
similar, in terms of 3D location, to the flat region (asphalt and grass) of
the camera’s point cloud $\mathbf{P}_{c}$ in Fig. LABEL:fig_vis_mdl_a, while
the white (uncertain) points have an inaccurate prediction. For more details
about the accuracy of the generated pointcloud, check [31]. It is worth
mentioning that the smoothness property of the GP denoises the pointclouds,
where individual points that do not belong to any point cluster are smoothed
out. Fig. LABEL:fig_vis_dpth_mdl_b shows $\hat{\mathbf{P}_{c}}$ with the color
representing the navigability $\hat{\iota}$ predicted by the N-SGP model;
$\hat{\mathbf{P}_{v}}=\\{(x_{i},y_{i},z_{i},\iota_{i})\\}_{i=1}^{N}$. The
prediction combination of D-SGP and N-SGP, $\hat{\mathbf{P}_{v}}$, is a
reconstruction of the navigability pointcloud $\mathbf{P}_{v}$ in Fig.
LABEL:fig_vis_mdl_e. Only the certain points whose $\hat{\rho}\leq\rho_{d}$
and $\sigma_{d}<V_{d_{th}}$, are considered for planning.

The combination of the D-SGP and the N-SGP models is considered as one entity,
called the visual SGP (V-SGP) model, see Fig. 3. The V-SGP model can predict
both the 3D location of the missing points in the camera’s FoV (points with a
‘nan’ return value) and their navigability status, check the difference
$\mathbf{P}_{v}$ in Fig. LABEL:fig_vis_mdl_e and $\hat{\mathbf{P}_{v}}$
predicted by the V-SGP model in Fig. LABEL:fig_vis_dpth_mdl_b. In the
geometry-based approach, G-LNPs are selected solely based on the geometry
information encoded in $\mathcal{S}_{g}$, Contrastingly, $\mathcal{S}_{v}$
introduces a more nuanced layer to the navigable space assessment by
incorporating the navigability status denoted as $\iota$. In the visual-based
approach, points falling within the safe elevation boundaries,
$\beta_{min_{s}}<\beta<\beta_{max_{s}}$, are considered as LNP candidates.
Then the navigability cost $\iota$ is estimated using the N-SGP model for all
candidates to form the Visual LNPs (V-LNPs), $\bm{\nu}_{lnp}$, see Fig.
LABEL:fig_vis_mdl_h.

(a)

(b)

Figure 4: D-SGP and D-SGP models of camera’s pointcloud: both (a) and (b)
shows the pointcloud generated from the D-SGP model, wherein the grey color
indicates the uncertainty in (a) and the predicted navigability class in (b).

###  II-C LNPs Selection For Goal-oriented Navigation

In general, to navigate the robot towards a given goal
$\mathbf{g}=(x_{f},y_{f})$ in $\mathcal{W}$, local mapless navigation
approaches use different criteria to select one LNP (i.e., an ideal G-LNP
($\mathbf{g}_{lnp}^{*}$) or V-LNP ($\mathbf{\nu}_{lnp}^{*}$)) from the
feasible LNPs. For example, the FGM [35] selects the ideal LNP based on the
area of the free space around it and its direction to the final goal
$\mathbf{g}$. While the admissible gap approach [36] selects the closest LNP
to $\mathbf{g}$. We utilize the cost function
$C_{g}\left(\mathbf{g}_{lnp_{i}}\right)$ introduced in [34], which accounts
for the distance to $\mathbf{g}$, the alignment relative to the robot’s
heading, and the elevation angle of the LNP,

$\displaystyle C_{g}\left(\mathbf{g}_{lnp_{i}}\right)$
$\displaystyle=k_{dst}d_{tg}+k_{dir}\|\alpha_{i}\|+k_{elv}\|\beta_{i}\|,$ (6)
$\displaystyle d_{tg}$
$\displaystyle=\rho_{i}+\sqrt{(x_{f}-x^{w}_{i})^{2}+(y_{f}-y^{w}_{i})^{2}},$
$\displaystyle\mathbf{g}_{lnp}^{*}$
$\displaystyle=\operatorname{arg}\min_{\mathbf{g}_{lnp_{i}}\in\text{LNPs}}\left(C_{g}\left(f_{i}\right)\right),$

where $k_{dst}$, $k_{dir},k_{elv}$ are weighting factors. The go-to-goal cost
$C_{g}\left(\mathbf{g}_{lnp_{i}}\right)$ is applied to the G-LNPs, where in
the case of the V-LNPs the navigability $\iota$ is used to mask the cost of
the visually non-navigable LNPs to the maximum as follow,

$\begin{gathered}C_{v}\left(\bm{\nu}_{lnp_{i}}\right)=\begin{cases}C_{g}\left(\mathbf{g}_{lnp_{i}}\right)&\text{if
}\iota_{i}=255,\\\ 1&\text{otherwise.}\end{cases}\end{gathered}$ (7)

Once the optimal LNP ($\mathbf{g}_{lnp}^{*}$ / $\mathbf{\nu}_{lnp}^{*}$) is
selected, a motion command with linear and angular velocities, $v$ and
$\omega$, is generated as follows: $v=k_{a}\rho_{*}-k_{b}\|\alpha_{*}\|$, and
$\omega=k_{c}\alpha_{*}$. $\rho_{*}$ and $\alpha_{*}$ are the attributes of
the optimal LNP ($\mathbf{g}_{lnp}^{*}$ / $\mathbf{\nu}_{lnp}^{*}$), where
$k_{a},\ k_{b}$ and $k_{c}$ are tunable coefficients.

###  II-D Joint Visual-Geometry Model

Our integration methodology excludes the D-SGP model from the V-SGP model,
thereby establishing two distinct SGP Models: the unaltered G-SGP model and
the N-SGP model. This process essentially substitutes the D-SGP model with a
segment of the G-SGP model, specifically the area where the FOVs of the camera
and LiDAR overlap, see Fig. LABEL:fig_jnt_mdl_b. The operational workflow
begins with the computation of G-LNPs derived from the G-SGP model. Following
this, the N-SGP model is engaged to determine the navigability of the
identified G-LNPs to define the Visual G-LNPs (VG-LNPs). For goal-oriented
navigation task, one VG-LNP is selected based on the same go-to-goal cost
function introduced in Eq. 6, however the navigability mask in Eq. 7 is
modified in such away that each G-LNP is allocated a navigability score as
follows,

$\begin{gathered}C_{vg}\left(\bm{\nu
g}_{lnp_{i}}\right)=\begin{cases}(1-k_{nav})\>C_{g}\left(\mathbf{g}_{lnp_{i}}\right)&\text{if
}\iota_{i}=255,\\\ 1&\text{if }\iota_{i}=0,\\\
k_{nav}\>C_{g}\left(\mathbf{g}_{lnp_{i}}\right)&\text{Out of camera's
FoV.}\end{cases}\end{gathered}$ (8)

where $0.5<k_{nav}<1$ is a parameter allowing the user to control the
preference between visually navigable VG-LNPs and G-LNPs outside the camera’s
FoV. When $k_{nav}=0.5$, there exists no preference between the two, granting
them an equal cost of $0.5\,C_{g}$. Conversely, setting $k_{nav}>0.5$ augments
the propensity to opt for VG-LNPs situated within the camera’s FoV.

##  III Experimental Design and Results

###  III-A Experimental Setup

A Clearpath Jackal robot equipped with a VLP-16 velodyne and Realsense D435
camera is used to validate our algorithm in both simulation and real-world
scenarios. For the simulation experiments, we used the grass-mud environment
[37] with ground-truth localization and an RGB-based segmentation to generate
the navigability image from the RGB image. While for real hardware
experiments, we used a microstrain GNSS/INS module for localization and the
mask2former [9] segmentation. The Velodyne pointcloud $\mathbf{P}_{l}$ is used
to construct the occupancy surface $\mathbf{S}_{g}$ to train the G-SGP model,
where the front Realsense D435 RGBD camera is used to generate the
navigability poincloud $\mathbf{P}_{v}$ and construct the visual surface
$\mathbf{S}_{v}$ to train the V-SGP model. Our algorithm runs in real-time
with a frequency of 4 Hz in real-world experiment and 10 Hz in simulation due
to the time difference required by RGB-Based and the mask2former segmentation.
Therefore, we limited the maximum velocity in the realworld experiment to 0.5
m/s. We assessed our algorithm’s effectiveness by comparing it to two
established methods. The first baseline is a purely visual-based navigation
named PovNav [38], that employs $I_{rgb}$ for generating $I_{nav}$, similar to
our proposed method. However, it uses image-based visual-servoing for motion
command generation. The second baseline is GPFrontiers framework [32, 34],
local mapless navigation which contrasts PovNav by focusing solely on
geometry-based navigation. For clarity, we refer to PovNav as V-Nav,
GPFrontiers as G-Nav, and our proposed approach as VG-Nav.

###  III-B Simulation Scenarios and Results

In the grass-mud environment, two visual classes are identified: grass and
mud, with grass assigned as navigable and mud as non-navigable. Furthermore,
the environment features two types of terrain: flat terrain and high slope
grass (HSG), with HSG being geometrically non-traversable. Our experimental
design request the robot to navigate from a start pose of $(15,-4,-\pi/2)$ to
a goal pose of $(-4,-16,\pi)$, avoiding mud and HGS regions. Initially facing
HSG, the robot must circumvent it, as well as a mud area presented along the
straight path to the goal. Successful task completion necessitates the
integration of both visual and geometric navigation capabilities to sidestep
HSG and mud obstacles. For each approach, we conducted 15 trials to examine
its navigation behavior. First, we executed the G-Nav, which, as anticipated,
avoided the HSG area and opted for a direct path to the goal, traversing the
mud area due to lack of a visual component. The paths generated by G-Nav are
depicted in red in Fig. 5. Subsequently, we deployed the V-Nav which failed to
circumvent the HSG, perceiving all grass as navigable due to missing the
geometric information. All 15 V-Nav trials ended with the robot stuck trying
to climb the HSG, depicted as black paths in Fig. 5. To evaluate the V-Nav
performance on flat terrain, we reoriented the robot to face flat areas
instead of HSG. In this revised configuration (V-Nav:Flat), V-Nav successfully
navigated through flat grass, avoiding mud areas in 12 out of 15 trials.
However, in 3 trials, V-Nav failed to circumvent mud due to erroneous motion
commands. The critical observation was that V-Nav, when detecting only non-
navigable classes (mud) within the camera’s FoV, tended to reach a local
minimum, failing to avoid going into non-navigable areas, check the blue paths
in Fig. 5.

(a)

(b)

Figure 5: Simulation Experiments: VG-LNPs cost on the right figure. $K_{nav}>0.5$ giving less cost for visually VG-LNPs than G-LNPS outside the Camer’a FoV.  TABLE I: Navigation metrics. (M: Mud and S: Slope)  | Path[m] | $avg(v_{max})[m/s]$ | Avoid M | Avoid S
---|---|---|---|---
G-Nav | $22.54\pm 0.15$ | $1.06\pm 0.10$ | 0% | 100%
V-Nav:HSG | - | - | - | 0%
V-Nav:Flat | $33.01\pm 0.30$ | $0.30\pm 0.01$ | 80% | -
VG-Nav | $29.21\pm 0.46$ | $1.03\pm 0.06$ | 100% | 100%

This behavior is replicated in real-world experiments as well, see Fig.
LABEL:fig_exp3_a. Finally, we tested our proposed VG-Nav in the original setup
with the robot facing HSG. VG-Nav successfully guided the robot to the goal
while processing both geometric and visual information, effectively avoiding
both HSG and mud areas. The resulting paths are illustrated in cyan in Fig. 5.

V-Nav tries to follow the wide visually-navigable space on the navigability
image, resulting in a longer path length of $~{}33$ m, in contrast to VG-Nav,
which prioritizes the goal direction unless avoiding mud areas, achieving a
shorter average path length of $~{}29.2$ m. Moreover, the average of the
achieved maximum velocity over the 15 trials was $1.06m/s$ for VG-Nav and
$0.3m/s$ for V-Nav. VG-Nav, by utilizing a 3D navigability point cloud instead
of a 2D navigability image, is more robust in avoiding non-navigable classes,
as evidenced in Table I.

###  III-C Real World Scenarios and Results

During real-world experiments, mask2former [9] is used for segmentation. While
the RGB-based segmentation was robust in simulation for distinguishing mud
from grass, mask2former shows occasional confusion among similar classes:
earth, ground, and grass. This issue is addressed by assigning these four
classes identical navigability values. Notably, mask2former reliably
distinguished these classes from asphalt/road. Three experiments were
conducted to validate the outcomes of the simulation studies. In the first
setup, the robot is instructed to move towards a goal positioned in a straight
line ahead of it. Despite the straightforward path, opting for this route
would necessitate crossing over a grassy area, where grass is defined as non-
navigable and asphalt/road is defined as navigable. Subsequently, we conducted
3-trials for each approach. Fig. 6(left) illustrates the different paths
achieved under each algorithm. Notably, the G-Nav failed to avoid the grass,
lacking the component to recognize the grass as a non-navigable zone. In
contrast, both V-Nav and VG-Nav successfully navigated around the grassy
region, with VG-Nav opting for a shorter path closer to the grass boundary.
These results corroborate the findings from the simulation experiments.

(a)

(b)

Figure 6: First real-world experiment: G-Nav can not avoid mud; G-Nav(yellow),
V-Nav(green), VG-Nav(cyan). $K_{nav}=0.5$.

(a)

(b)

(c)

(d)

Figure 7: Second real-world experiment: V-Nav can not avoid HSG.

(a)

(b)

(c)

(d)

Figure 8: Third real-world experiment: V-Nav local minimum due to limited FoV;
V-Nav(cyan), VG-Nav(green).

Fig.6(right) visualizes the cost associated with VG-LNPs. VG-LNPs situated
within the non-navigable part of the visual point cloud are assigned the
highest cost whereas the remaining VG-LNPs are assigned costs based on the go-
to-goal cost function according to Eq.8 with $K_{nav}=0.5$.

In the second experiment, the robot was positioned in an HSG area with grass
classified as navigable. The robot was directed towards a goal straight ahead,
yet its path passes through an area of grass with the steepest slope. We
conducted 3 trials using V-Nav and VG-Nav. VG-Nav successfully identified HSG
and avoided it. As illustrated in Fig.LABEL:fig_exp2_b, which corresponds to
the local observation shown in Fig.LABEL:fig_exp2_a, no VG-LNPs are present on
the robot’s right side that leads to HSG, instead, VG-LNPs are positioned
forward, backward, and to the left of the robot, in zones with elevation
angles bounded by $\beta_{max}$. Conversely, V-Nav attempted to guide the
robot through the HSG along the straight path. Fig.LABEL:fig_exp2_d depicts
the visual path planned based on the local observation in
Fig.LABEL:fig_exp2_c. For V-Nav, we ended the experiment before the robot goes
to the steepest slope area. These outcomes further validate the insights from
the simulation experiments. The third experiment replicated the local minimum
behavior observed with V-Nav. In this setup, the robot heading was initiated
to face a corner formed by a patch of grass (assigned as non-navgiable),
deliberately testing the constraints imposed by a limited visual FoV. Further
complicating the path, a table was placed on the direct path to the goal, to
introduce a geometric obstacle. Three trials were conducted, where V-Nav
failed to get out of the grass corner in all trials due to the local minimum
it encounters when the FoV contains only the non-navigable class, check the
cyan path in Fig LABEL:fig_exp3_b. In contrast, the VG-Nav demonstrated
superior adaptability, managing to identify geometrically feasible G-LNPs
outside the camera’s FoV. All VG-LNPs in the camera’s FoV were masked to the
maximum cost, see Fig LABEL:fig_exp3_b, therefore the VG-Nav command the robot
to follow the G-LNP with the minimum cost which lies outside the camera’s FOV,
this is interpreted from the high angular and low linear velocities values on
the first few seconds which led to a sharp rightward rotation, see Fig
LABEL:fig_exp3_c.

We also analyzed the computational costs associated with implementing our VG-
Nav by computing the time cost for each individual observation. Fig.
LABEL:fig_exp3_d shows the time costs for each component and their sum as the
the total time. The average training time of both the G-SGP and the V-SGP
models is around $19msec$, while the average prediction time of the G-SGP
($24msec$) is higher than the average navigability prediction time of the
V-SGP model ($17.5msec$) because V-SGP only predicts the navigability of the
feasible G-LNPs identified by the G-SGP model. The average total time for each
observation is $81msec$ for the VG-Nav and $43msec$ for G-Nav, while it is
$9msec$ for V-Nav. It is worth mentioning that the time cost for the VG-Nav
does not explode over the operation time (number of observations) because we
consider each observation individually. A video highlighting the real-world
experiments results is available at: https://youtu.be/0s6VSj5Z1dg

Contrary to the end-to-end approach, which does not allow the analysis of the
individual contributions from geometric and visual information, VG-Nav
distinctly employs G-SGP to invalidate all geometrically non-visible G-LNPs
and V-SGP to allocate visual traversability costs to the remaining VG-LNPs.
Moreover, our methodology is capable of integrating future advanced semantic
segmentation techniques to immediately improve navigation behavior without
necessitating a training phase. In the current implementation of VG-Nav,
traversability for each terrain type is determined manually, reflecting either
the robot’s capabilities or preferred navigation behavior, similar to [10, 11,
12]. Future developments aim to enhance the V-SGP model’s capacity to learn
terrain traversability, facilitating adaptation across diverse terrain types
in a similar way as [13, 14, 15, 16, 17].

##  IV Conclusion

We present the Visual Geometry Combined Spaces (VG-SGP) model along with a
corresponding navigation strategy (VG-Nav) designed to adeptly guide a robot
to its destination. This method leverages the environment analysis of two
separate SGP models, G-SGP and V-SGP, to pinpoint areas that are navigable
based on both visual and geometric characteristics in the robot’s surrounding
environment. By integrating visual and geometric information, our approach
facilitates more reliable and adaptive navigation. Simulation and real-world
experiments demonstrate that our VG-SGP model and its integrated navigation
strategy outperform systems reliant solely on either visual or geometric
navigation algorithms, showcasing superior adaptive behavoir required for
accomplishing flexible tasks.

## References

  * [1] K. Weerakoon, A. J. Sathyamoorthy, U. Patel, and D. Manocha, “Terp: Reliable planning in uneven outdoor environments using deep reinforcement learning,” in 2022 International Conference on Robotics and Automation (ICRA), pp. 9447–9453, IEEE, 2022.
  * [2] J. Liang, K. Weerakoon, T. Guan, N. Karapetyan, and D. Manocha, “Adaptiveon: Adaptive outdoor local navigation method for stable and reliable actions,” IEEE Robotics and Automation Letters, vol. 8, no. 2, pp. 648–655, 2022.
  * [3] S. Thrun, “Learning occupancy grid maps with forward sensor models,” Autonomous robots, vol. 15, pp. 111–127, 2003.
  * [4] P. Fankhauser and M. Hutter, “A Universal Grid Map Library: Implementation and Use Case for Rough Terrain Navigation,” in Robot Operating System (ROS) – The Complete Reference (Volume 1) (A. Koubaa, ed.), ch. 5, Springer, 2016.
  * [5] A. Hornung, K. M. Wurm, M. Bennewitz, C. Stachniss, and W. Burgard, “Octomap: An efficient probabilistic 3d mapping framework based on octrees,” Autonomous robots, vol. 34, pp. 189–206, 2013.
  * [6] J. Shin, D. Kwak, and K. Kwak, “Model predictive path planning for an autonomous ground vehicle in rough terrain,” International Journal of Control, Automation and Systems, vol. 19, pp. 2224–2237, 2021.
  * [7] Y. Zhu, R. Mottaghi, E. Kolve, J. J. Lim, A. Gupta, L. Fei-Fei, and A. Farhadi, “Target-driven visual navigation in indoor scenes using deep reinforcement learning,” in 2017 IEEE international conference on robotics and automation (ICRA), pp. 3357–3364, IEEE, 2017.
  * [8] W. Yang, X. Wang, A. Farhadi, A. Gupta, and R. Mottaghi, “Visual semantic navigation using scene priors,” arXiv preprint arXiv:1810.06543, 2018.
  * [9] B. Cheng, I. Misra, A. G. Schwing, A. Kirillov, and R. Girdhar, “Masked-attention mask transformer for universal image segmentation,” in Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pp. 1290–1299, 2022.
  * [10] T. H. Y. Leung, D. Ignatyev, and A. Zolotas, “Hybrid terrain traversability analysis in off-road environments,” in 2022 8th International Conference on Automation, Robotics and Applications (ICARA), pp. 50–56, IEEE, 2022.
  * [11] D. Maturana, P.-W. Chou, M. Uenoyama, and S. Scherer, “Real-time semantic mapping for autonomous off-road navigation,” in Field and Service Robotics: Results of the 11th International Conference, pp. 335–350, Springer, 2018.
  * [12] T. Guan, Z. He, R. Song, D. Manocha, and L. Zhang, “Tns: Terrain traversability mapping and navigation system for autonomous excavators,” arXiv preprint arXiv:2109.06250, 2021.
  * [13] G. Kahn, P. Abbeel, and S. Levine, “Badgr: An autonomous self-supervised learning-based navigation system,” IEEE Robotics and Automation Letters, vol. 6, no. 2, pp. 1312–1319, 2021.
  * [14] X. Cai, M. Everett, J. Fink, and J. P. How, “Risk-aware off-road navigation via a learned speed distribution map,” in 2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 2931–2937, IEEE, 2022.
  * [15] M. V. Gasparino, A. N. Sivakumar, Y. Liu, A. E. Velasquez, V. A. Higuti, J. Rogers, H. Tran, and G. Chowdhary, “Wayfast: Navigation with predictive traversability in the field,” IEEE Robotics and Automation Letters, vol. 7, no. 4, pp. 10651–10658, 2022.
  * [16] L. Nardi and C. Stachniss, “Actively improving robot navigation on different terrains using gaussian process mixture models,” in 2019 International Conference on Robotics and Automation (ICRA), pp. 4104–4110, IEEE, 2019.
  * [17] G. G. Waibel, T. Löw, M. Nass, D. Howard, T. Bandyopadhyay, and P. V. K. Borges, “How rough is the path? terrain traversability estimation for local and global path planning,” IEEE Trans. on Intelligent Transportation Sys., vol. 23, no. 9, pp. 16462–16473, 2022.
  * [18] A. Bernardino and J. Santos-Victor, “Visual behaviours for binocular tracking,” in The Second EUROMICRO Workshop on Advanced Mobile Robots, pp. 2–7, 1997.
  * [19] C. Boretti, P. Bich, Y. Zhang, and J. Baillieul, “Visual navigation using sparse optical flow and time-to-transit,” in 2022 Int. Conf. on Robotics and Automation (ICRA), pp. 9397–9403, 2022.
  * [20] S. R. Bista, P. R. Giordano, and F. Chaumette, “Appearance-based indoor navigation by ibvs using line segments,” IEEE Robotics and Automation Letters, vol. 1, no. 1, pp. 423–430, 2016.
  * [21] A. Dame and E. Marchand, “A new information theoretic approach for appearance-based navigation of non-holonomic vehicle,” in 2011 IEEE Int. Conf. on Robotics and Automation, pp. 2459–2464, 2011.
  * [22] D. Kim, D. Lee, H. Myung, and H. Choi, “Object detection and tracking for autonomous underwater robots using weighted template matching,” in 2012 Oceans - Yeosu, pp. 1–5, 2012.
  * [23] I. S. Mohamed, G. Allibert, and P. Martinet, “Model predictive path integral control framework for partially observable navigation: A quadrotor case study,” in 2020 16th Int. Conf. on Control, Automation, Robotics and Vision (ICARCV), pp. 196–203, IEEE, 2020.
  * [24] M. Zhang, Y. Li, and J. Yang, “Autonomous visual navigation guided by path boundaries for mobile robot,” in 2008 International Conference on Computer Science and Software Engineering, vol. 6, pp. 344–348, 2008.
  * [25] M. G. Jadidi, J. V. Miro, and G. Dissanayake, “Warped gaussian processes occupancy mapping with uncertain inputs,” IEEE Robotics and Automation Letters, vol. 2, no. 2, pp. 680–687, 2017.
  * [26] M. G. Jadidi, L. Gan, S. A. Parkison, J. Li, and R. M. Eustice, “Gaussian processes semantic map representation,” arXiv preprint arXiv:1707.01532, 2017.
  * [27] E. Snelson and Z. Ghahramani, “Sparse gaussian processes using pseudo-inputs,” Advances in neural information processing systems, vol. 18, p. 1257, 2006.
  * [28] R. Sheth, Y. Wang, and R. Khardon, “Sparse variational inference for generalized gp models,” in International Conference on Machine Learning, pp. 1302–1311, PMLR, 2015.
  * [29] M. Titsias, “Variational learning of inducing variables in sparse gaussian processes,” in Artificial intell. and statist., pp. 567–574, PMLR, 2009.
  * [30] H. Lee, J. Kwon, and C. Kwon, “Learning-based uncertainty-aware navigation in 3d off-road terrains,” in 2023 IEEE Int. Conf. on Robotics and Automation (ICRA), pp. 10061–10068, IEEE, 2023.
  * [31] M. Ali and L. Liu, “Light-weight pointcloud representation with sparse gaussian process,” in 2023 IEEE International Conference on Robotics and Automation (ICRA), pp. 4931–4937, 2023.
  * [32] M. Ali and L. Liu, “Gp-frontier for local mapless navigation,” in 2023 IEEE International Conference on Robotics and Automation (ICRA), pp. 10047–10053, 2023.
  * [33] M. Ali, H. Jardali, N. Roy, and L. Liu, “Autonomous navigation, mapping and exploration with gaussian processes.,” in Robotics: Science and Systems, 2023.
  * [34] H. Jardali, M. Ali, and L. Liu, “Autonomous mapless navigation on uneven terrains,” arXiv preprint arXiv:2402.13443, 2024.
  * [35] V. Sezer and M. Gokasan, “A novel obstacle avoidance algorithm:“follow the gap method”,” Robotics and Autonomous Systems, vol. 60, no. 9, pp. 1123–1134, 2012.
  * [36] M. Mujahed, D. Fischer, and B. Mertsching, “Admissible gap navigation: A new collision avoidance approach,” Robotics and autonomous systems, vol. 103, pp. 93–110, 2018.
  * [37] H. Lee, J. Kwon, and C. Kwon, “Learning-based uncertainty-aware navigation in 3d off-road terrains,” in 2023 IEEE Int. Conf. on Robotics and Automation (ICRA), pp. 10061–10068, IEEE, 2023.
  * [38] D. Pushp, Z. Chen, C. Luo, J. M. Gregory, and L. Liu, “Povnav: A pareto-optimal mapless visual navigator,” arXiv preprint arXiv:2310.14065, 2023.


# Resolvent analysis of a swimming foil

J. M. O. Massey1,2,3<EMAIL_ADDRESS>S. Symon 1 B. Ganapathisubramani1
G. D. Weymouth1,2  1Faculty of Engineering and Physical Sciences, University
of Southampton, UK 2Faculty of Mechanical Engineering, TU Delft, NL
3Department of Mechanical Engineering, Stanford University, Stanford, CA
94305, USA

###### Abstract

This study employs resolvent analysis to explore the dynamics and coherent
structures in the boundary layer of a foil that swims via a travelling wave
undulation. A modified NACA foil shape is used together with undulatory
kinematics to represent fish-like bodies at realistic Reynolds numbers
($\Rey=10,000$ and $\Rey=100,000$) in both thrust- and drag-producing
propulsion regimes. We introduce a novel coordinate transformation that
enables the implementation of the data-driven resolvent analysis (Herrmann et
al., 2021) to dissect the stability of the boundary layer of the swimming
foil. This is the first study to implement resolvent analysis on deforming
bodies with non-zero thickness and at realistic swimming Reynolds numbers. The
analysis reinforces the notion that swimming kinematics drive the system’s
physics. In drag-producing regimes, it reveals breakdown mechanisms of the
propulsive wave, while thrust-producing regimes show a uniform wave
amplification across the foil’s back half. The key thrust and drag mechanisms
scale with the boundary-layer thickness, implying geometric self-similarity in
this $\Rey$ regime. In addition, we identify a mechanism that is less strongly
coupled to the body motion. We offer a comparison to a rough foil that reduces
the amplification of this mechanism, demonstrating the potential of roughness
to control the amplification of key mechanisms in the flow. The results
provide valuable insights into the dynamics of swimming bodies and highlight
avenues for developing opposition control strategies.

###### keywords:

##  1 Introduction

The study of the boundary layer of swimming bodies has been addressed in the
literature by Anderson et al. (2001), who investigated the boundary layer of a
waving plate under various swimming conditions. They were unable to provide
information on the dynamics or coherent structures in the boundary layer,
information crucial to understanding the system (Rempfer & Fasel, 1994).

Swimming involves creating a propulsive wave through a time-periodic body
deformation. Resolvent analysis (McKeon & Sharma, 2010) has had vast success
in uncovering driving mechanisms in complex flows. The analysis treats the
non-linearity in the fluctuating part of the Navier-Stokes equation as an
unknown harmonic forcing. However, the traditional formulation linearises
around the time-averaged mean flow to study the transfer of energy from the
mean to dominant fluctuations. As such, only the first harmonic is identified
in oscillator-type flows, whereas higher harmonics are not, since they receive
energy from the first harmonic and not the mean (Symon et al., 2019). To
improve predictions at higher harmonics, the linearisation can be performed
around an alternative base flow. For example, Padovan et al. 2020 linearised
around a Fourier basis to offer insights into the unstable dynamics of the
flow over a NACA0012 aerofoil at $\Rey=200$ and an angle of attack of
$20^{\circ}$.

Herrmann et al. (2021) adopts a data-driven approach to obtain the resolvent
operator, that uses the link between Dynamic Mode Decomposition (DMD) and
Koopman eigenfunctions (Williams et al., 2015) to construct the linear basis
for analysis. DMD, introduced by Schmid (2010), enables the extraction of
spatio-temporal patterns from time-resolved data. The challenge of a moving
and deformable body is that the domain position containing the fluid and the
body depends on time. Goza & Colonius (2018) successfully coupled the fluid
and body domains to identify the driving mechanisms for the onset of chaotic
flapping, however, their body was an infinitesimally thin membrane, so they
did not have to consider the temporal domain switch between the fluid and body
regions. Menon & Mittal (2020) looked at the pitching and plunging of an
aerofoil by rotating and translating the reference frame, so the body remained
centred within the reference frame. This is limited in that it deals with a
rigid body.

In this paper, we investigate the dominant mechanisms in the boundary layer of
a swimming foil using the data-driven resolvent analysis (Herrmann et al.,
2021). We use a modified NACA foil shape to represent a fish-like body and a
generalised swimming motion to represent the kinematics. The foil swims at
both $\Rey=10,000$ and $\Rey=100,000$. We consider thrust- and drag-producing
propulsive regimes. We propose a novel coordinate transform to study the
boundary layer of a swimming foil in body coordinates. Finally, we provide a
comparison to a rough foil that reduces the amplification of a mechanism at
the third harmonic of the swimming frequency. To our knowledge, this is the
first study of its kind to implement resolvent analysis on deforming bodies
with nonzero thickness. In addition, we uncover key physical differences
between thrust- and drag-producing foils and highlight coherent structures
that scale with the boundary-layer thickness.

##  2 Methodology

###  2.1 Geometry and kinematics

Figure 1: Illustration of the computational set up and validation. (a)
Comparison of the NACA0012 profile in grey and the modified profile plotted on
top (black). (b) Schematic of the domain and grid. The inner box shows the
region where the grid is uniformly rectilinear; from there the grid stretches
toward the domain extent. The inset shows the grid resolution around the tail
of the foil. (c) Vortex structures in the wake visualised by the $1\times
10^{-5}$ contour of the Q-criterion, coloured by the y-vorticity.

We use a NACA0012 profile, modified with a sharpened leading edge, to
represent a fish-like body. Figure 1 compares the modified profile with the
NACA0012 excursion data.

Scaling all lengths by $L$, using the swimming speed $U$ to scale the velocity
and $L/U$ to scale the time, we define the excursion of the body from the
centre line as

$y(x)=A(x)\sin\big{(}2\pi\,[ft-x/\zeta]\big{)},$ (1)

where $f$ is the frequency, $\zeta$ is the phase speed of the travelling wave,
and $A(x)$ is the amplitude envelope. The Strouhal number is set to a value of
$St=0.3$, which determines the scaled frequency as $f=St/2A_{1}$, where
$A_{1}=A(x=1)$ is the trailing edge amplitude and defines the wake width of
the system. We use the result from Di Santo et al. (2021) for the envelope

$A(x)=\frac{A_{1}(a_{2}x^{2}+a_{1}x+a_{0})}{\sum_{i=0}^{2}a_{i}},$ (2)

using $a_{0,1,2}=(0.05,-0.13,0.28)$. $A_{1}=0.1$, which was found to be
optimal in Saadat et al. (2017). Finally, we set the wave speed to place the
foil in a thrust or drag-producing regime (figure 2 a). For the thrust-
producing foil $\zeta=1.25$ and for the drag $\zeta=1.75$. For the convergence
analysis, we use a representative simulation using $\zeta=1.42$, which we find
to result in a cycle averaged zero net thrust.

###  2.2 Numerical method

We simulate incompressible fluid flow with the dimensionless Navier-Stokes
combined with the continuity equation

$\displaystyle\frac{\partial\mathbf{u}}{\partial
t}+(\mathbf{u}\cdot\mathbf{\nabla})\mathbf{u}=$
$\displaystyle-\mathbf{\nabla}p+\frac{1}{\Rey}\nabla^{2}\mathbf{u},$ (3)
$\displaystyle\mathbf{\nabla}\cdot\mathbf{u}=$ $\displaystyle 0\text{ ,}$ (4)

where $\mathbf{u}(\mathbf{x},t)=(u,v,w)$ is the scaled velocity of the flow
and $p(\mathbf{x},t)$ is the scaled pressure. We use an in-house developed
solver based on the Boundary Data Immersion Method (BDIM) (Maertens &
Weymouth, 2015). The solver has been extensively validated in swimming studies
(Maertens et al., 2017; Zurman-Nasution et al., 2020; Massey et al., 2023;
Lauber et al., 2023) and converges at second order in both time and space
(Lauber et al., 2022).

The boundary conditions we enforce on the body is the no-slip condition.
Symmetry conditions are enforced on the upper and lower domain extents, and
there is a periodic condition in the spanwise direction. Figure 1 details the
grid and domain configuration. We use the smallest domain for which the forces
on the body remain invariant when the size increases. We use the thrust force
based on the surface pressure, $C_{T}=\oint\mathbf{f_{x}}ds/0.5A_{1}$, where
$\mathbf{f}=-p\hat{n}$ is the normal pressure stress on the body.

Figure 2: (a) Cycle normalised thrust coefficients for the 4 test cases. The
red lines are the thrust producing cases and the purple are the drag producing
cases. $\Rey$s are differentiated by the intensity of the colour, as labelled.
(b) Convergence of the standard deviation of the thrust coefficient for a
representative simulation with zero net thrust. The error of both 2-D and 3-D
simulations are relative to a 2-D simulation with $\Delta y=1/8192$.

We use a rectilinear grid in the domain $x\in[-1.5,5.5]$, $y\in[-1.5,1.5]$,
and $z\in[0,1/64]$ (figure 1a). For the area that contains the body motion and
the immediate wake, we use a uniform grid and then implement hyperbolic
stretching of the grid cells away from this area (figure 1). The grid is
refined in $y$ such that $\Delta y$ is a quarter of $\Delta x$ and $\Delta z$
(figure 1 inset). The total number of grid cells is $N=(4096,4096,64)$
totalling $N=1.07\times 10^{9}$ grid cells.

The grid resolution is shown to be sufficient by using representative two-
dimensional and three-dimensional simulations. Using 2-D simulations allows us
to increase the simulation comparison point one extra refinement step and
check the resolution at the final working resolution is sufficient. Figure 1c
shows that the flow has little variation in the spanwise direction similar to
the results shown in Zurman-Nasution et al. (2020); Massey et al. (2023).
Therefore, we use a highly refined two-dimensional simulation as a final
comparison point for both the 2-D and 3-D cases, avoiding the need for the
extremely large fine grid in 3-D. Figure 2 b shows that the standard deviation
of the thrust coefficient converges for both 2-D and 3-D. The domain size is
validated by comparing the thrust with the phase angle for the working domain
and a much larger domain with dimensions $x\in[-5,20]$ and $y\in[-5,5]$; the
thrust is invariant with the phase angle (not shown), indicating that the
domain size is sufficient. Finally, we save 800 snapshots of the velocity
field over 4 cycles for the subsequent analysis.

###  2.3 Coordinate transform

Figure 3: The velocity magnitude in (a) Cartesian coordinates from the
computational domain and (b) the local normal-tangential coordinates on the
upper section of the foil.

DMD produces inaccurate modes around a body with motion (Menon & Mittal,
2020). This is because the temporal switch between the body and the fluid
domain causes the DMD to scatter the coherent structures in the boundary layer
across the spatial domain. The analysis can still be performed on the wake but
misses valuable insights into the boundary-layer dynamics. Instead, we
transform the flow into local normal-tangential coordinates to investigate the
stability of the boundary layer of the swimming foil. This enables us to
identify spatially and temporally coherent structures through DMD.

Addressing the limitations of Dynamic Mode Decomposition (DMD) around moving
bodies, we employ a transformation to local normal-tangential coordinates, as
illustrated in figure 3. For simplicity, we spanwise average the flow. This
approach results in a loss of information about spanwise mechanisms; however,
we consider this loss to be minimal due to the small spanwise variability on
the foil (figure 1b), and the fact that spanwise fluctuations occur
predominantly in the wake of the foil. The global Cartesian coordinates are
transformed into a local normal-tangential coordinate system through a mapping
$\mathcal{T}$ such that

$\mathcal{T}:\mathbf{x}\rightarrow\mathbf{\xi}\text{ ,}$ (5)

where the Cartesian and local normal-tangential coordinates are

$\mathbf{x}=\begin{pmatrix}x\\\
y\end{pmatrix},\quad\mathbf{\xi}=\begin{pmatrix}\xi\\\ \eta\end{pmatrix}\text{
,}$ (6)

with $\xi$ aligned tangentially and $\eta$ aligned normally to the foil
surface. In this analysis, we will only consider the top boundary layer of the
foil, so the transformation can be determined by considering a superposition
between the polynomial describing the shape of the foil and the equation 1
describing the motion. The position of the body boundary is then defined by
$\mathbf{d}(t)$. The unit normal vector
($\hat{n}=\nabla\mathbf{d}/\left\|\nabla\mathbf{d}\right\|$) is calculated
with respect to spatial coordinates, establishing a normal direction at each
point along the body’s surface, thus facilitating the mapping
($\mathcal{T}(\mathbf{d},\eta)$) between Cartesian and local normal-tangential
frames

$\mathcal{T}=\mathbf{d}+\hat{n}\eta\text{ ,}$ (7)

which satisfies $\xi=x|_{\eta=0}$ on the body. The normal component ($\eta$)
is uniformly sampled, simplifying the interpolation process without
significantly impacting the accuracy of the analysis. The velocity field in
the local normal-tangential coordinates ($\mathbf{u}_{\mathcal{B}}$) is
obtained by an interpolation of the velocity in the Cartesian coordinates
($\mathbf{u}$) at the local normal-tangential sample points

$\mathbf{u}_{\mathcal{B}}=\mathcal{I}\left(\mathbf{u}_{\text{nn}},(\eta,\xi)\right)\text{
,}$ (8)

where $\mathcal{I}$, a linear interpolation operator chosen for its
computational efficiency and accuracy, maps the nearest-neighbour Cartesian
velocity field ($\mathbf{u}_{\text{nn}}$) to the local normal-tangential
coordinates. This transformed velocity field then underpins the resolvent
analysis, enabling focused investigation into boundary layer dynamics and the
identification of coherent structures. Throughout this study, we will refer to
the transformed velocity field without the subscript $\mathcal{B}$ for the
sake of brevity.

###  2.4 Resolvent analysis

We use the method of Herrmann et al. 2021 for the resolvent analysis. Briefly,
we consider a forced dynamical system represented by a divergence free
velocity subspace such that

$\dot{\mathbf{u}}=\mathcal{L}\mathbf{u}+\mathbf{f}\text{ ,}$ (9)

where $\mathbf{u}$ is the velocity field, $\mathcal{L}$ is a dynamical
operator driving the evolution of the system, and $\mathbf{f}$ is the forcing
term. To determine the linear operator, we consider the evolution of the
velocity field in discrete time such that

$\mathbf{U}_{k+1}=A\mathbf{U}_{k}\text{ ,}$ (10)

where $\mathbf{U}_{k}$ is the state at time $t_{k}=k\Delta t$. For
convenience, we use $\Delta t=0.005$, which is the flow field sampling
frequency–this rescales time to normalise the definition in section 2.2 by the
flapping frequency. To approximate the operator $A$, we implement the standard
DMD algorithm from Schmid (2010) using the open source software package Demo
et al. (2018). This algorithm is sufficient for our dataset and is able to
handle the number of snapshots we have. We construct flat overlapping matrices
$X$ and $Y$ such that

$\displaystyle X$
$\displaystyle=\begin{bmatrix}u_{0}&u_{1}&\cdots&u_{m-1}&v_{0}&v_{1}&\cdots&v_{m-1}\end{bmatrix}$
(11) $\displaystyle Y$
$\displaystyle=\begin{bmatrix}u_{1}&u_{2}&\cdots&u_{m}&v_{1}&v_{2}&\cdots&v_{m}\end{bmatrix}\text{
,}$ (12)

and employ a dimensionally efficient projection using the pseudoinverse to
approximate the linear map. We take the truncated SVD of the matrix $X$,
producing $X=U\Sigma V^{*}$ and estimate the linear map such that
$A=U^{*}YV\Sigma^{-1}$. For this study, we truncate the SVD at 6 as we find
that this is the largest value that does not introduce decaying modes into the
bases. The DMD modes are then computed by taking the eigendecomposition of
$A$, such that $AV=V\Gamma$ and $A^{+}W=W\Gamma^{*}$, where $A^{+}$ is the
adjoint of $A$.

For the resolvent analysis, we consider equation 9 in terms of eigenvector
coordinates such that $x(t)=Va(t)$ and $f(t)=Vb(t)$, with $a,b$ being
expansion coefficients. In continuous time, $\Lambda=\ln(\Gamma)/\Delta t$, so
the equation 9 can be written as

$\dot{a}=\Lambda a+b\text{ .}$ (13)

To maintain the 2-norm’s physical meaning, we adjust the inner product via a
weighting matrix defined from the Cholesky factorisation of
$V^{*}V=\widetilde{F}^{*}\widetilde{F}$. The resolvent operator is then
defined as

$\widetilde{F}(-i\omega
I-\Lambda)^{-1}\widetilde{F}^{-1}=\Psi_{\widetilde{F}}\Sigma\Phi^{*}_{\widetilde{F}}\text{
,}$ (14)

where we have decomposed the resolvent operator into the singular value
decomposition, $\Sigma$, and the forcing and response modes,
$\Phi^{*}_{\widetilde{F}}$ and $\Psi_{\widetilde{F}}$, respectively. By
synthesising resolvent modes into physical coordinates,
$\Phi=V\widetilde{F}^{-1}\Phi_{\widetilde{F}}$ and
$\Psi=V\widetilde{F}^{-1}\Psi_{\widetilde{F}}$ represent the forcing and
response modes, respectively.

##  3 Results

###  3.1 Maximum gain

Figure 4: The gain spectra of the four test cases. The upper row is the
spectra for the $\Rey=100,000$ case and the bottom row is the spectra for
$\Rey=10,000$. The two columns represent the propulsion conditions; (a) drag
producing and (b) thrust producing. The annotations indicate the peak values
for each of the cases.

The dominant time scales of the system dynamics are identified by the forcing
frequency resulting in the maximum gain of the resolvent operator. Figure 4
shows the gain spectra of the resolvent operator for the four test cases where
the top and bottom rows differentiate between $\Rey=10,000$ and $100,000$, and
the two columns are the propulsion conditions; (a) drag producing and (b)
thrust producing. For ease of comparison, we have labelled the values of the
local maxima. The two peaks at the first and second harmonics remain largely
constant across all the cases, although there is minimal variation between the
thrust and the drag cases, with the thrust exhibiting slightly more
sensitivity at the first harmonic. The differences in stability between the
cases are present in the third harmonic; there is also a significant
difference between the thrust and the drag cases in $\Rey$ values. The
dynamics at $f^{*}=3$ are less sensitive to forcing at $\Rey=100,000$ than at
$\Rey=10,000$. The amplification of this mechanism is greater than in the
second harmonic and, for the $\Rey=10,000$ drag case, more than in the first
harmonic.

###  3.2 Resolvent modes

Figure 5: The response modes represented by the normalised tangential velocity
($u_{\xi}/||u_{\xi}||_{\infty}$), for the first three harmonics of the
resolvent operator. Each row represents the four distinct test cases and is
labelled as appropriate. The colour palette ranges from -1 (blue) to 1 (red).
The grey dashed line depicts the boundary-layer thickness.

The response modes of the resolvent operator are the mechanisms driving the
system, as depicted in figure 5. These modes correspond to the peaks shown in
figure 4, representing the first three harmonics of the resolvent operator.
The first harmonic, with the exception of the $\Rey=100,000$ drag production
case, is the most energetic mode. It exhibits coherent structures in both the
thrust and drag regimes, reinforcing the idea that kinematics are the dominant
factor. The drag cases differ from the thrust cases by the presence of large
fluctuations ($k_{\xi}\approx 3$) associated with flow reversal in the back
half of the foil, which disrupt the propulsive wave. In the thrust-producing
cases, these fluctuations are absent, and the energy amplification across the
back two-thirds of the foil is more uniform. The time averaged boundary-layer
thickness is indicated by the grey dashed line in figure 5. This was
identified via the shear-threshold method (Uzun & Malik, 2021), using a
threshold of 0.001. The unique thrust- and drag-producing structures in the
first harmonic scale with the boundary-layer thickness across the viscous
regimes implying a geometric self-similarity of the energetic structures.

The second harmonic predominantly exhibits waves that do not depend on $\eta$
or the boundary-layer thickness and increase in amplitude towards the tail,
matching the increasing motion amplitude envelope along the body length. The
streamwise wavenumber of these structures is twice the prescribed locomotive
wave for each case; for the drag-producing cases ($k_{\xi}\approx 3.5$) and
for the thrust-producing cases ($k_{\xi}\approx 2.5$). At the nose of the
foil, the $\Rey=100,000$ cases display fluctuations associated with shedding
that are absent in the $\Rey=10,000$ cases. Finally, the third harmonic shows
positive and negative fluctuations that intensify towards the tail in both
regimes, though they are more pronounced in the drag regime. The amplitude of
these fluctuations is significantly higher within the boundary-layer region.
This implies that the third harmonic is less correlated with body motion and
is driven by boundary-layer dynamics. Furthermore, the wavenumber of the
structures does not seem to correlate with any integer multiple of the
propulsive wave. Given the high amplification of the third harmonic, its
modification offers an avenue to control that does not depend on changing the
kinematics.

###  3.3 Roughness comparison

Figure 6: Study on the effect roughness has on the third harmonic of the
$\Rey=100,000$ thrust-producing case. (a) is a Q-criterion contour of
$10^{-4}$ coloured by $\omega_{z}$ for the flow over the thin foil with egg-
carton roughness. (b) shows the absolute value of the normalised vorticity
flux for the smooth foil forcing mode at the third harmonic. (c) is the
difference in vorticity flux between the smooth and the rough foil with the
boundary-layer thicknesses and a black, dash-dot line marking the height of
the bumps. (d) is the gain spectra for the rough foil resolvent operator.

As a final example, we demonstrate the ability of the new dynamic resolvent
analysis to study the impact of roughness on the amplification of the dominant
mechanisms in the flow. Specifically, we look at the change in amplification
of the mechanism at the third harmonic which we have identified as one
decorrelated to the kinematics (section 3.2). We make the link to the
vorticity flux as the potential control mechanism because vorticity flux-based
control has been shown to reduce turbulent fluctuations in a turbulent channel
by $40\%$ (Koumoutsakos, 1999).

In figure 6b, we show the absolute value of the vorticity flux
$\partial(\mathbf{\nabla}\times\mathbf{u})/\partial\eta$ in the forcing mode
at the third harmonic for the $\Rey=100,000$ thrust-producing case. The
forcing mode consists of significant areas of high vorticity flux near the
wall. We add an egg-carton roughness by adding a surface excursion such that

$y(x,z)=\begin{cases}h\sin\Big{(}\frac{2\pi
x}{\lambda}\Big{)}\cos\Big{(}\frac{2\pi z}{\lambda}\Big{)},&\text{for }y>0\\\
h\sin\Big{(}\frac{2\pi x}{\lambda}-\pi\Big{)}\cos\Big{(}\frac{2\pi
z}{\lambda}\Big{)},&\text{for }y<0\end{cases}$ (15)

where h is the bump amplitude, set at $0.002$, and $\lambda$ is the bump
wavelength, set at $1/128$ to ensure the roughness is small but still resolved
by the grid. We include 16 roughness wavelengths across the span of the foil.
Figure 6a shows the flow of the rough foil, with the Q-criterion coloured by
the spanwise vorticity. The flow remains similar to the smooth foil, although
the structures are offset by the roughness amplitude and break down before
being ejected from the tail.

The vorticity flux in the roughness sublayer is greatly reduced by the
roughness: Figure 6 c shows the difference between rough and smooth, and the
region within the roughness height is marked with a black dashed dotted line.
The reduction in flux directly opposes the mechanisms present in the forcing
mode of the smooth foil (figure 6a). In figure 6c, we also mark the boundary-
layer thicknesses for the smooth ($\delta$) and the rough ($\delta_{r}$) foil
which is shifted in the rough case. Near the nose, there is a region where the
vorticity flux is greatly increased; this is a side effect of the control
method. The reduction in vorticity flux corresponds to a significant reduction
in the amplification of the mechanism at the third harmonic, as shown by the
gain spectra in figure 6d. The amplification of the first and second harmonics
is also slightly reduced. The comparison shows that roughness can be used to
control the amplification of key mechanisms in the flow.

##  4 Conclusions

In this study, we applied data-driven resolvent analysis to explore the
boundary-layer dynamics of a foil undergoing undulatory propulsion, across
thrust- and drag-producing regimes, and at realistic Reynolds numbers. Using a
novel coordinate transformation, we were able to perform the data-driven
resolvent analysis technique. The analysis revealed that the kinematics of the
system are the dominant factor driving the physics. The thrust- and drag-
producing cases exhibit distinct structures that scale with the boundary-layer
thickness. The drag producing cases show large fluctuations associated with
flow reversal in the back half of the foil, disrupting the propulsive wave.
Finally, we provide a comparison to a foil with egg-carton roughness. The
roughness reduces the amplification of the mechanism at the third harmonic,
which is less correlated with body motion. This demonstrates the potential for
roughness to control the amplification of key mechanisms in the flow. By
extending resolvent analysis to deforming bodies with non-zero thickness, this
work opens new avenues for future research in understanding propulsion
mechanisms.

## References

  * Anderson et al. (2001) Anderson, E. J., McGillis, W. R. & Grosenbaugh, M. A. 2001 The boundary layer of swimming fish. J. Exp. Biol. 204 (1), 81–102.
  * Demo et al. (2018) Demo, N., Tezzele, M. & Rozza, G. 2018 Pydmd: Python dynamic mode decomposition. JOSS 3 (22), 530.
  * Di Santo et al. (2021) Di Santo, V., Goerig, E., Wainwright, D. K., Akanyeti, O., Liao, J. C., Castro-Santos, T. & Lauder, G. V. 2021 Convergence of undulatory swimming kinematics across a diversity of fishes. PNAS 118 (49).
  * Goza & Colonius (2018) Goza, A. & Colonius, T. 2018 Modal decomposition of fluid-structure interaction with application to flag flapping. J. Fluids Struct. 81, 728–737.
  * Herrmann et al. (2021) Herrmann, B., Baddoo, P. J., Semaan, R., Brunton, S. L. & McKeon, B. J. 2021 Data-driven resolvent analysis. J. Fluid Mech. 918, A10.
  * Koumoutsakos (1999) Koumoutsakos, P. 1999 Vorticity ﬂux control for a turbulent channel ﬂow. Phys. Fluids 11 (2), 248–250.
  * Lauber et al. (2022) Lauber, M., Weymouth, G. D. & Limbert, G. 2022 Immersed boundary simulations of flows driven by moving thin membranes. J. Comput. Phys. 457, 111076.
  * Lauber et al. (2023) Lauber, M., Weymouth, G. D. & Limbert, G. 2023 Rapid flapping and fibre-reinforced membrane wings are key to high-performance bat flight. J. R. Soc. Interface. 20 (208), 20230466.
  * Maertens et al. (2017) Maertens, A. P., Gao, A. & Triantafyllou, M. S. 2017 Optimal undulatory swimming for a single fish-like body and for a pair of interactingÂ swimmers. J. Fluid Mech. 813, 301–345.
  * Maertens & Weymouth (2015) Maertens, A. P. & Weymouth, G. D. 2015 Accurate Cartesian-grid simulations of near-body flows at intermediate Reynolds numbers. Comput. Methods Appl. Mech. Eng. 283, 106–129.
  * Massey et al. (2023) Massey, J. M. O., Ganapathisubramani, B. & Weymouth, G. D. 2023 A systematic investigation into the effect of roughness on self-propelled swimming plates. J. Fluid Mech. 971, A39.
  * McKeon & Sharma (2010) McKeon, B. J. & Sharma, A. S. 2010 A critical-layer framework for turbulent pipe flow. J. Fluid Mech. 658, 336–382.
  * Menon & Mittal (2020) Menon, K. & Mittal, R. 2020 Dynamic mode decomposition based analysis of flow over a sinusoidally pitching airfoil. J. Fluids Struct. 94.
  * Padovan et al. (2020) Padovan, A., Otto, S. E. & Rowley, C. W. 2020 Analysis of amplification mechanisms and cross-frequency interactions in nonlinear flows via the harmonic resolvent. J. Fluid Mech. 900, A14.
  * Rempfer & Fasel (1994) Rempfer, D. & Fasel, H. F. 1994 Evolution of three-dimensional coherent structures in a flat-plate boundary layer. J. Fluid Mech. 260, 351–375.
  * Saadat et al. (2017) Saadat, M., Fish, F. E., Domel, A. G., Di Santo, V., Lauder, G. V. & Haj-Hariri, H. 2017 On the rules for aquatic locomotion. Phys. Rev. Fluids 2 (8), 083102.
  * Schmid (2010) Schmid, P. J. 2010 Dynamic mode decomposition of numerical and experimental data. J. Fluid Mech. 656, 5–28.
  * Symon et al. (2019) Symon, S., Sipp, D. & McKeon, B. J. 2019 A tale of two airfoils: resolvent-based modelling of an oscillator versus an amplifier from an experimental mean. J. Fluid Mech. 881, 51–83.
  * Uzun & Malik (2021) Uzun, A. & Malik, M. R. 2021 Simulation of a turbulent flow subjected to favorable and adverse pressure gradients. Theor. Comput. Fluid Dyn. 35 (3), 293–329.
  * Williams et al. (2015) Williams, M. O., Kevrekidis, I. G. & Rowley, C. W. 2015 A data-driven approximation of the koopman operator: Extending dynamic mode decomposition. J. Nonlinear Sci. 25 (6), 1307–1346.
  * Zurman-Nasution et al. (2020) Zurman-Nasution, A. N., Ganapathisubramani, B. & Weymouth, G. D. 2020 Influence of three-dimensionality on propulsive flapping. J. Fluid Mech. 886.