# Deploying a Reliable UAV-Aided Communication Service in Disaster Areas

Authors: Vicente Mayor; Rafael Estepa; Antonio Estepa; German Madinabeitia
Journal: Wireless Communications and Mobile Computing
 (2019)

---

## Abstract

When telecommunication infrastructure is damaged by natural disasters, creating a network that can handle voice channels can be vital for search and rescue missions. Unmanned Aerial Vehicles (UAV) equipped with WiFi access points could be rapidly deployed to provide wireless coverage to ground users. This WiFi access network can in turn be used to provide a reliable communication service to be used in search and rescue missions. We formulate a new problem for UAVs optimal deployment which considers not only WiFi coverage but also the mac sublayer (i.e., quality of service). Our goal is to dispatch the minimum number of UAVs for provisioning a WiFi network that enables reliable VoIP communications in disaster scenarios. Among valid solutions, we choose the one that minimizes energy expenditure at the user’s WiFi interface card in order to extend ground user’s smartphone battery life as much as possible. Solutions are found using well-known heuristics such as K-means clusterization and genetic algorithms. Via numerical results, we show that the IEEE 802.11 standard revision has a decisive impact on the number of UAVs required to cover large areas, and that the user’s average energy expenditure (attributable to communications) can be reduced by limiting the maximum altitude for drones or by increasing the VoIP speech quality.


---

## Body

## 1. Introduction

The use of UAVs in natural disasters has become popular in recent years [1, 2]. During the critical first 72 hours, UAVs can be used for tasks such as situational awareness [3], deploying communication systems [4–6], or search and rescue (SAR) missions [7]. In this work, we aim to dispatch UAVs to deploy a reliable communication system that can be used in SAR missions.There are many aspects to consider in the communication that takes place during SAR missions. In [8] the authors identify the main building blocks (e.g., commanding, surveying, relaying, etc.) and its communication requirements in terms of tolerance to delay, jitter and minimum throughput. Such requirements should be carefully considered when selecting wireless link technologies for either UAV-to-UAV or UAV-to-Infrastructure communication [9]. Similarly, in [10] the authors identify requirements for a flexible, secure, robust, and QoS-aware emergency response communication system, analyzing current wireless technologies (e.g., WiFi, WiMAX, cellular, TETRA, and Satellite) in the light of their capability to carry voice, video, and push-to-talk emergency communications. In the aforementioned works however, the fitness criteria are based on link properties such as throughput, delay or coverage versus the requirements of a single communication. Nevertheless, QoS performance in shared access networks such as WiFi does not only depend on signal coverage but also on the traffic generated/consumed by ground users sharing the medium, and if speech quality degrades, the communication service cannot be used. Unfortunately, this fact is frequently overlooked.In this work we take a first step into the problem of deploying a UAV-aided WiFi network that can support real-time communications in SAR missions such as the one illustrated in Figure1. We consider a reference SAR communication service composed of a QoS-guaranteed bidirectional voice channel for each ground user as well as a broadcast audio channel for announcements. Then, we define a new optimization problem for the deployment of a set of UAVs-mounted Access Points (APs) which consists of finding the minimum number of UAVs (and their position) to provide the aforementioned communication service to a set of known ground targets. Among alternative solutions (i.e., same number of UAVs but in different positions) we choose the one that exhibits the lowest energy expenditure in the network interface card (NIC) of terminals so ground users may extend their battery life. To the best of our knowledge, a similar problem has never been considered before. Note that we assume that the performance of the WiFi access network (see Figure 1) is critical for the feasibility of the QoS-guaranteed service while UAV-to-UAV or UAV-to-Infrastructure does not significantly impair QoS, or its impact can be added to our model as a constant extra delay and/or loss. Therefore, the scope of this paper is restricted to the IEEE 802.11 access network to be deployed.Figure 1
Reliable communication service scenario.Our optimization problem could be directly applied to cover known areas (e.g., meeting points or facilities) when the number of users is known or to provide service to semistatic ground users. Besides, we also analyze the limits of the applicability of this problem when ground users are moving.The contributions and originality of this paper are as follows:(i)
We mathematically formulate a new optimal drone location problem that considers coverage, quality of service, and energy expenditure attributable to the communication.(ii)
We provide a mathematical model to predict speech quality for a set of heterogeneous VoIP traffic sources and for different IEEE 802.11 standard amendments.(iii)
We propose a mathematical model to predict energy expenditure for a set of users and for different IEEE 802.11 standard amendments.(iv)
We provide a preliminary analysis of the applicability of our problem to ground users on the move.The findings of this paper can be useful not only for deploying a communication network, but also as a first step in the field of media coverage (audio), which is an unexplored area of application in disaster management according to [1].


## 2. Related Works

The deployment of standalone communication systems in disaster scenarios is not new [1]. UAVs-aided wireless communication systems have been previously proposed to provide ubiquitous coverage (e.g., off-load of existing Base Stations, or functional replacement of damaged ones [4, 9]), relaying by interconnecting groups of distant users [11], or to disseminate or collect information to/from ground terminals [12, 13]. All previous applications can be useful in disaster scenarios. However, the comparison of existing proposals is difficult as each one sets the focus on different aspects of the communication system (e.g., networking schemes and architectures [10, 14, 15], network formation in UAV-to-UAV communication [16], use of hybrid networks with various wireless link technologies [10, 17]) and, generally, does not specify the traffic generated by users, or else such traffic is not voice or video [15].Most works addressing the optimization of UAV deployment seek to find the minimum number of UAVs required to provide signal coverage (i.e., received signal strength) to a set of ground users at known positions [18, 19] disregarding the traffic generated by ground users. In some cases, the optimization problem seeks to cover all users while minimizing the energy expenditure in the UAV as a result of propulsion [20] or the communication [21].UAV-mounted IEEE 802.11 Access Points (APs) have also been explored in the past [8, 22–25]. A comprehensive summary of the challenges and implications of the various IEEE 802.11 amendments in long-range outdoor WiFi deployments can be found in [25]. Studies, however, only analyze the performance of IEEE 802.11 links in terms of delay, jitter, or throughput achievable by one traffic flow. While results confirm the technical feasibility of deploying an UAV-aided WiFi network for real-time voice communications, they have a common weakness: only one communication flow is considered, disregarding the negative effect of the traffic generated/consumed by WiFi users in the network performance. Moreover, it is well known [26] that there is an upper bound in the number of simultaneous conversations that each AP can take so the speech quality perceived by VoIP users is acceptable. For this reason, deploying a WiFi network usable in SAR missions should consider this restriction in the number and position of UAVs deployed. However, to the best of our knowledge, no previous work has taken this into consideration before.


## 3. Guaranteed QoS VoIP Service and Energy Expenditure in WiFi Networks

Speech quality perceived by users largely depends on end-to-end delay and packet loss experienced by VoIP packets, which is related to the performance of the WiFi network. This relation has been investigated over the last decades in its own research field [26–28] but its findings have been seldom applied to AP positioning. This section elaborates on the implications of WiFi network performance in both speech quality and NIC energy expenditure experienced by ground users in our context.Speech quality can be estimated at the planning stage by using the E-Model [29], which provides a quality score termed R factor (0-100). It is generally agreed that conversations are acceptable when its R factor is above a threshold (Rmin), commonly Rmin=65. In its simplest form, R can be expressed as [30](1)R=94.2-Id-Ie,effwhere, for now, it suffices to know that Id represents all impairments due to delay and Ie,eff is a factor that accounts for the impairments caused by low bit-rate coding and packet loss [31]. Both delay and packet loss are certainly affected by the performance of the IEEE 802.11 access network.Figure2 illustrates a system composed of a set of IEEE 802.11 stations associated with one UAV-mounted AP. Each station is represented by its physical layer and its medium access control (MAC) sublayer. While all stations share a common MAC protocol, each user station auto-configure its Modulation and Coding Schema (MCS) according to the Received Signal Strength Indicator (RSSI) from the AP. An example of different MCSs defined for OFDM modulations (e.g., 802.11a/g) is shown in Table 1. Thus, depending on the received signal power, each user terminal will auto-configure its MCS to the greatest bit-rate possible.Table 1
802.11 sensitivity relations for OFDM modulations with 20 MHz channels.
Modulation Coding rate Data rate Sensitivity (Mb/s) (dBm)  BPSK 1/2 6 -82 BPSK 3/4 9 -81 QPSK 1/2 12 -79 QPSK 3/4 18 -77 16-QAM 1/2 24 -74 16-QAM 3/4 36 -70 64-QAM 2/3 48 -66 64-QAM 3/4 54 -65Figure 2
Traffic in the VoWiFi system.At the application level, this work assumes that each ground user handles two real-time VoIP communication channels: one downlink broadcast audio channel to listen to announces from the rescue team and one bidirectional channel for a regular conversation. We assume that the broadcast traffic is always sent at the lowest bit-rate possible.In this scenario, it is well known that VoIP speech quality is mainly determined by the delay and loss experienced at the AP’s downlink transmission buffer [26], which in turn depends on the speech codec used by the VoIP application and the number of simultaneous calls at a given moment (i.e., traffic load). Some studies quantify this effect and determine the maximum number of conversations that can be simultaneously held so that delay and/or loss in the AP does not exceed a certain threshold. For example, in [26] it was found that a maximum of 15 calls for 64kb/s CBR traffic, or 38 calls for VBR VoIP traffic, can be held for IEEE 802.11b.Thus, UAV placement will be strongly influenced by the maximum number of simultaneous calls that each AP can take so that minimum level of speech quality is met. However, finding this number is a complex task which falls into the WiFi QoS research field. There is a plethora of analytical models of the IEEE 802.11 mac behaviour [32–35]. In most models, the central variable is the probability that an observed station attempts to transmit in a random time slot (τ). However, different models apply different assumptions to derive their analytical expression for τ. In this paper, the following assumptions are made: (i)
Heterogeneous traffic sources (users may have different physical data bit-rate according to their received signal strength).(ii)
Nonsaturated stations (transmission buffer could be empty).(iii)
Noisy channel (packets can be corrupted due to channel noise).Regarding the energy consumption attributable to the wireless card of IEEE 802.11 stations, one could break it down into the time spent on each possible power state (e.g., transmission, reception, idle) as(2)E=ρtx·Ttx+ρrx·Trx+ρidle·Tidlewhere Ttx, Trx, and Tidle stand for the time spent by the network interface on transmission, reception, and idle states and ρtx, ρrx, and ρidle stand for the power consumption of network interface during those states (the WiFi standard offers the possibility of using Power Saving Mode (PSM), which would add a new factor ρsleep·ρsleep, namely, sleep mode. Nevertheless, time traffic such as VoIP rarely uses the original PSM mode because it can introduce extra delays, so in this work it will not be considered.). According to [36], the average NIC consumption attributable to the VoIP application is 0.15W.Clearly, aspects like retransmissions due to noise or collisions, or the MCS used by a station have a great influence inTtx, Trx, and Tidle. For instance, the station’s data bit-rate (which depends on its MCS) will determine the duration of the transmission or reception of a packet. So we can conclude that the distance between drones and users (the greater distance the lower data rate) and the traffic supported by the WiFi network (which again depends on τ) will influence the NIC energy consumption.Then, the following two aspects are distinctive of our UAV placement problem:(i)
Each AP can only cover as many users as a minimum guaranteed speech quality level allows to (i.e.,R>Rmin).(ii)
From solutions with the same number of UAVs, the one that exhibits the lowest average power consumption will be preferred.


## 4. Problem Statement

In the scenario illustrated in Figure1, we want to deploy a set of drones to create a WiFi access network that can be used to provide the service described in Section 1. The following assumptions are made for the sake of tractability:(i)
The position of ground targets (i.e., users confined in a known area) is known.(ii)
Users have a smartphone running an app which handles the VoIP communications described previously using a known codec.(iii)
Channelisation between APs is done in such a manner that interferences are negligible.


### 4.1. Terminology

