# 3D Trajectory Design for Energy-constrained Aerial CRNs Under Probabilistic
LoS Channel ††thanks: Manuscript received.

Hongjiang Lei, Xiaqiu Wu, Ki-Hong Park, and Gaofeng Pan Hongjiang Lei and
Xiaqiu Wu are with the School of Communications and Information Engineering,
Chongqing University of Posts and Telecommunications, Chongqing 400065, China
(e-mail<EMAIL_ADDRESS>cquptwxq@163.com).Ki-Hong Park is with the CEMSE
Division, King Abdullah University of Science and Technology (KAUST), Thuwal
23955-6900, Saudi Arabia (e-mail: kihong.park@kaust.edu.sa).Gaofeng Pan is
with the School of Cyberspace Science and Technology, Beijing Institute of
Technology, Beijing 100081, China (e-mail: gaofeng.pan.cn@ieee.org).

###### Abstract

Unmanned aerial vehicles (UAVs) have been attracting significant attention
because there is a high probability of line-of-sight links being obtained
between them and terrestrial nodes in high-rise urban areas. In this work, we
investigate cognitive radio networks (CRNs) by jointly designing three-
dimensional (3D) trajectory, the transmit power of the UAV, and user
scheduling. Considering the UAV’s onboard energy consumption, an optimization
problem is formulated in which the average achievable rate of the considered
system is maximized by jointly optimizing the UAV’s 3D trajectory,
transmission power, and user scheduling. Due to the non-convex optimization
problem, a lower bound on the average achievable rate is utilized to reduce
the complexity of the solution. Subsequently, the original optimization
problem is decoupled into four subproblems by using block coordinate descent,
and each subproblem is transformed into manageable convex optimization
problems by introducing slack variables and successive convex approximation.
Numerical results validate the effectiveness of our proposed algorithm and
demonstrate that the 3D trajectories of UAVs can enhance the average
achievable rate of aerial CRNs.

###### Index Terms:

Cognitive radio networks, unmanned aerial vehicle, 3D trajectory design,
probabilistic LoS channel, power control.

##  I Introduction

###  I-A Background and Related Works

Unmanned aerial vehicles (UAVs) have become a revolutionary technology with
extensive applications in various fields. UAVs were initially developed for
military purposes and are now widely used in civilian and commercial sectors
such as agriculture, disaster management, surveillance, and aerial
photography, including surveillance and monitoring, aerial imaging, precision
agriculture, intelligent logistics, law enforcement, disaster response, and
prehospital emergency [1]. The rapid development of UAV technology, including
the improvement of autonomy, payload capacity, and flight endurance, has
promoted the widespread application of UAVs in various industries worldwide
[2, 3]. The continuous development of UAV technology has brought new
opportunities and challenges, requiring further exploration and research in
areas such as regulation, security, and integration into existing
infrastructure. Modifying the positioning or mapping out the flight path of
UAVs can ensure a stable line-of-sight (LoS) connection to terrestrial nodes
(TNs) with optimal reliability. Consequently, the altitude and lateral
placement of UAVs play a crucial role in enhancing the efficiency of UAV
communication systems. The trajectory designing has emerged as a critical
challenge to address in the development of UAV-assisted communication systems
[4].

To enhance the performance of the Internet of Thing (IoT) systems, various
schemes have been proposed in different scenarios. UAVs can play different
roles in communication systems, such as serving as aerial base stations (BSs)
to transmit information to ground users [5]-[8]. In [5], a multi-purpose UAV-
enabled wireless network was considered. Specifically, user scheduling, UAV
trajectories, and transmission power were jointly optimized to maximize the
minimum average rate among all users. The results indicated that, compared to
traditional static base stations, the mobility of UAVs offers advantages in
achieving better air-to-ground (A2G) channel conditions and providing
additional flexibility for interference mitigation, thereby enhancing system
throughput. Research on secure UAV-assisted downlink transmission in the
presence of colluding eavesdroppers was presented in [6]. A single-antenna UAV
served multiple ground users with imperfect location information, utilizing a
power-splitting scheme to transmit confidential information and artificial
noise. In [7], a scenario was considered where UAVs provide services to a
group of TNs and maximize the minimum secrecy rate (SR) to ensure fairness
between TNs. The results indicated that the minimum SR was significantly
improved. A UAV-to-ground communication system operating in the presence of
multiple potential Eves was investigated in [8] with incomplete information
about the Eves’ locations. effectively.

Zeng et al. derived the closed-form expression of propulsion energy
consumption for the fixed-wing UAVs in [9]. Their results demonstrated that
UAVs must operate within specific speed and acceleration limits to conserve
energy. The propulsion energy consumption for the fixed-wing UAVs is a
function of the flying velocity and acceleration. Moreover, for level flight
with fixed altitude, the UAV’s energy consumption only depends on the velocity
and acceleration rather than its actual location. Considering the propulsion
energy constraint, the authors in [10] investigated the minimum average rate
maximization and the energy efficiency (EE) maximization problems by jointly
designing the trajectory, velocity, and acceleration of the UAV and the
transmit power of the TNs. Considering the UAV’s mobility constraints, the
energy efficiency of the UAV-enabled communication system with multiple
terrestrial jammer was maximized by designing the trajectory of the UAV in
[11]. The EE of the UAV communication systems was also considered in [12]. To
reduce the computational complexity, a new method based on receding horizon
optimization was introduced to solve the formulated problem. The authors in
[13] studied the max-min fairness problem of a system consisting of multiple
UAVs and multiple TNs. The minimum achievable rate was maximized by jointly
designing the UAVs’ trajectories, power allocation, and user scheduling.

Rotary-wing UAVs have the advantage of being able to take off and land
vertically, as well as hover in a stationary position. Zeng et al. derived the
closed-form expression of propulsion energy consumption models for the rotary-
wing UAVs in [14]. Their results have shown that the propulsion power
consumption of rotary-wing UAVs consists of three components: blade profile,
induced, and parasite power. The first two parts increase with increasing
speed, and the third part decreases with increasing speed. The authors in [15]
minimized the maximum energy consumption of the rotary-wing UAV-enabled IoT
system by jointly designing the UAV’s trajectory, user scheduling, and power
allocation. Considering the average and peak power constraints, the authors in
[16] maximized the secrecy EE (SEE) of the full-duplex UAV communication
systems by jointly designing the trajectory and the power allocation. Their
results show that the SEE gains depended on the capability of the self-
interference cancelation at the full-duplex UAV. In [17], the worst-case
average SR and secrecy energy efficiency of the dual-UAV communication system
were maximized, respectively. Their results demonstrated there is a tradeoff
between maximizing the total information bits and minimizing the total
propulsion energy consumption. The authors in [18] investigated the passive
eavesdropping scenario without the eavesdropper’s instantaneous channel state
information. Considering the constraints of connection outage probability,
secrecy outage probability, and securely collected bits, the SEE was maximized
by jointly optimizing the user scheduling and transmit power, UAV’s transmit
power, trajectory, codeword rate, and redundancy rate. All the previous works
utilized the deterministic LoS channel model which is usually a valid
assumption for the scenarios without high and dense obstacles and the two-
dimensional (2D) UAV trajectory was designed with a fixed altitude.

When the UAV flies at a relatively low altitude, the shadowing effect becomes
more significant because obstacles will occasionally block the signal
propagation between the UAV and the TN. The probabilistic LoS (PLoS) channel
model was proposed in [19], and their results demonstrated that the
probability of the LoS are functions of the elevation angle. Specifically,
there are two states for the A2G channel, LoS and non-LoS (NLoS), and the
probabilities of LoS/NLoS states depend on the relative position between the
UAV and TN and the distributions of building density and height. Intuitively,
the LoS probability increases with the elevation angle, by either moving the
UAV horizontally closer to the ground node or increasing its altitude. In
[20], a UAV was utilized to transmit jamming signals to enhance the
terrestrial communication link. The expected SR was maximized by jointly
optimizing the transmit power at the BS, the jamming power of the UAV, and the
2D trajectory of the UAV. In [21], the minimum average data collection rate
was maximized by jointly optimizing the three-dimensional (3D) trajectory, the
flying speed, and user scheduling. In [22], the authors considered an aerial
wireless sensor network with a malicious terrestrial jammer and maximized the
minimum expected rate by jointly optimizing the transmission scheduling and
the 3D trajectory. In [23], a UAV was utilized as a relay to transmit signals
to the TN; the EE was maximized by optimizing the 3D trajectory of the UAV
considering the energy consumption of the UAV. In [24], the authors considered
the fairness among the TNs and maximized the minimum expected sum throughput
of the TNs by jointly designing the 3D trajectory of the UAV and the TN
scheduling. In [25], the authors studied the problem of secure data collection
problem in aerial communication systems with multiple location-uncertain
terrestrial eavesdroppers. The bandwidth allocation and the 3D trajectory of
the UAV were jointly designed to maximize the system’s overall fair SR subject
to flight energy consumption, user fairness, and secure transmission
constraints.

Cognitive radio technology is considered an advanced technology that can solve
the problem of spectrum scarcity in UAV-assisted communication systems by
using dynamic spectrum access techniques [26]. In [27], the authors considered
the aerial underlay IoT systems with a single location-uncertainty
eavesdropper. The cognitive UAV’s 2D trajectory, transmit power, and user
scheduling were jointly designed to maximize the average SR of the cognitive
users (CUs). The authors in [28] considered the aerial IoT systems with
multiple location-uncertainty full-duplex eavesdroppers worked in colluding
mode, and the cognitive UAV’s 2D trajectory and transmit power were jointly
designed to maximize the worst average SR of the CUs. In [29], a UAV was
utilized as a friendly jammer to transmit artificial noise for interrupting
the eavesdropper in both scenarios with perfect and imperfect locations of the
eavesdropper and TN. The average SR of the CU was maximized by jointly
optimizing the transmit power and UAV’s 3D trajectory. In [30], a UAV was
utilized as a cognitive relay to forward the signals to CUs and the sum
throughput of the CUs was maximized by jointly optimizing the power allocation
and 3D trajectory. It should be noted that LoS channel model was utilized in
[29] and [30], although the 3D trajectory was designed. Based on the PLoS
channel model, the authors in [31] considered both the quasi-stationary UAV
scenario and mobile UAV scenario. The CU’s achievable rate was maximized by
jointly designing the UAV’s 3D trajectory and power control, subject to the
UAV’s altitude and power constraints and the interference temperature (IT)
constraint for the primary user (PU). In [32], the SR of the UAV-enabled
cognitive radio network (CRN) was maximized by optimizing UAV’s 3D trajectory,
velocity, and acceleration considering the requirements location constraint,
the speed constraint of UAV, and the IT constraint for all the PU.

