# Multi-Antenna UAV Data Harvesting: Joint Trajectory and Communication
Optimization

Jingwei Zhang, Yong Zeng, and Rui Zhang  J. Zhang and R. Zhang are with the
Department of Electrical and Computer Engineering, National University of
Singapore, Singapore 117583 (e-mail<EMAIL_ADDRESS>elezhang@nus.edu.sg).Y. Zeng is with the National Mobile Communications
Research Laboratory, Southeast University, Nanjing 210096, China. He is also
with the Purple Mountain Laboratories, Nanjing 211111, China (e-mail:
yong_zeng@seu.edu.cn).

###### Abstract

Unmanned aerial vehicle (UAV)-enabled communication is a promising technology
to extend coverage and enhance throughput for traditional terrestrial wireless
communication systems. In this paper, we consider a UAV-enabled wireless
sensor network (WSN), where a multi-antenna UAV is dispatched to collect data
from a group of sensor nodes (SNs). The objective is to maximize the minimum
data collection rate from all SNs via jointly optimizing their transmission
scheduling and power allocations as well as the trajectory of the UAV, subject
to the practical constraints on the maximum transmit power of the SNs and the
maximum speed of the UAV. The formulated optimization problem is challenging
to solve as it involves non-convex constraints and discrete-value variables.
To draw useful insight, we first consider the special case of the formulated
problem by ignoring the UAV speed constraint and optimally solve it based on
the Lagrange duality method. It is shown that for this relaxed problem, the
UAV should hover above a finite number of optimal locations with different
durations in general. Next, we address the general case of the formulated
problem where the UAV speed constraint is considered and propose a traveling
salesman problem (TSP)-based trajectory initialization, where the UAV
sequentially visits the locations obtained in the relaxed problem with minimum
flying time. Given this initial trajectory, we then find the corresponding
transmission scheduling and power allocations of the SNs and further optimize
the UAV trajectory by applying the block coordinate descent (BCD) and
successive convex approximation (SCA) techniques. Finally, numerical results
are provided to illustrate the spectrum and energy efficiency gains of the
proposed scheme for multi-antenna UAV data harvesting, as compared to
benchmark schemes.

###### Index Terms:

Wireless sensor network (WSN), multi-antenna communication, unmanned aerial
vehicle (UAV), rate maximization, trajectory optimization.

##  I Introduction

As an essential part of the Internet of Things (IoT), wireless sensor networks
(WSNs) have found a proliferation of applications in many fields, such as
surveillance and monitoring, automated and cyber-physical systems, and so on
[1]. An important function of WSNs is data harvesting from a set of
distributed sensor nodes (SNs). Various techniques such as clustering,
multihop data relaying, and in-network data aggregation have been proposed to
improve the data collection efficiency and reduce the energy consumption of
the SNs in WSNs [2]. However, due to the static network infrastructure and the
short transmission range of low-power SNs, it is generally difficult for WSNs
to sustain a high throughput in the long term [3].

On the other hand, thanks to various advantages such as flexible and on-demand
deployment and high probability of having line-of-sight (LoS) communication
links with the ground terminals, unmanned aerial vehicle (UAV) has been
envisioned as a promising technique for the future wireless communication
systems to support massive IoT devices [4]. There are various UAV applications
in wireless communication systems such as UAV-aided ubiquitous coverage, UAV-
aided relaying, and UAV-aided information dissemination and data collection.
In particular, by leveraging the high and controllable mobility of UAVs, UAV-
mounted access point (AP) is regarded as a promising new solution to collect
data from widely spread SNs in a geographically large area. By properly
designing the UAV trajectory to move closer to the SNs, not only the coverage
and throughput performance of the WSN can be significantly improved as
compared to the traditional fixed APs on the ground, but also the energy
consumption of SNs can be considerably reduced, thus prolonging the lifetime
of WSNs. However, UAV-enabled data harvesting in WSNs also faces new design
challenges. Firstly, the trajectory of the UAV AP needs to be jointly designed
with the conventional transmission scheduling and power control of the SNs to
maximize the spectrum efficiency for data collection [5]. Secondly, besides
the limited energy of SNs, UAVs are usually battery-powered with limited
endurance, which needs to be taken into account in the above joint design to
minimize the UAV’s propulsion energy consumption [6], [7].

In recent years, significant research efforts have been devoted to designing
and optimizing the performance of UAV-enabled data harvesting systems. To
reduce the energy consumption of SNs, their wake-up scheduling was jointly
designed with the UAV’s trajectory in [8] to minimize the maximum energy
consumption of all SNs, subject to the data rate requirement for each SN. In
[9], the performance of data harvesting from massive IoT devices to a group of
UAVs was analyzed and evaluated in terms of coverage and rate. In [10], the
UAV trajectory was optimized in three-dimensional (3D) space under the angle-
dependent Rician fading channel model between the UAV and its served SNs on
the ground. The projection-based compressive data gathering was investigated
in [11] for energy-efficient UAV data collection from distributed SNs. In
[12], a low-complexity trajectory design algorithm was proposed based on the
principle of receding horizon. However, the above prior works considered the
scenarios that the UAV is equipped with a single antenna and thus only one SN
can be scheduled for transmission at each time instant. As a result, the UAV
usually needs to move sufficiently close to each SN for data collection, which
leads to not only high energy consumption of the UAV, but also limited
throughput for data collection from SNs.

To tackle the above issue, in this paper we propose the use of multiple
antennas at the UAV for simultaneously harvesting data from multiple SNs by
exploiting the multi-antenna beamforming and spatial multiplexing gains. This
will help improve the spectrum efficiency for data collection as compared to
the conventional single-antenna UAVs, and also reduce the UAV flying distance
and hovering time and hence its propulsion energy consumption. In [13], the
authors reported measurement results that show significant performance gains
of an airborne multiple-input multiple-output (MIMO) system over the single-
input single-output (SISO) system in terms of coverage and throughput. In
[14], the heading direction of a multi-antenna UAV was optimized for
maximizing the sum-rate from a set of ground nodes in their uplink
communication. An antenna array composed of multiple single-antenna UAVs was
proposed in [15] to provide services to ground users in a collaborative
manner, where UAVs are only allowed to communicate with ground users at a set
of hovering locations, but no communication occurs when UAVs are moving. In
[16], a virtual MIMO link was formed where a multi-antenna UAV was deployed to
serve a cluster of IoT devices with the aim of maximizing the data collection
efficiency. However, the UAV was assumed to follow a circular trajectory and
as such, the flexible mobility of the UAV was not fully exploited for
performance optimization. To our best knowledge, the problem of jointly
optimizing SNs’ transmission scheduling and multi-antenna UAV’s trajectory has
not been rigorously studied in the literature yet, to fully exploit the
spatial multiplexing gain as well as the UAV mobility gain. This thus
motivates our current work to study this problem from an optimization
perspective.

Specifically, in this paper, we consider a UAV-enabled WSN, where a group of
single-antenna SNs send their independently sensed data to a multi-antenna UAV
in a periodic manner. To eliminate the inter-user interference among
transmitting SNs at each time instant, zero-forcing (ZF)-based receive
beamforming is adopted at the UAV. The main contributions of this paper are
summarized below.

  * •

First, we formulate an optimization problem to maximize the minimum data
collection rate from all SNs via jointly optimizing the UAV trajectory as well
as the transmission scheduling and power allocations of the SNs, subject to
the practical constraints on the maximum transmit power of the SNs and the
maximum flying speed of the UAV. The formulated problem is a mixed-integer
non-convex optimization problem, which is difficult to be optimally solved in
general.

  * •

Next, to tackle this problem and draw useful insight, we consider the special
case of this problem by ignoring the UAV’s maximum speed constraint. For this
relaxed problem, it is shown that the strong duality holds and thus it can be
optimally solved by employing the Lagrange duality method. It is revealed that
the UAV should hover at a finite number of locations with optimal hovering
durations, and this solution becomes asymptotically optimal when the UAV’s
flight duration and/or maximum speed becomes sufficiently large such that its
flying time is negligible as compared to its hovering time.

  * •

Then, we address the minimum-rate maximization problem for the general case
with the UAV speed constraint considered. First, we propose an efficient
initial trajectory for the UAV by solving an equivalent traveling salesman
problem (TSP), i.e., minimizing the UAV’s flying time to sequentially visit
all those hovering locations obtained via solving the previous relaxed
problem. Given the initial trajectory, we then propose a suboptimal solution
to the general problem by applying the block coordinate descent (BCD) and
successive convex approximation (SCA) techniques, through iteratively updating
the transmission scheduling/power allocations of the SNs and the trajectory of
the UAV.

  * •

Finally, extensive simulation results are provided to evaluate the performance
of the proposed design for multi-antenna UAV-enabled data harvesting. It is
shown that, as compared to the maximal ratio combining (MRC) scheme (that
exploits the beamforming gain only) and the single-antenna scheme, the
proposed scheme requires fewer hovering locations and significantly improves
the max-min rate. Furthermore, the UAV flight time required by the proposed
scheme to meet the same throughput requirement of all SNs is also drastically
reduced, thus greatly saving the UAV energy consumption.

The remainder of this paper is organized as follows. Section II introduces the
system model and problem formulation. Section III considers the special case
without UAV speed constraint and obtains the optimal solution to the
corresponding relaxed problem. Section IV proposes an efficient suboptimal
solution to the general problem with UAV speed constraint included. Finally,
numerical results are given in Section V, followed by the conclusions in
Section VI.

In this paper, scalars and vectors are represented by italic letters and
boldface lower-case letters, respectively. I and 0 denote an identity matrix
and an all-zero matrix, respectively, with appropriate dimensions. For a
square matrix S, $[\textbf{S}]_{k,k}$ denotes its $k$th diagonal element. For
a matrix M of arbitrary size, $\textbf{M}^{T}$ and $\textbf{M}^{H}$ denote its
transpose and conjugate transpose, respectively. $\mathbb{E}[\cdot]$ denotes
the statistical expectation. $||\cdot||$ denotes the $l_{2}$ norm. The
distribution of a circularly symmetric complex Gaussian (CSCG) random vector
with mean $\bm{x}$ and covariance matrix $\bm{\Sigma}$ is denoted by
$\mathcal{CN}(\bm{x},\bm{\Sigma})$, and $\sim$ stands for “distributed as”.

##  II System Model and Problem Formulation

###  II-A System Model

