# 3-D Positioning and Resource Allocation for Multi-UAV Base Stations Under
Blockage-Aware Channel Model

Pengfei Yi, Lipeng Zhu, Zhenyu Xiao, Rui Zhang, Zhu Han, and Xiang-Gen Xia,
This work was supported in part by the National Key Research and Development
Program under grant number 2020YFB1806800, the National Natural Science
Foundation of China (NSFC) under grant numbers U22A2007, 62171010 and
61827901, and the Beijing Natural Science Foundation under grant number
L212003. The corresponding author is Dr. Zhenyu Xiao with Email
<EMAIL_ADDRESS>Yi and Z. Xiao are with the School of Electronic and
Information Engineering, Beihang University, Beijing 100191, China.
<EMAIL_ADDRESS>xiaozy@buaa.edu.cn).L. Zhu is with the Department of
Electrical and Computer Engineering, National University of Singapore,
Singapore 117583, Singapore. (zhulp@nus.edu.sg).R. Zhang is with The Chinese
University of Hong Kong, Shenzhen, and Shenzhen Research Institute of Big
Data, Shenzhen, China 518172 (rzhang@cuhk.edu.cn). He is also with the
Department of Electrical and Computer Engineering, National University of
Singapore, Singapore 117583 (elezhang@nus.edu.sg).Z. Han is with the
Department of Electrical and Computer Engineering in the University of
Houston, Houston, TX 77004 USA, and also with the Department of Computer
Science and Engineering, Kyung Hee University, Seoul, South Korea, 446-701.
(hanzhu22@gmail.com).X.-G. Xia is with the Department of Electrical and
Computer Engineering, University of Delaware, Newark, DE 19716, USA.
(xianggen@udel.edu).

###### Abstract

In this paper, we propose to deploy multiple unmanned aerial vehicle (UAV)
mounted base stations to serve ground users in outdoor environments with
obstacles. In particular, the geographic information is employed to capture
the blockage effects for air-to-ground (A2G) links caused by buildings, and a
realistic blockage-aware A2G channel model is proposed to characterize the
continuous variation of the channels at different locations. Based on the
proposed channel model, we formulate the joint optimization problem of UAV
three-dimensional (3-D) positioning and resource allocation, by power
allocation, user association, and subcarrier allocation, to maximize the
minimum achievable rate among users. To solve this non-convex combinatorial
programming problem, we introduce a penalty term to relax it and develop a
suboptimal solution via a penalty-based double-loop iterative optimization
framework. The inner loop solves the penalized problem by employing the block
successive convex approximation (BSCA) technique, where the UAV positioning
and resource allocation are alternately optimized in each iteration. The outer
loop aims to obtain proper penalty multipliers to ensure the solution of the
penalized problem converges to that of the original problem. Simulation
results demonstrate the superiority of the proposed algorithm over other
benchmark schemes in terms of the minimum achievable rate.

###### Index Terms:

UAV communication, geographic information, 3-D positioning, resource
allocation, blockage.

##  I Introduction

In recent years, unmanned aerial vehicle (UAV)-assisted communication systems
have attracted increasing attention for supporting the seamless coverage in
the beyond fifth-generation (B5G) and sixth-generation (6G) networks [1, 2, 3,
4]. Owing to their controllable three-dimensional (3-D) mobility and low cost,
UAVs can serve as aerial base stations (BSs), relays, or access points for
coverage enhancement [5, 6, 7], communication relaying [8, 9, 10], and data
broadcast/collection [11, 12, 13]. Compared to conventional terrestrial
communications with typically fixed infrastructures, UAV-assisted systems
offer new degrees of freedom in the spatial domain to further improve
communication performance by exploiting the flexible 3-D mobility of UAVs.

In this regard, several interesting topics arise in the study of UAV-assisted
communications, such as UAV placement, trajectory design, and resource
allocation. In particular, with the increasing number of users, their wider
distribution and various communication requirements make it necessary to
employ multiple UAVs forming a cooperative network to improve the access
capability, enlarge the coverage area, and enhance the communication
reliability. In [14], a successive UAV placement strategy was proposed to
minimize the number of required UAVs while satisfying the communication
requirement of ground users. The authors in [15] studied the joint
optimization for multi-UAV placement, user association, and resource
allocation, to maximize the downlink sum rate. Multi-UAV placement and user
association were also considered in [16], with a particular consideration on
constrained backhaul links. Aiming to support energy-efficient Internet of
Things (IoT) communications, multiple UAVs were deployed and the communication
resources were jointly optimized, to minimize the transmit power of IoT
devices in [11]. Combining millimeter-wave (mmWave) communications with UAV-
assisted systems, the authors in [6] investigated the joint optimization of
UAV placement, user clustering, and transmit/receive beamforming. In addition
to the placement optimization, there are also many works on the trajectory
design for UAV communication systems, aiming at ubiquitous coverage [5],
secure communications [17], interference coordination [18], energy-efficient
content coverage [19], and mobile Internet of vehicles [20].

Note that the aforementioned works significantly rely on the
simplified/statistical channel models for air-to-ground (A2G) communication
links between UAVs and ground users. For example, the A2G channels are assumed
to be dominated by line-of-sight (LoS) paths in [14, 5, 15, 17, 18, 20]. While
the existence of LoS links is probabilistically modeled as a function of the
elevation angle of the A2G link in [11, 19, 16], known as _probabilistic LoS
channels_ [2]. The simplified LoS channel models can make the positioning
optimization more tractable and are suitable for average performance analysis
in UAV communications. However, in practice, the terrain conditions, such as
buildings and other obstacles, may cause severe blockage to A2G links and
sharply weaken the strength of the received signals, especially for dense
urban areas [21]. In such cases, the communication design based on a
simplified LoS channel and/or statistical channel cannot guarantee the
performance under site-specific environments and may not be suitable for
practical UAV-assisted communications.

To overcome the drawback of the over-simplified channel models and capture
practical propagation conditions, there are emerging research directions that
exploit two kinds of information, namely _radio map_ and _geographic
information_ , for UAV-assisted communications. Constructed by a large number
of real-life channel measurements, a radio map can precisely describe the
average signal strength for all combinations of UAV-user locations [22]. Based
on the radio map, the joint optimization of UAV positioning, user association,
and wireless backhaul capacity allocation for a multi-UAV relay network was
studied in [7]. In [10], the UAV trajectory and resource allocation were
jointly optimized to ensure the fairness among users for a single UAV relay
system. Besides, a radio map was utilized to evaluate the A2G link quality for
the UAV trajectory design in a cellular-connected UAV system [23] and UAV
anti-jamming communications [24]. Radio map is theoretically appealing to
provide precise channel quality, but it encounters difficulty in obtaining
sufficient real-life channel measurements for radio map construction. Besides,
a large dataset is required for reconstructing a radio map, which results in
high overhead on storage and computation. In addition to radio map, geographic
information is also helpful for capturing practical propagation conditions.
With the available location and size information of the buildings, the
existence of a LoS channel can be inferred by evaluating whether the A2G link
is blocked by buildings, instead of being modeled as a random event [2]. Based
on geographic information, the authors in [25] proposed a geometric analysis
method to detect the blockage in a multi-UAV mmWave communication system, and
then developed a greedy user scheduling algorithm to decrease the probability
of the blockage. In [12], building blockage was considered for the trajectory
design and resource allocation for a UAV-enabled data collection system. In
[8], the blocked regions of the ground users with respect to (w.r.t.)
buildings were modeled as polyhedrons. By restricting the UAV to be deployed
outside all the blocked regions, the LoS links can be guaranteed for a UAV
relay system.

Motivated by the above works, in this paper, we study the joint optimization
of UAV 3-D positioning, power allocation, user association, and subcarrier
allocation to maximize the minimum achievable rate among multiple users in the
downlink of multi-UAV orthogonal frequency division multiple access (OFDMA)
communication systems, under the blockage-aware A2G channel model with the aid
of geographic information. Different from the works based on
deterministic/probabilistic LoS channel [14, 15, 11, 16], the geographic
information utilized in this paper can predict the channel condition with a
high precision. Compared to the works based on radio map [7, 10], this work
only needs the geographic information, which is easier to acquire in reality.
For example, geographic information can be derived offline from digital maps
[26], or constructed online by using photogrammetry techniques [27]. The main
contributions of this paper are summarized as follows:

  1. 1.

We propose to deploy multiple UAV BSs to serve multiple ground users via
OFDMA. Geographic information is utilized to capture the realistic propagation
environment. Specifically, we develop a blockage-aware A2G channel model,
where the LoS and Non-LoS (NLoS) channels can be uniformly expressed as a
continuous function of the normalized distance between the UAV and the ground
region. Then, to ensure fairness, we formulate an optimization problem to
maximize the minimum achievable rate among the users by jointly designing the
3-D positioning and resource allocation, including power allocation, user
association, and subcarrier allocation.

  2. 2.

The formulated optimization problem is non-convex and involves combinatorial
programming variables, and thus it is difficult to obtain the globally optimal
solution. Therefore, we develop a penalty-based double-loop iterative
optimization (PDLIO) algorithm to solve the problem suboptimally.
Specifically, the original problem is transformed to a penalized problem by
relaxing the binary association variables into continuous ones and introducing
a penalty component to the objective function. The inner-loop partitions the
penalized problem into a UAV positioning sub-problem and a resource allocation
sub-problem, which are alternately solved by employing the block successive
convex approximation (BSCA) technique [28, 29]. The outer-loop updates the
penalty multipliers to ensure the solution of the penalized problem converges
to that of the original problem.

  3. 3.

The performance of the proposed solution for the joint 3-D positioning and
resource allocation problem in geographic information-aided multi-UAV systems
is evaluated with different settings. Simulation results illustrate the
convergence of the developed PDLIO algorithm and reveal performance
superiority over other benchmark schemes in terms of the minimum achievable
rate.

The rest of this paper is organized as follows. In Section II, we introduce
the system model, propose the blockage-aware A2G channel model based on
geographic information, and formulate the joint positioning and resource
allocation problem. The problem transformation and the proposed PDLIO
algorithm are given in Section III. Section IV presents the simulation
results. Finally, the paper is concluded in Section V.

Notation: $a$, $\mathbf{a}$, $\mathbf{A}$, and $\mathcal{A}$ denote a scalar,
a vector, a matrix, and a set, respectively. $[\mathbf{a}]_{k}$ denotes the
$k$-th entry of vector $\mathbf{a}$. $\|\mathbf{a}\|$ represents the Euclidean
norm of vector $\mathbf{a}$. $(\cdot)^{\rm{T}}$ denotes transpose.
$\mathbb{R}^{M}$ denotes the space of the $M$-dimensional real vector.
$\mathcal{A}_{1}\bigcup\mathcal{A}_{2}$ represents the union of sets
$\mathcal{A}_{1}$ and $\mathcal{A}_{2}$.
$\mathcal{A}_{2}\setminus\mathcal{A}_{1}$ represents the elements of
$\mathcal{A}_{2}$ that are not included in $\mathcal{A}_{1}$. $|\mathcal{A}|$
denotes the cardinality of set $\mathcal{A}$. $\overrightarrow{AB}$ denotes
the vector from point $A$ to point $B$.
$\overrightarrow{AB}\cdot\overrightarrow{CD}$ and
$\overrightarrow{AB}\times\overrightarrow{CD}$ denote the inner product and
outer product between vector $\overrightarrow{AB}$ and vector
$\overrightarrow{CD}$, respectively. For a multivariate function
$f(\mathbf{x},\mathbf{y})$, $\nabla_{\mathbf{x}}f(\mathbf{x},\mathbf{y})$
denotes its gradient w.r.t. $\mathbf{x}$. $f(\mathbf{x};\mathbf{x}^{l})$
denotes a surrogate function of $f(\mathbf{x})$ constructed at local point
$\mathbf{x}^{{l}}$.

##  II System Model and Problem Formulation

As shown in Fig. 1, we consider a multi-UAV OFDMA downlink communication
network in outdoor environments with obstacles, such as buildings, where $M$
UAV BSs are deployed to serve $K$ ground users via $N$ orthogonal subcarriers
for each UAV. The sets of UAVs, users, and subcarriers are denoted as
$\mathcal{M}\triangleq\\{1,...,M\\}$, $\mathcal{K}\triangleq\\{1,...,K\\}$,
and $\mathcal{N}\triangleq\\{1,...,N\\}$, respectively. We assume that each
UAV allocates orthogonal subcarriers to its served ground users. As a result,
the interference among these users served by the same UAV can be eliminated.
However, since all the UAVs share the same $N$ subcarriers, the users served
by different UAVs via the same subcarrier may suffer from mutual interference.

Figure 1:  Illustration of the considered multi-UAV OFDMA downlink
communication network.

Without loss of generality, we employ a 3-D Cartesian coordinate system. For
user $k$, its coordinates are denoted by
$\mathbf{u}_{k}\in\mathbb{R}^{3},k\in\mathcal{K}$. The coordinates of UAV $m$
are given by $\mathbf{x}_{m}\in\mathbb{R}^{3},m\in\mathcal{M}$. We denote
$\mathbf{X}=\left\\{\mathbf{x}_{m}|m\in\mathcal{M}\right\\}$. $Q$ buildings
indexed by $\mathcal{Q}\triangleq\\{1,...,Q\\}$ are randomly distributed in
the considered area, whose 3-D locations and sizes are assumed to be available
with the aid of geographic information111A possible way to extract building
information is to use an open source geographic database _OpenStreetMap_ ,
where the contour and height of a building are given by its raw data tagged
with _geometry_ and _height_ , respectively [30]..

###  II-A Channel Model

