# UAV-Assisted Multi-Cluster Over-the-Air Computation

Min Fu, , Yong Zhou, , Yuanming Shi, , Chunxiao Jiang, , and Wei Zhang M. Fu,
Y. Zhou, and Y. Shi are with School of Information Science and Technology,
ShanghaiTech University, Shanghai 201210, China (e-mail: {fumin, zhouyong,
shiym}@shanghaitech.edu.cn).C. Jiang is with the Tsinghua Space Center and the
Beijing National Research Center for Information Science and Technology,
Tsinghua University, Beijing 100084, China (e-mail: jchx@tsinghua.edu.cn).W.
Zhang is with the School of Electrical Engineering and Telecommunications, The
University of New South Wales, Sydney, NSW 2052, Australia (e-mail:
w.zhang@unsw.edu.au).

###### Abstract

In this paper, we study unmanned aerial vehicles (UAVs) assisted wireless data
aggregation (WDA) in multi-cluster networks, where multiple UAVs
simultaneously perform different WDA tasks via over-the-air computation
(AirComp) without terrestrial base stations. This work focuses on maximizing
the minimum amount of WDA tasks performed among all clusters by optimizing the
UAV’s trajectory and transceiver design as well as cluster scheduling and
association, while considering the WDA accuracy requirement. Such a joint
design is critical for interference management in multi-cluster AirComp
networks, via enhancing the signal quality between each UAV and its associated
cluster for signal alignment and meanwhile reducing the inter-cluster
interference between each UAV and its non-associated clusters. Although it is
generally challenging to optimally solve the formulated non-convex mixed-
integer nonlinear programming, an efficient iterative algorithm as a
compromise approach is developed by exploiting bisection and block coordinate
descent methods, yielding an optimal transceiver solution in each iteration.
The optimal binary variables and a suboptimal trajectory are obtained by using
the dual method and successive convex approximation, respectively. Simulations
show the considerable performance gains of the proposed design over benchmarks
and the superiority of deploying multiple UAVs in increasing the number of
performed tasks while reducing access delays.

###### Index Terms:

Over-the-air computation, UAV communications, wireless data aggregation,
multi-cluster cooperation, interference management.

##  I Introduction

Machine-type communication (MTC) is one of the disruptive technologies
promised by 5G and beyond wireless networks [1]. Therein, it is crucial to
collect and leverage Big Data effectively for decision-making to automate
various intelligent applications. However, collecting data generated by an
enormous number of devices in the future Internet of Things (IoT) networks is
critically challenging due to limited spectrum resource [2], [3]. Meanwhile,
many emerging IoT applications (e.g., environmental monitoring) only aim to
collect a particular function of these massive data rather than to reconstruct
each individual data, which is referred to as wireless data aggregation (WDA)
[4]. To meet these demands, over-the-air computation (AirComp) has recently
been considered as an attractive technique to enable fast WDA among massive
devices by seamlessly integrating communication and computation processes [5].
The principle of AirComp is to exploit the waveform/signal superposition
property of multiple-access channels (MAC) such that an edge server directly
receives a function of concurrently transmitted data. This results in low
transmission delays regardless of the amount of devices, and makes AirComp
particularly appealing to data-intensive and/or latency-critical applications
such as consensus control[6], distributed sensing [7], and distributed machine
learning[8],[9].

So far, AirComp has been studied from various aspects in single-cell networks,
such as single-input-single-output (SISO) AirComp [7, 10], multiple-input-
single-output (MISO) AirComp [11, 12], and multiple-input-multiple-output
(MIMO) AirComp [13]. In single-cell networks, to achieve accurate computing,
AirComp requires the phase and magnitude of all signals to be aligned at the
receiver side. However, channel heterogeneity across devices makes signal
alignment challenging. To cope with this issue, different transceiver designs
have been proposed to compensate for the non-uniform channel fading and
suppress the noise. Specifically, for SISO AirComp, the authors in [7, 10]
proposed the optimal transmit power control and receive normalizing factor
design. For multi-antenna AirComp systems, beamforming vectors at receiver
and/or transmitter were designed to minimize the computation error in [11,
13]. However, the AirComp performance may deteriorate when one or more power-
constrained devices are in deep fading. To mitigate the communication
bottleneck, the authors in [12] employed a promising reconfigurable
intelligent surface (RIS) [14, 15] to jointly design the passive beamforming
at the RIS and transceiver. In addition, we in [16], [17] proposed to deploy
an unmanned aerial vehicle (UAV) as a mobile server and to exploit its
mobility to avoid any device being in deep fade, thereby enhancing the
performance of AirComp.

Meanwhile, the IoT networks generally involve different WDA tasks, each of
which is characterized by their applications (e.g., classification FL task
[18], regression FL task, and sensing tasks [19]), data types (e.g., model
parameters in machine learning and velocity in connected car platooning
applications), and computing functions (e.g., sum and mean). Hence,
researchers recently advocated the study of AirComp in a multi-cluster network
to simultaneously complete multiple WDA tasks, which is referred to as multi-
cluster AirComp [4]. In multi-cluster AirComp, besides signal alignment and
noise suppression, inter-cluster interference management arises as a new
challenge. Different from the orthogonal multiple access (OMA) based multi-
cluster networks [20], each server in the multi-cluster AirComp network not
only harnesses intra-cluster interference for function computation, but also
suppresses inter-cluster interference for computation error reduction. As an
initial study, the authors in [19] proposed a signal-and-interference
alignment scheme to simultaneously eliminate inter-cell interference and align
intra-cell signals in a two-cell MIMO AirComp network. Subsequently, the
authors in [21] studied the weighted sum mean-square error (MSE) minimization
problem in multi-cell SISO AirComp networks by optimizing the transmit power
of devices. Results in [21] showed that, the performance of AirComp is limited
by the inter-cell interference since the transmit power of the device needs to
balance the tradeoff between combating inter-cell interference and enhancing
signal alignment within the cell. These results are established when the
static terrestrial BS is available. However, in remote and under-developed
areas, the terrestrial BSs are usually sparsely deployed or not available. In
these harsh circumstances, it is critical to deploy more flexible BSs to
unleash the potential of multi-cluster AirComp.

As a parallel but complementary study, in this paper, we investigate a novel
multi-cluster AirComp framework with multiple UAVs dispatched as flying BSs to
cooperatively perform diverse AirComp tasks, where no terrestrial BS is
available. In fact, multiple UAVs have been deployed as BSs to assist the
terrestrial networks for rate-oriented communications [22], [23], [24]. For
example, the authors in [22] studied multi-UAV cooperative communications to
achieve higher data rate and lower access delay. The authors in [23]
considered a multi-UAV enabled wireless network for data collection. This
motivates us to study the use of multiple UAVs in multi-cluster AirComp by
taking into account UAV trajectory planning as well as cluster scheduling and
association for interference management, which has the following advantages.
First, multi-UAV cooperation exploits spectrum sharing to allow multiple
clusters to be simultaneously served, thereby increasing the amount of
performed task within the given time duration. Second, scheduling the clusters
that are far from each other can avoid strong co-channel interference and thus
improve computation accuracy. Furthermore, the joint design not only shortens
the communication distances between the UAV and its associated cluster to
enhance intra-cluster signal alignment but also enlarges the communication
distances between the UAV and its non-associated clusters to rein in the
inter-cluster interference. This joint design is crucial for interference
management in multi-cluster AirComp networks, but has not yet been studied in
the literature to the best of the authors’ knowledge.

A widely adopted performance metric for quantifying the AirComp computation
error is MSE between the estimated function value and the ground truth[8, 13,
10]. To promote fairness among clusters, under the given target MSE
requirements, we aim to maximize the minimum amount of performed WDA tasks
among all clusters by jointly designing cluster scheduling and association,
UAVs’ trajectories, and transceiver design in a given time duration. However,
such a non-convex max-min fairness problem presents unique challenges due to
their discontinuous objective functions (as a result of the binary cluster
scheduling and association variables), non-convex MSE constraints (because of
the coupling between all optimization variables), and non-convex trajectory
design. It is generally challenging to optimally solve such a mixed-integer
non-convex optimization problem.

###  I-A Contributions

The main contributions of this paper are summarized as follows.

  * •

This paper is one of the early attempts to study the multi-UAV enabled AirComp
in multi-cluster wireless networks, where multiple UAVs are deployed to
cooperatively perform different AirComp tasks. To ensure fairness among
clusters, we aim to maximize the minimum amount of performed WDA tasks among
all clusters by jointly optimizing scheduling and association for clusters,
UAVs’ trajectory planning, and transceiver design, taking into account the
target accuracy of AirComp, practical constraints on UAVs, inter-cluster
interference, as well as the total power budgets at devices.

  * •

To render the resulting non-convex mixed-integer nonlinear programming (MINLP)
tractable, by leveraging the bisection method, the original problem reduces to
a sequence of the minimum ratio maximization problems, which enables the
development of an iterative algorithm. After adopting block coordinate descent
(BCD) [25], besides both convex normalizing factors and power optimization
subproblems being optimally solved, the optimal cluster scheduling and
association are also obtained by applying the low-complexity Lagrange duality
method. For the non-convex UAV trajectory optimization problem, a suboptimal
solution is obtained by using the successive convex approximation (SCA)
method.

  * •

Simulation results are presented to show the effectiveness and superiority of
the proposed design and developed algorithm. In multi-cluster AirComp with a
single UAV, the performance achieved by the joint design outperforms other
benchmarks and can reach the upper bound. It is also shown that the max-min
task amount of the considered UAV network increases with the mission duration,
revealing a performance-access delay tradeoff in multi-cluster AirComp.
Compared to the single-UAV case, the use of multiple UAVs with effective
cooperative interference management can considerably increase the amount of
tasks performed by each cluster while reducing access delays.

###  I-B Organization

The rest of this paper is organized as follows. We presents the system model
and the problem formulation for multi-cluster AirComp assisted by UAVs in
Section II. We develop an iterative algorithm yielding high-quality solutions
to solve the formulated problem based on the Bisection method, BCD, and SCA in
Section III. In addition, numerical results are presented in Section IV to
evaluate the performance of the proposed design. In Section V, we draw the
conclusions.

Notations: Scalars, column vectors, and matrices are written in italic
letters, boldfaced lower-case letters, and boldfaced upper-case letters
respectively, e.g., $a$, $\mathbf{a}$, $\mathbf{A}$. $\mathbb{R}^{M\times N}$
denotes the space of a real-valued matrix with with $M$ rows and $N$ columns.
$\left\|\mathbf{a}\right\|_{2}$ denotes the Euclidean norm of vector
$\mathbf{a}$ and $\mathbf{a}^{T}$ represents its transpose. $|\mathcal{S}|$
denotes the cardinality of the set $\mathcal{S}$.

##  II System Model And Problem Formulation

Figure 1: Illustration of multi-UAV assisted AirComp in multi-cluster wireless
networks.

###  II-A System Model

As shown in Fig. 1, we consider a multi-cluster wireless network, where each
cluster consists of multiple ground devices and multiple UAVs111We consider
the scenario that the terrestrial BS is not available, e.g., in wild areas.
Under this circumstance, the UAVs are adopted as an alternative to provide
wireless services for ground devices [22].  are deployed to support multiple
WDA tasks using AirComp. Due to the limited physical size, we assume that UAVs
and ground devices are all equipped with one antenna. Let
$\mathcal{M}\triangleq\\{1,\ldots,M\\}$ and
$\mathcal{L}\triangleq\\{1,2,\ldots,L\\}$ denote the sets of $M$ UAVs and $L$
clusters, respectively. We denote
$\mathcal{K}_{l}\triangleq\\{1,\ldots,K_{l}\\}$ as the set of $K_{l}$ ground
devices in cluster $l$ with
$\mathcal{K}_{i}\cap\mathcal{K}_{j}=\emptyset,\forall i\neq
j,i,j\in\mathcal{L}$, and
$\mathcal{K}\triangleq\cup_{l\in\mathcal{L}}\mathcal{K}_{l}$ as the set of all
ground devices. Therein, the ground devices in set $\mathcal{K}_{l}$ need to
collaboratively perform the type-$l$ WDA task, i.e., computing a specific
nomographic function (e.g., arithmetic mean) with respect to their measured
data, and transmit the aggregated data to one of the UAVs. To avoid excessive
communication latency, we adopt the AirComp technique [8] to achieve fast WDA
by enabling concurrent transmission among multiple ground devices. Compared to
the OMA scheme, the unique feature of AirComp is the utilization of the
waveform superposition property of a MAC to harness the intra-cluster
interference for function computation, thereby reducing the communication
latency. Hereinafter, we refer to this type of WDA task by using the AirComp
technique as an AirComp task.