TABLE I:  Related Works to Trajectory Design of UAV Communication Systems. Reference | Role of UAVs |  LoS/PLoS |  2D/3D trajectory |  Propulsion Limitation | Rotary/Fixed |  Power Control |  Transmission Scheduling | CRNs
---|---|---|---|---|---|---|---|---
[5] | Transmitter |  LoS | 2D |  |  | $\checkmark$ | $\checkmark$ |
[6] | Transmitter |  LoS | 2D |  |  | $\checkmark$ | $\checkmark$ |
[7] | Transmitter |  LoS | 2D |  |  | $\checkmark$ | $\checkmark$ |
[8] | Transmitter |  LoS | 2D |  |  | $\checkmark$ | $\checkmark$ |
[9] | Receiver |  LoS | 2D | $\checkmark$ | Fixed |  |  |
[10] | Receiver |  LoS | 2D | $\checkmark$ | Fixed | $\checkmark$ | $\checkmark$ |
[11] | Receiver |  LoS | 2D | $\checkmark$ | Fixed |  |  |
[12] | Receiver |  LoS | 2D | $\checkmark$ | Fixed |  | $\checkmark$ |
[13] | Transmitter |  LoS | 2D | $\checkmark$ | Fixed | $\checkmark$ | $\checkmark$ |
[14] | Receiver |  LoS | 2D | $\checkmark$ | Rotary |  | $\checkmark$ |
[15] | Receiver |  LoS | 2D | $\checkmark$ | Rotary | $\checkmark$ | $\checkmark$ |
[16] | Receiver |  LoS | 2D | $\checkmark$ | Rotary | $\checkmark$ |  |
[17] | Receiver and Jammer |  LoS | 2D | $\checkmark$ | Rotary | $\checkmark$ | $\checkmark$ |
[18] | Receiver |  LoS | 2D | $\checkmark$ | Rotary | $\checkmark$ | $\checkmark$ |
[20] | Jammer |  PLoS | 2D |  |  | $\checkmark$ |  |
[21] | Receiver |  PLoS | 3D |  |  |  | $\checkmark$ |
[22] | Receiver |  PLoS | 3D |  |  |  | $\checkmark$ |
[23] | Relay |  PLoS | 3D |  |  |  |  |
[24] | Receiver |  PLoS | 3D | $\checkmark$ | Rotary |  | $\checkmark$ |
[25] | Transmitter |  PLoS | 3D | $\checkmark$ | Rotary | $\checkmark$ | $\checkmark$ |
[27] | Receiver and Jammer |  LoS | 2D |  |  | $\checkmark$ | $\checkmark$ | $\checkmark$
[28] | Transmitter |  LoS | 2D |  |  | $\checkmark$ | $\checkmark$ | $\checkmark$
[29] | Jammer |  LoS | 3D |  |  | $\checkmark$ |  | $\checkmark$
[30] | Relay |  LoS | 3D |  |  | $\checkmark$ | $\checkmark$ | $\checkmark$
[31] | Cognitive Transmitter |  LoS | 3D |  |  | $\checkmark$ | $\checkmark$ | $\checkmark$
[32] | Cognitive Transmitter |  PLoS | 3D |  |  | $\checkmark$ | $\checkmark$ | $\checkmark$
Our work | Transmitter |  PLoS | 3D | $\checkmark$ | Rotary | $\checkmark$ | $\checkmark$ | $\checkmark$
|  |  |  |  |  |  |  |

###  I-B Motivation and Contributions

The discussed works proved that the performance of aerial CRNs was
significantly improved by designing the UAV’s trajectory and transmit power.
However, these outstanding works on aerial CRNs did not consider the
propulsion energy constraint and PLoS channel model in designing the cognitive
UAV’s trajectory. In this work, based on the PLoS channel model, we consider
the limited onboard energy and jointly design the cognitive UAV’s 3D
trajectory and transmission power and user scheduling to maximize the
performance of the underlay IoT systems. The main contributions of this work
are summarized as follows:

  1. 1.

We consider an underlay aerial IoT system with multiple terrestrial cognitive
users and a primary user on the ground. Based on the PLoS channel model,
considering the UAV’s onboard energy consumption, the average achievable rate
of the considered system is maximized by jointly optimizing the UAV’s 3D
trajectory, transmission power, and user scheduling. The lower bound on the
average achievable rate is utilized to reduce the complexity of the solution
and the original optimization problem is decoupled into four subproblems by
using the block coordinate descent (BCD). Each subproblem is transformed into
manageable convex optimization problems by introducing slack variables and the
successive convex approximation (SCA).

  2. 2.

The simulation results of the proposed scheme are compared with benchmark
schemes, demonstrating that the proposed trajectory scheme effectively
improves the average achievable rate of CRN IoT systems. Moreover, in CRNs, IT
threshold, flight altitude, and UAV energy consumption are very crucial for
trajectory design.

  3. 3.

Although the 3D trajectory design of UAVs was investigated in some outstanding
works, such as [21] and [22], the energy consumption was ignored. Considering
the energy consumption, IT threshold, and transmitting power over 3D
trajectory design in this work makes the formulated optimization problem more
challenging to solve.

  4. 4.

Relative to [29]-[31], wherein the 3D trajectory design in the CRNs was
investigated based on the LoS channel model, the vertical position of the UAV
only depends on the IT constraint. In the CRNs based on the PLoS channel
model, optimizing the vertical position of the UAV must consider both the IT
constraints and the loss probability. Although the 3D trajectory design based
on the PLoS model was considered in [32], the elevation angle was assumed to
be fixed to simplify the trajectory design. This work considers both the PLoS
channel model and the energy consumption of the rotary-wing UAV in designing
the 3D trajectory.

##  II System Model

Figure 1: UAV-enabled CRN consists of an aerial base station (${B}$), multiple
cognitive users (${U_{r}}$), and a primary user (${D}$). The communication
between UAV and TNs has corresponding LoS and NLoS probabilities based on
their elevation angles.

As shown in Fig. 1, we consider an underlay IoT system where an aerial base
station $\left(B\right)$ transmits information to ground users. There are
multiple cognitive users $\left(U_{r},r=1,\cdots,R\right)$ and a primary user
$\left(D\right)$ on the ground. All the nodes are equipped with a single
antenna. The flight period of the UAV is $T$, and the period is decomposed
into $N$ time slots. The entire time is discredited as
$n\in\mathcal{N}\triangleq\\{1,\ldots,N\\}$, where the time interval for each
time slot is $\delta_{t}=\frac{T}{N}$. When $\delta_{t}$ is small enough, the
position of the UAV can be considered stationary at each time slot [5, 27]. In
this work, a 3D Cartesian coordinate system is utilized and the horizontal
positions of $U_{r}$ and $D$ are expressed as
$\mathbf{w}_{U_{r}}=\left[x_{r},y_{r}\right]^{T}\in\mathbb{R}^{2\times 1}$ and
$\mathbf{w}_{D}=\left[x_{D},y_{D}\right]^{T}\in\mathbb{R}^{2\times 1}$,
respectively.

The PLoS model is utilized and the LoS probability between $B$ and the ground
node $X\left({X\in\left\\{{{U_{r}},D}\right\\}}\right)$ in the $n$th time slot
is expressed as
$P_{BX}^{\mathrm{L}}\left[n\right]={\left({1+a{e^{\left({-b\left[{{\theta_{BX}}\left[n\right]-a}\right]}\right)}}}\right)^{-1}}$,
where $a>0$ and $b>0$ are constants specified by the actual environment and
the elevation angles between $B$ and $X$ [19] are expressed as

$\displaystyle{\theta_{B{U_{r}}}}\left[n\right]$
$\displaystyle=\frac{{180}}{\pi}\arctan\left({\frac{{{z_{B}}\left[{n}\right]}}{{\parallel{{\mathbf{q}}_{B}}\left[n\right]-{{\mathbf{w}}_{{U_{r}}}}\parallel}}}\right),\forall
n,r,$ (1a) $\displaystyle{\theta_{{BD}}}\left[n\right]$
$\displaystyle=\frac{{180}}{\pi}\arctan\left({\frac{{{z_{B}}\left[{n}\right]}}{{\parallel{{\mathbf{q}}_{B}}\left[n\right]-{{\mathbf{w}}_{D}}\parallel}}}\right),\forall
n,$ (1b)

where ${{\bf{q}}_{B}}$ and $z_{B}$ denote the horizontal and vertical position
of $B$, respectively.

The channel coefficient between $B$ and $X$ in the $n$th time slot is
expressed as