Note that the blockage effects caused by buildings may significantly
deteriorate the link quality, which should be considered in the channel
modeling. Since the UAV positioning is designed in a relatively large
timescale compared to the channel small-scale variation, we mainly focus on
the large-scale channel characteristics. The channel gain between user $k$ and
UAV $m$ is modeled as a function of UAV position $\mathbf{x}_{m}$, i.e.,

$g_{k}(\mathbf{x}_{m})=\frac{\beta_{k}(\mathbf{x}_{m})}{\|\mathbf{x}_{m}-\mathbf{u}_{k}\|^{\alpha_{k}(\mathbf{x}_{m})}},$
(1)

where $\alpha_{k}(\mathbf{x}_{m})$ is the path loss exponent, and
$\beta_{k}(\mathbf{x}_{m})$ is the channel gain at reference distance of 1
meter (m). Both $\alpha_{k}(\mathbf{x}_{m})$ and $\beta_{k}(\mathbf{x}_{m})$
are functions of UAV position $\mathbf{x}_{m}$ and user position
$\mathbf{u}_{k}$, which characterize the specific propagation conditions for
A2G links (LoS and/or NLoS channels) and are defined as

$\displaystyle\left(\alpha_{k}(\mathbf{x}_{m}),\beta_{k}(\mathbf{x}_{m})\right)=\begin{cases}\left(\alpha_{1},\beta_{1}\right),\text{for
LoS channel},\\\ \left(\alpha_{2},\beta_{2}\right),\text{for NLoS
channel}.\end{cases}$ (2)

Figure 2: Illustration of the blockage caused by building.

To distinguish LoS and NLoS channels for A2G links, the key is to model the
blockage effect between the UAV and user. The blocked region of user $k$
w.r.t. building $q$, denoted by $\mathcal{D}_{k,q}$, can be modeled as a
polyhedron

$\mathcal{D}_{k,q}=\\{\mathbf{x}\in\mathbb{R}^{3}|\mathbf{a}_{k,q,i}^{\rm{T}}\mathbf{x}-b_{k,q,i}\leq
0,i\in\mathcal{I}_{k,q}\\},$ (3)

where $\mathcal{I}_{k,q}$ is the set of indices of hyperplanes (boundaries)
for blocked region $\mathcal{D}_{k,q}$. $\mathbf{a}_{k,q,i}\in\mathbb{R}^{3}$
and $b_{k,q,i}\in\mathbb{R}$ are the outward normal vector and offset of the
$i$-th hyperplane, respectively, which are determined by the position of the
user and the visible flank surfaces of the building. The detailed procedures
to identify visible flank surfaces of the building and determine the
boundaries of the blocked region can be referred to [8, Algorithm 1]. Here, we
provide a toy example of an area with two buildings for illustration. As shown
in Fig. 2, the flank surface $A_{1}B_{1}B_{2}A_{2}$ of the first building is
visible for user $k$, and thus three hyperplanes, $SA_{1}A_{2}$,
$SA_{2}B_{2}$, and $SB_{2}B_{1}$, form the boundaries of blocked region
$\mathcal{D}_{k,1}$. The outward normal vector is obtained as
$\mathbf{a}_{k,1,1}=\frac{\overrightarrow{SA_{2}}\times\overrightarrow{SA_{1}}}{\|\overrightarrow{SA_{2}}\times\overrightarrow{SA_{1}}\|}$,
and the corresponding offset w.r.t. the origin $O$ is given by
$b_{k,1,1}=\mathbf{a}_{k,1,1}^{\rm{T}}\cdot\overrightarrow{OS}$, where $O$
denotes the origin of the coordinate system. The blocked region
$\mathcal{D}_{k,2}$ can be obtained in a similar manner, which is surrounded
by four hyperplanes.

With the union of the obtained blocked regions,
$\mathop{\bigcup}_{q\in\mathcal{Q}}\mathcal{D}_{k,q}$, of user $k$, the
existence of LoS path between user $k$ and UAV $m$ can be inferred directly by
evaluating whether the UAV is located in the blocked regions,
$\mathop{\bigcup}_{q\in\mathcal{Q}}\mathcal{D}_{k,q}$. If UAV $m$ is deployed
in any of the blocked regions of user $k$, i.e.,
$\mathbf{x}_{m}\in\mathop{\bigcup}_{q\in\mathcal{Q}}\mathcal{D}_{k,q}$, the
LoS path between them is blocked. Otherwise, for
$\mathbf{x}_{m}\notin\mathop{\bigcup}_{q\in\mathcal{Q}}\mathcal{D}_{k,q}$, the
LoS channel is guaranteed without any building blockage. Define

$d_{k,q}(\mathbf{x}_{m})\triangleq\max\limits_{i\in\mathcal{I}_{k,q}}\left\\{\mathbf{a}_{k,q,i}^{\rm{T}}\mathbf{x}_{m}-b_{k,q,i}\right\\},$
(4)

and then $\mathbf{x}_{m}\in\mathcal{D}_{k,q}$ is equivalent to
$d_{k,q}(\mathbf{x}_{m})\leq 0$. Therefore,
$\mathbf{x}_{m}\in\mathop{\bigcup}_{q\in\mathcal{Q}}\mathcal{D}_{k,q}$ (NLoS
condition) is equivalent to

$\min\limits_{q\in\mathcal{Q}}\left\\{d_{k,q}(\mathbf{x}_{m})\right\\}\leq 0,$

while
$\mathbf{x}_{m}\notin\mathop{\bigcup}_{q\in\mathcal{Q}}\mathcal{D}_{k,q}$ (LoS
condition) is equivalent to

$\min\limits_{q\in\mathcal{Q}}\left\\{d_{k,q}(\mathbf{x}_{m})\right\\}>0.$

As can be observed, the LoS/NLoS condition can be described as a _step
function_ w.r.t. $\min_{q\in\mathcal{Q}}$
$\left\\{d_{k,q}(\mathbf{x}_{m})\right\\}$. For the example shown in Fig. 2,
we have $\min_{q\in\\{1,2\\}}\left\\{d_{k,q}(\mathbf{x}_{m})\right\\}>0$,
because of
$d_{k,1}(\mathbf{x}_{m})=\max\limits_{i\in\mathcal{I}_{k,1}}\left\\{\mathbf{a}_{k,1,i}^{\rm{T}}\mathbf{x}_{m}-b_{k,1,i}\right\\}=\mathbf{a}_{k,1,1}^{\rm{T}}\mathbf{x}_{m}-b_{k,1,1}>0$
and
$d_{k,2}(\mathbf{x}_{m})=\max\limits_{i\in\mathcal{I}_{k,2}}\left\\{\mathbf{a}_{k,2,i}^{\rm{T}}\mathbf{x}_{m}-b_{k,2,i}\right\\}$
$=\mathbf{a}_{k,2,3}^{\rm{T}}\mathbf{x}_{m}-b_{k,2,3}>0$. Therefore, the A2G
link between UAV $m$ and user $k$ satisfies the LoS condition.

In order to facilitate the subsequent optimization procedure, the channel
model provided in (1) and (2) needs to be formulated as a continuous and
differentiable function w.r.t. the position of the UAV. Therefore, we employ
the _sigmoid function_ to approximate the step function, which is given by

$\displaystyle
s(\mathbf{x}_{m},\mathbf{u}_{k})=\frac{1}{1+\exp\left(-\eta\frac{\min_{q\in\mathcal{Q}}\left\\{d_{k,q}(\mathbf{x}_{m})\right\\}}{\|\mathbf{x}_{m}-\mathbf{u}_{k}\|}\right)},$
(5)

where the smooth parameter $\eta$ can be used to control the scale of the
approximation. From a physical point of view, $\eta$ characterizes the change
rate between different degrees of channel blockage, including LoS, obstructed
LoS, and NLoS [31]. The distance of the A2G link,
$\|\mathbf{x}_{m}-\mathbf{u}_{k}\|$, is used for normalization, which
guarantees an identical channel blockage condition along the same direction
relative to the user position. Thus, the channel parameters
$\alpha_{k}(\mathbf{x}_{m})$ and $\beta_{k}(\mathbf{x}_{m})$ in (2) can be
approximated as