We assume that the AirComp tasks in different clusters are distinct, i.e.,
different classes of IoT applications supported by different clusters (e.g.,
distributed machine learning [8], [9] and control consensus [6]). Among them,
their data types (e.g., model parameters in machine learning and velocity in
connected car platooning applications) and the intended computing functions
(e.g., sum and mean) can be different. For UAV deployment, it is not practical
to dispatch one or more UAVs for each cluster, especially when the number of
clusters served is relatively large. Thus, in this paper, we assume that
$M\leq L$ and each UAV can select an appropriate cluster for association.
Different UAVs can simultaneously deal with different AirComp tasks over the
same frequency band, thereby reducing the access delay of ground devices in
all clusters. As a result, this leads to inter-cluster interference among $M$
clusters and calls for interference management to balance the computing errors
among different clusters.

We aim to optimize the trajectories of UAVs to improve the WDA accuracy of all
clusters. This can be achieved by jointly designing the UAVs’ trajectories and
UAV-cluster association to shorten the communication distances between the UAV
and its associated cluster to enhance intra-cluster signal alignment, and
enlarge the communication distances between the UAV and its non-associated
clusters to rein in the inter-cluster interference. To facilitate trajectory
design, in this paper, time discretization technique is adopted to divide the
mission duration $T$ (s) into $N$ equal time-slots with time step-size
$\delta=\frac{T}{N}$, as in most of the previous studies on UAV-enabled
wireless networks [22, 23, 26]. The set of $N$ time slots is denoted as
$\mathcal{N}\triangleq\\{1,\ldots,N\\}$. Note that the time interval $\delta$
needs to be appropriately set, so that in each time slot, each UAV can
complete one AirComp task while the changes in the communication distances
between the UAVs and the ground devices are negligible.

In a three dimensional (3D) Cartesian coordinate system, we assume that all
UAVs has a fixed altitude of $H$ (m) over the horizontal plane, where $H$ can
be set to the minimum altitude to ensure that the UAVs can avoid obstacles
(e.g., building or terrain) without the need for frequent ascent and descent.
We denote, at time slot $n$, the horizontal coordinate of UAV $m$ as
$\bm{q}_{m}[n]\triangleq[x_{m}[n],y_{m}[n]]^{\sf T}\in\mathbb{R}^{2\times 1}$.
Moreover, the fixed horizontal coordinate of each ground device
$k\in\mathcal{K}$ is denoted as
$\mathbf{w}_{k}=[x_{k},y_{k}]^{T}\in\mathbb{R}^{2\times 1}$, which is assumed
to be known at the UAVs.

####  II-A1 Channel Model

A large body of studies on air-to-ground channel modeling [27], [28]. And
recent field experiments by Qualcomm [29] has verified that when UAVs fly
above a certain altitude, there is a high probability for the UAV-to-ground
channel to be dominated by a line-of-sight (LoS) link. Therefore, the uplink
channels from devices to UAVs are assumed to be dominated by LoS links in this
paper. The channel between UAV $m\in\mathcal{M}$ and device $k\in\mathcal{K}$
at time slot $n\in\mathcal{N}$ is modeled as

$\displaystyle
h_{k,m}[n]=\sqrt{\beta_{0}d^{-\gamma}_{k,m}[n]}e^{-j\theta_{k,m}[n]},$ (1)

where $\beta_{0}$ denotes the channel power gain at the reference distance of
$d_{0}=1$ m, $\gamma\geq 2$ is the path loss exponent, and $\theta_{k,m}[n]$
is phase shift component. We assume that all devices can perfectly estimate
its own phase component with its associated UAV. And $d_{k,m}[n]$ denotes the
distance from UAV $m$ to ground device $k$, which is represented as

$\displaystyle
d_{k,m}[n]=\sqrt{H^{2}+\left\|\bm{q}_{m}[n]-\mathbf{w}_{k}\right\|_{2}^{2}}.$
(2)

Note that the Doppler effect is assumed to be perfectly compensated [30].

####  II-A2 Multi-cluster AirComp

We define a set of binary variables
$\\{a_{l,m}[n]|l\in\mathcal{L},m\in\mathcal{M},n\in\mathcal{N}\\}$ to
represent UAV-cluster scheduling and association over different time slots. We
set $a_{l,m}[n]=1$ if there exist uplink transmissions from ground devices in
cluster $l$ to UAV $m$ at time slot $n$, and $a_{l,m}[n]=0$ otherwise.
Therefore, $a_{l,m}[n]$ not only indicates the association status between UAV
$m$ and cluster $l$ at slot $n$, but also determines the communication
scheduling for cluster $l$ at slot $m$. Suppose that, in each time slot, the
ground devices in each cluster can communicate with at most one UAV, and each
UAV can serve at most one cluster. In summary, the binary variables impose the
following constraints.

$\displaystyle\sum_{m=1}^{M}a_{l,m}[n]\leq 1,\forall
l\in\mathcal{L},n\in\mathcal{N},$ (3)
$\displaystyle\sum_{l=1}^{L}a_{l,m}[n]\leq 1,\forall
m\in\mathcal{M},n\in\mathcal{N},$ (4) $\displaystyle
a_{l,m}[n]\in\\{0,1\\},\forall
l\in\mathcal{L},m\in\mathcal{M},n\in\mathcal{N}.$ (5)

We consider the case that each cluster aims to compute the average of
distributed data measured by the ground devices [8], [10]. In the following,
we take cluster $l$ associated with UAV $m$ as an example. Let
$z_{k}[n]\in\mathbb{C}$ denote the measured data from device
$k\in\mathcal{K}_{l}$ in cluster $l$ at time slot $n$. We denote the target
function (e.g., arithmetic mean) of $K_{l}$ variables as
$f_{l}[n]:\mathbb{C}^{K_{l}}\rightarrow\mathbb{C}$. Therefore, the target
average function of type-$l$ data at time slot $n$ is expressed as

$\displaystyle{f}_{l}[n]=\frac{1}{K_{l}}\sum_{k\in\mathcal{K}_{l}}s_{k}[n],$
(6)

where $s_{k}[n]\triangleq\psi_{k}(z_{k}[n]),\forall k\in\mathcal{K}_{l}$ with
$\psi_{k}$ denoting the pre-processing function at device
$k\in\mathcal{K}_{l}$. We assume that $\\{s_{k}[n]\\}_{k\in\mathcal{K}_{l}}$
are independent and have zero mean and unit variance as in [10],[7], i.e.,
$\mathbb{E}(s_{k}[n])=0$, $\mathbb{E}(s_{k}[n]s_{k}^{\sf H}[n])=1$, and
$\mathbb{E}[s_{i}[n]s_{j}[n]^{\sf H}]=0,\forall i\neq j$.

In the sequel, we describe how to estimate the desired function $f_{l}[n]$ in
(6) at UAV $m$ via AirComp. Let
$\mathcal{J}(n)=\\{l|\sum_{m\in\mathcal{M}}a_{l,m}[n]=1,\forall
l\in\mathcal{L}\\}$ denote the index set of active clusters at time slot $n$.
After all devices in set $\mathcal{K}_{l}$ simultaneously send their pre-
processed signals $\\{s_{k}[n]\\}$ to UAV $m$ over the same radio channel, the
received signal at UAV $m$ is written as

$\displaystyle
y_{m}[n]=\sum_{k\in\mathcal{K}_{l}}b_{k}[n]h_{k,m}[n]s_{k}[n]+\sum_{j\in\mathcal{J}(n)\backslash\\{l\\}}\sum_{i\in\mathcal{K}_{j}}h_{i,m}[n]b_{i}[n]s_{i}[n]+e_{m}[n],$
(7)

where $b_{k}[n]\in\mathbb{C}$ denotes the transmit precoding coefficient at
device $k$ for channel-fading compensation at time slot $n$, and $e_{m}[n]$
denotes the additive white Gaussian noise (AWGN) at UAV $m$, i.e.,
$e[n]\thicksim\mathcal{C}\mathcal{N}$$(0,\sigma^{2})$. The ground devices’
transmissions are assumed to be synchronized [7, 10]. It is worth mentioning
that the synchronization techniques reported in the literature can be applied
in our work, such as the AirShare technique [31] and the timing advance
mechanism used in the long term evolution systems[4].

Upon receiving signal $y_{m}[n]$, after post-processing and scaling at UAV,
the estimated average function $m$ is given by

$\displaystyle\hat{f}_{l,m}[n]=\frac{\eta_{l,m}[n]y_{m}[n]}{K_{l}},$ (8)

where factor $1/K_{l}$ is employed for averaging purpose and
$\eta_{l,m}[n]\in\mathbb{C}$ is a normalizing factor at UAV $m$ applied to
scale received signal $y_{m}[n]$ for compensating channel fading and
suppressing noise, thereby accurately estimating the target function
$f_{l}[n]$. It can be observed from (7) and (8) that UAV $m$ harnesses the
intra-cluster interference for computing the target function but suffers from
inter-cluster interference due to the communication resource reuse among UAVs.
Each device $k$ has a total transmit power budget $P_{k}$, which is given by

$\displaystyle\sum_{n=1}^{N}|b_{k}[n]|^{2}\leq P_{k},\forall k,$ (9)

where the total power constraint at each ground device is considered for
adaptive power allocation over different time slots.

###  II-B Performance Metrics

To quantify the performance of AirComp, the computation distortion of the
ground true average function $f_{l}[n]$ is measured by its MSE. Specifically,
when $a_{l,m}[m]=1$, the corresponding instantaneous MSE of type-$l$ AirComp
task at UAV $m$ at time slot $n$ is given by

$\displaystyle{\sf{mse}}_{l,m,n}$ $\displaystyle=$
$\displaystyle\mathbb{E}[|\hat{f}_{l,m}[n]-f_{l}[n]|^{2}]=\frac{1}{K_{l}^{2}}\mathbb{E}\Big{[}\Big{|}{\eta_{l,m}[n]}{y_{m}[n]}-\sum_{k\in\mathcal{K}_{l}}s_{k}[n]\Big{|}^{2}\Big{]}$
(10) $\displaystyle=$
$\displaystyle\dfrac{1}{K_{l}^{2}}\Big{[}\sum_{k\in\mathcal{K}_{l}}\Big{(}\eta_{l,m}[n]b_{k}[n]h_{k,m}[n]\\!-1\Big{)}^{2}+\eta_{l,m}^{2}[n]\Big{(}\hat{I}_{l,m}[n]+\sigma^{2}\Big{)}\Big{]},$

where
$\hat{I}_{l,m}[n]\triangleq\sum_{j\in\mathcal{J}(n)\backslash\\{l\\}}\sum_{i\in\mathcal{K}_{j}}|b_{i}[n]|^{2}|h_{i,m}[n]|^{2}$
represents the inter-cluster interference received at UAV $m$ when it
associates with cluster $l$ at time slot $n$. The expectation in (10) is taken
over the distributions of signals $\\{s_{k}[n]\\}$ and noise $e_{m}[n]$. The
instantaneous MSE of type-$l$ AirComp task at time slot $n$ is represented as

$\displaystyle{\sf{MSE}}_{l,n}=\sum_{m=1}^{M}a_{l,m}[n]{\sf mse}_{l,m,n}.$
(11)

