# Stochastic geometric modelling and simulation of cellular systems for
coverage probability characterization

Hamed Nassar  Hamed Nassar Computer Science Department, Suez Canal University,
Ismailia 41522, Egypt<EMAIL_ADDRESS>,  Gehad Taher  Gehad Taher
Computer Science Department, Suez Canal University, Ismailia 41522, Egypt
<EMAIL_ADDRESS>and  El-Sayed El-Hady  El-Sayed El-Hady Basic
Sciences Department, Suez Canal University, Ismailia 41522, Egypt
<EMAIL_ADDRESS>
###### Abstract.

Stochastic geometry (SG) has been successfully used as a modelling tool for
cellular networks to characterize the coverage probability in both the
downlink (DL) and uplink (UL) systems, under the assumption that the base
stations (BS) are deployed as a Poisson point process. In the present article,
we extend this use and provide further results for interference limited and
Rayleigh fading networks, culminating in a multifaceted contribution. First,
we compactly model the two systems at once, allowing parallels to be drawn and
contrast to be created. Also, for DL we manage to obtain two closed form
expressions for two special cases. Moreover, for UL, notorious for being
difficult, we develop a clever approximation that overcomes the difficulty,
yielding excellent results. Additionally, we present two efficient Monte Carlo
simulation algorithms, designed primarily to validate the models, but can be
of great use for SG modelling of communications systems in general. Finally,
we prove two theorems at odds with popular belief in cellular communications
research. Specifically, we prove that under the SG model, the coverage
probability in both DL and UL is independent of BS density. Based on this
revelation, a plethora of results in the literature have to be re-examined to
rid them of a parameter that has been proven superfluous.

###### Key words and phrases:

Stochastic geometry; IoT; Cellular network; Downlink; Uplink; Coverage
probability

##  1\. Introduction

Stochastic geometry was initially stimulated by applications to biology,
astronomy and material sciences, then introduced seriously in the late
nineties to the field of wireless communications [1]. It is particularly
suited for modelling large scale wireless communication networks, where a
network is treated as a realization (snapshot) of a spatial point process in
the entire Euclidean plane [2]. In particular, it is a natural approach to
describe node locations in randomly formed networks, e.g. ad hoc networks, or
irregular and regular networks, e.g. cellular networks. It provides a natural
way of computing macroscopic properties, by averaging over the potential
geometrical patterns for all network nodes to obtain key performance
characteristics, such as connectivity, stability, and capacity, as functions
of a relatively small number of parameters, e.g. the intensity of the
underlying point process and the operational parameters. By a spatial average,
it is meant an average calculated over many points in the considered domain.
These locations in our context are the network elements at the time when the
snapshot is taken. Example spatial averages are the fraction of nodes which
transmit, the fraction of space which is covered or connected, the fraction of
nodes that transmit their data successfully, and the average geographic
progress achieved by a node forwarding a packet towards some destination [3].

The most common function of stochastic geometry in wireless communications is
to characterize the signal to interference and noise ratio (SINR), which can
be used to calculate many cellular performance metrics, such as outage
probability, coverage probability, spatial opportunity, spatial throughput,
network throughput, medium access probability and spectral efficiency [4], in
both downlink (DL) and uplink (UL) directions. It should be noted, however,
that in light the huge of influx of wireless emissions in recent years, the
impact of noise now pales in comparison with interference [5]. As such, there
is a growing modelling trend (see, for example, [6, 7, 8, 9]) to replace SINR
by SIR, and we will follow this trend in the present article. In all cases,
the most relevant assumption made when stochastic geometry is used in
modelling cellular networks is that base stations (BS) are deployed in the
Euclidean space as a Poisson Point Process (PPP), although many variant
processes, e.g. Thomas cluster processes or Matérn cluster processes, have
been suggested as well [10].

The stochastic geometric model of cellular networks, portrayed in Figure 1,
has established itself confidently in the past decade as a replacement to the
once popular hexagonal grid model, in which the base stations were placed at
the centers of the hexagonal lattices. It has been shown [11] that the PPP
approach provides much more accurate results than the hexagonal grid model
when both are used to model real world cellular installation. It has also been
shown [12] that the PPP model gives lower bounds, whereas the hexagonal model
gives upper bounds, of the coverage probability, which means that the former
is safer to rely on.

UEBS Figure 1. A Poisson-Voronoi tessellation model of a cellular network. The
BSs, induced by a PPP, partition the plane into Voronoi cells, each having one
BS and a number of UEs, with the property that a UE in some cell is closer to
the BS of that cell than to the BS of any other cell. This model is so natural
and realistic for cellular networks that it has superseded the once popular
model of hexagonal cells, which was artificial and idealistic.

This spatial consideration afforded by stochastic geometry is a key paradigm
shift in communications systems performance evaluation, where time averages
have traditionally reigned supreme. Before the advent of stochastic geometry,
performance analysis of wireless networks was performed using either damaging
mathematical simplifications or exhaustive simulations with the aim to average
out the many sources of randomness, such as BS and user equipment (UE)
locations and fading phenomena. The mathematical simplifications gave rise to
poor analytical results.

An important technique used in cellular networks, and considered in the
present article, is fractional power control (FPC), where each UE adjusts its
power level in the UL direction under control of its serving BS [8]. More
specifically, each UE controls its transmit power such that the received
signal power at its serving BS is equal to a predefined threshold. This
control benefits the UE by saving its energy and benefits the BS by decreasing
its SIR. The optimal levels of transmit power in a network depend on path
loss, shadowing, and multipath fading, as well as the network configuration.
It should be noted, however, that while FPC in homogeneous cellular networks,
it is less effective in heterogeneous networks (HetNets), where there are more
than one tier of BSs, because far UEs, especially those at the edge of the
cell, will use more transmit power, causing more interference to neighbor
cells.

We prove in this article two theorems that expose a prevalent flaw prevalent
in the stochastic geometric characterization of coverage probability, namely
that the latter is dependent of BS density. The Theorems, one for DL and one
for UL, prove that the opposite is true; namely that the coverage probability
is independent of BS density. This fact has only been alluded to in some
articles previously [13, 14]. The authors of [13] in the context of analyzing
DL coverage probability of millimeter-wave cellular networks noted that
“coverage does not scale with BS density.” Also, the authors of [14], while
analyzing UL FPC, noted that coverage is “ invariant to the density of
deployment of BSs when the shadowing is mild and power control is fractional.”
To the best of our knowledge, the present article is the first to provide a
rigorous mathematical proof that coverage probability is independent of BS
density.

The rest of the article is organized as follows. Section 2 reviews the
stochastic geometric models of cellular networks in the recent literature.
Section 3 provides the modelling process of a cellular system in both
directions, downlink (DL) and uplink (UL). The simulation algorithms developed
to validate the mathematical models are presented in Section 4. Some numerical
results for sample systems are presented in Section 5. Section 6 has the
conclusions.

##  2\. Related work

Perhaps the precursor of using stochastic geometry as a modelling vehicle for
wireless systems was the work in [15], approximately half a century ago.
However, the real momentum, especially for cellular networks, has been gained
over the past decade which witnessed an outbreak of hectic research in the
field. This research can be categorized in many ways, but the direction of
transmission, i.e. DL, UL or both, seems suitable enough.

Much work has focused on DL coverage probability. For example, in [16], the
problem of DL data transmission scheduling in wireless networks is studied,
proposing two algorithms concerning power assignment. In [13], the authors
study DL coverage probability in the context of millimeter-wave (mmWave)
cellular networks and note that mmWave coverage does not scale with BS
density. In [17], DL coverage probability is studied with composite
$\kappa-\mu$ shadowed and lognormal shadowed fading, also known as doubly-
selective (DS) fading. In [18], the authors derive an expression for DL
coverage probability, taking into account a queue setup. The static properties
of the physical layer of the network are studied by stochastic geometry and
the dynamic properties of the queue setup are studied with a discrete time
Markov chain. Upper and lower bounds on the dynamic coverage probability are
derived. In [7], the authors study DL in a setting of an ultra-dense network,
a promising technology for 5G cellular networks, where lots of low power Small
Base Stations (SBSs) overlapping with Macro Base Stations (MBSs) are deployed.
The authors derive the coverage probability.

Many authors have analyzed the DL coverage probability in the context of
HetNets. For example, in [19], the authors characterize the DL coverage
probability in the context of a 2-tier HetNet, under Line of Site (LOS) and
Non Line of Site (NLOS) Path loss models with different path loss exponents,
using stochastic geometry. Also working in HetNets are the authors of [20],
who consider interference to be a major performance bottleneck. They consider
$K$-tier transmission modeled by factorial moment and stochastic geometry to
obtain expressions for the coverage probability. They also compare the model
with a single-tier, traditional hexagonal grid model, concluding the
superiority of the former. Also, in [21], DL is studied using the stochastic
geometry for a two-tier ultradense HetNet with small-cell base stations (SBSs)
and UEs densely deployed in a traditional macrocell network. Performance is
studied in terms of the association probability, average link spectral
efficiency (SE), and average DL throughput. In [22], the authors study DL for
HetNets, where in each typical macrocell, many femtocells and few picocells
are deployed. They reveal that the coverage probability of the DL direction
decreases when the FPC increases in the UL direction, mentioning nothing of
the deployment density. In [23], the authors note that in current cellular
networks, cell association is heavily based on the DL signal power and that
all devices are associated with the same BS in both DL and UL, noting that
this technique is inappropriate in HetNets as transmission levels
significantly vary from BS to BS. They suggest decoupling DL and UL as a
solution, where the UL cell association is not necessarily based on the same
criteria as the DL association. The authors of [24] consider also the same
decoupling solution to alleviate the load imbalance problem in HetNets. The
BSs in each tier are modeled by an independent, homogeneous PPP, whereas the
spatially clustered UEs are modeled by a Matérn cluster process. Analytical
expressions for the coverage probability are derived, showing the merit of the
solution.

Like DL, a huge body of research work has focused on UL coverage probability
in different scenarios. For example, in [25], the authors study UL for a radio
frequency (RF) wake-up solution for Internet of Things (IoT) devices over
cellular networks. When the IoT device has no data to transmit, it turns its
main circuitry OFF and switches to a sleep mode. The transition back to the
active mode is only achieved upon receiving enough power at the device’s front
end. They derive an expression for the UL coverage probability after
successful wake up. In [5], the authors study UL in a densely deployed
cellular network in an interference limited scenario, meaning that the thermal
noise power is very small if compared with the interference. The authors
characterize the SIR, noting it is much more tractable than SINR in a Poisson
cellular network. In [26], UL is studied in the context of mmWave cellular
networks. The coverage probability is characterized under the assumption that
the BSs and UEs are modeled as two independent homogeneous PPPs.
Interestingly, the authors find no clear and consistent correlation between
the coverage probability and BS density, noting that the probability increases
with increasing BS density at first, then decreases as the network becomes
denser and eventually collapse as the density is further increased. In [27],
the authors employ stochastic geometry to study the SINR for UL Poisson
cellular network, where truncated fractional transmit power control is
performed. They obtains the UL SINR distribution under channel-adaptive user
scheduling, including cases in which edge users are both allowed and not
allowed to transmit at maximum transmit power. In [14], the authors analyze
three UL transmit power control schemes, assuming composite Rayleigh-lognormal
fading. Using stochastic geometry, they derive network coverage probability,
noting that it is highly dependent on the severity of shadowing, the power
control scheme, but _invariant_ of the density of deployment of base stations
when the shadowing is mild and power control is fractional. This note is
relevant to and consistent with the theorems proved in the present work.

