# Coverage Analysis and Trajectory Optimization for Aerial Users with
Dedicated Cellular Infrastructure

Yujie Qin, Mustafa A. Kishk, Member, IEEE, and Mohamed-Slim Alouini, Fellow,
IEEE Yujie Qin and Mohamed-Slim Alouini are with Computer, Electrical and
Mathematical Sciences and Engineering (CEMSE) Division, King Abdullah
University of Science and Technology (KAUST), Thuwal, 23955-6900, Saudi
Arabia. Mustafa Kishk is with the Department of Electronic Engineering,
National University of Ireland, Maynooth, W23 F2H6, Ireland. (e-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS>slim.alouini@kaust.edu.sa).

###### Abstract

In this paper, we consider a novel cellular network for aerial users, which is
composed of dedicated base stations (BSs), whose antennas are directed towards
aerial users, and traditional terrestrial BSs (TBSs). Besides, the dedicated
BSs are deployed on roadside furniture, such as lampposts and traffic lights,
to achieve multiple features while occupying less space. Therefore, the
locations of dedicated BSs and TBSs are modeled by a Poisson-line-Cox-process
(PLCP) and Poisson point process (PPP), respectively. For the proposed
network, we first compute the aerial coverage probability and show that the
deployment of dedicated BSs improves the coverage probability in both high
dense areas and rural areas. We then consider a cellular-connected UAV that
has a flying mission and optimize its trajectory to maximize the minimal
achievable signal-to-interference-plus-noise ratio (SINR) (Max-Min SINR). To
obtain the Max-Min SINR and minimal time trajectory that satisfies the Max-Min
SINR, we proposed two algorithms that are practical in large-scale networks.
Finally, our results show that the optimal density of dedicated BSs which
maximizes Max-Min SINR decreases with the increase of the road densities.

###### Index Terms:

Stochastic geometry, aerial transportation, trajectory planning, Max-Min SINR,
PPP, PLCP

##  I Introduction

Featured by high flexibility and mobility, unmanned aerial vehicles (UAVs),
also known as drones, have drawn great attention in a broad range of civilian
applications, such as transportation, package delivery, communications, and
surveillance. For instance, UAVs are expected to help in delivery of goods and
even passengers [1], security monitoring [2], traffic control, and even
finishing multiple tasks simultaneously [3, 4]. Among all these appealing
applications, it is of vital importance to ensure a reliable connection
between UAVs and ground BSs, and only in this way, the safety of the public
and residents can be guaranteed. Observing that the currently existing
infrastructure is mainly designed for serving the ground users while the
aerial users experience a sidelobe gain, which may be difficult to meet the
goal of full coverage of future aerial networks. The idea of dedicated
infrastructure for aerial users is expected to be a potential solution to
improve the connectivity of aerial users [5].

However, building new infrastructure for aerial users is tricky and may cost a
lot of money and energy. Generally, BSs are on-grid and connected by physical
lines, to cooperate and share energy [6], and both the installation and
maintenance costs are high. Besides, deploying a high density of new BSs in
space-limited urban areas is inefficient. Hence, we are seeking an alternative
solution such as integrating the aerial dedicated antennas into existing
infrastructure, such as roadside furniture. Considering the lighting
infrastructure and traffic signs as a commonly used source, deploying the
dedicated antennas on such structures to form aerial dedicated BSs may be a
suitable way to go: instead of being connected by fiber, these dedicated BSs
can be connected by wireless backhaul links and form a dense mesh network [7,
8]. For instance, radio access nodes can be deployed on lampposts or rooftops
and wirelessly connect with each other by using mmWave. Instead of connecting
each node to the fiber point of presence (PoP), most of the nodes
(distribution nodes) connect to the Fiber PoP indirectly via LoS links with
other nodes. To ensure reliable connectivity and large capacity, high-dense
deployment over large geographic areas is required, e.g., if some nodes are
blocked due to some unpredictable reasons, data can be delivered through other
LoS paths. Such a hybrid system provides multiple features with less space and
without the need for new grid infrastructure, which saves a huge amount of
money. Besides, such a system is not something unrealistic and it actually has
been tested and implemented (for ground users only). For instance,
Facebook/Meta terragraph technology [9], deploys highly-dense transmitters and
receivers on the rooftops and light posts to achieve high-speed and low-cost
transmission. Therefore, in this work, we seek the feasibility of integrating
aerial dedicated BSs on roadside furniture for aerial users. However, due to
the wireless connection feature, this network has lower capacity, higher
latency, and higher cumulative traffic compared to the fiber backhaul.

Motivated by the concept of smart city, which aims to integrate information
and communication technology into existing systems and achieves efficient
networks, we study the performance of a novel cellular network which is
composed of traditional TBSs and dedicated BSs that are deployed on the
roadside furniture. To get system-level insights of the proposed network, we
analyze the system performance from the perspective of static UAVs and UAVs
with mobility, respectively. For static UAVs, we study the coverage
probabilities, and for a moving UAV, e.g., a UAV that has a mission such as
delivering a package, we optimize the trajectory and study the Max-Min SINR
during traveling. Specifically, we use tools from stochastic geometry to model
the locations of TBSs and aerial dedicated BSs and optimization tools to
design the UAVs’ trajectories given the random locations.

###  I-A Related Work

Literature related to this work can be categorized into: (i) incorporating
aerial users into existing networks, (ii) optimization of UAV communication
networks , and (iii) stochastic geometry-based analysis of UAV networks. A
brief discussion on related works is provided in the following lines.

Incorporating aerial users into existing networks. An overview of UAV
communication enabled by cellular networks was provided in [10] and some
discussions about network design and potential technologies were elaborated.
Authors in [11] discussed designing and implementing of future aerial networks
and drew a conclusion that the technologies which are suitable for aerial
networks are highly restricted by regulations and mechanical limitations.
Authors in [12] investigated the performance a typical rural aerial network
and showed that when the load is high, aerial users experience a poor
connection. An altitude-dependent path loss exponent and fading function was
analyzed in [13], in which the authors extended the models for G2G links to
A2G channels. Authors in [14] investigated the communication model for
terrestrial-aerial communication channels.

Optimization-based UAV communication networks. Authors in [15] proposed an
energy-efficient trajectory design for UAVs given the constraints such as
initial/final locations, velocities, and accelerations. In [16], authors
optimized the placement of UAVs to ensure coverage for a group of ground
users. A jointly optimization of UAV trajectory and transmit power was
investigated in [17] and results showed a lower outage probability. A high
mobility UAV communication relay was analyzed in [18]. The authors studied the
throughput maximization problems by designing UAV trajectory and power
allocations. Authors in [19, 20] investigated a Max-Min SNR based UAV
trajectory optimization problem. Authors in [21] jointly designed multiple
UAVs’ trajectories and transmit powers to maximize the throughput of a group
of ground users. In [22], the authors considered a UAV that collects data from
multiple wireless sensor networks and jointly optimized UAV’s trajectory and
communication scheduling to maximize the minimum average data collection rate.
In [23], authors studied the reliability of command and control channel
between UAVs and a traditional cellular network with massive MIMO. Besides, a
deep reinforcement learning-based approach of optimizing multiple cellular-
connected UAVs to minimize their interference with ground users was analyzed
in [24]. Authors in [25] studied the uplink inter-cell interference
coordination design for a cellular network simultaneously serving UAVs and
ground users. Typically, the maximum weighted sum-rate of UAVs and users is
jointly optimized based on power allocations and uplink cell associations.
Authors in [26] maximized the sum rate of UAV uplink transmission by
optimizing the precoding vectors at the multi-antenna UAVs. A deep
reinforcement learning-based intelligent navigation task of a cellular-
connected UAV network was investigated in [27] which aims to minimize the
weighted sum of time cost and expected outage duration alongside UAV
trajectory.

Stochastic geometry-based analysis of UAV networks. Stochastic geometry is
widely used to model large-scale communication networks [28, 29] and modeling
the locations of UAVs as a Poisson point process (PPP) is commonly used in the
literature [30, 31, 32]. Authors in [33, 34, 35] studied the downlink coverage
probability and averaged data rate by modeling the locations of TBSs and UAVs
by two independent PPPs. Beside PPP, Matern cluster process (MCP) is another
commonly used point process. Authors in [36, 37, 38] use this model to jointly
model the locations of UAVs and users, where the UAVs are deployed above
cluster centers and users are uniformly distributed in the cluster. Authors in
[39, 40] studied the laser powered UAVs and studied the density of the laser
beam directors to ensure energy coverage probability. Authors in [41] studied
a 3D two-hop cellular network with coexisting TBSs and UAVs to serve users. To
guarantee line-of-sight (LoS) coverage to all users, authors in [42]
investigated a determined 3D placement and orientations problem of UAVs.

Unlike the existing literature which improves UAV networks’ performance by
optimizing the locations and trajectories, the proposed system improves the
connectivity by deploying dedicated aerial BSs; and unlike the existing work
which converts the current BSs into new BSs or simply deploys new BSs, the
proposed system integrates the new BSs to existing infrastructure to save
space and be energy efficient. Besides, we study the performance of moving
UAVs based on stochastic geometry models, which is more practical due to the
large-scale feature. Besides, we would like to clarify that the coverage
probability of ground users decreases with the increasing density of aerial
dedicated users due to the increasing interference, as discussed in our
previous work [5]. The differences between this work and our previous work are
(i) in this work we consider a more complex system model in which dedicated
BSs are modeled by a PLCP while in previous work the dedicated BSs are
deployed by converting existing TBSs and (ii) since we optimize the UAV
trajectory in this work, we simplify the antenna radiation pattern model,
which can be considered as a special case in [5], to be feasible for
optimization.

###  I-B Contribution

This paper systematically investigates the feasibility and performance of
integrating dedicated BSs into existing networks. Our main contributions of
this work can be summarized below.

Novel system model. To provide reliable connectivity for aerial users, we
propose a novel wireless cellular network that is composed of TBSs and aerial
dedicated BSs. To integrate this new infrastructure into existing networks, we
consider these dedicated BSs to be deployed on the roadside furniture (install
on the top and share the energy infrastructure, grid, battery, or solar
panels), such as lampposts and traffic lights, to save space and be energy
efficient. By doing so, aerial users, like UAVs, are able to experience
mainlobe signals, hence, have better communication channels.

Stochastic geometry-based UAV trajectory optimization. Besides static UAVs, we
also analyze the improvement of implementing the proposed network on the
performance of a flying UAV. We consider a connectivity-constrained UAV which
flies from the initial location to the destination and optimizes its
trajectory to achieve the Max-Min SINR during traveling. Besides, we propose
an algorithm to obtain the minimal time trajectory that satisfies the Max-Min
SINR constraint. Given the nature of stochastic geometry, which aims to
analyze the performance of large-scale networks with a large number of BSs,
this minimal time trajectory is a sub-optimal solution to the real case (the
optimal solution is tricky to obtain due to the large number of BSs which
results into infinite possible routes), however, it is more practical since
the proposed system is more generic and focuses on the spatial statistics of
the system rather than a particular system setup.

System-level insights. We use tools from stochastic geometry to model the
locations of TBSs and dedicated BSs. Specifically, we model the locations of
dedicated BSs by a Poisson-line-Cox-process (PLCP). Therefore, we are able to
obtain the averaged system performance (averaged over the locations) and the
proposed algorithms fit different scenarios. Our results show that the
deployment of dedicated BSs improves the aerial coverage probability in both
urban areas and rural areas, and the optimal value of point density varies
with road densities. For moving UAVs, we show that the dedicated BSs improve
the Max-Min SINR for different road densities (both rural and urban areas). We
also observe an interesting trend in averaged minimal time of the trajectory:
the averaged minimal time will increase first and then decrease with the
increasing of the point densities.

##  II System Model

We consider a cellular-enabled UAV that can be connected to two types of
ground BSs, traditional TBSs which mainly provide service to ground users, and
dedicated BSs which are dedicated to aerial users [5], such as UAVs. Since the
traditional TBSs have existed for decades while the new infrastructure,
dedicated BSs, can be added to the existing street furniture [43], the
locations of the traditional TBSs are modeled by a PPP, $\Phi_{tb}$ with
density $\lambda_{tb}$, and the locations of the dedicated BSs are modeled by
a PLCP $\Phi_{db}$ with road density $\lambda_{l}$ and point density
$\lambda_{p}$, since PLCP is commonly used in modeling the random nodes
(dedicated BSs) on random lines (roads) [44], as shown in Fig. 1. Besides, to
characterize the SINR, we let $\Phi_{l}$ denote the set of lines (roads) with
density $\lambda_{l}$. Firstly, for static UAVs hovering at fixed locations to
provide a particular service, we are interested in the aerial coverage
probability of this novel setup of PPP-modeled TBSs and PLCP-modeled dedicated
aerial BSs. Secondly, for traveling UAVs, we are interested in the maximum
achievable minimum SINR (Max-Min SINR, $\gamma^{*}$) of the UAV trajectory and
the minimum time of finishing the traveling mission. We assume that the
reference UAV has a mission of flying from its initial location ${\rm S}$ to
the final destination ${\rm D}$ which is $L$ km away and in a random
direction. During the UAV’s route, it flies at a fixed altitude $h_{u}$ and we
want it to achieve a reliable connection with ground BSs, either TBSs or
dedicated BSs. To ensure that, the minimum SINR of the whole trajectory should
be above a predefined threshold: at any time of the UAV trajectory, the SINR
should be greater than $\gamma^{*}$.

