# Fair Throughput Optimization with
a Dynamic Network of Drone Relays⋆††thanks: ⋆A preliminary 6-page version of
this manuscript appeared in the proceedings of IEEE INFOCOM WKSHPS – MiSARN
2019 [1].

Edgar Arribas, Vincenzo Mancuso, Vicent Cholvi  E. Arribas is with Universidad
CEU San Pablo, Madrid, Spain. E-mail<EMAIL_ADDRESS>V. Mancuso is
with IMDEA Networks Institute, Madrid, Spain. E-mail:
<EMAIL_ADDRESS>V. Cholvi is with Universitat Jaume I (UJI),
Castelló, Spain. E-mail<EMAIL_ADDRESS>The work of E. Arribas is partially
supported by the FPU15/02051 grant from the Spanish Ministry of Education,
Culture and Sports (MECD). The work of V. Mancuso is supported by the
RYC-2014-16285 grant from the Spanish Ministry of Economy and Competitiveness
and by the Spanish Ministry of Science and Innovation grant
PID2019-109805RB-I00 (ECID).

###### Abstract

Aiding the ground cellular network with aerial base stations carried by drones
has experienced an intensive raise of interest in the past years.
Reconfigurable air-to-ground channels enable aerial stations to enhance users
access links by means of seeking good line-of-sight connectivity while
hovering in the air. In this paper, we propose an analytical framework for the
3D placement of a fleet of coordinated drone relays. This framework optimizes
network performance in terms of user throughput fairness, expressed through
the $\alpha$-fairness metric. The optimization problem is formulated as a
mixed-integer non-convex program, which is intractable. Hence, we propose an
extremal-optimization-based algorithm, Parallelized Alpha-fair Drone
Deployment, that solves the problem online, in low-degree polynomial time. We
evaluate our proposal by means of numerical simulations over the real topology
of a dense city. We discuss the advantages of integrating drone relay stations
in current networks and test several resource scheduling approaches in both
static and dynamic scenarios, including with progressively larger and denser
crowds.

###### Index Terms:

Aerial networks; Relay; UAV; Mobile networks; Optimization; $\alpha$-Fairness.

##  1 Introduction

The infrastructure of cellular networks is evolving towards flexible and
reconfigurable solutions, able to cope with the highly variable densities of
users. Specifically, the current generation of cellular communications, namely
5G [2], embeds new transmission techniques as well as novel communication
paradigms, including smart and flexible relaying [3]. Besides, wireless
relaying with mobility of relays is possible thanks to precise beamforming and
highly efficient cooperative transmission techniques, which makes it possible
to operate broadband wireless backhaul links [4]. Without such mature
technological tools, many attempts toward mobile relaying have failed in the
past, since the advent of broadband wireless data networks [5].

It is therefore currently possible to mount mobile relays on, e.g., transport
vehicles and drones, which brings the possibility of moving the network with
the users and position relays where the fixed infrastructure cannot sustain
the user demand [6, 7]. However, there exist serious concerns on the
practicality of mobile relaying, due to interference management problems. For
instance, Guo and O’Farrel [8] have derived the capacity of OFDMA cellular
networks (like LTE-A, but also in 5G) in the presence of relays reusing
cellular frequencies, and showed that relays need to be operated onto
orthogonal frequencies. Besides, it is known that using over orthogonal
frequencies gives additional advantages in terms of simplified resource
allocation control [9]. Thus, we study the case of relay drones and base
stations (BSs) transmitting on orthogonal frequencies. Specifically, we tackle
the optimization of 3D positions for a fleet of coordinated drone relays
aiding a set of ground BSs, as depicted in Fig. 1, aiming at a fair throughput
distribution among users.

We base the optimization on the throughput as evaluated by means of an
$\alpha$-fairness metric [10], which is a high-level generalization of
fairness metrics. The parameter $\alpha$ can be tuned to analytically target,
e.g., maximum throughput, proportional fairness or max-min fairness with a
single framework. In the analysis, we account for the behavior of standard
3GPP networks and user access protocols, and model transmission technology
details like the random variations of signal quality received by users over
time, the interference caused by relays and BSs, the use of slotted time-
frequency resources, the wireless backhaul attachment, and cell selection and
resource allocation procedures. Specifically, we adopt stochastic models for
path-loss and availability of line-of-sight (LoS) and non-LoS (NLoS) links,
and cast our problem into an OFDMA-like resource allocation scheme with
several constraints.

Figure 1: Reference scenario: multi-drone-aided network with people moving
from a train station to a football stadium.

The problem of finding the exact optimal drone positions is NP-complete. In
particular, our analysis unveils that the role of interference caused by
drones and the stochastic characterization of LoS between drones and ground
users make the optimization problem intractable. However, is the problem
really intractable in practice? To answer this question we analyze the drone-
positioning problem under the perspective of a fair maximization of the
network throughput provided to mobile users, and unveil its structure, which
allows us to find efficient paths towards the computation of near-optimal
solutions with low time complexity. Moreover, we show that the most critical
part of the problem can be efficiently addressed by leveraging Extremal
Optimization (EO) algorithms, which are a class of algorithms specifically
designed for polynomial time optimization with intertwined variables [11]. The
EO operation is based on picking the “least fit” element of a finite set and
change its configuration parameters to improve a global utility function. EO
offers a lightweight solution to problems otherwise intractable, and for which
researchers proposed genetic algorithms. It is also more reliable than machine
learning approaches–e.g., reinforcement learning schemes–since it does not
require any training, which would instead be quite heavy and which does not
scale with the number of drones to control (the state space in which to train
the machine learning algorithm would quickly become unsustainable). We
therefore formulate a suitable utility function, targeting $\alpha$-fair user
throughput across the network, and design PADD: a Parallelized Alfa-fair Drone
Deployment Algorithm. By using an EO algorithm, PADD iteratively updates the
position of the least fit drone, i.e., the drone relay that contributes the
least to the utility function. We validate our algorithm and evaluate its
performance by means of simulations of realistic static and dynamic scenarios.
As an illustration of dynamic cases, we evaluate the performance of our
algorithm when customers move towards a stadium before a sport event, so that
their density grows over time. Beside illustrating the advantages offered by
PADD over state of the art algorithms, our numerical results show that
optimizing network throughput without considering fairness is not beneficial
in dense environments since drones serving the same area generate too much
interference.

The contributions of this paper are summarized:

  * •

We propose an analytical framework for the 3D placement of a fleet of
coordinated drone relays in cellular networks.

  * •

We propose an aerial cellular infrastructure that combines standard 3GPP
tools, beamforming techniques, wireless backhauling and fair resource
allocation.

  * •

We show that the analyzed problem is NP-Complete.

  * •

We formulate the 3D placement optimization problem in terms of the
$\alpha$-fairness metric to guarantee fair distribution of resources and
access link throughputs.

  * •

We propose PADD as an extremal-optimization-based algorithm to solve the
problem near-optimally.

  * •

We deterministically solve optimally in polynomial time the $\alpha$-fair
resource allocation problem for ground-plus-aerial cellular networks.

  * •

We assess and validate our proposal by means of comprehensive numerical
results.

The rest of the paper is structured as follows. Section 2 discusses related
work. Section 3 presents the system model, while Section 4 derives the
framework for optimizing drone positions under the $\alpha$-fairness metric.
Section 5 describes the design of our optimization algorithm. Section 6
provides numerical results. Section 7 discusses the findings of this paper and
possible practical implementation issues. Section 8 summarizes and concludes
the paper.

##  2 Related work

###  2.1 Relay Alternatives

Mobile and non-terrestrial relays have been investigated in several forms and
considering several technologies. For instance, the usage of satellite
networks [12] emerged several years ago. However, while satellites serve huge
areas, users can only achieve relatively low rates. Moreover, satellites
cannot adjust to the user’s topology, and service incurs high costs. In
contrast, drone relays may move dynamically at low altitudes and serve smaller
target regions on demand, where the ground network cannot sustain the high
demand from dense spots, so that a swarm of drones is able to rapidly act for
aerial connectivity assistance.

Drone relays are also different from fixed relays and Device-to-Device
(D2D)-based approaches [13, 14]. In fact, unlike those cases, drone
communications are neither fixed nor opportunistic, and the channel
propagation is impacted by the probability of communicating with LoS, as we
account for in our analysis. That behavior is properly described in [15],
where the authors highlight a number of factors that affect the air-to-ground
communications.

Thus, in general, satellites, balloons or terrestrial relays cannot face
scenarios as the ones studied in this paper. In fact, connectivity
requirements used to design protocols for satellite, balloon and D2D
communications, as well as technology constraints and signal propagation, are
radically different from our case. An overview of the recent advances in drone
communications, with an emphasis on integrating drones into the fifth-
generation can be found in [16].

###  2.2 Drone Position Optimization

In the recent years there have been various studies that optimize drone relay
placement mainly focusing on coverage in static and oversimplified
assumptions, as for instance neglecting inter-drone interference [17] or
ignoring fairness issues in resource allocation [18]. With that, the resulting
problem formulation is simple enough, typically quadratic, yet less realistic
than what we derive in this paper. In our previous work, we have presented
OnDrone, a lightweight multi-drone coverage framework that maximizes the
number of users covered. In there, we consider Quality-of-Service thresholds
and show how to design flying paths for fast repositioning and augmented
coverage over time [19]. A similar approach has been followed in [20], where
the authors propose the deployment of drones for video surveillance, obtaining
a suitable topology (in terms of achieving zero blind spots) with minimum
number of drones. However, they assume a region that lacks cellular
infrastructure. A different approach is taken in [21], where a Q-learning
solution is used to solve the deployment and movement of the drones that
maximizes the Quality of Experience [22], although it requires explicit
feedback from ground users and protocol changes in the way cells form. Galkin
et al. [23] use a stochastic model to analyze the coverage when drones act as
access points for users in urban areas. They demonstrate how the density of
the drones determines whether drones should position closer to user hotspots
to improve the received signal strength, or further away from one another to
mitigate interference. By using stochastic and fractal geometry as a
macromodel for both the vehicular networks and for the fixed infrastructures,
in [24] the authors derived analytic bounds in terms of connectivity. They
discovered two regimes, in function of the existing infrastructure, and
identified in which regime the deployment of nodes is desirable. In [25], the
authors propose a communications framework for HetNets composed by macrocells,
microcells and communication drones. Although the application scenario is
similar to the one of this paper, we find remarkable and fundamental
differences: they consider only a finite and small set of aerial locations for
drones, ignore key communication features such as throughput or coverage and
focus the decisions on minimizing the network energy consumption.

Although all these works provide valuable contribution and foundational
results, they do not shed light on problems like capacity optimization with
fairness targets when a fleet of drones is deployed to assist a cellular
network, which is our goal in this paper. Neither they provide a realistic
framework to integrate drone relay stations into current cellular networks,
whereas we show that our practical approach accounts for integration. In the
conference version of this paper [1], we address this problem from a quite
simple approach: in there, we relax critical constraints like wireless
backhaul and backbone bottlenecks, and air-space constraints where drones can
be located. Such simplifications allow for the design of a simpler system
model and less complex algorithms and optimizations, yet bring the risk of
being inaccurate. In this article, instead, we provide a fully general,
realistic, precise and tractable model for an aerial network integrated in a
cellular infrastructure.

Besides, here we use EO in a novel way. EO algorithms took relevance on fields
as biology or physics, and have been also applied to network science, e.g., to
optimize transport in complex networks [26]. EO has been thought and used so
far as a form of centralized static optimization. Instead, in this paper we
use EO to design an algorithm that works in parallel threads and dynamically,
as the system evolves.

Figure 2: System Model

##  3 System Model

Our goal is to derive an analytical framework that finds optimal 3D locations
of drone relays, given the position of users and ground stations, to optimize
downlink communications. We target $\alpha$-fair instantaneous user
throughput. To this aim, jointly with path-loss and interference in both
access and backhaul, we address how users perform cell selection and get
resources allocated. The modelled system is a plain 3GPP cellular system
comprising access and backhaul networks, with drones acting as standard 5G
relays [27].

In order to illustrate the system model and parameters involved in the
network, we show the main components in Fig. 2 and describe the meaning of the
parameters in Table I.

TABLE I: System model parameters Parameter | Description
---|---
$\alpha$ | Fairness level
$A_{g}$ | Maximum number of $a\\!B\\!S$s that $g\\!B\\!S$ $g$ can serve
$\mathcal{A}$ | Set of $a\\!B\\!S$s
$\mathcal{B}$ | Set of BSs, i.e., set of $a\\!B\\!S$s and $g\\!B\\!S$s jointly
$\mathcal{G}$ | Set of $g\\!B\\!S$s
$\gamma_{a,u}^{\mathcal{A}}$ | SINR of aerial access link $(a,u)$
$\gamma_{g,a}^{\mathcal{B}}$ | SINR of aerial backhaul link $(g,a)$
$\gamma_{g,u}^{\mathcal{G}}$ | SINR of ground access link $(g,u)$
$I_{a,u}^{\mathcal{A}}$ | Interference from $a\\!B\\!S$ $a$ to UE $u$
$I_{g,a}^{\mathcal{B}}$ | Interference from $g\\!B\\!S$ $g$ to $a\\!B\\!S$ $a$
$I_{g,u}^{\mathcal{G}}$ | Interference from $g\\!B\\!S$ $g$ to UE $u$
$L_{\mathcal{A}}(a,u)$ | Loss of air-to-ground channel from $a\\!B\\!S$ $a$ to UE $u$
$L_{\mathcal{B}}(g,a)$ | Loss of ground-to-air channel from $g\\!B\\!S$ $g$ to $a\\!B\\!S$ $a$
$L_{\mathcal{G}}(g,u)$ | Loss of ground-to-ground channel from $g\\!B\\!S$ $g$ to UE $u$
$\phi_{g,a}$ | Angle from main lobe of $g\\!B\\!S$ $g$ to $a\\!B\\!S$ $a$
$P_{L\\!o\\!S}(a,u)$ | LoS-likelihood of link $(a,u)$
$U_{\max}$ | Maximum number of UEs a BS can serve
$W_{\mathcal{A}}$ | Available bandwidth of aerial access channels
$W_{\mathcal{A}}^{\min}$ | Minimum bandwidth guaranteed to aerial access links
$W_{\mathcal{B}}$ | Available bandwidth of aerial backhaul channels
$W_{\mathcal{B}}^{\min}$ | Minimum bandwidth guaranteed to aerial backhaul links
$W_{\mathcal{G}}$ | Available bandwidth of ground access channels
$W_{\mathcal{G}}^{\min}$ | Minimum bandwidth guaranteed to ground access links
$\tau_{g}$ | Backbone capacity of $g\\!B\\!S$ $g$

###  3.1 Reference Scenario

We consider a flat ground surface $\mathcal{S}$ where a set $\mathcal{G}$ of
$G$ ground base stations ($g\\!B\\!S$) provide cellular service. The position
of each $g\\!B\\!S$ $g$ is denoted as $\Pi^{g}=(X^{g},Y^{g})$. We assume that
every $g\\!B\\!S$ $g$ is wired to the internet with a backbone capacity
$\tau_{g}$. We consider that a set $\mathcal{U}$ of $U$ user equipments (UEs)
is on the ground, requesting cellular service. We denote the position of each
user $u$ as $\pi_{u}=(x_{u},y_{u})$. The network disposes of a fleet
$\mathcal{A}$ of $A$ aerial base stations ($a\\!B\\!S$) that act as mobile
relays. Each $a\\!B\\!S$ is mounted on a drone (for ease of readability, we
may refer just as “drone” to an $a\\!B\\!S$). We assume that the system
operator disposes of two frequency bands to offer user access: the first band
is used by $g\\!B\\!S$s to provide both cellular access to ground users and
backhaul connectivity to $a\\!B\\!S$s; and the second band is used by
$a\\!B\\!S$s to provide cellular access to ground users. This assumption
mimics the common scenario in which an operator deploys different types of
base stations and equip them with different frequencies to simplify network
management. Indeed, it is common that an operator disposes of multiple
licenses to use more than one cellular band. It is also common practice to
coordinate transmissions from different base stations, especially in case of
mixed deployments with macro and small cells. Frequency reuse is one of the
possible and simplest schemes to use for coordination, which is what we adopt
for simplicity. Nonetheless, please note that the analysis carried out in this
article can be extended to the case of using a single frequency for all
transmissions or any arbitrary number of frequencies (e.g., one per aerial and
ground base station, in the most extreme case).

Under our assumptions, users connected to drones suffer no interference from
$g\\!B\\!S$s. Similarly, users connected to ground base stations are not
interfered by drone transmissions. Drones fly in the air space, so that we
denote as $\Pi_{a}=(X_{a},Y_{a},h_{a})$ the 3D position of drone $a$. We
denote as $\mathcal{B}=\mathcal{G}\cup\mathcal{A}$ the set of all the base
stations ($g\\!B\\!S$s jointly with $a\\!B\\!S$s) that form the whole network.

