# A Continuum Approach for Collaborative Task Processing in UAV MEC Networks

Lorson Blair, Carlos A. Varela, Stacy Patterson  Department of Computer
Science
Rensselaer Polytechnic Institute
Troy, NY, USA
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Unmanned aerial vehicles (UAVs) are becoming a viable platform for sensing and
estimation in a wide variety of applications including disaster response,
search and rescue, and security monitoring. These sensing UAVs have limited
battery and computational capabilities, and thus must offload their data so it
can be processed to provide actionable intelligence. We consider a compute
platform consisting of a limited number of highly-resourced UAVs that act as
mobile edge computing (MEC) servers to process the workload on premises. We
propose a novel distributed solution to the collaborative processing problem
that adaptively positions the MEC UAVs in response to the changing workload
that arises both from the sensing UAVs’ mobility and the task generation. Our
solution consists of two key building blocks: (1) an efficient workload
estimation process by which the UAVs estimate the task field—a continuous
approximation of the number of tasks to be processed at each location in the
airspace, and (2) a distributed optimization method by which the UAVs
partition the task field so as to maximize the system throughput. We evaluate
our proposed solution using realistic models of surveillance UAV mobility and
show that our method achieves up to 28% improvement in throughput over a non-
adaptive baseline approach.

###### Index Terms:

Mobile edge computing, task offloading, Voronoi partitioning.

##  I Introduction

In its 2021 climate change report, the United Nations’ International Panel on
Climate Change (IPCC) gave a grim forecast regarding the frequency, severity,
and devastating consequences of natural disasters—hurricanes, wild fires,
floods, etc.—as a result of climate change [1]. With more catastrophic
disasters, there is an urgent need for more sophisticated prediction models
and enhanced disaster management and search and rescue techniques. Due to
their flexible mobility, relatively cheap cost, and ease of deployment,
unmanned aerial vehicles (UAVs) are becoming increasingly popular for use in
these applications [2, 3]. UAVs can be quickly and efficiently deployed in
disaster areas to assess the level of damage, monitor the progression of
wildfires and floods, draw digital maps of the disaster area, and identify
rescue candidates [4].

These applications require the collection of large amounts of data, for
example, video, infrared imaging, or air quality data, which must be processed
to generate actionable intelligence. UAVs have limited physical space, which
is reserved for data collection hardware, as well as a limited battery
capacity that must be dedicated to this data collection. Thus, the UAVs, or
_mobile sensing agents_ (MSAs), must offload data processing tasks to external
resources. While a traditional solution is to offload this computation to the
cloud, the urgency of the situation and potentially compromised infrastructure
require offloading to nearby compute resources for immediate processing. A
promising approach is to utilize a compute platform consisting of additional
vehicles equipped with larger batteries and dedicated Mobile Edge Computing
(MEC) servers [5, 6, 7, 8, 9, 10]. These _mobile compute agents_ (MCAs) can be
UAVs with larger batteries and specialized computing hardware, or they can be
ground vehicles with their own power sources and servers. The MCAs have the
freedom and flexibility to move in close proximity to the MSAs and to change
locations in response to demand that changes both in location and quantity.

The successful deployment of such a platform requires addressing several
research challenges. First, both the locations and quantities of computing
demand are unpredictable and dynamic; however, this information is crucial for
determining the best placement of the MCAs. Second, even if the demand is
known, the problem of optimally adapting the locations of the MCAs is
intractable on its face due to the complexity of wireless communication and
the binary nature of associating tasks to MCAs [11, 12]. Third, unlike
offloading to the cloud, MCAs have limited processing capacity, and so the
assignment of tasks must also take these capacities into account to ensure
balanced task processing. Finally, in disaster scenarios, communication
infrastructure may be severely limited due to damage or the remoteness of the
location; thus, the MCAs must be able to autonomously coordinate to provide
the task processing services, reserving external communication for important
results.

Several early works proposed solutions for an MEC platform consisting of a
single UAV [5, 13, 14, 15, 16, 7]; however, a single UAV may not be sufficient
to meet the computing demands of a large number of MSAs. More recent works
have considered the deployment of multi-UAV MEC systems and address the
problem of optimizing the MEC UAV trajectories for performance measures like
latency and energy. Many of these works assume that the demand is stationary
[17, 11, 8, 12]. A few works support dynamic demand [18, 9, 10]. However, all
of these works propose centralized solutions, making them unsuitable for our
setting.

We propose a novel distributed solution to maximize processing throughput in a
multi-UAV MEC system. Our solution adaptively positions a limited number of
MCAs in response to the changing workload that arises both from the mobility
and task generation of the MSAs. We adopt a _continuum_ approach to model the
distribution of demand over the surveillance region. This unique approach
enables us to formulate a continuous optimization problem that admits a
distributed solution. Our solution consists of two key building blocks: (1) an
efficient demand estimation algorithm in which the MCAs collaborate to
generate an estimate of the demand as a continuous time-varying field, and (2)
a theoretically-verified two-phase distributed optimization algorithm by which
the MCAs optimize their locations to maximize the task transmission rate while
ensuring that the assignments of tasks is proportional to the compute
resources at each agent. We illustrate the effectiveness of our method through
realistic numerical simulations, which show that our method provides up to 28%
improvement processing in throughput over a non-adaptive baseline method.

The rest of the paper is organized as follows. In Section II, we describe the
system model and problem formulation. Our proposed method is presented in
Section III. Simulation results are presented in Section IV. We summarize
related work in Section V. Finally, we conclude in Section VI.

##  II System Model and Problem Formulation

Figure 1: System model.

As shown in Fig. 1, our system consists of a set of MSAs that are responsible
for surveilling a region of interest. The MSAs do not have the compute
resources or battery power for onboard processing, so they must offload the
computation tasks to more highly-resourced MCAs. The goal is to offload the
MSA tasks to the MCAs and process them as quickly as possible to inform system
operators and emergency personnel.

To achieve this goal, the MCAs should be positioned close to MSAs so that the
transmission rate is maximized. Since the transmission rate is inversely
proportional to the transmit time, maximizing the rate also minimizes the
energy expended to transmit each unit of data, which extends the surveillance
lifetime of the MSAs. We must also ensure that the amount of work assigned to
an MCA at any time does not exceed its computing capacity to avoid queuing
delays at the MCAs. Thus, the challenge is to dynamically update the positions
and task assignments of the MCAs in response to the changing demands of the
MSAs (both in terms of quantity and location), while respecting the capacities
of the MCAs.

In the remainder of this section, we formalize our system model and the
distributed task processing problem.

###  II-A System Model

There are $N$ MSAs that fly at a fixed mean sea-level (MSL) altitude of
$H_{s}$ within a 3-dimensional region $\Omega$. Using a fixed altitude, rather
than frequently ascending and descending, allows the MSAs to save energy [6,
14]. We denote the position of MSA $n$ at time $t$ by
${\mathbf{w}_{n}(t)=[x_{n}(t),y_{n}(t),H_{s}]^{T}}$.

The $M$ MCAs operate in the same region $\Omega$ at a fixed MSL altitude of
$H_{c}\neq H_{s}$ to avoid collisions with the MSAs. The MCAs update their
locations collaboratively to provide processing services to the MSAs. We
denote the position of MCA $m$ at time $t$ by
$\mathbf{u}_{m}(t)=[x_{m}(t),y_{m}(t),H_{c}]^{T}$. We drop the $t$ from
$\mathbf{w}_{n}(t)$ and $\mathbf{u}_{m}(t)$ when the context is clear. We
assume that the MCAs do not have control over, nor a priori knowledge of, the
MSA trajectories.

####  II-A1 Communication Model

We make the common assumptions that the MCAs and MSAs are equipped with omni-
directional antennas and that they communicate using an Orthogonal Frequency-
Division Multiple Access (OFDMA) transmission mechanism (as in [19, 20]). The
communication channel between MCA $m$ and MSA $n$ is dominated by the line-of-
sight link transmission with a distance dependent path-loss model. The maximum
transmission rate between MCA $m$ and MSA $n$ in bits per second (bps) at time
$t$ depends on the distance between them and is given by

$\displaystyle r(\mathbf{u}_{m},\mathbf{w}_{n})$
$\displaystyle=B\log_{2}\left(1+\frac{\beta
P}{\sigma^{2}\|\mathbf{u}_{m}-\mathbf{w}_{n}\|^{2}}\right)$ (1)

where $B$ is the bandwidth allocated to each MSA, $\beta$ is the channel gain
at a reference distance of 1m, $\sigma^{2}$ is the power of the white Gaussian
noise, and $P$ is the MSA’s transmit power.

A small portion of the OFDMA subcarriers are used as control channels for
lightweight communication among MSAs and MCAs. As is common, we assume that
the effects of this control channel communication on the maximum transmission
rate is negligible.

####  II-A2 Workload and Offloading Model

Each MSA $n$ stores a queue of its unprocessed tasks. Tasks are continuously
added to the queue as the MSA collects data and are removed from the queue
when the MSA offloads them to an MCA. Let $d_{n}(t)$ denote the quantity of
unprocessed tasks at MSA $n$ at time $t$. As in previous works (cf. [5, 13,
14, 16, 7, 18]), we assume that the tasks are arbitrarily divisible. At each
time $t$, each MSA $n$ is assigned to offload to a single MCA $m$. Each MCA
$m$ has a processing capacity $c_{m}$, which represents the amount of tasks it
can process per second (bps), i.e., its maximum throughput. The MSA offloads
to its assigned MCA up to the rate defined in (1), provided the MCA has
remaining capacity. If the MCA does not have remaining capacity, the MSA
stores the tasks until it is assigned to an MCA that has capacity to process
them. We note that multiple MSAs can be assigned to the same MCA
simultaneously. In this case, all MSAs offload in parallel until the MCA’s
capacity is saturated, at which point, the MSAs queue any remaining tasks.

###  II-B Problem Formulation

Our goal is to develop a distributed solution by which the MCAs
collaboratively maximize the system throughput by maximizing the total task
transmission rate, subject to the constraint that the amount of tasks assigned
to each MCA $m$ does not exceed its processing capacity $c_{m}$. To achieve
this, at each time $t$, the MCAs must identify both their optimal locations as
well as the assignment of MSAs. The MCAs can query the MSAs to learn the
current quantities of queued tasks $d_{n}(t)$ and the current MSA locations
$\mathbf{w}_{n}(t)$.

We define the optimization variables ${\tau(n,t)}$, $n=1\ldots N$; the
variable ${\tau(n,t)}$ gives the rate at which MSA $n$ transmits its data at
time $t$. We also define the set of binary assignment variables $a_{m,n}(t)$,
$n=1\ldots N,m=1\ldots M$, where $a_{m,n}(t)=1$ if MSA $n$ is assigned to MCA
$m$ at time $t$, and is 0 otherwise. Let $\Omega_{c}$ denote the two-
dimensional plane in $\Omega$ at MSL height $H_{c}$. We formalize our problem
as follows:

$\displaystyle\underset{{\begin{subarray}{c}\mathbf{u}_{m}\in\Omega_{c},m=1\ldots
M\\\ a_{m,n},m=1\ldots M,n=1\ldots
N\end{subarray}}}{\text{maximize}}~{}\sum_{m=1}^{M}\sum_{n=1}^{N}a_{m,n}(t){\tau(n,t)}$
(2) subject to $\displaystyle\sum_{n=1}^{N}a_{m,n}(t){\tau(n,t)}\leq
c_{m},~{}~{}{m=1\ldots M}$ (3) $\displaystyle{\tau(n,t)}\leq
d_{n}(t),~{}~{}{n=1\ldots N}$ (4) $\displaystyle a_{m,n}(t){\tau(n,t)}\leq
r(\mathbf{u}_{m},\mathbf{w}_{n}),~{}m=1\ldots M,n=1\ldots N$ (5)
$\displaystyle\sum_{m=1}^{M}a_{m,n}(t)=1,~{}~{}n=1\ldots N$ (6) $\displaystyle
a_{m,n}(t)\in\\{0,1\\},~{}~{}m=1\ldots M,n=1\ldots N.$ (7)

Here, the objective (2) is to maximize the total transmission rate. The
constraints (3) ensure that the total transmission rate to each MCA does not
exceed its processing capacity. The constraints (4) require that an MSA cannot
offload more data than it has in its queue, and the constraints (5) require
that the transmission rate between an MSA and an MCA does not exceed the
maximum transmission rate defined in (1). Finally, the constraints (6) and (7)
ensure that each MSA is assigned to exactly one MCA.