Before formulating the optimization problem, we first simplify the transmit
precoding coefficients and binary variable design. Note that the binary
variables in (11) are coupled across different clusters and UAVs due to the
interference term. The interference terms in ${\sf mse}_{l,m,n}$ are
simplified as

$\displaystyle
I_{l,m}[n]\triangleq\sum_{j\in\mathcal{L}\backslash\\{l\\}}\sum_{i\in\mathcal{K}_{j}}|b_{i}[n]|^{2}|h_{i,m}[n]|^{2}.$
(12)

This is because if there is no UAV associated with cluster $j$, i.e,
$\sum_{m\in\mathcal{M}}a_{j,m}[n]=0$, then the transmission power of device
$i\in\mathcal{K}_{j}$ must be zero (i.e., $b_{i}[n]=0$). Furthermore, it can
also be observed from (10) that the phases of $b_{i}[n]$ and $h_{i,m}[n]$ for
$\forall j\in\mathcal{L}\backslash\\{l\\},\forall i\in\mathcal{K}_{j}$ do not
affect the inter-cluster interference induced error in ${\sf{mse}}_{l,m,n}$
while the essential condition for ${\sf{mse}}_{l,m,n}$ to reach the minimum
value is the terms $\eta_{l,m}[n]b_{k}[n]h_{k,m}[n],\forall
k\in\mathcal{K}_{l}$ in (10) must be real and non-negative. Therefore, with
$\eqref{Cons:A UAV}$, when $a_{l,m}[n]=1$, we have

$\displaystyle{\sf{MSE}}_{l,n}={\sf
mse}_{l,m,n}\geq\dfrac{1}{K_{l}^{2}}\Big{[}\sum_{k\in\mathcal{K}_{l}}\Big{(}\eta_{l,m}[n]\sqrt{p_{k}[n]}|h_{k,m}[n]|-1\Big{)}^{2}+\eta_{l,m}^{2}[n]\Big{(}I_{l,m}[n]+\sigma^{2}\Big{)}\Big{]},$
(13)

where the equality holds when $\eta_{l,m}[n]\in\mathbb{R}^{+}$ and
$b_{k}[n]\triangleq\frac{\sqrt{p_{k}[n]}h_{k,m}^{\dagger}[n]}{|h_{k,m}[n]|}$
with $p_{k}[n]\in\mathbb{R^{+}}$ denoting the transmit power of device $k$ at
time slot $n$. please refer to [16] for its proof. Meanwhile, if
$\sum_{m=1}^{M}a_{l,m}[n]=0$, we have ${\sf{MSE}}_{l,n}=0$ and
$b_{k}[n]=0,\forall k\in\mathcal{K}_{l}$. Hence, at time slot $n$, the
achievable instantaneous MSE of type-$l$ data can be represented as

$\displaystyle{\sf{MSE}}_{l,n}$ $\displaystyle=$
$\displaystyle\sum_{m=1}^{M}a_{l,m}[n]{\sf mse}_{l,m,n}$ (14) $\displaystyle=$
$\displaystyle\sum_{m=1}^{M}\frac{a_{l,m}[n]}{K_{l}^{2}}\Big{[}\sum_{k\in\mathcal{K}_{l}}\Big{(}\eta_{l,m}[n]|h_{k,m}[n]|\sqrt{p_{k}[n]}-1\Big{)}^{2}+\eta_{l,m}^{2}[n]\Big{(}I_{l,m}[n]+\sigma^{2}\Big{)}\Big{]},$

which is a linear function of $a_{l,m}[n],\forall m$ but not coupled over
different clusters, and is related to real transmit power $p_{k}[n]$ instead
of complex transmit precoding coefficient $b_{k}[n]$.

######  Remark 1.

One can observe that ${\sf{MSE}}_{l,n}$ is composed of three components,
including the signal misalignment error (i.e.,
$\sum_{m=1}^{M}\frac{a_{l,m}[n]}{K_{l}^{2}}\sum_{k\in\mathcal{K}_{l}}\big{(}\eta_{l,m}[n]|h_{k,m}[n]|\sqrt{p_{k}[n]}-1\big{)}^{2}$),
the inter-cluster interference-induced error (i.e.,
$\sum_{m=1}^{M}\frac{a_{l,m}[n]}{K_{l}^{2}}\eta_{l,m}^{2}[n]\sum_{j\in\mathcal{L}\backslash\\{l\\}}\sum_{i\in\mathcal{K}_{j}}p_{i}[n]|h_{i,m}[n]|^{2}$
), and the noise-induced error (i.e.,
$\sum_{m=1}^{M}\frac{a_{l,m}[n]}{K_{l}^{2}}\eta_{l,m}^{2}[n]\sigma^{2}$). The
device can increase its power $p_{k}[n]$ to compensate for channel fading for
signal alignment whereas it leads to an increased interference for other co-
channel clusters. Meanwhile, enlarging $\eta_{l,m}[n]$ can reduce the signal-
misalignment error but magnify the negative impact of noise and inter-cluster
interference. Fortunately, the trajectory design with cluster scheduling and
association can be exploited to balance the tradeoff between signal
misalignment error reduction and co-channel interference reduction.
Specifically, a UAV moves closer to the associated cluster $l$ to construct
strong desired links $|h_{k,m}[n]|,\forall k\in\mathcal{K}_{l}$ for signal
alignment and keeps away from other scheduled clusters to construct weak co-
channel interference links $|h_{i,m}[n]|^{2},i\in\mathcal{K}_{j},j\neq l$.

###  II-C Problem Formulation

In the sequel, let $\bm{A}=\\{a_{l,m}[n],\forall l\in\mathcal{L},\forall
m\in\mathcal{M},\forall n\in\mathcal{N}\\}$, $\bm{P}=\\{p_{k}[n],\forall
k\in\mathcal{K},\forall n\in\mathcal{N}\\}$,
$\bm{\eta}=\\{\eta_{l,m}[n],\forall l\in\mathcal{L},\forall
m\in\mathcal{M},\forall n\in\mathcal{N}\\}$, and
$\bm{Q}=\\{\bm{q}_{m}[n],\forall m\in\mathcal{M},\forall n\\}$. Given mission
duration $T$ and certain instantaneous MSE requirements for different AirComp
tasks, we study the total amount of AirComp tasks maximization problem while
ensuring fairness among clusters. Therefore, we maximize the minimum amount of
AirComp tasks among all clusters, i.e.,
$\min_{l}\sum_{n=1}^{N}\sum_{m=1}^{M}a_{l,m}[n]$, by jointly optimizing the
cluster scheduling and association indicator variables $\bm{A}$, transmit
power $\bm{P}$ at devices, signal normalizing factors $\bm{\eta}$ at the UAVs,
and UAVs’ trajectories $\bm{Q}$. The optimization problem is formulated as

$\displaystyle\mathscr{P}:\mathop{\text{maximize}}\limits_{\bm{A},\bm{P},\bm{\eta},\bm{Q}}\mathop{\text{min}}\limits_{l}$
$\displaystyle\sum_{n=1}^{N}\sum_{m=1}^{M}a_{l,m}[n]$ subject to
$\displaystyle{\sf{MSE}}_{l,n}\leq\varepsilon_{l},\forall n,\forall l,$
$\displaystyle 0\leq\sum_{n=1}^{N}p_{k}[n]\leq P_{k},\forall k,$
$\displaystyle\bm{q}_{m}[0]=\bm{q}_{m}[N],\forall m\in\mathcal{M},$
$\displaystyle\left\|\bm{q}_{m}[n]-\bm{q}_{m}[n-1]\right\|_{2}\leq
V_{\text{max}}\delta,\forall n\in\mathcal{N},\forall m\in\mathcal{M},$
$\displaystyle\left\|\bm{q}_{m}[n]-\bm{q}_{i}[n]\right\|_{2}\geq
d_{\text{min}},\forall m\neq i,i,m\in\mathcal{M},\forall n,$ $\displaystyle
0\leq\eta_{m}[n],\forall m,\forall n,$ $\displaystyle\text{Constraints\
}\eqref{Cons:A UAV},\eqref{Cons:A cluster},\eqref{Cons:A binary}.$

Constraints (15) describe the minimum instantaneous MSE requirement
$\varepsilon_{l}$ for data aggregation from any cluster $l$ in each time slot.
To avoid trivial problems, it is assumed that $\varepsilon_{l}<1,\forall l$.
Otherwise, a trivial solution of problem (15), i.e, $p_{k}[n]=0,\forall
k\in\mathcal{K},\forall n\in\mathcal{N}$ and $\eta_{l,m}[n]=0,\forall
l\in\mathcal{L},\forall m\in\mathcal{M},\forall n\in\mathcal{N}$ achieves the
certain instantaneous MSE requirements for any given $\bm{Q}$ and $\bm{A}$.
Constraints (15) are the maximum total transmit power constraints of all
devices. Constraints (15) represent the UAVs’ initial and final locations
constraints. Constraints (15) represent that the UAVs are constrained by the
maximum speed $V_{\max}$ in meter/second (m/s). Constraints (15) limit the
minimum safe distance between UAVs for collision avoidance with $d_{\min}$
denoting the minimum inter-UAV distance in meters. Constraints (15) are
derived from (13). It can be easily verified that for problem $\mathscr{P}$,
if the obtained binary variables meet condition $\sum_{m=1}^{M}a_{l,m}[n]=0$,
the power solution for all devices in cluster $l$ must satisfy
$p_{k}[n]=0,\forall k\in\mathcal{K}_{l}$.

The main challenges in solving problem (15) arise from the following aspects.
First, due to the binary cluster scheduling and association variables
$\bm{A}$, max-min problem (15) involves a discontinuous objective function and
integer constraints (15), (4), and (3). Second, all continuous variables
$\\{\bm{P},\bm{\eta},\bm{Q}\\}$ are coupled with binary variables $\bm{A}$ in
constraints (15), which are non-convex. In addition, collision avoidance
constraints (15) are also non-convex. As such, problem (15) is a nonconvex
max-min problem with MINLP, which is generally challenging to be optimally
solved. As a compromise approach, in the following subsection, we transform
problem (15) into an equivalent MINLP problem to facilitate the development of
a low-complexity algorithm.

###  II-D Problem Reformulation

To address the challenges of the max-min problem with the discontinuous
objective function, an auxiliary variable
$\Gamma=\mathop{\text{max}}\limits_{(l,n)}\frac{{\sf{MSE}}_{l,n}}{\varepsilon_{l}}$
is defined as the maximum achievable AirComp MSE ratio among all clusters.
Furthermore, we denote the minimum amount of AirComp tasks among all clusters
as $D\in\mathbb{Z}$. Given any $D$, we first introduce the following
optimization problem:

$\displaystyle\mathscr{P}_{1}:\mathop{\text{minimize}}\limits_{\Gamma,\bm{A},\bm{P},\bm{\eta},\bm{Q}}$
$\displaystyle\Gamma(D)$ subject to
$\displaystyle\frac{{\sf{MSE}}_{l,n}}{\varepsilon_{l}}\leq\Gamma(D),\forall
n,\forall l,$ $\displaystyle\sum_{n=1}^{N}\sum_{m=1}^{M}a_{l,m}[n]\geq
D,\forall l,$ $\displaystyle\text{Constraints\ }\eqref{Cons:A
UAV},\eqref{Cons:A cluster},\eqref{Cons:A
binary},\eqref{Cons:power},\eqref{Cons:normalizing factor},\eqref{Cons:Q
position},\eqref{Cons:Q speed},\eqref{Cons:Q collision}.$

Problem $\mathscr{P}_{1}$ aims to minimize the maximum ratio $\Gamma$, which
is a function of $D$. By denoting the optimal objective value of
$\mathscr{P}_{1}$ as $\Gamma^{\star}(D)$, we have the following proposition.

######  Proposition 1.

The optimal objective value $\Gamma^{\star}(D)$ of problem $\mathscr{P}_{1}$
is a non-decreasing of $D$.

###### Proof.