We assume that $g\\!B\\!S$s reuse the downlink spectrum used for the
$g\\!B\\!S$–UE access links in order to provide backhaul wireless service to
$a\\!B\\!S$s. Usually, $g\\!B\\!S$s dispose of three antenna sectors pointing
mainly to the ground in order to provide cellular coverage to users. In
addition to these antennas, we assume that $g\\!B\\!S$s dispose of an
additional full dimensional antenna array that performs 3D-beamforming over
clear LoS links in order to establish $g\\!B\\!S$–$a\\!B\\!S$ links for
backhauling, as suggested and studied in [28]. Hence, $g\\!B\\!S$–UE and
$g\\!B\\!S$–$a\\!B\\!S$ links do not practically interfere. Moreover, 3D
antenna arrays allow $g\\!B\\!S$s to relay traffic to several drones, by
alternating transmissions every few milliseconds. Therefore, we assume that
each $g\\!B\\!S$ can set backhaul wireless links with more than one
$a\\!B\\!S$, and we also assume that a minimum backhaul bandwidth
$W_{\mathcal{B}}^{\min}$ is guaranteed, which implies that only a limited
number of drones can stay connected to a single ground base station.

In order to analyze a realistic network configuration, differently from what
we assume for BSs talking to each other, we do not assume that transmissions
to users can benefit from beamforming, therefore, we account for inter-cell
interference for the groups of users connected to the same type of BS, be it
$a\\!B\\!S$ or $g\\!B\\!S$. We do so because beamforming to users requires
non-negligible overheads to track the position of erratically moving users,
and to fine tune the beam within potentially narrow spaces crowded with users.
Instead, $a\\!B\\!S$s and $g\\!B\\!S$ can easily have precise information
about their positions at any time, and do not risk to stay packed in small
spaces. Having said that, our analysis can be straightforwardly simplified to
apply also to the ideal case in which all transmissions can be interference-
free due to beamforming, should that possibility become reality.

Another important aspect to consider in drone-assisted missions is the energy
consumption management of drones, generally by means of a recharge strategy
with additional drones. In this paper, we assume that $A$ drones are
guaranteed to be operational in a desired position when the optimization is
run. Hence, we assume that there is a parallel recharge scheduling strategy
that uses a minimum number of extra drones to guarantee a persistent
communications service with $A$ $a\\!B\\!S$s at any moment. This parallel
mechanism can certainly be achieved, as we have already proposed and
demonstrated in [29].

###  3.2 Channel Modelling

We assume that the network operator disposes of two orthogonal frequency
bands. One band is assigned to $g\\!B\\!S$s to provide access service to
ground users as well as aerial backhaul service to $a\\!B\\!S$s, with a fixed
bandwidth $W_{\mathcal{G}}$ and $W_{\mathcal{B}}$, respectively. The other
band is assigned to $a\\!B\\!S$s for aerial user access and it has a fixed
bandwidth $W_{\\!\mathcal{A}}$. Hence, we model three different kinds of
channels: (i) air-to-ground and (ii) ground-to-ground channels in the access
network, and (iii) ground-to-air channels in the backhaul network.

To optimize the network and compute link throughputs we need to compute the
signal strength of each wireless channel, measured as the SINR. In order to
compute the SINR, we study path-loss and fading of channels, as well as
suffered interference. Then, throughputs will be measured with the well known
Shannon formula.

The notation of each parameter of the channel model can be found in Table I,
and the mathematical details of the modelling of each kind of channel can be
found in separate sections of the online supplemental material, within Section
A.

Air-to-ground access channels. The loss of these channels differs notably
depending on whether links are free of obstacles [30]. Hence, this loss
depends on the LoS-likelihood [31] of links. As $a\\!B\\!S$s serve users in
OFDMA channels, these users do not interfere among themselves. However, users
suffer aggregated interference from other $a\\!B\\!S$s serving other users. We
denote the loss from $a\\!B\\!S$ $a$ to UE $u$ as $L_{\mathcal{A}}(a,u)$.

Ground-to-ground access channels. The loss of these channels follows the well
known path-loss model with slow fading [32]. As $g\\!B\\!S$s serve users in
OFDMA channels, there is no intra-cell interference. However, we account for
inter-cell interference from neighbour ground cells. We denote the loss from
$g\\!B\\!S$ $g$ to UE $u$ as $L_{\mathcal{G}}(g,u)$.

Ground-to-air backhaul channels. Since backhaul links follow directional beams
in LoS, the loss of these channels also follows the well known path-loss model
with slow fading [32]. Backhaul links reuse the spectrum from the ground
access channels. However, $g\\!B\\!S$s perform 3D beamforming for backhauling
and hence there is no interference between ground access and aerial backhaul.
Still, other $g\\!B\\!S$s might serve other $a\\!B\\!S$s in a beamforming
direction such that served $a\\!B\\!S$s suffer interference (attenuated
depending on the radiating angle). We account for such interference. We denote
the loss from $g\\!B\\!S$ $g$ to $a\\!B\\!S$ $a$ as $L_{\mathcal{B}}(g,a)$.

###  3.3 Cell Selection and Resource Allocation

BSs cannot provide service to unlimited numbers of users because $(i)$
available radio resources are limited and $(ii)$ it is necessary to guarantee
a minimum set of radio resources to each connected user, to guarantee
signaling exchange with the BS; this is needed to schedule data transmissions
to and from the BS. Of course, the number of devices is also limited by the
minimum bandwidth that the system aims to guarantee to each user. Therefore,
in general, the maximum number of users that can be simultaneously served is
limited, and we denote by $U_{\max}$ such number.

Cell Selection. We assume that users perform cell selection as in currently
operational 3GPP networks: first, UEs select the BS with strongest Signal-to-
Noise Ratio (SNR); if the request is rejected because channel conditions
deteriorate or the BS runs at maximum capacity, then the UE performs cell re-
selection, and tries to attach to the BS with next strongest SNR, and so on
until the user gets attached [33].

Resource Allocation. We assume that $g\\!B\\!S$s and $a\\!B\\!S$s schedule
cellular users according to an OFDMA system. Today’s BSs use an OFDMA system
and dispose of a finite set of physical resource blocks organized in
subframes, which repeat to form frames lasting a few milliseconds (1 to 10 ms
in 3GPP-compliant networks). A physical resource block is the smallest unit of
time-frequency resources that can be allocated to a user. Thus, we assume that
the minimum bandwidth allocated to a user is the bandwidth corresponding to
one resource block and the scheduler guarantees that each user receives, on
average, at least one block per subframe in each OFDMA frame. We denote as
$W_{\mathcal{G}}^{\min}$ and $W_{\\!\mathcal{A}}^{\min}$ the minimum bandwidth
that a $g\\!B\\!S$ or an $a\\!B\\!S$ can allocate to a single user.

Backhaul links also use an OFDMA system, although $a\\!B\\!S$s select a
$g\\!B\\!S$ according a the global network optimization criterion rather than
based on SNR. Moreover, $g\\!B\\!S$s can serve $A_{g}$ $a\\!B\\!S$s at most
and each backhaul link $(g,a)$ disposes of a minimum bandwidth
$W_{\mathcal{B}}^{\min}$ to relay traffic.

##  4 Optimization

Here we derive an analytic framework for the 3D positions of $a\\!B\\!S$s, to
optimize throughputs based on $\alpha$-fairness [10]. We seek optimal
$a\\!B\\!S$ positions and BS–UE associations, and also optimal backhaul
association and optimal allocation of physical resources. We formally present
the optimization problem addressed in this paper:

Throughput Problem $\boldsymbol{\mathcal{T}}$: Given a set $\mathcal{G}$ of
$G$ fixed $g\\!B\\!S$s, a fleet $\mathcal{A}$ of $A$ relay $a\\!B\\!S$s
hovering at heights in the range $[h_{\min},h_{\max}]$, a set $\mathcal{U}$ of
$U$ ground UEs that may connect to either a $g\\!B\\!S$ or an $a\\!B\\!S$s,
each of which can serve $U_{\max}$ UEs at most, find the optimal position of
each $a\\!B\\!S$ $a\\!\in\\!\mathcal{A}$, the optimal user association, the
optimal backhaul association and the optimal user resource allocation so to
maximize the $\alpha$-fair throughput utility function.

We formulate the drone positioning problem as a Mixed-Integer Non-Convex
Program (MINCP) in Fig. 3. The meaning of all the optimization variables is
shown in Table II.

On the access network side, we denote as $C_{a,u}\\!\in\\!\\{0,1\\}$ the
decision variable that tells whether $u$ connects to $a\\!B\\!S$ $a$.
Similarly, $C_{u}^{g_{u}}\\!\in\\!\\{0,1\\}$ tells whether $u$ connects to
$g\\!B\\!S$ $g_{u}$. Decision variables $W_{b,u}$ and $T_{b,u}$ denote
bandwidth and throughput allocated to link $(b,u)$. Throughput is a decision
variable and not directly computed with the Shannon formula, since, in
addition to bandwidth limitations we must account for access, backhaul and
backbone bottlenecks.

On the backhaul network side, we denote as $B^{g,a}\\!\\!\in\\!\\!\\{0,1\\}$
the decision variable that tells whether $a\\!B\\!S$ $a$ is attached to
$g\\!B\\!S$ $g$. Variables $W^{g,a}$ and $T^{g,a}$ denote bandwidth and
throughput of the backhaul link $(g,a)$, respectively.

#### Network utility and constraints

The utility function of the optimization problem in Fig. 3 corresponds to the
$\alpha$-fairness metric. This utility function is additive in terms of
utilities conveyed by single BSs. Depending on the value of
$\alpha\\!\geq\\!0$, known as the $\alpha$-fairness level, the metric captures
different fairness criteria such as weighted proportional fairness
($\alpha\\!=\\!1$), $\max$-$\min$ fairness
($\alpha\\!\rightarrow\\!+\\!\infty$) or the maximum capacity ($\alpha=0$).
Constraints in Fig. 3 correspond to the following restrictions:

$\begin{array}[]{ll}\\!\quad\quad\quad\quad\quad\quad\quad\qquad\qquad\quad\mathclap{\max\limits_{\begin{subarray}{c}\Pi_{a},C_{a,u},W_{b,u},\\\
B^{g,a},W^{g,a}\end{subarray}}\quad\\!\\!U^{\alpha}_{thr}\\!=\\!\begin{cases}\sum\limits_{u\in\mathcal{U}}\\!\\!\left(\sum\limits_{b\in\mathcal{B}}T_{b,u}\\!\right)^{1-\alpha}\\!\\!\\!\cdot\\!\frac{1}{1-\alpha},&\alpha\\!\neq\\!1;\\\
\sum\limits_{u\in\mathcal{U}}\\!\log\\!\left(\sum\limits_{b\in\mathcal{B}}T_{b,u}\\!\right),&\alpha\\!=\\!1;\end{cases}}\\\
\text{s.t.:}&\\\ \textit{$g\\!B\\!S$--$a\\!B\\!S$ association
constraints:}&\\\
\sum\limits_{g\in\mathcal{G}}B^{g,a}=1,\quad\sum\limits_{a\in\mathcal{A}}B^{g,a}\leq
A_{g},&\forall g\in\mathcal{G},\forall a\in\mathcal{A};\\\ \\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\textit{$g\\!B\\!S$--$a\\!B\\!S$
capacity constraints:}}&\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W_{\mathcal{B}}^{\min}\cdot
B^{g,a}\leq W^{g,a}\leq W_{\mathcal{B}}\cdot
B^{g,a},}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\forall
g\in\mathcal{G},\forall a\in\mathcal{A};}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum\limits_{a\in\mathcal{A}}W^{g,a}\leq
W_{\mathcal{B}},}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\forall
g\in\mathcal{G};}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}T^{g,a}\leq
W^{g,a}\log_{2}\left(1+\gamma_{g,a}^{\mathcal{B}}\right),}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\forall
g\in\mathcal{G},\forall a\in\mathcal{A};}\\\ \\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\textit{$g\\!B\\!S$
backbone constraint:}}&\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum\limits_{a\in\mathcal{A}}T^{g,a}+\sum\limits_{\begin{subarray}{c}u\in\mathcal{U}:\\\
g=g_{u}\end{subarray}}T_{g,u}\leq\tau_{g},}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\forall
g\in\mathcal{G};}\\\ \\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\textit{BS
--UE association constraints:}}&\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}C_{u}^{g_{u}}+\sum\limits_{a\in\mathcal{A}}C_{a,u}=1,}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\forall
u\in\mathcal{U};}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum\limits_{\begin{subarray}{c}u\in\mathcal{U}:\\\
g=g_{u}\end{subarray}}C_{u}^{g}\leq
U_{\max},\quad\sum\limits_{u\in\mathcal{U}}C_{a,u}\leq
U_{\max},}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\forall
g\in\mathcal{G},\forall a\in\mathcal{A};}\\\ \\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\textit{$g\\!B\\!S$--UE
capacity costraints:}}&\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W_{\mathcal{G}}^{\min}\cdot
C^{g_{u}}_{u}\leq W_{g_{u},u}\leq W_{\mathcal{G}}\cdot
C^{g_{u}}_{u},}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\forall
u\in\mathcal{U};}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum\limits_{u\in\mathcal{U}}W_{g,u}\leq
W_{\mathcal{G}},}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\forall
g\in\mathcal{G};}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}T_{g_{u},u}\leq
W_{g_{u},u}\log_{2}\left(1+\gamma^{\mathcal{G}}_{g_{u},u}\right),}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\forall
u\in\mathcal{U};}\\\ \\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\textit{$a\\!B\\!S$--UE
capacity constraints:}}&\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W_{\\!\mathcal{A}}^{\min}\cdot
C_{a,u}\leq W_{a,u}\leq W_{\\!\mathcal{A}}\cdot
C_{a,u},}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\forall
a\in\mathcal{A},\forall u\in\mathcal{U};}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum\limits_{u\in\mathcal{U}}W_{a,u}\leq
W_{\\!\mathcal{A}},}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\forall
a\in\mathcal{A};}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}T_{a,u}\leq
W_{a,u}\log_{2}\left(1+\gamma^{\mathcal{A}}_{a,u}\right),}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\forall
a\in\mathcal{A},\forall u\in\mathcal{U};}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum\limits_{u\in\mathcal{U}}T_{a,u}\leq\sum\limits_{g\in\mathcal{G}}T^{g,a},}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\forall
a\in\mathcal{A};}\\\ \\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\textit{Air
space
constraint:}\,\,\,\,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Pi_{a}\in\mathcal{S}_{a},}}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\forall
a\in\mathcal{A}.}\end{array}$

Figure 3: Fair Throughput Optimization Program.

$g\\!B\\!S$–$a\\!B\\!S$ association constraints state that every $a\\!B\\!S$
must associate with one $g\\!B\\!S$ to set a wireless backhaul link, and that
every $g\\!B\\!S$ can serve at most $A_{g}$ $a\\!B\\!S$s.

$g\\!B\\!S$–$a\\!B\\!S$ capacity constraints impose guarantees on bandwidth
and throughput allocation on backhaul links.

The $g\\!B\\!S$ backbone constraint restricts every $g\\!B\\!S$ to provide
connected $a\\!B\\!S$s and connected users with an aggregated throughput not
higher than the capacity $\tau_{g}$ of the backbone link that serves that
$g\\!B\\!S$.

BS–UE association constraints state that each user can associate only to one
BS, either a ground station or a drone, and that the maximum allowed number of
UEs served by each BS is limited to $U_{\max}$.

$g\\!B\\!S$–UE capacity constraints impose guarantees on bandwidth and
throughput allocation on $g\\!B\\!S$–UE links.

$a\\!B\\!S$–UE capacity constraints impose guarantees on bandwidth and
throughput allocation on $a\\!B\\!S$–UE links, while at the same time the
throughput of a user cannot exceed the backhaul link capacity and the
aggregate user throughput cannot exceed the backhaul throughput.

The Air space constraint delimits the 3D air space in within which an
$a\\!B\\!S$ can be moved, $\mathcal{S}_{a}$, and which is a ball centered in
the current position of the drone with a radius equal to the distance that the
drone can fly within a fixed time (i.e., time itself is the real constraint).

Modeling air-to-ground connections brings unavoidable non-convex functions, so
that the formulated problem is not tractable with currently available
optimizers, which are able to deal only with problems that are convex.

#### Feasibility and NP-Completeness

Here we remark that the formulated optimization problem in Fig. 3 always
admits a feasible solution. We just have to place each $a\\!B\\!S$
$a\\!\in\\!\mathcal{A}$ randomly in each $\mathcal{S}_{a}$, attach each UE and
$a\\!B\\!S$ to a random $g\\!B\\!S$, allocate to each UE and $a\\!B\\!S$ the
minimum allocatable bandwidth and assign null throughputs to every link in the
network. This naive network setting is always a feasible solution of the
problem.

Moreover, the problem of finding the exact optimal drone positions is NP-
Complete. Indeed, the NP-Complete Minimum-Geometric Disk-Cover (MGDC) problem
[34] can be reduced, in polynomial time, to a special instance of the problem
where users get 1 bps if a drone serves them and 0 bps otherwise. This result
is a direct consequence of the NP-Completeness proof of the Coverage Problem
that we have presented in [19], because coverage can be seen as a particularly
simplified throughput problem using an on/off, SINR-threshold-based,
throughput function.