The above optimization problem presents multiple challenges. First the
trajectories and the workloads of the MSAs are not known a priori, nor can
they be controlled by the MCAs. Second, even if these trajectories and
locations were known, this optimization problem is a mixed-integer nonlinear
programming problem, which is NP-Hard, in general. Further, both the objective
and constraints are coupled across the MCA locations and assignments. Thus,
even if the problem were tractable, it is not straightforward to devise a
distributed solution for it.

To address these challenges, we model the workload as a continuous spatial
field, a _task field_ , rather than considering the discrete locations of the
MSAs. This approach allows us to devise a distributed solution by which the
MCAs predict the evolution of the task field and adaptively position
themselves to maximize the transmission rates while maintaining their
processing capacity constraints. We describe our continuum approach and
proposed solution in the next section.

##  III Continuum Approach

As a first step towards making the problem tractable, we develop a continuous
representation of the tasks held by the MSAs. We define a continuous task
field over the sensing airspace $\Omega_{s}$, the two-dimensional plane in
$\Omega$ at height $H_{s}$. This field is described by a continuous function
$\rho$ where ${\rho({\mathbf{x}},t)\in\mathbb{R}_{\geq 0}}$ quantifies the
amount of tasks, in bits, that needs to be offloaded for processing from
location $\mathbf{x}\in\Omega_{s}$ at time $t$.

Let $\mathbf{u}(t)=[\mathbf{u}_{1}(t)^{T},\ldots,\mathbf{u}_{M}(t)^{T}]^{T}$
denote the positions of the MCAs at time $t$. To maximize the total
transmission rate, the tasks $\rho(\mathbf{x},t)$ should be assigned to the
MCA $m$ with maximal transmission rate $r(\mathbf{u}_{m},\mathbf{x})$. We
assign a cost to offloading one bit from location $\mathbf{x}$ to
$\mathbf{y}$, defined by

$\displaystyle{\omega}\left(\mathbf{x},\mathbf{y}\right)$
$\displaystyle=\frac{1}{r(\mathbf{x},\mathbf{y})}$ (8)

and note that minimizing ${\omega}\left(\mathbf{x},\mathbf{y}\right)$ will
maximize $r(\mathbf{x},\mathbf{y})$.

We create a Voronoi partitioning of $\Omega_{s}$ based on the locations of the
MCAs. The Voronoi region of MCA $m$ is defined as

$\displaystyle\mathcal{V}_{m}(\mathbf{u}(t))=\\{\mathbf{x}\in\Omega_{s}~{}|~{}{\omega}\left(\mathbf{u}_{m}(t),\mathbf{x}\right)\leq{\omega}\left(\mathbf{u}_{k}(t),\mathbf{x}\right)$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}k=1\ldots M,k\neq
m\\}.$ (9)

Note that the Voronoi region of MCA $m$ depends on the positions of all MCAs.
We denote the Voronoi partitioning of the entire region by:

$\displaystyle\mathcal{V}(\mathbf{u}(t))=\\{\mathcal{V}_{1}(\mathbf{u}_{1}(t)),\ldots,\mathcal{V}_{M}(\mathbf{u}_{m}(t))\\}.$
(10)

For a given Voronoi partitioning, the amount of tasks assigned to MCA $m$ at
time $t$ is given by the volume of its Voronoi region:

$\displaystyle|\mathcal{V}_{m}(\mathbf{u}_{m}(t))|$
$\displaystyle=\int_{\mathcal{V}_{m}}\rho(\mathbf{x},t)~{}\mathrm{d}\mathbf{x}.$
(11)

If $|\mathcal{V}_{m}|\leq c_{m}$, then the amount of tasks assigned to MCA $m$
is less than or equal to its processing capacity. The total cost of the tasks
assigned to MCA $m$ is:

$\displaystyle
f_{m}(\mathbf{u}(t))=\int_{\mathcal{V}_{m}}\rho(\mathbf{x},t){\omega}\left(\mathbf{u}_{m}(t),\mathbf{x}\right)\mathrm{d}\mathbf{x}$
(12)

and the total system cost is:

$\displaystyle F(\mathbf{u})=\sum_{m=1}^{M}f_{m}(\mathbf{u}(t)).$ (13)

We now reformulate problem (2) – (7) using this continuum task field model and
the Voronoi partitioning.

$\displaystyle\underset{\mathbf{u}_{m}\in\Omega_{c},m=1\ldots
M}{\text{minimize}}~{}~{}$ $\displaystyle F(\mathbf{u}(t))$ (14) subject to
$\displaystyle|\mathcal{V}_{m}|\leq c_{m},~{}~{}~{}~{}m=1\ldots M.$ (15)

The solution to this problem gives the locations of the MCAs as well as the
assignment of tasks to MCAs. By minimizing $F$, we maximize the total
transmission rate, while the constraints (15) ensure the capacity constraints
of the MCAs are respected. Provided the Voronoi region volumes do not violate
the capacity constraints, the tasks can all be offloaded and processed at
their respective maximal transmission rates.

The objective (14) and the constraints (15) are non-convex, and as far as we
are aware, there is no single distributed algorithm to solve such a problem.
Thus, we decompose the problem into tractable components that each admit a
distributed solution. We then combine these components to form our full
solution. We describe our proposed solution and each of its components below.

###  III-A Proposed Solution

We propose a distributed solution to problem (14) - (15).

Our solution is divided into four phases:

  1. 1.

Task field estimation. In this phase, the MCAs collect information about the
current state of the task queues of the MSAs and generate an estimate of
$\rho$.

  2. 2.

Transmission rate maximization. Next, the MCAs use a distributed algorithm to
determine the positions that minimize $F(\mathbf{u})$ over this estimated
field $\rho$.

  3. 3.

Capacity balancing. Then, the MCAs use a second distributed algorithm to
collaboratively update these positions so that the resulting Voronoi region
volumes are proportional to their processing capacities.

  4. 4.

Transmission and processing. The MCAs move to the new positions identified in
Phase 3, and the MSAs are assigned to MCAs using the Voronoi partitioning
defined in (10). The transmission rates are determined in a round robin
fashion. MCA $m$ considers each MSA $n$ one at a time and sets ${\tau(n,t)}$
to the maximum of $r(\mathbf{u}_{m},\mathbf{w}_{n})$ and its remaining
capacity.

Rather than continuously executing the first three phases, the MCAs execute
them every $\Delta$ seconds. A smaller value of $\Delta$ allows the MCAs to
more accurately estimate $\rho$ and adapt their positions as it changes over
time. This costs the MCAs energy to update their positions, which decreases
the energy available for task processing. A larger value of $\Delta$ leads to
less accurate estimates of $\rho$ and potentially lower transmission rates,
but it saves in the motion energy. We explore different values of $\Delta$ in
the experiments.

For the first $\Delta$ seconds, the MCAs hover in their initial positions.
After $\Delta$ seconds, they query the MSAs, which respond with the number of
tasks they each generated in the first $\Delta$ seconds. The MCAs use these
values as $d_{n}(t)$ in phase 1. The MCAs execute phases 2 and 3 to determine
their next locations, and then update their positions accordingly by flying
along straight paths to their destinations. They hover in these positions
until $2\Delta$ seconds total have elapsed, and the entire process is
repeated. Tasks are continuously assigned and processed as described in phase
4 using the Voronoi partitioning corresponding to the current MCA positions.

Next, we present the details of the first three phases.

###  III-B Task Field Estimation

We consider $\rho(x)$ as a spatial Gaussian process,
$GP(\mu,{\mathcal{K}}(\mathbf{x},\mathbf{y}))$ for
$\mathbf{x},\mathbf{y}\in\Omega_{s}$. The constant $\mu$ is the mean of the
Gaussian process, and ${\mathcal{K}}$ is the kernel function that gives the
spatial correlation of $\rho$. We use the squared exponential Gaussian kernel

$\displaystyle{\mathcal{K}}(\mathbf{x},\mathbf{y})=\varphi^{2}e^{\big{(}{-\frac{{\|\mathbf{x}-\mathbf{y}\|}^{2}}{L^{2}}}\big{)}}\hfill$
(16)

where $\varphi^{2}$ is the variance of the Gaussian process and $A_{1}$ is the
length scale. Gaussian processes are used in many applications to model
spatial processes [21].

At the beginning of each window, at time $t_{0}$, each MCA collects the value
$d_{n}(t_{0})$ and location $\mathbf{w}_{n}(t_{0})$ from each MSA in its
Voronoi region. The MCAs share these values with each other over the control
channels. Then, independently, each MCA uses the observed values of
$d_{n}(t_{0})$ and $\mathbf{w}_{n}(t_{0})$, ${n=1\ldots N}$, to construct an
estimate of $\rho$ using Gaussian process regression, where $\mu$,
$\varphi^{2}$, and $L$ are estimated from the observed values.

###  III-C Transmission Rate Maximization Algorithm

In this phase, the MCAs collaborate to solve the unconstrained optimization
problem:

$\displaystyle\underset{\mathbf{u}_{m}\in\Omega_{c},m=1,\ldots,M}{\text{minimize}}~{}~{}$
$\displaystyle F(\mathbf{u})=\sum_{m=1}^{M}f_{m}(\mathbf{u}).$ (17)

By doing so, the MCAs identify the positions and Voronoi partitioning that
maximizes the total _maximal_ transmission rate. However, this partitioning
does not necessarily respect the capacities of the MCAs.

Each MCA communicates with nearby MCAs to solve this problem. We say that two
MCAs are _neighbors_ for window $b$ if their Voronoi regions are adjacent in
window $b-1$. We define the communication graph for window $b$ as ${{\cal
G}_{b}=(V,E_{b})}$, where $V$ is the set of MCAs, and $(m,k)\in E_{b}$ if and
only if MCAs $m$ and $k$ are neighbors in window $b-1$. We denote the set of
neighbors of MCA $m$ in window $b$ by ${\cal N}_{b}(m)$.

Each MCA $m$ stores an estimate of the positions of all $M$ MCAs. We let
$\mathbf{u}^{m}_{k}$ denote MCA $m$’s estimate of the position of MCA $k$, and
$\mathbf{u}^{m}$ denotes the concatenation of all the estimates held by MCA
$m$, i.e.,

$\mathbf{u}^{m}=[(\mathbf{u}^{m}_{1})^{T},\ldots,(\mathbf{u}^{m}_{M})^{T}]^{T}.$

We can then rewrite (17) in the following decomposable form:

$\displaystyle\underset{\mathbf{u}^{1},\ldots,\mathbf{u}^{m}}{\text{minimize}}~{}$
$\displaystyle\hat{F}([(\mathbf{u}^{1})^{T},\ldots,(\mathbf{u}^{M})^{T}]^{T})=\sum_{m=1}^{M}f_{m}(\mathbf{u}^{m})$
(18) subject to $\displaystyle\mathbf{u}^{m}=\mathbf{u}^{k},~{}\text{for all
}(m,k)\in E_{b}.$ (19)

Algorithm 1 Transmission rate maximization algorithm.

1:Initialize: $\mathbf{u}^{m}(0)$ to current MCA positions

2:for $\ell=0,1,2,\ldots,T_{1}$ do

3: for each MCA $m$ do

4:
$\mathbf{z}^{m}(\ell)\leftarrow\text{Proj}_{\Omega_{c}}\left(\mathbf{u}^{m}(\ell)-\eta_{\ell}\nabla
f_{m}(\mathbf{u}^{m}(\ell))\right)$

5:
$\mathbf{u}^{m}(\ell+1)\leftarrow\xi\sum_{j\in\mathcal{N}_{b}(m)}\mathbf{z}^{j}(\ell)$

6:
$~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+(1-\xi|{\cal
N}_{b}(m)|)\mathbf{z}^{m}(\ell)$

7: end for

8:end for

