# Energy-Efficient UAV Communications in the Presence of Wind: 3D Modeling and
Trajectory Design

Xinhong Dai, Bin Duo, , Xiaojun Yuan, , Marco Di Renzo X. Dai, B. Duo and X.
Yuan are with the National Key Laboratory on Wireless Communications, the
University of Electronic Science and Technology of China, Chengdu 611731,
China (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>xjyuan@uestc.edu.cn). The corresponding author is Xiaojun Yuan.M. Di Renzo is
with Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire des Signaux
et Systèmes, 3 Rue Joliot-Curie, 91192 Gif-sur-Yvette, France. (marco.di-
renzo@universite-paris-saclay.fr).

###### Abstract

The rapid development of unmanned aerial vehicle (UAV) technology provides
flexible communication services to terrestrial nodes. Energy efficiency is
crucial to the deployment of UAVs, especially rotary-wing UAVs whose
propulsion power is sensitive to the wind effect. In this paper, we first
derive a three-dimensional (3D) generalised propulsion energy consumption
model (GPECM) for rotary-wing UAVs under the consideration of stochastic wind
modeling and 3D force analysis. Based on the GPECM, we study a UAV-enabled
downlink communication system, where a rotary-wing UAV flies subject to
stochastic wind disturbance and provides communication services for ground
users (GUs). We aim to maximize the energy efficiency (EE) of the UAV by
jointly optimizing the 3D trajectory and user scheduling among the GUs based
on the GPECM. We formulate the problem as stochastic optimization, which is
difficult to solve due to the lack of real-time wind information. To address
this issue, we propose an offline-based online adaptive (OBOA) design with two
phases, namely, an offline phase and an online phase. In the offline phase, we
average the wind effect on the UAV by leveraging stochastic programming (SP)
based on wind statistics; then, in the online phase, we further optimize the
instantaneous velocity to adapt the real-time wind. Simulation results show
that the optimized trajectories of the UAV in both two phases can better adapt
to the wind in changing speed and direction, and achieves a higher EE compared
with the windless scheme. In particular, our proposed OBOA design can be
applied in the scenario with dramatic wind changes, and makes the UAV adjust
its velocity dynamically to achieve a better performance in terms of EE.

###### Index Terms:

Rotary-wing UAV, wind effect, generalised propulsion energy consumption model,
UAV communications, energy efficiency, 3D trajectory design, stochastic
programming.

##  I Introduction

Unmanned aerial vehicles (UAVs) have attracted considerable attention on
providing cost-effective and diversified services in future wireless
communication systems [1, 2, 3]. Thanks to their unique features, such as
versatility, high maneuverability, and swift deployment capabilities, UAVs can
easily establish line-of-sight (LoS) links with communication nodes deployed
on the ground for high quality data delivery[4], and can adapt to diversely
complicated and time-varying environments[5, 6, 7]. UAVs can be used in many
applications, e.g., data collection for Internet of Things (IoTs), mobile edge
computing, secure communication, and emergency communication, etc.[8, 9, 10,
11, 12, 13, 14].

Prior works on UAV-enabled wireless communications mostly emphasize on the
placements and trajectories optimization of UAVs to improve communication
performance [15, 16, 17, 18, 19]. However, with the diversification of
missions, energy saving becomes a more and more important issue on the design
of UAV-enabled communication systems. Specifically, due to the size and weight
constraints, limited energy on-board significantly confines their long-term
propulsion energy consumption and other prolonged operations of UAVs[20].
Prior research on UAV-enabled wireless communications employed simple energy
consumption constraints[21, 22, 23] that cannot effectively characterize the
propulsion energy consumption of a UAV in a real environment. As such,
establishing a generalised propulsion energy consumption model (GPECM) to
quantify the amount of UAV energy consumption more accurately is of critical
importance to the design of UAV-enabled communication systems.

The existing studies introduced various UAV propulsion energy consumption
models from two or three-dimensional (2D or 3D) perspectives and proposed the
corresponding energy-efficient designs for UAV-enabled communication systems.
For rotary-wing UAVs, the authors of [24] derived a propulsion energy
consumption model in 2D space by considering UAV velocity, and optimized the
trajectory to minimize the total energy consumption. In [25], the authors
further improved the model in [24] to include the velocity, acceleration and
direction change of a UAV, and proposed a smooth and energy-efficient
trajectory design. The authors of [26] provided a 3D propulsion energy
consumption model, with a linear term accounting for the vertical component of
UAV velocity. For fixed-wing UAVs, the authors of [27] proposed a 2D
propulsion energy consumption model taking into account UAV velocity and
acceleration, where the energy efficiency (EE) was maximized for a given
flight duration. The authors of [28] expanded the model in [27] by analyzing
the dynamics of the fuselage axes and derived a 3D propulsion energy
consumption model that captures the impact of UAV velocity, acceleration, and
pitch angle, and minimized the total energy consumption of a UAV subject to
secret data collection under eavesdropping attack.

In general, the impact of stochastic winds cannot be neglected during UAV
flight [29]. No matter flying or hovering, a UAV has to withstand the
atmospheric drag caused by wind[30], which causes additional propulsion energy
consumption. Thus, the wind disturbance presents a real challenge to the
design of UAV-enabled communication systems, especially for rotary-wing UAVs
that are sensitive to the wind effect. Considering the wind effect, the
authors of [31] proposed an energy-saving path design of a UAV to utilize the
wind. In [32], the authors described the wind effect on a UAV by the
horizontal angle deviation between its ground velocity and the wind direction,
and minimized the propulsion energy consumption subject to the given data
collection constraints, by optimizing the 2D trajectories of the fixed-wing
UAV.

As a key factor affecting UAV flight and propulsion power, the wind
disturbance on UAV cannot be neglected for the energy-efficient design in UAV-
enabled communication systems. However, the existing approaches either did not
consider the wind effect, or did not accurately describe the impact of the
wind on aerodynamics parameters (e.g.,drag and thrust) and propulsion power of
UAV. Besides, real-time winds are usually random, such that traditional
energy-efficient trajectory designs make limited sense in practical
applications. Therefore, it is crucial to establish a GPECM to characterize
the 3D movement of a UAV in wind more accurately. Based on this model, energy-
efficient trajectory designs can be made possible in UAV-enabled communication
systems under the stochastic wind disturbance.

Motivated by the above, we consider a UAV-enabled downlink communication
system subject to the stochastic wind in an urban area, and propose an energy-
efficient trajectory design by jointly optimizing the 3D UAV trajectory and
the user scheduling. To the best of our knowledge, this is the first attempt
to the 3D energy-efficient trajectory design of the UAV-enabled communication
system that considers a 3D GPECM of UAV in wind. The main contributions of
this paper are summarized as follows:

  * •

We propose a 3D GPECM for a rotary-wing UAV in wind based on the wind modeling
and force analysis of the UAV. This model characterizes the real-time wind
effect on the propulsion power and 3D movement including UAV velocity and
acceleration.

  * •

With the proposed GPECM, we consider a UAV-enabled urban downlink
communication system. To maximize the EE of the UAV, we formulate a stochastic
optimization problem by jointly optimizing 3D trajectory and user scheduling
of the UAV. This problem is challenging to solve due to the difficulty in the
acquirement of real-time wind information. To address this difficulty, we
propose a two-phase offline-based online adaptive (OBOA) design to achieve an
energy-efficient design for the above system.

  * •

In the offline phase, we assume that the statistical distribution of the wind
speed and direction is known in prior. Then, we propose an offline design to
maximize the EE of the UAV, by alternately optimizing the UAV’s horizontal
trajectory, vertical trajectory and user scheduling based on the stochastic
programming (SP) technique. This offline solution is used as a reference to
the optimization in the online phase.

  * •

In the online phase, to make the UAV adapt better to the changeable real-time
wind and reduce the impact of the wind on its propulsion energy consumption,
an OBOA design is proposed to minimize the propulsion energy consumption in
each time slot, by optimizing UAV instantaneous velocity slot by slot with
offline trajectory taken as a reference. The OBOA design for each time slot
can be implemented with relatively low complexity.

Extensive simulations are conducted to validate the effectiveness of our
proposed design. Specifically, simulation results show that the the optimized
trajectories of both offline and online phases can adapt better to the wind in
changing speed and direction, and can achieve a higher EE compared to the
trajectory design in a windless environment. Furthermore, with online adaptive
optimization, the proposed OBOA design makes the UAV adapt better to the
changeable real-time wind by dynamically adjusting its instantaneous velocity,
yielding a reduced angular deviation between the UAV velocity and the real-
time wind. Numerical results show that even in the scenario with dramatic wind
changes, our proposed OBOA design shows a better performance in terms of EE.

The rest of this paper has the following structure: In Section ii@, the 3D
GPECM in wind of the UAV is derived. In Section iii@, we describe the
communication system model. In Section iv@, we formulate a problem of
maximizing the EE of the UAV and propose the OBOA design. The offline and
online phases are described in v@ and vi@, respectively. We discuss the
simulation results in Section vii@ and conclude this paper in Section<EMAIL_ADDRESS>
##  II 3D Generalised Propulsion Energy Consumption Model in Wind

In this section, we establish a GPECM of the UAV in wind in the 3D space. The
GPECM is formulated as a function of velocity $\mathbf{v}$, acceleration
$\mathbf{a}$, thrust $\mathbf{T}$ of the UAV, and wind velocity
$\mathbf{v}_{w}$, defined as the vectors in the 3D space.

###  II-A Wind Modeling and Force Analysis

We model the wind as a random vector, and establish the relationship between
thrust $\mathbf{T}$ and other parameters by the 3D force analysis on the UAV
in real-time wind. Since the wind speed in vertical direction is commonly
small [33], we only consider the wind in the horizontal direction for
simplicity, and the wind vector can be expressed as $\mathbf{v}_{w}\triangleq\
\|\mathbf{v}_{w}\|\angle\beta$, where $\|\mathbf{v}_{w}\|$ and $\beta$ denote
the speed and direction angle of the wind, respectively.

We model $\|\mathbf{v}_{w}\|$ as a random variable drawn from a Weibull
distribution[34], by considering the wind speed changes with respect to
(w.r.t.) altitude. Specifically, we denote $h_{\rm ref}$ and $v_{\rm ref}$ as
the reference altitude and the corresponding average wind speed. Then, from
[35] and [36], the wind speed at any altitude can be expressed as
$\|\mathbf{v}_{w}\|=v_{\rm ref}\left(\frac{z}{h_{\rm ref}}\right)^{p}$, where
$z$ is the measurement height of the wind, $p$ is an empirical exponent that
depends on atmospheric conditions and $p\in[0.4,0.6]$ in the city area as
shown in the following system model[35], and $v_{\rm ref}$ follows the Weibull
distribution with the probability density function (PDF)

$\displaystyle f(v_{\rm ref}|\lambda,c)=\frac{c}{\lambda}\left(\frac{v_{\rm
ref}}{\lambda}\right)^{c-1}e^{-\left(\frac{v_{\rm
ref}}{\lambda}\right)^{c}},\quad\forall v_{\rm ref}>0,$ (1)

where $\lambda$ is the scaling parameter and $c$ is the shape parameter111From
[34], the mean of the Weibull distribution is dominated by $\lambda$, and the
variance is dominated by $c$. Therefore, to simplify the analysis, we assume
that the increase of $\lambda$ represents the increase of the average wind
speed, and the increase of $c$ represents the decrease of the variance of the
wind speed change.. It is worth noting that, based on our modeling, the wind
effect on the UAV in the 3D space is closely related to the UAV flight height,
and becomes stronger with the increasing flight height. For the wind
direction, we model the direction angle $\beta$ as a random variable with a
Von-Mises distribution [37], and the PDF of $\beta$ can be written as

$\displaystyle f(\beta|\mu,\kappa)=\frac{e^{\kappa\cos{(\beta-\mu)}}}{2\pi\rm
I_{0}(\kappa)},$ (2)

where $\mu$ is the expectation angle, $\kappa$ is the concentration parameter
and $\rm I_{0}(\cdot)$ denotes the zero order Bessel function222From [37] the
concentration parameter $\kappa$ denotes the variance of the Von-Mise
distribution, i.e., the larger $\kappa$ means a more concentrated angle
distribution with a smaller variance of the wind direction change.. The wind
at the height $z$ can be expressed as a vector:

$\displaystyle\mathbf{v}_{w}(v_{\rm ref},\beta,z)=\left(\frac{z}{h_{\rm
ref}}\right)^{p}(v_{\rm ref}\cos{\beta},v_{\rm ref}\sin{\beta}).$ (3)

Figure 1: Instantaneous 3D force analysis on a rotary-wing UAV in real-time
wind

