Coverage Optimization with a Dynamic Network of Drone Relays

The integration of aerial base stations carried by drones in cellular networks offers promising opportunities to enhance the connectivity enjoyed by ground users. In this paper, we propose an optimization framework for the 3–D placement and repositioning of a fleet of drones with a realistic inter-drone interference model and drone connectivity constraints. We show how to maximize network coverage by means of an extremal-optimization algorithm. The design of our algorithm is based on a mixed-integer non-convex program formulation for a coverage problem that is NP-Complete, as we prove in the paper. We not only optimize drone positions in a 3–D space in polynomial time, but also assign flight routes solving an assignment problem and using a strong geometrical tool, namely Bézier curves, which are extremely useful for non-uniform and realistic topologies. Specifically, we propose to fly drones following Bézier curves to seek the chance of approaching to clusters of ground users. This enhances coverage over time while users and drones move. We assess the performance of our proposal for synthetic scenarios as well as realistic maps extracted from the topology of a capital city. We demonstrate that our framework is near-optimal and using Bézier curves increases coverage up to 47 percent while drones move.


INTRODUCTION
C ELLULAR networks are experiencing deep changes due to the advent of 5G technologies [1]. The need of flexible and adaptive management solutions, to address a highly mutable density of users, has allowed novel communication paradigms to emerge, e.g., device-to-device (D2D) communications, smart relay and the use of reconfigurable backhaul links controlled by Software-Defined Networking (SDN) tools [2]. Since the interest for mobile and topology-adaptive relays is now reviving-due to the techniques that make it doable in operational networks rather than just in theoretic speculations-in this paper we focus on the analysis of a key use-case: a fleet of coordinated drones carrying aerial relays meant to optimize cellular coverage.
The availability of broadband backhaul links allows manned and unmanned vehicles (e.g., drones) to carry mobile relays. The use of mobile relays brings unique opportunities to deploy adaptive and flexible networks that provide connectivity where fixed infrastructures lack operational connectivity [3]. Basic advantages from relays, like the reuse of cellular bands, have been studied by several authors, e.g., by Guo and O'Farrel [4], whose results illustrate how wireless networks are mainly impaired by interference. Other studies confirm that interference is commonly the main limiting factor also for the performance of relays operated over D2D sidelinks [5]. Hence, relays need to use dedicated frequencies in practice, which has the further advantage of simplifying resource allocation schemes [6].
Thus, we study the case of drones as relay stations transmitting on orthogonal frequencies with respect to ground base stations. However, due to limited spectrum resources, drones do interfere with each other and interference cannot be neglected both in users access and backhaul links, as we study in our work. These aspects have not been studied so far because existing works on drone relays focused on the optimal placement of a single drone in isolation, mainly under static scenarios, and neglect connectivity issues in the backhaul links between drones and ground base stations.
Although relay optimization can be used to support quality of service (QoS) in many different aspects, here we focus on providing support to user applications that require high probability to find the network available. This is important to guarantee essential communication means for voice and real time data communications. It is even more important for addressing the case of safety/emergency events and whenever notifying the status of users and having basic real-time communication tools is vital. These scenarios require to maximize the number of users that can be served, rather than maximizing on the allocated data rate or on the aggregate network rate. Therefore, we refer to QoS in terms of network availability with basic guaranteed rate.
Unlike other proposals, in this paper we focus on the optimization of 3-D hovering positions and flight routes for a fleet of drone relays aiding a ground cellular network, as depicted in Fig. 1. Drones are coordinated yet they mutually interfere. We optimize coverage based on the QoS offered by drones under realistic path-loss models for line-of-sight (LoS) and non-LoS (NLoS) communications and interference. Considering interference is key because it results in radically different coverage footprints and imposes strict constraints on the position of drones with respect to the position of ground base stations. We use extremal-optimization [7] and propose the On-demand Drone Coverage (OnDrone) algorithm, an extremal-optimization algorithm that computes near-optimally joint positions for drones, based on realistic assumptions on previous drone positions and interference, which is otherwise an intractable NP-Complete problem. We also propose for the first time the use of B ezier curves [8] for flight routes aiming to enhance communications over time. We assess the benefits of our optimization framework by (i) comparing OnDrone against the optimal solutions and state-of-the-art approaches for tractable cases; (ii) performing numerical simulations for larger networks with realistic topologies and environmental constraints; (iii) evaluating fleet repositioning using either B ezier curves or straight paths as drone routes. Our numerical results show that OnDrone is nearly optimal and outperforms state-of-the-art coverage solutions as proposed in [9] and [10]. Also, we show that the use of B ezier curves is key to boost coverage when studying drone repositioning in dense urban scenarios, and shows remarkable advantages over straight paths, as adopted in [11].
The main contributions of this paper are summarized: We propose a dynamic drone relay-aided network in which we maximize the coverage of ground users by means of aerial base stations with an interferenceaware on-demand multi-drone coverage framework that accounts for both user access and backhaul links. We prove that the problem is NP-Complete. We propose OnDrone, a light-weight algorithm based on extremal-optimization that solves the problem on-demand. We propose the use of a strong geometrical tool to design the flight paths of drones: the B ezier Scheme. We assess our proposals in realistic scenarios and topologies in comparison with state-of-the art solutions and show the gain of our proposals. The rest of the paper is structured as follows: Section 2 provides a discussion on related work regarding aerial networks and extremal-optimization. Section 3 presents the system model assumptions for the reference scenario and wireless channels, while Section 4 states and formulates the coverage problem, and shows that it is intractable. Section 5 details the overall optimization framework. Section 6 reports numerical results. Section 7 concludes the paper.

RELATED WORK
The usage of relays operating in the air space through mobile and non-terrestrial devices has been studied for several purposes, over different technologies. For instance, satellite networks [12] have already been deployed for several years. However, satellites aim to provide service to huge areas, typically at relatively low transmission rates. Moreover, satellites located hundreds of kilometers high are not able to adjust to ground users' topology, neither track the movement of small masses of users. In contrast, drone relay stations are able to dynamically change their position in real time at low altitudes, as the system evolves.
Google and Facebook carry the Loon [13] and Aquila [14] projects, respectively. They use aerial base stations mounted on high-altitude platforms-e.g., balloons-that fly several kilometers high. Balloons can drift slowly and provide basic network services to remote and rural areas. Instead, swarms of low-altitude drone relays can provide broadband connectivity and can reposition effectively on small time scales.
In [15], Zeng et al. study the use of aerial relays to relay traffic from two ground nodes whose links have been disrupted (due to big obstacles, environment or loss of infrastructures). Authors maximize throughput service and relay trajectories with an efficient algorithm that applies successive convex optimizations. In [16], Chen et al. extend this problem to the case where the relay network is offered by a swarm of multiple drones, and compare the effects of sending the traffic over one multi-hop link with using several two-hop links. In [17], Zhang et al. further focus on the multi-hop link of the aerial network to optimize the trajectories that maximize the end-to-end throughput and minimize transmit powers. These relay problems unveil the potentials of mounting relays on drones, and show clear use-cases for the applicability of such scenario. However, these problems are different to our approach, since they focus on the relay of traffic between two ground nodes. In contrast, we propose a drone relay-aided network that overlays the cellular network. Thus, we optimize the coverage service of users in a big region where a connected cellular infrastructure of ground base stations and mobile users exists, although it might become temporarily insufficient.
Moving drones to optimal locations to maximize various metrics, as coverage, is actively being investigated. Al-Hourani et al. [18] analytically model optimal altitude for one drone, to maximize coverage. Also, Hayajneh et al. [19] derive optimum drone altitude to minimize outages and biterror rate. Mei et al. [20] propose a decentralized inter-cell interference coordination scheme to maximize the weighted sum-rate of one aerial station and all users, as well as the uplink cell association over multiple resources. Mozaffari et al. [21] found the optimal location for multiple non-interfering drones to minimize the total transmission power. The same authors also analyze the performance of a singledrone-aided cell in the presence of underlaid D2D users, using stochastic geometry [22]. Chen et al. [23] optimize the location of drone relays to provide aerial caching for mobile users that connect to drones by means of millimeter wave links. Andryeyev et al. [10] estimate drone positions to increase cellular capacity by means of a self-organization scheme based on repulsion from base stations and drones, and attraction by users. They use conventional ground pathloss models for wireless channels, which differ substantially from the actual-and more complex-air-to-ground signal propagation model we use here.
Optimal location of multiple interfering drones was studied by Mozaffari et al. [24] and Rohde et al. [25]. However, differently from us, the former characterizes the network in the presence of only two drones, while the last obviates the effects of LoS channel conditions. Although they shed light on the nature of the problem, using such strong assumptions is limiting for real drone deployments. Conversely, in a companion work [26], we propose an optimization framework for cellular networks aided by a fleet of interfering coordinated drones, although the target optimization and use cases differ from the approach used here in important aspects, and so do problem formulation and heuristics. In [26], the network is optimized in terms of resource allocation, so to maximize throughput and fairness in the presence of significant crowds. That work does not focus on path optimization and does not consider backhaul limitations. Instead, here we maximize the amount of users covered for use cases in which providing network availability is the key performance indicator. Differently from [26] here we do not need to optimize resource allocation in detail, and can simplify the drone relocation problem by imposing bounds on the guaranteed user data rate. Nonetheless, we show that the problem remains NP-Complete. Similarly, Kalantari et al. [27] also optimize user association and fair bandwidth allocation in one cell with multiple drones. They show that interference does not allow to run more than three drones (per cell), and propose to obviate the problem by using advanced interference cancellation. This approach could be combined with our schemes, although it would increase complexity significantly.
Fotouhi et al. [11] propose a distributed algorithm for autonomous control of drones, and analyze the benefits of repositioning for spectral efficiency using straight paths. However, the literature does not offer yet any clever scheme to design drone paths to assist communication networks. B ezier curves [8] have been used to plan drone routes in [28] for military purposes, to smooth drone routes that have to fly over several check-points. However, B ezier curves have not been proposed yet to improve communications performance, as we approach for the first time in this paper.
A complete network architecture that could support a coordinated fleet of drone relays is still under design. Petrolo et al. [29] have designed ASTRO, a software-defined network for tetherless coordination of autonomous drones. Sundaresan et al. [30] have designed SkyCore, a network module integrated into an end-to-end network architecture called SkyLite. That is a complete network architecture for autonomous drone relay coordination, which demonstrates that operating an aerial network of drone relays is feasible provided the correct optimizations, like coverage maximization (which they do not approach). The maximum coverage framework that we discuss in this manuscript would fit in and be an asset for the above mentioned architectures.
We show in Section 4 that finding exact optimal drone coverage locations is an NP-Complete problem, which is non-linear and non-convex, thus intractable with solvers that provide strictly optimal results. However, we adopt an extremal-optimization scheme [7] that let us find near-optimal solutions requiring a much lower number of iterations.
Roughly speaking, such schemes focus on, at each time step, picking the "least fit" element of a discrete set and change its value to the "best fit" in order to improve a given utility function. In our case, the utility function that we want to maximize is the coverage, and the elements that we use to achieve such maximization are the positions of the drones. Since the position of a given drone affects which users other drones will cover, and also influences how the QoS is impaired by interference, extremal-optimization represents a suitable choice for our optimum drone placement and repositioning framework. At this point, we remark that although extremal-optimization is a form of optimization originally introduced to be used in a static manner, in our work we use it to design a deployable algorithm for the location of drones that works in a dynamic manner.