Also, like DL, many authors have studied UL coverage probability for HetNets.
For example, in [8], the authors study SIR for HetNets leveraging FPC. They
propose nonuniform SBS deployment (NU-SBSD) to reduce interference, where SBS
deployment (SBSD) near macro base station (MBS) is avoided, and MBS coverage
edge area is enriched with ultra-dense SBSD. Also, in [28], the authors
characterize the UL outage probability, which is the complement of coverage
probability, in multitier millimeter wave cellular networks, using stochastic
geometry. In each tier, the BSs are assumed to have their spatial density,
antenna gain, receiver sensitivity, blockage parameter, and pathloss
exponents. The results show that imposing a maximum power constraint on the
user significantly affects the SINR and outage probability. Moreover, in [29],
the authors indicate that in the conventionally coupled UL and DL association
(CUDA), the UL performance is limited greatly by the DL parameters such as the
density and power of base stations (BSs). To overcome this issue, they focus
on a three-tier HetNet as well as global performance evaluation, where the
cross-tier dual-connectivity (DC) and decoupled UL and DL association (DUDA)
are integrated. In [30], UL is studied for a two-tier dense HetNet, where
decoupled UL coverage probability is computed using multi-slope path loss
model. The authors find that the decoupled UL coverage probability is higher
when incorporating multi-slope path loss model as compared to single-slope
path loss model, while the decoupled UL spectral efficiency is observed to be
lower when incorporating dual slope path loss model.

Finally, some research work, such as ours, focuses on both DL and UL,
targeting comprehensive treatment. For example, in [31], an authoritative
tutorial on cellular network analysis using stochastic geometry is provided,
characterizing DL and UL coverage probabilities after a lavish introduction to
stochastic geometry. Also, in [32], both DL and UL are investigated in the
context of mmWave communication, known for being sensitive to blockage. The
authors consider a Time-Division Duplex (TDD) mode able to provide dynamic UL-
DL configurations, and compute the coverage probability based on stochastic
geometry. In [33], both DL and UL are studied in the context of full duplex
(FD) for a UE capable of transmitting and receiving data simultaneously in the
same frequency resource. The authors obtain the coverage probabilities,
revealing that while FD improves DL performance, it severely hurts UL
performance. Furthermore, in [34], the authors study both DL and UL for
massive multiple-input multiple-output (M-MIMO), promising for increasing the
spectral efficiency. They analyze an asynchronous DL M-MIMO system in terms of
the coverage probability by means of stochastic geometry.

Some authors investigate coverage probability for both DL and UL, but in the
context of HetNets. An interesting example is in [35] where the authors
consider DL simultaneous information transmission and power transfer and UL
information transmission for an unmanned aerial vehicle (UAV)-assisted mmWave
cellular network. Distinguishing features of mmWave communications, such as
different path loss models and directional transmissions are taken into
account. DL association probability and energy coverage of the two tiers of
UAVs and ground base stations (GBSs) are investigated.

Still, some other authors study both DL and UL for the purpose of decoupling,
aiming to alleviate BS load imbalance and also ease interference. For example,
in [23], the authors note that in current cellular networks, cell association
is heavily based on the DL signal power and that all devices are associated
with the same BS in both DL and UL. While this practice is adequate in
homogeneous networks, where all BSs have similar transmission levels, it can
fail in dense HetNets as transmission levels significantly vary from BS to BS.
They propose decoupling DL and UL as a solution, where the UL cell association
is not necessarily based on the same criteria as the DL association. Also, the
authors of [24] consider DL and UL decoupling is investigated to alleviate the
load imbalance problem in HetNets, aiming to increase the coverage
probabilities and data rates. The BSs in each tier are modeled by an
independent, homogeneous PPP, whereas the spatially clustered UEs are modeled
by a Matérn process. Analytical expressions for the coverage probability, are
derived.

##  3\. modelling

The subject of stochastic geometry and its use in modelling wireless systems
is vast [10]. The key aspect of the present study is that all the BSs are
located according to a PPP $\Phi$, which effectively means they are randomly
scattered in the Euclidean plane with independent locations. At any given
time, only one UE can be active communicating with the BS on any
time/frequency resource. Before delving even more into the modelling process,
some definitions used throughout the article are in order.

:

Definition 1 (BS-UE association): BS-UE association is the assignment of a UE
to a BS, for both to establish a communications session.

:

Definition 2 (Serving BS): Once a UE is associated with a BS, the latter is
said to be the serving BS of the UE.

:

Definition 3 (Typical receiver): The typical receiver is the receiving device
(UE or BS) where the SIR is to be assessed. It is always placed at the origin
of the Euclidean plane in the model, or the origin of the simulation window in
the simulation.

:

Definition 4 (Tagged transmitter): The tagged transmitter is transmitting
device (UE or BS) associated with the typical receiver.

:

Definition 5 (Typical circle): The typical circle is the circle centered at
the typical receiver and having the tagged transmitter on its circumference.

:

Definition 6 (Signal): A signal is the transmission arriving at, and intended
for, the typical receiver.

:

Definition 7 (Interference): An interference is the transmission arriving at,
but not intended for, the typical receiver.

:

Definition 8 (Interferer): An interferer is a transmitter causing interference
at the typical receiver. That is, it is any transmitter in the network other
than the tagged transmitter.

:

Definition 9 (Signal to Interference Ratio (SIR)): The quotient of the signal
at the typical receiver and the sum of all interferences at the typical
receiver.

:

Association rule: A UE will associate with the BS nearest it.

As per Definitions 3 and 4, in DL the typical receiver is a UE and the tagged
transmitter is a BS, whereas in UL the typical receiver is a BS and the tagged
transmitter is a UE.

As per the association rule, a UE will be nearer to its serving BS than to any
other BS in the cellular network. Strictly speaking, however, the UE
associates with the BS that provides the highest average SIR. Surprisingly,
depending on the fading and shadowing conditions, a more distant BS could
provide a higher instantaneous SIR, hence gets associated with the UE. That is
why, we rely on the average, which invariably translates to association with
the nearest BS. We will denote the distance between the associated BS-UE pair
throughout by $R$. If the BS density is $\lambda$, it can be shown that $R$ is
a random variable (RV) with the Rayleigh distribution

$f_{R}(r)=2\lambda\pi re^{-\lambda\pi r^{2}},\qquad r\geq 0$ (1)

Along the same line, the notation used throughout the article is provided in
Table 1. The notation of the PPP, pivotal in our work, depends on its setting
and the field it is being applied in and on the interpretation of the process.
For example, a simple PPP $\Phi$ may be considered as a random set, which
suggests the notation $x\in\Phi$, implying that $x$ is a random point
belonging to or being an element of the Poisson random set $N$. Another, more
general, interpretation is to consider a Poisson or any other point process as
a random counting measure, so one can write the number of points of a PPP
$\Phi$ being found or located in some (Borel measurable) region $B$ as
$\Phi(B)$, which is an RV (having a Poisson distribution, in the case of a
PPP) [10].

Table 1. Notation used in the model and simulation. Parameter | Description
---|---
BS | Base station
UE | User Equipment (can be a mobile phone, tablet, laptop, etc.)
$\Phi$ | Poisson point process (PPP) of BSs
$\Psi$ | Point process of UEs (not Poisson)
$\lambda$ | Density of BS (per m2), i.e. intensity of PPP $\Phi$
$\alpha$ | Path-loss exponent (per m)
SIR | Signal to interference ratio (dB)
$\xi$ | SIR threshold (dB)
$G$ | Rayleigh channel gain of tagged transmitter ($G\sim Exp(1)$)
p | Transmit power (Watts)
$p_{d}$ | DL coverage probability
$p_{u}$ | UL coverage probability
$\mathcal{N}$ | Number of simulation runs (PPP realizations)
$S$ | Side of the simulation square window (m)
$N$ | Number of BSs in a simulation run, $N\sim Pois(\lambda S^{2})$).
($x,y$) | Location of a BS in a simulation run
($u,v$) | Location of a UE in a simulation run
Covered | No. of times typical receiver is covered by tagged transmitter

Figure 2 shows a stochastic geometric construction to characterize the SIR at
the typical receiver of a cellular network, which with the tagged transmitter,
defines the typical circle. In part 2(a), we can see the DL model, where the
typical receiver is a UE, and the tagged transmitter is a BS at distance $R$.
All the BSs outside the typical circle are interferers to the typical UE. The
typical circle defines an exclusion zone, in the sense that there can be no BS
inside of it (or else the typical UE would associate with it.) In part 2(b),
we can see the UL model, where the typical receiver is a BS, and the tagged
transmitter is a UE at distance $R$. All the UEs except the tagged are
interferers to the typical BS. The typical circle does not define an exclusion
zone, as there can be UEs inside of it as we can see (the association problem
of the DL model is not present here.)

(a) (a) $z$ $R$ $D_{z}$ $R_{z}$

(b) (b) $\mathfrak{z}$ $R$ $U_{\mathfrak{z}}$ $R_{\mathfrak{z}}$

Figure 2. Stochastic geometric models for assessing the SIR at a typical
receiver at the origin of the cellular network. (a) The DL model, where the
typical receiver is a UE. The typical circle defines an exclusion zone, as it
cannot contain a BS inside. The BSs outside the typical circle, which cause
interference at the typical UE, form a PPP $\Phi$. (b) The UL model, where the
typical receiver is a BS. The typical circle does not define an exclusion
zone, as it can contain UEs inside. All the UEs form a PP $\Psi$. Since the
UEs are _s_ atellites to their serving BSs, as per the association rule,
$\Psi$ is not Poisson, a major challenge to the analysis. We mitigate this
challenge by relocating each UE, except the tagged, to the position of its
serving BS. That is, we relocate all the interferers.

The RV $R_{z}$ are identically distributed but not independent in general. The
dependence is induced by the structure of Poisson-Voronoi tessellation and the
restriction that only one BS can lie in each cell. To visualize this
dependence, recall that the presence of a BS in a particular Voronoi cell
forbids the presence of any other BS in that cell. However, as discussed in
detail in the next section, this dependence is weak, motivating the
approximation we will make in the UL model below. On the other hand, the RV
$R_{z}$ is upper bounded by the RV $D_{z}$, even for the pairs inside the
typical circle (or else the BS at $z$ would be closer to the typical UE,
violating the association rule.). The UEs are satellites to the points of the
PPP $\Phi$, based on the association rule, so they are not themselves a PPP,
making the UL model a challenge. This challenge is mitigated in the present
work through an elegant approximation.

We can relate Figure 1 and Figure 2 by recognizing that in the present
article, we consider orthogonal communications. This means that in each cell
there can be only one active UE on any time/frequency resource. Accordingly,
Figure 1 is a snapshot of the UEs that are active in all the cells at the same
time. Thus, we see a random number of UEs in each cell. Figure 2, on the other
hand, is a snapshot of the UEs that are active on the same frequency in all
the cells at the same time. Thus, we see exactly one UE per cell. We focus in
the model on one such resource. That is, every BS-UE associated pair in Figure
2 is operating on the same resource, hence is the interference that we are
going to characterize.

Given a homogeneous PPP with intensity $\lambda$, the number of process points
in the an area of size $A$, denoted by $\Phi(A)$, is Poisson distributed

$\mathbb{P}[\Phi(A)=k]=\frac{e^{-\lambda|A|}(\lambda|A|)^{k}}{k!}$

