# Fluid Antenna-Assisted Simultaneous Wireless Information and Power Transfer
Systems

Liaoshi Zhou, Junteng Yao, Tuo Wu, Ming Jin, Chau Yuen, _Fellow, IEEE_ ,
and Fumiyuki Adachi, _Life Fellow_ , _IEEE_ L. Zhou, J. Yao, and M. Jin are
with the faculty of Electrical Engineering and Computer Science, Ningbo
University, Ningbo 315211, China (E-mail: {2311100202, yaojunteng,
jinming}@nbu.edu.cn).T. Wu and C. Yuen are with the School of Electrical and
Electronic Engineering, Nanyang Technological University, 639798, Singapore
(E-mail: $\rm\\{tuo.wu,chau.yuen\\}@ntu.edu.sg$).F. Adachi is with the
International Research Institute of Disaster Science (IRIDeS), Tohoku
University, Sendai, Japan (E-mail: $\rm adachi@ecei.tohoku.ac.jp$).

###### Abstract

This paper examines a fluid antenna (FA)-assisted simultaneous wireless
information and power transfer (SWIPT) system. Unlike traditional SWIPT
systems with fixed-position antennas (FPAs), our FA-assisted system enables
dynamic reconfiguration of the radio propagation environment by adjusting the
positions of FAs. This capability enhances both energy harvesting and
communication performance. The system comprises a base station (BS) equipped
with multiple FAs that transmit signals to an energy receiver (ER) and an
information receiver (IR), both equipped with a single FA. Our objective is to
maximize the communication rate between the BS and the IR while satisfying the
harvested power requirement of the ER. This involves jointly optimizing the
BS’s transmit beamforming and the positions of all FAs. To address this
complex convex optimization problem, we employ an alternating optimization
(AO) approach, decomposing it into three sub-problems and solving them
iteratively using first and second-order Taylor expansions. Simulation results
validate the effectiveness of our proposed FA-assisted SWIPT system,
demonstrating significant performance improvements over traditional FPA-based
systems.

###### Index Terms:

Fluid antenna (FA), simultaneous wireless information and power transfer
(SWIPT), alternating optimization (AO).

##  I Introduction

Recently, fluid antenna (FA)-assisted wireless communication has garnered
substantial research attention, offering the ability to dynamically adjust
antenna positions and thus achieve higher data transmission rates for next-
generation communication systems [1, 2, 3, 4, 5]. Unlike traditional systems
that utilize fixed-position antennas (FPAs) and are limited in fully
exploiting spatial resources due to their immobility, FA-assisted systems
represent a new paradigm. These systems can actively reshape the wireless
propagation channel by maneuvering antennas within a designated area,
providing increased spatial degrees of freedom (DoFs) and achieving enhanced
communication rates [6, 7]. Consequently, FAs have been integrated into
various advanced wireless systems, such as FA-assisted multi-user systems [8],
FA-assisted over-the-air computation (AirComp) systems [9], and FA-assisted
non-orthogonal multiple access (NOMA) systems [10].

On the other hand, wireless communication networks are increasingly focused on
achieving energy-efficient communication. Information transmission enabled
simultaneous wireless information and power transfer (SWIPT) has emerged as a
promising energy harvesting (EH) technique to address the energy scarcity
problem in future energy-hungry wireless communication systems [11, 13]. The
fundamental concept of SWIPT involves a transmitter with a dedicated power
supply sending wireless signals to information receivers (IRs) and energy
receivers (ERs). While IRs decode the information, ERs harvest energy from
these signals. Given its advantages, SWIPT has been extensively explored
across various contexts, including NOMA systems [12], cognitive radio systems
[13], and secure systems [14].

Inspired by these benefits, integrating FA into SWIPT systems naturally
enhances communication rates and energy efficiency within wireless networks.
Unlike traditional SWIPT systems equipped with FPAs, where channel
characteristics between the IRs/ERs are static and performance gains are
limited to transmit beamforming, FAs offer the possibility to actively reshape
these channels. By strategically optimizing the positions of FAs at both the
transmitter and the receivers, improvements can be achieved in signal power,
harvested power, and beamforming flexibility. To the best of our knowledge, a
comprehensive optimization investigation into FA-assisted SWIPT systems has
not yet been conducted.

To address this research gap, we introduce a FA-assisted SWIPT system,
featuring a base station (BS) equipped with multiple FAs transmitting signals
to an IR with a single FA and an ER with a single FA. Our objective is to
maximize the IR’s communication rate, considering the harvested power
requirement of the ER and the power constraints of the BS. This involves
jointly optimizing the BS’s transmit beamforming, the positions of the BS’s
transmit FAs, and the positions of the IR’s and ER’s receive FAs. To tackle
this highly non-convex problem, we employ an alternating optimization (AO)
algorithm, utilizing first-order and second-order Taylor expansions to
precisely adjust the positions of all FAs. Simulation results confirm the
significant performance benefits of our proposed scheme over benchmarks,
particularly in terms of energy harvesting and communication efficacy.

Notations: $\mathrm{Tr}(\mathbf{A})$, $\mathbf{A}^{H}$, and $\mathbf{A}^{T}$
denote the spectral norm, trace, conjugate transpose, and transpose,
respectively; $\left\|\mathbf{a}\right\|_{2}$ denotes the 2-norm of vector
$\mathbf{a}$; $\mathbf{A}\succeq\mathbf{0}$ indicates that $\mathbf{A}$ is
positive semidefinite; $\Re\\{x\\}$ means the real part of $x$;
$\mathcal{CN}(0,\sigma^{2})$ denotes the distribution of a circularly
symmetric complex Gaussian variable with mean $0$ and variance $\sigma^{2}$.

##  II  System Model and Problem Formulation

Figure 1: The system model of FA-assisted SWIPT.

As shown in Fig. 1, we consider a FA-assisted SWIPT system, which contains a
BS, an IR, and an ER. We assume that the BS is equipped with $M(M\geq 2)$ FAs,
which can move freely within a certain range denoted as $\zeta_{t}$.The IR and
ER are both equipped with a single FA with corresponding certain movable
ranges, which are represented as $\zeta_{r,I}$ and $\zeta_{r,E}$,
respectively. We introduce a two-dimensional Cartesian coordinate system to
describe the locations of the FAs. Specifically, we use
$\mathbf{t}_{m}=[x_{m}^{t},y_{m}^{t}]^{T},m\in\mathcal{M}=\\{1,\cdots,M\\}$ to
denote the location of the $m$-th FA at the BS. Similarly, we use
$\mathbf{r}_{i}=[x^{r}_{i},y^{r}_{i}]^{T},i\in\\{I,E\\}$ to denote the
location of the FA at the IR or the ER.

During the downlink transmission, the signal received at the IR and the ER are
given by

$\mathbf{y}_{i}=\mathbf{h}_{i}\mathbf{w}x+n_{i},i\in\\{I,E\\}.$ (1)

Without loss of generality, we assume the data streams
$x\sim\mathcal{CN}(0,1)$ and $\mathbf{w}\in\mathbb{C}^{M\times 1}$ is the
transmit beamforming vector at the BS. $\mathbf{h}_{I}\in\mathbb{C}^{1\times
M}$ and $\mathbf{h}_{E}\in\mathbb{C}^{1\times M}$ are the channels from the BS
to the IR and the ER, respectively. $n_{I}\sim\mathcal{CN}(0,\sigma^{2}_{I})$
and $n_{E}\sim\mathcal{CN}(0,\sigma^{2}_{E})$ are the additive Gaussian noise
at the IR and the ER, respectively.

In this paper, we consider a far-field wireless channel model where the size
of the transmit/receive movable ranges are significantly smaller than the
signal propagation distance. Consequently, we assume a planar wavefront, and
the angles of arrival (AoA) and angles of departure (AoD) remain constant for
each channel path component [1, 2]. We denote the number of transmit paths for
the BS-IR and BS-ER links as $q_{t,I}$ and $q_{t,E}$, respectively. For the
transmit FAs of the BS, the signal propagation distance difference from the
origin $O^{t}$ on its $k$-th transmit path ($k\in{1,\cdots,q_{t,i}}$) is
defined as
$\rho_{t,i}^{k}(\mathbf{t})=x_{t}\sin\phi_{t,i}^{k}\cos\psi_{t,i}^{k}+y_{t}\cos\phi_{t,i}^{k}$,
where $i\in\\{I,E\\}$. Here, $\phi_{t,I}^{k}$ and $\psi_{t,I}^{k}$ denote the
elevation and azimuth AoDs of the $k$-th transmit path from the BS to the IR,
while $\phi_{t,E}^{k}$ and $\psi_{t,E}^{k}$ represent those from the BS to the
ER. Accordingly, the BS transmit field response vectors are then expressed as
follows [16]

$\bm{\xi}_{i}\mathbf{(t)}=\left[e^{j\frac{2\pi}{\lambda}\rho_{t,i}^{1}(\mathbf{t})},\cdots,e^{j\frac{2\pi}{\lambda}\rho_{t,i}^{q_{t,i}}(\mathbf{t})}\right]^{T},i\in\\{I,E\\},$
(2)

where $\lambda$ is the carrier wavelength. Thus, the transmission field
response matrices from the BS to the IR and the ER are given by

$\bm{\Xi}_{i}\mathbf{(\overline{t})}=\left[\bm{\xi}_{i}(\mathbf{t}_{1}),\bm{\xi}_{i}(\mathbf{t}_{2}),\cdots,\bm{\xi}_{i}(\mathbf{t}_{M})\right],i\in\\{I,E\\}.$
(3)

Similarly, let us denote ${q_{r,I}}$ and ${q_{r,E}}$ as the number of receive
paths for the BS-IR and BS-ER links, respectively. For the receive fluid
antennas (FAs) of the IR and ER, the signal propagation difference from their
respective reference points $O^{r}_{i}$ for the $s$-th receive path
$(s\in\\{1,\cdots,{q_{r,i}}\\})$ is given by
$\rho_{r,i}^{s}(\mathbf{r}_{i})=x_{r,i}\sin\phi_{r,i}^{s}\cos\psi_{r,i}^{s}+y_{r,i}\cos\phi_{r,i}^{s}$,
where $i\in\\{I,E\\}$. Here, $\phi_{r,i}^{s}$ and $\psi_{r,i}^{s}$ denote the
receive elevation and azimuth AoAs of the $s$-th path of the BS-IR and BS-ER
links. Consequently, the receive field response vectors at the IR and ER are
expressed as

$\displaystyle\mathbf{f}(\mathbf{r}_{i})=\left[e^{j\frac{2\pi}{\lambda}\rho_{r,i}^{1}(\mathbf{r}{{}_{i}})},\cdots,e^{j\frac{2\pi}{\lambda}\rho_{r,i}^{q_{r,i}}(\mathbf{r}{{}_{i}})}\right]^{T},i\in\\{I,E\\}.$
(4)

Moreover, let us define the path response matrix
$\bm{\Sigma}_{i}\in\mathbb{C}^{q_{r,i}\times q_{t}},i\in\\{I,E\\}$ as the
response coefficients of all paths from the origin $O^{t}$ of the BS to the
$O^{r}_{i},i\in\\{I,E\\}$ of the IR/ER. Then, the channel vectors of the BS-IR
and BS-ER links are given as

$\displaystyle\mathbf{h}_{I}=$
$\displaystyle\mathbf{f}^{H}(\mathbf{r}_{I})\bm{\Sigma}_{I}\bm{\Xi}_{I}\mathbf{(\overline{t})},$
(5) $\displaystyle\mathbf{h}_{E}=$
$\displaystyle\mathbf{f}^{H}(\mathbf{r}_{E})\bm{\Sigma}_{E}\bm{\Xi}_{E}\mathbf{(\overline{t})},$
(6)