Reference Scenario
We consider a ground surface S administrated by the ground network consisting of a set G of ground base stations, as shown in Fig. 1. 1 We refer to ground base stations as gNBs, using the new 3GPP jargon for next generation base stations (BSs). In the region S, the ground network provides service to a set U of user equipments (UEs), i.e., mobile users. In order to increase coverage, we consider that coverage assistance is provided by the presence of a finite set D consisting of D drone relay stations. Each drone is equipped with a mobile relay that gives access to UEs on an orthogonal downlink bandwidth with respect to the gNBs band. We refer to drones as aerial base stations (ABSs).
We assume that gNBs provide backhaul connectivity to ABSs over the reuse of the downlink spectrum used for gNB-UEs access links. Current gNBs provide cellular coverage through three sectors pointing mainly to the ground. We assume that in order to set backhaul gNB-ABS links, gNBs have an additional full dimensional antenna array that performs 3D-beamforming over clear LoS links, 2 as suggested and studied in [31]. Therefore, access links gNB-UEs and backhaul links gNB-ABS do not interfere. Furthermore, gNBs equipped with this kind of antenna array are able to perform 3D-beamforming to several relays, and alternate transmissions over time slots on a millisecond scale. Hence, each gNB g can provide backhaul service to a limited number of ABSs, namely D g .
The coverage of each ABS is an irregular ground area that depends on the drone height, cell environment and interference from other ABSs. The interference among ABSs directly affects the signal-to-interference-plus-noise ratio (SINR) that the ground users receive. The SINR 1. New SDN tools designed to manage large networks are able to coordinate ground base stations to perform any network optimization [2].
2. Based on the receiver location or instantaneous channel state information, 3D-beamforming allows to build directional beampatterns that focus the transmission energy on the direction where the receiver is. This flexible technique helps to mitigate interference so to provide higher rates. 3D-beamforming is very useful for backhaul wireless links from one source to few relays, as in an aerial backhaul network.
depends on the air-to-ground path-loss model, which is based on the link LoS probability between drones and users. As described later, such path-loss model clearly differs from the conventional attenuation models used in ground cellular networks. Indeed, such a LoS-based path-loss and interference model for the communications channel provides a multi-ABS coverage framework for aerial networks, which is radically different from conventional frameworks for ground networks-as e.g., ground D2D networks [5]-and whose characteristics we study. In this framework, we consider that a gNB g and an ABS d can realistically serve a limited number of users, namely U g and U d , respectively. We further assume that channel bandwidth is equally split among the users that a BS serves, although more sophisticated schedulers could be easily adopted in the analysis.
With the above, we aim to find optimal locations for D drones, so to maximize the number of users covered by gNBs and ABSs with a guaranteed bandwidth. Besides, we identify two additional problems to support fleet repositioning: (i) deciding which drone flies to which position upon an optimization update and (ii) designing flight routes.

Air & Ground channels: Path-Loss and Interference
We assume that the surface S is flat, so that the position of a UE u 2U is taken as an input and denoted by p u ¼ ðx u ; y u ; 0Þ. The position of a gNB g 2G is denoted by P g ¼ X g ; Y g ; h g À Á . The positions of all ABSs d 2D are the decision variables of the coverage problem, and denoted by For all drone d 2 D, and for all user u 2U, the horizontal distance between u and the ground projection of d is r d;u ¼ kðX d ; Y d Þ À ðx u ; y u Þk. The elevation of d is h d . Due to the low altitude of drones-few hundreds of meters at most-the channel conditions of communications between a serving drone and an end-user are much affected by the LoS. Depending on whether the access link is free of obstacles (like buildings, traffic, etc.), the attenuation differs considerably [32]. Thus, the air-to-ground path-loss among ABSs and UEs depends on the probability of LoS, which is a complex function of the elevation angle between user u and drone d, according to the following expression: where a, b are parameters depending on the environment, i.e., number of buildings and large signal obstructions per unit area, building's height distribution, ratio of built-up area and clean surfaces, etc., as it has been derived in [18], based on the ITU recommendations [33]. In Eq. (1), the elevation angle (in radians) appears as u d;u ¼ arctanðh d =r d;u Þ. As u d;u approaches p 2 , i.e., when the drone d hovers just above the user u, the probability of LoS reaches its maximum value. In Fig. 2, we see that the positions of the drones directly affect blockage conditions of the ABS-UE acess links. Thus, the higher a drone hovers, the more likely is to have LoS. However, the strength of the signal gets also attenuated with the distance. For single-drone missions, there is an optimal altitude that maximizes coverage [18]. However, in a multidrone scenario as the one we discuss in this paper, the effects on interference from other drones are a key additional issue to consider, one that makes the optimal drone hovering altitude depend on the positions and elevations of the rest of the drones. This also precludes the possibility to straightforwardly apply single-drone mission approaches to multidrone scenarios, since the former are not designed to account for interference, as in [9].

Ground-to-Ground Channel
Ground cellular links, i.e., gNB-UE access links, operate over an OFMDA channel with access bandwidth W A . Hence, users scheduled by the same gNB do not suffer intra-cell interference. However, ground users may enjoy poor QoS due to the presence of inter-cell interference, from other close gNBs. The path-loss of these channels is based on large-and smallscale fading, as widely studied in literature [34] and shown in Table 1. We denote the SINR of ground access links ðg; uÞ as G A g;u and impose that its minimum user access rate, i.e., 1 Ug W A log 2 ð1þG A g;u Þ, is above the guaranteed rate R A min .

Ground-to-Air Channel
In order to provide backhaul connectivity to ABSs, gNBs perform 3D-beamforming over clear LoS links. Hence, the attenuation that a signal from gNB g to ABS d suffers is where h B % 2 is the path-loss exponent in free-space LoS links, f B is the frequency of backhaul wireless links in Hz, c l is the speed of light in m/s and N s B is a random gaussian variable with zero mean and s B standard deviation, modelling the effects of slow fading and shadowing. 3D-beamforming builds antenna patterns that radiate much of the energy over a main lobe with a half-power  beam-width (HPBW) that may be wide, hence incurring high interference to other ABSs in LoS. Also, the formation of directional beam-patterns comes with the presence of sidelobes with non-negligible radiating power, that also incur (low) interference. Thus, depending on the radiating angle of other gNBs, a backhaul link may enjoy better or worse QoS, due to the presence of interference. 3 We denote the backhaul SINR of link ðg; dÞ as G B g;d , which is equal to where P g Tx is the transmission power of g, G g is the antenna gain over the main lobe of the beam-pattern of g, N d is the thermal noise, and most importantly, I B g;d is the actual interference that ABS d suffers from any other gNB, depending on the angle of their beam-patterns, i.e.
where f g 0 ;d is the angle between the direction of the main lobe of the antenna of g 0 and the position of ABS d. We impose that the minimum backhaul rate, 1 D g W B log 2 ð1þG B g;d Þ, is above a rate R B min . W B is the backhaul channel bandwidth.

Air-to-Ground Channel
While ground cellular links suffer from conventional attenuation based on fast and slow fading, the path-loss of an ABS-UE link ðd; uÞ differs notably and is affected by an excess attenuation, depending on the LoS likelihood presented in Eq. (1). The average attenuation is derived in [18] as where LoS ; NLoS are the excess path-loss components in LoS and NLoS connections respectively, h A ¼ 2 is the path-loss exponent and f A is the carrier frequency in Hz. As surveyed in [35], we have based the air-to-ground path-loss on the average large-scale fading in order to perform the ABS-UE association in the optimization. Nevertheless, we further consider fast and slow fading, modelled as log-normal and Rayleigh distributions, respectively, in the implementation of the framework when we analyze the results. We denote the access link SINR as G A d;u , which is equal to where P d Tx is the transmission power of the antenna integrated in the ABS d, N u is the thermal noise, and most importantly, I A d;u is the actual interference that user u suffers from any other ABSs. Thus, the interference depends on the 3-D position of the rest of the ABSs, i.e.
We impose that the minimum rate, 1 where W A is the access channel bandwidth. In our system model, ABSs operate in orthogonal bandwidth with the cellular band. Thus, there is no interference between cellular users and drone-served users, which is commonly the main limiting factor in aided cellular networks, as e.g., inband D2D networks [5].
In Table 1 we gather the path-loss model used for each kind of channel, where f G is the band used for gNB-UE access links and is equal to the gNB-ABS backhaul links, i.e., f G ¼ f A , and R & A is a random Rayleigh variable with scale parameter & A . In Table 2 we summarize the interference suffered in each of the channels, as it has been described along this section. We have shadowed the table cells that imply presence of interference.