That is, a homogeneous PPP is completely characterized by a single number,
$\lambda$. It is both stationary (translation invariant) and isotropic
(rotation invariant).

Similarly, if $F$ is a continuous RV, we can also write its Laplace transform
as a (continuous) expectation, namely

$\mathcal{L}_{F}(s)=\mathbb{E}\left[e^{-sF}\right]$ (2)

The probability generating function of a non-negative integer-valued RV leads
to the probability generating functional (PGFL) with respect to any non-
negative bounded function $f$ for a point process defined on $\mathbb{R}^{d}$,
with $0\leq f(x)\leq 1$. The need for the PGFL, which we will invoke twice in
this article, arises naturally when one wishes to find the Laplace transform
of a sum of functions. As the Laplace transform of a RV is already an
expectation of an exponential function, with the RV its exponent, then if the
RV is a sum of independent and identically distributed (iid) RVs, the
exponential becomes a product lending itself readily to the invocation of
PGFL. Due to this intimate relationship, some authors call the PGFL in this
context a Laplace functional.

We will also use in the present work Slivnyak’s theorem [10]. It states that
for a PPP, because of the independence between its points, conditioning on a
point at some location $x$ does not change the distribution of the rest of the
process. That is, removing an infinitesimally small area from the underlying
plane of the PPP, leaves still a PPP, since the distribution of points in all
nonoverlapping areas are independent for the PPP. Consequently, any property
seen from some location $x$ remains the same whether or not there is a point
at $x$. Though simple, this theorem allows us to add and remove points to/from
the PPP liberally. In our work, for example, it allows us in UL to place a BS
from those imposed by the PPP $\Phi$ and place it at the origin, treating the
remaining points, $\Phi\backslash x$, still as a PPP (for example, when
invoking a PGFL.)

Before starting the modelling process, a word on random channels is in order.
In our work, random channel effects are incorporated by multiplicative RVs,
namely $G$ for the signal and $G_{i}$ for each interferer $i$. For simplicity
we assume these all correspond to Rayleigh fading with mean $1$, denoted
$\exp(1)$, i.e. $\mu=1$ in $f_{G_{i}}(x)=\mu e^{-\mu x}$. We consider small-
scale Rayleigh fading, and assume $G$ and the $G_{i}$ iid RVs following an
exponential distribution with mean $1$.

###  3.1. Downlink model

The key assumptions of the DL system model are as follows:

  * •

BSs are located according to a homogeneous PPP $\Phi$ of intensity $\lambda$
in the Euclidean plane.

  * •

A UE wishing to start a communications session associates with the BS that is
closer to it than any other BS in the cellular network (association rule.)

  * •

The BS transmits at a fixed power $p$ to a designated UE on a particular time-
frequency resource, i.e. orthogonal multiple access communications within the
cell take place. The consequence of orthogonality is that the UE sees
interference from all other BSs in the plane using the same resource.

  * •

Signals attenuate with distance according to the standard power-law path loss
propagation model, with path loss exponent $\alpha>2$. That is, the average
received power at distance $r$ from a transmitter of power $p$ is
$pr^{-\alpha}$.

  * •

Random channel effects are incorporated by a multiplicative RVs $G$ for the
signal and $G_{z}$ for every interferer at $z$. For simplicity we assume these
all to correspond to Rayleigh fading with mean $1$, so $G$ and the $G_{z}$ are
iid RVs, having exponential distribution with mean $1$.

Let $I_{d}$ denote the interference experienced at the typical UE. The
interference is due to every BS $z$ in the plane, except the tagged BS,
denoted by $\mathfrak{b}$, at distance $D_{z}$ from the typical UE, as shown
in Figure 2(a),. That is, the interference $I_{d}$ is created by a PPP with
intensity $\lambda$ outside the typical circle, and is given by

$I_{d}=\sum_{z\in\Phi\backslash\\{\mathfrak{b}\\}}pG_{z}D_{z}^{-\alpha},$ (3)

where $p$ is the power of the BS at point $z$, and $D_{z}$ is a RV
representing the distance from the BS at $z$a and the typical UE. The SIR at
the typical UE is then given by

$\mathtt{SIR}_{\text{UE}}=\frac{pGR^{-\alpha}}{I_{d}}$ (4)

Note that the distance $R_{z}$ is not used in anything in DL, but it will be
used in the UL model. Note also that the typical circle defines an exclusion
zone. That is, there can never be a BS closer to the typical UE than the
tagged BS. For if there were such a BS, it would be the tagged BS associated
with the typical UE. Later in the UL model, we will see that we cannot define
an exclusion zone by any BS-UE pair.

The goal now is to derive the DL coverage probability $p_{d}$ in the cellular
network, which is exactly the complementary cumulative distribution function
(CCDF) of SIR over the entire network (Recall that the CDF gives
$\mathbb{P}[$SIR$\leq\xi]$). The coverage probability can also be visualized
as the probability that a randomly chosen user can achieve a target SIR $\xi$,
the average fraction of users who at any time achieve SIR $\xi$, or the
average fraction of network area that is in “coverage” at any time. We have
two routes now for obtaining $p_{d}$ as far as the manipulation of RVs is
concerned. Either use RVs throughout and keep resolving them via the
$\mathbb{E}$ operator vehicle, or condition the SIR expression at the very
beginning, using particular values for those RVs, then keep on deconditioning.
We will take the former route. At the outset, We will invoke the concept of
total probability, using (3) and (4), to get

$\displaystyle p_{d|R}$ $\displaystyle=$
$\displaystyle\mathbb{P}[\mathtt{SIR}_{\text{UE}}>\xi]$ $\displaystyle=$
$\displaystyle\mathbb{P}\left[\frac{pGR^{-\alpha}}{I_{d}}>\xi\right]$
$\displaystyle=$
$\displaystyle\mathbb{P}\left[G>\frac{\xi}{p}R^{\alpha}I_{d}\right]$
$\displaystyle\overset{(\text{a})}{=}$
$\displaystyle\mathbb{E}_{I_{d}}\left[\mathbb{P}\left[G>\frac{\xi}{p}R^{\alpha}I_{d}\right]\right]$
$\displaystyle\overset{(\text{b})}{=}$
$\displaystyle\mathbb{E}_{I_{d}}\left[e^{-\frac{\xi}{p}R^{\alpha}I_{d}}\right]$
$\displaystyle\overset{(\text{c})}{=}$
$\displaystyle\mathcal{L}_{I_{d}}(\frac{\xi}{p}R^{\alpha})$ (5)

where

$\mathcal{L}_{A}(s)=\int_{0}^{\infty}e^{-st}f_{A}(t)dt=\mathbb{E}\left[e^{-sA}\right]$
(6)

is the Laplace transform of the RV $I_{d}$ conditioned on the RV $R$ between
the typical UE and the tagged BS. In (a) we utilized the fact that we can
write a probability $\mathbb{P}\left[A>B)\right]$ as
$\mathbb{E}_{B}\left[\mathbb{P}\left[A>B\right]\right]$ (or
$\mathbb{E}_{A}\left[\mathbb{P}\left[A>B\right]\right]$), in (b) we benefited
from the fact that $G\sim\exp(1)$, i.e. $f_{G}(r)=e^{-r}$, and in (c) we used
the Laplace transform definition (6).

Clearly, the DL coverage probability $p_{d}$ in (5) is conditioned on $R$, the
distance between the typical UE and the tagged BS. We will now embark on
deconditioning $p_{d}$ (basically by obtaining its expectation with respect to
$R$.) Since $R$ is the distance between the typical UE and the closest BS (the
tagged BS), it is Rayleigh distributed, i.e. and $f_{R}(r)=2\lambda\pi
re^{-\lambda\pi r^{2}}$ from (1). Further, $R$ ranges from an arbitrarily
small positive real number greater than $0$ (to exclude the typical UE) to
$\infty$. Thus, the conditional coverage probability

$\displaystyle p_{d|R}$ $\displaystyle=$
$\displaystyle\mathbb{E}_{R}\left[p_{d|R}\right]$ (7) $\displaystyle=$
$\displaystyle\mathbb{E}_{R}\left[\mathcal{L}_{I_{d}}(\frac{\xi}{p}R^{\alpha})\right]$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\mathcal{L}_{I_{d}}(\frac{\xi}{p}r^{\alpha})f_{R}(r)dr$
$\displaystyle=$ $\displaystyle 2\lambda\pi\int_{0}^{\infty}e^{-\lambda\pi
r^{2}}\mathcal{L}_{I_{d}}(\frac{\xi}{p}r^{\alpha})rdr$

Next, we will embark on finding the Laplace transform $\mathcal{L}_{I_{d}}$ of
the DL interference $I_{d}$. Using (3) and (6), we get

$\displaystyle\mathcal{L}_{I_{d}}(s)$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[e^{-sI_{d}}\right]$ $\displaystyle=$
$\displaystyle\mathbb{E}_{\Phi,G_{z}}\left[e^{-s\sum_{z\in\Phi\backslash\\{\mathfrak{b}\\}}pG_{z}D_{z}^{-\alpha}}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\Phi,G_{z}}\left[\prod\limits_{z\in\Phi\backslash\\{\mathfrak{b}\\}}e^{-spG_{z}D_{z}^{-\alpha}}\right]$
$\displaystyle\overset{(\text{a})}{=}$
$\displaystyle\mathbb{E}_{\Phi}\left[\prod\limits_{z\in\Phi\backslash\\{\mathfrak{b}\\}}\mathbb{E}_{G_{z}}\left[e^{-spG_{z}D_{z}^{-\alpha}}\right]\right]$
$\displaystyle\overset{(\text{b})}{=}$
$\displaystyle\mathbb{E}_{\Phi}\left[\prod\limits_{z\in\Phi\backslash\\{\mathfrak{b}\\}}\mathcal{L}_{G_{z}}(spD_{z}^{-\alpha})\right]$
$\displaystyle\overset{(\text{c})}{=}$
$\displaystyle\exp\left(-\lambda\int_{\mathbb{R}^{2}\backslash
D(o,r)}(1-\mathcal{L}_{G_{z}}(spD_{z}^{-\alpha}))\right)$ (8)

where $D(o,r)$ is a disc centered at the origin and has a radius $r$. In (a)
we benefited from the independence of the $G_{z}$, which are iid and in (b) we
used the definition (6) of the Laplace transform. In (c), to decondition on
$D_{z}$ which is distributed _differently_ for each point $z$ of the PPP, we
invoked the PGFL $\mathbb{E}_{\Phi}\left[\prod\limits_{z\in\Phi}f(x)\right]$,
with $f(x)=\mathcal{L}_{G_{z}}(spx^{-\alpha})$, of the PPP $\Phi$. Note that
$\mathbb{E}_{G_{z}}$ will be taken at the level of a single point $z$, whereas
$\mathbb{E}_{\Phi}$ will be taken at the level of the PPP $\Phi$. However, for
the _expectation with respect to_ $\Phi$, we will apply it to
$D_{z}^{-\alpha}$, as the latter involves all the PPP points. Note finally
that some authors call the PGFL the Laplace functional, when the former is
invoked to find a Laplace transform as is the case here.

Note that $D_{z}$ ranges from $r^{+}$ (to exclude the tagged BS) to $\infty$.
That is, the integral is carried out outside the typical circle, namely from
$r$ to $\infty$. Switching to polar coordinates, with the interferer now at
$(x,\theta)\in\mathbb{R}^{2}$, then using the fact that $G_{z}\sim\exp(1)$,
i.e. $f_{G_{z}}(t)=e^{-t}$, then Using (8) yields