The MCAs use a distributed projected gradient algorithm to solve this problem
[22, 23]. We now describe the details of this algorithm. The pseudocode is
shown in Alg. 1. The algorithm executes in discrete iterations
$\ell=0,1,2,\ldots,T_{1}$. Each MCA initializes its estimate $\mathbf{u}^{m}$
to be the current MCA positions, i.e., their positions at time $t_{0}$. In
each iteration $\ell$, each MCA $m$ first computes the derivative of $f_{m}$
with respect to its estimate $\mathbf{u}^{m}$. It then performs a gradient
step on this estimate with step size $\eta_{\ell}$ and projects the result of
this gradient step onto $\Omega_{c}$, where
$\text{Proj}_{\Omega_{c}}\left(\mathbf{x}\right)$ denotes the Euclidean
projection onto the region $\Omega_{c}$. To complete the iteration, the MCA
updates its estimate by taking a weighted average of its own estimate and
those of its neighbors; MCA $m$ applies weight $\xi$ to each of its neighbor’s
estimates and weight $(1-\xi|{\cal N}_{b}(m)|)$ to its own estimate.

The MCAs repeat this iteration until a desired convergence is achieved
(typically within a hundred iterations in our experiments). Intuitively,
through the projected gradient steps of Alg. 1, each MCA updates its estimates
of the positions to maximize the transmission rate in its own Voronoi region,
while the averaging steps drive the estimates to agreement.

We now summarize the convergence behavior of Alg. 1. Define the set of
stationary points on $(\Omega_{c})^{M^{2}}$ as
${\mathcal{L}=\\{\mathbf{x}\in(\Omega_{c})^{M^{2}}~{}|~{}\nabla\hat{F}(\mathbf{x})\in\mathcal{C}(\mathbf{x})}\\}$,
where $\mathcal{C}(\mathbf{x})$ is the normal cone, i.e.,
${\mathcal{C}(\mathbf{x})=\\{\mathbf{v}\in\mathbb{R}^{2M}~{}|~{}\forall\mathbf{u}^{\prime}\in(\Omega_{c})^{M^{2}},\mathbf{v}^{T}(\mathbf{x}-\mathbf{u}^{\prime})^{\prime}\geq
0\\}}$. With this definition, we present the following theorem.

###### Theorem 1

Let $\Omega_{c}$ be a nonempty compact convex set, and assume that
$\hat{F}(\mathcal{L})$ has an empty interior and that
$\xi\leq\frac{1}{\delta}$, where $\delta$ is the maximum vertex degree of
${\cal G}_{b}$111If ${\cal G}_{b}$ is bipartite, we require
$\xi<\frac{1}{\delta}$.. Let the step size $\\{\eta_{\ell}\\}_{\ell\geq 0}$ be
such that $\sum_{\ell=0}^{\infty}\eta_{\ell}=\infty$ and
$\sum_{\ell=0}^{\infty}\eta_{\ell}^{2}<\infty$. Then, the sequence
$\\{[\mathbf{u}^{1}(\ell)^{T},\ldots,\mathbf{u}^{M}(\ell)^{T}]^{T}\\}_{\ell\geq
0}$ converges to the set
${\\{\mathbf{1}\otimes\mathbf{x}~{}|~{}\mathbf{x}\in\mathcal{L}\\}}$222The
symbol $\otimes$ denotes the Kronecker product.. Moreover,
$\frac{1}{M}\sum_{m=1}^{M}\mathbf{u}^{m}(\ell)$ converges to a connected
component of $\mathcal{L}$.

The proof of this theorem follows from Theorem 1 in [22] and the fact that
problem (18)-(19) satisfies certain technical conditions. The full proof is
given in the appendix.

Theorem 1 shows that Alg. 1 converges to a stationary point of the
transmission rate maximization problem, and further, that the MCAs position
vectors $\mathbf{u}^{m}$, $m=1\ldots M$, converge to agreement at this fixed
point.

Algorithm 2 Capacity balancing algorithm.

1:Initialize: $\mathbf{u}_{m}(0)$ to result of Alg. 1

2:for $\ell=0,1,2,\ldots,T_{2}$ do

3: for each MCA $m$ do

4: $\mathbf{d}_{m}(\ell)\leftarrow\sum_{j\in{\cal
N}_{b}(m)}\left(\left(\frac{\int_{\mathcal{V}_{m}}\rho(\mathbf{x})d\mathbf{x}}{c_{m}}-\frac{\int_{V_{j}}\rho(\mathbf{x})d\mathbf{x}}{c_{j}}\right)\times\right.$

5:
$~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\left.n_{mj}\int_{\mathcal{V}_{m}\cap\mathcal{V}_{j}}\rho(\gamma)d\gamma\right)$

6: $\mathbf{u}_{m}(\ell+1)=\mathbf{u}_{m}(\ell)-\alpha\mathbf{d}_{m}(\ell)$

7: end for

8:end for

###  III-D Capacity Balancing Algorithm

While Alg. 1 aims to optimize the Voronoi partitioning for transmission rate,
it may result in some Voronoi regions with workloads that exceed the
capacities of their MCAs. To optimize the system throughput, the MCAs must
have the capacity to process their assigned tasks. We next present a
distributed algorithm to adjust the MCA positions and thus adjust the Voronoi
partitioning, so that the region volumes are proportional to the processing
capacities of the MCAs.

We define the functions $g_{m}(\mathbf{u})=|\mathcal{V}_{m}|^{2}$, $m=1\ldots
M$, and the function $G$:

$\displaystyle G(\mathbf{u})=\sum_{m=1}^{M}\frac{1}{c_{m}}g_{m}(\mathbf{u}).$
(20)

We define the optimization problem:

$\displaystyle\underset{\mathbf{u}_{m}\in\Omega_{c},m=1\ldots
M}{\text{minimize}}~{}~{}G(\mathbf{u}).$ (21)

As shown in [24], the solution to problem (21) creates a Voronoi partitioning
where the total workload assigned to each MCA is proportional to its capacity.
This is formalized in the following lemma.

######  Lemma 1 (Lemma 1 in [24])

Let $x_{m}$, $m=1\ldots M$ be real variables, subject to the constraint that
$\sum_{m=1}^{M}x_{m}=X$ where $X$ is a constant. The function
$\sum_{m=1}^{M}(x_{m}^{2}/c_{m})$, with constants $c_{m}$, attains its minimum
when $x_{m}/c_{m}=X/\sum_{m=1}^{M}c_{m}$, for $i=1\ldots M$.

We present a discrete-time distributed gradient algorithm, based on the
continuous-time solution presented in [24], to solve problem (21). The
pseudocode is shown in Alg. 2. The MCAs initialize their positions to the
positions determined by the transmission rate maximization algorithm. In each
iteration $\ell$, each MCA communicates with its neighbors in the graph ${\cal
G}_{b}$ using the control channels to exchange their values
$\mathbf{u}_{m}(\ell)$. Each MCA $m$ then updates $\mathbf{u}_{m}(\ell+1)$ to
be a weighted average of its $\mathbf{u}_{m}(\ell)$ and that of its neighbor
MCAs:

$\displaystyle\mathbf{u}_{m}(\ell+1)$
$\displaystyle=\mathbf{u}_{m}(\ell)-\alpha\mathbf{d}_{m}(\ell)$ (22)

where $\alpha$ is the step size, and

$\displaystyle\mathbf{d}_{m}(\ell)$ $\displaystyle=\sum_{j\in{\cal
N}_{b}(m)}\left(\left(\frac{\int_{\mathcal{V}_{m}}\rho(\mathbf{x})d\mathbf{x}}{c_{m}}-\frac{\int_{V_{j}}\rho(\mathbf{x})d\mathbf{x}}{c_{j}}\right)\times\right.$
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\left.n_{mj}\int_{\mathcal{V}_{m}\cap\mathcal{V}_{j}}\rho(\gamma)d\gamma\right).$
(23)

The weights used in the averaging depend on the volumes of the Voronoi
regions. Here, $n_{mj}$ denotes the unit normal of the shared boundary of the
Voronoi regions of MCAs $m$ and $j$. Each term in the sum in (23) exerts a
“push” or “pull” on MCA $m$’s position proportional to how balanced $m$’s
capacity is with respect to its neighbor. The MCAs execute this algorithm
until it converges to a balanced partitioning (typically within a hundred
iterations in our experiments). The convergence behavior of Alg. 2 is given in
the following theorem.

###### Theorem 2

For an appropriately chosen step size $\alpha$, Alg. 2 converges
asymptotically to an equilibrium where the MCA positions are such that the
volume of each Voronoi region is proportional to the MCA’s capacity.

Theorem 2 follows directly from Theorem 1 in [24].

By initializing Alg. 2 with positions that optimize for transmission rate, and
then using Alg. 2 to adjust the positions to balance the workload to match the
capacities, the resulting Voronoi partitioning incorporates both transmission
rate maximization and processing capacity to optimize the total system
throughput.

##  IV Simulation Results

We evaluate our collaborative task processing method through numerical
simulations. We first detail our simulation setup, followed by experimental
results on the task field estimation and distributed task processing.

TABLE I: Simulation parameters. Parameters | Value | Parameter | Value
---|---|---|---
Region size | $(5000~{}\text{m})^{2}$ | Sub-region size | $(100~{}\text{m})^{2}$
$H_{c}$ |  $100$ m | $H_{s}$ |  $50$ m
MCA speed |  $25$ m/s | MSA speed |  $10-20$ m/s
Sub-region duration |  $50-60$ s | Sub-region pause |  $2-5$ s
$B$ |  $0.2$ Mbs | $P$ |  $40$ dBm
$\beta$ |  $-50$ dB | $\sigma^{2}$ |  $-60$ dBm
$T_{1}$ | $100$ | $\eta_{\ell}$ | $1e^{-4}$
$T_{2}$ | $200$ | $\alpha$ | $1e^{-2}$
$\xi$ | $1/\delta$ |  |

###  IV-A Simulation Setup

We implemented our simulations using Matlab. To simulate continuous time, we
discretize it into $0.1$s time steps. For Gaussian process regression, we use
the Matlab Statistics and Machine Learning toolbox.

We consider a scenario with $N=50$ MSAs that operate in a
$(5000~{}\text{m})^{2}$ area at a fixed MSL altitude of $H_{s}=50$ m. The MSAs
surveil the region using a modified random waypoint model based on [25]. Each
MSA starts at a randomly initialized location in the region and chooses a
destination at random within the region. The MSA flies to this destination
along the shortest path and then flies within a bounded sub-region of
$(100~{}\text{m})^{2}$ using a random waypoint trajectory for a period of time
chosen uniformly at random between 50s and 60s. The MSA pauses for a random
time between 2s and 5s, and repeats the entire process. The speed of each MSA,
in m/s, is randomly selected from the interval $[10,20]$.

As the MSAs fly, they collect sensing data. Rather than collecting data at a
constant rate, the collection rate increases as the MSAs are close to points
of interests, e.g., representing wildfire fronts. For the sake of comparison,
we standardize the total number of tasks generated per second to $6\times
10^{6}$ bps, and each MSA generates a number of tasks proportional to its
distance from the point(s) of interest. We consider two scenarios:

  1. 1.

Fixed points of interest. We use two points of interest at locations
$\mathbf{p}_{1}$ and $\mathbf{p}_{2}$. The point locations $\mathbf{p}_{1}$
and $\mathbf{p}_{2}$ are selected uniformly at random and remained fixed for
the duration of the simulation. The quantity of tasks generated at MSA $n$ at
time $t$ is proportional to
$\|\mathbf{w}_{n}(t)-\mathbf{p}_{1}\|^{-1.5}+\|\mathbf{w}_{n}(t)-\mathbf{p}_{2}\|^{-1.5}$.

  2. 2.

Moving point of interest. We use a single point of interest, initially at
location ${\mathbf{p}(0)=[4500,4500,0]^{T}}$, which moves $500$m every $15$s
at an angle of $-45^{\circ}$ across the region. The quantity of tasks
generated at MSA $n$ at time $t$ is proportional to
$\|\mathbf{w}_{n}(t)-\mathbf{p}(t)\|^{-1.5}$.

For the task processing experiments, there are ${M=6}$ MCAs. The MCAs fly at
fixed MSL altitude of $H_{c}=100$ m. At the beginning of each window, starting
with window 2, the MCAs pause for 0.1s to execute our algorithm and determine
their new positions. They then fly to their new positions at a constant speed
of $25$ m/s and hover at their new positions until the end of the time window.
The simulation parameters are given in Table I.

###  IV-B Task Field Estimation

We first explore the accuracy of the task field estimation phase. The MSAs
collect data for time window $(0,\Delta]$. We then sample $d_{n}(\Delta)$, the
amount of tasks they collected over these $\Delta$ seconds, from all 50 MSAs.
We use these observations to generate the estimate $\rho$ via Gaussian process
regression. In our solution, this estimate of $\rho$ is used to determine the
MCA positions for the time window $(\Delta,2\Delta]$. We therefore compare the
estimate $\rho$ against the distribution of tasks collected in time window
$(\Delta,2\Delta]$. Since we are interested in studying the estimation
accuracy for a single time window, we use fixed points of interest in these
experiments.