With the modeling of the wind, we further take a 3D force analysis of the UAV,
as shown in Fig. 1. Specifically, we denote the weight of the UAV and
gravitational acceleration as $m$ and $\mathbf{g}$, respectively. The UAV is
mainly affected by gravity $m\mathbf{g}$, thrust $\mathbf{T}$ from the rotors,
air drag $\mathbf{D}_{a}$ due to UAV flight and the drag $\mathbf{D}_{w}$
caused by the wind. To simplify the analysis, we neglect the viscous effect on
UAV movement of the air. Thus, $\mathbf{D}_{a}$ and $\mathbf{D}_{w}$ can be
handled as two independent forces. Without loss of generality, we assume the
drag caused by fluid of the UAV is isotropic. Then, the relative velocity of
the UAV to the air can be denoted as $\mathbf{v}_{w}-\mathbf{v}$. Based on the
definition of the drag in fluid [38], the resultant drag caused by the UAV
flight and wind can be expressed as

$\displaystyle\|\mathbf{D}\|=\|\mathbf{D}_{a}+\mathbf{D}_{w}\|=\frac{1}{2}\rho
S_{\rm FP}\|\mathbf{v}-\mathbf{v}_{w}\|^{2},$ (4)

where $\rho$ and $S_{\rm FP}$ represents the air density and fuselage
equivalent flat area of the UAV, respectively. Furthermore, from Newton’s
second law, we have $m\mathbf{a}=\mathbf{T}+\mathbf{D}+m\mathbf{g}$.
Therefore, the thrust of the UAV can be written as

$\displaystyle\|\mathbf{T}\|=\|m\mathbf{a}+\frac{1}{2}\rho S_{\rm
FP}\|\mathbf{v}-\mathbf{v}_{w}\|(\mathbf{v}-\mathbf{v}_{w})-m\mathbf{g}\|.$
(5)

###  II-B 3D GPECM Derivation

With the above wind modeling and force analysis, we derive the GPECM of the
UAV. Specifically, from [24] and [38], the propulsion power of a rotary-wing
UAV can be expressed as

$\displaystyle P=q_{c}\rho sAv_{\rm tip}^{3},$ (6)

where the rotor solidity $s$ and disc area $A$ are both constant parameters,
and $v_{\rm tip}$ denotes the tip speed of a rotor blade. From [38,(4.20)],
the torque coefficient $q_{c}$ can be written as:

$\displaystyle
q_{c}=\frac{\delta}{8}(1+3\hat{v}^{2})+(1+c_{f})\lambda_{i}t_{c}+c_{w}\hat{v}\sin{\tau_{c}}+\frac{1}{2}d_{0}\hat{v}^{3}.$
(7)

In addition, we have the following equations from [38]: 1) Forward speed
normalized on the tip speed $\hat{v}=\frac{\|\mathbf{v}\|}{v_{\rm tip}}$; 2)
thrust coefficient based on the disc axes
$t_{c}\approx\frac{\|\mathbf{T}\|}{\rho sAv_{\rm tip}^{2}}$; 3) mean induced
velocity
$\lambda_{i}=(\sqrt{v_{o}^{4}+\frac{\|\mathbf{v}\|^{4}}{4}}-\frac{\|\mathbf{v}\|^{2}}{2})^{\frac{1}{2}}$,
where $v_{o}=\sqrt{\frac{\|\mathbf{T}\|}{2\rho A}}$ is the thrust velocity; 4)
weight coefficient $c_{w}=\frac{\|m\mathbf{g}\|}{\rho sAv_{\rm tip}^{2}}$; 5)
fuselage drag ratio $d_{0}=\frac{S_{\rm FP}}{sA}$. According to the thrust
coefficient based on the disc area $c_{T}$, we have $\|\mathbf{T}\|=c_{T}\rho
Av_{\rm tip}^{2}$, and other parameters are detailed in [24, Table i@] and
[38]. With (4) and (5), we substitute the above equations into (6), yielding
the UAV’s 3D GPECM in wind as

$\displaystyle
P=P_{b}+P_{i}+m\|\mathbf{g}\|\|\mathbf{v}\|\sin{\tau_{c}}+\frac{1}{2}\rho
S_{\rm FP}\|\mathbf{v}-\mathbf{v}_{w}\|^{3},$ (8)

where $P_{b}=\frac{\delta}{8}(\frac{\|\mathbf{T}\|}{c_{T}\rho
A}+3\|\mathbf{v}\|^{2})\sqrt{\frac{\rho s^{2}A\|\mathbf{T}\|}{c_{T}}}$ and
$P_{i}=(1+c_{f})\|\mathbf{T}\|\left(\sqrt{\frac{\|\mathbf{T}\|^{2}}{(2\rho
A)^{2}}+\frac{\|\mathbf{v}\|^{4}}{4}}-\frac{\|\mathbf{v}\|^{2}}{2}\right)^{\frac{1}{2}}$
represent the blade profile power and the induced power, respectively. The
last two terms in (8) represent the climbing power and the power for
overcoming the air drag, respectively. This GPECM characterizes the practical
flight status of a rotary-wing UAV by the 3D aerodynamic analysis under the
consideration of the wind, which helps to provide a more accurate energy-
efficient design in the following UAV-enabled communication system. Note that
the derived GPECM (8) is a more general form of existing propulsion energy
consumption models for rotary-wing UAVs e.g.,[24, 25, 26], and it can be
reduced to these models without considering the wind effect under different
assumptions. The simplifications from the GPECM to the existing models in [24,
25, 26] are detailed in Appendix A.

##  III System Model for UAV-Enabled Communications

Figure 2: Downlink communication system aided by a rotary-wing UAV in wind

As shown in Fig. 2, we consider a general UAV-enabled downlink communication
system, where a rotary-wing UAV is dispatched to fly in an urban area and
provide the communication services to $K$ ground users (GUs) indexed by the
set $\mathcal{K}=\\{1,\dots,K\\}$ in the 3D space. The UAV has to withstand
strong wind with a given flight duration $T_{0}$. Assume that the location of
each GU _k_ is represented by $(\mathbf{g}_{k}^{\mathrm{T}},0)$,
$k\in\mathcal{K}$, where
$\mathbf{g}_{k}=[x_{k},y_{k}]^{\mathrm{T}}\in\mathbb{R}^{2\times 1}$.

###  III-A Discretized UAV Mobility Model

For simplicity, we discretize $T_{0}$ into $N$ equal time slots, indexed by
the $\mathcal{N}=\\{1,\dots,N\\}$. Let $\Delta=T_{0}/N$ denote the each time
slot, which is set sufficiently small. Denote the UAV trajectory by
$\mathbf{Q}=(\mathbf{q}^{\mathrm{T}},z)$, where
$\mathbf{q}=(x,y)^{\mathrm{T}}$ and $z$ represent the horizontal location and
the flight altitude, respectively. As such, the UAV trajectory can be
approximated by the $N$-length 3D position sequence
$\\{(\mathbf{q}^{\mathrm{T}}[n],z[n])\\}_{n=1}^{N}$, where the choice of $N$
largely determines the complexity of 3D UAV trajectory optimization. Then, UAV
velocity and acceleration can be defined as:

$\displaystyle\mathbf{v}[n]\triangleq\frac{\mathbf{Q}[n+1]-\mathbf{Q}[n]}{\Delta},\quad\mathbf{a}[n]\triangleq\frac{\mathbf{Q}[n+2]+\mathbf{Q}[n]-2\mathbf{Q}[n+1]}{\Delta^{2}},\quad\forall
n.$ (9)

Within each time slot $n$, we assume that the UAV location remains unchanged,
and the distance between the UAV and GU $k$ is given by
$d_{k}[n]=\sqrt{z[n]^{2}+\|\mathbf{q}[n]-\mathbf{g}_{k}\|^{2}}$, where
$\|\cdot\|$ denotes the Euclidean norm. Moreover, $\mathbf{Q}_{\rm
I}=(\mathbf{q}^{\mathrm{T}}[1],z[1])$ and $\mathbf{Q}_{\rm
F}=(\mathbf{q}^{\mathrm{T}}[N],z[N])$ denote the initial and final locations
of the UAV, respectively, and the UAV is launched and landed at certain
locations. We assume that the horizontal and vertical speeds of the UAV are
respectively limited by $V_{H}^{\rm max}$ and $V_{V}^{\rm max}$. The maximum
UAV flight distance in the horizontal and vertical directions within each time
slot are denoted by $D_{H}^{\rm max}=V_{H}^{\rm max}\Delta$ and $D_{V}^{\rm
max}=V_{V}^{\rm max}\Delta$, respectively. Then, the UAV trajectory satisfies
the following constraints:

$\displaystyle\|\mathbf{q}[n+1]-\mathbf{q}[n]\|\leq D_{H}^{\rm
max},\quad|z[n+1]-z[n]|\leq D_{V}^{\rm max},\quad\forall n.$ (10)

The flight height of the UAV is subject to the following constraints333In
general, most of small drones (e.g., rotary-wing UAVs) must avoid collisions
with obstacles when flying over cities, due to the highly random and
changeable buildings in the city. In addition, the UAV must fly at an altitude
within a legitimate range to not violate relevant aerial regulations.

$\displaystyle H_{\rm min}\leq z[n]\leq H_{\rm max},\quad\forall n,$ (11)

where $H_{\rm min}$ and $H_{\rm max}$ are the maximum and minimum altitudes of
the UAV. In addition, based on the above discretized model, the GPECM in (8)
can be rewritten as

$\displaystyle
P[n]=P_{b}[n]+P_{i}[n]+m\|\mathbf{g}\|\|\mathbf{v}[n]\|\sin{\tau_{c}[n]}+\frac{1}{2}\rho
S_{\rm FP}\|\mathbf{v}[n]-\mathbf{v}_{w}[n]\|^{3},$ (12)

where $\|\mathbf{T}[n]\|=\|m\mathbf{a}[n]+\frac{1}{2}\rho S_{\rm
FP}(\mathbf{v}[n]-\mathbf{v}_{w}[n])-m\mathbf{g}\|$,
$P_{b}[n]=\frac{\delta}{8}(\frac{\|\mathbf{T}[n]\|}{c_{T}\rho
A}+3\|\mathbf{v}[n]\|^{2})\sqrt{\frac{\rho s^{2}A\|\mathbf{T}[n]\|}{c_{T}}}$,
and
$P_{i}[n]=(1+c_{f})\|\mathbf{T}[n]\|\left(\sqrt{\frac{\|\mathbf{T}[n]\|^{2}}{(2\rho
A)^{2}}+\frac{\|\mathbf{v}[n]\|^{4}}{4}}-\frac{\|\mathbf{v}[n]\|^{2}}{2}\right)^{\frac{1}{2}}$.

###  III-B UAV-GUs Communication Model

We assume that the UAV and all GUs are each equipped with a single isotropic
antenna of unit gain to simplify the analysis. When the UAV communicates with
the GUs in an urban area, the communication links can be occasionally blocked
by buildings. Thus, the LoS link based on the free-space path loss model
cannot be directly applied in our proposed system. To characterize the UAV-GUs
communication channel model practically, we consider an improved probabilistic
LoS channel model in [39], where the typical parameters of the model are
obtained from a Manhattan-type city model444Note that the channel model is
divided into two categories, i.e., LoS and NLoS, and the UAV-GUS channel
depends on the real environment in practice. We consider this probabilistic
LoS channel model to obtain the maximum EE of the UAV from a statistical
perspective in the offline phase, and consider the practical channel model in
the online phase, which are detailed in the following sections.. Specifically,
the LoS probability between the UAV and GU $k$ in time slot $n$, denoted by
$P_{k}^{L}[n]$, can be rewritten as a generalized logistic function of the
UAV-GU elevation angle:

$\displaystyle
P_{k}^{L}[n]=A_{3}+\frac{A_{4}}{1+e^{-(A_{1}+A_{2}\theta_{k}[n])}},\quad\forall
n,k,$ (13)

where $A_{1}<0$, $A_{2}>0$, $A_{4}>0$ and $A_{3}$ are all constants, and
$A_{3}+A_{4}=1$. $\theta_{k}[n]$ is defined as the elevation angle between the
UAV and GU $k$ in time slot $n$, and

$\displaystyle\theta_{k}[n]=\frac{180}{\pi}\arctan\left(\frac{z[n]}{\|\mathbf{q}[n]-\mathbf{g}_{k}\|}\right).$
(14)

The NLoS probability can be obtained as $P_{k}^{N}[n]=1-P_{k}^{L}[n]$.
Considering the path loss and shadowing, the channel power gains in LoS and
NLoS states are expressed as