MULTI-DRONE COVERAGE FRAMEWORK
We aim to find optimal 3-D positions for a fleet of D drones in which the number of UEs under network coverage is maximum. The optimization is run at regular time intervals, considering every time the users as static, so that static drone positions solve the coverage problem. The coverage maximization provides a set of 3-D coordinates where drones have to fly during the time interval, provided a drone d has to fly towards a reachable destination, i.e., a point within the ball S d of radius given by the drone speed times the duration of the optimization update interval. However, the output of the optimization does not necessarily coincide with the assignment of fleet destinations that also minimizes flight time. Hence, we will solve the assignment of fleet destinations as a secondary problem, given the optimal coordinates found by the primary problem.
Although it is possible to formulate the coverage maximization and the minimum flight time assignment in one optimization problem, we believe that it is more clear to decouple both problems and solve them separately, while having the same optimal solution. This is due to the fact that the set of optimal drone positions is not changed by solving the secondary problem and at least one feasible solution

Ground-to-Ground
Ground-to-Air Air-to-Ground

Ground-to-Ground
Inter-cell interference Directional re-use Orthogonal bands Ground-to-Air Directional re-use Low interference: 3D-beamforming Orthogonal bands Air-to-Ground Orthogonal bands Orthogonal bands Inter-drone interference 3. In general, since beamforming builds antenna patterns with directional main lobes, interference remains low for non-aligned ABSs.
exists for the secondary problem, which is the output of the primary problem. Thus, we first present the optimal aerial coverage (Section 4.1) and then the assignment of fleet destinations (Section 4.2).

Optimal Aerial Coverage
Coverage Problem C. Given a fleet D of drone relays in a cellular network managed by a centralized orchestrator, U ground users, a height range ½h min ; h max for the ABSs, guaranteed coverage rates R A min , R B min for access and backhaul channels respectively, and a maximum number of users U g and U d that gNBs and ABSs can cover, find the optimal 3-D positions P d ¼ðX d ; Y d ; h d Þ of drones so to maximize the amount of users covered by ground-and drone-cells.
Since the positions of drones, including their heights, affect the shape of the covered regions, we can mathematically formulate the Coverage Problem C to search for optimal values of the continuous decision variables P d ¼ðX d ; Y d ; h d Þ 2 R 3 that maximize the number of users under network coverage. Denoting by C b;u the binary variable that takes value 1 if BS b 2G[A covers user u and 0 otherwise, and by B g;d the binary variable that takes value 1 if gNB g provides backhaul service to drone d and 0 otherwise, the formulation of the Coverage Problem C is Backhaul network constraints : 1 Access À backhaul constraints : Drone air À space constraint : (8) Access Network Constraints. The first constraint guarantees that the access link rate between BS b and user u is above R A min as soon as u accesses the network via b; the second constraint tells that a user cannot be covered by more than 1 BS; the third constraint accounts for the number of users that each BS (either gNB or ABS) can serve.
Backhaul Network Constraints. The forth constraint guarantees that the backhaul link rate between gNB g and ABS d is above R B min as soon as d connects to the network via g; the fifth constraint tells that each ABS cannot connect to more than 1 gNB; the sixth constraint limits the number of drones that each gNB g can serve to a maximum of D g drones.
Access-Backhaul Constraints. The seventh constraint states that a user u can connect to a drone d only if d is under the coverage of some gNB. Hence, each drone that provides network access to at least one ground user is connected to the network via one backhaul link, so that, every ground user served by a drone is indeed attached to the cellular network.
The eighth constraint states that the number of users covered by those drones attached to the same gNB g is limited by the maximum capacity of users U g in g.
Drone Air-Space Constraint. Finally, the air location of a drone d has to be within a 3-D region S d SÂ½h min ; h max which can be reached in the time interval used for optimization, depending on flight speed and current drone position.
Feasibility. The optimization problem in Eq. (8) is always feasible. For instance, consider a general instance of the problem. Take a random position for each ABS d inside their reachable regions S d . Now set all binary variables fC b;u }, fB g;d g equal to 0. Since the SINR functions G B g;d , G A b;u are always positive (see Eqs. (3), (6), respectively), the solution satisfies all the constraints. This solution provides a utility function of 0 users covered, but it is feasible.
Complexity. Problem C is NP-Complete because, as shown in Appendix a, the well-known NP-Complete minimum geometric disk cover (MGDC) problem [36] can be reduced, in polynomial time, to a particular case of Problem C.
The first constraint on the user access rate, when b 2 D, is non-linear and very complex (also the forth constraint), since it depends on the air-to-ground path-loss shown in Eq. (5) for one link, but also for the interfering links from other drones. To make the constraint more visual and remark its non-linearity and complexity, we develop its expression for a drone d as follows: (9) where the continuous variable r d;u ¼ kðx u ; y u ÞÀðX d ; Y d Þk is the distance between user u and the ground projection of drone d, and are constant. In Eq. (9) we see that this constraint depends on the position decision variables ðX d ; Y d ; h d Þ not only as an attenuation from the distance, but they also affect the LoS probability, as shown in Eq. (1).
Unlike previous works like [3], [9], [37], in which the drone-service condition is based only on the attenuation or the Signal-to-Noise Ratio (SNR), in (8) we have formulated a novel 3-D drone placement optimization that accounts for the actual inter-drone interference suffered by ground users.
Eq. (8) represents a Mixed-Integer Non-Convex Program (MINCP), which is not tractable with currently available optimizers dealing with problems that are, at least, convex. Since problem (8) presents a non-convex formulation mainly because of the attenuation depending on the LoS probability, we cannot apply any off-the-shelf optimizer to optimally solve this problem. In addition, the problem itself is NP-Complete, so we resort to a heuristic, as detailed in Section 5.

Assignment of Fleet Destinations
The second problem to solve when users move and drones have to be repositioned is an assignment problem. Since it does not matter which ABS goes to which destination (as long as such destination is reachable), we enforce each drone to fly towards the positions that minimize the aggregated flight-time by the fleet. Formally, we dispose of a fleet of D drones that must fly from source positions fp d g d2D and reach target positions fP d 0 g D d 0 ¼1 . Thus, we formulate the following assignment problem: where we have introduced the binary variable F d;d 0 2 f0; 1g to denote whether drone d flies from p d to P d 0 or not. Tðp d ; P d 0 Þ is the assignment weight, and depends on the time that drone d needs to fly from source position p d to the destination P d 0 . The equality constraints ensure that each destination P d 0 is reached by only one drone d, and that each drone d reaches only one destination P d 0 . The utility U fly is used to minimize the flight time of the fleet of drones.
For simplicity, we assume that drones d fly at a constant speeds of v d . Thus, the time needed to fly from p d to P d 0 is However, each drone d can reach only those destinations with coordinates within their reachable region S d . Hence, we impose infinite required time for each drone d to reach all destinations P d 0 that fall outside S d Hence, the optimization in Eq. (10) assigns to each drone d one destination P d 0 such that d is able to reach P d 0 and the aggregated flight time is minimized. Feasibility. The optimization problem presented in Eq. (10) is always feasible and finite. For instance, assigning to each drone d the destination P d ¼ ðX d ; Y d ; h d Þ obtained as an output of the optimization problem of Eq. (8) provides a (finite) feasible solution.
Complexity. The optimization problem in Eq. (10) is a special case of mixed-integer linear program (MILP), one that can be solved efficiently in polynomial time through the Hungarian method [38], with complexity OðD 3 Þ. Therefore, this second problem related to dynamic networks is easy to address optimally in low-degree polynomial time.

DYNAMIC DRONE REPOSITIONING ALGORITHMS
So far we have discussed the optimization of the placement of a fleet of ABSs hovering over a ground cellular network (see Section 4.1), and the optimal flight assignment that minimizes the flight time (see Section 4.2), thus overlaying a legacy cellular network managed by an orchestrator. Nevertheless, since users move, the optimization is reconsidered periodically, with updated drone air-space constraints. Repositioning has to be run frequently, so that we need efficient heuristics. Next, we describe how to practically implement the optimization framework described so far.

OnDrone: An Algorithm Suit for On-Demand Drone Coverage Optimization
Coverage Problem C is NP-Complete, thus optima cannot be reliably solved on-demand for fast placement of drones, which is key for dynamic repositioning cases. Thus, even if the problem was optimally solvable, the need of having an efficient heuristic would remain. To this aim, we propose here an On-demand Drone Coverage (OnDrone) algorithm, based on an extremal-optimization algorithm (EOA) that runs in polynomial time. For benchmarking purposes, we further consider state-of-the-art proposals from [9] and [10].