For the purposes of comparison, we discretize the region into
$(50~{}\text{m})^{2}$ cells. For each MSA $n$, we consider the cells that it
visits in window $(\Delta,2\Delta]$, and we assign an equal portion of
$d_{n}(t)$ to each of these cells. The resulting discrete field
${\mathbf{T}\in\mathbb{R}^{100\times 100}}$ is the _discretized task field_
for window $(\Delta,2\Delta]$. We generate the discretized task field estimate
$\mathbf{E}\in\mathbb{R}^{100\times 100}$ by computing the value of $\rho$ at
the center of each cell. To compare $\mathbf{T}$ and $\mathbf{E}$, we use the
normalized mean-square error (NMSE)
${\|\mathbf{E}-\mathbf{T}\|_{F}/(\text{max}(\mathbf{T})-\text{min}(\mathbf{T}))}$.
Here, $\|\cdot\|_{F}$ denotes the Frobenius norm.

Figure 2: Change in NMSE with increasing $\Delta$.

Fig. 2 shows the NMSE for various values of $\Delta$. For each $\Delta$, the
NMSE is the average of three experiments. As can be seen, the NMSE increases
as $\Delta$ increases. This is as anticipated since the smaller the value of
$\Delta$, the closer our approximation is to a continuous time model.

(a) Observations for $\Delta=10$s.

(b) Estimated field for $\Delta=10$s.

(c) Task field for $\Delta=10$s.

(d) Observations for $\Delta=50$s.

(e) Estimated field for ${\Delta=50\text{s}}$.

(f) Task field for $\Delta=50$s.

Figure 3: Examples of the observed task queues, estimated task fields from
these observations $\mathbf{E}$ , and the discretized task fields $\mathbf{T}$
for the window, for $\Delta=10$s and $\Delta=50$s.

This trend is further corroborated in Fig. 3. In Figs. 3(a) and 3(d), we show
the observed values $d_{n}(\Delta)$, each in the cell where its corresponding
MSA $n$ was located at time $\Delta$, for ${\Delta=10\text{s}}$ and
${\Delta=50\text{s}}$, respectively. Figs. 3(b) and 3(e) show the
corresponding estimates of the task fields. We can see that for both values of
$\Delta$, the estimate captures the locations and magnitudes of the
observations. The estimated field is smooth; the benefit of this smoothing is
that it can, in some part, account for the movement of the MSAs throughout the
subsequent window. Figs. 3(c) and 3(f) show the discretized task fields for
the subsequent window. We observe that the estimates reflect the locations of
higher volumes of tasks, however, the magnitudes of the estimate are much
larger than those of the discretized task fields. This is due to the fact that
in our discretized task field, the MSA’s tasks are split over all cells that
it visits in the window. The magnitude discrepancy is larger for $\Delta=50$s
than $\Delta=10$s because the MSAs move for a longer period of time and thus
cover more cells. These results indicate that while the estimates differ in
magnitude from the discretized task fields, they are still valuable for
selecting the locations for MCAs to optimize the throughput.

###  IV-C Collaborative Task Processing

We next explore the performance of our full proposed solution. Based on the
results of the previous section, we select two smaller values of $\Delta$,
$\Delta=10$s and $\Delta=20$s. Each experiment lasts $120$s.

We evaluate three approaches:

  1. 1.

Baseline. The MCA positions remained fixed for the duration of the experiment.

  2. 2.

Transmission rate maximization. Every $\Delta$ seconds, the MCAs generate an
estimate of $\rho$ and then execute the transmission rate maximization
algorithm and update their positions accordingly.

  3. 3.

Full solution. Every $\Delta$ seconds, the MCAs generate an estimate of $\rho$
and then execute the transmission rate maximization algorithm followed by the
capacity balancing algorithm. They then update their positions accordingly.

For all three approaches, the MCAs are initially positioned in two rows of
three, so that each MCA has the same size Voronoi region. For the baseline
approach, each MSA continuously offloads to its closest MCA, provided that the
MCA has capacity to process the workload. If the closest MCA does not have
capacity, the MSA queues its data until the closest MCA has available
capacity. Note that since the MSA moves throughout the experiment, the closest
MCA may be different at different times. For the transmission rate
maximization and full solution approaches, the MCAs start executing the
algorithms after the first $\Delta$ seconds. At the start of each subsequent
window, the MCAs pause for $0.1$s to generate $\rho$ and compute their new
positions. They then move to those new positions. The MSAs continuously
offload to the MCAs assigned to them by the Voronoi partitioning. If an MCA
does not have sufficient capacity, the MSA queues any remaining tasks until it
is assigned to an MCA that does have capacity as in the baseline approach.

For ease of comparison, we set the total processing capacity of the system to
be $6\times 10^{6}$ bps so that it matches the total task generation rate. We
study two scenarios, a _homogenous capacity_ scenario, where each MCA has the
same capacity of $10^{6}$ bps, and a _heterogeneous capacity_ scenario, where
MCA 1 has capacity $2\times 10^{6}$ bps, MCAs 2, 3, and 4 have capacities
$10^{6}$ bps, and MCAs 5 and 6 have capacities $0.5\times 10^{6}$ bps.

For each experiment, we calculate the total amount of tasks that are processed
in two ways. We show the _cold start total_ , where we compute the total
number of tasks processed over the entire $120$s. In all three approaches, the
MCAs have the same positions, and therefore the same amount of tasks processed
in the first window, and so we also calculate the _warm start total_ , which
excludes the first $\Delta$ seconds, so as to only include the time in which
the MCAs adapt their positions and task assignments. We also calculate the
percentage improvement of our full solution over the baseline. All results are
the averages over three experiments.

####  IV-C1 Fixed points of interest

TABLE II: Total tasks processed ($\times 10^{8}$) and warm start percent increase over baseline, with homogeneous capacities and fixed points of interest. Approach | $\mathbf{\Delta=10\text{s}}$ | $\mathbf{\Delta=20\text{s}}$
---|---|---
cold | warm | % inc. | cold | warm | % inc.
Baseline | 4.54 | 4.14 |  | 4.54 | 3.75 |
Rate Max. | 4.56 | 4.17 | 1 | 4.54 | 3.75 | 0
Full Solution | 5.29 | 4.90 | 18 | 5.01 | 4.22 | 13
TABLE III: Total tasks processed ($\times 10^{8}$) and warm start percent increase over baseline, with heterogeneous capacities and fixed points of interest. Approach | $\mathbf{\Delta=10\text{s}}$ | $\mathbf{\Delta=20\text{s}}$
---|---|---
cold | warm | % inc. | cold | warm | % inc.
Baseline | 3.88 | 3.54 |  | 3.88 | 3.21 |
Rate Max. | 3.87 | 3.53 | 0 | 3.86 | 3.19 | -1
Full Solution | 4.78 | 4.45 | 26 | 4.45 | 3.78 | 18

Tables II and III show the results for the homogeneous and heterogeneous
capacities, respectively, for the fixed points of interest scenario. The
transmission rate maximization approach does not show much benefit over the
baseline, indicating that optimizing for transmission rate alone is not
sufficient to maximize the system throughput. We observe that the full
solution yields significant improvement in total tasks processed over the
baseline. Further, the benefit is greater for the smaller value of $\Delta$,
with an $18\%$ improvement over the baseline for the warm start total with
homogeneous capacities, and a $26\%$ improvement for heterogeneous capacities.
These results show the importance of aligning the task assignment with the MCA
capacities—a novel feature of our proposed solution. We also observe that in
all cases, the total amount of tasks processed is less than the total amount
of tasks generated ($7.2\times 10^{8}$ bits for cold start; $6.6\times 10^{8}$
bits and $6.0\times 10^{8}$ bits for warm start for $\Delta=10$s and $20$s,
respectively). We believe this is due to two factors. First, the MCAs spend
part of each window moving to their new positions. During this transition
period, the Voronoi partitioning is sub-optimal with respect to transmission
rates and capacities. Second, any inaccuracies in the estimate of $\rho$ lead
the MCAs to select positions that are less compatible with the true
distribution of tasks. Thus, there is opportunity to improve on our solution
with more sophisticated task field estimation techniques.

(a) Window 1.

(b) Window 4.

(c) Window 8.

Figure 4: Examples of homogeneous MCAs adjusting their positions to the task
field for fixed points of interest for $\Delta=10$s. The triangles represent
the MCAs. The background color represents the intensity of the estimated task
field.

Fig. 4 shows the evolution of the MCA locations and the Voronoi partitioning
over multiple time windows for the homogeneous capacity scenario. The MCA
locations and region boundaries are superimposed on a heat map of $\rho$ for
that window. Initially (window 1), each MCA has a similar-sized Voronoi
region; however, the volume of tasks in each cell is very different. In window
4, the Voronoi partitioning more accurately reflects the task distribution,
and in window 8, the Voronoi partitioning is further adapted to minor changes
in the tasks distribution that result from the MSAs’ mobility and residual
queued tasks.

####  IV-C2 Moving point of interest

TABLE IV: Total tasks processed ($\times 10^{8}$) and warm start percent increase over baseline, with homogeneous capacities and moving point of interest. Approach | $\mathbf{\Delta=10\text{s}}$ | $\mathbf{\Delta=20\text{s}}$
---|---|---
cold | warm | % inc. | cold | warm | % inc.
Baseline | 4.09 | 3.75 |  | 4.09 | 3.42 |
Rate Max. | 4.25 | 3.91 | 4 | 4.25 | 3.58 | 5
Full Solution | 4.88 | 4.55 | 21 | 4.81 | 4.14 | 21
TABLE V: Total tasks processed ($\times 10^{8}$) and warm start percent increase over baseline, with heterogeneous capacities and moving point of interest. Approach | $\mathbf{\Delta=10\text{s}}$ | $\mathbf{\Delta=20\text{s}}$
---|---|---
cold | warm | % inc. | cold | warm | % inc.
Baseline | 3.33 | 3.06 |  | 3.33 | 2.80 |
Rate Max. | 3.52 | 3.25 | 6 | 3.43 | 2.90 | 3
Full Solution | 4.20 | 3.93 | 28 | 4.06 | 3.53 | 26

Tables IV and V show the results for the homogeneous and heterogeneous
capacities, respectively for the moving point of interest. There is a small
increase in the total number of tasks processed for the transmission rate
maximization approach for both capacity types. Since the task field changes as
the point of interest moves, it is intuitive that adapting the MCA positions
should improve the overall transmission rate, thus improving the system
throughput. Our full solution yields significant improvement over the
baseline. For MCAs with homogeneous capacities, our full solution achieved a
$21\%$ improvement in the warm start total for both $\Delta=10$s and
$\Delta=20$s. For the heterogeneous capacities, we observe improvements of
$28\%$ and $26\%$ for $\Delta=10$s and $\Delta=20$s, respectively. This is a
greater improvement than observed for the fixed points of interest, which
shows a key benefit of our method in that it adapts to both small and large
changes in the task field. As in the previous set of experiments, no approach
processes all of the generated tasks; we believe the same reasons apply to
this scenario. We also observe that in most cases, our full solution performs
better with $\Delta=10$s than with $\Delta=20$s. This is expected because the
estimate of $\rho$ is more accurate for the smaller $\Delta$. The improvement
over $\Delta=20$s is not very large though, indicating that it is possible to
save on the motion energy of the MCAs without much decay in the system
throughput.

(a) Window 1.

(b) Window 4.

(c) Window 8.

Figure 5: Examples of heterogeneous MCAs adjusting their positions to the task
field for a moving point of interest for $\Delta=10$s. The triangles represent
the MCAs. The background color represents the intensity of the estimated task
field.

Fig. 5 shows the evolution of the MCA locations and the Voronoi partitioning
over multiple time windows for the heterogeneous capacity scenario with a
moving point of interest. As the point of interest moves, the task field
distribution changes accordingly, and in response, the MCAs adapt their
positions and Voronoi regions. As expected, the MCA position changes are more
pronounced than in the fixed points of interest setting since the changes in
the task field are more significant. These figures illustrate the power of our
method to autonomously adapt to changes in the task demand.

##  V Related Work