$\displaystyle h_{k}^{L}[n]=\rho_{0}d_{k}^{-\alpha_{L}}[n],\quad
h_{k}^{N}[n]=\mu_{0}\rho_{0}d_{k}^{-\alpha_{N}}[n],\quad\forall n,k,$ (15)

where $\rho_{0}$ denotes the channel power gain at the unit reference distance
in LoS state, $\mu_{0}$ is the additional signal attenuation factor due to the
more complex electromagnetic propagation environment in NLoS state [40].
$\alpha_{L}$ and $\alpha_{N}$ denote the path loss exponents for LoS and NLoS
states, respectively.

Assume that the UAV communicates with all GUs via the time-division multiple
access (TDMA) in the same bandwidth B, and transmits the signal with power
$P_{0}$. We define the binary transmission scheduling variable as
$\mathbf{A}=\\{a_{k}[n],\forall n\in\mathcal{N},k\in\mathcal{K}\\}$, where the
GU $k$ is served by the UAV in time slot $n$ if $a_{k}[n]=1$; otherwise,
$a_{k}[n]=0$. Without loss of generality, we assume that the UAV can only
communicate with one user in each time slot, leading to the following
constraints:

$\displaystyle a_{k}[n]\in{\\{0,1\\}},\quad\forall n,k,$ (16)
$\displaystyle\sum_{k=1}^{K}a_{k}[n]\leq 1,\quad\forall n.$ (17)

Thus, if GU $k$ is scheduled, an achievable rate of GU $k$ is given by:

$\displaystyle
R_{k}[n]=\log_{2}\left(1+\frac{h_{k}[n]P_{0}}{\sigma^{2}\Gamma}\right),\quad\forall
n,$ (18)

where $\sigma^{2}$ is the power of the additive white Gaussian noise (AWGN) at
the receiver, and $\Gamma$ is the signal-to-noise ratio (SNR) gap due to the
practical modulation-and-coding loss. Then, the achievable rates of GU $k$
conditioned on the LoS and NLoS states are respectively given by

$\displaystyle R_{k}^{L}[n]=\log_{2}(1+\gamma_{0}d_{k}[n]^{-\alpha_{L}}),\quad
R_{k}^{N}[n]=\log_{2}(1+\mu_{0}\gamma_{0}d_{k}[n]^{-\alpha_{N}}),\quad\forall
n,$ (19)

where $\gamma_{0}=\frac{\rho_{o}P_{o}}{\sigma^{2}\Gamma}$. The expected rate
of GU $k$ can be expressed as

$\displaystyle\mathbb{E}(R_{k}[n])=P^{L}_{k}[n]R^{L}_{k}[n]+(1-P^{L}_{k}[n])R^{N}_{k}[n],\quad\forall
n,$ (20)

where $\mathbb{E}(\cdot)$ denotes the expectation.

##  IV Problem Formulation

Considering the UAV-enabled downlink communication system where $K$ GUs are
randomly assigned to the fixed locations, we aim to maximize the EE of the UAV
by jointly optimizing the 3D trajectory $\mathbf{Q}$ and user scheduling
$\mathbf{A}$. We assume that all GUs’ locations and the probabilistic LoS
channel model are known a priori, and the UAV can estimate the instantaneous
channel state information (CSI) perfectly with individual GU in real-time
along its flight555Practically, the UAV can only estimate the CSI by receiving
signals from the GUs within its communication coverage. However, based on
proposed system model, the UAV will not transmit information to the GUs far
away from it. Thus such assumption does not compromise the practicability.[1].
To ensure the fairness of the UAV communication to each GU, we introduce a
relaxation variable $R_{\rm min}$, defined as the minimum threshold to the
expected achievable rate of each user over $N$ time slots. Then, the original
optimization problem can be formulated as:

$\displaystyle\underset{\mathbf{Q},\mathbf{A},R_{\rm min}}{\max}\quad$
$\displaystyle\frac{R_{\rm min}}{\sum_{n=1}^{N}P[n]}$ (21a) $\displaystyle
s.t.\quad$ $\displaystyle\mathbf{Q}[1]=\mathbf{Q}_{\rm
I},\quad\mathbf{Q}[N]=\mathbf{Q}_{\rm F},$ (21b) $\displaystyle R_{\rm
min}\leq\sum_{n=1}^{N}a_{k}[n]\mathbb{E}(R_{k}[n]),\quad\forall k,$ (21c)
$\displaystyle\eqref{constaint_v},\eqref{constaint_H},\eqref{constraint_a_0,1},\eqref{constraint_a<=1},$

since $P[n]$ and (21c) are non-convex w.r.t. $\mathbf{Q}$, and (16) gives
integer constraints, problem (21) is a mixed-integer non-convex fractional
program, which is challenging to solve in general. Additionally, problem (21)
is difficult to handle due to the lack of real-time wind information in
$P[n]$, i.e., speeds and directions in all time slots prior to flight. As
such, we propose an OBOA design with offline and online phases, shown in Fig.
3666From [1], the proposed algorithm in the offline phase can be executed on a
ground control center, which has enough computing capability and can exchange
information and control signal to the UAV. Then, the proposed algorithm in the
online phase can be executed in a UAV-mounted computer., to reduce the real-
time wind impact on UAV flight and obtain an energy-efficient solution to
problem (21). Specifically, the general descriptions of the two phases are as
follows:

  * •

Offline phase: To solve problem (21), before the flight, we assume that the
statistical distributions of both the wind speed and direction are
known777Before the flight, the UAV can be launched to hover at a fixed
altitude to collect wind information based on the UAV-mounted
electronic/mechanical anemometer and wind direction indicator or other high
precision sensors. Then, the parameters of the proposed distributions can be
estimated by existing estimation algorithms e.g., maximum likelihood
estimation. Similarly, the UAV can also collect real-time wind information
based on the mounted sensors in the online phase. The proposed design can be
extended to address other distributions of the wind.. With given wind
statistics, we propose an offline design to average the wind effect on the
UAV, by jointly optimizing the 3D trajectory $\mathbf{Q}$ and user scheduling
$\mathbf{A}$. The proposed offline design based on SP yields a maximum EE of
the UAV from a statistically favorable perspective. Note that the offline
design is obtained based on the wind statistics instead of the real-time wind
information. Therefore, the offline design leads to a suboptimal solution to
problem (21) which is used as a reference for the further online adaptive
optimization.

Figure 3: The proposed OBOA design for the UAV-enabled communication in wind

  * •

Online phase: During the flight, to make the UAV adapt to the changeable real-
time wind better and achieve a higher EE, we assume that the instantaneous
wind information can be perfectly measured and collected by the UAV. Then, an
OBOA design is proposed to minimize propulsion energy consumption slot by
slot, by optimizing UAV instantaneous velocity with the real-time wind
information and taking the obtained offline trajectory and user scheduling as
references. Based on the OBOA design, the UAV can dynamically adjusts its
velocity to further reduce the impact of the real-time wind and save more
propulsion energy. In this case, the proposed OBOA design leads to a more
energy-efficient solution to problem (21) compared with the offline design.

In the following sections, we describe the two phases of the proposed OBOA
design, including the offline design based on the given wind distribution and
the OBOA design with real-time wind information and offline references.

##  V Offline Design

In this section, to obtain a statistically suboptimal solution, we approximate
the original problem (21) into a tractable one, which is formulated as a
stochastic program. Further, we decompose this problem into three subproblems
and alternately optimize them to obtain an offline solution.

###  V-A Problem Approximation

To maximize the EE of the UAV, we need to relax the objective function in (21)
as the ratio of the minimum average rate $R_{\rm min}$ from each GU to the
upper bound of $P[n]$, denoted as $P^{\rm ub}[n]$. To this end, we tackle the
non-convexity of $P_{b}[n]$, $P_{i}[n]$, and
$m\|\mathbf{g}\|\|\mathbf{v}[n]\|\sin{\tau_{c}[n]}$ in $P[n]$. Specifically,
to deal with $P_{b}[n]$, we introduce the slack variables
$\mathbf{M}_{1}=\\{M_{1}[n]\\}_{n=1}^{N}$ and
$\mathbf{S}_{1}=\\{S_{1}[n]\\}_{n=1}^{N}$ such that $M_{1}[n]=\frac{1}{2\rho
A}(\|m\mathbf{a}[n]-m\mathbf{g}\|+\frac{1}{2}\rho S_{\rm
FP}\|(\mathbf{v}[n]-\mathbf{v}_{w}[n])\|^{2})\geq\frac{\|\mathbf{T}[n]\|}{2\rho
A}$ and
$S_{1}[n]\geq(\frac{2}{c_{T}}\eta[n]+3\|\mathbf{v}[n]\|^{2})\sqrt{\frac{2\eta[n]}{c_{T}}}$,
yielding

$\displaystyle\frac{c_{T}}{2}\frac{S_{1}^{2}[n]}{M_{1}[n]}\geq\left(\frac{2}{c_{T}}M_{1}[n]+3\|\mathbf{v}[n]\|^{2}\right)^{2},\quad\forall
n.$ (22)

Similarly, for $P_{i}[n]$, we introduce the slack variables
$\mathbf{M}_{2}=\\{M_{2}[n]\\}_{n=1}^{N}$ and
$\mathbf{S}_{2}=\\{S_{2}[n]\\}_{n=1}^{N}$ satisfying
$M_{2}[n]\geq\sqrt{M_{1}^{2}[n]+\frac{\|\mathbf{v}[n]\|^{4}}{4}}-\frac{\|\mathbf{v}[n]\|^{2}}{2}$
and $S_{2}[n]\geq M_{1}[n]\sqrt{M_{2}[n]}$. Hence we have:

$\displaystyle\frac{S_{2}^{2}[n]}{M_{2}[n]}\geq M_{1}^{2}[n],\quad\forall n,$
(23) $\displaystyle
M_{2}[n]+\|\mathbf{v}[n]\|^{2}\geq\frac{M_{1}^{2}[n]}{M_{2}[n]},\quad\forall
n.$ (24)

Note that the left-hand side (LHS) of (22), (23) and (24) are all convex and
differentiable functions w.r.t. the slack variables, which must hold with
equalities to obtain the optimal solution to problem (21). Otherwise
$M_{1}[n]$ and $M_{2}[n]$ can be increased and so can $S_{1}[n]$ and
$S_{2}[n]$ to further reduce the objective value of problem (21). To convert
the inequalities (22)-(24) into convex constraints, we apply the first-order
Taylor expansion to obtain their linear lower bounds, at given local points
denoted by $\mathbf{M}_{1}^{t}=\\{M_{1}^{t}[n]\\}_{n=1}^{N}$,
$\mathbf{M}_{2}^{t}=\\{M_{2}^{t}[n]\\}_{n=1}^{N}$,
$\mathbf{S}_{1}^{t}=\\{S_{1}^{t}[n]\\}_{n=1}^{N}$, and
$\mathbf{S}_{2}^{t}=\\{S_{2}^{t}[n]\\}_{n=1}^{N}$ in the _t_ -th iteration,
respectively. Then, we define
$\mathbf{v}^{t}[n]=\frac{\mathbf{Q}^{t}[n+1]-\mathbf{Q}^{t}[n]}{\Delta}$, and
have the following convex constraints:

$\displaystyle\left(\frac{2}{c_{T}}M_{1}[n]+3\|\mathbf{v}[n]\|^{2}\right)^{2}\leq$
$\displaystyle\frac{c_{T}}{2}\left(\frac{(S_{1}^{t}[n])^{2}}{M_{1}^{t}[n]}+\frac{2S_{1}^{t}[n](S_{1}[n]-S^{t}_{1}[n])}{M_{1}^{t}[n]}\right.$
$\displaystyle\left.-\left(\frac{S_{1}^{t}[n]}{M_{1}^{t}[n]}\right)^{2}(M_{1}[n]-M_{1}^{t}[n])\right),\quad\forall
n,$ (25) $\displaystyle\frac{M_{1}^{2}[n]}{M_{2}[n]}\leq
M_{2}[n]-\|\mathbf{v}^{t}[n]\|^{2}+2(\mathbf{v}^{t}[n])^{\mathrm{T}}\mathbf{v}[n],\quad\forall
n,$ (26) $\displaystyle
M_{1}^{2}[n]\leq\frac{S_{2}^{t}[n]}{M_{2}^{t}[n]}+\frac{2S_{2}^{t}[n](S_{2}[n]-S_{2}^{t}[n])}{M_{2}^{t}[n]}-\left(\frac{S_{2}^{t}[n]}{M_{2}^{t}[n]}\right)^{2}(M_{2}[n]-M_{2}^{t}[n]),\quad\forall
n.$ (27)