OnDrone for Multi-Drone Coverage
EOAs are evolutionary algorithms that restrict the search space and achieve near-optimal results in polynomial time [7]. EOAs are based on a fitness metric and, at each step, try to improve the configuration of the element of the system that yields the least contribution to the fitness metric. Therefore, EOA's principles perfectly match the nature of the coverage problem addressed. Specifically, the fitness function is the number of covered users and the least significant contribution comes from the drone that covers the least number of users. Such drone may be either far from users, where its transmissions are severely affected by the interference coming from the rest of ABSs, or in a position where the backhaul service is low. Thus, it is convenient to reposition such drone and increase the coverage. The search space for drone locations is restricted to a lattice, as shown in Fig. 3. We derive a cylindrical lattice that contains the ground network surface S composed by the positions over which OnDrone moves the ABSs to get the best coverage utility. Those lattice points that fall outside the ground region S are discarded, so the design of the lattice applies to any shape of S. We split the base into a grid of equal areas (see Fig. 4), and the height into equidistant altitudes. This leads to a cylindrical lattice of equal volume subspaces. Specifically, we divide the square of the radius  and the angle of the base grid in N r and M u equal pieces. In polar coordinates, the resulting points and base grid are where Dð0; R C Þ is the closed disk in R 2 centered at the origin and of radius R C (where R C is big enough to make Dð0; R C Þ contain S). In this way, the base area is divided into N r ÂM u regions with the same area Fig. 4). The height of the cylinder is divided into H equidistant segments in the interval between minimum and maximum drone hovering height, h min and h max . In cylindrical coordinates, the resulting lattice is The proposed EOA-OnDrone-begins with an initial feasible and suboptimal (random) implementation of the system. Then, it updates the positions of the ABSs providing worst individual contribution to the full performance, i.e., the drone covering less users. At each iteration, the "least fit" ABS is selected and moved to a reachable position where the system utility increases as much as possible, considering the coverage by ground-and drone-cells. Also, OnDrone provides a new position where the drone can be attached to a gNB that provides backhaul service with the guaranteed QoS. To decrease the probability of finding only local optima, we consider that if the worst performing drone cannot be moved to improve coverage, then we try to reposition the next worst-performing drone. OnDrone keeps moving the ABSs with lowest contribution until it does not find any better location for any ABS, or it reaches a maximum number of iterations i 0 2 N. The optimality of this algorithm is studied in Section 6.
Complexity. Algorithm 1 reports the pseudocode of the proposed OnDrone, in order to target maximum ground-plus-air users coverage, thus approximating the optimal solution of Coverage Problem C. The complexity of OnDrone can be evaluated as follows. At each iteration, one drone d 0 2 D is selected and repositioned. Such drone d 0 selects the 3-D position P d 0 in the lattice L (see Eq. (13)) at which gNBs and ABSs cover more users as long as there exists g 2G providing the guaranteed backhaul QoS in that position. A user u is covered by an ABS d only if the user rate experienced is greater than minimal user access rate R A min and d is covering at most U d users, so that u enjoys the guaranteed QoS. Thus, the signal strength from d 0 and from the rest of the drones in Dnfd 0 g must be checked. This means that the complexity of each iteration is OðjLj Á D Á UÞ, where L is the size of the lattice. Since i 0 is the maximum number of iterations needed for the algorithm to converge to a solution, the complexity of OnDrone is Oði 0 Á jLj Á D Á UÞ, where i 0 is constant and can be omitted.
We remark that, unlike the NP-Complete problem presented in the previous section, OnDrone requires few iterations. Indeed, we have performed all the experiments in this paper in few seconds in a personal computer, on average. As a matter of fact, OnDrone is intended to be used in an on-demand fashion, dynamically repositioning drones to adapt to user's moves over time. OnDrone is then practical and can be executed at the network orchestrator.
If the same coverage remains when placing d 0 at P d 0 , go back to step 3 ignoring unsuccessful d 0 's. 6: Place d 0 at P d 0 and set P fP d g d2D . 7: Go back to step 2 until: Ground-plus-air coverage is no longer improved. Maximum number of iterations i 0 is reached.

Seq: Sequential Multi-Placement
In addition to OnDrone, we have also developed a simple heuristic that finds feasible solutions to the Coverage Problem C in polynomial time. We base this algorithm on the Efficient 3-D Placement-hence the name of the algorithm-scheme derived in [9], and thus we adapt it to Sequential Multi-Placement (Seq) in order to support aerial networks with more than one ABS. In [9], the authors model the presence of one drone providing coverage in one single cell, i.e., they maximize coverage for single drone missions. They model a circular drone footprint, which allows to formulate a convex optimization problem due to the presence of only one drone. We adapt the proposal in [9] to place ABSs one by one, according to realistic interference metrics not originally considered in that work. Seq will be used as a benchmark for OnDrone.
We build the Seq heuristic by induction as follows: Since the framework proposed in [9] only considers drone-coverage, we first compute the amount of users covered by groundcells. Second, we select one ABS and maximize coverage as described in [9] (i.e., using Efficient 3-D Placement) for the remaining non-covered users, namely U 0 . Here, the placement space for drones is restricted to those positions where at least one gNB provides backhaul connectivity with the guaranteed QoS. We denote as U 1 to the set of users that are covered by the first selected ABS. For this first ABS, there are no interference issues. Let i > 1 and assume that we have located iÀ1 ABSs and that we want to locate the ith ABS. Assume that U i k are the sets of UEs that each previous i k th ABS covers at the moment of its placement. Then, Seq finds a 3-D position at which the ith ABS covers more users from the set U 0 n S 1 i k < i U i k and at least one gNB provides the requested backhaul QoS. Thus, the ith ABS aims to cover the maximum number of users that are not covered yet.
Seq ends when the Dth ABS is located. After this, the algorithm computes the actual number of users served according to interference (in both the backhaul and access network). Hence, Seq has the same objective as OnDrone: covering the maximum number of users according to QoS guarantees. We report the pseudocode in Algorithm 2. At each iteration i > 1, the ith ABS sequentially selects the best position for it based on Efficient 3-D Placement from [9] maximizing its own coverage, no matter how its position affects the coverage of the remaining ABSs or the already set backhaul links. Since placing the new ith ABS adds interference to the system, the previous iÀ1 dronecells shrink, and cover less UEs than the ones originally intended by the Seq choice.

Algorithm 2. Seq: Sequential Multi-Placement
Complexity. The complexity of Seq is evaluated as follows. Seq makes D steps, one per each drone. At each step d, Seq defines the final position of drone d 2D. In [9], the authors do not propose any algorithm for placing the drone, but they solve a convex mixed-integer non-linear program with a convex optimizer. Such optimizer does not run an algorithm with polynomial-time complexity. Instead, it performs a combination of interior-point methods with a branch&bound search [39]. Hence, we opt for approximating their problem through a search on the lattice L. As in OnDrone, checking whether a user is covered requires to check the signal strength of the serving drone along with the interfering drones, i.e., D signal strengths. The complexity of this process is OðjLjÁDÁUÞ. After the last iteration, the algorithm checks which users are actually covered, since those users that at some iteration i were covered, may no longer enjoy the guaranteed QoS because of the repositioning of other drones in successive iterations. This check has a complexity of OðD 2 ÁUÞ. Thus, the complexity of the Sequential Multi-Placement algorithm is OðjLjÁDÁU þ D 2 ÁUÞ, that is similar to the complexity of OnDrone because jLj ! D, since drones cannot be co-located and so the number of possible distinct drone positions cannot be smaller than the number of drones (indeed, in a well designed system, jLj ) D).

RA: a Repulsion-Attraction Scheme
Here we briefly describe the Repulsion-Attraction (RA) scheme derived in [10]. RA is a multi-drone placement scheme in which several self-organized ABSs dynamically change their position to track clusters of users. The approach is based on the assumption that ABSs will be attracted by the presence of users in the ground, and will be repulsed by gNBs and other ABSs in order to avoid interference.
Complexity. RA consists into maximizing a non-integer metric without any constraints. Hence, one can apply standard methods like line-search or trust-regions methods for unconstrained optimization. Such methods have low complexity and their performance depends on the target tolerance on the error. Moreover they converge really quickly [39], with a linear dependance on the number of iterations i RA , i.e., the complexity is Oði RA Þ.

B ezier Flight Routes
The last problem to solve consists in designing drone's trajectories. The output of OnDrone (Section 5.1), and the Hungarian method (Section 4.2), provide the source and destination for each drone carrying an ABS. Therefore, we now design a route optimization scheme. While drones fly, both the backhaul and users association change. Initially, each ABS attaches to a gNB with the required QoS and starts the flight. Once the backhaul QoS level is low due to long distance or interference impairment, the ABS sets a backhaul link with a new gNB with the guaranteed QoS. Similarly, while a drone flies, UEs attach to that ABS in case that the minimum access data link rate (i.e., QoS) is guaranteed. Upon arrival to the destination, the association is already optimal in terms of coverage.
On the one hand, drones have high aerial mobility and fly over a ground cellular network in a 3-D space, without many restrictions of walls, streets or vehicles. On the other hand, drones hovering over regions with good QoS from some gNBs or under-populated regions with a big surface may lead to under-utilized ABSs and low coverage, depending on the topology of users. Furthermore, if a drone is not fast enough, it might occur that when the drone arrives at the destination the users topology has changed too much so the destination is no longer optimal, and the network needs to be reoptimized. To avoid such undesired effects, and knowing that drones may have to be redirected while flying towards a destination, we propose drone paths following B ezier curves [8], instead of commonly assumed straight lines, as adopted in [11]. Indeed, using B ezier curves allows to deflect drone trajectories towards areas with higher user density, so to enhance drone coverage and enable unique coverage opportunities while drones seek their optimal position. Since we leverage on the notion of B ezier curve, Appendix B provides some background on the subject.
In our proposal, we define the set of anchor points for our B ezier-based flight path and use the standard de Casteljau algorithm [8] to derive the B ezier curve corresponding to the selected anchor points (see Appendix B for details).
We obtain the set of anchor points inductively, as detailed next. Let p d and P d be source and destination of a drone d 2 D, let v > 0 be the width for the two-sided offset region of the curves and let B > 1 be the maximum number of anchor points for the B ezier curve. B is determined as the density of users covered per drone-cell, since in case a drone cannot cover more than B users, it does not make sense that such drone wishes to deflect its path attracted by more than B users. We take v ¼ 2R d , where R d is the maximum range at which drone d can provide coverage in its optimal position. We define as the initial set of anchor points P 0 both nodes: P 0 ¼ fp d ; P d g. Thus, we define the B ezier curve b P 0 ðtÞ for P 0 , which is the segment joining p d and P d , computed with the de Casteljau algorithm [8]. Now, we iteratively modify the current B ezier curve until we derive the final B ezier path.