$\displaystyle
h_{BX}\left[n\right]=\left\\{\begin{array}[]{ll}\rho_{0}d_{BX}\left[n\right]^{-\alpha_{L}},&\mathrm{LoS},\\\
\mu\rho_{0}d_{BX}\left[n\right]^{-\alpha_{N}},&\mathrm{NLoS},\end{array}\right.$
(4)

where
$d_{BX}\left[n\right]=\sqrt{z_{B}\left[n\right]^{2}+\left\|\mathbf{q}_{B}\left[n\right]-\mathbf{w}_{X}\right\|^{2}}$,
$\rho_{0}$ denotes the channel gain with the unit reference distance in the
LoS environment, $\mu$ is the additional signal attenuation factor for the
NLoS environment, and $\alpha_{L}$ and $\alpha_{N}$ are the path loss
exponents for LoS and NLoS environments, typically satisfying
$\alpha_{L}$$\leq$$\alpha_{N}$ [24]. Similar to [22, 27], the small-scale
fading is not considered in this work. Therefore, the expected rate of $U_{r}$
in the $n$th time slot is expressed as [20]

$\displaystyle\bar{R}_{r}\left[n\right]=P_{BU_{r}}^{\mathrm{L}}\left[n\right]R_{r}^{\mathrm{L}}\left[n\right]+P_{BU_{r}}^{\mathrm{N}}\left[n\right]R_{r}^{\mathrm{N}}\left[n\right],$
(5)

where
$R_{r}^{\mathrm{L}}\left[n\right]=\log_{2}\left(1+\frac{P_{B}\left[n\right]\gamma}{d_{BU_{r}}^{\alpha_{L}}\left[n\right]}\right)$,
$R_{r}^{\mathrm{N}}\left[n\right]=\log_{2}\left(1+\frac{P_{B}\left[n\right]\gamma\mu}{d_{BU_{r}}^{\alpha_{N}}\left[n\right]}\right)$,
$\gamma=\frac{\rho_{0}}{\sigma^{2}}$, where $P_{B}\left[n\right]$ represents
the transmission power of $B$, and $\sigma^{2}$ represents the variance of the
additive white Gaussian noise (AWGN). Because the rate in the NLoS scenario is
much lower than that in the LoS scenario [24], the average rate of the system
is approximated as

$\displaystyle\bar{R}_{r}^{\mathrm{lb}}\left[n\right]=P_{BU_{r}}^{\mathrm{L}}\left[n\right]\log_{2}\left(1+\frac{P_{B}\left[n\right]\gamma}{d_{BU_{r}}^{\alpha_{L}}\left[n\right]}\right).$
(6)

As can be seen from (4) and (5), ${z_{B}}\left[n\right]$ increases results in
$d_{BX}\left[n\right]$ increases, the path loss increases and the rate
decreases. However, according to (1a) and (1b), an increase of
${z_{B}}\left[n\right]$ also leads to an increase of ${\theta_{BX}}$, thus
increasing $P_{BX}^{\mathrm{L}}\left[n\right]$. Thus, there is a tradeoff
between the LoS probability and the rate, and there is an optimal altitude to
maximize the expected rate given in (6). To illustrate this tradeoff, a simple
scenario is illustrated in Fig. 2, in which the UAV flies from (-200,0) to
(200,0) with different-vertical line trajectory shown in Fig. 2(a), and the
destination node is located at (0,0). As can be seen from Fig. 2(b), in the
beginning, the LoS probability of the UAV located at 100 m is more significant
than that at 30 m, and the probability of LoS increases when it approaches the
user until it flies over the user. Plan A has the lowest LoS probability, and
Plan B has the highest LoS probability. It can be observed from Fig. 2(c) that
Plan B has the lowest rate and Plan A has the highest rate at each slot due to
the path loss. Fig. 2(d) shows the expected rate with different schemes. We
find that the results of (5) and (6) are almost equal, and the expected value
of Plan c outperforms that of other schemes, which is consistent with the
results given in Sec. V.

((a)) The different schemes.

((b)) The LoS probability.

((c)) The rate.

((d)) The expected rate.

Figure 2: The comparison among different schemes.

The horizontal energy consumption of $B$ is expressed as [14]

$\displaystyle
P_{\mathrm{hor}}\left[n\right]=\underbrace{P_{0}\left(1+\frac{3\left\|\mathbf{v}_{xy}\left[n\right]\right\|^{2}}{U_{\mathrm{tip}}^{2}}\right)}_{\mathrm{bladeprofile}}+\underbrace{\frac{1}{2}d_{0}\rho
sA\left\|\mathbf{v}_{xy}\left[n\right]\right\|^{3}}_{\mathrm{parasite}}$
$\displaystyle+\underbrace{P_{1}\left(\sqrt{1+\frac{\left\|\mathbf{v}_{xy}\left[n\right]\right\|^{4}}{4v_{0}^{4}}}-\frac{\left\|\mathbf{v}_{xy}\left[n\right]\right\|^{2}}{2v_{0}^{2}}\right)^{1/2}}_{\mathrm{induced}},$
(7)

where $\mathbf{v}_{xy}\left[n\right]$ denotes the flying speed of $B$ in the
hovering state, $P_{0}$ and $P_{1}$ are two constants representing the
inherent blade surface power and induced power, respectively,
$U_{\mathrm{tip}}$ denotes the tip speed of the rotor blade, $d_{0}$ and
$\rho$ represent the body resistance ratio and air density, and the average
rotor induced velocity is represented as $v_{0}$, $s$ and $A$ represent the
rotor stiffness and rotor disc area, respectively.

The energy consumption of $B$ in the vertical direction is as expressed as
[24, 39]

$\displaystyle
P_{\mathrm{ver}}\left[n\right]=Wv_{z}\left[n\right],\quad\forall
v_{z}\left[n\right]>0,$ (8)

where $W$ represents the weight of $B$. UAVs consume no power during vertical
descent [39]. Therefore, when $v_{z}\left[n\right]<0$, we have
$P_{\mathrm{ver}}\left[n\right]=0$.

A binary variable $\alpha_{r}\left[n\right]\in\\{0,1\\}$ is utilized to
represent whether $B$ sends information to $U_{r}$ in the $n$th time slot.
Specifically, $\alpha_{r}\left[n\right]=1$ denotes that $U_{r}$ is scheduled,
otherwise, $\alpha_{r}\left[n\right]=0$. It is assumed that $B$ can only
communicate with one user in each time slot, then we have the following
constraint

$\displaystyle\sum\limits_{r=1}^{R}\alpha_{r}\left[n\right]\leq 1,\forall n,$
(9a) $\displaystyle\alpha_{r}\left[n\right]\in\\{0,1\\},\forall r,n.$ (9b)

The average rate of the considered system is expressed as

$\displaystyle{\bar{\Phi}}=\frac{1}{N}\sum\limits_{n=1}^{N}{\Phi\left[n\right]},$
(10)

where
$\Phi\left[n\right]=\sum\limits_{r=1}^{R}{{\alpha_{r}}\left[n\right]\bar{R}_{r}^{{\rm{lb}}}\left[n\right]}$
signifies the sum rate of the considered system in each slot.

In underlay CRNs, the expected interference power constraint on $D$ is
expressed as

$\displaystyle
P_{BD}^{\mathrm{L}}\left[n\right]\frac{P_{B}\left[n\right]\gamma}{d_{BD}^{\alpha_{L}}\left[n\right]}+P_{BD}^{\mathrm{N}}\left[n\right]\frac{P_{B}\left[n\right]\gamma\mu}{d_{BD}^{\alpha_{N}}\left[n\right]}\leq\Gamma,$
(11)

where $\Gamma$ denotes the peak interference power threshold value at $D$.

##  III Problem Formulation

In this work, the average rate of the system is optimized, which is related to
user scheduling, the transmission power and 3D trajectory, the horizontal and
vertical velocities of $B$. Then the following optimization problem is
formulated

$\displaystyle\mathcal{P}_{0}:$
$\displaystyle\max_{\mathbf{A},\mathbf{P},\mathbf{\Theta},\mathbf{Q},\mathbf{H}}{\bar{\Phi}}$
(12a) $\displaystyle\mathrm{s.t.}$ $\displaystyle 0\leq
P_{B}\left[n\right]\leq P_{B}^{\max},\forall n,$ (12b) $\displaystyle
0\leq\frac{1}{N}\sum\limits_{n=1}^{N}P_{B}\left[n\right]\leq{P_{B}^{\rm{ave}}},$
(12c)
$\displaystyle\frac{1}{N}\sum\limits_{n=1}^{N}P_{\mathrm{hor}}\left[n\right]\leq{P_{{\rm{hor}}}^{{\rm{ave}}}},\forall
n,$ (12d)
$\displaystyle\frac{1}{N}\sum\limits_{n=1}^{N}P_{\mathrm{ver}}\left[n\right]\leq{P_{{\rm{ver}}}^{{\rm{ave}}}},\forall
n,$ (12e)
$\displaystyle{\mathbf{q}}_{B}\left[N\right]={\mathbf{q}}_{B}\left[1\right],$
(12f)
$\displaystyle\left\|\mathbf{q}_{B}\left[{n+1}\right]-\mathbf{q}_{B}\left[n\right]\right\|=\delta_{t}v_{xy}\left[n\right],\forall
n,$ (12g) $\displaystyle\left\|\mathbf{v}_{xy}\left[n\right]\right\|\leq
V_{\max},\forall n,$ (12h)
$\displaystyle\left\|\mathbf{v}_{xy}\left[{n+1}\right]-\mathbf{v}_{xy}\left[n\right]\right\|\leq\delta_{t}a_{\max},\forall
n,$ (12i) $\displaystyle{z_{B}}\left[N\right]={z_{B}}\left[1\right],$ (12j)
$\displaystyle\left\|z_{B}\left[{n+1}\right]-z_{B}\left[n\right]\right\|=\delta_{t}v_{z}\left[n\right],\forall
n,$ (12k)
$\displaystyle\left\|{v}_{z}\left[n\right]\right\|\leq\hat{V}_{\max},\forall
n,$ (12l) $\displaystyle H_{\min}\leq z_{B}\left[n\right]\leq H_{\max},\forall
n,$ (12m)
$\displaystyle(\rm{\ref{eq1athetaBU}}),(\rm{\ref{eq1bthetaBD}}),(\rm{\ref{diaodu1}}),(\rm{\ref{diaodu2}}),(\rm{\ref{underlayCon}}).$

where $\mathbf{A}=\left\\{\alpha_{r}\left[n\right],\forall r,n\right\\}$,
$\mathbf{P}=\left\\{P_{B}\left[n\right],\forall n\right\\}$,
$\mathbf{Q}=\\{\mathbf{q}_{B}\left[n\right],v_{xy}\left[n\right],\forall
n\\}$, $\mathbf{H}=\\{z_{B}\left[n\right],v_{z}\left[n\right],\forall n\\}$,
$\mathbf{\Theta}=\left\\{\theta_{BU_{r}}\left[{n}\right],\theta_{BD}\left[{n}\right],\forall
n,r\right\\}$, $V^{max}$ and $\hat{V}^{max}$ represent the maximum horizontal
and vertical velocity, respectively, $a_{\max}$ and $\hat{a}_{\max}$ represent
the maximum horizontal and vertical acceleration, respectively, $H_{\min}$ and
$H_{\max}$ represent the minimum and maximum vertical altitude of $B$,
respectively, and $P_{B}^{\max}$ and $P_{B}^{\rm{ave}}$ represent the maximum
and average transmission power of $B$, respectively. (12d) and (12e) represent
the propulsion energy constraint for $B$, (12f) and (12j) are the constraints
on the starting and final positions of $B$, that is, $B$ must fly back to the
starting position at the end of each flight cycle, (12g)-(12i) and (12k)-(12m)
represent the velocity and acceleration constraints, respectively, (12m)
denotes the flight altitude constraint of $B$, (1a) and (1b) are the
constraints on elevation angle between $B$ and TN, and (9a) and (9b) signify
the user scheduling constraints, and (11) denotes the interference power
constraint of $B$ on $D$,

Some factors make the optimization problem $\mathcal{P}_{0}$ difficult to
solve. Firstly, the average rate expression in (6) is related to the UAV’s 3D
trajectory and also to the probability of LoS/NLoS communication. Therefore,
the objective function (12a) in optimization problem $\mathcal{P}_{0}$ is
highly coupled. The second constraint brought by user scheduling (9b) is
binary integer constraint, which is non-convex. Also, consider that the
elevation constraints (1a) and (1b) caused by PLoS modeling are nonlinear
constraints. Therefore, the optimization problem $\mathcal{P}_{0}$ is a
nonlinear mixed integer non-convex optimization problem, which is usually
troublesome to solve.

##  IV Proposed Algorithm for Problem $\mathcal{P}_{0}$

To solve $\mathcal{P}_{0}$, we utilize the BCD technology to decompose the
original problem into multiple subproblems. Specifically, for the given other
variables, $\mathbf{A}$, $\mathbf{P}$, $\mathbf{H}$, and $\mathbf{Q}$ are
optimized in each subproblem respectively. In addition, the SCA technology is
utilized to transform the non-convex constraints into convex constraints.

###  IV-A Subproblem 1: Optimizing User Scheduling Variable

In this subsection, the user scheduling variable $\mathbf{A}$ is optimized
with given {$\mathbf{P}$,$\mathbf{H}$,$\mathbf{Q}$,$\mathbf{\Theta}$}. The
binary variable constraint $\alpha_{r}\left[n\right]\in\\{0,1\\}$ is relaxed
as a linear variable $0\leq\alpha_{r}\left[n\right]\leq 1$ and
$\mathcal{P}_{0}$ is reformulated as

$\displaystyle\mathcal{P}_{1.1}:$ $\displaystyle\max_{\mathbf{A},\eta}\eta$
(13a) $\displaystyle\mathrm{s.t.}$ $\displaystyle{\bar{\Phi}}\geq\eta,$ (13b)
$\displaystyle 0\leq\alpha_{r}\left[n\right]\leq 1,\forall n,r,$ (13c)
$\displaystyle(\rm{\ref{diaodu1}}).$

where $\eta$ is a slack variable. $\mathcal{P}_{1.1}$ is a strictly linear
optimization problem, which can be effectively solved by using the CVX toolbox
in [38].

###  IV-B Subproblem 2: Optimizing Transmit Power of $B$

For given {$\mathbf{A}$,$\mathbf{Q}$,$\mathbf{H}$,$\mathbf{\Theta}$}, we have
an optimization problem as follows

$\displaystyle\mathcal{P}_{2.1}:$ $\displaystyle\max_{\mathbf{P},\eta}\eta$
(14a) $\displaystyle\mathrm{s.t.}$ $\displaystyle{\bar{\Phi}}\geq\eta,$ (14b)
$\displaystyle(\rm{\ref{underlayCon}}),(\rm{\ref{eq9b}}),(\rm{\ref{eq9c}}).$

The objective function in $\mathcal{P}_{2.1}$ is linear, the left-hand side
(LHS) of (14b) is a concave function in the form of $\log_{2}(1+x)$, which is
a convex constraint, and (11), (12b), and (12c) are a linear constraints, so
$\mathcal{P}_{2.1}$ is a convex optimization problem and can be solved through
the CVX toolbox.

###  IV-C Subproblem 3: Optimizing Horizontal Trajectory and Velocity of $B$

In this subsection, the horizontal trajectory and velocity of $B$ is optimized
for provided {$\mathbf{A}$,$\mathbf{P}$,$\mathbf{H}$}. The original
optimization problem is rewritten as

$\displaystyle\mathcal{P}_{3.1}:$
$\displaystyle\max_{\mathbf{Q},\mathbf{\Theta},\eta}\eta$ (15a)
$\displaystyle\mathrm{s.t.}$ $\displaystyle{\bar{\Phi}}\geq\eta,$ (15b)
$\displaystyle(\rm{\ref{underlayCon}}),(\rm{\ref{eq9d}}),(\rm{\ref{eq9e}})-(\rm{\ref{eq9h}}),$
$\displaystyle(\rm{\ref{eq1athetaBU}}),(\rm{\ref{eq1bthetaBD}}).$

$\mathcal{P}_{3.1}$ is non-convex because (12d) and (11) are non-convex
constraints, (1a) and (1b) are nonlinear constraints, and
$\bar{R}_{r}^{\mathrm{lb}}\left[n\right]$ in constraint (15b) not only depends
on $\mathbf{Q}$, but also on $P_{BU_{r}}^{\mathrm{L}}\left[n\right]$.

$\displaystyle\frac{1}{N}\sum\limits_{n=1}^{N}{\left[{{P_{0}}\left({1+\frac{{3{{\left\|{{{\bf{v}}_{xy}}\left[n\right]}\right\|}^{2}}}}{{U_{{\rm{tip}}}^{2}}}}\right)+\frac{1}{2}{d_{0}}\rho
sA{{\left\|{{{\bf{v}}_{xy}}\left[n\right]}\right\|}^{3}}+{P_{1}}{\lambda_{B}}\left[n\right]}\right]}\leq
P_{{\rm{hor}}}^{{\rm{ave}}}$ (16a)
$\displaystyle\frac{1}{\lambda_{{B}}\left[n\right]^{2}}\leq\lambda_{{B}}\left[n\right]^{2}+\frac{\left\|\mathbf{v}_{xy}\left[n\right]\right\|^{2}}{v_{0}^{2}}$
(16b)

By introducing a relaxation variable $\left\\{\lambda_{{B}}\left[n\right]\geq
0\right\\}$, (12d) is rewritten as (16a) and (16b), shown at the top of this
page, It should be noted that the right-hand side (RHS) of (16b) is a convex
function for both $\lambda_{{B}}\left[n\right]$ and
$\left\|\mathbf{v}_{xy}\left[n\right]\right\|$, which means it is not a convex
constraint. By utilizing SCA and replacing the RHS with its convex lower
bound, (16b) is rewritten as

$\displaystyle\frac{1}{\lambda_{{B}}\left[n\right]^{2}}\leq{\lambda}^{{\left(\kappa\right)}}_{{B}}\left[n\right]^{2}+2{\lambda}^{{\left(\kappa\right)}}_{{B}}\left[n\right]\left(\lambda_{{B}}\left[n\right]-{\lambda}^{{\left(\kappa\right)}}_{{B}}\left[n\right]\right)$
$\displaystyle+\frac{\left\|{\mathbf{v}}^{{\left(\kappa\right)}}_{xy}\left[n\right]\right\|^{2}}{v_{0}^{2}}+\frac{2}{v_{0}^{2}}\left({\bf{v}}^{{\left(\kappa\right)}}_{xy}\left[n\right]\right)^{T}\left({\bf{v}}_{xy}\left[n\right]-{{\bf{v}}}^{{\left(\kappa\right)}}_{xy}\left[n\right]\right),$
(15)

where ${\bf{v}}_{xy}^{{\left(\kappa\right)}}\left[n\right]$ represents the
horizontal velocity obtained in the $k$th iteration.

Similarly, (11) is reformulated as

$\displaystyle\Gamma$
$\displaystyle\geq\frac{P_{B}\left[n\right]\gamma}{p_{\mathrm{L}}\left[n\right]t_{D,{\rm
L}}\left[n\right]}+\frac{P_{B}\left[n\right]\gamma\mu}{p_{\mathrm{N}}\left[n\right]t_{D,{\rm
N}}\left[n\right]},$ (16a) $\displaystyle p_{\mathrm{L}}\left[n\right]$
$\displaystyle\leq 1+ae^{ab}e^{-b\theta_{BD}\left[n\right]},$ (16b)
$\displaystyle\theta_{BD}\left[n\right]$
$\displaystyle\geq\frac{180}{\pi}\arctan\left(\frac{z_{B}\left[n\right]}{\left\|\mathbf{q}_{B}\left[n\right]-\mathbf{w}_{D}\right\|}\right),$
(16c) $\displaystyle p_{\mathrm{N}}\left[n\right]$ $\displaystyle\leq
1+\frac{1}{a}e^{-ab}e^{b\phi_{BD}\left[n\right]},$ (16d)
$\displaystyle\phi_{BD}\left[n\right]$
$\displaystyle\leq\frac{180}{\pi}\arctan\left(\frac{z_{B}\left[n\right]}{\left\|\mathbf{q}_{B}\left[n\right]-\mathbf{w}_{D}\right\|}\right),$
(16e) $\displaystyle t_{D,h}\left[n\right]$
$\displaystyle\leq\left(\left\|\mathbf{q}_{B}\left[n\right]-\mathbf{w}_{D}\right\|^{2}+z_{B}^{2}\left[n\right]\right)^{\alpha_{h}/2},$
(16f)

where $h\in\\{\mathrm{L},\mathrm{N}\\}$,
$\mathbf{p}_{h}=\left\\{p_{h}\left[n\right],\forall n\right\\}$, and
$\mathbf{t}_{D,h}=\left\\{t_{D,h}\left[n\right],\forall n,r\right\\}$ are
relaxation variables, (16c) and (16e) relax (1b) into inequality constraints.
It should be noted that the RHS of (16b), (16d), (16e) and (16f) are convex
functions, which make them non-convex constraints. Then, the SCA technology is
utilized to convert the terms $e^{-b\theta_{BD}\left[n\right]}$,
$e^{b\phi_{BD}\left[n\right]}$,
$\arctan(\frac{z_{B}\left[n\right]}{\|\mathbf{q}_{B}\left[n\right]-\mathbf{w}_{D}\|})$
and
$(\|\mathbf{q}_{B}\left[n\right]-\mathbf{w}_{D}\|^{2}+z_{B}^{2}\left[n\right])^{\alpha_{h}/2}$
to linear functions. By replacing them with their respective convex lower
bounds, we obtain

$\displaystyle e^{-b\theta_{BD}\left[n\right]}$ $\displaystyle\geq
F_{1}^{{\left(\kappa\right)}}\left[n\right],$ (17a) $\displaystyle
e^{b\phi_{BD}\left[n\right]}$ $\displaystyle\geq
F_{2}^{{\left(\kappa\right)}}\left[n\right],$ (17b)
$\displaystyle\left(\left\|\mathbf{q}_{B}\left[n\right]-\mathbf{w}_{D}\right\|^{2}+z_{B}^{2}\left[n\right]\right)^{\alpha_{h}/2}$
$\displaystyle\geq F_{3,h}^{{\left(\kappa\right)}}\left[n\right],$ (17c)
$\displaystyle\arctan\left(\frac{z_{B}\left[n\right]}{\left\|\mathbf{q}_{B}\left[n\right]-\mathbf{w}_{D}\right\|}\right)$
$\displaystyle\geq F_{5}^{{\left(\kappa\right)}}\left[n\right],$ (17d)

where
$F_{1}^{\left(\kappa\right)}\left[n\right]={e^{-b\theta_{BD}^{\left(\kappa\right)}\left[n\right]}}-b{e^{-b\theta_{BD}^{\left(\kappa\right)}\left[n\right]}}\left({{\theta_{BD}}\left[n\right]-\theta_{BD}^{\left(\kappa\right)}\left[n\right]}\right)$,
$F_{2}^{\left(\kappa\right)}\left[n\right]={e^{b\phi_{BD}^{\left(\kappa\right)}\left[n\right]}}+b{e^{-b\phi_{BD}^{\left(\kappa\right)}\left[n\right]}}\left({{\phi_{BD}}\left[n\right]-\phi_{BD}^{\left(\kappa\right)}\left[n\right]}\right)$,
$F_{3,h}^{\left(\kappa\right)}\left[n\right]={\left({d_{{B_{D}}}^{\left(\kappa\right)}\left[n\right]}\right)^{{\alpha_{h}}}}+{\alpha_{h}}{\left({d_{{B_{D}}}^{\left(\kappa\right)}\left[n\right]}\right)^{{\alpha_{h}}-2}}+{\left({{\mathbf{q}}_{B}^{\left(\kappa\right)}\left[n\right]-{{\mathbf{w}}_{D}}}\right)^{T}}\left({{{\mathbf{q}}_{B}}\left[n\right]-{\mathbf{q}}_{B}^{\left(\kappa\right)}\left[n\right]}\right)$,
and
$F_{5}^{\left(\kappa\right)}\left[n\right]=\arctan\left({\frac{{{z_{B}}\left[n\right]}}{{\left\|{{\mathbf{q}}_{B}^{\left(\kappa\right)}\left[n\right]-{{\mathbf{w}}_{D}}}\right\|}}}\right)-\frac{{{z_{B}}\left[n\right]\left({\left\|{{{\mathbf{q}}_{B}}\left[n\right]-{{\mathbf{w}}_{D}}}\right\|-\left\|{{\mathbf{q}}_{B}^{{\left(\kappa\right)}}\left[n\right]-{{\mathbf{w}}_{D}}}\right\|}\right)}}{{{{\left\|{{\mathbf{q}}_{B}^{\left(\kappa\right)}\left[n\right]-{{\mathbf{w}}_{D}}}\right\|}^{2}}+z_{B}^{2}\left[n\right]}}$,
${\theta_{BD}^{\left(\kappa\right)}\left[n\right]}$,
${\phi_{BD}^{\left(\kappa\right)}\left[n\right]}$,
${d_{{B_{D}}}^{\left(\kappa\right)}\left[n\right]}$, and
$\mathbf{q}_{B}^{{\left(\kappa\right)}}\left[n\right]$ denote
${\theta_{BD}\left[n\right]}$, ${\phi_{BD}\left[n\right]}$,
${d_{{B_{D}}}\left[n\right]}$, and $\mathbf{q}_{B}\left[n\right]$ in the $k$th
iteration, respectively. Therefore, (16b), (16d), (16e), and (16f) are
reformulated as

$\displaystyle p_{\mathrm{L}}\left[n\right]$ $\displaystyle\leq
1+ae^{ab}F_{1}^{{\left(\kappa\right)}}\left[n\right],$ (18a) $\displaystyle
p_{\mathrm{N}}\left[n\right]$ $\displaystyle\leq
1+\frac{1}{a}e^{-ab}F_{2}^{{\left(\kappa\right)}}\left[n\right],$ (18b)
$\displaystyle\phi_{BD}\left[n\right]$
$\displaystyle\leq\frac{180}{\pi}F_{5}^{{\left(\kappa\right)}}\left[n\right],$
(18c) $\displaystyle t_{D,h}\left[n\right]$ $\displaystyle\leq
F_{3,h}^{{\left(\kappa\right)}}\left[n\right].$ (18d)

By introducing slack variables
$\mathbf{x}_{L}=\left\\{x_{L}\left[n\right],\forall n\right\\}$ and
$\mathbf{t}_{r,L}=\left\\{t_{r,L}\left[n\right],\forall n,r\right\\}$, (15b)
is reformulated as

$\displaystyle{\bar{\Phi}_{1}}\geq\eta,$ (19a) $\displaystyle
x_{L}\left[n\right]\geq 1+ae^{ab}e^{-b\theta_{BU_{r}}\left[n\right]},$ (19b)
$\displaystyle\theta_{BU_{r}}\left[n\right]\leq\frac{180}{\pi}\arctan\left(\frac{z_{B}\left[n\right]}{\left\|\mathbf{q}_{B}\left[n\right]-\mathbf{w}_{U_{r}}\right\|}\right),$
(19c) $\displaystyle
t_{r,L}\left[n\right]\geq\left(\left\|\mathbf{q}_{B}\left[n\right]-\mathbf{w}_{U_{r}}\right\|^{2}+z_{B}^{2}\left[n\right]\right)^{\alpha_{L}/2},$
(19d)

where
${\bar{\Phi}_{1}}=\frac{1}{N}\sum\limits_{n=1}^{N}\sum\limits_{r=1}^{R}\alpha_{r}\left[n\right]\bar{R}_{r,1}\left[n\right]$,
$\bar{R}_{r,1}\left[n\right]=\frac{1}{x_{L}\left[n\right]}\log_{2}\left(1+\frac{P_{B}\left[n\right]\gamma}{t_{r,L}\left[n\right]}\right)$.
It should be noted that (19c) relaxes (1a) into an inequality constraint. With
the same method as (16b), we have

$\displaystyle\theta_{BU_{r}}\left[n\right]\leq\frac{180}{\pi}F_{4}^{{\left(\kappa\right)}}\left[n\right],$
(20)

where
$F_{4}^{\left(\kappa\right)}\left[n\right]=\arctan\left({\frac{{{z_{B}}\left[n\right]}}{{\left\|{{\mathbf{q}}_{B}^{\left(\kappa\right)}\left[n\right]-{{\mathbf{w}}_{{U_{r}}}}}\right\|}}}\right)-\frac{{{z_{B}}\left[n\right]\left({\left\|{{{\mathbf{q}}_{B}}\left[n\right]-{{\mathbf{w}}_{{U_{r}}}}}\right\|-\left\|{{\mathbf{q}}_{B}^{\left(\kappa\right)}\left[n\right]-{{\mathbf{w}}_{{U_{r}}}}}\right\|}\right)}}{{{{\left\|{{\mathbf{q}}_{B}^{\left(\kappa\right)}\left[n\right]-{{\mathbf{w}}_{{U_{r}}}}}\right\|}^{2}}+z_{B}^{2}\left[n\right]}}$.

To address the non-convexity in (19a), Lemma 1 is introduced.

_Lemma 1_ : Given any $A\geq 0$, the function $f\left(x,y\right)=$
$\frac{1}{x}\log_{2}\left(1+\frac{A}{y}\right)$ is jointly convex with respect
to $x>0$ and $y>0$.

_Proof_ : See the Appendix A.

Based on Lemma 1, the first-order Taylor expansion of
$\bar{R}_{r,1}\left[n\right]$ is obtained as

$\displaystyle\bar{R}_{r,1}^{\mathrm{lb}}\left[n\right]$
$\displaystyle=A_{r}^{{\left(\kappa\right)}}\left[n\right]+B_{r,L}^{{\left(\kappa\right)}}\left[n\right]\left(t_{r,L}\left[n\right]-t_{r,L}^{{\left(\kappa\right)}}\left[n\right]\right)$
$\displaystyle\,\,\,\,\,\,\,+C_{r,L}^{{\left(\kappa\right)}}\left[n\right]\left(x_{L}\left[n\right]-x_{L}^{{\left(\kappa\right)}}\left[n\right]\right),$
(21)

where
$A_{r}^{{\left(\kappa\right)}}\left[{n}\right]=\frac{{{{\log}_{2}}\left({1+\frac{{{P_{B}}\left[{n}\right]\gamma}}{{t_{r,{\mathrm{L}}}^{{\left(\kappa\right)}}\left[{n}\right]}}}\right)}}{{x_{\mathrm{L}}^{{\left(\kappa\right)}}\left[{n}\right]}}$,
$B_{r,L}^{{\left(\kappa\right)}}\left[n\right]=\frac{-P_{B}\left[n\right]\gamma}{\ln(2)x_{L}^{{\left(\kappa\right)}}\left[n\right]\left(t_{r,L}^{{\left(\kappa\right)}}\left[n\right]\right)^{2}\left(1+\frac{P_{B}\left[n\right]\gamma}{t_{r,L}^{{\left(\kappa\right)}}\left[n\right]}\right)}$,
$C_{r,L}^{{\left(\kappa\right)}}\left[n\right]=-\frac{\log_{2}\left(1+\frac{P_{B}\left[n\right]\gamma}{t_{r,L}^{{\left(\kappa\right)}}\left[n\right]}\right)}{\left(x_{L}^{{\left(\kappa\right)}}\left[n\right]\right)^{2}}$,
and $x_{L}^{{\left(\kappa\right)}}\left[n\right]$,
$t_{r,L}^{{\left(\kappa\right)}}\left[n\right]$ and
$t_{D,h}^{{\left(\kappa\right)}}\left[n\right]$ are
$x_{\mathrm{L}}\left[n\right]$, $t_{r,L}\left[n\right]$ and
$t_{D,h}\left[n\right]$ at the $k$th iteration, respectively.

Finally, $\mathcal{P}_{3.1}$ is rewritten as

$\displaystyle\mathcal{P}_{3.2}$
$\displaystyle\max_{\mathbf{Q},\boldsymbol{\Theta},\eta,\lambda_{B},{x}_{\mathrm{L}},\mathbf{p}_{h},\mathbf{t}_{r,L},\mathbf{t}_{D,h}}\eta$
(22a) $\displaystyle\mathrm{s.t.}$ $\displaystyle{\bar{\Phi}_{2}}\geq\eta,$
(22b)
$\displaystyle(\rm{\ref{eq9e}})-(\rm{\ref{eq9h}}),(\rm{\ref{shuipingslack}}),(\rm{\ref{lambdasca}}),(\rm{\ref{eq15agama2}}),(\rm{\ref{eq15cthetabd1}}),$
$\displaystyle(\rm{\ref{17aF1}})-(\rm{\ref{17dF3}}),(\rm{\ref{p32b}}),(\rm{\ref{p32d}}),(\rm{\ref{F4}}),$

where
${\bar{\Phi}_{2}}=\frac{1}{N}\sum\limits_{n=1}^{N}\sum\limits_{r=1}^{R}\alpha_{r}\left[n\right]\bar{R}_{r,1}^{\mathrm{lb}}\left[n\right]$.
$\mathcal{P}_{3.2}$ is a convex problem that can be solved using existing
optimization tools such as CVX.

###  IV-D Subproblem 4: Optimizing Horizontal Trajectory and Velocity of $B$

In this subsection, for given {$\mathbf{A}$,$\mathbf{P}$,$\mathbf{Q}$}, the
vertical trajectory $\mathbf{H}$ of $B$ is optimized. The optimization problem
is expressed as

$\displaystyle\mathcal{P}_{4.1}$
$\displaystyle:\max_{\mathbf{H},\mathbf{\Theta},\eta}\eta$ (23a)
$\displaystyle\mathrm{s.t.}$ $\displaystyle{\bar{\Phi}}\geq\eta,$ (23b)
$\displaystyle(\rm{\ref{underlayCon}}),(\rm{\ref{eq9d2}}),(\rm{\ref{eq9i}})-(\rm{\ref{eq9m}}),(\rm{\ref{eq1athetaBU}}),(\rm{\ref{eq1bthetaBD}}).$

It should be noted that $\mathcal{P}_{4.1}$ has the same problems as
$\mathcal{P}_{3.1}$. Specifically, (23b) not only depends on $\mathbf{H}$, but
also on $P_{BU_{r}}^{\mathrm{L}}\left[n\right]$, (11) is non-convex
constraint, and (1a) and (1b) are nonlinear constraints.

With the same method as (15b), (23b) is reformulated as (19a)-(19d) and (1a)
and (1b) are reformulated as (16c), (16e), and (19c). With the same method in
Subproblem 3, (11) in this subsection is reformulated as (16a)-(16f) wherein
(16b) and (16d) are reformulated as (18a) and (18b), respectively.

The SCA is utilized to deal with (16c) and (16f) with respect to
$(\|\mathbf{q}_{B}\left[n\right]-\mathbf{w}_{D}\|^{2}+z_{B}^{2}\left[n\right])^{\alpha_{h}/2}$
and
$\arctan(\frac{z_{B}\left[n\right]}{\|\mathbf{q}_{B}\left[n\right]-\mathbf{w}_{D}\|})$,
respectively. Similarly, we obtain

$\displaystyle{\left({{{\left\|{{{\mathbf{q}}_{B}}\left[{n}\right]-{{\mathbf{w}}_{D}}}\right\|}^{2}}+z_{B}^{2}\left[{n}\right]}\right)^{{\alpha_{h}}/2}}\geqslant
F_{6,h}^{(\kappa)}\left[{n}\right],$ (24a)
$\displaystyle\arctan\left({\frac{{{z_{B}}\left[{n}\right]}}{{\left\|{{{\mathbf{q}}_{B}}\left[{n}\right]-{{\mathbf{w}}_{D}}}\right\|}}}\right)\geqslant
F_{7}^{(\kappa)}\left[{n}\right],$ (24b)

where
$F_{6,h}^{(\kappa)}\left[{n}\right]={\left({d_{BD}^{(\kappa)}\left[{n}\right]}\right)^{{\alpha_{h}}}}+{\alpha_{h}}{\left({d_{BD}^{(\kappa)}\left[{n}\right]}\right)^{{\alpha_{h}}-2}}{\left({z_{B}^{(\kappa)}\left[{n}\right]}\right)^{T}}\left({{z_{B}}\left[{n}\right]-z_{B}^{(\kappa)}\left[{n}\right]}\right)$,
$F_{7}^{(\kappa)}\left[{n}\right]=$
$\arctan\left({\frac{{z_{B}^{(\kappa)}\left[{n}\right]}}{{\left\|{{{\mathbf{q}}_{B}}\left[{n}\right]-{{\mathbf{w}}_{D}}}\right\|}}}\right)+\frac{{\left\|{{{\mathbf{q}}_{B}}\left[{n}\right]-{{\mathbf{w}}_{D}}}\right\|\left({{z_{B}}\left[{n}\right]-z_{B}^{(\kappa)}\left[{n}\right]}\right)}}{{{{\left\|{{{\mathbf{q}}_{B}}\left[{n}\right]-{{\mathbf{w}}_{D}}}\right\|}^{2}}+z_{B}^{(\kappa)}{{\left[{n}\right]}^{2}}}}$,
and $z_{B}^{(\kappa)}\left[n\right]$ denotes ${z_{B}}\left[n\right]$ in the
$k$th iteration. Finally, (16c) and (16f) are reformulated as

$\displaystyle t_{D,h}\left[n\right]$ $\displaystyle\leq
F_{6,h}^{{\left(\kappa\right)}}\left[n\right],$ (25a)
$\displaystyle\theta_{BD}\left[n\right]$
$\displaystyle\geq\frac{180}{\pi}F_{7}^{{\left(\kappa\right)}}\left[n\right].$
(25b)

Then $\mathcal{P}_{4.1}$ is rewritten as

$\displaystyle\mathcal{P}_{4.2}$
$\displaystyle\max_{\mathbf{H},\mathbf{\Theta},\eta,\mathbf{x}_{L},\mathbf{p}_{h},\mathbf{t}_{r,L},\mathbf{t}_{D,h}}\eta$
(26a) $\displaystyle\mathrm{s.t.}$ $\displaystyle{\bar{\Phi}_{2}}\geq\eta,$
(26b)
$\displaystyle(\rm{\ref{eq9d2}}),(\rm{\ref{eq9i}})-(\rm{\ref{eq9m}}),(\rm{\ref{eq15agama2}}),(\rm{\ref{eq15ephibd}}),(\rm{\ref{17aF1}}),$
$\displaystyle(\rm{\ref{17bF2}}),(\rm{\ref{p32b}})-(\rm{\ref{p32d}}),(\rm{\ref{F6}}),(\rm{\ref{F7}}).$

$\mathcal{P}_{4.2}$ is a convex optimization problem that can be solved using
existing optimization tools such as CVX.

###  IV-E Convergence Analysis of Algorithm 1

$\mathcal{P}_{0}$ is a non-convex optimization problem, which decouples the
variables $\mathbf{A}$, $\mathbf{P}$, $\mathbf{\Theta}$, $\mathbf{H}$, and
$\mathbf{Q}$ in the original problem. Through BCD, solving $\mathcal{P}_{0}$
converts to solving $\mathcal{P}_{1.1}$, $\mathcal{P}_{2.1}$,
$\mathcal{P}_{3.1}$ and $\mathcal{P}_{4.1}$ in turn. By using SCA and
introducing relaxation variables, $\mathcal{P}_{3.1}$ and $\mathcal{P}_{4.1}$
that have non-convex constraints are rewritten as $\mathcal{P}_{3.2}$ and
$\mathcal{P}_{4.2}$, which are all convex optimization problems. The value of
the objective function for the original problem $\mathcal{P}_{0}$ is defined
as
$R\left(\mathbf{A}^{\kappa},\mathbf{P}^{\kappa},\mathbf{\Theta}^{\kappa},\mathbf{H}^{\kappa},\mathbf{Q}^{\kappa}\right)$
at the $k$th iteration, Algorithm 1 summarizes the details of the overall
iteration for solving problem $\mathcal{P}_{0}$, where $\varepsilon$
represents the exact convergence value.

The obtained suboptimal solution of the transformed subproblem is also the
suboptimal solution of the original non-convex subproblem, and each subproblem
is solved using SCA convex transformation iteration. Finally, all suboptimal
solutions of the subproblems that satisfy the threshold $\varepsilon$
constitute the suboptimal solution of the original problem. Therefore, our
algorithm is to alternately solve the subproblem $\mathcal{P}_{1.1}$,
$\mathcal{P}_{2.1}$, $\mathcal{P}_{3.2}$ and $\mathcal{P}_{4.2}$ to obtain the
suboptimal solution of the original problem until a solution that satisfies
the threshold $\varepsilon$ is obtained.

It is worth noting that in the classic BCD, to ensure the convergence of the
algorithm, it is necessary to accurately solve and update the subproblems of
each variable block with optimality in each iteration. But when we solve
$\mathcal{P}_{3.1}$ and $\mathcal{P}_{4.1}$ , we can only optimally solve
their approximation problem $\mathcal{P}_{3.2}$ and $\mathcal{P}_{4.2}$.
Therefore, we cannot directly apply the convergence analysis of the classical
BCD, and further proof of the convergence of Algorithm 1 is needed, as shown
below.

Input: Initialization of feasible points.

_1_ while
_$R(\mathbf{A}^{\kappa},\mathbf{P}^{\kappa},\mathbf{\Theta}^{\kappa},\mathbf{H}^{\kappa},\mathbf{Q}^{\kappa})-$_

2 __
$R(\mathbf{A}^{\kappa-1},\mathbf{P}^{\kappa-1},\mathbf{\Theta}^{\kappa-1},\mathbf{H}^{\kappa-1},\mathbf{Q}^{\kappa-1})\succ\varepsilon$
do

3 Solve $\left(\mathcal{P}_{1.1}\right)$ for given
$\left\\{\mathbf{P}^{\kappa},\mathbf{H}^{\kappa},\mathbf{Q}^{\kappa},\mathbf{\Theta}^{\kappa}\right\\}$

4 and obtain the solution $\mathbf{A}^{\kappa+1}$;

5 Solve $\left(\mathcal{P}_{2.1}\right)$ for given
$\left\\{\mathbf{A}^{\kappa+1},\mathbf{H}^{\kappa},\mathbf{Q}^{\kappa},\mathbf{\Theta}^{\kappa}\right\\}$

6 and obtain the solution $\mathbf{P}^{\kappa+1}$;

7 Solve $\left(\mathcal{P}_{3.2}\right)$ for given
$\left\\{\mathbf{A}^{\kappa+1},\mathbf{H}^{\kappa},\mathbf{P}^{\kappa+1}\right\\}$
and

8 obtain the solution $\mathbf{Q}^{\kappa+1}$;

9 Solve $\left(\mathcal{P}_{4.2}\right)$ for given
$\left\\{\mathbf{A}^{\kappa+1},\mathbf{Q}^{\kappa+1},\mathbf{P}^{\kappa+1}\right\\}$
and

10 obtain the solution $\mathbf{H}^{\kappa+1}$;

11 $\kappa=\kappa+1$;

12 Compute the objective value

13
$R(\mathbf{A}^{\kappa},\mathbf{P}^{\kappa},\mathbf{\Theta}^{\kappa},\mathbf{H}^{\kappa},\mathbf{Q}^{\kappa})$.

14 end while

Output:
$R(\mathbf{A}^{\kappa},\mathbf{P}^{\kappa},\mathbf{\Theta}^{\kappa},\mathbf{H}^{\kappa},\mathbf{Q}^{\kappa})$
with $\mathbf{A}^{*}=\mathbf{A}^{\kappa}$,

$\mathbf{P}^{*}=\mathbf{P}^{\kappa}$,$\mathbf{\Theta}^{*}=\mathbf{\Theta}^{\kappa}$,
$\mathbf{H}^{*}=\mathbf{H}^{\kappa}$,$\mathbf{Q}^{*}=\mathbf{Q}^{\kappa}$.

Algorithm 1 Iterative Procedure of Problem $\mathcal{P}_{0}$

Firstly, for step 1 of our proposed Algorithm 1, we obtain the optimal
solution to user scheduling problem $\mathcal{P}_{1.1}$ under the condition of
given other feasible variables $\mathbf{P}^{\kappa}$, $\mathbf{H}^{\kappa}$,
$\mathbf{Q}^{\kappa}$ and $\mathbf{\Theta}^{\kappa}$, therefore we can obtain

$\displaystyle
R\left(\mathbf{A}^{\kappa},\mathbf{P}^{\kappa},\mathbf{\Theta}^{\kappa},\mathbf{H}^{\kappa},\mathbf{Q}^{\kappa}\right)\leq
R\left(\mathbf{A}^{\kappa+1},\mathbf{P}^{\kappa},\mathbf{\Theta}^{\kappa},\mathbf{H}^{\kappa},\mathbf{Q}^{\kappa}\right).$
(27)

Similarly, step 2 in Algorithm 1 is to solve the power optimization problem
$\mathcal{P}_{2.1}$, given other feasible variables $\mathbf{A}^{\kappa+1}$,
$\mathbf{H}^{\kappa}$, $\mathbf{Q}^{\kappa}$ and $\mathbf{\Theta}^{\kappa}$,
this problem is a convex optimization problem, so we can also obtain the
optimal solution for $\mathcal{P}_{2.1}$, so we have

$\displaystyle
R\left(\mathbf{A}^{\kappa+1},\mathbf{P}^{\kappa},\mathbf{\Theta}^{\kappa},\mathbf{H}^{\kappa},\mathbf{Q}^{\kappa}\right)\leq
R\left(\mathbf{A}^{\kappa+1},\mathbf{P}^{\kappa+1},\mathbf{\Theta}^{\kappa},\mathbf{H}^{\kappa},\mathbf{Q}^{\kappa}\right).$
(28)

Then, in step 3 of Algorithm 1, the approximate problem $\mathcal{P}_{3.2}$ of
Horizontal Trajectory Optimization Problem $\mathcal{P}_{3.2}$ is solved by
fixing other feasible variables $\mathbf{A}^{\kappa+1}$,
$\mathbf{P}^{\kappa+1}$ and $\mathbf{H}^{\kappa}$. We obtain

$\displaystyle
R\left(\mathbf{A}^{\kappa+1},\mathbf{P}^{\kappa+1},\mathbf{\Theta}^{\kappa},\mathbf{H}^{\kappa},\mathbf{Q}^{\kappa}\right)$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}R^{lb}\left(\mathbf{A}^{\kappa+1},\mathbf{P}^{\kappa+1},\mathbf{\Theta}^{\kappa},\mathbf{H}^{\kappa},\mathbf{Q}^{\kappa}\right)$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}R^{lb}\left(\mathbf{A}^{\kappa+1},\mathbf{P}^{\kappa+1},\mathbf{\Theta}^{\kappa},\mathbf{H}^{\kappa},\mathbf{Q}^{\kappa+1}\right)$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{\leq}}R\left(\mathbf{A}^{\kappa+1},\mathbf{P}^{\kappa+1},\mathbf{\Theta}^{\kappa},\mathbf{H}^{\kappa},\mathbf{Q}^{\kappa+1}\right).$
(29)