As shown in Fig. 1, we consider a UAV-enabled WSN, where a UAV equipped with
$M$ antennas is dispatched to collect data from $K$ distributed single-antenna
SNs. The set of SNs is denoted by $\mathcal{K}=\\{1,\cdots,K\\}$. We consider
the uplink transmission from SNs to the UAV for data collection, where the
results can be similarly extended to the downlink transmission as well.

We consider a 3D cartesian coordinate system, where the location of SN
$k\in\mathcal{K}$ is denoted as
$[\textbf{d}_{k}^{T},0]^{T}\in\mathbb{R}^{3\times 1}$ with
$\textbf{d}_{k}=[x_{k},y_{k}]^{T}$ denoting the horizontal coordinate. For
ease of exposition, the time discretization technique is applied, where the
time horizon $T$ for each periodic data collection operation of the UAV is
divided into $N$ time slots, each with equal length $\delta$, i.e.,
$t_{n}=n\delta$, $n=1,\cdots,N$ [4]. As such, the UAV trajectory is
approximated by a finite number of line segments with endpoints
$[\textbf{q}[n]^{T},z[n]]^{T}\in\mathbb{R}^{3\times 1}$, with
$\textbf{q}[n]=[x[n],y[n]]^{T}$ and $z[n]$ representing the horizontal and
vertical coordinates, respectively. Thus, the time-dependent distance between
the UAV and SN $k$ is expressed as

$\displaystyle d_{k}[n]=\sqrt{z[n]^{2}+||\textbf{q}[n]-\textbf{d}_{k}||^{2}},\
\ k\in\mathcal{K}.$ (1)

Let the baseband equivalent complex channel between SN $k$ and the UAV be
modelled as

$\displaystyle\textbf{h}_{k}[n]=\sqrt{\beta_{k}[n]}\textbf{g}_{k}[n],$ (2)

where $\beta_{k}[n]$ denotes the large-scale channel power gain over time due
to the distance-dependent path loss and shadowing, and $\textbf{g}_{k}[n]$
models the small-scale channel fading. Specifically, $\beta_{k}[n]$ is
modelled as

$\displaystyle\beta_{k}[n]=\beta_{0}d_{k}^{-\alpha}[n],$ (3)

where $\beta_{0}$ represents the channel power gain at the reference distance
of $d_{0}=1$ m, and $\alpha\geq 2$ is the path loss exponent. Furthermore,
$\textbf{g}_{k}[n]$ is modelled as the Rician fading with
$\mathbb{E}[||\textbf{g}_{k}[n]||^{2}]=1$ and

$\displaystyle\textbf{g}_{k}[n]=\sqrt{\frac{G}{G+1}}\bar{\textbf{g}}_{k}[n]+\sqrt{\frac{1}{G+1}}\tilde{\textbf{g}}_{k}[n],$
(4)

where $G$ is the Rician factor;
$\bar{\textbf{g}}_{k}[n]=[e^{j\theta_{k,1}}[n],\cdots,e^{j\theta_{k,M}}[n]]^{T}$
denotes the LoS channel component with $\theta_{k,m}[n]$ representing the
phase of the LoS path between SN $k$ and the $m$th antenna of the UAV;
$\tilde{\textbf{g}}_{k}[n]$ $\sim$ $\mathcal{CN}(\textbf{0},\textbf{I}_{M})$
denotes the Rayleigh fading channel component.

Moreover, define a binary variable $a_{k}[n]$, which indicates that SN $k$ is
scheduled for transmission to the UAV in time slot $n$ if $a_{k}[n]=1$ and
otherwise if $a_{k}[n]=0$. We then have the following constraints,

$\displaystyle a_{k}[n]\in\\{0,1\\},\ \ k\in\mathcal{K},\forall n.$ (5)

Define $\mathcal{K}_{n}=\\{k\in\mathcal{K}:a_{k}[n]=1\\}$ as the set of
transmitting SNs in time slot $n$, and $K_{n}=|\mathcal{K}_{n}|$. Then the
corresponding channel vectors are denoted by
$\textbf{H}[n]=[\textbf{h}_{\mathcal{K}_{n}(1)}[n],\cdots,\textbf{h}_{\mathcal{K}_{n}(K_{n})}[n]]$.
By denoting the transmit power of SN $k$ in time slot $n$ as $p_{k}[n]$, the
received signal at the UAV can be expressed as

$\displaystyle\textbf{y}[n]=\sum_{k\in\mathcal{K}_{n}}\sqrt{p_{k}[n]}\textbf{h}_{k}[n]s_{k}[n]+\textbf{z}[n]$
$\displaystyle=\textbf{H}[n]\textbf{P}[n]+\textbf{z}[n],~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
(6)

where $s_{k}[n]\sim\mathcal{CN}(0,1)$ denotes the signal sent by SN
$k\in\mathcal{K}_{n}$ in time slot $n$, and
$\textbf{z}[n]\sim\mathcal{CN}(\textbf{0},\sigma^{2}\textbf{I}_{M})$
represents the additive white Gaussian noise (AWGN) at the UAV receiver with
$\sigma^{2}$ denoting the noise power, and
$\textbf{P}[n]=[a_{\mathcal{K}_{n}(1)}[n]\sqrt{p_{\mathcal{K}_{n}(1)}[n]}s_{\mathcal{K}_{n}(1)}[n],\cdots,$
$a_{\mathcal{K}_{n}(K_{n})}[n]\sqrt{p_{\mathcal{K}_{n}(K_{n})}[n]}s_{\mathcal{K}_{n}(K_{n})}[n]]^{T}$.

With linear receive beamforming applied at the UAV, the processed signal is
given by

$\displaystyle\tilde{\textbf{s}}[n]=\textbf{W}[n]^{H}\textbf{y}[n],$ (7)

where
$\textbf{W}[n]=[\textbf{w}_{\mathcal{K}_{n}(1)}[n],\cdots,\textbf{w}_{\mathcal{K}_{n}(K_{n})}[n]]$,
with $\textbf{w}_{k}[n]\in\mathcal{C}^{M\times 1}$ denoting the beamforming
vector for extracting the signal of SN $k$ and $||\textbf{w}_{k}[n]||=1$,
$k\in\mathcal{K}_{n}$.

In particular, assuming the practical ZF beamforming, the inter-user
interference among SNs can be completely eliminated, i.e.,
$\textbf{w}_{k}[n]^{H}\textbf{h}_{l}[n]=0$, $l\neq k$,
$k,l\in\mathcal{K}_{n}$. Let

$\displaystyle\bar{\textbf{W}}[n]\triangleq[\bar{\textbf{w}}_{1}[n],\cdots,\bar{\textbf{w}}_{K_{n}}[n]]=\textbf{H}[n]\left(\textbf{H}[n]^{H}\textbf{H}[n]\right)^{-1}.$
(8)

Then the ZF beamforming vector for SN $k\in\mathcal{K}_{n}$ is given by

$\displaystyle\textbf{w}_{k}[n]=\frac{\bar{\textbf{w}}_{k}[n]}{||\bar{\textbf{w}}_{k}[n]||},\
\ k\in\mathcal{K}_{n},\forall n.$ (9)

Figure 1:  Multi-antenna UAV-enabled WSN.

As a result, the receive signal-to-noise ratio (SNR) for decoding the signal
from SN $k$ in time slot $n$ is given by

$\displaystyle\gamma_{k}[n]=\frac{p_{k}[n]|\textbf{w}_{k}[n]^{H}\textbf{h}_{k}[n]|^{2}}{\sigma^{2}},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle=\frac{p_{k}[n]}{\left[(\textbf{H}[n]^{H}\textbf{H}[n])^{-1}\right]_{k,k}\sigma^{2}},\
\ k\in\mathcal{K}_{n},\forall n.\\!\\!\\!\\!$ (10)

By averaging over the random small-scale channel fading, the achievable rate
of SN $k$ over the $n$th time slot in bits/second/Hertz (bps/Hz) is given by

$\displaystyle
R_{k}[n]=a_{k}[n]\mathbb{E}\left[\log_{2}\left(1+\frac{p_{k}[n]|\textbf{w}_{k}[n]^{H}\textbf{h}_{k}[n]|^{2}}{\sigma^{2}}\right)\right]$
$\displaystyle=a_{k}[n]\mathbb{E}\left[\log_{2}\left(1+\frac{p_{k}[n]}{\left[(\textbf{H}[n]^{H}\textbf{H}[n])^{-1}\right]_{k,k}\sigma^{2}}\right)\right],\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$
$\displaystyle k\in\mathcal{K}_{n},\forall n.$ (11)

The average achievable rate for SN $k$ over the entire time horizon is then
obtained as

$\displaystyle\bar{R}_{k}=\frac{1}{N}\sum_{n=1}^{N}R_{k}[n],\ \
k\in\mathcal{K}.$ (12)

Note that $R_{k}[n]$ in general does not admit a closed-form expression due to
the difficulty in handling the expectation operation in (II-A). Thus, in the
following, we derive a lower bound of $R_{k}[n]$ based on the results in [17].
Specifically, since both the LoS component and Rayleigh fading component in
the Rician channel model are independent over different SNs within each time
slot, a tight lower bound of $R_{k}[n]$ can be derived as

$\displaystyle\\!\\!\\!\\!\\!R_{k}[n]\geq
a_{k}[n]\log_{2}\left(1+\frac{p_{k}[n]}{\mathbb{E}\left[\left[(\textbf{H}[n]^{H}\textbf{H}[n])^{-1}\right]_{k,k}\right]\sigma^{2}}\right)~{}~{}~{}~{}~{}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$
$\displaystyle=a_{k}[n]\log_{2}\left(1+\frac{p_{k}[n]\beta_{0}d_{k}[n]^{-\alpha}}{\frac{\sigma^{2}}{M-K_{n}}}\right)~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$
$\displaystyle\\!\\!=a_{k}[n]\log_{2}\left(1+\frac{(M-K_{n})p_{k}[n]\gamma_{0}}{(z[n]^{2}+||\textbf{q}[n]-\textbf{d}_{k}||^{2})^{\alpha/2}}\right),$
(13)

where $\gamma_{0}\triangleq\beta_{0}/\sigma^{2}$. It should be noted that the
above rate lower bound of $R_{k}[n]$ only applies for $K_{n}\geq 2$ for the
feasibility of ZF beamforming assumed. While if $K_{n}=1$, we can simply
employ the optimal MRC-based receive beamforming at the UAV such that the
achievable rate can also be approximated similarly as (II-A) with the term
$(M-K_{n})$ replaced by $M$. Hence, for consistency, we combine the above two
cases and approximate $R_{k}[n]$ as