respectively 111 The information of $\mathbf{h}_{I}$ and $\mathbf{h}_{E}$ can
be determined by using some channel estimation methods, e.g., compressed
sensing [15]..

Building upon (5), the communication rate at the IR is

$R=\log_{2}\left(1+\frac{\mathbf{h}_{I}\mathbf{W}\mathbf{h}^{H}_{I}}{\sigma^{2}_{I}}\right),$
(7)

where $\mathbf{W}=\mathbf{ww}^{H}$ is rank-one matrix.

For the ER, we consider a linear energy harvesting model, and the total
harvested power at the ER is given by [11], [14]

$Q=\tau\mathrm{Tr}\left(\mathbf{h}_{E}\mathbf{W}\mathbf{h}_{E}^{H}\right),$
(8)

where $0<\tau\leq 1$ is the energy harvesting efficiency.

In this paper, we aim to maximize the communication rate between the BS and IR
under the requirement of the harvested power of the ER by jointly optimizing
the transmit beamforming of the BS, the positions of FAs of the BS, the IR,
and the ER. The corresponding optimization problem can be expressed as

$\displaystyle\max\limits_{\mathbf{\overline{t}},\mathbf{r}_{i},\mathbf{W}\succeq\mathbf{0}}\quad$
$\displaystyle R$ (9a) $\displaystyle\mathrm{s.t.}\quad\ $
$\displaystyle\mathbf{\overline{t}}\in\zeta_{t},$ (9b)
$\displaystyle\mathbf{r}_{i}\in\zeta_{r,i},i\in\\{I,E\\},$ (9c)
$\displaystyle||\mathbf{t}_{m}-\mathbf{t}_{v}||_{2}\geq
D,~{}m,v\in\mathcal{M},~{}m\neq v,$ (9d)
$\displaystyle\mathrm{Tr}(\mathbf{W})\leq P_{\max},$ (9e) $\displaystyle
Q\geq\overline{Q},$ (9f) $\displaystyle\mathrm{rank}(\mathbf{W})=1,$ (9g)

where (9b) and (9c) are the movable range constraints of the FAs of the BS,
the IR, and the ER; (9d) represents the minimum distance constraint for the
transmitting FAs to prevent coupling; (9e) denotes the maximum transmit power
constraint at the BS. (9f) means the power received at the ER should be larger
than a certain value $\overline{Q}$, where $\overline{Q}$ is the harvested
power threshold at the ER. Since the optimization variables are coupled in
(9a), (9d) and (9f), Problem (9) is non-convex, which is difficult to obtain
an optimal solution.

##  III Proposed Algorithm

In this section, we decompose Problem (9) into three sub-problems about
transmit beamforming, transmit FAs’ locations, and receive FA’s locations, and
use the AO algorithm to iteratively optimize these sub-problems to obtain the
locally optimal solution.

###  III-A Optimization of Transmit Beamforming of BS

Given $\mathbf{\overline{t}}$, $\mathbf{r}_{I}$, and $\mathbf{r}_{E}$, Problem
(9) is non-convex due to the rank-one constraint (9g). By omitting the
constraint (9g), Problem (9) can be rewritten as

$\displaystyle\max\limits_{\mathbf{W}\succeq\mathbf{0}}\quad$ $\displaystyle
R$ (10a) $\displaystyle\mathrm{s.t.}\quad\ $
$\displaystyle\eqref{02e},\eqref{02f}.$ (10b)

It is readily to find that the objective function is concave with respect to
$\mathbf{W}$, and the constraints (9e) and (9f) are both linear. Thus, Problem
(10) is convex, which can be solved by the CVX toolbox [17]. After obtaining
the optimal $\mathbf{W}$, we can reconstruct the rank-one solution by using
Gaussian randomization [18].

###  III-B Optimization of Transmit FAs’ Locations of BS

Given $\mathbf{W}$, $\mathbf{r}_{I}$, and $\mathbf{r}_{E}$, our aim is to
optimize the position of the transmit FAs $\mathbf{\overline{t}}$. Since
$\log_{2}(1+x)$ is an increasing function with respect to $x$, maximizing $R$
is equivalent to maximizing
$\mathrm{Tr}\left(\mathbf{h}_{I}\mathbf{W}\mathbf{h}^{H}_{I}\right)$. Thus,
Problem (9) can be reformulated as

$\displaystyle\max\limits_{\mathbf{\overline{t}}}\quad$
$\displaystyle\mathrm{Tr}\left(\mathbf{h}_{I}\mathbf{W}\mathbf{h}^{H}_{I}\right)$
(11a) $\displaystyle\mathrm{s.t.}\quad$
$\displaystyle~{}~{}~{}\eqref{02b},\eqref{02d},\eqref{02f}.$ (11b)

We can see that the objective function
$\mathrm{Tr}\left(\mathbf{h}_{I}\mathbf{W}\mathbf{h}^{H}_{I}\right)$ and
$\mathrm{Tr}\left(\mathbf{h}_{E}\mathbf{W}\mathbf{h}_{E}^{H}\right)$ in (9f)
have the same form. Therefore, we can handle them by using the same approach.
Specifically, utilizing the property of the trace of a matrix, we can rewrite
them as

$\displaystyle\mathrm{Tr}\left(\mathbf{h}_{i}^{H}\mathbf{h}_{i}\mathbf{W}\right)=$
$\displaystyle\mathrm{Tr}\left(\sum_{m=1}^{M}{\gamma}_{i}(\mathbf{t}_{m}){\gamma}_{i}^{H}(\mathbf{t}_{m})\mathbf{W}\right)$
$\displaystyle=$
$\displaystyle\left(\eta_{i}+\gamma_{i}(\mathbf{t}_{m})\gamma_{i}^{H}(\mathbf{t}_{m})\right)\mathrm{Tr}\left(\mathbf{W}\right)$
$\displaystyle=$
$\displaystyle\left(\eta_{i}+\underbrace{\bm{\xi}_{i}^{H}(\mathbf{t}_{m})\bm{\omega}_{i}\bm{\xi}_{i}(\mathbf{t}_{m})}_{G_{i}(\mathbf{t}_{m})}\right)\mathrm{Tr}\left(\mathbf{W}\right),i\in\\{I,E\\},$
(12)

where

$\displaystyle\eta_{i}=$ $\displaystyle\sum_{j\neq
m}^{M}{\gamma}_{i}(\mathbf{t}_{j}){\gamma}_{i}^{H}(\mathbf{t}_{j}),i\in\\{I,E\\},$
(13) $\displaystyle\gamma_{i}(\mathbf{t}_{m})=$
$\displaystyle\mathbf{f}^{H}(\mathbf{r}_{i})\mathbf{\Sigma}_{i}\bm{\xi}_{i}(\mathbf{t}_{m}),i\in\\{I,E\\},$
(14) $\displaystyle\bm{\omega}_{i}=$
$\displaystyle\bm{\Sigma}_{i}^{H}\mathbf{f}(\mathbf{r}_{i})\mathbf{f}^{H}(\mathbf{r}_{i})\bm{\Sigma}_{i},i\in\\{I,E\\},$
(15)

Since $\\{\mathbf{t}_{j},j\neq m\\}_{j=1}^{M}$ is known, $\eta_{i}$ in (13) is
a scalar that is not related with $\mathbf{t}_{m}$. Therefore, we only need to
focus on $G_{i}(\mathbf{t}_{m})$. It can be observed that
$G_{i}(\mathbf{t}_{m})$ is a convex function of
$\bm{\xi}_{i}(\mathbf{t}_{m})$, we can obtain a lower bound for
$G_{i}(\mathbf{t}_{m})$ by employing a first-order Taylor expansion at point
$\widetilde{\mathbf{t}}_{m}$, which is given by

$\displaystyle
G_{i}(\mathbf{t}_{m};\widetilde{\mathbf{t}}_{m})=2\underbrace{\Re\left\\{\bm{\xi}_{i}^{H}(\widetilde{\mathbf{t}}_{m})\bm{\omega}_{i}\bm{\xi}_{i}(\mathbf{t}_{m})\right\\}}_{\widehat{G}_{i}(\mathbf{t}_{m})}-\mu_{i},i\in\\{I,E\\},$
(16)

where
$\mu_{i}=\bm{\xi}_{i}^{H}(\widetilde{\mathbf{t}}_{m})\bm{\omega}_{i}\bm{\xi}_{i}(\widetilde{\mathbf{t}}_{m})$.
Since $\widehat{G}_{i}(\mathbf{t}_{m})$ is neither convex nor concave with
respect to $\mathbf{t}_{m}$, we utilize a second-order Taylor expansion to
construct a surrogate function for
$\widehat{G}_{i}(\mathbf{t}_{m}),i\in\\{I,E\\}$ [2]. Then, we obtain the
following equation [1]

$\displaystyle\widehat{G}_{i}(\mathbf{t}_{m};\widetilde{\mathbf{t}}_{m})=$
$\displaystyle\widehat{G}_{i}(\widetilde{\mathbf{t}}_{m})+\nabla\widehat{G}_{i}(\widetilde{\mathbf{t}}_{m})^{T}(\mathbf{t}_{m}-\widetilde{\mathbf{t}}_{m})$
$\displaystyle-\frac{{\kappa}_{m}^{i}}{2}(\mathbf{t}_{m}-\widetilde{\mathbf{t}}_{m})^{T}(\mathbf{t}_{m}-\widetilde{\mathbf{t}}_{m}),$
(17)

where

$\displaystyle\nabla\widehat{G}_{i}(\widetilde{\mathbf{t}}_{m})=$
$\displaystyle\left[\frac{\partial\widehat{G}_{i}(\widetilde{\mathbf{t}}_{m})}{\partial
x^{t}_{m}},\frac{\partial\widehat{G}_{i}(\widetilde{\mathbf{t}}_{m})}{\partial
y^{t}_{m}}\right],$ (18) $\displaystyle\kappa_{m}^{i}=$
$\displaystyle\frac{8\pi^{2}}{\lambda^{2}}\sum_{k=1}^{q_{t,i}}\lvert\varrho^{k}_{i}\rvert,$
(19)
$\displaystyle\frac{\partial\widehat{G}_{i}(\widetilde{\mathbf{t}}_{m})}{\partial
x^{t}_{m}}=$
$\displaystyle-\frac{2\pi}{\lambda}\sum_{k=1}^{q_{t,i}}\lvert\varrho_{i}^{k}\rvert\sin\phi_{t,i}^{k}\cos\varphi_{t,i}^{k}\sin(\varsigma^{k}_{i}(\mathbf{\widetilde{t}}_{m})),$
(20)
$\displaystyle\frac{\partial\widehat{G}_{i}(\widetilde{\mathbf{t}}_{m})}{\partial
y^{t}_{m}}=$
$\displaystyle-\frac{2\pi}{\lambda}\sum_{k=1}^{q_{t,i}}\lvert\varrho_{i}^{k}\rvert\cos\phi_{t,i}^{k}\sin(\varsigma^{k}_{i}(\mathbf{\widetilde{t}}_{m})),$
(21) $\displaystyle\varsigma^{k}_{i}(\mathbf{\widetilde{t}}_{m})=$
$\displaystyle
2\pi\rho_{t,i}^{k}(\mathbf{\widetilde{t}}_{m})/\lambda-\angle\varrho_{i}^{k},$
(22)

$\lvert\varrho_{i}^{k}\rvert$ is amplitude for the $k$-th entry of
$\bm{\xi}_{i}^{H}(\widetilde{\mathbf{t}}_{m})\bm{\omega}_{i}$, and
$\angle\varrho_{i}^{k}$ is the phase.