Among them, (a) holds because given $\mathbf{Q}^{\kappa}$, $\mathcal{P}_{3.1}$
and $\mathcal{P}_{3.2}$ has the same objective function value; (b) is bounded
because given $\mathbf{A}^{\kappa+1}$, $\mathbf{H}^{\kappa}$ and
$\mathbf{P}^{\kappa+1}$, $\mathbf{Q}^{\kappa+1}$ is maximizing $R^{\rm{lb}}$
obtained from $\mathcal{P}_{3.2}$; (c) is established because the target value
of problem $\mathcal{P}_{3.2}$ is the lower bound of the target value of its
original problem $\mathcal{P}_{3.2}$ at $\mathbf{Q}^{\kappa+1}$. The
inequality in (IV-E) indicates that although the UAV trajectory was obtained
by solving approximate optimization problem, the objective value of problem
$\mathcal{P}_{3.2}$ still does not decrease after each iteration. Finally, for
step 4 in Algorithm 1, it is also achieved by fixing other feasible variables
$\mathbf{A}^{\kappa+1}$, $\mathbf{Q}^{\kappa+1}$ and $\mathbf{P}^{\kappa+1}$,
similar to step 3, we obtain

$\displaystyle
R\left(\mathbf{A}^{\kappa+1},\mathbf{P}^{\kappa+1},\mathbf{\Theta}^{\kappa},\mathbf{H}^{\kappa},\mathbf{Q}^{\kappa+1}\right)$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}R^{lb}\left(\mathbf{A}^{\kappa+1},\mathbf{P}^{\kappa+1},\mathbf{\Theta}^{\kappa},\mathbf{H}^{\kappa},\mathbf{Q}^{\kappa+1}\right)$
$\displaystyle\stackrel{{\scriptstyle(b)}}{{\leq}}R^{lb}\left(\mathbf{A}^{\kappa+1},\mathbf{P}^{\kappa+1},\boldsymbol{\Theta}^{\kappa+1},\mathbf{H}^{\kappa+1},\mathbf{Q}^{\kappa+1}\right)$
$\displaystyle\stackrel{{\scriptstyle(c)}}{{\leq}}R\left(\mathbf{A}^{\kappa+1},\mathbf{P}^{\kappa+1},\mathbf{\Theta}^{\kappa+1},\mathbf{H}^{\kappa+1},\mathbf{Q}^{\kappa+1}\right).$
(30)