TABLE II: Optimization variables & parameters. Variable | Description
---|---
$B^{g,a}$ | Binary variable for the existence of backhaul link $(g,a)$
$C_{a,u}$ | Binary variable for the existence of access link $(a,u)$
$C_{u}^{g_{u}}$ | Binary variable for the existence of access link $(g_{u},u)$
$\Pi_{a}$ | Aerial position of $a\\!B\\!S$ $a$
$\mathcal{S}_{a}$ | Air space restricted to $a\\!B\\!S$ $a$
$T_{b,u}$ | Throughput allocated to access link $(b,u)$
$T^{g,a}$ | Throughput allocated to backhaul link $(g,a)$
$W_{b,u}$ | Bandwidth allocated to access link $(b,u)$
$W^{g,a}$ | Bandwidth allocated to backhaul link $(g,a)$

##  5 Extremal Optimization

The optimization framework proposed in Section 4 is non-convex and mixed-
integer, hence not solvable with any off-the-shelf optimizer [35]. The problem
is hard to solve because any change in a decision variable (e.g., a position
of a drone) affects backhaul and user association, as well as interference and
resource allocation for all users. Due to the intertwined nature of the
decision variables of the problem, this is the kind of problems EO has been
thought for.

To find time-efficient and near-optimal solutions, we propose a Parallelized
Alfa-fair Drone Deployment (PADD) algorithm, based on a typical EO operation.
We base the design of PADD on decoupling the problem in the four main
decisions that the optimization framework must make: $(i)$ the 3D positions of
the fleet of $a\\!B\\!S$s; $(ii)$ the sets of users attached to each
$g\\!B\\!S$ and $a\\!B\\!S$; $(iii)$ $g\\!B\\!S$–$a\\!B\\!S$ backhaul
association; $(iv)$ bandwidth allocation to backhaul links
$g\\!B\\!S$–$a\\!B\\!S$ as well as to access links from each $g\\!B\\!S$ or
$a\\!B\\!S$ to their attached users.

In what follows, we formally describe the PADD operation and provide details
of each step that PADD takes. The algorithm iteratively solves four steps, as
pictured in Fig. 4 as a flow-chart: after deriving an initial feasible system
setting, the least fit $a\\!B\\!S$ $a_{0}$ is selected in order to locally
probe non-searched positions that improve the current network performance.
While probed positions do not provide a relative improvement of $\delta\geq 0$
over the current network performance, new local non-searched positions are
probed. In case a probed position improves performance, $a\\!B\\!S$ $a_{0}$ is
moved to such position and the current system performance is updated. In case
this new performance provides a relative improvement higher than
$\varepsilon\geq\delta$, the new least fit $a\\!B\\!S$ is selected and the
process begins again. Otherwise, the algorithm converges and outputs the
decided system setting. We next formalize the operation of PADD more in
detail.

Initial system setting. Performance: $\mathcal{U}_{\mathrm{thr},0}^{\alpha}$
Select $a_{0}\\!\\!\in\\!\\!\mathcal{A}$ as the least fit $a\\!B\\!S$ (see
§5.5) Set:
$\mathcal{U}_{\mathrm{thr},0}^{\alpha}\\!\leftarrow\\!\mathcal{U}_{\mathrm{thr}\\!,\mathrm{new}}^{\alpha}$
Move $a\\!B\\!S$ $a_{0}$ locally to a non-search position BS–UE association:
best-signal policy (see §5.2)
$\frac{\mathcal{U}_{\mathrm{thr},\mathrm{new}}^{\alpha}\\!-\\!\mathcal{U}_{\mathrm{thr},0}^{\alpha}}{\mathcal{U}_{\mathrm{thr},0}^{\alpha}}\\!\\!<\\!\\!\varepsilon$
$g\\!B\\!S$–$a\\!B\\!S$ backhaul association: knap heuristic (see §5.3)
Optimal bandwidth allocation (CP), $\forall g\\!\in\\!\mathcal{G}$. Output:
$\mathcal{U}_{\mathrm{thr},\mathrm{new}}^{\alpha}$ (see §5.4) Allocate to
$a_{0}$ this new position
$\frac{\mathcal{U}_{\mathrm{thr},\mathrm{new}}^{\alpha}\\!-\\!\mathcal{U}_{\mathrm{thr},0}^{\alpha}}{\mathcal{U}_{\mathrm{thr},0}^{\alpha}}\\!\\!>\\!\\!\delta$
Decided system setting. Performance: $\mathcal{U}_{\mathrm{thr}}^{\alpha}$
yesnonoyes Figure 4: Flow diagram of PADD operation.

###  5.1 Initial System Setting

Initially, we consider a naïve drone positioning. For instance, any random
placement of drones provides a feasible solution that could be iteratively
improved. However, we more efficiently select those locations that are closer
to $g\\!B\\!S$s–so to guarantee good backhaul links–and those locations that
are above densely populated regions–which are regions that potentially need
drone relay assistance. This initial drone setting yields utility
$U_{\mathrm{thr},0}^{\alpha}$. However, in order to compute
$U_{\mathrm{thr},0}^{\alpha}$, we need to know also user association, backhaul
association and resource allocation. These decisions are made as described in
the following subsections.

###  5.2 BS–UE Association: Best-Signal Policy

For fixed positions of the drones, the BS–UE association is performed
individually by each user according to the best-signal policy described in
Section 3.3, which is at most linear in the number of base stations
$|\mathcal{B}|$, with $\mathcal{B}=\mathcal{G}\cup\mathcal{A}$. However,
before association, we need to compute $|\mathcal{B}|$ SNR values for each
UEs, and sort them in decreasing order, which goes with
$|\mathcal{B}|\log|\mathcal{B}|$. The complexity of user association is
therefore $\mathcal{O}(|\mathcal{B}|\log|\mathcal{B}|)$ comparisons, for each
user.
In the following, the set of users attached to a $g\\!B\\!S$ $g\in\mathcal{G}$
and the set of users attached to an $a\\!B\\!S$ $a\in\mathcal{A}$ are denoted
as $\mathcal{U}_{g}$ and $\mathcal{U}_{a}$, respectively.

###  5.3 $g\\!B\\!S$–$a\\!B\\!S$ Backhaul Association: GAP-Knap Heuristic

For fixed drone positions and given BS–UE association, the backhaul
association is solved by assuming that the bandwidth $W_{\mathcal{G}}$ is
proportionally shared among users connected to a $g\\!B\\!S$ and to the
$a\\!B\\!S$s to be connected to that $g\\!B\\!S$.

Specifically, the backhaul throughput is computed with the Shannon formula,
using a fraction of bandwidth proportional to the number of users attached to
the drone, and with the SINR resulting from current position of drone and the
fixed position of $g\\!B\\!S$. Maximizing the $\alpha$-fairness of such
throughput values for all drones, translates into a generalized assignment
problem (GAP) [36]. With the notation used in Fig. 3, in which
$\gamma_{g,a}^{\mathcal{B}}$ is the backhaul SINR of link $(b,a)$ and
$B^{g,a}$ denotes the binary decision variable that tells wether $g\\!B\\!S$
$g$ associates with $a\\!B\\!S$ $a$, we show the backhaul association
optimization in Fig. 5.

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\begin{array}[]{ll}\\!\max\limits_{B^{g,a}}\,U_{\mathrm{bhl}}^{\alpha}\\!=\\!\\!\begin{cases}\\!\sum\limits_{g\in\mathcal{G}}\sum\limits_{a\in\\!\mathcal{A}}\\!\\!\\!\left(\\!\frac{\left|\mathcal{U}_{a}\right|}{\left|\mathcal{U}_{g}\bigcup\mathcal{U}_{a}\right|}\\!\log_{2}\\!\\!\left(1\\!+\\!\gamma_{g,a}^{\mathcal{B}}\\!\right)\\!\\!\right)^{\\!\\!1-\alpha}\\!\\!\\!\\!\cdot\\!\frac{B^{g,a}}{1\\!-\\!\alpha},&\\!\\!\alpha\\!\neq\\!1;\\\
\\!\sum\limits_{g\in\mathcal{G}}\sum\limits_{a\in\\!\mathcal{A}}\\!\log\\!\left(\\!\frac{\left|\mathcal{U}_{a}\right|}{\left|\mathcal{U}_{g}\bigcup\mathcal{U}_{a}\right|}\\!\log_{2}\\!\\!\left(1\\!+\\!\gamma_{g,a}^{\mathcal{B}}\right)\\!\\!\right)\\!\\!\cdot\\!\\!B^{g,a}\\!\\!,&\\!\\!\alpha\\!=\\!1;\\\
\vspace{-3mm}&\end{cases}\\\ \text{s.t.:}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum\limits_{a\in\mathcal{A}}B^{g,a}\leq
A_{g},\,\quad\quad\forall g\\!\in\\!\mathcal{G};}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum\limits_{g\in\mathcal{G}}B^{g,a}=1,\qquad\quad\forall
a\\!\in\\!\mathcal{A}.}\end{array}}$

Figure 5: Backhaul Association Optimization Program.

The first constraint states that each $g\\!B\\!S$ can provide backhaul service
to at most $A_{g}$ $a\\!B\\!S$s. The second constraint states that each
$a\\!B\\!S$ $a$ must be associated with one $g\\!B\\!S$.

GAP is an NP-hard MILP [37]. Hence, although the GAPs that PADD needs to solve
have a small size and could be optimally solved by means of standard methods
as a Branch&Bound search [38], such an approach would not lead to a
polynomial-time algorithm. Hence, multiple heuristics have been proposed in
literature to find near-optimal solutions to the GAP [39, 37]. In particular,
we perform a simple heuristic based on an approximation to the $0$-$1$
knapsack problem [40] by means of dynamic programming [41]. We name this
heuristic as GAP-knap, which has a polynomial complexity of
$\mathcal{O}(|\mathcal{G}|\\!\cdot\\!|\mathcal{A}|\\!\cdot\\!A_{g})$. The
complexity of GAP-knap is linear with the sizes of $\mathcal{G}$ and
$\mathcal{A}$ and the maximum number of $a\\!B\\!S$s allowed to attached to
$g\\!B\\!S$s, $A_{g}$.

###  5.4 Optimal Bandwidth Allocation: Convex Program

For fixed drone positions, BS–UE association and $g\\!B\\!S$–$a\\!B\\!S$
backhaul association, the optimal bandwidth allocation is solved by each
$g\\!B\\!S$ in parallel, with a convex program.

Each $g\\!B\\!S$ must perform bandwidth allocation to split backhaul resources
among the served $a\\!B\\!S$s, split $g\\!B\\!S$–UE access resources among the
attached users $\mathcal{U}_{g}$, and the attached $a\\!B\\!S$s $a$ must split
$a\\!B\\!S$–UE access resources among the users attached to each $a\\!B\\!S$
$a$, $\mathcal{U}_{a}$. Since all these bandwidth allocations are intertwined
($a\\!B\\!S$–UE allocation depends on the backhaul bottleneck,
$g\\!B\\!S$–$a\\!B\\!S$ allocation depends on the number of attached
$a\\!B\\!S$s to the same $g\\!B\\!S$ and the number of final users of each
$a\\!B\\!S$, and $g\\!B\\!S$–UE access resources depend also on the backbone
service that the $g\\!B\\!S$ gets from the internet), we formulate the
bandwidth allocation optimization as a convex program in Fig. 6.

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\begin{array}[]{ll}\,\,\quad\qquad\qquad\qquad\qquad\qquad\mathclap{\max\limits_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w^{a},}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w_{u},}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}T_{u}}}U_{\mathrm{cvx},g}^{\alpha}=\begin{cases}\sum\limits_{u\in\\!\\!\bigcup\limits_{b\in\mathcal{B}_{g}}\\!\\!\\!\mathcal{U}_{b}}\left(T_{u}\right)^{1-\alpha}\cdot\frac{1}{1-\alpha},&\alpha\neq
1;\\\
\sum\limits_{u\in\\!\\!\bigcup\limits_{b\in\mathcal{B}_{g}}\\!\\!\\!\mathcal{U}_{b}}\log\left(T_{u}\right),&\alpha=1;\end{cases}}\\\
\text{s.t.:}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w^{a}\geq
W^{\min}_{\mathcal{B}},}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\\!\\!\forall
a\\!\in\\!\mathcal{A}_{g};}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum\limits_{a\in\mathcal{A}_{g}}\\!w^{a}=W_{\mathcal{B}};}&\\!\\!\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}T^{a}\\!\leq\\!w^{a}\log_{2}\\!\left(1\\!+\\!\gamma_{g,a}^{\mathcal{B}}\right)\\!,}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\\!\\!\forall
a\in\mathcal{A}_{g};}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w_{u}\geq
W^{\min}_{\mathcal{G}},}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\\!\\!\forall
u\\!\in\\!\mathcal{U}_{g};}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum\limits_{u\in\mathcal{U}_{g}}\\!w_{u}=W_{\mathcal{G}};}&\\!\\!\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w_{u}\geq
W^{\min}_{\mathcal{A}},}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\\!\\!\forall
u\\!\in\\!\bigcup\limits_{a\in\mathcal{A}_{g}}\\!\\!\mathcal{U}_{a};}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum\limits_{u\in\mathcal{U}_{a}}\\!w_{u}=W_{\\!\mathcal{A}},}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\\!\\!\forall
a\\!\in\\!\mathcal{A}_{g};}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}T_{u}\\!\leq\\!w_{u}\log_{2}\\!\left(1\\!+\\!\gamma_{b,u}\right)\\!,}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\\!\\!\forall(b,u)\\!\in\\!\mathcal{B}_{g}\\!\times\\!\\!\bigcup\limits_{b^{\prime}\in\mathcal{B}_{g}}\\!\\!\mathcal{U}_{b^{\prime}}\\!\,\mid\,\\!u\\!\in\\!\mathcal{U}_{b};}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum\limits_{u\in\mathcal{U}_{a}}\\!T_{u}\\!\leq
T^{a},}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\\!\\!\forall
a\\!\in\\!\mathcal{A}_{g};}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum\limits_{u\in\mathcal{U}_{g}}\\!T_{u}\\!+\\!\sum\limits_{a\in\mathcal{A}_{g}}\\!T^{a}\leq\tau_{g}.}&\\!\\!\vspace{-0.15mm}\end{array}}$

Figure 6: Bandwidth Allocation Optimization Program.

In the program, $w^{a}$ is the share of total bandwidth that $g\\!B\\!S$ $g$
allocates to $a\\!B\\!S$ $a$, $T^{a}$ is the backhaul throughput for
$a\\!B\\!S$ $a$, $w_{u}$ is the share bandwidth that BS $b$ allocates to user
$u$, and $T_{u}$ is the access service throughput for user $u$.
$\mathcal{A}_{g}$ is the set of $a\\!B\\!S$s attached to $g\\!B\\!S$ $g$, and
$\mathcal{B}_{g}=\\{g\\}\\!\cup\\!\mathcal{A}_{g}$ is the set of $g\\!B\\!S$
$g$ jointly with $\mathcal{A}_{g}$. $U_{\mathrm{cvx},g}^{\alpha}$ is the
utility function of $g\\!B\\!S$ $g$, which is based on the $\alpha$-fairness
metric.

The problem is convex, hence it is optimally solvable in polynomial time by
means of standard interior-point methods [42]. However, these methods run in
cubic time with respect to the number of users, i.e., the complexity is
$\mathcal{O}(|\mathcal{U}|^{3})$ [43]. This cubic complexity might soon become
prohibitive for real-time applications such as drone-aiding and fast
repositioning in wireless networks (specially for big populations), as
addressed in this paper. However, we have derived KKT conditions [44] for
problem of Fig. 6 that allow us to find the exact solution analytically (see
Section A in the online supplemental material). The complexity of finding the
exact solution for each $g\\!B\\!S$ $g$ is linear with respect to the number
of users and number of $a\\!B\\!S$s attached to $g\\!B\\!S$ $g$, i.e.,
$\mathcal{O}(U_{\max}\cdot A_{g})$, in the worst case.

###  5.5 Least Fit Drone Selection

PADD iteratively improves the utility until convergence. The algorithm uses
the idea behind EO algorithms. It selects the least fit element and re-sets
its parameters in order to improve system performance. In our case, an
$a\\!B\\!S$ is selected and a new location is probed.

The choice of which $a\\!B\\!S$ is the least fit is made based on the
consideration that there are two factors that cause sub-optimality of
$a\\!B\\!S$ positions: (i) the $a\\!B\\!S$ has a bad backhaul connectivity,
hence it provides a worse service to users than what the access channel
conditions allow, i.e., access resources are wasted; or (ii) the $a\\!B\\!S$
offers bad access connectivity to users due to inter-drone interference, even
though it has a good backhaul service, so that backhaul resources are wasted.
Accordingly, we derive two indicators of sub-optimality as the relative
difference between the aggregate utility from $a\\!B\\!S$-served users (namely
$U_{\mathrm{thr},a}^{\alpha}$) and the following quantities: (i) the utility
of the $a\\!B\\!S$ assuming infinite backhaul capacity, namely
$U_{\mathrm{thr},a}^{\alpha,B_{\infty}}$; and (ii) the utility of the
$a\\!B\\!S$ assuming no inter-drone interference, namely
$U_{\mathrm{thr},a}^{\alpha,S\\!N\\!R}\\!$. Eventually, we pick as least fit
$a\\!B\\!S$ the drone $a_{0}$ that provides the higher value of all sub-
optimality indicators:

$\displaystyle
a_{0}\\!=\\!\arg\max\limits_{a\in\mathcal{A}}\left(\\!\\!\max\\!\left(\left|\frac{U_{\mathrm{thr},a}^{\alpha,B_{\\!\infty}}\\!-\\!U_{\mathrm{thr},a}^{\alpha}}{U_{\mathrm{thr},a}^{\alpha,B_{\\!\infty}}}\right|\\!,\left|\frac{U_{\mathrm{thr},a}^{\alpha,S\\!N\\!R}\\!-\\!U_{\mathrm{thr},a}^{\alpha}}{U_{\mathrm{thr},a}^{\alpha,S\\!N\\!R}}\right|\right)\\!\\!\right)\\!.$
(1)

Computing utilities for $a\\!B\\!S$ $a$ needs $\mathcal{O}(U_{\max})$ sums and
powers (logarithms for $\alpha\\!=\\!1$). Hence, finding the maximum shown in
Eq. (1) has complexity $\mathcal{O}(|\mathcal{A}|\\!\cdot\\!U_{\max})$ powers
(or logarithms), the sums incurring negligible extra complexity.

Having identified $a_{0}$, PADD selects a new random position within the
allowed 3D ball space around the current position of the drone.

###  5.6 Complexity of PADD

In this section we provide the complexity of PADD. The algorithm consists of
sequential steps, involving: user association with BSs; the backhaul
association solved with the GAP problem of Fig. 5; the resource allocation of
each $g\\!B\\!S$ solved with the convex optimization of Fig. 6; the system
utility evaluation with the $\alpha$-fairness metric; and the least fit drone
selection of the EO operation. The overall PADD’s complexity is
$\mathcal{O}\left(N\left(\frac{1}{2\log\\!2}|\mathcal{B}|\log|\mathcal{B}|+2A_{g}|\mathcal{G}||\mathcal{A}|+4A_{g}U_{\max}+|\mathcal{U}||\mathcal{B}|\right)\right)$
simple operations (e.g., sums or comparisons) and
$\mathcal{O}\left(N|\mathcal{B}|+4|\mathcal{A}|U_{\max}\right)$ complex
operations (e.g., powers or logarithms). This means that PADD’s complexity is
low-order polynomial, linear with respect to the number of mobile users and
log-linear with respect to the number of BSs. More details are available in
Section B in the online supplemental material.

Moreover, PADD’s iterations converge very fast (we have observed no more than
a hundred iterations) and our algorithm implemented on Matlab took at most
about 800 billion floating-point operations to optimize drone positions, which
requires, e.g., about 10 s on a linux machine with an Intel Xeon E5-2670
processor at $2.6$ GHz with a processor performing 32 flops/cycle with 64-bit
operations.

##  6 Numerical Simulations

Here we present the numerical evaluation of the $\alpha$-fair cellular
capacity optimization in both static and dynamic scenarios. All our
simulations are performed over the real topology of an operational network
deployed in a dense city: 150,000 inhabitants in a 10 km2 area. (Leganés,
south of Madrid, Spain). The area is covered by 10 $g\\!B\\!S$s using the same
LTE band, as shown in Fig. 7.

As there are several operators in Spain providing LTE service over multiple
bands, and we only consider one operator in one band, we use a number of 1000
UEs that request service at the same time, unless otherwise specified.

First, we numerically validate our proposed PADD in comparison with optimal
results approximated by means of Monte-Carlo (MC) simulations in static
networks with limited population: We perform $10^{7}$ MC runs per instance
where $a\\!B\\!S$s are placed in feasible random positions. Then, we optimize
the network for each run and save the best setting as the approximated optima.
Second, we analyze the robustness of the system model in order to prove that
assuming perfect knowledge about user positions has limited impact on
performance. Eventually, we study three significative static and dynamic
scenarios:

  * •

Poisson Point Process (PPP): we statically place UEs on the city map according
to a Poisson point process [45].

  * •

Stadium: we statically place 60% of the UEs in the surroundings of a stadium,
and the rest like in PPP.

  * •

Event: 40% of the population of UEs moves according to the random way-point
model, whereas other UEs keep arriving at an official scheduled rate of a
train station and move towards the stadium.

For each scenario, we study drone placement, utility and throughput achieved,
and system fairness measured with the Jain’s index [46].

Table III reports the evaluation parameters used in our simulations. We take a
carrier bandwidth of $20$ MHz for both $g\\!B\\!S$ and $a\\!B\\!S$ channels
(out of which, we consider that $10\%$ is for guard bands, so we only use $18$
MHz [33]). The carrier frequency of cellular links is $1815.1$ MHz, while for
air-to-ground links we use $2630$ MHz. These are two commonly used LTE bands.
The transmission power of $a\\!B\\!S$s is 25 dBm, notably lower than the 44
dBm power transmission of the ground $g\\!B\\!S$s. In addition, probed aerial
positions arise from a lattice that spans equal-volume subspaces.

To position our proposal, we consider the performance without drones as
baseline (Ground in the figures). Moreover, to compare PADD to state of the
art drone-position optimization frameworks, we have implemented and tested the
Repulsion-Attraction scheme (RA) [47].

Figure 7: Topology of Leganés (Spain) and $g\\!B\\!S$s placement.

Although our analysis holds for generic values of $\alpha$, we present results
for three specific and most interesting cases:

  1. (i)

For $\alpha=0$, we obtain the maximum throughput achievable (MaxThr in the
figures).

  2. (ii)

For $\alpha=1$, we optimize the proportional fairness network metric (PropFair
in the figures).

  3. (iii)

For $\alpha\rightarrow\infty$, we optimize the max-min fairness network metric
(MaxMin in the figures).

We study these three metrics for the following reasons. The MaxThr metric
provides the maximum achievable network capacity, without fairness. The
PropFair metric takes into account the network capacity but it also does not
let fairness decay. The mathematical design of this metric, which optimizes
the aggregation of logarithms of user throughputs, allows for finding a trade-
off between a high system capacity and high fairness values. The MaxMin metric
only targets fairness of the weakest customer, which comes at the cost of
providing lower aggregate throughput. Both the PropFair and MaxMin have been
adopted in the implementation of real telecommunication systems [48] as well
as in many research works [49, 50]. Complementary, in Section D in the online
supplemental material, we provide more theoretical results for the PADD scheme
for a generic value of $\alpha$.

We have simulated every analyzed use-case $1000$ times using Matlab R2020a and
show average results. Error bars in the figures are 95% confidence intervals.

TABLE III: Evaluation Parameters Parameter | Value
---|---
$\xi_{LoS}$, $\xi_{NLoS}$, $\beta_{1}$, $\beta_{2}$ | 1.6 dB, 23 dB, 12.08, 0.11
Carrier frequencies, $f_{\mathcal{G}}$, $f_{\mathcal{A}}$ |  $1815.1$ MHz, $2.63$ GHz
Bandwidths, $W_{\mathcal{G}}$, $W_{\\!\mathcal{A}}$ |  $18$ MHz, $18$ MHz
Tx power, $P_{Tx}^{g}$, $P_{Tx}^{a}$ |  $44$ dBm, $25$ dBm
Thermal Noise Power | -174 dBm/Hz
Ground path loss exponent, $\eta_{BS}$ | $3$
Height range, $[h_{\min},h_{\max}]$ |  $[40,300]$ m
Urban area, $|\mathcal{S}|$ |  $10$ km2
Average walking speed |  $2.5$ m/s
Monte-Carlo runs per instance | $10^{7}$
Instances of simulations | $1000$

Figure 8: Utility validation for $\alpha~{}\in~{}\\{0,1,\infty\\}$. $G=10$,
$U=1000$. Scenario: PPP.

Figure 9: Network capacity validation for $\alpha\in\\{0,1,\infty\\}$. $G=10$,
$U=1000$. Scenario: PPP.

Figure 10: Robustness validation for $\alpha=1$. $a\\!B\\!S$s placement error.
$G=10$, $A=5$, $U=1000$. Scenario: PPP.

###  6.1 Validation of PADD Operation

Here we compare PADD results and optima (approximated by means of MC, since
there are no computationally feasible alternatives) in the PPP scenario, and
comment on the basic properties of our approach.

In Fig. 10, we present a utility comparison between the Ground scheme, optima
results from MC and our proposal PADD with the three main metrics studied in
this paper (MaxThr, PropFair and MaxMin). Note that different values of
$\alpha$ make utility values very different, so we cannot compare utilities
across fairness metrics. So, we will compare PADD and the existing schemes on
a per-metric basis.

We observe that the difference between utilities achieved with PADD and MC is
negligible, below $1\%$ in all cases. This shows that our proposal PADD is
able to achieve close-to-optimal results with a much lower complexity.
Furthermore, PADD matches very well also the throughput achievable in the
optimal case, as illustrated in Fig. 10. This shows that PADD is able to
pursue the optimal system configuration, not just an operational configuration
that is near-optimal according to the chosen utility metric.

The figures show that, in the PPP scenario, utility and throughput increase
with fleet size. This is because user’s density is homogeneously spread over
the entire area, so that drones are spaced apart, and inter-drone interference
constraints do not kick in. Note also that the gain in terms of utility and
throughput with respect to the Ground configuration is remarkable under all
considered metrics, which confirms that coordinated drone relays have huge
potential.

###  6.2 Robustness of PADD

Figure 11: Utility of all users for $\alpha\in\\{0,1,\infty\\}$. $G=10$,
$U=1000$. Scenario: Stadium with $U_{d}=600$.

Figure 12: Utility of stadium users for $\alpha\in\\{0,1,\infty\\}$. $G=10$,
$U=1000$. Scenario: Stadium with $U_{d}=600$.

Figure 13: Throughput of stadium users for $\alpha\in\\{0,1,\infty\\}$.
$G=10$, $U=1000$. Scenario: Stadium with $U_{d}=600$.

In the system model described in Section 3, we have assumed that the positions
of users are known. However, it might be not realistic to estimate user
locations with negligible error unless GPS is enabled or many drones are
available [51]. For instance, using a trilateration on the signal strength at
the base station, the error on the position is normally below $50$ m [52].
Hence, here we consider the PPP scenario to numerically analyze the robustness
of PADD by introducing uncertainty in the position of the ground users,
uniformly at random, within $50$ m.

In Fig. 10, we observe that the position of drones selected by PADD is similar
with and without localization errors. Indeed, we have quantified an average
utility relative loss below $5\%$ (see Section C in the online supplemental
material for further details). The reason of such robust behavior stands in
the fact that the optimization of drone positions is done based on many users
and in relatively large areas, so that multiple errors in user positions are
not so important, whereas the presence of a distributed mass of users in a
given area is what actually catalyzes the presence of a drone.

Next, we analyze non-homogeneous user topologies in static and dynamic cases.

###  6.3 Performance Evaluation in the Static Stadium Case

The Stadium scenario allows to study network performance when the ground
network cannot sustain the UE’s demand.

Figure 14: Network state at $t=25$ min. $G=10$, $A=5$, $U=640$. Scenario:
PropFair, Event.

Figure 15: Network state at $t=60$ min. $G=10$, $A=5$, $U=880$. Scenario:
PropFair, Event.

Figure 16: Drone trajectories during 75 minutes. $G\\!=\\!10$, $A\\!=\\!5$,
$U\\!=\\![400,$ $\ldots,1000]$. Scenario: PropFair, Event.

Figure 17: Attendance utility for $\alpha\in\\{0,1,\infty\\}$. $G=10$, $A=5$,
$U=[400,\ldots,1000]$. Scenario: Event.

In Fig. 13 we show the average performance experienced by all users in the
scenario, while in Fig. 13 we focus on the performance of users by the
stadium. We observe that PADD with MaxThr benefits of the presence of drones
(see Fig. 13, top), although adding drones is negative for users by the
stadium (see Fig. 13, top). In fact, adding $a\\!B\\!S$s introduces more
capacity and connectivity opportunities, hence the benefit at aggregate city
level. However, PADD with MaxThr seeks for drone positions where they are less
impaired by interference, so to generate a few good-quality channels. This
behavior results into having only one or two $a\\!B\\!S$s at most located near
the stadium. More $a\\!B\\!S$s would interfere too much. Hence, as long as
$a\\!B\\!S$s are added, PADD with MaxThr positions them apart from the
stadium, yet they generate some interference, which progressively worsens the
performance of users by the stadium.

PADD with PropFair behaves completely different. It brings drones where it
favors users otherwise served below the average, so to increase the system’s
log utility. This results in deploying $a\\!B\\!S$s by the stadium, where the
density of users in presence of limited radio resources hinders performance
more than radio channel quality issues.

With MaxMin, PADD positions drones where users would otherwise suffer the
worst connection quality, irrespective of their closeness to densely populated
areas. Therefore the aggregated throughput is lower than with PropFair, which
in turn is much lower than with MaxThr.

If we now consider only users by the stadium, Fig. 13 illustrates how using
PropFair can clearly outperform MaxThr in terms of throughput. This is due to
the fact that the optimization of throughput requires positioning drones where
they interfere less, which is not necessarily by the stadium. Indeed, due to
interference, users by the stadium could incur a loss by increasing the number
of drones even in the case of using MaxMin, as shown in the figure.

As a result of the previous numerical analysis, we observe that using
$a\\!B\\!S$s can be beneficial to help the ground cellular network in dense
spots, except it cannot help in purely maximizing throughput (e.g., with
MaxThr). Drones cannot either “rescue” all users with bad channel conditions,
as PADD would seek with MaxMin and PropFair. However, PADD can always provide
fairness and large throughput gain with respect to the Ground case. Besides,
the version of PADD with PropFair results to be quite effective in case of
dense masses of users.

###  6.4 Performance Evaluation in the Dynamic Event Case

The last scenario we consider is dynamic and allows us to study the evolution
of network performance while the density of users increases. Moreover, it
shows the importance of designing a fast and reactive algorithm to re-position
drones as the user topology changes.

In the Event scenario, small masses of 40 users arrive periodically to the
train station of the city every 5 minutes. The initial population is 400 users
and keeps growing during 75 minutes up to 1000 users. There are 5 $a\\!B\\!S$s
hovering the area in this example. Upon a train arrival, the new mobile users
walk towards the stadium, located 1.5 km away from the train station. The
fleet is repositioned every 5 simulated minutes, using as initial condition
the positions of the 5 drones in the previous optimization epoch.

We first illustrate how the network evolves over time in Figs. 16 and 16,
using PADD with PropFair. After 25 minutes (Fig. 16) some people are already
at the stadium, while many others keep arriving and are walk towards it. At
that point in time, 3 $a\\!B\\!S$s are getting prepared to serve the users
nearby the stadium and also the smaller groups on their way from the train
station to the stadium. We see how drones adapt their positions after other 35
minutes, in Fig. 16, when much more users have reached the stadium. By that
time, one more drone has been dispatched as well, to serve the stadium. The
trajectories of the drones for a 75-minute simulation instance are shown in
Fig. 16, where the red drones represent the source position of $a\\!B\\!S$s.
In the figure, we can see that $a\\!B\\!S$ 5 is not required to assist the
dense stadium spot, while $a\\!B\\!S$s 1 and 2 are always hovering between the
train station and the stadium. Also, $a\\!B\\!S$s 3 and 4 keep moving back and
forth within different regions, to fairly supply the demand of users.

As concerns performance, Fig. 17 illustrates utility as they evolve over time.
Here, in addition to PADD with the three selected $\alpha$’s, we also compare
the Repulsion-Attraction (RA) scheme. With RA, $a\\!B\\!S$s are attracted by
UE’s inverse SNR, and repulsed by proximity to $g\\!B\\!S$s to avoid
interference. RA does not target any fairness metric, so we quantify its
impact with the same utility functions used for PADD computed for users
arriving over time (the attendance).

As time passes by, and more people reach the stadium, Fig. 17 shows a
significant utility raise under the adoption of either MaxThr or PropFair
schemes. With MaxThr, the gain of PADD over RA and Ground schemes is high,
although it saturates quickly. Instead, with PropFair, PADD exhibits a
smoother behavior, as its gain keeps increasing.

Under the MaxMin scheme, PADD performs better than RA, although now we observe
a decay of performance over time, for all schemes. This is due to the fact
that, with more attendance, the minimum per-user achieved rate will have
decrease, unless more drones where deployed.

Clearly, the RA scheme is not able to opportunistically take advantage of
user’s diversity and improve utility because it does not target a throughput-
based metric, unlike PADD. This permits PADD to optimize the network with much
better guarantees in terms of throughput and fairness.