For the left side of constraint (9d), we can relax it using a first-order
Taylor expansion at point $\widetilde{\mathbf{t}}_{m}$ as the following form

$\displaystyle\beta(\mathbf{t}_{m};\widetilde{\mathbf{t}}_{m})=\frac{1}{||\widetilde{\mathbf{t}}_{m}-\mathbf{t}_{v}||_{2}}(\widetilde{\mathbf{t}}_{m}-\mathbf{t}_{v})^{T}(\mathbf{t}_{m}-\mathbf{t}_{v}).$
(23)

By removing the terms that are not related with $\mathbf{t}_{m}$, we can
reformulate Problem (11) as

$\displaystyle\max\limits_{\mathbf{t}_{m}}$
$\displaystyle~{}~{}\widehat{G}_{I}(\mathbf{t}_{m};\mathbf{t}_{m}^{(q)})$
(24a) $\displaystyle\mathrm{s.t.}$
$\displaystyle~{}~{}\mathbf{\overline{t}}\in\zeta_{t},$ (24b)
$\displaystyle~{}~{}\beta(\mathbf{t}_{m};\mathbf{t}_{m}^{(q)})\geq{D},m,v\in\mathcal{M},~{}m\neq
v,$ (24c)
$\displaystyle~{}~{}\widehat{G}_{E}(\mathbf{t}_{m};\mathbf{t}_{m}^{(q)})\geq\frac{\overline{Q}}{\tau\mathrm{Tr}(\mathbf{W})}+\mu_{E}-\eta_{E},$
(24d)

where $\mathbf{t}_{m}^{(q)}$ is the optimal $\mathbf{t}_{m}$ of the $q$-th
iteration. It can be seen that the objective function (24a) is a concave
quadratic function of $\mathbf{t}_{m}$, constraints (24b) and (24c) are
linear, and (24d) is a convex constraint. Therefore, Problem (24) is convex.

###  III-C Optimization of Receive FA’s Locations of IR and ER

Given $\mathbf{W}$ and $\mathbf{\overline{t}}$, our aim is to optimize the
FA’s location at the IR and the FA’s location at the ER. Because the
optimization variables $\mathbf{r}_{I}$ and $\mathbf{r}_{E}$ are decoupled in
Problem (9), Problem (9) can be decomposed into two optimization problems

$\displaystyle\max\limits_{\mathbf{r}_{I}}\quad$
$\displaystyle\mathrm{Tr}\left(\mathbf{h}_{I}\mathbf{W}\mathbf{h}_{I}^{H}\right)$
(25a) $\displaystyle\mathrm{s.t.}\quad$
$\displaystyle\mathbf{r}_{I}\in\zeta_{r,I},$ (25b)

and

$\displaystyle\max\limits_{\mathbf{r}_{E}}\quad$
$\displaystyle\mathrm{Tr}\left(\mathbf{h}_{E}\mathbf{W}\mathbf{h}_{E}^{H}\right)$
(26a) $\displaystyle\mathrm{s.t.}\quad$
$\displaystyle\mathbf{r}_{E}\in\zeta_{r,E}.$ (26b)

Denote the eigenvalue decomposition (EVD) of $\mathbf{W}$ by
$\mathbf{W}=\mathbf{U}_{W}\mathbf{\Lambda}_{W}\mathbf{U}_{W}^{H}$, we can get
$\mathbf{Z}(\mathbf{r}_{i})=\mathbf{h}_{i}\mathbf{U}_{W}\mathbf{\Lambda}_{W}^{\frac{1}{2}}\in\mathbb{C}^{1\times
M}$. Furthermore, we let
$\mathbf{z}(\mathbf{r}_{i})=\mathbf{\Lambda}_{W}^{\frac{1}{2}}\mathbf{U}_{W}^{H}\bm{\xi}_{i}^{H}\mathbf{(\overline{t})}\bm{\Sigma}^{H}\mathbf{f}(\mathbf{r}_{i})\in\mathbb{C}^{M\times
1}$. Rewrite the objective functions in Problem (25) and Problem (26) as

$\displaystyle\mathrm{Tr}\left(\mathbf{h}_{i}\mathbf{W}\mathbf{h}_{i}^{H}\right)=$
$\displaystyle\mathrm{Tr}\left(\mathbf{Z}(\mathbf{r}_{i})\mathbf{Z}^{H}(\mathbf{r}_{i})\right)$
$\displaystyle=$
$\displaystyle\mathbf{z}^{H}(\mathbf{r}_{i})\mathbf{z}(\mathbf{r}_{i})$
$\displaystyle=$
$\displaystyle\mathbf{f}^{H}(\mathbf{r}_{i})\bm{\delta}_{i}\mathbf{f}(\mathbf{r}_{i}),i\in\\{I,E\\},$
(27)

where
$\bm{\delta}_{i}=\bm{\Sigma}_{i}\bm{\xi}_{i}\mathbf{(\overline{t})}\mathbf{U}_{W}\bm{\Lambda}_{W}^{\frac{1}{2}}\bm{\Lambda}_{W}^{\frac{1}{2}}\mathbf{U}_{W}^{H}\bm{\xi}_{i}^{H}\mathbf{(\overline{t})}\bm{\Sigma}_{i}^{H}$.
Therefore, we can rewrite Problem (25) and Problem (26) as follow

$\displaystyle\max\limits_{\mathbf{r}_{i}}\quad$
$\displaystyle\mathbf{f}^{H}(\mathbf{r}_{i})\bm{\delta}_{i}\mathbf{f}(\mathbf{r}_{i})$
(28a) $\displaystyle\mathrm{s.t.}\quad$
$\displaystyle\mathbf{r}_{i}\in\zeta_{r,i},i\in\\{I,E\\}.$ (28b)

It can be seen that the objective function in (28a) is similar to
$G_{i}(\mathbf{t}_{m})$ in (III-B). Therefore, we can also use the second-
order Taylor expansion method to solve (28a), hence we do not repeat it here.

The procedures for solving Problem (9) is summarized in Algorithm 1, where the
optimal $\mathbf{\overline{t}}$, $\mathbf{{r}}_{I}$, $\mathbf{{r}}_{E}$, and
$\mathbf{W}$ in the $l$th iteration are denoted as
$\mathbf{\overline{t}}^{(l)}$, $\mathbf{{r}}_{I}^{(l)}$,
$\mathbf{{r}}_{E}^{(l)}$, and $\mathbf{W}^{(l)}$, respectively.

Algorithm 1 The Proposed Alternating Optimization Algorithm

1: Initialize: $l=0$, $\mathbf{\overline{t}}^{(0)}$, $\mathbf{{r}}_{I}^{(0)}$,
and $\mathbf{{r}}_{E}^{(0)}$;

2: Repeat

3: $l:=l+1$;

4: Update $\mathbf{W}^{(l)}$ by solving Problem (10) ;

5: Repeat

6: $q:=q+1$;

7: for $m=1\rightarrow M$

8: Update $\mathbf{\overline{t}}_{m}^{(q)}$ by solving Problem (24) ;

9: end

10: Until: Convergence.

11: Update $\mathbf{\overline{t}}^{(l)}$;

12: Update $\mathbf{{r}}_{I}^{(l)}$ and $\mathbf{{r}}_{E}^{(l)}$ by solving
Problem (28);

13: Until: Convergence.

##  IV Numerical Result

In our simulation experiments, we assume that the energy harvesting efficiency
$\tau$ is set to $0.5$. The number of transmit and receive paths is set as
$q_{t,I}=q_{t,E}=q_{r,I}=q_{r,E}=3$, and the elevation and azimuth angles
$\phi_{t,i}^{k},\psi_{t,i}^{k},\phi_{r,i}^{s},\psi_{r,i}^{s},i\in\\{I,E\\}$
are all independent and identically distributed variables randomly distributed
in $[0,\pi]$. The minimum distance between two FAs is set to be $D=\lambda/2$
with $\lambda=1$ m [16], and the movable region for FAs at the BS is
$\zeta_{t}=\left[-{\frac{A}{2}},{\frac{A}{2}}\right]\times\left[-{\frac{A}{2}},{\frac{A}{2}}\right]$,
while the movable regions for the FA at the ER and the IR are both
$\zeta_{r,I}=\zeta_{r,E}=\left[-{\frac{A}{4}},{\frac{A}{4}}\right]\times\left[-{\frac{A}{4}},{\frac{A}{4}}\right]$,
where $A=4\lambda$ [1]. Furthermore, we use the geometric channel model to
describe the channels of the BS-IR and BS-ER links [2]. The path response
matrixs $\bm{\Sigma}_{I}$ and $\bm{\Sigma}_{E}$ are assumed to be diagonal
with $\bm{\Sigma}_{i}[1,1]\sim\mathcal{CN}({0},{\nu/(\nu+1)})$ and
$\bm{\Sigma}_{i}[p,p]\sim\mathcal{CN}({0},{1/(\nu+1)(q_{r,i}-1)})$ for
$p=2,3,\cdots,q_{r,i},i\in\\{I,E\\}$, where $\nu$ represents the ratio of the
average power of the line-of-sight (LoS) path to that of the non-line-of-sight
(NLoS) path and $\nu=1$ is assumed for the simulation [1]. We assume that the
noise powers of the IR and the ER are $\sigma_{I}^{2}=\sigma_{E}^{2}=1$, the
maximum transmit power to noise ratio is $P_{\max}/\sigma_{I}^{2}$ = 5 dB, and
the harvested power threshold to noise ratio is $\overline{Q}/\sigma_{E}^{2}$
= 0 dB.

Figure 2: Convergence of the proposed AO algorithm.

In Fig. 2, we demonstrate the convergence behavior of our proposed AO
algorithm for $M=2,3,4$. As observed from Fig. 2, the AO algorithm converges
after approximately 11 iterations. Additionally, the communication rate
increases with the number of FAs, denoted by $M$. This increase is attributed
to the greater DoFs provided by the FAs as $M$ increases, enhancing the
system’s capacity to optimize signal paths.

In Fig. 3, Fig. 4, and Fig. 5, we compare the proposed scheme with the
following benchmarks:

$\mathbf{transmit\ FA\ (TFA)}$: The BS is equipped with $M$ FAs, while the IR
and the ER are both equipped with a fixed single antenna.

$\mathbf{receive\ FA\ (RFA)}$: The BS is equipped with $M$ fixed antennas,
while the IR and the ER are both equipped with a single FA.

$\mathbf{FPA}$: The BS is equipped with $M$ fixed antennas, while the IR and
the ER are bothe equipped with a single fixed antenna.

Figure 3: $P_{\max}/\sigma^{2}_{I}$ versus $R$, where $M=4$.

In Fig. 3, we compare the communication rate $R$ of various schemes under
different $P_{\max}/\sigma^{2}_{I}$ ratios, where $M=4$. It is observed that
as the $P_{\max}/\sigma^{2}_{I}$ increases, the communication rate of all
schemes increases. Additionally, our proposed scheme consistently outperforms
the other schemes. Specifically, at $P_{\max}/\sigma^{2}_{I}=10$ dB, our
proposed scheme achieves improvements of $13.2\%$, $40.6\%$, and $73.1\%$ over
the “TFA” scheme, the “RFA” scheme, and the “FPA” scheme, respectively.

Figure 4: $A/\lambda$ versus $R$, where $M=4$.

In Fig. 4, we investigate the impact of the normalized FAs’ movable range
$A/\lambda$ on the communication rate $R$, where $M=4$. The results indicate
that, as $A/\lambda$ increases, the communication rates of all schemes except
for the “FPA” scheme increase. This stagnation in the “FPA” scheme occurs
because it employs fixed-position antennas (FPA). Notably, the communication
rates of all schemes except for the “FPA” scheme converge when $A/\lambda$
exceeds 3. This suggests that the communication rate $R$ reaches its maximum
within a limited movable range. Furthermore, our proposed scheme consistently
outperforms the other benchmark schemes, demonstrating the advantages of
utilizing FAs.