$\displaystyle\\!\\!r_{k}[n]\triangleq
a_{k}[n]\log_{2}\left(1+\frac{\kappa_{n}p_{k}[n]\gamma_{0}}{\left(z[n]^{2}+||\textbf{q}[n]-\textbf{d}_{k}||^{2}\right)^{\alpha/2}}\right),$
(14)

where

$\displaystyle\kappa_{n}=\left\\{\begin{array}[]{ll}M-K_{n},&K_{n}\geq 2,\\\
M,&K_{n}=1.\end{array}\right.$ (17)

Then, the average achievable rate is approximated as
$\bar{r}_{k}\triangleq\frac{1}{N}\sum_{n=1}^{N}r_{k}[n]$, $k\in\mathcal{K}$.
Numerical results in Section V will show that the above rate approximation
achieves good accuracy, especially for large $M$.

From (14), it is observed that the achievable rate of each transmitting SN in
each time slot (i.e., $a_{k}[n]=1$) depends on not only its distance from the
UAV, but also the spatial multiplexing gain of the UAV (i.e., $\kappa_{n}$).
Moreover, it is noted that to achieve the maximum transmission rate, the UAV
should always fly at the minimum possible altitude, i.e.,
$z[n]=H_{\mathrm{min}}$, $\forall n$, to minimize the distances to all SNs,
where $H_{\mathrm{min}}$ corresponds to the minimum altitude allowed to ensure
safety in practice.

###  II-B Problem Formulation

Our objective is to maximize the minimum average rate from all SNs by jointly
optimizing the UAV trajectory $\\{\textbf{q}[n]\\}$, SNs’ transmission
scheduling $\\{a_{k}[n]\\}$ and power allocations $\\{p_{k}[n]\\}$. By
defining $r\triangleq\min_{k\in\mathcal{K}}\bar{r}_{k}$, this optimization
problem is formulated as

$\displaystyle\mathrm{(P1)}~{}~{}\max_{\begin{subarray}{c}\\{\textbf{q}[n]\\},\\{a_{k}[n]\\},\\\
\\{p_{k}[n]\\},r\end{subarray}}r~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\mathrm{s.t.}~{}~{}\frac{1}{N}\sum_{n=1}^{N}r_{k}[n]\geq r,\ \
k\in\mathcal{K},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$ (18a)
$\displaystyle\sum_{k=1}^{K}a_{k}[n]\leq M,\ \ \forall
n,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$ (18b)
$\displaystyle a_{k}[n]\in\\{0,1\\},\ \ k\in\mathcal{K},\forall
n,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$ (18c)
$\displaystyle\frac{1}{N}\sum_{n=1}^{N}p_{k}[n]\leq\bar{P},\ \
k\in\mathcal{K},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$ (18d)
$\displaystyle p_{k}[n]\geq 0,\ \ k\in\mathcal{K},\forall
n,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$ (18e)
$\displaystyle||\textbf{q}[n+1]-\textbf{q}[n]||\leq V_{\mathrm{h}},\ \
n=1,\cdots,N-1,\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$ (18f)
$\displaystyle\textbf{q}[1]=\textbf{q}_{\mathrm{I}},\ \
\textbf{q}[N]=\textbf{q}_{\mathrm{F}},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\!$
(18g)

where $\bar{P}$ is the average power limit at each SN, and
$V_{\mathrm{h}}\triangleq v_{\mathrm{h}}\delta$ with $v_{\mathrm{h}}$
representing the maximum horizontal speed of the UAV,
$\textbf{q}_{\mathrm{I}}$ and $\textbf{q}_{\mathrm{F}}$ are the initial and
final locations of the UAV, respectively.

Note that problem (P1) is challenging to solve due to the following reseasons.
Firstly, $r_{k}[n]$ in constraint (18a) is not jointly concave with respect to
the optimization variables $\\{\textbf{q}[n]\\}$, $\\{a_{k}[n]\\}$ and
$\\{p_{k}[n]\\}$. Secondly, even with fixed trajectory $\\{\textbf{q}[n]\\}$,
problem (P1) is still a mixed-integer non-linear programming as the binary
variable $a_{k}[n]$ is in general coupled with $\kappa_{n}$ and $p_{k}[n]$ in
$r_{k}[n]$ given by (14). Consequently, problem (P1) is a mixed-integer non-
convex optimization problem, which is difficult to solve in general. To tackle
this problem, we first consider the special case of (P1) by ignoring the UAV
maximum speed constraint in (18f) as well as the initial/final location
constraints in (18g). The relaxed problem of (P1) is thus given by

$\displaystyle\mathrm{(P2)}~{}~{}\max_{\begin{subarray}{c}\\{\textbf{q}[n]\\},\\{a_{k}[n]\\},\\\
\\{p_{k}[n]\\},r\end{subarray}}r~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\mathrm{s.t.}~{}~{}\eqref{p1001}-\eqref{p1005}.~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$

In the following, we first find the optimal solution to (P2) in Section III.
Then based on the optimal solution obtained, we propose an efficient solution
to solve the general problem (P1) sub-optimally in Section IV.

##  III Optimal Solution to (P2)

In this section, we solve the relaxed problem (P2), which is still challenging
due to the non-convex constraints and binary variables. Fortunately, it can be
verified that (P2) satisfies the so-called time-sharing condition in [18], so
that the strong duality holds between (P2) and its Lagrange duality problem.
Therefore, (P2) can be optimally solved by the Lagrange duality method.

Specifically, the partial Lagrangian of (P2) is given by

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mathcal{L}_{1}(\\{\textbf{q}[n]\\},\\{p_{k}[n]\\},\\{a_{k}[n]\\},r,\bm{\lambda},\bm{\mu})=~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\left(1-\sum_{k=1}^{K}\lambda_{k}\right)r+\sum_{k=1}^{K}\lambda_{k}\bar{r}_{k}+\sum_{k=1}^{K}\mu_{k}\left(N\bar{P}-\sum_{n=1}^{N}p_{k}[n]\right),\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$

where $\bm{\lambda}\triangleq\\{\lambda_{k}\\}$ and
$\bm{\mu}\triangleq\\{\mu_{k}\\}$ are the dual variables associated with the
constraints in (18a) and (18d), respectively. The Lagrange dual function of
(P2) is then given by

$\displaystyle\\!\\!\\!\\!\\!\\!g_{1}(\bm{\lambda},\bm{\mu})=~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\\!\\!\\!\left\\{\begin{aligned}
\max_{\begin{subarray}{c}\\{\textbf{q}[n]\\},\\{p_{k}[n]\\},\\\
\\{a_{k}[n]\\},r\end{subarray}}\mathcal{L}_{1}\left(\\{\textbf{q}[n]\\},\\{p_{k}[n]\\},\\{a_{k}[n]\\},r,\bm{\lambda},\bm{\mu}\right)\\\
\mathrm{s.t.}~{}~{}~{}~{}~{}\eqref{p1002},~{}\eqref{p1003},~{}\eqref{p1005}.~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\end{aligned}\right.\\!\\!\\!\\!\\!\\!\\!\\!$
(21)

For $g_{1}(\bm{\lambda},\bm{\mu})$ in (III) to be bounded, we should have
$1-\sum_{k=1}^{K}\lambda_{k}=0$. Therefore, the dual problem of (P2) is given
by

$\displaystyle\mathrm{(D2)}~{}~{}\min_{\bm{\lambda},\bm{\mu}}~{}~{}~{}g_{1}(\bm{\lambda},\bm{\mu})~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\mathrm{s.t.}~{}~{}~{}1-\sum_{k=1}^{K}\lambda_{k}=0,~{}~{}~{}~{}~{}~{}~{}~{}\\!~{}~{}~{}~{}~{}$
(22a) $\displaystyle\lambda_{k}\geq 0,\ \ \forall
k\in\mathcal{K},~{}~{}~{}\\!~{}~{}~{}~{}~{}~{}~{}$ (22b)
$\displaystyle\mu_{k}\geq 0,\ \ \forall
k\in\mathcal{K}.~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\!\\!$ (22c)

Since the strong duality holds between (P2) and (D2), we can solve (P2) by
equivalently solving (D2). Let the feasible set of $\bm{\lambda}$ and
$\bm{\mu}$ characterized by the constraints in (22a)-(22c) be denoted as
$\mathcal{X}_{1}$. In the following, we first obtain
$g_{1}(\bm{\lambda},\bm{\mu})$ by solving problem (III) under any given
$(\bm{\lambda},\bm{\mu})\in\mathcal{X}_{1}$, and then solve (D2) to find the
optimal $(\bm{\lambda},\bm{\mu})$ to minimize $g_{1}(\bm{\lambda},\bm{\mu})$,
and finally construct the optimal solution to (P2).

###  III-A Obtaining $g_{1}($$\bm{\lambda}$$,~{}\\!\\!$$\bm{\mu}$$)$ by
solving problem (III) for given $($$\bm{\lambda}$$,~{}\\!\\!$$\bm{\mu}$$)$
$\in\mathcal{X}_{1}$

For any given $(\bm{\lambda},\bm{\mu})\in\mathcal{X}_{1}$, problem (III) is
still a non-convex optimization problem. In order to solve this problem, we
first focus on the case with given SNs’ scheduling $\\{a_{k}[n]\\}$ and
trajectory $\\{\textbf{q}[n]\\}$ to optimize the power allocations
$\\{p_{k}[n]\\}$ only, which is reformulated as

$\displaystyle\max_{\\{p_{k}[n]\\}}~{}~{}\sum_{k=1}^{K}\lambda_{k}\bar{r}_{k}+\sum_{k=1}^{K}\mu_{k}\left(N\bar{P}-\sum_{n=1}^{N}p_{k}[n]\right)$
(23) $\displaystyle\mathrm{s.t.}~{}~{}~{}p_{k}[n]\geq 0,\ \
k\in\mathcal{K},\forall n.~{}~{}~{}~{}~{}~{}~{}~{}~{}\\!\\!~{}~{}~{}~{}~{}~{}$
(23a)

As it can be verified that the objective of problem (23) is a concave function
with respect to $p_{k}[n]$, problem (23) is a convex optimization problem.
Then by applying the standard Lagrange duality method and the Karush-Kuhn-
Tucker (KKT) conditions, it can be shown that the optimal solution to (23)
follows the classic water-filling (WF) structure given by

$\displaystyle
p_{k}^{*}[n]=\left[\frac{\lambda_{k}a_{k}[n]}{N\mu_{k}\log(2)}-\frac{d_{k}[n]^{\alpha}}{\kappa_{n}\gamma_{0}}\right]^{+},~{}~{}~{}~{}~{}~{}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}a_{k}[n]\left[\frac{\lambda_{k}}{N\mu_{k}\log(2)}-\frac{d_{k}[n]^{\alpha}}{\kappa_{n}\gamma_{0}}\right]^{+},\\!$
(24)

where $(a)$ holds since $a_{k}[n]\in\\{0,1\\}$, and $[b]^{+}=\max\\{b,0\\}$.
As a result, the achievable rate is expressed as

$\displaystyle
r^{*}_{k}[n]=a_{k}[n]\left[\log_{2}\left(\frac{\lambda_{k}a_{k}[n]\kappa_{n}\gamma_{0}}{N\mu_{k}\log(2)d_{k}[n]^{\alpha}}\right)\right]^{+}~{}$
$\displaystyle=a_{k}[n]\left[\log_{2}\left(\frac{\lambda_{k}\kappa_{n}\gamma_{0}}{N\mu_{k}\log(2)d_{k}[n]^{\alpha}}\right)\right]^{+},$
(25)

and thus $\bar{r}^{*}_{k}=\frac{1}{N}\sum_{n=1}^{N}r^{*}_{k}[n]$. With the
obtained $\\{p_{k}^{*}[n]\\}$ and $\\{r^{*}_{k}[n]\\}$, problem (III) is
further recast to

$\displaystyle\max_{\\{\textbf{q}[n]\\},\\{a_{k}[n]\\}}\sum_{k=1}^{K}\lambda_{k}\bar{r}^{*}_{k}+\sum_{k=1}^{K}\mu_{k}\left(N\bar{P}-\sum_{n=1}^{N}p^{*}_{k}[n]\right)$
(26)
$\displaystyle\mathrm{s.t.}~{}~{}~{}~{}~{}\eqref{p1002},~{}\eqref{p1003}.~{}~{}~{}~{}~{}\\!~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$

It is worth noting that (26) consists of $N$ sub-problems, each corresponding
to one time slot. As all sub-problems are identical across different time
slots, we can drop the time slot index $n$ and re-express each problem as

$\displaystyle\max_{\textbf{q},\\{a_{k}\\},\kappa}~{}~{}\sum_{k=1}^{K}a_{k}f_{k}(\textbf{q},\kappa)+N\bar{P}\sum_{k=1}^{K}\mu_{k}~{}~{}~{}$
(27)
$\displaystyle\mathrm{s.t.}~{}~{}~{}~{}~{}~{}\eqref{p1002},\eqref{p1003},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$

where

$\displaystyle\\!\\!\\!f_{k}(\textbf{q},\kappa)\triangleq\lambda_{k}\left[\log_{2}\left(\frac{\lambda_{k}\kappa\gamma_{0}}{N\mu_{k}\log(2)d_{k}^{\alpha}}\right)\right]^{+}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle-N\mu_{k}\left[\frac{\lambda_{k}}{N\mu_{k}\log(2)}-\frac{d_{k}^{\alpha}}{\kappa\gamma_{0}}\right]^{+}.$
(28)

The constant term in the objective of (27) has been omitted for brevity. Note
that (27) is a mixed-integer optimization problem, which is non-convex and
difficult to solve in general. In the following, we optimally solve it via a
two-dimensional (2D) search.

To start with, we consider the box region
$[\underline{x},\overline{x}]^{T}\times[\underline{y},\overline{y}]^{T}$,
where $\underline{x}=\min_{k\in\mathcal{K}}x_{k}$,
$\overline{x}=\max_{k\in\mathcal{K}}x_{k}$,
$\underline{y}=\min_{k\in\mathcal{K}}y_{k}$,
$\overline{y}=\max_{k\in\mathcal{K}}y_{k}$, and discretize it with a finite
granularity $\Delta_{g}$. Note that the value of $\Delta_{g}$ is typically
chosen to be sufficiently small to ensure certain accuracy. Next, we focus on
finding $\\{a_{k}\\}$ and $\kappa$ at each discretized location q to maximize
the objective value in (27) and then compare them to find the maximum one.

Note that with fixed q and $\kappa$, the values of
$\\{f_{k}(\textbf{q},\kappa)\\}$ are also determined. Therefore, to maximize
the objective value in (27) with given q and $\kappa$, we only need to choose
the $(M-\kappa)$ largest values (or choose the largest value when $\kappa=M$)
among $\\{f_{k}(\textbf{q},\kappa)\\}$ and set the corresponding scheduling
variables $\\{a_{k}\\}=1$. Then by comparing the objective value achieved by
each possible value of $\kappa$, where
$\kappa=M,M-2,M-3,\cdots,M-K_{\mathrm{max}}$ with
$K_{\mathrm{max}}\triangleq\min\\{M-1,K\\}$, we are able to find the optimal
value of $\kappa$ and the corresponding $\\{a_{k}\\}$ under the given location
q. Finally, by comparing the objective values achieved at each discretized
location, we can find the optimal solution to problem (27), which is given by

$\displaystyle\textbf{q}^{\\{\bm{\lambda},\bm{\mu}\\}}[n]=\textbf{q}^{\\{\bm{\lambda},\bm{\mu}\\}},\
\
\\{a^{\\{\bm{\lambda},\bm{\mu}\\}}_{k}[n]\\}=\\{a^{\\{\bm{\lambda},\bm{\mu}\\}}_{k}\\},$
$\displaystyle\kappa^{\\{\bm{\lambda},\bm{\mu}\\}}_{n}=\kappa^{\\{\bm{\lambda},\bm{\mu}\\}},\
\ \forall n.$ (29)

It is worth pointing out that although the above optimal solution to (27) is
generally non-unique, we can arbitrarily choose one of them to obtain the dual
function $g_{1}(\bm{\lambda},\bm{\mu})$ in (D2). Note that as the optimal dual
solution is generally non-unique, the solution obtained for (D2) may not be
optimal for (P2) after solving the dual problem; thus, an additional step is
usually needed to construct the optimal solution to (P2) based on the dual
optimal solution, as will be shown in Section III-C.

###  III-B Finding optimal $\bm{\lambda}$ and $\bm{\mu}$ to solve (D2)

In the following, we search over $\bm{\lambda}$ and $\bm{\mu}$ to minimize
$g_{1}(\bm{\lambda},\bm{\mu})$ for solving (D2). Since the dual problem (D2)
is always convex but non-differentiable in general, we employ the subgradient
based method, such as the ellipsoid method [19], to obtain the optimal dual
solution $\bm{\lambda}$ and $\bm{\mu}$, which are denoted by
$\bm{\lambda}^{*}$ and $\bm{\mu}^{*}$. In each iteration, the dual variables
$\bm{\lambda}$ and $\bm{\mu}$ are updated based on the objective and
constraints in (D2). Specifically, the subgradients with respect to
$(\bm{\lambda},\bm{\mu})$ are given by

$\displaystyle\\!\Delta\lambda_{k}=r_{k}^{*},~{}~{}~{}~{}~{}\\!~{}~{}~{}~{}~{}~{}$
(30) $\displaystyle\Delta\mu_{k}=N\bar{P}-Np_{k}^{*},$ (31)

where

$\displaystyle
r_{k}^{*}=a_{k}\left[\log_{2}\left(\frac{\lambda_{k}^{*}\kappa\gamma_{0}}{N\mu^{*}_{k}\log(2)d_{k}^{\alpha}}\right)\right]^{+},$
(32) $\displaystyle
p_{k}^{*}=a_{k}\left[\frac{\lambda^{*}_{k}}{N\mu^{*}_{k}\log(2)}-\frac{d_{k}^{\alpha}}{\kappa\gamma_{0}}\right]^{+}.~{}~{}~{}~{}$
(33)

###  III-C Constructing optimal solution to (P2)

Based on the optimal dual solution $\bm{\lambda}^{*}$ and $\bm{\mu}^{*}$, it
remains to obtain the optimal (primal) solution to (P2), which is denoted as
$\\{\textbf{q}^{*}[n]\\}$, $\\{p^{*}_{k}[n]\\}$, and $\\{a^{*}_{k}[n]\\}$. It
is worth pointing out that when the Lagrange duality method is employed to
solve (P2) via (D2), the optimal solution to problem (III) under the optimal
dual solution $\bm{\lambda}^{*}$ and $\bm{\mu}^{*}$ (i.e.,
$\\{\textbf{q}^{\bm{\lambda}^{*},\bm{\mu}^{*}}[n]\\}$,
$\\{p^{\bm{\lambda}^{*},\bm{\mu}^{*}}_{k}[n]\\}$, and
$\\{a_{k}^{\bm{\lambda}^{*},\bm{\mu}^{*}}[n]\\}$) is the optimal solution to
(P2), if such a solution is feasible [20]. On the other hand, when the
obtained $\\{\textbf{q}^{\bm{\lambda}^{*},\bm{\mu}^{*}}[n]\\}$,
$\\{p^{\bm{\lambda}^{*},\bm{\mu}^{*}}_{k}[n]\\}$ and
$\\{a_{k}^{\bm{\lambda}^{*},\bm{\mu}^{*}}[n]\\}$ for problem (III) are non-
unique, they may not be feasible nor optimal to problem (P2) in general. In
this case, additional procedures are needed to construct the optimal solution
to (P2) by applying time-sharing over these non-unique solutions.

Specifically, with the optimal dual solution $\bm{\lambda}^{*}$ and
$\bm{\mu}^{*}$, suppose that problem (26) has a total number of $\Omega$
solutions, denoted by $\\{\textbf{q}^{*}_{\omega}\\}$,
$\\{p^{*}_{k,\omega}\\}$, and $\\{a^{*}_{k,\omega}\\}$,
$\omega=1,\cdots,\Omega$. The corresponding average rate is then expressed as
$\bar{r}^{*}_{k,\omega}$. Due to the zero duality gap between (P2) and (D2),
it is obvious that at the optimal solution for (P2), the UAV should choose
from the $\Omega$ locations to hover above in each time slot. Since the
achievable rate of SNs remains the same when the UAV stays at the same
hovering location in different time slots, the optimal solution to (P2) can be
constructed via allocating the flight duration $T$ over the $\Omega$ hovering
locations obtained. Let $\tau_{\omega}$ denote the hovering duration of the
UAV above the location $\textbf{q}^{*}_{\omega}$. The extra problem to
maximize the minimum average rate of SNs is formulated as

$\displaystyle\mathrm{(P2.1)}~{}~{}\max_{\\{\tau_{\omega}\\},r}~{}~{}r~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\!~{}~{}~{}~{}~{}$
$\displaystyle\mathrm{s.t.}~{}~{}\frac{1}{T}\sum_{\omega=1}^{\Omega}\tau_{\omega}\bar{r}^{*}_{k,\omega}\geq
r,\ \ k\in\mathcal{K},~{}$ (34a)
$\displaystyle\sum_{\omega=1}^{\Omega}\tau_{\omega}=T,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
(34b) $\displaystyle\tau_{\omega}\geq 0,\ \
\omega=1,\cdots,\Omega.~{}~{}~{}\\!~{}~{}~{}$ (34c)

Problem (P2.1) is a linear programming (LP), which can be efficiently solved
via the standard convex optimization technique [21]. It should be noted that
after solving (P2.1), there might exist some hovering locations with
$\tau_{\omega}=0$, which indicates that these hovering locations are not
required for the optimal solution to (P2). Therefore, we should only choose
those locations with the corresponding $\tau_{\omega}>0$ and the number of
such locations is denoted as $\Omega^{*}$.

In summary, the details to solve (P2) are presented in Algorithm 1. The
computational complexity of Algorithm 1 mainly consists of three parts. The
first one is the exhaustive search over the region in step 3 with complexity
$O(K_{\mathrm{max}}\Delta_{G})$, where
$\Delta_{G}\triangleq(\bar{x}-\underline{x})(\bar{y}-\underline{y})/\Delta_{g}^{2}$.
The second part is for updating the dual variables via the ellipsoid method in
step 4 with complexity $O(K^{2})$. Since the ellipsoid method takes $O(K^{2})$
to converge [19], the overall complexity from step 2 to step 5 is
$O((K_{\mathrm{max}}\Delta_{G}+K^{2})K^{2})$. The third part is due to step 7
for solving the LP with complexity $O(\Omega^{3})$. Due to the much smaller
value of $\Omega^{3}$ compared to $K^{4}$, the overall complexity of Algorithm
1 is approximately $O((K_{\mathrm{max}}\Delta_{G}+K^{2})K^{2})$.

Algorithm 1 Proposed Algorithm for Solving (P2).

1: Initialize an ellipsoid $\bm{\epsilon}((\bm{\lambda},\bm{\mu}),\textbf{A})$
containing $(\bm{\lambda}^{*},\bm{\mu}^{*})$,
where $(\bm{\lambda},\bm{\mu})$ is the center point of the ellipsoid, and A is
a positive definite matrix that characterizes its size.
2: Repeat
3:  Obtain $\\{\textbf{q}^{\bm{\lambda},\bm{\mu}}\\}$,
$\\{p_{k}^{\bm{\lambda},\bm{\mu}}\\}$, and
$\\{a_{k}^{\bm{\lambda},\bm{\mu}}\\}$ via 2D search over
the region
$[\underline{x},\overline{x}]^{T}\times[\underline{y},\overline{y}]^{T}$.
4: Compute the subgradients of $g_{1}(\bm{\lambda},\bm{\mu})$, and update
$(\bm{\lambda},\bm{\mu})$
by the ellipsoid method.
5: Until $\bm{\lambda}$ and $\bm{\mu}$ converge within a given accuracy.
6: Set $(\bm{\lambda}^{*},\bm{\mu}^{*})\leftarrow(\bm{\lambda},\bm{\mu})$.
7: Obtain the optimal solution to (P2) via solving (P2.1).

Note that the above results demonstrate that the UAV should hover above a
finite number of locations with optimal allocations of the hovering duration
to maximize the minimum average rate of all SNs by solving problem (P2)
optimally without considering the UAV speed and initial/final location
constraints111Note that the initial/final UAV location constraints can be
added back to (P2) without changing the above optimal solution, since if they
do not belong to the optimal hovering locations, we can simply assign zero
hovering time for them without loss of generality. In practice, even with
finite speed constraint for the UAV, when the given time horizon $T$ or the
speed of the UAV $v_{\mathrm{h}}$ is sufficiently large such that the flying
time of the UAV becomes negligible as compared its hovering time, the above
solution can be shown to be asymptotically optimal to problem (P1) [22].

##  IV Proposed Solution to (P1)

In this section, we consider the general problem (P1) with the UAV speed and
initial/final location constraints included, i.e., the flying time of the UAV
as well as the communications with the SNs during its flying cannot be
ignored. Although the results obtained in Section III are not directly
applicable to the general case of (P1), the max-min rate achieved by solving
problem (P2) can serve as an upper bound for that achievable by solving
problem (P1).

In the following, we propose an efficient suboptimal solution to (P1) based on
the optimal solution obtained for (P2). Specifically, we firstly design an
efficient initial trajectory for the UAV based on the classic TSP, where the
UAV sequentially visits the $\Omega^{*}$ locations obtained by solving (P2)
with the minimum flying time. Then, we find the optimal transmission
scheduling and power allocations of SNs with the given UAV trajectory and
iteratively optimize them with the UAV trajectory by applying the BCD and SCA
techniques, until all of them get converged.

###  IV-A TSP-based initial trajectory

In this subsection, we design a TSP-based initial trajectory, which minimizes
the flying time of the UAV to sequentially visit those $\Omega^{*}$ hovering
locations obtained in Section III. Note that the proposed TSP-based trajectory
is similar to the successive hover-and-fly trajectory proposed in [7], [23,
24, 25].

To maximize the hovering time of the UAV above the obtained hovering
locations, the UAV should always fly with the maximum speed $v_{\mathrm{h}}$
among these locations, and the flying time, or equivalently the flying
distance, needs to be minimized. By including the initial/final location
constraints in (18g), the problem of minimizing the UAV flying distance is
reminiscent of the classic TSP, which can be efficiently solved via the
techniques outlined in [25], [26]. After solving the corresponding TSP, we
obtain the minimum flying time required, denoted by $T_{\mathrm{tsp}}$, to
visit all the hovering locations, and the permutation order
$\bm{\pi}\triangleq[\pi(1),\cdots,\pi(\Omega^{*})]$, with $\pi(\omega)$
representing the index of the $\omega$th hovering location to be visited. In
the following, for any given flight duration $T$, we design the initial
trajectory depending on whether $T$ is greater than $T_{\mathrm{tsp}}$ or not.

####  IV-A1 $T\geq T_{\mathrm{tsp}}$

In this case, the UAV is able to reach all the hovering locations. With the
permutation order $\bm{\pi}$ and flying time $T_{\mathrm{tsp}}$, the
trajectory design reduces to allocating the remaining time
$T-T_{\mathrm{tsp}}$ among all hovering locations, which can be efficiently
obtained via proportionally allocating the time according to
$\\{\tau_{\omega}\\}$ obtained in (P2.1). We denote the above obtained
trajectory as $\\{\textbf{q}^{0}[n]\\}_{n=1}^{N}$.

####  IV-A2 $T<T_{\mathrm{tsp}}$

In this case, it is infeasible for the UAV to visit all the hovering locations
within the given time $T$. To design a feasible trajectory, we specify a disk-
shaped region centered at each hovering location $\\{\textbf{q}_{\omega}\\}$
with the same radius $R$. The main idea here is to find the UAV trajectory and
radius $R$ such that the UAV is able to reach each disk region. This problem
is equivalent to the traveling salesman problem with neighborhoods (TSPN), for
which an efficient suboptimal solution can be found in [25]. For notational
convenience, we also denote the above trajectory obtained as
$\\{\textbf{q}^{0}[n]\\}_{n=1}^{N}$.

Given the initial UAV trajectory, the transmission scheduling and power
allocations of SNs can be optimized. Then the UAV trajectory can be further
optimized jointly with the transmission scheduling and power allocations. In
the following, we consider two sub-problems of (P1), namely, optimizing SNs’
transmission scheduling and power allocations with fixed UAV trajectory, and
optimizing UAV trajectory with fixed SNs’ scheduling and power allocations.
Finally, based on the solutions obtained for these two sub-problems, an
iterative algorithm is proposed to alternately optimize these two blocks of
variables until they both converge.

###  IV-B Transmission scheduling and power allocations optimization with
fixed trajectory

In this subsection, we consider the sub-problem to optimize the SNs’
transmission scheduling $\\{a_{k}[n]\\}$ and power allocations
$\\{p_{k}[n]\\}$ with fixed UAV trajectory $\\{\textbf{q}[n]\\}$. This problem
is given by

$\displaystyle\mathrm{(P3)}~{}~{}\max_{\\{p_{k}[n]\\},\\{a_{k}[n]\\},r}~{}~{}r~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\\!\\!\\!\\!\mathrm{s.t.}~{}~{}\frac{1}{N}\sum_{n=1}^{N}a_{k}[n]\log_{2}\left(1+\frac{\kappa_{n}p_{k}[n]\gamma_{0}}{d_{k}[n]^{\alpha}}\right)\geq
r,\ \ k\in\mathcal{K},\\!\\!\\!\\!\\!$ (35a)
$\displaystyle\sum_{k=1}^{K}a_{k}[n]\leq M,\ \ \forall
n,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
(35b) $\displaystyle a_{k}[n]\in\\{0,1\\},\ \ k\in\mathcal{K},\forall
n,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$ (35c)
$\displaystyle\frac{1}{N}\sum_{n=1}^{N}p_{k}[n]\leq\bar{P},\ \
k\in\mathcal{K},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
(35d) $\displaystyle p_{k}[n]\geq 0,\ \ k\in\mathcal{K},\forall
n.~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
(35e)

Although (P3) is also a mixed-integer non-convex optimization problem, it can
be verified that it satisfies the time-sharing condition in [18] such that it
can be optimally solved using the Lagrange duality method, which is similar to
problem (P2). The details for solving (P3) are given in Appendix.

###  IV-C Trajectory optimization with fixed transmission scheduling and power
allocations

In this subsection, we consider the other sub-problem to further optimize the
UAV trajectory $\\{\textbf{q}[n]\\}$ given the SNs’ transmission scheduling
$\\{a_{k}[n]\\}$ and power allocations $\\{p_{k}[n]\\}$ obtained by solving
(P3). The problem is formulated as

$\displaystyle\mathrm{(P4)}~{}~{}\max_{\\{\textbf{q}[n]\\},r}~{}~{}r~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mathrm{s.t.}~{}~{}\frac{1}{N}\sum_{n=1}^{N}\log_{2}\left(1+\frac{\varepsilon_{k}[n]}{(z[n]^{2}+||\textbf{q}[n]-\textbf{d}_{k}||^{2})^{\alpha/2}}\right)\geq
r,\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!~{}$
$\displaystyle k\in\mathcal{K},$ (36a)
$\displaystyle||\textbf{q}[n+1]-\textbf{q}[n]||\leq V_{\mathrm{h}},\ \
n=1,\cdots,N-1,\\!\\!\\!$ (36b)
$\displaystyle\textbf{q}[1]=\textbf{q}_{\mathrm{I}},\ \
\textbf{q}[N]=\textbf{q}_{\mathrm{F}},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\!~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
(36c)

where $\varepsilon_{k}[n]\triangleq\kappa_{n}p_{k}[n]\gamma_{0}$. Although
problem (P4) is a non-convex optimization problem due to the non-convex
constraints in (36), an efficient suboptimal solution can be found by applying
the SCA technique with any given initial UAV trajectory (e.g.,
$\\{\textbf{q}^{0}[n]\\}$ obtained in Section IV-A). To this end, we need the
following result.

######  Lemma 1.

For any given trajectory $\\{\textbf{q}^{j}[n]\\}$, we have

$\displaystyle\\!\\!\\!\\!\\!\\!r_{k}[n]\geq
r^{lb}_{k}[n]~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\triangleq
a_{k}[n]\log_{2}\left(1+\frac{\varepsilon_{k}[n]}{(z[n]^{2}+||\textbf{q}^{j}[n]-\textbf{d}_{k}||^{2})^{\alpha/2}}\right)\\!\\!\\!\\!\\!\\!$
$\displaystyle-\vartheta^{j}_{k}[n]\left(||\textbf{q}[n]-\textbf{d}_{k}||^{2}-||\textbf{q}^{j}[n]-\textbf{d}_{k}||^{2}\right),~{}~{}$
(37)

where

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\vartheta^{j}_{k}[n]\triangleq~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\frac{a_{k}[n](\log_{2}e)\varepsilon_{k}[n](\alpha/2)}{\left(z[n]^{2}+||\textbf{q}^{j}[n]-\textbf{d}_{k}||^{2}\right)\left(\left(z[n]^{2}+||\textbf{q}^{j}[n]-\textbf{d}_{k}||^{2}\right)^{\alpha/2}+\varepsilon_{k}[n]\right)}.\\!\\!\\!\\!\\!\\!\\!\\!\\!$
(38)

###### Proof:

The proof is similar to that given in [4], [5], and thus omitted for brevity.
∎

As a result, a lower bound of the optimal value to problem (P4) with fixed
SNs’ transmission scheduling and power allocations can be obtained via solving
the following problem

$\displaystyle\mathrm{(P4.1)}~{}~{}\max_{\\{\textbf{q}[n]\\},r}~{}~{}r~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\mathrm{s.t.}~{}~{}\frac{1}{N}\sum_{n=1}^{N}r^{lb}_{k}[n]\geq
r,\ \ k\in\mathcal{K},~{}~{}~{}~{}$ (39)
$\displaystyle\eqref{p31002},~{}\eqref{p31003}.~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$

Problem (P4.1) can be verified to be convex, which can be efficiently solved
by CVX [21].

Finally, based on the solutions obtained to the above two sub-problems, the
overall algorithm to solve problem (P1) is summarized in Algorithm 2. Although
the TSP involved in step 2 is an NP-hard problem, it can be efficiently solved
with high accuracy by existing algorithms with polynomial complexity. In step
5, the complexity for obtaining the scheduling and power allocation solution
is similar to that of Algorithm 1, which is
$O((NK_{\mathrm{max}}\Delta_{G}+K^{2})K^{2})$. The complexity for updating the
trajectory in step 6 is $O(N^{3.5})$ [25]. By denoting the number of
iterations needed for convergence as $J$, the total complexity of Algorithm 2
is $O(J((NK_{\mathrm{max}}\Delta_{G}+K^{2})K^{2}+N^{3.5}))$.

Algorithm 2 Proposed Algorithm for Solving (P1).

1: Solve problem (P2) with Algorithm 1 to obtain the hovering
locations $\\{\textbf{q}^{*}_{\omega}\\}$, scheduling
$\\{a^{*}_{k,\omega}\\}$, and power allocations
$\\{p^{*}_{k,\omega}\\}$, $\omega=1,\cdots,\Omega^{*}$.
2: Solve the TSP to obtain the permutation order $\bm{\pi}$ and
flying time $T_{\mathrm{tsp}}$ to visit all the above hovering locations.
3: Construct the initial trajectory $\\{\textbf{q}^{0}[n]\\}$ by comparing $T$
and $T_{\mathrm{tsp}}$ according to Section IV-A. Let $j$$=$$0$.
4: Repeat
5: For given trajectory $\\{\textbf{q}^{j}[n]\\}$, obtain the optimal
transmission scheduling $\\{a^{j+1}_{k}[n]\\}$ and power allocations
$\\{p^{j+1}_{k}[n]\\}$ based on the algorithm given in Appendix.
6:  For given $\\{a^{j+1}_{k}[n]\\}$ and $\\{p^{j+1}_{k}[n]\\}$, update the
trajectory $\\{\textbf{q}^{j+1}[n]\\}$ via solving (P4.1).
7: Update $j=j+1$.
8: Until the objective value $r$ converges within a given
accuracy.

##  V Simulation Results

In this section, numerical results are provided to evaluate the performance of
the proposed design for multi-antenna UAV data harvesting. We consider a WSN
with $K=8$ SNs that are distributed in a square area with side length equal to
1000 m. The simulation results are based on one realization of SNs’ locations,
as shown in Fig. 3. The channel power gain at the reference distance $d_{0}=1$
m is set as $\beta_{0}=-60$ dB. The total available bandwidth is $B=0.1$ MHz
and the noise power spectrum density is $N_{0}=-154$ dBm/Hz such that the
noise power $\sigma^{2}=BN_{0}=-104$ dBm. The average power limit at all SNs
is set as $\bar{P}=0.01$ W. The path loss exponent is set as $\alpha=2$ and
the Rician factor is $G=0.94$ [27]. The UAV flying altitude is set as
$H_{\mathrm{min}}=130$ m, and the initial and final locations are
$\textbf{q}_{\mathrm{I}}=[400,0]^{T}$ m and
$\textbf{q}_{\mathrm{F}}=[1000,500]^{T}$ m, respectively. Furthermore, the
maximum UAV speed is $v_{\mathrm{h}}=20$ m/s, and the length of time slot is
set as $\delta=0.5$ s. We assume that a uniform rectangular array (URA) with 4
antennas per row is mounted on the UAV, where the adjacent antennas are
equally separated both horizontally and vertically. The number of antennas at
the UAV is $M=4$, $12$ or $20$.

For comparison, two benchmark schemes are considered: 1) MRC receive
beamforming scheme, where at most one SN can be scheduled for transmission in
each time slot, i.e., $K_{n}=1$ and $\kappa_{n}=M$, $\forall n$; 2) single-
antenna scheme, where only one antenna is mounted on the UAV, i.e., $M=1$.
Specifically, for the single-antenna scheme, only one SN can transmit in each
time slot and the achievable rate can be approximated as (14) with
$\kappa_{n}=1$, $\forall n$.