$\displaystyle\mathcal{L}_{I_{d}}(s)$ $\displaystyle=$
$\displaystyle\exp\left(-\lambda\int_{0}^{2\pi}\int_{r}^{\infty}(1-\mathcal{L}_{G_{z}}(spx^{-\alpha}))xdxd\theta\right)\text{
}$ (9) $\displaystyle=$
$\displaystyle\exp\left(-2\pi\lambda\int_{r}^{\infty}(\frac{spx^{-\alpha}}{1+spx^{-\alpha}})xdx\right)$

For use in (7), we write this result as

$\displaystyle\mathcal{L}_{I_{d}}(\frac{\xi}{p}r^{\alpha})$ $\displaystyle=$
$\displaystyle\exp\left(-2\pi\lambda\int_{r}^{\infty}(\frac{(\frac{\xi}{p}r^{\alpha})px^{-\alpha}}{1+(\frac{\xi}{p}r^{\alpha})px^{-\alpha}})xdx\right)$
$\displaystyle=$ $\displaystyle\exp\left(-\pi\lambda
r^{2}\xi^{\frac{2}{\alpha}}\int_{\xi^{-2/\alpha}}^{\infty}\frac{1}{1+u^{\alpha/2}}du\right)$

where $u=\left(x/r\right)^{2}\xi^{-\frac{2}{\alpha}}$. Substituting this in
(7), we get

$p_{d}=2\widetilde{\lambda}\int_{0}^{\infty}e^{-\widetilde{\lambda}r^{2}}e^{-\widetilde{\lambda}r^{2}\sqrt[\kappa]{\xi}\int_{\frac{1}{\sqrt[\kappa]{\xi}}}^{\infty}\frac{1}{1+u^{\kappa}}du}rdr$
(10)

where $\widetilde{\lambda}=\lambda\pi$ and $\kappa=\alpha/2$.

From (10), it _appears_ that the DL coverage probability $p_{d}$ is dependent
on the BS density parameter $\lambda$. Indeed, this parameter exists in a
plethora of coverage probability (and derivatives) expressions all over the
literature on stochastic geometric modelling of DL cellular systems (See, for
example, [3, 5, 7, 9, 13, 17, 18, 19, 21, 22, 32, 30, 31, 33, 34, 35]).
However, this dependence is false, as we will show in the next Theorem.

:

Theorem 1: Under the stochastic geometric model of the cellular DL system, the
DL coverage probability $p_{d}$ is independent of the BS density $\lambda$.

Proof: The proof is attained through two changes of variables. Starting with
(10), use the substitution $x=r^{2}$ to get

$\displaystyle p_{d}$ $\displaystyle=$ $\displaystyle
2\widetilde{\lambda}\int_{0}^{\infty}e^{-\widetilde{\lambda}r^{2}}e^{-\widetilde{\lambda}r^{2}\xi^{\frac{2}{\alpha}}\int_{\xi^{-\frac{2}{\alpha}}}^{\infty}\frac{1}{1+u^{\kappa}}du}rdr$
$\displaystyle=$
$\displaystyle\widetilde{\lambda}\int_{0}^{\infty}e^{-\widetilde{\lambda}x\left(1+\sqrt[\kappa]{\xi}\int_{\frac{1}{\sqrt[\kappa]{\xi}}}^{\infty}\frac{1}{1+u^{\kappa}}du\right)}dx$

Now use the substitution $z=\widetilde{\lambda}x$ to get

$\displaystyle p_{d}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}e^{-z\left(1+\sqrt[\kappa]{\xi}\int_{\frac{1}{\sqrt[\kappa]{\xi}}}^{\infty}\frac{1}{1+u^{\kappa}}du\right)}dz$
(11) $\displaystyle=$
$\displaystyle\frac{1}{1+\sqrt[\kappa]{\xi}\int_{\frac{1}{\sqrt[\kappa]{\xi}}}^{\infty}\frac{1}{1+u^{\kappa}}du}$

where $\lambda$ has totally disappeared, proving the theorem. $\blacksquare$

Theorem 1 then calls for a revisit to the DL coverage probability expressions
(and their derivatives) that include $\lambda$. Those expressions exist in a
huge number of articles, such as those cited just above the Theorem. It is
worth examining them for possible elimination of a parameter that may exist
superfluously.

DL Special cases:

First, we will consider the special case of $\alpha=4$, i.e. $\kappa=2$. From
(11), we can obtain (using a computer algebra system, e.g. Maxima (R)), a
simple closed form expression for the DL coverage probability, namely

$\displaystyle p_{d}$ $\displaystyle=$
$\displaystyle\frac{1}{1+\sqrt{\xi}\int_{\frac{1}{\sqrt{\xi}}}^{\infty}\frac{1}{1+u^{2}}du}$
$\displaystyle p_{d}$ $\displaystyle=$
$\displaystyle\frac{1}{1+\sqrt{\xi}\left(\frac{\pi}{2}-\arctan(\frac{1}{\sqrt{\xi}})\right)}$
(12)

Second, we will consider the special case of $\alpha=6$, i.e. $\kappa=3$. From
(11), we can obtain (using a computer algebra system, e.g. Maxima (R)), a not-
so-simple closed form expression for the DL coverage probability, namely

$\displaystyle p_{d}$ $\displaystyle=$
$\displaystyle\frac{1}{1+\sqrt[3]{\xi}\int_{\frac{1}{\sqrt[3]{\xi}}}^{\infty}\frac{1}{1+u^{3}}du}$
$\displaystyle p_{d}$ $\displaystyle=$
$\displaystyle\frac{6c}{6c+\ln\left(c^{2}-c+1\right)+\sqrt{3}\left({\pi}-2\arctan\left(\frac{2c-1}{\sqrt{3}}\right)\right)-2\ln\left(c+1\right)}$
(13)

where $c=1/\sqrt[3]{\xi}$. To further attest to the accuracy of the last two
results, they will be plotted later together with the corresponding simulation
counterparts.

###  3.2. Uplink model

The net interference at the typical BS is the sum of the received
transmissions from all the UEs except the tagged, including those inside the
typical circle. Referring to Figure 2(b), for each UE $\mathfrak{z}\in\Psi$,
we denote its distance to its serving BS by $R_{\mathfrak{z}}$. Although the
RVs $R_{\mathfrak{z}}$ are identically distributed, they are not independent.
However, as shown later in this section, this dependence is weak and we will
henceforth assume the $R_{\mathfrak{z}}$ iid. Under this assumption, we will
first derive the coverage probability for the general distribution of
$R_{\mathfrak{z}}$.

The set of interferers are the points of $\Psi$, which is not a PPP. The
reason is that they are not scattered uniformly in $R^{2}$, but rather
associated to a PPP as satellites. This poses difficulty, as we will not have
the luxury of using such useful tools as the PGFL here. However, we can get
around this difficutly as follows. Note that each point of $\Psi$ is
associated, by being closest, to a point in the PPP $\Phi$ of BSs, which we
used above in the down link analysis. More importantly, as each point in
$\Psi$ is located somewhere around the corresponding point in $\Phi$, we can
approximate the ”spatial” average of the former to be the latter.
Consequently, we can approximate the locations of the interfering UEs $\Psi$
by the points of $\Phi$. Accordingly, below we will carry out the sums and
products involved in the stochastic geometry analysis over $\Phi$, rather than
$\Psi$. Specifically, for the sake of calculating the interference, we will
consider that each interfering UE is placed exactly at its serving BS’s
location. Consequently, referring to Figure 2(b), we will employ the distance
$R_{\mathfrak{z}}$ between this UE and its serving BS to calculate its emitted
power. We will then consider this as interference at the typical BS at
distance $D_{z}$ away, not $U_{\mathfrak{z}}$.

In UL, FPC leads to amplifying the transmit power $p$ at the UE based on its
distance to the serving BS. If the distance is $R$ and the FPC factor is
$\epsilon$, with values in $[0,1]$, then $p$ is amplified by
$R^{\epsilon\alpha}$ to offset the path loss, which is $R^{-\alpha}$, where
$\alpha$ is the path loss exponent, with values greater than $2$. Reducing the
power of a UE as it gets closer to its serving BS is useful for two reasons.
First it saves the UE’s battery. Second, it makes the UE less of an interferer
to UEs of neighbor cells. Combining the effects of FPC, power loss and fading,
the amount of power reaching the serving BS from the UE will be
$pGR^{-\alpha(1-\epsilon)}$.

Referring to Figure 2(b), the RV $R_{\mathfrak{z}}$ is upper bounded by
$U_{\mathfrak{z}}$, otherwise the sample UE at $\mathfrak{z}$ would associate
with the typical BS. Moreover, as can be seen, the typical circle does not
define an exclusion zone, as UEs can be located inside of it without breaking
the associatiation rule. Let $I_{u}$ denote the interference caused by all the
UEs, except the tagged, at the typical BS. Accordingly,

$I_{u}=\sum_{\mathfrak{z}\in\Psi\backslash\\{\mathfrak{u}\\}}pG_{\mathfrak{z}}R_{\mathfrak{z}}^{\alpha\epsilon}U_{\mathfrak{z}}^{-\alpha}$
(14)

In the UL, the SIR of the typical BS, at distance $R$ from the tagged UE, is

$\mathtt{SIR}_{\text{BS}}=\frac{pGR^{-\alpha(1-\epsilon)}}{I_{u}}$ (15)

Consequently, the probability $p_{u}$ of UL coverage is

$p_{u}=\mathbb{P}[\mathtt{SIR}_{\text{BS}}>\xi]\text{.}$

This probability can be visualized as being the average area or the average
fraction of users in coverage. As noted earlier, we perform analysis on a
randomly chosen BS assumed to be located at the origin associated with the
closest UE. Under assumption 1, the distribution of the distance $R$ of the
closest mobile from the randomly chosen BS can be assumed Rayleigh given by
(1). This chosen BS is placed at the origin, and is henceforth called the
typical BS.

Referring to Figure 2(b), both $R$ and $R_{\mathfrak{z}}$ are Rayleigh
distributed, i.e. $f_{R}(r)=f_{R_{\mathfrak{z}}}(r)=2\lambda\pi
re^{-\lambda\pi r^{2}}$. Thus, $R_{\mathfrak{z}}$ ranges from $0$ to $R$, with
$R$ ranging from an arbitrarily small positive real number to $\infty$. Now,
the conditional UL coverage probability is defined as

$\displaystyle p_{u|R}$ $\displaystyle=$
$\displaystyle\mathbb{P}[\mathtt{SIR}_{\text{BS}}>\xi]$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[\mathbb{P}\left[\frac{pGR^{-\alpha(1-\epsilon)}}{I_{u}}>\xi\right]\right]$
$\displaystyle=$ $\displaystyle\mathbb{E}\left[\mathbb{P}\left[G>\xi
p^{-1}R^{\alpha(1-\epsilon)}I_{u}\right]\right]$
$\displaystyle\overset{(\text{a})}{=}$ $\displaystyle\mathbb{E}\left[e^{-\xi
p^{-1}I_{u}R^{\alpha(1-\epsilon)}}\right]$ $\displaystyle=$
$\displaystyle\mathcal{L}_{I_{u}}(\xi p^{-1}R^{\alpha(1-\epsilon)})$