Figure 5: $\overline{Q}/\sigma^{2}_{E}$ versus $R$, where $M=4$.

In Fig. 5, we examine the impact of $\overline{Q}/\sigma^{2}_{E}$ on the
communication rate $R$, where $M=4$. Observations from Fig. 5 reveal that as
$\overline{Q}/\sigma^{2}_{E}$ increases, the communication rate $R$ decreases
for all schemes. This decline is attributable to the higher energy harvesting
requirement from the ER, which necessitates the BS to focus its transmission
towards the ER, thus compromising the achievable rate $R$. Moreover, our
proposed scheme demonstrates its capability to dynamically adjust the FAs’
locations at the BS, the IR, and the ER to mitigate this impact, consistently
achieving higher $R$ than other benchmark schemes under the same
$\overline{Q}/\sigma^{2}_{E}$ conditions.

##  V Conclusion

In this paper, we proposed a FA-assisted SWIPT system, where the BS, the IR,
and the ER are both equipped with FA. By jointly designing the transmit
beamforming of the BS and the positions of all FAs, we addressed the
communication rate of the IR maximization problem under the constraint of the
harvested energy of the ER. Due to the non-convex nature of the problem, we
employ the AO algorithm to iteratively optimize the transmit beamforming and
the positions of the FAs to obtain a sub-optimal solution. Numerical results
show that the proposed FA-assisted SWIPT system achieves a higher
communication rate compared to systems based on FPAs under the same harvested
energy requirement of the ER.

## References

  * [1] W. Ma, L. Zhu, and R. Zhang, “MIMO capacity characterization for movable antenna systems,” IEEE Trans. Wireless Commun., vol. 23, no. 4, pp. 3392–3407, Apr. 2023.
  * [2] L. Zhu et al., “Movable-antenna enhanced multiuser communication via antenna position optimization,” IEEE Trans. Wireless Commun., vol. 23, no. 7, pp. 7214–7229, Jul. 2024.
  * [3] J. Yao et al., “Proactive monitoring via jamming in fluid antenna systems,” IEEE Commun. Lett., vol. 28, no. 7, pp. 1698–1702, Jul. 2024.
  * [4] K. K. Wong, A. Shojaeifard, K.-F. Tong, and Y. Zhang, “Performance limits of fluid antenna systems,” IEEE Commun. Lett., vol. 24, no. 11, pp. 2469–2472, Nov. 2020.
  * [5] K. K. Wong et al., “Fluid antenna systems,” IEEE Trans. Wireless Commun., vol. 20, no. 3, pp. 1950–1962, Mar. 2021.
  * [6] L. Zhu and K. K. Wong, “Historical review of fluid antennas and movable antennas,” arXiv preprint arXiv:2401.02362v2, Jan. 2024.
  * [7] X. Lai et al., “On performance of fluid antenna system using maximum ratio combining,” IEEE Commun. Lett., vol. 28, no. 2, pp. 402–406, Feb. 2024.
  * [8] H. Qin, W. Chen, Z. Li, Q. Wu, N. Cheng, and F. Chen, “Antenna positioning and beamforming design for fluid antenna-assisted multi-user downlink communications,” IEEE Wireless Commun. Lett., vol. 13, no. 4, pp. 1073–1077, Apr. 2024.
  * [9] D. Zhang et al., “Fluid antenna array enhanced over-the-air computation,” IEEE Wireless Commun. Lett., vol. 13, no. 6, pp. 1541–1545, Jun. 2024.
  * [10] J. Zheng et al., “FAS-assisted NOMA short-packet communication systems,” IEEE Trans. Veh. Technol., early access, doi:10.1109/TVT.2024.3363115, 2024.
  * [11] C. Pan et al.,“Intelligent reflecting surface aided MIMO broadcasting for simultaneous wireless information and power transfer,” IEEE J. Sel. Areas Commun., vol. 38, no. 8, pp. 1719–1734, Aug. 2020.
  * [12] J. Ren et al., “RIS-assisted cooperative NOMA with SWIPT,” IEEE Wireless Commun. Lett., vol. 12, no. 3, pp. 446–450, Mar. 2023.
  * [13] C. Hu, Q. Li, Q. Zhang, and J. Qin, “Security optimization for an AF MIMO two-way relay-assisted cognitive radio nonorthogonal multiple access networks with SWIPT,” IEEE Trans. Inf. Forensics Security, vol. 17, pp. 1481–1496, 2022.
  * [14] Z. Zhu et al., “Robust beamforming design for STAR-RIS-aided secure SWIPT system with bounded CSI error,” IEEE Trans. Green Commun. Netw., vol. 8, no. 1, pp. 345–361, Mar. 2024.
  * [15] W. Ma et al., “Compressed sensing based channel estimation for movable antenna communications,” IEEE Commun. Lett., vol. 27, no. 10, pp. 2747–2751, Oct. 2023.
  * [16] Y. Ye et al., “Fluid antenna-assisted MIMO transmission exploiting statistical CSI,” IEEE Commun. Lett., vol. 28, no. 1, pp. 223–227, Jan. 2024.
  * [17] M. Grant and S. Boyd, “CVX: MATLAB software for disciplined convex programming,” [Online]. Available: http://cvxr.com/cvx.
  * [18] B. Wei, P. Zhang, and Q. Zhang, “Active reconfigurable intelligent surface-aided over-the-air computation networks,” _IEEE Wireless Commun. Lett._ vol. 13, no. 4, pp. 1148-1152, Apr. 2024.


#  Eclares: Energy-Aware Clarity-Driven Ergodic Search

Kaleb Ben Naveed1, Devansh Agrawal2, Christopher Vermillion3 and Dimitra
Panagou1,2 ∗The authors would like to acknowledge the support of the National
Science Foundation (NSF) under grant no. 2223845 and grant no.
1942907.1Department of Robotics, University of Michigan, Ann Arbor, MI, 48109
USA<EMAIL_ADDRESS>Department of Aerospace Engineering, University of
Michigan, Ann Arbor, MI, 48109 USA. 3 Department of Mechanical Engineering,
University of Michigan, Ann Arbor, MI, 48109 USA.

###### Abstract

Planning informative trajectories while considering the spatial distribution
of the information over the environment, as well as constraints such as the
robot’s limited battery capacity, makes the long-time horizon persistent
coverage problem complex. Ergodic search methods consider the spatial
distribution of environmental information while optimizing robot trajectories;
however, current methods lack the ability to construct the target information
spatial distribution for environments that vary stochastically across space
and time. Moreover, current coverage methods dealing with battery capacity
constraints either assume simple robot and battery models, or are
computationally expensive. To address these problems, we propose a framework
called Eclares, in which our contribution is two-fold. 1) First, we propose a
method to construct the target information spatial distribution for ergodic
trajectory optimization using clarity, an information measure bounded between
$[0,1]$. The clarity dynamics allows us to capture information decay due to
lack of measurements and to quantify the maximum attainable information in
stochastic spatiotemporal environments. 2) Second, instead of directly
tracking the ergodic trajectory, we introduce the energy-aware (eware) filter,
which iteratively validates the ergodic trajectory to ensure that the robot
has enough energy to return to the charging station when needed. The proposed
eware filter is applicable to nonlinear robot models and is computationally
lightweight. We demonstrate the working of the framework through a simulation
case study. [Code]aaaCode:
https://github.com/kalebbennaveed/Eclares.git[Video]bbbVideo:
https://youtu.be/1ZCgxlHitzk

##  I Introduction

Autonomous robots are widely used in tasks that involve data acquisition over
long time horizons, such as search-and-rescue [1, 2], characterization of
ocean currents such as the Gulf Stream to study cross-shelf exchanges [3, 4],
water body exploration [5, 6] and wildfire monitoring [7]. Coverage planning
for data acquisition poses diverse challenges, including considering the
spatial distribution of environmental information while optimizing search
trajectories [8, 9, 10]. This is important because some regions have higher
density of information as compared to others. Other challenges involve
persistently monitoring a spatiotemporal environment that varies across space
and time [9, 11, 12, 7], and ensuring task persistence by enabling the robot
to return to the charging station and recharge when/as needed during
persistent coverage tasks [13, 14, 15, 16]. In this paper, we consider these
problems and propose a solution framework Eclares. Moreover, we also consider
the case when a spatiotemporal environment has process uncertainty, i.e., we
consider a class of stochastic spatiotemporal environments. Probabilistic
methods can be used to compute estimates of the environmental variation.
However, in general the stochastic variations in the environment can result in
information decay without continuous monitoring.

Recent advances in ergodic search [8, 9, 17, 18, 10] consider the spatial
distribution of the information over the environment while generating coverage
trajectories. This is acheieved by constructing a _target information spatial
distribution_ (TISD) over the environment and then adjusting the time spent in
specific regions according to the TISD. However, these methods either assume
spatiostatic environment [8, 10], which only varies through space, or
spatiotemporal with known variation model [9, 19, 20]. Therefore, these
methods can not incorporate potential information decay, which will result
from the environment’s stochastic variation across space and time. To address
this problem, we use clarity [21], an information measure bounded between
$[0,1]$, to construct the TISD. The proposed method to construct the TISD
captures information decay due to a lack of measurements, and quantifies the
maximum attainable information in stochastic spatiotemporal environments.

In order to ensure persistent monitoring of a stochastic spatiotemporal
environment for a long time horizon, the robot’s limited battery capacity must
be taken into consideration so that it can return to a charging station when
needed. Most of the recent work on task persistence mainly either uses control
barrier function (CBF)-based methods [22, 13, 14, 15] or Hamilton-Jacobi
reachability based methods [23, 16, 21]. CBF-based methods are computationally
efficient; however, they only assume simple robot and battery models.
Alternatively, reachability-based methods can handle nonlinear systems and are
robust with respect to external disturbances, but they are mostly used
offline, suffer from curse of dimensionality, and are computationally
expensive. The proposed online filter eware inspired by [24] addresses these
problems. eware iteratively validates the ergodic trajectory based on the
robot’s current battery levels, ensuring that the robot is aware of its
remaining battery energy during the coverage task and returns to the charging
station before the energy is depleted. In summary, the main contributions of
the proposed framework Eclares include:

  * •

Clarity-driven method to generate TISD, Which incorporates potential
information decay and quantifies maximum attainable information in a
stochastic spatiotemporal environment;

  * •

eware filter which iteratively validates the generated ergodic trajectory to
ensure that the robot has enough energy to return to the charging station when
needed.

The paper is organized as follows: Section II introduce preliminaries, Section
III formulates the problem statement, Section IV present the proposed solution
framework Eclares and Section V discuss simulation setup and results.

##  II Preliminaries

In this section, we outline the notation used throughout the paper, provide an
overview of the ergodic search methods, and introduce clarity as an
information measure.

####  II-1 Notation

Let $\mathbb{N}=\\{0,1,2,...\\}$. Let $\mathbb{R}$, $\mathbb{R}_{\geq 0}$,
$\mathbb{R}_{>0}$ be the set of reals, non-negative reals, and positive reals
respectively. Let $\mathcal{S}^{n}_{++}$ denote set of symmetric positive-
definite in $\mathbb{R}^{n\times n}$. Let $\mathcal{N}(\mu,\Sigma)$ denote a
normal distribution with mean $\mu$ and covariance
$\Sigma\in\mathcal{S}^{n}_{++}$. The space of continuous functions
$f:\mathcal{A}\to\mathcal{B}$ is denoted $C(\mathcal{A},\mathcal{B})$. The
$Q\in\mathbb{S}_{++}^{n}$ norm of a vector $x\in\mathbb{R}^{n}$ is denoted
$\left\|x\right\|_{Q}=\sqrt{x^{T}Qx}.$