Algorithm 3. B ezier Scheme
Given a curve bðtÞ, we denote its length as bðtÞ ð Þ, and take ðb P 0 ðtÞÞ as a reference length for the final B ezier curve. Indeed, we build a B ezier Scheme such that the obtained curve is not longer than t ¼ ð1þaÞ Á ðb P 0 ðtÞÞ, for a given a > 0. a is determined as the fraction of the time interval in which a drone would have already arrived to its destination in case of following a straight path. Hence, we make sure that a drone reaches its destination but has the flexibility to deflect its path to improve coverage. Let S v ðbðtÞÞ be the two-sided offset region of width v that results from stroking a curve bðtÞ, and let U P 0 ¼ S v ðb P 0 ðtÞÞ \ U be the set of users that are inside the offset region of b P 0 ðtÞ. We denote by g P 0 u the gravity of a user u 2U P 0 and define it as the density of users in a disk centered in u with radius v=2 The first user u 1 selected as anchor point is the one with highest gravity in U P 0 such that the resulting B ezier curve is not longer than t ¼ ð1þaÞ Á b P 0 ðtÞ À Á , i.e.
Now, P 1 ¼ P 0 [ fu 1 g defines a new B ezier curve, b P 1 ðtÞ. The sorting of positions in the sets P k is important, since each order defines a different curve. Thus, we sort the points in increasing order according to the distance to the source, p d . Finally, U P 1 ¼U P 0 \ S v ðb P 1 ðtÞÞ is the updated set of users for next iteration. Then, we keep selecting anchor points for the B ezier curve while they exist, until the length of the B ezier curve exceeds t or the maximum number B of anchor points is reached. At each iteration k, we build a new curve from P kÀ1 , U P kÀ1 and b P kÀ1 ðtÞ.
Complexity. We show the derived B ezier Scheme in Algorithm 3. Its complexity can be evaluated as follows. There is an iteration of the de Casteljau algorithm to derive each B ezier curve b P k ðtÞ. The complexity of the de Casteljau algorithm is quadratic with the number of anchor points of the B ezier curve that it builds [8]. The B ezier Scheme runs at most B iterations, and at each iteration k it uses k þ 2 ¼ jP k j B anchor points.
Thus, the complexity of the B ezier Scheme is OðB 3 Þ, where B is the maximum number of anchor points allowed.

Overall Complexity
In Fig. 5 we show the flow chart of the proposed optimization framework for drone-aided dynamic cellular networks. As we have fully detailed in the paper, first we solve, through OnDrone, the framework for maximizing network coverage in polynomial time. OnDrone outputs the optimal 3-D positions at which drones must locate, given their previous positions. However, drones do not move instantaneously. The Hungarian method solves then the assignment problem that minimizes the distance flown in polynomial time, and outputs the source and destination of each drone concretely. Finally, our B ezier Scheme designs the flight routes of drones in order to deflect their paths towards clusters of users and increase the coverage efficiency over the time interval. This is another polynomial time algorithm. Therefore, the overall optimization framework runs in polynomial time. The overall complexity of our framework for fast repositioning of drone-aided cells is the result of summing the complexities of the schemes that compose such framework, which is shown in Fig. 5. When the 3-D placement optimization is performed by OnDrone, the complexity is OðjLj Á D Á UÞ þ OðD 3 Þ þ OðB 3 Þ ¼ OðjLj Á D Á UÞ þ OðB 3 Þ, since jLj ) U ) D.
In Table 3 we summarize the complexities of the detailed algorithms, including the overall optimization framework.

Orchestration of the Optimization Framework
The possibility to integrate the above described framework into a communication system raises several alternatives regarding where the different algorithms can be executed. First, there is the issue of deciding where the drones must be located at each time instant. This task (i.e., to execute OnDrone for drone placement) should be assigned to a device on which all drones have direct communication and therefore can easily know what must be their destination. Therefore, the gNBs seem to be an appropriate place for the orchestration of drones' positions. In alternative, and considering current trends in 5G networks architecture design, the orchestrator can be a software slice in the MEC, which is the edge-cloud computing platform of future 5G cellular networks and which resides just next to base stations [40].
The second issue consists in designing drone's trajectories, provided it must change its current location. However, the B ezier curves used for the flight routes have been already  designed in a discretized manner with the de Casteljau algorithm, using short straight segments to build a B ezier curve. At this point, we note that currently, most drones are already capable of autonomously travel to concrete positions following straight lines (see for instance [41]). So, drones can follow such short segments (using their originally integrated traveling mechanisms) without significantly deviating from the flight route provided by the B ezier scheme. In that case, those responsible for such tasks (i.e., obtain the discretized B ezier curves and follow the corresponding straight segments) are the drones themselves.

EXPERIMENTAL RESULTS
Here we numerically assess the performance of our multidrone optimization framework. We assess the coverage offered by the network, with the assistance of a fleet of drones, in a circular surface. We compare optimal placement results yielded by OnDrone (presented in Section 5.1.1) with the ones obtained with Seq (based on [9] and described in Section 5.1.2) and with the RA scheme (from [10], and described in Section 5.1.3). We also compare the results with the optimal achieved by means of Monte Carlo simulations, since computing exact optima is not doable in networks with as few as a fistful of UEs, gNBs, and drones. Besides, we compare our scheme to a modified OnDrone that neglects interference (referred as "iNeg" in the figures). Hence, we assess the importance of introducing interference in the analysis. We also compare coverage results while drones reposition following the B ezier Scheme (presented in Section 5.2) with a simpler Straight Scheme, as adopted in [11], which consists in moving drones over straight paths towards a certain destination identified with the Hungarian method (see Section 4.2). Since RA has been designed in [10] to dynamically track UEs, we compare also the B ezier and Straight Schemes when the path planning is based on RA.
We mainly study the placement and repositioning of drones in synthetic and realistic scenarios, and the effects of interference and LoS on user coverage over time, under four classes of environmental scenarios: suburban, urban, dense and high-rise. These four environments correspond to different densities of elements (e.g., buildings) that affect the LoS probability. Moreover we study three distinct cases of deployment scenarios, in which the location of users follows different distributions: PPP: We place UEs in a circular surface, according to a Poisson point process. Cheese: We define a surface that includes certain regions which are not accessible to UEs, and locate UEs uniformly random in the rest of the network.
Then, we have a surface with empty areas (like in a Swiss cheese). This users' distribution is typical of public gardens or areas with restricted zones. Capital: We also run our numerical evaluation over a simplified map of the center of a dense capital city, Madrid, considering main zones of people affluence and no users in indoor installations. In the PPP and Cheese scenarios, we locate gNBs according to the same distribution. Heterogeneous synthetic distributions would be of interest for traffic demand-based optimizations [42], which is out of the scope of this paper, and we therefore do not consider them. However, in the Capital scenario, we consider the actual locations of those gNBs that belong to the main network operator in the city. In all experiments, the number of iterations required by OnDrone was very small (see Appendix C.1). Table 4 gathers the parameters for the air-to-ground pathloss model and the probability of LoS described in Eqs. (1) and (5), depending on the density of the environment that we consider in the numerical simulations. These parameters are obtained based on the number of buildings and large signal obstructions per unit area, building's height distribution, ratio of built-up area and clean surfaces, etc., as it has been derived in [18], based on the ITU recommendations [33]. Such parameters allow to differentiate the main four environmental conditions. Table 5 gathers the parameters that, unless otherwise specified, we have used for the network model, regardless of the simulation environment. We take a circular surface of R C ¼ 1:5 km of radius where there are 10 gNBs, and a height range between h min ¼ 60 m and h max ¼ 600 m for drones. This is a generally doable height range, 4 since lower values would be too close to ground (and, e.g., vehicles or even people) and higher elevation would be affected by high-speed winds which are unsafe for an aerial network of simple drones. However, in our numerical evaluation the actual maximum  drone altitude is rather determined by the environment density. For instance, in a high-rise environment, although high altitudes increase the probability of LoS, the attenuation is much stronger, so that drones need to fly closer to the ground. In contrast, the suburban or urban environments do not suffer strong attenuation, so that drones can fly higher. However, a high altitude turns into links with higher LoS probability, thus yielding more interference for far ground users. The power transmission from ABSs is 10 dBm, as adopted in [43], [44], [45], which is notably lower than the usual 44 dBm used for gNBs in the ground network (as we adopt). This is because ABSs have much higher LoS probability than ground base stations and do not use omnidirectional antennas, and hence require much less power. Using higher ABS transmission power, as 25 or 44 dBm, may provide better coverage due to better signal strength, although also provides less resilience to interference impairment, as we discuss in our results. Moreover, ABSs are carried by flying drones, which spend most of their energy into hovering, flying towards desired positions at a given speed, and carrying the weight of the communication equipment. This poses serious constraints on transmission power, as evaluated in [43]. Hence, we have chosen to use a 10 dBm of power transmission for analyzing our framework and algorithms. The guaranteed user data rate is 0.72 Mbps, which guarantees video streaming with 240p resolution [46], and allows 360p and 480p resolution in many cases over the MPEG-4 standard. It also allows adequate videoconferencing quality using video compression [47]. Such guaranteed rate, with customary 20 MHz bands and assuming that no more than 100 users can attach to a BS, corresponds to guaranteeing that the SINR is higher than g A ¼ 10:9 dB, according to the Shannon capacity formula (see Appendix C.2 for a discussion on the minimum data rate experienced under coverage). Besides, such SINR value allows for a 16QAM modulation and a coding rate of 1=2 in LTE communications, as derived in [48], although we also study the impact of using other guaranteed rates. The carrier frequency used by gNBs (either for backhaul channels between gNBs and ABSs or access channels between gNBs and UEs) and for the access channel between ABSs and users are 1815.1 MHz and 2630 MHz, respectively. These channel parameters are as in the LTE/ LTE-A network provided by Movistar in Madrid. As it has been discussed in [31], the antenna patterns of gNBs for backhaul connectivity are directional with a HPBW of 65 degrees, according to the 3GPP technical specifications [49]. For dynamic experiments, we slot time into intervals of length T ¼ 60 s. This means that every minute, we reoptimize the network by means of OnDrone (or with any of the benchmark schemes, as Seq or RA schemes) and immediately re-direct drone flights for dynamic repositioning using the solution of our assignment problem and a flight route computed with either the B ezier or the basic Straight Scheme from [11]. Users move according to the random way-point model (RWP) 5 (see [50], [51]) with an average speed of 2 m/s. We update user's positions every second, while drones fly at a constant speed of 15 m/s over a continuous path. Moreover, such flight speed allows for low drone energy consumption, according to the energy consumption model for aerial aircrafts derived in [43].
We use MATLAB R2018a to simulate channel conditions and mobility of drones and UEs. For small networks, we have performed Monte Carlo simulations to search for the optimal drone placement. For each network instance, we have run 10 7 random positions for the fleet and taken the settings for the best coverage. Such simulations are very timeconsuming even for small networks, yet we have observed that such number of Monte Carlo runs per instance provides a coverage output that is hard to improve, since the output placement provided remains invariable for a large number of additional runs, as we have preliminary tested. We have simulated every scenario 1000 times in order to derive the statistics shown in the figures.
As mentioned earlier, the high complexity of the exact solution of the Coverage Problem C, allows us to find optima only for small instances of the problem, with a reduced number of gNBs, ABSs and UEs. Instead, as further discussed in Appendix C.1, OnDrone only requires a few iterations to converge. Thus, for realistically larger deployments, we only show the results obtained with OnDrone, and compare the results to what achieved by Seq based on [9], the RA scheme [10], and a modified OnDrone scheme that neglects inter-ABS interference ("iNeg" in the figures). We have evaluated the coverage performance with denser lattices in OnDrone and observed that the results cannot be significantly improved (since they are already close to optimal, as shown in the comparison with the optimal placements), while imposing more computational complexity. We also evaluate the impact of B ezier routes versus straight flying routes [11]. The flight assignment problem is solved optimally and efficiently by the Hungarian method [38], so we do not comment on its performance. At the end of the section we also provide a summary of the lessons learnt from our performance evaluation study.