Several works have proposed solutions using a single UAV as an MEC server for
users on the ground. The majority of these [5, 13, 14, 15] assume that user
demand is static and known to the UAV, and they propose methods to design a
UAV trajectory that minimizes the motion energy cost while also optimizing for
offloading latency or rate. [16] proposes a solution to optimize the UAV
trajectory when the user demand distribution changes but in a constrained and
known way along straight road segments. In contrast, in our setting, the
demand distribution is continuously changing, with no constraints. Finally,
[7] proposes a solution for a single UAV system that can adapt to changes in
the user demand, however, the limitation to a single MEC UAV simplifies the
problem setting over the one we study.

Several multi-UAV MEC systems for ground users have also been proposed. These
works focus on designing the UAV trajectories and task assignments to optimize
for similar performance measures as the single UAV systems. One line of work
assumes that the user demand is fixed and known. Due to the intractability of
the problem, these works utilize various relaxation approaches and meta-
optimization methods such convex relaxation [17, 11], Deep Reinforcement
Learning [8], and Differential Evolution [12]. The solutions proposed in these
works are all centralized, and this centralization, along with the assumption
about fixed demand, make them unsuitable for our setting. A few works support
dynamic user demand, e.g., [18], which also proposes a solution based on
Differential Evolution, [9], which transforms the problem to a maximum clique
problem, and [10], which formulates the problem as an NP-Hard bin packing
problem and provides an online heuristic solution. The algorithms in these
three works are centralized, which limits their applicability in disaster
scenarios that do not have the infrastructure to support centralized
computation.

The recent work by Chen et al. [6] studies the problem of task offloading from
UAVs to MEC servers on ground vehicles. The work assumes that the trajectories
of the UAVs and ground vehicles are known and cannot be altered, and the UAVs
pay the ground vehicles for use of the MEC servers. They present a centralized
solution based on maximal matching that optimizes for the UAV’s satisfaction
and the ground vehicle’s payoff. The problem considered in this work is
distinctly different from ours in that the MEC server trajectories cannot be
controlled.

Finally, a work that is similar in spirit to ours is [26], which considers a
setting where UAVs serve as aerial base stations for a time-varying continuous
distribution of ground user demand. The authors present a distributed
algorithm that adapts the UAV positions so as to maximize the network coverage
of the ground users. While this work also uses Gaussian process regression to
estimate user demand, it assumes the demand itself is a continuous
distribution, whereas the demand in our problem is generated by individual
UAVs, and we demonstrate that this demand can be approximated using Gaussian
process regression. Further, in the coverage problem, the UAVs are assumed to
have unlimited capacity and can cover any user demand within a certain range.
This is in contrast with our setting, where the MCAs have throughput
capacities that must be respected when determining the MCA positions. Thus,
our objective is mathematically distinct from [26], and further, the inclusion
of the capacity constraints makes the problem even more challenging,
necessitating our multi-phase solution.

##  VI Conclusion

We have presented a novel distributed solution for collaborative task
processing in UAV MEC networks. One key contribution of this work is to show
that the workload can be efficiently approximated by a continuous task field.
By adopting this continuum model, we are able to devise a distributed approach
to maximize the system throughput; the MCAs adaptively update their positions
and the task assignments to maximize the task transmission rates while meeting
their processing capacities. Our solution shows up to 28% improvement over a
non-adaptive baseline strategy in experiments.

This work demonstrates the potential of utilizing continuous models for
problems like task offloading, which are traditionally treated as discrete
optimization problems. As discussed in Section IV-B, the accuracy of our task
field estimation decreases as the time window $\Delta$ increases. To improve
this estimation accuracy, in future work, we plan to investigate more
sophisticated estimation techniques that incorporate MSA mobility prediction.
In addition, we plan to extend our approach to optimize for the MCA motion
energy and use techniques like work stealing to balance workload. Finally, we
aim to relax the model to a continuous time approach, rather than a tumbling
window approach, so that the MCAs move at independently determined times to
optimize for the combination of energy and system throughput.

## Acknowledgment

This research was partially supported by National Science Foundation Grant
CNS-1816307 and Air Force Office of Scientific Research DDDAS Grant
FA9550-19-1-0054.

## Proof of Theorem 1

To prove Theorem 1, we must show that our system model and Alg. 1 satisfy the
assumptions needed for convergence of the projected stochastic gradient
descent algorithm as required by Theorem 1 in [22]. The projected stochastic
gradient descent algorithm addresses problems of the form:

$\displaystyle\underset{\theta\in C}{\text{minimize}}~{}$
$\displaystyle~{}\sum_{m=1}^{M}f_{m}(\theta).$ (24)

It is required that $C$ is a non-empty compact convex set and that the
functions $f_{m}$ are (possibly) non-convex and continuously differentiable.
In our setting, $C$ corresponds to $\Omega_{c}$, which satisfies the
requirements. For example, the square region used for $\Omega_{c}$ in the
experiments is non-empty, compact, and convex. Recall the definition of $f$:

$\displaystyle
f_{m}(\mathbf{u})=\int_{\mathcal{V}_{m}}\rho(\mathbf{x},t){\omega}\left(\mathbf{u}_{m},\mathbf{x}\right)\mathrm{d}\mathbf{x}.$
(25)

We note that $\omega$ is non-decreasing and continuously differentiable. Thus,
as shown in [27, 28], each $f_{m}$ is continuously differentiable.

The final requirement to satisfy Theorem 1 in [22] relates to the weight $\xi$
that the MCAs use when averaging their position estimates. Let $\mathbf{I}$ be
the $M\times M$ identity matrix. Further, let $\mathbf{A}$ be the unweighted
adjacency matrix of ${\cal G}_{b}$, and let $\mathbf{D}$ be the $M\times M$
diagonal matrix where the $(m,m)$th entry is the vertex degree of MCA $m$ in
${\cal G}_{b}$. We define the weight matrix
$\mathbf{W}=\mathbf{I}-\xi(\mathbf{D}-\mathbf{A})$, and note that in Alg. 1,
when MCA $m$ computes a weighted average of the position estimates, it gives
weight $\mathbf{W}_{mj}$ to the estimate of MCA $j$, for each of its neighbors
$j$, and it gives weight $\mathbf{W}_{mm}$ to its own estimate. Theorem 1 in
[22] requires that the spectral radius of $\mathbf{W}$ is strictly less than
1. We first observe that $\mathbf{W}$ is symmetric and so all of its
eigenvalues are real valued. If $\xi$ is such that $\xi\leq\frac{1}{\delta}$
where $\delta$ is the maximum vertex degree of ${\cal G}_{b}$, and ${\cal
G}_{b}$ is not bipartite, then it holds that the eigenvalues of $\mathbf{W}$
lie strictly between $-1$ and $1$ [29], and thus, the requirement on the
spectral radius of $\mathbf{W}$ is satisfied.

## References

  * [1] S. Seneviratne, X. Zhang, M. Adnan, W. Badi, C. Dereczynski, A. Di Luca, S. Ghosh, I. Iskandar, J. Kossin, S. Lewis, F. Otto, I. Pinto, S. Vicente-Serrano, M. Wehner, and B. Zhou, _2021: Weather and Climate Extreme Events in a Changing Climate._ Cambridge University Press, 2021.
  * [2] S. Hayat, E. Yanmaz, and R. Muzaffar, “Survey on unmanned aerial vehicle networks for civil applications: A communications viewpoint,” _IEEE Commun. Surveys Tuts._ , vol. 18, no. 4, pp. 2624–2661, 2016.
  * [3] M. Mozaffari, W. Saad, M. Bennis, Y.-H. Nam, and M. Debbah, “A tutorial on UAVs for wireless networks: Applications, challenges, and open problems,” _IEEE Commun. Surveys Tuts._ , vol. 21, no. 3, pp. 2334–2360, 2019.
  * [4] S. Verykokou, A. Doulamis, G. Athanasiou, C. Ioannidis, and A. Amditis, “UAV-based 3D modeling of disaster scenes for urban search and rescue,” in _Proc. IEEE Int. Conf. Imaging Syst. Tech._ , 2016, pp. 106–111.
  * [5] M. Li, N. Cheng, J. Gao, Y. Wang, L. Zhao, and X. Shen, “Energy-efficient UAV-assisted mobile edge computing: Resource allocation and trajectory optimization,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 3, pp. 3424–3438, 2020.
  * [6] W. Chen, Z. Su, Q. Xu, T. H. Luan, and R. Li, “VFC-based cooperative UAV computation task offloading for post-disaster rescue,” in _Proc. IEEE Conf. Comput. Commun._ , 2020, pp. 228–236.
  * [7] Y. Liu, K. Xiong, Q. Ni, P. Fan, and K. B. Letaief, “UAV-assisted wireless powered cooperative mobile edge computing: Joint offloading, CPU control, and trajectory optimization,” _IEEE Internet Things J._ , vol. 7, no. 4, pp. 2777–2790, 2020.
  * [8] N. Zhao, Z. Ye, Y. Pei, Y.-C. Liang, and D. Niyato, “Multi-agent deep reinforcement learning for task offloading in UAV-assisted mobile edge computing,” _IEEE Trans. Wireless Commun._ , pp. 1–1, 2022.
  * [9] Z. Liao, Y. Ma, J. Huang, J. Wang, and J. Wang, “HOTSPOT: A UAV-assisted dynamic mobility-aware offloading for mobile-edge computing in 3-D space,” _IEEE Internet Things J._ , vol. 8, no. 13, pp. 10 940–10 952, 2021.
  * [10] J. Wang, K. Liu, and J. Pan, “Online UAV-mounted edge server dispatching for mobile-to-mobile edge computing,” _IEEE Internet Things J._ , vol. 7, no. 2, pp. 1375–1386, 2020.
  * [11] X. Huang, X. Yang, Q. Chen, and J. Zhang, “Task offloading optimization for UAV-assisted fog-enabled Internet of Things networks,” _IEEE Internet Things J._ , vol. 9, no. 2, pp. 1082–1094, 2022.
  * [12] L. Yang, H. Yao, J. Wang, C. Jiang, A. Benslimane, and Y. Liu, “Multi-UAV-enabled load-balance mobile-edge computing for IoT networks,” _IEEE Internet Things J._ , vol. 7, no. 8, pp. 6898–6908, 2020\.
  * [13] Q. Hu, Y. Cai, G. Yu, Z. Qin, M. Zhao, and G. Y. Li, “Joint offloading and trajectory design for UAV-enabled mobile edge computing systems,” _IEEE Internet Things J._ , vol. 6, no. 2, pp. 1879–1892, 2019\.
  * [14] F. Zhou, Y. Wu, R. Q. Hu, and Y. Qian, “Computation rate maximization in UAV-enabled wireless-powered mobile-edge computing systems,” _IEEE J. Sel. Areas Commun._ , vol. 36, no. 9, pp. 1927–1941, 2018.
  * [15] Z. Na, Y. Liu, J. Shi, C. Liu, and Z. Gao, “UAV-supported clustered NOMA for 6G-enabled Internet of Things: Trajectory planning and resource allocation,” _IEEE Internet Things J._ , vol. 8, no. 20, pp. 15 041–15 048, 2021.
  * [16] L. Zhang, Z. Zhao, Q. Wu, H. Zhao, H. Xu, and X. Wu, “Energy-aware dynamic resource allocation in UAV assisted mobile edge computing over social Internet of Vehicles,” _IEEE Access_ , vol. 6, pp. 56 700–56 715, 2018\.
  * [17] H. Wang, H. Ke, and W. Sun, “Unmanned-aerial-vehicle-assisted computation offloading for mobile edge computing based on deep reinforcement learning,” _IEEE Access_ , vol. 8, pp. 180 784–180 798, 2020.
  * [18] H. Yan, W. Bao, X. Zhu, J. Wang, and L. Liu, “Data offloading enabled by heterogeneous UAVs for IoT applications under uncertain environments,” _IEEE Internet Things J._ , pp. 1–1, 2022.
  * [19] T. X. Tran and D. Pompili, “Joint task offloading and resource allocation for multi-server mobile-edge computing networks,” _IEEE Trans. Veh. Technol._ , vol. 68, no. 1, pp. 856–868, 2019.
  * [20] K. Cheng, Y. Teng, W. Sun, A. Liu, and X. Wang, “Energy-efficient joint offloading and wireless resource allocation strategy in multi-MEC server systems,” in _Proc. IEEE Int. Conf. Commun._ , 2018, pp. 1–6.
  * [21] E. Schulz, M. Speekenbrink, and A. Krause, “A tutorial on Gaussian process regression: Modelling, exploring, and exploiting functions,” _J. Math. Psychol._ , vol. 85, pp. 1–16, 2018.
  * [22] P. Bianchi and J. Jakubowicz, “Convergence of a multi-agent projected stochastic gradient algorithm for non-convex optimization,” _IEEE Trans. Autom. Control_ , vol. 58, no. 2, pp. 391–405, 2012.
  * [23] J. Zeng and W. Yin, “On nonconvex decentralized gradient descent,” _IEEE Trans. Signal Process._ , vol. 66, no. 11, pp. 2834–2848, 2018.
  * [24] H. Sayyaadi and M. Moarref, “A distributed algorithm for proportional task allocation in networks of mobile agents,” _IEEE Trans. Autom. Control_ , vol. 56, no. 2, pp. 405–410, 2011.
  * [25] A. Bujari, C. T. Calafate, J.-C. Cano, P. Manzoni, C. E. Palazzi, and D. Ronzani, “Flying ad-hoc network application scenarios and mobility models,” _Int. J. Distrib. Sens. Netw._ , vol. 13, no. 10, p. 1550147717738192, 2017. [Online]. Available: https://doi.org/10.1177/1550147717738192
  * [26] T. Kimura and M. Ogura, “Distributed collaborative 3D-deployment of UAV base stations for on-demand coverage,” in _Proc. IEEE Conf. Comput. Commun._ , 2020, pp. 1748–1757.
  * [27] Q. Du, V. Faber, and M. Gunzburger, “Centroidal Voronoi tessellations: Applications and algorithms,” _SIAM review_ , vol. 41, no. 4, pp. 637–676, 1999.
  * [28] J. Cortes, S. Martinez, T. Karatas, and F. Bullo, “Coverage control for mobile sensing networks,” _IEEE trans. Robot. Autom._ , vol. 20, no. 2, pp. 243–255, 2004.
  * [29] W. Anderson and T. Morley, “Eigenvalues of the Laplacian of a graph,” _Lin. Multilin. Algebra_ , vol. 18, no. 2, pp. 141–145, 1985.