Figure 2: Max-min rate approximations versus simulation results.

(a) MRC scheme with $M=12$.

(b) Proposed scheme with $M=12$.

(c) Proposed scheme with $M=20$.

Figure 3:  Hovering locations and trajectories of different schemes under
$T=100$ s.

(a) MRC scheme with $M=12$.

(b) Proposed scheme with $M=12$.

(c) Proposed scheme with $M=20$.

Figure 4:  Hovering locations and trajectories of different schemes under
$T=200$ s.

First, in order to evaluate the accuracy of the proposed rate approximation
given in (14), Fig. 2 shows the achievable max-min rate for all SNs computed
via numerical simulations compared with the closed-form rate approximation in
(14) with different number of antennas, where the rate is averaged over
$10^{3}$ random channel realizations at each UAV location. It is observed that
the approximation (14) achieves good accuracy, especially for large number of
antennas, which is also in accordance with [17]. In particular, there exists a
small gap between the numerical simulation and the approximation for small
number of antennas, e.g., $M=4$, while such gap can be practically ignored for
sufficiently large number of antennas, e.g., $M=12,20$.

Next, we evaluate the performance of the proposed design for the minimum-rate
maximization problem. Besides the optimal hovering locations obtained by
solving (P2) without considering the UAV speed constraint, in Figs. 3 and 4,
the converged UAV trajectories obtained by solving (P1) with the speed and
initial/final location constraints are also shown for different schemes under
$T=100$ s and 200 s, respectively. First, it is observed that for the MRC
scheme in the case of $T=100$ s as shown in Fig. 3, there are 8 UAV hovering
locations, each above a different SN for receiving its data. When the number
of antennas for the MRC scheme is $M=4$ or $M=20$, the optimal hovering
locations of the UAV are the same as that for $M=12$ and thus are not shown
for brevity. The max-min rates achieved by the MRC scheme after solving (P2)
with $M=4$, 12, and 20 are 1.04 bps/Hz, 1.27 bps/Hz, and 1.37 bps/Hz,
respectively; thus, the rate gain with the increasing number of antennas is
only marginal. The reason is that under the MRC scheme, only the beamforming
gain is achieved and the UAV needs to receive data for each SN over orthogonal
time slots, thus resulting in low spectrum efficiency. As for the single-
antenna scheme, it is obvious that at the optimal solution to (P2), the UAV
should hover above each SN with the same duration, similar to the case of the
MRC scheme. Thus, the optimal hovering locations and converged trajectories of
this scheme are omitted for brevity. The max-min rate achieved by the single-
antenna scheme after solving (P2) is 0.74 bps/Hz, which is lower than those
obtained by the MRC scheme, as expected.