We discretize the flying zone as illustrated in Figure3. The set P represents the set of 3D coordinates of each edge of the grid (i.e., potential locations of a drone). The following terms and definitions will be used for the remainder of this paper:(i)
Users are denoted by the setU={1,2,…,U} and at known locations given by {wk∣k∈U}, where wk∈R3 represents the 3D coordinates of user k. There are U=12 ground users illustrated in Figure 3.(ii)
Drones are denoted by the setD={1,2,…,D} and at known locations given by the set X={xi∣i∈D,xi∈P}, where xi represents the 3D coordinates of drone i. There are D=3 drones in Figure 3.(iii)
C(i)⊂U represents the set of users associated with the AP installed at drone i∈D. Figure 3 illustrates the sets C(1), C(2), and C(3) with 4, 1, and 7 users, respectively.(iv)
Letδij be 1 if user i∈C(j) and 0 otherwise. Then, C=∑j=1D‍∑i=1U‍δij=∑i=1D‍|C(i)| is the number of ground users associated with the WiFi network (C=12=U in Figure 3).(v)
R(i) represents the speech quality level for the users associated with the AP at drone i∈D.(vi)
Eij represents the average energy consumption per second of the WiFi network card installed at user i’s terminal when it is associated with drone j∈D.(vii)
E-(j) represents the NIC average energy consumption of stations associated with the AP at drone j. Then E-(j)=∑i=1U‍Eijδij/|C(j)|.Figure 3
UAVs placement scenario.


### 4.2. Problem Definition