By Cauchy-Schwarz inequality and since $\left\vert Z \right\vert \leq \alpha$, we obtain

$$ \alpha^n + \beta^n \leq \left(\alpha^{n-2} + \beta^{n-2}\right) \alpha \beta $$

so that since $2\alpha \beta \leq \alpha^2 + \beta^2$,

$$ \alpha^n + \beta^n \leq \left( \alpha^2 \beta^{n-2} + \beta^2 \alpha^{n-2} \right)$$

which implies that

$$ \left(\alpha^{n-2} - \beta^{n-2}\right)\left(\alpha^2 -\beta^2\right) \leq 0 $$

and that all the previous inequalities are in fact equalities: $\alpha = \beta = \left\vert Z \right\vert $ and $\left\langle X,Z \right\rangle = \alpha^2$. In particular $Z = X$ and $X_k$ strongly converges to $Z = X$. Since the limit is independant from the subsequence, Step 2 is proved.

STEP 3:  Conclusion: since $\left\vert \nabla u_k -\nabla u \right\vert$ is bounded in $W^{1,n}$ and converges almost everywhere to $0$ in $(M,g)$, then, $\left\vert \nabla u_k -\nabla u \right\vert$ strongly converges to $0$ in $L^q(M,g)$ for any $q < n$.

§.§ Point removability

Let $\Psi : \mathbb{B}\setminus\{0\} \to \mathbb{R}^{\mathbb{N}}$ be a $\mathcal{C}^1$ $n$-harmonic map into an infinite dimensional sphere such that

$$\int_{\mathbb{B}}  \left\vert \nabla \Psi \right\vert_g^n dv_g < +\infty \text{ and }  \left\vert \nabla \Psi \right\vert(x) \leq \frac{C}{\left\vert x \right\vert} $$

then $\Psi : \mathbb{B}\setminus\{0\} \to \mathbb{R}^{\mathbb{N}}$ extends to a $\mathcal{C}^{1,\alpha}$ function on $\mathbb{B}$ which is $n$-harmonic.

We can follow and generalize the proof by [31] to infinite dimensional target spheres.

§.§ Regular n-harmonic maps are locally $(\tau,n)$-approximated harmonic maps

The aim of the section is to prove the following for $\mathcal{C}^1$ harmonic maps into a possibly infinite dimensional target sphere.

There is $\delta_{n,g}>0$ such that for any weak $n$-harmonic map $\Phi : M \to \mathbb{R}^{\mathbb{N}}$ into a possibly infinite dimensional target sphere, such that $\Phi \in \mathcal{C}^{1}$ and for any ball $\mathbb{B}_r(p)$ such that

$$ \int_{\mathbb{B}_{2r}(p)} \left\vert \nabla \Phi \right\vert^n_g dv_g \leq \delta_{n,g} $$

for any  $\left(\Psi_{\tau,m}\right)_{\tau>0, n\in \mathbb{N}}$ minimizers of

$$\Psi \mapsto \int_{\mathbb{B}_r(p)} \left( \left\vert \nabla \Psi\right\vert_g^2 + \tau \right)^{\frac{n}{2}} dv_g   $$

under the $W^{1,n}$ maps $\Psi : M \to \mathbb{S}^m$ such that $\left\vert \Psi \right\vert^2 =_{a.e} 1$ and $\Psi = \frac{(\phi_0,\cdots,\phi_m)}{\left(\sum_{i=0}^m (\phi_i)^2 \right)^{\frac{1}{2}}}$ on $\partial \mathbb{B}_{r}(p)$ we have that $\Psi_{\tau,m}$ converges to $\Phi$ as $(\tau,m) \to (0,+\infty)$.

Step 1:  Since $\Phi$ is a $\mathcal{C}^1$ function, it satisfies the a priori estimate on the gradient for $x \in \mathbb{B}_r(p)$

$$\left\vert \nabla \Phi(x) \right\vert_g^2 \leq C \int_{\mathbb{B}_{2r}(p)} \left\vert \nabla \Phi \right\vert^n_g dv_g \leq C \delta_{n,g} $$

We have for $x \in \mathbb{B}_r(p)$ that

$$ \left\vert \phi_1(0) - \phi_1(x) \right\vert^2 \leq r^2 C_g \left\vert \nabla \Phi(z) \right\vert_g^2 \leq  C_g C \delta_{n,g}$$

so that taking  $\delta_{n,g}$ small enough, and up to a rotation of coordinates so that $\phi_1(0) = 1$ we can assume that

$$ \forall x \in \mathbb{B}_r(p) \left\vert\phi_1(x) \right\vert \geq \frac{3}{4}. $$

Moreover, up to reduce $\delta_{n,g}$, and for $\tau$ small enough $\Psi_{\tau}$ has to satisfy

$$ \forall x \in \mathbb{B}_r(p) \left\vert\left(\psi_{\tau,m}\right)_1(x) \right\vert \geq \frac{1}{2} $$

thanks to the same trick as in the proof of Step 1 of Claim <ref> based on Courant-Lebesgue lemma and Morrey-Sobolev embeddings.

Step 2:  We set $\tilde \Phi_m = \frac{\Phi_m}{\left\vert \Phi_m \right\vert}$ where $ \Phi_m = (\phi_0,\cdots,\phi_m, 0,\cdots)  $. We denote $\tilde{\Psi}_{\tau,m} : M \to \mathbb{R}^{\mathbb{N}}$ the minimizer of

$$\Psi \mapsto \int_{\mathbb{B}_r(p)} \left( \left\vert \nabla \Psi\right\vert_g^2 + \tau \right)^{\frac{n}{2}} dv_g   $$

under the $W^{1,n}$ maps $\Psi : M \to \mathbb{S}^m$ such that $\left\vert \Psi \right\vert^2 =_{a.e} 1$ and $\Psi =  \tilde \Phi_m$ on $\partial \mathbb{B}_{r}(p)$. we also set

$$\Psi_{\tau,m} = \left\vert \Phi_m \right\vert \tilde{\Psi}_{\tau,m} + \Phi- \Phi_m \text{ so that } \Psi_{\tau,m} - \Phi =  \left\vert \Phi_m \right\vert \left( \tilde{\Psi}_{\tau,m} - \tilde\Phi_m  \right).$$

We have:

\begin{equation*}
\begin{split}\int_{\mathbb{B}_r(p)} \left\vert \nabla\left( \Phi -  \Psi_{\tau,m}\right) \right\vert_g^2 \left\vert \nabla \Phi \right\vert_g^{n-2} - \left(  \int_{\mathbb{B}_r(p)} \left\vert \nabla \Psi_{\tau,m} \right\vert_g^{2} \left\vert \nabla \Phi \right\vert_g^{n-2}  - \int_{\mathbb{B}_r(p)} \left\vert \nabla \Phi \right\vert_g^n \right) \\
= 2 \int_{\mathbb{B}_r(p)} \left\vert \nabla \Phi \right\vert_g^{n-2} \nabla \Phi . \nabla\left( \Phi -  \Psi_{\tau,m}\right)
= 2 \int_{\mathbb{B}_r(p)} \left\vert \nabla \Phi \right\vert_g^{n}  \Phi . \left( \Phi -  \Psi_{\tau,m} \right) \\
=  \int_{\mathbb{B}_r(p)} \left\vert \nabla \Phi \right\vert_g^{n}  \left\vert \Phi -  \Psi_{\tau,m} \right\vert^2
\end{split}
\end{equation*}

and we obtain

\begin{equation*}
\begin{split}\int_{\mathbb{B}_r(p)} \left\vert \nabla\left( \Phi -  \Psi_{\tau,m}\right) - \nabla \ln \phi_1\left( \Phi -  \Psi_{\tau,m}\right)  \right\vert_g^2 \left\vert \nabla \Phi \right\vert_g^{n-2}  \\
= \int_{\mathbb{B}_r(p)} \left\vert \nabla\left( \Phi -  \Psi_{\tau,m}\right) \right\vert_g^2 \left\vert \nabla \Phi \right\vert_g^{n-2} - \int_{\mathbb{B}_r(p)} \left\vert \nabla \Phi \right\vert_g^{n}  \left\vert \Phi -  \Psi_{\tau,m} \right\vert^2
\end{split}
\end{equation*}

so that

\begin{equation*}
\begin{split}
\int_{\mathbb{B}_r(p)} \left\vert \nabla\left( \Phi -  \Psi_{\tau,m}\right) - \nabla \ln \phi_1\left( \Phi -  \Psi_{\tau,m}\right)  \right\vert_g^2 \left\vert \nabla \Phi \right\vert_g^{n-2}  = \int_{\mathbb{B}_r(p)} \left\vert \nabla \Psi_{\tau,m} \right\vert_g^{2} \left\vert \nabla \Phi \right\vert_g^{n-2}  \\
- \int_{\mathbb{B}_r(p)} \left\vert \nabla \Phi \right\vert_g^n \\
= \int_{\mathbb{B}_r(p)} \left\vert \nabla \Phi \right\vert_g^{n-2} \left\vert \Phi_m \right\vert^2 \left(  \left\vert \nabla \tilde{\Psi}_{\tau,m} \right\vert_g^{2} - \left\vert \nabla \tilde{\Phi}_m \right\vert_g^{2} \right)
\end{split}
\end{equation*}

Moreover, we have by similar computations that

\begin{equation*}
\begin{split}
\int_{\mathbb{B}_r(p)} & \left\vert \nabla\left( \tilde{\Phi}_m -  \tilde\Psi_{\tau,m}\right) - \nabla \ln \left(\tilde\Psi_{\tau,m}\right)_1\left( \tilde\Phi_m -  \tilde\Psi_{\tau,m}\right)  \right\vert_g^2  \left(\left\vert \nabla \tilde\Psi_{\tau,m} \right\vert_g^2 + \tau\right)^{\frac{n-2}{2}} \\
= & \int_{\mathbb{B}_r(p)} \left\vert \nabla \tilde\Phi_m \right\vert_g^{2}  \left(\left\vert \nabla \tilde\Psi_{\tau,m} \right\vert_g^2 + \tau\right)^{\frac{n-2}{2}}  - \int_{\mathbb{B}_r(p)} \left\vert \nabla \tilde\Psi_{\tau,m} \right\vert_g^2 \left(\left\vert \nabla \tilde\Psi_{\tau,m} \right\vert_g^2 + \tau\right)^{\frac{n-2}{2}} \\
\end{split}
\end{equation*}