Furthermore, for $P_{c}[n]=m\|\mathbf{g}\|\|\mathbf{v}[n]\|\sin{\tau_{c}[n]}$,
we have

$\displaystyle m\|\mathbf{g}\|\|\mathbf{v}[n]\|\sin{\tau_{c}[n]}\leq
m\|\mathbf{g}\|\|\mathbf{v}[n]\|.$ (28)

Therefore, $P^{\rm ub}[n]$ can be rewritten as:

$\displaystyle P^{\rm ub}[n]=\frac{\delta\rho sA}{8}S_{1}[n]+2(1+c_{f})\rho
AS_{2}[n]+m\|\mathbf{g}\|\|\mathbf{v}[n]\|+\frac{1}{2}\rho S_{\rm
FP}\|\mathbf{v}[n]-\mathbf{v}_{w}[n]\|^{3}.$ (29)

We now tackle constraint (21c). From [39], with given GUs’ locations, the
achievable rate in LoS state is typically much higher than that in NLoS state.
To express the minimum expected rate, we take an approximation as:

$\displaystyle\mathbb{E}(R_{k}[n])$
$\displaystyle=P^{L}_{k}[n]R^{L}_{k}[n]+(1-P^{L}_{k}[n])R^{N}_{k}[n]\geq
P^{L}_{k}[n]R^{L}_{k}[n]$
$\displaystyle=\left(A_{3}+\frac{A_{4}}{1+e^{(A_{1}+A_{2}\theta_{k}[n])}}\right)\log_{2}(1+\gamma_{0}d_{k}^{-\alpha_{L}}[n])\triangleq\overline{R}_{k}[n],\quad\forall
n,k.$ (30)

Then, to address the non-affine relationship between $\mathbf{Q}$ and
$\overline{R}_{k}[n]$, combining (V-A)-(V-A), we reformulate problem (21) into
the following approximate form:

$\displaystyle\underset{\mathbf{Q},\mathbf{A},\mathbf{\Theta},\mathbf{M}_{1},\mathbf{M}_{2},\mathbf{S}_{1},\mathbf{S}_{2},R_{\rm
min}}{\max}\quad$ $\displaystyle\frac{R_{\rm min}}{\sum_{n=1}^{N}P^{\rm
ub}[n]}$ (31a) $\displaystyle s.t.\quad$ $\displaystyle
M_{1}[n]=\frac{1}{2\rho A}(\|m\mathbf{a}[n]-m\mathbf{g}\|+\frac{1}{2}\rho
S_{\rm FP}\|(\mathbf{v}[n]-\mathbf{v}_{w}[n])\|^{2}),\quad\forall n,$ (31b)
$\displaystyle\theta_{k}[n]\leq\frac{180}{\pi}\arctan\left(\frac{z[n]}{\|\mathbf{q}[n]-\mathbf{g}_{k}\|}\right),\quad\forall
k,n,$ (31c) $\displaystyle R_{\rm
min}\leq\sum_{n=1}^{N}a_{k}[n]\overline{R}_{k}[n],\quad\forall k,$ (31d)
$\displaystyle\eqref{constaint_v},\eqref{constaint_H},\eqref{constraint_a_0,1},\eqref{constraint_a<=1},\eqref{contraint_start_final},\eqref{constraint_S1M1_cvx},\eqref{constraint_M1M2_cvx},\eqref{constraint_S2M2_cvx},$

where $\mathbf{\Theta}=\\{\theta_{k}[n]\\}_{n=1}^{N},\forall k$. Note that
(31d) is still a non-convex constraint due to the coupled horizontal and
vertical trajectory variables of $\theta_{k}[n]$ in $\overline{R}_{k}[n]$.

###  V-B Proposed Offline Design on SP

To average the effect of $\mathbf{v}_{w}$, a stochastic fractional program
yielding a statistical maximum EE of the UAV is formulated as:

$\displaystyle\underset{\mathbf{Q},\mathbf{A},\mathbf{\Theta},\mathbf{M}_{1},\mathbf{M}_{2},\mathbf{S}_{1},\mathbf{S}_{2},R_{\rm
min}}{\max}\quad$ $\displaystyle\frac{R_{\rm min}}{\sum_{n=1}^{N}(P^{\rm
ub}_{b}[n]+P^{\rm ub}_{i}[n]+P^{\rm ub}_{c}[n]+\frac{1}{2}\rho S_{\rm
FP}\mathbb{E}(\|\mathbf{v}[n]-\mathbf{v}_{w}[n]\|^{3}))}$ (32a) $\displaystyle
s.t.\quad$ $\displaystyle M_{1}[n]=\frac{1}{2\rho
A}(\|m\mathbf{a}[n]-m\mathbf{g}\|+\frac{1}{2}\rho S_{\rm
FP}\mathbb{E}\|(\mathbf{v}[n]-\mathbf{v}_{w}[n])\|^{2}),\quad\forall n,$ (32b)
$\displaystyle\eqref{constaint_v},\eqref{constaint_H},\eqref{constraint_a_0,1},\eqref{constraint_a<=1},\eqref{contraint_start_final},\eqref{constraint_S1M1_cvx},\eqref{constraint_M1M2_cvx},\eqref{constraint_S2M2_cvx},\eqref{constraint_theta},\eqref{constraint_fairness_lb},$

where $P^{\rm ub}_{b}[n]=\frac{\delta\rho sA}{8}S_{1}[n]$, $P^{\rm
ub}_{i}[n]=2(1+c_{f})\rho AS_{2}[n]$ and $P^{\rm
ub}_{c}[n]=m\|\mathbf{g}\|\|\mathbf{v}[n]\|$. Note that $\mathbf{v}_{w}[n]$ is
a function of $v_{\rm ref}[n]$, $\beta[n]$ and $z[n]$ as shown in (3),
$\mathbf{v}[n]-\mathbf{v}_{w}[n]$ is convex w.r.t. $\mathbf{z}[n]$ and also
$\mathbf{Q}[n]$ when $p\leq 1$. Since the expectation preserves convexity,
$\mathbb{E}(\|\mathbf{v}[n]-\mathbf{v}_{w}[n]\|^{3})$ in $P^{\rm ub}[n]$ and
constraint (32b) are still convex.

We use the Monte Carlo sample average approximation (MCSAA), i.e., a method of
SP [41] to deal with the expectations in problem (32). Specifically, we take
an approximation as
$\mathbb{E}\\{f(\mathbf{v}[n],\mathbf{v}_{w}[n])\\}\approx\frac{1}{S}\sum_{i=1}^{S}f(\mathbf{v}[n],\mathbf{v}^{i}_{w}[n])$,
where $f$ is a general function of the optimization variable $\mathbf{v}[n]$
and the random variable $\mathbf{v}_{w}[n]$; $\mathbf{v}^{i}_{w}[n]$ is the
$i$-th sample generated from the given random distribution of the wind; $S$
denotes a finite set of random samples. Clearly, the above approximation
becomes accurate for a sufficiently large $S$.

To tackle the non-convexity of binary $\mathbf{A}$ and (31d), we next
decompose problem (32) into three subproblems, where the user scheduling
$\mathbf{A}$, the horizontal trajectory $\mathbf{q}$ and the vertical
trajectory $\mathbf{z}$ are optimized in an alternate fashion.

####  V-B1 User Scheduling Optimization

For given
$\hat{\mathbf{Q}}=\\{(\hat{\mathbf{q}}[n]^{\mathrm{T}},\hat{z}[n])\\}_{n=1}^{N}$,
the subproblem is equivalent to maximize $R_{\rm min}$ in (32a). To make the
subproblem tractable, we relax the integer variables in (16) into continuous
variables by considering that all $a_{k}[n]$ are continuous between 0 and 1,
yielding the following subproblem:

$\displaystyle\underset{\mathbf{A},R_{\rm min}}{\max}\quad$ $\displaystyle
R_{\rm min}$ (33a) $\displaystyle s.t.\quad$ $\displaystyle 0\leq a_{k}[n]\leq
1,\quad\forall n,k,$ (33b)
$\displaystyle\eqref{constraint_a<=1},\eqref{constraint_fairness_lb}.$

Since subproblem (33) is a linear program, it can be solved efficiently by
CVX.

####  V-B2 UAV Horizontal Trajectory Optimization

Given $\hat{\mathbf{A}}$ and $\hat{\mathbf{z}}$, the UAV horizontal trajectory
$\mathbf{q}$ can be optimized by solving the following subproblem:

$\displaystyle\underset{\mathbf{q},\mathbf{\Theta},\mathbf{M}_{1},\mathbf{M}_{2},\mathbf{S}_{1},\mathbf{S}_{2},R_{\rm
min}}{\max}\quad$ $\displaystyle\frac{R_{\rm min}}{\sum_{n=1}^{N}(P^{\rm
ub}_{b}[n]+P^{\rm ub}_{i}[n]+P^{\rm ub}_{c}[n]+\frac{1}{2}\rho S_{\rm
FP}\mathbb{E}(\|\mathbf{v}[n]-\mathbf{v}_{w}[n]\|^{3}))}$ (34a) $\displaystyle
s.t.\quad$
$\displaystyle\mathbf{Q[n]}=(\mathbf{q}^{\mathrm{T}}[n],\hat{z}[n]),\quad\forall
n,$ (34b)
$\displaystyle\theta_{k}[n]\leq\frac{180}{\pi}\arctan\left(\frac{\hat{z}[n]}{\|\mathbf{q}[n]-\mathbf{g}_{k}\|}\right),\quad\forall
k,n,$ (34c) $\displaystyle R_{\rm
min}\leq\sum_{n=1}^{N}\hat{a}_{k}[n]\overline{R}_{k}[n],\quad\forall k,$ (34d)
$\displaystyle\eqref{constaint_v},\eqref{contraint_start_final},\eqref{constraint_S1M1_cvx},\eqref{constraint_M1M2_cvx},\eqref{constraint_S2M2_cvx},\eqref{constraint_M1_SP}.$

To solve this subproblem, from [39, Lemma 1], we see that
$\overline{R}_{k}[n]$ is a convex function w.r.t.
$(1+e^{-(A_{1}+A_{2}\theta_{k}[n])})$ and
$(z[n]^{2}+\|\mathbf{q}[n]-\mathbf{g}_{k}\|^{2})$. Thus, we obtain a lower
bound of $\overline{R}_{k}[n]$, by using the first-order Taylor expansion for
any local UAV horizontal trajectory
$\mathbf{Q}^{t}[n]=((\mathbf{q}^{t}[n])^{\mathrm{T}},\hat{z}[n])$ in the
$t$-th iteration:

$\displaystyle\overline{R}_{k}[n]\geq$
$\displaystyle(A_{3}+\frac{A_{4}}{\xi^{t}_{k}[n]})\log_{2}\left(1+\frac{\gamma_{0}}{(\phi^{t}_{k}[n])^{\frac{1}{2}\alpha_{L}}}\right)-\frac{A_{4}\log_{2}e}{(\xi^{t}_{k}[n])^{2}}\ln\left(1+\frac{\gamma_{0}}{(\phi^{t}_{k}[n])^{\frac{1}{2}\alpha_{L}}}\right)$
$\displaystyle\times(e^{-u_{k}[n]}-e^{-u^{t}_{k}[n]})-(A_{3}+\frac{A_{4}}{\xi^{t}_{k}[n]})\left(\frac{\gamma_{0}\alpha_{L}}{2((\phi^{t}_{k}[n])^{\frac{1}{2}\alpha_{L}}+\gamma_{0})\phi^{t}_{k}[n]}\right)$
$\displaystyle\times(\|\mathbf{q}[n]-\mathbf{g}_{k}\|^{2}-\|\mathbf{q}^{t}[n]-\mathbf{g}_{k}\|^{2})\triangleq\overline{R}^{\rm
lb,\mathbf{q}}_{k}[n],$ (35)

where $u_{k}[n]\triangleq A_{1}+A_{2}\theta_{k}[n]$, $u^{t}_{k}[n]\triangleq
A_{1}+A_{2}\theta^{t}_{k}[n]$,
$\xi^{t}_{k}[n]=1+e^{-(A_{1}+A_{2}\theta^{t}_{k}[n])}$ and
$\phi^{t}_{k}[n]=\hat{z}[n]^{2}+\|\mathbf{q}^{t}[n]-\mathbf{g}_{k}\|^{2}$.
Similarly, for non-concave term
$\arctan\left(\frac{\hat{z}[n]}{\|\mathbf{q}[n]-\mathbf{g}_{k}\|}\right)$
w.r.t. $\mathbf{q}[n]$, since $\arctan(1/x)$ is convex for $x>0$, we obtain a
lower bound as