This is similar to the representation in (IV-E), and from (27) to (IV-E), we
obtain

$\displaystyle
R\left(\mathbf{A}^{\kappa},\mathbf{P}^{\kappa},\mathbf{\Theta}^{\kappa},\mathbf{H}^{\kappa},\mathbf{Q}^{\kappa}\right)$
$\displaystyle\leq
R\left(\mathbf{A}^{\kappa+1},\mathbf{P}^{\kappa+1},\mathbf{\Theta}^{\kappa+1},\mathbf{H}^{\kappa+1},\mathbf{Q}^{\kappa+1}\right).$
(31)

The above analysis indicates that the target value of $\mathcal{P}_{0}$ does
not decrease after each iteration of Algorithm 1. Due to the objective value
of $\mathcal{P}_{0}$ is a finite upper bound, therefore the proposed Algorithm
1 ensures convergence. The simulation results in the next section indicate
that the proposed BCD-based method converges rapidly for the setting we are
considering. In addition, since only convex optimization problems need to be
solved in each iteration of Algorithm 1, which have polynomial complexity,
Algorithm 1 can actually converge quickly for wireless networks with a
moderate number of users.

In summary, by alternately solving the four subproblems $\mathcal{P}_{1.1}$,
$\mathcal{P}_{2.1}$, $\mathcal{P}_{3.2}$ and $\mathcal{P}_{4.2}$, until the
algorithm converges to the predetermined accuracy $\varepsilon>0$, we can
obtain the suboptimal solution to problem $\mathcal{P}_{0}$. In addition, the
computational complexity of this algorithm is
$O\left(N^{3.5}\log(1/\varepsilon)\right)$.