and we obtain

\begin{equation*}
\begin{split}
A:= & \int_{\mathbb{B}_r(p)}  \left\vert \nabla\left( \tilde\Phi_m -  \tilde\Psi_{\tau,m}\right) - \nabla \ln \left(\tilde\Psi_{\tau,m}\right)_1\left( \tilde\Phi_m -  \tilde\Psi_{\tau,m}\right)   \right\vert_g^2   \left( \left\vert \nabla \tilde\Psi_{\tau,m} \right\vert_g^2 + \tau\right)^{\frac{n-2}{2}}   \\
& + \int_{\mathbb{B}_r(p)} \left\vert \nabla\left( \Phi -  \Psi_{\tau,m}\right) - \nabla \ln \phi_1\left( \Phi -  \Psi_{\tau,m}\right)   \right\vert_g^2 \left\vert \nabla \Phi \right\vert_g^{n-2}÷   \\
\leq & \int_{\mathbb{B}_r(p)}\left(  \left\vert \nabla \tilde{\Psi}_{\tau,m} \right\vert_g^{2} - \left\vert \nabla \tilde{\Phi}_m \right\vert_g^{2}  \right)\left( \left\vert \nabla \Phi \right\vert_g^{n-2} \left\vert \Phi_m \right\vert^2- \left(\left\vert \nabla \tilde\Psi_{\tau,m} \right\vert_g^2 + \tau\right)^{\frac{n-2}{2}}  \right) \\
\end{split}
\end{equation*}

so that splitting the right-hand side into three terms, we have

\begin{equation*}
\begin{split}
A \leq & \int_{\mathbb{B}_r(p)}\left(  \left\vert \nabla \tilde{\Psi}_{\tau,m} \right\vert_g^{2} - \left\vert \nabla \tilde{\Phi}_m \right\vert_g^{2}  \right)\left( \left\vert \nabla \tilde\Phi_m \right\vert_g^{n-2} -\left\vert \nabla \tilde\Psi_{\tau,m} \right\vert_g^{n-2}   \right) \\
& + \int_{\mathbb{B}_r(p)}\left(  \left\vert \nabla \tilde{\Psi}_{\tau,m} \right\vert_g^{2} - \left\vert \nabla \tilde{\Phi}_m \right\vert_g^{2}  \right)\left( \left\vert \nabla \Phi \right\vert_g^{n-2} \left\vert \Phi_m \right\vert^2 -  \left\vert \nabla \tilde{\Phi}_m \right\vert_g^{n-2} \right) \\
& + \int_{\mathbb{B}_r(p)}\left(  \left\vert \nabla \tilde{\Psi}_{\tau,m} \right\vert_g^{2} - \left\vert \nabla \tilde{\Phi}_m \right\vert_g^{2}  \right) \left(  \left\vert \nabla \tilde\Psi_{\tau,m} \right\vert_g^{n-2} - \left(\left\vert \nabla \tilde\Psi_{\tau,m} \right\vert_g^2 + \tau\right)^{\frac{n-2}{2}}  \right) \end{split}
\end{equation*}

and since the first right-hand term is non positive, we obtain

\begin{equation*}
\begin{split}
A \leq & \int_{\mathbb{B}_r(p)}\left(  \left\vert \nabla \tilde{\Psi}_{\tau,m} \right\vert_g^{2} + \left\vert \nabla \tilde{\Phi}_m \right\vert_g^{2}  \right)\left\vert \left\vert \nabla \Phi \right\vert_g^{n-2} \left\vert \Phi_m \right\vert^2 -  \left\vert \nabla \tilde{\Phi}_m \right\vert_g^{n-2} \right\vert \\
& + \int_{\mathbb{B}_r(p)}\left(  \left\vert \nabla \tilde{\Psi}_{\tau,m} \right\vert_g^{2} + \left\vert \nabla \tilde{\Phi}_m \right\vert_g^{2}  \right) \left\vert  \left\vert \nabla \tilde\Psi_{\tau,m} \right\vert_g^{n-2} - \left(\left\vert \nabla \tilde\Psi_{\tau,m} \right\vert_g^2 + \tau\right)^{\frac{n-2}{2}}  \right\vert .
\end{split}
\end{equation*}

Knowing that

$$ \left\vert \nabla \Phi \right\vert^2 = \left\vert \Phi_m \right\vert^2\left\vert \nabla \tilde\Phi_m \right\vert^2 + \left\vert \nabla\left\vert \Phi_m \right\vert \right\vert^2 + \left\vert \nabla \left(\Phi-\Phi_m\right) \right\vert_g^2  $$

and that $\Phi_m$ and $\nabla \Phi_m$ converge a.e to $\Phi$ and $\nabla \Phi$, it is clear by the dominated convergence theorem that the first integral converges to $0$ as $m\to +\infty$. Thanks to the a priori estimates of Claim <ref>, up to a subsequence, $\Psi$ is the limit as $(\tau,m) \to (0,+\infty)$ of $\tilde\Psi_{\tau,m}$. We easily obtain that $\Psi$ is also the limit of $\Psi_{\tau,m}$ and that the second integral in the previous inequality converges to $0$ as $(\tau,m)\to (0,+\infty)$. Moreover, if we set

$$ Z = \{ x \in \mathbb{B}_r(p) ; \left\vert \nabla \Psi \right\vert_g^2 + \left\vert \nabla \Phi \right\vert_g^2 = 0  \}, $$

we have

$$ \left\vert \nabla\left( \Phi -  \Psi \right) - \nabla \ln \psi_1\left( \Phi -  \Psi \right)   \right\vert_g^2 + \left\vert \nabla\left( \Phi -  \Psi \right) - \nabla \ln \phi_1\left( \Phi -  \Psi \right)   \right\vert_g^2 = 0 \text{ in } \mathbb{B}_r(p)\setminus Z. $$

so that

\begin{equation*}
\begin{split}
\int_{\mathbb{B}_r(p)}  \left\vert \nabla\left( \Phi -  \Psi \right)  \right\vert_g^2 = \int_{\mathbb{B}_r(p)\setminus Z}   \left\vert \nabla\left( \Phi -  \Psi \right)  \right\vert_g^2 \leq 2 \int_{\mathbb{B}_r(p)\setminus Z}  \left(   \left\vert \nabla \ln \psi_1 \right\vert^2 + \left\vert \nabla \ln \phi_1 \right\vert^2 \right) \left\vert \Phi-\Psi \right\vert^2 \\
\leq 2 C \delta_{n,g}^{\frac{2}{n}} \int_{\mathbb{B}_r(p)} \frac{\left\vert \Phi-\Psi \right\vert^2}{\left(r-\left\vert x \right\vert\right)^2} \leq  \tilde{C} \delta_{n,g}^{\frac{2}{n}}   \int_{\mathbb{B}_r(p)}  \left\vert \nabla\left( \Phi -  \Psi \right)  \right\vert_g^2
\end{split}
\end{equation*}

where the second inequality comes from classical weak estimates on the gradient when it satisfies a $L^{\infty}$ $\eps$-regularity property and the third one comes from a classical Hardy inequality. Letting $\delta_{n,g}$ be small enough gives that $\Phi= \Psi$.