####  II-2 Dynamics Model

Consider the continuous-time system dynamics comprising the robot and battery
discharge dynamics:

$\dot{\chi}=\begin{bmatrix}\dot{x}\\\ \dot{e}\end{bmatrix}=f(\chi,u)\\\
=\begin{bmatrix}f_{r}(x,u)\\\ f_{e}(e,x,u)\end{bmatrix}$ (1)

where
$\chi=\begin{bmatrix}{x}^{T},&e\end{bmatrix}^{T}\in\mathcal{Z}\subset\mathbb{R}^{n+1}$
is the system state consisting of the robot’s state
$x\in\mathcal{X}\subset\mathbb{R}^{n}$ and robot’s State-of-Charge (SoC)
$e\in\mathbb{R}_{\geq 0}$. $u\in\mathcal{U}\subset\mathbb{R}^{m}$ is the
control input, $f:\mathcal{Z}\times\mathcal{U}\rightarrow\mathbb{R}^{n+1}$
defines the continuous-time system dynamics,
$f_{r}:\mathcal{X}\times\mathcal{U}\rightarrow\mathbb{R}^{n}$ define robot
dynamics and $f_{e}:\mathbb{R}_{\geq
0}\times\mathcal{X}\times\mathcal{U}\to\mathbb{R}$ define battery discharge
dynamics.

####  II-3 Ergodic Search

Ergodic search [8, 9] is a technique to generate trajectories
$x:[t_{0},T]\to\mathcal{X}$ that cover a rectangular domain
$\mathcal{P}=[0,L_{1}]\times\cdots[0,L_{s}]\subset\mathbb{R}^{s}$, matching a
specified _target information spatial distribution_ (TISD)
$\phi:\mathcal{P}\to\mathbb{R}$, where $\phi(p)$ is the density at
$p\in\mathcal{P}$. Moreover, the spatial distribution of the trajectory $x(t)$
is defined as

$c(x(t),p)=\frac{1}{T-t_{0}}\int_{t_{0}}^{T}\delta(p-\Psi(x(\tau)))d\tau$ (2)

where $\delta:\mathcal{P}\to\mathbb{R}$ is the Dirac delta function and
$\Psi:\mathcal{X}\to\mathcal{P}$ is a mapping such that $\Psi(x(\tau))$ is the
position of the robot at time $\tau\in[t_{0},T]$. In other words, given a
trajectory $x(t)$, $c(x(t),p)$ represents the fraction of time the robot
spends at a point $p\in\mathcal{P}$ over the interval $[t_{0},T]$. Then, the
_ergodicity_ of $x(t)$ w.r.t to a TISD $\phi$ is

$\displaystyle\Phi(x(t),\phi)=\left\|c-\phi\right\|_{H^{-(s+1)/2}}$ (3)

where $\left\|\cdot\right\|_{H^{-(s+1)/2}}$ is the Sobolev space norm defined
in [8], i.e., $\Phi$ is a function space norm measuring the difference between
the TISD $\phi$ and the spatial distribution of the trajectory $c$. In [9] a
Projection-based Trajectory Optimization (PTO) method is proposed to generate
ergodic trajectories. It optimizes (over the space of trajectories $x(t)\in
C([t_{0},T],\mathcal{X}),u(t)\in C([t_{0},T],\mathcal{U})$)

$\displaystyle\min_{x(t),u(t)}$
$\displaystyle\Bigl{[}\Phi(x(t),\phi)+J_{b}(x(t))+\int^{T}_{t_{0}}\left\|u(\tau)\right\|^{2}d\tau\Bigr{]}$
(4) s.t. $\displaystyle\dot{x}=f(x,u)$ $\displaystyle x_{0}=x(t_{0})$

where $J_{b}$ penalizes trajectories that leave the domain, and $x_{0}$ is the
robot’s state at time $t_{0}$. The algorithm uses Fourier decomposition of
$c,\phi$ to numerically evaluate $\Phi$ and its gradients with respect to
$x(t)$ and $u(t)$. This enables an efficient gradient descent algorithm to
optimize (4). In this paper, we do not modify the algorithm in [9]. Instead,
we focus on developing a strategy to systematically construct the TISD $\phi$
given the environment model and the robot’s sensing model, as discussed later.

####  II-4 Clarity

Consider a stochastic variable (quantity of interest) $h\in\mathbb{R}$
governed by the process and output (measurement) models:

$\displaystyle\dot{h}$ $\displaystyle=w(t),$ $\displaystyle
w(t)\sim\mathcal{N}(0,Q)$ (5a) $\displaystyle y$ $\displaystyle=C(x)h+v(t),$
$\displaystyle v(t)\sim\mathcal{N}(0,R)$ (5b)

where $Q\in\mathbb{R}_{\geq 0}$ is the known variance associated with the
process noise, $y\in\mathbb{R}$ is the measurement,
$C:\mathcal{X}\to\mathbb{R}$ is the mapping between robot state and sensor
statecccC is set to 1 when sensed by the robot and 0 otherwise., and
$R\in\mathbb{R}$ is the known variance of the measurement noise.

In [21], we introduced the notion of clarity $q$ of the random quantity $h$,
which lies between $[0,1]$ and is defined such that $q=0$ represents $h$ being
unknown, and $q=1$ corresponds to $h$ being completely known. The clarity
dynamics for the subsystem (5a), (5b) are given by [21]

$\dot{q}=\frac{C(x)^{2}}{R}(1-q)^{2}-Qq^{2}$ (6)

If $C(x)=C$ is constant, (6) admits a closed-form solution for the initial
condition $q(0)=q_{0}$:

$q(t;q_{0})=q_{\infty}\left(1+\frac{2\gamma_{1}}{\gamma_{2}+\gamma_{3}e^{2kQt}}\right)$
(7)

where $q_{\infty}=k/(k+1)$, $k=C/\sqrt{QR}$, $\gamma_{1}=q_{\infty}-q_{0}$,
$\gamma_{2}=\gamma_{1}(k-1)$, $\gamma_{3}=(k-1)q_{0}-k$.

Notice that as $t\to\infty$, $q(t;q_{0})\to q_{\infty}\leq 1$ monotonically.
Thus $q_{\infty}$ defines the maximum attainable clarity. Equation (7) can be
inverted to determine the time required to increase the clarity from $q_{0}$
to some $q_{1}$. This time is denoted $\Delta T:[0,1]^{2}\to\mathbb{R}_{\geq
0}$:

$\displaystyle\Delta T(q_{0},q_{1})=t\text{ s.t.
}q(t,q_{0})=q_{1}\quad\text{for }q_{1}\in[q_{0},q_{\infty})$ (8)

For $q_{1}<q_{0}$, we set $\Delta T(q_{0},q_{1})=0$ while $\Delta
T(q_{0},q_{1})$ is undefined for $q_{1}\geq q_{\infty}$.

##  III Problem Formulation

###  III-A Environment Specification

Figure 1: The high-level overview of the proposed framework. Blocks (a) and
(c) are the proposed modules of the framework, while block (b) is borrowed
from [9].

Consider the coverage space $\mathcal{P}$. We discretized the domain into a
set of $N_{c}$ cells each with size $V$.dddSize is width in 1D, area in 2D,
and volume in 3D. Let $m_{c}:[t_{0},\infty)\to R$ be the (time-varying)
quantity of interest at each cell $c\in\mathcal{C}=\\{1,...,N_{c}\\}$. We
model the quantities of interest as independent stochastic processes:

$\displaystyle\dot{m}_{c}$ $\displaystyle=w_{c}(t),$ $\displaystyle
w_{c}(t)\sim\mathcal{N}(0,Q_{c})$ (9a) $\displaystyle y_{c}$
$\displaystyle=C_{c}(x)m_{c}+v_{c}(t),$ $\displaystyle
v_{c}(t)\sim\mathcal{N}(0,R)$ (9b)

where $y_{c}\in\mathbb{R}$ is the output corresponding to cell $c$. $R$ is the
measurement noise, and $Q_{c}\in\mathbb{R}_{>0}$ is the process noise variance
at each cell $c$. Since $m_{c}$ varies spatially and temporally under process
noise $Q_{c}$ for each cell $c\in\mathcal{C}$, the environment becomes a
_stochastic spatiotemporal environment_.

###  III-B Problem Statement

Consider a robotic system performing persistent coverage of a stochastic
spatiotemporal environment (9), i.e., over a time interval $[0,\infty)$.
Assume the desired quality of information at each cell has been specified.
This is encoded using a _target clarity_ $\overline{q}_{c}<q_{\infty,c}$ for
each cell $c\in\mathcal{C}$. The target clarity can be different at each cell,
indicating a different desired quality of information at each cell, but must
be less than $q_{\infty,c}$, the maximum attainable clarity of the cell. If
$\overline{q}_{c}\geq q_{\infty,c}$ for a cell $c\in\mathcal{C}$, then the
robot would try to spend infinite amount of time at a cell $c$, which is
undesirable.

We use clarity as our information metric since it is particularly effective
for stochastic spatiotemporal environments:

  * •

The clarity decay rate in cell $c$, i.e. $-Q_{c}q_{c}^{2}$, is explicitly
dependent on the stochasticity of the environment $Q_{c}$ in (9). This allows
the information decay rate to be determined from the environment model, and
not set heuristically [21]. Furthermore, spatiostatic environments are a
special case: by setting $Q_{c}=0$, clarity cannot decay.

  * •

While taking measurements of cell $c$, clarity $q_{c}$ monotonically
approaches $q_{\infty,c}<1$ for $Q_{c},R>0$. This indicates that maximum
attainable information is upper bounded.

In this persistent task, the robot replans a trajectory every
$T_{H}\in\mathbb{R}_{>0}$ seconds, i.e., at times $\\{t_{0},t_{1},...\\}$ for
$t_{k}=kT_{H}$, $k\in\mathbb{N}$. At the $k$-th iteration, the objective is to
minimize the mean clarity deficit $q_{d}(t_{k}+T_{H})$, which is defined as

$q_{d}(t_{k}+T_{H})=\frac{1}{N_{c}}\sum_{c=1}^{N_{c}}\max(0,\overline{q}_{c}-q_{c}(t_{k}+T_{H}))$
(10)

where $q_{c}(t_{k}+T_{H})$ is the clarity at time $t_{k}+T_{H}$ of cell
$c\in\mathcal{C}$. However, in order to persistently monitor a stochastic
spatiotemporal environment over a long time horizon, the robot’s energy
constraints must be taken into consideration. Thus, the overall problem can be
posed as follows:

$\displaystyle\min_{\chi(t),u(t)}\quad$ $\displaystyle q_{d}(t_{k}+T_{H})$
(11a) s.t. $\displaystyle\chi(t_{k})=\chi_{k}$ (11b)
$\displaystyle\dot{\chi}=f(\chi,u)$ (11c)
$\displaystyle\dot{q}_{c}=g(x,q_{c})\ \forall c\in\mathcal{C}$ (11d)
$\displaystyle e(t)\geq{e_{min}}$ (11e)

where $q_{d}(t_{k}+T_{H})$ is the mean clarity deficit at the end of system
trajectory $\chi(t;t_{k},\chi_{k})$, $\forall t\in[t_{k},t_{k}+T_{H}]$ given
by (10), $g:\mathcal{X}\times[0,1]\rightarrow\mathbb{R}_{\geq 0}$ define the
clarity dynamics (6), and $e_{min}$ is the minimum energy level allowed for
the robot.