The optimal hovering locations obtained by the proposed scheme for solving
(P2) are shown in Figs. 3 and 3, respectively, for the case of $T=100$ s. It
is observed that the number of UAV hovering locations is less than that in the
two benchmark schemes, and it decreases with the increasing number of
antennas, which is 5 for $M=12$ and 4 for $M=20$. Specifically, 3 SNs in the
upper-left corner or 2 SNs in the lower-left corner need to be served by the
UAV at one single hovering location under the proposed scheme, rather than one
hovering location for each SN under the MRC/single-antenna scheme. Moreover,
for the proposed scheme, different from the case with $M=12$ where one
hovering location is needed to cover the SN around the center of the area of
interest, it is no longer needed when the number of antennas increases to
$M=20$. The max-min rate achieved with $M=12$ and $M=20$ under the proposed
scheme after solving (P2) is 2.57 bps/Hz and 3.41 bps/Hz, respectively, which
are significantly higher than those of the MRC scheme or the single-antenna
scheme. The reason is that different from the two benchmark schemes, where the
UAV only serves one SN in each time slot, the UAV is able to serve multiple
SNs simultaneously under the proposed scheme by exploiting the spatial
multiplexing gain with ZF receive beamforming; as a result, fewer hovering
locations are needed for the UAV. Accordingly, the flying distance of the UAV
can also be reduced after solving (P1) with the speed constraint, such that
the UAV has more time to stay at the hovering locations to achieve higher
rate, as shown next. Note that similar results can be observed for the case of
$T=200$ s in Fig. 4.