$\displaystyle\frac{180}{\pi}\arctan\left(\frac{\hat{z}[n]}{\|\mathbf{q}[n]-\mathbf{g}_{k}\|}\right)$
$\displaystyle\geq\frac{180}{\pi}\theta_{k}^{\rm lb}[n]\geq\theta_{k}[n],$
(36)

where

$\displaystyle\theta_{k}^{\rm
lb}[n]\triangleq\arctan\left(\frac{\hat{z}[n]}{\|\mathbf{q}^{t}[n]-\mathbf{g}_{k}\|}\right)-\frac{\hat{z}[n]}{\phi^{t}_{k}[n]}(\|\mathbf{q}[n]-\mathbf{g}_{k}\|-\|\mathbf{q}^{t}[n]-\mathbf{g}_{k}\|).$
(37)

Therefore, subproblem (34) can be recast as

$\displaystyle\underset{\mathbf{q},\mathbf{\Theta},\mathbf{U},\mathbf{M}_{1},\mathbf{M}_{2},\mathbf{S}_{1},\mathbf{S}_{2},R_{\rm
min}}{\max}\quad$ $\displaystyle\frac{R_{\rm min}}{\sum_{n=1}^{N}(P^{\rm
ub}_{b}[n]+P^{\rm ub}_{i}[n]+P^{\rm ub}_{c}[n]+\frac{1}{2}\rho S_{\rm
FP}\mathbb{E}(\|\mathbf{v}[n]-\mathbf{v}_{w}[n]\|^{3}))}$ (38a) $\displaystyle
s.t.\quad$ $\displaystyle R_{\rm
min}\leq\sum_{n=1}^{N}\hat{a}_{k}[n]\overline{R}^{\rm
lb,\mathbf{q}}_{k}[n],\quad\forall k,$ (38b)
$\displaystyle\theta_{k}[n]\leq\frac{180}{\pi}\theta^{\rm
lb}_{k}[n],\quad\forall k,n,$ (38c) $\displaystyle
u_{k}[n]=A_{1}+A_{2}\theta_{k}[n],\quad\forall k,n,$ (38d)
$\displaystyle\eqref{constaint_v},\eqref{contraint_start_final},\eqref{constraint_S1M1_cvx},\eqref{constraint_M1M2_cvx},\eqref{constraint_S2M2_cvx},\eqref{constraint_M1_SP},\eqref{constraint_Q_q},$

where $\mathbf{U}=\\{u_{k}[n]\\}_{n=1}^{N},\forall k$. Then, subproblem (38)
is a quasi-convex optimization problem, which can be solved efficiently via
successive convex approximation (SCA) and Dinkelbach’s method [42].

####  V-B3 UAV Vertical Trajectory Optimization

Given $\hat{\mathbf{A}}$ and $\hat{\mathbf{q}}$, the UAV vertical trajectory
$\mathbf{z}$ can be optimized by solving the following subproblem:

$\displaystyle\underset{\mathbf{z},\mathbf{\Theta},\mathbf{M}_{1},\mathbf{M}_{2},\mathbf{S}_{1},\mathbf{S}_{2},R_{\rm
min}}{\max}\quad$ $\displaystyle\frac{R_{\rm min}}{\sum_{n=1}^{N}(P^{\rm
ub}_{b}[n]+P^{\rm ub}_{i}[n]+P^{\rm ub}_{c}[n]+\frac{1}{2}\rho S_{\rm
FP}\mathbb{E}(\|\mathbf{v}[n]-\mathbf{v}_{w}[n]\|^{3}))}$ (39a) $\displaystyle
s.t.\quad$
$\displaystyle\mathbf{Q[n]}=(\hat{\mathbf{q}}^{\mathrm{T}}[n],z[n]),\quad\forall
n,$ (39b)
$\displaystyle\theta_{k}[n]\leq\frac{180}{\pi}\arctan\left(\frac{z[n]}{\|\hat{\mathbf{q}}[n]-\mathbf{g}_{k}\|}\right),\quad\forall
k,n,$ (39c) $\displaystyle R_{\rm
min}\leq\sum_{n=1}^{N}\hat{a}_{k}[n]\overline{R}_{k}[n],\quad\forall k,$ (39d)
$\displaystyle\eqref{constaint_v},\eqref{constaint_H},\eqref{contraint_start_final},\eqref{constraint_S1M1_cvx},\eqref{constraint_M1M2_cvx},\eqref{constraint_S2M2_cvx},\eqref{constraint_M1_SP}.$

Note that subproblem (39) is similar to (34), yielding the following
subproblem:

$\displaystyle\underset{\mathbf{z},\mathbf{U},\mathbf{\Theta},\mathbf{M}_{1},\mathbf{M}_{2},\mathbf{S}_{1},\mathbf{S}_{2},R_{\rm
min}}{\max}\quad$ $\displaystyle\frac{R_{\rm min}}{\sum_{n=1}^{N}(P^{\rm
ub}_{b}[n]+P^{\rm ub}_{i}[n]+P^{\rm ub}_{c}[n]+\frac{1}{2}\rho S_{\rm
FP}\mathbb{E}(\|\mathbf{v}[n]-\mathbf{v}_{w}[n]\|^{3}))}$ (40a) $\displaystyle
s.t.\quad$ $\displaystyle R_{\rm
min}\leq\sum_{n=1}^{N}\hat{a}_{k}[n]\overline{R}^{\rm
lb,\mathbf{z}}_{k}[n],\quad\forall k,$ (40b)
$\displaystyle\eqref{constaint_v},\eqref{constaint_H},\eqref{contraint_start_final},\eqref{constraint_S1M1_cvx},\eqref{constraint_M1M2_cvx},\eqref{constraint_S2M2_cvx},\eqref{constraint_M1_SP},\eqref{constraint_u},\eqref{constraint_Q_z},\eqref{constraint_theta_z},$

where $\overline{R}^{\rm lb,\mathbf{z}}_{k}[n]$ is the lower bound of
$\overline{R}_{k}[n]$ for any local UAV vertical trajectory
$\mathbf{Q}^{t}[n]=(\hat{\mathbf{q}}[n])^{\mathrm{T}},z^{t}[n])$ in $t$-th
iteration, which is similar to $\overline{R}^{\rm lb,\mathbf{q}}_{k}[n]$. Now
all constraints of subproblem (40a) are convex, and it can be solved similarly
as in (38).

Algorithm 1 Proposed Offline Design on SP

0: $\mathbf{g}_{k},k\in\mathcal{K}$, wind samples
$\mathbf{v}^{i}_{w}[n],n\in\mathcal{N},i\in\mathcal{S}$, convergence threshold
$\epsilon_{1}$ and $\epsilon_{2}$.

1: Initialization: Let $\mathbf{q}^{0},\mathbf{z}^{0}$, $\mathbf{A}^{0}$ be
feasible for problem (32); Set iteration index $r=0$;

2: repeat

3: Solve problem (33) with given $\mathbf{q}^{r}$ and $\mathbf{z}^{r}$, and
denote $\mathbf{A}^{r+1}$ by the optimal solution;

4: Solve problem (38) based on SCA and Dinkelbach’s method with given
$\mathbf{A}^{r+1}$ and $\mathbf{z}^{r}$, until its computed objective value
converges within $\epsilon_{1}>0$, and denote $\mathbf{q}^{r+1}$ by the
optimal solution;

5: Solve problem (40) based on the same methods in step 4 with given
$\mathbf{A}^{r+1}$ and $\mathbf{q}^{r+1}$, until its computed objective value
converges within $\epsilon_{1}>0$, and denote $\mathbf{z}^{r+1}$ by the
optimal solution;

6: Update $r=r+1$.

7: until The computed objective value of problem (32) converges within
$\epsilon_{2}>0$.

7: $\mathbf{Q[n]}=(\mathbf{q}^{\mathrm{T}}[n],z[n]),\forall n,\mathbf{A}$.

###  V-C Overall Algorithm and Computational Complexity

From the analysis in preceding subsections, the proposed offline design on SP
(summarized in Algorithm 1) provides the suboptimal offline trajectory and
user scheduling policy of problem (32) by alternately optimizing $\mathbf{A}$,
$\mathbf{q}$ and $\mathbf{z}$. Note that the user scheduling obtained by
solving subproblem (33) is continuous. To compute the objective value of
problem (32) accurately, we rounding the optimized user scheduling into a
binary solution after Algorithm 1 is completed [43]. Note that the decomposed
three subproblems satisfy that their optimal values all serve as a lower bound
of the optimal value of problem (32), and the values are non-decreasing after
each iteration. In addition, it is clear that the maximum objective value of
problem (32) is upper bounded by a finite value, which guarantees the
convergence of Algorithm 1 as analyzed in [44]. For the horizontal and
vertical trajectory optimizations, both subproblems involve $2N$ $l_{2}$-norm
terms with samples size $S$. Then, since all subproblems above are solved
based on the standard interior-point method, and based on the complexity
analysis method in [45], the required computational complexity of Algorithm 1
is thus
$\mathcal{O}\left(((KN+8N)^{3.5}+K^{1.5}S^{2}N^{2.5})\log(1/\epsilon_{1})\log(1/\epsilon_{2})\right)$.

##  VI Offline-Based Online Adaptive Design

The proposed offline design achieves an expected EE by jointly optimizing 3D
UAV trajectory and user scheduling over the whole flight period, based on the
given wind statistics. Although we obtain a statistically favorable suboptimal
solution to problem (21), the practical wind affecting UAV flight may be
significantly different from the predicted statistics, e.g., the speed and
direction of the real-time wind have large standard deviations, and change
dramatically. In this case, the obtained offline solution is not energy-
efficient enough.

In this section, the OBOA design with low complexity is proposed to make the
UAV better adapt to changeable real-time wind. Although the offline solution
is only statistically favorable, it still provides a suboptimal trajectory
design and user scheduling policy considering the real-time wind. Thus, in the
online phase, to take full advantage of the offline design and further improve
the EE, we denote the obtained trajectory and user scheduling policy in the
offline phase as $\mathbf{Q}^{\rm off}$ and $\mathbf{A}^{\rm off}$,
respectively, which are used as references for the OBOA design. In practice,
the complexity of the online adaptive design must be relatively low that we
can execute the algorithm within the given time slot. From the analysis of the
offline design, (38b) and other communication related constraints lead to high
computational complexity so that the proposed design is difficult to implement
in practice. Since the expected achievable rate is only related to the
distance and scheduling between the UAV and each GU, we make the following
assumptions: i) Under reasonable trajectory deviation, expected achievable
rate of the UAV with the online design is slightly lower than the rate with
the offline design. ii) The UAV can measure and collect the wind data in real
time, including the speed and direction, and we ignore the collection delay of
sensors deployed on the UAV.

Based on the assumptions above, the OBOA design can be implemented by only
optimizing the real-time velocity of the UAV to minimize its propulsion energy
consumption in each slot, regardless of the constraints related to the
communication. As a result, an online adaptation problem at each slot n can be
modified as

$\displaystyle\underset{\mathbf{v[n]}}{\min}\quad$ $\displaystyle{P[n]}$ (41a)
$\displaystyle s.t.\quad$
$\displaystyle\mathbf{Q}[n+1]=\mathbf{Q}[n]+\mathbf{v}[n]\Delta,$ (41b)
$\displaystyle\|\mathbf{Q}[n+1]-\mathbf{Q}^{\rm
off}[n+1]\|\leq\varepsilon_{Q},$ (41c)
$\displaystyle\|\mathbf{v}[n]-\mathbf{v}^{\rm off}[n]\|\leq\varepsilon_{v},$
(41d) $\displaystyle\|\mathbf{Q}[n+1]-\mathbf{Q}_{\rm
F}\|\leq\|\mathbf{Q}^{\rm off}[n+1]-\mathbf{Q}_{\rm F}\|,$ (41e)

where $\hat{\mathbf{v}}_{w}[n]$ denotes the wind measured by the UAV in real-
time in time slot $n$. $\mathbf{Q}[n]$ represents the current position of the
UAV, i.e., 3D position obtained by the online optimization in the last time
slot. $\varepsilon_{Q}$ and $\varepsilon_{v}$ are two tolerance parameters to
restrict the extents of deviation from the offline trajectory $\mathbf{Q}^{\rm
off}$ and velocity $\mathbf{v}^{\rm off}$, respectively, and the larger
tolerances mean that the UAV has more space to adjust. (41e) ensures that in
the whole online adaptation process, the UAV can dynamically adjust flight
direction and control its speed $\|\mathbf{v}[n]\|$ properly to avoid flying
at the maximum-endurance speed defined in [24] in each slot, so that
$\|\mathbf{v}[n]\|$ does not deviate much from $\|\mathbf{v}^{\rm off}[n]\|$.
Problem (41) is difficult to solve due to the non-convexity of $P[n]$. To
tackle this, we minimize its upper bound by the approximation analyzed in
Section v@-A. Then, problem (41) is rewritten as:

$\displaystyle\underset{\mathbf{v[n]},M_{1}[n],M_{2}[n],S_{1}[n],S_{2}[n]}{\min}\quad$
$\displaystyle{P^{\rm ub}_{b}[n]+P^{\rm ub}_{i}[n]+P^{\rm
ub}_{c}[n]+\frac{1}{2}\rho S_{\rm
FP}\|\mathbf{v}[n]-\hat{\mathbf{v}}_{w}[n]\|^{3}}$ (42a) $\displaystyle
s.t.\quad$ $\displaystyle M_{1}[n]=\frac{1}{2\rho
A}(\|m\mathbf{a}[n]-m\mathbf{g}\|+\frac{1}{2}\rho S_{\rm
FP}\|(\mathbf{v}[n]-\hat{\mathbf{v}}_{w}[n])\|^{2}),$ (42b)
$\displaystyle\eqref{constraint_S1M1_cvx},\eqref{constraint_M1M2_cvx},\eqref{constraint_S2M2_cvx},\eqref{constraint_online_Q},\eqref{constraint_online_distance},\eqref{constraint_online_v},\eqref{constraint_online_final},$

which can be solved efficiently by CVX. Note that the proposed adaptive design
requires to perform just once the SCA iteration in each given time slot $n$
(summarized in Algorithm 2), and the optimization variable dimension of
problem (42) is fixed at 7 (i.e., a three-dimensional vector $\mathbf{v}[n]$
and four one-dimensional slack variables). Since solving problem (42) in each
slot by using the CVX solver is based on the standard interior-point method,
its computational complexity scale as
$\mathcal{O}\left((7)^{3.5}\log(1/\epsilon_{1})\right)$, which is practically
affordable for common UAV-mounted computer. Then, the totally required
computational complexity of Algorithm 2 in the whole flight duration is
$\mathcal{O}\left((N-1)(7)^{3.5}\log(1/\epsilon_{1})\right)$, which do not
increase with the number of GUs $K$. Therefore, the proposed OBOA design can
be implemented practically within a given time slot888The average running time
of the OBOA design is about 0.5 s for each time slot. The result is obtained
by using CVX and running with 3.6 GHz CPU and 8 GB RAM memory..

Algorithm 2 Proposed OBOA Design

0: $\varepsilon_{Q}$, $\varepsilon_{v}$, offline references $\mathbf{Q}^{\rm
off}$, $\mathbf{v}^{\rm off}$ and $\mathbf{A}^{\rm off}$, convergence
threshold $\epsilon_{3}$.

1: Initialization: Set $\mathbf{v}^{0}[n]=\mathbf{v}^{\rm off}[n]$,
$\mathbf{Q}[1]=\mathbf{Q}_{{\rm I}}$, $\mathbf{Q}[N]=\mathbf{Q}_{{\rm F}}$ and
SCA iteration index $t=0$;

2: for $n=1$ to $N-1$ do

3: Collect the instantaneous wind information $\hat{\mathbf{v}}_{w}[n]$;

4: Solve problem (42) based on SCA with $\mathbf{Q}[n]$, $\mathbf{Q}^{\rm
off}[n+1]$, $\mathbf{v}^{\rm off}[n]$ and $\hat{\mathbf{v}}_{w}[n]$, until its
computed objective value converges within $\epsilon_{3}>0$, and denote
$\mathbf{v}^{\rm opt}[n]$ by the optimal solution;

5: Obtain $\mathbf{Q}[n+1]=\mathbf{Q}[n]+\mathbf{v}^{\rm opt}[n]\Delta$.

6: end for

6: $\mathbf{Q},\mathbf{v}^{\rm opt}$.

7: UAV flies following the obtained $\mathbf{Q}$ and $\mathbf{v}^{\rm opt}$,
and communicates to each GU with the scheduling policy $\mathbf{A}^{\rm off}$.

It is worth mentioning that, as shown in step 7 of Algorithm 2, the UAV flies
following the obtained $\mathbf{Q}$ and $\mathbf{v}^{\rm opt}$, and
communicates to each GU with the scheduling policy $\mathbf{A}^{\rm off}$.
Since we consider the real-time wind affecting UAV flight in the online
adaptive optimization, we also consider the practical channel model depending
on the real environment instead of the probability, to obtain the achievable
rate and also the EE of the UAV more accurately in the online phase.
Specifically, the achievable rate $R_{k}[n]$ is computed by $R^{L}_{k}[n]$
when the UAV-GU link is conditioned on the LoS state, and computed by
$R^{N}_{k}[n]$ when buildings block the UAV-GU link. The practical channel in
the LoS/NLoS states is based on a concrete Manhattan city model[46], as shown
in the following numerical results.

TABLE I: Simulation Parameters Parameter |  Description |  Value
---|---|---
$T_{0}$, $\Delta$ |  Flight duration and time slot |  150 s, 1 s
$\mathbf{Q}_{\rm I},\mathbf{Q}_{\rm F}$ |  Initial and final position |  $(0,500,100)$ m, $(1000,500,100)$ m
$V^{\rm max}_{H},V^{\rm max}_{V}$ |  Maximum horizontal and vertical flight speeds |  $40$ m/s, $20$ m/s
$H_{\rm max}$, $H_{\rm min}$ |  Maximum and minimum flight altitudes |  $300$ m, $50$ m
$K$, $\mathbf{g}_{k}$ |  Number of GUs and corresponding horizontal coordinates |  4, $(100,300)$ m, $(500,800)$ m, $(500,200)$ m, $(900,600)$ m
$A_{1}-A_{4}$ |  Probabilistic LoS channel parameters |  -1, 0.05, 0.1, 0.9
$\rho_{0}$, $\mu_{0}$, $\Gamma$, $\alpha_{L}$, $\alpha_{N}$, $\sigma^{2}$, $B$, $P_{0}$ |  Communication parameters |  -60 dB, -20 dB, 8.2 dB, 2.5, 5, -110 dBm, 1 MHz, 0.1 W
$m$, $\|\mathbf{g}\|$ |  Weight of the UAV and gravitational acceleration |  2 $\rm kg$, 9.8 $\rm m/s^{2}$
$\rho$, $A$, $\delta$, $s$, $S_{FP}$, $c_{T}$, $c_{f}$ |  UAV aerodynamics parameters |  1.225 $\rm kg/m^{3}$, 0.79 $\rm m^{2}$, 0.012, 0.1, 0.01 $\rm m^{2}$, 0.3, 0.13
$h_{\rm ref}$ |  Reference altitude |  50 m
$p$ |  Empirical exponent of the wind modeling |  0.5
$\varepsilon_{v}$ |  Speed tolerance |  20 m/s
$\epsilon_{1}$, $\epsilon_{2}$, $\epsilon_{3}$, $S$ |  Convergence thresholds and size of random wind samples set |  $10^{-3}$, 300

##  VII Numerical Results

In this section, numerical simulations are conducted to verify the
effectiveness of the proposed offline and OBOA designs compared with the
benchmark scheme (denoted as Windless). The trajectory and user scheduling of
windless scheme are optimized by following the similar steps in Algorithm 1
considering a windless environment, i.e., solving the corresponding problems
in Algorithm 1 based on the GPECM with $\|\mathbf{v}_{w}\|=0$. The parameter
setting for the simulations are summarized in Table I. Note that the
parameters about communication and probabilistic LoS channel are similar to
those in [39, 47], and the UAV aerodynamics parameters are set similarly as in
[24, 25, 26].

###  VII-A Performance Analysis for Offline Design

For the proposed offline design, we use MCSAA to average the wind effect.
Thus, the optimized trajectories (i.e., Fig. 4 and 5) are obtained from a
statistical perspective based on corresponding wind distributions.
Additionally, the expected rates of the UAV are computed by the probabilistic
LoS channel model as seen in (V-A), and the total propulsion energy of offline
design and windless scheme are computed by repeating $S=300$ independent
experiments for randomized wind and averaging their numerical results. Then,
we have the EE comparisons (i.e., Fig. 7 and 7) as follows.

(a)

(b)

Figure 4:  UAV trajectories based on different $\mu$ ($\kappa=20$,
$\lambda=10$, $c=10$)

Fig. 4 shows the optimized UAV trajectories based on different values of
$\mu$. From the 3D view, we see that, the UAV tends to communicate with the
GUs at a higher flight altitude in a windless environment. This is because a
higher flight altitude leads to a higher LoS probability, and hence a higher
EE. As approaching each GU, the UAV lowers its flight height. Because a
shorter horizontal distance between the UAV and each GU leads to a higher LoS
probability, and also a higher EE. To probe into the impact of wind direction,
we notice the UAV trajectories from GU 2, i.e., $(500,800)$ m to GU 3, i.e.,
$(500,200)$ m in the 3D view. Along this direction, the UAV flight direction
is approximately equal to 90 degrees as shown in the 2D view in Fig. 4. From
the definition in Section ii@, $\mu$ is the expectation angle of the wind. At
this point, the winds with $\mu=90$ and $\mu=270$ are respectively the
tailwind and the headwind for the UAV. Therefore, the UAV flies lower (seen in
the red trajectory) to withstand the headwind, and save more propulsion energy
to achieve a higher EE. On the contrary, the UAV increases its flight height
(seen in the black trajectory) to achieve a higher expected achievable rate
and a higher EE with the tailwind. From the 2D view in Fig. 4, we see that all
trajectories are slightly different horizontally. Specifically, as shown in
the magnified part of the 2D view, the green trajectory is more biased in the
direction of 180 degrees but the blue is opposite, due to the wind effect in
the corresponding direction.

(a)

(b)

Figure 5: UAV trajectories based on different $\lambda$ ($\mu=270$,
$\kappa=20$, $c=10$)

Fig. 5 shows the optimized UAV trajectories based on different $\lambda$.
Similar to the 3D view in Fig. 4, without the wind impact, the UAV tends to
fly higher. Additionally, we see that the average height of the UAV decreases
with a larger $\lambda$, because the UAV has to fly at a lower altitude to
reduce the wind impact, which becomes stronger at a higher altitude, and save
more propulsion energy. For the 2D view in Fig. 5, all trajectories are also
slightly different similar to Fig. 4.

Figure 6:  EE comparison based on different $\mu$ ($\kappa=20$, $\lambda=10$,
$c=10$)

Figure 7:  EE comparison based on different $\lambda$ ($\mu=270$, $\kappa=20$,
$c=10$)

Fig. 7 and Fig. 7 show the EE comparisons with varying $\mu$ and $\lambda$. As
expected, under the impact of the wind, the proposed offline design based on
SP greatly improves the EE compared with the windless scheme. Specifically,
the EEs vary greatly over $\mu$, due to the wind impact in different
directions. In addition, as $\lambda$ increases, the EEs decrease, and the
offline design outperforms significantly compared with the windless scheme.
Therefore, the above results demonstrate the validity of the proposed offline
design.

###  VII-B Performance Analysis for OBOA Design

For the proposed OBOA design, we consider the real environment with real-time
wind effect and practical channel model. Clearly, the results of the OBOA
design are highly correlated with particular parameters settings. Thus, to
compare the performances of the proposed three schemes more accurately in the
real environment, the achievable rates for all schemes (e.g., Fig. 11) are
obtained by repeating $S=300$ independent experiments for the randomized
Manhattan city model and averaging their numerical results. Similarly, the
total propulsion energy of all schemes (e.g., Fig. 11) is computed by
repeating $S=300$ independent experiments for randomized wind and averaging
their numerical results. Then, we have the EE comparisons (i.e., Fig. 13 and
Fig. 13) as follows. Also, the optimized trajectory with the OBOA design and
Manhattan city model in Fig. 8, and Fig. 9 are generated by an independent
experiment for the randomized Manhattan city model and wind. The corresponding
conclusion does not rely on any specific random parameters settings.

(a)

(b)

Figure 8: UAV trajectories based on different schemes ($\varepsilon_{Q}=50$ m,
$\mu=180$, $\kappa=5$, $\lambda=10$, $c=5$)

Fig. 8 shows the UAV trajectories comparison. Specifically, from the 3D view,
we see that the flight heights of the UAV with the OBOA and offline designs
are lower than the height with the windless scheme as the analysis above.
Furthermore, the optimized trajectory with the OBOA design is different from
the offline design in each time slot, due to the adaptation to the real-time
wind changing dramatically. To better verify this, we continue to analyze the
2D view in Fig. 8 and Fig. 9 in the following and observe the trajectory and
velocity changes of the UAV based on the online adaptation.