[1]
B. Ammann, E. Humbert,  The second Yamabe invariant, J. Funct. Anal., 235, 2006, no.2, 377–412
[2]
B. Colbois, A. El Soufi,
Extremal eigenvalues of the Laplacian in a conformal class of metrics: the `conformal spectrum',
Ann. Global Anal. Geom. 24 (2003), no.4, 337–349.
[3]
B. Colbois, J. Dodziuk, Riemannian metrics with large $\lambda_1$, Proc. Amer. Math. Soc, 122, 1994, no.3, 905–906
[4]
Y. Colin de Verdière, Sur la multiplicité de la première valeur propre non nulle du laplacien, Comment. Math. Helv., 61, 1986, no.2, 254–270
[5]
P. Courilleau, A compactness result for $p$-harmonic maps, Differential Integral Equations,  14, 2001, no.1, 75–84
[6]
A. El Soufi, S. Ilias, Le volume conforme et ses applications d'apres Li et Yau, Sém. Théorie Spectrale et Géométrie, Institut Fourier, 1983?1984, No.VII, (1984).

[7]
A. El Soufi, S. Ilias, Immersions minimales, première valeur propre du laplacien et volume conforme, Mathematische Annalen, 1986, 275, 257-267
[8]
A. El Soufi, S. Ilias, Riemannian manifolds admitting isometric immersions by their first eigenfunctions, Pacific J. Math. 195, 2000, 91–99.

[9]
A. El Soufi and S. Ilias, Extremal metrics for the first eigenvalue of the Laplacian in a conformal class, Proc. Amer. Math. Soc. 131, 2003, 1611-1618.
[10]
A. El Soufi, S. Ilias, Laplacian eigenvalue functionals and metric deformations on compact manifolds, Journal of Geometry and Physics, 58, Issue 1, January 2008, 89-104
[11]
A. Fraser, R. Schoen, Sharp eigenvalue bounds and minimal surfaces in the ball, Invent. Math. 203, 2016, 823–890.
[12]
A. Fraser, R. Schoen, Some results on higher eigenvalue optimization,  Calc. Var. Partial Differential Equations, 59, 2020, 5, Paper No. 151, 22
[13]
L. Friedlander, N. Nadirashvili, A differential invariant related to the first eigenvalue of
the Laplacian, Internat. Math. Res. Notices, 1999, 17, 939–952,

[14]
M.J. Gursky, S. Pérez-Ayala,  Variational properties of the second eigenvalue of the conformal Laplacian, J. Funct. Anal., 282, 2022, no.8, Paper No. 109371, 60
[15]
Q. Han, F. Lin , Elliptic partial differential equations, Courant Lecture Notes in Mathematics, 1, New York University, Courant Institute of Mathematical
Sciences, New York; American Mathematical Society, Providence,
RI, 1997, x+144,
[16]
A. Hassannezhad, Conformal upper bounds for the eigenvalues of the Laplacian and Steklov problem, J. Funct. Anal., 261, 2011, no 12, 3419–3436
[17]
J. Hersch, Quatre propriétés isopérimétriques de
membranes sphériques homogènes, C.R. Acad. Sci. Paris Sér. A-B 270, 1970, A1645–A1648.
[18]
M. Karpukhin,
Index of minimal spheres and isoperimetric eigenvalue inequalities, Invent. Math., 223, 2021, no. 1, 335–377

[19]
M. Karpukhin,  V. Medvedev, On the Friedlander-Nadirashvili invariants of surfaces, Math. Ann., 379, 2021, 3-4, 1767–1805,
[20]
M. Karpukhin, A. Métras, Laplace and Steklov extremal metrics via $n$-harmonic maps, J. Geom. Anal., 32, 2022, no.5, Paper No. 154, 36
[21]
M.A. Karpukhin, N. Nadirashvili, A. Penskoi, I.  Polterovich, An isoperimetric inequality for Laplace eigenvalues on the sphere. arXiv: 1706.05713.
To appear in J. Differential Geom..

[22]
M.A. Karpukhin, N. Nadirashvili, A. Penskoi, I.  Polterovich, Conformally maximal metrics for Laplace eigenvalues on surfaces, arXiv:2003.02871v2

[23]
M.A. Karpukhin, D. L. Stern Min-max harmonic maps and a new characterization of conformal eigenvalues
arXiv preprint (2020),	arXiv:2004.04086, 69 pp
[24]
M.A. Karpukhin, D. L. Stern Existence of harmonic maps and eigenvalue optimization in higher dimensions
arXiv preprint (2022), arXiv:2207.13635, 60 pp

[25]
G. Kokarev, Variational aspects of Laplace eigenvalues on Riemannian surfaces, Adv. Math. 258, 2014, 191–239.

[26]
H.N. Kim, Maximization of the second Laplacian eigenvalue on the
sphere, Proc. Amer. Math. Soc., Proceedings of the American Mathematical Society, 150, 2022, no.8, 3501–3512
[27]
N. Korevaar, Upper bounds for eigenvalues of conformal metrics, J. Differential Geom., 37, 1993, 1, 73–93

[28]
R. Hardt, F-H Lin, Mappings minimizing the $L^p$ norm of the gradient, Comm. Pure Appl. Math., Communications on Pure and Applied Mathematics, 40, 1987, no.5, 555–588
[29]
P. Li, S.T. Yau, A new conformal invariant and its applications to the Willmore conjecture and the first eigenvalue of compact surfaces, Invent.. Math. 69, 1982, 269–291.
[30]
H. Matthiesen, A. Siffert,
Handle attachement and the normalized first eigenvalue, preprint.
[31]
L. Mou, P. Yang, Regularity for $n$-harmonic maps, J. Geom. Anal., 6, 1996, no.1, 91–112
[32]
N. Nadirashvili, Berger's isoperimetric problem and minimal immersions of surfaces, Geom. Func. Anal. 6, 1996, 877-897.

[33]
N. Nadirashvili, Isoperimetric inequality for the second eigenvalue of a sphere, J. Differential Geom. 61, 2002, 335–340.
[34]
R. Petrides, Maximization of the second conformal eigenvalue of spheres, Proc. Amer. Math. Soc., 142, 2014, 7 , 2385–2394

[35]
R. Petrides, On a rigidity result for the first conformal eigenvalue of the
Laplacian, J. Spectr. Theory, 5, 2015, no.1, 227–234
[36]
R. Petrides, Existence and regularity of maximal metrics for the first Laplace eigenvalue on surfaces, Geom. Funct. Anal. 24, 2014, 1336–1376.

[37]
R. Petrides, On the existence of metrics which maximize Laplace eigenvalues on surfaces, Int. Math. Res. Not. , 14, 2018, 4261–4355.
[38]
R. Petrides, A variational method for functionals depending on eigenvalues, submitted
[39]
R. Petrides, D. Tewodrose, Subdifferentials and critical points of eigenvalue functionals
[40]
S. Sarsa, Note on an elementary inequality and its application to the
regularity of $p$-harmonic functions, Ann. Fenn. Math., 47, 2022, no.1, 139–153
[41]
Strzelecki, Paweł, Regularity of $p$-harmonic maps from the $p$-dimensional ball into a sphere, Manuscripta Math., 82, 1994, no.3-4, 407–415
[42]
K. Uhlenbeck, K., Regularity for a class of non-linear elliptic systems, Acta Math., 138, 1977, no.3-4, 219–240
[43]
P.C. Yang, S.T. Yau, Eigenvalues of the Laplacian of compact Riemannian surfaces and minimal submanifolds, Ann. Scuola Norm. Sup. Pisa Cl. Sci. (4) 7, 1980, no.1, 53–63.


$\operatorname{\underline{Isom}}^{\otimes}(\omega_{1},\omega_{2})\times\operatorname{\underline{Aut}}^{\otimes}(\omega_{1})\operatorname{\longrightarrow}\operatorname{\underline{Isom}}^{\otimes}(\omega_{1},\omega_{2})$

of $\operatorname{\underline{Aut}}^{\otimes}(\omega_{1})$ on
$\operatorname{\underline{Isom}}^{\otimes}(\omega_{1},\omega_{2})$, given by
pre-composition. We note that
$\operatorname{\underline{Aut}}^{\otimes}(\omega_{1})=\operatorname{\underline{Aut}}^{\otimes}(\omega^{G}\otimes
K^{\operatorname{nr}})\cong G_{K^{\operatorname{nr}}}$,222For this
isomorphism, we refer to the discussion above Proposition 3.11. so
$\operatorname{\underline{Isom}}^{\otimes}(\omega_{1},\omega_{2})$ may be
viewed as a $G_{K^{\operatorname{nr}}}$-torsor over $K^{\operatorname{nr}}$.
By Lemma 4.20, we have
$H^{1}(K^{\operatorname{nr}},G_{K^{\operatorname{nr}}})=1$. Thus,
$\operatorname{\underline{Isom}}^{\otimes}(\omega_{1},\omega_{2})$ is
isomorphic to the trivial $G_{K^{\operatorname{nr}}}$-torsor over
$K^{\operatorname{nr}}$, i.e., we have
$\operatorname{\underline{Isom}}^{\otimes}(\omega_{1},\omega_{2})_{K^{\operatorname{nr}}}\cong
G_{K^{\operatorname{nr}}}$.

On the other hand, we have an isomorphism
$\gamma\colon\omega_{2}\otimes\mathcal{B}_{0}\to\omega_{1}\otimes\mathcal{B}_{0}$
of tensor functors, induced by the $\mathcal{B}_{0}$-linear extension of the
inclusion

${(V_{\mathcal{B}_{0}})^{\Theta_{X_{0}}=0}}$${V_{\mathcal{B}_{0}}}$

for all $(V,\rho)\in\operatorname{\mathbf{Rep}}_{F}(G)$. We now fix
$\beta\in\operatorname{\underline{Isom}}^{\otimes}(\omega_{1},\omega_{2})(K^{\operatorname{nr}})$,
we then have an element
$\tilde{\beta}\coloneqq\gamma\circ\beta_{\mathcal{B}_{0}}\in\operatorname{\underline{Aut}}^{\otimes}(\omega_{1}\otimes\mathcal{B}_{0})(\mathcal{B}_{0})=G_{\mathcal{B}_{0}}$.
Let $b\in G(\mathcal{B}_{0})$ be the inverse of the image of $\tilde{\beta}$
under the following isomorphism

$\operatorname{\underline{Aut}}^{\otimes}(\omega_{1}\otimes\mathcal{B}_{0})(\mathcal{B}_{0})\operatorname{\longrightarrow}G_{\mathcal{B}_{0}}(\mathcal{B}_{0})=G(\mathcal{B}_{0}).$

Since $F[G]$ is finitely presented over $F$, the functor
$\operatorname{Hom}_{\operatorname{\mathbf{Alg}}_{F}}(F[G],\rule{5.69054pt}{0.42677pt})$
commutes with colimits. We have

$G(\mathcal{B}_{0})=G(\varinjlim_{L}\mathcal{R}_{L})=\varinjlim_{L}G(\mathcal{R}_{L}),$

where $L$ runs over all finite separable extensions of $\kappa(\\!(t)\\!)$, we
thus find a finite separable extension $L$ of $\kappa(\\!(t)\\!)$ such that
$b\in G(\mathcal{R}_{L})$.

For any $(V,\rho)\in\operatorname{\mathbf{Rep}}_{F}(G)$, it follows from the
construction of $b$ that the automorphism $\rho(b^{-1})\colon
V_{\mathcal{B}_{0}}\to V_{\mathcal{B}_{0}}$ factors through
$(V_{\mathcal{B}_{0}})^{\Theta_{X_{0}}=0}\otimes\mathcal{B}_{0}$. Notice that
$\Theta_{X_{0}}$ and $X_{0}$ agree on
$\omega_{1}(V)=V_{K^{\operatorname{nr}}}$. Therefore, we have

$\displaystyle\rho(b)X_{0}\rho(b^{-1})-\partial(\rho(b))\rho(b^{-1})=0.$ (10)

We now fix a faithful representation $(V,\rho)$. The equality (10) then
implies

$\Gamma_{b}(X_{0})=0.$

Put $X_{1}\coloneqq X-X_{0}\in\mathfrak{g}_{\mathcal{R}}$, we then have

$\displaystyle\Gamma_{b}(X)$
$\displaystyle=\operatorname{Ad}(b)(X_{0}+X_{1})-\operatorname{dlog}(b)$
$\displaystyle=\operatorname{Ad}(b)(X_{0})-\operatorname{dlog}(b)+\operatorname{Ad}(b)(X_{1})$
$\displaystyle=\Gamma_{b}(X_{0})+\operatorname{Ad}(b)(X_{1})$
$\displaystyle=\operatorname{Ad}(b)(X_{1}).$

Conserving the notation as in the second paragraph,
$\Theta_{X}=\rho(b)X_{1}\rho(b^{-1})$ acts in the basis
$\mathbf{w}_{1}^{(1)},\cdots,\mathbf{w}_{d_{1}}^{(1)},\cdots,\mathbf{w}_{1}^{(l)},\cdots,\mathbf{w}_{d_{l}}^{(l)}$
via a matrix of the form

$\footnotesize\begin{pmatrix}0&&\\\
&&0&~{}~{}~{}~{}\mathop{\scalebox{1.9}{\raisebox{3.68558pt}{$\ast$}}}\\\
&&&\ddots\\\ &&&&0\end{pmatrix}\in\operatorname{Mat}_{d,d}(\mathcal{R}_{L}).$

Here, the $i$-th $0$ in the diagonal denotes the zero matrix of size
$d_{i}\times d_{i}$. Hence,
$\Gamma_{b}(X)\in\operatorname{Lie}\big{(}U_{G_{\mathcal{R}_{L}}}(-\lambda_{g,\mathcal{R}_{L}})\big{)}=\operatorname{Lie}\big{(}U_{G_{\mathcal{R}}}(-\lambda_{g})_{\mathcal{R}_{L}}\big{)}=\operatorname{Lie}\big{(}U_{G_{\mathcal{R}}}(-\lambda_{g})\big{)}_{\mathcal{R}_{L}}$,
as desired.

∎

## References

  * [1] Y. André.  Filtrations de type Hasse-Arf et monodromie $p$-adique.  Invent. Math., 148:285–317, 2002.
  * [2] J. Anschütz.  Reductive group schemes over the fargues-fontaine curve.  Math. Ann., 374:1219–1260, 2019.
  * [3] L. Berger.  Représentations $p$-adiques et équations différentielles.  Invent. Math., 148:219–284, 2002.
  * [4] B. Conrad, O. Gabber, and G. Prasad.  Pseudo-reductive groups, volume 26 of New mathematical monographs.  Cambridge University Press, 2 edition, 2015.
  * [5] J.-F. Dat, S. Orlik, and M. Rapoport.  Period domains over finite and $p$-adic fields, volume 183 of Cambridge Tracts in Mathematics.  Cambridge University Press, 2010.
  * [6] P. Deligne and J. Milne.  Tannakian categories.  In Hodge cycles, motives, and Shimura varieties, volume 900 of Lecture Notes in Mathematics, pages 101–228. Springer-Verlag, 1982.
  * [7] M. Demazure and P. Gabriel.  Groupes algébriques. Tome I: Géométrie algébrique, généralités, groupes commutatifs.  Masson, 1970.
  * [8] K. S. Kedlaya.  A $p$-adic local monodromy theorem.  Ann. Math., 160:93–184, 2004.
  * [9] K. S. Kedlaya.  Slope filtrations revisited.  Doc. Math., 10:447–525, 2005.
  * [10] K. S. Kedlaya.  Slope filtrations for relative Frobenius.  Astérisque, 317:259–301, 2008.
  * [11] K. S. Kedlaya.  $p$-adic differential equations, volume 125 of Cambridge Studies in Advanced Mathematics.  Cambridge University Press, 2010.
  * [12] K. S. Kedlaya.  Notes on isocrystals.  2016\.  arXiv: 1606.01321.
  * [13] R. Kottwitz.  Isocrystals with additional structures.  Compos. Math., 56:201–220, 1985.
  * [14] R. Liu.  Slope filtrations in families.  J. Inst. Math. Jussieu, 12:249–296, 2013.
  * [15] A. Marmora.  Facteurs epsilon $p$-adiques.  Compos. Math., 144:439–483, 2008.
  * [16] Z. Mebkhout.  Analogue $p$-adique du théorème de turrittin et le théorème de la monodromie $p$-adique.  Invent. Math., 148:319–351, 2002.
  * [17] J. Milne.  Algebraic groups, volume 170 of Cambridge Studies in Advanced Mathematics.  Cambridge University Press, 2017.
  * [18] J.-P. Serre.  Galois cohomology.  Springer Monographs in Mathematics. Springer-Verlag, corrected second printing edition, 2002.
  * [19] T. A. Springer.  Linear algebraic groups.  Modern Birkhäuser Classics. Springer-Verlag, 2 edition, 1998.
  * [20] R. Steinberg.  Regular elements of semisimple algebraic groups.  Inst. Hautes Études Sci. Publ. Math., 25:49–80, 1965.
  * [21] N. Tsuzuki.  Finite local monodromy of overconvergent unit-root $F$-isocrystals on a curve.  Amer. J. Math., 120:1165–1190, 1998.
  * [22] P. Ziegler.  Graded and filtered fiber functors on tannakian categories.  J. Inst. Math. Jussieu, 14:87–130, 2015.

Department of Mathematics, East China Normal University,
---
500 Dongchuan Road, Shanghai, 200241 P.R. China
E-mail address<EMAIL_ADDRESS>