Figure 5: Max-min rate comparison between different schemes.

Figure 6: Minimum flight time required versus common throughput
requirement.

In addition, the initial and converged trajectories of the UAV for solving
(P1) with the speed constraint under $T=100$ s and $T=200$ s are also shown in
Figs. 3 and 4, respectively. The minimum time required for the UAV to visit
all hovering locations for the MRC scheme and that for the proposed scheme
with $M=12$ and $M=20$ is 181 s, 154 s, and 144 s, respectively. Hence, in the
initial trajectory design, when the given time $T=100$ s is insufficient for
the UAV to visit all hovering locations, the UAV tries to fly as close as
possible to each hovering location to shorten the link distance with the SNs,
as shown in Fig. 3. While when the given time increases to $T=200$ s, the UAV
can successively visit all hovering locations, as shown in Fig. 4. Besides,
with shorter flying distance and thus time by the proposed scheme as shown in
Figs. 4 and 4, the UAV has more time to stay at the hovering locations to
achieve higher max-min rate, which is 2.48 bps/Hz for $M=12$ and 3.32 bps/Hz
for $M=20$, as compared to 1.09 bps/Hz for the MRC scheme shown in Fig. 4.

The max-min rate achieved by solving (P1) under different $T$ by the
considered schemes is compared in Fig. 6, where that achieved by solving (P2)
without the speed constraint is also shown as the rate upper bound. It is
observed that for all the three schemes considered, the max-min rate gradually
converges to the rate upper bound with the increase of $T$. This is expected
since with a longer flight duration, the UAV has more time to stay at the
hovering locations obtained in (P2) and thus the rate achieved when it is
flying becomes more negligible. Furthermore, the achievable rate of the
proposed scheme is significantly improved with the increasing number of
antennas at the UAV and also greatly outperforms those of the benchmark
schemes.