$\displaystyle\left\\{\begin{matrix}\alpha_{k}(\mathbf{x}_{m})=(\alpha_{1}-\alpha_{2})\cdot
s(\mathbf{x}_{m},\mathbf{u}_{k})+\alpha_{2},\\\
\beta_{k}(\mathbf{x}_{m})=(\beta_{1}-\beta_{2})\cdot
s(\mathbf{x}_{m},\mathbf{u}_{k})+\beta_{2}.\end{matrix}\right.$ (6)

By substituting (6) into (1), we then obtain a blockage-aware A2G channel
model, which involves position-dependent channel parameters to distinguish in
different propagation environments. In particular, for $\eta\to\infty$, we
have

$\lim\limits_{{}^{\eta\to\infty}}g_{k}(\mathbf{x}_{m})=\begin{cases}\frac{\beta_{1}}{\|\mathbf{x}_{m}-\mathbf{u}_{k}\|^{\alpha_{1}}},\min\limits_{q\in\mathcal{Q}}\left\\{d_{k,q}(\mathbf{x}_{m})\right\\}\geq\Delta,\\\
\frac{\beta_{2}}{\|\mathbf{x}_{m}-\mathbf{u}_{k}\|^{\alpha_{2}}},\min\limits_{q\in\mathcal{Q}}\left\\{d_{k,q}(\mathbf{x}_{m})\right\\}\leq-\Delta,\end{cases}$
(7)

where $\Delta>0$ is an arbitrarily small positive constant. The model in (7)
is consistent with the traditional segmented channel models which classify the
propagation into LoS and NLoS conditions [22, 2].

###  II-B Problem Formulation

Denote $\mathbf{X}=\left\\{\mathbf{x}_{m},m\in\mathcal{M}\right\\}$ as the
collection of the positions of UAVs, and
$\mathbf{P}=\left\\{p_{m,n},m\in\right.$
$\left.\mathcal{M},n\in\mathcal{N}\right\\}$ as the collection of the transmit
powers of UAVs, in which $p_{m,n}$ is the transmit power allocated on
subcarrier $n$ at UAV $m$. Assuming that user $k$ is served by UAV $m$ over
subcarrier $n$, then the received signal-to-interference-plus-noise-ratio
(SINR) of user $k$ on this subcarrier can be expressed as

$\mathrm{SINR}_{k,m,n}\left(\mathbf{X},\mathbf{P}\right)=\frac{p_{m,n}g_{k}(\mathbf{x}_{m})}{\sum\limits_{j\in\mathcal{M}\setminus\\{m\\}}p_{j,n}g_{k}(\mathbf{x}_{j})+\sigma^{2}},$
(8)

where $\sigma^{2}$ is the power of the additive white Gaussian noise (AWGN) at
user $k$.

Denote
$\mathbf{C}=\left\\{c_{k,m,n},k\in\mathcal{K},m\in\mathcal{M},n\in\mathcal{N}\right\\}$
as the collection of user-UAV-subcarrier association indicators. If user $k$
is served by UAV $m$ over subcarrier $n$, we have $c_{k,m,n}=1$; otherwise,
$c_{k,m,n}=0$. Then, the achievable rate (in bits/sec/Hz) of user $k$ served
by UAV $m$ over subcarrier $n$ is expressed as a function of $\mathbf{X}$,
$\mathbf{P}$, and $\mathbf{C}$, i.e.,

$\displaystyle R_{k,m,n}\left(\mathbf{X},\mathbf{P},\mathbf{C}\right)$
$\displaystyle=c_{k,m,n}\log\left(1+\mathrm{SINR}_{k,m,n}\left(\mathbf{X},\mathbf{P}\right)\right)$
(9)
$\displaystyle=c_{k,m,n}\log\left(1+\frac{p_{m,n}g_{k}(\mathbf{x}_{m})}{\sum\limits_{j\in\mathcal{M}\setminus\\{m\\}}p_{j,n}g_{k}(\mathbf{x}_{j})+\sigma^{2}}\right),$

with the base of $\log$ equal to 2. Then, the achievable rate of user $k$ is
given by

$\displaystyle R_{k}(\mathbf{X},\mathbf{P},\mathbf{C})$
$\displaystyle=\sum_{m\in\mathcal{M}}\sum_{n\in\mathcal{N}}R_{k,m,n}\left(\mathbf{X},\mathbf{P},\mathbf{C}\right).$
(10)

To maximize the minimum achievable rate among all the ground users, we
formulate the following problem for optimizing the UAV positioning
$\mathbf{X}$, power allocation $\mathbf{P}$, and user association along with
subcarrier allocation $\mathbf{C}$:

$\displaystyle\max\limits_{\mathbf{X},\mathbf{P},\mathbf{C}}~{}~{}$
$\displaystyle\min\limits_{k\in\mathcal{K}}R_{k}(\mathbf{X},\mathbf{P},\mathbf{C})$
(11) s.t. $\displaystyle
c_{k,m,n}\in\\{0,1\\},k\in\mathcal{K},m\in\mathcal{M},n\in\mathcal{N},$ (11a)
$\displaystyle\sum_{k\in\mathcal{K}}c_{k,m,n}\leq
1,m\in\mathcal{M},n\in\mathcal{N},$ (11b)
$\displaystyle\sum_{m\in\mathcal{M}}\sum_{n\in\mathcal{N}}c_{k,m,n}=1,k\in\mathcal{K},$
(11c) $\displaystyle p_{m,n}\geq 0,m\in\mathcal{M},n\in\mathcal{N},$ (11d)
$\displaystyle\sum_{n\in\mathcal{N}}p_{m,n}\leq
P_{\mathrm{max}},m\in\mathcal{M},$ (11e)
$\displaystyle\mathbf{x}_{m}\in\mathcal{D},m\in\mathcal{M},$ (11f)
$\displaystyle\|\mathbf{x}_{m}-\mathbf{x}_{j}\|^{2}\geq
d_{\mathrm{min}}^{2},m,j\in\mathcal{M},m\neq j,$ (11g)

where $\min_{k\in\mathcal{K}}R_{k}(\mathbf{X},\mathbf{P},\mathbf{C})$ denotes
the minimum achievable rate among all the users. Constraint (11a) indicates
that $c_{k,m,n}$ is binary. Constraint (11b) ensures that on each subcarrier,
each UAV can serve one user at most. Constraint (11c) guarantees that each
user can be connected to one UAV over one subcarrier. Constraints (11d) and
(11e) indicate that the transmit power of each UAV on each subcarrier is
nonnegative and the total transmit power does not exceed a maximum value
$P_{\mathrm{max}}$. Constraint (11f) confines the region that each UAV can be
positioned, where
$\mathcal{D}=\\{\mathbf{x}\in\mathbb{R}^{3}|[\mathbf{x}]_{1}\in[0,x_{\mathrm{D}}],[\mathbf{x}]_{2}\in[0,y_{\mathrm{D}}],[\mathbf{x}]_{3}\geq
h_{\mathrm{min}}\\}$ denotes the whole considered region. Finally, in
practice, the positions of UAVs are subject to the collision avoidance
constraint (11g), where $d_{\mathrm{min}}$ denotes the minimum inter-UAV
distance to ensure collision avoidance. As can be seen, problem (11) involves
combinatorial programming variables, and the variables are highly coupled. It
is challenging to obtain the globally optimal solution for this problem. To
address this issue, we propose a suboptimal solution for problem (11) in the
following section.

##  III Proposed Solution

In this section, we propose to employ the penalty method and BSCA technique to
solve the optimization problem (11). In Section III-A, we first transform
problem (11) to a penalized problem, where the binary variables $\mathbf{C}$
are relaxed to continuous ones and a penalty term is introduced to the
objective function. Then, a double-loop optimization framework is developed.
The inner-loop solves the penalized problem for given penalty multipliers by
alternately optimizing UAV positioning $\mathbf{X}$ and resource allocation
$\\{\mathbf{P},\mathbf{C}\\}$, which are introduced in Sections III-B and
III-C, respectively. The outer-loop updates the penalty multipliers to
decrease the violation of the relaxed constraints, as detailed in Section
III-D. We finally present the overall PDLIO algorithm in Section III-E.

###  III-A Problem Transformation

To make the binary constraint in (11a) more tractable, we replace it
equivalently by

$\displaystyle 0\leq c_{k,m,n}\leq
1,k\in\mathcal{K},m\in\mathcal{M},n\in\mathcal{N},$ (12) $\displaystyle
c_{k,m,n}(1-c_{k,m,n})\leq 0,k\in\mathcal{K},m\in\mathcal{M},n\in\mathcal{N}.$
(13)

In this way, $c_{k,m,n}$ becomes a continuous optimization variable between 0
and 1. By penalizing constraint (13) into the objective function with
multipliers
$\mathbf{\Lambda}=\\{\lambda_{k,m,n},k\in\mathcal{K},m\in\mathcal{M},n\in\mathcal{N}\\}$,
we obtain a penalized problem

$\displaystyle\max\limits_{\mathbf{X},\mathbf{P},\mathbf{C}}~{}~{}$
$\displaystyle\min\limits_{k\in\mathcal{K}}R_{k}(\mathbf{X},\mathbf{P},\mathbf{C})+\rho(\mathbf{\Lambda},\mathbf{C})$
(14) s.t. $\displaystyle\lambda_{k,m,n}\geq
0,k\in\mathcal{K},m\in\mathcal{M},n\in\mathcal{N},$ (14a)
$\displaystyle\eqref{c:orthogonal},\eqref{c:onetoone},\eqref{c:power0},\eqref{c:power_max},\eqref{c:region},\eqref{c:secure_dis},\eqref{c:c_kmn_continuous},$

where
$\rho(\mathbf{\Lambda},\mathbf{C})\triangleq-\sum\limits_{k\in\mathcal{K}}\sum\limits_{m\in\mathcal{M}}\sum\limits_{n\in\mathcal{N}}\lambda_{k,m,n}c_{k,m,n}(1-c_{k,m,n})$.

Note that multiplier $\lambda_{k,m,n}$ avoids the violation of constraint
(13). In particular, for
$\lambda_{k,m,n}\to\infty,k\in\mathcal{K},m\in\mathcal{M},n\in\mathcal{N}$,
problem (14) is equivalent to the original problem (11) [32], and thus they
have the same optimal solution. However, problem (14) is still non-convex and
may be suboptimally solved during the iterations. In such case, it is not wise
to initialize $\lambda_{k,m,n}$ to be too large, since the objective will be
dominated by the penalty term $\rho(\mathbf{\Lambda},\mathbf{C})$ and the
minimum achievable rate term
$\min_{k\in\mathcal{K}}R_{k}(\mathbf{X},\mathbf{P},\mathbf{C})$ will be
diminished. Therefore, we initialize $\lambda_{k,m,n}$ as a small value to
provide enough degrees of freedom for UAV positioning and resource allocation
to obtain a good solution. Then, by gradually increasing the value of
$\lambda_{k,m,n}$, the violation of constraint (13) can be gradually decreased
until it is strictly satisfied, which ensures that the solution for penalized
problem (14) converges to a feasible solution for the original problem (11).

In the next two subsections, for the $(L+1)$-th outer loop with given
multipliers $\Lambda^{L}$, we introduce a method to solve the penalized
problem (14) through inner-loop iterations. Define

$\displaystyle\hat{R}_{k,m,n}\left(\mathbf{X},\mathbf{P},\mathbf{C}\right)$
$\displaystyle=c_{k,m,n}\log\left(1+\sum_{j\in\mathcal{M}}\frac{p_{j,n}}{\sigma^{2}}g_{k}(\mathbf{x}_{j})\right),$
$\displaystyle\bar{R}_{k,m,n}\left(\mathbf{X},\mathbf{P},\mathbf{C}\right)$
$\displaystyle=c_{k,m,n}\log\left(1+\sum_{j\in\mathcal{M}\setminus\\{m\\}}\frac{p_{j,n}}{\sigma^{2}}g_{k}(\mathbf{x}_{j})\right),$

and then (9) is equivalent to

$\displaystyle
R_{k,m,n}\left(\mathbf{X},\mathbf{P},\mathbf{C}\right)=\hat{R}_{k,m,n}\left(\mathbf{X},\mathbf{P},\mathbf{C}\right)-\bar{R}_{k,m,n}\left(\mathbf{X},\mathbf{P},\mathbf{C}\right).$

For notational simplicity, we define

$Z(\mathbf{X},\mathbf{P},\mathbf{C})=\min\limits_{k\in\mathcal{K}}R_{k}\left(\mathbf{X},\mathbf{P},\mathbf{C}\right)+\rho(\mathbf{\Lambda}^{L},\mathbf{C}),$

which is the objective of problem (14) with given $\Lambda^{L}$.

###  III-B UAV Positioning

For the $(l+1)$-th iteration of the inner-loop, given the power allocation
$\mathbf{P}^{{l}}=\\{p_{m,n}^{{l}}\\}$ and association
$\mathbf{C}^{{l}}=\\{c_{k,m,n}^{{l}}\\}$, problem (14) is transformed to the
following UAV positioning problem:

$\displaystyle\max\limits_{\mathbf{X}}~{}~{}$ $\displaystyle
Z(\mathbf{X},\mathbf{P}^{{l}},\mathbf{C}^{{l}})$ (15) s.t.
$\displaystyle\eqref{c:region},\eqref{c:secure_dis}.$

Problem (15) is a non-convex problem because of the non-concave objective
function and non-convex constraint (11g). To tackle this problem, we solve a
local approximation of problem (15) based on the BSCA technique [28, 29]. This
method mainly involves three steps as follows:

  * •

Construct a concave approximation of
$Z(\mathbf{X},\mathbf{P}^{{l}},\mathbf{C}^{{l}})$ at local point
$\mathbf{X}^{{l}}$ in the $l$-th iteration of $\mathbf{P}^{{l}}$ and
$\mathbf{C}^{{l}}$;

  * •

Solve the local approximation problem to obtain an ascent direction of the
original objective function;

  * •

Select a proper stepsize to update positioning variables $\mathbf{X}$ for
yielding an increase of the objective value.

####  III-B1 Construction of Approximation Functions

We aim to design a surrogate function
$Z^{\mathrm{Pos}}\left(\mathbf{X};\right.$
$\left.\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$, which is a
concave function approximation of the original objective function
$Z(\mathbf{X},\mathbf{P}^{{l}},$ $\mathbf{C}^{{l}})$ at local point
$\mathbf{X}^{{l}}=\\{\mathbf{x}_{m,n}^{{l}}\\}$. The surrogate function
$Z^{\mathrm{Pos}}\left(\mathbf{X};\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$
is chosen to be a strongly concave function such that problem (15) can be
relaxed to a convex problem. Moreover, it should satisfy the following
conditions to ensure the local equivalence property [28]:

$\displaystyle
Z^{\mathrm{Pos}}\left(\mathbf{X}^{{l}};\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$
$\displaystyle=Z\left(\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right),$
(16)
$\displaystyle\nabla_{\mathbf{X}}Z^{\mathrm{Pos}}\left(\mathbf{X}^{{l}};\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$
$\displaystyle=\nabla_{\mathbf{X}}Z\left(\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right).$
(17)

To this end, we first design a concave function approximation of the channel
gain $g_{k}(\mathbf{x}_{m})$. Note that in practice, it is unlikely that the
large-scale propagation environment rapidly changes over a small area.
Therefore, the propagation parameters $\alpha_{k}(\mathbf{x}_{m})$ and
$\beta_{k}(\mathbf{x}_{m})$ are approximately the same over a local region
around $\mathbf{x}_{m}^{{l}}$. Then, an approximation on the channel gain
$g_{k}(\mathbf{x}_{m})$ in (1) can be given by

$\displaystyle
g_{k}(\mathbf{x}_{m})\approx\underbrace{\tilde{g}_{k}(\mathbf{x}_{m};\mathbf{x}_{m}^{{l}})}_{\text{Predicted
channel
gain}}+\underbrace{\delta_{k}(\mathbf{x}_{m}^{{l}})\left(\mathbf{x}_{m}-\mathbf{x}_{m}^{{l}}\right)}_{\text{Correction
term}},$ (18)

where
$\tilde{g}_{k}(\mathbf{x}_{m};\mathbf{x}_{m}^{{l}})\triangleq\beta_{k}(\mathbf{x}_{m}^{{l}})\|\mathbf{x}_{m}-\mathbf{u}_{k}\|^{-\alpha_{k}(\mathbf{x}_{m}^{{l}})}$
represents the predicted channel gain between UAV $m$ located at
$\mathbf{x}_{m}$ and user $k$, based on the current channel state with fixed
channel parameters $\alpha_{k}(\mathbf{x}_{m}^{{l}})$ and
$\beta_{k}(\mathbf{x}_{m}^{{l}})$.
$\delta_{k}(\mathbf{x}_{m}^{{l}})\left(\mathbf{x}_{m}-\mathbf{x}_{m}^{{l}}\right)$
acts as a correction term to rectify the prediction errors, as well as ensures
the gradient consistency, with
$\delta_{k}(\mathbf{x}_{m}^{{l}})=\nabla_{\mathbf{x}_{m}}^{\mathrm{T}}g_{k}\left(\mathbf{x}_{m}^{{l}}\right)-\nabla_{\mathbf{x}_{m}}^{\mathrm{T}}\tilde{g}_{k}\left(\mathbf{x}_{m}^{{l}};\mathbf{x}_{m}^{{l}}\right)$.
Besides, $\nabla_{\mathbf{x}_{m}}g_{k}\left(\mathbf{x}_{m}^{{l}}\right)$ and
$\nabla_{\mathbf{x}_{m}}\tilde{g}_{k}\left(\mathbf{x}_{m}^{{l}};\mathbf{x}_{m}^{{l}}\right)$
can be derived as

$\displaystyle\nabla_{\mathbf{x}_{m}}g_{k}\left(\mathbf{x}_{m}^{{l}}\right)=$
$\displaystyle-
g_{k}\left(\mathbf{x}_{m}^{{l}}\right)\cdot\nabla_{\mathbf{x}_{m}}\alpha_{k}\left(\mathbf{x}_{m}^{{l}}\right)\cdot\log(\|\mathbf{x}_{m}^{{l}}-\mathbf{u}_{k}\|)$
$\displaystyle-
g_{k}\left(\mathbf{x}_{m}^{{l}}\right)\cdot\frac{\alpha_{k}\left(\mathbf{x}_{m}^{{l}}\right)\cdot(\mathbf{x}_{m}^{{l}}-\mathbf{u}_{k})}{\|\mathbf{x}_{m}^{{l}}-\mathbf{u}_{k}\|^{2}}+\frac{\nabla_{\mathbf{x}_{m}}\beta_{k}\left(\mathbf{x}_{m}^{{l}}\right)}{\|\mathbf{x}_{m}^{{l}}-\mathbf{u}_{k}\|^{\alpha_{k}\left(\mathbf{x}_{m}^{{l}}\right)}},$
$\displaystyle\nabla_{\mathbf{x}_{m}}\tilde{g}_{k}\left(\mathbf{x}_{m}^{{l}};\mathbf{x}_{m}^{{l}}\right)=-\tilde{g}_{k}\left(\mathbf{x}_{m}^{{l}};\mathbf{x}_{m}^{{l}}\right)\cdot\frac{\alpha_{k}(\mathbf{x}_{m}^{{l}})\cdot(\mathbf{x}_{m}^{{l}}-\mathbf{u}_{k})}{\|\mathbf{x}_{m}^{{l}}-\mathbf{u}_{k}\|^{2}}.$

The detailed procedures to derive
$\nabla_{\mathbf{x}_{m}}g_{k}\left(\mathbf{x}_{m}\right)$,
$\nabla_{\mathbf{x}_{m}}\alpha_{k}\left(\mathbf{x}_{m}\right)$, and
$\nabla_{\mathbf{x}_{m}}\beta_{k}\left(\mathbf{x}_{m}\right)$ are provided in
Appendix A.

Note that $\tilde{g}_{k}(\mathbf{x}_{m};\mathbf{x}_{m}^{{l}})$ is still not
concave w.r.t. $\mathbf{x}_{m}$. By taking the first-order Taylor expansion of
$\tilde{g}_{k}(\mathbf{x}_{m};\mathbf{x}_{m}^{{l}})$ w.r.t.
$\|\mathbf{x}_{m}-\mathbf{u}_{k}\|^{2}$, we can obtain a concave function
approximation of $\tilde{g}_{k}(\mathbf{x}_{m};\mathbf{x}_{m}^{{l}})$ as

$\displaystyle\tilde{g}_{k}(\mathbf{x}_{m};\mathbf{x}_{m}^{{l}})$
$\displaystyle\approx
A_{k,m}^{{l}}\left(\|\mathbf{x}_{m}^{{l}}-\mathbf{u}_{k}\|^{2}-\|\mathbf{x}_{m}-\mathbf{u}_{k}\|^{2}\right)+\tilde{g}_{k}(\mathbf{x}_{m}^{{l}};\mathbf{x}_{m}^{{l}}),$
(19)

with
$A_{k,m}^{{l}}=\alpha_{k}(\mathbf{x}_{m}^{{l}})\beta_{k}(\mathbf{x}_{m}^{{l}})/\left(2\|\mathbf{x}_{m}^{{l}}-\mathbf{u}_{k}\|^{2+\alpha_{k}(\mathbf{x}_{m}^{{l}})}\right)$.

Based on the above approximation of the channel gain, a concave function
approximation of
$\hat{R}_{k,m,n}\left(\mathbf{X},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$
around $\mathbf{X}^{{l}}$ can be obtained as

$\displaystyle\hat{R}_{k,m,n}\left(\mathbf{X},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$
(20)
$\displaystyle\overset{(a)}{\approx}B_{k,m,n}^{{l}}\sum_{j\in\mathcal{M}}\frac{p_{j,n}^{{l}}}{\sigma^{2}}\left(g_{k}(\mathbf{x}_{j})-g_{k}(\mathbf{x}_{j}^{{l}})\right)+\hat{R}_{k,m,n}\left(\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$
$\displaystyle\overset{(b)}{\approx}B_{k,m,n}^{{l}}\sum_{j\in\mathcal{M}}\frac{p_{j,n}^{{l}}}{\sigma^{2}}\left(A_{k,j}^{{l}}\left(\|\mathbf{x}_{j}^{{l}}-\mathbf{u}_{k}\|^{2}-\|\mathbf{x}_{j}-\mathbf{u}_{k}\|^{2}\right)+\delta_{k}(\mathbf{x}_{j}^{{l}})\left(\mathbf{x}_{j}-\mathbf{x}_{j}^{{l}}\right)\right)+\hat{R}_{k,m,n}\left(\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$
$\displaystyle\triangleq\hat{R}_{k,m,n}^{\mathrm{appr}}\left(\mathbf{X};\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$

with
$B_{k,m,n}^{{l}}=c_{k,m,n}^{{l}}/\left(1+\sum_{j\in\mathcal{M}}\frac{p_{j,n}^{{l}}}{\sigma^{2}}g_{k}(\mathbf{x}_{j}^{{l}})\right)$.
The approximation in step $(a)$ is obtained according to the first-order
Taylor expansion of $\log(1+y)$, and the approximation in step $(b)$ is
obtained by replacing $g_{k}(\mathbf{x}_{j})$ with its approximation in (18),
in which $\tilde{g}_{k}(\mathbf{x}_{j};\mathbf{x}_{j}^{{l}})$ is replaced by
its concave function approximation in (19).

On the other hand, a linear function approximation of
$\bar{R}_{k,m,n}\left(\mathbf{X},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$ can
be obtained as

$\displaystyle\bar{R}_{k,m,n}\left(\mathbf{X},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$
$\displaystyle\overset{(c)}{\approx}B_{k,-m,n}^{{l}}\sum_{j\in\mathcal{M}\setminus\\{m\\}}\frac{p_{j,n}^{{l}}}{\sigma^{2}}\left(g_{k}(\mathbf{x}_{j})-g_{k}(\mathbf{x}_{j}^{{l}})\right)+\bar{R}_{k,m,n}\left(\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$
(21)
$\displaystyle\overset{(d)}{\approx}B_{k,-m,n}^{{l}}\sum_{j\in\mathcal{M}\setminus\\{m\\}}\frac{p_{j,n}^{{l}}}{\sigma^{2}}\cdot\nabla_{\mathbf{x}_{j}}^{\mathrm{T}}g_{k}\left(\mathbf{x}_{j}^{{l}}\right)\cdot\left(\mathbf{x}_{j}-\mathbf{x}_{j}^{{l}}\right)+\bar{R}_{k,m,n}\left(\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$
$\displaystyle\triangleq\bar{R}_{k,m,n}^{\mathrm{appr}}\left(\mathbf{X};\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$

with
$B_{k,-m,n}^{{l}}=c_{k,m,n}^{{l}}/\left(1+\sum_{j\in\mathcal{M}\setminus\\{m\\}}\frac{p_{j,n}^{{l}}}{\sigma^{2}}g_{k}(\mathbf{x}_{j}^{{l}})\right)$.
The approximation in step $(c)$ is obtained according to the first-order
Taylor expansion of $\log(1+y)$, and the approximation in step $(d)$ is
obtained by replacing $g_{k}(\mathbf{x}_{j})$ with its first-order Taylor
expansion.

Therefore, we can obtain a concave function approximation of
$R_{k}\left(\mathbf{X},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$ as

$\displaystyle
R_{k}^{\mathrm{Pos}}\left(\mathbf{X};\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$
$\displaystyle=\sum_{m\in\mathcal{M}}\sum_{n\in\mathcal{N}}\hat{R}_{k,m,n}^{\mathrm{appr}}\left(\mathbf{X};\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)-\sum_{m\in\mathcal{M}}\sum_{n\in\mathcal{N}}\bar{R}_{k,m,n}^{\mathrm{appr}}\left(\mathbf{X};\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right).$

We then obtain an approximation function of
$Z(\mathbf{X},\mathbf{P}^{{l}},\mathbf{C}^{{l}})$ as

$\displaystyle
Z^{\mathrm{Pos}}\left(\mathbf{X};\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)=\min\limits_{k\in\mathcal{K}}R_{k}^{\mathrm{Pos}}\left(\mathbf{X};\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)+\rho(\mathbf{\Lambda}^{L},\mathbf{C}^{l})$
(22)

From (20) and (21), it is not hard to verify that
$Z^{\mathrm{Pos}}\left(\mathbf{X};\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$
is concave in terms of $\mathbf{X}$. Also, (16) and (17) hold because the
functions in (20) and (21) are from the first-order Taylor approximations of
that in $Z(\mathbf{X},\mathbf{P}^{l},\mathbf{C}^{l})$ in terms of
$\mathbf{X}$.

Finally, the left-hand-side (LHS) of constraint (11g) is lower-bounded by its
first-order Taylor expansion w.r.t. $\|\mathbf{x}_{m}-\mathbf{x}_{j}\|$, i.e.,
$\|\mathbf{x}_{m}-\mathbf{x}_{j}\|^{2}\geq
2\|\mathbf{x}_{m}^{l}-\mathbf{x}_{j}^{l}\|^{\mathrm{T}}\|\mathbf{x}_{m}-\mathbf{x}_{j}\|-\|\mathbf{x}_{m}^{l}-\mathbf{x}_{j}^{l}\|^{2}.$
Then, constraint (11g) can be relaxed as the following convex constraint

$\displaystyle
2\|\mathbf{x}_{m}^{l}-\mathbf{x}_{j}^{l}\|^{\mathrm{T}}\|\mathbf{x}_{m}-\mathbf{x}_{j}\|-\|\mathbf{x}_{m}^{l}-\mathbf{x}_{j}^{l}\|^{2}\geq
d_{\mathrm{min}}^{2},m,j\in\mathcal{M},m\neq j.$ (23)

####  III-B2 Solution for Relaxed Problem

Given a local point $\mathbf{X}^{l}$ in the $(l+1)$-th iteration, by replacing
$Z\left(\mathbf{X},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$ with
$Z^{\mathrm{Pos}}\left(\mathbf{X};\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$
and replacing constraint (11g) with (23), problem (15) is relaxed as

$\displaystyle\max\limits_{\mathbf{X}}~{}~{}$ $\displaystyle
Z^{\mathrm{Pos}}\left(\mathbf{X};\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$
(24) s.t. $\displaystyle\eqref{c:region},\eqref{c:secure_dis_2}.$

Problem (24) is to maximize a concave objective function subject to convex
constraints on $\mathbf{X}$. Therefore, it is a convex problem and can be
readily solved by standard convex program solvers, such as CVXPY [33]. We
denote the optimal solution for problem (24) as
$\tilde{\mathbf{X}}^{l\star}=\\{\tilde{\mathbf{x}}_{m}^{{l\star}}\\}$.

####  III-B3 UAV Positioning Update

Since $\tilde{\mathbf{X}}^{l\star}$ is an optimal solution for problem (24),
we have

$\displaystyle 0$
$\displaystyle\overset{(e)}{\leq}Z^{\mathrm{Pos}}\left(\tilde{\mathbf{X}}^{l\star};\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)-Z^{\mathrm{Pos}}\left(\mathbf{X}^{{l}};\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$
(25)
$\displaystyle\overset{(f)}{\leq}(\tilde{\mathbf{X}}^{l\star}-\mathbf{X}^{{l}})^{\mathrm{T}}\nabla_{\mathbf{X}}Z^{\mathrm{Pos}}\left(\mathbf{X}^{{l}};\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)\overset{(g)}{=}(\tilde{\mathbf{X}}^{l\star}-\mathbf{X}^{{l}})^{\mathrm{T}}\nabla_{\mathbf{X}}Z\left(\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right),$

where $(e)$, $(f)$, and $(g)$ hold, respectively, due to the optimality of
$\tilde{\mathbf{X}}^{l\star}$, the concavity of
$Z^{\mathrm{Pos}}\left(\mathbf{X};\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$
w.r.t. $\mathbf{X}$, and the consistency conditions (16) and (17). Therefore,
$\tilde{\mathbf{X}}^{l\star}-\mathbf{X}^{{l}}$ is an ascent direction of
$Z\left(\mathbf{X},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$ at
$\mathbf{X}=\mathbf{X}^{l}$, along which the objective value can be further
increased compared with that at the current local point $\mathbf{X}^{l}$.

Given the ascent direction $\tilde{\mathbf{X}}^{l\star}-\mathbf{X}^{{l}}$, the
UAV positioning variables are updated as

$\displaystyle\mathbf{X}^{l+1}=\mathbf{X}^{{l}}+\gamma_{1}^{l}(\tilde{\mathbf{X}}^{l\star}-\mathbf{X}^{{l}}),$
(26)

where $\gamma_{1}^{l}$ is the stepsize that needs to be selected properly.
Backtracking line search can be adopted to choose $\gamma_{1}^{l}$ to
efficiently increase
$Z\left(\mathbf{X},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$ and yield fast
convergence [34]. Given constants $\zeta,\tau\in(0,1)$, $\gamma_{1}^{l}$ is
set to be $\gamma_{1}^{l}=\zeta^{t_{l}}$ with $t_{l}$ being the smallest
nonnegative integer that ensures constraint (11g), i.e.,
$\|\mathbf{x}_{m}^{l+1}-\mathbf{x}_{j}^{l+1}\|^{2}\geq
d_{\mathrm{min}}^{2},m,j\in\mathcal{M},m\neq j$, as well as the following
inequality:222Constant $\tau$ can be interpreted as the fraction of the
increment in $Z(\mathbf{X},\mathbf{P}^{l},\mathbf{C}^{l})$ predicted by linear
extrapolation that we will accept. Starting from $\gamma_{1}^{l}=\zeta^{0}=1$
and gradually decreasing $\gamma_{1}^{l}$, there always exists a
$\gamma_{1}^{l}$ (at least we have $\gamma_{1}^{l}=\zeta^{\infty}=0$)
satisfying both (11g) and (27). For numerical evaluation, if $\gamma_{1}^{l}$
falls below a small positive threshold, we set $\gamma_{1}^{l}=0$, and update
UAV positioning variables as $\mathbf{X}^{l+1}=\mathbf{X}^{{l}}$.

$\displaystyle
Z\left(\mathbf{X}^{{l+1}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)-Z\left(\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)\geq\tau\gamma_{1}^{l}(\tilde{\mathbf{X}}^{l\star}-\mathbf{X}^{{l}})^{\mathrm{T}}\nabla_{\mathbf{X}}Z\left(\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right),$
(27)

with

$\displaystyle(\tilde{\mathbf{X}}^{l\star}-\mathbf{X}^{{l}})^{\mathrm{T}}\nabla_{\mathbf{X}}Z\left(\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)=\sum_{m\in\mathcal{M}}\sum_{n\in\mathcal{N}}(\tilde{\mathbf{X}}^{l\star}-\mathbf{X}^{{l}})^{\mathrm{T}}\nabla_{\mathbf{X}}R_{k^{\prime},m,n}\left(\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$
$\displaystyle=\sum_{m\in\mathcal{M}}\sum_{n\in\mathcal{N}}\left(B_{k^{\prime},m,n}^{{l}}\sum_{j\in\mathcal{M}}\frac{p_{j,n}^{{l}}}{\sigma^{2}}\left(\tilde{\mathbf{x}}_{j}^{{l\star}}-\mathbf{x}_{j}^{{l}}\right)^{\mathrm{T}}\nabla_{\mathbf{x}_{j}}g_{k^{\prime}}\left(\mathbf{x}_{j}^{{l}}\right)\right.$
$\displaystyle\left.~{}~{}~{}~{}-B_{k^{\prime},-m,n}^{{l}}\sum_{j\in\mathcal{M}\setminus\\{m\\}}\frac{p_{j,n}^{{l}}}{\sigma^{2}}\left(\tilde{\mathbf{x}}_{j}^{{l\star}}-\mathbf{x}_{j}^{{l}}\right)^{\mathrm{T}}\nabla_{\mathbf{x}_{j}}g_{k^{\prime}}\left(\mathbf{x}_{j}^{{l}}\right)\right),$

and
$k^{\prime}=\arg\min\limits_{k}R_{k}\left(\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$
such that
$Z\left(\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)=R_{k^{\prime}}\left(\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)+\rho(\mathbf{\Lambda}^{L},\mathbf{C}^{l})$.
The corresponding objective value of problem (14) is updated as
$Z(\mathbf{X}^{l+1},\mathbf{P}^{l},\mathbf{C}^{l})$.

###  III-C Resource Allocation

Given UAVs’ positions $\mathbf{X}^{l+1}=\\{\mathbf{x}_{m}^{{l+1}}\\}$, problem
(11) is transformed to the following resource allocation problem

$\displaystyle\max\limits_{\mathbf{P},\mathbf{C}}~{}$ $\displaystyle
Z(\mathbf{X}^{{l+1}},\mathbf{P},\mathbf{C})$ (28) s.t.
$\displaystyle\eqref{c:orthogonal},\eqref{c:onetoone},\eqref{c:power0},\eqref{c:power_max},\eqref{c:c_kmn_continuous}.$

Problem (28) is a non-convex problem because of the non-concave objective
function. Similar to the method for optimizing UAV positioning, we aim to
design a concave function approximation for
$Z(\mathbf{X}^{{l+1}},\mathbf{P},\mathbf{C})$, solve the relaxed counterpart
of problem (28), and update the resource allocation variables.

####  III-C1 Construction of Approximation Functions

We aim to design a surrogate function
$Z^{\mathrm{RA}}\left(\mathbf{P},\mathbf{C};\right.$
$\left.\mathbf{X}^{{l+1}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$, which is
a concave function approximation of
$Z(\mathbf{X}^{{l+1}},\mathbf{P},\mathbf{C})$ at local point
$(\mathbf{P}^{{l}},\mathbf{C}^{{l}})$. Besides, the following conditions
should be satisfied

$\displaystyle
Z^{\mathrm{RA}}\left(\mathbf{P}^{{l}},\mathbf{C}^{{l}};\mathbf{X}^{{l+1}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)=Z(\mathbf{X}^{{l+1}},\mathbf{P}^{l},\mathbf{C}^{l}),$
(29)
$\displaystyle\nabla_{\mathbf{P}}Z^{\mathrm{RA}}\left(\mathbf{P}^{{l}},\mathbf{C}^{{l}};\mathbf{X}^{{l+1}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)=\nabla_{\mathbf{P}}Z\left(\mathbf{X}^{{l+1}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right),$
(30)
$\displaystyle\nabla_{\mathbf{C}}Z^{\mathrm{RA}}\left(\mathbf{P}^{{l}},\mathbf{C}^{{l}};\mathbf{X}^{{l+1}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)=\nabla_{\mathbf{C}}Z\left(\mathbf{X}^{{l+1}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right).$
(31)

To this end, we adopt the approximation of
$R_{k,m,n}\left(\mathbf{X}^{{l+1}},\mathbf{P},\mathbf{C}\right)$ in the
$(l+1)$-th iteration as

$\displaystyle
R_{k,m,n}\left(\mathbf{X}^{{l+1}},\mathbf{P},\mathbf{C}\right)\approx
R_{k,m,n}\left(\mathbf{X}^{{l+1}},\mathbf{P},\mathbf{C}^{l}\right)+R_{k,m,n}\left(\mathbf{X}^{{l+1}},\mathbf{P}^{l},\mathbf{C}\right)-R_{k,m,n}\left(\mathbf{X}^{{l+1}},\mathbf{P}^{l},\mathbf{C}^{l}\right).$
(32)

Note that
$R_{k,m,n}\left(\mathbf{X}^{{l+1}},\mathbf{P},\mathbf{C}^{l}\right)=\hat{R}_{k,m,n}\left(\mathbf{X}^{{l+1}},\mathbf{P},\mathbf{C}^{l}\right)-\bar{R}_{k,m,n}\left(\mathbf{X}^{{l+1}},\mathbf{P},\mathbf{C}^{l}\right)$
is the difference of two concave functions w.r.t. $\mathbf{P}$. Since any
concave function is globally upper-bounded by its first-order Taylor expansion
at any point, we have

$\displaystyle\bar{R}_{k,m,n}\left(\mathbf{X}^{{l+1}},\mathbf{P},\mathbf{C}^{l}\right)=c_{k,m,n}^{l}\log\left(1+\sum_{j\in\mathcal{M}\setminus\\{m\\}}\frac{g_{k}(\mathbf{x}_{j}^{{l+1}})}{\sigma^{2}}p_{j,n}\right)$
(33) $\displaystyle\leq
B_{k,-m,n}^{{}^{\prime}{l}}\sum_{j\in\mathcal{M}/\\{m\\}}\frac{g_{k}(\mathbf{x}_{j}^{{l+1}})}{\sigma^{2}}\left(p_{j,n}-p_{j,n}^{{l}}\right)+\bar{R}_{k,m,n}\left(\mathbf{X}^{{l+1}},\mathbf{P}^{l},\mathbf{C}^{l}\right)$
$\displaystyle\triangleq\bar{R}_{k,m,n}^{\mathrm{ub}}\left(\mathbf{P};\mathbf{X}^{{l+1}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right),$

with
$B_{k,-m,n}^{{}^{\prime}{l}}=c_{k,m,n}^{{l}}/\left(1+\sum_{j\in\mathcal{M}\setminus\\{m\\}}\frac{p_{j,n}^{{l}}}{\sigma^{2}}g_{k}(\mathbf{x}_{j}^{{l+1}})\right)$.
Therefore,
$R_{k,m,n}\left(\mathbf{X}^{{l+1}},\mathbf{P},\mathbf{C}^{l}\right)$ is lower-
bounded by

$\displaystyle
R_{k,m,n}\left(\mathbf{X}^{{l+1}},\mathbf{P},\mathbf{C}^{l}\right)\geq\hat{R}_{k,m,n}\left(\mathbf{X}^{{l+1}},\mathbf{P},\mathbf{C}^{l}\right)-\bar{R}_{k,m,n}^{\mathrm{ub}}\left(\mathbf{P};\mathbf{X}^{{l+1}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right).$
(34)

By substituting (34) into (32),
$R_{k,m,n}\left(\mathbf{X}^{{l+1}},\mathbf{P},\mathbf{C}\right)$ is further
approximated as

$\displaystyle R_{k,m,n}\left(\mathbf{X}^{{l+1}},\mathbf{P},\mathbf{C}\right)$
$\displaystyle\approx\hat{R}_{k,m,n}\left(\mathbf{X}^{{l+1}},\mathbf{P},\mathbf{C}^{l}\right)-\bar{R}_{k,m,n}^{\mathrm{ub}}\left(\mathbf{P};\mathbf{X}^{{l+1}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$
(35)
$\displaystyle~{}~{}~{}+R_{k,m,n}\left(\mathbf{X}^{{l+1}},\mathbf{P}^{l},\mathbf{C}\right)-R_{k,m,n}\left(\mathbf{X}^{{l+1}},\mathbf{P}^{l},\mathbf{C}^{l}\right)$
$\displaystyle\triangleq
R_{k,m,n}^{\mathrm{appr}}\left(\mathbf{P},\mathbf{C};\mathbf{X}^{{l+1}},\mathbf{P}^{l},\mathbf{C}^{l}\right).$

In addition, since $c_{k,m,n}^{2}\geq
2c_{k,m,n}^{l}c_{k,m,n}-(c_{k,m,n}^{l})^{2}$ for a given local point
$c_{k,m,n}^{l}$, we have the following lower bound on
$\rho(\mathbf{\Lambda}^{L},\mathbf{C})$:

$\displaystyle\rho(\mathbf{\Lambda}^{L},\mathbf{C})$
$\displaystyle=-\sum_{k\in\mathcal{K}}\sum_{m\in\mathcal{M}}\sum_{n\in\mathcal{N}}\lambda_{k,m,n}^{L}c_{k,m,n}(1-c_{k,m,n})$
(36)
$\displaystyle\geq\sum_{k\in\mathcal{K}}\sum_{m\in\mathcal{M}}\sum_{n\in\mathcal{N}}\lambda_{k,m,n}^{L}\left(\left(2c_{k,m,n}^{l}-1\right)c_{k,m,n}-\left(c_{k,m,n}^{l}\right)^{2}\right)$
$\displaystyle\triangleq\rho^{\mathrm{lb}}(\mathbf{C};\mathbf{\Lambda}^{L},\mathbf{C}^{l}).$

According to (35) and (36), we obtain an approximation function of
$Z\left(\mathbf{X}^{{l+1}},\mathbf{P},\mathbf{C}\right)$ as

$\displaystyle
Z^{\mathrm{RA}}\left(\mathbf{P},\mathbf{C};\mathbf{X}^{{l+1}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)=\min\limits_{k\in\mathcal{K}}\left(\sum_{m\in\mathcal{M}}\sum_{n\in\mathcal{N}}R_{k,m,n}^{\mathrm{appr}}\left(\mathbf{P},\mathbf{C};\mathbf{X}^{{l+1}},\mathbf{P}^{l},\mathbf{C}^{l}\right)\right)+\rho^{\mathrm{lb}}(\mathbf{C};\mathbf{\Lambda}^{L},\mathbf{C}^{l}).$
(37)

It can be verified that conditions (29), (30), and (31) are satisfied.

####  III-C2 Solution for Relaxed Problem

Problem (28) can be relaxed as the following convex problem

$\displaystyle\max\limits_{\mathbf{P},\mathbf{C}}~{}$ $\displaystyle
Z^{\mathrm{RA}}\left(\mathbf{P},\mathbf{C};\mathbf{X}^{{l+1}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$
(38) s.t.
$\displaystyle\eqref{c:orthogonal},\eqref{c:onetoone},\eqref{c:power0},\eqref{c:power_max},\eqref{c:c_kmn_continuous},$

which can be readily solved by standard convex program solvers. We denote the
optimal solution for problem (38) as
$(\tilde{\mathbf{P}}^{l\star}=\\{\tilde{p}_{m,n}^{l\star}\\},\tilde{\mathbf{C}}^{l\star}=\\{\tilde{c}_{k,m,n}^{l\star}\\})$,
which is an ascent direction of
$Z\left(\mathbf{X}^{l+1},\mathbf{P},\mathbf{C}\right)$ at
$(\mathbf{P}^{l},\mathbf{C}^{l})$.

####  III-C3 Resource Allocation Update

We update the resource allocation variables as

$\displaystyle\left\\{\begin{matrix}\mathbf{P}^{l+1}=\mathbf{P}^{l}+\gamma_{2}^{l}(\tilde{\mathbf{P}}^{l\star}-\mathbf{P}^{l}),\\\
\mathbf{C}^{l+1}=\mathbf{C}^{l}+\gamma_{2}^{l}(\tilde{\mathbf{C}}^{l\star}-\mathbf{C}^{l}),\end{matrix}\right.$
(39)

where $\gamma_{2}^{l}$ is the stepsize determined by adopting the backtracking
line search, with the following condition satisfied:

$\displaystyle
Z\left(\mathbf{X}^{{l+1}},\mathbf{P}^{{l+1}},\mathbf{C}^{{l+1}}\right)-Z\left(\mathbf{X}^{{l+1}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)\geq\tau\gamma_{2}^{l}\begin{bmatrix}(\tilde{\mathbf{P}}^{l\star}-\mathbf{P}^{{l}})^{\mathrm{T}}\\\
(\tilde{\mathbf{C}}^{l\star}-\mathbf{C}^{{l}})^{\mathrm{T}}\end{bmatrix}^{\mathrm{T}}\begin{bmatrix}\nabla_{\mathbf{P}}Z\left(\mathbf{X}^{{l+1}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)\\\
\nabla_{\mathbf{C}}Z\left(\mathbf{X}^{{l+1}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)\end{bmatrix},$
(40)

with

$\displaystyle(\tilde{\mathbf{P}}^{l\star}-\mathbf{P}^{{l}})^{\mathrm{T}}\nabla_{\mathbf{P}}Z\left(\mathbf{X}^{{l+1}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)=\sum_{m\in\mathcal{M}}\sum_{n\in\mathcal{N}}(\tilde{\mathbf{P}}^{l\star}-\mathbf{P}^{{l}})^{\mathrm{T}}\nabla_{\mathbf{P}}R_{k^{\prime},m,n}\left(\mathbf{X}^{{l+1}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$
$\displaystyle=\sum_{m\in\mathcal{M}}\sum_{n\in\mathcal{N}}\left(B_{k^{\prime},m,n}^{{}^{\prime}{l}}\sum_{j\in\mathcal{M}}\frac{g_{k^{\prime}}(\mathbf{x}_{j}^{{l+1}})}{\sigma^{2}}\left(\tilde{p}_{j,n}^{{l\star}}-p_{j,n}^{{l}}\right)-B_{k^{\prime},-m,n}^{{}^{\prime}{l}}\sum_{j\in\mathcal{M}/\\{m\\}}\frac{g_{k^{\prime}}(\mathbf{x}_{j}^{{l+1}})}{\sigma^{2}}\left(\tilde{p}_{j,n}^{{l\star}}-p_{j,n}^{{l}}\right)\right),$
$\displaystyle(\tilde{\mathbf{C}}^{l\star}-\mathbf{C}^{{l}})^{\mathrm{T}}\nabla_{\mathbf{C}}Z\left(\mathbf{X}^{{l+1}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)=$
$\displaystyle\sum_{m\in\mathcal{M}}\sum_{n\in\mathcal{N}}(\tilde{c}_{k^{\prime},m,n}^{l\star}-c_{k^{\prime},m,n}^{l})\frac{R_{k^{\prime},m,n}\left(\mathbf{X}^{{l+1}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)}{c_{k^{\prime},m,n}^{l}}$
$\displaystyle-\sum_{m\in\mathcal{M}}\sum_{n\in\mathcal{N}}\sum\limits_{k\in\mathcal{K}}(\tilde{c}_{k,m,n}^{l\star}-c_{k,m,n}^{l})\lambda_{k,m,n}^{L}(1-2c_{k,m,n}^{l}),$

and
$k^{\prime}=\arg\min\limits_{k}R_{k}\left(\mathbf{X}^{{l+1}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)$,
$B_{k^{\prime},m,n}^{{}^{\prime}{l}}=c_{k^{\prime},m,n}^{{l}}/\left(1+\sum_{j\in\mathcal{M}}\frac{p_{j,n}^{{l}}}{\sigma^{2}}g_{k^{\prime}}\left(\mathbf{x}_{j}^{{l+1}}\right)\right)$.
The corresponding objective value of problem (14) is updated as
$Z(\mathbf{X}^{l+1},\mathbf{P}^{l+1},\mathbf{C}^{l+1})$.

###  III-D Penalty Multiplier Update

For the $(L+1)$-th outer-loop with given $\Lambda^{L}$, we have obtained a
suboptimal solution for problem (14) by inner-loop iterations (as described in
Sections III-B and III-C), denoted by
$\bar{\mathbf{X}}^{L+1}=\\{\bar{\mathbf{x}}_{m}^{L+1}\\}$,
$\bar{\mathbf{P}}^{L+1}=\\{\bar{p}_{m,n}^{L+1}\\}$, and
$\bar{\mathbf{C}}^{L+1}=\\{\bar{c}_{k,m,n}^{L+1}\\}$. To decrease the
violation of constraint (13), the multipliers in $\mathbf{\Lambda}$ are
updated by using the following strategy [35]

$\displaystyle\lambda_{k,m,n}^{L+1}=\lambda_{k,m,n}^{L}+\gamma^{L}\bar{c}_{k,m,n}^{L+1}(1-\bar{c}_{k,m,n}^{L+1}),$
(41)

which basically optimizes a dual problem w.r.t. $\mathbf{\Lambda}$ via the
gradient method [32]. $\gamma^{L}$ is the stepsize given by

$\displaystyle\gamma^{L}=\frac{\mu^{L}}{\sum\limits_{k\in\mathcal{K}}\sum\limits_{m\in\mathcal{M}}\sum\limits_{n\in\mathcal{N}}\left(\bar{c}_{k,m,n}^{L+1}(1-\bar{c}_{k,m,n}^{L+1})\right)^{2}},$
(42)

where $\mu^{L}$ is an adaption parameter which is set as $\mu^{0}=2$ and
$\mu^{L}\leftarrow 2\mu^{L}$ for the case when the maximum constraint
violation $\max\limits_{k,m,n}\bar{c}_{k,m,n}^{L}(1-\bar{c}_{k,m,n}^{L})$ does
not decrease in the $L$-th iteration.

It can be seen that (41) guarantees an increasing update for
$\lambda_{k,m,n}$. By gradually increasing the value of $\lambda_{k,m,n}$, the
violation of constraint (13) can decrease. For sufficiently large
$\lambda_{k,m,n}$, constraint (13) is strictly satisfied such that the
obtained solution for problem (14) is also a feasible solution for the
original problem (11).

###  III-E Overall Solution

Hereto, we are ready to show the overall PDLIO algorithm for solving the joint
positioning and resource allocation problem (11) for a multi-UAV OFDMA
communication network aided by geographic information. As summarized in
Algorithm 1, in line 1, we invoke [8, Algorithm 1] to calculate the blocked
regions $\\{\mathcal{D}_{k,q},k\in\mathcal{K},q\in\mathcal{Q}\\}$ for $K$
users caused by $Q$ buildings based on geographic information. Then, in lines
3-13, we employ the double-loop framework to iteratively solve the penalized
problem (14) and update the multipliers $\mathbf{\Lambda}$. Lines 5-9 solve
the penalized problem (14) with given $\mathbf{\Lambda}^{L}$, where the UAV
positioning and resource allocation are optimized in an alternate manner. The
inner-loop terminates when the increase of the objective value of problem (14)
from one iteration to the next falls bellow a positive threshold
$\epsilon_{\mathrm{l}}$. The outer-loop terminates when the maximum constraint
violation $\max\limits_{k,m,n}\bar{c}_{k,m,n}^{L}(1-\bar{c}_{k,m,n}^{L})$
falls below a positive threshold $\epsilon_{\mathrm{L}}$.

0: The coordinates of the vertices of the buildings,
$\\{\mathbf{x}_{k},k\in\mathcal{K}\\}$, $P_{\mathrm{max}}$, $\sigma^{2}$,
$\alpha_{1}$, $\beta_{1}$, $\alpha_{2}$, $\beta_{2}$, $\eta$,
$x_{\mathrm{D}}$, $y_{\mathrm{D}}$, $h_{\mathrm{min}}$, $d_{\mathrm{min}}$,
$\zeta$, $\tau$, $\epsilon_{\mathrm{L}}$, $\epsilon_{\mathrm{l}}$.  1:
Calculate the blocked regions
$\\{\mathcal{D}_{k,q},k\in\mathcal{K},q\in\mathcal{Q}\\}$ according to [8,
Algorithm 1].  2: Set the iteration index of outer-loops as $L=0$, and
initialize $\mathbf{\Lambda}^{0}$, $\bar{\mathbf{X}}^{0}$,
$\bar{\mathbf{P}}^{0}$, $\bar{\mathbf{C}}^{0}$.  3: repeat 4: Set the
iteration index of inner-loops as $l=0$, and initialize
$\mathbf{X}^{0}\leftarrow\bar{\mathbf{X}}^{L}$,
$\mathbf{P}^{0}\leftarrow\bar{\mathbf{P}}^{L}$,
$\mathbf{C}^{0}\leftarrow\bar{\mathbf{C}}^{L}$.  5: repeat 6: $\bullet$ Solve
problem (24) to obtain the optimal solution $\tilde{\mathbf{X}}^{l\star}$, and
update the UAV position $\mathbf{X}^{l+1}$ according to (26).  7: $\bullet$
Solve problem (38) to obtain the optimal solution
$(\tilde{\mathbf{P}}^{l\star},\tilde{\mathbf{C}}^{l\star})$, and update the
power allocation $\mathbf{P}^{l+1}$ and association $\mathbf{C}^{l+1}$
according to (39).  8: $\bullet$ Update $l\leftarrow l+1$.  9: until The
increase of objective value falls bellow $\epsilon_{\mathrm{l}}$.  10: Update
$\bar{\mathbf{X}}^{L+1}\leftarrow\mathbf{X}^{l}$,
$\bar{\mathbf{P}}^{L+1}\leftarrow\mathbf{P}^{l}$,
$\bar{\mathbf{C}}^{L+1}\leftarrow\mathbf{C}^{l}$.  11: Update multiplier
$\mathbf{\Lambda}^{L+1}$ based on (41).  12: Update $L\leftarrow L+1$.  13:
until The maximum constraint violation
$\max\limits_{k,m,n}\bar{c}_{k,m,n}^{L}(1-\bar{c}_{k,m,n}^{L})$ falls below
$\epsilon_{\mathrm{L}}$.  13: $\bar{\mathbf{X}}^{L}$, $\bar{\mathbf{P}}^{L}$,
$\bar{\mathbf{C}}^{L}$,
$Z\left(\bar{\mathbf{X}}^{L},\bar{\mathbf{P}}^{L},\bar{\mathbf{C}}^{L}\right)$.

Algorithm 1 Joint positioning and resource allocation for multi-UAV OFDMA
communication network.

For the proposed algorithm, the calculation of blocked regions in line 1
entails a computational complexity of $\mathcal{O}\left(KQ\right)$. The
computational complexity of lines 5-9 is dominated by solving the 3-D
positioning sub-problem (24) and resource allocation sub-problem (38), whose
complexities are $\mathcal{O}\left((3M)^{3.5}\right)$ and
$\mathcal{O}\left((KMN+MN)^{3.5}\right)$, respectively, by using the interior
point method [34]. The complexity of updating penalty multipliers in line 11
is $\mathcal{O}\left(KMN\right)$. Therefore, the worst-case computational
complexity of Algorithm 1 is
$\mathcal{O}\left(KQ+L_{\mathrm{out}}L_{\mathrm{in}}((3M)^{3.5}+(KMN+MN)^{3.5})\right)$,
where $L_{\mathrm{out}}$ and $L_{\mathrm{in}}$ denote the numbers of the
outer-loop and inner-loop iterations, respectively.

The inner-loop solves the penalized problem (14) by using the BSCA technique.
Due to the use of backtracking line search, we have

$\displaystyle
Z(\mathbf{X}^{{l+1}},\mathbf{P}^{{l}},\mathbf{C}^{{l}})-Z(\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}})\geq\tau\gamma_{1}^{l}(\tilde{\mathbf{X}}^{l\star}-\mathbf{X}^{{l}})^{\mathrm{T}}\nabla_{\mathbf{X}}Z\left(\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)\geq
0,$ (43)

and

$\displaystyle
Z\left(\mathbf{X}^{{l+1}},\mathbf{P}^{{l+1}},\mathbf{C}^{{l+1}}\right)-Z\left(\mathbf{X}^{{l+1}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)\geq\tau\gamma_{2}^{l}\begin{bmatrix}(\tilde{\mathbf{P}}^{l\star}-\mathbf{P}^{{l}})^{\mathrm{T}}\\\
(\tilde{\mathbf{C}}^{l\star}-\mathbf{C}^{{l}})^{\mathrm{T}}\end{bmatrix}^{\mathrm{T}}\begin{bmatrix}\nabla_{\mathbf{P}}Z\left(\mathbf{X}^{{l+1}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)\\\
\nabla_{\mathbf{C}}Z\left(\mathbf{X}^{{l+1}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)\end{bmatrix}\geq
0.$ (44)

Based on (43) and (44), we have

$\displaystyle
Z\left(\mathbf{X}^{{l+1}},\mathbf{P}^{{l+1}},\mathbf{C}^{{l+1}}\right)\geq
Z\left(\mathbf{X}^{{l+1}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right)\geq
Z\left(\mathbf{X}^{{l}},\mathbf{P}^{{l}},\mathbf{C}^{{l}}\right),$

which indicates that the objective value of problem (14) is non-decreasing
over the iteration. Since the objective value is upper-bounded, it always
converges to a finite value. For the outer-loop, the increasing update for
$\mathbf{\Lambda}$ leads to a sufficiently large multiplier which ensures that
the solution for the penalized problem (14) converges to a feasible solution
for the original problem in (11). Therefore, Algorithm 1 is guaranteed to
obtain a suboptimal solution for problem (11). The convergence of Algorithm 1
will be further evaluated by simulations in Section IV-B.

##  IV Performance Evaluation

In this section, we provide simulation results to evaluate the performance of
the proposed joint 3-D positioning and resource allocation scheme for multi-
UAV communication networks aided by geographic information.

###  IV-A Simulation Setup and Benchmark Schemes

As shown in Fig. 3, we consider a dense urban area of campus at Beihang
University with size $1,500\times 1,500\text{ m}^{2}$, i.e.,
$x_{\mathrm{D}}=1,500$, $y_{\mathrm{D}}=1,500$. The geometries of buildings
are obtained from the _OpenStreetMap_
database333https://www.openstreetmap.org/ and then processed into a number of
cubes. In this area, the maximum building height is $96$ m. Therefore, the
minimum flight altitude of UAVs, $h_{\mathrm{min}}$, is set to $100$ m such
that no collision will occur. The adopted simulation parameter settings are
listed in Table I [8, 36], unless specified otherwise. The users are randomly
generated on the ground, and each point in the simulation figures is the
average performance over 500 user location realizations.

Figure 3: (Left) An orthoimagery of the $1,500\times 1,500\text{ m}^{2}$ campus area at Beihang University, Beijing, China. (Right) The corresponding elevation map of buildings. TABLE I: Simulation Parameters Parameter | Description | Value
---|---|---
$P_{\mathrm{max}}$ | Maximum transmit power of each UAV | 30 dBm
$\sigma^{2}$ | Power of the noise at users | -107 dBm
$f_{c}$ | Carrier frequency | 5 GHz
$\alpha_{\mathrm{1}}$ | Channel gain exponent for LoS path | 2
$\alpha_{\mathrm{2}}$ | Channel gain exponent for NLoS path | 3.3
$\beta_{\mathrm{1}}$ |  | Channel gain at the reference distance of 1 m for LoS path
---
-46.43 dB
$\beta_{\mathrm{2}}$ |  | Channel gain at the reference distance of 1 m for NLoS path
---
-56.43 dB
$d_{\mathrm{min}}$ |  | Minimum inter-UAV distance
---
25 m
$\eta$ |  | Smooth parameter of the channel model
---
1000
$\zeta$ |  | Scaling factor in backtracking line search
---
0.9
$\tau$ |  | Coefficient in backtracking line search
---
0.01
$\epsilon_{\mathrm{l}}$ |  | Threshold for convergence of inner-loop
---
$10^{-3}$
$\epsilon_{\mathrm{L}}$ |  | Threshold for convergence of outer-loop
---
$10^{-4}$
$\lambda_{k,m,n}^{0}$ |  | Initial value of the multiplier
---
$0.2K/(MN)$

The initial state of the system, i.e.,
$\\{\bar{\mathbf{X}}^{0},\bar{\mathbf{P}}^{0},\bar{\mathbf{C}}^{0}\\}$, is
determined via the following steps. First, $M$ UAVs are deployed right above
$M$ out of $K$ users at an altitude of 500 m, where the users that are the
closest to the four endpoints of the area are sequentially selected. Then,
each user selects the one with the highest channel gain among the UAVs with
idle subcarriers for UAV association, and occupies the least used idle
subcarrier to minimize the interference. Finally, each UAV evenly allocates
the transmit power on its occupied subcarriers.

We label the proposed method as “Proposed”, and define three benchmark schemes
for performance comparison, namely “Fixed Association”, “K-means Position”,
and “No GoeInfo”, respectively, explained as follows.

  * •

_Fixed Association_ : In this scheme, the association variables $\mathbf{C}$
are fixed as initial value $\bar{\mathbf{C}}^{0}$, while the UAV positioning
and power allocation are jointly optimized following lines 5-9 in Algorithm 1.
Note that there is no need to perform any outer-loop iteration since
constraint (13) is always satisfied with $\mathbf{C}=\bar{\mathbf{C}}^{0}$.

  * •

_K-means Position_ : This scheme partitions the users into $M$ groups based on
their horizontal coordinates by employing the K-means algorithm, and places
the UAVs right above the cluster centers with an altitude of 500 m. Then, all
the steps in Algorithm 1 except line 6, are executed to optimize resource
allocation.

  * •

_No GeoInfo_ : This scheme assumes that geographic information is unavailable
for the system. Therefore, the joint positioning and resource allocation is
performed with the assumption of LoS A2G channels. The achievable rate is
calculated according to the actual LoS/NLoS channels at the obtained positions
of UAVs.

###  IV-B Simulation Results

First, we provide a demonstration of the proposed solution for the UAV
positioning and resource allocation in Fig. 4. The users are marked by
‘$\blacksquare$’, and the UAVs are marked by ‘$\bullet$’. Each UAV and its
associated users are marked with the same color, and are connected by a solid
(dashed) line for LoS (NLoS) channel condition. Links with different colors
represent that they occupy different subcarriers, and the linewidths are
proportional to the allocated transmit powers. As shown in Fig. 4 (a), 4 UAVs
are initially placed right above 4 out of 12 users at an altitude of 500 m.
Since there is one user who cannot establish LoS links with the UAVs, the
minimum achievable rate of the system is only 0.042 bits/sec/Hz. The final
state of the system by employing the proposed solution is shown in Fig. 4 (b),
where the optimization process of the UAV positioning is demonstrated by a
dash-dotted line. As can be observed, the UAVs tend to decrease their
altitudes and adjust their horizontal positions to decrease path losses as
well as get rid of building blockages. Besides, users farther away from their
associated UAVs are allocated with more power to compensate for higher path
loss. User association and subcarrier allocation have also adapted to
coordinate mutual interference along with UAV positioning and power
allocation. Compared to the initial state, the minimum achievable rate after
optimization increases from 0.042 bits/sec/Hz to 5.344 bits/sec/Hz, which
demonstrates the effectiveness of our proposed Algorithm 1.

(a) Initial state. (b) Final state.

Figure 4: Demonstration of the proposed solution for the UAV positioning and
resource allocation in Algorithm 1.

Fig. 5 presents the convergence of the proposed PDLIO algorithm for different
numbers of users with $M=4$ and $N=4$. In Fig. 5 (a), we evaluate the
convergence of the outer-loop iteration. As can be observed, the maximum
constraint violation
$\max\limits_{k,m,n}\bar{c}_{k,m,n}^{L}(1-\bar{c}_{k,m,n}^{L})$ decreases with
the iteration and falls below $10^{-4}$ within 30 iterations for all settings.
As the number of users increases, more iterations are needed to tackle more
blockage constraints and optimization variables. The convergence of the outer-
loop iterations indicates that appropriate values of the multipliers are
obtained to ensure that the obtained user-UAV-subcarrier association variables
in $\bar{\mathbf{C}}^{L}$ are binary. Moreover, the convergence of the inner-
loop iteration is shown in Fig. 5 (b). It can be observed that the inner-loop
converges within 20 iterations for all settings. The convergence of the outer-
loop iteration and inner-loop iteration guarantees the convergence of the
proposed algorithm, and ensures the feasibility of the obtained suboptimal
solution.

(a) Maximum constraint violation versus the number of outer-loop iterations.
(b) Objective value of problem (14) versus the number of inner-loop
iterations.

Figure 5: Convergence of Algorithm 1 for different numbers of users with 4
UAVs and 4 subcarriers.

Fig. 6 compares the minimum achievable rates for different schemes versus the
number of users $K$ with $M=4$ and $N=4$. It can be observed that the proposed
solution outperforms all other benchmark schemes. As the number of users
increases, the minimum achievable rate decreases. The reason is as follows.
First, as $K$ increases, the transmit power that can be potentially allocated
to each user is reduced, leading to lower power of the received signals.
Second, a limited number of subcarriers may lead to more severe mutual
interference with more users. Last, as $K$ increases, more blocked regions are
involved, and the UAVs tend to be deployed at higher altitudes to avoid signal
blockage, leading to higher path loss. Besides, the proposed solution has a
performance similar to “Fixed Association” and “No GeoInfo” for $K=M=N=4$,
because each user can be served by an independent UAV employing a unique
orthogonal subcarrier without mutual interference, and each UAV is deployed
right above its served user without any blockage. In addition, when the number
of UAVs and the available subcarrier resources are limited compared to the
number of users $K$, the proposed solution achieves a performance similar to
“Fixed Association”, since there is little freedom for user-UAV-subcarrier
association optimization. Finally, without geographic information, “No
GeoInfo” scheme gets the worst results for all settings. Since the blockage
effect is not properly considered during the UAV positioning and resource
allocation, this scheme cannot guarantee practical communication performance.

Figure 6: Minimum achievable rates for different schemes versus the number of
users with 4 UAVs and 4 subcarriers.

Fig. 8 compares the minimum achievable rates for different schemes versus the
number of UAVs with $K=8$ and $N=4$. As can be observed again, the proposed
scheme outperforms all other benchmark schemes. As $M$ increases, the minimum
achievable rate increases because there is more freedom for user association
and UAV positioning. Besides, compared to “Fixed Association” and “K-means
Position”, the growth rate of the proposed scheme is even greater, which
demonstrates the significance of user-UAV-subcarrier association optimization
and positioning optimization. Finally, although the performance of “No
GeoInfo” scheme increases with $M$, it shows the worst practical performance,
since the blockage effect is not properly addressed.

Figure 7: Minimum achievable rates for different schemes versus the number of
UAVs with $8$ users and $4$ subcarriers.

Figure 8: Minimum achievable rates for different schemes versus the number of
subcarriers with $8$ users and $4$ UAVs.

Finally, in Fig. 8 we compare the minimum achievable rates for different
schemes versus the number of subcarriers with $K=8$ and $M=4$. The proposed
method still outperforms all the other benchmark schemes. Besides, we can
observe that as $N$ increases, the performance gap between the proposed method
and “Fixed Association” scheme first increases and then decreases. The reason
is as follows. When $N$ is small, such as $N=K/M=2$, there is little freedom
for association optimization, thus leading to a small performance gap. For a
large $N$, such as $N=K=8$, each user can employ a unique orthogonal
subcarrier with little mutual interference. In such a case, the UAV-user-
subcarrier association optimization yields little improvement in rate
performance.

##  V Conclusion

In this paper, we proposed to use geographic information to characterize the
LoS blockage caused by buildings, for a multi-UAV OFDMA communication system.
Assisted by geographic information, a realistic channel model with blockage-
aware parameters was introduced. The UAV positioning and resource allocation
were then optimized to fully exploit the more favorable LoS channel conditions
between UAVs and users, such that the minimum achievable rate among all the
users is maximized. A penalty-based double-loop iterative algorithm was
proposed to solve the challenging optimization problem. The inner-loop is to
solve a penalized problem by optimizing UAV positioning sub-problem and
resource allocation sub-problem in an alternating way. The outer-loop is to
update multipliers to gradually decrease the violation of relaxed constraints
and finally obtain a feasible solution for the original problem. Simulation
results demonstrated that the proposed scheme can guarantee practical
communication performance compared to conventional LoS channel-based UAV
positioning design, and achieves higher a minimum achievable rate compared to
the stationary positioning and fixed association cases.

##  Appendix A Derivation of
$\nabla_{\mathbf{x}_{m}}g_{k}\left(\mathbf{x}_{m}\right)$,
$\nabla_{\mathbf{x}_{m}}\alpha_{k}\left(\mathbf{x}_{m}\right)$, and
$\nabla_{\mathbf{x}_{m}}\beta_{k}\left(\mathbf{x}_{m}\right)$

First,
$\nabla_{\mathbf{x}_{m}}\left(\frac{\min_{q\in\mathcal{Q}}\left\\{d_{k,q}(\mathbf{x}_{m})\right\\}}{\|\mathbf{x}_{m}-\mathbf{u}_{k}\|}\right)$
is derived as

$\displaystyle\nabla_{\mathbf{x}_{m}}\left(\frac{\min_{q\in\mathcal{Q}}\left\\{d_{k,q}(\mathbf{x}_{m})\right\\}}{\|\mathbf{x}_{m}-\mathbf{u}_{k}\|}\right)$
$\displaystyle=\frac{\nabla_{\mathbf{x}_{m}}\min\limits_{q\in\mathcal{Q}}\left\\{d_{k,q}(\mathbf{x}_{m})\right\\}}{\|\mathbf{x}_{m}-\mathbf{u}_{k}\|}-\frac{\min\limits_{q\in\mathcal{Q}}\left\\{d_{k,q}(\mathbf{x}_{m})\right\\}\nabla_{\mathbf{x}_{m}}\|\mathbf{x}_{m}-\mathbf{u}_{k}\|}{\|\mathbf{x}_{m}-\mathbf{u}_{k}\|^{2}}$
$\displaystyle=\frac{\mathbf{a}_{k}(\mathbf{x}_{m})}{\|\mathbf{x}_{m}-\mathbf{u}_{k}\|}-\frac{\min\limits_{q\in\mathcal{Q}}\left\\{d_{k,q}(\mathbf{x}_{m})\right\\}\nabla_{\mathbf{x}_{m}}\|\mathbf{x}_{m}-\mathbf{u}_{k}\|}{\|\mathbf{x}_{m}-\mathbf{u}_{k}\|^{2}}$
$\displaystyle=\frac{\mathbf{a}_{k}(\mathbf{x}_{m})}{\|\mathbf{x}_{m}-\mathbf{u}_{k}\|}-\frac{(\mathbf{a}_{k}^{\mathrm{T}}(\mathbf{x}_{m})\mathbf{x}_{m}-b_{k}(\mathbf{x}_{m}))(\mathbf{x}_{m}-\mathbf{u}_{k})}{\|\mathbf{x}_{m}-\mathbf{u}_{k}\|^{3}},$

with
$\left(\mathbf{a}_{k}(\mathbf{x}_{m}),b_{k}(\mathbf{x}_{m})\right)=\arg\min\limits_{\mathbf{a}_{k,q,i},b_{k,q,i}}\left\\{d_{k,q}(\mathbf{x}_{m})|q\in\mathcal{Q}\right\\}$.

Furthermore, $\nabla_{\mathbf{x}_{m}}s(\mathbf{x}_{m},\mathbf{u}_{k})$ is
derived as

$\displaystyle\nabla_{\mathbf{x}_{m}}s(\mathbf{x}_{m},\mathbf{u}_{k})$
$\displaystyle=\nabla_{\mathbf{x}_{m}}\left(\frac{1}{1+\exp\left(-\eta\frac{\min_{q\in\mathcal{Q}}\left\\{d_{k,q}(\mathbf{x}_{m})\right\\}}{\|\mathbf{x}_{m}-\mathbf{u}_{k}\|}\right)}\right)=-\frac{\nabla_{\mathbf{x}_{m}}\exp\left(-\eta\frac{\min_{q\in\mathcal{Q}}\left\\{d_{k,q}(\mathbf{x}_{m})\right\\}}{\|\mathbf{x}_{m}-\mathbf{u}_{k}\|}\right)}{\left(1+\exp\left(-\eta\frac{\min_{q\in\mathcal{Q}}\left\\{d_{k,q}(\mathbf{x}_{m})\right\\}}{\|\mathbf{x}_{m}-\mathbf{u}_{k}\|}\right)\right)^{2}}$
$\displaystyle=\eta
s(\mathbf{x}_{m},\mathbf{u}_{k})(1-s(\mathbf{x}_{m},\mathbf{u}_{k}))\nabla_{\mathbf{x}_{m}}\left(\frac{\min_{q\in\mathcal{Q}}\left\\{d_{k,q}(\mathbf{x}_{m})\right\\}}{\|\mathbf{x}_{m}-\mathbf{u}_{k}\|}\right).$

Therefore, $\nabla_{\mathbf{x}_{m}}\alpha_{k}\left(\mathbf{x}_{m}\right)$ and
$\nabla_{\mathbf{x}_{m}}\beta_{k}\left(\mathbf{x}_{m}\right)$ are given by

$\displaystyle\nabla_{\mathbf{x}_{m}}\alpha_{k}\left(\mathbf{x}_{m}\right)=(\alpha_{1}-\alpha_{2})\nabla_{\mathbf{x}_{m}}s(\mathbf{x}_{m}^{{l}},\mathbf{u}_{k})=\eta(\alpha_{1}-\alpha_{2})s(\mathbf{x}_{m},\mathbf{u}_{k})\left(1-s(\mathbf{x}_{m},\mathbf{u}_{k})\right)\times$
(45)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\frac{\left(\|\mathbf{x}_{m}-\mathbf{u}_{k}\|^{2}\mathbf{a}_{k}(\mathbf{x}_{m})-(\mathbf{a}_{k}^{\mathrm{T}}(\mathbf{x}_{m})\mathbf{x}_{m}-b_{k}(\mathbf{x}_{m}))(\mathbf{x}_{m}-\mathbf{u}_{k})\right)}{\|\mathbf{x}_{m}-\mathbf{u}_{k}\|^{3}},$
$\displaystyle\nabla_{\mathbf{x}_{m}}\beta_{k}\left(\mathbf{x}_{m}\right)=(\beta_{1}-\beta_{2})\nabla_{\mathbf{x}_{m}}s(\mathbf{x}_{m}^{{l}},\mathbf{u}_{k})=\eta(\beta_{1}-\beta_{2})s(\mathbf{x}_{m},\mathbf{u}_{k})\left(1-s(\mathbf{x}_{m},\mathbf{u}_{k})\right)\times$
(46)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\frac{\left(\|\mathbf{x}_{m}-\mathbf{u}_{k}\|^{2}\mathbf{a}_{k}(\mathbf{x}_{m})-(\mathbf{a}_{k}^{\mathrm{T}}(\mathbf{x}_{m})\mathbf{x}_{m}-b_{k}(\mathbf{x}_{m}))(\mathbf{x}_{m}-\mathbf{u}_{k})\right)}{\|\mathbf{x}_{m}-\mathbf{u}_{k}\|^{3}}.$

Finally, the derivation of
$\nabla_{\mathbf{x}_{m}}g_{k}\left(\mathbf{x}_{m}\right)$ is given as follows:

$\displaystyle\nabla_{\mathbf{x}_{m}}g_{k}\left(\mathbf{x}_{m}\right)=\nabla_{\mathbf{x}_{m}}\left(\frac{\beta_{k}(\mathbf{x}_{m})}{\|\mathbf{x}_{m}-\mathbf{u}_{k}\|^{\alpha_{k}(\mathbf{x}_{m})}}\right)=\frac{\nabla_{\mathbf{x}_{m}}\beta_{k}(\mathbf{x}_{m})}{\|\mathbf{x}_{m}-\mathbf{u}_{k}\|^{\alpha_{k}(\mathbf{x}_{m})}}-\frac{\beta_{k}(\mathbf{x}_{m})\nabla_{\mathbf{x}_{m}}\|\mathbf{x}_{m}-\mathbf{u}_{k}\|^{\alpha_{k}(\mathbf{x}_{m})}}{\|\mathbf{x}_{m}-\mathbf{u}_{k}\|^{2\alpha_{k}(\mathbf{x}_{m})}},$

in which

$\displaystyle\nabla_{\mathbf{x}_{m}}\|\mathbf{x}_{m}-\mathbf{u}_{k}\|^{\alpha_{k}(\mathbf{x}_{m})}=\nabla_{\mathbf{x}_{m}}\exp\left(\alpha_{k}(\mathbf{x}_{m})\log(\|\mathbf{x}_{m}-\mathbf{u}_{k}\|)\right)$
$\displaystyle=\|\mathbf{x}_{m}-\mathbf{u}_{k}\|^{\alpha_{k}(\mathbf{x}_{m})}\left(\nabla_{\mathbf{x}_{m}}\alpha_{k}(\mathbf{x}_{m})\cdot\log(\|\mathbf{x}_{m}-\mathbf{u}_{k}\|)+\frac{\alpha_{k}(\mathbf{x}_{m})\cdot\left(\mathbf{x}_{m}-\mathbf{u}_{k}\right)}{\|\mathbf{x}_{m}-\mathbf{u}_{k}\|^{2}}\right).$

Therefore, we have

$\displaystyle\nabla_{\mathbf{x}_{m}}g_{k}\left(\mathbf{x}_{m}\right)=$
$\displaystyle-
g_{k}\left(\mathbf{x}_{m}\right)\cdot\nabla_{\mathbf{x}_{m}}\alpha_{k}\left(\mathbf{x}_{m}\right)\cdot\log(\|\mathbf{x}_{m}-\mathbf{u}_{k}\|)$
(47) $\displaystyle-
g_{k}\left(\mathbf{x}_{m}\right)\cdot\frac{\alpha_{k}\left(\mathbf{x}_{m}\right)\cdot(\mathbf{x}_{m}-\mathbf{u}_{k})}{\|\mathbf{x}_{m}-\mathbf{u}_{k}\|^{2}}+\frac{\nabla_{\mathbf{x}_{m}}\beta_{k}\left(\mathbf{x}_{m}\right)}{\|\mathbf{x}_{m}-\mathbf{u}_{k}\|^{\alpha_{k}\left(\mathbf{x}_{m}\right)}}.$

## References

  * [1] 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, Third quarter 2019.
  * [2] 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.
  * [3] Z. Xiao, L. Zhu, Y. Liu, P. Yi, R. Zhang, X.-G. Xia, and R. Schober, “A survey on millimeter-wave beamforming enabled uav communications and networking,” _IEEE Commun. Surveys Tuts._ , vol. 24, no. 1, pp. 557–610, First quarter 2022.
  * [4] G. Geraci, A. Garcia-Rodriguez, M. M. Azari, A. Lozano, M. Mezzavilla, S. Chatzinotas, Y. Chen, S. Rangan, and M. D. 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, Third quarter 2022.
  * [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] L. Zhu, J. Zhang, Z. Xiao, X.-G. Xia, and R. Zhang, “Multi-UAV aided millimeter-wave networks: Positioning, clustering, and beamforming,” _IEEE Trans. Wireless Commun._ , vol. 21, no. 7, pp. 4637–4653, Jul. 2022\.
  * [7] A. Liu and V. K. N. Lau, “Optimization of multi-UAV-aided wireless networking over a ray-tracing channel model,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 9, pp. 4518–4530, Sep. 2019.
  * [8] P. Yi, L. Zhu, L. Zhu, Z. Xiao, Z. Han, and X.-G. Xia, “Joint 3-D positioning and power allocation for UAV relay aided by geographic information,” _IEEE Trans. Wireless Commun._ , vol. 21, no. 10, pp. 8148–8162, Oct. 2022.
  * [9] L. Zhu, J. Zhang, Z. Xiao, X. Cao, X.-G. Xia, and R. Schober, “Millimeter-wave full-duplex UAV relay: Joint positioning, beamforming, and power control,” _IEEE J. Select. Areas Commun._ , vol. 38, no. 9, pp. 2057–2073, Sep. 2020.
  * [10] Q. Hu, Y. Cai, A. Liu, G. Yu, and G. Y. Li, “Low-complexity joint resource allocation and trajectory design for UAV-aided relay networks with the segmented ray-tracing channel model,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 9, pp. 6179–6195, Sep. 2020.
  * [11] Y. Liu, K. Liu, J. Han, L. Zhu, Z. Xiao, and X.-G. Xia, “Resource allocation and 3-D placement for UAV-enabled energy-efficient IoT communications,” _IEEE Internet Thing J._ , vol. 8, no. 3, pp. 1322–1333, Feb. 2021.
  * [12] Y. Cai, W. Yuan, Z. Wei, C. Liu, S. Hu, and D. W. Kwan Ng, “Trajectory design and resource allocation for UAV-enabled data collection in wireless sensor networks with 3D blockages,” in _Proc. Int. Conf. 6G Networking (6GNet)_ , Paris, France, Jul. 2022.
  * [13] Z. Xiao, L. Zhu, and X.-G. Xia, “UAV communications with millimeter-wave beamforming: Potentials, scenarios, and challenges,” _China Communications_ , vol. 17, no. 9, pp. 147–166, Sep. 2020.
  * [14] 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, Mar. 2017.
  * [15] S. Yin, L. Li, and F. R. Yu, “Resource allocation and base station placement in downlink cellular networks assisted by multiple wireless powered UAVs,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 2, pp. 2171–2184, Feb. 2020\.
  * [16] C. Qiu, Z. Wei, X. Yuan, Z. Feng, and P. Zhang, “Multiple UAV-mounted base station placement and user association with joint fronthaul and backhaul optimization,” _IEEE Trans. Commun._ , vol. 68, no. 9, pp. 5864–5877, Sep. 2020.
  * [17] Y. Cai, F. Cui, Q. Shi, M. Zhao, and G. Y. Li, “Dual-UAV-enabled secure communications: Joint trajectory design and user scheduling,” _IEEE J. Select. Areas Commun._ , vol. 36, no. 9, pp. 1972–1985, Sep. 2018.
  * [18] C. Shen, T.-H. Chang, J. Gong, Y. Zeng, and R. Zhang, “Multi-UAV interference coordination via joint trajectory and power control,” _IEEE Trans. Signal Processing_ , vol. 68, pp. 843–858, Jan. 2020.
  * [19] C. Zhao, J. Liu, M. Sheng, W. Teng, Y. Zheng, and J. Li, “Multi-UAV trajectory planning for energy-efficient content coverage: A decentralized learning-based approach,” _IEEE J. Select. Areas Commun._ , vol. 39, no. 10, pp. 3193–3207, Oct. 2021.
  * [20] X. Liu, B. Lai, B. Lin, and V. C. M. Leung, “Joint communication and trajectory optimization for multi-UAV enabled mobile internet of vehicles,” _IEEE Trans. Intell. Transport. Syst._ , 2022 (Early Access).
  * [21] T. Bai, R. Vaze, and R. W. Heath, “Analysis of blockage effects on urban cellular networks,” _IEEE Trans. Wireless Commun._ , vol. 13, no. 9, pp. 5070–5083, Sep. 2014.
  * [22] J. Chen, U. Yatnalli, and D. Gesbert, “Learning radio maps for UAV-aided wireless networks: A segmented regression approach,” in _Proc. IEEE Int. Conf. Commun. (ICC)_ , Chengdu, China, Jul. 2017.
  * [23] S. Zhang and R. Zhang, “Radio map-based 3D path planning for cellular-connected UAV,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 3, pp. 1975–1989, Mar. 2021.
  * [24] Y. Dong, C. He, Z. Wang, and L. Zhang, “Radio map assisted path planning for UAV anti-jamming communications,” _IEEE Signal Processing Lett._ , vol. 29, pp. 607–611, Feb. 2022.
  * [25] J. Zhao, J. Liu, J. Jiang, and F. Gao, “Efficient deployment with geometric analysis for mmWave UAV communications,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 7, pp. 1115–1119, Jul. 2020.
  * [26] H. Kim and S. Han, “Interactive 3D building modeling method using panoramic image sequences and digital map,” _Multimedia Tools and Applications_ , vol. 77, no. 20, pp. 27 387–27 404, Oct. 2018.
  * [27] Q.-Y. Zhou and U. Neumann, “Complete residential urban area reconstruction from dense aerial LiDAR point clouds,” _Graphical Models_ , vol. 75, no. 3, pp. 118–125, May 2013.
  * [28] M. Razaviyayn, M. Hong, and Z.-Q. Luo, “A unified convergence analysis of block successive minimization methods for nonsmooth optimization,” _SIAM Journal on Optimization_ , vol. 23, no. 2, pp. 1126–1153, Jun. 2013\.
  * [29] Y. Yang, M. Pesavento, Z.-Q. Luo, and B. Ottersten, “Inexact block coordinate descent algorithms for nonsmooth nonconvex optimization,” _IEEE Trans. Signal Processing_ , vol. 68, pp. 947–961, Feb. 2020.
  * [30] L. Kang, Q. Wang, and H. W. Yan, “Building extraction based on Openstreetmap tags and very high spatial resolution image in urban area,” _International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences_ , vol. 42.3, pp. 715–718, Apr. 2018.
  * [31] J. Chen, U. Mitra, and D. Gesbert, “3D urban UAV relay placement: Linear complexity algorithm and analysis,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 8, pp. 5243–5257, Aug. 2021.
  * [32] D. P. Bertsekas, _Nonlinear Programming_ , 3rd ed. Belmont, MA: Athena Scientific, 1999.
  * [33] S. Diamond and S. Boyd, “CVXPY: A Python-embedded modeling language for convex optimization,” _Journal of Machine Learning Research_ , vol. 17, no. 83, pp. 1–5, Apr. 2016.
  * [34] S. Boyd and L. Vandenberghe, _Convex Optimization_. Cambridge, UK: Cambridge Univ. Press, Mar. 2004.
  * [35] M. L. Fisher, “The Lagrangian relaxation method for solving integer programming problems,” _Management Science_ , vol. 27, no. 1, pp. 1–18, 1981\.
  * [36] A. Al-Hourani, S. Kandeepan, and A. Jamalipour, “Modeling air-to-ground path loss for low altitude platforms in urban environments,” in _Proc. IEEE Glob. Commun. Conf._ , Austin, Tx, Dec. 2014, pp. 2898–2904.