Figure 1: Illustration of the system model.

As shown in Fig. 1, to ensure a reliable connectivity with aerial users,
dedicated antennas are deployed on roadside furniture, such as lampposts.
Instead of flying in a straight line between the initial point to the
destination, the reference UAV adjusts its trajectory to achieve a target
SINR. In the following part, we illustrate the proposed system from the
perspective of aerial coverage probability and optimal trajectory,
respectively.

###  II-A Communication Channel Model

The detailed communication channel model of the dedicated BSs is defined in
[5]. To simplify the selection of notation and have a feasible zone in the
trajectory optimization part, the antenna gain when the aerial user associates
with a dedicated BS in this work is

$\displaystyle g_{a}(d)=\left\\{\begin{aligned} g_{m},&\quad d<z_{db}\\\
g_{s},&\quad d>z_{db},\\\ \end{aligned}\right.$ (1)

where $d$ is the horizontal distance between the UAV and the dedicated BS and
$z_{db}$ denotes the horizontal coverage range of the mainlobe of dedicated
BSs, and $g_{m}$ and $g_{s}$ denote the mainlobe and sidelobe gain,
respectively.

Let $\rho_{tb}$ and $\rho_{db}$ be the transmission power of TBSs and
dedicated BSs, respectively. Given the horizontal distance between the UAV and
the serving TBS/dedicated BS, $R_{tb}$ and $R_{db}$, respectively, the
received power of the UAV is

$\displaystyle p_{tb}(R_{tb})=\left\\{\begin{aligned}
p_{tb,l}(R_{tb,l})=\eta_{l}\rho_{tb}g_{s}G_{l}D^{-\alpha_{l}}_{tb,l},&\quad\text{in
the case of LoS},\\\
p_{tb,n}(R_{tb,n})=\eta_{n}\rho_{tb}g_{s}G_{n}D^{-\alpha_{n}}_{tb,n},&\quad\text{in
the case of NLoS},\\\ \end{aligned}\right.$ $\displaystyle
p_{db}(R_{db})=\left\\{\begin{aligned}
p_{db,l}(R_{db,l})=g_{a}(R_{db,l})\eta_{l}\rho_{db}G_{l}D^{-\alpha_{l}}_{db,l},&\quad\text{in
the case of LoS},\\\
p_{db,n}(R_{db,n})=g_{a}(R_{db,n})\eta_{n}\rho_{db}G_{n}D^{-\alpha_{n}}_{db,n},&\quad\text{in
the case of NLoS},\\\ \end{aligned}\right.$ (2)

where $D_{tb,\\{l,n\\}}$ and $D_{db,\\{l,n\\}}$ are the Euclidean distances
and obtained by $D_{tb,\\{l,n\\}}=\sqrt{\Delta
h_{2}^{2}+R_{tb,\\{l,n\\}}^{2}}$ and $D_{db,\\{l,n\\}}=\sqrt{\Delta
h_{1}^{2}+R_{db,\\{l,n\\}}^{2}}$, respectively, in which $R_{tb,l}$,
$R_{tb,n}$, $R_{db,l}$ and $R_{db,n}$ denote the distance between the UAV and
the serving LoS TBS, NLoS TBS, LoS dedicated BS and NLoS dedicated BS,
respectively. $\Delta h_{1}=h_{u}-h_{db}$ and $\Delta h_{2}=h_{u}-h_{tb}$,
where $h_{tb}$ and $h_{db}$ are the heights of the TBSs and dedicated BSs and
$h_{u}$ is the altitude of UAVs, $\alpha_{l}$, $\eta_{l}$ and $G_{l}$ denote
the path loss, additional loss and channel fading of LoS transmission and
$\alpha_{n}$, $\eta_{n}$ and $G_{n}$ denote the path loss, additional loss and
channel fading of NLoS transmission.

Besides, the probability of establishing a LoS/NLoS link between the UAV and
the serving BS is given in [45] as

$\displaystyle P_{l}(r)$
$\displaystyle=\frac{1}{1+a\exp(-b(\frac{180}{\pi}\arctan(\frac{\Delta
h}{r})-a))},$ $\displaystyle P_{n}(r)$ $\displaystyle=1-P_{l}(r),$ (3)

where $r\in\\{R_{tb},R_{db}\\}$, $\Delta h\in\\{\Delta
h_{1}=h_{u}-h_{tb},\Delta h_{2}=h_{u}-h_{db}\\}$, and $a$ and $b$ are two
environment variables.

Assuming that the reference UAV associates with the BS that provides the
strongest average received power. We define the association probability with
TBSs and dedicated BSs as follows.

######  Definition 1 (Association Probability).

Let $\mathcal{A}_{db,l}(d,\theta_{d})$, $\mathcal{A}_{db,n}(d,\theta_{d})$,
$\mathcal{A}_{tb,l}(d,\theta_{d})$ and $\mathcal{A}_{tb,n}(d,\theta_{d})$ be
association probabilities between the reference UAV and the nearest LoS
dedicated BS, NLoS dedicated BS, nearest LoS TBS, and NLoS TBS, respectively,

$\displaystyle\mathcal{A}_{db,l}(d,\theta_{d})$
$\displaystyle=\mathcal{A}_{db,l-db,n}(d,\theta_{d})\mathcal{A}_{db,l-tb,l}(d)\mathcal{A}_{db,l-tb,n}(d),$
$\displaystyle\mathcal{A}_{db,n}(d,\theta_{d})$
$\displaystyle=\mathcal{A}_{db,n-db,l}(d,\theta_{d})\mathcal{A}_{db,n-tb,l}(d)\mathcal{A}_{db,n-tb,n}(d),$
$\displaystyle\mathcal{A}_{tb,l}(d,\theta_{d})$
$\displaystyle=\mathcal{A}_{tb,l-db,n}(d,\theta_{d})\mathcal{A}_{tb,l-db,l}(d,\theta_{d})\mathcal{A}_{tb,l-tb,n}(d),$
$\displaystyle\mathcal{A}_{tb,n}(d,\theta_{d})$
$\displaystyle=\mathcal{A}_{tb,n-db,l}(d,\theta_{d})\mathcal{A}_{tb,n-tb,l}(d)\mathcal{A}_{tb,n-db,n}(d,\theta_{d}),$
(4)

where $d$ denotes the horizontal communication distance and
$\theta_{d}=\arcsin(\frac{y_{l}}{d})$, in which $y_{l}$ is the distance to the
line from the PLCP that contains the serving or interfering dedicated BS (more
details are provided in Lemma 1 and Lemma 2 and Fig. 2), and
$\mathcal{A}_{a-b}(d,\theta_{d})$ denotes the probability that the average
received power from BS $a$ is stronger than the average received power from BS
$b$.

As mentioned, we are interested in the aerial coverage probability of the
proposed network. The coverage probability is defined as the probability that
a typical user is successfully served by the BS, which refers to the
probability that its SINR is greater than some predefined threshold.

######  Definition 2 (Coverage Probability).

Let $\tau$ be the SINR threshold, the coverage probability is defined as

$\displaystyle P_{\rm cov}=P_{\rm cov,tb}+P_{\rm cov,db},$ (5)

where $P_{\rm cov,tb}$ and $P_{\rm cov,db}$ denote the coverage probability
when the reference UAV associates with TBS and dedicated BS, respectively,

$\displaystyle P_{\rm cov,tb}$
$\displaystyle=\mathbb{E}_{R_{tb,l},\theta_{d}}[\mathcal{A}_{tb,l}(R_{tb,l},\theta_{d})\mathbb{P}({\rm
SINR_{tb,l}}(R_{tb,l})>\tau)]+\mathbb{E}_{R_{tb,n},\theta_{d}}[\mathcal{A}_{tb,n}(R_{tb,n},\theta_{d})\mathbb{P}({\rm
SINR_{tb,n}}(R_{tb,n})>\tau)],$ $\displaystyle P_{\rm cov,db}$
$\displaystyle=\mathbb{E}_{R_{db,l},\theta_{d}}[\mathcal{A}_{db,l}(R_{db,l},\theta_{d})\mathbb{P}({\rm
SINR_{db,l}}(R_{db,l})>\tau)]$
$\displaystyle\quad+\mathbb{E}_{R_{db,n},\theta_{d}}[\mathcal{A}_{db,n}(R_{db,n},\theta_{d})\mathbb{P}({\rm
SINR_{db,n}}(R_{db,n})>\tau)],$ (6)

in which SINR and the aggregated interference are, respectively, given by

$\displaystyle{\rm SINR}_{\\{tb,db\\}}(\\{R_{tb},R_{db}\\})$
$\displaystyle=\frac{p_{\\{tb,db\\}}(R_{\\{tb,db\\}})}{\sigma^{2}+I_{\\{tb,db\\}}},$
(7) $\displaystyle I_{\\{tb,db\\}}$
$\displaystyle=\sum_{tb_{l},i\in\Phi_{tb,l}\setminus
t_{0}}\eta_{l}\rho_{tb}g_{s}G_{l,tb_{l},i}D_{tb_{l},i}^{-\alpha_{l}}+\sum_{tb_{n},j\in\Phi_{tb,n}\setminus
t_{0}}\eta_{n}\rho_{tb}g_{s}G_{n,tb_{n},k}D_{tb_{n},k}^{-\alpha_{n}}$
$\displaystyle+\sum_{db_{l},k\in\Phi_{db}\setminus
t_{0}}\eta_{l}g_{a}(\sqrt{D_{db_{l},k}^{2}-\Delta
h_{1}^{2}})\rho_{db}G_{l,db_{l},k}D_{db_{l},k}^{-\alpha_{l}}$
$\displaystyle+\sum_{db_{n},m\in\Phi_{db,n}\setminus
t_{0}}\eta_{n}g_{a}(\sqrt{D_{db_{n},m}^{2}-\Delta
h_{1}^{2}})\rho_{db}G_{n,db_{n},m}D_{db_{n},m}^{-\alpha_{n}},$ (8)

where $t_{0}$ denotes the location of the serving BS, $\Phi_{tb,l}$ and
$\Phi_{tb,n}$ are subsets of $\Phi_{tb}$ denote the locations of the LoS TBS
and NLoS TBSs, $\Phi_{db,l}$ and $\Phi_{db,n}$ are subsets of $\Phi_{db}$
denote the locations of the LoS and NLoS dedicated BSs, respectively, and the
subscript $\\{l,n\\}$ denotes the LoS or NLoS transmission, as defined in
(II-A).

###  II-B Trajectory Design

As mentioned, we are also interested in the Max-Min SINR $\gamma^{*}$ that the
reference UAV can achieve when it flies from a starting point to a destination
and the minimal time given the target $\gamma^{*}$. Let $\mathbf{u}(t)$ be the
trajectory of UAV with S and D denoting the locations of the initial location
and the destination of the UAV. Let $T$ be defined as the time of finishing
the traveling mission from ${\rm S}$ to ${\rm D}$ and $\gamma^{t}$ denotes the
mean SINR achieved by UAV at time $t$, in which $0<t\leq T$.

We first observe that the mean SINR $\gamma^{t}$ at time $t$ is a function of
the serving distance (distance to the serving BS). Conditioned on the serving
BS being a dedicated BS or TBS, the mean SINR given the horizontal
communication distance $d$ are, respectively, given by

$\displaystyle{\gamma}_{db}(d)$
$\displaystyle=\mathbb{E}_{\theta_{d}}\bigg{[}\int_{0}^{\infty}P_{l}(d)P_{s,db,l}(\tau,d,\theta_{d})+P_{n}(d)P_{s,db,n}(\tau,d,\theta_{d}){\rm
d}\tau\bigg{]},$ $\displaystyle{\gamma}_{tb}(d)$
$\displaystyle=\mathbb{E}_{\theta_{d}}\bigg{[}\int_{0}^{\infty}P_{l}(d)P_{s,tb,l}(\tau,d,\theta_{d})+P_{n}(d)P_{s,tb,n}(\tau,d,\theta_{d}){\rm
d}\tau\bigg{]},$ (9)

where $P_{s,\\{\cdot\\}}(\tau,d,\theta_{d})$ is defined as the conditional
success probability and obtained by

$\displaystyle P_{s,\\{\cdot\\}}(\tau,d,\theta_{d})=\mathbb{P}({\rm
SINR_{\\{\cdot\\}}}(d)>\tau),$ (10)

in which ${\rm SINR_{\\{\cdot\\}}}(d)$ is given in (7) and the subscript
$\\{l,n\\}$ in $P_{s,\\{\cdot\\}}(\tau,d)$ denotes the LoS or NLoS
transmission, as defined in (II-A).

######  Definition 3 (Max-Min SINR).

Let $\bar{\gamma}$ be the minimum SINR during a UAV trajectory and
$\gamma^{*}$ is the Max-Min SINR of the UAV trajectory given an arbitrary
realization of locations of TBSs and dedicated BSs,

$\displaystyle\mathcal{P}_{1}:\quad$
$\displaystyle\gamma^{*}=\max_{\mathbf{u}(t)}\bar{\gamma},$ $\displaystyle{\rm
s.t.}\quad$ $\displaystyle\gamma^{t}=\max(\gamma_{tb}(d),\gamma_{db}(d)),$
$\displaystyle\gamma^{t}\geq\bar{\gamma},\quad\forall\quad 0<t\leq T,$ (11)

where $\gamma^{t}$ denotes the mean SINR achieved by UAV at time $t$.

Let $r_{tb}({\gamma})$ and $r_{db}(\gamma)$ be the maximal communication
distances between the reference UAV and the serving TBS and dedicated BS that
achieve a target $\gamma$, respectively, where $\gamma$ denotes a value of
mean SINR. The $r_{tb}({\gamma})$ and $r_{db}(\gamma)$ can be obtained by
solving the inverse function of (II-B).

Consequently, the optimization problem of the minimal time trajectory given
$\gamma^{*}$ is formulated as

$\displaystyle\mathcal{P}_{2}:$ $\displaystyle\min_{\mathbf{u}(t)}T$
$\displaystyle{\rm s.t.}\quad\mathbf{u}(0)={\rm S},\quad\mathbf{u}(T)={\rm
D},$ (12) $\displaystyle\qquad\mathbf{R}_{tb}\leq
r_{tb}({\gamma}^{*}),\quad\mathbf{R}_{db}\leq r_{db}(\gamma^{*}),$ (13)

where $\mathbf{R}_{tb}$ and $\mathbf{R}_{db}$ denote the distances between the
reference UAV and the serving TBSs and dedicated BSs along its trajectory,
respectively. Besides, we would like to mention that the initial location and
the destination of UAVs are totally random, which means that the UAVs are not
necessarily to start exactly above a BS and ends exactly above a BS.

##  III Coverage Probability

In this part, we provide the analysis for the aerial coverage probability.
Before deriving the equations for the coverage probability, we first need to
compute the association probability and some distance distributions, such as
the PDFs of $R_{tb}$ and $R_{db}$.

###  III-A Distance Distribution

Recall that $R_{tb,l}$, $R_{tb,n}$, $R_{db,l}$ and $R_{db,n}$ are the
distances to the nearest LoS TBS, NLoS TBS, LoS dedicated BS and NLoS
dedicated BS, respectively. The probability density functions (PDFs) and
cumulative density functions (CDFs) of the above mentioned distances are
provided in the following lemma.

######  Lemma 1 (PDFs and CDFs of the $R_{tb,l}$, $R_{tb,n}$, $R_{db,l}$ and
$R_{db,n}$).

The PDFs of the $R_{tb,l}$, $R_{tb,n}$, $R_{db,l}$ and $R_{db,n}$ are,
respectively, given by

$\displaystyle f_{R_{db,\\{l,n\\}}}(r)=$ $\displaystyle\bigg{[}2\pi
P_{\\{l,n\\}}(r)\cos\theta_{d}$
$\displaystyle+2\pi\lambda_{l}\int_{0}^{r}2\lambda_{p}\frac{r}{\sqrt{r^{2}-\rho^{2}}}P_{\\{l,n\\}}(r)\exp\bigg{(}-2\lambda_{p}\int_{0}^{\sqrt{r^{2}-\rho^{2}}}P_{\\{l,n\\}}(\sqrt{\rho^{2}+z^{2}}){\rm
d}z\bigg{)}{\rm d}\rho\bigg{]}$
$\displaystyle\exp\bigg{(}-2\lambda_{p}\int_{0}^{r\cos(\theta_{d})}P_{\\{l,n\\}}(\sqrt{z^{2}+r^{2}\sin^{2}(\theta_{d})}){\rm
d}z$
$\displaystyle\qquad-2\pi\lambda_{l}\int_{0}^{r}1-\exp\bigg{(}-2\lambda_{p}\int_{0}^{\sqrt{r^{2}-\rho^{2}}}P_{\\{l,n\\}}(\sqrt{\rho^{2}+z^{2}}){\rm
d}z\bigg{)}{\rm d}\rho\bigg{)},$ $\displaystyle f_{R_{tb,\\{l,n\\}}}(r)=$
$\displaystyle
2\pi\lambda_{tb}P_{\\{l,n\\}}(r)\exp\bigg{(}-2\pi\lambda_{tb}\int_{0}^{r}zP_{\\{l,n\\}}(z){\rm
d}z\bigg{)},$ (14)

consequently, the corresponding CDFs are obtained by taking the integration
and given by

$\displaystyle F_{R_{db,\\{l,n\\}}}(r,\theta_{d})$
$\displaystyle=1-\exp\bigg{(}-2\lambda_{p}\int_{0}^{r\cos(\theta_{d})}P_{\\{l,n\\}}(\sqrt{z^{2}+r^{2}\sin^{2}(\theta_{d})}){\rm
d}z$
$\displaystyle\qquad-2\pi\lambda_{l}\int_{0}^{r}1-\exp\bigg{(}-2\lambda_{p}\int_{0}^{\sqrt{r^{2}-\rho^{2}}}P_{\\{l,n\\}}(\sqrt{\rho^{2}+z^{2}}){\rm
d}z\bigg{)}{\rm d}\rho\bigg{)},$ $\displaystyle F_{R_{tb,\\{l,n\\}}}(r)$
$\displaystyle=1-\exp\bigg{(}-2\pi\lambda_{tb}\int_{0}^{r}zP_{\\{l,n\\}}(z){\rm
d}z\bigg{)},$ (15)

in which $\theta_{d}=\arcsin(\frac{y_{l}}{r})$ as shown in Fig. 11 and we need
this variable to compute the lower bound of integration in the Laplace
transform of interference in (3).

Recall that we assume that the reference UAV associates with the BS that
provides the strongest average received power and
$\mathcal{A}_{db,l}(d,\theta_{d})$, $\mathcal{A}_{db,n}(d,\theta_{d})$,
$\mathcal{A}_{tb,l}(d,\theta_{d})$ and $\mathcal{A}_{tb,n}(d,\theta_{d})$ are
association probabilities that the reference UAV associates with the nearest
LoS dedicated BS, NLoS dedicated BS, nearest LoS TBS and NLoS TBS at distance
$d$ away, respectively.

######  Lemma 2 (Association Probability).

The association probabilities given the horizontal communication distance, are
given by

$\displaystyle\mathcal{A}_{db,l}(d,\theta_{d})$
$\displaystyle=\mathcal{A}_{db,l-db,n}(d,\theta_{d})\mathcal{A}_{db,l-tb,l}(d)\mathcal{A}_{db,l-tb,n}(d)$
$\displaystyle=\bar{F}_{R_{db,n}}(d_{db,db,ln}(d),\theta_{d})\bar{F}_{R_{tb,l}}(d_{db,tb,ll}(d))\bar{F}_{R_{tb,n}}(d_{db,tb,ln}(d)),$
$\displaystyle\mathcal{A}_{db,n}(d,\theta_{d})$
$\displaystyle=\mathcal{A}_{db,n-db,l}(d,\theta_{d})\mathcal{A}_{db,n-tb,l}(d)\mathcal{A}_{db,n-tb,n}(d)$
$\displaystyle=\bar{F}_{R_{db,l}}(d_{db,db,nl}(d),\theta_{d})\bar{F}_{R_{tb,l}}(d_{db,tb,nl}(d))\bar{F}_{R_{tb,n}}(d_{db,tb,nn}(d)),$
$\displaystyle\mathcal{A}_{tb,l}(d,\theta_{d})$
$\displaystyle=\mathcal{A}_{tb,l-db,n}(d)\mathcal{A}_{tb,l-db,l}(d,\theta_{d})\mathcal{A}_{tb,l-tb,n}(d)$
$\displaystyle=\bar{F}_{R_{db,n}}(d_{tb,db,ln}(d),\theta_{d})\bar{F}_{R_{db,l}}(d_{tb,db,ll}(d),\theta_{d})\bar{F}_{R_{tb,n}}(d_{tb,tb,ln}(d)),$
$\displaystyle\mathcal{A}_{tb,n}(d,\theta_{d})$
$\displaystyle=\mathcal{A}_{tb,n-db,l}(d,\theta_{d})\mathcal{A}_{tb,n-tb,l}(d)\mathcal{A}_{tb,n-db,n}(d,\theta_{d})$
$\displaystyle=\bar{F}_{R_{db,l}}(d_{tb,db,nl}(d),\theta_{d})\bar{F}_{R_{db,n}}(d_{tb,db,nn}(d),\theta_{d})\bar{F}_{R_{tb,l}}(d_{tb,tb,nl}(d)),$
(16)

where $\bar{F}_{\\{\cdot\\}}(\cdot)$ denotes the complementary cumulative
density function (CCDF) which equals to $1-F_{\\{\cdot\\}}(\cdot)$, where
$F_{\\{\cdot\\}}(\cdot)$ is the CDF defined in (1) in Lemma 1 and the
horizontal distances $d_{\\{\cdot\\}}$ are provided in Appendix -A.

###  III-B Laplace Transform

Now we derive the Laplace transform for the aggregate interference in Lemma 3
later in this subsection. Recall that we model the locations of dedicated BSs
as a PLCP, which results in a slight difference in deriving the Laplace
transform of the aggregate interference compared to the interference of a PPP.

When the reference UAV associates with a dedicated BS, a typical line from the
PLP which passes through the serving dedicated BS exists. While when the
reference UAV associates with a TBS, this typical line does not exist, as
shown in Fig. 2. Hence, we discuss these two scenarios separately.

Figure 2: Illustration of interfers when the reference UAV associates with (a)
a dedicated BS, (b) a TBS.

As shown in Fig. 2 (a), when the reference UAV associates with a dedicated BS,
the aggregate interference is composed of three sets, (i) the interference
$I_{db,1}$ from the typical line $L_{0}$, (ii) the interference $I_{db,2}$
from the rest of the lines $\Phi_{l}^{{}^{\prime}}$, in which
$\Phi_{l}^{{}^{\prime}}=\Phi_{l}\backslash L_{0}$ denotes the set of lines,
and (iii) the interference from TBSs. Therefore, the aggregate interference
$I_{tb}$ when associating with a dedicated BS is

$\displaystyle I_{db}$ $\displaystyle=I_{db,1}+I_{db,2}+I_{tb}$
$\displaystyle=I_{db,1,l}+I_{db,1,n}+I_{db,2,l}+I_{db,2,n}+I_{tb,l}+I_{tb,n}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}\sum_{db_{l},i\in
L_{0}\setminus
t_{0}}p_{db}g_{a}(Z_{db_{l},i})\eta_{l}G_{l,db_{l},i}(Z_{db_{l},i}^{2}+\Delta
h_{1}^{2})^{-\alpha_{l}/2}$ $\displaystyle+\sum_{db_{n},j\in L_{0}\setminus
t_{0}}p_{db}g_{a}(Z_{db_{n},j})\eta_{n}G_{n,db_{n},j}(Z_{db_{n},j}^{2}+\Delta
h_{1}^{2})^{-\alpha_{n}/2}$
$\displaystyle+\sum_{db_{l},i\in\Phi_{l}^{{}^{\prime}}}p_{db}g_{a}(Z_{db_{l},i})\eta_{l}G_{l,db_{l},i}(Z_{db_{l},i}^{2}+\Delta
h_{1}^{2})^{-\alpha_{l}/2}+\sum_{db_{n},j\in\Phi_{l}^{{}^{\prime}}}p_{db}g_{a}(Z_{db_{n},j})\eta_{n}G_{n,db_{n},j}(Z_{db_{n},j}^{2}+\Delta
h_{1}^{2})^{-\alpha_{n}/2}$
$\displaystyle+\sum_{tb_{l,i}\in\Phi_{tb,l}}p_{tb}\eta_{l}g_{s}G_{l,tb_{l},i}D_{tb_{l},i}^{-\alpha_{l}}+\sum_{tb_{n},j\in\Phi_{tb,n}}p_{tb}\eta_{n}g_{s}G_{n,tb_{n},j}D_{tb_{n},j}^{-\alpha_{n}},$
(17)

where step (a) follows from the fact that the line process can be represented
by $\Phi_{l}=\Phi_{l}^{{}^{\prime}}\cup L_{0}$, $Z_{db,\\{\cdot\\}}$ denotes
the horizontal distance between the reference UAV and the interfering
dedicated BSs, and $t_{0}$ denotes the location of the serving BS.

As shown in Fig. 2 (b), when the reference UAV associates with a TBS, the
aggregate interference is given by

$\displaystyle I_{tb}$ $\displaystyle=I_{db,2}+I_{tb}$
$\displaystyle=I_{db,2,l}+I_{db,2,n}+I_{tb,l}+I_{tb,n}$
$\displaystyle=\sum_{db_{l},i\in\Phi_{l}}p_{db}g_{a}(Z_{db_{l},i})\eta_{l}G_{l,db_{l},i}(Z_{db_{l},i}^{2}+\Delta
h_{1}^{2})^{-\alpha_{l}/2}+\sum_{db_{n},j\in\Phi_{l}}p_{db}g_{a}(Z_{db_{n},j})\eta_{n}G_{n,db_{n},j}(Z_{db_{n},j}^{2}+\Delta
h_{1}^{2})^{-\alpha_{n}/2}$
$\displaystyle+\sum_{tb_{l},i\in\Phi_{tb,l}\setminus
t_{0}}p_{tb}\eta_{l}g_{s}G_{l,tb_{l},i}D_{tb_{l},i}^{-\alpha_{l}}+\sum_{tb_{n},j\in\Phi_{tb,n}\setminus
t_{0}}p_{tb}\eta_{n}g_{s}G_{n,tb_{n},j}D_{tb_{n},j}^{-\alpha_{n}},$ (18)

where $t_{0}$ denotes the location of the serving BS.

######  Lemma 3 (Laplace Transform of The Interference).

The Laplace transform the the aggregate interference is given by

$\displaystyle\mathcal{L}_{I_{db}}(s,\theta_{d})$
$\displaystyle=\mathcal{L}_{\rm gen}(s)\mathcal{L}_{\rm typ}(s,\theta_{d}),$
$\displaystyle\mathcal{L}_{I_{tb}}(s)$ $\displaystyle=\mathcal{L}_{\rm
gen}(s),$ (19)

where $\mathcal{L}_{\rm gen}(s)$ and $\mathcal{L}_{\rm typ}(s,\theta_{d})$
denote the Laplace transform of the interference from general BSs, including
TBSs and dedicated BSs on general lines, and typical dedicated BSs (dedicated
BSs on $L_{0}$), respectively,

$\displaystyle\mathcal{L}_{\rm
gen}(s)=\exp\bigg{(}-2\pi\lambda_{tb}\sum_{c=\\{l,n\\}}\int_{v_{tb,c}(d)}^{\infty}(1-\kappa_{1}(c,s,z))zP_{c}(z){\rm
d}z\bigg{)}$
$\displaystyle\times\prod_{c=\\{l,n\\}}\exp\bigg{(}-2\pi\lambda_{l}\int_{0}^{v_{1,c}(d)}1-\exp\bigg{(}-2\lambda_{p}\int_{\sqrt{v_{1,c}^{2}(d)-\rho^{2}}}^{\infty}(1-\kappa_{2}(c,s,z))P_{c}(\sqrt{u^{2}+\rho^{2}}){\rm
d}u\bigg{)}{\rm d}\rho\bigg{)}$
$\displaystyle\times\prod_{c=\\{l,n\\}}\exp\bigg{(}-2\pi\lambda_{l}\int_{v_{1,c}(d)}^{\infty}1-\exp\bigg{(}-2\lambda_{p}\int_{0}^{\infty}(1-\kappa_{2}(c,s,z))P_{c}(\sqrt{u^{2}+\rho^{2}}){\rm
d}u\bigg{)}{\rm d}\rho\bigg{)},$ (20) $\displaystyle\mathcal{L}_{\rm
typ}(s,\theta_{d})=\exp\bigg{(}-2\lambda_{p}\sum_{c=\\{l,n\\}}\int_{v_{0,c}(d)}^{\infty}(1-\kappa_{3}(c,s,z))P_{c}(x){\rm
d}x\bigg{)},$ (21)

where

$\displaystyle\kappa_{1}(c,s,z)$
$\displaystyle=\bigg{(}\frac{m_{c}}{m_{c}+sp_{tb}g_{s}\eta_{c}(z^{2}+\Delta
h_{2}^{2})^{-\alpha_{c}/2}}\bigg{)}^{m_{c}},$ $\displaystyle\kappa_{2}(c,s,z)$
$\displaystyle=\bigg{(}\frac{m_{c}}{m_{c}+sp_{db}g_{a}(\sqrt{u^{2}+\rho^{2}})\eta_{c}(u^{2}+\rho^{2}+\Delta
h_{1}^{2})^{-\alpha_{c}/2}}\bigg{)}^{m_{c}}$ $\displaystyle\kappa_{3}(c,s,z)$
$\displaystyle=\bigg{[}\frac{m_{c}}{m_{c}+sp_{db}g_{a}(x)\eta_{c}(x^{2}+\Delta
h_{1}^{2})^{-\alpha_{c}/2}}\bigg{]}^{m_{c}},$ (22)

and,

$\displaystyle v_{0,l}(d)$ $\displaystyle=\left\\{\begin{aligned}
d\cos\theta_{d},&\quad t_{0}\in\Phi_{db,l},\\\
\sqrt{d_{db,db,nl}^{2}(d)-(d\sin\theta_{d})^{2}},&\quad
t_{0}\in\Phi_{db,n},\\\ \end{aligned}\right.$ $\displaystyle v_{0,n}(d)$
$\displaystyle=\left\\{\begin{aligned}
\sqrt{(d_{db,db,ln}^{2}(d)-(d\sin\theta_{d})^{2})^{+}},&\quad
t_{0}\in\Phi_{db,l},\\\ d\cos\theta_{d}&\quad t_{0}\in\Phi_{db,n},\\\
\end{aligned}\right.$ $\displaystyle v_{tb,\\{l,n\\}}(d)$
$\displaystyle=\left\\{\begin{aligned} d_{db,tb,l\\{l,n\\}}(d),&\quad
t_{0}\in\Phi_{db,l},\\\ d_{db,tb,n\\{l,n\\}}(d),&\quad t_{0}\in\Phi_{db,n},\\\
d_{tb,tb,l\\{l,n\\}}(d),&\quad t_{0}\in\Phi_{tb,l},\\\
d_{tb,tb,n\\{l,n\\}}(d),&\quad t_{0}\in\Phi_{tb,n},\\\
\end{aligned}\right.v_{1,\\{l,n\\}}(d)=\left\\{\begin{aligned}
d_{db,db,l\\{l,n\\}}(d),&\quad t_{0}\in\Phi_{db,l},\\\
d_{db,db,n\\{l,n\\}}(d),&\quad t_{0}\in\Phi_{db,n},\\\
d_{tb,db,l\\{l,n\\}}(d),&\quad t_{0}\in\Phi_{tb,l},\\\
d_{tb,db,n\\{l,n\\}}(d),&\quad t_{0}\in\Phi_{tb,n},\\\ \end{aligned}\right.$
(23)

where the distances $d_{\\{\cdot\\}}$ are defined in Lemma 2.

###### Proof:

See Appendix -B. ∎

Now we are able to obtain the coverage probability of this network. Observe
from (2) that coverage probability can be represented using the success
probability,

$\displaystyle P_{\rm cov,db}$
$\displaystyle=\mathbb{E}_{R_{db,l}}[\mathcal{A}_{db,l}(R_{db,l},\theta_{d})\mathbb{P}({\rm
SINR_{db,l}}(R_{db,l})>\tau)]+\mathbb{E}_{R_{db,n}}[\mathcal{A}_{db,n}(R_{db,n},\theta_{d})\mathbb{P}({\rm
SINR_{db,n}}(R_{db,n})>\tau)],$
$\displaystyle=\mathbb{E}_{R_{db,l}}[\mathcal{A}_{db,l}(R_{db,l},\theta_{d})P_{s,db_{l}}(\tau,d,\theta_{d})]+\mathbb{E}_{R_{db,n}}[\mathcal{A}_{db,n}(R_{db,n},\theta_{d})P_{s,db_{n}}(\tau,d,\theta_{d})],$
(24) $\displaystyle P_{\rm cov,tb}$
$\displaystyle=\mathbb{E}_{R_{tb,l}}[\mathcal{A}_{tb,l}(R_{tb,l},\theta_{d})\mathbb{P}({\rm
SINR_{tb,l}}(R_{tb,l})>\tau)]+\mathbb{E}_{R_{tb,n}}[\mathcal{A}_{tb,n}(R_{tb,n})\mathbb{P}({\rm
SINR_{tb,n}}(R_{tb,n})>\tau)]$
$\displaystyle=\mathbb{E}_{R_{tb,l}}[\mathcal{A}_{tb,l}(R_{tb,l},\theta_{d})P_{s,t_{l}}(\tau,d)]+\mathbb{E}_{R_{tb,n}}[\mathcal{A}_{tb,n}(R_{tb,n},\theta_{d})P_{s,t_{n}}(\tau,d)],$
(25)

and we need the success probability to compute the mean SINR, as defined in
(II-B). Hence, we first derive the success probability in the following lemma.

######  Lemma 4 (Conditional Success Probability).

The conditional success probabilities given the reference UAV associates with
a LoS dedicated BS, NLoS dedicated BS, LoS TBS and NLoS TBS, respectively, are
given by

$\displaystyle P_{s,db,{\\{l,n\\}}}(\tau,d,\theta_{d})$
$\displaystyle=\sum_{k=1}^{m_{\\{l,n\\}}}\binom{m_{\\{l,n\\}}}{k}(-1)^{k+1}\mathcal{L}_{I_{db,\\{l,n\\}}+\sigma^{2}}(\xi_{db,\\{l,n\\}}(\tau,d),\theta_{d}),$
$\displaystyle P_{s,tb,{\\{l,n\\}}}(\tau,d)$
$\displaystyle=\sum_{k=1}^{m_{\\{l,n\\}}}\binom{m_{\\{l,n\\}}}{k}(-1)^{k+1}\mathcal{L}_{I_{tb,\\{l,n\\}}+\sigma^{2}}(\xi_{tb,\\{l,n\\}}(\tau,d)),$
(26)

where

$\displaystyle\xi_{db,\\{l,n\\}}(\tau,d)$
$\displaystyle=k\beta_{2}(m_{\\{l,n\\}})m_{\\{l,n\\}}\tau(p_{db}\eta_{\\{l,n\\}}g_{a}(d))^{-1}(d^{2}+\Delta
h_{1}^{2})^{\frac{\alpha_{\\{l,n\\}}}{2}},$
$\displaystyle\xi_{tb,\\{l,n\\}}(\tau,d)$
$\displaystyle=k\beta_{2}(m_{\\{l,n\\}})m_{\\{l,n\\}}\tau(p_{tb}g_{s}\eta_{\\{l,n\\}})^{-1}(d^{2}+\Delta
h_{2}^{2})^{\frac{\alpha_{\\{l,n\\}}}{2}}.$ (27)

in which the subscript $l,n$ in the Laplace transform of the interference term
denotes the LoS or NLoS transmission link between the UAV and serving BS,
which influences the lower bound of the integration of the interference, as
defined in (3), and $\beta_{2}(m)=(m!)^{-1/m}$ when $m\geq 1$ and this term
appears due to the use of the upper bound approximation of Gamma distribution,
which can be found in [46, 33].

Now we obtain the final expression of the coverage probability of the proposed
network.

######  Theorem 1 (Coverage Probability).

The coverage probability is given by

$\displaystyle P_{\rm cov}$ $\displaystyle=P_{\rm cov,tb}+P_{\rm cov,db}$
$\displaystyle=P_{\rm cov,tb,l}+P_{\rm cov,tb,n}+P_{\rm cov,db,l}+P_{\rm
cov,db,n},$ (28)

where,

$\displaystyle P_{\rm cov,db,\\{l,n\\}}$
$\displaystyle=\mathbb{E}_{\theta_{d}}\bigg{[}\int_{0}^{\infty}\mathcal{A}_{db,\\{l,n\\}}(d,\theta_{d})P_{s,db,\\{l,n\\}}(\tau,d,\theta_{d})f_{R_{db,\\{l,n\\}}}(r,\theta_{d}){\rm
d}r\bigg{]},$ $\displaystyle P_{\rm cov,tb,\\{l,n\\}}$
$\displaystyle=\mathbb{E}_{\theta_{d}}\bigg{[}\int_{0}^{\infty}\mathcal{A}_{tb,\\{l,n\\}}(d,\theta_{d})P_{s,tb,\\{l,n\\}}(\tau,d)f_{R_{tb,\\{l,n\\}}}(r){\rm
d}r\bigg{]},$ (29)

in which $\theta_{d}$ is uniformly distributed in $(0,\frac{\pi}{2})$, the
association probability $\mathcal{A}_{\\{\cdot\\}}$ is defined in (2) in Lemma
2, success probability $P_{s,\\{\cdot\\}}$ is defined in (4) in Lemma 4, and
the distance distribution is defined in (1) in Lemma 1.

######  Remark 1.

Observing from (2) that UAVs are mainly associated with LoS TBSs and dedicated
BSs since association probabilities are functions of CCDFs with horizontal
distances $d_{\\{\cdot\\}}$. Taking
$\mathcal{A}_{db,l}(d,\theta_{d})=\bar{F}_{R_{db,n}}(d_{db,db,ln}(d),\theta_{d})$
$\bar{F}_{R_{tb,l}}(d_{db,tb,ll}(d))$ $\bar{F}_{R_{tb,n}}(d_{db,tb,ln}(d))$
for example, the UAV associates with a LoS dedicated BS at a distance $d$ away
and $d_{db,db,ln}(d)$ denotes the distance to the nearest NLoS dedicated BS
and $d_{db,tb,ln}(d)$ denotes the nearest NLoS TBS. Generally, $\eta_{n}$ is
much lower than $\eta_{l}$ and $\alpha_{l}$ is less than $\alpha_{n}$,
therefore, $d_{db,db,ln}(d)\approx d_{db,tb,ln}(d)\approx\Delta h_{\\{1,2\\}}$
which means that
$\bar{F}_{R_{db,n}}(d_{db,db,ln}(d),\theta_{d})\approx\bar{F}_{R_{tb,n}}(d_{db,tb,ln}(d),\theta_{d})\approx
1$. That is, UAVs have a very high probability of associating with LoS TBSs or
dedicated BSs. Therefore, coverage probability is dominated by $P_{\rm
cov,db,l}+P_{\rm cov,tb,l}$.

##  IV Trajectory Optimization

In this section, we optimize the trajectory of a flying UAV to maintain a
reliable connection with BSs. To do so, we first need to compute the mean SINR
given the transmission distance, which is given in the following lemma.

######  Lemma 5 (Mean SINR).

Given the transmission distance $d$, the averaged SINR conditioned on
association with a dedicated BS or a TBS are, respectively, given by

$\displaystyle{\gamma}_{db}(d)=$
$\displaystyle\mathbb{E}_{\theta_{d}}\bigg{[}\sum_{c=\\{l,n\\}}P_{c}(d)\int_{0}^{\infty}\sum_{k=1}^{m_{c}}\binom{m_{c}}{k}(-1)^{k+1}\mathcal{L}_{I_{db,c}+\sigma^{2}}(\xi_{db,c}(\tau,d),\theta_{d}){\rm
d}\tau\bigg{]},$ $\displaystyle{\gamma}_{tb}(d)=$
$\displaystyle\mathbb{E}_{\theta_{d}}\bigg{[}\sum_{c=\\{l,n\\}}P_{c}(d)\int_{0}^{\infty}\sum_{k=1}^{m_{c}}\binom{m_{c}}{k}(-1)^{k+1}\mathcal{L}_{I_{tb,c}+\sigma^{2}}(\xi_{tb,c}(\tau,d)){\rm
d}\tau\bigg{]}.$ (30)

######  Remark 2.

We capture the correlation in the trajectory by considering the time
correlation in the serving distance. For instance, we optimize the UAV
trajectory while the locations of BSs are fixed and the distances to the
serving BSs are correlated on time. For the interference correlation, authors
in [47] showed that the interference correlation coefficient decreases
dramatically with slight increases in distances between two locations.
Therefore, we investigate it in an average sense. For instance, we compute the
mean SINR, which is a function of distance to the serving BSs while the
interference, as well as channel fadings, are averaged. Then, the trajectory
optimization is based on the mean SINR.

Now that the mean SINR is represented by a function of transmission distance,
and we can start the trajectory optimization and to maximize $\bar{\gamma}$.
Without loss of generality, we assume that the reference UAV travels from
${\rm S}(-L/2,0)$ to the destination located at ${\rm D}(L/2,0)$. This is due
to the stationarity of PPP and PLCP [48, 44]. Given the initial location and
destination, we first obtain the ground BSs that UAV associates with at the
beginning and end of the trajectory, which are denoted by $\mathbf{g}_{s}$ and
$\mathbf{g}_{d}$, which are both 2 by 1 vectors. Note that $\mathbf{g}_{s}$
and $\mathbf{g}_{d}$ are not necessarily to be the nearest BS to the ${\rm S}$
or the destination, e.g. $\mathbf{g}_{s}$ may not be the nearest BS to the
${\rm S}$, but the BS that provides the highest mean SINR. Besides, the
maximal achievable mean SINR with $\mathbf{g}_{s}$ and $\mathbf{g}_{d}$ is an
upper bound of the Max-Min SINR of the trajectory given an arbitrary but fixed
realization, otherwise, the reference UAV cannot leave or reach the initial
location or the destination. Let $\gamma_{u}$ denotes this upper bound.

Let $\mathbf{g}=\\{\mathbf{g}_{1},\mathbf{g}_{2},\cdots,\mathbf{g}_{n}\\}$ be
a sequence of the locations of ground BSs, either dedicated BSs or TBSs, that
the reference UAV associates with along its trajectory, and
$\mathbf{g}_{1}=\mathbf{g}_{s}$ and $\mathbf{g}_{n}=\mathbf{g}_{d}$. Note that
$\mathbf{g}$ is a $2\times n$ matrix, where the first row denotes the $x$
coordinates and the second row denotes the $y$ coordinates. Let
$\mathbf{r}(\gamma)=\\{r_{1}(\gamma),\cdots,r_{n}(\gamma)\\}$ be the sequence
of the maximum transmission distance given a target mean SINR $\gamma$, in
which $r_{\\{\cdot\\}}(\gamma)\in\\{r_{tb}(\gamma),r_{db}(\gamma)\\}$. For a
fixed realization, the constraint (13) in $\mathcal{P}_{2}$ is satisfied if
and only if a ground BS sequence existed which satisfies the following
constraint,

$\displaystyle|\mathbf{g}_{i}-\mathbf{g}_{i-1}|\leq\mathbf{r}_{i}(\gamma)+\mathbf{r}_{i+1}(\gamma),\quad\forall
i\in(2,n),$ (31)

therefore, the constraints for the locations of the association ground BSs are
equivalent to,

$\displaystyle|\mathbf{g}_{1}|$ $\displaystyle\leq
r_{db}(\gamma)\quad\text{or}\quad r_{tb}(\gamma),$
$\displaystyle|\mathbf{g}_{i}-\mathbf{g}_{i-1}|$
$\displaystyle\leq\mathbf{r}_{i}(\gamma)+\mathbf{r}_{i+1}(\gamma),\quad\forall
i\in(2,n),$ $\displaystyle|\mathbf{g}_{n}-{\rm D}(L,0)|$ $\displaystyle\leq
r_{db}(\gamma)\quad\text{or}\quad r_{tb}(\gamma),$ (32)

where $|\cdot|$ denotes the horizontal distance.

We now provide the details of solving $\mathcal{P}_{1}$ and obtaining the Max-
Min SINR $\gamma^{*}$. In this step, we first construct an association BS
graph based on the distances and $r_{tb}(\gamma)$ and $r_{db}(\gamma)$. If
(31) is satisfied, UAV can travel from these two BSs while maintaining a
reliable connection. Hence, these two BSs are connected on the graph. After
computing all the horizontal distances between the BSs, we obtain the graph
showing the connection of BSs. The target SINR can be achieved if there exists
at least one association BS sequence. Then, the Max-Min SINR is obtained by
iterating over different values of $\gamma$.

Figure 3: Illustration of (a) covered area formed by ground BSs that achieves
$\gamma^{*}$, (b) BS connection graph constructed by (a) and $\gamma^{*}$, (c)
covered area formed by $\mathbf{g}^{*}$.

We plot an illustration of the system model in Fig. 3 in which $\gamma^{*}$ is
achieved. In Fig. 3(a), the shadowed area is formed by the ground BSs with
radii $r_{tb}(\gamma^{*})$ and $r_{db}(\gamma^{*})$. If the UAV flies within
the shadowed area, $\gamma^{t}\geq\gamma^{*}$ at any time $t$ can always be
achieved, e.g., the minimum $\gamma^{t}=\gamma^{*}$ achieved at the boundary.
Hence, if we can find a trajectory from ${\rm S}$ to ${\rm D}$ within the
covered area, the target SINR can be satisfied. In Fig. 3(b), BSs are
connected if (31) is satisfied, and as shown, some BS association sequences
can be found to connect ${\rm S}$ and ${\rm D}$. Among all the BS association
sequence, let $\mathbf{g}^{*}$ denote the BS sequence that has the shortest
distance, which is shown in Fig. 3 (c).

Due to the nature of large-scale deployment of BSs, which results in a huge
number of possible association BS sequences, we choose the the BS sequence
that has the shortest horizontal distance, which is obtained by Dijkstra
algorithm and denoted by $\mathbf{g}^{*}$, and optimize the UAV trajectory
based on $\mathbf{g}^{*}$. The algorithm for the first step is proposed in
Algorithm 1. The complexity of the proposed algorithm is
$\mathcal{O}(\bar{M}!)$, where $\bar{M}=(\lambda_{tb}+\lambda_{db})S$ denotes
the average number of BSs in the simulation area and $S$ is the simulation
area.

Input: $\Phi_{\\{tb,db\\}}$: Locations of the ground BSs
${\rm S},{\rm D}$: Locations of the inital point and destination

Output: $\gamma^{*}$: Max Min SINR, $\mathbf{g}$, $\mathbf{r}(\gamma^{*})$

Initialization: $s_{m}=\gamma_{m}$, $s_{u}=\gamma_{u}$, $\mathbf{g}$,
$\mathbf{r}$, $i=0$, $s^{0}=s_{m}$
Function  _MaxMinSINR(_$\Phi_{\\{tb,db\\}}$ ,${\rm S}$ ,${\rm D}$_)_:

repeat

Update $r_{1}=r_{tb}(s^{i})$, $r_{2}=r_{db}(s^{i})$

Compute the horizontal distance between any two ground BSs, if the distance is
shorter than the sum of their maximal transmission distances, either $r_{1}$
or $r_{2}$, save the value, otherwise, save as infinity

Construct the graph based on the obtained distance and find the shortest path
between $g_{s}$ and $g_{d}$ via Dijkstra algorithm

if  _If the distance of the shortest path is finite_ then

Update $s^{i+1}=\frac{s_{u}+s^{i}}{2}$, $s_{m}=s^{i}$, $\mathbf{g}$,
$\mathbf{r}(\gamma^{*})$

else

Update $s^{i+1}=\frac{s_{m}+s^{i}}{2}$, $s_{u}=s^{i}$

end if

Update $i=i+1$

until  _$|s^{i}-s^{i-1}| <\epsilon$_

Update $\gamma^{*}=s^{i}$

return $\gamma^{*}$, $\mathbf{g}^{*}$, $\mathbf{r}(s^{i})$

End Function

Algorithm 1 Algorithm for Max Min SINR and $\mathbf{g}^{*}$

In what follows, we optimize the UAV trajectory to minimize the traveling time
$T_{\rm min}$ based on the association BS sequence $\mathbf{g}^{*}$. We obtain
the optimal trajectory by exhaustive search which is similar to Method II
introduced in [20], but we modify this method by adding one more step to
achieve a better performance. However, this trajectory is still a sub-optimal
solution and slightly worse than the real shortest time path, denoted by
$\mathbf{u}^{*}(t)$, which is obtained by Method I in [20], since we only
consider a route within the covered area formed by the BS sequence
$\mathbf{g}^{*}$. The reason that the Method I in [20] is not practical in
this work is that we consider a large-scale network and a large number of BSs,
which enables us to analyze the interference but cause an infinite number of
possible BS association sequences. Therefore, obtaining all the possible BS
association sequences and all the possible paths that can achieve the Max-Min
SINR is tricky.

As mentioned in [20] Method II, an approximation solution of the UAV
trajectory is obtained by graph constructed by the intersections of the
circles formed by the ground BSs and their maximal transmission distances.
Observing that UAVs are not necessarily to travel to the intersections of each
BS one by one, it can skip some intersections. For instance, in Fig. 4 (b), if
the reference UAV travels from $\mathbf{a}_{i}$ to
$\mathbf{a}_{i+n}^{{}^{\prime\prime}}$, it can travel along a straight line
without passing though any intersections. It implies that if we plot a
connection graph, the points $\mathbf{a}_{i}$ and
$\mathbf{a}_{i+n}^{{}^{\prime\prime}}$ should be connected directly, and the
points $\mathbf{a}_{i}$ and $\mathbf{a}_{i+n}^{{}^{\prime}}$ should not be
connected since its segment is not totally covered. Therefore, instead of
boundary quantization, we use the following lemma to check the connection
between intersections, ${\rm S}$ and ${\rm D}$.

Figure 4: Illustration of Lemma 6. (a) Instead of traveling from
$\mathbf{m}_{i}(1)$ to $\mathbf{m}_{i+1}(2)$ to $\mathbf{m}_{i+2}(2)$, the UAV
can fly from $\mathbf{m}_{i}(1)$ to $\mathbf{m}_{i+2}(2)$ directly. (b)
General case.

######  Lemma 6 (Possible Traveling Path of UAVs).

Let $\mathbf{m}_{i}(1),\mathbf{m}_{i}(2)$ be the locations of two
intersections of $B(\mathbf{g}_{i},r_{i})$ and $B(\mathbf{g}_{i+1},r_{i+1})$,
and $\mathbf{m}_{i+1}(1),\mathbf{m}_{i+1}(2)$ be the locations of two
intersections of $B(\mathbf{g}_{i+1},r_{i+1})$ and
$B(\mathbf{g}_{i+2},r_{i+2})$, as shown in Fig. 4 (a). If the segment
$\overrightarrow{\mathbf{m}_{i}(1)\mathbf{m}_{i+1}(1)}$ is covered by the
$B(\mathbf{g}_{i},r_{i})$ and $B(\mathbf{g}_{i+1},r_{i+1})$, it should satisfy
the following constraints,

$\displaystyle k_{\mathbf{m}_{i}(1)\mathbf{m}_{i+1}(1)}\geq
k_{\mathbf{m}_{i}(1)\mathbf{m}_{i+2}(1)}\geq
k_{\mathbf{m}_{i}(1)\mathbf{m}_{i+1}(2)},$ (33)

where $k_{ab}$ denotes the slope of the segment $\overrightarrow{ab}$.
Consequently, for a general case, if a path
$\overrightarrow{\mathbf{a}_{i}\mathbf{a}_{i+n}}$ is covered by the circles
formed by the BS $\mathbf{g}_{i},\cdots,\mathbf{g}_{i+n}$, it should satisfy,

$\displaystyle k_{\mathbf{a}_{i}\mathbf{a}_{i+n}}$
$\displaystyle\leq\min(k_{\mathbf{a}_{i}\mathbf{m}_{i+n}(1)},k_{\mathbf{a}_{i}\mathbf{m}_{i+n-1}(1)}\cdots,k_{\mathbf{a}_{i}\mathbf{m}_{i}(1)}),$
$\displaystyle k_{\mathbf{a}_{i}\mathbf{a}_{i+n}}$
$\displaystyle\geq\max(k_{\mathbf{a}_{i}\mathbf{m}_{i+n}(2)},k_{\mathbf{a}_{i}\mathbf{m}_{i+n-1}(2)}\cdots,k_{\mathbf{a}_{i}\mathbf{m}_{i}(2)}).$
(34)

###### Proof:

The geometric illustration is shown in Fig. 4 (b):
$\overrightarrow{\mathbf{a}_{i}\mathbf{a}_{i+n}^{{}^{\prime\prime}}}$ is
covered by the circles since
$k_{\mathbf{a}_{i}\mathbf{a}_{i+n}^{{}^{\prime\prime}}}$ is in the range,
$\overrightarrow{\mathbf{a}_{i}\mathbf{a}_{i+n}^{{}^{\prime}}}$ satisfies the
constraints due to its slope is greater than
$\min(k_{\mathbf{a}_{i}\mathbf{m}_{i+n}(1)},k_{\mathbf{a}_{i}\mathbf{m}_{i+n-1}(1)}\cdots,k_{\mathbf{a}_{i}\mathbf{m}_{i}(1)})$.
∎

######  Remark 3.

This step is a simplification of the discrete boundary step in [20] method II.
If the segment of the initial location and the destination is covered, say
$\mathbf{a}_{i},\mathbf{a}_{i+n}^{{}^{\prime\prime}}$, then the minimal
distance trajectory is achieved by the line
$\overrightarrow{\mathbf{a}_{i}\mathbf{a}_{i+n}^{{}^{\prime\prime}}}$. In this
case, discretizing the boundary over-complicates the analysis since we do not
compute the handover locations. In the case that the initial location and the
destination are not covered, $\mathbf{a}_{i},\mathbf{a}_{i+n}^{{}^{\prime}}$,
the minimal distance trajectory is achieved at the boundary point, e.g.,
$|\overrightarrow{\mathbf{a}_{i}\mathbf{m}_{i+2}(1)}|+|\overrightarrow{\mathbf{m}_{i+2}(1)\mathbf{a}_{i+n}^{{}^{\prime}}}|$.

Now we are able to construct a graph which is composed of all the possible
paths that the reference UAV can travel within the area covered by the circles
formed by BS sequence $\mathbf{g}^{*}$. The algorithm for this step is
provided in Algorithm 2.

Input: $\mathbf{g}^{*}$: Sequence of the locations of the ground BSs
$\mathbf{r}(\gamma^{*})$: Sequence of the maximal transmission distance

Output: $T_{\rm min}$: Minimal traveling time, $\mathbf{u}$: UAV trajectory

Initialization: $\mathbf{u}$, $T_{\rm min}$
Function  _MinTime(_$\mathbf{g}^{*}$ , $\mathbf{r}({\rm\gamma^{*}})$, S, D_)_:

Compute the intersections given the BS sequence $\mathbf{g}^{*}$ and the radii
$\mathbf{r}(\gamma^{*})$

Use Lemma 6 to obtain all the segments between the intersections, ${\rm S}$
and ${\rm D}$, which are covered by the area formed by $\mathbf{g}^{*}$ with
radii $\mathbf{r}(\gamma^{*})$

Compute the shortest path $\mathbf{u}$ and the minimal time $T_{\rm min}$ via
Dijkstra algorithm return $\mathbf{u}$, $T_{\rm min}$

End Function

Algorithm 2 Algorithm for minimal time trajectory

However, we would like to point out that the $T_{\rm min}$ obtained by
Algorithm 2 is a sub-optimal solution even though we use exhaustive search
since we only consider the routes among the BS sequence $\mathbf{g}^{*}$, but
this do not influence the obtained Max-Min SINR since $\gamma^{*}$ is obtained
by considering all the BSs. For instance, in Fig. 5(a) we plot the
$\mathbf{u}(t)$ obtained by Algorithm 2. Since we only consider the covered
area formed by $\mathbf{g}^{*}$, which has a smaller covered area than
considering all the BSs, such as the yellow area in Fig. 5(b), the obtained
$\mathbf{u}(t)$ is slightly longer than the optimal solution, e.g., in Fig.
5(a) the reference UAV travels to the intersection while it travels in a
straight line in Fig. 5(b).

Figure 5: The proposed trajectory which is a sub-optimal solution to the
minimal time trajectory.

##  V Numerical Results

In this section, we validate our analytical results with simulations and
evaluate the impact of various system parameters on the network performance.
Unless stated otherwise, we use the simulation parameters as listed herein
Table I.

TABLE I: Table of Parameters

Parameter | Symbol | Simulation Value
---|---|---
Density of TBSs | $\lambda_{tb}$ | 1 km-2
Height of TBSs, dedicated BSs and altitude of UAVs |  $h_{tb}$, $h_{db}$, $h_{u}$ |  $30$, $10$, $100$ m
Transmission power of TBS and dedicated BSs |  $\rho_{tb}$, $\rho_{db}$ |  $1$, $1$ w
Mainlobe, sidelobe gain and coverage range |  $g_{m}$, $g_{s}$, $z_{db}$ |  $10$, $0$ dB, $534$ m
UAV traveling distance and velocity |  $L$, $v$ |  $5$ km, $18$ m/s
SINR threshold | $\tau$ |  $0$ dB
N/LoS environment variable | $a,b$ | 12, 0.11
Noise power | $\sigma^{2}$ |  $10^{-9}$ W
N/LoS and TBS path-loss exponent | $\alpha_{n},\alpha_{l},\alpha_{tb}$ | $4,2.1,4$
N/LoS fading gain | $m_{n},m_{l}$ | $1,3$
N/LoS additional loss | $\eta_{n},\eta_{l}$ |  $-20,0$ dB

Figure 6: Coverage probability in the case of (a) high dense urban area
($a=12,b=0.11$) and $\lambda_{tb}=1$ km-2, (b) rural area ($a=4.88,b=0.43$,
$\eta_{\\{l,n\\}}=0$ dB) and $\lambda_{tb}=0.1$ /km-2.

For the simulation of the considered system setup, we compute the coverage
probability and the Max-Min SINR separately and run a large number of
iterations, $10^{3}$ for Max-Min SINR and $10^{6}$ for coverage probability,
to ensure accuracy. For each iteration of computing coverage and Max-Min SINR,
we first generate a PPP realization for the locations of TBSs and a PLCP
realization for the locations of the dedicated BSs. For the coverage
probability simulation, assuming the typical UAV is located at the origin, we
compute the SINR at the origin. Besides, channel fading is also generated for
each iteration. For the simulation of the Max-Min SINR, we first compute the
SINR given the transmission distance in the case of different road and point
densities. Next, we use Algorithm 1 and Algorithm 2 to compute the
$\gamma^{*}$ and $T_{\rm min}$. Note that the randomness of BS locations is
considered since PPP and PLCP realizations vary with iterations. To avoid UAVs
traveling a very long distance, we define a simulation area $S=25\times 25$
km2 to restrict the traveling area of UAVs. For the road density, we choose
the value $\lambda_{l}=(2,4,8)/\pi$ km/km2 in which low value is used for
rural areas and high value is used for urban areas [49].

In Fig. 6, we plot the coverage probability of the considered network. In high
dense urban area, with the increase of the dedicated BS density, the coverage
probability increases first due to a better communication channel and higher
received signal power. However, with further increase of the dedicated BS
density, the coverage probability decreases, which is because of higher
interference. Our results show that optimal value of point density exists for
different values of road density. With the decrease of the road density, the
optimal value of point density increases and the optimal dedicated BS density
is about $\lambda_{l}\lambda_{p}\pi\approx 4$ km-2. In rural areas (in which
$\lambda_{tb}=0.1$ km-2), the network also benefits from deploying dedicated
BSs but the density of dedicated BSs is much lower than high dense urban area,
which is about $\lambda_{l}\lambda_{p}\pi\approx 0.75$ km-2. Optimal dedicated
BS density decreases with the increases of the probability of establishing LoS
links due to the interference. We would like to point out that even if the
dedicated BSs are modeled by a PPP, aerial users coverage performance would
still improve. However, as stated in the motivation, it is more practical to
model them using PLCP to accurately capture the scenario of deploying them on
roadside furniture, which we believe is a most more cost and space efficient
deployment strategy.

Figure 7: Optimal point density under different TBS density in the case of (a)
high dense urban area ($a=12,b=0.11$), (b) rural area ($a=4.88,b=0.43$,
$\eta_{\\{l,n\\}}=0$ dB).

To further study the influence of TBS density on the dedicated BS density, we
plot optimal $\lambda_{p}$ under different $\lambda_{tb}$ and road density in
Fig. 7, in which optimal $\lambda_{p}$ is obtained by maximizing the coverage
probability. We show that optimal $\lambda_{p}$ decreases with the increase of
$\lambda_{tb}$ and rural areas require less dedicated BSs compared Fig. 7 (a)
with Fig. 7 (b).

Figure 8: Illustration of minimal time UAV trajectory which maximizes the
minimal SINR when $\lambda_{l}=4/\pi$ km/km2 and $\lambda_{p}=0.4$ km-1.

Fig. 8 shows the final minimal time trajectory of the reference UAV when
$\lambda_{l}=4/\pi$ km/km2 and $\lambda_{p}=0.4$ /km and the red dash circle
denotes the maximal transmission distance of the BSs that achieves
$\gamma^{*}=17.8$ dB. Since we consider that the BS on the street furniture is
dedicated to aerial users, the maximal transmission distance of dedicated BSs
is much larger than traditional TBSs. Besides, this trajectory shows the use
of Lemma 6 to skip some intersections, hence, achieving a better performance
without exhaustive search among the whole covered area.

Figure 9: (a) Max-Min SINR under different line densities and dedicated BS
density (point density). (b) $T_{\rm min}$ of the UAV trajectory achieves
$\gamma^{*}$ and $\lambda_{l}$ is in km/km2.

Fig. 9 shows the results of UAV trajectory optimization, in which the analysis
results are obtained based on (5) and the two proposed algorithms. In Fig. 9
(a), we show the Max-Min SINR under different road densities. While the
deployment of dedicated BSs improves the Max-Min SINR compared with the TBS-
only network, we show the optimal point densities for different line densities
and have the same shift trend as Fig. 6 (a). For a given line density, with
the increase of the point densities, the distance between BSs decreases,
however, the mean SINR given a transmission distance also decreases due to
higher interference. Therefore, both sides in (31) decreases. At low point
density, the distances decrease dramatically, while at high point density the
decrease of $\gamma$ has a higher impact.

In Fig. 9 (b), we show the $T_{\rm min}$ under different road densities.
Interestingly, $T_{\rm min}$ increases first with the increase of point
density and then decreases with the further increase of $\lambda_{p}$. The
reasons are as follows. At low point densities (low BS density), ${\rm S}$ and
${\rm D}$ are connected by a BS sequence with a small number of elements. To
connect each element, the corresponding radii are large, hence, the UAV
trajectory is short: as shown in Fig. 10 (a), ${\rm S}$ and ${\rm D}$ are
connected by three BSs with large radii. With the increase of the point
density, $T_{\rm min}$ increases first due to some association BS sequence
existing but the UAV needs to travel farther, and then decreases owing to more
nearby BSs existing and connected, as shown in Fig. 10 (b) and (c).

Figure 10: UAV trajectories in three different point densities,
$\lambda_{p,1}<\lambda_{p,2}<\lambda_{p,3}$, and achieve three different Max-
Min SINR $\gamma_{1}^{*}<\gamma_{3}^{*}<\gamma_{2}^{*}$, while the road
density keeps the same.

##  VI Conclusion

In this paper, we consider a novel cellular network for aerial users, which is
composed of dedicated BSs and traditional TBS, where the dedicated BSs are
deployed on roadside furniture. For the proposed network, we first compute the
aerial coverage probability and show that the deployment of dedicated BSs
improves the coverage probability in high dense urban areas and rural areas.
We then consider a cellular-connected UAV that has a flying mission and
optimize its trajectory to maximize the minimum SINR. To obtain the Max-Min
SINR and minimal time trajectory, we proposed two algorithms, which are more
practical in large-scale networks. Finally, our results show that the optimal
densities of dedicated BSs that maximize Max-Min SINR are different for
different road densities.

###  -A Horizontal Distances in Lemma 2

The horizontal distances in Lemma 2 are given by:

$\displaystyle d_{db,tb,ll}(d)$
$\displaystyle=\sqrt{\bigg{(}g_{a}(d)^{-\frac{2}{\alpha_{l}}}(d^{2}+\Delta
h_{1}^{2})-\Delta h_{2}^{2}\bigg{)}^{+}},$ $\displaystyle d_{db,tb,ln}(d)$
$\displaystyle=\sqrt{\bigg{(}\bigg{(}\frac{\eta_{n}}{\eta_{l}g_{a}(d)}\bigg{)}^{\frac{2}{\alpha_{n}}}(d^{2}+\Delta
h_{1}^{2})^{\frac{\alpha_{l}}{\alpha_{n}}}-\Delta h_{2}^{2}\bigg{)}^{+}},$
$\displaystyle d_{db,db,ll}(d)$ $\displaystyle=d,$ $\displaystyle
d_{db,db,ln}(d)$ $\displaystyle=\left\\{\begin{aligned}
\sqrt{\bigg{(}\bigg{(}\frac{\eta_{n}}{\eta_{l}}\bigg{)}^{\frac{2}{\alpha_{n}}}(d^{2}+\Delta
h_{1}^{2})^{\frac{\alpha_{l}}{\alpha_{n}}}-\Delta
h_{1}^{2}\bigg{)}^{+}},&\quad d<z_{db},\\\
\max\bigg{[}\sqrt{\bigg{(}\bigg{(}\frac{\eta_{n}}{\eta_{l}}\bigg{)}^{\frac{2}{\alpha_{n}}}(d^{2}+\Delta
h_{1}^{2})^{\frac{\alpha_{l}}{\alpha_{n}}}-\Delta h_{1}^{2}\bigg{)}^{+}},&\\\
\min\bigg{[}z_{db},\sqrt{\bigg{(}\bigg{(}\frac{g_{a}(d)\eta_{n}}{\eta_{l}}\bigg{)}^{\frac{2}{\alpha_{n}}}(d^{2}+\Delta
h_{1}^{2})^{\frac{\alpha_{l}}{\alpha_{n}}}-\Delta
h_{1}^{2}\bigg{)}^{+}}\bigg{]}\bigg{]},&\quad d>z_{db},\\\
\end{aligned}\right.$ (35) $\displaystyle d_{db,tb,nl}(d)$
$\displaystyle=\sqrt{\bigg{(}\bigg{(}\frac{\eta_{l}}{\eta_{n}g_{a}(d)}\bigg{)}^{\frac{2}{\alpha_{l}}}(d^{2}+\Delta
h_{1}^{2})^{\frac{\alpha_{n}}{\alpha_{l}}}-\Delta h_{2}^{2}\bigg{)}^{+}},$
$\displaystyle d_{db,tb,nn}(d)$
$\displaystyle=\sqrt{\bigg{(}g_{a}(d)^{-\frac{2}{\alpha_{n}}}(d^{2}+\Delta
h_{1}^{2})-\Delta h_{2}^{2}\bigg{)}^{+}},$ $\displaystyle d_{db,db,nl}(d)$
$\displaystyle=\left\\{\begin{aligned}
\min\bigg{[}\sqrt{\bigg{(}\bigg{(}\frac{\eta_{l}}{\eta_{n}}\bigg{)}^{\frac{2}{\alpha_{l}}}(d^{2}+\Delta
h_{1}^{2})^{\frac{\alpha_{n}}{\alpha_{l}}}-\Delta h_{1}^{2}\bigg{)}^{+}},&\\\
\max\bigg{[}\sqrt{\bigg{(}\bigg{(}\frac{\eta_{l}}{g_{a}(d)\eta_{n}}\bigg{)}^{\frac{2}{\alpha_{l}}}(d^{2}+\Delta
h_{1}^{2})^{\frac{\alpha_{n}}{\alpha_{l}}}-\Delta
h_{1}^{2}\bigg{)}^{+}},z_{db}\bigg{]}\bigg{]},&\quad d<z_{db}\\\
\sqrt{\bigg{(}\bigg{(}\frac{\eta_{l}}{\eta_{n}}\bigg{)}^{\frac{2}{\alpha_{l}}}(d^{2}+\Delta
h_{1}^{2})^{\frac{\alpha_{n}}{\alpha_{l}}}-\Delta
h_{1}^{2}\bigg{)}^{+}},&\quad d>z_{db},\\\ \end{aligned}\right.$
$\displaystyle d_{db,db,nn}(d)$ $\displaystyle=d,$ (36) $\displaystyle
d_{tb,tb,ll}(d)$ $\displaystyle=d,$ $\displaystyle d_{tb,tb,ln}(d)$
$\displaystyle=\sqrt{\bigg{(}\bigg{(}\frac{\eta_{n}}{\eta_{l}}\bigg{)}^{\frac{2}{\alpha_{n}}}(d^{2}+\Delta
h_{2}^{2})^{\frac{\alpha_{l}}{\alpha_{n}}}-\Delta h_{2}^{2}\bigg{)}^{+}},$
$\displaystyle d_{tb,db,ll}(d)$ $\displaystyle=\left\\{\begin{aligned}
\min\bigg{[}\sqrt{\bigg{(}g_{a}(d)^{\frac{2}{\alpha_{l}}}(d^{2}+\Delta
h_{2}^{2})-\Delta h_{1}^{2}\bigg{)}^{+}},z_{db}\bigg{]},&\quad d<z_{db},\\\
d,&\quad d>z_{db},\\\ \end{aligned}\right.$ $\displaystyle d_{tb,db,ln}(d)$
$\displaystyle=\min\bigg{[}\sqrt{\bigg{(}\bigg{(}\frac{g_{a}(d)\eta_{n}}{\eta_{l}}\bigg{)}^{\frac{2}{\alpha_{n}}}(d^{2}+\Delta
h_{2}^{2})^{\frac{\alpha_{l}}{\alpha_{n}}}-\Delta h_{1}^{2}\bigg{)}^{+}},$
$\displaystyle\min\bigg{[}z_{db},\sqrt{\bigg{(}\bigg{(}\frac{\eta_{n}}{\eta_{l}}\bigg{)}^{\frac{2}{\alpha_{n}}}(d^{2}+\Delta
h_{2}^{2})^{\frac{\alpha_{l}}{\alpha_{n}}}-\Delta
h_{1}^{2}\bigg{)}^{+}}\bigg{]}\bigg{]},$ (37) $\displaystyle d_{tb,tb,nl}(d)$
$\displaystyle=\sqrt{\bigg{(}\bigg{(}\frac{\eta_{l}}{\eta_{n}}\bigg{)}^{\frac{2}{\alpha_{l}}}(d^{2}-\Delta
h_{2}^{2})^{\frac{\alpha_{n}}{\alpha_{l}}}-\Delta h_{1}^{2}\bigg{)}^{+}},$
$\displaystyle d_{tb,tb,nn}(d)$ $\displaystyle=d,$ $\displaystyle
d_{tb,db,nl}(d)$
$\displaystyle=\min\bigg{[}\sqrt{\bigg{(}\bigg{(}\frac{g_{a}(d)\eta_{l}}{\eta_{n}}\bigg{)}^{\frac{2}{\alpha_{l}}}(d^{2}+\Delta
h_{2}^{2})^{\frac{\alpha_{n}}{\alpha_{l}}}-\Delta h_{1}^{2}\bigg{)}^{+}},$
$\displaystyle\max\bigg{[}z_{db},\sqrt{\bigg{(}\bigg{(}\frac{\eta_{l}}{\eta_{n}}\bigg{)}^{\frac{2}{\alpha_{l}}}(d^{2}+\Delta
h_{2}^{2})^{\frac{\alpha_{n}}{\alpha_{l}}}-\Delta
h_{1}^{2}\bigg{)}^{+}}\bigg{]}\bigg{]},$ $\displaystyle d_{tb,db,nn}(d)$
$\displaystyle=\min\bigg{[}\sqrt{\bigg{(}g_{a}(d)^{\frac{2}{\alpha_{n}}}(d^{2}+\Delta
h_{2}^{2})-\Delta
h_{1}^{2}\bigg{)}^{+}},\max\bigg{[}z_{db},\sqrt{(d^{2}+\Delta
h_{2}^{2})-\Delta h_{1}^{2}}\bigg{]}\bigg{]},$ (38)

where $(a)^{+}$ denotes $\max(0,a)$.

###  -B Proof of Lemma 3

Figure 11: Illustration of the lower bound of integration in the case of
associating with a dedicated BS.

The Laplace transform of the interference when associated with a dedicated BS
is obtained by

$\displaystyle\mathcal{L}_{I_{db}}(s)=$
$\displaystyle\mathbb{E}_{I_{db}}[\exp(-s(I_{db}))]$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}$
$\displaystyle\mathbb{E}_{\Phi_{db}}\bigg{[}\exp\bigg{(}-s\bigg{(}\sum_{db_{l,i}\in\Phi_{l}^{{}^{\prime}}\cup
L_{0}}p_{db}g_{a}(Z_{l})\eta_{l}G_{l}(Z_{l,i}^{2}+\Delta
h_{1}^{2})^{-\alpha_{l}/2}$
$\displaystyle+\sum_{db_{n,j}\Phi_{l}^{{}^{\prime}}\cup
L_{0}}p_{db}g_{a}(Z_{n})\eta_{n}G_{n}(Z_{n,j}^{2}+\Delta
h_{1}^{2})^{-\alpha_{n}/2}\bigg{)}\bigg{)}$
$\displaystyle\times\mathbb{E}_{\Phi_{tb}}\bigg{[}\exp\bigg{(}-s\bigg{(}\sum_{t_{l,i}\in\Phi_{tb}}p_{tb}g_{s}\eta_{l}G_{l}D_{t_{l,i}}^{-\alpha_{l}}+\sum_{t_{n,j}\in\Phi_{tb}}p_{tb}g_{s}\eta_{n}G_{n}D_{t_{n,j}}^{-\alpha_{n}}\bigg{)}\bigg{)}\bigg{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\Phi_{db}}\bigg{[}\prod_{db_{l,i}\in\Phi_{l}^{{}^{\prime}}\cup
L_{0}}\exp\bigg{(}-sp_{db}g_{a}(Z_{l})\eta_{l}G_{l}(Z_{l}^{2}+\Delta
h_{1}^{2})^{-\alpha_{l}/2}\bigg{)}$
$\displaystyle\times\prod_{db_{n,j}\Phi_{l}^{{}^{\prime}}\cup
L_{0}}\exp\bigg{(}-sp_{db}g_{s}g_{a}(Z_{n})\eta_{n}G_{n}(Z_{n}^{2}+\Delta
h_{1}^{2})^{-\alpha_{n}/2}\bigg{)}\bigg{]}$
$\displaystyle\times\mathbb{E}_{\Phi_{tb}}\bigg{[}\prod_{t_{l,i}\in\Phi_{tb}}g_{s}\exp\bigg{(}-sp_{tb}g_{s}\eta_{l}G_{l}D_{t_{l}}^{-\alpha_{l}}\bigg{)}\prod_{t_{n,j}\in\Phi_{tb}}\exp\bigg{(}-sp_{tb}g_{s}\eta_{n}G_{n}D_{t_{n}}^{-\alpha_{n}}\bigg{)}\bigg{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\Phi_{db}}\bigg{[}\prod_{db_{l,i}\Phi_{l}^{{}^{\prime}}\cup
L_{0}}\bigg{(}\frac{m_{l}}{m_{l}+sp_{db}g_{a}(Z_{l})\eta_{l}(Z_{l}^{2}+\Delta
h_{1}^{2})^{-\alpha_{l}/2}}\bigg{)}^{m_{l}}$
$\displaystyle\times\prod_{db_{n,j}\Phi_{l}^{{}^{\prime}}\cup
L_{0}}\bigg{(}\frac{m_{n}}{m_{n}+sp_{db}g_{a}(Z_{n})\eta_{n}G_{n}(Z_{n}^{2}+\Delta
h_{1}^{2})^{-\alpha_{n}/2}}\bigg{)}^{m_{n}}\bigg{]}$
$\displaystyle\times\mathbb{E}_{\Phi_{tb}}\bigg{[}\prod_{t_{l,i}\in\Phi_{tb}}\bigg{(}\frac{m_{l}}{m_{l}+sp_{tb}g_{s}\eta_{l}D_{t_{l}}^{-\alpha_{l}}}\bigg{)}^{m_{l}}\prod_{t_{n,j}\in\Phi_{tb}}\bigg{(}\frac{m_{n}}{m_{n}+sp_{tb}g_{s}\eta_{n}D_{t_{n}}^{-\alpha_{n}}}\bigg{)}^{m_{n}}\bigg{]},$
(39)

where step (a) follows from we write the typical line with general lines
together, and proof completes by using the probability generating functional
of inhomogeneous PPP and PLCP. Just notice that for the interference from
general lines, we need to compute the lines intersect with $B({\rm
UAV,d_{db,db,nn}})$, ($\rho<d_{db,db,nn}$), and lines do not intersect with
$B({\rm UAV,d_{db,db,nl}})$, ($\rho>d_{db,db,nn}$), separately for the NLoS
interference, and the lines intersect with $B({\rm UAV,d_{db,db,nn}})$,
($\rho<d_{db,db,nn}$), and lines do not intersect with $B({\rm
UAV,d_{db,db,nl}})$, ($\rho>d_{db,db,nn}$), separately for the LoS
interference, where $B(a,b)$ denotes the ball centered at $a$ with radius $b$,
as shown in Fig. 11.

The Laplace transform $\mathcal{L}_{I_{tb}}$ follows a similar method and
omits the the interference from the typical line, thus, omitted here.

## References

  * [1] D. Cvitanić, “Drone applications in transportation,” in _2020 5th International Conference on Smart and Sustainable Technologies (SpliTech)_. IEEE, 2020, pp. 1–4.
  * [2] H. Shakhatreh, A. H. Sawalmeh, A. Al-Fuqaha, Z. Dou, E. Almaita, I. Khalil, N. S. Othman, A. Khreishah, and M. Guizani, “Unmanned aerial vehicles (UAVs): A survey on civil applications and key research challenges,” _IEEE Access_ , vol. 7, pp. 48 572–48 634, 2019.
  * [3] M. Khosravi, S. Enayati, H. Saeedi, and H. Pishro-Nik, “Multi-purpose drones for coverage and transport applications,” _IEEE Transactions on Wireless Communications_ , vol. 20, no. 6, pp. 3974–3987, 2021.
  * [4] Y. Qin, M. A. Kishk, and M.-S. Alouini, “Stochastic geometry-based analysis of multipurpose UAVs for package and data delivery,” _IEEE Internet of Things Journal_ , vol. 10, no. 5, pp. 4664–4676, 2023.
  * [5] L. Chen, M. A. Kishk, and M.-S. Alouini, “Dedicating cellular infrastructure for aerial users: Advantages and potential impact on ground users,” _IEEE Transactions on Wireless Communications_ , to appear.
  * [6] K.-J. Jung, K.-H. Park, Y.-C. Ko, and M.-S. Alouini, “Renewable energy-enabled cellular networks,” _Available at SSRN 3967953_ , 2021.
  * [7] “Terragraph virtual fiber for high speed fixed broadband,” Available online: https://terragraph.com/.
  * [8] S. Lalwani, “Economics of Terragraph backhaul for AT&T’s 5G network in San Jose,” Available online: https://rb.gy/elbmkx.
  * [9] A. Nordrum, “Facebook pushes networking tech: The company’s Terragraph technology will soon be available in commercial gear-[news],” _IEEE Spectrum_ , vol. 56, no. 4, pp. 8–9, 2019.
  * [10] Y. Zeng, J. Lyu, and R. Zhang, “Cellular-connected UAV: Potential, challenges, and promising technologies,” _IEEE Wireless Communications_ , vol. 26, no. 1, pp. 120–127, 2019.
  * [11] S. Chandrasekharan, K. Gomez, A. Al-Hourani, S. Kandeepan, T. Rasheed, L. Goratti, L. Reynaud, D. Grace, I. Bucaille, T. Wirth, and S. Allsopp, “Designing and implementing future aerial communication networks,” _IEEE Communications Magazine_ , vol. 54, no. 5, pp. 26–34, 2016.
  * [12] H. C. Nguyen, R. Amorim, J. Wigard, I. Z. KovaCs, T. B. Sorensen, and P. E. Mogensen, “How to ensure reliable connectivity for aerial vehicles over cellular networks,” _IEEE Access_ , vol. 6, pp. 12 304–12 317, 2018.
  * [13] M. M. Azari, F. Rosas, K.-C. Chen, and S. Pollin, “Ultra reliable UAV communication using altitude and cooperation diversity,” _IEEE Transactions on Communications_ , vol. 66, no. 1, pp. 330–344, 2017.
  * [14] ——, “Optimal UAV positioning for terrestrial-aerial communication in presence of fading,” in _2016 IEEE Global Communications Conference (GLOBECOM)_ , 2016, pp. 1–7.
  * [15] Y. Zeng and R. Zhang, “Energy-efficient UAV communication with trajectory optimization,” _IEEE Transactions on Wireless Communications_ , vol. 16, no. 6, pp. 3747–3760, 2017.
  * [16] J. Lyu, Y. Zeng, R. Zhang, and T. J. Lim, “Placement optimization of UAV-mounted mobile base stations,” _IEEE Communications Letters_ , vol. 21, no. 3, pp. 604–607, 2017.
  * [17] S. Zhang, H. Zhang, Q. He, K. Bian, and L. Song, “Joint trajectory and power optimization for UAV relay networks,” _IEEE Communications Letters_ , vol. 22, no. 1, pp. 161–164, 2018.
  * [18] Y. Zeng, R. Zhang, and T. J. Lim, “Throughput maximization for UAV-enabled mobile relaying systems,” _IEEE Transactions on Communications_ , vol. 64, no. 12, pp. 4983–4996, 2016.
  * [19] S. Zhang, Y. Zeng, and R. Zhang, “Cellular-enabled UAV communication: Trajectory optimization under connectivity constraint,” in _2018 IEEE International Conference on Communications (ICC)_ , 2018, pp. 1–6.
  * [20] ——, “Cellular-enabled UAV communication: A connectivity-constrained trajectory optimization perspective,” _IEEE Transactions on Communications_ , vol. 67, no. 3, pp. 2580–2604, 2019.
  * [21] Q. Wu, Y. Zeng, and R. Zhang, “Joint trajectory and communication design for multi-UAV enabled wireless networks,” _IEEE Transactions on Wireless Communications_ , vol. 17, no. 3, pp. 2109–2121, 2018.
  * [22] C. You and R. Zhang, “3D trajectory optimization in Rician fading for UAV-enabled data harvesting,” _IEEE Transactions on Wireless Communications_ , vol. 18, no. 6, pp. 3192–3207, 2019.
  * [23] G. Geraci, A. Garcia-Rodriguez, L. Galati Giordano, D. Lopez-Perez, and E. Bjornson, “Understanding UAV cellular communications: From existing networks to massive MIMO,” _IEEE Access_ , vol. 6, pp. 67 853–67 865, 2018.
  * [24] U. Challita, W. Saad, and C. Bettstetter, “Interference management for cellular-connected UAVs: A deep reinforcement learning approach,” _IEEE Transactions on Wireless Communications_ , vol. 18, no. 4, pp. 2125–2140, 2019.
  * [25] W. Mei, Q. Wu, and R. Zhang, “Cellular-connected UAV: Uplink association, power control and interference coordination,” _IEEE Transactions on wireless communications_ , vol. 18, no. 11, pp. 5380–5393, 2019.
  * [26] X. Pang, G. Gui, N. Zhao, W. Zhang, Y. Chen, Z. Ding, and F. Adachi, “Uplink precoding optimization for NOMA cellular-connected UAV networks,” _IEEE Transactions on Communications_ , vol. 68, no. 2, pp. 1271–1283, 2019\.
  * [27] Y. Li, A. H. Aghvami, and D. Dong, “Path planning for cellular-connected UAV: A DRL solution with quantum-inspired experience replay,” _IEEE Transactions on Wireless Communications_ , vol. 21, no. 10, pp. 7897–7912, 2022\.
  * [28] H. ElSawy, A. Sultan-Salem, M.-S. Alouini, and M. Z. Win, “Modeling and analysis of cellular networks using stochastic geometry: A tutorial,” _IEEE Communications Surveys Tutorials_ , vol. 19, no. 1, pp. 167–203, Firstquarter 2017.
  * [29] H. ElSawy, E. Hossain, and M. Haenggi, “Stochastic geometry for modeling, analysis, and design of multi-tier and cognitive cellular wireless networks: A survey,” _IEEE Communications Surveys Tutorials_ , vol. 15, no. 3, pp. 996–1019, Third 2013.
  * [30] Y. Zeng, R. Zhang, and T. J. Lim, “Wireless communications with unmanned aerial vehicles: Opportunities and challenges,” _IEEE Communications Magazine_ , vol. 54, no. 5, pp. 36–42, 2016.
  * [31] 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 Communications Surveys & Tutorials_, vol. 21, no. 3, pp. 2334–2360, 2019.
  * [32] Y. Qin, M. A. Kishk, and M.-S. Alouini, “Performance evaluation of UAV-enabled cellular networks with battery-limited drones,” _IEEE Communications Letters_ , vol. 24, no. 12, pp. 2664–2668, 2020.
  * [33] M. Alzenad and H. Yanikomeroglu, “Coverage and rate analysis for vertical heterogeneous networks (VHetNets),” _IEEE Transactions on Wireless Communications_ , vol. 18, no. 12, pp. 5643–5657, Dec. 2019.
  * [34] B. Galkin, J. Kibilda, and L. A. DaSilva, “A stochastic model for UAV networks positioned above demand hotspots in urban environments,” _IEEE Transactions on Vehicular Technology_ , vol. 68, no. 7, pp. 6985–6996, 2019.
  * [35] M. Alzenad, A. El-Keyi, and H. Yanikomeroglu, “3-D placement of an unmanned aerial vehicle base station for maximum coverage of users with different QoS requirements,” _IEEE Wireless Communications Letters_ , vol. 7, no. 1, pp. 38–41, 2017.
  * [36] Y. Qin, M. A. Kishk, and M.-S. Alouini, “Drone charging stations deployment in rural areas for better wireless coverage: Challenges and solutions,” _IEEE Internet of Things Magazine_ , vol. 5, no. 1, pp. 148–153, 2022.
  * [37] C. Saha, M. Afshang, and H. S. Dhillon, “Enriched $k$-tier HetNet model to enable the analysis of user-centric small cell deployments,” _IEEE Transactions on Wireless Communications_ , vol. 16, no. 3, pp. 1593–1608, March 2017.
  * [38] Y. Qin, M. A. Kishk, and M.-S. Alouini, “On the influence of charging stations spatial distribution on aerial wireless networks,” _IEEE Transactions on Green Communications and Networking_ , vol. 5, no. 3, pp. 1395–1409, 2021.
  * [39] M.-A. Lahmeri, M. A. Kishk, and M.-S. Alouini, “Laser-powered UAVs for wireless communication coverage: A large-scale deployment strategy,” _IEEE Transactions on Wireless Communications_ , pp. 1–1, 2022.
  * [40] ——, “Stochastic geometry-based analysis of airborne base stations with laser-powered UAVs,” _IEEE Communications Letters_ , vol. 24, no. 1, pp. 173–177, 2020.
  * [41] M. Banagar and H. S. Dhillon, “3D two-hop cellular networks with wireless backhauled UAVs: Modeling and fundamentals,” _IEEE Transactions on Wireless Communications_ , vol. 21, no. 8, pp. 6417–6433, 2022.
  * [42] J. Sabzehali, V. K. Shah, H. S. Dhillon, and J. H. Reed, “3D placement and orientation of mmWave-based UAVs for guaranteed LoS coverage,” _IEEE Wireless Communications Letters_ , vol. 10, no. 8, pp. 1662–1666, 2021\.
  * [43] Y. Qin, M. A. Kishk, and M.-S. Alouini, “Everything on lamp posts: Towards multi-functional street furniture,” Available online: shorturl.at/bglmR.
  * [44] H. S. Dhillon and V. V. Chetlur, “Poisson line Cox process: Foundations and applications to vehicular networks,” _Synthesis Lectures on Learning, Networks, and Algorithms_ , vol. 1, no. 1, pp. 1–149, 2020.
  * [45] A. Al-Hourani, S. Kandeepan, and S. Lardner, “Optimal LAP altitude for maximum coverage,” _IEEE Wireless Communications Letters_ , vol. 3, no. 6, pp. 569–572, 2014.
  * [46] H. Alzer, “On some inequalities for the incomplete Gamma function,” _Mathematics of Computation_ , vol. 66, no. 218, pp. 771–778, 1997.
  * [47] S. Krishnan and H. S. Dhillon, “Spatio-temporal interference correlation and joint coverage in cellular networks,” _IEEE Transactions on Wireless Communications_ , vol. 16, no. 9, pp. 5659–5672, 2017.
  * [48] M. Haenggi, _Stochastic geometry for wireless networks_. Cambridge University Press, 2012.
  * [49] H. Lu, P. Zhao, H. Hu, L. Zeng, K. S. Wu, and D. Lv, “Transport infrastructure and urban-rural income disparity: A municipal-level analysis in China,” _Journal of Transport Geography_ , vol. 99, p. 103292, 2022.