Coverage Optimization
Here we numerically evaluate the coverage performance at a precise time instant. We pictorially show the drone footprints and indicate the altitude, computed with the different schemes that target coverage maximization, namely with the optimum (i.e., with its approximation obtained by means of Monte Carlo simulations), with OnDrone, and with Seq. With OnDrone, we clearly see in Fig. 6 that ABSs are positioned very close to the optimal positions, while in Fig. 7 we see that the Seq scheme locates drones in much distinct positions (as well as the RA scheme, not shown here to keep the figure clear). Indeed, Seq makes a greedy decision at each step for a given drone. Hence, Seq finds a good position for such drone knowing the interference incurred by the previously located drones. However, the additional interference from drones located afterwards is ignored, which results on shrunk coverage areas as the final performance. Also, the altitudes provided by Monte Carlo simulations and OnDrone are only slightly different. In general, ABSs avoid locations already covered by gNBs.
Another important aspect to pay attention to is the shape of the drone footprints. Unlike in currently used models [3], 5. The RWP model is one of the most studied and used mobility models to assess mobile networks. It is simple and has wide availability. [21], [22], [27], the area served by a drone is not circular, due to inter-drone interference, which is not considered in the mentioned works. In Fig. 8, we show the average number of covered users in a network with 100 users distributed in the ground surface according to a PPP, and different fleet sizes. Solid bars represent the total amount of covered users, while stripped bars represent users covered by ABSs. The figure shows that OnDrone approximates well the optimal solution (within a mere 1 percent from reaching the optimal) both to total coverage and ABS coverage, while neglecting interference leads to very inaccurate results, getting worse as the fleet size increases (and hence backhaul and inter-drone interference increases). In the figure, we also see that Seq only covers around 80 percent of the optimal coverage, depending on the fleet size. RA is not even able to outperform the coverage results from Seq, and its performance soon decays due to interference issues in the presence of as few as 5 ABSs. This shows that although it is practical, the intuition behind RA is not accurate enough for optimizing coverage.
In Fig. 9 we further compare the coverage achieved when considering the four reference environments of Table 4, for the same PPP case discussed above. The LoS likelihood between ABSs and UEs decreases in denser environments. Thus, drone-cells shrink, and ABSs cover less users. Indeed, the figure shows a factor $6 between the ABS coverage achievable in suburban and high-rise environments, and also that the ground network handles higher total coverage percentages as the environment grows denser, at least for small fleet sizes. Here, we again see the high accuracy of OnDrone in comparison to the optimum searched by means of Monte Carlo simulations, in both total and ABS coverage. As before, Seq and RA provide coverage noticeably lower. We also see that neglecting interference is very counter-productive for ABS coverage in low dense environments, since there is higher LoS probability and links are attenuated easily from interfering signals. The opposite behavior is observed with RA, which provides reasonable results for suburban scenarios (better than Seq) but then its performance decays with denser environments (and RA becomes the worst scheme). Results derived in the Cheese and Capital scenarios are qualitatively similar to those discussed above for the PPP scenario. Thus, we omit the results here.
To give a performance sample of OnDrone for larger fleets, Figs. 10 and 11 show the placement of 8 and 6 ABSs in a dense Cheese and Capital scenario, respectively. Here, drones tend to follow UEs distribution, avoiding empty areas and regions covered by gNBs. Indeed, OnDrone avoids also overlapping drone-cells, thus incurring low inter-ABS interference and being able to cover more users.
In Fig. 12 we further study the impact of the fleet size on coverage, with U ¼ 1000 users, for each scenario. Here, solid lines represent total coverage and dashed lines report the portion of users that would be served by ABSs. We also show coverage performance in absence of drones, labelled as Ground in the figures. First, for the PPP scenario in Fig. 12a, adding more drones increases total and ABS coverage because there are more ABSs that can cover larger areas with limited interference. Each analyzed scheme behaves significantly different. When the fleet size becomes larger than    D ¼ 7, the coverage remains stable or even slightly decays for OnDrone. Here, we notice that the larger the fleet is, the more interference issues appear in both, the backhaul and the access network. OnDrone is able to maintain a stable coverage by reallocating positions to drones to have good backhaul connectivity while providing stable and wide coverage to users, thanks to the design based on extremal-optimization. However, Seq does not have a design that allows a reconfiguration of the aerial network. Hence, it suffers more from interference in the access and backhaul sides. The figure also shows that neglecting interference leads to very poor coverage, as well as with the RA scheme. In both cases, adding drones becomes soon counter-productive. While we have adapted in Seq the scheme proposed in [9] to account for interference, the RA scheme proposed in [10] does not target any specific interference metric in the computation of its repulsion component. This explains why RA behaves poorly with as few as 4 or more ABSs in a dense scenario. Instead, as we have checked numerically, when the interference is neglected-or approximated by a constant value-the coverage apparently never stops increasing with the fleet size. In fact, without interference, having more drones implies covering more non-interfering drone-cells. However, in reality, it happens that there is an exact number of drones that maximizes coverage, depending on the environment.
In Figs. 12b and 12c we show the same type of graph for the Cheese and Capital deployment, respectively. The Cheese case confirms that OnDrone is a good option for irregular deployments, in which the area to cover in the ground surface is smaller than the complete circular surface, so drone positions are more packed. Here, we observe more clearly the effect that increasing the fleet size is not beneficial for coverage purposes, since the surface S is irregular, so that users are more packed and hence it is harder for ABSs to fully avoid interference. Indeed, Seq is not able to avoid interference as good as OnDrone, and suffers larger performance drops with more than 5 drones, as well as the scheme neglecting interference. RA appears to be able to manage interference issues, although with no good coverage results overall. In the Capital case, no matter the adopted scheme, once the optimum number of ABSs is reached, it is hard to add more drones without incurring interference, although in this scenario Ondrone substantially overcomes the rest of schemes. This is due to the fact that users are concentrated in relatively small areas and nearby drones can interfere large masses of users. In any case, OnDrone largely outperforms any other benchmarking schemes also in realistic networks as the one extracted from a dense capital city (Madrid).
So far we analyzed our framework in comparison with significant state-of-the-art proposals, and shown that our approach provides significant gain. We now show, in Figs. 13a, 13b, and 13c, users' coverage obtained with OnDrone when tuning a few key parameters. Specifically, in a dense Cheese deployment, we consider that case in which the guaranteed user data rate varies based on the SINR value g A , the drone height becomes fixed, and the ABS power transmission increases, respectively. We show in any case both total and ABS coverage. In Fig. 13a, we consider fleet sizes of D ¼ 2; 4; 6 drones with different user data rate guarantees. Such SINR values correspond to the MCS (Modulation and Coding Scheme) values marked in the figure, according to [48]. Here, we observe that g A ¼ 10:9 dB is a good election: on the one side, lower g A values provide   higher total and ABS coverage since the QoS requirement is less strict, but the MCS is only QPSK, which renders a considerably lower final user throughput; on the other side, the highest QoS requirements lead to much better MCS and coding rate ("c.r." in the figure), but here the QoS requirement gets too strict so that coverage performance falls notably. In Fig. 13b, we fix the drone height to altitudes between 60 and 600 meters, for fleet sizes of D ¼ 2; 4; 6 drones. Moreover, in the left side of the figure we show a histogram with the total and ABS coverage results when the height is left as an adaptive choice of OnDrone, as intended by the framework proposed in this paper. The results show that, depending on the fleet size, there is an optimal fixed height for the fleet where signal strength is good and the interference impairment remains stable, hence providing the best coverage. However, the difference with respect to the case in which the height is adaptive is around 20 percent for ABS coverage, which supports the idea that flexible non-uniform altitudes are convenient for drone-aided networks. In Fig. 13c, we analyze the impact of transmission power. We have selected three typical values for P d Tx : 10 dBm [43], [44], 25 dBm [16], [52], and 44 dBm [25] (besides, 44 dBm is the power transmission used for gNBs in cellular networks). Here, we clearly see that with 25 and 44 dBm the coverage is significantly increased due to better signal strength from the serving drone. However, fleet size increases, the framework is not able to keep the coverage stable and the interference quickly impairs coverage performance. Conversely, a 10 dBm power transmission allows the framework to combat the interference from multiple ABSs and keep the coverage stable. This shows that lower power transmissions make the framework more resilient. Moreover, since aerial networks are very energy-limited [43], using high power transmission needs to guarantee that the performance is more efficient. Note also that the difference in Watts from 10 to 25 dBm is more than 96 percent, while the corresponding coverage improvement is only 30 percent in the best case.