##  IV Energy Aware Clarity Driven Ergodic Search (Eclares)

In this section, we introduce Eclares as a solution to the optimization
problem (11). We first provide the solution overview and discuss its
motivation. Next, we provide details on the method to construct the TISD and
lastly, we introduce the eware filter.

###  IV-A Method Motivation & Overview

In order to solve problem (11), we take inspiration from ergodic search [8,
9]. As discussed in Section II-3, ergodic search can generate trajectories by
solving problem (4). If the construction of TISD $\phi$ is based on the
current clarity $q_{c}(t)$ and the target clarity $\overline{q}_{c}$ at each
cell, then ergodic search can be used to minimize the mean clarity deficit
(10). Therefore, in this work, we propose a method to construct $\phi$ using
clarity. However, the ergodic trajectory optimization in (4) does not consider
the energy constraint (11e). One approach would be to add the energy
constraint in (4), but the non-convexity of (4) implies that guaranteeing
convergence and feasibility is challenging. Therefore, we propose the
framework Eclares, shown in Figure 1, as an approximate solution to the
problem (11).eeeIt is also important to note $q_{d}(T)$ is not differentiable
so computing gradients for optimization problem (11) can be challenging.
However, the problem (4) is differentiable and can be approximately solved
using PTO [9].

Algorithm 1 The genTISD algorithm

1:function genTISD ($q_{c}$, $\overline{q}_{c}$, environment model (9))

2: for $c\in\\{1,...,N_{c}\\}$ do

3: $k\leftarrow 1/\sqrt{Q_{c}R}$

4: $q_{\infty}\leftarrow k/(k+1)$

5: $\overline{q}\leftarrow\min(\overline{q}_{c},q_{\infty,c}-\epsilon)$

6: $\phi_{c}\leftarrow\Delta T(\overline{q},q_{c})$

7: end for

8: $\phi_{c}\leftarrow\phi_{c}/({\sum_{c=1}^{N_{c}}\phi_{c}}),\quad\forall
c\in\\{1,...,N_{c}\\}$

9: return $\phi_{c}\ \forall c\in\\{1,...,N_{c}\\}$

10:end function

Algorithm 2 The eware algorithm

1:function eware($\tau_{j},x_{k}^{ergo},x_{j-1}^{com}$)

2: Solve back-to-base trajectory (12a)

3: Solve initial value problem (13a)

4: $\tau_{f}\leftarrow\tau_{j}+T_{N}+T_{B}$

5: if $e(t)\geq 0\ \forall t\in[\tau_{j},\tau_{f}]$ and $x(\tau_{f})=x_{c}$
then

6: $x_{j}^{com}\leftarrow x_{j}^{can}$

7: else

8: $x_{j}^{com}\leftarrow x_{j-1}^{com}$

9: end if

10: return $x_{j}^{com}$

11:end function

Our solution involves decoupling (11) into two sub-problems: (A) we design a
trajectory that maximizes information collection without the energy constraint
(referred to as the _ergodic trajectory_), and (B) we construct a trajectory
that tracks part of the ergodic trajectory but also reaches the charging
station before the robot’s energy depletes (referred to as the _committed
trajectory_). The low level controller of the robot always tracks the last
committed trajectory. In this way, we ensure that the robot persistently
explores the domain without violating the energy constraint.

We run these steps at different rates. The ergodic trajectory is replanned
every $T_{H}$ seconds, while the committed trajectory is updated every
$T_{E}<T_{H}$ seconds.fff$T_{E},T_{H}\in\mathbb{R}_{>0}$ are user-defined
parameters. That is,

  * •

at each time $t_{k}=kT_{H},k\in\mathbb{N}$,

    * –

recompute the TISD $\phi$ using genTISD.

    * –

recompute an ergodic trajectory using PTO [9].

  * •

at each time $\tau_{j}=jT_{E},j\in\mathbb{N}$,

    * –

update the committed trajectory using eware.

In this work, we adopt PTO directly from [9].

Before describing genTISD and eware, we establish notation for trajectories.
Let $x([t_{k},t_{k}+T_{H}];t_{k},x_{k})$ be the ergodic trajectory generated
at time $t_{k}$ starting at state $x_{k}$ and defined over a time horizon
$T_{H}$. For compactness, we refer to $x([t_{k},t_{k}+T_{H}];t_{k},x_{k})$ as
$x^{ergo}_{k}$. We use similar notation to describe other trajectories later
on.

###  IV-B Generate Target Spatial Distribution (genTISD)

The genTISD algorithm is described in Algorithm 1. Let $\phi_{c}$ denote the
target information density evaluated for cell $c$. At the $k$-th iteration
(i.e, at time $t_{k}=kT_{H}$), we set $\phi_{c}$ to be the time that the robot
would need to increase the clarity from $q_{c}(t_{k})$ to the target
$\overline{q}_{c}$ by observing cell $c$ (Lines 3-6). This is determined using
(8). The small positive constant $\epsilon>0$ in Line 5 ensures that target
clarity is always less than the maximum attainable clarity, i.e.,
$\overline{q}_{c}<q_{\infty,c}$. Finally, we normalize $\phi_{c}$ such that
the sum of $\sum_{c\in\mathcal{C}}\phi_{c}=1$ (Line 8). Once $\phi$ is
constructed, PTO method [9] is used to generate the ergodic trajectory.

Figure 2: This figure illustrates the generation of candidate trajectory at
time $\tau_{j}$, combining a small portion of the ergodic trajectory (shown in
red) and the b2b trajectory (shown in green). Figure 3: The still frames from
quadrotor simulation are shown where the quadrotor with a downward-facing
camera and a gimbal is persistently monitoring a stochastic spatiotemporal
environment. The TISD is shown as a point cloud in each frame and the colored
bar in each frame shows the battery level. Frame (a) shows the quadrotor
starting the coverage mission with full energy. Moving across frames (b) and
(c), the quadrotor energy depletes and returns to the charging station
(modeled as the origin of the world frame). After battery swapping, the
quadrotor resumes the coverage mission as shown in the frame (d). The same
trend follows in the subsequent frames. Moreover, it can be seen that as the
quadrotor explores the environment, the TISD becomes more sparse, which
reflects increase in clarity in the region. For full simulation video, please
see https://youtu.be/i2yD3A2b5nA

###  IV-C Energy-Aware (eware)

Inspired by [24], eware is a lightweight filter to ensure that the robot
reaches the charging station before running out of energy. To do so, it
constructs a candidate trajectory that follows the ergodic trajectory for a
short time horizon $T_{N}$ before attempting to return to the charging
station. If the candidate trajectory is valid (defined below), it becomes a
committed trajectory. The low-level tracking controller always tracks the last
committed trajectory. Algorithm 2 describes the eware algorithm.

The $j$-th iteration of eware starts at time $\tau_{j}=jT_{E}$. Let the robot
state at time $\tau_{j}$ be $x_{j}\in\mathcal{X}$ and the system system at
time $\tau_{j}$ be $\chi_{j}\in\mathcal{Z}$. Suppose the last ergodic
trajectory generated is $x_{k}^{ergo}$ (defined over the interval
$[t_{k},t_{k}+T_{H}]$).

(Line 2) We construct a back-to-base (b2b) trajectory $x_{j}^{b2b}$, defined
over $[\tau_{j}+T_{N},\tau_{j}+T_{N}+T_{B}]$, by solving the optimal control
problem:

$\displaystyle\min_{x(t),u(t)}\ $
$\displaystyle\int_{\tau_{j}+T_{N}}^{\tau_{j}+T_{N}+T_{B}}\left\|x(\tau)-x_{c}\right\|_{\mathbf{Q}}^{2}+\left\|u(\tau)\right\|_{\mathbf{R}}^{2}d\tau$
(12a) s.t. $\displaystyle x(\tau_{j}+T_{N})=x_{j}^{ergo}(\tau_{j}+T_{N}),$
(12b) $\displaystyle\dot{x}=f_{r}(x,u)$ (12c)

where $x_{c}\in\mathcal{X}$ is the state of the charging station,
$\mathbf{Q}\in\mathbb{S}_{++}^{n}$ weights state cost, and
$\mathbf{R}\in\mathbb{S}_{++}^{m}$ weights control cost.

(Line 3): Once b2b trajectory $x_{j}^{b2b}$ is generated, we numerically
construct the candidate trajectory by solving the initial value problem

$\displaystyle\dot{\chi}$ $\displaystyle=f(\chi,u(t)),$ (13a)
$\displaystyle\chi(\tau_{j})$ $\displaystyle=\chi_{j}$ (13b) $\displaystyle
u(t)$
$\displaystyle=\begin{cases}\pi(\chi,x_{k}^{ergo}(t)),&t\in[\tau_{j},\tau_{j}+T_{N})\\\
\pi(\chi,x_{j}^{b2b}(t)),&t\in[\tau_{j}+T_{N},\tau_{j}+T_{N}+T_{B}]\end{cases}$
(13c)

where $\pi:\mathcal{Z}\times\mathcal{X}\rightarrow\mathcal{U}$ is a control
policy to track the portion of the ergodic trajectory and the b2b trajectory.
Figure 2 illustrates the generation of a candidate trajectory.

(Lines 5-9): Once the candidate trajectory is constructed, we check whether
the robot can reach the charging station before running out of energy. If it
can, the candidate replaces the committed trajectory.

##  V Simulation Case Study

###  V-A Simulation Setup

In this section, we evaluate Eclares using a simulation case study. We
synthesize a hypothetical scenario where a quadrotor modeled with nonlinear
dynamics described in [25, Eq. (10)] is persistently monitoring a stochastic
spatiotemporal environment modeled according to the dynamics described in (9).
The coverage domain is 20$\times$20 meter in area with grid cells of dimension
0.20$\times$0.20. We use the PTO method [9] with double integrator dynamics to
solve the optimization problem (4) for the time horizon $T_{H}=30.0$s. To
generate b2b trajectory, we solve the problem (12) using model predictive
control (MPC) with reduced linear quadrotor dynamics proposed by [25]. Another
MPC problem with the same reduced linear quadrotor dynamics was posed to track
the ergodic trajectory for shorter time horizon $T_{N}=2.0$s. Then the
committed trajectory is generated using the eware algorithm. The committed
trajectories are generated in the receding-horizon manner at $0.50$ Hz (i.e.,
$T_{E}=2.0$s) and tracked at $20.0$ Hz with zero-order hold.gggFor more
details on algorithm implementation, trajectory generation, and on
hypothetical stochastic environment, please see
https://github.com/kalebbennaveed/Eclares.git. If one of the committed
trajectory returns quadrotor back to the charging station, it stays there to
simulate recharging. Note that no assumptions are made about the charging
method; it may either involve stationary charging or battery swapping. We use
the battery dynamics described in [14, Eq. (6)] as it directly incorporates
robot control input as opposed to the worst-case approximation used in other
works [13, 15]. Throughout the simulation, we use the Runge-Kutta 4th order
(RK4) integration scheme. Figure 3 shows still frames from the light UAV
simulator where a quadrotor with downward-facing camera is exploring a
stochastic spatiotemporal environment.