((a)) The average rate with different schemes.

((b)) The average rate with different power limitation.

((c)) The user scheduling.

Figure 3: The average rate and user scheduling.

((a)) 3D trajectories of $B$ with different schemes.

((b)) 3D trajectories of $B$ under different scenarios.

Figure 4: 3D trajectories of $B$ under different schemes and scenarios.

((a)) The horizontal trajectory of $B$.

((b)) The vertical attitude of $B$.

((c)) The transmit power of $B$.

((d)) The achievable rate of $B$.

((e)) The average rate versus varying $T$ with different schemes.

((f)) The average rate versus varying $\Gamma$ with different schemes.

Figure 5: Simulation results with different schemes.

((a)) The horizontal trajectory of $B$.

((b)) The vertical attitude of $B$.

((c)) The transmit power of $B$.

((d)) The horizontal speed of $B$.

((e)) The vertical speed of $B$.

((f)) The achievable rate of $B$.

Figure 6: Simulation results of the proposed scheme for varying $P_{\lim}$ and
$\Gamma$.

##  V Numerical Results

In this section, extensive simulations are conducted to validate the
effectiveness and the convergence of the proposed algorithm. To illustrate the
advantages of the proposed algorithm, the following three benchmark schemes
are compared.

  1. 1.

Benchmark I: $B$ transmit signals with the average power and both 3D
trajectory of $B$ and user scheduling are jointly optimized based on the PLoS
model, which is denoted by ‘NPC’.

  2. 2.