Continuous Repositioning
Having assessed the basic properties and performance figures of OnDrone for drone placement in Section 6.1, we now consider flight routes. Specifically, we study the performance of OnDrone and RA (specifically proposed for dynamic cases) with repositioning routes computed with either the B ezier Scheme or with the Straight Scheme every T ¼ 60 s, in two practical and realistic scenarios: Cheese and Capital. Irregular topologies like Cheese and Capital are more interesting to study with respect to regular PPP cases because they clearly provide visual fact of the importance of our B ezier Scheme or alike schemes using deflected routes when several regions have really low densities of users.
In order to assess dynamic topologies, we consider a random initial position of drones in the network. In the successive time intervals, the source position of each drone is the last location it was occupying in the previous interval.
In order to compute the B ezier curves used as flight routes, we keep adding anchor points and run iterations of the de Casteljau algorithm in the B ezier Scheme until the segments of the piece-wise curve approximating the B ezier curve are all shorter than 3 m, which corresponds to flight segments of 0.2 s. Hence we can consider the network as practically static in each of such segments, and solve the corresponding coverage problem.
In Figs. 14 and 15 we analyze the performance of continuous ABS repositioning with the B ezier Scheme, in comparison to the Straight Scheme in the dense Cheese and high-rise Capital scenarios, respectively. In both cases, we show an illustration of drone trajectories (subfigure (a)) with the corresponding ABS coverage over time (subfigure (b)) and the average ABS coverage for the scenario obtained with longer and repeated simulation runs (subfigure (c)).
In Fig. 14a we show an instance of a network surface where four empty regions have no user (or few users) demanding for coverage, and 1 drone provides coverage assistance. While a drone following straight lines can be barely used during the full time interval since (i) it traverses part of the empty regions and (ii) it does not avoid ground cells covered by fixed gNBs, the B ezier Scheme deflects drone routes towards regions with denser population, avoiding empty regions and gNB-served users. The resulting gain of a simple example is quantified in Fig. 14b. The figure depicts the count of ABS-covered users over time for 60 s (the count is updated at each segment of the piece-wise curve approximating the B ezier path, every 0.2 s). In contrast with the Straight Scheme, in this example the B ezier Scheme provides a coverage gain from the beginning of the interval, and shows a gain of 26 percent after about 45 seconds. The gain remains high until the drone ends its path. The figure also shows that OnDrone largely outperforms RA. Indeed, B ezier and Straight Schemes used with RA most of the time perform worse than OnDrone with the B ezier Scheme. These coverage performance gaps vary depending on the instance. Hence, to fairly compare each scheme, in Fig. 14c we report the average cumulated ABS coverage with each scheme, calculated over the entire numerical simulation of the scenario and over the tested initial positions. On average, by using B ezier flight paths in combination with OnDrone, at the end of an interval T , ABS coverage increases by a remarkable 18 percent compared to Straight Scheme. Using B ezier paths is beneficial also when using the RA placement algorithm. However, since the coverage guarantees of RA are much weaker than with OnDrone, the achievable coverage is lower.
In Fig. 15a we show a snapshot of drone trajectories optimized with the B ezier Scheme over high-rise environmental conditions in the Capital scenario with 1000 UEs located uniformly at random over the city, plus five masses of 200 UEs. Again, to make the presentation visual and simpler, we use a case with a single drone, although the results for more drones are similar. In the figure, we also report the corresponding flight path computed with the Straight Scheme. We observe that the ABS flies from the bottom towards the right side of the city, where there are more users. Moreover, the B ezier path avoids two gNBs in order not to overlap ABS coverage with gNB coverage. In addition, although not shown to keep the figure clear, the B ezier path meets two masses of 200 users on its way to the final destination, where another mass of 200 users is targeted also by the Straight Scheme. However, the path followed under the Straight Scheme does not avoid coverage overlapping with gNBs and includes regions with very low density of users, hence missing key coverage opportunities. Fig. 15b quantifies the gain due to deflected drone paths for an example over these conditions. The figure shows that the B ezier Scheme considerably increases the cumulated number of users covered since the beginning of the interval, reaching a coverage efficiency increase of 33 percent with respect to the Straight Scheme, and with minimal route adjustments. As done in the dense Cheese scenario, to fairly compare the performance obtained with either B ezier or straight paths, in Fig. 15c we report the average gain in terms of the count of users covered since the beginning of a time interval T . We see that, on average, by using B ezier flight paths in combination with the OnDrone placement algorithm, at the end of an interval T the ABS coverage is increased by a remarkable 47 percent. Using B ezier paths is beneficial also when using the RA placement algorithm. However, since the coverage guarantees of RA are much weaker than with OnDrone, the achievable coverage is much lower.
To show the repositioning operation of multiple drones with OnDrone and the B ezier Scheme, in Fig. 16a we show the trajectories of 4 drones in the high-rise Capital scenario, during an interval of duration 10 T . Here, we capture the main behaviours of our repositioning schemes. The figure reflects backhaul gNB associations by showing multi-color paths, each color corresponding to a different gNB indicated on the map with a marker of the same color. Drones at the bottom and left side start at a suboptimal position. On their path, they deflect over more populated zones and follow  streets that connect the origin and final destinations. Once they get to optimal positions, they remain hovering the zone, adjusting according to changes in the presence of users. The drone starting at the top behaves similarly, although it does not benefit from following any street from source to destination because (i) this street is partially covered by the yellow gNB and (ii) user density is low. The remaining drone has the longest path and traverses several regions. First, it flies upwards to comply with the optimal position of the first stages, but soon deflects its path towards the center of the city where a new optimal position is identified, and remains hovering this central region over the rest of the time. We also see that when this drone is very close to the red gNB, it does not switch its backhaul association because the yellow gNB provides enough connectivity while the drone coming from the top has no alternative but to attach to the red gNB. Pictorially showing correlation between drone trajectories and user movements is not simple with the scenarios presented throughout this paper. However, for a simple case with a dozen of users and one ABS, we refer the reader to the illustrative example described in Appendix C.3.
Finally, in Figs. 16b and 16c we analyze, for the same experiment with four repositioning drones, the ABS coverage count (i.e., the average number of users covered by drones in a one-minute time window) and the average time during which an ABS-covered user receives service within the same one-minute time window. We use the Capital scenario in all possible environmental density cases, with 2000 users in total. Fig. 16b shows that the B ezier Scheme provides a coverage gain around 25 percent with respect to using straight paths. More interestingly, the denser environments present higher coverage gains. While the urban environment provides a gain of 25 percent, the dense environment increases such gain up to a remarkable 33 percent. Therefore, deflected paths are more important in denser scenarios, like in historical districts of old cities and in modern downtowns.
The number of seconds spent under ABS coverage, reported in Fig. 16c, tells the quality of coverage opportunities offered by drones on the move. The figure shows little coverage time differences between the B ezier and Straight Schemes, with some improvements observed with the former. Therefore, we can conclude that the increased coverage count offered by using deflected flight paths is not obtained at the expenses of the time spent by users under coverage. In general, we remark that although the B ezier Scheme needs more time to reach the optimal placement identified with OnDrone, it outperforms significantly approaches based on straight flights for ABSs.

Discussion and Lesson Learnt
Our analysis has shown that optimizing the position of a fleet of drones in a coordinated manner is unfeasible with standard solvers. However, using extremal-optimization techniques, we have seen that it is possible to achieve nearlyoptimal results in polynomial time. Thus, it is doable to run drone repositioning on a minute time scale. Indeed, our results show that with realistic topologies, drones are able to follow mobile masses of people over time, or either react quickly to changes, thanks to OnDrone.
The performance of OnDrone that we have proposed and validated depends on the thickness of the lattice used to reduce the space of options used to seek (near-)optimal drone positions. Denser lattices would result in improved accuracy, although at the expense of computational cost. With as much as 10 drones and a few thousands of users, it is possible to achieve accurate near-optimal results on a subsecond time scale using commodity hardware, as we have done for our numerical evaluation. For larger fleets, more powerful hardware is required, which is not a big problem for base stations and definitely not a problem for cellular back-offices, e.g., in the MEC of future 5G networks [53].
Our numerical results have quantitatively shown the importance of integrating a fleet of drones in a cellular network. The presence of inter-drone interference is key in the optimization of drone positions, so that it cannot be neglected, contrary to what so far used in the literature. The presence of interference makes the optimal number of drones finite, with no advantages coming from deploying unnecessarily dense fleets.
We have also seen that repositioning is a key component of the overall drone-aided network framework. Repositioning requires solving not only for 3-D drone positions, but also finding a flight assignment and deflect flight routes. We have used the Hungarian method to solve the assignment problem optimally and B ezier curves to obtain very efficient and dynamic trajectories that offer coverage opportunities to many users without reducing their time under drone coverage. Such dynamic behavior is key for network surfaces in which some regions cannot host users that can benefit of the presence of ABSs, e.g., in forests or in indoor installations in residential and commercial areas, and also to avoid flying over those ground regions already served by gNBs. Indeed, our dynamic B ezier Scheme presents remarkable results, and we have observed that it enhances a lot coverage experienced in realistic topologies, especially in densely populated cities with dense or high-rise profiles. Where main avenues and landmarks attract users, our B ezier Scheme allows the drones to easily follow masses of users on selected paths and areas.
In general, we have shown that it is feasible to have an autonomous dynamic aerial network that reorganizes itself to optimize network coverage.