##  7 Lessons learnt and discussion

The performance assessment carried in this paper shows the importance of
integrating a fleet of drone relay stations in a cellular network. It also
unveils that optimizing drone positions to maximize throughput, without taking
into account fairness, has little relevance in the presence of dense spots or
ground users. Instead, a fair metric like PropFair provides notable throughput
and utility improvements.

The fact that PADD is fast, allows to design an almost continuous
reconfiguration of positions in realistic networks, as shown in this paper. In
turn, our scheme is fast because we have designed it by segmenting the problem
into a few phases: we use EO for optimizing drone positions, while wireless
backhaul attachment, BS selection and distributed resource allocation are
sequentially solved optimally and analytically for each drone deployment
topology.

The actual implementation of PADD requires the exchange of signaling
information between drones and a centralized orchestrator, as well as the
implementation of a mechanism to track users. Signaling incurs some limited
overhead to instruct drones and to gather user positions and interference
reports (which are however already collected by current BSs), depending on the
frequency of reconfigurations. However, we have seen that PADD is robust to
imprecise tracking of user’s positions, which can be therefore strongly
simplified, so that the additional overhead due to PADD will mostly be due to
controlling drones.

Finally, we comment on drone mechanicals. Current commercial drones can carry
small BSs and access points, although they cannot fly for long time due to
battery limitations (around 30 minutes at most). They can easily move at the
reasonable speed of 15 m/s. Therefore, it is possible to derive a
repositioning scheme that accounts for replacing drones that go back to the
charging station if the drones do not fly too far from it. For the case
considered in Section 6.4, the routes flown by drones in 5 minutes are short
enough, and it is easy to hover a city district in a few minutes. Thus,
notwithstanding the intricacies of the analysis, the performance evaluation
discussed in this paper is quite relevant for realistic systems.

##  8 Conclusions

In this paper we have proposed an analytic framework to optimize drone-aided
cellular networks in terms of an $\alpha$-fair throughput utility function
under realistic stochastic models. Specifically, we have studied the
integration of a coordinated fleet of aerial base stations carried by drones
and relaying traffic for ground base stations by means of 3D-beamforming
wireless backhaul connections. Due to the complexity of the studied problem,
we have proposed PADD, an optimization algorithm based on extremal
optimization that works in parallel threads. PADD provides near-optimal
solutions in low-degree polynomial time, with a linear dependency on all
parameters but for the number of base stations, which causes a sub-quadratic
dependency. This makes PADD suitable for implementation in dynamically
changing environments. The performance evaluation presented in the paper shows
that PADD brings significant gain and outperforms existing approaches. It also
unveils that using fairness is key to get benefit from coordinated yet
interfering drone relay stations.

## References

  * [1] E. Arribas, V. Mancuso, and V. Cholvi, “Fair Cellular Throughput Optimization with the Aid of Coordinated Drones,” in _IEEE INFOCOM WKSHPS – MiSARN_ , pp. 1–6, 2019. Link to paper.
  * [2] W. Mohr, “5G empowering vertical industries,” Tech. Rep., 2016.
  * [3] M. Maternia _et al._ , “5G PPP use cases and performance evaluation models,” 5G-PPP, Tech. Rep., Apr. 2016, White Paper.
  * [4] C. Liu et al., “The case for re-configurable backhaul in cloud-RAN based small cell networks,” in _IEEE INFOCOM_ , 2013.
  * [5] V. Genc, S. Murphy, and J. Murphy, “Performance analysis of transparent relays in 802.16j mmr networks,” in _WiOPT_ , 2008.
  * [6] I. Bor-Yaliniz and H. Yanikomeroglu, “The new frontier in RAN heterogeneity: Multi-tier drone-cells,” _IEEE Communications Magazine_ , vol. 54, no. 11, 2016.
  * [7] W. Guo, C. Devine, and S. Wang, “Performance analysis of micro unmanned airborne communication relays for cellular networks,” in _IEEE CSNDSP_. IEEE, 2014.
  * [8] W. Guo and T. O’Farrell, “Relay deployment in cellular networks: Planning and optimization,” _IEEE JSAC_ , vol. 31, no. 8, 2013.
  * [9] M. Mozaffari, W. Saad, M. Bennis, and M. Debbah, “Wireless communication using unmanned aerial vehicles (UAVs): Optimal transport theory for hover time optimization,” _IEEE Transactions on Wireless Communications_ , vol. 16, no. 12, pp. 8052–8066, 2017.
  * [10] J. Mo and J. Walrand, “Fair end-to-end window-based congestion control,” _IEEE/ACM Transactions on networking_ , vol. 8, no. 5, pp. 556–567, 2000\.
  * [11] S. Boettcher and A. Percus, “Optimization with extremal dynamics,” _Phys. Rev. Lett._ , vol. 86, pp. 5211–5214, 2001.
  * [12] R. Ferrús et al. , “SDN/NFV-enabled satellite communications networks: Opportunities, scenarios and challenges,” _Physical Communication_ , vol. 18, pp. 95–112, 2016.
  * [13] E. Arribas and V. Mancuso, “Multi-path D2D leads to satisfaction,” in _WoWMoM_. IEEE, 2017, pp. 1–7.
  * [14] ——, “Achieving per-flow satisfaction with multi-path D2D,” _Ad Hoc Networks_ , vol. 106, p. 102162, 2020.
  * [15] P. S. Bithas, V. Nikolaidis, A. G. Kanatas, and G. K. Karagiannidis, “Uav-to-ground communications: Channel modeling and uav selection,” _IEEE Trans. Commun._ , vol. 68, no. 8, pp. 5135–5144, 2020.
  * [16] M. Mozaffari, W. Saad, M. Bennis, Y. Nam, and M. Debbah, “A tutorial on uavs for wireless networks: Applications, challenges, and open problems,” _IEEE Commun. Surv. Tutorials_ , vol. 21, no. 3, pp. 2334–2360, 2019.
  * [17] I. Strumberger, M. Sarac, D. Markovic, and N. Bacanin, “Moth search algorithm for drone placement problem,” _International Journal of Computers_ , vol. 3, 2018.
  * [18] M. Mozaffari, W. Saad, M. Bennis, and M. Debbah, “Unmanned aerial vehicle with underlaid device-to-device communications: Performance and tradeoffs,” _IEEE Transactions on Wireless Communications_ , vol. 15, no. 6, pp. 3949–3963, 2016.
  * [19] E. Arribas, V. Mancuso, and V. Cholvi, “Coverage Optimization with a Dynamic Network of Drone Relays,” _IEEE Transactions on Mobile Computing_ , vol. 19, no. 10, pp. 2278–2298, 2020.
  * [20] N. Kumar, M. Ghosh, and C. Singhal, “UAV network for surveillance of inaccessible regions with zero blind spots,” in _IEEE INFOCOM Workshops_ , 2020, pp. 1213–1218.
  * [21] X. Liu, Y. Liu, and Y. Chen, “Reinforcement learning in multiple-uav networks: Deployment and movement design,” _IEEE Trans. Veh. Technol._ , vol. 68, no. 8, pp. 8036–8049, 2019.
  * [22] J. Cui, Y. Liu, Z. Ding, P. Fan, and A. Nallanathan, “Qoe-based resource allocation for multi-cell NOMA networks,” _IEEE Trans. Wirel. Commun._ , vol. 17, no. 9, pp. 6160–6176, 2018.
  * [23] B. Galkin, J. Kibilda, and L. DaSilva, “A stochastic model for UAV networks positioned above demand hotspots in urban environments,” _IEEE Trans. Veh. Technol._ , vol. 68, no. 7, pp. 6985–6996, 2019\.
  * [24] P. Jacquet, D. Popescu, and B. Mans, “Connecting flying backhauls of drones to enhance vehicular networks with fixed 5g NR infrastructure,” in _IEEE INFOCOM Workshops_ , 2020, pp. 472–477.
  * [25] A. Alsharoa, H. Ghazzai, A. Kadri, and A. E. Kamal, “Spatial and temporal management of cellular hetnets with multiple solar powered drones,” _IEEE Transactions on Mobile Computing_ , vol. 19, no. 4, pp. 954–968, 2019\.
  * [26] B. Danila, Y. Yu, J. Marsh, and K. Bassler, “Optimal transport on complex networks,” _Physical Review E_ , vol. 74, no. 4.
  * [27] A. BenMimoune and M. Kadoch, _Relay Technology for 5G Networks and IoT Applications_. Cham: Springer International Publishing, 2017, pp. 3–26.
  * [28] Y. Zeng, J. Lyu, and R. Zhang, “Cellular-connected UAV: potential, challenges, and promising technologies,” _IEEE Wirel. Commun._ , vol. 26, no. 1, pp. 120–127, 2019.
  * [29] E. Arribas, V. Cholvi, V. Mancuso _et al._ , “An optimal scheme to recharge communication drones,” in _IEEE GLOBECOM_ , 2021.
  * [30] A. Al-Hourani, S. Kandeepan, and A. Jamalipour, “Modeling air-to-ground path loss for low altitude platforms in urban environments,” in _GLOBECOM_. IEEE, 2014, pp. 2898–2904.
  * [31] A. Al-Hourani, S. Kandeepan, and S. Lardner, “Optimal LAP altitude for maximum coverage,” _IEEE Wireless Communications Letters_ , vol. 3, no. 6, pp. 569–572, 2014.
  * [32] J. Goldhirsh and W. J. Vogel, “Handbook of propagation effects for vehicular and personal mobile satellite systems,” _NASA Reference Publication_ , vol. 1274, pp. 40–67, 1998.
  * [33] S. Sesia, M. Baker, and I. Toufik, _LTE-the UMTS long term evolution: from theory to practice_. John Wiley & Sons, 2011.
  * [34] G. Das, R. Fraser, A. Lóopez-Ortiz, and B. Nickerson, “On the discrete unit disk cover problem,” _International Journal of Comp. Geometry & Applications_, vol. 22, no. 05, pp. 407–419, 2012.
  * [35] S. Burer and A. N. Letchford, “Non-convex mixed-integer nonlinear programming: A survey,” _Surveys in Operations Research and Management Science_ , vol. 17, no. 2, pp. 97–106, 2012.
  * [36] M. L. Fisher, R. Jaikumar, and L. N. Van Wassenhove, “A multiplier adjustment method for the generalized assignment problem,” _Management Science_ , vol. 32, no. 9, pp. 1095–1103, 1986.
  * [37] R. Cohen, L. Katzir, and D. Raz, “An efficient approximation for the generalized assignment problem,” _Information Processing Letters_ , vol. 100, no. 4, pp. 162–166, 2006.
  * [38] W. Sun and Y.-X. Yuan, _Optimization theory and methods: non linear programming_. Springer Sci. & Business Media, 2006, vol. 1.
  * [39] D. G. Cattrysse and L. N. Van Wassenhove, “A survey of algorithms for the generalized assignment problem,” _European journal of operational research_ , vol. 60, no. 3, pp. 260–272, 1992.
  * [40] A. Fréville, “The multidimensional 0–1 knapsack problem: An overview,” _European Journal of Operational Research_ , vol. 155, no. 1, pp. 1–21, 2004\.
  * [41] R. Andonov, V. Poirriez, and S. Rajopadhye, “Unbounded knapsack problem: Dynamic programming revisited,” _European Journal of Operational Research_ , vol. 123, no. 2, pp. 394–407, 2000.
  * [42] Y. Nesterov and A. Nemirovskii, _Interior-point polynomial algorithms in convex programming_. Siam, 1994, vol. 13.
  * [43] S. Bubeck _et al._ , “Convex optimization: Algorithms and complexity,” _Foundations and Trends® in Machine Learning_ , vol. 8, no. 3-4, pp. 231–357, 2015.
  * [44] H. Kuhn and A. Tucker, “Nonlinear programming,” in _Traces and emergence of nonlin. programming_. Springer, 2014, pp. 247–258.
  * [45] G. Last and M. Penrose, _Lectures on the Poisson process_. Cambridge University Press, 2017, vol. 7.
  * [46] R. K. Jain, D.-M. W. Chiu, and W. R. Hawe, “A quantitative measure of fairness and discrimination,” _Eastern Research Laboratory, Digital Equipment Corporation, Hudson, MA_ , 1984.
  * [47] O. Andryeyev and A. Mitschele-Thiel, “Increasing the cellular network capacity using self-organized aerial base stations,” in _DroNet_. ACM, 2017, pp. 37–42.
  * [48] G. Miao, J. Zander, K. W. Sung, and S. B. Slimane, _Fundamentals of mobile data networks_. Cambr. Univ. Press, 2016.
  * [49] S. Mosleh, L. Liu, and J. Zhang, “Proportional-fair resource allocation for coordinated multi-point transmission in LTE-advanced,” _IEEE Trans. on Wireless Comm._ , vol. 15, no. 8, pp. 5355–5367, 2016.
  * [50] R. Li and P. Patras, “Max-Min Fair Resource Allocation in Millimetre-Wave Backhauls,” _IEEE Trans. on Mob. Comput._ , 2019.
  * [51] R. Petrolo, Y. Lin, and E. Knightly, “ASTRO: Autonomous, sensing, and tetherless networked drones,” in _DroNet_. ACM, 2018, pp. 1–6.
  * [52] J. del Peral-Rosado et al., “Survey of cellular mobile radio localization methods: from 1G to 5G,” _IEEE Communications Surveys & Tutorials_, vol. 20, no. 2, pp. 1124–1148, 2017.
  * [53] ITU-R, “Propagation data and prediction methods required for the design of terrestrial line-of-sight systems,” _ITU-R_ , 2015.

|  Edgar Arribas obtained his PhD in Telematic Engineering in 2020 at IMDEA
Networks Institute and Universidad Carlos III de Madrid (UC3M), funded by the
MECD FPU15/02051 grant. He is currently a lecturer and researcher at the
Applied Mathematics and Statistics Department of Universidad CEU San Pablo
(Spain). He works on optimization of dynamic relay in wireless networks.
---|---
|  Vincenzo Mancuso is Research Associate Professor at IMDEA Networks, Madrid,
Spain, and recipient of a Ramon y Cajal research grant of the Spanish Ministry
of Science and Innovation. Previously, he was with INRIA (France), Rice
University (USA) and University of Palermo (Italy), from where he obtained his
Ph.D. in 2005. His research focus is on analysis, design, and experimental
evaluation of opportunistic wireless architectures and mobile broadband
services.
---|---
|  Vicent Cholvi graduated in Physics from the University of Valencia, Spain
and received his doctorate in Computer Science in 1994 from the Polytechnic
University of Valencia. In 1995, he joined the Jaume I University in
Castellón, Spain where he is currently a Professor. His interests are in
distributed and communication systems.
---|---

## Online Supplemental Material

##  A Channel Modelling

We assume that the network operator disposes of two orthogonal frequency
bands. One band is assigned to $g\\!B\\!S$s to provide access service to
ground users as well as aerial backhaul service to $a\\!B\\!S$s. The other
band is assigned to $a\\!B\\!S$s for aerial user access. Hence, we model three
different channel types: (i) air-to-ground and (ii) ground-to-ground channels
in the access network, and (iii) ground-to-air channels in the backhaul
network.

Indeed, the access network communication channels between BSs and UEs differ
much depending on whether users connect to a $g\\!B\\!S$ or to an $a\\!B\\!S$.
While the ground-to-ground channel attenuation for $g\\!B\\!S$–UE links
follows conventional path-loss modeling based on slow and fast fading, air-to-
ground channels ($a\\!B\\!S$–UE links) suffer additional attenuation depending
on the LoS—or NLoS—state of the channel. Such additional attenuation is
referred to in the literature as an excess attenuation [30]. Moreover,
antennas used for the access network differ from backhaul network antennas
performing 3D-beamforming, which directly affects the interference suffered in
each case. In the following sections, we detail these features for each type
of modelled channel.

###  A.1 Air-to-Ground Channels

Depending on whether links between $a\\!B\\!S$s and UEs are free of obstacles
(e.g., buildings, traffic, etc.), the attenuation differs notably [30]. The
LoS-likelihood is a complex function of the elevation angle between UE
$u\in\mathcal{U}$ and $a\\!B\\!S$ $a\in\mathcal{A}$:

$\displaystyle
P_{LoS}(a,u)=\frac{1}{1+\beta_{1}\cdot\exp\left(-\beta_{2}\left(\frac{180}{\pi}\arctan\left(\frac{h_{a}}{r_{a,u}}\right)-\beta_{1}\right)\right)},$
(2)