Benchmark II: The 2D trajectory of $B$, the transmission power of $B$, and
user scheduling are jointly optimized based on the LoS model, which is denoted
by ‘2D-LoS’.

  3. 3.

Benchmark III: $B$ works with fixed vertical trajectory and its horizontal and
transmission power of $B$, and user scheduling are jointly optimized based on
the PLoS model, which is denoted by ‘2D-PLoS’.

TABLE II: List of Simulation Parameters. Parameter | Value
---|---
$\left[{{{\bf{q}}_{B}}\left[1\right],{z_{B}}\left[1\right]}\right]$ | ${\left[{{\rm{336,187,30}}}\right]^{T}}$
${{\mathbf{w}}_{{U_{r}}}}$ |  ${\left[{162,23;112,301;332,50;298,381}\right]^{T}}$
${{\mathbf{w}}_{{D}}}$ |  ${\left[{200,360}\right]^{T}}$
$H_{\min},H_{\max}$ | 30 m, 100 m
${\Gamma}$ | -121 dB
${\sigma^{2}}$ | -100 dBm
$P_{B}^{\max}$ | 0.4 W
${P_{B}^{\text{ave }}}$ | 0.1 W
${\alpha_{L}},{\alpha_{N}}$ | 2.2, 3.5
${a},{b}$ | 11.95, 0.14
$V_{\max}$ | 25 m/s
${\hat{V}_{\max}}$ | 10 m/s
$a_{\max}$ | 6 ${\rm{m/}}{{\rm{s}}^{2}}$
${\rho_{0}}$ | -60 dB
${\mu}$ | -20 dB
${W}$ | 100 N
${T}$ | 70 s
${\delta_{t}}$ | 1 s
$\varepsilon$ | 0.0001
${P_{\text{hor }}^{\text{ave }}}$ | 880 W
${P_{\text{ver }}^{\text{ave }}}$ | 265 W

Similar to [5], the initial flight trajectory of $B$ is set as a circular
trajectory with the center at the geometric center of $U_{r}$, which is
expressed as $\mathbf{C}=\frac{1}{R}\sum\limits_{r=1}^{R}\mathbf{w}_{U_{r}}$.
To access all CUs as much as possible, the initial flight radius of the UAV is
$R_{B}=\min\left(\frac{V_{\max}T}{2\pi},\left\|\mathbf{C}-\mathbf{w}_{U_{r}}\right\|\right)$
[5, 27]. Details of parameter settings of the UAV communication system are
shown in Table II [14], [15], [20], [21], [22].

Fig. 3 demonstrates the convergence of all the schemes and user scheduling
results. Fig. 3(a) indicates that all the schemes converge fast, and the
performance with 2D-LoS and 2D-PLOS have the best and the worse performance,
respectively. Fig. 3(b) plots the rate with varying
$P_{{\rm{hor}}}^{{\rm{ave}}}$, $P_{{\rm{ver}}}^{{\rm{ave}}}$, and $\Gamma$ and
demonstrates that the impact of the power limitation and IT threshold on the
average rate of the considered system. Fig. 3(c) indicates that each user is
scheduled in turn.

Figs. 4(a) and 4(b) provide the 3D trajectory of $B$ with the different
schemes and in the scenarios with different propulsion power limitations,
respectively. From Fig. 4(a), it can be observed that compared with the
scenarios with the 2D-PLoS scheme, $B$ in the 2D-LoS scheme must keep away
from $D$ to meet IT constraints. This requirement for $B$ to keep away from
$D$ in the 2D-LoS scheme is due to the exaggerated probability of LoS for the
A2G link. The results in Fig. 4(b) show that different propulsion power
limitation results in different vertical trajectories.

Fig. 5 demonstrates the simulation results under different schemes with
$P_{{\rm{hor}}}^{{\rm{ave}}}=880$ W, $P_{{\rm{ver}}}^{{\rm{ave}}}=265$ W, and
$\Gamma=-121$ dB. Figs. 5(a), 5(b), 5(c), and 5(d) plot the horizontal
trajectory, vertical trajectory, transmit power, and instantaneous achievable
rate of different schemes, respectively. It can be observed from Fig. 5(a)
that different schemes achieve the same horizontal trajectory in the region
away from $D$. However, in the area near $D$, the results of the proposed
scheme, NPC, and 2D-PLoS are different from those of the 2D-LoS scheme. This
is because the A2G link is treated as LoS in the 2D-LoS scheme, $B$ must move
horizontally away from $D$ to satisfy the IT constraint. As can be seen in
Fig. 5(b), when $B$ moves away from the TNs, it first raises its altitude and
then dives down to get a better elevation angle, which also can be found from
Figs. 4(a) and 4(b) and stated in [25]. Fig. 5(c) shows the the transmit power
is optimized to meet the IT constraint and maximize performance. In
particular, in the 2D-PLoS scheme, the higher power is used when the UAV is
near TNs, while the lower power is used when the UAV is far away from TNs.
That is, performance is maximized by allocating different power at different
locations. Moreover, in the 2D-LoS scheme, since A2G links are treated as LoS,
the transmit power is constrained by the IT constraint. Specifically, the
transmit power in 2D-LoS decreases as the distance between $B$ and PU
increases. As can be seen from Fig. 5(d), the considered system with the
2D-LoS has the lowest performance in the areas near $D$ because the lowest
transmit power is used subject to the IT constraints. The proposed scheme’s
performance in the regions near the TNs is between that of the NPC and the
2D-PLoS schemes because both the transimit power and the vertical trajectory
are optimized to obtain a larger elevation angle (the larger the probability
of LoS). Figs. 5(e) and 5(f) show the relationship between the average rate
and varying $T$ and IT thresholds, respectively. From Fig. 5(e), one can
observe that, same as in Fig. 3(a), the performance improves as $T$ increases
because $B$ can spend more time hovering over TNs. The 2D-LoS scheme has the
best performance. The reason is, in the areas away from PU, more power can be
used to maximize the achievable rate. In addition, in the 2D-LoS scheme, the
A2G link is always considered LoS $\left(P_{BX}^{\mathrm{L}}=1\right)$. Fig.
5(f) demonstrates the effect of the IT thresholds on the average rate with
different schemes. One can observe that the average rate increases as the IT
threshold becomes more relaxed. This is because, with the higher IT threshold,
$B$ can use higher transmit power and thus obtain higher performance. When
$\Gamma$ increases to a specific range, the system enters the no-cognitive
state, and the system performance is independent of the IT threshold. In the
low-$\Gamma$ region, the requirement for the CUs is very harsh and the
transmit power is limited to the IT constraint. Specifically, the NPC scheme
does not converge when $\Gamma=-128$ dB. Limited by the IT constraint, the
performance of the 2D-LoS scheme is the worst in the low-$\Gamma$ region.
However, in the high-$\Gamma$ region, the system enters the no-cognitive state
and higher transmission power can be used, thus, the 2D-LoS has the best
performance.

Fig. 6 describes the simulation results of the proposed scheme with varying
$P_{{\rm{hor}}}^{{\rm{ave}}}$, $P_{{\rm{ver}}}^{{\rm{ave}}}$, and $\Gamma$.
Fig. 6(a) show that, in the scenarios with lower propulsion power and lower IT
thresholds, $B$ is needed away from $D$ in the region near $D$ to meet the IT
constraint. This is because the lower propulsion power responds to the lower
horizontal speed, shown at Fig. 6(d), and the lower IT threshold denotes that
the requirements for the cognitive users are more strict, so $B$ needs to be
away from $D$. With the same reason, in the area close to $D$, the elevation
angle cannot reach the maximum to ensure that IT constraints are met,
presented in Fig. 6(b). In areas away from the PU, IT constraints are not
limited and vertical trajectory is maximized. Fig. 6(c) shows that in the
scenarios with small $\Gamma$, the UAV uses lower power in the area near the
PU. In areas away from PU, higher power is used. Figs. 6(d) and 6(e) plot the
optimized horizontal and vertical speeds for the different power and IT
thresholds, respectively. It can be observed that the larger propulsion power
budget, the longer hovering time and the longer fall time, which result in
higher performance.

##  VI Conclusion

This work investigated the achievable rate of the underlay IoT system with an
energy-constrained UAV under PLoS channels. The achievable rate of the
considered systems was maximized by jointly considering the UAV’s 3D
trajectory, transmission power, and user scheduling, which is a nonlinear
mixed-integer non-convex problem. The lower bound of the average achievable
rate was utilized and the original non-convex problem was transformed into
several solvable convex subproblems by using BCD and SCA techniques, and an
efficient iterative algorithm was proposed. The numerical results not only
verified the convergence and effectiveness of the algorithm but also
illustrated the impact of propulsion power and interference thresholds on the
average achievable rate.

##  Appendix A Proof of Lemma 1

For simplicity, we define a bivariate function
$f(x,y)=\frac{1}{x}\log_{2}\left(1+\frac{A}{y}\right)$, where $A\geq 0$, $x$
and $y$ are both positive variables of the function. The Jacobian of $f(x,y)$
is

$\displaystyle\nabla f(x,y)$
$\displaystyle=\left[\begin{array}[]{ll}\frac{\partial f(x,y)}{\partial
x}&\frac{\partial f(x,y)}{\partial y}\end{array}\right]$ (33)
$\displaystyle=\left[\begin{array}[]{ll}-\frac{1}{x^{2}}\log_{2}(1+\frac{A}{y})&-\frac{A\log_{2}e}{xy^{2}(1+\frac{A}{y})}\end{array}\right].$
(35)