where $\mathcal{L}_{I_{u}}$ is the Laplace transform of the distribution of
the $I_{u}$ RV. In (a), we used the fact that $G\sim\exp(1)$, i.e.
$f_{G}(x)=e^{-x}$, which implies that $\mathbb{P}[G>x]=e^{-x}$. Now, we
decondition on $R$, getting

$\displaystyle p_{u}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\left.\mathcal{L}_{I_{u}}(\xi
p^{-1}R^{\alpha(1-\epsilon)})\right|_{R=r}f_{R}(r)dr$ (16) $\displaystyle=$
$\displaystyle\int_{0}^{\infty}2\widetilde{\lambda}re^{-\widetilde{\lambda}r^{2}}\mathcal{L}_{I_{u}}(\xi
p^{-1}r^{\alpha(1-\epsilon)})dr$

where $\widetilde{\lambda}=\pi\lambda$. We integrate from a point just outside
the origin, to skip the typical BS that resides there, to $\infty$ where the
closest UE can possibly exist.

Next, we will embark on finding $\mathcal{L}_{I_{u}}$, the Laplace transform
of the distribution of the RV $I_{u}$. Substituting for $I_{u}$ from (14),
gives

$\displaystyle\mathcal{L}_{I_{u}}(s)$ $\displaystyle=$
$\displaystyle\mathbb{E}[e^{-sI_{u}}]$ (17) $\displaystyle=$
$\displaystyle\mathbb{E}\left[\exp\left(s\sum_{\mathfrak{z}\in\Psi}-pG_{\mathfrak{z}}R_{\mathfrak{z}}^{\alpha\epsilon}U_{\mathfrak{z}}^{-\alpha}\right)\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\prod\limits_{\mathfrak{z}\in\Psi}\exp\left(-spG_{\mathfrak{z}}R_{\mathfrak{z}}^{\alpha\epsilon}U_{\mathfrak{z}}^{-\alpha}\right)\right]$

In (17), for each point $\mathfrak{z}\in\Psi$ there are three RVs:
$G_{\mathfrak{z}},R_{\mathfrak{z}},U_{\mathfrak{z}}$. The $G_{\mathfrak{z}}$
are independent of the $R_{\mathfrak{z}}$ and of the $U_{\mathfrak{z}}$.
However, $U_{\mathfrak{z}}$ and $R_{\mathfrak{z}}$ are dependant in that
$R_{\mathfrak{z}}<U_{\mathfrak{z}}$ (Recall that $R_{z}$ is the distance
between an interfering UE $\mathfrak{z}$ and its typical BS, and
$U_{\mathfrak{z}}$ is the distance between the same interfering UE
$\mathfrak{z}$ and the typical BS at the origin). That is
$\mathbb{P}[R_{\mathfrak{z}}<x|U_{\mathfrak{z}}=x]=1$, since if
$U_{\mathfrak{z}}<R_{\mathfrak{z}}$ the interfering UE $\mathfrak{z}$ would
associate with the typical BS at the origin.

We will next apply the expectation operator $\mathbb{E}$ to the three RVs in
(17), one at a time, starting with the one that is independent of the other
two, namely $G_{\mathfrak{z}}$. Only _one_ of these expectations, the one with
respect to the PPP $\Psi$, will be evaluated using a PGFL and the others using
the standard definition of an expectation. Note that
$\mathbb{E}_{G_{\mathfrak{z}}}$, $G_{\mathfrak{z}}\sim{Exp}(1)$, i.e.
$f_{G_{\mathfrak{z}}}(x)=e^{-x}$.

$\displaystyle\mathcal{L}_{I_{u}}(s)$ $\displaystyle=$
$\displaystyle\mathbb{E}_{\Psi,R_{\mathfrak{z}},G_{\mathfrak{z}}}\left[\prod\limits_{\mathfrak{z}\in\Psi}e^{-spG_{\mathfrak{z}}R_{\mathfrak{z}}^{\alpha\epsilon}U_{\mathfrak{z}}^{-\alpha}}\right]$
$\displaystyle\overset{(\text{a})}{=}$
$\displaystyle\mathbb{E}_{\Psi,R_{\mathfrak{z}}}\left[\prod\limits_{\mathfrak{z}\in\Psi}\mathbb{E}_{G_{\mathfrak{z}}}\left[e^{-spG_{\mathfrak{z}}R_{\mathfrak{z}}^{\alpha\epsilon}U_{\mathfrak{z}}^{-\alpha}}\right]\right]$
$\displaystyle\overset{(\text{b})}{=}$
$\displaystyle\mathbb{E}_{\Psi,R_{\mathfrak{z}}}\left[\prod\limits_{\mathfrak{z}\in\Psi}\int_{0}^{\infty}e^{-\left(1+spR_{\mathfrak{z}}^{\alpha\epsilon}U_{\mathfrak{z}}^{-\alpha}\right)x}dx\right]$
(18) $\displaystyle=$
$\displaystyle\mathbb{E}_{\Psi,R_{\mathfrak{z}}}\left[\prod\limits_{\mathfrak{z}\in\Psi}\frac{1}{1+spR_{\mathfrak{z}}^{\alpha\epsilon}U_{\mathfrak{z}}^{-\alpha}}\right]$

In (a) we used the fact that the $G_{\mathfrak{z}}$ are iid and in (b) we used
the fact that $f_{G_{\mathfrak{z}}}(x)=e^{-x}$.

Now, we consider the expectation with respect to $\Psi$, to uncondition on
$U_{\mathfrak{z}}$, the distance between every point $\mathfrak{z}\in\Psi$ and
the origin. We will use for this expectation a PGFL, since $U_{\mathfrak{z}}$
is distributed differently for each point $\mathfrak{z}\in\Psi$. Based on the
approximation we are proposing, the points of the non-Poisson $\Psi$ PP are
now relocated to the points of the (homogeneous) Poisson $\Phi$ PP. In
particular, each point $\mathfrak{z}\in\Psi$ will be relocated to the position
of the associated point $z\in\Phi$. That is, we will consider each interfering
UE at point $z\in\Phi$ emitting power $pR_{\mathfrak{z}}^{\alpha\epsilon}$,
but causing interference with this same power at the typical BS, at a distance
$D_{z}$ based on the UE relocation. This allows us to write

$\mathbb{E}_{\Psi}\left[\prod\limits_{\mathfrak{z}\in\Psi}f(\mathfrak{z})\right]\approx\mathbb{E}_{\Phi}\left[\prod\limits_{z\in\Phi}f(z)\right]=e^{-\lambda\int_{\mathbb{R}^{2}}(1-f(x))}\text{.}$

Substituting for $f(y)$ from (18), converting to polar coordinates, and
substituting for the angle integral by $2\pi$, then

$\displaystyle\mathcal{L}_{I_{u}}(s)$ $\displaystyle=$
$\displaystyle\mathbb{E}_{R_{\mathfrak{z}}}\left[\mathbb{E}_{\Phi}\left[\prod\limits_{z\in\Phi}\frac{1}{1+spR_{\mathfrak{z}}^{\alpha\epsilon}D_{z}^{-\alpha}}\right]\right]$
(19) $\displaystyle=$
$\displaystyle\mathbb{E}_{R_{\mathfrak{z}}}\left[e^{-2\pi\lambda\int_{0}^{\infty}\frac{1}{1+\left(sp\right)^{-1}R_{\mathfrak{z}}^{-\alpha\epsilon}x^{\alpha}}xdx}\right]$

The integral includes the entire Euclidean space, so $r$ goes from $0$, ”very
near” the origin to avoid having a BS at the origin where interference is
assessed, to infinity. However, we approximate this nearness by $0$, not
fearing a singularity in the denominator of the integrand since there is a
$1$.

Note that $R_{\mathfrak{z}}$ is lower bounded by $U_{\mathfrak{z}}$, for if
$R_{\mathfrak{z}}<U_{\mathfrak{z}}$, the UE at $\mathfrak{z}$ would associate
with the typical BS. But note that $U_{\mathfrak{z}}$ has been replaced now,
through PP relocation, by $D_{\mathfrak{z}}$, which after instantiation at the
deconditioning step above becomes $x$. Therefore, when we decondition on
$R_{\mathfrak{z}}$ below, the integral will be from $0$ to $x$.

We will apply the last expectation, $\mathbb{E}_{R_{\mathfrak{z}}}$, using the
property that the expectation of an integral equals the integral of the
expectation of the integrand. The integral will extend from $R=r$ as a minimum
to $\infty$. Note that the distribution of $R_{\mathfrak{z}}$ is Rayleigh, as
it is an association distance (based on the closest distance). In light of
(19), using the Rayleigh distribution $f_{R_{\mathfrak{z}}}(y)=2\lambda\pi
ye^{-\lambda\pi y^{2}}$, we have

$\displaystyle\mathcal{L}_{I_{u}}(s)$ $\displaystyle=$ $\displaystyle
e^{-2\pi\lambda\int_{0}^{\infty}\left(\mathbb{E}_{R_{\mathfrak{z}}}\left[\frac{1}{1+\left(sp\right)^{-1}R_{\mathfrak{z}}^{-\alpha\epsilon}x^{\alpha}}\right]\right)xdx}$
(20) $\displaystyle=$ $\displaystyle
e^{-2\pi^{2}\lambda^{2}\int_{0}^{\infty}x\int_{0}^{x^{2}}\frac{e^{-\lambda\pi
u}}{1+\left(sp\right)^{-1}u^{-\alpha\epsilon/2}x^{\alpha}}dudx}$

where $u=y^{2}$. Recall that $D_{z}^{-\alpha}$ and $R_{\mathfrak{z}}$ are
dependent in that if $D_{z}=x$ then $R_{\mathfrak{z}}<x$. That is, the
distance $R_{\mathfrak{z}}$ between an interferer and the typical BS is upper
bounded by the distance $D_{z}$ between the (relocated) interfering UE and the
typical BS. From (14) and (16), it follows that

$\displaystyle\mathcal{L}_{I_{u}}(\xi p^{-1}r^{\alpha(1-\epsilon)})$
$\displaystyle=$ $\displaystyle
e^{-2\pi^{2}\lambda^{2}\int_{0}^{\infty}x\int_{0}^{x^{2}}\frac{e^{-\lambda\pi
u}}{1+\left(sp\right)^{-1}u^{-\alpha\epsilon/2}x^{\alpha}}dudx}$ (21)
$\displaystyle=$ $\displaystyle
e^{-2\pi^{2}\lambda^{2}\int_{0}^{\infty}x\int_{0}^{x^{2}}\frac{\xi
r^{\alpha(1-\epsilon)}e^{-\lambda\pi u}}{\xi
r^{\alpha(1-\epsilon)}+u^{-\alpha\epsilon/2}x^{\alpha}}dudx}$

From (16) and (21), we get

$p_{u}=2\widetilde{\lambda}\int_{0}^{\infty}re^{-\widetilde{\lambda}r^{2}}e^{-2\widetilde{\lambda}^{2}\xi
r^{2\kappa(1-\epsilon)}\int_{0}^{\infty}x\int_{0}^{x^{2}}\frac{e^{-\widetilde{\lambda}u}}{\xi
r^{2\kappa(1-\epsilon)}+u^{-\epsilon\kappa}x^{2\kappa}}dudx}dr$ (22)

where $\widetilde{\lambda}=\pi\lambda$ and $\kappa=\alpha/2$.