For any given task amount, i.e., $D$ and $D^{\prime}$, such that
$D^{\prime}>D$, we should prove that
$\Gamma^{\star}(D^{\prime})\geq\Gamma^{\star}(D)$. Given $D^{\prime}$, the
optimal solution to $\mathscr{P}_{1}$ is denoted as
$\\{\Gamma^{\star},\bm{A}^{\star},\bm{P}^{\star},\bm{\eta}^{\star},\bm{Q}^{\star}\\}$.
As the minimum task amount decreases to $D$, a feasible solution to
$\mathscr{P}_{1}$ with the given $D$ denoting as
$\\{\bar{\bm{A}},\bar{\bm{P}},\bar{\bm{\eta}},\bar{\bm{Q}}\\}$ can be
constructed by letting $\bar{\bm{Q}}=\bm{Q}^{\star}$,
$\bar{\bm{\eta}}=\bm{\eta}^{\star}$, $\bar{\bm{A}}=\bm{A}^{\star}$, and
$\bar{\bm{P}}=\bm{P}^{\star}$. In this case, the objective value
$\Gamma^{\star}(D^{\prime})=\bar{\Gamma}(D)$. $\bar{\Gamma}(D)$ is the
objective value for a feasible solution
$\\{\bar{\bm{A}},\bar{\bm{P}},\bar{\bm{\eta}},\bar{\bm{Q}}\\}$ to
$\mathscr{P}_{1}$, which is obviously no less than the objective function for
the optimal solution $\Gamma^{\star}(D)$. Thus, we have
$\Gamma^{\star}(D^{\prime})=\bar{\Gamma}(D)\geq\Gamma^{\star}(D)$, which
completes the proof. ∎

For any given $D$, the MSE requirements $\\{\varepsilon_{l}\\}$ in each
cluster are achievable if and only if $\Gamma^{\star}(D)\leq 1$. Thus, it can
be verified that problem $\mathscr{P}$ is equivalent to

$\displaystyle\mathscr{P}_{2}:\mathop{\text{maximize}}\limits_{D}$
$\displaystyle D$ subject to $\displaystyle\Gamma^{\star}(D)\leq 1.$ (17a)

With such a transformation, we can concentrate on addressing problem
$\mathscr{P}_{2}$ in the remainder of the paper. Fortunately, based on
Proposition 1, problem $\mathscr{P}_{2}$ can be efficiently solved by applying
the bisection search over the minimum amount of AirComp tasks $D$ until the
equality in (17a) holds. As a result, the main difficulty of solving problem
$\mathscr{P}_{2}$ lies in finding an efficient algorithm to solve problem
$\mathscr{P}_{1}$ with any given amount $D$. Compared to the original problem
$\mathscr{P}$, problem $\mathscr{P}_{1}$ eliminates the max-min programming
and discontinuous objective function. Although it still involves coupled
optimization variables and non-convex constraints, such as (5), (16), and
(15), problem $\mathscr{P}_{1}$ promotes the development of an efficient
suboptimal algorithm, which will be elaborated in section III.

######  Remark 2.

The equations of all the constraints in (16) hold for the optimal solution of
problem $\mathscr{P}_{1}$. If any equality of constraints in (16) is not met,
we can always reduce the value of $\sum_{n=1}^{N}\sum_{m=1}^{M}a_{l,m}[n]$
when other variables are fixed, while all other constraints are still
satisfied, and the objective value of $\mathscr{P}_{1}$ remains unchanged.

##  III Proposed Algorithm

In this section, we develop an efficient iterative algorithm to solve problem
$\mathscr{P}_{1}$ based on the principle of BCD and SCA techniques.
Specifically, to tackle the coupling between optimization variables, we adopt
the well-known BCD technique [25] to decouple problem $\mathscr{P}_{1}$ into
four subproblems, which optimize cluster association $\bm{A}$, transmit power
$\bm{P}$, normalizing factors $\bm{\eta}$, and UAVs’ trajectories $\bm{Q}$ in
an iterative manner. Besides the convex transmit power $\bm{P}$ and
normalizing factors $\bm{\eta}$ optimization subproblems, we propose a low-
complexity dual ascent method to obtain a closed-form binary solution for
cluster association optimization subproblem and exploit the SCA to deal with
the non-convexity of UAVs’ trajectories optimization subproblem.

###  III-A Cluster Scheduling and Association Optimization

1: Input: $M$, $N$, $L$, $D$, $\\{{\sf mse}_{l,m,n}/\varepsilon_{l},\forall
l,\forall m,\forall n\\}$.

2: Initialize dual variables $\\{\lambda_{l,n}=1/(KN)\\}$,
$\\{\beta_{l,n}=0\\}$, and $\\{\nu_{m,n}=0\\}$.

3: repeat

4: Update the primal variables $\bm{A}$ and $\bm{\Gamma}$ according to (23)
and (25).

5: Update the dual variables $\bm{\lambda}$, $\bm{\beta}$, and $\bm{\nu}$
according to (26)-(27).

6: until $\bm{\lambda}$, $\bm{\beta}$, and $\bm{\nu}$ converge within a
prescribed accuracy

7: Set $\bm{a}^{\star}$ = $\bm{a}$ and $\Gamma^{\star}$ = $\Gamma$.

8: Output: $\bm{a}^{\star}$ and $\Gamma^{\star}$

Algorithm 1 Dual Method for Problem (18)

For given $\\{\bm{\eta},\bm{Q},\bm{P}\\}$, the cluster scheduling and
association optimization subproblem is reduced to as follows.

$\displaystyle\mathop{\text{minimize}}\limits_{\Gamma,\bm{A}}$
$\displaystyle\Gamma$ subject to
$\displaystyle\sum_{m=1}^{M}a_{l,m}[n]\frac{{\sf{mse}}_{l,m,n}}{\varepsilon_{l}}\leq\Gamma,\forall
l,\forall n,$ $\displaystyle\text{Constraints~{}}\eqref{Cons:A
UAV},\eqref{Cons:A cluster},\eqref{Cons:A binary},\eqref{Cons:target amount}.$

Problem (18) is a binary linear programming optimization problem. Although
such a problem can be directly solved via CVX and solver (e.g., Mosek), it
requires prohibitively high computational complexity due to exhausted search,
which is not practical even for the moderate sizes of $N$ and $L$.
Alternatively, to gain more design insights and reduce the computational cost,
by relaxing integer constraints (5) with $a_{l,m}[n]\in[0,1]$, we exploit the
Lagrange duality to solve problem (18) with relaxed constraints. It will show
that the obtained closed-form optimal solution ensures not only the dual
optimality but also the feasibility of problem (18). Specifically, the partial
Lagrangian of problem (18) with relaxed constraints is given by

$\displaystyle\mathcal{L}(\eta,\bm{A},\bm{\lambda},\bm{\beta},\bm{\nu})$
$\displaystyle=\Big{(}1-\sum_{l=1}^{L}\sum_{n=1}^{N}\lambda_{l,n}\Big{)}\Gamma+\sum_{l=1}^{L}\sum_{n=1}^{N}\sum_{m=1}^{M}a_{l,m}[n]\Big{(}\lambda_{l,n}\frac{{\sf{mse}}_{l,m,n}}{\varepsilon_{l}}+\beta_{l,n}+\nu_{m,n}\Big{)}$
(19)
$\displaystyle-\sum_{l=1}^{L}\sum_{n=1}^{N}\beta_{l,n}-\sum_{m=1}^{M}\sum_{n=1}^{N}\nu_{m,n},$

where $\bm{\lambda}=\\{\lambda_{l,n},\forall l,\forall n\\}$,
$\bm{\beta}=\\{\beta_{l,n},\forall l,\forall n\\}$, and
$\bm{\nu}=\\{\nu_{m,n},\forall m,\forall n\\}$ are the Lagrange multipliers
associated with constraints (18), (3), and (4). Accordingly, the Lagrange dual
function of problem (18) with relaxed constraints is expressed as

$\displaystyle g(\bm{\lambda},\bm{\beta},\bm{\nu})=$
$\displaystyle\min\limits_{\Gamma,\bm{A}}\mathcal{L}(\eta,\bm{A},\bm{\lambda},\bm{\beta},\bm{\nu})$
(20b) $\displaystyle\sum_{n=1}^{N}\sum_{m=1}^{M}a_{l,m}[n]\geq D,\forall l,$
$\displaystyle 0\leq a_{l,m}[n]\leq 1,\forall l,\forall m,\forall n.$

which leads to the following lemma.

######  Lemma 1.

To make $g(\bm{\lambda},\bm{\beta},\bm{\nu})$ in problem (20) bounded, it must
follows that $(1-\sum_{l=1}^{L}\sum_{n=1}^{N}\lambda_{l,n})=0$.

###### Proof.

This can be proved by contradiction. If
$(1-\sum_{l=1}^{L}\sum_{n=1}^{N}\lambda_{l,n})>0$ or
$(1-\sum_{l=1}^{L}\sum_{n=1}^{N}\lambda_{l,n})<0$, then
$g(\bm{\lambda},\bm{\beta},\bm{\nu})\rightarrow-\infty$ by setting
$\Gamma\rightarrow-\infty$ or $\Gamma\rightarrow\infty$. Therefore, both of
the two above inequalities cannot be true. This completes the proof. ∎

From Lemma 1, dual variables $\bm{\lambda}$ are subject to an additional
constraint. Accordingly, the dual problem of (18) with relaxed constraints is
expressed as

$\displaystyle\mathop{\text{maximize}}\limits_{\bm{\lambda},\bm{\beta},\bm{\nu}}$
$\displaystyle g(\bm{\lambda},\bm{\beta},\bm{\nu})$ subject to
$\displaystyle\sum_{l=1}^{L}\sum_{n=1}^{N}\lambda_{l,n}=1,\bm{\beta}\succeq
0,\bm{\nu}\succeq 0.$ (21a)

The next step is to apply the Lagrange duality to find the primal optimal
solution.

####  III-A1 Obtaining $g(\bm{\lambda},\bm{\beta},\bm{\nu})$

First, we obtain the dual function $g(\bm{\lambda},\bm{\beta},\bm{\nu})$ under
given $\bm{\lambda},\bm{\beta},\bm{\nu}$ by solving problem (20). It is
observed that we can decompose Problem (20) into $L+1$ subproblems, each of
which can be solved in parallel. Particularly, the $L$ subproblems are for
optimizing $\bm{A}$, each of which is given by

$\displaystyle\mathop{\text{minimize}}\limits_{a_{l,m}[n],\forall m,\forall
n}$
$\displaystyle\sum_{n=1}^{N}\sum_{m=1}^{M}a_{l,m}[n]\Big{(}\lambda_{l,n}\frac{{\sf{mse}}_{l,m,n}}{\varepsilon_{l}}+\beta_{l,n}+\nu_{m,n}\Big{)}$
subject to $\displaystyle\sum_{n=1}^{N}\sum_{m=1}^{M}a_{l,m}[n]\geq D,$ (22b)
$\displaystyle 0\leq a_{l,m}[n]\leq 1,\forall m,\forall n,$

To minimize the objective function in (22), which is a linear combination of
$a_{l,m}[n]$, we should let the association coefficient corresponding to the
UAV with the top $D$ smallest
$\big{(}\lambda_{l,n}\frac{{\sf{mse}}_{l,m,n}}{\varepsilon_{l}}+\beta_{l,n}+\nu_{m,n}\big{)}$
be $1$ for any $l$. The optimal solution of problem (22) is given by the
following theorem.

######  Theorem 1.

For problem (22), by denoting
$x_{l,m}[n]=\lambda_{l,n}\frac{{\sf{mse}}_{l,m,n}}{\varepsilon_{l}}+\beta_{l,n}+\nu_{m,n}$,
the optimal cluster scheduling and association $\bm{A}$ can be expressed as