where the elevation of $a$ is $h_{a}$, while $\beta_{1}$ and $\beta_{2}$ are
parameters depending on the number of large signal obstructions per unit area,
building’s height distribution, ratio of built-up area and clean surfaces,
etc., as derived in [31], based on ITU recommendations [53]. In Eq. (2),
$\theta_{a,u}\\!=\\!\arctan(h_{a}/r_{a,u})$ is the elevation angle.
$\theta_{a,u}$ approaches $\frac{\pi}{2}$ when the $a\\!B\\!S$ $a$ hovers just
above the user $u$, i.e., when the LoS likelihood reaches its maximum. The
elevation angle $\theta_{a,u}$ is characterized by the $a\\!B\\!S$ height and
the ground distance between the user and the $a\\!B\\!S$, that is
$r_{a,u}=\|(X_{a},Y_{a})-(x_{u},y_{u})\|$. Since safety, legislation and
technology constraints allow drones to fly only at low altitudes—a few
hundreds of meters at most—the ground distance can be of the same order of the
drone’s elevation, so that the argument of the $\arctan$ function in Eq. (2)
can vary sensibly with user and drone positions, and so can LoS likelihood.
This means that drone positions, including their elevation, play an important
role on air-to-ground channel conditions.

In particular, the average attenuation (in dB units) of an air-to-ground
channel between drone $a$ and user $u$ depends on the LoS likelihood, with the
following expression [31]:

$\displaystyle L_{\mathcal{A}}(a,u)$ $\displaystyle=$ $\displaystyle
20\log_{10}\left(\frac{4\pi
f_{\mathcal{A}}}{c}\cdot\sqrt{h_{a}^{2}+r_{a,u}^{2}}\right)+$ (3)
$\displaystyle+$ $\displaystyle
P_{LoS}(a,u)\cdot\left(\xi_{LoS}-\xi_{N\\!LoS}\right)+\xi_{N\\!LoS},$

where $\xi_{LoS},\xi_{N\\!LoS}$ are the excess attenuation components in
LoS/NLoS conditions; $f_{\\!\mathcal{A}}$ is the carrier frequency in Hz; and
$c$ is the speed of light in m/s.

Since $g\\!B\\!S$s and $a\\!B\\!S$s operate onto orthogonal bands, there is no
interference between drone-served users and cellular users, which is commonly
the main limiting factor in aided cellular networks, as, e.g., inband D2D
networks [13]. With the above, the experienced SINR for air-to-ground access
links $(a,u)$ is:

$\displaystyle\gamma^{\mathcal{A}}_{a,u}=\frac{P^{a}_{Tx}\\!\cdot\\!10^{-L_{\mathcal{A}}(a,u)/10}}{N_{a,u}+I_{a,u}^{\mathcal{A}}},$
(4)

where $P^{a}_{Tx}$ is the transmission power of an omnidirectional antenna in
the $a\\!B\\!S$s $a\in\mathcal{A}$; $N_{a,u}$ is thermal noise according to
the allocated bandwidth; and $I_{a,u}^{\mathcal{A}}$ is the interference level
that user $u$ suffers from other $a\\!B\\!S$s. However, note that the 3D
position of an $a\\!B\\!S$ is a decision parameter that directly affects
interfering signals received by user $u$, i.e.:

$\displaystyle
I_{a,u}^{\mathcal{A}}=\sum\limits_{a^{\prime}\in\mathcal{A}\setminus\\{a\\}}P^{a}_{Tx}\\!\cdot\\!10^{-L_{\mathcal{A}}(a^{\prime},u)/10},\quad\forall
a\in\mathcal{A},$ (5)

where $L_{\mathcal{A}}(a^{\prime},u)$ depends on the 3D position of
$a\\!B\\!S$s $a^{\prime}\\!\in\\!\mathcal{A}$, as shown in Eq. (3).

###  A.2 Ground-to-Ground Channels.

Connections in the access network between $g\\!B\\!S$s and users experience an
attenuation based on a well known path-loss model with slow fading (in dB
units):

$\displaystyle
L_{\mathcal{G}}(g,u)=10\eta_{\mathcal{G}}\log_{10}\left(\frac{4\pi
f_{\mathcal{G}}}{c_{l}}\cdot\|\Pi^{g}-\pi_{u}\|\right)+\mathcal{N}(0,\sigma_{\mathcal{G}}^{2}),$
(6)

where $\eta_{\mathcal{G}}>2$ is the path-loss exponent in ground
communications; $f_{\mathcal{G}}$ is the operating carrier frequency of the
$g\\!B\\!S$s; and $\sigma_{\mathcal{G}}$ is the standard deviation of the
gaussian random variable $\mathcal{N}(0,\sigma^{2}_{\mathcal{G}})$, modelling
the effects of shadowing.

As mentioned above, since there is no interference between cellular users and
drone-served users, the SINR for access links $(g,u)$ is:

$\displaystyle\gamma^{\mathcal{G}}_{g,u}=\frac{P^{g}_{Tx}\\!\cdot\\!10^{-L_{\mathcal{G}}(g,u)/10}}{N_{g,u}+I_{g,u}^{\mathcal{G}}},$
(7)

where $P^{g}_{Tx}$ is the transmission power of an omnidirectional antenna
integrated in the $g\\!B\\!S$s $g\\!\in\\!\mathcal{G}$; $N_{g,u}$ represents
thermal noise according to the allocated bandwidth; and most importantly,
$I_{g,u}^{\mathcal{G}}$ is the interference level that user $u$ suffers from
other $g\\!B\\!S$s.

###  A.3 Ground-to-Air Channels

The aerial network relays traffic from the $g\\!B\\!S$s by means of LoS
backhaul wireless links. Hence, the attenuation of a $g\\!B\\!S$–$a\\!B\\!S$
link $(g,a)$ is the following:

$\displaystyle
L_{\mathcal{B}}(g,a)=10\eta_{\mathcal{B}}\log_{10}\left(\frac{4\pi
f_{\mathcal{B}}}{c}\\!\cdot\\!\left\|\Pi^{g}\\!-\\!\Pi_{a}\right\|\right)\\!+\\!\mathcal{N}\left(0,\sigma_{\mathcal{B}}^{2}\right),$
(8)

where $\eta_{\mathcal{B}}\approx 2$ is the path-loss exponent in LoS;
$f_{\mathcal{B}}$ is the operating carrier frequency of the backhaul wireless
links; and $\sigma_{\mathcal{B}}^{2}$ is the standard deviation of the
gaussian random variable $\mathcal{N}\left(0,\sigma_{\mathcal{B}}^{2}\right)$,
modeling the effects of shadowing.

Backhaul links operate on the bandwidth shared with user access to
$g\\!B\\!S$s. However, as backhaul links perform 3D-beamforming pointing to
the air (where $a\\!B\\!S$s hover), the interference between
$g\\!B\\!S$-served users and backhaul-served $a\\!B\\!S$s is very limited.
Although the majority of the $g\\!B\\!S$ radiating power is focused in one
direction towards the air thanks to the adoption of 3D-beamforming, non-ideal
beam-patterns also radiate energy in other directions. Therefore, the SINR
experienced by an $a\\!B\\!S$ $a\\!\in\\!\mathcal{A}$ depends also on the
direction in which other $g\\!B\\!S$s transmit to other $a\\!B\\!S$s. The SINR
experienced by a $g\\!B\\!S$–$a\\!B\\!S$ link $(g,a)$ is:

$\displaystyle\gamma^{\mathcal{B}}_{g,a}=\frac{P_{Tx}^{g}\\!\cdot\\!G_{g}\\!\cdot\\!10^{-L_{\mathcal{B}}(g,a)/10}}{N^{g,a}+I^{\mathcal{B}}_{g,a}},$
(9)

where $P_{Tx}^{g}$ is the transmission power of the $g\\!B\\!S$ $g$; $G_{g}$
is the antenna gain over the main lobe of the beam-pattern of $g\\!B\\!S$ $g$;
$N^{g,a}$ is the thermal noise; and $I_{g,a}^{\mathcal{B}}$ is the
interference coming from the remaining backhaul links of the network.

Backhaul links reuse the spectrum used for ground cellular connections,
although using beam-patterns pointing to the air, while antennas that provide
service to ground users are pointing mainly to the ground. Hence, we assume
that the interference suffered by a backhaul link $(g,a)$ is dominated by the
interference from other backhaul links. Hence, the interference suffered by a
backhaul link $(g,a)$ is:

$\displaystyle
I_{g,a}^{\mathcal{B}}=\sum\limits_{g^{\prime}\in\mathcal{G}\setminus\\{g\\}}P_{Tx}^{g^{\prime}}\\!\cdot\\!G_{g^{\prime}}(\phi_{g^{\prime},a})\\!\cdot\\!10^{-L_{\mathcal{B}}(g^{\prime},a)/10},$
(10)

where $\phi_{g^{\prime},a}$ is the angle between the main lobe direction of
the antenna of $g^{\prime}$ and the position of $a\\!B\\!S$ $a$. In case a
$g\\!B\\!S$ $g^{\prime}$ does no set any backhaul wireless link, this
$g\\!B\\!S$ will not affect interference, and $P_{Tx}^{g^{\prime}}$ will be
considered as zero.

##  B Overall Complexity of PADD

Our proposed algorithm consists of sequential steps, some of which involve
operations that can be parallelized and can run on a centralized or
distributed network orchestrator. Specifically, at each iteration, i.e., for
fixed positions of drones and having identified the least fit drone, we have:

  1. Step 1

User association. This can be implemented on parallel threads: one thread per
UE ranks the candidate list of BSs to attach to, then a separate thread
computes the association in at most as many rounds as the number of BSs, as
seen in Section 5.2. With $|\mathcal{U}|$ parallel threads, the complexity
required for this step to rank the list of BSs is
$\mathcal{O}\\!\left(\frac{1}{2\log\\!2}|\mathcal{B}|\log|\mathcal{B}|\right)$
sums and comparisons.

  2. Step 2

Backhaul association. The next step consists in solving the GAP-knap problem
for backhaul association, which must be done with a single thread, as shown in
Section 5.3. The complexity of combining dynamic programming and the knapsack
problem to solve this step is analyzed in [36] and it is
$\mathcal{O}\left(2A_{g}|\mathcal{G}||\mathcal{A}|\right)$ sums or
comparisons.

  3. Step 3

Resource allocation. This requires a thread per each $g\\!B\\!S$ for which we
need to solve the convex problem of Fig. 6, as shown in Section 5.4. The time
needed to complete this step is therefore the time needed to solve a single
problem. The complexity of the costliest case is
$\mathcal{O}\left(4A_{g}U_{\max}\right)$ simple operations (e.g., sums), as
discussed in Section D.

  4. Step 4

Utility evaluation. The current configuration is evaluated in terms of
$\alpha$-fair utility, which has the cost of $|\mathcal{U}||\mathcal{B}|$ sums
or comparisons (for $\alpha\\!\rightarrow\\!\infty$) plus $|\mathcal{B}|$
power operations (for $\alpha\neq 0,1$) or logarithms (for $\alpha=1$). As
discussed at the beginning of this section and depicted in Fig. 4, the current
configuration can be discarded and a new position is probed for the current
least fit drone (back to Step 1).111When the optimization problem is
initialized, there is no least fit drone yet; thus, the first iteration
executes Step 1-Step 4 only once, then it moves to Step 5 to make the first
least fit selection.

  5. Step 5

Least fit selection. Eventually, the algorithm computes the least fit drone
with a single thread, with complexity $\mathcal{O}(4|\mathcal{A}|U_{\max})$
operations, as two kinds of additional different utilities need to be computed
in Eq. (1) of Section 5.5.

With discretized drone positions, and indicating $N$ the number of points that
can be probed at each iteration for each drone, the complexity of one
iteration is $N$ times the complexity of Step 1-Step 4, plus the complexity of
Step 5. The overall complexity is therefore
$\mathcal{O}\left(N\left(\frac{1}{2\log\\!2}|\mathcal{B}|\log|\mathcal{B}|+2A_{g}|\mathcal{G}||\mathcal{A}|+4A_{g}U_{\max}+|\mathcal{U}||\mathcal{B}|\right)\right)$
sums or comparisons and
$\mathcal{O}\left(N|\mathcal{B}|+4|\mathcal{A}|U_{\max}\right)$ powers (or
logarithms) in the worst case. This complexity is low-degree polynomial, with
small constant factors, and linear in most of the parameters. The number of
iterations required by EO is not bounded unless a hard limit is imposed,
however, the approach is designed to quickly approach a local optimum.
In our experiments we have observed no more than a hundred of iterations,
without imposing any limit to the number of iterations. The number of possible
testing positions for the least fit drone is $N=400$ but we have observed
numbers below $5$. Indeed, on average, the initial position can be improved
with probability 0.5, which means that, if we assume that positions are
selected at random with no memory, the iteration stops after testing, on
average, 2 positions in the ball around the least fit drone (it would be a
Bernoulli process). If we consider memory, the Bernoulli approximation upper
bounds the average number of probing attempts, because the probability that
the next position will be better than the current one, and the iteration will
hence stop, will grow attempt after attempt.
Besides, our algorithm implemented on Matlab and running on a linux machine
with an Intel Xeon E5-2670 processor at $2.6$ GHz took 10 s in the worst case
(i.e., it required about 800 floating-point operations, at 32 flops/cycle
involving 64-bit operations), consuming only $1$ GB of RAM, to optimize drone
positions.

##  C Robustness of PADD

Fig. 18 depicts the CDF of the relative loss due to erroneous user position
estimation, i.e., the relative loss of utility due to optimizing drone
positions according to erroneous user positions. The figure also shows the
error in terms of throughput and fairness separately. The loss is below $10\%$
in all cases, and below $3\%$ for utility and throughput in more than half of
the cases, while the average loss is below $5\%$.

Figure 18: Robustness validation for $\alpha\in\\{0,1,\infty\\}$. CDF of the
relative loss. $G=10$, $A=5$, $U=1000$. Scenario: PPP.

##  D Optimal Solution for Bandwidth Allocation

Bandwidth allocation is done by solving the convex optimization program
presented in Fig. 6.We apply KKT conditions [44] so to find solutions without
having to resort to a solver. To simplify the notation, we consider that,
given a set of generic entities $\mathcal{E}$, there is a bijective function
$\sigma\\!:\\!\\{1,\ldots,|\mathcal{E}|\\}\\!\longrightarrow\\!\mathcal{E}$
that maps entities onto integer numbers. With that, we can denote, e.g., the
throughput of a user $u$ indistinctly as $T_{u}$ (with $u$ being a user) or
$T_{i}$ (with $i$ being an integer number), where $\sigma(i)=u$.

We next show how to solve the bandwidth allocation problem by considering two
separate cases: with and without drones.

###  D.1 Without Drones

In this scenario, there are no drone. CP of Fig. 6 simplifies considerably
since $g\\!B\\!S$ $g$ only needs to manage resources to be split among
$g\\!B\\!S$-served users, i.e., the set $\mathcal{U}_{g}$. Hence, we solve the
following CP:

${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\begin{cases}\qquad\qquad\qquad\qquad\qquad\mathclap{\max\limits_{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w_{u},}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}T_{u}}}U_{\mathrm{cvx},g}^{\alpha}=\begin{cases}\sum\limits_{u\in\mathcal{U}_{g}}\left(T_{u}\right)^{1-\alpha}\cdot\frac{1}{1-\alpha},&\alpha\neq
1;\\\
\sum\limits_{u\in\mathcal{U}_{g}}\log\left(T_{u}\right),&\alpha=1;\end{cases}}\\\
\text{s.t.:}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w_{u}\geq
W^{\min}_{\mathcal{G}},}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\\!\\!\forall
u\\!\in\\!\mathcal{U}_{g};}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum\limits_{u\in\mathcal{U}_{g}}\\!w_{u}=W_{\mathcal{G}};}&\\!\\!\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}T_{u}\\!\leq\\!w_{u}\log_{2}\\!\left(1\\!+\\!\gamma_{b,u}\right)\\!,}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\\!\\!\forall
u\\!\in\\!\mathcal{U}_{g};}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum\limits_{u\in\mathcal{U}_{g}}\\!T_{u}\\!\leq\tau_{g}.}&\\!\\!\end{cases}}$
(11)

Denote as
$x=[\\{w_{u}\\}_{u\in\mathcal{U}_{g}},\\{T_{u}\\}_{u\in\mathcal{U}_{g}}]\in\mathds{R}^{2|\mathcal{U}_{g}|}$
the vector of variables. Then, depending on the value of $\alpha$, we define
the KKT functions as follows:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!f^{\alpha}(x)=$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-\sum\limits_{u\in\mathcal{U}_{g}}T_{u}^{1-\alpha}$
$\displaystyle\text{if }\alpha\neq 1;$ (12)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!f^{\log}(x)=$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!-\sum\limits_{u\in\mathcal{U}_{g}}\log
T_{u}$ $\displaystyle\text{if }\alpha=1;$ (13)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!g_{u}(x)=$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!W_{\mathcal{G}}^{\min}-w_{u},$
$\displaystyle\forall u\in\mathcal{U}_{g};$ (14)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!g_{|\mathcal{U}_{g}|+u}(x)=$
$\displaystyle T_{u}-w_{u}\log_{2}\\!\left(1\\!+\\!\gamma_{g,u}\right),$
$\displaystyle\forall u\in\mathcal{U}_{g};$ (15)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!g_{2|\mathcal{U}_{g}|+1}(x)=$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\sum\limits_{u\in\mathcal{U}_{g}}T_{u}-\tau_{g};$
(16) $\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!h(x)=$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\sum\limits_{u\in\mathcal{U}_{g}}w_{u}-W_{\mathcal{G}}.$
(17)