Our goal is to minimize the number of drones deployed to provide service to ground users and, among solutions with the same number of drones, minimize the NICs average power consumption attributable to VoWiFi. This can be formulated as(3)minXD+1D∑i=1DE-iEmaxsubject  to∑j=1D∑i=1Uδij=URi≥Rmin,∀i∈Dxi∈P,∀i∈DD≤DmaxwhereD represents the number of drones (which is bounded by Dmax), Emax is a constant that represents the maximum NIC energy consumption possible (i.e., a station in the highest consumption state, e.g., transmission), and E-(i) is the average NIC energy consumption of stations associated with drone i. The constant Rmin represents the minimum acceptable speech quality (e.g., 65). Observe that the objective function is composed of a principal part (number of UAVs, which is an integer between 1 and Dmax, and a subordinated one (which is <1), so its integer part always indicates the minimum number of UAVs that meet all constraints.By minimizing the number of drones launched, the deployment cost is reduced. Besides, choosing the configuration that minimizes the average energy expenditure in stations also benefits the SAR mission as it prolongs the battery life of ground users’ terminals.


## 4.1. Terminology

We discretize the flying zone as illustrated in Figure3. The set P represents the set of 3D coordinates of each edge of the grid (i.e., potential locations of a drone). The following terms and definitions will be used for the remainder of this paper:(i)
Users are denoted by the setU={1,2,…,U} and at known locations given by {wk∣k∈U}, where wk∈R3 represents the 3D coordinates of user k. There are U=12 ground users illustrated in Figure 3.(ii)
Drones are denoted by the setD={1,2,…,D} and at known locations given by the set X={xi∣i∈D,xi∈P}, where xi represents the 3D coordinates of drone i. There are D=3 drones in Figure 3.(iii)
C(i)⊂U represents the set of users associated with the AP installed at drone i∈D. Figure 3 illustrates the sets C(1), C(2), and C(3) with 4, 1, and 7 users, respectively.(iv)
Letδij be 1 if user i∈C(j) and 0 otherwise. Then, C=∑j=1D‍∑i=1U‍δij=∑i=1D‍|C(i)| is the number of ground users associated with the WiFi network (C=12=U in Figure 3).(v)
R(i) represents the speech quality level for the users associated with the AP at drone i∈D.(vi)
Eij represents the average energy consumption per second of the WiFi network card installed at user i’s terminal when it is associated with drone j∈D.(vii)
E-(j) represents the NIC average energy consumption of stations associated with the AP at drone j. Then E-(j)=∑i=1U‍Eijδij/|C(j)|.Figure 3
UAVs placement scenario.


## 4.2. Problem Definition

Our goal is to minimize the number of drones deployed to provide service to ground users and, among solutions with the same number of drones, minimize the NICs average power consumption attributable to VoWiFi. This can be formulated as(3)minXD+1D∑i=1DE-iEmaxsubject  to∑j=1D∑i=1Uδij=URi≥Rmin,∀i∈Dxi∈P,∀i∈DD≤DmaxwhereD represents the number of drones (which is bounded by Dmax), Emax is a constant that represents the maximum NIC energy consumption possible (i.e., a station in the highest consumption state, e.g., transmission), and E-(i) is the average NIC energy consumption of stations associated with drone i. The constant Rmin represents the minimum acceptable speech quality (e.g., 65). Observe that the objective function is composed of a principal part (number of UAVs, which is an integer between 1 and Dmax, and a subordinated one (which is <1), so its integer part always indicates the minimum number of UAVs that meet all constraints.By minimizing the number of drones launched, the deployment cost is reduced. Besides, choosing the configuration that minimizes the average energy expenditure in stations also benefits the SAR mission as it prolongs the battery life of ground users’ terminals.


## 5. Solving by Exhaustive Search

The objective function can be evaluated for all edges of the gridP with one drone (D=1), if no solution is found, the number of drones will be increased (D=2) and every possible combination of positions {x1,x2∣x1∈P,x2∈P,x1≠x2} will be evaluated, and so on, until an optimal solution is found or a maximum number of UAVs is reached (Dmax).Algorithm1 indicates the main steps followed. The input of the algorithm is the set of ground users and their location, the set of edges P, Dmax, and the level of speech quality required Rmin. For an incremental number of drones (D), all possible combinations of UAVs locations X are checked. Firstly, the number of users associated with each drone in the WiFi network is evaluated. If the first constraint is met, the QoS constraint is assessed (considering the lowest R among all UAVs). If both constraints are met, a potential solution is found and the objective function is evaluated. Among potential solutions with the same value of D, the one with the lowest objective function is selected. The algorithm ends after finding the optimal solution, or after trying unsuccessfully with Dmax drones. The output is the set of optimal UAVs locations and the value of the objective function. A null location will be returned if no solution is found.Algorithm 1:Exhaustive search algorithm.
Input: U, {wk}, P, Dmax, Rmin Output: L (location), ofmin (obj. function) Initialization: D = 0, L=⌀, o f m i n = D m a x + 1 (1)while (L=⌀orD≤Dmax) do (2)D++; // increase UAVs (3)forx∈Xdo (4)forj=1toDdo ( 5 )  C ( j ) = associate(U,{wk},X); (6)end (7)if∑i=1D‍|C(i)|==Uthen /∗    firstconstaintmet∗/ (8)forj=1toDdo (9)R(j)=QoSEval(C(j),xj); (10)E-(j)=EnergyEval(C(j),xj); (11)end (12)R=mink=1..D{R(k)}; (13)ifR≥Rminthen /∗    candidate solution∗/ (14)of = D + 1D⁢∑i=1DE-(i)Emax; (15)ifof<ofminthen (16)ofmin=of; (17)L=X; (18)end (19)end (20)end (21)end (22)endNext, we elaborate on the functions used for the assessment of signal coverage, speech quality, and power consumption.


### 5.1. Signal Coverage Evaluation: Associate()

A WiFi station requires Signal to Noise Ratio (SNR) and RSSI to be above minimum levels to properly demodulate the signal from an AP. In case of multiple APs (i.e.,D>1), a station associates to the one with greatest RSSI.LetRSSIij be the power received by user i located at wi from the UAV j located at xj. We assume an open area and Line of Sight. Although more complex path loss models could be used (e.g., [37]) we use a free space propagation model:(4)RSSIij=PTX+147,55+Gij-n·10log10wi-xj-20log10f(5)Gij=10log1010Gmax/20·cos2θijwherePTX (dBm) stands for the power delivered by the transmitter antennas; Gij represents the gain of the antenna between user i and UAV j as indicated in (5); n is the propagation exponent (e.g., approximately 3.3 in outdoor); wi-xj is the Euclidean distance between user i and drone j; f (hertz) is the channel frequency. In (5), θij (radians) accounts for the elevation angle between user’s position wi and UAV’s position xj and Gmax is a constant indicating the maximum power.GivenRSSIij, theSignal to Noise Ratio (SNRij)(dB) can be readily obtained by subtracting the receiver’s noise figure (NF) and thermal noise (N) as indicated in (6) and (7). (6)SNRij=RSSIij-NF-N(7)N=-174+10log10CBWwhere CBW is signal bandwidth as specified by the IEEE 802.11 standard in use.For alli∈U and for all j∈D, let us define γij as a Boolean variable that represents whether user i satisfies minimum thresholds SNRij≥SNRmin and RSSIij≥RSSImin with respect to UAV j or not. Then, the functionassociate() returns, for a specific UAV k, the set of users that meet minimum thresholds and whose RSSI with drone k is greater than with any other drone:(8)Ck=i∈U∣γik=1,RSSIik>RSSIij∀j≠k


### 5.2. VoWiFi Service QoS Evaluation: QoS.eval()

This function uses an analytical model to estimate the delay and packet loss experienced by the set of users associated with a specific AP. This is then used as input to the E-Model to assess the speech quality level of the set of users associated with a particular UAV.Remember from Section3 that each user will be executing an app which handles a broadcast VoIP traffic flow (for announcements) and a bidirectional VoIP conversation. Traffic from all users will be aggregated in the AP-mounted UAV to which users are associated which will act as transparent bridge. Each station is autoconfigured with a specific MCS which determines its data bit-rate and bit error rate (BER). This is illustrated in Figure 4.Figure 4
WiFi system in a generic setC(k).All stations are assumed to run the sameDistributed Coordination Function (DCF) of the 802.11 MAC sublayer, which uses CSMA/CA (Carrier Sense Multiple Access with Collision Avoidance) for medium access control. In summary CSMA/CA works as follows: each contending station must sense the medium during a period of time in order to ensure that it is idle before transmission. If the channel is busy, the station waits a random backoff interval before trying again. The backoff process is based on theBinary Exponential algorithm. Time is discretized by defining a time slot duration, and the algorithm picks a random number of time slots between 0 and 2iWo, where Wo accounts for the minimum contention window value, and i increases by one in each failed attempt up to a ceiling m. After a maximum number of retransmissions is reached (M), the packet is discarded.


#### 5.2.1. Delay and Packet Loss Estimation

The set of stations inC(k) and its AP (see Figure 4) constitute a system whose performance has been largely studied in scientific literature [38–41]. We have used a mix of different analytical models from [32–35] suited for the conditions assumed in Section 3: heterogeneous traffic sources, stations are nonsaturated, and noisy channel. For the sake of clarity, we keep this section as simple as possible using a homogeneous notation. However, the reader is encouraged to read more elaborated information in the corresponding references.The probability that an observed station attempts to transmit in a randomly chosen time slot (τ) can be modelled as in [32]:(9)τ=1η11-qr2Wo1-p1-1-rWo-qr1-p(10)η=1-r+r2WoWo+121-1-rWo+Wo+121-qr2qWo1-1-rWo+rp1-q-rq1-p2+p21-q1-pr2Wo1-1-rWo+qr1-p22Wo1-p-p2pm-11-2p+1where η is defined in (10), q is the probability of having at least one packet queued at the transmission buffer after an average MAC service time, p is the probability that a packet suffers any transmission errors, and r is the probability that at least one packet arrives during an idle state. Assuming nonsaturated stations and Poisson packet arrivals to the queue (with rate λ), r and q can be expressed as (11)r=1-e-λET(12)q=1-e-λETEBwhere E[T] represents the expected average slot duration and E[B] is expected average number of backoff slots that a packet waits before transmission. Due to its complexity, a closed-form of E[T] is deduced in Appendix A. E[B] can, however, be expressed as in [35]: (13)EB=Wo21-p1-p-2pm1-2p-2mpM+1Letτ(j), r(j), q(j), E[B(j)] and λ(j) denote τ, r, q, E[B] and λ for the station j in the system under consideration (i.e., ground users and AP such as in Figure 4), where j refers to either one user station j={1,2,..|C(k)|}, or the AP (j=|C(k)|+1=AP). Let p(j) be the probability of packet transmission error p for station j. Then, p(j) can be broken down as (14)pj=1-Pij+FERj-1-Pij·FERjwhere Pi(j) accounts for the probability that thej-th station finds the channel idle and FER(j) stands for the Frame Error Rate due to channel noise. The probability of finding the channel idle can be further expressed as (15)Pij=∏i=1,i≠jS1-τiwhere S stands for the number of stations in the system, including the AP (i.e., |C(k)|+1).Assuming that frames have a constant size ofL bits (L=preamble+header+data), it is possible to obtain the FER of station j as(16)FERj=1-1-BERjLwhere BER(j) represents the bit error rate of station j, which can be readily calculated (see BER equations in [42] and [43] for DSS and OFDM modulations respectively) if one knows user’s modulation (i.e., MCSj).Finally, solving the nonlinear equation system, thepacket loss of station j can be expressed as(17)PLj=1-1-FERjτj∏i=1,i≠jS1-τiλjETSince we assume a very small buffer size, the queueing delay can be neglected and the only delay component will be the channel access delay. Then, the delay of a stationj can be expressed as (18)DELj=EBjETAs justified in [26], the AP (j=AP) is the most saturated station, leading the packet loss and delay in the system (observe that the AP is also part of the system and its MCS and FER change dynamically according to its communication partner. Our approach is to consider average values of its data bit-rate and FER.). Thus, we can take its packet loss and delay as representative of the worst case. Then, the output of this step is (19)PL=PLAP(20)DEL=DELAP


#### 5.2.2. Speech Quality Estimation

As stated in Section3, the E-Model rates the conversation quality R factor, which can be calculated using (1) [44], whose terms were(i)
Ie,eff is the effective impairment equipment parameter, which is a combination between the impairment equipment parameter at zero packet loss (Ie), and a function of Ie that is dependent on packet loss rate and packet loss behaviour. It can be expressed as(21)Ie,eff=Ie+95-IePplPpl/BurstR+BplwhereIe is a codec-dependent constant associated with codec compression degradation (a list of values from ITU-T codecs were presented in ITUT-T Rec. G.113 Appendix I), Ppl represents the packet loss rate, BurstR accounts for the burst ratio (i.e., equals 1 if packet loss if random and greater otherwise), and Bpl represents the codec packet loss robustness, which also has a specific value for each codec (listed in ITU-T Rec. G117 Appendix I).(ii)
Id accounts for all impairments due to delay of communication chain. A widely accepted approximation for Id can be obtained from one-way delay in communication path (d) as follows: (22)Id=0.024d+0.11d-177.3Hd-177.3whereH is the heavy side function (i.e., H(x)=0 for x<0 and H(x)=1 for x>0). This shows that in practical terms, small delays (e.g. <100 ms) can be disregarded.In this paper we will use the G.711 codec (Ie=0), with BurstR=1 and Bpl=25.1. So, assuming that Ppl=PL (from Eq. (19)) and d=DEL (from Eq. (20)) + 20 ms (from the VoIP codec packetization), R for the set of users associated with a drone can be expressed (the 20 ms of packetization delay are included) as (23)R=93.72-95PLPL+25.1-0.024DEL+0.11DEL-157.3HDEL-157.3


### 5.3. Energy at Terminals: EnergyEval()

As stated in Section3 the energy consumption depends on how long the NIC spends on each one of the possible energy states (TX, RX or IDLE) as shown in (2). Considering all disjoint events that can happen during a slot, the energy consumption at the NIC of station j can be further decomposed as(24)EJj=JσjPi︸idleinterval+Jtx,sjPtx,sj+Jtx,ejPtx,ej+Jtx,cjPtx,cj︸tranmissioninterval+Jrx,sjPrx,sj+Jrx,ejPrx,ej+Jrx,cjPrx,cj︸receptioninterval+Jo,sjPo,sj︸otherstationswhere Jtx,s(j), Jtx,e(j), and Jtx,c(j) represent the expected energy consumption during a successful, erroneous, or collided transmission, respectively, Jrx,s(j), Jrx,e(j), and Jrx,c(j) represent the expected energy consumption during a successful, erroneous, or collided reception respectively, and Jo,s(j) represents the expected energy consumption when listening to a packet successfully sent to other station. The terms Ptx,s(j), Ptx,e(j), Ptx,c(j), Prx,s(j), Prx,e(j), Prx,c(j), Po,s(j) indicate the probability of such events during a time slot (E[T]). Finally, Jσ(j) represents the energy spent during an idle interval (i.e., none of the other events happened) and Pi is the probability of that.As an example, let us analyze the first factor of (24): Jσ(j) and Pi. The expected energy consumption during an idle interval of duration σ will be (25)Jσ=ρidle·σand the probability that stations do not transmit in an interval is (26)Pi=∏j=1S1-τjA detailed expression of the remainder factors in (24) is deduced in Appendix B.Given thatE[J(j)] is the expected energy consumption for each slot, it is possible to carry out the average energy consumption per second of station j as (27)Ej=EJjETThen, the average energy consumption per second of all ground user stations that belong toC(k) will be (28)E-=∑i∈CkEiCkwhich is the output of this function.


## 5.1. Signal Coverage Evaluation: Associate()

A WiFi station requires Signal to Noise Ratio (SNR) and RSSI to be above minimum levels to properly demodulate the signal from an AP. In case of multiple APs (i.e.,D>1), a station associates to the one with greatest RSSI.LetRSSIij be the power received by user i located at wi from the UAV j located at xj. We assume an open area and Line of Sight. Although more complex path loss models could be used (e.g., [37]) we use a free space propagation model:(4)RSSIij=PTX+147,55+Gij-n·10log10wi-xj-20log10f(5)Gij=10log1010Gmax/20·cos2θijwherePTX (dBm) stands for the power delivered by the transmitter antennas; Gij represents the gain of the antenna between user i and UAV j as indicated in (5); n is the propagation exponent (e.g., approximately 3.3 in outdoor); wi-xj is the Euclidean distance between user i and drone j; f (hertz) is the channel frequency. In (5), θij (radians) accounts for the elevation angle between user’s position wi and UAV’s position xj and Gmax is a constant indicating the maximum power.GivenRSSIij, theSignal to Noise Ratio (SNRij)(dB) can be readily obtained by subtracting the receiver’s noise figure (NF) and thermal noise (N) as indicated in (6) and (7). (6)SNRij=RSSIij-NF-N(7)N=-174+10log10CBWwhere CBW is signal bandwidth as specified by the IEEE 802.11 standard in use.For alli∈U and for all j∈D, let us define γij as a Boolean variable that represents whether user i satisfies minimum thresholds SNRij≥SNRmin and RSSIij≥RSSImin with respect to UAV j or not. Then, the functionassociate() returns, for a specific UAV k, the set of users that meet minimum thresholds and whose RSSI with drone k is greater than with any other drone:(8)Ck=i∈U∣γik=1,RSSIik>RSSIij∀j≠k


## 5.2. VoWiFi Service QoS Evaluation: QoS.eval()

This function uses an analytical model to estimate the delay and packet loss experienced by the set of users associated with a specific AP. This is then used as input to the E-Model to assess the speech quality level of the set of users associated with a particular UAV.Remember from Section3 that each user will be executing an app which handles a broadcast VoIP traffic flow (for announcements) and a bidirectional VoIP conversation. Traffic from all users will be aggregated in the AP-mounted UAV to which users are associated which will act as transparent bridge. Each station is autoconfigured with a specific MCS which determines its data bit-rate and bit error rate (BER). This is illustrated in Figure 4.Figure 4
WiFi system in a generic setC(k).All stations are assumed to run the sameDistributed Coordination Function (DCF) of the 802.11 MAC sublayer, which uses CSMA/CA (Carrier Sense Multiple Access with Collision Avoidance) for medium access control. In summary CSMA/CA works as follows: each contending station must sense the medium during a period of time in order to ensure that it is idle before transmission. If the channel is busy, the station waits a random backoff interval before trying again. The backoff process is based on theBinary Exponential algorithm. Time is discretized by defining a time slot duration, and the algorithm picks a random number of time slots between 0 and 2iWo, where Wo accounts for the minimum contention window value, and i increases by one in each failed attempt up to a ceiling m. After a maximum number of retransmissions is reached (M), the packet is discarded.


### 5.2.1. Delay and Packet Loss Estimation

The set of stations inC(k) and its AP (see Figure 4) constitute a system whose performance has been largely studied in scientific literature [38–41]. We have used a mix of different analytical models from [32–35] suited for the conditions assumed in Section 3: heterogeneous traffic sources, stations are nonsaturated, and noisy channel. For the sake of clarity, we keep this section as simple as possible using a homogeneous notation. However, the reader is encouraged to read more elaborated information in the corresponding references.The probability that an observed station attempts to transmit in a randomly chosen time slot (τ) can be modelled as in [32]:(9)τ=1η11-qr2Wo1-p1-1-rWo-qr1-p(10)η=1-r+r2WoWo+121-1-rWo+Wo+121-qr2qWo1-1-rWo+rp1-q-rq1-p2+p21-q1-pr2Wo1-1-rWo+qr1-p22Wo1-p-p2pm-11-2p+1where η is defined in (10), q is the probability of having at least one packet queued at the transmission buffer after an average MAC service time, p is the probability that a packet suffers any transmission errors, and r is the probability that at least one packet arrives during an idle state. Assuming nonsaturated stations and Poisson packet arrivals to the queue (with rate λ), r and q can be expressed as (11)r=1-e-λET(12)q=1-e-λETEBwhere E[T] represents the expected average slot duration and E[B] is expected average number of backoff slots that a packet waits before transmission. Due to its complexity, a closed-form of E[T] is deduced in Appendix A. E[B] can, however, be expressed as in [35]: (13)EB=Wo21-p1-p-2pm1-2p-2mpM+1Letτ(j), r(j), q(j), E[B(j)] and λ(j) denote τ, r, q, E[B] and λ for the station j in the system under consideration (i.e., ground users and AP such as in Figure 4), where j refers to either one user station j={1,2,..|C(k)|}, or the AP (j=|C(k)|+1=AP). Let p(j) be the probability of packet transmission error p for station j. Then, p(j) can be broken down as (14)pj=1-Pij+FERj-1-Pij·FERjwhere Pi(j) accounts for the probability that thej-th station finds the channel idle and FER(j) stands for the Frame Error Rate due to channel noise. The probability of finding the channel idle can be further expressed as (15)Pij=∏i=1,i≠jS1-τiwhere S stands for the number of stations in the system, including the AP (i.e., |C(k)|+1).Assuming that frames have a constant size ofL bits (L=preamble+header+data), it is possible to obtain the FER of station j as(16)FERj=1-1-BERjLwhere BER(j) represents the bit error rate of station j, which can be readily calculated (see BER equations in [42] and [43] for DSS and OFDM modulations respectively) if one knows user’s modulation (i.e., MCSj).Finally, solving the nonlinear equation system, thepacket loss of station j can be expressed as(17)PLj=1-1-FERjτj∏i=1,i≠jS1-τiλjETSince we assume a very small buffer size, the queueing delay can be neglected and the only delay component will be the channel access delay. Then, the delay of a stationj can be expressed as (18)DELj=EBjETAs justified in [26], the AP (j=AP) is the most saturated station, leading the packet loss and delay in the system (observe that the AP is also part of the system and its MCS and FER change dynamically according to its communication partner. Our approach is to consider average values of its data bit-rate and FER.). Thus, we can take its packet loss and delay as representative of the worst case. Then, the output of this step is (19)PL=PLAP(20)DEL=DELAP


### 5.2.2. Speech Quality Estimation

As stated in Section3, the E-Model rates the conversation quality R factor, which can be calculated using (1) [44], whose terms were(i)
Ie,eff is the effective impairment equipment parameter, which is a combination between the impairment equipment parameter at zero packet loss (Ie), and a function of Ie that is dependent on packet loss rate and packet loss behaviour. It can be expressed as(21)Ie,eff=Ie+95-IePplPpl/BurstR+BplwhereIe is a codec-dependent constant associated with codec compression degradation (a list of values from ITU-T codecs were presented in ITUT-T Rec. G.113 Appendix I), Ppl represents the packet loss rate, BurstR accounts for the burst ratio (i.e., equals 1 if packet loss if random and greater otherwise), and Bpl represents the codec packet loss robustness, which also has a specific value for each codec (listed in ITU-T Rec. G117 Appendix I).(ii)
Id accounts for all impairments due to delay of communication chain. A widely accepted approximation for Id can be obtained from one-way delay in communication path (d) as follows: (22)Id=0.024d+0.11d-177.3Hd-177.3whereH is the heavy side function (i.e., H(x)=0 for x<0 and H(x)=1 for x>0). This shows that in practical terms, small delays (e.g. <100 ms) can be disregarded.In this paper we will use the G.711 codec (Ie=0), with BurstR=1 and Bpl=25.1. So, assuming that Ppl=PL (from Eq. (19)) and d=DEL (from Eq. (20)) + 20 ms (from the VoIP codec packetization), R for the set of users associated with a drone can be expressed (the 20 ms of packetization delay are included) as (23)R=93.72-95PLPL+25.1-0.024DEL+0.11DEL-157.3HDEL-157.3


## 5.2.1. Delay and Packet Loss Estimation

The set of stations inC(k) and its AP (see Figure 4) constitute a system whose performance has been largely studied in scientific literature [38–41]. We have used a mix of different analytical models from [32–35] suited for the conditions assumed in Section 3: heterogeneous traffic sources, stations are nonsaturated, and noisy channel. For the sake of clarity, we keep this section as simple as possible using a homogeneous notation. However, the reader is encouraged to read more elaborated information in the corresponding references.The probability that an observed station attempts to transmit in a randomly chosen time slot (τ) can be modelled as in [32]:(9)τ=1η11-qr2Wo1-p1-1-rWo-qr1-p(10)η=1-r+r2WoWo+121-1-rWo+Wo+121-qr2qWo1-1-rWo+rp1-q-rq1-p2+p21-q1-pr2Wo1-1-rWo+qr1-p22Wo1-p-p2pm-11-2p+1where η is defined in (10), q is the probability of having at least one packet queued at the transmission buffer after an average MAC service time, p is the probability that a packet suffers any transmission errors, and r is the probability that at least one packet arrives during an idle state. Assuming nonsaturated stations and Poisson packet arrivals to the queue (with rate λ), r and q can be expressed as (11)r=1-e-λET(12)q=1-e-λETEBwhere E[T] represents the expected average slot duration and E[B] is expected average number of backoff slots that a packet waits before transmission. Due to its complexity, a closed-form of E[T] is deduced in Appendix A. E[B] can, however, be expressed as in [35]: (13)EB=Wo21-p1-p-2pm1-2p-2mpM+1Letτ(j), r(j), q(j), E[B(j)] and λ(j) denote τ, r, q, E[B] and λ for the station j in the system under consideration (i.e., ground users and AP such as in Figure 4), where j refers to either one user station j={1,2,..|C(k)|}, or the AP (j=|C(k)|+1=AP). Let p(j) be the probability of packet transmission error p for station j. Then, p(j) can be broken down as (14)pj=1-Pij+FERj-1-Pij·FERjwhere Pi(j) accounts for the probability that thej-th station finds the channel idle and FER(j) stands for the Frame Error Rate due to channel noise. The probability of finding the channel idle can be further expressed as (15)Pij=∏i=1,i≠jS1-τiwhere S stands for the number of stations in the system, including the AP (i.e., |C(k)|+1).Assuming that frames have a constant size ofL bits (L=preamble+header+data), it is possible to obtain the FER of station j as(16)FERj=1-1-BERjLwhere BER(j) represents the bit error rate of station j, which can be readily calculated (see BER equations in [42] and [43] for DSS and OFDM modulations respectively) if one knows user’s modulation (i.e., MCSj).Finally, solving the nonlinear equation system, thepacket loss of station j can be expressed as(17)PLj=1-1-FERjτj∏i=1,i≠jS1-τiλjETSince we assume a very small buffer size, the queueing delay can be neglected and the only delay component will be the channel access delay. Then, the delay of a stationj can be expressed as (18)DELj=EBjETAs justified in [26], the AP (j=AP) is the most saturated station, leading the packet loss and delay in the system (observe that the AP is also part of the system and its MCS and FER change dynamically according to its communication partner. Our approach is to consider average values of its data bit-rate and FER.). Thus, we can take its packet loss and delay as representative of the worst case. Then, the output of this step is (19)PL=PLAP(20)DEL=DELAP


## 5.2.2. Speech Quality Estimation

As stated in Section3, the E-Model rates the conversation quality R factor, which can be calculated using (1) [44], whose terms were(i)
Ie,eff is the effective impairment equipment parameter, which is a combination between the impairment equipment parameter at zero packet loss (Ie), and a function of Ie that is dependent on packet loss rate and packet loss behaviour. It can be expressed as(21)Ie,eff=Ie+95-IePplPpl/BurstR+BplwhereIe is a codec-dependent constant associated with codec compression degradation (a list of values from ITU-T codecs were presented in ITUT-T Rec. G.113 Appendix I), Ppl represents the packet loss rate, BurstR accounts for the burst ratio (i.e., equals 1 if packet loss if random and greater otherwise), and Bpl represents the codec packet loss robustness, which also has a specific value for each codec (listed in ITU-T Rec. G117 Appendix I).(ii)
Id accounts for all impairments due to delay of communication chain. A widely accepted approximation for Id can be obtained from one-way delay in communication path (d) as follows: (22)Id=0.024d+0.11d-177.3Hd-177.3whereH is the heavy side function (i.e., H(x)=0 for x<0 and H(x)=1 for x>0). This shows that in practical terms, small delays (e.g. <100 ms) can be disregarded.In this paper we will use the G.711 codec (Ie=0), with BurstR=1 and Bpl=25.1. So, assuming that Ppl=PL (from Eq. (19)) and d=DEL (from Eq. (20)) + 20 ms (from the VoIP codec packetization), R for the set of users associated with a drone can be expressed (the 20 ms of packetization delay are included) as (23)R=93.72-95PLPL+25.1-0.024DEL+0.11DEL-157.3HDEL-157.3


## 5.3. Energy at Terminals: EnergyEval()

As stated in Section3 the energy consumption depends on how long the NIC spends on each one of the possible energy states (TX, RX or IDLE) as shown in (2). Considering all disjoint events that can happen during a slot, the energy consumption at the NIC of station j can be further decomposed as(24)EJj=JσjPi︸idleinterval+Jtx,sjPtx,sj+Jtx,ejPtx,ej+Jtx,cjPtx,cj︸tranmissioninterval+Jrx,sjPrx,sj+Jrx,ejPrx,ej+Jrx,cjPrx,cj︸receptioninterval+Jo,sjPo,sj︸otherstationswhere Jtx,s(j), Jtx,e(j), and Jtx,c(j) represent the expected energy consumption during a successful, erroneous, or collided transmission, respectively, Jrx,s(j), Jrx,e(j), and Jrx,c(j) represent the expected energy consumption during a successful, erroneous, or collided reception respectively, and Jo,s(j) represents the expected energy consumption when listening to a packet successfully sent to other station. The terms Ptx,s(j), Ptx,e(j), Ptx,c(j), Prx,s(j), Prx,e(j), Prx,c(j), Po,s(j) indicate the probability of such events during a time slot (E[T]). Finally, Jσ(j) represents the energy spent during an idle interval (i.e., none of the other events happened) and Pi is the probability of that.As an example, let us analyze the first factor of (24): Jσ(j) and Pi. The expected energy consumption during an idle interval of duration σ will be (25)Jσ=ρidle·σand the probability that stations do not transmit in an interval is (26)Pi=∏j=1S1-τjA detailed expression of the remainder factors in (24) is deduced in Appendix B.Given thatE[J(j)] is the expected energy consumption for each slot, it is possible to carry out the average energy consumption per second of station j as (27)Ej=EJjETThen, the average energy consumption per second of all ground user stations that belong toC(k) will be (28)E-=∑i∈CkEiCkwhich is the output of this function.


## 6. Example Solution

This section introduces a first example of the results obtained after implementing Algorithm1 in Matlab®. The scenario includes 50 users randomly distributed among an area of 100 m × 100 m. The X-Y step used was 10 m. The parameters used are listed in Table 3. The consumption values for ρtx, ρrx, and ρidle were 2.5, 0.9, and 0.11, respectively. Unless otherwise specified, such parameters are common to all experiments in this paper.Figure5 provides a graphical representation of the solution obtained, showing users’ distribution among UAVs by color. In this case, 4 drones were necessary to service ground users with the standard IEEE 802.11n.Figure 5
Solution obtained.The solution obtained with our exhaustive search algorithm has been used as input scenario to the network simulator ns-3 in order to validate the IEEE 802.11 analytical model provided in Section5.2. The ns-3 modelYansWifiPhy with a transmission buffer size of one packet has been used in simulation. Table 2 summarizes the results obtained, comparing R and the normalized energy expenditure from our analytical models with that of the simulation. Results presented are the average of 30 runs.Table 2
Solution obtained versus simulation.
Drone (i) x i C ( i ) Algorithm / Simulation R ( i ) [1 - 100] E ¯ ( i ) E m a x  1 ( 20,20,20 ) 11 86 / 87 0.10 / 0.09 2 ( 30,70,20 ) 11 85 / 87 0.11 / 0.10 3 ( 80,30,20 ) 15 80 / 83 0.12 / 0.11 4 ( 80,80,20 ) 13 83 / 85 0.11 / 0.10Table 3
Example solution input parameters.
IEEE Standard Scenario Traffic Constraints  Revision 802.11n Users 50 Call length 180 s R S S I m i n -82 dBm GI 800 ns Size 100 m × 100 m VoIP codec G.711 S N R m i n 20 dB Preamble Greenfield X-Y step 10 m On/Off times CBR R m i n 65 Bandwidth 20 MHz Altitude layers { 20, 30} m Packet interval 20 ms Retries (M) 7 Prop. Exponent 3.3 Broadcast channels 1Looking at Table2 one can see acceptable differences in the output of the E-Model, which suggests that the methods and models proposed in Section 5.2 are accurate enough.While exhaustive search always finds the optimal solution, it exhibits poor scalability since assessingD drones requires to evaluate ∑d=1D‍|P|d possibilities (e.g., ~7·107 points in the scenario under consideration). For that reason, exhaustive search is not always feasible and heuristic search methods are commonly used as a general way to find a (semi) optimal solution.


## 7. Heuristic Solutions

Genetic algorithms (GAs) have demonstrated considerable success in providing good solutions to a wide variety of optimization problems, including UAVs deployment [45, 46]. Other metaheuristic methods, including particle swarms optimization, artificial immune system, and simulated annealing, can be used instead or in hybridization with GAs. However, we opt for using GAs for the sake of simplicity since the performance obtained (which is analyzed in the following subsection) is fair enough for the purpose of providing a numerical analysis of our new problem which is our main contribution.We propose a heuristic search method that replaces lines(2)-(21) in Algorithm 1 by a function call (GAsearch) that runs a genetic algorithm to find the optimal solution for a specific number of drones D as shown in Algorithm 2. D is increased until a valid solution is found up to a maximum value of Dmax drones.Algorithm 2:Heuristic search pseudocode.
Input: U, {wk}, P, Dmax, Rmin Output: L (location), ofmin (obj. function) Initialization: D = 1, L=⌀ (1)while (L=⌀orD≤Dmax) do (2)L, of = GAsearch(D, {wk}, Rmin); /∗    call GA∗/ (3)D++; (4)end


### 7.1. Genetic Algorithm:GAsearch()

We have used the Matlab® R2017A Global Optimization Toolbox, whose main operators are defined in [47]. The basic concepts and steps followed are detailed bellow.


#### 7.1.1. Individuals

An individual is defined as a possible solution to the problem. Thus, each individual is a location ofD drones X={x1,x2,…,xD∣xi∈P,i∈D,xi≠xj|∀j≠i}, where xi∈P represents the 3D coordinates of drone i. For example, for D=2, an individual {x1,x2} is composed by a sequence of genes (e.g., {x1,y1,z1,x2,y2,z2}).


#### 7.1.2. Algorithm Steps

The main steps followed by our GA can be summarized as follows:(1)
Aninitial population is generated.(2)
Each individual from the generation isevaluated and ranked by assessing araw fitness score (e.g., objective function of the problem (3)).(3)
Some individuals areselected to be parents according to their position in the ranking.(4)
A new generation is created as follows:
(i)
A 5% of the new individuals are a copy of the top 5% of the previous generation (elite individuals).(ii)
Of the remaining individuals:
(a)
80% is created byCrossover-and-Mutation (CM) combining the genes of two selected parents (crossover) and applying a mutation to these new individuals with a very low probability pmCM.(b)
20% is created by an operation of mutation of some parents (termed self-reproduction and mutation, SRM) with probabilitypmSRM.(5)
Theexit criteria are checked after creating the new generation. The algorithm finishes when the lowest raw fitness score found cannot be improved (i.e., lowered) after 50 consecutive generations by at least 10-3. If this condition is not met, go to step (2). However, if within the first MAXG generations no individual meets problem constraints, the algorithm returns L=⌀ to indicate that no solutions can be found with that number of drones.For example, for a population size of 200 individuals, each generation would be composed of 10 elite individuals (5%), 152 individuals formed by CM, and 38 individuals formed by SRM. Although this differs slightly from the traditional sequential application of operations done by the canonical GA, it improves convergence according to [48, 49].


#### 7.1.3. Initial Population

An initial population ofp individuals is created. Experimentally, we found that a population size of p=200 individuals provides results that cannot be improved in the scenarios tested.According to [50], choosing an initial population which is not completely random, but containing individuals prone to be fit, improves the performance of the GA. For this reason, we use thek-means clusterization method illustrated in Figure 6. In particular, inside an area similar to the grid X-Y dimensions, U points—which represent 2D user locations from the set {wk∣k∈U}—are taken as input to create D clusters so that the mean distance from each point to its cluster centroid is minimized. Once thek-means method is applied, we add to each 2D centroid the average altitude (i.e., vertical axis) of the grid, obtaining a set of 3D centroids. The first individual of the initial population will locate its drones in the grid edges closest to these centroids.Figure 6
Example of 2D spatial regions where drones of initial population are confined (D=4).The other(p-1) individuals of the initial population are created by distributing each individual’s drones randomly among the spatial regions created around each centroid. In particular, we define cubes with a volume equal to the d-th part of the grid volume. Observe that since the centroid of each cube was determined by k-means clusterization, there could be overlapping regions or regions out of the grid bounds such as those illustrated (in 2D) in Figure 6. Note that drones are always located at grid edges.To guarantee convergence to the optimal, after creating the initial population, these spatial regions are no longer taken into consideration for next generations. This implies that genes carried by individuals from next generations could be any edge∈P without any restriction.


#### 7.1.4. Operators

The following operators have been used:(1)
Evaluation and ranking: each individual is assigned a raw fitness score by assessing the problem objective function as detailled in Appendix C. Since the genetic algorithm is called to search for solutions with D drones, individuals that do not meet the problem constraints are assigned a higher score by adding a penalty (a number between 1 and 2, according to the ratio of ground users that meet QoS constraints) to the value of their objective function. Individuals are then sorted according to their raw fitness score. Finally, if an individual is in the n position in the ranking, it is assigned a new scored termed expectation value of 1/n.(2)
ParentsSelection: we use a stochastic uniform selection process among individuals according to their expectation value. In our case, 342 parents (152×2CM+38SRM) are selected for crossover and/or mutation operations. Therefore, individuals in the top positions are chosen multiple times to be parents.(3)
Crossover: the crossover operator combines the characteristics of two parents to create a new individual. We use a uniform crossover operation by generating a random binary vector which determines for each gene of the child whether it comes from one parent or the other.(4)
Mutation: We apply an exchange-type uniform mutation that consists of changing one gene of an individual with a given probability (pmCM for those individuals generated after crossover, or pmSRM for parents) for a random coordinate within the grid edges (we set pmSRM=1/(3·D) and pmCM=1/2·pmSRM).


### 7.2. Heuristics Performance

This section validates the heuristic method proposed above in terms of accuracy and convergence speed. To do so, we propose a series of experiments using the parameters listed in Table4. The number of users has been increased from 10 to 40 in a fixed area of 50 m × 50 m and as a result, the number of required drones has also increased. Each experiment has been repeated 30 times (users’ position is generated randomly on each run). Results shown represent averaged values.Table 4
Heuristics validation input parameters.
IEEE Standard Scenario Traffic Constraints  Revision 802.11n Users 10 , … , 40 Call length 180 s R S S I m i n -82 dBm GI 800 ns Size 50 m × 50 m VoIP codec G.711 S N R m i n 20 dB Preamble Greenfield X-Y step 1 m On/Off times CBR R m i n 65 Bandwidth 20 MHz Altitude layers 20 m Packet interval 20 ms Retries 7 Prop. Exponent 3.3 Broadcast channels 1


#### 7.2.1. Quality of the Solutions

The optimal solution found with our heuristic method is compared with that from exhaustive search and from simply applying the k-means clusterization method (i.e., the first individual in our initial population). Table5 shows for each method the number of drones D, the average energy consumption at the stations as E, and the value of the problem objective function. Results also include the standard deviation between parentheses when greater than 10-2. It can be observed that, while results obtained with our heuristic method are very close to the optimal ones (i.e., provided by exhaustive search), the k-means clusterization method always overestimates the number of required drones.Table 5
Comparison of exhaustive, heuristic, and k-means.
U 10 20 30 40  exhaustive D 1 1 2 2 E 0.24 0.34 0.28 0.33 o f 1.10 1.14 2.12 2.13  genetic D 1 1 2 2 E 0.24 0.34 0.29 0.33 o f 1.10 1.14 2.12 2.13  k-means D 3 (0.62) 4 (0.88) 5 (0.96) 5 (1.07) E 0.17 (0.01) 0.20 (0.02) 0.21 (0.02) 0.22 (0.02) o f 3.47 (0.62) 4.38 (0.87) 5.18 (0.95) 5.49 (1.06)Results have been extended from 50 to 150 ground users for our heuristic method and k-means but not for the exhaustive search one (such extension poses a computational burden too high for exhaustive search method with the grid under consideration. For example, forD=3 and |P|= 2 500 edges, more than 15·109 possible UAV location combinations would have to be evaluated). Figure 7 plots the value of the objective function in our extended comparison. Results suggest that the number of drones obtained with our heuristics tends to grow linearly with the density of users in the scenario under consideration, outperforming the solutions found with k-means in the studied cases. Observe that since the energy term in (3) is normalized and typically small, the objective function value is almost entirely determined by the number of drones (which explains the steps in the plot).Figure 7
Comparing heuristics, exhaustive search, and clusterization.


#### 7.2.2. Convergence Speed and Complexity

The number of generations until our heuristic algorithm finds the solution is plotted in Figure8. Notice that if no individual satisfies the problem constraints during MAXG consecutive generations, the number of drones D is increased and theGAsearch procedure is called again (this value was determined experimentally: initially, we tried MAXG=200 generations in an scenario of 10 000 m with 1 m X-Y steps but, since solutions were always found during the first five generations, we decided to reduce it to 50.). Results show that the best individual is always found within the first five generations in all scenarios tested. This explains that solutions with D drones (see Figures 7 and 8) commonly require evaluating around ~50D generations (see our exit criteria).Figure 8
Generations evaluated until exit criteria are met.The fast convergence speed obtained can be explained as follows:(i)
Both the initial population selection based on applying k-means and the large population size (p=200) create a first generation of individuals prone to be excellent candidates.(ii)
The number of grid edges (|P|) (which depends on the discretization step and the terrain size) also impacts on convergence speed. Observe that the number of grid edges that belong to each cluster (i.e., |P|/D) is not too large when compared to the population size p. In general, scenarios with a greater number of grid edges would need a larger population to attain fast convergence.In terms of computational complexity, the most demanding procedure in our heuristic method is the assessment of the raw fitness score (see Appendix C), whose computation time can be expressed asO(U·D). However, if one considers that D is upper bounded by the constant Dmax and that the number of evaluations of the fitness function is also upper bounded in our heuristic, then the computational complexity of our heuristic method would be O(U) (if the upper limit of D was not bounded by the problem constraint but by the number of ground users U, then the computational complexity would be O(U3). At any rate, solutions can be computed in polynomial time.). Figure 9 shows the execution time obtained in the previous experiments with a computer with 16 GB RAM and CPU Intel Broadwell x86 (8 cores and 4 threads per core) at 2.2 GHz. Results confirm computation time monotonically increases with U in less than polynomial time for our heuristic method. The exhaustive search method, however, shows a stepped rise in computation time when the number of drones in the solution is increased from 1 to 2 at U=30 (see Figure 7).Figure 9
CPU execution time obtained for the methods used.From the previous results it can be concluded that the quality of the solutions found with the proposed heuristic is acceptable and that computational complexity and computation time are also acceptable for practical purposes.


## 7.1. Genetic Algorithm:GAsearch()

We have used the Matlab® R2017A Global Optimization Toolbox, whose main operators are defined in [47]. The basic concepts and steps followed are detailed bellow.


### 7.1.1. Individuals

An individual is defined as a possible solution to the problem. Thus, each individual is a location ofD drones X={x1,x2,…,xD∣xi∈P,i∈D,xi≠xj|∀j≠i}, where xi∈P represents the 3D coordinates of drone i. For example, for D=2, an individual {x1,x2} is composed by a sequence of genes (e.g., {x1,y1,z1,x2,y2,z2}).


### 7.1.2. Algorithm Steps

The main steps followed by our GA can be summarized as follows:(1)
Aninitial population is generated.(2)
Each individual from the generation isevaluated and ranked by assessing araw fitness score (e.g., objective function of the problem (3)).(3)
Some individuals areselected to be parents according to their position in the ranking.(4)
A new generation is created as follows:
(i)
A 5% of the new individuals are a copy of the top 5% of the previous generation (elite individuals).(ii)
Of the remaining individuals:
(a)
80% is created byCrossover-and-Mutation (CM) combining the genes of two selected parents (crossover) and applying a mutation to these new individuals with a very low probability pmCM.(b)
20% is created by an operation of mutation of some parents (termed self-reproduction and mutation, SRM) with probabilitypmSRM.(5)
Theexit criteria are checked after creating the new generation. The algorithm finishes when the lowest raw fitness score found cannot be improved (i.e., lowered) after 50 consecutive generations by at least 10-3. If this condition is not met, go to step (2). However, if within the first MAXG generations no individual meets problem constraints, the algorithm returns L=⌀ to indicate that no solutions can be found with that number of drones.For example, for a population size of 200 individuals, each generation would be composed of 10 elite individuals (5%), 152 individuals formed by CM, and 38 individuals formed by SRM. Although this differs slightly from the traditional sequential application of operations done by the canonical GA, it improves convergence according to [48, 49].


### 7.1.3. Initial Population

An initial population ofp individuals is created. Experimentally, we found that a population size of p=200 individuals provides results that cannot be improved in the scenarios tested.According to [50], choosing an initial population which is not completely random, but containing individuals prone to be fit, improves the performance of the GA. For this reason, we use thek-means clusterization method illustrated in Figure 6. In particular, inside an area similar to the grid X-Y dimensions, U points—which represent 2D user locations from the set {wk∣k∈U}—are taken as input to create D clusters so that the mean distance from each point to its cluster centroid is minimized. Once thek-means method is applied, we add to each 2D centroid the average altitude (i.e., vertical axis) of the grid, obtaining a set of 3D centroids. The first individual of the initial population will locate its drones in the grid edges closest to these centroids.Figure 6
Example of 2D spatial regions where drones of initial population are confined (D=4).The other(p-1) individuals of the initial population are created by distributing each individual’s drones randomly among the spatial regions created around each centroid. In particular, we define cubes with a volume equal to the d-th part of the grid volume. Observe that since the centroid of each cube was determined by k-means clusterization, there could be overlapping regions or regions out of the grid bounds such as those illustrated (in 2D) in Figure 6. Note that drones are always located at grid edges.To guarantee convergence to the optimal, after creating the initial population, these spatial regions are no longer taken into consideration for next generations. This implies that genes carried by individuals from next generations could be any edge∈P without any restriction.


### 7.1.4. Operators

The following operators have been used:(1)
Evaluation and ranking: each individual is assigned a raw fitness score by assessing the problem objective function as detailled in Appendix C. Since the genetic algorithm is called to search for solutions with D drones, individuals that do not meet the problem constraints are assigned a higher score by adding a penalty (a number between 1 and 2, according to the ratio of ground users that meet QoS constraints) to the value of their objective function. Individuals are then sorted according to their raw fitness score. Finally, if an individual is in the n position in the ranking, it is assigned a new scored termed expectation value of 1/n.(2)
ParentsSelection: we use a stochastic uniform selection process among individuals according to their expectation value. In our case, 342 parents (152×2CM+38SRM) are selected for crossover and/or mutation operations. Therefore, individuals in the top positions are chosen multiple times to be parents.(3)
Crossover: the crossover operator combines the characteristics of two parents to create a new individual. We use a uniform crossover operation by generating a random binary vector which determines for each gene of the child whether it comes from one parent or the other.(4)
Mutation: We apply an exchange-type uniform mutation that consists of changing one gene of an individual with a given probability (pmCM for those individuals generated after crossover, or pmSRM for parents) for a random coordinate within the grid edges (we set pmSRM=1/(3·D) and pmCM=1/2·pmSRM).


## 7.1.1. Individuals

An individual is defined as a possible solution to the problem. Thus, each individual is a location ofD drones X={x1,x2,…,xD∣xi∈P,i∈D,xi≠xj|∀j≠i}, where xi∈P represents the 3D coordinates of drone i. For example, for D=2, an individual {x1,x2} is composed by a sequence of genes (e.g., {x1,y1,z1,x2,y2,z2}).


## 7.1.2. Algorithm Steps

The main steps followed by our GA can be summarized as follows:(1)
Aninitial population is generated.(2)
Each individual from the generation isevaluated and ranked by assessing araw fitness score (e.g., objective function of the problem (3)).(3)
Some individuals areselected to be parents according to their position in the ranking.(4)
A new generation is created as follows:
(i)
A 5% of the new individuals are a copy of the top 5% of the previous generation (elite individuals).(ii)
Of the remaining individuals:
(a)
80% is created byCrossover-and-Mutation (CM) combining the genes of two selected parents (crossover) and applying a mutation to these new individuals with a very low probability pmCM.(b)
20% is created by an operation of mutation of some parents (termed self-reproduction and mutation, SRM) with probabilitypmSRM.(5)
Theexit criteria are checked after creating the new generation. The algorithm finishes when the lowest raw fitness score found cannot be improved (i.e., lowered) after 50 consecutive generations by at least 10-3. If this condition is not met, go to step (2). However, if within the first MAXG generations no individual meets problem constraints, the algorithm returns L=⌀ to indicate that no solutions can be found with that number of drones.For example, for a population size of 200 individuals, each generation would be composed of 10 elite individuals (5%), 152 individuals formed by CM, and 38 individuals formed by SRM. Although this differs slightly from the traditional sequential application of operations done by the canonical GA, it improves convergence according to [48, 49].


## 7.1.3. Initial Population

An initial population ofp individuals is created. Experimentally, we found that a population size of p=200 individuals provides results that cannot be improved in the scenarios tested.According to [50], choosing an initial population which is not completely random, but containing individuals prone to be fit, improves the performance of the GA. For this reason, we use thek-means clusterization method illustrated in Figure 6. In particular, inside an area similar to the grid X-Y dimensions, U points—which represent 2D user locations from the set {wk∣k∈U}—are taken as input to create D clusters so that the mean distance from each point to its cluster centroid is minimized. Once thek-means method is applied, we add to each 2D centroid the average altitude (i.e., vertical axis) of the grid, obtaining a set of 3D centroids. The first individual of the initial population will locate its drones in the grid edges closest to these centroids.Figure 6
Example of 2D spatial regions where drones of initial population are confined (D=4).The other(p-1) individuals of the initial population are created by distributing each individual’s drones randomly among the spatial regions created around each centroid. In particular, we define cubes with a volume equal to the d-th part of the grid volume. Observe that since the centroid of each cube was determined by k-means clusterization, there could be overlapping regions or regions out of the grid bounds such as those illustrated (in 2D) in Figure 6. Note that drones are always located at grid edges.To guarantee convergence to the optimal, after creating the initial population, these spatial regions are no longer taken into consideration for next generations. This implies that genes carried by individuals from next generations could be any edge∈P without any restriction.


## 7.1.4. Operators

The following operators have been used:(1)
Evaluation and ranking: each individual is assigned a raw fitness score by assessing the problem objective function as detailled in Appendix C. Since the genetic algorithm is called to search for solutions with D drones, individuals that do not meet the problem constraints are assigned a higher score by adding a penalty (a number between 1 and 2, according to the ratio of ground users that meet QoS constraints) to the value of their objective function. Individuals are then sorted according to their raw fitness score. Finally, if an individual is in the n position in the ranking, it is assigned a new scored termed expectation value of 1/n.(2)
ParentsSelection: we use a stochastic uniform selection process among individuals according to their expectation value. In our case, 342 parents (152×2CM+38SRM) are selected for crossover and/or mutation operations. Therefore, individuals in the top positions are chosen multiple times to be parents.(3)
Crossover: the crossover operator combines the characteristics of two parents to create a new individual. We use a uniform crossover operation by generating a random binary vector which determines for each gene of the child whether it comes from one parent or the other.(4)
Mutation: We apply an exchange-type uniform mutation that consists of changing one gene of an individual with a given probability (pmCM for those individuals generated after crossover, or pmSRM for parents) for a random coordinate within the grid edges (we set pmSRM=1/(3·D) and pmCM=1/2·pmSRM).


## 7.2. Heuristics Performance

This section validates the heuristic method proposed above in terms of accuracy and convergence speed. To do so, we propose a series of experiments using the parameters listed in Table4. The number of users has been increased from 10 to 40 in a fixed area of 50 m × 50 m and as a result, the number of required drones has also increased. Each experiment has been repeated 30 times (users’ position is generated randomly on each run). Results shown represent averaged values.Table 4
Heuristics validation input parameters.
IEEE Standard Scenario Traffic Constraints  Revision 802.11n Users 10 , … , 40 Call length 180 s R S S I m i n -82 dBm GI 800 ns Size 50 m × 50 m VoIP codec G.711 S N R m i n 20 dB Preamble Greenfield X-Y step 1 m On/Off times CBR R m i n 65 Bandwidth 20 MHz Altitude layers 20 m Packet interval 20 ms Retries 7 Prop. Exponent 3.3 Broadcast channels 1


### 7.2.1. Quality of the Solutions

The optimal solution found with our heuristic method is compared with that from exhaustive search and from simply applying the k-means clusterization method (i.e., the first individual in our initial population). Table5 shows for each method the number of drones D, the average energy consumption at the stations as E, and the value of the problem objective function. Results also include the standard deviation between parentheses when greater than 10-2. It can be observed that, while results obtained with our heuristic method are very close to the optimal ones (i.e., provided by exhaustive search), the k-means clusterization method always overestimates the number of required drones.Table 5
Comparison of exhaustive, heuristic, and k-means.
U 10 20 30 40  exhaustive D 1 1 2 2 E 0.24 0.34 0.28 0.33 o f 1.10 1.14 2.12 2.13  genetic D 1 1 2 2 E 0.24 0.34 0.29 0.33 o f 1.10 1.14 2.12 2.13  k-means D 3 (0.62) 4 (0.88) 5 (0.96) 5 (1.07) E 0.17 (0.01) 0.20 (0.02) 0.21 (0.02) 0.22 (0.02) o f 3.47 (0.62) 4.38 (0.87) 5.18 (0.95) 5.49 (1.06)Results have been extended from 50 to 150 ground users for our heuristic method and k-means but not for the exhaustive search one (such extension poses a computational burden too high for exhaustive search method with the grid under consideration. For example, forD=3 and |P|= 2 500 edges, more than 15·109 possible UAV location combinations would have to be evaluated). Figure 7 plots the value of the objective function in our extended comparison. Results suggest that the number of drones obtained with our heuristics tends to grow linearly with the density of users in the scenario under consideration, outperforming the solutions found with k-means in the studied cases. Observe that since the energy term in (3) is normalized and typically small, the objective function value is almost entirely determined by the number of drones (which explains the steps in the plot).Figure 7
Comparing heuristics, exhaustive search, and clusterization.


### 7.2.2. Convergence Speed and Complexity

The number of generations until our heuristic algorithm finds the solution is plotted in Figure8. Notice that if no individual satisfies the problem constraints during MAXG consecutive generations, the number of drones D is increased and theGAsearch procedure is called again (this value was determined experimentally: initially, we tried MAXG=200 generations in an scenario of 10 000 m with 1 m X-Y steps but, since solutions were always found during the first five generations, we decided to reduce it to 50.). Results show that the best individual is always found within the first five generations in all scenarios tested. This explains that solutions with D drones (see Figures 7 and 8) commonly require evaluating around ~50D generations (see our exit criteria).Figure 8
Generations evaluated until exit criteria are met.The fast convergence speed obtained can be explained as follows:(i)
Both the initial population selection based on applying k-means and the large population size (p=200) create a first generation of individuals prone to be excellent candidates.(ii)
The number of grid edges (|P|) (which depends on the discretization step and the terrain size) also impacts on convergence speed. Observe that the number of grid edges that belong to each cluster (i.e., |P|/D) is not too large when compared to the population size p. In general, scenarios with a greater number of grid edges would need a larger population to attain fast convergence.In terms of computational complexity, the most demanding procedure in our heuristic method is the assessment of the raw fitness score (see Appendix C), whose computation time can be expressed asO(U·D). However, if one considers that D is upper bounded by the constant Dmax and that the number of evaluations of the fitness function is also upper bounded in our heuristic, then the computational complexity of our heuristic method would be O(U) (if the upper limit of D was not bounded by the problem constraint but by the number of ground users U, then the computational complexity would be O(U3). At any rate, solutions can be computed in polynomial time.). Figure 9 shows the execution time obtained in the previous experiments with a computer with 16 GB RAM and CPU Intel Broadwell x86 (8 cores and 4 threads per core) at 2.2 GHz. Results confirm computation time monotonically increases with U in less than polynomial time for our heuristic method. The exhaustive search method, however, shows a stepped rise in computation time when the number of drones in the solution is increased from 1 to 2 at U=30 (see Figure 7).Figure 9
CPU execution time obtained for the methods used.From the previous results it can be concluded that the quality of the solutions found with the proposed heuristic is acceptable and that computational complexity and computation time are also acceptable for practical purposes.


## 7.2.1. Quality of the Solutions

The optimal solution found with our heuristic method is compared with that from exhaustive search and from simply applying the k-means clusterization method (i.e., the first individual in our initial population). Table5 shows for each method the number of drones D, the average energy consumption at the stations as E, and the value of the problem objective function. Results also include the standard deviation between parentheses when greater than 10-2. It can be observed that, while results obtained with our heuristic method are very close to the optimal ones (i.e., provided by exhaustive search), the k-means clusterization method always overestimates the number of required drones.Table 5
Comparison of exhaustive, heuristic, and k-means.
U 10 20 30 40  exhaustive D 1 1 2 2 E 0.24 0.34 0.28 0.33 o f 1.10 1.14 2.12 2.13  genetic D 1 1 2 2 E 0.24 0.34 0.29 0.33 o f 1.10 1.14 2.12 2.13  k-means D 3 (0.62) 4 (0.88) 5 (0.96) 5 (1.07) E 0.17 (0.01) 0.20 (0.02) 0.21 (0.02) 0.22 (0.02) o f 3.47 (0.62) 4.38 (0.87) 5.18 (0.95) 5.49 (1.06)Results have been extended from 50 to 150 ground users for our heuristic method and k-means but not for the exhaustive search one (such extension poses a computational burden too high for exhaustive search method with the grid under consideration. For example, forD=3 and |P|= 2 500 edges, more than 15·109 possible UAV location combinations would have to be evaluated). Figure 7 plots the value of the objective function in our extended comparison. Results suggest that the number of drones obtained with our heuristics tends to grow linearly with the density of users in the scenario under consideration, outperforming the solutions found with k-means in the studied cases. Observe that since the energy term in (3) is normalized and typically small, the objective function value is almost entirely determined by the number of drones (which explains the steps in the plot).Figure 7
Comparing heuristics, exhaustive search, and clusterization.


## 7.2.2. Convergence Speed and Complexity

The number of generations until our heuristic algorithm finds the solution is plotted in Figure8. Notice that if no individual satisfies the problem constraints during MAXG consecutive generations, the number of drones D is increased and theGAsearch procedure is called again (this value was determined experimentally: initially, we tried MAXG=200 generations in an scenario of 10 000 m with 1 m X-Y steps but, since solutions were always found during the first five generations, we decided to reduce it to 50.). Results show that the best individual is always found within the first five generations in all scenarios tested. This explains that solutions with D drones (see Figures 7 and 8) commonly require evaluating around ~50D generations (see our exit criteria).Figure 8
Generations evaluated until exit criteria are met.The fast convergence speed obtained can be explained as follows:(i)
Both the initial population selection based on applying k-means and the large population size (p=200) create a first generation of individuals prone to be excellent candidates.(ii)
The number of grid edges (|P|) (which depends on the discretization step and the terrain size) also impacts on convergence speed. Observe that the number of grid edges that belong to each cluster (i.e., |P|/D) is not too large when compared to the population size p. In general, scenarios with a greater number of grid edges would need a larger population to attain fast convergence.In terms of computational complexity, the most demanding procedure in our heuristic method is the assessment of the raw fitness score (see Appendix C), whose computation time can be expressed asO(U·D). However, if one considers that D is upper bounded by the constant Dmax and that the number of evaluations of the fitness function is also upper bounded in our heuristic, then the computational complexity of our heuristic method would be O(U) (if the upper limit of D was not bounded by the problem constraint but by the number of ground users U, then the computational complexity would be O(U3). At any rate, solutions can be computed in polynomial time.). Figure 9 shows the execution time obtained in the previous experiments with a computer with 16 GB RAM and CPU Intel Broadwell x86 (8 cores and 4 threads per core) at 2.2 GHz. Results confirm computation time monotonically increases with U in less than polynomial time for our heuristic method. The exhaustive search method, however, shows a stepped rise in computation time when the number of drones in the solution is increased from 1 to 2 at U=30 (see Figure 7).Figure 9
CPU execution time obtained for the methods used.From the previous results it can be concluded that the quality of the solutions found with the proposed heuristic is acceptable and that computational complexity and computation time are also acceptable for practical purposes.


## 8. Numerical Results

In this section we provide a numerical analysis using the heuristics proposed in Section7. All results shown represent the average value after repeating each scenario 30 times (each time the user disposition was random). We set the focus on three different aspects: the impact of the IEEE 802.11 standard in use in the number of UAVs used, the energy consumption and its dependence of traffic load, and the impact of speech quality and grid bounds constraints in the energy consumption of ground users.


### 8.1. Influence of the IEEE 802.11 Standard Revision

In this section, we compare the performance of the widely adopted IEEE 802.11 g/n/ac standard amendments with the configuration listed in Table6. In particular we would like to find out if certain amendments are more fit than others for the two following scenarios of application:(i)
Low user density scenario: for a fixed number of users (U=100) we change the terrain size from 100 to 10 000 m2, obtaining a range of user sparsity from 1 to 100 m2/user. In this type of scenarios the number of UAVs required is expected to increase mainly due to the signal coverage requirement.(ii)
High user density scenario, (sparsity ≤ 1m2/user): in a small area of 25 m2 we change the number of users from 25 to 100 which yields a user sparsity range from 1 to 0.25 m2/user. In these circumstances the number of drones required is mainly determined by the need to satisfy the Speech Quality constraint (Rmin).Table 6
IEEE 802.11 revisions considered.
Frequency Channel Bandwidth Revision  2.4 GHz 20 MHz 802.11g 20 MHz 802.11n 40 MHz  5 GHz 20 MHz 802.11ac 80 MHzResults are shown in Figures10(a) and 10(b). In the low user density, large area scenario (Figure 10(a)), it can be observed that those standard revisions using the 2.4GHz band (g and n) outperform (for large areas) the revision ac which uses 5GHz as a result of lower path loss. Regarding the small scenario with high-density shown in Figure 10(b), no significant differences can be observed among different revisions.Figure 10
Comparison of various revisions of standard IEEE 802.11.
(a)
Low user density, large area scenario (b)
High user density, small area scenario


### 8.2. NIC Energy Consumption

We study how the area size and user density impact on the average energy consumption of ground users. In our study we define two terrain sizes of 10 m× 10 m and 100 m × 100 m (100 m2 and 10 000 m2 respectively), changing the number of ground users from 10 to 100 in steps of 10. The energy is expected to increase with the number of users associated per drone (i.e., |C(k)|) as a result of an increment in collisions.Results shown in Figure11 represent the average energy consumption of ground users and the UAV altitude from all drones. The following points can be made after examining these results:(i)
In the small size scenario, at some points the number of users increase but the average energy decreases (e.g., between 20 and 30 users in Figure11(b)) which seems anti-intuitive. This is attributable to an increment in the number of UAVs (see Figure 11(a)), which is consistent with our hypothesis since the number of users associated with each AP decreases due to the increment in the number of UAVs.(ii)
Ground users in the smaller size terrain exhibit higher energy expenditure than those in the bigger size up to a number of 60 users. This can also be explained as a result of less drones being deployed and a higher number of users associated with each drone, and therefore more traffic load per drone is generated.(iii)
Both curves get slightly inverted after 70 users. In this case the number of deployed UAVs is the same in both scenarios, so energy consumption in both cases tend to be more similar although a little higher in the bigger terrain. This can be explained looking at the average drone altitude in Figure11(c). The larger area produces UAVs in higher altitude to provide coverage. As a consequence, path loss is higher and more users receive a weaker signal, setting a MCS that sets a lower data bit-rate, which translates into more transmission time.Figure 11
Energy consumption analysis.
(a)
Objective function (b)
Average energy consumption (c)
Average UAV altitude


### 8.3. Relation between Energy and Constraints:Rmin and Grid Bounds

We would like to show how the minimum speech quality (Rmin) and the maximum height of the grid (hmax) used in the problem have a direct impact in the energy consumption of ground users. We perform a different experiment for each constraint:(i)
Limited altitude scenario: we create a large area (10 000m2) with a range ground users from 10 to 100 and limit the maximum altitude of the grid to (hmax=10,25,40 meters). In this scenario it is expected that the low upper bound results into more UAVs in the solutions and, in turn, less users per drone which should translate into less energy consumption.(ii)
Higher speech quality scenario: we create a small area of 100 m2 with a range of ground users from 10 to 100 so that coverage is not demanding. Then we solve our deployment problem for different values of the speech quality constraint (Rmin=70,75,80). It is expected that higher quality turns into more UAVs and consequently, as previously stated, should translate into less energy consumption.The results obtained in both scenarios are shown in Figures12 and 13 respectively. Results suggest that our rationale is confirmed: the greater maximum height is, the less UAVs are launched and as a consequence the energy expenditure in ground user stations tends to decrease. The same logic applies to decreasing the minimum speech quality Rmin. Choosing one or the other constraint to regulate the energy consumption has various implications. Limiting the altitude in the grid bound will be more effective in large areas where users are disperse, while increasing the minimum speech quality will have more impact in scenarios with higher density of ground users. However, note than choosing very strict bounds (i.e., too high quality or too low altitude) can severely overestimate the number of UAVs to be deployed.Figure 12
Experiment: limited altitude scenario.
(a)
Objective function (b)
Average energy consumption (c)
Average UAV altitudeFigure 13
Experiment: higher speech quality scenario.
(a)
Objective function (b)
Average energy consumption (c)
Average UAV altitude


## 8.1. Influence of the IEEE 802.11 Standard Revision

In this section, we compare the performance of the widely adopted IEEE 802.11 g/n/ac standard amendments with the configuration listed in Table6. In particular we would like to find out if certain amendments are more fit than others for the two following scenarios of application:(i)
Low user density scenario: for a fixed number of users (U=100) we change the terrain size from 100 to 10 000 m2, obtaining a range of user sparsity from 1 to 100 m2/user. In this type of scenarios the number of UAVs required is expected to increase mainly due to the signal coverage requirement.(ii)
High user density scenario, (sparsity ≤ 1m2/user): in a small area of 25 m2 we change the number of users from 25 to 100 which yields a user sparsity range from 1 to 0.25 m2/user. In these circumstances the number of drones required is mainly determined by the need to satisfy the Speech Quality constraint (Rmin).Table 6
IEEE 802.11 revisions considered.
Frequency Channel Bandwidth Revision  2.4 GHz 20 MHz 802.11g 20 MHz 802.11n 40 MHz  5 GHz 20 MHz 802.11ac 80 MHzResults are shown in Figures10(a) and 10(b). In the low user density, large area scenario (Figure 10(a)), it can be observed that those standard revisions using the 2.4GHz band (g and n) outperform (for large areas) the revision ac which uses 5GHz as a result of lower path loss. Regarding the small scenario with high-density shown in Figure 10(b), no significant differences can be observed among different revisions.Figure 10
Comparison of various revisions of standard IEEE 802.11.
(a)
Low user density, large area scenario (b)
High user density, small area scenario


## 8.2. NIC Energy Consumption

We study how the area size and user density impact on the average energy consumption of ground users. In our study we define two terrain sizes of 10 m× 10 m and 100 m × 100 m (100 m2 and 10 000 m2 respectively), changing the number of ground users from 10 to 100 in steps of 10. The energy is expected to increase with the number of users associated per drone (i.e., |C(k)|) as a result of an increment in collisions.Results shown in Figure11 represent the average energy consumption of ground users and the UAV altitude from all drones. The following points can be made after examining these results:(i)
In the small size scenario, at some points the number of users increase but the average energy decreases (e.g., between 20 and 30 users in Figure11(b)) which seems anti-intuitive. This is attributable to an increment in the number of UAVs (see Figure 11(a)), which is consistent with our hypothesis since the number of users associated with each AP decreases due to the increment in the number of UAVs.(ii)
Ground users in the smaller size terrain exhibit higher energy expenditure than those in the bigger size up to a number of 60 users. This can also be explained as a result of less drones being deployed and a higher number of users associated with each drone, and therefore more traffic load per drone is generated.(iii)
Both curves get slightly inverted after 70 users. In this case the number of deployed UAVs is the same in both scenarios, so energy consumption in both cases tend to be more similar although a little higher in the bigger terrain. This can be explained looking at the average drone altitude in Figure11(c). The larger area produces UAVs in higher altitude to provide coverage. As a consequence, path loss is higher and more users receive a weaker signal, setting a MCS that sets a lower data bit-rate, which translates into more transmission time.Figure 11
Energy consumption analysis.
(a)
Objective function (b)
Average energy consumption (c)
Average UAV altitude


## 8.3. Relation between Energy and Constraints:Rmin and Grid Bounds

We would like to show how the minimum speech quality (Rmin) and the maximum height of the grid (hmax) used in the problem have a direct impact in the energy consumption of ground users. We perform a different experiment for each constraint:(i)
Limited altitude scenario: we create a large area (10 000m2) with a range ground users from 10 to 100 and limit the maximum altitude of the grid to (hmax=10,25,40 meters). In this scenario it is expected that the low upper bound results into more UAVs in the solutions and, in turn, less users per drone which should translate into less energy consumption.(ii)
Higher speech quality scenario: we create a small area of 100 m2 with a range of ground users from 10 to 100 so that coverage is not demanding. Then we solve our deployment problem for different values of the speech quality constraint (Rmin=70,75,80). It is expected that higher quality turns into more UAVs and consequently, as previously stated, should translate into less energy consumption.The results obtained in both scenarios are shown in Figures12 and 13 respectively. Results suggest that our rationale is confirmed: the greater maximum height is, the less UAVs are launched and as a consequence the energy expenditure in ground user stations tends to decrease. The same logic applies to decreasing the minimum speech quality Rmin. Choosing one or the other constraint to regulate the energy consumption has various implications. Limiting the altitude in the grid bound will be more effective in large areas where users are disperse, while increasing the minimum speech quality will have more impact in scenarios with higher density of ground users. However, note than choosing very strict bounds (i.e., too high quality or too low altitude) can severely overestimate the number of UAVs to be deployed.Figure 12
Experiment: limited altitude scenario.
(a)
Objective function (b)
Average energy consumption (c)
Average UAV altitudeFigure 13
Experiment: higher speech quality scenario.
(a)
Objective function (b)
Average energy consumption (c)
Average UAV altitude


## 9. Discussion of the applicability to SAR Missions and Open Issues

As stated in Section1, this work is a first step into the provision of a WiFi network that enables a reliable communication service for SAR missions. The problem defined in Equation (3) is aimed at solving the initial deployment of UAVs, but it can also be applied to provide service to ground users gathered in specific areas such as meeting points or rescue areas, or in situations when ground users remain semistatic (e.g., users are conducted to remain still by the communication broadcast channel). In this section we provide a brief analysis of its applicability when ground users are moving and identify key issues to be addressed in further research.


### 9.1. Conditions for the Validity of the Solutions

Given a solution obtained at timet1, the conditions for its validity in t2 (being t1<t2) are(i)
Ci(t1)=Ci(t2), ∀i∈D (i.e., wireless clients do not change their AP)(ii)
RSSIij(t1)≤RSSIij(t2), ∀i,∀j∈U (i.e., ground station’s data rate, energy, and speech quality are not negatively affected)When ground users can freely move the previous conditions are unlikely to be satisfied after a while. In such case, a new solution search has to be carried out, arising new issues to be discussed next.


### 9.2. Preliminary Analysis of the Applicability When Users Are Moving

Assuming that ground users can be tracked, their position could be periodically checked in a process such as the one illustrated in Figure14 where T represents the period between consecutive observations (an alternative scheme to Figure 14 is to ignore validity conditions and recalculate new solutions every T seconds.). In such process, the following key points have to be considered:(i)
The observation period (T): it should be set according to how ground users are expected to move (e.g., speed and direction). Note that the difference between the theoretical set of user locations {wk} at the beginning of each observation period and the actual position of ground users at the end of such period will likely grow with T. As such, long values of T will produce less accurate speech quality prediction and higher chances of not satisfying the validity conditions. Although in Figure 14 we assume that T is long enough to relocate UAVs, its strict lower bound simply requires T≥Tsol, where Tsol stands for the time required to compute a new solution. Obviously, the shorter computation delay the better, which in turn depends on the computational resources available as well as the heuristic method used to find solutions (our heuristic method can be sped up in various ways and other alternative heuristics could also be explored).(ii)
Relocating UAVs: relocating drones may result in transient periods where the problem constraints are not met for some users. The actual strategy for UAVs displacement constitutes a problem itself that deserves further analysis as it should minimize not only the relocation delay (e.g., see [51]) but also service disruption to ground users. For example, initial strategies could include minimizing the sum of all UAVs displacements, the longest UAV displacement, or the number of users affected by service outages.(iii)
Problem objective function: after the previous discussion it seems reasonable to define a new problem targeted to the relocation of UAVs during service provision. This new problem could be an extension of the one defined in this paper but including a cost that accounts periods of disrupted communication, which allows one to balance the pros and cons of relocating UAVs. For example, it may be reasonable to keep the last solution as long as the problem constraints are satisfied if the new solution implies a transient period where some users lose their connectivity.Figure 14
Continuous process for the provision of the communication service.For illustrative purposes and as a prospective exercise, we have performed the following experiment: using the scenario and parameters shown in Table4 we have placed users who move according to a correlated random walk model (each ground user moves at a walking speed of 5.3 km/h with a probability of 0.8 and, in such case, the direction remains unchanged with probability 0.8. Users rotate 180 degrees when they reach the area bounds.). Then, we perform an off-line computation (i.e., Tsol=0) of the optimal solution every T seconds. Algorithm 2 is used for the initial deployment and the number of drones deployed D remains unchanged for the remainder of the simulation (30 min). When recalculating new solutions we simply get the best solution possible with D drones found by our GA (even though problem constraints were not satisfied). With every new solution, drones are relocated (considering a speed of 60 Km/h) so that the sum of the distances travelled by all drones is the minimum possible. Every five seconds the position of ground users is sampled and the problem constraints (i.e., coverage and speech quality) are verified. For each user, we define a service disruption event when the speech quality constraint is not satisfied for two consecutive samples.Figure15(a) shows the percentage of the simulation time (i.e., 30 minutes) that a ground user on average experiences a service free of disruption events. As expected, results show that longer observation periods provide worse performance, obtaining uneventful service during more than 90% of the time for T≤60 seconds. In our simulations, the longest distance ever travelled by any drone was always shorter than 20 meters, which results in a maximum relocation delay of 1.2 seconds.Figure 15
Experiment: percentage of time with uneventful communication per average user.
(a)
LaunchingD drones using the optimal solution (b)
Launching optimal solution plus one extra droneA final experiment was done similar to the previous one but launching one extra drone (D+1) in the initial deployment (i.e., one extra iteration in the loop in Algorithm 2). The results obtained are plotted in Figure 15(b) and suggest that launching more drones than strictly necessary (D+1 in our case) can be an effective method to avoid service degradation during transient periods.


## 9.1. Conditions for the Validity of the Solutions

Given a solution obtained at timet1, the conditions for its validity in t2 (being t1<t2) are(i)
Ci(t1)=Ci(t2), ∀i∈D (i.e., wireless clients do not change their AP)(ii)
RSSIij(t1)≤RSSIij(t2), ∀i,∀j∈U (i.e., ground station’s data rate, energy, and speech quality are not negatively affected)When ground users can freely move the previous conditions are unlikely to be satisfied after a while. In such case, a new solution search has to be carried out, arising new issues to be discussed next.


## 9.2. Preliminary Analysis of the Applicability When Users Are Moving

Assuming that ground users can be tracked, their position could be periodically checked in a process such as the one illustrated in Figure14 where T represents the period between consecutive observations (an alternative scheme to Figure 14 is to ignore validity conditions and recalculate new solutions every T seconds.). In such process, the following key points have to be considered:(i)
The observation period (T): it should be set according to how ground users are expected to move (e.g., speed and direction). Note that the difference between the theoretical set of user locations {wk} at the beginning of each observation period and the actual position of ground users at the end of such period will likely grow with T. As such, long values of T will produce less accurate speech quality prediction and higher chances of not satisfying the validity conditions. Although in Figure 14 we assume that T is long enough to relocate UAVs, its strict lower bound simply requires T≥Tsol, where Tsol stands for the time required to compute a new solution. Obviously, the shorter computation delay the better, which in turn depends on the computational resources available as well as the heuristic method used to find solutions (our heuristic method can be sped up in various ways and other alternative heuristics could also be explored).(ii)
Relocating UAVs: relocating drones may result in transient periods where the problem constraints are not met for some users. The actual strategy for UAVs displacement constitutes a problem itself that deserves further analysis as it should minimize not only the relocation delay (e.g., see [51]) but also service disruption to ground users. For example, initial strategies could include minimizing the sum of all UAVs displacements, the longest UAV displacement, or the number of users affected by service outages.(iii)
Problem objective function: after the previous discussion it seems reasonable to define a new problem targeted to the relocation of UAVs during service provision. This new problem could be an extension of the one defined in this paper but including a cost that accounts periods of disrupted communication, which allows one to balance the pros and cons of relocating UAVs. For example, it may be reasonable to keep the last solution as long as the problem constraints are satisfied if the new solution implies a transient period where some users lose their connectivity.Figure 14
Continuous process for the provision of the communication service.For illustrative purposes and as a prospective exercise, we have performed the following experiment: using the scenario and parameters shown in Table4 we have placed users who move according to a correlated random walk model (each ground user moves at a walking speed of 5.3 km/h with a probability of 0.8 and, in such case, the direction remains unchanged with probability 0.8. Users rotate 180 degrees when they reach the area bounds.). Then, we perform an off-line computation (i.e., Tsol=0) of the optimal solution every T seconds. Algorithm 2 is used for the initial deployment and the number of drones deployed D remains unchanged for the remainder of the simulation (30 min). When recalculating new solutions we simply get the best solution possible with D drones found by our GA (even though problem constraints were not satisfied). With every new solution, drones are relocated (considering a speed of 60 Km/h) so that the sum of the distances travelled by all drones is the minimum possible. Every five seconds the position of ground users is sampled and the problem constraints (i.e., coverage and speech quality) are verified. For each user, we define a service disruption event when the speech quality constraint is not satisfied for two consecutive samples.Figure15(a) shows the percentage of the simulation time (i.e., 30 minutes) that a ground user on average experiences a service free of disruption events. As expected, results show that longer observation periods provide worse performance, obtaining uneventful service during more than 90% of the time for T≤60 seconds. In our simulations, the longest distance ever travelled by any drone was always shorter than 20 meters, which results in a maximum relocation delay of 1.2 seconds.Figure 15
Experiment: percentage of time with uneventful communication per average user.
(a)
LaunchingD drones using the optimal solution (b)
Launching optimal solution plus one extra droneA final experiment was done similar to the previous one but launching one extra drone (D+1) in the initial deployment (i.e., one extra iteration in the loop in Algorithm 2). The results obtained are plotted in Figure 15(b) and suggest that launching more drones than strictly necessary (D+1 in our case) can be an effective method to avoid service degradation during transient periods.


## 10. Conclusions and Further Work

In this paper we have addressed the problem of UAVs deployment to provide a reliable emergency communication service usable in SAR missions. We have formulated a new problem suited for the initial deployment of UAVs which can also be applied to provide service to ground users gathered in specific areas such as meeting points or rescue areas or in situations when ground users remain semistatic. We have also analyzed the limits of its applicability to moving users. Our optimization problem minimizes the number of drones required to provide a reliable communication service. Among equivalent solutions we choose the one that minimizes energy expenditure at ground user’s stations due to communication. Optimal solutions are found using well-known metaheuristic such as GA since exhaustive search is not generally feasible due to its high computational complexity.Numerical results suggest large area scenarios with highly disperse users benefit from IEEE standards using 2.4 GHz and the lowest bandwidth possible, as less drones will be needed to provide the service. The average number of ground users associated with each UAV will also have impact in the energy consumption and, therefore, in the ground users battery life. For that reason, problem constraints such as the maximum altitude allowed to UAVs or the minimum speech quality required can be used to increase energy expenditure at ground user stations.A number of open issues linked to users’ mobility have been pointed for further research. In addition to including the cost of relocating drones, other optimization problems can also be studied in the context of SAR missions (e.g., optimizing drones’ energy expenditure, maximizing the number of covered users given a fixed number of drones, etc.). Relocation strategies should also be explored, considering not only the minimization of the deployment delay but also the number of users affected and the level of degradation of the service.



---

*Source: 7521513-2019-04-08.xml*