(a) $\mu=90$ (b) $\mu=270$

Figure 9:  Angular deviation of the velocity with two designs w.r.t. real-time
wind in each slot ($\varepsilon_{Q}=100$ m, $\kappa=5$, $\lambda=10$, $c=5$)

From the magnified part in the 2D view in Fig. 8, we see that the optimized
trajectories with the offline design and the windless scheme are both smooth,
but the UAV with the OBOA design continuously adjusts its velocity in each
slot. The main reason is that the trajectories with the offline design and the
windless scheme are globally optimized based on MCSAA from a statistical
perspective, and the wind effect on the UAV is averaged in this case. On the
contrary, the trajectory with the OBOA design is constructed by a series of
independent optimizations in each time slot. Furthermore, from Fig. 9, we see
that the angular deviation of the OBOA design is less than that of the offline
design in most time slots, yielding a smaller wind impact on UAV flight. Thus,
with the OBOA design, the UAV adjusts its current flight velocity in
adaptation to the real-time wind under the corresponding flight constraints,
such that the OBOA design can better deal with the impact of the wind
randomness, and save more propulsion energy of the UAV to achieve a higher EE.

Figure 10:  Achievable rate comparison with different $\varepsilon_{Q}$
($\mu=180$, $\kappa=5$, $\lambda=10$, $c=5$)

Figure 11:  Propulsion energy consumption comparison with different
$\varepsilon_{Q}$ ($\mu=180$, $\kappa=5$, $\lambda=10$, $c=5$)

Fig. 11 and Fig. 11 show the expected achievable rate and propulsion energy
consumption comparisons of the UAV under varying distance constraint. Fig. 11
demonstrates the assumption in the online phase. The expected achievable rate
with the OBOA design is lower than the rate with the offline design. However,
the difference of the expected achievable rate between the OBOA and offline
designs is no more than $10\%$ at best at $\varepsilon_{Q}=100$ $\rm m$. The
OBOA design in Fig. 11 saves propulsion energy of the UAV significantly
compared with the offline design, and the difference of energy consumption
goes up to $20\%$ at best at $\varepsilon_{Q}=100$ $\rm m$. As we pay more
attention to the improvement of the EE, Fig. 11 and Fig. 11 verify the
effectiveness of the OBOA design.

Figure 12:  EE comparison with different $\kappa$ ($\varepsilon_{Q}=100\rm m$
, $\mu=180$, $\lambda=10$, $c=5$)

Figure 13:  EE comparison with different $c$ ($\varepsilon_{Q}=100\rm m$ ,
$\mu=180$, $\lambda=10$, $\kappa=5$)

Fig. 13 and Fig. 13 show the EE comparisons with different levels of wind
randomness by varying $\kappa$ and $c$. As is shown in Fig. 13, we see that
the OBOA and offline designs outperforms the windless scheme in terms of EE,
and the OBOA design is more energy-efficient than the offline design.
Specifically, when $\kappa\geq 3$, all three EEs gradually go down with
decreasing $\kappa$. Because the wind varying more frequently in direction
makes the UAV difficult to make timely adjustments for the flight direction.
However, when $\kappa\leq 3$, the EEs increase as $\kappa$ decreases. The main
reason is that, a smaller $\kappa$, i.e., the direction of the wind changes
more frequently, leads to a higher tailwind probability for the UAV. In this
case, the UAV can save more propulsion energy compared with the condition that
the direction of the wind is more concentrated. For Fig. 13, the general
conclusion about the EE comparison is similar to Fig. 13. In addition, we see
that the EEs go down with decreasing $c$ in Fig. 13. This is because frequent
changes in wind speed, especially when the wind speed becomes higher, will
cause more propulsion energy consumption for the UAV. To sum up, the above
results corroborate that the proposed OBOA design is more energy-efficient,
which makes the UAV better adapt to the real-time wind changing more
dramatically.

##  VIII Conclusion

In this paper, we investigated the effect of the stochastic wind on the UAV-
enabled communication system. We first derived a 3D GPECM of a rotary-wing UAV
affected by the stochastic wind. Then, we formulated a stochastic optimization
problem to maximize the EE of the UAV, by jointly optimizing the 3D trajectory
and user scheduling. To tackle this problem, we proposed an OBOA design
yielding a statistical solution based on wind statistics and SP technique in
the offline phase, and further optimized UAV instantaneous velocity in
adaptation to the real-time wind in each slot based on the offline design in
the online phase. Numerical results show the optimized UAV trajectories of
both two phases can better adapt to the wind in changing speed and direction,
and achieve a higher EE compared to the benchmark scheme. Moreover, the
proposed OBOA design can be applied in the scenario with dramatic wind changes
and achieves a better performance in terms of EE.

##  Appendix A Simplifications For the GPECM

Since the existing models in [24, 25, 26] all lack the analysis of the impact
of wind on the UAV propulsion power, we first simplify the GPECM into a form
without considering the wind effect. Specifically, the GPECM is derived based
on the 3D force analysis on a rotary-wing UAV, and winds exert impact on the
UAV as shown in (4) and (5). Hence, if we neglect the wind effect, i.e., by
setting $\|\mathbf{v}_{w}\|=0$ in (4) and (5), the GPECM without considering
the wind effect which has a similar form as (8), can be expressed as

$\displaystyle P=$
$\displaystyle\frac{\delta}{8}\left(\frac{\|\mathbf{T}\|}{c_{T}\rho
A}+3\|\mathbf{v}\|^{2}\right)\sqrt{\frac{\rho
s^{2}A\|\mathbf{T}\|}{c_{T}}}+(1+c_{f})\|\mathbf{T}\|\left(\sqrt{\frac{\|\mathbf{T}\|^{2}}{(2\rho
A)^{2}}+\frac{\|\mathbf{v}\|^{4}}{4}}-\frac{\|\mathbf{v}\|^{2}}{2}\right)^{\frac{1}{2}}$
$\displaystyle+m\|\mathbf{g}\|\|\mathbf{v}\|\sin{\tau_{c}}+\frac{1}{2}\rho
S_{\rm FP}\|\mathbf{v}\|^{3},$ (43)

where $\|\mathbf{T}\|=\|m\mathbf{a}+\frac{1}{2}\rho S_{\rm
FP}\|\mathbf{v}\|\mathbf{v}-m\mathbf{g}\|$. We now clarify how (A) reduces to
the existing models under different assumptions.

####  A-1 Simplification for the GPECM to the model in [24]

In [24], a 2D propulsion energy consumption model is derived, which a function
of UAV velocity and other constant fuselage parameters. Specifically, the
authors of [24] only considered the impact of the horizontal velocity and the
flight drag on UAV propulsion power, and assumed that the thrust always
equaled to the weight considering the horizontal flight. Thus, we have
$\mathbf{T}=m\mathbf{g}$, and can omit the terms $m\mathbf{a}$ and climbing
power $m\|\mathbf{g}\|\|\mathbf{v}\|\sin{\tau_{c}}$ in (A). Then, based on
$v_{o}=\sqrt{\frac{\|\mathbf{T}\|}{2\rho
A}}=\sqrt{\frac{\|m\mathbf{g}\|}{2\rho A}}$ and $d_{0}=\frac{S_{\rm FP}}{sA}$,
(A) is recast as

$\displaystyle P=$
$\displaystyle\frac{\delta}{8}\left(\frac{\|\mathbf{T}\|}{c_{T}\rho
A}+3\|\mathbf{v}\|^{2}\right)\sqrt{\frac{\rho
s^{2}A\|\mathbf{T}\|}{c_{T}}}+(1+c_{f})\frac{\|m\mathbf{g}\|^{\frac{3}{2}}}{\sqrt{2\rho
A}}\kappa\left(\sqrt{\kappa^{2}+\frac{\|\mathbf{v}\|^{4}}{4v_{o}^{4}}}-\frac{\|\mathbf{v}\|^{2}}{2v_{o}^{2}}\right)^{\frac{1}{2}}+\frac{1}{2}d_{o}\rho
sA\|\mathbf{v}\|^{3},$ (44)

where $\kappa=\|\mathbf{T}\|/\|m\mathbf{g}\|$ denotes the thrust-to-weight
ratio defined in [24], and we have $\kappa=1$ based the above assumption.
Then, based on $\|\mathbf{T}\|=c_{T}\rho Av_{\rm tip}^{2}$ and
$\hat{v}=\frac{\|\mathbf{v}\|}{v_{\rm tip}}$, (44) can be rewritten as

$\displaystyle P=$
$\displaystyle\frac{\delta}{8}\left(1+3\frac{\|\mathbf{v}\|^{2}}{v^{2}_{\rm
tip}}\right)+(1+c_{f})\frac{\|m\mathbf{g}\|^{\frac{3}{2}}}{\sqrt{2\rho
A}}\left(\sqrt{1+\frac{\|\mathbf{v}\|^{4}}{4v_{o}^{4}}}-\frac{\|\mathbf{v}\|^{2}}{2v_{o}^{2}}\right)^{\frac{1}{2}}+\frac{1}{2}d_{o}\rho
sA\|\mathbf{v}\|^{3}.$ (45)

This is exactly the same as model (12) in [24].

####  A-2 Simplification for the GPECM to the model in [25]

In [25], an improved 2D propulsion energy consumption model based on (45) is
derived. The authors of [25] considered the impact of thrust-to-weight ratio
$\kappa$ and obtain a model w.r.t. the velocity, acceleration and direction
change in the horizontal direction of the UAV. Note the assumptions of the
model in [25] are the same as in [R15] besides considering $\kappa\neq 1$.
Thus (A) is reduced to

$\displaystyle P=$
$\displaystyle\frac{\delta}{8}\left(1+3\frac{\|\mathbf{v}\|^{2}}{v^{2}_{\rm
tip}}\right)+(1+c_{f})\frac{\|m\mathbf{g}\|^{\frac{3}{2}}}{\sqrt{2\rho
A}}\kappa\left(\sqrt{\kappa^{2}+\frac{\|\mathbf{v}\|^{4}}{4v_{o}^{4}}}-\frac{\|\mathbf{v}\|^{2}}{2v_{o}^{2}}\right)^{\frac{1}{2}}+\frac{1}{2}d_{o}\rho
sA\|\mathbf{v}\|^{3}.$ (46)

From Fig. 1 in [25], we similarly denote $\phi$ as the angle between
$\mathbf{T}$ and $Z$-axis, and denote $\varphi$ as the angle between the drag
$\mathbf{D}$ and $\mathbf{F}$ (i.e., the resultant force on the UAV). Note
that $\mathbf{D}$ and $\mathbf{F}$ both are forces in the horizontal
direction. Then, based on $\|\mathbf{D}\|=\frac{1}{2}\rho S_{\rm
FP}\|\mathbf{v}\|^{3}$ and the law of cosines, we take the force analysis on
the UAV in both horizontal and vertical direction, yielding

$\displaystyle\|\mathbf{T}\||\sin\phi|=\sqrt{\|\mathbf{F}\|^{2}+\|\mathbf{D}\|^{2}-2\|\mathbf{F}\|\|\mathbf{D}\|\cos\varphi},\quad\|\mathbf{T}\||\cos\phi|=\|m\mathbf{g}\|.$
(47)

Then, based on $|\sin\phi|^{2}+|\cos\phi|^{2}=1$ and
$\kappa=\|\mathbf{T}\|/\|m\mathbf{g}\|$, we have

$\displaystyle\kappa=\sqrt{1+\frac{4m^{2}\|\mathbf{a}\|^{2}+\rho^{2}S_{FP}^{2}\|\mathbf{v}\|^{4}+4m\rho
S_{FP}\mathbf{a}\mathbf{v}\|\mathbf{v}\|}{4\|m\mathbf{g}\|^{2}}}.$ (48)

Therefore, (A) is reduced to equation (5) in [25] by substituting (48) into
(46).

####  A-3 Simplification for the GPECM to the model in [26]

In [26], the authors proposed a 3D propulsion energy consumption model by
considering the propulsion energy consumption in horizontal and vertical
direction as two independent components. In addition, the authors of [26] also
neglected the impact of the vertical flight drag and acceleration on the UAV
propulsion power, similarly as in [R15]. Then from [26], we similarly
decompose the velocity $\mathbf{v}$ into three components in $X$, $Y$ and
$Z$-axis, respectively, i.e., $\mathbf{v}=(v_{x},v_{y},v_{z})$. Since
$\tau_{c}$ denotes the climbing angle and $\|\mathbf{v}\|$ in (45) is the
horizontal velocity of the UAV, we have $\|\mathbf{v}\|\sin{\tau_{c}}=v_{z}$.
Thus, (A) can be recast as