Figure 4: The clarity deficit Map for the spatiostatic environment at time T =
0 s, T = 210 s, and T = 780 s is shown comparing two methods, uniform
distribution, and proposed method, with the same clarity $q$ and target
clarity $\overline{q}$ at T = 0 s. Figure 5: (a) Plot of the mean clarity
deficit $(q_{d})$ against time is shown for the spatiostatic environment for
three methods. (b) The plot of the $q_{d}$ against time is shown for the
stochastic spatiotemporal environment for three methods. It can be seen from
both plots that the proposed method outperforms other methods: Ergodic search
with uniform target spatial distribution and classic lawnmower coverage
method. Moreover, it is important to note that in the stochastic
spatiotemporal environment, it is not possible to acquire zero clarity error
as information is always decaying due to term $-Qq^{2}$ in (6). (c) Plot of
the quadrotor’s distance to charging station $(\|x-x_{c}\|)$ against battery’s
State of Charge (SoC) is shown for one quadrotor flight, where $x_{c}$ is the
position of the charging station. In case of eware, quadrotor safely returns
to the charging station as energy depletes with almost 3 % SoC. However,
without eware, quadrotor unaware of its remaining energy keep on following the
nominal path, leading to a crash.

###  V-B Results and Discussion

####  V-B1 Performance Comparison between Spatiostatic and Stochastic
Spatiotemporal Environment

Figure 5(a) and Figure 5(b) show the evolution of mean clarity deficit $q_{d}$
for the spatiostatic and stochastic spatiotemporal environments respectively.
For the proposed method, in the case of the spatiostatic environment, the
$q_{d}$ approached almost zero at time T = 780s; however, for the stochastic
spatiotemporal environment, the $q_{d}$ reduced to a non-zero value. This is
expected due to the inherent stochasticity present in the environment as the
information is decaying constantly due to the term $-Qq^{2}$ in (6), and it is
not possible to achieve zero clarity deficit. However, the proposed method
does provide an ability to set realistic expectations on the maximum possible
information that can be collected in a stochastic environment by computing
$q_{\infty,c}$, $\forall c\in C$.

####  V-B2 Comparison to Baseline methods

We compare the performance of the proposed methodology genTISD to the uniform
TISD method used in most of the ergodic literature [8, 10]. The uniform TISD
assigns equal importance to all cells $c\in\mathcal{C}$ in the domain. As
evident from Figure 5(a) and Figure 5(b), the proposed method shows a faster
convergence rate as compared to the uniform method. This is explained by the
robot trajectory shown in Figure 4 as the robot spends more time in the region
with higher clarity deficit as compared to the uniform method. It is important
to note that the past ergodic literature does not consider the stochastic
nature of the environment in their construction of target distribution. We
also evaluate the performance of the proposed coverage method against the
classic lawnmower coverage method. In the lawnmower coverage method, the robot
follows the prescribed path which covers the whole domain. The lawnmower
coverage pattern results in a repetitive paths as the robot covers same areas
multiple times while leaving gaps in some areas. This causes clarity deficit
to stabilize to a non-zero value as shown in Figure 5(a). In case of using
lawnmower coverage method for stochastic spatiotemporal environment, the
clarity deficit eventually starts increasing as the information constantly
decays due to the term $-Qq^{2}$ in (6).

####  V-B3 Evaluation of eware Filter

Figure 5(c) compares the quadrotor flight with and without the eware filter.
Evidently, with the eware filter, the quadrotor safely returns to the charging
station. However, in the absence of eware filter, quadrotor unaware of its
remaining battery energy keeps on following the nominal ergodic trajectory,
leading to a crash. Moreover, we observed that although the minimum SoC
specified was $0\%$, the quadrotor returns to the charging station with almost
$3\pm 1$ % remaining SoC. This number can be brought down by increasing the
frequency of the candidate trajectory generation from $0.5$ Hz to a larger
number; however, this comes at the expense of overall computational time. We
also compared the compute time of proposed method eware with that of the CBF-
based methods used for robotic task persistence since both are online methods.
The eware filter, updated at 0.5 Hz, took average time of 30.2 ms to
numerically forward integrate the candidate trajectory through quadrotor and
battery’s nonlinear dynamics. Alternatively, the CBF-based method used in
[13], updated at 20 Hz, took average time of 0.0551 ms to solve the
optimization problem described in [13, Eq. (2)] with single integrator
dynamics. While considering the update frequency, the CBF-based method [13] is
almost 13 times faster than the proposed method; however, it assumes simple
robot and battery models. Conversely, the proposed eware filter achieved real-
time performance with nonlinear quadrotor and battery dynamics. This allows
for more flexibility with robot and battery models and thus increases
generalizability of the proposed method .

##  VI Conclusion

In this work, we introduced a framework Eclares, in which we proposed a method
to update target information spatial distribution for the ergodic trajectory
optimization, which considers the stochastic spatiotemporal environments. We
also integrated a lightweight filter called eware between the high-level
ergodic planner and low-level tracking controller, which iteratively ensures
that the robot has enough battery energy to continue the coverage task.
Through simulation case study, we demonstrated that the proposed framework
incorporates potential information decay and quantifies maximum attainable
information in a stochastic environment while ensuring robot’s safe return to
the charging station before the battery energy depletes. However, We note that
the decoupling used in Eclars can lead to a suboptimal solution, since the b2b
trajectory is not optimized for information collection, which we hope to
address in future work. Other future directions include considering
environmental disturbances and extension of the proposed method to multi-agent
exploration scenarios.

## References

  * [1] S. Mayer, L. Lischke, and P. W. Woźniak, “Drones for search and rescue,” in _1st International Workshop on Human-Drone Interaction_ , 2019.
  * [2] S. Waharte and N. Trigoni, “Supporting search and rescue operations with uavs,” in _2010 international conference on emerging security technologies_. IEEE, 2010, pp. 142–147.
  * [3] G. Gawarkiewicz, R. E. Todd, W. Zhang, J. Partida, A. Gangopadhyay, M.-U.-H. Monim, P. Fratantoni, A. M. Mercer, and M. Dent, “The changing nature of shelf-break exchange revealed by the ooi pioneer array,” _Oceanography_ , vol. 31, no. 1, pp. 60–70, 2018.
  * [4] R. E. Todd, “Export of middle atlantic bight shelf waters near cape hatteras from two years of underwater glider observations,” _Journal of Geophysical Research: Oceans_ , vol. 125, no. 4, p. e2019JC016006, 2020.
  * [5] P. Sujit, J. Sousa, and F. L. Pereira, “Uav and auvs coordination for ocean exploration,” in _Oceans 2009-Europe_. IEEE, 2009, pp. 1–7.
  * [6] S. Manjanna, A. Q. Li, R. N. Smith, I. Rekleitis, and G. Dudek, “Heterogeneous multi-robot system for exploration and strategic water sampling,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2018, pp. 4873–4880.
  * [7] K. D. Julian and M. J. Kochenderfer, “Distributed wildfire surveillance with autonomous aircraft using deep reinforcement learning,” _Journal of Guidance, Control, and Dynamics_ , vol. 42, no. 8, pp. 1768–1778, 2019.
  * [8] G. Mathew and I. Mezić, “Metrics for ergodicity and design of ergodic dynamics for multi-agent systems,” _Physica D: Nonlinear Phenomena_ , vol. 240, no. 4, pp. 432–442, 2011. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S016727891000285X
  * [9] L. Dressel and M. J. Kochenderfer, “Tutorial on the generation of ergodic trajectories with projection-based gradient descent,” _IET Cyper-Phys. Syst.: Theory & Appl._, vol. 4, no. 2, pp. 89–100, 2019. [Online]. Available: https://doi.org/10.1049/iet-cps.2018.5032
  * [10] D. Dong, H. Berger, and I. Abraham, “Time optimal ergodic search,” _arXiv preprint arXiv:2305.11643_ , 2023.
  * [11] B. Haydon, K. D. Mishra, P. Keyantuo, D. Panagou, F. Chow, S. Moura, and C. Vermillion, “Dynamic coverage meets regret: Unifying two control performance measures for mobile agents in spatiotemporally varying environments,” in _2021 60th IEEE Conference on Decision and Control (CDC)_ , 2021, pp. 521–526.
  * [12] Y. Sung, D. Dixit, and P. Tokekar, “Environmental hotspot identification in limited time with a UAV equipped with a downward-facing camera,” _CoRR_ , vol. abs/1909.08483, 2019. [Online]. Available: http://arxiv.org/abs/1909.08483
  * [13] G. Notomista, S. F. Ruf, and M. Egerstedt, “Persistification of robotic tasks using control barrier functions,” _IEEE Robotics and Automation Letters_ , vol. 3, no. 2, pp. 758–763, 2018.
  * [14] G. Notomista, C. Pacchierotti, and P. R. Giordano, “Online robot trajectory optimization for persistent environmental monitoring,” _IEEE Control Systems Letters_ , vol. 6, pp. 1472–1477, 2022.
  * [15] H. Fouad and G. Beltrame, “Energy autonomy for robot systems with constrained resources,” _IEEE Transactions on Robotics_ , vol. 38, no. 6, pp. 3675–3693, 2022.
  * [16] W. Bentz and D. Panagou, “Energy-aware persistent coverage and intruder interception in 3d dynamic environments,” in _2018 Annual American Control Conference (ACC)_ , 2018, pp. 4426–4433.
  * [17] I. Abraham, A. Prabhakar, and T. D. Murphey, “An ergodic measure for active learning from equilibrium,” _IEEE Transactions on Automation Science and Engineering_ , vol. 18, no. 3, pp. 917–931, 2021.
  * [18] H. Coffin, I. Abraham, G. Sartoretti, T. Dillstrom, and H. Choset, “Multi-agent dynamic ergodic search with low-information sensors,” in _2022 International Conference on Robotics and Automation (ICRA)_. IEEE, 2022, pp. 11 480–11 486.
  * [19] A. Rao, A. Breitfeld, A. Candela, B. Jensen, D. Wettergreen, and H. Choset, “Multi-objective ergodic search for dynamic information maps,” in _2023 IEEE International Conference on Robotics and Automation (ICRA)_ , 2023, pp. 10 197–10 204.
  * [20] A. C. Garza, “Bayesian models for science-driven robotic exploration,” Ph.D. dissertation, Carnegie Mellon University, Pittsburgh, PA, September 2021.
  * [21] D. R. Agrawal and D. Panagou, “Sensor-based planning and control for robotic systems: Introducing clarity and perceivability,” _IEEE Control Systems Letters_ , vol. 7, pp. 2623–2628, 2023.
  * [22] A. D. Ames, X. Xu, J. W. Grizzle, and P. Tabuada, “Control barrier function based quadratic programs for safety critical systems,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 8, pp. 3861–3876, 2016.
  * [23] S. Bansal, M. Chen, S. Herbert, and C. J. Tomlin, “Hamilton-jacobi reachability: A brief overview and recent advances,” in _2017 IEEE 56th Annual Conference on Decision and Control (CDC)_. IEEE, 2017, pp. 2242–2253.
  * [24] D. Agrawal, R. Chen, and D. Panagou, “gatekeeper: Safety verification and control for nonlinear systems in unknown and dynamic environments,” 2023.
  * [25] B. E. Jackson, K. Tracy, and Z. Manchester, “Planning with attitude,” _IEEE Robotics and Automation Letters_ , vol. 6, no. 3, pp. 5658–5664, 2021\.
  * [26] T. Lee, M. Leok, and N. H. McClamroch, “Geometric tracking control of a quadrotor uav on se (3),” in _49th IEEE conference on decision and control (CDC)_. IEEE, 2010, pp. 5420–5425.
  * [27] “Ncss for grids (grid as point dataset),” accessed: 2023-07-13. [Online]. Available: https://ncss.hycom.org/thredds/ncss/grid/GOMl0.04/expt_32.5/2019/dataset.html

###  -A System Dynamics Model

####  -A1 3D Quadrotor Dynamics

We modeled the quadrotor using the dynamics described in [25, equation (10)].
The dynamics have the form:

$\dot{x}=\begin{bmatrix}\dot{r}\\\ \dot{q}\\\ \dot{v}\\\
\dot{\omega}\end{bmatrix}=\begin{bmatrix}v\\\ (\frac{1}{2})L(q)H\omega\\\
(\frac{1}{m})^{W}F(x,u)\\\
J^{-1}(^{B}\tau(x,u)-\hat{\omega}J\omega)\end{bmatrix}$ (14)

where $x\in\mathbb{R}^{13}$ and $u\in\mathbb{R}^{4}$ are the state and control
vectors, $r\in\mathbb{R}^{3}$ is the position vector, $q\in\mathbb{H}$ is
attitude where $\mathbb{H}$ is space of quaternions, $v\in\mathbb{R}^{3}$ is
the linear velocity, $\omega\in\mathbb{R}^{3}$ is the angular velocity, $m$ is
the mass, $J\in\mathbb{R}^{3\times 3}$ is the inertia matrix
,${}^{W}F(x,u)\in\mathbb{R}^{3}$ are the forces in the world frame,
${}^{B}\tau(x,u)$ are the moments in the body frame, $H\in\mathbb{R}^{4\times
3}$ is the matrix used to convert the 3-dim vector into 4-dim quaternion with
zero scalar part,

$\ H=\begin{bmatrix}0.0&0.0&0.0\\\ 1.0&0.0&0.0\\\ 0.0&1.0&0.0\\\
0.0&0.0&1.0\end{bmatrix}$ (15)

$L(q)\in\mathbb{R}^{4\times 4}$ is the matrix used for left multiplication of
quaternion,

$\ L(q)=\begin{bmatrix}q_{s}&-q_{v}^{T}\\\ q_{v}&q_{s}I+\hat{q_{v}}\\\
\end{bmatrix}$ (16)

where $q_{s}\in\mathbb{R}$ is the scalar part of the quaternion,
$q_{v}\in\mathbb{R}^{3}$ is the vector part of the quaternion, and $\hat{x}$
is the skew-symmetric matrix operator

$\hat{x}=\begin{bmatrix}0&-x_{3}&x_{2}\\\ x_{3}&0&-x_{1}\\\ -x_{2}&x_{1}&0\\\
\end{bmatrix}$ (17)

####  -A2 Battery Dynamics

We modeled the battery dynamics using the dynamics described in [14, equation
(6)]. The dynamics have the form:

$\dot{e}=-\eta C^{-1}\alpha(\|u\|^{2})$ (18)

where $C$ is battery’s rated capacity, $\eta$ is its Columbic efficiency,
$\alpha$ is a class $\mathcal{K}$ function, and $\|u\|$ is the magnitude of
the control input.

###  -B Ergodic Trajectory Generation

We use the Projection-based Trajectory Optimization (PTO) proposed by [9] for
ergodic trajectory optimization. For implementation sake, discrete
approximation of system dynamics (1) is used. The optimization problem to
generate the $T_{H}$ horizon discrete trajectory $\tilde{x}_{1:N_{h}}$ is
defined as

$\displaystyle\min_{\tilde{x}_{1:N_{h}},\tilde{u}_{1:N_{h}-1}}$
$\displaystyle\\{\Phi(\tilde{x}_{1:N_{h}},\phi)+J_{b}(\tilde{x}_{1:N_{h}})+\frac{\Delta
t}{2}\sum_{k=1}^{N_{h}-1}\tilde{u}_{k}^{T}\mathbf{R}\tilde{u}_{k}\\}$ (19)
s.t. $\displaystyle\tilde{x}_{1}=\tilde{x}_{ic},$
$\displaystyle\tilde{x}_{k+1}=\tilde{x}_{k}+f(\tilde{x}_{k},\tilde{u}_{k})\Delta
t$

where $\mathbf{R}\in\mathbb{S}_{++}^{m}$ weights control cost,
$\tilde{x}_{ic}$ is the initial robot state, and $J_{b}(\tilde{x})$ is the
boundary penalty to keep the robot within the boundary of the rectangular
coverage space $\mathcal{P}$. The $J_{b}(\tilde{x})$ is given as

$J_{b}(\tilde{x}_{1:N_{h}})=c_{b}\sum^{N_{h}}_{k=1}\sum^{s}_{i=1}(\max(\tilde{x}_{k,i}-L_{i},0)^{2}+\min(\tilde{x}_{k,i},0)^{2})$
(20)

where $c_{b}>0$ weights the boundary penalty, $L_{i}$ is the $i^{th}$ length
of the rectangular domain, and $\tilde{x}_{k,i}$ is the $i^{th}$ element of
state $\tilde{x}_{k}$. We used the double integrator dynamics while solving
the problem (19). Comprehensive details on the PTO method can be found in [9].
We iteratively generate ergodic trajectories of time horizon $T_{H}=30$ sec
with $dt=0.2$ at a frequency of $\frac{1}{30}$ Hz. If the quadrotor returns to
the charging station while following the nominal trajectories due to energy
depletion, the new ergodic trajectory is generated based on the current
clarity level after recharging.

###  -C Candidate Trajectory Generation

Figure 6: This figure illustrates the process of generating candidate
trajectories iteratively. The green trajectories are the valid candidate
trajectories or the committed trajectories. These trajectories follow part of
the ergodic trajectory (shown in yellow) and then go to the charging station.
If the candidate trajectory cannot ensure the quadrotor’s safe return to the
charging station, it is considered invalid (shown in red) and the quadrotor
takes the last committed trajectory back to the charging station.

The candidate trajectories are constructed by using the method described in
Section IV-C. The candidate trajectories are generated at $0.5$ Hz and tracked
with zero-order hold at $20.0$ Hz. Essentially, the candidate trajectories are
generated by stitching the shorter $T_{N}$ horizon ergodic trajectory with the
$T_{B}$ horizon b2b trajectory. Here, we describe the method to generate the
$T_{N}$ horizon ergodic trajectory where $T_{N}<T_{H}$ and the $T_{B}$ horizon
b2b trajectory.

####  -C1 $T_{N}$ horizon Ergodic Trajectory

The shorter horizon ergodic state trajectory $\tilde{x}_{1:N_{n}}$ and the
control trajectory $\tilde{u}_{1:N_{n}-1}$ are generated using the convex
Model Predictive Controller (MPC) formulation given as

$\displaystyle\min_{\tilde{x}_{1:N_{n}},\tilde{u}_{1:N_{n}-1}}$
$\displaystyle\sum_{k=1}^{N_{n}-1}\Bigl{[}(\tilde{x}_{k}-\tilde{x}^{ref}_{k})^{T}\mathbf{Q}(\tilde{x}_{k}-\tilde{x}_{k}^{ref})$
(21) $\displaystyle+$
$\displaystyle(\tilde{u}_{k}-\tilde{u}_{k}^{ref})^{T}\mathbf{R}(\tilde{u}_{k}-\tilde{u}_{k}^{ref})\Bigr{]}$
s.t. $\displaystyle\tilde{x}_{1}=\tilde{x}_{ic},$
$\displaystyle\tilde{x}_{k+1}=f(\bar{x}_{k},\bar{u}_{k})+\Biggl{[}\frac{\partial
f}{\partial
x}\bigg{|}_{\bar{x}_{k},\bar{u}_{k}}\Biggr{]}(\tilde{x}_{k}-\bar{x}_{k})$
$\displaystyle\hskip 28.45274pt+\Biggl{[}\frac{\partial f}{\partial
u}\bigg{|}_{\bar{x}_{k},\bar{u}_{k}}\Biggr{]}(\tilde{u}_{k}-\bar{u}_{k})$

where $(\bar{x}_{k},\bar{u}_{k})$ is the $k^{th}$ equilibrium point,
$\mathbf{Q}\in S_{++}^{m}$ weights state cost, $\mathbf{R}\in S_{++}^{m}$
weights control effort, $\tilde{x}_{k}^{ref}$ is the $k^{th}$ reference state
from the reference $T_{H}$ horizon ergodic state trajectory
$\tilde{x}^{ref}_{1:N_{n}}$, $\tilde{u}_{k}^{ref}$ is the $k^{th}$ reference
control input from $T_{H}$ horizon ergodic control trajectory
$\tilde{u}^{ref}_{1:N_{n}}$, and $N_{n}=\frac{T_{N}}{dt}$. For the convex MPC
problem, we use the linearized version of the nonlinear quadrotor dynamics
described in (14). It is important to note that the linearized quadrotor
system with $A\in\mathbb{R}^{13\times 13}$ and $B\in\mathbb{R}^{13\times 4}$
is not controllable with rank 12. Therefore, we reduce the system to
$\bar{A}\in\mathbb{R}^{12\times 12}$ and $\bar{B}\in\mathbb{R}^{12\times 4}$
using the attitude Jacobin for error dynamics proposed by [25]. For more
details, please see [25]. The reference state trajectory and reference control
trajectory used in the MPC problem (21) are formed using the double integrator
ergodic trajectory with unit quaternion. Apart from this method, differential
flatness property can be used to track the double integrator trajectory using
the nonlinear geometric controller [26]. We use the values of $T_{N}=2.0$ sec
with sampling time $dt=0.05$.

####  -C2 back-to-base(b2b) Trajectory

The b2b robot state trajectory $\tilde{x}_{1:N_{b}}$ and control trajectory
$,\tilde{u}_{1:N_{b}-1}$ are generated using the convex MPC formulation given
as

$\displaystyle\min_{\tilde{x}_{1:N_{b}},\tilde{u}_{1:N_{b}-1}}$
$\displaystyle\sum_{k=1}^{N_{b}-1}\Bigl{[}(\tilde{x}_{k}-x_{c})^{T}\mathbf{Q}(\tilde{x}_{k}-x_{c})+(\tilde{u}_{k})^{T}\mathbf{R}(\tilde{u}_{k})\Bigr{]}$
(22) s.t. $\displaystyle\tilde{x}_{1}=\tilde{x}_{N_{h}},$
$\displaystyle\tilde{x}_{k+1}=f(\bar{x}_{k},\bar{u}_{k})+\Biggl{[}\frac{\partial
f}{\partial
x}\bigg{|}_{\bar{x}_{k},\bar{u}_{k}}\Biggr{]}(\tilde{x}_{k}-\bar{x}_{k})$
$\displaystyle\hskip 28.45274pt+\Biggl{[}\frac{\partial f}{\partial
u}\bigg{|}_{\bar{x}_{k},\bar{u}_{k}}\Biggr{]}(\tilde{u}_{k}-\bar{u}_{k})$

where $x_{c}$ is the state of the charging station and the initial state
$\tilde{x}_{1}$ is the last state of the $T_{N}$ horizon ergodic trajectory.
We use the same reduced linearized quadrotor dynamics to generate a b2b
trajectory. We use the values of $T_{B}=10.0$ sec with sampling time
$dt=0.05$.

Figure 6 illustrates the process of generating candidate trajectories
iteratively. The green trajectories are the valid candidate trajectories or
the committed trajectories. These trajectories follow part of the ergodic
trajectory (shown in yellow) and then go to the charging station. If the
candidate trajectory cannot ensure the quadrotor’s safe return to the charging
station, it is considered invalid (shown in red) and then the quadrotor takes
the last committed trajectory back to the charging station.

###  -D Construction of the hypothetical stochastic spatiotemporally varying
environment

We constructed a hypothetical stochastic spatiotemporally varying environment
using the Salinity data from the Gulf of Mexico [27]. Salinity, one of the
oceanographic quantities, refers to the concentration of dissolved salts in
water. In the oceanographic community, salinity evolution in a water body can
be used to study cross-shelf exchanges. Cross-shelf exchanges refer to the
movement of nutrients from the shallow waters to the open ocean. We extracted
the salinity mean values along with the standard deviation for a specific
region and then scaled down the environment spatially and temporally to form a
hypothetical environment.