$\displaystyle\bm{a}_{l}^{\star}(\bm{\pi}_{l}(i))=\left\\{\begin{aligned}
&1,i=1,\ldots D,\\\ &0,\ \text{otherwise},\end{aligned}\right.$ (23)

where vector
$\bm{a}_{l}^{\star}\triangleq\big{[}a_{l,1}[1],\ldots,a_{l,1}[n],\ldots
a_{l,1}[N],a_{l,2}[1],\ldots a_{l,M}[N]\big{]}^{\sf T}\in\mathbb{R}^{MN\times
1}$, vector $\bm{x}_{l}\triangleq\big{[}x_{l,1}[1],\ldots,x_{l,1}[n],\ldots
x_{l,1}[N],\ldots,x_{l,M}[N]\big{]}^{\sf T}\in\mathbb{R}^{MN\times 1}$, and
permutation $\bm{\pi}_{l}=[\pi_{l}(1),\ldots,\pi_{l}(MN)]$ corresponds to the
ascend order such that
$\bm{x}_{l}(\pi_{l}(1))\leq\cdots\leq\bm{x}_{l}(\pi_{l}(MN))$.

Theorem 1 states that cluster $l$ should be associated with UAV $m$ at time
slot $n$ when coefficient $x_{l,m}[n]$ is smaller. According to the expression
of coefficient $x_{l,m}[n]$, the term
$\lambda_{l,n}\frac{{\sf{mse}}_{l,m,n}}{\varepsilon_{l}}$ is the effect of
introducing the ratio between the instantaneous MSE ${\sf{mse}}_{l,m,n}$ that
can be achieved for type-$l$ AirComp task and the target MSE threshold
$\varepsilon_{l}$ that should be achieved if cluster $l$ is associated with
UAV $m$ at time slot $n$. Furthermore, the remaining term represents other
effects due to the problem constraints. When $x_{l,m}[n]$ is smaller, it means
that it brings a smaller achievable ratio if cluster $l$ is associated with
UAV $m$ at time slot $n$.

And one subproblem is for optimizing $\Gamma$, i.e.,

$\displaystyle\mathop{\text{minimize}}\limits_{\Gamma}$
$\displaystyle(1-\sum_{l=1}^{L}\sum_{n=1}^{N}\lambda_{l,n})\Gamma$ (24)

Since $\sum_{l=1}^{L}\sum_{n=1}^{N}\lambda_{l,n}=1$, the optimal solution of
$\Gamma^{\star}$ of problem (24) can be any arbitrary real number. Without
loss of generality, we simply set

$\displaystyle\Gamma^{\star}=\max\\{a_{l,m}[n]\frac{{\sf{mse}}_{l,m,n}}{\varepsilon_{l}},\forall
l,\forall n,\forall m\\}.$ (25)

####  III-A2 Updating $\bm{\lambda},\bm{\beta},\bm{\nu}$

After obtaining $\bm{A}^{\star}$ and $\Gamma^{\star}$ for given
$\bm{\lambda},\bm{\beta}$, and $\bm{\nu}$, the optimal dual variables are
obtained by solvingthe dual problem (21). Since the dual function
$g(\bm{\lambda},\bm{\beta},\bm{\nu})$ is concave but non-differentiable,
subgradient-based methods such as the subgradient projection method are
adopted to update the dual variables $(\bm{\lambda},\bm{\beta},\bm{\nu})$.
Specifically, in the $t+1$-th iteration, the update of
$(\bm{\lambda},\bm{\beta},\bm{\nu})$ is given by

$\displaystyle\beta_{l,n}^{t+1}=\Big{[}\beta_{l,n}^{t}+\gamma(\sum_{m=1}^{M}a_{l,m}[n]-1)\Big{]}^{+},\nu_{m,n}^{t+1}=\Big{[}\nu_{m,n}^{t}+\gamma(\sum_{L=1}^{L}a_{l,m}[n]-1)\Big{]}^{+},$
(26)
$\displaystyle\lambda_{l,n}^{t+\frac{1}{2}}=\Big{[}\lambda_{l,n}^{t}+\gamma(\sum_{m=1}^{M}a_{l,m}[n](\frac{{\sf{mse}}_{l,m,n}}{\varepsilon_{l}}-\Gamma)\Big{]}^{+},\lambda_{l,n}^{t+1}=\lambda_{l,n}^{t+\frac{1}{2}}/(\sum_{l=1}^{L}\sum_{n=1}^{N}\lambda_{l,n}^{t+\frac{1}{2}}),\
$ (27)

where $\gamma$ is a dynamically chosen step-size sequence and
$[a]^{+}=\max(a,0)$.

####  III-A3 Constructing the optimal solution $\Gamma$

With the updated dual variables $\bm{\lambda}^{\star}$, $\bm{\beta}^{\star}$,
and $\bm{\nu}^{\star}$, we need to construct the primal solutions
$\Gamma^{\star}$ and $\bm{A}^{\star}$ to problem (18) with relaxed
constraints. The key observation is that with given $\bm{\lambda}^{\star}$,
$\bm{\beta}^{\star}$, and $\bm{\nu}^{\star}$, the optimal solution
$\bm{A}^{\star}$ can be uniquely obtained from (23). By substituting
$\bm{A}^{\star}$ into the primal problem (18), it is evident that
$\Gamma^{\star}=\max\\{a_{l,m}[n]\frac{{\sf{mse}}_{l,m,n}}{\varepsilon_{l}},\forall
l,\forall n,\forall m\\}$.

By iteratively optimizing primal variables $(\bm{A},\Gamma)$ and dual
variables $(\bm{\lambda},\bm{\beta},\bm{\nu})$, the optimal cluster scheduling
and association variables are obtained. Algorithm 1 summarizes how to solve
problem (18) via the dual method. It is important to note that the obtained
solution $a^{\star}_{l,m}[n]$ is either 1 or 0 according to (23), even though
we relax $a_{l,m}[n]$ as (20b). Consequently, the optimal solution to problem
(18) is obtained by using the dual method.

###  III-B Transmit Power Optimization

By substituting $\bm{A}$ into constraints (16), ${\sf MSE}_{l,n}$ is rewritten
as

$\displaystyle{\sf{MSE}}_{l,n}=\left\\{\begin{aligned} &{\sf
mse}_{l,m,n},(l,m,n)\in\mathcal{A},\\\ &0,\
\text{otherwise},\end{aligned}\right.$ (28)

where $\mathcal{A}\triangleq\\{(l,m,n)|a_{l,m}[n]=1,\forall l,\forall
m,\forall n\\}$. Thus, for any given $\\{\bm{\eta},\bm{Q},\bm{A}\\}$, the
transmit power optimization problem reduces to the following subproblem:

$\displaystyle\mathop{\text{minimize}}_{\Gamma,\bm{P}}$ $\displaystyle\Gamma$
subject to
$\displaystyle\sum_{k\in\mathcal{K}_{l}}\\!\\!\Big{(}{\theta_{l,k,m}[n]\sqrt{p_{k}[n]}}\\!-\\!1\Big{)}^{2}\\!\\!+\sum_{j\neq
l}\sum_{i\in\mathcal{K}_{j}}\theta_{l,i,m}^{2}[n]p_{i}[n]+\phi_{l,m}[n]\leq\Gamma\varepsilon_{l}K_{l}^{2},\forall(l,m,n)\\!\in\\!\mathcal{A},$
$\displaystyle\text{Constraints\ }\eqref{Cons:power},$

where $\theta_{l,k,m}[n]\triangleq\eta_{l,m}[n]|h_{k,m}[n]|,\forall l,\forall
m,\forall k$ and $\phi_{l,m}[n]\triangleq\eta_{l,m}^{2}[n]\sigma^{2},\forall
l,\forall m,\forall k$. Note that the number of constraints in (16) is $LN$,
whereas the number of constraints in (29) is $LD$. Problem (29) is a
quadratically constrained quadratic programming (QCQP), which can be solved
via CVX with interior-point solvers (e.g., Mosek).

######  Remark 3.

From constraints (29), to minimize $\Gamma$, if there is no UAV associated
with cluster $j$, i.e, $\sum_{m\in\mathcal{M}}a_{j,m}[n]=0$, then the
transmission power of device $i\in\mathcal{K}_{j}$ must be zero (i.e.,
$p_{i}[n]=0$). Otherwise, it will introduce interference for UAVs and thus
lead to a larger $\Gamma$.

###  III-C Normalizing Factors Optimization

For given $\\{\bm{A},\bm{Q},\bm{P}\\}$, the normalizing factor optimization
subproblem is reduce to the following formulation:

$\displaystyle\mathop{\text{minimize}}_{\Gamma,\bm{\eta}}$
$\displaystyle\Gamma$ subject to
$\displaystyle\sum_{m=1}^{M}\frac{a_{l,m}[n]}{K_{l}^{2}}\big{(}{\eta_{l,m}[n]\varphi_{k,m}[n]}-1\big{)}^{2}+\eta_{l,m}^{2}[n]\big{(}I_{l,m}[n]+\sigma^{2}\big{)}\leq\Gamma\varepsilon_{l}K_{l}^{2},\forall
l,\forall n,$ $\displaystyle\text{Constraints\ }\eqref{Cons:normalizing
factor},$

where $\varphi_{k,m}[n]=|h_{k,m}[n]|\sqrt{p_{k}[n]}$. Problem (30) is also a
convex QCQP. Furthermore, we can observe that problem (30) can be decoupled
into $LD$ subproblems, each of which optimizes $\eta_{l,m}[n]$ with
$a_{l,m}[n]=1$ to minimize ${\sf{mse}_{l,m,n}}$. For any $(l,m,n)$ with
$a_{l,m}[n]=1$, the subproblem is given by

$\displaystyle\mathop{\text{minimize}}\limits_{\eta_{l,m}[n]\geq
0}\sum_{k\in\mathcal{K}_{l}}\big{(}{\eta_{l,m}[n]\varphi_{k,m}[n]}-1\big{)}^{2}+\eta_{l,m}^{2}[n]\big{(}I_{l,m}[n]+\sigma^{2}\big{)}.$
(31)

And the optimal solution of problem (30) is given by the following
proposition.

######  Proposition 2.

Given by $\\{\bm{A},\bm{Q},\bm{P}\\}$, by setting the first derivative of the
objective function to zero, the optimal solution $\bm{\eta}$ to problem (31)
is expressed as

$\displaystyle\eta_{l,m}^{\star}[n]=\left\\{\begin{aligned}
&\frac{\sum\limits_{k\in\mathcal{K}_{l}}\sqrt{p_{k}[n]}|h_{k,m}[n]|}{\sum\limits_{k\in\mathcal{K}_{l}}p_{k}[n]|h_{k,m}[n]|^{2}+\sum_{j\in\mathcal{L}\backslash\\{l\\}}\sum_{i\in\mathcal{K}_{j}}p_{i}[n]|h_{i,m}[n]|^{2}+\sigma^{2}},\text{if}\
a_{l,m}[n]=1,\\\ &0,\ \text{otherwise}.\end{aligned}\right.$ (32)

By substituting the solution $\eta_{l,m}^{\star}[n]$ into problem (30), the
optimal solution $\Gamma$ is given by

$\displaystyle\Gamma^{\star}=\max\limits_{\forall l,\forall
n}\sum_{m=1}^{M}\frac{a_{l,m}[n]}{K_{l}^{2}}\Big{[}\sum_{k\in\mathcal{K}_{l}}\Big{(}\eta_{l,m}[n]\sqrt{p_{k}[n]}|h_{k,m}[n]|-1\Big{)}^{2}+\eta_{l,m}^{2}[n]\Big{(}I_{l,m}[n]+\sigma^{2}\Big{)}\Big{]}.$
(33)

######  Remark 4.

Note that with $a_{l,m}[n]=1$, the normalizing factor $\eta_{l,m}^{\star}[n]$
monotonically decreases with respect to the noise power $\sigma^{2}$ and the
inter-cluster interference power from other devices associated with other
UAVs, i.e.,
$\sum_{j\in\mathcal{L}\backslash\\{l\\}}\sum_{i\in\mathcal{K}_{j}}p_{i}[n]|h_{i,m}[n]|^{2}$.
From (33), we observe that reducing $\eta_{l,m}^{\star}[n]$ can suppress the
noise-induced error and inter-cluster interference induced error components
but increase the signal misalignment error.

###  III-D UAV Trajectory Optimization

For given $\\{\bm{A},\bm{P},\bm{\eta}\\}$, the UAVs’ trajectories subproblem
is reduced to the following subproblem:

$\displaystyle\mathop{\text{minimize}}\limits_{\Gamma,\bm{Q}}$
$\displaystyle\Gamma$ subject to $\displaystyle\text{Constraints\
}\eqref{Cons:Q position},\eqref{Cons:Q speed},\eqref{Cons:Q
collision},\eqref{Cons:target MSE1}.$ (34)

Due to the nonconvexity of constraints (15) and (16), problem (III-D) is
nonconvex. Generally speaking, there is no efficient method that can be used
to attain the optimal solution. To address their non-convexity, we exploit the
SCA technique in the sequel. We first transform constraints (16) into a
tractable form, which facilitates the development of the SCA technique. We
define

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!G_{l,m}[n]$ $\displaystyle=$
$\displaystyle\eta^{2}_{l,m}[n]\sum_{l^{\prime}\in\mathcal{L}}\sum\limits_{k\in\mathcal{K}_{l^{\prime}}}\frac{p_{k}[n]\beta_{0}}{(H^{2}+\left\|\bm{q}_{m}[n]-\mathbf{w}_{k}\right\|_{2}^{2})^{\frac{\gamma}{2}}},$
(35) $\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!F_{l,m}[n]$ $\displaystyle=$
$\displaystyle
2\eta_{l,m}[n]\sum\limits_{k\in\mathcal{K}_{l}}\dfrac{\sqrt{\beta_{0}p_{k}[n]}}{(H^{2}+\left\|\bm{q}_{m}[n]-\mathbf{w}_{k}\right\|_{2}^{2})^{\frac{\gamma}{4}}},$
(36) $\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!C_{l,m}[n]$ $\displaystyle=$
$\displaystyle\left\\{\begin{aligned}
&K+\eta^{2}_{l,m}[n]\sigma^{2},~{}\text{if~{}}a_{l,m}[n]=1,\\\ &0,\
\text{otherwise}.\end{aligned}\right.$ (37)

Based on the above introduced functions, we substitute $\bm{A}$ into ${\sf
MSE}_{l,n}$ in constraints (16), which is given by

$\displaystyle{\sf{MSE}}_{l,n}=\left\\{\begin{aligned}
&G_{l,m}[n]+C_{l,m}[n]-F_{l,m}[n],(l,m,n)\in\mathcal{A},\\\ &0,\
\text{otherwise}.\end{aligned}\right.$ (38)

With (38), constraints (16) are transformed into

$\displaystyle
G_{l,m}[n]+C_{l,m}[n]-F_{l,m}[n]\leq\Gamma\varepsilon_{l}K_{l}^{2},\forall(l,m,n)\in\mathcal{A}.$
(39)

Although $G_{l,m}[n]$ and $F_{l,m}[n]$ in constraints (39) are non-convex
w.r.t $\bm{q}_{m}[n]$, they are convex w.r.t
$\|\bm{q}_{m}[n]-\mathbf{w}_{k}\|_{2}^{2}$. Base on the key observation, we
first introduce slack variables
$\mathbf{S}=\\{s_{k,m}[n]\triangleq\left\|\bm{q}_{m}[n]-\mathbf{w}_{k}\right\|_{2}^{2}\
|\ \forall k,\forall m,\forall n\\}$ to tackle the non-convexity of
$G_{l,m}[n]$ in constraints (39). The reformulated problem (III-D) is given by

$\displaystyle\mathop{\text{minimize}}\limits_{\Gamma,\textbf{Q}}$
$\displaystyle\Gamma$ subject to
$\displaystyle\\!\\!\\!\hat{G}_{l,m}[n]+C_{l,m}[n]-F_{l,m}[n]\leq\varepsilon_{l}\Gamma,\forall(l,m,n)\in\mathcal{A},$
$\displaystyle
s_{k,m}[n]\leq\left\|\bm{q}_{m}[n]-\mathbf{w}_{k}\right\|_{2}^{2},\forall
k,\forall m,\forall n,$ $\displaystyle\text{Constraints\ }\eqref{Cons:Q
position},\eqref{Cons:Q speed},\eqref{Cons:Q collision},$

where

$\displaystyle\hat{G}_{l,m}[n]$ $\displaystyle=$
$\displaystyle\eta^{2}_{l,m}[n]\sum_{l^{\prime}\in\mathcal{L}}\sum\limits_{k\in\mathcal{K}_{l^{\prime}}}\frac{p_{k}[n]\beta_{0}}{(H^{2}+s_{k,m}[n])^{\frac{\gamma}{2}}}.$
(41)

Note that for problem (40), it can be easily verified that all constraints in
(40) can be met with equality. Otherwise, when other variables are fixed, we
can increase the value of $\bm{S}$ to further decrease the value of
$\\{\hat{G}_{l,m}[n]\\}$ without violating all constraints in problem (40),
thereby reducing the objective function. However, problem (40) is still a
nonconvex optimization problem because of the nonconvexity of constraints
(40), (40), and (15).

1: Input: $D$, tolerance $\epsilon>0$.

2: Initialize: $\bm{P}^{0}$, $\bm{\eta}^{0}$, and $\bm{Q}^{0}$, let $r=0$.

3: repeat

4: Given $\\!\\{\\!\bm{P}^{r},\bm{\eta}^{r},\bm{Q}^{r}\\!\\}$, solve problem
(18) with Algo. 1, and denote the solution as $\bm{A}^{r+1}$.

5: Given $\\!\\{\\!\bm{A}^{r+1},\bm{\eta}^{r},\bm{Q}^{r}\\!\\}$, solve convex
problem (29), and the solution is denoted as $\bm{P}^{r+1}$.

6: Given $\\!\\{\\!\bm{A}^{r+1},\bm{P}^{r+1},\bm{Q}^{r}\\!\\}$, solve problem
(30) based on (32), denote the solution as $\bm{\eta}^{r+1}$.

7: Given $\\!\\{\\!\bm{A}^{r+1},\bm{P}^{r+1},\bm{\eta}^{r+1}\\}$, solve
problem (46), and denote the solution as $\bm{Q}^{r+1}$.

8: Set $r=r+1$.

9: until The fractional decrease of the objective value is below a tolerance
$\epsilon$.

10: Output: $\Gamma$, $\bm{A}^{r}$, $\bm{P}^{r}$, $\bm{\eta}^{r}$, and
$\bm{Q}^{r}$.

Algorithm 2 BCD-SCA Algorithm for $\mathscr{P}_{1}$

Fortunately, $F_{l,m}[n]$ in constraints (40) is convex w.r.t
$\|\bm{q}_{m}[n]-\mathbf{w}_{k}\|_{2}^{2}$. And
$\|\bm{q}_{m}[n]-\mathbf{w}_{k}\|_{2}^{2}$ in constraints (40) and
$\left\|\bm{q}_{m}[n]-\bm{q}_{i}[n]\right\|_{2}$ in constraints (15) are
convex w.r.t $\bm{Q}$. It is important to recall that the first-order Taylor
expansion of any function at any point severs as a lower bound [32]. Thus,
similarly to [22], the SCA technique is adopted to tackle their non-convexity,
thus obtaining a suboptimal solution. Specifically, defining $\bm{Q}^{r}$ as
the output of the $r$-th iteration, we obtain that $F_{l,m}[n]$ is lower
bounded by the following expression, i.e.,

$\displaystyle F_{l,m}[n]$ $\displaystyle\geq$ $\displaystyle
F^{r}_{l,m}[n]\\!+\\!\nabla_{\|\bm{q}_{m}[n]-\mathbf{w}_{k}\|^{2}}F_{l,m}[n]\big{|}{\bm{q}_{m}[n]=\bm{q}_{m}^{r}[n]}\big{(}\left\|\bm{q}_{m}[n]-\mathbf{w}_{k}\right\|_{2}^{2}-\left\|\bm{q}_{m}^{r}[n]\\!-\\!\mathbf{w}_{k}\right\|_{2}^{2}\big{)}$
(42) $\displaystyle\triangleq$ $\displaystyle F^{\text{lb}}_{l,m}[n],$

where

$\displaystyle\nabla_{\|\bm{q}_{m}[n]-\mathbf{w}_{k}\|^{2}}F_{l,m}[n]\big{|}_{\bm{q}_{m}[n]=\bm{q}_{m}^{r}[n]}=\sum\limits_{k\in\mathcal{K}_{l}}-\dfrac{{\eta_{l,m}[n]}\gamma\sqrt{p_{k}[n]\beta_{0}}}{2(H^{2}+\left\|\bm{q}_{m}^{r}[n]-\mathbf{w}_{k}\right\|_{2}^{2})^{\frac{\gamma+4}{4}}}.$
(43)

It is easily verified that $F^{\text{lb}}_{l,m}[n]$ is concave with regard to
$\bm{q}_{m}[n]$. Hence, $\hat{G}_{l,m}[n]+C_{l,m}[n]-F^{\text{lb}}_{l,m}[n]$
is a convex function. Likewise, with the first-order Taylor expansion at point
$\bm{q}_{m}^{r}[n]$, $\left\|\bm{q}_{m}[n]-\mathbf{w}_{k}\right\|_{2}^{2}$ in
constraints (40) are subject to

$\displaystyle\left\|\bm{q}_{m}[n]-\mathbf{w}_{k}\right\|_{2}^{2}$
$\displaystyle\geq$
$\displaystyle\left\|\bm{q}^{r}_{m}[n]-\mathbf{w}_{k}\right\|_{2}^{2}+2({\bm{q}}^{r}_{m}[n]-\mathbf{w}_{k})^{T}(\bm{q}_{m}[n]-\bm{q}^{r}_{m}[n])\triangleq
D^{\text{lb}}_{k,m}[n].$ (44)

Note that the lower bound function $D^{\text{lb}}_{k,m}[n]$ is a linear
function w.r.t $\bm{q}_{m}[n]$. Likewise,
$\left\|\bm{q}_{m}[n]-\bm{q}_{i}[n]\right\|_{2}^{2}$ in constraint (15) is
lower-bounded by

$\displaystyle\\!\\!\left\|\bm{q}_{m}[n]-\bm{q}_{i}[n]\right\|^{2}$
$\displaystyle\geq$
$\displaystyle-\left\|\bm{q}^{r}_{m}[n]-\bm{q}_{i}^{r}[n]\right\|^{2}+2(\bm{q}^{r}_{m}[n]-\bm{q}_{i}^{r}[n])^{T}(\bm{q}_{m}[n]-\bm{q}_{i}[n])\triangleq
d^{\text{lb}}_{m,i}[n].$ (45)

Given point $\bm{Q}^{r}$, based on expressions (42), (44), and (45), problem
(40) approximates as follows:

$\displaystyle\mathop{\text{minimize}}\limits_{\Gamma,\textbf{Q}}$
$\displaystyle\Gamma$ subject to
$\displaystyle\\!\\!\hat{G}_{l,m}[n]\\!+\\!C_{l,m}[n]\\!-\\!F^{\text{lb}}_{l,m}[n]\\!\leq\\!\varepsilon_{l}\Gamma,\\!\forall(l,m,n)\in\mathcal{A},$
$\displaystyle 0\leq s_{k,m}[n]\leq D^{\text{lb}}_{k,m}[n],\forall k,\forall
m,\forall n,$ $\displaystyle d^{\text{lb}}_{m,i}[n]\geq d_{\text{min}},\forall
n,\forall m,m\neq i,$ $\displaystyle\text{Constraints~{}}\eqref{Cons:Q
position},\eqref{Cons:Q speed}.$

Since constraints in (46), (46), and (46) are jointly convex w.r.t $\bm{Q}$
and $\bm{S}$ now, problem (46) is convex, which can be solved via modeling
framework CVX and interior-point solvers (e.g., Mosek). Based on the lower
bounds in (46), we can conclude that any feasible solution to problem (46) is
also feasible to problem (III-D). Additionally, $-F^{\rm lb}_{l,m}[n]$ is an
upper bound of $-F_{l,m}[n]$. Therefore, the optimal objective value of the
approximate problem (46) is generally an upper bound to problem (III-D).

1: Input: $D_{\min}=0$, $D_{\max}=\lfloor\frac{TM}{\delta L}\rfloor$.

2: Initialize: $\bm{P}^{0}$ , $\bm{\eta}^{0}$, and $\bm{Q}^{0}$, let $t=1$.

3: repeat

4: Update $D=\lceil\frac{D_{\min}+D_{\max}}{2}\rceil$. Initialize $\bm{P}^{0}$
and $\bm{\eta}^{0}$.

5: Obtain $\Gamma^{t}$, $\bm{A}^{t}$, $\bm{P}^{t}$, $\bm{\eta}^{t}$, and
$\bm{Q}^{t}$ using Algorithm 2.

6: if $\Gamma^{t}\leq 1$ then

7: Set $D_{\min}=D$. And let $D^{\star}=D$, $\bm{A}^{\star}=\bm{A}^{t}$,
$\bm{P}^{\star}=\bm{P}^{t}$, $\bm{\eta}^{\star}=\bm{\eta}^{t}$, and
$\bm{Q}^{\star}=\bm{Q}^{t}$.

8: else

9: Set $D_{\max}=D$.

10: end if

11: Update $t=t+1$.

12: until $D_{\max}-D_{\min}\leq 1$.

13: Output: $D^{\star}$, $\bm{A}^{\star}$, $\bm{P}^{\star}$ ,
$\bm{\eta}^{\star}$, and $\bm{Q}^{\star}$.

Algorithm 3 Overall Bisection Algorithm for $\mathscr{P}_{2}$

###  III-E Proposed Overall Algorithm

The BCD-SCA method for solving $\mathscr{P}_{1}$ is summarized in Algorithm 2,
where the cluster association, transmit power, normalizing factors, and UAV
trajectories are successively optimized while keeping the other variables
fixed until convergence. In addition, the obtained solution at current
iteration will be applied to be the input of the next iteration. And, the
computational complexity of solving problem (18) is $\mathcal{O}(LMN)$. The
complexity of solving problem (29) is $\mathcal{O}(K^{3.5}L^{1.5}D^{1.5})$.
The complexity of solving problem (30) is $\mathcal{O}(KMN)$. The complexity
of solving problem (III-D) is $\mathcal{O}(K^{1.5}M^{3.5}N^{3.5})$. According
to the aforementioned complexity analysis of each subproblem, the overall
complexity of Algorithm 2 is
$\mathcal{O}((K^{3.5}L^{1.5}D^{1.5}+K^{1.5}M^{3.5}N^{3.5})log(1/\epsilon))$.

Recall that problem $\mathscr{P}_{2}$ can be efficiently solved by applying
the bisection search over the minimum amount of AirComp tasks $D$, then the
overall bisection algorithm for $\mathscr{P}_{2}$ is summarized in Algorithm
3. Although the obtained solution is generally suboptimal, we show the
effectiveness of our proposed algorithm in increasing the number of AirComp
tasks among all clusters via numerical simulations, and compare it to other
benchmarks in Section IV.

Since our proposed algorithm is an iterative algorithm, in the following, we
present the initial procedures for transmit power, normalizing factors, and
trajectories. 1) Transmit power initialization: The transmit powers of the
ground devices are initialized by the equal transmit power in all time slots,
i.e., $p_{k}[n]=P_{k}/D,\forall k,\forall n$. 2) Trajectory initialization:
The trajectory is initialized by a simple circular trajectory scheme, which is
detailed in [22]. 3) Normalizing factor initialization: Given the initial
transmit power and trajectories, the initial normalizing factors can be
obtained by computing (32).

##  IV Numerical Results

In this section, numerical results are provided to verify performance gain of
the proposed design in terms of the max-min AirComp task amount and the
effectiveness of Algorithm 3. We consider $L=6$ clusters. Therein, each
cluster has $|K_{\ell}|=20,\forall\ell$ ground devices, each of which randomly
distributes in a circular area with a radius of $r=50$ meters, where the
centers of circles are set to $(x_{C_{1}},y_{C_{1}},z_{C_{1}})=(100,50,0)$
meters, $(x_{C_{2}},y_{C_{2}},z_{C_{2}})=(200,200,0)$ meters,
$(x_{C_{3}},y_{C_{3}},z_{C_{3}})=(-100,100,0)$ meters,
$(x_{C_{4}},y_{C_{4}},z_{C_{4}})=(-400,150,0)$ meters,
$(x_{C_{5}},y_{C_{5}},z_{C_{5}})=(-200,-200,0)$ meters, and
$(x_{C_{6}},y_{C_{6}},z_{C_{6}})=(-250,-100,0)$ meters. The UAVs fly at a
fixed altitude $H=100$ meters to comply with the rule that all commercial UAVs
should not fly over $400$ feet ($122$ meters) [33]. In addition, the minimum
distance between any two UAVs is set as $d_{\text{min}}=100$ meters [22]. The
maximum speed of all UAVs is assumed to be the same and set as
$V_{\text{max}}=30$ m/s. Time step size $\delta$ is set as $0.5$ s, which is
small enough to satisfy $\delta V_{\text{max}}\ll H$. The noise power at the
receiver and the channel power gain at the reference distance of $d_{0}=1$ m
are set as $\sigma^{2}=-80$ dBm and $\beta_{0}=-50$ dB, respectively. The
target AirComp MSE threshold for all clusters is set as
$\varepsilon_{l}=2\times 10^{-3},\forall l$. Let the threshold $\epsilon$ in
Algorithm 2 be $10^{-3}$.

###  IV-A Multi-Cluster AirComp with a Single UAV

In this subsection, to show the superiority of our proposed design and the
effectiveness of Algorithm 3 for solving the minimum amount of AirComp tasks
maximization problem, we first consider a special case that there is only one
UAV serving ground devices, i.e., $M=1$. In this case, the system is free of
inter-cluster interference.

Figure 2: Optimized UAV trajectories in a single UAV case.

Figure 3: Cluster scheduling schemes over timeline.

Fig. 3 illustrates the obtained UAV trajectories optimized by using Algorithm
3 with varying mission duration $T$ when $E=0.8$ W. Each trajectory is sampled
every three seconds marked with “$\triangle$” by using the same colors. The
locations of devices are marked by dark blue “$\square$”. We observe that,
when $T$ is large, the UAV makes use of its maneuverability and adjusts its
trajectory to approach its associated cluster as much as possible. In
particular, the UAV hovers over its associated cluster for a certain amount of
time, except for the minimum time spent in flight between clusters. Under the
given mission duration $T=150$ s, the UAV can move sufficiently close to each
cluster in sequence. To this end, the optimized UAV trajectory contains line
segments, which connect the points on the top of all clusters.

Fig. 3 shows the corresponding cluster scheduling along the timeline. The
different colored rectangles represent different clusters, and their lengths
represent how long these clusters are scheduled. It is observed that for both
$T=50$ s and $T=150$ s, the UAV visits each cluster in turn and associates
with each cluster within an equal time period to complete the corresponding
AirComp task under the target MSE requirements, which is in accordance with
the observed results in Fig. 3.

To show the performance gain of our proposed joint design, we consider the
following benchmarks.

  * •

Static UAV: This scheme fixes the location of UAV on the geometric point of
all devices and only optimizes $\bm{A}$, $\bm{\eta}$, and $\bm{P}$ by using
Algorithm 3.

  * •

Equal power transmission: At each scheduled time slot, each device is
allocated equal power (i.e., $p_{k}[n]=P/D$). Only three variables (i.e.,
$\bm{A}$, $\bm{\eta}$, and $\bm{Q}$) are optimized by using Algorithm 3.

  * •

Upper bound: In this scheme, we assume that all devices have sufficient power
budgets and the system is interference-free. Therefore, the maximum task
amount for each cluster can be achieved, i.e.,

$\displaystyle D^{\text{ub}}=\lfloor\frac{MT}{\delta L}\rfloor.$ (47)

Figure 4: Max-min AirComp task amount versus power budget $P$ in a single UAV
case.

Figure 5: Max-min AirComp task amount versus mission duration $T$ in a single
UAV case.

Fig. 5 shows the average max-min task amount versus the total power budget
when $T=200$ s. First, with the fixed mission duration, observing from (47),
the upper bound of max-min AirComp task amount is fixed. With the increase of
power budget $P$, the max-min task amount achieved by all three schemes
increases. This is because, with a larger $P$, the ground devices can afford
more energy to compensate channel fading to reduce the achievable AirComp MSE.
However, as long-distance transmission may lead to deep fading, the
performance gap between the static UAV scheme and upper bound is large even in
the relatively large $P$ (i.e., $P=0.8$ W). Conversely, by utilizing UAV
mobility, the proposed joint design scheme and the equal power transmission
scheme significantly outperform the static UAV scheme. This is because the
mobile UAV dynamically constructs favorable communication channels to avoid
deep attenuation. In addition, the performance of the proposed joint design is
the best and reaches the upper bound when $P$ is relatively small (i.e.,
$P=0.3$ W). To achieve uniform powers of received signals, higher power is
required to compensate for worse channel conditions. Using identical power
transmission results in several misaligned signals since all devices’ channel
conditions change due to UAV mobility, especially when the power budget is
small but ineffectively allocated. Fortunately, the proposed joint design
scheme exploits the synergy of trajectory design and power control to improve
signal alignment. All the above results illustrate the importance and
necessity of the joint design in maximizing the max-min AirComp task amount.
This demonstrates the effectiveness of Algorithm 3 for solving the joint
design optimization problem $\mathscr{P}_{2}$. It shows that in the single-UAV
case without inter-cluster interference, the max-min task amount can achieve
the upper bound by increasing the transmit power.

Fig. 5 illustrates the average max-min task amount versus the mission duration
$T$ when $P=0.8$ W. In the case of the static UAV, performance is independent
of mission duration because of the time-invariant channel conditions between
the UAV and devices. Conversely, with optimized trajectory design, the max-min
AirComp task amount achieved by the proposed joint design scheme and the equal
power transmission scheme increase with $T$ increasing. There are two main
reasons. First, the UAV with the optimized flight trajectory can establish
good enough channel conditions with its associated cluster. Second, as $T$
increases, the UAV can spend more time performing tasks from its associated
cluster under the established favorable channel conditions. We can further
observe that the performance gap between the upper bound and the equal power
transmission scheme increases as $T$ increases. However, thanks to effectively
allocating power over all time slots, the proposed joint design scheme can
reach the upper bound from $T=160$ s to $T=360$ s. Nevertheless, we can
observe that there is a fundamental tradeoff between performance and delay. In
the case of a single UAV, an increased max-min AirComp task amount means an
increased access delay since the UAV sequentially visits each cluster and thus
the cluster needs to wait longer for communication until the previous
scheduled clusters to complete task amount.

Figure 6:  UAV trajectories optimized by the proposed joint design in two-UAVs
case.

Figure 7: The corresponding cluster-UAV scheduling and association schemes
over timeline in two-UAVs case.

###  IV-B Multi-Cluster AirComp with Multiple UAVs

Next, we study the multiple UAV case to illustrate its further performance
gain. To validate the effectiveness of our proposed cooperative interference
management, here we take double UAVs as a simple example. Additionally, our
proposed algorithm is also valid for more than two UAVs. We compare the
performance of the proposed design to three benchmarks. In the static UAV
scheme, the UAVs are located at the geometric points of devices in clusters 1,
2, 6 and 3, 4, 5, respectively. The following orthogonal UAV transmission
scheme is included for comparison.

  * •

Orthogonal transmission: The UAVs provide services to clusters via TDMA
without inter-cluster interference. In this scheme, the binary variables are
subject to

$\displaystyle\sum\limits_{l=1}^{L}\sum\limits_{m=1}^{M}a_{l,m}[n]\leq
1,n\in\mathcal{N};a_{l,m}[n]\in\\{0,1\\},\forall
l\in\mathcal{L},m\in\mathcal{M},n\in\mathcal{N}.$ (48)

This problem is a special instance of problem (15), which can be solved by
using a similar Algorithm 3.

Fig. 7 shows the trajectories of double UAVs optimized by the proposed design
when $T=80$ s and $P=0.8$ W. To better observe phenomena, we also plot the
corresponding cluster scheduling and association schemes along the timeline in
Fig. 7. From Fig. 7 and Fig. 7, it can be observed that each cluster is served
by the same UAV and each UAV only needs to be responsible for half of
clusters. From Fig. 7, we can observe that there is at least one UAV to
complete one AirComp task at each time slot. Additionally, at some time
intervals (e.g., from $t=0$ s to $t=3$ s and from $t=10$ s to $t=30$ s), there
are two scheduled clusters, which means two UAVs actually work simultaneously
to deal with different AirComp tasks. Furthermore, observed from Fig. 7 and
Fig. 7, with the cluster-UAV scheduling and association design, to avoid
strong co-channel interference, in most time slots, UAVs tend to choose to
serve two clusters that are geographically far away from each other. As a
result, the two UAVs can properly enlarge their trajectories to move closer to
their associated cluster for better channel conditions establishment while
achieving weak co-channel interference for the other cluster. For example, as
shown in Fig. 7, from $t=70$ s to $t=80$ s (i.e., $t=0$ s), UAV $2$ associates
with cluster $C_{4}$ while UAV $1$ associates with cluster $C_{1}$. During
this time interval, from Fig. 7, we can observe that UAV $2$ and UAV $1$ are
respectively closer to cluster $C_{4}$ and cluster $C_{1}$. In order to
complete as many tasks as possible for each cluster, two UAVs will at some
point have to serve two clusters simultaneously when flying. For example, as
shown in Fig. 7, from $t=35$ s to $t=40$ s, UAV $2$ associates with cluster
$C_{3}$ while UAV $1$ associates with cluster $C_{6}$. In these time instants,
the optimized trajectories show that UAVs tend to stay away from each other to
reduce co-channel inter-cluster interference. This is why the UAVs fly along
arc paths with opposite opening directions during that time interval.

Figure 8: Max-min task amount versus power budget $P$ in two-UAVs case.

Figure 9: Max-min task amount versus mission duration $T$ in two-UAVs case.

Fig. 9 shows the average max-min task amount versus device power budget when
$T=80$ s. First, due to effective interference management, the proposed joint
design scheme outperforms all benchmarks with all regimes considered. Second,
at the low power budget regime, the performance gaps between the proposed
joint design and other benchmarks are relatively large. This implies the
effectiveness of power control optimization in suppressing the noise-induced
error that is dominant for the achievable MSE in the low power budget regime.
Moreover, the proposed joint design and the equal-power transmission scheme
outperform the orthogonal transmission scheme, demonstrating the benefit of
cooperative transmission over clusters. And their performance gaps become
large as the power budget increases. In addition, it is observed that all
schemes suffering from inter-cluster interference become saturated in the high
power budget regime, meaning that the achievable MSE performance is limited by
the inter-cell interference and cannot be reduced further by simply increasing
the transmit power.

Fig. 9 shows the average max-min task amount versus the mission duration $T$
when $P=0.8$ W in the case of two UAVs. It can be observed that the
performance of static UAV is the worst due to long distance transmission.
Compared to the curves of the proposed joint design scheme in Fig. 5, it can
be observed that with the same given $T$, the deployment of more than one UAV
with cooperative interference management achieves a larger max-min task
amount, thereby improving the performance-delay tradeoff. For example, when
$T=160$ s, max-min AirComp task amount achieved by the proposed joint design
in the single-UAV system is 53, whereas that in the two-UAV system exceeds 92.

##  V Conclusion

We study the minimum task amount maximization problem in a UAV-aided AirComp
system with multiple clusters, taking into account cluster scheduling and
association, the UAV trajectory design, and AirComp transceiver design. The
joint design aims to align signals within each cluster while mitigating the
detrimental inter-cluster interference, thereby improving WDA performance. By
applying the Bisection technique, the formulated mixed-integer non-convex
optimization problem is reduced to a sequence of the maximum ratio
minimization problems. This resulting problem is further solved by developing
an efficient algorithm by applying BCD, Lagrange dual ascent, and SCA
techniques. Simulation results demonstrate significant performance gains of
our proposed joint design over the benchmark schemes.

## References

  * [1] K. B. Letaief, Y. Shi, J. Lu, and J. Lu, “Edge artificial intelligence for 6G: Vision, enabling technologies, and applications,” _IEEE J. Sel. Areas Commun._ , vol. 40, no. 1, pp. 5–36, Jan. 2022.
  * [2] M. Agiwal, A. Roy, and N. Saxena, “Next generation 5G wireless networks: A comprehensive survey,” _IEEE Commun. Surveys Tut._ , vol. 18, no. 3, pp. 1617–1655, 3rd Quart., 2016., 2016.
  * [3] Y. Shi, K. Yang, T. Jiang, J. Zhang, and K. B. Letaief, “Communication-efficient edge AI: Algorithms and systems,” _IEEE Commun. Surv. Tutorials_ , vol. 22, no. 4, pp. 2167–2191, 2020.
  * [4] G. Zhu, J. Xu, K. Huang, and S. Cui, “Over-the-air computing for wireless data aggregation in massive IoT,” _IEEE Wirel. Commun._ , vol. 28, no. 4, pp. 57–65, Aug. 2021.
  * [5] M. Goldenbaum, H. Boche, and S. Stańczak, “Harnessing interference for analog function computation in wireless sensor networks,” _IEEE Trans. Signal Process._ , vol. 61, no. 20, pp. 4893–4906, Oct. 2013.
  * [6] F. Molinari, N. Agrawal, S. Stanczak, and J. Raisch, “Max-consensus over fading wireless channels,” _IEEE Trans. Control Netw. Syst._ , vol. 8, no. 2, pp. 791–802, Jun. 2021.
  * [7] W. Liu, X. Zang, Y. Li, and B. Vucetic, “Over-the-air computation systems: Optimization, analysis and scaling laws,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 8, pp. 5488–5502, Dec. 2020.
  * [8] K. Yang, T. Jiang, Y. Shi, and Z. Ding, “Federated learning via over-the-air computation,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 3, pp. 2022–2035, Mar. 2020.
  * [9] Z. Wang, J. Qiu, Y. Zhou, Y. Shi, L. Fu, W. Chen, and K. B. Letaief, “Federated learning via intelligent reflecting surface,” _IEEE Trans. Wireless Commun._ , vol. 21, no. 2, pp. 808–822, Feb. 2022.
  * [10] X. Cao, G. Zhu, J. Xu, and K. Huang, “Optimized power control for over-the-air computation in fading channels,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 11, pp. 7498–7513, Nov. 2020.
  * [11] L. Chen, X. Qin, and G. Wei, “A uniform-forcing transceiver design for over-the-air function computation,” _IEEE Wireless Commun. Lett._ , vol. 7, no. 6, pp. 942–945, Dec. 2018.
  * [12] W. Fang, Y. Jiang, Y. Shi, Y. Zhou, W. Chen, and K. B. Letaief, “Over-the-air computation via reconfigurable intelligent surface,” _IEEE Trans. Commun._ , vol. 69, no. 12, pp. 8612–8626, Dec. 2021.
  * [13] G. Zhu and K. Huang, “MIMO over-the-air computation for high-mobility multimodal sensing,” _IEEE Internet Things J._ , vol. 6, no. 4, pp. 6089–6103, Aug. 2019.
  * [14] X. Yuan, Y. J. Angela Zhang, Y. Shi, W. Yan, and H. Liu, “Reconfigurable-intelligent-surface empowered wireless communications: Challenges and opportunities,” _IEEE Wireless Commun._ , vol. 28, no. 2, pp. 136–143, Apr. 2021.
  * [15] M. Fu, Y. Zhou, Y. Shi, and K. B. Letaief, “Reconfigurable intelligent surface empowered downlink non-orthogonal multiple access,” _IEEE Trans. Commun._ , vol. 69, no. 6, pp. 3802–3817, Jun. 2021.
  * [16] M. Fu, Y. Zhou, Y. Shi, T. Wang, and W. Chen, “UAV-assisted over-the-air computation,” in _Proc. IEEE Int. Conf. Commun. (ICC)_ , 2021.
  * [17] M. Fu, Y. Zhou, Y. Shi, W. Chen, and R. Zhang, “UAV aided over-the-air computation,” _IEEE Trans. Wireless Commun._ , 2021, doi: 10.1109/TWC.2021.3134327.
  * [18] W. Wen, Y. Jia, and W. Xia, “Joint scheduling and resource allocation for federated learning in SWIPT-enabled micro UAV swarm networks,” _China Communications_ , vol. 19, no. 1, pp. 119–135, Jan. 2022.
  * [19] Q. Lan, H. S. Kang, and K. Huang, “Simultaneous signal-and-interference alignment for two-cell over-the-air computation,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 9, pp. 1342–1345, Sept. 2020.
  * [20] Y. Shi, J. Zhang, and K. B. Letaief, “Group sparse beamforming for green cloud-RAN,” _IEEE Trans. Wireless Commun._ , vol. 13, no. 5, pp. 2809–2823, May 2014.
  * [21] X. Cao, G. Zhu, J. Xu, and K. Huang, “Cooperative interference management for over-the-air computation networks,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 4, pp. 2634–2651, Apr. 2021.
  * [22] Q. Wu, Y. Zeng, and R. Zhang, “Joint trajectory and communication design for multi-UAV enabled wireless networks,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 3, pp. 2109–2121, Mar. 2018.
  * [23] C. Zhan and Y. Zeng, “Completion time minimization for multi-UAV-enabled data collection,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 10, pp. 4859–4872, Oct. 2019.
  * [24] C. Shen, T.-H. H. Chang, J. Gong, Y. Zeng, and R. Zhang, “Multi-UAV interference coordination via joint trajectory and power control,” _IEEE Trans. Signal Process._ , vol. 68, pp. 843–858, Jan. 2020.
  * [25] Y. Xu and W. Yin, “A block coordinate descent method for regularized multiconvex optimization with applications to nonnegative tensor factorization and completion,” _SIAM J. Imag. Sci._ , vol. 6, no. 3, pp. 1758–1789, 2013.
  * [26] T. Ma, H. Zhou, B. Qian, N. Cheng, X. Shen, X. Chen, and B. Bai, “UAV-LEO integrated backbone: A ubiquitous data collection approach for B5G internet of remote things networks,” _IEEE J. Sel. Areas Commun._ , vol. 39, no. 11, pp. 3491–3505, Nov. 2021.
  * [27] D. W. Matolak and R. Sun, “Air-ground channel characterization for unmanned aircraft systems—part III: The suburban and near-urban environments,” _IEEE Trans. Veh. Techn._ , vol. 66, no. 8, pp. 6607–6618, Aug. 2017.
  * [28] _Study on Channel Model for Frequencies From 0.5 to 100 GHz_ , document 3GPP-TR-38.901, 3GPP, 2017.
  * [29] I. Qualcomm Technologies, _LTE unmanned aircraft systems_ , San Diego, CA, USA, Trial report v.1.0.1, 2017.
  * [30] U. Mengali and A. N. D’Andrea, _Synchronization Techniques for Digital Receivers_ , New York, NY, USA: Springer, 1997.
  * [31] O. Abari, H. Rahul, D. Katabi, and M. Pant, “Airshare: Distributed coherent transmission made seamless,” in _Proc. IEEE Conf. Comput. Commun. (INFOCOM)_ , 2015, pp. 1742–1750.
  * [32] S. Boyd and L. Vandenberghe, _Convex Optimization_. Cambridge university press, 2004.
  * [33] FAA, _Summary of small unmanned aircraft rule_ , Washington, DC: Federal Aviation Administration, 2016.