From (22), it _appears_ that the UL coverage probability $p_{u}$ is dependent
on the BS density $\lambda$. Indeed, this belief is popular research work on
stochastic geometric models of the cellular UL system (See, for example, [3,
5, 8, 14, 25, 26, 31, 29, 32, 30, 33, 34, 35]). However, we will show in
Theorem 2 below that this belief is flawed. Specifically, we will show that
the presence of $\lambda$ in (22), as in a large number of similar expressions
in the literature of stochastic geometric modelling of cellular networks, is
superfluous.

:

Theorem 2: Under the stochastic geometric model of the cellular UL system, the
UL coverage probability $p_{u}$ is independent of the BS density $\lambda$.

Proof: The proof is attained through a sequence of changes of variables.
Starting with (22), use the substitution $v=r^{2}$ to get

$p_{u}=\widetilde{\lambda}\int_{0}^{\infty}e^{-\widetilde{\lambda}v}e^{-2\widetilde{\lambda}^{2}\xi
v^{\kappa(1-\epsilon)}\int_{0}^{\infty}x\int_{0}^{x^{2}}\frac{e^{-\widetilde{\lambda}u}}{\xi
v^{\kappa(1-\epsilon)}+u^{-\epsilon\kappa}x^{2\kappa}}dudx}dv$

Use $y=x^{2}$ to get

$p_{u}=\widetilde{\lambda}\int_{0}^{\infty}e^{-\widetilde{\lambda}v}e^{-\widetilde{\lambda}^{2}\xi
v^{\kappa(1-\epsilon)}\int_{0}^{\infty}\int_{0}^{y}\frac{e^{-\widetilde{\lambda}u}}{\xi
v^{\kappa(1-\epsilon)}+u^{-\epsilon\kappa}y^{\kappa}}dudy}dv$

Use $x=\widetilde{\lambda}u$ to get

$p_{u}=\widetilde{\lambda}\int_{0}^{\infty}e^{-\widetilde{\lambda}v}e^{-\widetilde{\lambda}\xi
v^{\kappa(1-\epsilon)}\int_{0}^{\infty}\int_{0}^{\widetilde{\lambda}y}\frac{e^{-x}}{\xi
v^{\kappa(1-\epsilon)}+\left(\frac{x}{\widetilde{\lambda}}\right)^{-\epsilon\kappa}y^{\kappa}}dxdy}dv$

Use $z=\widetilde{\lambda}v$ to get

$p_{u}=\int_{0}^{\infty}e^{-z}e^{-\widetilde{\lambda}\xi
z^{\kappa(1-\epsilon)}\int_{0}^{\infty}\int_{0}^{\widetilde{\lambda}y}\frac{e^{-x}}{\xi
z^{\kappa(1-\epsilon)}+x^{-\epsilon\kappa}\left(\widetilde{\lambda}y\right)^{\kappa}}dxdy}dz$

Finally, use $u=\widetilde{\lambda}y$ to get

$p_{u}=\int_{0}^{\infty}e^{-z\left(1+\xi
z^{\kappa(1-\epsilon)-1}\int_{0}^{\infty}\int_{0}^{u}\frac{e^{-x}}{\xi
z^{\kappa(1-\epsilon)}+x^{-\epsilon\kappa}u^{\kappa}}dxdu\right)}dz$ (23)

where $\lambda$ has totally disappeared, proving the theorem. $\blacksquare$

This is the final expression for the UL coverage probability. The striking
observation about the above $p_{u}$ is that for all $\alpha$, the UL coverage
probability is independent of the density $\lambda$ of the BSs. A possible
interpretation is that as the BS density gets higher, the Voronoi cells get
smaller, and the UEs get closer to their serving BSs, and vice versa. This
means that the increase in interference due to the increase of BS density,
will be offset by an increase in signal due to the nearness of the tagged UE
to the typical BS, keeping SIR, hence the coverage probability, unchanged.

Just like Theorem 1, Theorem 2 calls for a revisit to the UL coverage
probability expressions (and their derivatives) that include $\lambda$. Those
expressions exist in a huge number of articles, such as those cited just above
Theorem 2. It is worth examining them for possible elimination of a parameter
that may exist misleadingly.

UL Special cases:

First, we will consider the special case of $\epsilon=0$, i.e. there is no FPC
by the UE. For this case, we get from (23) that

$\displaystyle p_{u}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}e^{-z\left(1+\xi
z^{\kappa-1}\int_{0}^{\infty}\int_{0}^{u}\frac{e^{-x}}{\xi
z^{\kappa}+u^{\kappa}}dxdu\right)}dz$ (24) $\displaystyle=$
$\displaystyle\int_{0}^{\infty}e^{-z\left(1+\xi
z^{\kappa-1}\int_{0}^{\infty}\frac{1-e^{-u}}{\xi
z^{\kappa}+u^{\kappa}}du\right)}dz$

Using a computer algebra system (e.g., Maxima (R)) or a table of integrals
(e.g., [36]: p. 325, #3.241], one can find

$\int_{0}^{\infty}\frac{1}{\xi
z^{\kappa}+u^{\kappa}}du=\frac{1}{\kappa\xi^{1-1/\kappa}z^{\kappa-1}}B\left(1-\frac{1}{\kappa},\frac{1}{\kappa}\right),\qquad
0<\frac{1}{\kappa}<1$

where $B(x,y)$ is the Bessel function of $x$ and $y$. Substituting in (24), we
get

$p_{u}=\int_{0}^{\infty}e^{-z\left(1+\left(\frac{\sqrt[\kappa]{\xi}}{\kappa}B\left(1-\frac{1}{\kappa},\frac{1}{\kappa}\right)-\xi
z^{\kappa-1}\int_{0}^{\infty}\frac{e^{-u}}{\xi
z^{\kappa}+u^{\kappa}}du\right)\right)}dz$

For $\alpha=4$ ($\kappa=2$), one can find

$\int_{0}^{\infty}\frac{1}{\xi
y^{2}+x^{2}}dx=\frac{1}{2y\sqrt{\xi}}B\left(0.5,0.5\right)=\frac{\pi}{2y\sqrt{\xi}}$

Substituting in (24), we get

$p_{u}=\int_{0}^{\infty}e^{-z\left(1+\left(\frac{\sqrt{\xi}}{2}B\left(\frac{1}{2},\frac{1}{2}\right)-\xi
z\int_{0}^{\infty}\frac{e^{-u}}{\xi z^{2}+u^{2}}du\right)\right)}dz$ (25)

Using a computer algebra system (e.g., Maxima (R)) or a table of integrals
(e.g., [36]: p. 343, #3.354], one can find the integral at the exponent

$\int_{0}^{\infty}\frac{e^{-u}}{\xi
z^{2}+u^{2}}dx=\frac{1}{\sqrt{\xi}z}\left[\mathop{\rm
ci}(\sqrt{\xi}z)\sin(\sqrt{\xi}z)-\mathop{\rm
si}(\sqrt{\xi}z)\cos(\sqrt{\xi}z)\right]$ (26)

where

$\mathop{\rm si}(x)=\mathop{\rm
Si}(x)-\frac{\pi}{2}=-\int_{x}^{\infty}\frac{\sin t}{t}dt$

and

$\mathop{\rm ci}(x)=\mathop{\rm Ci}(x)=-\int_{x}^{\infty}\frac{\cos t}{t}dt$

are the sine and cosine integrals, respectively. Using (26) and the fact that
$B\left(\frac{1}{2},\frac{1}{2}\right)=\pi$ in (25), the coverage probability
for $\alpha=4$ becomes

$p_{u}=\int_{0}^{\infty}e^{-z\left(1+\sqrt{\xi}\left(\frac{\pi}{2}-\mathop{\rm
ci}(\sqrt{\xi}z)\sin(\sqrt{\xi}z)+\mathop{\rm
si}(\sqrt{\xi}z)\cos(\sqrt{\xi}z)\right)\right)}dz$ (27)

Second, we will consider the special case of $\epsilon=1$, i.e. there is
channel inversion by the UE. For this case, we get from (23) that

$\displaystyle p_{u}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}e^{-z\left(1+\xi
z^{-1}\int_{0}^{\infty}\int_{0}^{u}\frac{x^{\kappa}e^{-x}}{\xi
x^{\kappa}+u^{\kappa}}dxdu\right)}dz$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}e^{-z}e^{-\xi\int_{0}^{\infty}\int_{0}^{u}\frac{x^{\kappa}e^{-x}}{\xi
x^{\kappa}+u^{\kappa}}dxdu}dz$ $\displaystyle=$
$\displaystyle\exp\left(-\xi\int_{0}^{\infty}\int_{0}^{u}\frac{x^{\kappa}e^{-x}}{\xi
x^{\kappa}+u^{\kappa}}dxdu\right)$

At $\kappa=2$:

$p_{u}=\exp\left(-\xi\int_{0}^{\infty}\int_{0}^{u}\frac{x^{2}e^{-x}}{\xi
x^{2}+u^{2}}dxdu\right)$

which could not be simplified further because the inner integral could not be
evaluated.

##  4\. Simulation

Despite diligent search in the literature, we have not spotted any formal
stochastic geometric simulation algorithms. What is available is only code,
and even that is rare, poorly documented, and sporadic. In this section, we
provide the simulation algorithms we have developed for the present article.
Although the algorithms were primarily intended to validate our analytical
results, we feel they can be of value in similar stochastic geometric
research. Since DL and UL are different, as we have seen in the modelling,
their simulation algorithms are different, and therefore will be introduced
separately. However, we will give below some points common to both.

Although we have proved that the BS density $\lambda$ has no bearing on the
coverage probability in either DL or UL, we have included it in the two
simulation algorithms. The reason for that is to allow the interested reader
to verify practically what we have proven theoretically. Indeed, one can use
the simulation algorithms below to ascertain that for any value of $\lambda$,
the coverage probability remains the same for the same path loss exponent
$\alpha$ and SIR threshold $\xi$. Actually, it was experimentation of this
sort that led us to discover the irrelevance of $\lambda$ and motivated the
theoretical proof presented earlier.

###  4.1. Downlink

The modelling of the DL mode was undoubtedly simple. Its simulation is just as
simple. After scattering the BSs of a realization in the simulation window, we
select the closest BS to the origin and make it the tagged BS, considering
that the typical UE is at the origin already. That is, we do not need to
generate any UEs. The transmission received at the origin from the tagged BS
will be then considered signal and every other transmission received there
considered interference. It is a simple matter then to calculate the SIR at
the typical UE and use it to determine if the typical UE is covered by the
tagged BS for this realization, for the given $\xi$. The coverage probability
is evaluated after all $\mathcal{N}$ realizations have been generated, by
dividing the number of realizations where the typical UE was covered by the
total number $\mathcal{N}$ of realizations. The algorithm then proceeds to
calculate the SIR at the typical UE, using (4).

###  4.2. Uplink

UL simulation is more involved than DL simulation, as was the case with their
modelling. A key element for a successful UL simulation is the association
table, $A$, an example of which is shown in Table 2. For a network of $N$ BSs,
this $N\times 2$ table has one column for the BSs and another for the UEs. In
particular, the table stores the coordinates of each BS and those of its
associated UE (or conversely, each UE and its serving BS) in one row. As such,
the table not only tells which UE is associated to which BS, but also allows
calculating the distance between any BS or UE and any other BS or UE. These
distances are crucial to calculate the interference at a given receiver. For
example, referring to Table 2, the first BS is at the origin and its
associated UE is $734.2$ m west and $628.4$ m north, i.e. the BS and UE are
$\sqrt{734.2^{2}+628.4^{2}}=966.4$ m apart. We note that the third BS has no
associated UE yet.