$\displaystyle P=$
$\displaystyle\frac{\delta}{8}\left(1+\frac{3(v_{x}^{2}+v_{y}^{2})}{v^{2}_{\rm
tip}}\right)+(1+c_{f})\frac{\|m\mathbf{g}\|^{\frac{3}{2}}}{\sqrt{2\rho
A}}\left(\sqrt{1+\frac{(v_{x}^{2}+v_{y}^{2})^{2}}{4v_{o}^{4}}}-\frac{v_{x}^{2}+v_{y}^{2}}{2v_{o}^{2}}\right)^{\frac{1}{2}}$
$\displaystyle+m\|\mathbf{g}\|v_{z}+\frac{1}{2}d_{o}\rho
sA(v_{x}^{2}+v_{y}^{2})^{\frac{3}{2}}.$ (49)

Therefore, (A) is reduced to the form as equation (17) in [26].

## References

  * [1] Y. Zeng, Q. Wu, and R. Zhang, “Accessing from the sky: A tutorial on UAV communications for 5G and beyond,” _Proceedings of the IEEE_ , vol. 107, no. 12, pp. 2327–2375, 2019.
  * [2] G. Geraci, A. Garcia-Rodriguez, M. M. Azari, A. Lozano, M. Mezzavilla, S. Chatzinotas, Y. Chen, S. Rangan, and M. Di Renzo, “What will the future of UAV cellular communications be? a flight from 5G to 6G,” _IEEE Commun. Surveys Tuts._ , vol. 24, no. 3, pp. 1304–1335, 2022.
  * [3] C. Pan, H. Ren, K. Wang, J. F. Kolb, M. Elkashlan, M. Chen, M. Di Renzo, Y. Hao, J. Wang, A. L. Swindlehurst _et al._ , “Reconfigurable intelligent surfaces for 6G systems: Principles, applications, and research directions,” _IEEE Commun. Mag._ , vol. 59, no. 6, pp. 14–20, 2021.
  * [4] Y. Zeng, J. Lyu, and R. Zhang, “Cellular-connected UAV: Potential, challenges, and promising technologies,” _IEEE Wireless Commun._ , vol. 26, no. 1, pp. 120–127, 2018.
  * [5] A. A. Khuwaja, Y. Chen, N. Zhao, M.-S. Alouini, and P. Dobbins, “A survey of channel modeling for UAV communications,” _IEEE Commun. Surveys Tuts._ , vol. 20, no. 4, pp. 2804–2821, 2018.
  * [6] M. Mozaffari, W. Saad, M. Bennis, Y.-H. Nam, and M. Debbah, “A tutorial on UAVs for wireless networks: Applications, challenges, and open problems,” _IEEE Commun. Surveys Tuts._ , vol. 21, no. 3, pp. 2334–2360, 2019.
  * [7] M. Di Renzo, A. Zappone, M. Debbah, M.-S. Alouini, C. Yuen, J. De Rosny, and S. Tretyakov, “Smart radio environments empowered by reconfigurable intelligent surfaces: How it works, state of research, and the road ahead,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 11, pp. 2450–2525, 2020.
  * [8] Y. Wang, Z. Gao, J. Zhang, X. Cao, D. Zheng, Y. Gao, D. W. K. Ng, and M. Di Renzo, “Trajectory design for UAV-based internet of things data collection: A deep reinforcement learning approach,” _IEEE Internet Things J._ , vol. 9, no. 5, pp. 3899–3912, 2021.
  * [9] S. Hayat, E. Yanmaz, and R. Muzaffar, “Survey on unmanned aerial vehicle networks for civil applications: A communications viewpoint,” _IEEE Commun. Surveys Tuts._ , vol. 18, no. 4, pp. 2624–2661, 2016.
  * [10] S. Li, B. Duo, X. Yuan, Y.-C. Liang, and M. Di Renzo, “Reconfigurable intelligent surface assisted UAV communication: Joint trajectory design and passive beamforming,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 5, pp. 716–720, 2020.
  * [11] B. Duo, Q. Wu, X. Yuan, and R. Zhang, “Energy efficiency maximization for full-duplex UAV secrecy communication,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 4, pp. 4590–4595, 2020.
  * [12] Z. Zhai, X. Dai, B. Duo, X. Wang, and X. Yuan, “Energy-efficient UAV-mounted RIS assisted mobile edge computing,” _IEEE Wireless Commun. Lett._ , vol. 11, no. 12, pp. 2507–2511, 2022.
  * [13] S. Li, B. Duo, M. Di Renzo, M. Tao, and X. Yuan, “Robust secure UAV communications with the aid of reconfigurable intelligent surfaces,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 10, pp. 6402–6417, 2021.
  * [14] M. Erdelj, E. Natalizio, K. R. Chowdhury, and I. F. Akyildiz, “Help from the sky: Leveraging UAVs for disaster management,” _IEEE Pervasive Comput._ , vol. 16, no. 1, pp. 24–32, 2017.
  * [15] Y. Lei, Y. Liu, Q. Wu, X. Yuan, J. Ning, and B. Duo, “Enhancing UAV-enabled communications via multiple intelligent omni-surfaces,” _IEEE Commun. Lett._ , 2022.
  * [16] J. Lyu, Y. Zeng, R. Zhang, and T. J. Lim, “Placement optimization of UAV-mounted mobile base stations,” _IEEE Commun. Lett._ , vol. 21, no. 3, pp. 604–607, 2016.
  * [17] 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, 2016.
  * [18] A.-A. A. Boulogeorgos, A. Alexiou, and M. Di Renzo, “Outage performance analysis of RIS-assisted UAV wireless systems under disorientation and misalignment,” _arXiv preprint arXiv:2201.12056_ , 2022.
  * [19] L. Yang, F. Meng, J. Zhang, M. O. Hasna, and M. Di Renzo, “On the performance of RIS-assisted dual-hop UAV communication systems,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 9, pp. 10 385–10 390, 2020.
  * [20] Y. Zeng and R. Zhang, “Wireless communications with unmanned aerial vehicles: Opportunities and challenges,” _IEEE Commun. Mag._ , vol. 54, no. 5, pp. 36–42, 2016.
  * [21] D. Yang, Q. Wu, Y. Zeng, and R. Zhang, “Energy tradeoff in ground-to-UAV communication via trajectory design,” _IEEE Transactions on Vehicular Technology_ , vol. 67, no. 7, pp. 6721–6726, 2018.
  * [22] X. Pang, J. Tang, N. Zhao, X. Zhang, and Y. Qian, “Energy-efficient design for mmwave-enabled NOMA-UAV networks,” _Sci. China Inf. Sci._ , vol. 64, pp. 1–14, 2021.
  * [23] Y. Su, X. Pang, S. Chen, X. Jiang, N. Zhao, and F. R. Yu, “Spectrum and energy efficiency optimization in IRS-assisted UAV networks,” _IEEE Trans. Commun._ , vol. 70, no. 10, pp. 6489–6502, 2022.
  * [24] 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, 2019.
  * [25] X. Dai, B. Duo, X. Yuan, and W. Tang, “Energy-efficient UAV communications: A generalized propulsion energy consumption model,” _IEEE Wireless Commun. Lett._ , vol. 11, no. 10, pp. 2150–2154, 2022.
  * [26] Y. Cai, Z. Wei, S. Hu, C. Liu, D. W. K. Ng, and J. Yuan, “Resource allocation and 3D trajectory design for power-efficient IRS-assisted UAV-NOMA communications,” _IEEE Trans. Wireless Commun._ , vol. 21, no. 12, pp. 10 315–10 334, 2022.
  * [27] Y. Zeng and R. Zhang, “Energy-efficient UAV communication with trajectory optimization,” _IEEE Trans. Wireless Commun._ , vol. 16, no. 6, pp. 3747–3760, 2017.
  * [28] X. Xiong, C. Sun, W. Ni, and X. Wang, “Three-dimensional trajectory design for unmanned aerial vehicle-based secure and energy-efficient data collection,” _IEEE Trans. Veh. Technol._ , 2022.
  * [29] T. Stathopoulos, H. Alrawashdeh, A. Al-Quraan, B. Blocken, A. Dilimulati, M. Paraschivoiu, and P. Pilay, “Urban wind energy: Some views on potential and challenges,” _J. Wind. Eng. Ind. Aerodyn._ , vol. 179, pp. 146–157, 2018\.
  * [30] P. Abichandani, D. Lobo, G. Ford, D. Bucci, and M. Kam, “Wind measurement and simulation techniques in multi-rotor small unmanned aerial vehicles,” _IEEE Access_ , vol. 8, pp. 54 910–54 927, 2020.
  * [31] G. Nachmani, “Minimum-energy flight paths for UAVs using mesoscale wind forecasts and approximate dynamic programming,” Master’s thesis, Monterey California Naval Postgraduate School, Tech. Rep., 2007.
  * [32] Y. Zhang, J. Lyu, and L. Fu, “Energy-efficient trajectory design for UAV-aided maritime data collection in wind,” _IEEE Trans. Wireless Commun._ , vol. 21, no. 12, pp. 10 871–10 886, 2022.
  * [33] R. Calmer, G. C. Roberts, J. Preissler, K. J. Sanchez, S. Derrien, and C. O’Dowd, “Vertical wind velocity measurements using a five-hole probe with remotely piloted aircraft to study aerosol–cloud interactions,” _Atmos Meas Tech_ , vol. 11, no. 5, pp. 2583–2599, 2018.
  * [34] G. Bowden, P. Barker, V. Shestopal, and J. Twidell, “The weibull distribution function and wind power statistics,” _Wind. Eng._ , pp. 85–98, 1983.
  * [35] T.-H. Yeh and L. Wang, “A study on generator capacity for wind turbines under various tower heights and rated wind speeds using weibull distribution,” _IEEE Trans. Energy Convers._ , vol. 23, no. 2, pp. 592–602, 2008.
  * [36] K. Tar, “Some statistical characteristics of monthly average wind speed at various heights,” _Renew. Sust. Energ. Rev._ , vol. 12, no. 6, pp. 1712–1724, 2008.
  * [37] J. A. Carta, C. Bueno, and P. Ramírez, “Statistical modelling of directional wind speeds using mixtures of von mises distributions: Case study,” _Energy Convers. Manag._ , vol. 49, no. 5, pp. 897–907, 2008.
  * [38] A. R. S. Bramwell, D. Balmford, and G. Done, _Bramwell’s Helicopter Dynamics_. Elsevier, 2001.
  * [39] C. You and R. Zhang, “Hybrid offline-online design for UAV-enabled data harvesting in probabilistic los channels,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 6, pp. 3753–3768, 2020.
  * [40] X. Qian, M. Di Renzo, J. Liu, A. Kammoun, and M.-S. Alouini, “Beamforming through reconfigurable intelligent surfaces in single-user MIMO systems: SNR distribution and scaling laws in the presence of channel fading and phase noise,” _IEEE Wireless Commun. Lett._ , vol. 10, no. 1, pp. 77–81, 2020.
  * [41] A. Shapiro, “Stochastic programming approach to optimization under uncertainty,” _Math. Program._ , vol. 112, no. 1, pp. 183–220, 2008.
  * [42] K. Shen and W. Yu, “Fractional programming for communication systems—Part I: Power control and beamforming,” _IEEE Trans. Signal Processing_ , vol. 66, no. 10, pp. 2616–2630, 2018.
  * [43] 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, 2018.
  * [44] A. Ben-Tal and A. Nemirovski, _Lectures on modern convex optimization: analysis, algorithms, and engineering applications_. SIAM, 2001.
  * [45] K.-Y. Wang, A. M.-C. So, T.-H. Chang, W.-K. Ma, and C.-Y. Chi, “Outage constrained robust transmit optimization for multiuser MISO downlinks: Tractable approximations by conic optimization,” _IEEE Trans. Signal Processing_ , vol. 62, no. 21, pp. 5690–5705, 2014.
  * [46] C. You, X. Peng, and R. Zhang, “3D trajectory design for UAV-enabled data harvesting in probabilistic LoS channel,” in _2019 IEEE Global Commun. Conf. (GLOBECOM)_. IEEE, 2019, pp. 1–6.
  * [47] B. Duo, Q. Wu, X. Yuan, and R. Zhang, “Anti-jamming 3D trajectory design for UAV-enabled wireless sensor networks under probabilistic LoS channel,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 12, pp. 16 288–16 293, 2020\.