The corresponding KKT gradients are as follows:

$\displaystyle\\!\\!\nabla\\!f^{\alpha}(\\!x\\!)=$
$\displaystyle\left[0,\ldots,\left\\{(\alpha\\!-\\!1)T_{u}^{-\alpha}\right\\}_{u\in\mathcal{U}_{g}}\right]\\!;$
(18) $\displaystyle\\!\\!\nabla\\!f^{\log}(\\!x\\!)=$
$\displaystyle\left[0,\ldots,0,\left\\{\frac{-1}{T_{u}}\right\\}_{u\in\mathcal{U}_{g}}\right]\\!;$
(19) $\displaystyle\\!\\!\nabla\\!g_{u}(\\!x\\!)=$
$\displaystyle\left[0,\ldots,0,-1|_{u},0,\ldots,0\right]\\!,$
$\displaystyle\forall u\\!\in\\!\mathcal{U}_{g};$ (20)
$\displaystyle\\!\\!\nabla\\!g_{|\mathcal{U}_{g}|\\!+\\!u}(\\!x\\!)=$
$\displaystyle\left[0,...,0,-\\!\log_{2}\\!(\\!1\\!\\!+\\!\\!\gamma_{g,u})|_{u},0,...,0,1|_{|\mathcal{U}_{g}|\\!+\\!u},0,...,0\right]\\!,$
$\displaystyle\forall u\\!\in\\!\mathcal{U}_{g};$ (21)
$\displaystyle\\!\\!\nabla\\!g_{2|\mathcal{U}_{g}|\\!+\\!1}(\\!x\\!)=$
$\displaystyle\left[0,\ldots,0,1|_{|\mathcal{U}_{g}|+1},1,\ldots,1\right]\\!;$
(22) $\displaystyle\\!\\!\nabla\\!h(\\!x\\!)=$
$\displaystyle\left[1,\ldots,1|_{|\mathcal{U}_{g}|},0,\ldots,0\right]\\!.$
(23)

Now, the KKT conditions state that if we find a vector
$x^{*}\in\mathds{R}^{2|\mathcal{U}_{g}|}$ that is feasible (i.e., it satisfy
the original constraints of the problem), and multipliers $\mu_{i}\geq 0$
$\forall 1\leq i\leq 2|\mathcal{U}_{g}|+1$ and $\nu$, and not all of them
null, then the following equations hold and $x^{*}$ is the optimal solution:

$\displaystyle\overrightarrow{0}$ $\displaystyle=$ $\displaystyle\nabla
f(x^{*})+\sum\limits_{i=1}^{2|\mathcal{U}_{g}|+1}\mu_{i}\nabla
g_{i}(x^{*})+\nu\nabla h(x^{*});$ (24) $\displaystyle\mu_{i}g_{i}(x^{*})$
$\displaystyle=$ $\displaystyle 0,\qquad\forall 1\leq i\leq
2|\mathcal{U}_{g}|+1;$ (25) $\displaystyle h(x^{*})$ $\displaystyle=$
$\displaystyle 0,$ (26)

where $f\in\\{\\{f^{\alpha}\\}_{\alpha\in[0,1[},f^{\log}\\}$.

Hence, in what follows we describe how we have found the optimal solution for
each particular case. We analyze separately the MaxThr case (when $\alpha=0$),
the $\alpha$Fair case (when $\alpha\in]0,1[$), the PropFair case (when
$\alpha=1$) and the MaxMin case (when $\alpha\rightarrow\infty$).

####  D.1.1 $\alpha$Fair optimum ($\alpha\in]0,1[$)

In this case,
$f(x)=f^{\alpha}(x)=-\\!\\!\sum\limits_{u\in\mathcal{U}_{g}}\\!T_{u}^{1-\alpha}$,
and the KKT conditions are as follows:

$\displaystyle-\mu_{i}-\log_{2}\\!\left(1\\!+\\!\gamma_{g,i}\right)\\!\cdot\\!\mu_{|\mathcal{U}_{g}|+i}\\!+\\!\nu=0,$
$\displaystyle\forall 1\leq i\leq|\mathcal{U}_{g}|;$ (27)
$\displaystyle\mu_{|\mathcal{U}_{g}|+i}+\mu_{2|\mathcal{U}_{g}|+1}-(1-\alpha)T_{i}^{-\alpha}=0,$
$\displaystyle\forall 1\leq i\leq|\mathcal{U}_{g}|;$ (28)
$\displaystyle\mu_{i}\cdot\left(W_{\mathcal{G}}^{\min}-w_{i}\right)=0,$
$\displaystyle\forall 1\leq i\leq|\mathcal{U}_{g}|;$ (29)
$\displaystyle\mu_{|\mathcal{U}_{g}|+i}\\!\cdot\\!\left(T_{i}-w_{i}\log_{2}\\!\left(1\\!+\\!\gamma_{g,i}\right)\right)=0,$
$\displaystyle\forall 1\leq i\leq|\mathcal{U}_{g}|;$ (30)
$\displaystyle\mu_{2|\mathcal{U}_{g}|+1}\cdot\Big{(}\sum\limits_{u\in\mathcal{U}_{g}}T_{u}-\tau_{g}\Big{)}=0;$
(31) $\displaystyle\sum\limits_{u\in\mathcal{U}_{g}}w_{u}-W_{\mathcal{G}}=0.$
(32)

With unbounded backbone capacity $\boldsymbol{\tau_{g}}$. Assuming that
$\tau_{g}\rightarrow+\infty$, then
$\sum\limits_{u\in\mathcal{U}_{g}}T_{u}<\tau_{g}$ for all feasible $x$. Hence,
according to Eq. (31), $\mu_{2|\mathcal{U}_{g}|+1}=0$. Therefore, for $1\leq
i\leq|\mathcal{U}_{g}|$, Eq. (28) becomes
$\mu_{|\mathcal{U}_{g}|+i}=(1-\alpha)T_{i}^{-\alpha}$, and Eq. (27) becomes
$\nu=\mu_{i}+(1-\alpha)T_{i}^{-\alpha}\cdot\log_{2}\left(1+\gamma_{g,i}\right)$.
As $\nu$ is a constant, we search for expressions of $\mu_{i}$ and $T_{i}$
that make $\nu$ independent of index $i$. For convenience, we define set
$\mathcal{J}$ as those indices $i$ for which $w_{i}$ is the minimum guaranteed
value, i.e.:

$\displaystyle\mathcal{J}=\left\\{i\,\mid\,w_{i}=W_{\mathcal{G}}^{\min}\right\\}.$
(33)

We start by assuming that $\mu_{i}=0$, $\forall i\notin\mathcal{J}$, so that
(29) is satisfied $\forall i\notin\mathcal{J}$. Given the SINR values, this
assumption makes $T_{i}$ become a function of $\nu$ only, $\forall
i\notin\mathcal{J}$, i.e.,
$T_{i}=\left(\frac{1-\alpha}{\nu}\log_{2}\left(1+\gamma_{g,i}\right)\right)^{\frac{1}{\alpha}}$,
$\forall i\notin\mathcal{J}$. In addition, as $\tau_{g}$ is unbounded, we can
assume that $T_{i}$ must take the highest possible value, which is the Shannon
capacity $w_{i}\log_{2}\left(1+\gamma_{g,i}\right)$, $\forall 1\leq
i\leq|\mathcal{U}_{g}|$. Hence,
$w_{i}=\left(\frac{1-\alpha}{\nu}\right)^{\frac{1}{\alpha}}\left(\log_{2}\left(1+\gamma_{g,i}\right)\right)^{\frac{1-\alpha}{\alpha}}$
for all cases in which $w_{i}>W_{\mathcal{G}}^{\min}$ (i.e., $\forall
i\notin\mathcal{J}$) and $w_{i}=W_{\mathcal{G}}^{\min},\forall
i\in\mathcal{J}$. Replacing the above expression for $w_{i}$ in Eq. (32), we
obtain an equation in which $\nu$ is the only unknown to be derived. We
therefore have obtained the following solution:

Algorithm 1 Derivation of set $\mathcal{J}$ in $\alpha$Fair (unbounded
$\tau_{g}$).

0: $\log_{2}\left(1+\gamma_{g,i}\right)$, $\forall 1\leq
i\leq|\mathcal{U}_{g}|$.

1: Initialize: $\mathcal{J}\leftarrow\emptyset$;
$w_{i}\leftarrow\frac{W_{\mathcal{G}}-|\mathcal{J}|\cdot
W_{\mathcal{G}}^{\min}}{\sum\limits_{j\notin\mathcal{J}}\left(\frac{\log_{2}\left(1+\gamma_{g,j}\right)}{\log_{2}\left(1+\gamma_{g,i}\right)}\right)^{\frac{1-\alpha}{\alpha}}}$,
$\forall i\notin\mathcal{J}$.

2: while $\exists
i_{1}\notin\mathcal{J}\,\mid\,w_{i_{1}}<W_{\mathcal{G}}^{\min}$ do

3:
$i_{0}\leftarrow\arg\max_{i\notin\mathcal{J}}\log_{2}\left(1+\gamma_{g,i}\right)$.

4: $\mathcal{J}\leftarrow\mathcal{J}\cup\\{i_{0}\\}$.

5: $w_{i}\leftarrow\frac{W_{\mathcal{G}}-|\mathcal{J}|\cdot
W_{\mathcal{G}}^{\min}}{\sum\limits_{j\notin\mathcal{J}}\left(\frac{\log_{2}\left(1+\gamma_{g,j}\right)}{\log_{2}\left(1+\gamma_{g,i}\right)}\right)^{\frac{1-\alpha}{\alpha}}}$,
$\forall i\notin\mathcal{J}$.

6: end while

$\displaystyle w_{i}=W_{\mathcal{G}}^{\min},$ $\displaystyle\forall
i\in\mathcal{J};$ (34) $\displaystyle
w_{i}=\frac{W_{\mathcal{G}}-|\mathcal{J}|\cdot
W_{\mathcal{G}}^{\min}}{\sum\limits_{j\notin\mathcal{J}}\\!\\!\left(\frac{\log_{2}\\!\left(1\\!+\\!\gamma_{g,j}\right)}{\log_{2}\\!\left(1\\!+\\!\gamma_{g,i}\right)}\right)^{\frac{1\\!-\\!\alpha}{\alpha}}}\\!,$
$\displaystyle\forall i\notin\mathcal{J};$ (35) $\displaystyle
T_{i}=w_{i}\log_{2}(1\\!+\\!\gamma_{g,i})$ $\displaystyle\forall 1\leq
i\leq|\mathcal{U}_{g}|;$ (36)
$\displaystyle\mu_{i}=\nu-\log_{2}\left(1+\gamma_{g,i}\right)(1-\alpha)T_{i}^{-\alpha},$
$\displaystyle\forall i\in\mathcal{J};$ (37) $\displaystyle\mu_{i}=0,$
$\displaystyle\forall i\notin\mathcal{J};$ (38)
$\displaystyle\mu_{|\mathcal{U}_{g}|+i}=(1\\!-\\!\alpha)\,T_{i}^{-\alpha},$
$\displaystyle\forall 1\leq i\leq|\mathcal{U}_{g}|;$ (39)
$\displaystyle\mu_{2|\mathcal{U}_{g}|+1}=0;$ (40) $\displaystyle\nu=\hskip
71.13188pt\mathclap{(1\\!-\\!\alpha)\\!\\!\left(\frac{\sum\limits_{j\notin\mathcal{J}}\left(\log_{2}\left(1+\gamma_{g,j}\right)\right)^{\frac{1-\alpha}{\alpha}}}{W_{\mathcal{G}}-|\mathcal{J}|\cdot
W_{\mathcal{G}}^{\min}}\\!\right)^{\\!\\!\\!\alpha}\\!\\!\\!.}$ (41)

The solution we have built accomplishes all KKT equations and works also for
the case $\sum\limits_{u\in\mathcal{U}_{g}}T_{u}=\tau_{g}$. In particular,
note in Eq. (37) that $\mu_{i}\geq 0$, $\forall i\in\mathcal{J}$, which can be
seen by replacing the expressions for $T_{i}$, $w_{i}$ and $\nu$ in the
equation and consider that, $\forall i\in\mathcal{J}$,
$W_{\mathcal{G}}^{\min}\geq\frac{W_{\mathcal{G}}-|\mathcal{J}|\cdot
W_{\mathcal{G}}^{\min}}{\sum\limits_{j\notin\mathcal{J}}\\!\\!\left(\frac{\log_{2}\\!\left(1\\!+\\!\gamma_{g,j}\right)}{\log_{2}\\!\left(1\\!+\\!\gamma_{g,i}\right)}\right)^{\frac{1\\!-\\!\alpha}{\alpha}}}$.

A valid set $\mathcal{J}$ always exists and it is built as shown in Algorithm
1. The algorithm starts with an empty set and, while it is not true that
$w_{i}\geq W_{\mathcal{G}}^{\min}$, $\forall i\notin\mathcal{J}$, keeps adding
to $\mathcal{J}$ one user at a time, the one with the highest SINR, which is
the user that consumes less resources. The algorithm stops for sure,
eventually after moving all users, with a set $\mathcal{J}$ for which KKT and
feasibility conditions are satisfied, therefore the built solution is optimal.

With limited backbone capacity $\boldsymbol{\tau_{g}}$. In order to find the
optimal solution in the general case, we start from the solution for
$\tau_{g}\rightarrow+\infty$. In case this procedure yields
$\sum\limits_{u\in\mathcal{U}_{g}}T_{u}\leq\tau_{g}$, we have the optimal
solution described above. Otherwise, we have that
$\sum\limits_{u\in\mathcal{U}_{g}}T_{u}>\tau_{g}$, which is unfeasible.
However, we can build a feasible solution by decreasing some values of
$T_{i}$, motivated by the consideration that the utility function is strictly
concave, so that the least utility reduction is obtained by reducing the
highest throughput value, as shown in Algorithm 2 (for $n=|\mathcal{U}_{g}|$,
$\\{x_{i}\\}_{i=1}^{n}=\\{T_{u}\\}_{u\in\mathcal{U}_{g}}$ and $X=\tau_{g}$).

Algorithm 2 is iterative, and it is based on the following principle: Because
of concavity, if only one user has maximal throughput
$T^{\prime}\\!=\\!\max_{i}\\{T_{i}\\}$, and second best throughput is
$T^{\prime\prime}$, any throughput reduction from $T^{\prime}$ to
$T^{\prime}-y\geq T^{\prime\prime}$ maximizes utility for a reduction of $y$
in $\sum_{i}T_{i}$; if we have two users with maximal throughput, and we have
to reduce the sum of throughputs by $y$, then strict concavity assures that
$2\left(T^{\prime}-y/2\right)^{1-\alpha}>\left(T^{\prime}-(y/2-\epsilon)\right)^{1-\alpha}+\left(T^{\prime}-(y/2+\epsilon)\right)^{1-\alpha},\forall\epsilon\in]0;y/2],y/2\in]0,T^{\prime}-T^{\prime\prime}]$,
so the best aggregate utility is obtained by decreasing the two highest
throughput values both to $T^{\prime}-y/2$. As it is easy to see, in case of
tie between $n$ users, the best choice consists in reducing their throughputs
to $T^{\prime}-y/n,\forall y/n\in]0,T^{\prime}-T^{\prime\prime}]$.

####  D.1.2 MaxThr optimum ($\alpha=0$)