From the above Jacobian matrix, we can obtain the Hessian matrix of $f(x,y)$
as

$\displaystyle\nabla^{2}f(x,y)$
$\displaystyle=\left[\begin{array}[]{ll}\frac{\partial^{2}f(x,y)}{\partial
x^{2}}&\frac{\partial^{2}f(x,y)}{\partial x\partial y}\\\
\frac{\partial^{2}f(x,y)}{\partial y\partial
x}&\frac{\partial^{2}f(x,y)}{\partial y^{2}}\end{array}\right]$ (38)
$\displaystyle=\left[\begin{array}[]{cc}\frac{2}{x^{3}}\log_{2}(1+\frac{A}{y})&\frac{A\log_{2}e}{x^{2}y^{2}(1+\frac{A}{y})}\\\
\frac{A\log_{2}e}{x^{2}y^{2}(1+\frac{A}{y})}&\frac{A(2y+A)\log_{2}e}{xy^{4}(1+\frac{A}{y})^{2}}\end{array}\right].$
(41)

We need to scale $\log_{2}(1+\frac{1}{x})$ to take its lower bound. When
$x\geq 0$, we can use Lagrange’s mean value theorem to determine
$\log_{2}(1+\frac{1}{x})$$\geq$$\frac{\log_{2}e}{1+x}$. Therefore, we replace
it with the lower bound of $\log_{2}(1+\frac{A}{y})$ and find the determinant
of its Hessian matrix. Furthermore, the determinant of the Hessian matrix of
$f(x,y)$ is

$\displaystyle\left|\nabla^{2}f(x,y)\right|$
$\displaystyle=\frac{2}{x^{4}y^{4}}\frac{A(2y+A)\log_{2}e}{(1+\frac{A}{y})^{2}}\log_{2}\left(1+\frac{A}{y}\right)$
$\displaystyle-\frac{1}{x^{4}y^{4}}\frac{A^{2}(\log_{2}e)^{2}}{(1+\frac{A}{y})^{2}}$
$\displaystyle\geq\frac{1}{x^{4}y^{4}}\frac{A^{2}(2y+A)(\log_{2}e)^{2}}{(1+\frac{A}{y})^{2}(y+A)}$
$\displaystyle-\frac{1}{x^{4}y^{4}}\frac{A^{2}(\log_{2}e)^{2}}{(1+\frac{A}{y})^{2}}$
$\displaystyle=\frac{1}{x^{4}y^{4}}\frac{A^{2}y(\log_{2}e)^{2}}{(1+\frac{A}{y})^{2}(y+A)}\geq
0.$ (42)

According to (38) and (42), for any given $A\geq 0$, We can determine that the
determinant of the cofactor matrix of the Hessian matrix of function $f(x,y)$
is greater than or equal to 0. Therefore, the Hessian matrix of $f(x,y)$ is a
semi-positive definite matrix, so $f(x,y)$ is a convex function.

## References

  * [1] Q. Wu, J. Xu, Y. Zeng, D. W. K. Ng, N. Al-Dhahir, R. Schober, and A. L. Swindlehurst, “A comprehensive overview on 5G-and-beyond networks with UAVs: From communications to sensing and intelligence,” IEEE J. Sel. Areas Commun., vol. 39, no. 10, pp. 2912-2945, Oct. 2021.
  * [2] B. Li, Z. Fei, and Y. Zhang, “UAV communications for 5G and beyond: recent advances and future trends,” IEEE Internet Things J., vol. 6, no. 2, pp. 2241-2263, Apr. 2019.
  * [3] 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.
  * [4] Q. Wu, L. Liu, and R. Zhang, “Fundamental trade-offs in communication and trajectory design for UAV-enabled wireless network,” IEEE Wireless Commun., vol. 26, no. 1, pp. 36-44, Feb. 2019.
  * [5] Q. Wu, Y. Zeng, and R. Zhang, “Joint trajectory and communication design for multi-UAV enabled wireless networks,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 3, pp. 2109-2121, Mar. 2018.
  * [6] Z. Li, M. Chen, C. Pan, N. Huang, Z. Yang, and A. Nallanathan, “Joint trajectory and communication design for secure UAV networks,” IEEE Commun. Lett., vol. 23, no. 4, pp. 636-639, Apr. 2019.
  * [7] H. Fu, Z. Sheng, A. A. Nasir, A. H. Muqaibel, and L. Hanzo, “Securing the UAV-aided non-orthogonal downlink in the face of colluding eavesdroppers,” IEEE Trans. Veh. Technol., vol. 71, no. 6, pp. 6837-6842, Jun. 2022.
  * [8] M. Cui, G. Zhang, Q. Wu, and D. W. K. Ng, “Robust trajectory and transmit power design for secure UAV communications,” IEEE Trans. Veh. Technol., vol. 67, no. 9, pp. 9042-9046, Sep. 2018.
  * [9] 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.
  * [10] S. Eom, H. Lee, J. Park, and I. Lee, “UAV-aided wireless communication designs with propulsion energy limitations,” IEEE Trans. Veh. Technol., vol. 69, no. 1, pp. 651-662, Jan. 2020.
  * [11] Y. Wu, W. Yang, X. Guan, and Q. Wu, “Energy-efficient trajectory design for UAV-enabled communication under malicious jamming,” IEEE Wireless Commun. Lett., vol. 10, no. 2, pp. 206-210, Feb. 2021.
  * [12] J. Zhang, Y. Zeng, and R. Zhang, “Receding horizon optimization for energy-efficient UAV communication,” IEEE Wireless Commun. Lett., vol. 9, no. 4, pp. 490-494, Apr. 2020.
  * [13] A. Bejaoui, K.-H. Park, and M.-S. Alouini, “A QoS-oriented trajectory optimization in swarming unmanned-aerial-vehicles communications,” IEEE Wireless Commun. Lett., vol. 9, no. 6, pp. 791-794, Jun. 2020.
  * [14] 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.
  * [15] C. Zhan and H. Lai, “Energy minimization in internet-of-things system based on rotary-wing UAV,” IEEE Wireless Commun. Lett., vol. 8, no. 5, pp. 1341-1344, Oct. 2019.
  * [16] 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, Apr. 2020.
  * [17] R. Zhang, X. Pang, W. Lu, N. Zhao, Y. Chen, and D. Niyato, “Dual-UAV enabled secure data collection with propulsion limitation,” IEEE Trans. Wireless Commun., vol. 20, no. 11, pp. 7445-7459, Nov. 2021.
  * [18] M. Li, X. Tao, N. Li, H. Wu, and J. Xu, “Secrecy energy efficiency maximization in UAV-enabled wireless sensor networks without eavesdropper’s CSI,” IEEE Internet Things J., vol. 9, no. 5, pp. 3346-3358, Mar. 2022.
  * [19] A. Al-Hourani, S. Kandeepan, and S. Lardner, “Optimal LAP altitude for maximum coverage,” IEEE Wireless Commun. Lett., vol. 3, no. 6, pp. 569-572, Dec. 2014.
  * [20] B. Duo, H. Hu, Y. Li, and X. Yuan, “Joint trajectory and power design in probabilistic LoS channel for UAV-enabled cooperative jamming,” Proc. IEEE International Conference on Communications (ICC 2021), Montreal, QC, Canada, Jun. 2021, pp. 1-6.
  * [21] 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, Jun. 2020.
  * [22] 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. 16288-16293, Dec. 2020.
  * [23] Z. Guan, S. Wang, L. Gao, and W. Xu, “Energy-efficient UAV communication with 3D trajectory optimization,” Proc. 2021 7th International Conference on Computer and Communications (ICCC), Chengdu, China, Dec. 2021, pp. 312-317.
  * [24] A. Meng, X. Gao, Y. Zhao, and Z. Yang, “Three-dimensional trajectory optimization for energy-constrained UAV-enabled IoT system in probabilistic LoS channel,” IEEE Internet Things J., vol. 9, no. 2, pp. 1109-1121, Jan. 2022.
  * [25] H. Lei, H. Ran, I. S. Ansari, K.-H. Park, G. Pan, and M.-S. Alouini, “DDPG-based aerial secure data collection,” IEEE Trans. Commun., doi: 10.1109/TCOMM.2024.3379417, 2024.
  * [26] Y. Saleem, M. H. Rehmani, and S. Zeadally, “Integration of cognitive radio technology with unmanned aerial vehicles: Issues, opportunities, and future research challenges,” J. Netw. Comput. Appl., vol. 50, pp. 15-31, Apr. 2015.
  * [27] H. Lei, H. Yang, K.-H. Park, I. S. Ansari, J. Jiang, and M.-S. Alouini, “Joint trajectory design and user scheduling for secure aerial underlay IoT systems,” IEEE Internet Things J., vol. 10, no. 15, pp. 13637-13648, Aug. 2023.
  * [28] H. Lei, J. Jiang, H. Yang, K.-H. Park, I. S. Ansari, G. Pan, and M.-S. Alouini. “Trajectory and power design for aerial CRNs with colluding eavesdroppers,” arXiv:2310.13931, Oct. 2023, [Online]: https://arxiv.org/abs/2310.13931
  * [29] P. X. Nguyen, V.-D. Nguyen, H. V. Nguyen, and O.-S. Shin, “UAV-assisted secure communications in terrestrial cognitive radio networks: Joint power control and 3D trajectory optimization,” IEEE Trans. Veh. Technol., vol. 70, no. 4, pp. 3298-3313, Apr. 2021.
  * [30] Z. Wang, F. Zhou, Y. Wang, and Q. Wu, “Joint 3D trajectory and resource optimization for a UAV relay-assisted cognitive radio network,” China Commun., vol. 18, no. 6, pp. 184-200, Jun. 2021.
  * [31] Y. Huang, W. Mei, J. Xu, L. Qiu, and R. Zhang, “Cognitive UAV communication via joint maneuver and power control,” IEEE Trans. Commun., vol. 67, no. 11, pp. 7872-7888, Nov. 2019.
  * [32] Y. Jiang and J. Zhu, “Three-dimensional trajectory optimization for secure UAV-enabled cognitive communications,” China Commun., vol. 18, no. 12, pp. 285-296, Dec. 2021.
  * [33] Y. Sun, D. Xu, D. W. K. Ng, L. Dai, and R. Schober, “Optimal 3D-trajectory design and resource allocation for solar-powered UAV communication systems,” IEEE Trans. Commun., vol. 67, no. 6, pp. 4281-4298, Jun. 2019.
  * [34] J. Yao and J. Xu, “Joint 3D maneuver and power adaptation for secure UAV communication with CoMP reception,” IEEE Trans. Wireless Commun., vol. 19, no. 10, pp. 6992-7006, Oct. 2020.
  * [35] E. Tzoreff and A. J. Weiss, “Path design for best emitter location using two mobile sensors,” IEEE Trans. Signal Process., vol. 65, no. 19, pp. 5249-5261, Oct. 2017.
  * [36] H. Xing, L. Liu, and R. Zhang, “Secrecy wireless information and power transfer in fading wiretap channel,” IEEE Trans. Veh. Technol., vol. 65, no. 1, pp. 180-190, Jan. 2016.
  * [37] T. Liu, M. Cui, G. Zhang, Q. Wu, X. Chu, and J. Zhang, “3D trajectory and transmit power optimization for UAV-enabled multi-link relaying systems,” IEEE Trans. Green Commun. Netw., vol. 5, no. 1, pp. 392-405, Mar. 2021.
  * [38] S. Boyd and L. Vandenberghe, Convex Optimization: Cambridge Univ Press, 2004.
  * [39] A. Filippone, Flight Performance of Fixed and Rotary Wing Aircraft. Amsterdam, The Netherlands: Elsevier, 2006.