CONCLUSIONS
In this paper we have proposed a novel optimization framework for drone-aided cellular networks in terms of user coverage under guaranteed signal quality. The framework is novel because it addresses multiple coordinated drones as well as legacy gNBs and analytically accounts for complex interference expressions caused by drone transmissions under non-negligible and variable LoS probability. Specifically, we have studied the integration of a finite fleet of drones acting as aerial base stations connected to and aiding a ground network of gNBs, from/to which they are able to relay traffic. Since cellular users can move, we have shown that implementing a solution for our framework requires solving three important subproblems: finding optimal positions of drones at a given time instant, map drones onto coverage points, and plan deflected flight routes, to optimize network coverage also while moving the drones.
The coverage problem is a non-linear, non-convex and mixed-integer NP-Complete problem, so it is not possible to handle it with conventional off-the-shelf optimizers. The main issue in our novel formulation lies in the intertwined nature of interference, drone positioning, and LoS probability. Hence, we have proposed OnDrone, an extremaloptimization algorithm that performs near-optimally in loworder polynomial time for various user topologies and under different density of LoS obstructions, from suburban to highrise scenarios. OnDrone outperforms state of the art mechanisms because they do not address the root causes of interference. Interestingly, we have found out that unnecessarily large drone fleets would only have negative impacts on coverage, due to interference. Besides, we have unveiled that optimally mapping drones onto coverage targets is doable in negligible time and, most importantly, we have introduced and evaluated a novel and dynamic scheme to compute intelligent drone trajectories upon repositioning. With our B ezier Scheme, it is possible to deflect flight routes of drones so to dramatically improve coverage performance (up to 47 percent in our experiments) in a variety of scenarios, including in real, dense and high-populated cities.

APPENDIX A PROOF OF NP-COMPLETENESS
Proof. We claim that the coverage problem C is NP-Complete, since the minimum geometric disk cover problem can be reduced, in polynomial time, to a particular case of the Coverage Problem C. Since the MGDC problem is a very well known NP-Complete problem [36], then C is also NP-Complete. We prove such claim in what follows.
Let M be the MGDC problem. From the statement of problem M, we have a set of points U (users in our case) and a fixed radius R > 0. Thus we want to minimize the number of disks with radius R that cover all the points of U. We consider a specific set of instances of the Coverage Problem C in which the interference is negligible, all drones hover at the same elevation and are fast enough to reach any position in the surface in the time interval (i.e., S d ¼ S). In this way, the ground coverage area of each ABS is a circle with a radius that depends on the elevation. Let's use elevation h R at which the ground coverage area has the radius R of the MGDC. Thus, we show that this is a particular set of instances of the coverage problem that cannot be deterministically solved in polynomial time.
Assume, reductio ad absurdum, that we can solve problem C in polynomial time for all its instances. Then, we present in Algorithm 4 an algorithm that solves any instance of the MGDC problem in polynomial time. Specifically, the MGDC problem that finds the minimal number of disks of radius R that cover all points in a set of 2-D coordinates can be solved by a linear search of the minimum number of disks. The search proceeds by adding a disk per iteration. It verifies that a given number of disks covers all points by solving the Coverage Problem C for the special instances described above. If the solution of the Coverage Problem C covers all points, the algorithm stops since we have found the minimum number of disks that cover all points. Note that in Algorithm 4, the iteration counter d (which is also the number of disks to use in one iteration) is bounded by the number of points jUj. Thus, Algorithm 4 has a polynomial complexity, since we are assuming that problem C is solvable in polynomial time, which is a contradiction with the fact that the MGDC problem is NP-Complete.
Since MGDC is NP-Complete [36], the claim follows. t u APPENDIX B B EZIER CURVES FOR DRONE FLIGHT PATHS B ezier curves are an outstanding geometrical tool that deflects the trajectory of an ABS close to clusters of UEs while targeting the optimum positioning. B ezier curves are smooth, endpoint interpolators-i.e., the curve begins and ends in a provided source and destination-and are contained in the convex hull of a set of anchor points. The curve is attracted by anchor points, which in our case are specific users. Below, we mathematically define a B ezier curve.
Definition 1. Given a set of points P ¼fP k g n k¼0 & R m , the resulting B ezier curve is b P ðtÞ ¼ X n k¼0 n k Á P k Á t k Á ð1ÀtÞ nÀk ; t 2 ½0; 1: From the definition, we see that the B ezier curve b P ðtÞ begins at the source b P ðt¼0Þ¼P 0 , and ends at the destination b P ðt ¼ 1Þ ¼P n . The variable range for variable t can be adapted with a linear transformation to modify the speed of the curve. Since we are interested only in the flight path, i.e., the trajectory of the curve, we consider t 2½0; 1.
We need to obtain the offset region of the B ezier curve (a. k.a. stroking the curve). Since B ezier offset curves are not analytically obtainable [54], we use the de Casteljau algorithm [8]. This is the main recursive and numerically stable method to draw B ezier curves from a set of anchor points P, by approximating the curve with small segments. We use B ezier curves for the ground projection of the trajectories, while the elevation of drones varies linearly from source to destination heights, indicated as h src and h dest . The actual 3-D trajectory of a drone with anchor points P is # P ðtÞ ¼ b P ðtÞ; h src þ t Á ðh dest Àh src Þ À Á ; t 2 ½0; 1: Finally, we stroke b P ðtÞ by stroking the segments obtained by the de Casteljau algorithm. Notwithstanding the complexity of the problem formulation expressed in Eq. (8), our near-optimal heuristic, OnDrone, requires few iterations. We report in Fig. 17a typical example of OnDrone iterations for each of the three deployment types studied in the manuscript, with 1000 users and 6 drones. The figure plots coverage versus number of iterations. Coverage increases fast in the initial iterations and then flattens until the OnDrone stopping condition is reached. In each of the deployment scenarios of this case, we do not need more than 14 iterations to converge. With up to 10 drones and 2000 users, we have observed that convergence never required more than a few tens of iterations.

C.2 Guaranteed User Data Rate
The coverage optimization problem addressed in this paper aims to maximize the amount of ground users that are covered, either by gNBs or by ABSs, with a guaranteed service availability that allows for a minimum user data rate. Hence, as soon as we determine parameters as the maximum number of users U g and U d allowed to attach to a gNB g or ABS d, respectively, and the data rates R A min and R B min , there is a user data rate guaranteed for each covered user. In this paper, we mainly use an SINR value of 10.9 dB because it allows to use a 16QAM MCS and a reasonable coding rate of 1=2, so that covered users can enjoy a decent network access service. Indeed, we have evaluated in Fig. 13a the impact of varying the SINR values studied in [48], and observed that higher values let coverage decay considerably. In contrast, lower values do not present much relevant impact on coverage, albeit they guarantee worse data rates. Moreover, the maximum allowed number of users served by a gNB or an ABS is set here to 100 users, which reflects the capacity scale of current 3GPP-compliant cells in terms of active users per BS sector.
In Fig. 18, we show the minimum data rate experienced by covered users. The rate is computed with the Shannon formula, using simulated SINR values for a network of 1000 UEs with up to 10 ABSs. With the selected SINR thresholds (10.9 dB) and maximum number of users per base station (100), and with a bandwidth of 20 MHz, the guaranteed data rate is 0.72 Mbps, which is plotted as an horizontal dotted line in the figure.
On the one side, since in the PPP scenario the distribution of users is uniform, the footprint of gNBs and ABSs is not sufficient to cover the maximum allowed number of users. Hence, resources are split among less users, so that they experience rates well above the guaranteed one. On the other side, scenarios like Cheese and, more clearly, Capital present a non-uniform distribution of users, with groups that form spontaneously. Hence, as previously seen in Fig. 12, coverage is higher than for the PPP scenario. The result is that both gNBs and ABSs tend to cover almost as many users as they are allowed to cover, and minimum rates approach the guaranteed one. In all scenarios, increasing the number of  drones is beneficial because the resulting number of users per base station diminishes, on average. The analysis of other scenarios, thresholds, limits in the number of users per base stations and in the backhaul, and environmental conditions lead to very similar qualitative conclusions, so we do not include those results in the figure.
Overall, here we have shown that the good performance discussed for our coverage optimization schemes, are not obtained in change of poor or, what would be worst, uncontrollably low data rates. In contrast, the guaranteed rate computed based on a fistful of optimization constraints in our problem formulation, results in a realistically close approximation for the lower bound of users performance.

C.3 User Mobility and Drone Trajectories
In order to evaluate the relation between user mobility and drone trajectories, in Fig. 19 we picture a small scenario with 12 UEs and 1 ABS. We depict both the ABS movement and the users movement, during 10 minutes. Here, to further simply presentation, we have took out gNBs.
In the figure, trajectories are marked with positions sampled once per minute, although the actual movement followed by users is more erratic, according to the random way-point model simulated with a one-second time resolution. An X marks the starting point for each user, while a green diamond indicates the initial position of the drone.
If we analyze the trajectories minute by minute, we can observe the evolution of optimum coverage over time. For instance, at the beginning, the drone only covers users 4, 5 and 6, which is optimum at that stage (other users are far and there are no larger groups to cover). However, after one minute, these three users walk away from each other making not possible to cover them altogether, so that the ABS immediately reacts and flies where it can cover other three users. Note that, on its first-minute trajectory, the drone deflects its path towards user 4 to cover it momentarily. The effects caused by using the B ezier Scheme are evident throughout the entire trajectory of the drone. During the following minutes, the drone hovers on the left side of the map, until the users begin to gather in the area on the right part of the map. Hence, the drone quickly reacts to this change and gets repositioned on this region, where it stays for the last few minutes, eventually covering four users.
From what presented in this appendix, the proposed OnDrone heuristic, with the help of B ezier trajectories, is suitable for fast reconfiguration and therefore allows to dynamically adjust network coverage at the same time-scale at which the topology of users evolves. The presence of users plays the role of "attractor" for both optimizing the position of drones and for shaping the trajectories followed while repositioning.