In this case, $f(x)=f^{0}(x)=-\\!\\!\sum\limits_{u\in\mathcal{U}_{g}}T_{u}$,
and the KKT conditions are like for the case $\alpha\in]0,1[$, except
$\alpha=0$.

Algorithm 2 Optimum for strictly concave increasing utility

0: $\\{x_{i}\\}_{i=1}^{n}$, $X$.

1: while $X-\sum\limits_{i=1}^{n}x_{i}<0$ do

2: $X_{M}\leftarrow\\{x_{i}\,\mid\,x_{i}=\max_{1\leq j\leq n}x_{j}\\}$.

3: $x_{M_{2}}\leftarrow\arg\max\\{x_{i}\,\mid\,x_{i}\notin X_{M}\\}$ (assume
$|x_{M_{2}}|=1$).

4: $R\leftarrow-X+\sum\limits_{i=1}^{n}x_{i}>0$.

5: $R_{j}\leftarrow\frac{R}{|X_{M}|}$, $\forall 1\leq j\leq n\text{ with
}x_{j}\in X_{M}$.

6: $x_{j}\leftarrow
x_{j}\\!-\\!\min\\!\left\\{R_{j},x_{j}\\!-\\!x_{M_{2}}\right\\}$, $\\!\forall
1\\!\leq\\!j\\!\leq\\!n\text{ with }x_{j}\\!\in\\!X_{M}$.

7: end while

With unbounded backbone capacity $\boldsymbol{\tau_{g}}$. First, we look for a
solution assuming that $\tau_{g}$ is unbounded, i.e., assuming that
$\sum\limits_{u\in\mathcal{U}_{g}}T_{u}<\tau_{g}$ for all feasible $x$. Hence,
according to Eq. (31), $\mu_{2|\mathcal{U}_{g}|+1}=0$. Also, according to Eq.
(28), $\mu_{|\mathcal{U}_{g}|+i}=1$, $\forall 1\leq i\leq|\mathcal{U}_{g}|$,
and according to Eq. (30), $T_{i}=w_{i}\log_{2}\left(1+\gamma_{g,i}\right)$,
$\forall 1\leq i\leq|\mathcal{U}_{g}|$. Moreover, according to Eq. (27),
$\nu-\mu_{i}=\log_{2}\left(1+\gamma_{g,i}\right)$, $\forall 1\leq
i\leq|\mathcal{U}_{g}|$.

Now, let $i_{0}=\arg\max_{1\leq
i\leq|\mathcal{U}_{g}|}\log_{2}\left(1+\gamma_{g,i}\right)$ and let
$w_{i}=W_{\mathcal{G}}^{\min}$, $\forall 1\leq i\neq
i_{0}\leq|\mathcal{U}_{g}|$. Hence, according to Eq. (32),

$w_{i_{0}}=W_{\mathcal{G}}-\sum\limits_{\begin{subarray}{c}i=1\\\ i\neq
i_{0}\end{subarray}}^{|\mathcal{U}_{g}|}w_{i}=W_{\mathcal{G}}-\sum\limits_{\begin{subarray}{c}i=1\\\
i\neq
i_{0}\end{subarray}}^{|\mathcal{U}_{g}|}W_{\mathcal{G}}^{\min}=W_{\mathcal{G}}-(|\mathcal{U}_{g}|-1)W_{\mathcal{G}}^{\min}.$

Now, let $\mu_{i_{0}}=0$. Hence, according to Eq. (27),
$\nu=\log_{2}\left(1+\gamma_{g,i_{0}}\right)$ and
$\mu_{i}=\log_{2}\left(1+\gamma_{g,i_{0}}\right)-\log_{2}\left(1+\gamma_{g,i}\right)$,
$\forall 1\leq i\leq|\mathcal{U}_{g}|$ (note that $\mu_{i}\geq 0$ $\forall
1\leq i\leq|\mathcal{U}_{g}|$ by definition of $i_{0}$).

Since the found solution $x$ is feasible and accomplishes all KKT conditions,
this is necessary the optimal solution.

With limited backbone capacity $\boldsymbol{\tau_{g}}$. We can observe that
the optimal solution when $\sum\limits_{u\in\mathcal{U}_{g}}T_{u}<\tau_{g}$
corresponds to assigning the maximum possible amount of resources to the users
with highest SINR. If we now consider a bounded value for $\tau_{g}$, we can
build the solution following the same scheme to build Algorithm 3, which finds
the optimal solution, as explained in what follows.

Let $i_{0}=\arg\max_{1\leq
i\leq|\mathcal{U}_{g}|}\log_{2}\left(1+\gamma_{g,i}\right)$. Let
$w_{i}=W_{\mathcal{G}}^{\min}$, $\forall 1\leq i\neq
i_{0}\leq|\mathcal{U}_{g}|$ and
$w_{i_{0}}=W_{\mathcal{G}}-\sum\limits_{\begin{subarray}{c}i=1\\\ i\neq
i_{0}\end{subarray}}W_{\mathcal{G}}^{\min}=W_{\mathcal{G}}-(|\mathcal{U}_{g}|-1)W_{\mathcal{G}}^{\min}$.
In the algorithm, we initially assign
$T_{i_{0}}\\!=\\!\min\left(\tau_{g},w_{i_{0}}\log_{2}\\!\left(1\\!+\\!\gamma_{g,i_{0}}\right)\right)$,
$T_{i}=0$, $\forall 1\\!\leq\\!i\\!\neq\\!i_{0}\\!\leq\\!|\mathcal{U}_{g}|$,
and define the set
$\overline{\mathcal{U}_{g}}=\mathcal{U}_{g}\setminus\\{i_{0}\\}$ (step 1).
Next, while $\sum\limits_{i=1}^{|\mathcal{U}_{g}|}T_{i}<\tau_{g}$ and
$\overline{\mathcal{U}_{g}}\neq\emptyset$, we do the following three
assignments: (i)
$j_{0}=\arg\max_{j\in\overline{\mathcal{U}_{g}}}\log_{2}\left(1+\gamma_{g,j}\right)$
(step 3), (ii)
$T_{j_{0}}=\min\left(\tau_{g}-\sum\limits_{i\in\mathcal{U}_{g}}T_{i},w_{j_{0}}\log_{2}\left(1+\gamma_{g,j_{0}}\right)\right)$
(step 4), and (iii)
$\overline{\mathcal{U}_{g}}=\overline{\mathcal{U}_{g}}\setminus\\{j_{0}\\}$
(step 5). This procedure finds values of
$\\{T_{i}\\}_{i=1}^{|\mathcal{U}_{g}|}$ such that
$\sum\limits_{i=1}^{|\mathcal{U}_{g}|}T_{i}\leq\tau_{g}$. Besides, if
$\overline{\mathcal{U}_{g}}\neq\emptyset$ when the algorithm concludes, we
have that $\sum\limits_{i=1}^{|\mathcal{U}_{g}|}T_{i}=\tau_{g}$.

Note that in case $\sum\limits_{i=1}^{|\mathcal{U}_{g}|}T_{i}<\tau_{g}$, the
solution of Algorithm 3 corresponds with the solution for unbounded
$\tau_{g}$. Moreover, in case the output of the algorithm is such that
$\sum\limits_{i=1}^{|\mathcal{U}_{g}|}T_{i}=\tau_{g}$, the solution is also
optimal because, with $\alpha=0$, $\sum\limits_{i=1}^{|\mathcal{U}_{g}|}T_{i}$
is the utility function, and $\tau_{g}$ is the maximum that that sum can
achieve, according to the last constraint in (11). Therefore, the solution
computed with Algorithm 3 is always optimal for $\alpha=0$.

Algorithm 3 Optimum MaxThr in Backhaul-free scenario

1: Initialize:
$T_{i_{0}}=\min\left(\tau_{g},\;w_{i_{0}}\log_{2}\left(1+\gamma_{g,i_{0}}\right)\right)$,
$T_{i}=0$, $\forall 1\leq i\neq i_{0}\leq|\mathcal{U}_{g}|$,
$\overline{\mathcal{U}_{g}}=\mathcal{U}_{g}\setminus\\{i_{0}\\}$.

2: while $\sum\limits_{i=1}^{|\mathcal{U}_{g}|}T_{i}<\tau_{g}$ and
$\overline{\mathcal{U}_{g}}\neq\emptyset$ do

3:
$j_{0}=\arg\max_{j\in\overline{\mathcal{U}_{g}}}\log_{2}\left(1+\gamma_{g,j}\right)$.

4:
$T_{j_{0}}=\min\left(\tau_{g}-\sum\limits_{i\in\mathcal{U}_{g}}T_{i},w_{j_{0}}\log_{2}\left(1+\gamma_{g,j_{0}}\right)\right)$.

5:
$\overline{\mathcal{U}_{g}}=\overline{\mathcal{U}_{g}}\setminus\\{j_{0}\\}$.

6: end while

####  D.1.3 PropFair optimum ($\alpha=1$)

This case aims to maximize the network proportional fairness metric,
$f(x)=f^{\log}(x)=-\sum\limits_{u\in\mathcal{U}_{g}}\log\left(T_{u}\right)$.
The KKT conditions for this case are like for the case $\alpha\in]0,1[$,
except condition (28) is replaced by the following one:

$\displaystyle\mu_{|\mathcal{U}_{g}|+i}+\mu_{2|\mathcal{U}_{g}|+1}-\frac{1}{T_{i}}=0,$
$\displaystyle\forall 1\leq i\leq|\mathcal{U}_{g}|;$ (42)

With unbounded backbone capacity $\boldsymbol{\tau_{g}}$. If $\tau_{g}$ is
unbounded, Eq. (31) requires that $\mu_{2|\mathcal{U}_{g}|+1}=0$. Hence,
according to Eq. (42), $\mu_{|\mathcal{U}_{g}|+i}=\frac{1}{T_{i}}$, $\forall
1\leq i\leq|\mathcal{U}_{g}|$.

Now, $\forall 1\leq i\leq|\mathcal{U}_{g}|$, let $\mu_{i}=0$ as in the
previous cases, and $w_{i}=\frac{W_{\mathcal{G}}}{|\mathcal{U}_{g}|}$,
$T_{i}=\frac{W_{\mathcal{G}}}{|\mathcal{U}_{g}|}\log_{2}\left(1+\gamma_{g,i}\right)$.
Hence, according to Eq. (27), $\nu=\frac{|\mathcal{U}_{g}|}{W_{\mathcal{G}}}$.

Since the found feasible solution $x$ accomplishes all KKT conditions, this is
the optimal solution. Note that, as expected, in a PropFair without backbone
constraints, all users would obtain the same amount of resources.

With limited backbone capacity $\boldsymbol{\tau_{g}}$. To find the optimal
solution for the generic case, in which $\tau_{g}$ is limited, we start from
the solution for unbounded backbone capacity. If this results in
$\sum\limits_{u\in\mathcal{U}_{g}}T_{u}>\tau_{g}$, being $\log(x)$ is a
strictly concave increasing function, we can apply Algorithm 2 with
$n=|\mathcal{U}_{g}|$,
$\\{x_{i}\\}_{i=1}^{n}=\\{T_{u}\\}_{u\in\mathcal{U}_{g}}$ and $X=\tau_{g}$ in
order to find the optimal solution of the PropFair case.

####  D.1.4 MaxMin optimum ($\alpha\rightarrow\infty$)

The utility function $f(x)=\min_{u\in\mathcal{U}_{g}}T_{u}$ is not
differentiable, so we cannot directly apply KKT conditions. However, following
the spirit of MaxMin optimization, we can reformulate the problem by adding a
new decision variable $T$, changing the utility function and adding one extra
set of constraints. This results in the following convex program:

$\begin{cases}\max T;&\\\ \text{s.t.:}\\\ T_{u}\geq T,&\\!\\!\forall
u\\!\in\\!\mathcal{U}_{g};\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}w_{u}\geq
W^{\min}_{\mathcal{G}},}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\\!\\!\forall
u\\!\in\\!\mathcal{U}_{g};}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum\limits_{u\in\mathcal{U}_{g}}\\!w_{u}=W_{\mathcal{G}};}&\\!\\!\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}T_{u}\\!\leq\\!w_{u}\log_{2}\\!\left(1\\!+\\!\gamma_{b,u}\right)\\!,}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\\!\\!\forall
u\\!\in\\!\mathcal{U}_{g};}\\\
{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum\limits_{u\in\mathcal{U}_{g}}\\!T_{u}\\!\leq\tau_{g}.}&\\!\\!\end{cases}$
(43)

Since we have added a new decision variable $T$, we denote the solution as
$x=\left[\\{w_{u}\\}_{u\in\mathcal{U}_{g}}\\},\\{T_{u}\\}_{u\in\mathcal{U}_{g}},T\right]\in\mathds{R}^{2|\mathcal{U}_{g}|+1}$.
Hence, we add to the KKT functions of Eqs. (12)–(17) the following KKT
functions:

$\displaystyle f^{\min}(x)=$ $\displaystyle-T;$ (44) $\displaystyle
g_{2|\mathcal{U}_{g}|+1+u}(x)=$ $\displaystyle T-T_{u},$ $\displaystyle\forall
u\in\mathcal{U}_{g}.$ (45)

Also, we have to add to KKT gradients of Eqs. (18)–(23) the following KKT
gradients:

$\displaystyle\nabla\\!f^{\min}(x)=$ $\displaystyle[0,\ldots,0,-1];$ (46)
$\displaystyle\nabla\\!g_{2|\mathcal{U}_{g}\\!|\\!+\\!1\\!+\\!u}(x)=$
$\displaystyle[0,...,0,-\\!1|_{|\mathcal{U}_{g}\\!|\\!+\\!u},0,...,0,1],$
$\displaystyle\\!\\!\forall u\\!\in\\!\mathcal{U}_{g}.$ (47)

In this case, applying KKT conditions, i.e., finding $\nu\in\mathds{R}$ and
positive constants $\\{\mu_{i}\\}_{i=1}^{3|\mathcal{U}_{g}|+1}$ not all of
them null, and $x\in\mathds{R}^{2|\mathcal{U}_{g}|+1}$ such that
$\overrightarrow{0}=\nabla
f^{\min}(x)+\sum\limits_{i=1}^{3|\mathcal{U}_{g}|+1}\mu_{i}\nabla
g_{i}(x)+\nu\nabla h(x)$ results into the following set of KKT equations:

$\displaystyle-\mu_{i}-\log_{2}\\!\left(1\\!+\\!\gamma_{g,i}\right)\\!\cdot\\!\mu_{|\mathcal{U}_{g}|+i}\\!+\\!\nu=0,$
$\displaystyle\\!\forall 1\\!\leq\\!i\\!\leq\\!|\mathcal{U}_{g}|;$ (48)
$\displaystyle-\mu_{|\mathcal{U}_{g}|+i}+\mu_{2|\mathcal{U}_{g}|+1}-\mu_{2|\mathcal{U}_{g}|+1+i}=0,$
$\displaystyle\\!\forall 1\\!\leq\\!i\\!\leq\\!|\mathcal{U}_{g}|;$ (49)
$\displaystyle-1+\sum\limits_{i=1}^{|\mathcal{U}_{g}|}\mu_{2|\mathcal{U}_{g}|+1+i}=0;$
(50) $\displaystyle\mu_{i}\cdot\left(W_{\mathcal{G}}^{\min}-w_{i}\right)=0,$
$\displaystyle\\!\forall 1\\!\leq\\!i\\!\leq\\!|\mathcal{U}_{g}|;$ (51)
$\displaystyle\mu_{|\mathcal{U}_{g}|+i}\\!\cdot\\!\left(T_{i}-w_{i}\log_{2}\\!\left(1\\!+\\!\gamma_{g,i}\right)\right)=0,$
$\displaystyle\\!\forall 1\\!\leq\\!i\\!\leq\\!|\mathcal{U}_{g}|;$ (52)
$\displaystyle\mu_{2|\mathcal{U}_{g}|+1}\cdot\Big{(}\sum\limits_{u\in\mathcal{U}_{g}}T_{u}-\tau_{g}\Big{)}=0;$
(53) $\displaystyle\mu_{2|\mathcal{U}_{g}|+1+i}\cdot\left(T-T_{i}\right)=0,$
$\displaystyle\\!\forall 1\\!\leq\\!i\\!\leq\\!|\mathcal{U}_{g}|;$ (54)
$\displaystyle\sum\limits_{u\in\mathcal{U}_{g}}w_{u}-W_{\mathcal{G}}=0.$ (55)

With unbounded backbone capacity $\boldsymbol{\tau_{g}}$. First, we look for a
solution assuming that $\tau_{g}$ is unbounded. Hence, according to Eq. (53),
$\mu_{2|\mathcal{U}_{g}|+1}=0$. As $\\{T_{i}\\}_{i=1}^{|\mathcal{U}_{g}|}$ is
not limited by $\tau_{g}$, we take
$T_{i}=w_{i}\log_{2}\left(1+\gamma_{g,i}\right)$, $\forall 1\leq
i\leq|\mathcal{U}_{g}|$.

We initially assume that $w_{i}=W_{\mathcal{G}}^{\min}$, $\forall 1\leq
i\leq|\mathcal{U}_{g}|$ and build an algorithm that finds the optimal solution
of CP (43) (see Algorithm 4). Please note that now
$T_{i}=W_{\mathcal{G}}^{\min}\cdot\log_{2}\left(1+\gamma_{g,i}\right)$,
$\forall 1\leq i\leq|\mathcal{U}_{g}|$ (step 1). We assume, without loss of
generality, that $T_{i-1}\leq T_{i}$, $\forall 2\leq i\leq|\mathcal{U}_{g}|$
(step 2), and define the set of indices $\mathcal{J}$ as those indices $i$
such that $T_{i}$ is minimum (step 3), i.e.,

$\displaystyle\mathcal{J}=\left\\{i\in\\{1,\ldots,|\mathcal{U}_{g}|\\}\,\mid\,T_{i}=\min\limits_{1\leq
j\leq|\mathcal{U}_{g}|}T_{j}\right\\}.$ (56)

Algorithm 4 works in the following way: while there are still bandwidth
resources to be allocated, i.e.,
$\sum\limits_{i=1}^{|\mathcal{U}_{g}|}w_{i}<W_{\mathcal{G}}$, and
$|\mathcal{J}|\neq|\mathcal{U}_{g}|$, we select the index
$j_{0}=\arg\min\\{T_{i}\,\mid\,i\notin\mathcal{J}\\}$ such that $T_{j_{0}}$ is
the lowest user throughput rate not equal to the minimum of the throughput
rates (step 5). Now, we aim to increase $w_{i}$, $\forall i\in\mathcal{J}$ as