Fig. 6 compares the minimum flight time required for the UAV by different
schemes versus the common (minimum) throughput requirement of all SNs. Since
the proposed scheme with $M=20$ significantly outperforms other schemes as
shown in Fig. 6, we only consider the proposed scheme with $M=12$, MRC scheme
with $M=12$, and the single-antenna scheme for ease of comparison. The minimum
time required for different throughput requirements can be found by solving
(P1) with an additional bisection search. It is observed in Fig. 6 that the
proposed scheme always needs the least flight time to meet the throughput
requirement of the SNs, compared to the two benchmark schemes. On the other
hand, with less time to meet the throughput requirement, the energy
consumption of the UAV can also be saved, as shown next.

Finally, to show the benefits brought by multi-antenna UAV data harvesting in
terms of energy saving, we consider the trade-off between the average transmit
power limit of SNs and the energy consumption of the UAV. In general, the
energy consumption of the UAV consists of two parts for our considered
problem. The first part is the propulsion energy, while the second part is the
communication-related energy, which is much smaller than the former for
practical UAVs and thus is ignored for simplicity. Specifically, based on [7],
the propulsion power of the UAV can be modelled as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!P_{\mathrm{h}}[n]=P_{0}\left(1+\frac{3v_{\mathrm{h}}[n]^{2}}{U_{\mathrm{tip}}^{2}}\right)+P_{i}\left(\sqrt{1+\frac{v_{\mathrm{h}}[n]^{4}}{4v_{0}^{4}}}-\frac{v_{\mathrm{h}}[n]^{2}}{2v_{0}^{2}}\right)^{1/2}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$
$\displaystyle+\frac{1}{2}d_{1}\varrho sAv_{\mathrm{h}}[n]^{3},$ (40)

where $v_{\mathrm{h}}[n]=||(\textbf{q}[n+1]-\textbf{q}[n])||/\delta$ is the
horizontal speed of the UAV, $P_{0}=79.8563$ and $P_{i}=88.6279$ are two
constants, $U_{\mathrm{tip}}=120$ m/s represents the tip speed of the rotor
blade, $v_{0}=4.03$ is the mean rotor induced velocity in hover, $d_{1}=0.6$
and $s=0.05$ are the fuselage drag ratio and rotor solidity, respectively,
$\varrho=1.225$ $\mathrm{kg}/\mathrm{m}^{3}$ and $A=0.503$ $\mathrm{m}^{2}$
denote the air density and rotor disc area, respectively. Fig. 7 shows the
trade-off between the average transmit power limit of SNs and the (propulsion)
energy consumption of the UAV for different schemes under the same throughput
requirement of 4 Mbits for all SNs. It is observed that for all the three
schemes, as the average transmit power $\bar{P}$ of the SNs increases, the
energy consumption of the UAV decreases, which demonstrates the energy trade-
off between SNs and the UAV [28], and such a trade-off is more evident for
smaller $\bar{P}$. This is expected since with higher transmit power of SNs,
the UAV generally needs less flight time to meet their throughput requirement,
thus leading to lower energy consumption. On the other hand, it is observed
that the proposed scheme always outperforms both benchmark schemes in terms of
the SNs-UAV energy consumption trade-off, as expected.

Figure 7: Trade-off between average transmit power limit of SNs and
energy consumption of UAV.

##  VI Conclusion

This paper studies a UAV-enabled WSN, where a multi-antenna UAV is employed to
collect data from a group of single-antenna SNs. Our aim is to maximize the
minimum rate of all SNs via jointly optimizing the transmission scheduling and
power allocations of SNs as well as the UAV trajectory, subject to practical
SNs’ transmit power and UAV speed constraints. To tackle this challenging non-
convex optimization problem, we first consider its relaxed problem by ignoring
the UAV maximum speed constraint and derive the optimal solution to it. Based
on this solution, we further propose a suboptimal solution to the general
problem with the UAV speed and initial/final location constraints. Numerical
results are provided to show significant performance gains by the proposed
scheme, as compared to benchmark schemes, in terms of both spectrum efficiency
and SNs-UAV energy consumption trade-off.

## Appendix

## Optimal Solution to (P3)

In this appendix, we derive the optimal solution to problem (P3). Since the
problem can be optimally solved by the Lagrange duality method similar to
(P2), we first express the partial Lagrangian of (P3) as

$\displaystyle\mathcal{L}_{2}(\\{p_{k}[n]\\},\\{a_{k}[n]\\},r,\bm{\nu},\bm{\phi})=~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\left(1-\sum_{k=1}^{K}\nu_{k}\right)r+\sum_{k=1}^{K}\nu_{k}\bar{r}_{k}+\sum_{k=1}^{K}\phi_{k}\left(N\bar{P}-\sum_{n=1}^{N}p_{k}[n]\right),$
(41)

where $\bm{\nu}\triangleq\\{\nu_{k}\\}$ and
$\bm{\phi}\triangleq\\{\phi_{k}\\}$ are the dual variables associated with the
constraints (35a) and (35d), respectively. The Lagrange dual function of (P3)
is thus given by