Table 2. Example BS-UE Association table $A$ for UL simulation. BS | UE
---|---
(0,0) | (-734.2,628.4)
(1243.2,-221.4) | (973.2,1628.4)
(-345.2,928.4) | NULL
$\ldots$ | $\ldots$

The association table $A$ is built by the simulation algorithm as follows.
First, the BS column is filled at once by generating $N$ pairs of random
numbers (drawn from a uniform distribution) in the interval $[-S,S]$, where
$S$ is the simulation window (square) side. The first number in the pair is
taken as the $x$ coordinate, and the second as the $y$ coordinate. Then the
UEs are added incrementally, from top to bottom, one at a time based on the
BS-UE association rule, which says: a UE must be closer to its serving BS than
to any other BS in the cellular network. Adding the first UE to the table is
simplest. A pair of random numbers is first generated (drawn from a uniform
distribution) in the interval $[-S,S]$ to determine the location of the UE in
the simulation window. Then the distances between this UE and all the BSs in
the network, whose locations are in column $1$ of the table, are calculated.
The UE is assigned immediately to the closest BS.

1

2Input: $\mathcal{N}$, $\xi$, $p$, $\lambda$, $S$, $\alpha$

3 Output: $p_{d}$

4 Covered := 0

5 for  _$i:=1$ to $\mathcal{N}$_ do

/* Find number $N$ of BSs to deploy in simulation window: */

6 Generate a Possion distributed random number $N$, $N\sim$ Pois($\lambda
S^{2}$)

7

/* Scatter $N$ BSs over an $S\times S$ square window, storing the location at
a 1 column array $D$: */

8

9 for  _$i:=1$ to $N$_ do

10

11 Generate two uniformly distributed random numbers $(x,y)$, where $x,y\sim$
U($-S,S$), for BS $i$ location.

12 $D(i)=(x,y)$

13 end for

14

15 Identify BS $n$ nearest the origin and designate it tagged BS

16

/* Find SIR at typical UE and check if the latter is covered by tagged BS for
the given $\xi$: */

17

18 Generate a Poisson random number $G$, where $G\sim$ Pois(1)

19 Calculate distance $R$ between the tagged BS and origin

20 NUM $=pGR^{-\alpha}$

21 DENUM $=0$

22 for  _$i:=1$ to $N$_ do

23 if  _$i\neq n$_ then

24 Generate a Poisson random number $G_{i}$, where $G_{i}\sim$ Pois(1)

25 Use array $D$ to calculate the distance $D_{i}$ between BS $i$ and origin

26 DENUM $=$ DENUM $+pG_{i}{D_{i}}^{-\alpha}$

27

28 end if

29

30 end for

31 SIR := NUM/DENUM

32

33 if  _SIR $>\xi$_ then

34 Covered := Covered + 1

35 end if

36

37 end for

/* Calculate DL coverage probability $p_{d}$ for this threshold $\xi$: */

38 $p_{d}$ := Covered $/\mathcal{N}$

Algorithm 1 Simulation of cellular DL system

1

2Input: $\mathcal{N}$, $\xi$, $p$, $\lambda$, $S$, $\alpha$, $\epsilon$

3 Output: $p_{u}$

4 Covered := 0

5 for  _$i:=1$ to $\mathcal{N}$_ do

/* Find number $N$ of BSs to deploy in simulation window: */

6 Generate a Possion distributed random number $N$, $N\sim$ Pois($\lambda
S^{2}$)

7

/* Build table $A$ for a network of $N$ BSs and $N$ UEs: */

8

9 for  _$i:=1$ to $N$_ do

$A(i,1)=A(i,2)=$ NULL /* Empty the association table, $A$ */

10

11 end for

12

/* 1- Fill in column 1 of $A$ with BSs, making BS 1 typical. */

13

14 $A(1,1)=(0,0)$

15 for  _$i:=2$ to $N$_ do

16

17 Generate random location $(x,y)$, where $x,y\sim$ U($-S,S$), for BS $i$.

18 $A(i,1)=(x,y)$

19 end for

20

/* 2- Fill in column 2 of table $A$ with UEs. */

21

22 while  _( $\exists$ NULL $\in A$)_ do

23 Generate random location $(u,v)$, where $u,v\sim$ U($-S,S$), for a UE.

24

25 Calculate the $N$ distances between the UE and every BS $i$

26 Identify BS $n$, the BS nearest the UE

27 if  _$A(n,2)=$ NULL_ then

$A(n,2)=(u,v)$ /* Assocate UE with BS $n$ if available. */

28

29 end if

30

31 end while

32

/* Find SIR at typical BS, update ‘Covered’ & calculate $p_{u}$: */

33

34 Generate a Poisson random number $G$, where $G\sim$ Pois(1)

35 Calculate from table $A$ distance $R$ between typical BS and tagged UE

36 NUM $=pGR^{-\alpha(1-\epsilon)}$

37 DENUM $=0$

38 for  _$i:=2$ to $N$_ do

39 Generate a Poisson random number $G_{i}$, where $G_{i}\sim$ Pois(1)

40 Calculate from table $A$ distance $R_{i}$ between UE $i$ and serving BS

41 Calc. from table $A$ distance $U_{i}$ between UE $i$ and typical BS

42 DENUM $=$ DENUM $+pG_{i}{R_{i}}^{\alpha\epsilon}{U_{i}}^{-\alpha}$

43

44 end for

45 SIR := NUM/DENUM

46

47 if  _SIR $>\xi$_ then

48 Covered := Covered + 1

49 end if

50

51 end for

52$p_{u}$ := Covered $/\mathcal{N}$

Algorithm 2 Simulation of cellular UL system

The addition of the second UE gets a little harder. First, we generate the
location and calculate the distance as we did for the first UE. However, now
there is a possibility that the second UE is closest to the BS to which the
first UE was associated. Here the attempt fails. We keep making attempts till
we find a nearest BS that is available, to which we make the association by
recording the coordinates of the UE in the same row of that nearest BS. Of
course, adding UEs keeps getting harder and harder as the number of
associations increases since the chance of finding a BS that is both closest
to the UE and at the same time available becomes increasingly small. For this
probabilistic behavior, building the association table is the most time
consuming part of the algorithm. All the other parts are just fast
computations based on distances obtained from the numbers of this table.

The algorithm then proceeds to calculate the SIR at the typical BS, using
(15). Note that in the simulation we can calculate exactly the two distances:
$R_{\mathfrak{z}}$, the distance between the interfering UE and its serving
BS, and $U_{\mathfrak{z}}$, the distance between the interfering UE and the
typical BS. That is, we can do in an exact way in the simulation, what we did
in an approximate way in the modelling. Interestingly, however, the results of
both the modelling and simulation conform impressively, testifying to the
validity of the approximation made in the modelling.

Assessing the time complexity of both algorithms is difficult since the number
$N$ of BSs, i.e. PPP points, generated in each realization is not constant,
but a Poisson distributed RV ranging from $0$ to $\infty$. This difficulty
increases even more in the UL case, where the way the association table $A$ is
built is also random, as explained earlier.

##  5\. Experimental Work

We implemented the simulation algorithms in the Matlab (R) language as two
separate programs, one for DL and the other for UL. We coded the algorithms in
Matlab (R). The simulation time took around 1 hr on a Laptop of $8$ GB RAM and
Core $5$ CPU at $2.4$ GHz, for $3000$ runs, on a simulation window in the form
of a square of side $2000$ m. Matlab (R) was also used to perform the
computations of the analytical results, and proved particularly efficient with
numerical integration of formidable functions.

An important implementation tip is in order. The pseudo code shown in the two
simulation algorithms calculates the coverage probability for only one value
of the threshold, $\xi$. To sketch a smooth curve over a reasonable range of
the threshold, such as that of the coverage probability figures below, where
the range is from $-15$ dB to $15$ dB, we should then run the algorithm
$\mathcal{N}$ times, with $\mathcal{N}$ typically exceeding $3000$ for good
convergence, for each of the $31$ dB values. To save simulation time, however,
we did something interesting that yielded the same results nonetheless. In
each realization, we tested the resulting SIR with all the values of $\xi$ in
the range, starting from $-15$ dB and going upwards, incrementing the variable
‘Covered’ each time the former exceeds the latter for this particular value of
$\xi$, for which an array, rather than a single variable, is created. Once the
former stops exceeding the latter, we quit the comparisons for this
realization. This trick reduces the simulation time significantly, yet
produces the same results. The same thing, by the way, applies in UL when
using more than one value of the power control factor, $\epsilon$. Instead of
running $\mathcal{N}$ realizations for each $\epsilon$, we use the same
realization to get results for as many $\epsilon$ values as desired. Of
course, one would use arrays for these multiple results, instead of the
variables in the algorithms below which are intended for only one $\xi$ value
and one $\epsilon$ value (in case of UL.) We always calculate the SIR at the
typical receiver placed at the origin (which is a UE in DL and a BS in UL.)

From our experience, there are two tricks that can make the running of the
algorithms extremely fast. First, we use the same realization to check for all
required thresholds. For example, in our experiments we use 31 thresholds,
from $-15$ dB to $15$ dB, in increments of $1$ dB. Then, when we obtain the
SIR for a realization, we keep checking (using a loop) it against the
thresholds from the smaller going towards the bigger, within the same
realization. Once the SIR fails to exceed a threshold, we exit (or break) the
loop, since there is no point in checking the SIR against greater thresholds.
Moreover, in the UL system in particular, we use also the same realization for
all the values of the power control factor $\epsilon$ under consideration. In
some implementations we have seen, there is one realization per threshold
value per control factor value, totalling $93$ realizations if there are $31$
threshold values and $3$ power factor values. These $93$ realizations are
replaced by only $1$ realization, a huge savings in terms of simulation time.

Furthermore, in our UL algorithm we used a trick in building the association
table that accelerated the simulation immensely. Unlike popular
implementations which insist on building the association table strictly from
top down, starting by associating the first BS, then the second BS, and so on,
we associate with the BS that turns out to be closest to the UE just
generated, provided the BS is available. We note in passing that we start in
UL by placing one BS at the origin, designating it to be the typical BS. This
BS can in fact be at any row of the association table $A$, so it might as well
be at the first row, for easy reference.

$-15$$-10$$-5$$0$$5$$10$$15$$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$$\alpha=6$
($\kappa=3$)$\alpha=4$ ($\kappa=2$)SIR threshold, $\xi$ (dB)DL coverage
probability, $p_{d}$ AnalayticSimulation Figure 3. DL coverage probability,
$p_{d}$, for two values of $\alpha$
$-15$$-10$$-5$$0$$5$$10$$15$$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$$\alpha=6$
($\kappa=3$)$\alpha=4$ ($\kappa=2$)SIR threshold, $\xi$ (dB)UL coverage
probability, $p_{u}$ AnalayticSimulation Figure 4. UL coverage probability,
$p_{u}$, for two values of $\alpha$, and with no FPC, i.e. $\epsilon=0$.

