Performance Characterization of Canonical Mobility Models in Drone Cellular Networks

TL;DR

This paper rigorously analyzes the performance of canonical drone mobility models in cellular networks, providing exact expressions for key metrics and establishing the straight line model as a performance lower bound.

Contribution

It introduces a unified mathematical framework for modeling drone mobility in cellular networks and compares different models, highlighting the straight line model as a performance lower bound.

Findings

Exact expressions for average and session rates

Straight line model provides a lower performance bound

Unified approach for various mobility and service models

Abstract

In this paper, we characterize the performance of several canonical mobility models in a drone cellular network in which drone base stations (DBSs) serve user equipments (UEs) on the ground. In particular, we consider the following four mobility models: (i) straight line (SL), (ii) random stop (RS), (iii) random walk (RW), and (iv) random waypoint (RWP), among which the SL mobility model is inspired by the simulation models used by the third generation partnership project (3GPP) for the placement and trajectory of drones, while the other three are well-known canonical models (or their variants) that offer a useful balance between realism and tractability. Assuming the nearest-neighbor association policy, we consider two service models for the UEs: (i) UE independent model (UIM), and (ii) UE dependent model (UDM). While the serving DBS follows the same mobility model as the other DBSs in the UIM, it is assumed to fly towards the UE of interest in the UDM and hover above its location after reaching there. The main contribution of this paper is a unified approach to characterize the point process of DBSs for all the mobility and service models. Using this, we provide exact mathematical expressions for the average received rate and the session rate as seen by the typical UE. Further, using tools from calculus of variations, we concretely demonstrate that the simple SL mobility model provides a lower bound on the performance of other general mobility models (including the ones in which drones follow curved trajectories) as long as the movement of each drone in these models is independent and identically distributed (i.i.d.). To the best of our knowledge, this is the first work that provides a rigorous analysis of key canonical mobility models for an infinite drone cellular network and establishes useful connections between them.

Figures

Equations

$$\mathrm{SIR}(t)=\frac{P{h}_0(t)r_0(t{)}^{-\alpha }}{I(t)}.$$

$$\mathrm{SR}(T)=\frac{1}{T}{\int }_0^{T}R(t)\mathrm{d}t.$$