$\displaystyle\\!\\!\\!\\!\\!\\!g_{2}(\bm{\nu},\bm{\phi})=~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\\!\\!\\!\left\\{\begin{aligned}
\max_{\begin{subarray}{c}\\{p_{k}[n]\\},\\{a_{k}[n]\\},r\end{subarray}}\mathcal{L}_{2}\left(\\{p_{k}[n]\\},\\{a_{k}[n]\\},r,\bm{\nu},\bm{\phi}\right)~{}~{}~{}~{}~{}~{}\\\
\mathrm{s.t.}~{}~{}~{}~{}~{}~{}~{}~{}\eqref{p3002},~{}\eqref{p3003},~{}\eqref{p3005},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\end{aligned}\right.\\!\\!\\!\\!\\!\\!\\!\\!$
(42)

For (Optimal Solution to (P3)) to be bounded, we should have
$1-\sum_{k=1}^{K}\nu_{k}=0$. Therefore, the dual problem of (P3) is given by

$\displaystyle\mathrm{(D3)}~{}~{}\min_{\bm{\nu},\bm{\phi}}~{}~{}~{}g_{2}(\bm{\nu},\bm{\phi})~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\mathrm{s.t.}~{}~{}~{}1-\sum_{k=1}^{K}\nu_{k}=0,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
(43a) $\displaystyle\nu_{k}\geq 0,\ \ \forall
k\in\mathcal{K},~{}~{}~{}\\!~{}~{}~{}~{}~{}~{}$ (43b)
$\displaystyle\phi_{k}\geq 0,\ \ \forall
k\in\mathcal{K}.~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\!\\!$ (43c)

Let the feasible set of $\bm{\nu}$ and $\bm{\phi}$ characterized by the
constraints in (43a)-(43c) as $\chi_{2}$. In the following, we solve (P3)
based on the similar procedures as for (P2).

###  VI-A Obtaining $g_{2}($$\bm{\nu}$$,~{}\\!\\!$$\bm{\phi}$$)$ by solving
problem (Optimal Solution to (P3)) for given
$($$\bm{\nu}$$,~{}\\!\\!$$\bm{\phi}$$)$ $\in\mathcal{X}^{2}$

Similar to (III-A) and (III-A), the optimal transmit power and achievable rate
over time slots are respectively given by

$\displaystyle
p_{k}^{*}[n]=a_{k}[n]\left[\frac{\nu_{k}}{N\phi_{k}\log(2)}-\frac{d_{k}[n]^{\alpha}}{\kappa_{n}\gamma_{0}}\right]^{+},~{}~{}~{}~{}~{}$
(44) $\displaystyle
r_{k}^{*}[n]=a_{k}[n]\left[\log_{2}\left(\frac{\nu_{k}\kappa_{n}\gamma_{0}}{N\phi_{k}\log(2)d_{k}[n]^{\alpha}}\right)\right]^{+}.$
(45)

Then, in (Optimal Solution to (P3)) we have

$\displaystyle\mathcal{L}_{2}\left(\\{p_{k}[n]\\},\\{a_{k}[n]\\},r,\bm{\nu},\bm{\phi}\right)~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!=\sum_{k=1}^{K}\nu_{k}\left(\frac{1}{N}\sum_{n=1}^{N}r_{k}^{*}[n]\right)+\sum_{k=1}^{K}\phi_{k}\left(N\bar{P}-\sum_{n=1}^{N}p_{k}[n]^{*}\right)\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$
$\displaystyle=\sum_{n=1}^{N}\left(\sum_{k=1}^{K}a_{k}[n]\varphi_{k}[n]\right)+N\bar{P}\sum_{k=1}^{K}\phi_{k},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
(46)

where

$\displaystyle\varphi_{k}[n]\triangleq\frac{\nu_{k}}{N}\left[\log_{2}\left(\frac{\nu_{k}\kappa_{n}\gamma_{0}}{N\phi_{k}\log(2)d_{k}[n]^{\alpha}}\right)\right]^{+}$
$\displaystyle-\phi_{k}\left[\frac{\nu_{k}}{N\phi_{k}\log(2)}-\frac{d_{k}[n]^{\alpha}}{\kappa_{n}\gamma_{0}}\right]^{+}.$
(47)

Thus, problem (Optimal Solution to (P3)) is recast to

$\displaystyle\max_{\\{a_{k}[n]\\},\kappa_{n}}~{}~{}\sum_{n=1}^{N}\left(\sum_{k=1}^{K}a_{k}[n]\varphi_{k}[n]\right)+N\bar{P}\sum_{k=1}^{K}\phi_{k}$
(48)
$\displaystyle\mathrm{s.t.}~{}~{}~{}~{}\\!~{}~{}~{}\eqref{p3002},~{}\eqref{p3003}.~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$

Note that problem (48) can be decoupled into $N$ sub-problems, each
corresponding to one time slot. For each sub-problem, the optimal scheduling
can be easily found via an exhaustive search over all possible values of
$\kappa_{n}$, similar to (27). However, different from problem (26) where each
sub-problem is the same across different time slots, each sub-problem in
problem (48) is different and thus needs to be solved separately.

###  VI-B Finding optimal $\bm{\nu}$ and $\bm{\phi}$ to solve (D3)

Similar to Section III-B, we can use the ellipsoid method to obtain optimal
$\bm{\nu}$ and $\bm{\phi}$. The subgradients with respect to
$(\bm{\nu},\bm{\phi})$ are given by

$\displaystyle\\!\Delta\nu_{k}=\frac{1}{N}\sum_{n=1}^{N}r_{k}[n]^{*},~{}~{}~{}~{}\\!~{}$
(49) $\displaystyle\Delta\phi_{k}=N\bar{P}-\sum_{n=1}^{N}p_{k}[n]^{*},$ (50)

where $r_{k}[n]^{*}$ and $p_{k}[n]^{*}$ are given in (45) and (44),
respectively.

Finally, since the obtained $\\{a_{k}[n]\\}$ for each sub-problem in (48) is
unique over different time slots, there is no need to further construct the
optimal solution for (P3) by time-sharing as for (P2). The rest of the
algorithm for solving (P3) is similar to Algorithm 1 for (P2), and thus is
omitted for brevity.

## References

  * [1] M. Tubaishat and S. Madria, “Sensor networks: An overview,” _IEEE Potentials_ , vol. 22, no. 2, pp. 20–23, Apr. 2003.
  * [2] P. Rawat, K. D. Singh, H. Chaouchi, and J. M. Bonnin, “Wireless sensor networks: A survey on recent developments and potential synergies,” _J.. Supercomput._ , vol. 68, no. 1, pp. 1–48, 2014.
  * [3] F. Wang and J. Liu, “Networked wireless sensor data collection: Issues, challenges, and approaches,” _IEEE Commun. Surveys Tuts._ , vol. 13, no. 4, pp. 673–687, Fourth Quarter 2011.
  * [4] Y. Zeng, Q. Wu, and R. Zhang, “Accessing from the sky: A tutorial on UAV communications for 5G and beyond,” _Proc. IEEE_ , vol. 107, no. 12, pp. 2327–2375, Dec. 2019.
  * [5] Y. Zeng, R. Zhang, and T. J. Lim, “Throughput maximization for UAV-enabled mobile relaying systems,” _IEEE Trans. Commun._ , vol. 64, no. 12, pp. 4983–4996, Dec. 2016.
  * [6] Y. Zeng and R. Zhang, “Energy-efficient UAV communication with trajectory optimization,” _IEEE Trans. Wireless Commun._ , vol. 16, no. 6, pp. 3747–3760, Jun. 2017.
  * [7] Y. Zeng, J. Xu, and R. Zhang, “Energy minimization for wireless communication with rotary-wing UAV,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 4, pp. 2329–2345, Apr. 2019.
  * [8] C. Zhan, Y. Zeng, and R. Zhang, “Energy-efficient data collection in UAV enabled wireless sensor network,” _IEEE Wireless Commun. Lett._ , vol. 7, no. 3, pp. 328–331, Jun. 2018.
  * [9] C. Choi, F. Baccelli, and G. de Veciana, “Analysis of data harvesting by unmanned aerial vehicles,” _arXiv preprint arXiv:1902.06350_ , 2019.
  * [10] C. You and R. Zhang, “3D trajectory optimization in Rician fading for UAV-enabled data harvesting,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 6, pp. 3192–3207, Jun. 2019.
  * [11] D. Ebrahimi, S. Sharafeddine, P. Ho, and C. Assi, “UAV-aided projection-based compressive data gathering in wireless sensor networks,” _IEEE Internet Things J._ , vol. 6, no. 2, pp. 1893–1905, Apr. 2019.
  * [12] J. Zhang, Y. Zeng, and R. Zhang, “Receding horizon optimization for energy-efficient UAV communication,” _IEEE Wireless Commun. Lett._ , Early Access, 2019.
  * [13] J. Chen, B. Daneshrad, and Weijun Zhu, “MIMO performance evaluation for airborne wireless communication systems,” in _Mil. Commun. Conf. (MILCOM)_ , Nov. 2011, pp. 1827–1832.
  * [14] F. Jiang and A. L. Swindlehurst, “Optimization of UAV heading for the ground-to-air uplink,” _IEEE J. Sel. Areas Commun._ , vol. 30, no. 5, pp. 993–1005, Jun. 2012.
  * [15] M. Mozaffari, W. Saad, M. Bennis, and M. Debbah, “Communications and control for wireless drone-based antenna array,” _IEEE Trans. Commun._ , vol. 67, no. 1, pp. 820–834, Jan. 2019.
  * [16] W. Feng, J. Wang, Y. Chen, X. Wang, N. Ge, and J. Lu, “UAV-aided MIMO communications for 5G Internet of things,” _IEEE Internet Things J._ , vol. 6, no. 2, pp. 1731–1740, Apr. 2019.
  * [17] L. Liu, S. Zhang, and R. Zhang, “CoMP in the sky: UAV placement and movement optimization for multi-user communications,” _IEEE Trans. Commun._ , vol. 67, no. 8, pp. 5645–5658, Aug. 2019.
  * [18] W. Yu and R. Lui, “Dual methods for nonconvex spectrum optimization of multicarrier systems,” _IEEE Trans. Commun._ , vol. 54, no. 7, pp. 1310–1322, Jul. 2006.
  * [19] S. Boyd and J. Duchi, “EE364b: Convex optimization II,” _Course Notes_. [Online]. Available: http://www.stanford.edu/class/ee364b
  * [20] S. Boyd and L. Vandenberghe, _Convex optimization_. Cambridge university press, 2004.
  * [21] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex programming,” 2014.
  * [22] Q. Wu, J. Xu, and R. Zhang, “Capacity characterization of UAV-enabled two-user broadcast channel,” _IEEE J. Sel. Areas Commun._ , vol. 36, no. 9, pp. 1955–1971, Sep. 2018.
  * [23] J. Xu, Y. Zeng, and R. Zhang, “UAV-enabled wireless power transfer: Trajectory design and energy optimization,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 8, pp. 5092–5106, Aug. 2018.
  * [24] L. Xie, J. Xu, and R. Zhang, “Throughput maximization for UAV-enabled wireless powered communication networks,” _IEEE Internet Things J._ , vol. 6, no. 2, pp. 1690–1703, Apr. 2018.
  * [25] J. Zhang, Y. Zeng, and R. Zhang, “UAV-enabled radio access network: Multi-mode communication and trajectory design,” _IEEE Trans. Signal Process._ , vol. 66, no. 20, pp. 5269–5284, Oct. 2018.
  * [26] Y. Zeng, X. Xu, and R. Zhang, “Trajectory design for completion time minimization in UAV-enabled multicasting,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 4, pp. 2233–2246, Apr. 2018.
  * [27] Iskandar and S. Shimamoto, “The channel characterization and performance evaluation of mobile communication employing stratospheric platform,” in _IEEE/ACES Int. Conf. Wireless Commun. Appl. Computat. Electromagn. (ICWCACE)_ , Apr. 2005, pp. 828–831.
  * [28] D. Yang, Q. Wu, Y. Zeng, and R. Zhang, “Energy tradeoff in ground-to-UAV communication via trajectory design,” _IEEE Trans. Veh. Technol._ , vol. 67, no. 7, pp. 6721–6726, Jul. 2018.