We sketch in Figures 3 and 4 the analytical and simulation values of the DL
and UL coverage probabilities, respectively, for a range of a threshold,
$\xi$, extending from $-15$ dB to $15$ dB. For Figure 3, we used (12) and (13)
to sketch the analytical results for the cases of $\alpha=4$ and $\alpha=6$,
respectively, and Aglorithm 1 to sketch the simulation results. On the other
hand, for Figure 4, we used (23) and (27) to sketch the analytical results for
the cases of $\alpha=6$ and $\alpha=4$, respectively, and Algorithm 2 to
sketch the simulation results. The UL coverage probability is sketched for the
case of no FPC, i.e. $\epsilon=0$.

As can be seen the agreement between the two types of values is excellent.
From the actual values we used for graphig the curves, the difference between
the analytical and simulation results is within only $1\%$, out of as few as
$3000$ simulation runs (realizations) on a $2000\times 2000$ m simulation
window.

##  6\. Conclusions

In this article we have investigated thoroughly the modelling and simulation
of cellular DL and UL cellular channels using stochastic geometry. A number of
contributions have been made in the process. First, we introduced the DL and
UL models concisely together so as to expose their similarities and
differences vividly. In addition, we were able to derive closed form
expressions for the DL coverage probability in two special cases. Then, for
UL, known to be challenging, we introduced an approximation that circumvented
the challenge and yet yielded excellent results validated later by simulation.
For UL also, we obtained for the coverage probability an expression, though
not closed form, but simple enough to calculate easily. For the simulation
dimension of our study, we presented two effecient Monte Carlo simulation
algorithms designed to validate the models, but can be useful to anyone using
stochastic geometric in the communications field or otherwise. Finally, we
have proved two theorems that go against established belief in the cellular
modelling community. Namely, we have proved that under the stochastic
geometric model, the coverage probability in either DL and UL is independent
of the BS density. This finding calls for a revisit to a large body of results
published in the past decade with the BS density present superfluously in
them.

## References

  * [1] F. Baccelli, B. Blaszczyszyn, Stochastic geometry and wireless networks: Volume II Applications, Foundations and Trends in Networking, Boston, 2010. https://doi.org/10.1561/1300000026.
  * [2] M. Haenggi, Stochastic Geometry for Wireless Networks, Cambridge University Press, Cambridge, 2012. https://doi.org/10.1017/CBO9781139043816.
  * [3] H. ElSawy, A. Sultan, M. Alouini, M. Win, Modeling and Analysis of Cellular Networks Using Stochastic Geometry: A Tutorial, IEEE Communications Surveys and Tutorials, 19 (1) (2017) pp. 167-203. https://doi.org/10.1109/COMST.2016.2624939.
  * [4] S. Okegbile, B. Maharaj, A. Alfa, Stochastic geometry approach towards interference management and control in cognitive radio network: A survey, Computer Communications, 166 (15) (2021) pp. 174-195.
  * [5] C. Liu, Y. Shen, C. Lee, Energy-Efficient Activation and Uplink Transmission for Cellular IoT, IEEE Internet of Things Journal, 7 (2) (2020) pp. 906-921. https://doi.org/10.1109/JIOT.2019.2946331.
  * [6] X. Tang, X. Xu, M. Haenggi, Meta Distribution of the SIR in Moving Networks, IEEE Transactions on Communications, 68 (6) (2020) pp. 3614-3626. https://doi.org/10.1109/TCOMM.2020.2975996.
  * [7] Q. Liu, Z. Zhang, The analysis of coverage probability, ASE and EE in heterogeneous ultra-dense networks with power control, Digital Communications and Networks, 6 (4) (2020) pp. 524-533.
  * [8] M. Haroon, F. Muhammad, Z. Abbas, G. Abbas, N. Ahmed, S. Kim, Proactive Uplink Interference Management for Nonuniform Heterogeneous Cellular Networks, IEEE Access, 8 (2020) pp. 55501-55512. https://doi.org/10.1109/ACCESS.2020.2981631.
  * [9] N. Kouzayha, H. Elsawy, H. Dahrouj, T. Al-Naffouri, Meta Distribution of Downlink SIR for Binomial Point Processes, IEEE Wireless Communications Letters, 10 (7) (2021) pp. 1557-1561. https://doi.org/10.1109/LWC.2021.3074399.
  * [10] B. Blaszczyszyn, M. Haenggi, P. Keeler, S. Mukherjee, Stochastic Geometry Analysis of Cellular Networks, Cambridge University Press, Cambridge, 2018. https://doi.org/10.1017/9781316677339.
  * [11] C. Lee, C. Shih, Y. Chen, Stochastic geometry based models for modeling cellular networks in urban areas, Wireless networks, 19 (6) (2013) pp. 1063-1072. https://doi.org/10.1007/s11276-012-0518-0.
  * [12] J. Andrews, F. Baccelli, R. Ganti, A Tractable Approach to Coverage and Rate in Cellular Networks, IEEE Transactions on Communications, 59 (11) (2011) pp. 3122-3134. https://doi.org/10.1109/TCOMM.2011.100411.100541.
  * [13] T. Bai, R. Heath, Coverage and Rate Analysis for Millimeter-Wave Cellular Networks, IEEE Transactions on Wireless Communications, 14 (2) (2015) pp. 1100-1114. https://doi.org/10.1109/TWC.2014.2364267.
  * [14] P. Herath, C. Tellambura, W. Krzymień, Coverage probability analysis of three uplink power control schemes: Stochastic geometry approach. EURASIP journal on wireless communications and networking, 1 (2018) pp. 1-14.
  * [15] E. Gilbert, Random plane networks, Journal of the Society for Industrial and Applied Mathematics, 9 (4) (1961) pp. 533-543. https://doi.org/10.1137/0109045.
  * [16] K. Li, Downlink data transmission scheduling algorithms in wireless networks, Simulation Modelling Practice and Theory 19 (2011) pp. 1427-1444. https://doi.org/10.1016/j.simpat.2010.11.010.
  * [17] J. Chen, C. Yuan, Coverage probability and average rate of downlink user-centric wireless cellular networks with composite $\kappa$-$\mu$ shadowed and lognormal shadowed fading, IET Communications, 13 (17) (2019) pp. 2805-2813. https://doi.org/10.1049/iet-com.2019.0220.
  * [18] Q. Liu, J. Baudais, P. Mary, A Tractable Coverage Analysis in Dynamic Downlink Cellular Networks, 2020 IEEE 21st International Workshop on Signal Processing Advances in Wireless Communications (SPAWC), 2020, pp. 1-5, https://doi.org/10.1109/SPAWC48557.2020.9154321.
  * [19] M. Ouamri, M. Oteşteanu, A. Isar, M. Azni, Coverage, Handoff and cost optimization for 5G Heterogeneous Network, Physical Communication, 39 (5) (2020) 101037. https://doi.org/10.1016/j.phycom.2020.101037.
  * [20] M. Fadoul, Modeling multi-tier heterogeneous small cell networks: rate and coverage performance, Telecommunication Systems, 75 (2020) pp. 369-382. https://doi.org/10.1007/s11235-020-00680-y.
  * [21] J. Lei, H. Chen, F. Zhao, Stochastic geometry analysis of downlink spectral and energy efficiency in ultradense heterogeneous cellular networks, Mobile Information Systems, 2018 (2018) 1684128. https://doi.org/10.1155/2018/1684128.
  * [22] M. Fadoul, Rate and coverage analysis in multi-tier heterogeneous network using stochastic geometry approach, Ad Hoc Networks, 98 (2020) 102038. https//doi.org/10.1016/j.adhoc.2019.102038.
  * [23] C. Bouras, R. Kalogeropoulos, A QoS driven adaptive mechanism for downlink and uplink decoupling in 5G, Internet of Things, 11 (2020) 100217. https://doi.org/10.1016/j.iot.2020.100217.
  * [24] M. Arif, S. Wyne, J. Ahmed, Efficiency analysis of a K-tier clustered HCN using dual connectivity with DUDe access, AEU - International Journal of Electronics and Communications, 123 (2020) 153291. https://doi.org/10.1016/j.aeue.2020.153291.
  * [25] N. Kouzayha, Z. Dawy, J. Andrews, H. ElSawy, Joint Downlink/Uplink RF Wake-Up Solution for IoT Over Cellular Networks, IEEE Transactions on Wireless Communications, 17 (3) (2018) pp. 1574-1588, https://doi.org/10.1109/TWC.2017.2781696.
  * [26] H. Mariam, I. Ahmed, M. Aslam, Coverage probability of uplink millimeter wave cellular network with non-homogeneous interferers’ point process, Physical Communication, 45 (2021) 101274. https://doi.org/10.1016/j.phycom.2021.101274.
  * [27] S. Kamiya, K. Yamamoto, S. Kim, T. Nishio, M. Morikura, SINR Distribution and Scheduling Gain Analysis of Uplink Channel-Adaptive Scheduling, IEEE Transactions on Wireless Communications, 19 (4) (2020) pp. 2321-2335. https://doi.org/10.1109/twc.2019.2963866.
  * [28] O. Onireti, L. Zhang, A. Imran, M. Imran, Outage Probability in the Uplink of Multitier Millimeter Wave Cellular Networks, IEEE Systems Journal, 14 (2) (2020) pp. 2520-2531. https://doi.org/10.1109/jsyst.2020.2965001.
  * [29] X. Jia, Q. Fan, W. Xu, L. Yang, Cross-Tier Dual-Connectivity Designs of Three-Tier Hetnets With Decoupled Uplink/Downlink and Global Coverage Performance Evaluation, IEEE Access, 7 (2019) pp. 16816-16836. https://doi.org/10.1109/ACCESS.2019.2895389.
  * [30] S. Ali, M. Aslam, I. Ahmed, Uplink coverage probability and spectral efficiency for downlink uplink decoupled dense heterogeneous cellular network using multi-slope path loss model, Telecommunication systems, 72 (4) (2019) pp. 505-516.
  * [31] J. Andrews, A. Gupta, H. Dhillon, A primer on cellular network analysis using stochastic geometry, ArXiv: 1604.03183, 2016, https://arxiv.org/abs/1604.03183.
  * [32] Y. Gao, S. Yang, S. Wu, M. Wang, X. Song, Coverage Probability Analysis for mmWave Communication Network With ABSF-Based Interference Management by Stochastic Geometry, IEEE Access, 7 (2019) pp. 133572-133582. https://doi.org/10.1109/ACCESS.2019.2940537.
  * [33] A. Kundu, R. Pal, M. Kumar, T. Sreejith, Uplink and Downlink Performance Bounds for Full Duplex Cellular Networks, 2020 IEEE International Black Sea Conference on Communications and Networking (BlackSeaCom), Odessa, Ukraine, 26-29 May 2020, pp. 1-6. https://doi.org10.1109/BlackSeaCom48709.2020.9235024.
  * [34] E. Sadeghabadi, S. Abarghouyi, B. Makki, M. Kenari, T. Svensson, Asynchronous Downlink Massive MIMO Networks: A Stochastic Geometry Approach, IEEE Transactions on Wireless Communications, 19 (1) (2020) pp. 579-594, https:/doi.org/10.1109/TWC.2019.2946824.
  * [35] X. Wang, M. Gursoy, Coverage Analysis for Energy-Harvesting UAV-Assisted mmWave Cellular Networks, IEEE Journal on Selected Areas in Communications, 37 (12) (2019) pp. 2832-2850. https://doi.org/10.1109/JSAC.2019.2947929.
  * [36] I. Gradshteyn, I. Ryzhik, Table of Integrals, Series, and Products, 8th edn., ed. by D. Zwillinger, Elsevier Academic Press, Waltham, 2015. https://doi.org/10.1016/C2010-0-64839-5.