$$\lambda (t;u_{\mathbf{x}},u_0)=\{\begin{array}{cc}{\lambda }_0& u_{\mathbf{x}}>u_0(t)\\ 0& u_{\mathbf{x}}\le u_0(t)\end{array}.$$

$$\lambda (t;u_{\mathbf{x}},u_0)={\lambda }_0\{\begin{array}{ll}1& u_0+vt\le u_{\mathbf{x}}\\ \beta (t,u_{\mathbf{x}},u_0)& |u_0-vt|\le u_{\mathbf{x}}\le u_0+vt\\ \beta (t,u_{\mathbf{x}},u_0)1\left(t>\frac{u_0}{v}\right)& 0\le u_{\mathbf{x}}\le |u_0-vt|\end{array},$$

$$\beta (t,u_{\mathbf{x}},u_0)=1-F_L(u_0-u_{\mathbf{x}};t)-{\int }_{|u_{\mathbf{x}}-u_0|}^{\min \{vt,u_{\mathbf{x}}+u_0\}}{f}_L(l;t)\frac{1}{\pi }{\cos }^{-1}\left(\frac{l^2+u_{\mathbf{x}}^2-u_0^2}{2l{u}_{\mathbf{x}}}\right)\mathrm{d}l.$$

$$\lambda (t;u_{\mathbf{x}},u_0)={\lambda }_0\{\begin{array}{ll}1& u_0+vt\le u_{\mathbf{x}}\\ \frac{1}{\pi }{\cos }^{-1}\left(\frac{u_0^2-u_{\mathbf{x}}^2-v^2{t}^2}{2u_{\mathbf{x}}vt}\right)& |u_0-vt|\le u_{\mathbf{x}}\le u_0+vt\\ 1\left(t>\frac{u_0}{v}\right)& 0\le u_{\mathbf{x}}\le |u_0-vt|\end{array}.$$

$$\beta (t,u_{\mathbf{x}},u_0)=1-1(u_0-u_{\mathbf{x}}-vt)-{\int }_{|u_{\mathbf{x}}-u_0|}^{\min \{vt,u_{\mathbf{x}}+u_0\}}\delta (l-vt)\frac{1}{\pi }{\cos }^{-1}\left(\frac{l^2+u_{\mathbf{x}}^2-u_0^2}{2l{u}_{\mathbf{x}}}\right)\mathrm{d}l.$$

$$\beta (t,u_{\mathbf{x}},u_0)$$

$$F_L(l;t)$$

$$f_L(l;t)$$

$$\beta (t,u_{\mathbf{x}},u_0)=1-F_R(u_0-u_{\mathbf{x}})-{\int }_{|u_{\mathbf{x}}-u_0|}^{u_{\mathbf{x}}+u_0}{f}_R(l)\frac{1}{\pi }{\cos }^{-1}\left(\frac{l^2+u_{\mathbf{x}}^2-u_0^2}{2l{u}_{\mathbf{x}}}\right)\mathrm{d}l.$$

$$\displaystyle\scalebox{0.85}{$S_{n-1}\!=\!\sum_{i=1}^{n-1}\!R_{i}$},\hskip 12.80365pt\scalebox{0.85}{$Z_{n-1}\!=\!\sqrt{\!\left(\sum_{i=1}^{n-1}\!R_{i}\cos(\Theta_{i})\right)^{2}\!\!+\!\left(\sum_{i=1}^{n-1}\!R_{i}\sin(\Theta_{i})\right)^{2}}$},\hskip 12.80365pt\scalebox{0.85}{$\Psi_{n-1}\!=\!\tan^{-1}\!\left(\!\frac{\sum_{i=1}^{n-1}R_{i}\sin(\Theta_{i})}{\sum_{i=1}^{n-1}R_{i}\cos(\Theta_{i})}\!\right)$}.$$

$$f_{S_n,Z_n}(s,z)=\frac{4z}{(2\pi {)}^n}\underset{\mathcal{R}}{\int \cdots \int }\frac{f_R(s-J_n){\prod }_{i=1}^{n-1}{f}_R(x_i)}{\sqrt{\left(z^2-(s-J_n-K_n{)}^2\right)\left((s-J_n+K_n{)}^2-z^2\right)}}\mathrm{d}\mathbf{x}\mathrm{d}\mathbf{\xi },$$

$$f_{S_2,Z_2}(s,z)=\frac{2z}{\pi \sqrt{s^2-z^2}}{\int }_{\frac{s-z}{2}}^{\frac{s+z}{2}}\frac{f_R(x)f_R(s-x)}{\sqrt{z^2-(2x-s{)}^2}}\mathrm{d}x,$$

$$F_L(l;t)=$$

$$\displaystyle\sum_{n=2}^{\infty}\int_{\frac{vt-l}{2}}^{vt}\int_{|l-(vt-s)|}^{\min\{s,\,l+(vt-s)\}}\scalebox{1.0}{$f_{S_{n-1},Z_{n-1}}(s,z)(1\!-\!F_{R}(vt\!-\!s))\frac{1}{\pi}\cos^{-1}\left(\frac{z^{2}+(vt-s)^{2}-l^{2}}{2z(vt-s)}\right)\,{\rm d}z\,{\rm d}s,$}$$

$$f_L(l;t)=$$

$$F_L(l;t|n=2)$$

$$f_L(l;t|n=2)$$

$$F_L(l;t)=$$

$$\sum _{n=2}^{\infty }{\int }_0^{vt}{\int }_0^{\min \{s,l\}}{f}_{S_{n-1},Z_{n-1}}(s,z)\left(F_{W_{n-2}}(t-\frac{s}{v})-F_{W_{n-1}}(t-\frac{s}{v})\right)\mathrm{d}z\mathrm{d}s+$$

$$\sum _{n=2}^{\infty }{\int }_{vt-l}^{vt}{\int }_0^{l-(vt-y)}{\int }_0^{\frac{y}{v}}{f}_{W_{n-1}}(w)f_{S_{n-1},Z_{n-1}}(y-vw,z)(1-F_R(vt-y))\mathrm{d}w\mathrm{d}z\mathrm{d}y+$$

$$\sum _{n=2}^{\infty }{\int }_{\frac{vt-l}{2}}^{vt}{\int }_{|l-(vt-y)|}^{\min \{y,l+(vt-y)\}}{\int }_0^{\frac{y}{v}}{f}_{W_{n-1}}(w)f_{S_{n-1},Z_{n-1}}(y-vw,z)(1-F_R(vt-y))\times$$

$$\frac{1}{\pi }{\cos }^{-1}\left(\frac{z^2+(vt-y{)}^2-l^2}{2z(vt-y)}\right)\mathrm{d}w\mathrm{d}z\mathrm{d}y,$$

$$f_L(l;t)=$$

$$\displaystyle\frac{2l}{\pi}\sum_{n=2}^{\infty}\int_{\frac{vt-l}{2}}^{vt}\int_{|l-(vt-y)|}^{\min\{y,\,l+(vt-y)\}}\!\!\!\!\int_{0}^{\frac{y}{v}}\scalebox{1.03}{$\frac{f_{W_{n-1}}(w)f_{S_{n-1},Z_{n-1}}(y-vw,z)(1-F_{R}(vt-y))}{\sqrt{\left[l^{2}-(z-(vt-y))^{2}\right]\left[(z+(vt-y))^{2}-l^{2}\right]}}\,{\rm d}w\,{\rm d}z\,{\rm d}y,$}$$

$$F_L(l;t|n=2)=$$

$$\displaystyle\int_{0}^{t-\frac{l}{v}}\int_{\frac{vt-vw-l}{2}}^{\frac{vt-vw+l}{2}}\scalebox{1.02}{$f_{R}(r)f_{W_{1}}(w)(1-F_{R}(vt-vw-r))\frac{1}{\pi}\cos^{-1}\left(\frac{r^{2}+(vt-vw-r)^{2}-l^{2}}{2r(vt-vw-r)}\right)\,{\rm d}r\,{\rm d}w,$}$$

$$f_L(l;t|n=2)=$$

$$\frac{2l}{\pi }{\int }_0^{t-\frac{l}{v}}{\int }_{\frac{vt-vw-l}{2}}^{\frac{vt-vw+l}{2}}\frac{f_R(r)f_{W_1}(w)(1-F_R(vt-vw-r))}{\sqrt{\left[(vt-vw{)}^2-l^2\right]\left[l^2-(2r-(vt-vw){)}^2\right]}}\mathrm{d}r\mathrm{d}w.$$

$$R(t)$$

$$\displaystyle\mathcal{L}_{I(t)}(s\bigr{\rvert}\mathbf{x}_{0}(t))$$

$$\displaystyle\mathbb{P}\left[{\rm SIR}(t)\geq\gamma\bigr{\rvert}\mathbf{x}_{0}(t)\right]\!$$

$$\displaystyle\overset{(b)}{=}\mathbb{E}\left[\sum_{k=0}^{m_{0}-1}\frac{(m_{0}\gamma r_{0}^{\alpha}(t)I(t))^{k}}{k!}{\rm e}^{-m_{0}\gamma r_{0}^{\alpha}(t)I(t)}\middle|\mathbf{x}_{0}(t)\right]=\sum_{k=0}^{m_{0}-1}\frac{(-s)^{k}}{k!}\frac{\partial^{k}}{\partial s^{k}}\mathcal{L}_{I(t)}(s\bigr{\rvert}\mathbf{x}_{0}(t))\biggr{|}_{s=m_{0}\gamma r_{0}^{\alpha}(t)},$$

$$\displaystyle\mathcal{L}_{I(t)}(s\bigr{\rvert}\mathbf{x}_{0}(t))$$

$$\displaystyle\overset{(a)}{=}\mathbb{E}\!\!\left[\prod_{\mathbf{x}(t)\in\Phi_{\rm D}^{\prime}(t)}\!\!\!\!\scalebox{0.9}{$\left(1+\frac{sr_{\mathbf{x}}(t)^{-\alpha}}{m}\right)^{\!-m}$}\middle|u_{0}(t)\right]\!\overset{(b)}{=}\!\exp\left[-2\pi\!\!\int_{0}^{\infty}\!\!\!\!\!\scalebox{0.9}{$u_{\mathbf{x}}(t)\lambda(t;u_{\mathbf{x}},u_{0})\left(1\!-\!\left(1\!+\!\frac{sr_{\mathbf{x}}(t)^{\!-\alpha}}{m}\right)^{-m}\right){\rm d}u_{\mathbf{x}}(t)$}\right],$$

$$\displaystyle\int_{0}^{\infty}\log(1+\gamma)f_{\Gamma}(\gamma;t)\,{\rm d}\gamma=\int_{0}^{\infty}\int_{0}^{\infty}\frac{2\pi\lambda_{0}u_{0}{\rm e}^{-\pi\lambda_{0}u_{0}^{2}}}{1+\gamma}\mathbb{P}\left[{\rm SIR}(t)\geq\gamma\bigr{\rvert}\mathbf{x}_{0}(t)\right]\,{\rm d}u_{0}\,{\rm d}\gamma,$$

$${\lambda }_1(t;u_{\mathbf{x}},u_0)$$

Full text

I Introduction

Wireless networks are all set to undergo a major transformation from being predominantly terrestrial to the ones that will have an elaborate and dynamic aerial component in the form of drone networks [3]. Owing to their deployment flexibility, drones are currently being considered for a variety of use cases, such as acting as mobile relays to expand the coverage of cellular networks, providing network connectivity for public safety applications, and setting up temporary networks in the times of natural disasters or large social gatherings. While the flexibility offered by the mobility of drones is highly appealing to the network designers, it also adds an entirely new dimension to the system design that was not present in the traditional terrestrial networks. In particular, while mobility is known to have a fundamental impact on the system-level performance of wireless networks, e.g., see [4], the research focus in this direction has traditionally been on the setting in which the UEs are mobile while the base stations (BSs) are static [5, 6]. This is clearly not the case in drone-assisted communication networks where some mobile drones may act as BSs [7]. Not surprisingly, the support of mobile drones, either as UEs or BSs, has also been recently explored in 3GPP studies [8, 9]. In general, there has been an increasing interest in the community to characterize the effect of drone mobility on the design and performance of drone-assisted cellular networks. Due to the irregularity of the drone placements and their trajectories, it is expected that powerful tools from stochastic geometry and point process theory could be leveraged for accurate modeling and tractable analysis of drone networks. Inspired by this, the main focus of this paper is to present a unified performance analysis of a DBS network under various mobility models that involves a novel characterization of the underlying point processes as a function of time.

I-A Related Works

This paper builds on the following two key lines of research: (i) stochastic geometry for drone networks, and (ii) mobility models in wireless networks. Although sparse, there are some works that lie at their intersection and will be discussed below after establishing the prior art for each research direction separately.

Stochastic Geometry for Drone Networks. Owing to its ability to capture irregularity in the placement and movement of drones, stochastic geometry has recently found many applications in the performance analysis of drone networks [10, 11, 12, 13, 14, 15, 16]. In [10], the authors considered a finite network of DBSs distributed as a uniform binomial point process (BPP) [17] and derived the coverage probability of the network. In [12], the problem of spectrum sharing for a network of drone small cells as an underlay to a conventional cellular network has been studied. On the similar lines, the authors of [13] investigated the coexistence of BSs and DBSs using probabilistic line-of-sight (LoS) and non-line-of-sight (NLoS) propagation models [18], where the locations of BSs and DBSs are modeled as a superposition of a Poisson point process (PPP) and a BPP. The work presented in [14] considered a network of DBSs modeled as a PPP serving ground UEs. In particular, incorporating LoS and NLoS propagation models, the authors derived approximations for the coverage probability and the received rate in the network. In [15], the coexistence between a single DBS and an underlaid device-to-device (D2D) network has been analyzed in terms of coverage probability and rate. Motivated by the requirement to provide coverage in a post-disaster scenario, the authors in [16] analyzed the performance of a DBS network serving clustered UEs.

Mobility Models in Wireless Networks. Mobility modeling is a well-established area of research in wireless networks [19, 20, 21, 22]. Perhaps the simplest mobility model is the one in which nodes move along straight lines in random directions with a constant speed. Variants of this simple model have been used extensively in the literature, e.g., see [23, 24, 21, 25, 8]. Although this model may appear simplistic, it is known to provide performance bounds and useful insights in wireless networks [23]. In fact, this model has also been used recently to model drone mobility in 3GPP studies related to drone networks [8]. Among other mobility models that have been studied in wireless networks (such as RW, RWP, random direction, Brownian motion, Levy walk, and Gauss-Markov), RW and RWP mobility models have been more popular because of their tractability [26, 5, 27, 28, 19]. Hence, we also study these two models in this paper.

In the RW model, the mobile nodes can change their directions [19, 26, 24]. In particular, each node selects a uniformly random direction and a random speed in each time slot and moves along a straight line using the selected direction and speed. Upon arrival at the destination, it repeats this procedure. As a generalization of the RW model, the finite RWP model was proposed in [29] by adding a random pause time between direction and/or speed changes in the RW model. Specifically, the mobility process in the finite RWP model starts with a random pause time at the initial waypoint of a node in a finite region. Then, a destination waypoint is selected randomly in that region and the node moves to this destination waypoint along a straight line with a random speed. Upon its arrival, the node pauses for another random time and repeats this procedure. Due to some major drawbacks of this model, such as nonuniform distribution of the nodes [5, 30, 31], the authors in [6] extended the finite RWP to an infinite RWP model, where the nodes are allowed to move over the entire plane. In particular, each node first pauses for a random time at its initial waypoint and determines the next waypoint by choosing a uniformly random direction and a random transition distance. It then moves towards the chosen waypoint at a random speed, pauses at this new waypoint for another random time, and repeats this procedure.

Mobile Drone Networks. Among many works on the analysis of drone networks, only a handful of them considered mobile drones. Considering a mobile network of DBSs that are initially modeled as a BPP, the authors in [32] used results from [10] to design stochastic trajectory processes for the mobility of DBSs that provide the same coverage as given in [10]. In [33] and [34], the authors modeled the motion of DBSs in a 3D finite network by the finite RWP and RW models, and using the results of [5], derived the coverage probability of the network.

I-B Contributions

This paper develops a unified approach for the performance analysis of drone cellular networks under several key canonical models. In particular, we model the initial locations of a mobile network of DBSs as a homogeneous PPP operating at a constant height that serve UEs on the ground. The serving DBS is selected based on the nearest-neighbor association policy and all the other DBSs are regarded as interferers. We consider four canonical mobility models for the interfering DBSs, i.e., SL, RS, RW, and RWP, and two service models for the serving DBS: (i) UIM, in which the serving DBS follows the same mobility model as the interfering DBSs, and (ii) UDM, in which the serving DBS moves towards the typical UE at a constant height and keeps hovering above its location after reaching there until its transmission is completed. For this setup, our key contributions are described next.

I-B1 Distributional Properties of Canonical Mobility Models

We derive various distributional properties for all the mobility models. Most notably, for the RW and RWP, we derive the joint distribution of the net displacement of a DBS and its total traveled distance at each waypoint. Using this joint distribution, we compute the distribution of the net displacement of a DBS at any time $t$. We also provide insightful asymptotic results for these distributions.

I-B2 Unified Framework for Characterizing the Point Process of DBSs

We present a novel characterization of the point process of DBSs as seen by the typical UE at any time $t$ for all the service and mobility models. Since the displacement of each DBS is independent from the others, we apply displacement theorem from stochastic geometry along with the net displacement results described above to characterize these point processes. Using this, we characterize the aggregate interference and key performance metrics, such as average and session rates, as seen by the typical UE for all the service and mobility models.

I-B3 Establishing Connections among Mobility Models

Borrowing tools from the calculus of variations, we establish meaningful connections among our mobility models in terms of network performance. Specifically, we prove that the SL mobility model provides a lower bound on the average received rate over the space of all i.i.d. mobility models in which the drones can also follow curved trajectories (as long as they are i.i.d.).

II System Model

II-A Spatial Setup

We consider a network of mobile DBSs serving UEs on the ground. We assume that the DBSs are located at a constant height $h$ from the ground and the temporal evolution of the DBS locations is modeled as the sequence of point processes $\Phi_{\rm D}(t)\subset{\mathbb{R}}^{2}$, indexed by $t\in{\mathbb{R}}^{+}$. Further, we assume that the initial locations of the DBSs are distributed as a homogeneous PPP with density $\lambda_{0}$, i.e., $\Phi_{\rm D}(0)\sim{\rm PPP}(\lambda_{0})$. Terrestrial UEs are distributed as an independent homogeneous PPP $\Phi_{\rm U}$ on the ground. We also assume that the origin ${\mathbf{o}}=(0,0,0)$ of the 3D coordinate system is located on the ground, which is assumed to be aligned with the $xy$-plane. Throughout the paper, we refer to the $z=h$ plane as the DBS plane. In this setup, the projection of ${\mathbf{o}}$ onto the DBS plane is denoted by ${\mathbf{o}}^{\prime}=(0,0,h)$. Without loss of generality, the analysis will be performed for the typical UE placed at ${\mathbf{o}}$. As shown in Fig. 1, the distances of a DBS at time $t$ located at ${\mathbf{x}}(t)\in\Phi_{\rm D}(t)$ from $\bf{o^{\prime}}$ and $\bf{o}$ are denoted by $u_{\mathbf{x}}(t)=\|{\mathbf{x}}(t)-\bf{o^{\prime}}\|$ and $r_{\mathbf{x}}(t)=\|{\mathbf{x}}(t)-{\bf o}\|=\sqrt{u_{\mathbf{x}}(t)^{2}+h^{2}}$, respectively. Moreover, the location of the nearest DBS to $\bf{o^{\prime}}$ and its corresponding distance at time $t$ are denoted by ${\mathbf{x}}_{0}(t)$ and $u_{0}(t)$, respectively. Thus, the distance of the closest DBS to $\bf{o}$ at time $t$ is $r_{0}(t)=\sqrt{u_{0}(t)^{2}+h^{2}}$. For simplicity, we drop the time index $t$ for the distances defined at $t=0$, i.e., $u_{0}\triangleq u_{0}(0)$, $r_{0}\triangleq r_{0}(0)$, $u_{\mathbf{x}}\triangleq u_{\mathbf{x}}(0)$, and $r_{\mathbf{x}}\triangleq r_{\mathbf{x}}(0)$.

II-B Service and Mobility Models

In this paper, we assume that each UE connects to its nearest DBS. For the typical UE, this DBS is called the serving DBS while all the other DBSs are regarded as interfering DBSs. In order to motivate the drone mobility models studied in this paper, we first take cue from the mobility model used in the 3GPP studies related to drone networks [8]. In this model, the drones start their movement at randomly selected locations in the network and then move at a constant speed and height in straight lines and in uniformly random directions independently from each other for the entire duration of the simulation. It turns out that this simple enough model is already considered sufficient to capture key effects of drone mobility on the system-level performance. One can of course generalize this straight-line mobility model to arrive at classical canonical models, such as RW and RWP, which offer a useful balance between realism and tractability. Since all these models are important in their own right, we will develop a unified approach to analyze their performance jointly. Before that, we formally define the four mobility models considered in this paper next.

Definition 1**.**

(Mobility Models).* We define the mobility models used in this paper as follows. Without loss of generality, we choose the $x$-axis as the reference for measuring all angles.*

SL*: DBSs move constantly in random directions $\Theta\sim U[0,2\pi)$ on straight lines with a constant speed $v$, independently from each other.* 2. 2.

RS*: Each DBS moves a random distance $R\sim f_{R}(.)$ in a random direction $\Theta\sim U[0,2\pi)$ on a straight line with a constant speed $v$, independently from the other DBSs, and then stops and hovers over the stopping location.* 3. 3.

RW*: Each DBS selects a random direction $\Theta\!\sim\!U[0,2\pi)$, independently of the other DBSs, and moves a random distance $R\sim f_{R}(.)$ in this direction with a constant speed $v$. Upon its arrival, it selects another random direction $\Theta$ and distance $R$ and repeats this procedure.* 4. 4.

RWP*: In the beginning, each DBS hovers for a random time $T\!\sim\!f_{T}(.)$ at its initial location. It then selects a random direction $\Theta\!\sim\!U[0,2\pi)$, independently of the other DBSs, and moves a random distance $R\sim f_{R}(.)$ in this direction with a constant speed $v$. Upon its arrival, it hovers for another random time $T$, sampled independently from $f_{T}(.)$, selects another random direction $\Theta$ and distance $R$, and repeats this procedure.*

Remark 1**.**

All the mobility models considered in this paper lie in the space of i.i.d. mobility models, where the drone trajectories are chosen independently of each other from the same common distribution. As will be clear in the sequel, some of our initial results hold for all i.i.d. models (including the ones with curved trajectories).

Definition 2**.**

(Flight).* In the RW and RWP mobility models, a flight is defined as the distance traveled by a DBS between two consecutive stop points.*

For the serving DBS, we consider the following two service models.

UIM: The serving DBS follows the same mobility model as the interfering DBSs, independent of the typical UE location. 2. 2.

UDM: The serving DBS moves towards $\bf{o}^{\prime}$ in the DBS plane and keeps hovering at this location until its transmission to the typical UE is completed.

Note that in all the aforementioned models, the speeds of the serving and interfering DBSs are assumed to be the same. Moreover, since we are not considering any dependency across the user locations, DBS trajectories will be independent of each other in both the UIM and the UDM.

Remark 2**.**

In the UIM, since all DBSs move in different directions based on the defined mobility models, handover may occur. On the other hand, in the UDM, no matter what the mobility model is for the interfering DBSs, as long as the interfering DBSs move with the same speed as the serving DBS, the serving DBS will remain the closest DBS to the typical UE. Hence, based on our association policy, the serving DBS will not change in our mobility models, and thus, handover will not occur. Furthermore, UDM can be considered as a best-case model from the perspective of minimizing the distance between the typical UE and its serving DBS.

II-C Channel Model

We assume that all DBSs transmit with the same power level $P$ at all times. The received power at the typical UE from the serving DBS is assumed to be $Ph_{0}(t)r_{0}(t)^{-\alpha}$, where $h_{0}(t)$ represents the small-scale fading power gain between the typical UE and the serving DBS and $\alpha>2$ is the path-loss exponent. Likewise, the interference power is $I(t)=\sum_{{\mathbf{x}}(t)\in\Phi_{\rm D}^{\prime}(t)}Ph_{\mathbf{x}}(t)r_{\mathbf{x}}(t)^{-\alpha}$, where $\Phi_{\rm D}^{\prime}(t)\equiv\Phi_{\rm D}(t)\backslash{\mathbf{x}}_{0}(t)$ represents the point process of interfering DBSs and $h_{\mathbf{x}}(t)$ is the small-scale fading power gain between the typical UE and the interfering DBSs. Since the air-to-ground links may experience various fading scenarios, Nakagami-$m$ fading is used here to capture a large class of fading environments. The Nakagami-$m$ fading parameter is assumed to be $m_{0}$ and $m_{\mathbf{x}}$ for the serving and interfering links, respectively. Thus, the channel fading power gains $h_{0}(t)$ and $h_{\mathbf{x}}(t)$ follow gamma distributions with probability density function (pdf) $f_{H}(h)=\frac{m^{m}}{\Gamma(m)}h^{m-1}{\rm e}^{-mh}$, where $\Gamma(x)=\int_{0}^{\infty}t^{x-1}{\rm e}^{-t}\,{\rm d}t$ is the gamma function. For the serving and interfering links, we assume $m=m_{0}$ and $m=m_{\mathbf{x}}$, respectively, and we consider integer values for $m_{0}$ and $m_{\mathbf{x}}$ for mathematical tractability.

II-D Metrics of Interest

The received signal-to-interference ratio (${\rm SIR}$) at the typical UE is defined as

[TABLE]

We now define our ${\rm SIR}$-based performance metrics as follows.

Average rate: Average received rate is given as $R(t)={\mathbb{E}}[\log\left(1+{\rm SIR}(t)\right)]$, where the expectation is taken over the PPP $\Phi_{\rm D}(t)$ and the trajectories. This is essentially the average rate experienced by the typical UE at time $t$ across different network and trajectory realizations.

Session rate: This metric is defined as the average received rate by the typical UE at each session of duration $T$. Mathematically speaking, we have

[TABLE]

III Point Process of Interferers

In this section, we characterize the temporal evolution of the point process of interferers for all the mobility models described in Section II-B. We start our analysis by first considering the UIM in the following lemma.

Lemma 1**.**

Let $\Phi$ be a homogeneous PPP with density $\lambda_{0}$. If all the points of $\Phi$ are independently displaced based on the four mobility models mentioned in II-B, then the displaced points at every time $t$ form another homogeneous PPP $\Psi$ with the same density $\lambda_{0}$.

Proof:

Based on the displacement theorem [35], we need to argue that the displaced distances of DBSs at every time $t$ are i.i.d. and also independent of their original locations in $\Phi$. For the first two mobility models, this is clear from our model constructions. For the RW and RWP models, note that different flights are independent from the locations of DBSs and the overall displaced distance at time $t$ is a function of these flights. Hence, the overall displacement of DBSs are independent from their locations in the original PPP $\Phi$. Therefore, $\Psi$ is distributed as a homogeneous PPP with density $\lambda_{0}$. This completes the proof. ∎

In the UIM, the serving DBS and all the interfering DBSs are displaced in uniformly random directions based on our different mobility models. Thus, we can infer from Lemma 1 that the network of all DBSs at any time $t$ will remain a homogeneous PPP with density $\lambda_{0}$. As a result, the network of interfering DBSs will be distributed as an inhomogeneous PPP in the DBS plane with density

[TABLE]

Note that although the serving distance $u_{0}(t)$ varies over time, its distribution does not change. In the UDM, it is clear from our construction that $\Phi_{\rm D}^{\prime}(0)$ is an inhomogeneous PPP with density given by (3), which introduces an exclusion zone, ${\mathcal{X}}=b({\bf o^{\prime}},u_{0})$, for the interfering DBSs in the DBS plane, where $b({\bf o},r)$ is a disc of radius $r$ centered at ${\bf o}$. Using displacement theorem, we observe that $\Phi_{\rm D}^{\prime}(t)$ remains an inhomogeneous PPP for the UDM as well. In the next lemma, we provide a mathematical characterization of $\Phi_{\rm D}^{\prime}(t)$ in the UDM at every time $t$.

Lemma 2**.**

(Point Process of Interferers).* Consider the UDM with an i.i.d. mobility model. Let ${\mathbf{x}}(0)$ and ${\mathbf{x}}(t)$ denote the initial location and the location at time $t$ of an interfering DBS, respectively. Define $L(t)=\|{\mathbf{x}}(0)-{\mathbf{x}}(t)\|$ as the net displacement of this interfering DBS until time $t$ and denote its cumulative distribution function (cdf) and pdf by $F_{L}(l;t)$ and $f_{L}(l;t)$, respectively. Then the network of interfering DBSs will be an inhomogeneous PPP with density*

[TABLE]

where ${\bf 1}(.)$ is the indicator function and

[TABLE]

Proof:

See Appendix -A. ∎

Remark 3**.**

We observe the following directly from Lemma 2: (i) (4) is continuous at boundaries, i.e., at $u_{\mathbf{x}}=|u_{0}\pm vt|$, (ii) as $u_{0}\to 0$, we get $\beta(t,u_{\mathbf{x}},u_{0})\to 1$, i.e., the network of interfering DBSs becomes homogeneous, and (iii) as $u_{0}\to 0$ or $t\to 0$, (3) and (4) become identical.

The density derived in Lemma 2 is valid for all i.i.d. mobility models where DBSs move with the same constant speed. Hence, for mobility models in this paper, we only need to characterize the distribution of the net displacement of each DBS at every time $t$. For the SL model, we provide this distribution, and thus, the density of the network of interferers, in the next corollary.

Corollary 1**.**

(Point Process of Interferers for the SL Model).* When the interferers move based on the SL mobility model, the network of interferers will be an inhomogeneous PPP with density*

[TABLE]

Proof:

In the SL mobility model, we have $L(t)=vt$, and thus, $f_{L}(l;t)=\delta(l-vt)$ and $F_{L}(l;t)={\bf 1}(l-vt)$, where $\delta(.)$ is the Dirac delta function. Hence, $\beta(t,u_{\mathbf{x}},u_{0})$ can be written as

[TABLE]

We need to consider two cases:

• •

$0\leq u_{\mathbf{x}}\leq|u_{0}-vt|$

In this case, since $\lambda(t;u_{\mathbf{x}},u_{0})=0$ when $u_{0}\geq vt$, we need to only evaluate $\beta(t,u_{\mathbf{x}},u_{0})$ when $u_{0}<vt$. This gives $0\leq u_{\mathbf{x}}\leq vt-u_{0}$, which makes the integral term in (7) zero. Since $u_{0}-u_{\mathbf{x}}\leq u_{0}+u_{\mathbf{x}}\leq vt$, the term for the indicator function in (7) will become zero as well. Hence, $\beta(t,u_{\mathbf{x}},u_{0})={\bf 1}\left(t>\frac{u_{0}}{v}\right)$.

• •

$|u_{0}-vt|\leq u_{\mathbf{x}}\leq u_{0}+vt$

In this case, we have the triangle inequality for the triple $(u_{\mathbf{x}},u_{0},vt)$, which yields $|u_{0}-u_{\mathbf{x}}|\leq vt\leq u_{0}+u_{\mathbf{x}}$. Hence,

[TABLE]

This completes the proof. ∎

Remark 4**.**

We can make sense of (6) intuitively by considering the exclusion zone ${\mathcal{X}}$ and different regions illustrated in Fig. 2 as follows: (i) ${\mathcal{R}}_{1}=\{u_{\mathbf{x}}\geq u_{0}+vt\}$: No DBS that was initially in ${\mathcal{R}}_{1}$ can enter ${\mathcal{X}}$ within time $t$, and thus, $\lambda(t;u_{\mathbf{x}},u_{0})=\lambda_{0}$, (ii) ${\mathcal{R}}_{2}=\{0\leq u_{\mathbf{x}}\leq u_{0}-vt|vt\leq u_{0}\}$: The serving DBS is in motion towards $\mathbf{o}^{\prime}$ and since there is no interfering DBS in ${\mathcal{R}}_{2}$, we have $\lambda(t;u_{\mathbf{x}},u_{0})=0$, and (iii) ${\mathcal{R}}_{3}=\{0\leq u_{\mathbf{x}}\leq vt-u_{0}|vt>u_{0}\}$: In this region, we can use the same method introduced in the proof of Lemma 2 as follows. We first calculate the density contributed by points that would have originally fallen in ${\mathcal{X}}$ and then subtract it from $\lambda_{0}$ to get the density of the network of interfering DBSs. Since all the DBSs that are initially inside ${\mathcal{X}}$ will leave ${\mathcal{R}}_{3}$ after the displacement of $vt$, the density contributed by ${\mathcal{X}}$ in ${\mathcal{R}}_{3}$ is [math], which gives $\lambda(t;u_{\mathbf{x}},u_{0})=\lambda_{0}$. However, for non-linear i.i.d. mobility models, such as RW or RWP, where direction changes are allowed during flights, there is a possibility that some of the DBSs that are initially inside ${\mathcal{X}}$ do not leave ${\mathcal{R}}_{3}$ after the displacement of $vt$. This fact is highlighted in Fig. 2 (b), where a DBS that is initially at ${\mathbf{x}}$ could possibly fall into ${\mathcal{R}}_{3}$ after a displacement of $vt$ with a non-linear i.i.d. mobility model. Hence, the impact of ${\mathcal{X}}$ is not zero in this case, which implies that the density of interfering DBSs is less than $\lambda_{0}$ in ${\mathcal{R}}_{3}$.

It is in fact possible to make a more formal statement about the comparison of SL mobility model with the other models, which is done next.

Theorem 1**.**

The expected number of interferers in the disc ${\mathcal{B}}=b({\mathbf{o}}^{\prime},u_{0}+vt)$ at any time $t$ is maximized over the space of i.i.d. mobility models (including the ones where drones follow curved trajectories) when interferers follow the SL mobility model.

Proof:

See Appendix -B. ∎

Remark 5**.**

Theorem 1 demonstrates that the average number of interferers in any neighborhood of the typical UE is higher in the SL mobility model compared to the other i.i.d. mobility models. Consequently, the average received rate at the typical UE under the SL mobility model is lower compared to the other i.i.d. mobility models.

In the RS mobility model, the net displacement of DBSs until time $t$ is $L(t)=\min\{vt,R\}$, where $R$ is a random variable that determines the displaced distances of DBSs. The following corollary gives the network density for the RS mobility model.

Corollary 2**.**

(Point Process of Interferers for the RS Model).* When interfering DBSs move based on the RS mobility model, the network of interfering DBSs will be an inhomogeneous PPP with the same density as given in Lemma 2. Furthermore, the function $\beta(t,u_{\mathbf{x}},u_{0})$ can be written as*

[TABLE]

where $r=\min\{vt,u_{\mathbf{x}}+u_{0}\}$.

Proof:

Starting with $L(t)=\min\{vt,R\}$, we can write the cdf and pdf of $L(t)$ as

[TABLE]

Similar to the proof of Corollary 1, we consider two cases to derive $\beta(t,u_{\mathbf{x}},u_{0})$ as follows.

• •

$0\leq u_{\mathbf{x}}\leq|u_{0}-vt|$.

[TABLE]

• •

$|u_{0}-vt|\leq u_{\mathbf{x}}\leq u_{0}+vt$.

[TABLE]

Combining these equations into a single one, we end up with (8) and the proof is complete. ∎

Note that the SL mobility model is a special case of the RS mobility model when $R\to\infty$, which means that the DBSs never stop. Mathematically speaking, evaluating $\beta(t,u_{\mathbf{x}},u_{0})$ in Corollary 2 with $F_{R}(l)=0$ and $f_{R}(l)=\delta(l-\infty)=0$ for $l<\infty$, we end up with (6).

IV RW and RWP Mobility Models

In the previous section, we characterized the point process of interferers for both service models, which required the distribution of the net displacement of each DBS as a function of time. As it was shown in Corollaries 1 and 2, this distribution can be easily derived for the SL and RS mobility models. However, characterization of this distribution for the RW and RWP mobility models is not straightforward and is the main focus of this section. Note that the distributional results provided here are novel and may be useful in their own right.

As stated in Definition 1, we assume that the direction of the $i$-th movement of a DBS is $\Theta_{i}\sim U[0,2\pi)$, which is selected independently of the other DBSs, and the flights $R_{i}$ are i.i.d. with cdf and pdf of $F_{R}(.)$ and $f_{R}(.)$, respectively. Fig. 3 shows an example of the RW mobility model. Let ${\mathbf{x}}(0)$ and ${\mathbf{x}}(n-1)$ be the initial location of a DBS and its location after $(n-1)$ flights, respectively. Define $Z_{n-1}$ and $\Psi_{n-1}$ as the net displacement of a DBS between ${\mathbf{x}}(0)$ and ${\mathbf{x}}(n-1)$ and the angle between the $x$-axis and the line connecting ${\mathbf{x}}(0)$ and ${\mathbf{x}}(n-1)$, respectively. Furthermore, we define $S_{n-1}$ and $L(t)$ as the total distance traveled from ${\mathbf{x}}(0)$ to ${\mathbf{x}}(n-1)$ and the net displacement of a DBS until time $t$, respectively. In order to properly find the distribution of $L(t)$, we need the distributional characteristics of $S_{n-1}$, $Z_{n-1}$, and $\Psi_{n-1}$, which is done next.

From Fig. 3, we can derive the equations for $S_{n-1}$, $Z_{n-1}$, and $\Psi_{n-1}$ as

[TABLE]

In the next lemma, we derive the distribution of $\Psi_{n}$.

Lemma 3**.**

When interfering DBSs move based on the RW mobility model, then irrespective of the flight distribution $R$, the random variable $\Psi_{n}$ is distributed uniformly in $[0,2\pi)$.

Proof:

See Appendix -C. ∎

Using the same methodology as in the proof of Lemma 3, one can also find the distribution of $Z_{n}$. However, it turns out that the $(n-1)$-fold integral that arises in the derivation of the pdf of $Z_{n}$ does not have a closed-form solution in general. Nevertheless, when the flights are i.i.d. Rayleigh random variables, we can derive the distribution of $Z_{n}$ in closed form.

Lemma 4**.**

When interfering DBSs move based on the RW mobility model and the flights are Rayleigh distributed with parameter $\sigma$, then $Z_{n}$ is also Rayleigh distributed with parameter $\sigma\sqrt{n}$.

Proof:

See Appendix -D. ∎

Lemma 4 establishes a useful result for the distribution of $Z_{n}$ when the flights are Rayleigh distributed. On the other hand, if the flights have a general non-Rayleigh distribution, we can find an asymptotic distribution for $Z_{n}$ as $n\to\infty$. We present this result in the next lemma.

Lemma 5**.**

When interfering DBSs move based on the RW mobility model and the flights have general non-Rayleigh but i.i.d. distributions with a mean and variance of $\mu_{R}$ and $\sigma_{R}^{2}$, respectively, then $\frac{Z_{n}}{\sqrt{n}}$ will have a Rayleigh distribution with parameter $\sqrt{\frac{\mu_{R}^{2}+\sigma_{R}^{2}}{2}}$ as $n\to\infty$.

Proof:

Define $X=\sum_{i=1}^{n}\frac{1}{\sqrt{n}}R_{i}\cos(\Theta_{i})$ and $Y=\sum_{i=1}^{n}\frac{1}{\sqrt{n}}R_{i}\sin(\Theta_{i})$. Since $\Theta_{i}$’s are i.i.d. with uniform distribution in $[0,2\pi)$, the central limit theorem (CLT) asserts that as $n\to\infty$, $X$ and $Y$ will have Gaussian distributions. Since $R_{i}$ and $\Theta_{i}$ are independent of each other, the moments of $X$ can be computed as follows: ${\mathbb{E}}[X]=\sum_{i=1}^{n}\frac{1}{\sqrt{n}}{\mathbb{E}}[R_{i}]{\mathbb{E}}[\cos(\Theta_{i})]=0$ and ${\mathbb{E}}[X^{2}]={\mathbb{E}}\left[\sum_{i=1}^{n}\sum_{j=1}^{n}\frac{1}{n}R_{i}R_{j}\cos(\Theta_{i})\cos(\Theta_{j})\right]={\mathbb{E}}\left[\sum_{i=1}^{n}\frac{1}{n}R_{i}^{2}\cos^{2}(\Theta_{i})\right]=\frac{\mu_{R}^{2}+\sigma_{R}^{2}}{2}$. Note that the same is also true for $Y$. Hence, $X\sim{\mathcal{N}}(0,\frac{\mu_{R}^{2}+\sigma_{R}^{2}}{2})$ and $Y\sim{\mathcal{N}}(0,\frac{\mu_{R}^{2}+\sigma_{R}^{2}}{2})$. Now since ${\mathbb{E}}[XY]=0$, $X$ and $Y$ are uncorrelated, and thus, independent. Therefore, $\frac{Z_{n}}{\sqrt{n}}=\sqrt{X^{2}+Y^{2}}$ is Rayleigh distributed with parameter $\sqrt{\frac{\mu_{R}^{2}+\sigma_{R}^{2}}{2}}$ as $n\to\infty$. This completes the proof. ∎

After characterizing $Z_{n}$ and $\Psi_{n}$, we need to find the distribution of $S_{n}$ as well. Recall that $S_{n}$ is the sum of $n$ i.i.d. random variables and for a general distribution of $R_{i}$, the distribution of $S_{n}$ is not known. There have been some works that establish the distribution of $S_{n}$ for different distributions of $R_{i}$ (see [36] for a comprehensive review). When $R_{i}$’s are exponentially distributed, the distribution of $S_{n}$ is known to be Erlang. However, the exact distribution of $Z_{n}$ is not known for exponentially distributed flights. When $R_{i}$’s are Rayleigh, the distribution of $S_{n}$ has been investigated in the literature [37] and the result is given here without proof.

Lemma 6**.**

The pdf of the sum of $n$ i.i.d. Rayleigh random variables with parameter $\sigma$ can be approximated as

[TABLE]

*where $b=\sigma^{2}\frac{\sqrt[n]{(2n-1)!!}}{n}$ and the constants $a_{0}$, $a_{1}$, and $a_{2}$ are derived numerically using a nonlinear curve-fitting least square method based on the trust region reflective algorithm [38]. *

In order to better understand the distributional properties of the RW mobility model, we also need to find the joint distribution of $S_{n}$ and $Z_{n}$ for a given number of flights $n$. Given $S_{n}$, the range of the values for $Z_{n}$ is upper bounded by the value of $S_{n}$, which suggests a dependency between $S_{n}$ and $Z_{n}$. In the next proposition, we characterize this joint distribution.

Proposition 1**.**

The joint pdf of $S_{n}$ and $Z_{n}$ for a given number of flights $n\geq 2$ can be written as a $2(n-1)$-fold integral given as

[TABLE]

where $J_{n}=\sum_{i=1}^{n-1}x_{i}$, $K_{n}=\sqrt{\sum_{i,j=1}^{n-1}x_{i}x_{j}\cos(\xi_{i}-\xi_{j})}$, and ${\mathcal{R}}$ is defined as the region where $0\leq x_{i}<\infty$, $0\leq\xi_{i}<2\pi$ for $1\leq i\leq n-1$, and $|s-J_{n}-K_{n}|\leq z\leq s-J_{n}+K_{n}$.

Proof:

See Appendix -E. ∎

Note that for $n=1$, the random variables $S_{n}$ and $Z_{n}$ will become identical and equal to $R_{1}$, and thus, their joint distribution will be the same as the distribution of $R_{1}$, i.e., $f_{R}(.)$. For $n=2$, the result of Proposition 1 can be further simplified, which is given in the next corollary.

Corollary 3**.**

The joint pdf between $S_{2}$ and $Z_{2}$ can be written as

[TABLE]

when $s>z$ and zero otherwise.

Remark 6**.**

As $n$ gets larger, the dependency between $S_{n}$ and $Z_{n}$ will become less significant, and thus, we can ignore it for large enough $n$. Hence, we approximate the joint pdf of $S_{n}$ and $Z_{n}$ by assuming them to be independent for $n\geq 3$, thus giving $f_{S_{n},Z_{n}}(s,z)\approx f_{S_{n}}(s)f_{Z_{n}}(z)$. Note that $n\geq 3$ is large enough for this result to be reasonably accurate.

Using the results that we have derived for the RW mobility model so far, we can now compute the distribution of $L(t)$, i.e., the distance between the location of a DBS at time $t$ and its original location at time $t=0$. Fig. 3 shows a DBS that is flying in its $n$-th flight. Note that the distance traveled by the DBS until time $t$ is $vt$, and thus, the residual distance in the $n$-th flight will be $vt-S_{n-1}$. We now state the main result in the next proposition.

Proposition 2**.**

When interfering DBSs move based on the RW mobility model, the cdf and pdf of $L(t)$ are given as

[TABLE]

when $l\leq vt$, respectively. Otherwise, we have $F_{L}(l;t)=1$ and $f_{L}(l;t)=0$.

Proof:

See Appendix -F. ∎

For $n=2$, the double integrals in (13) and (14) can be written as single integrals since $S_{1}=Z_{1}=R_{1}$. Mathematically speaking, when $l<vt$ we have

[TABLE]

This result will be useful for our further approximations on the distribution of $L(t)$. Applying the results of Proposition 2 to Lemma 2, we can compute $\beta(t,u_{\mathbf{x}},u_{0})$ in (5) and get the density of the network of interfering DBSs for the RW mobility model.

Building on the RW mobility model, we now derive the distribution of $L(t)$ for the RWP mobility model as well. The RWP mobility model is defined by a sequence of quadruples at the $i$-th flight period: two waypoints ${\mathbf{p}}_{i-1}$ and ${\mathbf{p}}_{i}$ as the starting and destination waypoints, respectively, one transition length $R_{i}$, and one waiting time $T_{i}$ at the destination waypoint. In an infinite network, we can observe that the only difference between the RWP and the RW mobility models is the inclusion of the random variable $T_{i}$ as the waiting time at the end of each flight. Assuming that the random variables $\{T_{i}\}$ are i.i.d. with cdf and pdf of $F_{T}(.)$ and $f_{T}(.)$, respectively, and also independent from $R_{i}$, the following proposition extends the results of Proposition 2 to the RWP mobility model. Before stating the main result, we define $W_{n}$ as the aggregate waiting time until the end of the $n$-th flight, i.e., $W_{n}=\sum_{i=0}^{n}T_{i}$. Note that $W_{n}$ will be independent from $S_{n}$ and $Z_{n}$. Moreover, we assume that $W_{0}=T_{0}\neq 0$, which implies that there is an initial random waiting time before the DBSs start to move.

Proposition 3**.**

When interfering DBSs move based on the RWP mobility model, the cdf and pdf of $L(t)$ are given as

[TABLE]

when $l\leq vt$, respectively. Otherwise, we have $F_{L}(l;t)=1$ and $f_{L}(l;t)=0$.

Proof:

See Appendix -G. ∎

As in the RW mobility model, the results can be further simplified for $n=2$ and $l<vt$ as

[TABLE]

Similar to the RW scenario, we can apply the distribution of $L(t)$ derived in Proposition 3 to Lemma 2 and compute the density of the network of interferers for the RWP mobility model.

V Average and Session Rates

In this section, we compute the ${\rm SIR}$-based metrics defined in Section II-D for the typical UE in the network. Equipped with the density of the network of interfering DBSs for each mobility model under both the UIM and the UDM, we first derive the average received rate by the typical UE in the UDM. The following theorem provides this result.

Theorem 2**.**

In the UDM and using any i.i.d. mobility model described in the previous sections, the average received rate by the typical UE at time $t$ is given as

[TABLE]

where \mathcal{L}_{I(t)}(s\bigr{\rvert}\mathbf{x}_{0}(t)) is the conditional Laplace transform of interference given as

[TABLE]

Proof:

We start by writing the complementary cumulative distribution function (ccdf) of ${\rm SIR(t)}$ conditioned on $\mathbf{x}_{0}(t)$ as

[TABLE]

where in $(a)$ the Nakagami-$m$ fading assumption is used and in $(b)$ we used the definition of the incomplete gamma function for integer values of $m_{0}$. The conditional Laplace transform of interference at time $t$ can be computed as

[TABLE]

where (a) results from the moment generating function (MGF) of the gamma distribution and (b) follows from the probability generating functional (PGFL) of a general PPP. Now, the average rate at time $t$ can be written as $R(t)=\mathbb{E}\left[\log\left(1+{\rm SIR}(t)\right)\right]=$

[TABLE]

where $f_{\Gamma}(\gamma;t)$ is the pdf of ${\rm SIR}(t)$ and in the last equation, we used integration by parts and deconditioned the result on $u_{0}(t)$. Note that in the UDM, we have $u_{0}(t)=[u_{0}-vt]^{+}$, where $[x]^{+}=x$ if $x\geq 0$ and $[x]^{+}=0$ otherwise. This completes the proof. ∎

For the UIM, since the density of the network of interfering DBSs is given in (3), the received rate by the typical UE will be given as in (21) evaluated at $t=0$. Finally, the session rate at time $T$ is derived by integrating over the average rate with respect to $t$ as in (2).

VI Simulation Results

In this section, we provide numerical simulations to verify our analytical results and provide several useful insights about the system-level performance of this network. We assume that the network density is $\lambda_{0}=10^{-6}$ and DBSs move with a constant speed $v=45~{}{\rm km/h}$ using one of the four mobility models defined in Section II. We consider the low altitude platform (LAP) for the flight of DBSs and assume $h\in\{100,\,200\}~{}{\rm m}$ as typical values for the height in LAP. Furthermore, we assume the path loss exponent is $\alpha=3$. For the RS, RW, and RWP mobility models, we assume that the flight distances are distributed as i.i.d. Rayleigh random variables with mean $500~{}{\rm m}$. Moreover, in the RWP mobility model, the hovering times follow i.i.d. exponential random variables with mean $5~{}{\rm s}$. Note that the evaluations of the multiple integrals that occur in the analyses of the RW and RWP mobility models are carried out using the Monte Carlo (MC) integration method.

VI-A Distribution of $L(t)$ in RW and RWP

Characterizing the distribution of $L(t)$ for a given time $t$ is essential for the analysis of the RW and RWP mobility models. We derived the exact cdf and pdf of $L(t)$ for the RW and RWP mobility models in Propositions 2 and 3, respectively. Fig. 4 shows the pdf of $L(t)$ for both the RW and RWP mobility models at $t\in\{50,100,300\}~{}{\rm s}$. Note that as $t\to\infty$, the distance traveled during the $n$-th flight will be much smaller than $Z_{n-1}$, and thus, we have $L(t)\approx Z_{n-1}$. Now, according to Lemma 5, as $t\to\infty$, the distribution of $Z_{n-1}$ converges to a Rayleigh distribution. Hence, $L(t)$ is also Rayleigh distributed as $t\to\infty$. This trend can also be noticed in Fig. 4.

VI-B Point Process of Interferers

In Figs. 5 and 6, we plot the density of the network of interfering DBSs for all the mobility models considered in this paper at $t\in\{20,40,50,200\}~{}{\rm s}$. We assume that the serving DBS follows the UDM and the exclusion zone radius is $u_{0}=500~{}{\rm m}$. In the SL mobility model, as also highlighted in Remark 4, the density will be divided into two homogeneous parts and one bowl-shaped inhomogeneous part after $t=\frac{u_{0}}{v}$. Furthermore, the inhomogeneous part will become homogeneous as $t\to\infty$, which ultimately makes the point process of interferers homogeneous. This fact can also be directly inferred from Corollary 1 by taking the limit of (6) as $t\to\infty$ and $u_{\mathbf{x}}\to vt$. According to Fig. 6, this homogenization happens for the RW and RWP mobility models as well. However, this is not the case in the RS mobility model, since the DBSs “stop” moving after a period of time. Hence, the point process of interferers does not evolve with time, making it inhomogeneous for all time $t$. In Fig. 6, based on Remark 6, we also plot the density of the network of interfering DBSs for the RW and RWP mobility models by assuming that $S_{n}$ and $Z_{n}$ are independent for $n\geq 3$. Clearly, this approximation is fairly accurate.

VI-C Impact of Fading and Height

In order to show the effect of the Nakagami-$m$ fading parameter on the performance of the network, we plot the average rate in the UDM under the SL mobility model for $m=m_{0}\in\{1,2\}$ in Fig. 8. Since increasing $m$ and $m_{0}$ decreases the severity of the fading, the average rate will increase as well. Fig. 8 depicts the impact of DBS heights in the performance of the network for both service models under the SL mobility model. As is clear in this figure, the average rate increases as height decreases, which can also be observed directly from (21) and (22).

VI-D Impact of Mobility Models

Under the UDM, we compare different mobility models studied in this paper in terms of the average and session rates in Figs. 10 and 10, respectively. As can be seen in these figures, the SL mobility model acts as a lower bound on the performance to these models. This result is essentially the one that we proved in Theorem 1 by demonstrating that the average number of interferers in the vicinity of the typical UE is higher in the SL mobility model than the other i.i.d. mobility models. Hence, the performance of the network under the SL mobility model will be worse than the other mobility models.

VII Conclusion

In this paper, we presented an in-depth and unified analysis of a mobile drone cellular network operating at a constant height to serve the UEs on the ground. Specifically, we considered four mobility models for the DBSs, i.e., (i) SL, (ii) RS, (iii) RW, and (iv) RWP, and provided several fundamental distributional properties of these models. The use of the SL mobility model for drone networks was inspired by the simulation models used in the 3GPP studies of drone networks, while the others are useful canonical models (or their variants) that have been used extensively in wireless networks and provide a reasonable balance between realism and tractability. The serving DBS is selected based on the nearest-neighbor association policy and moves according to two service models: (i) based on the same mobility model as the interfering DBSs (UIM), and (ii) towards the typical UE at a constant height and keeps hovering above the location of the typical UE (UDM). We proposed a novel characterization of the point process of DBSs for both the UIM and the UDM, using which we analyzed the average received rate and the session rate at the typical UE. Borrowing tools from the calculus of variations, we mathematically showed that the SL mobility model acts as a lower bound on the system-level performance of our mobile drone network over the space of all i.i.d. mobility models. To the best of our understanding, this is the first work that offers a unified analysis of canonical mobility models for a drone cellular network in an infinite plane and establishes meaningful connections between them. While this work offers many useful insights in the canonical settings, it will be useful to extend some of these results to more realistic mobility models, including the ones developed from actual mobility traces of drones (as and when they become available).

-A Proof of Lemma 2

Since we have started with an inhomogeneous PPP with the density given as (3) for $t=0$ in the UDM and the displacements are independent of each other, the network of interfering DBSs will also be an inhomogeneous PPP at every time $t$ due to displacement theorem [35]. According to Lemma 1, if there was no exclusion zone and DBSs moved independently of each other, as in our mobility models, the network of all DBSs (including the serving DBS) would have remained a homogeneous PPP with density $\lambda_{0}$. Taking ${\mathcal{X}}$ into account, the resulting density of the network can be partitioned into two sets: (i) density contributed by ${\mathcal{X}}$ (denoted as $\lambda_{1}(t;u_{\mathbf{x}},u_{0})$), i.e., due to the points that are initially inside ${\mathcal{X}}$, and (ii) density of interferers (denoted as $\lambda(t;u_{\mathbf{x}},u_{0})$), i.e., due to the points that are initially outside ${\mathcal{X}}$. Since the resulting network density is $\lambda_{0}$, we have $\lambda(t;u_{\mathbf{x}},u_{0})=\lambda_{0}-\lambda_{1}(t;u_{\mathbf{x}},u_{0})$.

Define $N(t)$ as the average number of points that are initially inside ${\mathcal{X}}$ and after a displacement of $L(t)$ land on an infinitesimal annulus ${\mathcal{A}}$ with an inner and outer radii of $u_{\mathbf{x}}$ and $u_{\mathbf{x}}+{\rm d}u_{\mathbf{x}}$, respectively. Since the density in this case is rotation invariant, it is sufficient to consider an annulus centered at ${\mathbf{o}}^{\prime}$ for our analysis. By definition, we can write

[TABLE]

where the denominator denotes the area of ${\mathcal{A}}$. As shown in Fig. 11, let ${\mathbf{x}}$, ${\mathbf{y}}_{1}$, and ${\mathbf{y}}_{2}$ denote a randomly selected point inside ${\mathcal{X}}$ with distance $r$ from ${\bf o^{\prime}}$ and two intersection points of $b({\mathbf{x}},L(t))$ with ${\mathcal{A}}$, respectively. Writing the cosine law in triangles ${{\bf o^{\prime}}{\mathbf{x}}{\mathbf{y}}_{1}}$ and ${{\bf o^{\prime}}{\mathbf{x}}{\mathbf{y}}_{2}}$, we have

[TABLE]

where $\theta=\angle{\bf o^{\prime}}{\mathbf{x}}{\mathbf{y}}_{1}$ and $\varphi=\angle{\mathbf{y}}_{1}{\mathbf{x}}{\mathbf{y}}_{2}$. Eliminating $\theta$ in both of these equations gives

[TABLE]

where $a=\frac{L(t)^{2}+r^{2}-u_{\mathbf{x}}^{2}}{2u_{\mathbf{x}}},b=\frac{2rL(t)}{2u_{\mathbf{x}}},c=\sqrt{b^{2}-a^{2}}$. Note that the probability that a DBS at ${\mathbf{x}}$ lands on ${\mathcal{A}}$ after a displacement of $L(t)$ is $\frac{2\varphi}{2\pi}$. Hence, we write $N(t)$ by considering all such points ${\mathbf{x}}$ in ${\mathcal{X}}$ as $N(t)={\mathbb{E}}\left[\int_{0}^{u_{0}}\frac{\varphi}{\pi}2\pi r\lambda_{0}\,{\rm d}r\right]$, where the expectation is taken over $L(t)$. This gives

[TABLE]

where we have changed the order of limit with the expectation and the integration due to the continuity of the integrand. Representing ${\rm d}u_{\mathbf{x}}$ with $x$ for simplicity, we compute the limit as

[TABLE]

where in $(*)$ we used the Taylor series expansion $\cos^{-1}(1-x)=\sqrt{2x}+\Theta(x^{3/2})$ as $x\to 0$, where $f(x)=\Theta(g(x))$ implies that $f(x)$ is asymptotically bounded by $g(x)$ both from above and below. Note that since the triangle inequality holds for the triple $(u_{\mathbf{x}},r,L(t))$, the result in (-A) is real and positive, as expected. Plugging (-A) into (24), we have

[TABLE]

where ${\mathcal{I}}_{1}=\left\{|l-u_{\mathbf{x}}|\leq r\leq l+u_{\mathbf{x}}\right\}\bigcap\left\{0\leq r\leq u_{0}\right\}$. Since the net displacement of a DBS at time $t$ cannot exceed its aggregate traveled distance $vt$, we have $L(t)\leq vt$. Hence,

[TABLE]

where ${\mathcal{I}}_{2}=\left\{|vt-u_{\mathbf{x}}|\leq r\leq vt+u_{\mathbf{x}}\right\}\bigcap\left\{0\leq r\leq u_{0}\right\}$. Simplifying the last step requires tedious integrations and the details are skipped to maintain brevity. Finally, the density of the network of interfering DBSs is summarized as (4) and (5) in the lemma statement. ∎

-B Proof of Theorem 1

The intensity measure or the expected number of interfering DBSs in the Borel set ${\mathcal{B}}=b({\mathbf{o}}^{\prime},u_{0}+vt)$ is given as $\Lambda({\mathcal{B}})=2\pi\int_{0}^{u_{0}+vt}u_{\mathbf{x}}\lambda(t;u_{\mathbf{x}},u_{0})\,{\rm d}u_{\mathbf{x}}$. Defining $\Lambda_{1}({\mathcal{B}})$ and $\Lambda_{2}({\mathcal{B}})$ as the intensity measures for the SL mobility model and an i.i.d. mobility model, respectively, we need to show that $\Lambda_{1}({\mathcal{B}})\geq\Lambda_{2}({\mathcal{B}})$. Note that we only prove this for $t\leq\frac{u_{0}}{v}$ and the proof for $t>\frac{u_{0}}{v}$ follows on the similar lines. Using (4) and (5), we write $\Lambda_{1}({\mathcal{B}})-\Lambda_{2}({\mathcal{B}})$ as

[TABLE]

where $g(l,u_{\mathbf{x}})=\frac{1}{\pi}\cos^{-1}\left(\frac{l^{2}+u_{\mathbf{x}}^{2}-u_{0}^{2}}{2lu_{\mathbf{x}}}\right)$ and we denote the integrand by $\mathfrak{L}(u_{\mathbf{x}},F_{L}(u_{\mathbf{x}}),f_{L}(u_{\mathbf{x}}))$, which is a functional with $u_{\mathbf{x}}$ being the independent variable. Note that we also have $F_{L}(0)=0$ and $F_{L}(vt)=1$ by definition. Using the Euler-Lagrange equation, we show that $\Lambda_{1}({\mathcal{B}})-\Lambda_{2}({\mathcal{B}})$ attains its minimum at zero. From the calculus of variations, we know that $F_{L}$ is a critical (extremum) function for the functional $\mathfrak{L}$, if it satisfies the Euler-Lagrange equation: $\frac{\partial\mathfrak{L}}{\partial F_{L}}-\frac{{\rm d}}{{\rm d}u_{\mathbf{x}}}\left(\frac{\partial\mathfrak{L}}{\partial f_{L}}\right)=0$. We compute the first term as $\frac{\partial\mathfrak{L}}{\partial F_{L}}=u_{\mathbf{x}}\mathbf{1}(u_{0}-u_{\mathbf{x}})$ and the second term as

[TABLE]

where in $(a)$ we used the Leibniz integral rule along with $g(|u_{\mathbf{x}}-u_{0}|,u_{\mathbf{x}})=\mathbf{1}(u_{0}-u_{\mathbf{x}})$. Applying these derivatives to the Euler-Lagrange equation, we get $\int_{|u_{\mathbf{x}}-u_{0}|}^{vt}\frac{\partial}{\partial u_{\mathbf{x}}}\left(u_{\mathbf{x}}g(l,u_{\mathbf{x}})\right)\,{\rm d}l=0$, which does not include either $F_{L}$ or $f_{L}$. This means that the original function in (27) attains its extremum at the boundaries. We now evaluate $\mathfrak{L}(u_{\mathbf{x}},F_{L}(u_{\mathbf{x}}),f_{L}(u_{\mathbf{x}}))$ at $u_{\mathbf{x}}=u_{0}\pm vt$ as

[TABLE]

Since both of the boundary values are [math], we conclude that the minimum value in (27) is [math]. ∎

-C Proof of Lemma 3

In order to derive the distribution of $\Psi_{n}$, we introduce $(n-1)$ auxiliary random variables $\Xi_{i}=\Theta_{i}$, $1\leq i\leq n-1$, and find the joint pdf of $n$ random variables $\Psi_{n}$ and $\Xi_{i}$, $1\leq i\leq n-1$. We then integrate out these auxiliary random variables to find the pdf of $\Psi_{n}$. We start by solving the system of $n$ equations (one equation in (9) for the definition of $\Psi_{n}$ and $(n-1)$ equations introduced by the auxiliary random variables) to derive $\Theta_{i}$’s in terms of $\Psi_{n}$ and $\Xi_{i}$’s. The result can be written as two sets of solutions as follows:

[TABLE]

where $\Delta=\sum_{i=1}^{n-1}R_{i}\sin(\Xi_{i}-\Psi_{n})$ and $\Delta^{\prime}=\sqrt{R_{n}^{2}-\Delta^{2}}$. Note that $2\Psi_{n}=\Xi_{n_{1}}+\Xi_{n_{2}}$. Computing the determinant of the Jacobian matrix, we get $|{\mathcal{J}}|=|\frac{\partial\Theta_{n}}{\partial\Psi_{n}}|$ for both solution sets. Hence,

[TABLE]

By some algebraic manipulations, we find that |{\mathcal{J}}|\bigr{\rvert}_{\Theta_{n}=\Xi_{n_{1}}}+|{\mathcal{J}}|\bigr{\rvert}_{\Theta_{n}=\Xi_{n_{2}}}=2. Now, since $\Theta_{i}\sim U[0,2\pi)$, the joint distribution of $\Psi_{n}$ and $\Xi_{i}$’s can be written as

[TABLE]

where the boldface letters represent vector random variables. Integrating out $\Xi_{i}$, $1\leq i\leq n-1$, the distribution of $\Psi_{n}$ is derived as $f_{\Psi_{n}}(\psi_{n})=\frac{1}{\pi}$ for $\psi_{n}\in[-\frac{\pi}{2},\frac{\pi}{2})$, due to the range of the $\tan^{-1}$ function. Finally, since the true range of $\Psi_{n}$ is $[-\pi,\pi)$, we conclude that $f_{\Psi_{n}}(\psi_{n})=\frac{1}{2\pi}$ for $\psi_{n}\in[-\pi,\pi)$, and the proof is complete. ∎

-D Proof of Lemma 4

We consider an RW of two flights $F_{1}$ and $F_{2}$, where $F_{1}$ and $F_{2}$ are two independent Rayleigh random variables with parameters $\sigma_{1}$ and $\sigma_{2}$, respectively. Let $F$ be the distance from the start of the first flight to the end of the second flight. We show that $F$ is Rayleigh distributed with parameter $\sigma=\sqrt{\sigma_{1}^{2}+\sigma_{2}^{2}}$, which proves the lemma by induction. Since $F_{1}$ and $F_{2}$ are Rayleigh distributed, they can be written as $F_{1}=\sqrt{X_{1}^{2}+Y_{1}^{2}}$ and $F_{2}=\sqrt{X_{2}^{2}+Y_{2}^{2}}$, where $X_{1},Y_{1}\sim{\mathcal{N}}(0,\sigma_{1}^{2})$ are $X_{2},Y_{2}\sim{\mathcal{N}}(0,\sigma_{2}^{2})$ are independent Gaussian random variables. Note that due to the independence of $X_{1}$ and $X_{2}$, their sum will also be a Gaussian random variable. Therefore, $X=X_{1}+X_{2}\sim{\mathcal{N}}(0,\sigma_{1}^{2}+\sigma_{2}^{2})$. Likewise, $Y=Y_{1}+Y_{2}\sim{\mathcal{N}}(0,\sigma_{1}^{2}+\sigma_{2}^{2})$. We can now represent $F_{1}$ and $F_{2}$ in the abscissa and ordinate axes as $[X_{1},Y_{1}]$ and $[X_{2},Y_{2}]$, respectively. Hence, $F=\sqrt{X^{2}+Y^{2}}$ will have a Rayleigh distribution with parameter $\sigma=\sqrt{\sigma_{1}^{2}+\sigma_{2}^{2}}$. This completes the proof. ∎

-E Proof of Proposition 1

Similar to the proof of Lemma 3, we define $2(n-1)$ auxiliary random variables $X_{i}=R_{i}$ and $\Xi_{i}=\Theta_{i}$, $1\leq i\leq n-1$ and find the joint pdf of these $2n$ random variables. Solving this system of $2n$ equations, we can write $R_{i}$ and $\Theta_{i}$, $1\leq i\leq n$ in terms of $S_{n}$, $Z_{n}$, $X_{i}$, and $\Xi_{i}$, $1\leq i\leq n-1$. Note that $(2n-1)$ of these equations have trivial solutions, i.e., $R_{i}=X_{i}$ and $\Theta_{i}=\Xi_{i}$ for $1\leq i\leq n-1$ and $R_{n}\!=\!S_{n}\!-\!\sum_{i=1}^{n-1}X_{i}$. For $\Theta_{n}$, we write the governing equation as

[TABLE]

Note that the exact expression for $\Theta_{n}$ will not be required for the derivation of our joint pdf of interest. We then derive the determinant of the Jacobian matrix as

[TABLE]

Solving (28) and (29), we can eliminate $R_{n}$ and $\Theta_{n}$ to get $|{\mathcal{J}}|$ only in terms of $S_{n}$, $Z_{n}$, $X_{i}$, and $\Xi_{i}$, $1\leq i\leq n-1$. Hence, we can write the joint distribution of these $2n$ random variables as

[TABLE]

where $J_{n}=\sum_{i=1}^{n-1}x_{i}$, $K_{n}=\sqrt{\sum_{i,j=1}^{n-1}x_{i}x_{j}\cos(\xi_{i}-\xi_{j})}$. Integrating this pdf $2(n-1)$ times with respect to $x_{i}$ and $\xi_{i}$, $1\leq i\leq n-1$, we end up with (11) and the proof is complete. ∎

-F Proof of Proposition 2

Considering the RW mobility model at time $t$, a DBS is either in its first flight, or in its second flight, etc. Hence, initializing $S_{0}=0$, the event $\Omega=\bigcup_{n=1}^{\infty}\left(S_{n-1}\leq vt<S_{n}\right)$ has unit probability. Thus, we can write the cdf of $L(t)$ as

[TABLE]

where in $(a)$ we used the cosine law and in $(b)$ we conditioned the probability on knowing $S_{n-1}$ and $Z_{n-1}$. Note that $\Phi_{n}=\pi-\Theta_{n}+\Psi_{n-1}$ is the angle between the direction of $Z_{n-1}$ and the direction of the $n$-th flight. We now rewrite the probability in (-F) as ${\mathbb{P}}\left[\cos(\Phi_{n})\geq x\right]$, where $x=\frac{z^{2}+(vt-s)^{2}-l^{2}}{2z(vt-s)}$. In order for this probability to be non-zero, we have two cases: (i) $x\leq-1$ which gives $z+(vt-s)\leq l$, and (ii) $-1<x\leq 1$ which gives $|z-(vt-s)|\leq l<z+(vt-s)$, which is the triangle inequality. Note also that it is clear from our setup that $z\leq s\leq vt$. When $l\geq vt$, only the first case will hold and we can write the cdf as

[TABLE]

where in the last equality we used the convolution integral that arises in the derivation of the cdf of $S_{n}=S_{n-1}+R_{n}$. Hence, $F_{L}(l;t)=1$ for $l\geq vt$, as expected. On the other hand, both cases can occur when $l<vt$. For the first case, we have $0\leq z\leq l-(vt-s)$ and $vt-l\leq s\leq vt$. Similarly, for the second case we have $|l-(vt-s)|\leq z\leq\min\{s,\,l+(vt-s)\}$ and $\frac{vt-l}{2}\leq s\leq vt$ based on the triangle inequality. Now, according to Lemma 3, $\Psi_{n-1}$ is uniformly distributed in $[0,2\pi)$, and since $\Theta_{n}\sim[0,2\pi)$ is independent of $\Psi_{n-1}$, the random variable $\Phi_{n}$ will have a symmetric triangular distribution. However, since the range of values of $\Phi_{n}$ is between [math] and $2\pi$, we have $\Phi_{n}\sim[0,2\pi)$. Hence, the cdf of $L(t)$ for $l<vt$ can be written as in (13). Note that when $n=1$, we have $L(t)=vt$ and the cdf becomes ${\mathbb{P}}[vt\leq l,vt<R_{1}]=(1-F_{R}(vt)){\bf 1}(l-vt)$. Differentiating the derived cdf with respect to $l$ using the Leibniz integral rule for two-dimensional integrals, we end up with the pdf of $L(t)$ as written in (14). ∎

-G Proof of Proposition 3

Define $M_{n}=S_{n}+vW_{n-1}$ and $Y_{n}=S_{n}+vW_{n}$. Considering the RWP mobility model at time $t$, a DBS is either in one of its flight states, or in one of its waiting time states. Hence, initializing $S_{0}=0$ and $W_{-1}=0$, the event $\Omega^{\prime}=\bigcup_{n=1}^{\infty}\left(V_{n-1}\cup F_{n}\right)$ has unit probability, where $V_{n}=\{M_{n}\leq vt<Y_{n}\}$ and $F_{n}=\{Y_{n-1}\leq vt<M_{n}\}$ are the $n$-th waiting and flight periods, respectively. Thus, we can write the cdf of $L(t)$ as follows.

[TABLE]

where in $(a)$ we used the fact that $L(t)=Z_{n-1}$ when a DBS is in its $(n-1)$-th waiting state and $L(t)$ follows the cosine rule when a DBS is in its $n$-th flight state. In the first summation of $(b)$, we first conditioned the probability in the integrand on knowing $T_{n-1}$ and then used the identity $W_{n-1}=W_{n-2}+T_{n-1}$ to simplify the result. In the second summation of $(b)$, we used the definition of $Y_{n-1}$ to write the joint pdf of $Y_{n-1}$ and $Z_{n-1}$ in terms of the joint pdf of $S_{n-1}$ and $Z_{n-1}$. Now, with the same reasoning as in the proof of Proposition 2, we get the cdf and pdf of $L(t)$ when $l\leq vt$ as in (17) and (18), respectively. Note that for $n=1$, we have

[TABLE]

When $l>vt$, we can write the cdf as

[TABLE]

where in $(a)$ we used the definition of $Y_{n-1}$ to get $F_{Y_{n-1}}(vt)=\int_{0}^{vt}f_{S_{n-1}}(s)F_{W_{n-1}}(t-\frac{s}{v})\,{\rm d}s$ and in $(b)$ we used the definition of $M_{n-1}$ and the identity $M_{n}=Y_{n-1}+R_{n}$ to get the result. ∎